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