/*==============================================================================
|
|  NAME
|
|     ppArith.cpp
|
|  DESCRIPTION
|
|     Miscellaneous integer multiple precision math routines.
|
|  LEGAL
|
|     Primpoly Version 9.1 - A Program for Computing Primitive
|     Polynomials.  Copyright (C) 1999-2008 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;  version 2 of the
|     License.
|
|     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, write to the Free
|     Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
|     MA  02111-1307, USA.
|
|     The author's address is artifex@seanerikoconnor.freeservers.com.
|
==============================================================================*/

/*------------------------------------------------------------------------------
|                                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.h"		  // Global functions.
#include "ppArith.h"		  // Basic arithmetic functions.
#include "ppBigInt.h"       // Arbitrary precision integer arithmetic.
#include "ppStatistics.h"   // Statistics collection for factoring and poly finding.
#include "ppFactor.h"       // Prime factorization and Euler Phi.
#include "ppParser.h"	       // Parsing of polynomials and I/O services.
#include "ppPolynomial.h"	  // Polynomial operations and mod polynomial operations.
#include "ppUnitTest.h"     // Complete unit test.

/*==============================================================================
|                                  mod                                         |
================================================================================

DESCRIPTION

     A mod function for both positive and negative arguments.

INPUT

    n (int, -infinity < n < infinity)
    p (int, p > 0)

RETURNS

    n mod p (0 <= n mod p < p)

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, FORTRAN's MOD function does the same thing.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int ModP::operator()( int n )
{

int raw_mod = n % p_ ;

if (raw_mod == 0)

    return( 0 ) ;

else if (n >= 0)             /*  mod is not 0.  n >= 0. */

    return( raw_mod ) ;

else

    return( raw_mod + p_ ) ;  /* mod is not 0.  n < 0. */

} /* =========================== end of function mod ======================== */



/*==============================================================================
|                                power_mod                                     |
================================================================================

DESCRIPTION

     Raise a positive integer to a positive power modulo an integer.

INPUT

     a (int, a >= 0)   The base.
     n (int, n >= 0)   The exponent.
     p (int, p >= 2)   The modulus.

RETURNS

      n
     a  (modulo p)
    -1 if a < 0, n < 0, or p <= 1, or 0^0 case.

EXAMPLE CALLING SEQUENCE

     power_mod( 2, 3, 7 ) => 1

METHOD

     Multiplication by repeated squaring.


--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

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

//  Out of range conditions.
if (p_ <= 1u || (a == 0u && n == 0u))
{
    // TODO fix up the throw
    ostringstream os ;
    os << "PowerMod::operator() "
       << "out of range a = " << a << " n = " << n  << " p_ = " << p_
       << " p_ <= 1u is " << (p_ <= 1u)
       << " a == 0u is " <<  (a == 0u)
       << " n == 0u is " <<  (n == 0u) 
       << " at " << __FILE__ << ": line " << __LINE__ ;
    throw PolynomialRangeError( os.str() ) ;
}

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

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

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

// Find the number of the leading bit.
int bitNum = n.maxBitNumber() ; // Number of highest possible bit.

#ifdef DEBUG_PP_ARITH
cout << "initial max bitNum = " << bitNum << endl ;
cout << "a = " << a << endl ;
#endif

while (!n.testBit( bitNum ))
    --bitNum ;

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

if (bitNum == -1)
{
    // TODO fix up the throw
    ostringstream os ;
    os << "PowerMod::operator() " << "bitNum == -1 internal error in PowerMod"
       << " at " << __FILE__ << ": line " << __LINE__ ;
    throw PolynomialRangeError( os.str() ) ;
}

#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<BigInt>( p_ ) ; // Square mod p.

    if (n.testBit( bitNum ))
        a1 = (a1 * a) % static_cast<BigInt>( 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 ;
}

// Specialize for integer.
template<>
uint PowerMod<uint>::operator()( const uint & a, const uint & 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.  Signed or
//  unsigned type for bigint shouldn't matter here, since we just mask and
//	compare bits.

int mask = (1 << (8 * sizeof( int ) - 1)) ;
int bit_count = 0 ;  // Number of bits in exponent to the right of the leading bit.
int n1 = n ;

/*  Use extra precision since we need to square an integer within the
    calculations below. */
long int product = (int) a ;

/*
    Out of range conditions.
*/

if (a < 0 || n1 < 0 || p_ <= 1 || (a == 0 && n1 == 0))
    return 0 ;

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

if (n == 0)
    return 1 ;

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

//    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( int )) - bit_count ;


/*
    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 )
{
    /*  Expose the next bit. */
    n1 <<= 1 ;

	/* Square modulo p. */
	product = (product * product) % p_ ;

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

/*  We don't lose precision because product is always reduced modulo p. */
return (int) product ;
}



/*==============================================================================
|                             is_primitive_root                                |
================================================================================

DESCRIPTION

     Test whether the number a is a primitive root of the prime p.

INPUT

     a      (int, a >= 1)     The number to be tested.
     p      (int, p >= 2)     A prime number.

RETURNS

     true      If a is a primitive root of p.
     false     If a isn't a primitive root of p or if inputs are out of range.

EXAMPLE
                                                     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.

TODO

    We don't check if p is not a prime.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

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

/*  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)))
{
	return false ;
}

/*  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)) )
{
	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.
Factor<uint> factor( p_ - 1 ) ;


//            p-1
//  Check if a    != 1 (mod p) for all primes q | (p-1).
//          ---
//           q
//  If so, we have a primitive root of p, otherwise, we exit early.
//
for (int i = 0 ;  i < factor.num() ;  ++i)
{
	/*  The prime factors can be safely cast from bigint to int
	    because p-1 needs only integer precision.  */
    if (powermod( a, (p_ - 1) / factor[ i ]) == 1)
	{
	    return false ;
	}
}

return true ;

} /* =============== end of function is_primitive_root ===================== */



/*==============================================================================
|                                 inverse_mod_p                                |
================================================================================

DESCRIPTION
                                                            -1
     Find the inverse of u modulo p.  In other words, find u,   u in GF( p ).

INPUT

     u      0 <= u < p         Number for which we want to find the inverse.
     p      p is prime >= 2    The modulus.

RETURNS

	  -1
     u  (mod p)

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.

TODO

    We don't check if p is prime.
	If we return 0, need to check this in calling routines.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int InverseModP::operator()( int u )
{
    ModP mod( p_ ) ;

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

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

	int inv_v = 0 ;

	while( v3 != 0)
	{
        q = (int)(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)
	{
		return 0 ;
	}

	return inv_v ;
}


/*==============================================================================
|                               const_coeff_test                               |
================================================================================

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 order_r() test.

INPUT

    f (int *)  nth degree monic mod p polynomial f(x).
    n (int)    Its degree.
    p (int)    Modulus for coefficient arithmetic.
    a (int)

RETURNS

    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

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

bool ArithModP::const_coeff_test( int a, int a0, int n )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int constantCoeff = a0 ;

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

ModP mod( p_ ) ; // Initialize the functionoid.

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

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

} /* =================== end of function const_coeff_test ================ */



/*==============================================================================
|                         const_coeff_is_primitive_root                        |
================================================================================

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

INPUT

    f (int *)  nth degree monic mod p polynomial f(x).
    n (int)    Its degree.
    p (int)    Modulus for coefficient arithmetic.

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

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

bool ArithModP::const_coeff_is_primitive_root( int a0, int n )
{

int constantCoeff = a0 ;

ModP mod( p_ ) ;
IsPrimitiveRoot isroot( p_ ) ;

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

return isroot( mod( constantCoeff ) ) ;

} /* ============= of function const_coeff_is_primitive_root ================ */

/*

Mathematica 

    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 ]]

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

*/

// gcd( u, v ) = gcd( v, u - q v ) for any q.
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 ;
       cout << " r = " ; printNumber( r ) ;
       cout << " u2 = " ; printNumber( u2 ) ;
       cout << " v2 = " ; printNumber( v2 ) ;
       #endif
    }

    return u2 ;
}


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

// C++ doesn't automatically generate templated classes or functions unless
// we USE them on particular data types.

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

template         PowerMod<uint>::PowerMod( const uint & p ) ;

// We already specialized this function for uint.
template BigInt  PowerMod<BigInt>::operator()( const BigInt & a, const BigInt & n ) ;