1 /*============================================================================== 2 | 3 | NAME 4 | 5 | ppArith.cpp 6 | 7 | DESCRIPTION 8 | 9 | Miscellaneous integer multiple precision math routines. 10 | 11 | User manual and technical documentation are described in detail in my web page at 12 | http://seanerikoconnor.freeservers.com/Mathematics/AbstractAlgebra/PrimitivePolynomials/overview.html 13 | 14 | LEGAL 15 | 16 | Primpoly Version 16.1 - A Program for Computing Primitive Polynomials. 17 | Copyright (C) 1999-2021 by Sean Erik O'Connor. All Rights Reserved. 18 | 19 | This program is free software: you can redistribute it and/or modify 20 | it under the terms of the GNU General Public License as published by 21 | the Free Software Foundation, either version 3 of the License, or 22 | (at your option) any later version. 23 | 24 | This program is distributed in the hope that it will be useful, 25 | but WITHOUT ANY WARRANTY; without even the implied warranty of 26 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 27 | GNU General Public License for more details. 28 | 29 | You should have received a copy of the GNU General Public License 30 | along with this program. If not, see http://www.gnu.org/licenses/. 31 | 32 | The author's address is seanerikoconnor!AT!gmail!DOT!com 33 | with the !DOT! replaced by . and the !AT! replaced by @ 34 | 35 ==============================================================================*/ 36 37 /*------------------------------------------------------------------------------ 38 | Include Files | 39 ------------------------------------------------------------------------------*/ 40 41 #include <cstdlib> // abort() 42 #include <iostream> // Basic stream I/O. 43 #include <new> // set_new_handler() 44 #include <cmath> // Basic math functions e.g. sqrt() 45 #include <complex> // Complex data type and operations. 46 #include <fstream> // File stream I/O. 47 #include <sstream> // String stream I/O. 48 #include <vector> // STL vector class. 49 #include <string> // STL string class. 50 #include <algorithm> // Iterators. 51 #include <stdexcept> // Exceptions. 52 #include <cassert> // assert() 53 54 using namespace std ; 55 56 /*------------------------------------------------------------------------------ 57 | PP Include Files | 58 ------------------------------------------------------------------------------*/ 59 60 #include "Primpoly.h" // Global functions. 61 #include "ppArith.h" // Basic arithmetic functions. 62 #include "ppBigInt.h" // Arbitrary precision integer arithmetic. 63 #include "ppOperationCount.h" // OperationCount collection for factoring and poly finding. 64 #include "ppFactor.h" // Prime factorization and Euler Phi. 65 #include "ppPolynomial.h" // Polynomial operations and mod polynomial operations. 66 #include "ppParser.h" // Parsing of polynomials and I/O services. 67 #include "ppUnitTest.h" // Complete unit test. 68 69 /*============================================================================= 70 | 71 | NAME 72 | 73 | ModP::operator() 74 | 75 | DESCRIPTION 76 | 77 | Computes k = n mod p where 0 <= k < p for both positive and negative 78 | arguments n for p >= 1. 79 | Takes a surprising proportion of the total compute time according to the 80 | profiler, so worth optimizing. 81 | 82 | EXAMPLE 83 | 84 | The C language gives -5 % 3 = - ( -(-5) mod 3 ) = -( 5 mod 3 ) = -2. 85 | The result is nonzero, so we add p=3 to give 1. 86 | 87 | C computes -3 % 3 = - ( -(-3) mod 3 ) = 0, which we leave unchanged. 88 | 89 | METHOD 90 | 91 | For n >= 0, the C language defines r = n mod p by the equation 92 | 93 | n = kp + r 0 <= r < p 94 | 95 | but when n < 0, C returns the quantity 96 | 97 | r = - ( (-n) mod p ) 98 | 99 | To put the result into the correct range 0 to p-1, add p to r if 100 | r is non-zero. 101 | 102 | By the way, dear old FORTRAN's MOD function does the same thing. 103 | 104 +============================================================================*/ 105 106 template <typename UIntType, typename SIntType> 107 inline UIntType ModP<UIntType,SIntType>::operator()( SIntType n ) 108 { 109 if (p_ <= 0) 110 { 111 ostringstream os ; 112 os << "ModP::operator() The modulus p = " << p_ << " <= 0" ; 113 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ; 114 } 115 116 if (n >= 0) 117 { 118 // Avoid integer division in this special case. 119 if (p_ == 2) 120 return n - ((n >> 1) << 1) ; 121 else 122 return n % p_ ; 123 } 124 // Both types must be signed if we want n % p to compute correctly for n < 0. 125 else 126 return( (n % static_cast<SIntType>(p_)) + p_ ) ; 127 } 128 129 /*============================================================================= 130 | 131 | NAME 132 | 133 | PowerMod::operator() 134 | 135 | DESCRIPTION 136 | 137 | n 138 | Compute a (modulo p) for generic integer types. 139 | 0 140 | where a >= 0, p >= 2, n >= 0. Other cases including 0 141 | will throw an ArithModPError. 142 | 143 | We have two versions: one is generic for any integer type, and 144 | the other is specialized for ppuint. 145 | 146 | EXAMPLE 147 | 148 | BigInt a = 2 ; 149 | BigInt n = 3 ; 150 | BigInt p = 7 ; 151 | try 152 | { 153 | PowerMod<BigInt> powermod( p ) ; 154 | BigInt pwr = powermod( a, n ) ; 155 | } 156 | catch( ArithModPError e ) 157 | { 158 | cout << "should have gotten pwr = 1\n" ; 159 | } 160 | 161 | METHOD 162 | 163 | Multiplication by repeated squaring. 164 | 165 +============================================================================*/ 166 167 // Generic. 168 template <typename IntType> 169 IntType PowerMod<IntType>::operator()( const IntType & a, const IntType & n ) 170 { 171 IntType a1 = a ; 172 173 // Out of range conditions. 174 if (a < static_cast<IntType>( 0u ) || 175 n < static_cast<IntType>( 0u ) || 176 p_ <= static_cast<IntType>( 1u ) || 177 (a == static_cast<IntType>( 0u ) && n == static_cast<IntType>( 0u ))) 178 { 179 ostringstream os ; 180 os << "PowerMod::operator() One or more parameters out of range: a = " 181 << a << " n = " << n << " p_ = " << p_ ; 182 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ; 183 } 184 185 // Quick return for 0 ^ n, a^0 and a^1. 186 if (a == static_cast<IntType>( 0u )) 187 return static_cast<IntType>( 0u ) ; 188 189 if (n == static_cast<IntType>( 0u )) 190 return static_cast<IntType>( 1u ) ; 191 192 if (n == static_cast<IntType>( 1u )) 193 return a % static_cast<IntType>( p_ ) ; 194 195 int bitNum = n.maxBitNumber() ; // Index to highest bit. 196 197 #ifdef DEBUG_PP_ARITH 198 cout << "initial max bitNum = " << bitNum << endl ; 199 cout << "a = " << a << endl ; 200 #endif 201 202 // Find the index of the leading 1 bit. 203 while (!n.testBit( bitNum )) 204 --bitNum ; 205 206 #ifdef DEBUG_PP_ARITH 207 cout << "after skipping leading 0 bits, bitNum = " << bitNum << endl ; 208 #endif 209 210 if (bitNum == -1) 211 { 212 ostringstream os ; 213 os << "PowerMod::operator() " << "bitNum == -1" ; 214 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ; 215 } 216 217 #ifdef DEBUG_PP_ARITH 218 cout << "\nAfter skipping zero bits, bitNum = " << bitNum << endl ; 219 #endif 220 221 // Exponentiation by repeated squaring. Discard the leading 1 bit. 222 // Thereafter, square for every 0 bit; square and multiply by a for 223 // every 1 bit. 224 while ( --bitNum >= 0 ) 225 { 226 a1 = (a1 * a1) % static_cast<IntType>( p_ ) ; // Square mod p. 227 228 if (n.testBit( bitNum )) 229 a1 = (a1 * a) % static_cast<IntType>( p_ ) ; // Times a mod p. 230 231 #ifdef DEBUG_PP_ARITH 232 cout << "S " ; 233 if (n.testBit( bitNum )) 234 cout << "X " ; 235 cout << "Bit num = " << bitNum << " a1 = " << a1 << endl ; 236 #endif 237 } 238 239 #ifdef DEBUG_PP_ARITH 240 cout << "Out of the loop bitNum = " << bitNum << " a1 = " << a1 << endl ; 241 #endif 242 243 return a1 ; 244 } 245 246 247 248 /*============================================================================= 249 | 250 | NAME 251 | 252 | addMod() 253 | 254 | DESCRIPTION 255 | 256 | Computes (x + y) mod p for generic integer types, up to the maximum 257 | unsigned type. We extend the range to the limit by testing and compensating 258 | for carries internally. 259 | 260 | Example: 261 | a = 4294967290 262 | n = 4294967294 263 | (a + b) % n = 4294967286 264 | 265 +============================================================================*/ 266 267 template <typename IntType> 268 IntType addMod( IntType a, IntType b, IntType n ) 269 { 270 #ifdef DEBUG_PP_ARITH 271 cout << "addMod" << endl ; 272 cout << " sizeof IntType = " << 8 * sizeof( IntType ) << " bits" << endl ; 273 cout << " a = " << a << " b = " << b << " n = " << n << endl ; 274 #endif 275 276 // Assume the positive arguments might not already be in the range [0, n). 277 a %= n ; 278 b %= n ; 279 280 // Add with carry. Carry bit thrown away by the language. 281 IntType c = a + b ; 282 283 #ifdef DEBUG_PP_ARITH 284 cout << " a mod n = " << a << " b mod n = " << b << endl ; 285 cout << " c = a + b (discarding carry bit) = " << c << endl ; 286 #endif 287 288 // Test for a carry. Most computer languages don't give us access to 289 // the carry bit, so we must infer whether a carry occurred during addition. 290 // Logically, 291 // (no carry for a+b) => a+b <= a && a+b <= b 292 // The contrapositive is 293 // a+b < a || a+b < b => carry 294 if (c < a || c < b) 295 { 296 // Lemma 1. If a carry occured, n < c < 2n. 297 // Pf. 298 // m 299 // When a carry occurs, c = a + b > 2 > n. But 0 <= a, b < n. So n < c and a + b < n + n = 2n. 300 // 301 // Lemma 2. If a carry occured, we can simply compute c - n as usual, discarding the carry. 302 // Pf. 303 // Since a carry occurred, we an write 304 // m 305 // c = a + b = 2 + c' 306 // 307 // Let's assume the computer hardware does subtraction in two's complement. 308 // Now use two's complement arithmetic on c' - n and see the results. 309 // 310 // c' + TwosComplement( n ) <discarding carry bit> = 311 // 312 // c' + (~n + 1) <discarding carry bit> = 313 // m m 314 // (c' + (2 - n)) mod 2 = 315 // m m m m m m 316 // (c - 2 + 2 + (2 - n)) mod 2 = (c + (2 - n)) mod 2 = 317 // 318 // c + TwosComplement( n ) <discarding carry bit> = 319 // 320 // c - n done the usual way in computer hardware. 321 // 322 // Th. If a carry occurs, c mod n = c - n. 323 // Pf. See Lemma 1 and 2. 324 // 325 // e.g. Let m = 4 bits and 326 // 327 // a = 1 1 0 1 = 13 328 // b = 1 1 0 1 = 13 329 // n = 1 1 1 0 = 14 330 // 331 // c = 1 1 0 1 0 332 // c' = 1 0 1 0 333 // ~n = 0 0 0 1 + 1 = 334 // 0 0 1 0 335 // c' - n = c' + ~n = 336 // 1 0 1 0 337 // + 0 0 1 0 338 // = 1 1 0 0 = 12 = 26 mod 14 339 340 c -= n ; 341 342 #ifdef DEBUG_PP_ARITH 343 cout << " Carry!" << " c < a = " << (c < a ) << "|| c < b = " << (c < b) << endl ; 344 #endif 345 } 346 else 347 { 348 c %= n ; 349 350 #ifdef DEBUG_PP_ARITH 351 cout << " No carry" << endl ; 352 #endif 353 } 354 355 return c ; 356 } 357 358 359 360 /*============================================================================= 361 | 362 | NAME 363 | 364 | timesTwoMod() 365 | 366 | DESCRIPTION 367 | 368 | Computes (2 * x) mod p for generic integer types, up to the maximum 369 | unsigned type. We extend the range to the limit by testing and compensating 370 | for carries internally. 371 | 372 | Example: 373 | a = 10376293541461622782 374 | n = 18446744073709551610 375 | (2 * a) % n = 2305843009213693954 376 | 377 +============================================================================*/ 378 379 template <typename IntType> 380 IntType timesTwoMod( IntType a, IntType n ) 381 { 382 #ifdef DEBUG_PP_ARITH 383 cout << "timesTwoMod" << endl ; 384 cout << " sizeof IntType = " << 8 * sizeof( IntType ) << " bits" << endl ; 385 cout << " a = " << a << " n = " << n << endl ; 386 #endif 387 388 // Assume the positive argument might not already be in the range [0, n). 389 a %= n ; 390 391 #ifdef DEBUG_PP_ARITH 392 cout << " a mod n = " << a << endl ; 393 #endif 394 395 // High bit mask. 396 IntType mask = ((IntType)1 << (8 * sizeof( IntType ) - 1)) ; 397 398 IntType c = (a << 1) ; // Lose the carry bit. 399 400 #ifdef DEBUG_PP_ARITH 401 cout << " mask = " << hex << mask << dec << endl ; 402 cout << " c = 2 a (losing the carry bit) = " << c << endl ; 403 #endif 404 405 // High bit is set so (2 a) will have a carry. 406 if (mask & a) 407 { 408 // See notes in function addMod() above for why this works correctly. 409 c -= n ; 410 411 #ifdef DEBUG_PP_ARITH 412 cout << " Carry!" << endl ; 413 #endif 414 } 415 else 416 { 417 c %= n ; 418 419 #ifdef DEBUG_PP_ARITH 420 cout << " No carry" << endl ; 421 #endif 422 } 423 424 return c ; 425 } 426 427 428 429 /*============================================================================= 430 | 431 | NAME 432 | 433 | multiplyMod() 434 | 435 | DESCRIPTION 436 | 437 | Computes (x * y) mod p for generic integer types, up to the maximum 438 | unsigned type. We extend the range this far by testing and compensating 439 | for carries internally. 440 | 441 | Multiply a * b mod n 442 | 443 | m-1 m-2 444 | a b = a (b 2 + b 2 + ... + b 2 + b ) 445 | m-1 m-2 1 0 446 | 447 | Expanding using Horner's rule, 448 | 449 | a b = (0 2 + a b ) 2 + a b ) 2 + ... + a b ) 2 + a b 450 | m-1 m-2 1 0 451 | 452 | where b = 0 or 1 453 | i 454 | At each step call specialized functions to do r = 2 a mod n and r = r + a mod n 455 | which will work even when carries occur. 456 | 457 | The process takes O( ln n ) operations. An example, 458 | 459 | a = 4294967290 460 | b = 4294967290 461 | n = 4294967294 462 | (a * b) % n = 16 463 | 464 | Reference: "Comments on 'A Computer Algorithm for Calculating the Product 465 | A B Modulo M'" K. R. Sloane, IEEE Transactions on Computers, Vol C-34, No. 3, 466 | March 1985. 467 | 468 +============================================================================*/ 469 470 template<typename IntType> 471 IntType multiplyMod( const IntType a, const IntType b, const IntType n ) 472 { 473 // High bit mask. 474 const int numBits = 8 * sizeof( IntType ) ; 475 IntType mask = ((IntType)1 << (numBits - 1)) ; 476 477 IntType r { 0 } ; 478 479 for (int i = numBits-1 ; i >= 0 ; --i, mask >>= 1) 480 { 481 r = timesTwoMod( r, n ) ; 482 483 if (mask & b) 484 r = addMod( r, a, n ) ; 485 } 486 487 return r ; 488 } 489 490 491 492 /*============================================================================= 493 | 494 | NAME 495 | 496 | PowerMod::operator() 497 | 498 | DESCRIPTION 499 | 500 | Specialized for ppuint type. 501 | 502 +============================================================================*/ 503 504 template<> 505 ppuint PowerMod<ppuint>::operator()( const ppuint & a, const ppuint & n ) 506 { 507 // mask = 1000 ... 000. That is, all bits of mask are zero except for the 508 // most significant bit of the computer word holding its value. 509 ppuint mask{ static_cast<ppuint>(1) << (8 * sizeof( ppuint ) - 1) } ; 510 int bit_count{ 0 } ; // Number of bits in exponent to the right of the leading bit. 511 ppuint n1{ n } ; 512 ppuint product{ a } ; 513 514 // Out of range conditions. 515 if (p_ <= 1 || (a == 0 && n1 == 0)) 516 { 517 ostringstream os ; 518 os << "PowerMod<ppuint>::operator() parameters are out of range: a = " << a << " n = " << n << " p_ = " << p_ ; 519 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ; 520 } 521 522 // n 0 1 523 // Quick return for 0 , a and a 524 if (a == 0) 525 return 0 ; 526 527 if (n == 0) 528 return 1 ; 529 530 if (n == 1) 531 return a % p_ ; 532 533 #ifdef DEBUG_PP_ARITH 534 cout << "a = " << a << endl ; 535 cout << "n = " << n << endl ; 536 cout << "p = " << p_ << endl ; 537 cout << "n1 (before shifting) = " << n1 << endl ; 538 #endif 539 540 // Advance the leading bit of the exponent up to the word's left hand boundary. 541 // Count how many bits were to the right of the leading bit. 542 while (! (n1 & mask)) 543 { 544 n1 <<= 1 ; 545 ++bit_count ; 546 } 547 548 bit_count = (8 * sizeof( ppuint )) - bit_count ; 549 550 #ifdef DEBUG_PP_ARITH 551 cout << "n1 = " << n1 << endl ; 552 cout << "mask = " << mask << endl ; 553 cout << "bit_count = " << bit_count << endl ; 554 #endif 555 556 // Exponentiation by repeated squaring. Discard the leading 1 bit. 557 // Thereafter, square for every 0 bit; square and multiply by x for every 1 bit. 558 while ( --bit_count > 0 ) 559 { 560 #ifdef DEBUG_PP_ARITH 561 cout << "product (before squaring) = " << product << " n1 = " << n1 << endl ; 562 #endif 563 564 // Expose the next bit. 565 n1 <<= 1 ; 566 567 // Out of range conditions. 568 if (product > BigInt::getBase() || a > BigInt::getBase()) 569 { 570 // Square modulo p. 571 product = multiplyMod( product, product, p_ ) ; 572 573 // Leading bit is 1: multiply by a modulo p. 574 if (n1 & mask) 575 product = multiplyMod( a, product, p_ ) ; 576 } 577 else 578 { 579 // Square modulo p. 580 product = (product * product) % p_ ; 581 582 // Leading bit is 1: multiply by a modulo p. 583 if (n1 & mask) 584 product = (a * product) % p_ ; 585 } 586 587 #ifdef DEBUG_PP_ARITH 588 cout << "S " ; 589 if (n1 & mask) 590 cout << "X " ; 591 cout << "product = " << product << " n1 = " << n1 << endl ; 592 #endif 593 } 594 595 return product ; 596 } 597 598 599 600 /*============================================================================= 601 | 602 | NAME 603 | 604 | IsPrimitiveRoot::operator() 605 | 606 | DESCRIPTION 607 | 608 | Returns true if a is a primitive root of prime p, and false otherwise. 609 | Throws ArithModPError if a < 1 or p < 2 or p is even. 610 | p must be a prime number, but we only test that p is 2 or odd. 611 | 612 | EXAMPLE 613 | 614 | ppuint p{ 7 } ; ppuint a{ 3 } ; 615 | IsPrimitiveRoot isRoot( p ) ; 616 | try 617 | { 618 | bool isPrim = isRoot( a ) ; 619 | } 620 | catch( ArithModPError e ) 621 | { 622 | cout << "out of range" << endl ; 623 | } 624 | 1 2 3 625 | 3 is a primitive root of the prime p = 7 since 3 = 3, 3 = 2, 3 = 6, 626 | 4 5 p-1 6 627 | 3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3 = 3 = 1 (mod 7). 628 | 629 | So 3 is a primitive root of 7 because it has maximal period. a = 2 630 | 3 631 | isn't a primitive root of p = 7 because already 2 = 1 (mod 7). 632 | 633 | METHOD 634 | 635 | From number theory, a is a primitive root of the prime p iff 636 | (p-1)/q 637 | a != 1 (mod p) for all prime divisors q of (p-1). 638 | (p-1) 639 | Don't need to check if a = 1 (mod p): Fermat's little 640 | theorem guarantees it. 641 | 642 +============================================================================*/ 643 644 bool IsPrimitiveRoot::operator()( ppuint a ) 645 { 646 PowerMod<ppuint> powermod( p_ ) ; 647 648 // a = 0 (mod p); Zero can't be a primitive root of p. 649 if (a == 0) 650 return false ; 651 652 653 // Error for out of range inputs, including p 654 // being an even number greater than 2. 655 if (p_ < 2 || a < 1 || (p_ > 2 && (p_ % 2 == 0))) 656 { 657 ostringstream os ; 658 os << "IsPrimitiveRoot::operator() Inputs are out of range: p = " << p_ << " a = " << a ; 659 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ; 660 } 661 662 // Special cases: 663 // 1 is the only primitive root of 2; i.e. 1 generates the unit elements 664 // of GF( 2 ); 2 is the only primitive root of 3, and 2 and 3 are the only 665 // primitive roots of 5. 666 // 667 if ( (p_ == 2 && a == 1) || 668 (p_ == 3 && a == 2) || 669 (p_ == 5 && (a == 2 || a == 3)) || 670 (p_ == 7 && (a == 3 || a == 5)) || 671 (p_ == 11 && (a == 2 || a == 6 || a == 7 || a == 8)) || 672 (p_ == 13 && (a == 2 || a == 6 || a == 7 || a == 11))) 673 { 674 return true ; 675 } 676 677 // Reduce a down modulo p. 678 a = a % p_ ; 679 680 // a = 0 (mod p); Zero can't be a primitive root of p. 681 if (a == 0) 682 return false ; 683 684 // Factor p-1 into distinct primes. 685 Factorization<ppuint> factorization( p_ - 1 ) ; 686 687 // p-1 688 // a 689 // Test --- != 1 (mod p) for all primes q | (p-1). 690 // q 691 // If so, we have a primitive root of p, otherwise, we exit early. 692 // 693 for (unsigned int i = 0 ; i < factorization.numDistinctFactors() ; ++i) 694 { 695 if (powermod( a, (p_ - 1) / factorization.primeFactor( i )) == 1) 696 { 697 return false ; 698 } 699 } 700 701 return true ; 702 } 703 704 705 706 /*============================================================================= 707 | 708 | NAME 709 | 710 | inverse_mod_p 711 | 712 | DESCRIPTION 713 | 714 | -1 715 | Find the inverse of u modulo p. In other words, find u (mod p) 716 | p >= 2. 717 | 718 | EXAMPLE 719 | -1 720 | For p=7, and u = 2, u = 4 because 2 * 4 = 8 (mod 7) = 1. 721 | 722 | METHOD 723 | 724 | Do extended Euclid's algorithm on u and v to find u1, u2, and u3 such that 725 | 726 | u u1 + v u2 = u3 = gcd( u, v ). 727 | 728 | Now let v = p = prime number, so gcd = u3 = 1. Then we get 729 | 730 | u u1 + p ? = 1 731 | 732 | or u u1 = 0 (mod p) which makes u (mod p) the unique multiplicative 733 | inverse of u. 734 | 735 +============================================================================*/ 736 737 ppsint InverseModP::operator()( ppsint u ) 738 { 739 ModP<ppsint,ppsint> mod( p_ ) ; 740 741 // Do Euclid's algorithm to find the quantities. 742 ppsint t1 = 0 ; 743 ppsint t3 = 0 ; 744 ppsint q = 0 ; 745 746 ppsint u1 = 1 ; 747 ppsint u3 = u ; 748 ppsint v1 = 0 ; 749 ppsint v3 = p_ ; 750 751 ppsint inv_v = 0 ; 752 753 while( v3 != 0) 754 { 755 // Casting to do a floor function. 756 q = static_cast<ppsint>(u3 / v3) ; 757 758 t1 = u1 - v1 * q ; 759 t3 = u3 - v3 * q ; 760 761 u1 = v1 ; 762 u3 = v3 ; 763 764 v1 = t1 ; 765 v3 = t3 ; 766 } 767 768 inv_v = mod( u1 ) ; 769 770 // -1 771 // Self check: does u u = 1 (mod p)? 772 // 773 if ( mod( u * inv_v ) != 1) 774 { 775 ostringstream os ; 776 os << "InverseModP::operator(): inverse self check failed: u = " << u << " * u^(-1) = " << inv_v << " != 1" ; 777 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ; 778 } 779 780 return inv_v ; 781 } 782 783 784 /*============================================================================= 785 | 786 | NAME 787 | 788 | constCoeffTest 789 | 790 | DESCRIPTION 791 | 792 | n 793 | Test if a = (-1) a (mod p) where a is the constant coefficient 794 | 0 0 795 | of polynomial f(x) and a is the number from the orderR() test. 796 | Return true if we pass the test, false otherwise. 797 | 798 | EXAMPLE 799 | 11 3 800 | Let p = 5, n = 11, f( x ) = x + 2 x + 1, and a = 4. 801 | Thus a = 1. 802 | 0 803 | Then return true since 804 | 805 | 11 806 | (-1) * 1 (mod 5) = -1 (mod 5) = 4 (mod 5). 807 | 808 | METHOD 809 | n 810 | We test if (a - (-1) a ) mod p = 0. 811 | 0 812 | 813 +============================================================================*/ 814 815 bool ArithModP::constCoeffTest( ppsint a, ppsint a0, int n ) 816 { 817 ppsint constant_coeff = a0 ; 818 819 ModP<ppuint,ppsint> mod( p_ ) ; // Initialize the functionoid. 820 821 if (n % 2 != 0) 822 constant_coeff = -constant_coeff ; // (-1)^n < 0 for odd n. 823 824 return( (mod( a - constant_coeff ) == 0) ? true : false ) ; 825 } 826 827 828 829 /*============================================================================= 830 | 831 | NAME 832 | 833 | constCoeffIsPrimitiveRoot 834 | 835 | DESCRIPTION 836 | 837 | n 838 | Test if (-1) a (mod p) is a primitive root of the prime p where 839 | 0 840 | a is the constant term of the polynomial f(x). 841 | 0 842 | 843 | EXAMPLE 844 | 11 3 845 | Let p = 7, n = 11, f( x ) = x + 2 x + 4. Thus a = 4. 846 | 0 847 | Then return true since 848 | 849 | 11 850 | (-1) * 4 (mod 7) = -4 (mod 7) = 3 (mod 7), and 3 is a 851 | 1 2 3 852 | primitive root of the prime 7 since 3 = 3, 3 = 2, 3 = 6, 853 | 4 5 7-1 854 | 3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3 = 1 (mod 7). 855 | 856 | METHOD 857 | 858 | n 859 | We test if (-1) a mod p is a primitive root mod p. 860 | 0 861 | 862 +============================================================================*/ 863 864 bool ArithModP::constCoeffIsPrimitiveRoot( ppuint a0, int n ) 865 { 866 ppsint constant_coeff = a0 ; 867 868 ModP<ppuint,ppsint> mod( p_ ) ; 869 IsPrimitiveRoot isroot( p_ ) ; 870 871 // (-1)^n < 0 for odd n. 872 if (n % 2 != 0) 873 constant_coeff = -constant_coeff ; 874 875 return isroot( mod( constant_coeff ) ) ; 876 } 877 878 879 880 /*============================================================================= 881 | 882 | NAME 883 | 884 | gcd 885 | 886 | DESCRIPTION 887 | 888 | Euclid's method for computing greatest common divisor. 889 | 890 | METHOD 891 | 892 | Iterate the equation 893 | 894 | gcd( u, v ) = gcd( v, r ) 895 | 896 | where r = | u / v | = u - q v 897 | -- -- 898 | until v is zero, terminating in gcd( u, 0 ) = u. 899 | 900 | Let's use Mathematica to illustrate the method for large integers, 901 | 902 | gcd[ u_, v_] := 903 | (* Compute GCD using Euclid's method. *) 904 | Module[ {u2, v2, r}, u2 = u ; v2 = v ; 905 | Print[ "Computing GCD( u = ", u, " v = ", v , " )" ] ; 906 | While[ v2 != 0, 907 | r = Mod[u2, v2]; u2 = v2; v2 = r; 908 | Print[ "r = ", r, " u2 = ", u2, " v2 = ", v2]]; 909 | Return[ u2 ]] 910 | 911 | gcd[ 779953197883173551166308319545, 1282866356929526866866376009397 ] 912 | 913 | Computing GCD( u = 779953197883173551166308319545 v = 1282866356929526866866376009397 ) 914 | u = 779953197883173551166308319545 v = 1282866356929526866866376009397 915 | r = 779953197883173551166308319545 u2 = 1282866356929526866866376009397 v2 = 779953197883173551166308319545 916 | r = 502913159046353315700067689852 u2 = 779953197883173551166308319545 v2 = 502913159046353315700067689852 917 | r = 277040038836820235466240629693 u2 = 502913159046353315700067689852 v2 = 277040038836820235466240629693 918 | r = 225873120209533080233827060159 u2 = 277040038836820235466240629693 v2 = 225873120209533080233827060159 919 | r = 51166918627287155232413569534 u2 = 225873120209533080233827060159 v2 = 51166918627287155232413569534 920 | r = 21205445700384459304172782023 u2 = 51166918627287155232413569534 v2 = 21205445700384459304172782023 921 | r = 8756027226518236624068005488 u2 = 21205445700384459304172782023 v2 = 8756027226518236624068005488 922 | r = 3693391247347986056036771047 u2 = 8756027226518236624068005488 v2 = 3693391247347986056036771047 923 | r = 1369244731822264511994463394 u2 = 3693391247347986056036771047 v2 = 1369244731822264511994463394 924 | r = 954901783703457032047844259 u2 = 1369244731822264511994463394 v2 = 954901783703457032047844259 925 | r = 414342948118807479946619135 u2 = 954901783703457032047844259 v2 = 414342948118807479946619135 926 | r = 126215887465842072154605989 u2 = 414342948118807479946619135 v2 = 126215887465842072154605989 927 | r = 35695285721281263482801168 u2 = 126215887465842072154605989 v2 = 35695285721281263482801168 928 | r = 19130030301998281706202485 u2 = 35695285721281263482801168 v2 = 19130030301998281706202485 929 | r = 16565255419282981776598683 u2 = 19130030301998281706202485 v2 = 16565255419282981776598683 930 | r = 2564774882715299929603802 u2 = 16565255419282981776598683 v2 = 2564774882715299929603802 931 | r = 1176606122991182198975871 u2 = 2564774882715299929603802 v2 = 1176606122991182198975871 932 | r = 211562636732935531652060 u2 = 1176606122991182198975871 v2 = 211562636732935531652060 933 | r = 118792939326504540715571 u2 = 211562636732935531652060 v2 = 118792939326504540715571 934 | r = 92769697406430990936489 u2 = 118792939326504540715571 v2 = 92769697406430990936489 935 | r = 26023241920073549779082 u2 = 92769697406430990936489 v2 = 26023241920073549779082 936 | r = 14699971646210341599243 u2 = 26023241920073549779082 v2 = 14699971646210341599243 937 | r = 11323270273863208179839 u2 = 14699971646210341599243 v2 = 11323270273863208179839 938 | r = 3376701372347133419404 u2 = 11323270273863208179839 v2 = 3376701372347133419404 939 | r = 1193166156821807921627 u2 = 3376701372347133419404 v2 = 1193166156821807921627 940 | r = 990369058703517576150 u2 = 1193166156821807921627 v2 = 990369058703517576150 941 | r = 202797098118290345477 u2 = 990369058703517576150 v2 = 202797098118290345477 942 | r = 179180666230356194242 u2 = 202797098118290345477 v2 = 179180666230356194242 943 | r = 23616431887934151235 u2 = 179180666230356194242 v2 = 23616431887934151235 944 | r = 13865643014817135597 u2 = 23616431887934151235 v2 = 13865643014817135597 945 | r = 9750788873117015638 u2 = 13865643014817135597 v2 = 9750788873117015638 946 | r = 4114854141700119959 u2 = 9750788873117015638 v2 = 4114854141700119959 947 | r = 1521080589716775720 u2 = 4114854141700119959 v2 = 1521080589716775720 948 | r = 1072692962266568519 u2 = 1521080589716775720 v2 = 1072692962266568519 949 | r = 448387627450207201 u2 = 1072692962266568519 v2 = 448387627450207201 950 | r = 175917707366154117 u2 = 448387627450207201 v2 = 175917707366154117 951 | r = 96552212717898967 u2 = 175917707366154117 v2 = 96552212717898967 952 | r = 79365494648255150 u2 = 96552212717898967 v2 = 79365494648255150 953 | r = 17186718069643817 u2 = 79365494648255150 v2 = 17186718069643817 954 | r = 10618622369679882 u2 = 17186718069643817 v2 = 10618622369679882 955 | r = 6568095699963935 u2 = 10618622369679882 v2 = 6568095699963935 956 | r = 4050526669715947 u2 = 6568095699963935 v2 = 4050526669715947 957 | r = 2517569030247988 u2 = 4050526669715947 v2 = 2517569030247988 958 | r = 1532957639467959 u2 = 2517569030247988 v2 = 1532957639467959 959 | r = 984611390780029 u2 = 1532957639467959 v2 = 984611390780029 960 | r = 548346248687930 u2 = 984611390780029 v2 = 548346248687930 961 | r = 436265142092099 u2 = 548346248687930 v2 = 436265142092099 962 | r = 112081106595831 u2 = 436265142092099 v2 = 112081106595831 963 | r = 100021822304606 u2 = 112081106595831 v2 = 100021822304606 964 | r = 12059284291225 u2 = 100021822304606 v2 = 12059284291225 965 | r = 3547547974806 u2 = 12059284291225 v2 = 3547547974806 966 | r = 1416640366807 u2 = 3547547974806 v2 = 1416640366807 967 | r = 714267241192 u2 = 1416640366807 v2 = 714267241192 968 | r = 702373125615 u2 = 714267241192 v2 = 702373125615 969 | r = 11894115577 u2 = 702373125615 v2 = 11894115577 970 | r = 620306572 u2 = 11894115577 v2 = 620306572 971 | r = 108290709 u2 = 620306572 v2 = 108290709 972 | r = 78853027 u2 = 108290709 v2 = 78853027 973 | r = 29437682 u2 = 78853027 v2 = 29437682 974 | r = 19977663 u2 = 29437682 v2 = 19977663 975 | r = 9460019 u2 = 19977663 v2 = 9460019 976 | r = 1057625 u2 = 9460019 v2 = 1057625 977 | r = 999019 u2 = 1057625 v2 = 999019 978 | r = 58606 u2 = 999019 v2 = 58606 979 | r = 2717 u2 = 58606 v2 = 2717 980 | r = 1549 u2 = 2717 v2 = 1549 981 | r = 1168 u2 = 1549 v2 = 1168 982 | r = 381 u2 = 1168 v2 = 381 983 | r = 25 u2 = 381 v2 = 25 984 | r = 6 u2 = 25 v2 = 6 985 | r = 1 u2 = 6 v2 = 1 986 | r = 0 u2 = 1 v2 = 0 987 | ---> 1 988 | 989 +============================================================================*/ 990 991 template <typename IntType> 992 IntType gcd( const IntType & u, const IntType & v ) 993 { 994 IntType r ; 995 IntType u2 { u } ; 996 IntType v2 { v } ; 997 998 #ifdef DEBUG_PP_ARITH 999 cout << "gcd: u = " << u << " v = " << v << endl ; 1000 #endif 1001 1002 while (v2 != static_cast<IntType>(0u)) 1003 { 1004 r = u2 % v2 ; 1005 u2 = v2 ; 1006 v2 = r ; 1007 1008 #ifdef DEBUG_PP_ARITH 1009 cout << " r = " << r << " u2 = " << u2 << " v2 = " << v2 << endl ; 1010 #endif 1011 } 1012 1013 return u2 ; 1014 } 1015 1016 1017 /*============================================================================== 1018 | Forced Template Instantiations | 1019 ==============================================================================*/ 1020 1021 // C++ doesn't automatically generate templated classes or functions unless 1022 // we USE them on particular data types. 1023 template ppuint gcd( const ppuint &, const ppuint & ) ; 1024 template BigInt gcd( const BigInt &, const BigInt & ) ; 1025 1026 template ppuint addMod( ppuint a, ppuint b, ppuint n ) ; 1027 template ppuint timesTwoMod( ppuint a, ppuint n ) ; 1028 template ppuint multiplyMod( const ppuint a, const ppuint b, const ppuint n ) ; 1029 1030 // We already specialized this function for ppuint in the source code implementation above, so we can omit 1031 // the template instantiation: template ppuint PowerMod<ppuint>::operator()( const ppuint & a, const ppuint & n ) ; 1032 template BigInt PowerMod<BigInt>::operator()( const BigInt & a, const BigInt & n ) ; 1033 1034 template ModP<ppuint,ppsint>::ModP( ppuint p ) ; 1035 template ModP<ppuint,ppsint>::ModP( const ModP & ) ; 1036 template ppuint ModP<ppuint,ppsint>::operator()( ppsint ) ; 1037