/*==============================================================================
|
|  NAME
|
|     ppArith.cpp
|
|  DESCRIPTION
|
|     Miscellaneous integer multiple precision math routines.
|
|     User manual and technical documentation are described in detail in my web page at
|     http://seanerikoconnor.freeservers.com/Mathematics/AbstractAlgebra/PrimitivePolynomials/overview.html
|
|  LEGAL
|
|     Primpoly Version 16.4 - A Program for Computing Primitive Polynomials.
|     Copyright (C) 1999-2025 by Sean Erik O'Connor.  All Rights Reserved.
|
|     This program is free software: you can redistribute it and/or modify
|     it under the terms of the GNU General Public License as published by
|     the Free Software Foundation, either version 3 of the License, or
|     (at your option) any later version.
|
|     This program is distributed in the hope that it will be useful,
|     but WITHOUT ANY WARRANTY; without even the implied warranty of
|     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
|     GNU General Public License for more details.
|
|     You should have received a copy of the GNU General Public License
|     along with this program.  If not, see http://www.gnu.org/licenses/.
|     
|     The author's address is seanerikoconnor!AT!gmail!DOT!com
|     with the !DOT! replaced by . and the !AT! replaced by @
|
==============================================================================*/

/*------------------------------------------------------------------------------
|                                Include Files                                 |
------------------------------------------------------------------------------*/

#include <cstdlib>      // abort()
#include <iostream>     // Basic stream I/O.
#include <new>          // set_new_handler()
#include <cmath>        // Basic math functions e.g. sqrt()
#include <complex>      // Complex data type and operations.
#include <fstream>      // File stream I/O.
#include <sstream>      // String stream I/O.
#include <vector>       // STL vector class.
#include <string>       // STL string class.
#include <algorithm>    // Iterators.
#include <stdexcept>    // Exceptions.
#include <cassert>      // assert()

using namespace std ;

/*------------------------------------------------------------------------------
|                                PP Include Files                              |
------------------------------------------------------------------------------*/

#include "Primpoly.hpp"		  // Global functions.
#include "ppArith.hpp"		  // Basic arithmetic functions.
#include "ppBigInt.hpp"         // Arbitrary precision integer arithmetic.
#include "ppOperationCount.hpp" // OperationCount collection for factoring and poly finding.
#include "ppFactor.hpp"         // Prime factorization and Euler Phi.
#include "ppPolynomial.hpp"	  // Polynomial operations and mod polynomial operations.
#include "ppParser.hpp"	      // Parsing of polynomials and I/O services.
#include "ppUnitTest.hpp"       // Complete unit test.

/*=============================================================================
 | 
 | NAME
 | 
 |     ModP::operator()
 | 
 | DESCRIPTION
 | 
 |     Computes k = n mod p where 0 <= k < p for both positive and negative 
 |     arguments n for p >= 1.
 |     Takes a surprising proportion of the total compute time according to the 
 |     profiler, so worth optimizing.
 | 
 | EXAMPLE
 | 
 |     The C language gives -5 % 3 = - ( -(-5) mod 3 ) = -( 5 mod 3 ) = -2.
 |     The result is nonzero, so we add p=3 to give 1.
 | 
 |     C computes -3 % 3 = - ( -(-3) mod 3  ) = 0, which we leave unchanged.
 | 
 | METHOD
 | 
 |      For n >= 0, the C language defines r = n mod p by the equation
 | 
 |          n = kp + r    0 <= r < p
 | 
 |      but when n < 0, C returns the quantity
 | 
 |          r = - ( (-n) mod p )
 | 
 |      To put the result into the correct range 0 to p-1, add p to r if
 |      r is non-zero.
 | 
 |      By the way, dear old FORTRAN's MOD function does the same thing.
 | 
 +============================================================================*/

template <typename UIntType, typename SIntType>
inline UIntType ModP<UIntType,SIntType>::operator()( SIntType n )
{
    if (p_ <= 0)
    {
        ostringstream os ;
        os << "ModP::operator()  The modulus p = " << p_ << " <= 0" ;
        throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
    }
    
    if (n >= 0)
    {
        // Avoid integer division in this special case.
        if (p_ == 2)
            return n - ((n >> 1) << 1) ;
        else
            return n % p_ ;
    }
    // Both types must be signed if we want n % p to compute correctly for n < 0.
     else
        return( (n % static_cast<SIntType>(p_)) + p_ ) ;
}

/*=============================================================================
|
| NAME
|
|     PowerMod::operator()
|
| DESCRIPTION
|
|               n 
|     Compute  a  (modulo p)  for generic integer types.
|                                                           0
|     where a >= 0, p >= 2, n >= 0.  Other cases including 0
|     will throw an ArithModPError.
|
|     We have two versions:  one is generic for any integer type, and
|     the other is specialized for ppuint.
|
| EXAMPLE
|
|     BigInt a = 2 ;
|     BigInt n = 3 ;
|     BigInt p = 7 ;
|     try 
|     {
|         PowerMod<BigInt> powermod( p ) ;
|         BigInt pwr = powermod( a, n ) ;
|     } 
|     catch( ArithModPError e )
|     { 
|         cout << "should have gotten pwr = 1\n" ;
|     }
|
| METHOD
|
|     Multiplication by repeated squaring.
|
+============================================================================*/

// Generic.
template <typename IntType>
IntType PowerMod<IntType>::operator()( const IntType & a, const IntType & n )
{
    IntType a1 = a ;

    //  Out of range conditions.
    if (a  <  static_cast<IntType>( 0u ) || 
        n  <  static_cast<IntType>( 0u ) || 
        p_ <= static_cast<IntType>( 1u ) || 
       (a  == static_cast<IntType>( 0u ) && n == static_cast<IntType>( 0u )))
    {
        ostringstream os ;
        os << "PowerMod::operator()  One or more parameters out of range:  a = "
           << a << " n = " << n  << " p_ = " << p_ ;
        throw ArithModPError( os.str(),  __FILE__, __LINE__ ) ;
    }

    //  Quick return for 0 ^ n, a^0 and a^1.
    if (a == static_cast<IntType>( 0u ))
        return static_cast<IntType>( 0u ) ;

    if (n == static_cast<IntType>( 0u ))
        return static_cast<IntType>( 1u ) ;

    if (n == static_cast<IntType>( 1u ))
        return a % static_cast<IntType>( p_ ) ;

    int bitNum = n.maxBitNumber() ; // Index to highest bit.

    #ifdef DEBUG_PP_ARITH
    cout << "initial max bitNum = " << bitNum << endl ;
    cout << "a = " << a << endl ;
    #endif
    
    // Find the index of the leading 1 bit.
    while (!n.testBit( bitNum ))
        --bitNum ;

    #ifdef DEBUG_PP_ARITH
    cout << "after skipping leading 0 bits, bitNum = " << bitNum << endl ;
    #endif

    if (bitNum == -1)
    {
        ostringstream os ;
        os << "PowerMod::operator() " << "bitNum == -1" ;
        throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
    }

    #ifdef DEBUG_PP_ARITH
    cout << "\nAfter skipping zero bits, bitNum = " << bitNum << endl ;
    #endif

    //  Exponentiation by repeated squaring.  Discard the leading 1 bit.
    //  Thereafter, square for every 0 bit;  square and multiply by a for
    //  every 1 bit.
    while ( --bitNum >= 0 )
    {
        a1 = (a1 * a1) % static_cast<IntType>( p_ ) ; // Square mod p.

        if (n.testBit( bitNum ))
            a1 = (a1 * a) % static_cast<IntType>( p_ ) ; // Times a mod p.

        #ifdef DEBUG_PP_ARITH
        cout << "S " ;
        if (n.testBit( bitNum ))
            cout << "X " ;
        cout << "Bit num = " << bitNum << " a1 = " << a1 << endl ;
        #endif
    }

    #ifdef DEBUG_PP_ARITH
        cout << "Out of the loop bitNum = " << bitNum << " a1 = " << a1 << endl ;
    #endif

    return a1 ;
}



/*=============================================================================
 | 
 | NAME
 | 
 |     addMod()
 | 
 | DESCRIPTION
 | 
 |     Computes (x + y) mod p for generic integer types, up to the maximum
 |     unsigned type.  We extend the range to the limit by testing and compensating
 |     for carries internally.
 |
 |     Example:
 |         a           = 4294967290
 |         n           = 4294967294
 |         (a + b) % n = 4294967286
 | 
 +============================================================================*/

template <typename IntType>
IntType addMod( IntType a, IntType b, IntType n )
{
    #ifdef DEBUG_PP_ARITH
    cout << "addMod" << endl ;
    cout << "    sizeof IntType = " << 8 * sizeof( IntType ) << " bits" << endl ;
    cout << "    a = " << a << " b = " << b << " n = " << n << endl ;
    #endif

    // Assume the positive arguments might not already be in the range [0, n).
    a %= n ;
    b %= n ;

    // Add with carry.  Carry bit thrown away by the language.
    IntType c = a + b ;

    #ifdef DEBUG_PP_ARITH
    cout << "    a mod n = " << a << " b mod n = " << b << endl ;
    cout << "    c = a + b (discarding carry bit) = " << c << endl ;
    #endif

    // Test for a carry.  Most computer languages don't give us access to
    // the carry bit, so we must infer whether a carry occurred during addition.
    // Logically, 
    //     (no carry for a+b)  =>  a+b <= a && a+b <= b
    // The contrapositive is
    //      a+b < a || a+b < b => carry
    if (c < a || c < b)
    {
        // Lemma 1.  If a carry occured, n < c < 2n.
        // Pf.
        //                                       m
        //     When a carry occurs, c = a + b > 2  > n.  But 0 <= a, b < n.  So n < c and a + b < n + n = 2n.
        //
        // Lemma 2.  If a carry occured, we can simply compute c - n as usual, discarding the carry.
        // Pf.
        //     Since a carry occurred, we an write
        //                  m
        //     c = a + b = 2  + c'
        //
        //     Let's assume the computer hardware does subtraction in two's complement.
        //     Now use two's complement arithmetic on c' - n and see the results.
        //
        //     c' + TwosComplement( n ) <discarding carry bit> = 
        //
        //     c' + (~n + 1) <discarding carry bit> =
        //             m           m
        //     (c' + (2 - n)) mod 2  = 
        //           m   m     m           m          m           m
        //     (c - 2 + 2  + (2 - n)) mod 2  = (c + (2 - n)) mod 2 =
        //
        //     c + TwosComplement( n ) <discarding carry bit> =
        //
        //     c - n done the usual way in computer hardware.
        //
        // Th. If a carry occurs, c mod n = c - n.
        // Pf. See Lemma 1 and 2.
        //
        // e.g.  Let m = 4 bits and
        //
        //       a  =   1 1 0 1 = 13
        //       b  =   1 1 0 1 = 13
        //       n  =   1 1 1 0 = 14
        //
        //       c  = 1 1 0 1 0
        //       c' =   1 0 1 0
        //      ~n  =   0 0 0 1 + 1 = 
        //              0 0 1 0
        //       c' -  n = c' + ~n =
        //              1 0 1 0
        //           +  0 0 1 0
        //        =     1 1 0 0 = 12 = 26 mod 14
        
        c -= n ;

        #ifdef DEBUG_PP_ARITH
        cout << "    Carry!" << " c < a = " << (c < a ) << "|| c < b = " << (c < b) << endl ;
        #endif
    }
    else
    {
        c %= n ;

        #ifdef DEBUG_PP_ARITH
        cout << "    No carry" << endl ;
        #endif
    }

    return c ;
}



/*=============================================================================
 | 
 | NAME
 | 
 |     timesTwoMod()
 | 
 | DESCRIPTION
 | 
 |     Computes (2 * x) mod p for generic integer types, up to the maximum
 |     unsigned type.  We extend the range to the limit by testing and compensating
 |     for carries internally.
 |
 |     Example:
 |         a           = 10376293541461622782
 |         n           = 18446744073709551610
 |         (2 * a) % n =  2305843009213693954
 | 
 +============================================================================*/

template <typename IntType>
IntType timesTwoMod( IntType a, IntType n )
{
    #ifdef DEBUG_PP_ARITH
    cout << "timesTwoMod" << endl ;
    cout << "    sizeof IntType = " << 8 * sizeof( IntType ) << " bits" << endl ;
    cout << "    a = " << a << " n = " << n << endl ;
    #endif

    // Assume the positive argument might not already be in the range [0, n).
    a %= n ;

    #ifdef DEBUG_PP_ARITH
    cout << "    a mod n = " << a << endl ;
    #endif

    // High bit mask.
    IntType mask = ((IntType)1 << (8 * sizeof( IntType ) - 1)) ;

    IntType c = (a << 1) ; // Lose the carry bit.

    #ifdef DEBUG_PP_ARITH
    cout << "    mask = " << hex << mask << dec << endl ;
    cout << "    c = 2 a (losing the carry bit) = " << c << endl ;
    #endif

    // High bit is set so (2 a) will have a carry.
    if (mask & a)
    {
        // See notes in function addMod() above for why this works correctly.
        c -= n ;

        #ifdef DEBUG_PP_ARITH
        cout << "    Carry!" << endl ;
        #endif
    }
    else
    {
        c %= n ;

        #ifdef DEBUG_PP_ARITH
        cout << "    No carry" << endl ;
        #endif
    }

    return c ;
}



/*=============================================================================
 | 
 | NAME
 | 
 |     multiplyMod()
 | 
 | DESCRIPTION
 | 
 |     Computes (x * y) mod p for generic integer types, up to the maximum
 |     unsigned type.  We extend the range this far by testing and compensating
 |     for carries internally.
 |
 |     Multiply a * b mod n 
 |
 |                   m-1         m-2
 |     a b = a (b   2    + b    2    + ... + b  2 + b )
 |               m-1        m-2               1      0
 |
 |     Expanding using Horner's rule,
 |
 |     a b = (0 2 + a b   ) 2 + a b   ) 2 + ... + a b ) 2 + a b
 |                     m-1         m-2               1         0
 |
 |     where b  = 0 or 1
 |            i
 |     At each step call specialized functions to do r = 2 a mod n and r = r + a mod n
 |     which will work even when carries occur.
 |
 |     The process takes O( ln n ) operations.  An example,
 |
 |          a = 4294967290 
 |          b = 4294967290
 |          n = 4294967294
 |         (a * b) % n = 16
 | 
 |     Reference:  "Comments on 'A Computer Algorithm for Calculating the Product
 |     A B Modulo M'" K. R. Sloane, IEEE Transactions on Computers, Vol C-34, No. 3,
 |     March 1985.
 | 
 +============================================================================*/

template<typename IntType>
IntType multiplyMod( const IntType a, const IntType b, const IntType n )
{
    // High bit mask.
    const int numBits = 8 * sizeof( IntType ) ;
    IntType mask = ((IntType)1 << (numBits - 1)) ;

    IntType r { 0 } ;

    for (int i = numBits-1 ;  i >= 0 ;  --i, mask >>= 1)
    {
        r = timesTwoMod( r, n ) ;

        if (mask & b)
            r = addMod( r, a, n ) ;
    }

    return r ;
}



/*=============================================================================
 | 
 | NAME
 | 
 |     PowerMod::operator()
 | 
 | DESCRIPTION
 | 
 |     Specialized for ppuint type.
 | 
 +============================================================================*/

template<>
ppuint PowerMod<ppuint>::operator()( const ppuint & a, const ppuint & n )
{
    //  mask = 1000 ... 000.  That is, all bits of mask are zero except for the
    //  most significant bit of the computer word holding its value.
    ppuint mask{ static_cast<ppuint>(1) << (8 * sizeof( ppuint ) - 1) } ;
    int  bit_count{ 0 } ;  // Number of bits in exponent to the right of the leading bit.
    ppuint n1{ n } ;
    ppuint product{ a }  ;

    // Out of range conditions. 
    if (p_ <= 1 || (a == 0 && n1 == 0))
    {
        ostringstream os ;
        os << "PowerMod<ppuint>::operator() parameters are out of range:  a = " << a << " n = " << n  << " p_ = " << p_ ;
        throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
    }

    //                    n   0      1
    //  Quick return for 0 , a  and a
    if (a == 0)
        return 0 ;

    if (n == 0)
        return 1 ;

    if (n == 1)
        return a % p_ ;

    #ifdef DEBUG_PP_ARITH
    cout << "a                    = " << a << endl ;
    cout << "n                    = " << n << endl ;
    cout << "p                    = " << p_ << endl ;
    cout << "n1 (before shifting) = " << n1 << endl ;
    #endif

    // Advance the leading bit of the exponent up to the word's left hand boundary.  
    // Count how many bits were to the right of the leading bit.
    while (! (n1 & mask))
    {
        n1 <<= 1 ;
        ++bit_count ;
    }

    bit_count = (8 * sizeof( ppuint )) - bit_count ;

    #ifdef DEBUG_PP_ARITH
    cout << "n1        = " << n1 << endl ;
    cout << "mask      = " << mask << endl ;
    cout << "bit_count = " << bit_count << endl ;
    #endif

    // Exponentiation by repeated squaring.  Discard the leading 1 bit.
    // Thereafter, square for every 0 bit;  square and multiply by x for every 1 bit.
    while ( --bit_count > 0 )
    {
        #ifdef DEBUG_PP_ARITH
        cout << "product (before squaring) = " << product << " n1 = " << n1 << endl ;
        #endif

        // Expose the next bit.
        n1 <<= 1 ;

        // Out of range conditions. 
        if (product > BigInt::getBase() || a > BigInt::getBase())
        {
            // Square modulo p.
            product = multiplyMod( product, product, p_ ) ;

            //  Leading bit is 1: multiply by a modulo p.
            if (n1 & mask)
                product = multiplyMod( a, product, p_ ) ;
        }
        else
        {
            // Square modulo p.
            product = (product * product) % p_ ;

            //  Leading bit is 1: multiply by a modulo p.
            if (n1 & mask)
                product = (a * product) % p_ ;
        }

        #ifdef DEBUG_PP_ARITH
        cout << "S " ;
        if (n1 & mask)
            cout << "X " ;
        cout << "product = " << product << " n1 = " << n1 << endl ;
        #endif
    }

    return product ;
}



/*=============================================================================
|
| NAME
|
|     IsPrimitiveRoot::operator()
|
| DESCRIPTION
|
|     Returns true if a is a primitive root of prime p, and false otherwise.
|     Throws ArithModPError if a < 1 or p < 2 or p is even.
|     p must be a prime number, but we only test that p is 2 or odd.
|
| EXAMPLE
|
|     ppuint p{ 7 } ;  ppuint a{ 3 } ;
|     IsPrimitiveRoot isRoot( p ) ;
|     try 
|     {
|        bool isPrim = isRoot( a ) ;
|     } 
|     catch( ArithModPError e )
|     {  
|         cout << "out of range" << endl ; 
|     }
|                                                     1      2      3
|     3 is a primitive root of the prime p = 7 since 3 = 3, 3 = 2, 3 = 6,
|      4      5                                       p-1    6
|     3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3    = 3 = 1 (mod 7).
|
|     So 3 is a primitive root of 7 because it has maximal period.  a = 2
|                                                      3
|	  isn't a primitive root of p = 7 because already 2 = 1 (mod 7).
|
| METHOD
|
|    From number theory, a is a primitive root of the prime p iff
|     (p-1)/q
|    a        != 1 (mod p) for all prime divisors q of (p-1).
|                            (p-1)
|    Don't need to check if a     = 1 (mod p):  Fermat's little
|    theorem guarantees it.
|
+============================================================================*/

bool IsPrimitiveRoot::operator()( ppuint a )
{
    PowerMod<ppuint> powermod( p_ ) ;

    //  a = 0 (mod p);  Zero can't be a primitive root of p.
    if (a == 0)
        return false ;


    //  Error for out of range inputs, including p
    //  being an even number greater than 2.
    if (p_ < 2 || a < 1 || (p_ > 2 && (p_ % 2 == 0)))
    {
        ostringstream os ;
        os << "IsPrimitiveRoot::operator()  Inputs are out of range: p = " << p_ << " a = " << a ;
        throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
    }

    //  Special cases:
    //  1 is the only primitive root of 2;  i.e. 1 generates the unit elements
    //  of GF( 2 );  2 is the only primitive root of 3, and 2 and 3 are the only
    //  primitive roots of 5.
    //
    if ( (p_ == 2  &&  a == 1) ||
         (p_ == 3  &&  a == 2) ||
         (p_ == 5  && (a == 2  || a == 3)) ||
         (p_ == 7  && (a == 3  || a == 5)) ||
         (p_ == 11 && (a == 2  || a == 6   || a == 7 || a == 8)) ||
         (p_ == 13 && (a == 2  || a == 6   || a == 7 || a == 11)))
    {
        return true ;
    }

    //  Reduce a down modulo p.
    a = a % p_ ;

    //  a = 0 (mod p);  Zero can't be a primitive root of p.
    if (a == 0)
        return false ;

    //  Factor p-1 into distinct primes.
    Factorization<ppuint> factorization( p_ - 1 ) ;

    //            p-1
    //           a    
    //  Test    ---    != 1 (mod p) for all primes q | (p-1).
    //           q
    //  If so, we have a primitive root of p, otherwise, we exit early.
    //
    for (unsigned int i = 0 ;  i < factorization.numDistinctFactors() ;  ++i)
    {
        if (powermod( a, (p_ - 1) / factorization.primeFactor( i )) == 1)
        {
            return false ;
        }
    }

    return true ;
} 



/*=============================================================================
 | 
 | NAME
 | 
 |     inverse_mod_p
 | 
 | DESCRIPTION
 | 
 |                                                             -1
 |      Find the inverse of u modulo p.  In other words, find u   (mod p)
 |      p >= 2.
 | 
 | EXAMPLE
 |                           -1
 |      For p=7, and u = 2, u   = 4 because 2 * 4 = 8 (mod 7) = 1.
 | 
 | METHOD
 | 
 |     Do extended Euclid's algorithm on u and v to find u1, u2, and u3 such that
 | 
 |         u u1 + v u2 = u3 = gcd( u, v ).
 | 
 |     Now let v = p = prime number, so gcd = u3 = 1.  Then we get
 | 
 |         u u1 + p ? = 1
 | 
 |     or u u1 = 0 (mod p) which makes u (mod p) the unique multiplicative
 |     inverse of u.
 | 
 +============================================================================*/

ppsint InverseModP::operator()( ppsint u )
{
    ModP<ppsint,ppsint> mod( p_ ) ;

    //  Do Euclid's algorithm to find the quantities.
	ppsint t1 = 0 ;
	ppsint t3 = 0 ;
	ppsint q  = 0 ;

    ppsint u1 = 1 ;
    ppsint u3 = u ;
    ppsint v1 = 0 ;
    ppsint v3 = p_ ;

	ppsint inv_v = 0 ;

	while( v3 != 0)
	{
        // Casting to do a floor function.
        q = static_cast<ppsint>(u3 / v3) ;

		t1 = u1 - v1 * q ;
		t3 = u3 - v3 * q ;

		u1 = v1 ;
		u3 = v3 ;

		v1 = t1 ;
		v3 = t3 ;
	}

    inv_v = mod( u1 ) ;

    //                       -1
	// Self check:  does u  u   = 1 (mod p)?
    //
	if ( mod( u * inv_v ) != 1)
	{
        ostringstream os ;
        os << "InverseModP::operator():  inverse self check failed:  u = " << u << " * u^(-1) = " << inv_v  << " != 1" ;
        throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
	}

	return inv_v ;
}


/*=============================================================================
 | 
 | NAME
 | 
 |     constCoeffTest
 | 
 | DESCRIPTION
 | 
 |                   n
 |   Test if a = (-1)  a  (mod p) where a  is the constant coefficient
 |                      0                0
 |   of polynomial f(x) and a is the number from the orderR() test.
 |   Return true if we pass the test, false otherwise. 
 | 
 | EXAMPLE
 |                                  11     3
 |     Let p = 5, n = 11, f( x ) = x  + 2 x + 1, and a = 4.
 |     Thus a  = 1.
 |           0
 |     Then return true since
 | 
 |         11
 |     (-1)  * 1 (mod 5) = -1 (mod 5) = 4 (mod 5).
 | 
 | METHOD
 |                         n
 |     We test if (a - (-1) a ) mod p = 0.
 |                           0
 | 
 +============================================================================*/

bool ArithModP::constCoeffTest( ppsint a, ppsint a0, int n )
{
    ppsint constant_coeff = a0 ;

    ModP<ppuint,ppsint> mod( p_ ) ; // Initialize the functionoid.

    if (n % 2 != 0)
        constant_coeff = -constant_coeff ;    // (-1)^n < 0 for odd n.

    return( (mod( a - constant_coeff ) == 0) ? true : false ) ;
} 



/*=============================================================================
 | 
 | NAME
 | 
 |     constCoeffIsPrimitiveRoot
 | 
 | DESCRIPTION
 | 
 |               n
 |   Test if (-1)  a  (mod p) is a primitive root of the prime p where
 |                  0
 |   a  is the constant term of the polynomial f(x).
 |    0
 | 
 | EXAMPLE
 |                                  11     3
 |     Let p = 7, n = 11, f( x ) = x  + 2 x + 4.  Thus a  = 4.
 |                                                      0
 |     Then return true since
 | 
 |         11
 |     (-1)  * 4 (mod 7) = -4 (mod 7) = 3 (mod 7), and 3 is a
 |                                          1      2      3
 |     primitive root of the prime 7 since 3 = 3, 3 = 2, 3 = 6,
 |      4      5                                       7-1
 |     3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3   = 1 (mod 7).
 | 
 | METHOD
 | 
 |                    n
 |     We test if (-1) a  mod p is a primitive root mod p.
 |                      0
 |
 +============================================================================*/

bool ArithModP::constCoeffIsPrimitiveRoot( ppuint a0, int n )
{
    ppsint constant_coeff = a0 ;

    ModP<ppuint,ppsint> mod( p_ ) ;
    IsPrimitiveRoot isroot( p_ ) ;

    // (-1)^n < 0 for odd n.
    if (n % 2 != 0)
        constant_coeff = -constant_coeff ;

    return isroot( mod( constant_coeff ) ) ;
} 



/*=============================================================================
|
| NAME
|
|     gcd
|
| DESCRIPTION
|
|    Euclid's method for computing greatest common divisor.
|
| METHOD
|
|     Iterate the equation 
|
|        gcd( u, v ) = gcd( v, r )
|   
|     where r = | u / v | = u - q v
|              --     --
|     until v is zero, terminating in gcd( u, 0 ) = u.
|
|     Let's use Mathematica to illustrate the method for large integers,
|
|     gcd[ u_, v_] :=
|         (* Compute GCD using Euclid's method. *)
|         Module[ {u2, v2, r}, u2 = u ; v2 = v ; 
|                  Print[ "Computing GCD( u = ", u, " v = ", v , " )" ] ;
|                  While[ v2 != 0, 
|                         r = Mod[u2, v2]; u2 = v2; v2 = r;
|                         Print[ "r = ", r, " u2 = ", u2, " v2 = ", v2]];
|                         Return[ u2 ]]
|
|     gcd[ 779953197883173551166308319545, 1282866356929526866866376009397 ]
|
|     Computing GCD( u = 779953197883173551166308319545 v = 1282866356929526866866376009397 )
|        u = 779953197883173551166308319545  v = 1282866356929526866866376009397
|        r = 779953197883173551166308319545 u2 = 1282866356929526866866376009397  v2 = 779953197883173551166308319545
|        r = 502913159046353315700067689852 u2 =  779953197883173551166308319545  v2 = 502913159046353315700067689852
|        r = 277040038836820235466240629693 u2 =  502913159046353315700067689852  v2 = 277040038836820235466240629693
|        r = 225873120209533080233827060159 u2 =  277040038836820235466240629693  v2 = 225873120209533080233827060159
|        r =  51166918627287155232413569534 u2 =  225873120209533080233827060159  v2 =  51166918627287155232413569534
|        r =  21205445700384459304172782023 u2 =   51166918627287155232413569534  v2 =  21205445700384459304172782023
|        r =   8756027226518236624068005488 u2 =   21205445700384459304172782023  v2 =   8756027226518236624068005488
|        r =   3693391247347986056036771047 u2 =    8756027226518236624068005488  v2 =   3693391247347986056036771047
|        r =   1369244731822264511994463394 u2 =    3693391247347986056036771047  v2 =   1369244731822264511994463394
|        r =    954901783703457032047844259 u2 =    1369244731822264511994463394  v2 =    954901783703457032047844259
|        r = 414342948118807479946619135 u2 =  954901783703457032047844259 v2 = 414342948118807479946619135
|        r = 126215887465842072154605989 u2 =  414342948118807479946619135 v2 = 126215887465842072154605989
|        r = 35695285721281263482801168 u2 =  126215887465842072154605989 v2 = 35695285721281263482801168
|        r = 19130030301998281706202485 u2 =  35695285721281263482801168 v2 = 19130030301998281706202485
|        r = 16565255419282981776598683 u2 =  19130030301998281706202485 v2 = 16565255419282981776598683
|        r = 2564774882715299929603802 u2 =  16565255419282981776598683 v2 = 2564774882715299929603802
|        r = 1176606122991182198975871 u2 =  2564774882715299929603802 v2 = 1176606122991182198975871
|        r = 211562636732935531652060 u2 =  1176606122991182198975871 v2 = 211562636732935531652060
|        r = 118792939326504540715571 u2 =  211562636732935531652060 v2 = 118792939326504540715571
|        r = 92769697406430990936489 u2 =  118792939326504540715571 v2 = 92769697406430990936489
|        r = 26023241920073549779082 u2 =  92769697406430990936489 v2 = 26023241920073549779082
|        r = 14699971646210341599243 u2 =  26023241920073549779082 v2 = 14699971646210341599243
|        r = 11323270273863208179839 u2 =  14699971646210341599243 v2 = 11323270273863208179839
|        r = 3376701372347133419404 u2 =  11323270273863208179839 v2 = 3376701372347133419404
|        r = 1193166156821807921627 u2 =  3376701372347133419404 v2 = 1193166156821807921627
|        r = 990369058703517576150 u2 =  1193166156821807921627 v2 = 990369058703517576150
|        r = 202797098118290345477 u2 =  990369058703517576150 v2 = 202797098118290345477
|        r = 179180666230356194242 u2 =  202797098118290345477 v2 = 179180666230356194242
|        r = 23616431887934151235 u2 =  179180666230356194242 v2 = 23616431887934151235
|        r = 13865643014817135597 u2 =  23616431887934151235 v2 = 13865643014817135597
|        r = 9750788873117015638 u2 =  13865643014817135597 v2 = 9750788873117015638
|        r = 4114854141700119959 u2 =  9750788873117015638 v2 = 4114854141700119959
|        r = 1521080589716775720 u2 =  4114854141700119959 v2 = 1521080589716775720
|        r = 1072692962266568519 u2 =  1521080589716775720 v2 = 1072692962266568519
|        r = 448387627450207201 u2 =  1072692962266568519 v2 = 448387627450207201
|        r = 175917707366154117 u2 =  448387627450207201 v2 = 175917707366154117
|        r = 96552212717898967 u2 =  175917707366154117 v2 = 96552212717898967
|        r = 79365494648255150 u2 =  96552212717898967 v2 = 79365494648255150
|        r = 17186718069643817 u2 =  79365494648255150 v2 = 17186718069643817
|        r = 10618622369679882 u2 =  17186718069643817 v2 = 10618622369679882
|        r = 6568095699963935 u2 =  10618622369679882 v2 = 6568095699963935
|        r = 4050526669715947 u2 =  6568095699963935 v2 = 4050526669715947
|        r = 2517569030247988 u2 =  4050526669715947 v2 = 2517569030247988
|        r = 1532957639467959 u2 =  2517569030247988 v2 = 1532957639467959
|        r = 984611390780029 u2 =  1532957639467959 v2 = 984611390780029
|        r = 548346248687930 u2 =  984611390780029 v2 = 548346248687930
|        r = 436265142092099 u2 =  548346248687930 v2 = 436265142092099
|        r = 112081106595831 u2 =  436265142092099 v2 = 112081106595831
|        r = 100021822304606 u2 =  112081106595831 v2 = 100021822304606
|        r = 12059284291225 u2 = 100021822304606  v2 = 12059284291225
|        r = 3547547974806 u2 = 12059284291225  v2 = 3547547974806
|        r = 1416640366807 u2 = 3547547974806 v2  = 1416640366807
|        r = 714267241192 u2 = 1416640366807 v2  = 714267241192
|        r = 702373125615 u2 = 714267241192 v2 =  702373125615
|        r = 11894115577 u2 = 702373125615 v2 =  11894115577
|        r = 620306572 u2 = 11894115577 v2 =  620306572
|        r = 108290709 u2 = 620306572 v2 =  108290709
|        r = 78853027 u2 = 108290709 v2 =  78853027
|        r = 29437682 u2 = 78853027 v2 =  29437682
|        r = 19977663 u2 = 29437682 v2 =  19977663
|        r = 9460019 u2 = 19977663 v2 = 9460019
|        r = 1057625 u2 = 9460019 v2 = 1057625
|        r = 999019 u2 = 1057625 v2 = 999019
|        r = 58606 u2 = 999019 v2 = 58606
|        r = 2717 u2 = 58606 v2 = 2717
|        r = 1549 u2 = 2717 v2 = 1549
|        r = 1168 u2 = 1549 v2 = 1168
|        r = 381 u2 = 1168 v2 = 381
|        r = 25 u2 = 381 v2 = 25
|        r = 6 u2 = 25 v2 = 6
|        r = 1 u2 = 6 v2 = 1
|        r = 0 u2 = 1 v2 = 0
|        ---> 1
|
+============================================================================*/

template <typename IntType>
IntType gcd( const IntType & u, const IntType & v )
{
    IntType r ;
    IntType u2 { u } ;
    IntType v2 { v }  ;

    #ifdef DEBUG_PP_ARITH
    cout << "gcd:  u = " << u << " v = " << v << endl ;
    #endif

    while (v2 != static_cast<IntType>(0u))
    {
        r  = u2 % v2 ;
        u2 = v2 ;
        v2 = r ;

       #ifdef DEBUG_PP_ARITH
       cout << "  r = " << r << " u2 = " << u2 << " v2 = " << v2 << endl ;
       #endif
    }

    return u2 ;
}


/*==============================================================================
|                     Forced Template Instantiations                           |
==============================================================================*/

// C++ doesn't automatically generate templated classes or functions until they
// are first used.  So depending on the order of compilation, the linker may say
// the templated functions are undefined.
//
// Therefore, explicitly instantiate ALL templates here:

template ppuint   gcd( const ppuint &, const ppuint & ) ;
template BigInt   gcd( const BigInt &, const BigInt & ) ;

template ppuint addMod( ppuint a, ppuint b, ppuint n ) ;
template ppuint timesTwoMod( ppuint, ppuint ) ;
template ppuint multiplyMod( const ppuint, const ppuint, const ppuint ) ;

// We already specialized this function for ppuint in the source code implementation above, so we can omit
// the template instantiation:   template ppuint PowerMod<ppuint>::operator()( const ppuint & a, const ppuint & n ) ;
template BigInt  PowerMod<BigInt>::operator()( const BigInt &, const BigInt & ) ;

template ModP<ppuint,ppsint>::ModP( ppuint ) ;
template ModP<ppuint,ppsint>::ModP( const ModP & ) ;
template ppuint ModP<ppuint,ppsint>::operator()( ppsint ) ;