/*
    // Temporary copy for the power table.
    // If it throws bad_alloc, we are OK, since
    // nothing got constructed, and the temporary copy's destructor gets
    // called during stack unwind.
    vector< vector<int> > tempPowerTable ;

    for (int i = 0 ;  i < powerTable_.size() ;  ++i)
       	swap( powerTable_[ i ], tempPowerTable[ i ] ) ;
   	swap( powerTable_, tempPowerTable ) ;

//              n
//     Compute p  = maximum number of polynoimals to to test,
//     which is the largest number the multiple precision arithmetic
//     system will be allocated to handle,
//               n
//             p  - 1
//     and r = ------.
//             p - 1
//
BigInt::reserve( p, n ) ;

BigInt max_num_poly = pow( p, n ) ;
BigInt r            = max_num_poly / (p - 1) ;

//  Factor r into distinct primes.
cout << "\nFactoring r =  %s into\n    " << r.to_string() ;

Factor<BigInt> factor( r ) ;

int prime_count = factor( r, primes, count ) ;   	


It would also be nice to be able to do the same from strings back to
numbers, no? But the C++ standard disallows it:

int convert(char*);
unsigned long convert(char*);
float convert(char*);
// and so forth

because we can't include the return type as part of the signature.

We *can*, however, do this:

class NumericResultProxy
{
public:
        NumericResultProxy(char* str)
        { m_stringVersion = new char[strlen(str)+1]; strcpy(m_stringVersion,
str); }
        ~NumericResultProxy()
        { delete [] m_stringVersion; }

        operator int()
        { return atoi(m_stringVersion); }
        operator unsigned long()
        { return atol(m_stringVersion); }
        operator float()
        { return atof(m_stringVersion); }
        // and so forth, for each distinct type we care to convert to

};

NumericResultProxy convert(char* numStr)
{
        return NumericResultProxy(numStr);

}

The trick here, of course, is that the unnamed temporary
NumericResultProxy object, when it is tried to be applied as a int or
long to a function or as an rvalue in an assignment statement, silently
invokes the appropriate user-defined conversion operator, and performs
the translation at that point. We get all the benefits of
overloads-by-return-type, without having to change a single line of the
C++ standard. 
*/

/*==============================================================================
|
|  FILE NAME
|
|     ppPolynomial.cpp
|
|  DESCRIPTION
|
|     Polynomial arithmetic and exponentiation routines.
|
|  NOTES
|
|  LEGAL
|
|     Primpoly Version 7.7 - 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 <iomanip>      // I/O manipulators.
#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 "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.

Polynomial::Polynomial()
                throw( PolynomialError, bad_exception )
                : f_()
                , n_( 0 )
                , p_( 2 )
                , mod( p_ )
{
    // f(x) = 0
    f_.push_back( 0 ) ;
}



/*=============================================================================
|
| NAME
|
|     ~Polynomial (destructor)
|
| DESCRIPTION
|
|     Destructor.
|
| SYNTAX
|
| EXAMPLE
|
|     void add1( const Polynomial & f )
|     {
|         const Polynomial g = 1 ;
|         return f + g ;
|         // Destructor for g automatically called as f goes out of scope.
|     }
|
| NOTES
|
| BUGS
|
+============================================================================*/

Polynomial::~Polynomial()
{
    // vector f_ frees itself and mod_ calls its own destructor.
}



/*=============================================================================
|
| NAME
|
|     Polynomial (copy constructor)
|
| DESCRIPTION
|
|     Copy constructor.
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f ;
|         Polynomial f( g ) ;
|     }
|     catch( PolynomialError & e )
|     {
|         cout << "Error in constructing polynomial f(x) or g(x)" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/
 
Polynomial::Polynomial( const Polynomial & g )
                throw( PolynomialError, bad_exception )
                : f_( g.f_ )
                , n_( g.n_ )
                , p_( g.p_ )
                , mod( p_ )
{
    // The classes in the initializer above all copy themselves.
}



/*=============================================================================
|
| NAME
|
|     operator=
|
| DESCRIPTION
|
|     Assigment operator for polynomials, f( x ) = g( x )
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f ;
|         Polynomial g ;
|         f = g ;
|     }
|     catch( PolynomialError & e )
|     {
|         cout << "Error in constructing polynomial f(x)" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

Polynomial & Polynomial::operator=( const Polynomial & g )
                                  throw( PolynomialError, bad_exception )
{
	// If assigned to itself, return a reference to the object itself,
    // without doing anything further.
    if (this == &g)
        return *this ;

    // Assign the scalars first.
    n_ = g.n_ ;
    p_ = g.p_ ;

    // And the modulus functionoid.
    mod( g.p_ ) ;

    // Create a temporary copy of the coefficients.
    try
    {
        vector<int> tempCoeff( g.f_ ) ;

  	    // Exchange polynomial coefficients.  
        // Cannot throw any exceptions.
	    swap( f_, tempCoeff ) ;
    }
    // Failed to construct tempPoly.
    catch( bad_alloc & e )
    {
        throw PolynomialError( "memory failure in Polynomial.operator=()" ) ;
    }


    // Return a reference to the altered object.
    return *this ;
}



/*=============================================================================
|
| NAME
|
|     Polynomial, construct from string.
|
| DESCRIPTION
|
|     Construct a polynomial from a string.
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f( "x^2 + 2x + 1, 3" ) ;
|     }
|     catch( PolynomialRangeError & e )
|     {
|         cout << "Error in construction polynomial f(x) from string" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

Polynomial::Polynomial( string s, uint p )
                throw( PolynomialRangeError, bad_exception )
                : mod( 0 )
{
    //  The polynomial string must have at least one character in it.
    if (s.empty())
        throw PolynomialRangeError( "polynomial string is empty" ) ;

    try
    {
        // Initialize an LALR(1) parser for polynomials.
        Parser pp ;
        
        Value v = pp.parse( s ) ;
        
        // Get the modulus.
        p_ = v.scalar_ ;
		
		// If the modulus is explicitly specified, use that instead
		// of the polynomial's modulus.
		if (p > 0)
		    p_ = p ;

        // Modulus must be in range.
        if (p_ <= 0)
		{
		    ostringstream os ;
            os << "Error converting polynomial from string:" << s 
			   << "modulus <= 0"
               << " at " << __FILE__ << ": line " << __LINE__ ;
            throw PolynomialRangeError( os.str() ) ;
        }
		
        mod.set( p_ ) ;

        // Check the degree of the polynomial.
        n_ = v.f_.size() - 1 ;
        if (n_ < 0)
        {
            ostringstream os ;
            os << "Error converting polynomial from string: " << s
			   << "with n = " << n_ << " and p = " << p_ 
               << " at " << __FILE__ << ": line " << __LINE__ ;
            throw PolynomialRangeError( os.str() ) ;
        }
  
        // Reduce all the polynomial coefficients modulo p.
        vector< int >::iterator i ;
        for (i = v.f_.begin() ;  i != v.f_.end() ;  ++i)
            *i = mod( *i ) ;

        // Copy over the polynomial coefficients.
        f_ = v.f_ ;
    }
    catch( ParserError & e )
    {
            ostringstream os ;
            os << "Error in parser converting polynomial from string: " << s
			   << "with n = " << n_ << " and p = " << p_ 
               << " at " << __FILE__ << ": line " << __LINE__ ;
            throw PolynomialRangeError( os.str() ) ;
    }
}


/*=============================================================================
|
| NAME
|
|     string
|
| DESCRIPTION
|
|     Convert a polynomial to a string.
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f( "x^2 + 1,3" ) ;
|         string poly = f ;
|     }
|     catch( PolynomialRangeError & e )
|         cout << "Error in construction polynomial f(x) from string" << endl ;
|     catch( bad_exception & e )
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

// Operator casting to string type.
Polynomial::operator string() const
                 throw( PolynomialRangeError, bad_exception )
{
    // Set up a string stream for convenience.
    ostringstream os ;

    // Spin out the polynomial coefficients from high to low degree.
    // Special case of f(x) = const.
    if (f_.size() == 1)
    {
        if (!(os << f_[ 0 ])) 
		{
            ostringstream os ;
            os << "Error in converting polynomial to string: "
			   << "with n = " << n_ << " and p = " << p_ 
               << " at " << __FILE__ << ": line " << __LINE__ ;
            throw PolynomialRangeError( os.str() ) ;
		}
    }
    else
    {    
		int lowestDegreeTerm = -1 ;
		for (int deg = n_ ;  deg >= 0 ;  --deg)
			if (f_[ deg ] != 0)
			    lowestDegreeTerm = deg ;

        for (int deg = n_ ;  deg >= 0 ;  --deg)
        {
			if (f_[ deg ] != 0)
            {
                // x^n has a nonzero coefficient.
                int coeff = f_[ deg ] ;
                                 
                 // Print coeff of x^n unless it is 1.
				 // But print the constant term regardless.
                if (coeff != 1 || deg == 0)
                    if (!(os << coeff)) 
	             	{
						ostringstream os ;
						os << "Error in converting polynomial to string: "
						<< "with n = " << n_ << " and p = " << p_ 
						<< " at " << __FILE__ << ": line " << __LINE__ ;
						throw PolynomialRangeError( os.str() ) ;
	             	}
                
                // Print x (no exponent).
                if (deg == 1)
				{
                    if ( !(os << " x ") )
	             	{
						ostringstream os ;
						os << "Error in converting polynomial to string: "
						<< "with n = " << n_ << " and p = " << p_ 
						<< " at " << __FILE__ << ": line " << __LINE__ ;
						throw PolynomialRangeError( os.str() ) ;
	             	}
				}
                // Print x^n ... x^2
                else if (deg != 0)
				{
                    if (!(os << " x ^ " << deg << " "))
	             	{
						ostringstream os ;
						os << "Error in converting polynomial to string: "
						<< "with n = " << n_ << " and p = " << p_ 
						<< " at " << __FILE__ << ": line " << __LINE__ ;
						throw PolynomialRangeError( os.str() ) ;
	             	}
				}
                
				// Print +, but only when followed by a lower degree term or constant.
                // e.g. x^2 + 2 x, x + 3
                if (lowestDegreeTerm != -1 && deg > lowestDegreeTerm)
                    if (!(os << " + "))
	             	{
						ostringstream os ;
						os << "Error in converting polynomial to string: "
						<< "with n = " << n_ << " and p = " << p_ 
						<< " at " << __FILE__ << ": line " << __LINE__ ;
						throw PolynomialRangeError( os.str() ) ;
	             	}
            } // end coeff != 0.
        } // end for deg
    } // end f(x) = const.

    // Print out the modulus.
    if (!(os << ", " << p_  ))
	             	{
						ostringstream os ;
						os << "Error in converting polynomial to string: "
						<< "with n = " << n_ << " and p = " << p_ 
						<< " at " << __FILE__ << ": line " << __LINE__ ;
						throw PolynomialRangeError( os.str() ) ;
	             	}

    // Return the string from the stream.
    return os.str() ;
}



/*=============================================================================
|
| NAME
|
|     <<
|
| DESCRIPTION
|
|     Polynomial stream output.
|
| SYNTAX
|
|     try 
|     {
|         
|         Polynomial f( "x^2 + 1,3" ) ;
|         cout << f << endl ;
|     }
|     catch( PolynomialRangeError & e )
|     {
|         cout << "Error in outputting polynomial f(x)" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

ostream & operator<<( ostream & out, const Polynomial & p )
                  throw( PolynomialRangeError, bad_exception )
{
    // Cast to polynomial to a string first, then output as a string.  
    // May throw a PolynomialRangeError.
    out << static_cast<string>( p ) << endl ;

	return out ;
}



/*=============================================================================
|
| NAME
|
|     >>
|
| DESCRIPTION
|
|     Polynomial stream input.
|
| SYNTAX
|
|     try 
|     {
|         
|         Polynomial f ;
|         cin >> f ;
|     }
|     catch( PolynomialRangeError & e )
|     {
|         cout << "Error in inputting polynomial f(x)" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

istream & operator>>( istream & in, Polynomial & p )
                  throw( PolynomialRangeError, bad_exception )
{
    // Input as a string.
    string s ;
    in >> s ;

    // Copy into argument polynomial.
    // TODO handle exceptions.
    p = Polynomial( s ) ;

	return in ;
}



/*=============================================================================
|
| NAME
|
|     Polynomial operator[]
|
| DESCRIPTION
|
|     Polynomial indexing operator which allows an lvalue:  f[ 33 ] = 2 ;
|     If we don't have a coefficient of this degree, create it and backfill
|     earlier coefficients with zero.
|     
|     Throws an exception if out of bounds.
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f( "2x^2 + 1, 3" ) ;
|
|         f[ 5 ] = 2 ;
|
|         // Now f(x) = 2 x^5 + 0 x^4 + 0 x^3 + 2 x^2 + 0 x + 1
|         // f_.size() => 5 + 1 = 6
|
|     }
|     catch( PolynomialRangeError & e )
|     {
|         cout << "Error in assigning polynomial f(x) coefficient" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

// TODO don't need to throw a range error, unless resize fails.
// TODO use reserve in main()
int & Polynomial::operator[]( int i ) throw( PolynomialRangeError, bad_exception )
{
    // We attempt to access beyond the current degree.
    if (i > n_)
    {
        // Extend the vector size with zeros.
        f_.resize( i+1, 0 ) ;
        n_ = i ;
    }

    // Return a reference to the coefficient.
    return f_[ i ] ;
}



/*=============================================================================
|
| NAME
|
|     Polynomial operator[]
|
| DESCRIPTION
|
|     Polynomial indexing operator for read only access:  int coeff = f[ 33 ] ;
|     Throws an exception if out of bounds.
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f ;
|         int value = f[ 33 ] ;
|     }
|     catch( PolynomialRangeError & e )
|     {
|         cout << "Error in getting polynomial f(x) coefficient" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

const int Polynomial::operator[]( int i ) const throw( PolynomialRangeError, bad_exception )
{
    // We throw our own exception here.
    if (i > n_)
	             	{
						ostringstream os ;
						os << "Error accessing polynomial with coefficients p[0]...p[n] = (" ;
                        for (int j = 0 ;  j <= n_ ;  ++j)
                            os << f_[ j ] << " " ;
                        os << ")" << endl
                           << " at index i = " << i
						   << " of degree n = " << n_ << " modulo p = " << p_
						   << " at " << __FILE__ << ": line " << __LINE__ ;
						throw PolynomialRangeError( os.str() ) ;
	             	}

    return f_[ i ] ;
}


/*=============================================================================
|
| NAME
|
|     deg
|
| DESCRIPTION
|
|     Return the degree of f(x).
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f ;
|         cout << "Degree of f(x) = " << f.deg() << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

int Polynomial::deg() const
                throw()
{
    return n_ ;
}

       // Return the modulus p of f(x).
int Polynomial::modulus() const
                  throw()
{
                         return p_ ;
                         }

void Polynomial::setModulus( int p )
			throw( PolynomialRangeError, bad_exception )
			{
			p_ = p ;
	        // And the modulus functionoid.
            mod( p_ ) ;
			}

/*=============================================================================
|
| NAME
|
|     Polynomial operator+=
|
| DESCRIPTION
|
|     Polynomial sum f(x) += g(x)
|
|
| SYNTAX
|
|     try 
|     {
|         Polynomial f ;
|         Polynomial g ;
|         f += g ;
|     }
|     catch( PolynomialRangeError & e )
|     {
|         cout << "Error in polynomial sum f(x) += g(x)" << endl ;
|     }
|     catch( bad_exception & e )
|     {
|         cout << "Unexpected exception was thrown.  Please write the author" << endl ;
|     }
|
| EXAMPLE
|
| NOTES
|
| BUGS
|
+============================================================================*/

Polynomial & Polynomial::operator+=( const Polynomial & g )
                         throw( PolynomialRangeError, bad_exception )
{
    // f(x) = x^2 +   + 1
    // g(x) =       x + 3
	//
    // f(x) =       x + 3
    // g(x) = x^2 +   + 1
	//
    int minDeg = (n_ < g.n_) ? n_ : g.n_ ;

    // Add coefficients modulo p.
    for (int i = 0 ;  i <= minDeg ;  ++i)
        f_[ i ] = mod( f_[ i ] + g.f_[ i ] ) ;
		
	// If g(x) has larger degree, extend f(x) and copy the coefficients of g(x).
	if (g.n_ > n_)
	{
	    // Extend the vector size with zeros.
        f_.resize( g.n_ + 1, 0 ) ;

	    for (int i = n_ + 1 ;  i <= g.n_ ;  ++i)
            f_[ i ] = g.f_[ i ] ;
	}

// TODO
	// Trim leading zero coefficients, but leave a constant term of zero.
//    while( f_.back() == 0 && f_.size() > 0)
//	    f_.pop_back() ;
		
    // Return current object now containing the sum.
	return *this ;
}


const Polynomial operator+( const Polynomial & f, const Polynomial &g )
                         throw( PolynomialRangeError, bad_exception )
{
    // Do + in terms of += to maintain consistency.
	// Copy construct temporary copy, then add to it.
    // Return value optimization compiles away the copy constructor.
    // const on return type disallows doing (u+v) = w ;
	return Polynomial( f ) += g ;
}

// Scalar multiple.
Polynomial & Polynomial::operator*=( const uint k )
             throw( bad_exception ) 
{
    // Multiply coefficients modulo p.
    for (int i = 0 ;  i <= n_ ;  ++i)
        f_[ i ] = mod( f_[ i ] * k ) ;
		
    // Return current object now containing the scalar product.
    return *this ;
}

const Polynomial operator*( const Polynomial & f, const uint k )
                         throw( bad_exception )
{
    // Do * in terms of *= to maintain consistency.
	// Copy construct temporary copy, then add to it.
    // Return value optimization compiles away the copy constructor.
    // const on return type disallows doing (u*k) = w ;
	return Polynomial( f ) *= k ;
}


/*=============================================================================
|
| NAME
|
|     Polynomial operator()
|
| DESCRIPTION
|
|     Evaluate the monic polynomial f( x ) with modulo p arithmetic.
|
|               n         n-1
|     f( x ) = x  +  a   x  + ... + a    0 <= a  < p
|                     n-1            0         i
|
| EXAMPLE
|                                  4
|     Let n = 4, p = 5 and f(x) = x  + 3x + 3.
|
|	 By Horner's rule, f(x) = (((x + 0)x + 0)x + 3)x + 3.
|
|	 Then f(2) = (((2 + 0)2 + 0)2 + 3) = (8 + 3)2 + 3 = 1 + 2 + 3 (mod 5) = 0.
|    and f(3) = (((3 + 0)3 + 0)3 + 3)3 + 3 (mod 5) = 3
|
| METHOD
|
|     By Horner's rule, f(x) = ( ... ( (x + a   )x + ... )x + a .
|                                            n-1               0

|     We evaluate recursively, f := f * x + a (mod p), starting
|                                            i
|     with f = 1 and i = n-1.
|
| BUGS
|
|     None.
|
+============================================================================*/

int
Polynomial::operator()( int x )
                 throw( PolynomialRangeError, bad_exception )
{

int val = 1 ;

for (int degree = n_- 1 ;  degree >= 0 ;  --degree)
    val = mod( val * x + f_[ degree ]) ;

return( val ) ;

} /* ========================= end of function eval ========================= */





/*

DESCRIPTION

     Check if the polynomial f(x) has any linear factors.

INPUT

    f (int *)  nth degree monic mod p polynomial f(x).

    n (int)    Its degree.

    p (int)    Modulus for coefficient arithmetic.

RETURNS

    true    if f( a ) = 0 (mod p) for a = 1, 2, ... p-1.
    false   otherwise

    i.e. check if f(x) contains a linear factor (x - a).  We don't need to test
    for the root a = 0 because f(x) was chosen in main to have a non-zero
    constant term, hence no zero root.

EXAMPLE
                                 4                               2   2
    Let n = 4, p = 5 and f(x) = x  + 3x + 3.  Then f(x) = (x + 3)  (x + 4x + 2)
	
	Then f(0) = 3 (mod 5), f(1) = 2 (mod 5), but
    f(2) = 0 (mod 5), so we exit immediately with a true.
	
	                 4      2
	However, f(x) = x  + 3 x  + x + 1 is irreducible, so has no linear factors.

METHOD

    Evaluate f(x) at x = 0, ..., p-1 by Horner's rule.  Return instantly the
    moment f(x) evaluates to 0.

BUGS

    None.

*/

bool
Polynomial::hasLinearFactor()
                 throw( PolynomialRangeError, bad_exception )
{

for (int i = 0 ;  i <= p_ - 1 ;  ++i)

    if ((*this)( i ) == 0)

        return( true ) ;

return( false ) ;

}


/*

DESCRIPTION

     Test if a polynomial is a constant.

INPUT

    t (int *)  mod p polynomial t(x).
    n (int)    Its degree.

RETURNS

    true    if t(x) is a constant polynomial.
    false   otherwise

EXAMPLE
                               2
    Let n = 3.  Then t(x) = 2 x  is not constant but t(x) = 2 is.

METHOD

    A constant polynomial is zero in its first through n th degree
    terms.  Return immediately with false if any such term is non-zero.

BUGS

    None.
*/

bool
Polynomial::isInteger() const
                 throw( PolynomialRangeError, bad_exception )
{
// Degree 0 is constant.
if (n_ == 0)
    return true ;

// Not integer if any coefficients above zero degree term are non-zero.
for (int i = 1 ;  i <= n_ ;  ++i)
    if (f_[ i ] != 0)
        return( false ) ;

return( true ) ;

} /* ====================== end of function is_integer ====================== */



/*

DESCRIPTION

    Return the initial monic polynomial in the sequence for f(x).

INPUT

     f (int *)                Monic polynomial f(x).
     n      (int, n >= 1)     Degree of f(x).

RETURNS
                                             n
     f (int *)                Sets f( x ) = x  - 1

  EXAMPLE
                             4
     Let n = 4.  Set f(x) = x  - 1.

METHOD

BUGS

    None.
*/

void Polynomial::initial_trial_poly( const uint n, const uint p )
                 throw( PolynomialRangeError, bad_exception )
{
// Set the modulus.
(*this).setModulus(p);

(*this)[ n ] = 1 ;
(*this)[ 0 ] = -1 ;

} /* ================= end of function initial_trial_poly ==================== */



/*

DESCRIPTION

    Return the next monic polynomial in the sequence after f(x).

INPUT

    f (int *)           Monic polynomial f(x).
    n (int, n >= 1)     Degree of monic polynomial f(x).
    p (int, p >= 2)     Modulo p coefficient arithmetic.

RETURNS

     f (int *)                Overwrites f(x) with the next polynomial
                              after it in the sequence (explained below).
EXAMPLE
                                           3
     Let n = 3 and p = 5.  Suppose f(x) = x  + 4 x + 4.  As a mod p number,
     this is 1 0 4 4.  Adding 1 gives 1 0 4 5.  We reduce modulo
     5 and propagate the carry to get the number 1 0 5 0.  Propagating
     the carry again and reducing gives 1 1 0 0.  The next polynomial
                     3    2
     after f(x) is  x  + x .

METHOD

     Think of the polynomial coefficients as the digits of a number written
     in base p.  The "next" polynomial is the one you would get by adding 1
     to this number in multiple precision arithmetic.  Our intention is to
     run through all possible monic polynomials modulo p.

     Propagate carries in digits 1 through n-2 when any digit exceeds p.  No
     carries take place in the n-1 st digit because our polynomial is monic.

BUGS

    None.

*/
void Polynomial::next_trial_poly()
                 throw( PolynomialRangeError, bad_exception )
{
++f_[ 0 ] ;     /* Add 1, i.e. increment the coefficient of the x term. */

/*
     Sweep through the number from right to left, propagating carries.  Skip
     the constant and the nth degree terms.
*/

for (int digit_num = 0 ;  digit_num <= n_ - 2 ;  ++digit_num)
{
    if (f_[ digit_num ] == p_)   /*  Propagate carry to next digit. */
    {
        f_[ digit_num ] = 0 ;
        ++f_[ digit_num + 1 ] ;
    }
}

} /* ================= end of function next_trial_poly ====================== */






/*------------------------------------------------------------------------------
|                              PolyMod Implementation                          |
------------------------------------------------------------------------------*/

PolyMod::PolyMod()
           throw( PolynomialRangeError, bad_exception )
           : g_()
           , f_()
           , powerTable_()
           , mod( f_.modulus() )
{
    constructPowerTable() ;
	modf() ;
}

PolyMod::~PolyMod()
{
// Member fields will clean up themselves.
}

PolyMod::PolyMod( const string & g, const Polynomial & f )
         throw( PolynomialRangeError, bad_exception )
         : g_( g )
         , f_( f )
         , powerTable_()
         , mod( f.modulus() )
{
    constructPowerTable() ;
	modf() ;
}

PolyMod::PolyMod( const Polynomial & g, const Polynomial & f )
         throw( PolynomialRangeError, bad_exception )
         : g_( g )
         , f_( f )
         , powerTable_()
         , mod( f.modulus() )
{
    constructPowerTable() ;
	modf() ;
}
     
// Operator casting to string type.
PolyMod::operator string() const
                 throw( PolynomialRangeError, bad_exception )
{
return static_cast<string>( g_ ) ;
}

ostream & operator<<( ostream & out, const PolyMod & p )
                  throw( PolynomialRangeError, bad_exception )
{
    // Cast to polynomial to a string first, then output as a string.  
    // May throw a PolynomialRangeError.
    out << static_cast<string>( p.g_ ) << endl ;

	return out ;
}


const Polynomial PolyMod::getf() const
{
    return f_ ;
}
		
const int PolyMod::getModulus() const
{
    return f_.modulus() ;
}

PolyMod::PolyMod( const PolyMod & g2 )
         throw( PolynomialRangeError, bad_exception )
         : g_( g2.g_ )
         , f_( g2.f_ )
         , powerTable_( g2.powerTable_ )
         , mod( f_.modulus() )
{
}
                  

PolyMod & PolyMod::operator=( const PolyMod & g2 )
                   throw( PolynomialRangeError, bad_exception )
{
	// If assigned to itself, return a reference to the object itself,
    // without doing anything further.
    if (this == &g2)
        return *this ;

    g_ = g2.g_ ;
    g_ = g2.f_ ;

    powerTable_ = g2.powerTable_ ;
    mod = g2.mod ;

    // Return a reference to the altered object.
    return *this ;
}

// Bounds checked indexing operator for read only access:
// coeff = p[ i ] ;
const int PolyMod::operator[]( int i ) const
                          throw( PolynomialRangeError, bad_exception )
{
// Can throw an exception.                          
return g_[ i ] ;
}
                                
                                
                                
/*=============================================================================
|
| NAME
|
|     constructPowerTable
|
| DESCRIPTION
|
|     Construct a table of powers of x:
|
|      n                     2n-2
|     x  (mod f(x), p)  ... x    (mod f(x), p)
|
| SYNTAX
|
|    f (int *)   Coefficients of f(x), a monic polynomial of degree n.
|    n (int, -infinity < n < infinity)
|    p (int, p > 0)
|
|    powerTable_ (int *)  powerTable_[i][j] is the coefficient of
|     j       n+i
|    x   in  x   (mod f(x), p) where 0 <= i <= n-2 and 0 <= j <= n-1.
|
| EXAMPLE
|                                  4     2                     4
|     Let n = 4, p = 5 and f(x) = x  +  x  +  2x  +  3.  Then x  =
|
|         2                  2
|     -( x  +  2x  + 3) = 4 x  + 3 x + 2 (mod f(x), 5), and we get
|
|      4                    2
|     x  (mod f(x), 5) = 4 x  + 3 x + 2 = powerTable_[0].
|
|      5                       2                 3      2
|     x  (mod f(x), 5) = x( 4 x  + 3 x + 2) = 4 x  + 3 x  + 2x
|                      = powerTable_[1].
|
|      6                       3      2           4      3      2
|     x  (mod f(x), 5) = x( 4 x  + 3 x + 2 x) = 4x  + 3 x  + 2 x
|
|                              2                 3      2
|                      = 4 ( 4x  + 3 x + 2) + 3 x  + 2 x  =
|
|                           3     2
|                      = 3 x + 3 x + 2 x + 3 = powerTable_[2].
|
|                                    j
|     powerTable_[i][j]:       | 0  1  2  3
|                           ---+-------------
|                            0 | 2  3  4  0
|                        i   1 | 0  2  3  4
|                            2 | 3  2  3  3
|
| NOTES
|                              n-1
|     Beginning with t( x ) = x,    compute the next power of x from the last
|                                                         n
|     one by multiplying by x.  If necessary, substitute x  =
|             n-1
|     -( a   x    + ... + a ) to reduce the degree.  This formula comes from
|         n-1              0
|                               n         n-1
|     the observation f( x ) = x   + a   x    + ... + a    = 0 (mod f(x), p).
|                                     n-1              0
| BUGS
|
|     None.
|
+============================================================================*/

void
PolyMod::constructPowerTable()
{
// Get hold of the degree of f(x).
int n = f_.deg() ;

// Empty the power table.
powerTable_.clear() ;

//
//  t(x) is temporary storage for x ^ k (mod f(x),p)
//   n <= k <= 2n-2.  Its degree can go as high as
//   n before it is reduced again.
Polynomial t ;


//                         n-1
//    Initialize t( x ) = x    mod p.
t[ n-1 ] = 1 ;
t.setModulus( getModulus() ) ;



//                                      i+n
//  Fill the ith row of the table with x   (mod f(x), p)
//  for i = 0 ... n-2.
//
for (int i = 0 ;  i <= n - 2 ;  ++i)
{
    // Compute t(x) = x t(x) by shifting the coefficients
    // to the left and filling with zero.
    for (int j = n ;  j >= 1 ;  --j)
        t[ j ] = t[ j-1 ] ;

    t[ 0 ] = 0 ;

	//  Coefficient of the x ^ n degree term of t(x).  
	int coeff = 0 ;
    if ( (coeff = t[ n ]) != 0)
    {
        //  Zero out the x ^ n th term. 
        t[ n ] = 0 ;

        //          n       n                        n-1
        // Replace x  with x  (mod f(x), p) = -(a   x   + ... + a )
	    //	                                     n-1             0
        for (int j = 0 ;  j <= n-1 ;  ++j)

            t[ j ] = mod( t[ j ] +
                          mod( -coeff * f_[ j ]) ) ;
    }  // end if 

    //  Copy t(x) into row i of power_table. 
    powerTable_.push_back( t ) ;

} // end for 

#ifdef DEBUG_PP_POLYNOMIAL
    cout << "PowerTable of polynomials x^n ... x^2n-2 mod f(x), p" << endl ;
	cout << "f(x) = " << getf() << " n = " << n << " p = " << getModulus() << endl ;
    for  (int i = n ;  i <= 2*n-2 ;  ++i)
        cout << "powerTable[ x^" << i << " ] = " << powerTable_[ offset(i) ] << endl ;
#endif


// TODO check size of powerTable

// t will be automagically freed upon exit.
return ;

} /* ================== end of function constructPowerTable =============== */


/*=============================================================================
|
| NAME
|
|     modf
|
| DESCRIPTION
|
|     Reduce g(x) modulo f(x), and p.
|
| SYNTAX
|

| EXAMPLE
|                                  4     2                   
|     Let n = 4, p = 5 and f(x) = x  +  x  +  2x  +  3.  
|
|      6                       3      2           4      3      2
|     x  (mod f(x), 5) = x( 4 x  + 3 x + 2 x) = 4x  + 3 x  + 2 x
|
|
| NOTES
|
|
| BUGS
|
|     None.
|
+============================================================================*/

void
PolyMod::modf() throw( PolynomialRangeError, bad_exception )
{
// Get hold of the degree of f(x).
int n = f_.deg() ;
int m = g_.deg() ;

if (m > 2 * n - 2)
{
    ostringstream os ;
    os << "Error in PolyMod::modf():  degree of g(x) higher than power table can handle "
       << "with deg f = " << n << " deg g = " << m << " and p = " << getModulus() 
       << " at " << __FILE__ << ": line " << __LINE__ ;
	throw PolynomialRangeError( os.str() ) ;
}


//                                      i+n
//  Fill the ith row of the table with x   (mod f(x), p)
//  for i = 0 ... n-2.
//
for (int i = n ;  i <= m ;  ++i)
{
    #ifdef DEBUG_PP_POLYNOMIAL
     cout << "\nBefore converting, g( x ) = " << g_ << endl ;
    #endif
	
	//  Coefficient of the x ^ i degree term of g(x).  
	int coeff = 0 ;
    if ( (coeff = g_[ i ]) != 0)
    {
        //  Subtract (zero out) the x ^ i th term. 
        g_[ i ] = 0 ;

        //          i       i
        // Replace x  with x  (mod f(x), p) from the power table * coeff.
		g_ += (powerTable_[ offset(i) ] * coeff) ;
	 }
	 
     #ifdef DEBUG_PP_POLYNOMIAL
     cout << "\nAfter converting with coeff = " << coeff << " g( x ) = " << g_ << endl ;
     #endif

} // end for 

return ;

} /* ================== end of function constructPowerTable =============== */



/*==============================================================================
|                                 autoconvolve                                 |
================================================================================

DESCRIPTION

     Compute a convolution-like sum for use in function coeffOfSquare.

INPUT
                                             n-1
    t (int *)   Coefficients of t(x) = t    x    + ... + t x  + t
                                        n-1               1      0
    k (int), 0 <= k <= 2n - 2)

    lower (int, 0 <= lower <= n-1)

    uppper (int, 0 <= upper <= n-1)

    p (int, p > 0)

RETURNS

     upper
     ---
     \    t  t       But define the sum to be 0 when lower > upper to catch
     /     i  k-i    the special cases that come up in function coeffOfSquare.
     ---
     i=lower

EXAMPLE
                       3      2
     Suppose t(x) = 4 x  +  x  +  3 x  +  3, lower = 1, upper = 3, n = 3,

     and p = 5.  For k = 3, autoConvolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +

     t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3.  For lower = 0,

     upper = -1, or for lower = 3 and upper = 2, autoConvolve = 0, no matter what

     k is.

METHOD

    A "for" loop in the C language is not executed when its lower limit

    exceeds its upper limit, leaving sum = 0.

BUGS

     None.

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

int autoConvolve( const Polynomial & t, int k, int lower, int upper )
{
    ModP mod( t.modulus() ) ;
    int deg_t = t.deg() ;

    int sum = 0 ;
    for (int i = lower ;  i <= upper ;  ++i)
    {
        // Coeff is zero if higher or lower than degree of polynomial.
        uint coeff_ti   = 0u ;
        uint coeff_tkmi = 0u ;

        if (i <= deg_t && i >= 0)
            coeff_ti = t[ i ] ;

        if (k-i <= deg_t && k-i >= 0)
            coeff_tkmi = t[ k - i ] ;

        sum = mod( sum + mod( coeff_ti * coeff_tkmi )) ;
    }

    return( sum ) ;

} /* ==================== end of function autoConvolve ===================== */


/*==============================================================================
|                                 convolve                                     |
================================================================================

DESCRIPTION

     Compute a convolution-like sum.

INPUT
                                             n-1
    s (int *)   Coefficients of s(x) = s    x    + ... + s x  + s
                                        n-1               1      0
										     n-1
    t (int *)   Coefficients of t(x) = t    x    + ... + t x  + t
                                        n-1               1      0

    k (int), 0 <= k <= 2n - 2)

    lower (int, 0 <= lower <= n-1)

    uppper (int, 0 <= upper <= n-1)

    p (int, p > 0)

RETURNS

     upper
     ---
     \    s  t       But define the sum to be 0 when lower > upper to catch
     /     i  k-i    the special cases
     ---
     i=lower

EXAMPLE
                       3      2
     Suppose s(x) = 4 x  +  x  +  3 x  +  3,

	                   3     2
     Suppose t(x) = 4 x  +  x  +  3 x  +  3,


	 lower = 1, upper = 3, n = 3,

     and p = 5.  For k = 3, convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +

     t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3.  For lower = 0,

     upper = -1, or for lower = 3 and upper = 2, convolve = 0, no matter what

     k is.

METHOD

    A "for" loop in the C language is not executed when its lower limit

    exceeds its upper limit, leaving sum = 0.

BUGS

     None.

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

int convolve( const Polynomial & s, const Polynomial & t,
              const int k, const int lower, const int upper )
{

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

int
    sum = 0,      /* Convolution sum. */
    i ;           /* Loop counter.    */

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

ModP mod( s.modulus() ) ; // Initialize the functionoid.

int deg_s = s.deg() ;
int deg_t = t.deg() ;

for (i = lower ;  i <= upper ;  ++i)
{
    // Coeff is zero if higher or lower than degree of polynomial.
    uint coeff_s = 0u ;
    uint coeff_t = 0u ;

    if (i <= deg_s && i >= 0)
        coeff_s = s[ i ] ;

    if (k-i <= deg_t && k-i >= 0)
        coeff_t = t[ k- i ] ;

   sum = mod( sum + mod( coeff_s * coeff_t )) ;
}

return( sum ) ;

} /* ======================= end of function convolve ======================= */


/*

DESCRIPTION
                                    th                2
     Return the coefficient of the k   power of x in g ( x )  modulo p.

INPUT

    Coefficients of g(x), of degree <= n-1.

    k  0 <= k <= 2n-2

RETURNS
                        th                 2
    coefficient of the k   power of x in g (x) mod p.

EXAMPLE
                                    3     2                 2
    Let n = 4, p = 5, and g(x) = 4 x  +  x  + 3 x + 3.  g(x) =

     6      5
    x  + 3 x + 3 x + 3 x + 4, all modulo 5.

            k        |  0  1  2  3  4  5  6
     ----------------+---------------------
     coeffOfSquare   |  4  3  0  0  0  3  1

METHOD
                                                         2
    The formulas were gotten by writing out the product g (x) explicitly.

    The sum is 0 in two cases:

	    (1) when k = 0 and the limits of summation are 0 to -1

        (2) k = 2n - 2, when the limits of summation are n to n-1.

    To derive the formulas, let

                     n-1
    Let g(x) = g    x     +  ... + g x + g
                n-1                 1     0

    Look at the formulas in coeffOfProduct for each power of x,
	replacing s with t, and observe that half of the terms are
	duplicates, so we can save computation time.

    Inspection yields the formulas,

    for 0 <= k <= n-1, even k,

     k/2-1
      ---             2
   2  \   g  g     + g
	  /    i  k-i     k/2
 	  ---
	  i=0

    for 0 <= k <= n-1, odd k,

      (k-1)/2
      ---
    2 \   g  g
	  /    i  k-i
	  ---
	  i=0


    and for n <= k <= 2n-2, even k,

      n-1
      ---            2
   2  \   g  g    + g
	  /    i  k-i    k/2
	  ---
	  i=k/2+1

      and for n <= k <= 2n-2, odd k,

      n-1
      ---
   2  \   g  g
	  /    i  k-i
	  ---
	  i=(k+1)/2


BUGS

     None.
*/

int coeffOfSquare( const Polynomial & g, const int k, const int n )
{
ModP mod( g.modulus() ) ; // Initialize the functionoid.

                   //                          2
int sum = 0 ;      // kth coefficient of g( x )

// Coeff is zero if higher or lower than degree of polynomial.
uint coeff_gkd2 = 0 ;
if (k/2 <= g.deg() && k/2 >= 0)
    coeff_gkd2 = g[ k/2 ] * g[ k/2 ] ;

if (0 <= k && k <= n-1)
{
    if (k % 2 == 0)        // Even k 
        sum = mod( mod( 2 * autoConvolve( g, k, 0, k/2 - 1) ) + coeff_gkd2 ) ;

     else                  // Odd k 
         sum = mod( 2 * autoConvolve( g, k, 0, (k-1)/2)) ;
}
else if (n <= k && k <= 2 * n - 2)
{

    if (k % 2 == 0)        // Even k 
        sum = mod( mod( 2 * autoConvolve( g, k, k/2 + 1, n-1)) + coeff_gkd2 ) ;

     else                  // Odd k 
         sum = mod( 2 * autoConvolve( g, k, (k+1)/2, n-1)) ;
}

return( sum ) ;

} /* ================== end of function coeffOfSquare ===================== */


/*==============================================================================
|                               coeffOfProduct                               |
================================================================================

DESCRIPTION
                                    th
     Return the coefficient of the k   power of x in s( x ) t( x )  modulo p.

EXAMPLE
                              3     2                  2
  Let n = 4, p = 5, t(x) = 4 x  +  x  + 4, s( x ) = 3 x  + x + 2

                           5      4      3      2
  then s ( x ) t( x ) = 2 x  + 2 x  + 4 x  + 4 x  + 4 x + 3

  We'll do the case k=3,

  t3 s0 + t2 s1 + t1 s2 + t0 s3 = 4 * 2 + 1 * 1 + 0 * 3 + 4 * 0 = 9 = 4 (mod 5).

            k       |  0  1  2  3  4  5  6
     -----------------+---------------------
     coeffOfProduct |  3  4  4  4  2  2  0

METHOD

    The formulas were gotten by writing out the product s(x) t (x) explicitly.

    The sum is 0 in two cases:

	    (1) when k = 0 and the limits of summation are 0 to -1

        (2) k = 2n - 2, when the limits of summation are n to n-1.


    To derive the formulas, let

                      n-1
    Let s (x) = s    x     +  ... + s x + s  and
                 n-1                 1     0

                      n-1
        t (x) = t    x     +  ... + t x + t
                 n-1                 1     0

    and multiply out the terms, collecting like powers of x:


    Power of x     Coefficient
	==========================
     2n-2
    x              s    t
                    n-1  n-1

     2n-3
    x              s    t    +  s    t
                    n-2  n-1     n-1  n-2

     2n-4
    x              s    t    +  s    t    +  s    t
                    n-3  n-1     n-2  n-2     n-3  n-1

     2n-5
    x              s    t    +  s    t    +  s    t    +  s    t
                    n-4  n-1     n-3  n-2     n-2  n-3     n-1  n-4

     . . .

     n
    x              s  t    +  s  t    + ...  +  s    t
                    1  n-1     2  n-2            n-1  1

     n-1
    x              s  t    +  s  t    + ...  +  s    t
                    0  n-1     1  n-2            n-1  0

    . . .

     3
    x              s  t  +  s  t  +  s  t  +  s  t
                    0  3     1  2     2  1     3  0

     2
    x              s  t  +  s  t  +  s  t
                    0  2     1  1     2  0


    x              s  t  +  s  t
                    0  1     1  0

    1              s  t
                    0  0


    Inspection yields the formulas,


    for 0 <= k <= n-1,

      k
     ---
     \   s  t
	 /    i  k-i
	 ---
	 i=0


    and for n <= k <= 2n-2,

     n-1
     ---
     \   s  t
	 /    i  k-i
	 ---
	 i=k-n+1

BUGS

     None.

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

int coeffOfProduct( const Polynomial & s, const Polynomial & t, const int k, const int n )
{
// TODO check p is the same for s and t.
// TODO check degree of s and t < n.
    
int sum = 0 ;      /* kth coefficient of t(x) ^ 2. */

if (0 <= k && k <= n-1)
{
    sum = convolve( s, t, k, 0, k ) ;
}
else if (n <= k && k <= 2 * n - 2)
{
    sum = convolve( s, t, k, k - n + 1, n - 1 ) ;
}

return( sum ) ;

} 


/*==============================================================================
|                                  product                                     |
================================================================================

DESCRIPTION

     Compute s( x ) t( x ) (mod f(x), p)

     Uses a precomputed table of powers of x.

INPUT

    s (int *)              Coefficients of s(x), of degree <= n-1.

    t (int *)              Coefficients of t(x), of degree <= n-1.

    powerTable (int **)    x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.

    n (int, -infinity < n < infinity)  Degree of f(x).

    p (int, p > 0)         Mod p coefficient arithmetic.

OUTPUT

    s (int *)              Overwritten with s( x ) t( x ) (mod f(x), p)

EXAMPLE
                                     3    2                 2
    Let n = 4 and p = 5, t( x ) = 4 x  + x + 4, s( x ) = 3 x + x + 2

                            5      4     3     2
	Then s( x ) t( x ) = 2 x  + 2 x + 4 x + 4 x + 4 x + 3, modulo 5,

                                         4   2
    and after reduction modulo f( x ) = x + x + 2 x + 3, using the power

	                   4      2               5      3      2
	table entries for x  = 4 x + 3 x + 2 and x  = 4 x  + 3 x + 2 x, we get

	                                   3      2
    s( x ) t( x ) (mod f( x ), p) = 2 x  + 3 x  + 4 x + 2


METHOD

    Compute the coefficients using the function coeffOfProduct.

    The next step is to reduce s(x) t(x) modulo f(x) and p.  To do so, replace

                           k                                      k
    each non-zero term t  x,  n <= k <= 2n-2, by the term t * [ x   mod f(x), p)]
                        k                                  k

    which we get from the array powerTable.

BUGS

     None.

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

const PolyMod operator*( const PolyMod & s,
                         const PolyMod & t )
                         throw( PolynomialRangeError, bad_exception )
{
    // Do * in terms of *= to maintain consistency.
    // Return value optimization compiles away the copy constructor.
    // const on return type disallows doing (u*v) = w ;
    return PolyMod( s ) *= t ;
}


PolyMod &
PolyMod::operator*=( const PolyMod & t )
                   throw( PolynomialRangeError, bad_exception )
{
int
    i, j,    /* Loop counters. */
    coeff;   /* Coefficient of x ^ k term of t(x) ^2 */

// Temporary storage for the new t(x).  Can have degree up to n.
Polynomial temp ;

// Get hold of the degree of f(x).
int n = f_.deg() ;

/*
                                 0        n-1
    Compute the coefficients of x , ..., x.   These terms do not require

    reduction mod f(x) because their degree is less than n.
*/

for (i = 0 ;  i <= n ;  ++i)
    temp[ i ] = coeffOfProduct( g_, t.g_, i, n ) ;

/*
                                 n        2n-2             k
    Compute the coefficients of x , ..., x.    Replace t  x  with
                                                        k
            k
    t  * [ x  (mod f(x), p) ] from array powerTable when t is
     k                                                    k
    non-zero.
*/

for (i = n ;  i <= 2 * n - 2 ;  ++i)

    if ( (coeff = coeffOfProduct( g_, t.g_, i, n)) != 0 )

        for (j = 0 ;  j <= n - 1 ;  ++j)

            temp[ j ] = mod( temp[ j ] +
                             mod( coeff * powerTable_[ offset(i) ] [ j ])) ;

for (i = 0 ;  i <= n - 1 ;  ++i)

    g_[ i ] = temp[ i ] ;

// Return (reference to) the product.
return *this ;

} /* ======================== end of function product ======================= */




/*

DESCRIPTION

     Compute x g(x) (mod f(x), p)

SYNTAX

    g.timesX( t ) ;

EXAMPLE
                                    3       2
    Let n = 4, p = 5, and g(x) = 2 x  +  4 x  + 3 x.  Let f(x) =
     4    2                                  4      3      2
    x  + x  + 2 x + 3.  Then x t (x) = 2 x  + 4 x  + 3 x  and
                                     2                3     2
    x g(x) (mod f(x), p) = 2 * (4 x + 3 x + 2) + 4 x + 3 x  =
       3    2
    4 x  + x + x + 4.
         3     2
    = 4 x + 2 x + 3 x + 2.

METHOD

	Uses a precomputed table of powers of x.

                          n-1         n-2
    Multiply g(x) = g    x   +  g    x   + ... + g  by shifting the coefficients
                     n-1         n-2              0

                         n
    to the left.  If an x   term appears, eliminate it by

    substitution using powerTable.

BUGS

     None.

*/

void
PolyMod::timesX()
                 throw( PolynomialRangeError, bad_exception )
{
int n = f_.deg() ;

#ifdef DEBUG_PP_POLYNOMIAL
cout << "timesX:  g( x ) = " << g_ << endl ;
#endif

//  Multiply g(x) by x by shifting the coefficients left in the array, giving
//         n-1
//   g    x    + ... + g  x + 0
//    n-2               1
//                          
// but save the coefficient g    first before overwriting it.
//                           n-1 
int g_coeff = g_[ n - 1 ] ;

for (int i = n-1 ;  i >= 1 ;  --i)

    g_[ i ] = g_[ i-1 ] ;

g_[ 0 ] = 0 ;



//                 n                n
//    Replace g   x  with g    * [ x  (mod f(x), p) ] using
//             n-1         n-1
//     n
//    x  from powerTable

if (g_coeff != 0)
{
    for (int i = 0 ;  i <= n - 1 ;  ++i)

        g_[ i ] = mod( g_[ i ] +
                      mod( g_coeff * powerTable_[ offset(n) ] [ i ] )) ;
}

#ifdef DEBUG_PP_POLYNOMIAL
     cout << "timesX:  x g( x ) = " << g_ << endl ;
#endif


} /* ======================== end of function timesX ======================= */



/*==============================================================================
|                                  square                                      |
================================================================================

DESCRIPTION

              2
     Compute g (x) (mod f(x), p)

SYNTAX

    g.square() ;


EXAMPLE
                                    3     2
    Let n = 4, p = 5, and g(x) = 4 x  +  x  + 4.  Let f(x) =

     4    2                   2       6      5     4      3     2
    x  + x  + 2 x + 3.  Then t (x) = x  + 3 x  +  x  + 2 x + 3 x +  1

    Now subsituting powers of x modulo f(x) from the power table,

      2                         3     2
     t (x) (mod f(x), p) =  (3 x + 3 x + 2 x + 3) +

            3      2                 2                3     2
    3 * (4 x  + 3 x + 2 x) + 4 * (4 x + 3 x + 2) + 4 x + 4 x + 3 x + 1

         3     2
    = 2 x + 4 x + x + 1.


METHOD

    Uses a precomputed table of powers of x.


         2            2n-2              n         n-1
    Let g (x) = g    x     +  ... + g  x  +  g   x   +  ... + g .
                 2n-2                n        n-1              0

    Compute the coefficients g  using the function coeffOfSquare.
                              k

                                2
    The next step is to reduce g (x) modulo f(x).  To do so, replace

                           k                                      k
    each non-zero term g  x,  n <= k <= 2n-2, by the term g * [ x   mod f(x), p)]
                        k                                  k

    which we get from the array powerTable_.

BUGS

     None.

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

void
PolyMod::square()
                   throw( PolynomialRangeError, bad_exception )
{

// Get hold of the degree of f(x).
int n = f_.deg() ;

#ifdef DEBUG_PP_POLYNOMIAL
cout << "square:  g( x ) = " << g_ << endl ;
#endif


// TODO exceptions?
// Temporary storage for the new g(x).  Can have degree up to n.
Polynomial t ;

/*
                                 0        n-1
    Compute the coefficients of x , ..., x.   These terms do not require

    reduction mod f(x) because their degree is less than n.
*/

for (int i = 0 ;  i <= n ;  ++i)

    t[ i ] = coeffOfSquare( g_, i, n ) ;

/*
                                 n        2n-2             k
    Compute the coefficients of x , ..., x.    Replace g  x  with
            k                                           k
			
    g  * [ x  (mod f(x), p) ] from array powerTable_ when g is
     k                                                     k
    non-zero.

*/

for (int i = n ;  i <= 2 * n - 2 ;  ++i)
{
    uint coeff = 0 ;
	
    if ( (coeff = coeffOfSquare( g_, i, n )) != 0 )

        for (int j = 0 ;  j <= n- 1 ;  ++j)

            t[ j ] = mod( t[ j ] +
                             mod( coeff * powerTable_[ offset(i) ] [ j ])) ;
}

for (int i = 0 ;  i <= n - 1 ;  ++i)

    g_[ i ] = t[ i ] ;

#ifdef DEBUG_PP_POLYNOMIAL
cout << "square:  g( x ) ^ 2 = " << g_ << endl ;
#endif


} /* ======================== end of function square ======================== */


/*==============================================================================
|                                  power                                       |
================================================================================

DESCRIPTION
                      m
     Compute g(x) = x  (mod f(x), p).

EXAMPLE
                              4    2
    Let n = 4, p = 5, f(x) = x  + x  + 2 x + 3, and m = 156.

    156 = 0  . . . 0  1  0  0  1  1  1  0  0 (binary representation)
          |<- ignore ->| S  S SX SX SX  S  S (exponentiation rule,
		                                      S = square, X = multiply by x)
     m
    x  (mod f(x), p) =

         2     2
6   S   x  =  x

         4       2
5   S   x  =  4 x  + 3 x + 2


         8       3      2
4   S   x  =  4 x  + 4 x + 1

         9       3     2
4   X   x  =  4 x  +  x + 3 x + 3


         18       2
3   S   x  =  2 x  +  x + 2

         19      3     2
3   X   x  =  2 x  +  x  + 2 x


         38      3       2
2   S   x  =  2 x  +  4 x  +  3 x

         39      3       2
2   X   x  =  4 x  +  2 x  +  x + 4


         78      3       2
1   S   x  =  4 x  +  2 x  +  3 x + 2

         156
0   S   x    = 3

METHOD

    Exponentiation by repeated squaring, using precomputed table of
	powers.  See ART OF COMPUTER PROGRAMMING, vol. 2, 2nd Ed.,
	D. E. Knuth,  pgs 441-443.

     n         2n-2
    x,  ... , x    (mod f(x), p)

                    m
    to find g(x) = x   (mod f(x), p), expand m into binary,

           k        k-1
    m = a 2  + a   2    + . . . + a 2 + a
         k      k-1                1     0

                            m
	where a = 1, and split x   apart into
           k

            k      k
           2      2  a             2 a    a
     m                k-1             1    0
    x  =  x     x           . . .  x     x


    Then to raise x to the mth power, do


    t( x ) = x

    return if m = 1


    for i = k-1 downto 0 do

                   2
        t(x) = t(x)  (mod f(x), p)       // Square each time.

        if a = 1 then
		    i

            t(x) = x t(x) (mod f(x), p)  // Times x only if current bit is 1
        endif

    endfor
                                                        k
                                                       2
    The initial t(x) = x gets squared k times to give x  .  If a  = 1 for
	                                                            i
    0 <= i <= k-1, we multiply by x which then gets squared i times more

              i
             2
    to give x .  On a binary computer, we use bit shifting and masking to

    identify the k bits  { a    . . .  a  } to the right of the leading 1
	                        k-1         0

    bit.  There are log ( m ) - 1 squarings and (number of 1 bits) - 1
		                   2
	multiplies.

BUGS

     None.

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

const PolyMod power( const PolyMod & g1, const BigInt & m )
            throw( PolynomialRangeError, bad_exception )
{
// Return if g(x) != x
if (g1.f_.deg() == 1 && g1[ 0 ] == 0 && g1[ 1 ] == 1)
{
    ostringstream os ;
    os << "Error in PolyMod::power():  g( x ) != x "
       << "with deg g = " << g1.f_.deg() << " m = " << m  
       << " at " << __FILE__ << ": line " << __LINE__ ;
	throw PolynomialRangeError( os.str() ) ;
}

// Exit right away if m = 1 and return a copy of g(x).
PolyMod g( g1 ) ;

if (m == 1u)
    return g ;

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

#ifdef DEBUG_PP_POLYNOMIAL
cout << "initial max bitNum = " << bitNum << endl ;
cout << "g( x ) = " << g << endl ;
#endif 

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

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

if (bitNum == -1)
{
    ostringstream os ;
    os << "PolyMod::x_to_power " << "bitNum == -1 internal error in PolyMod"
       << " at " << __FILE__ << ": line " << __LINE__ ;
    throw PolynomialRangeError( os.str() ) ;
}

#ifdef DEBUG_PP_POLYNOMIAL
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 x for
    every 1 bit.
*/

while ( --bitNum >= 0 )
{
    g.square() ;

    if (m.testBit( bitNum ))
       g.timesX() ;

    #ifdef DEBUG_PP_POLYNOMIAL
    cout << "S " ;
    if (m.testBit( bitNum ))
        cout << "X " ;
    cout << "Bit num = " << bitNum << " g( x ) = " << g << endl ;
    #endif
}

#ifdef DEBUG_PP_POLYNOMIAL
    cout << "Out of the loop bitNum = " << bitNum << " g( x ) = " << g << endl ;
#endif

return g ;
 
}


bool PolyMod::isInteger() const
{
    return g_.isInteger() ;
}
        
// Set a new value of f(x) with same degree n and modulus p.
void PolyOrder::newPolynomial( const Polynomial & f )
                       throw( PolynomialRangeError, bad_exception )
{
    f_ = f ;
}

PolyOrder::PolyOrder( const Polynomial & f )
                 throw( PolynomialRangeError, bad_exception )
             : f_( f )
             , r_( 0 )
             , a_( 0 )
			 , factor_( 1 )
             , Q_( 0 )
			 , p_( f.modulus() )
			 , n_( f.deg() )
			 , mod( f.modulus() )

{
    //                                          n
    //	                                       p  - 1 
	//  Get everything ready by computing r = --------
	//                                         p - 1
	// and factoring it into the product of primes.
    BigInt pn = power( p_, n_ ) ;
    r_ = (pn - 1u) / (p_ - 1u) ;

    Factor<BigInt> fact( r_ ) ;
	factor_ = fact ;
	
	// TODO what if we can't allocate Q?
	// Prepare the Q matrix to the proper size.
    Q_.clear() ;
    Q_.resize( n_ ) ;

    for (int row = 0 ;  row < n_ ;  ++row)
    {
        Q_[ row ].resize( n_ ) ;
    }
}
             
   

/*
DESCRIPTION
                 m
    Check that x  (mod f(x), p) is not an integer for m = r / p  but skip
                                                               i
                                           n
                                          p  - 1           th
    this test if p  | (p-1).  Recall r = -------, and p = i   prime in 
				  i                       p - 1        i
				  
    the factorization of r.


EXAMPLE
                                           2
     Let n = 4 and p = 5.  Then r = 156 = 2 * 3 * 13, and p = 2, p = 3,
														   1      2
														   
     and p = 13.  m = { 156/2, 156/3, 156/13 } = { 78, 52, 12 }.  We can
          3
		  
     skip the test for m = 78 because p = 2 divides p-1 = 4.  Exponentiation
                                       1
									   
            52       3   2                                    12
     gives x    = 2 x + x + 4 x, which is not an integer and x   =
	 
        3       2
     4 x  +  2 x  +  4 x  + 3 which is not an integer either, so we return

     true.

METHOD

    Exponentiate x with x_to_power and test the result with is_integer.
    Return right away if the result is not an integer.

BUGS

    None.
*/

bool PolyOrder::order_m()
                throw( PolynomialRangeError, bad_exception )
{
int p = f_.modulus() ;
                
for (int i = 0 ;  i < factor_.num() ;  ++i)
{
    if (!factor_.skip_test( p, i ))
    {
        BigInt m = r_ / factor_[ i ] ;

	    Polynomial x1( "x" ) ;
		x1.setModulus( p ) ;
        PolyMod x( x1, f_ ) ;
		
        PolyMod x_to_m = power( x, m ) ;
		
        #ifdef DEBUG_PP_POLYNOMIAL
        cout << "Prime factor p[ " << i << " ] = " << factor_[ i ] << endl ;
        cout << "m = " << m << endl ;
		cout << "x^m = " << x_to_m << endl ;
        #endif 

        if (x_to_m.isInteger())
		{
 			return( false );
		}
    }
}

return( true ) ;

} /* ====================== end of function order_m ========================= */



/*==============================================================================
|                                order_r                                       |
================================================================================

DESCRIPTION
                                              n
             r                               p - 1
    Compute x  (mod f(x), p) = a, where r = -------
                                             p - 1

    If a is not an integer, return 0, else return a itself.

EXAMPLE
             4    2
     f(x) = x  + x  + 2 x + 3, n = 4 and p = 5.  Then r = 156 and
	 
      r    156
     x  = x    = 3 (mod f(x), 5) = 3, so we return 3.

                     4    
     But for f(x) = x  + x + 3, n = 4, p = 5,

      r    156      3
     x  = x    = 3 x + 2 x + 1 (mod f(x), 5) so we return 0.

METHOD
                          r
    First compute g(x) = x (mod f(x), p).
    Then test if g(x) is a constant polynomial.

BUGS

    None.

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

uint PolyOrder::order_r()
                 throw( PolynomialRangeError, bad_exception )
{
    Polynomial x1( "x", p_ ) ;
    PolyMod x( x1, f_ ) ;
		
    PolyMod x_to_r = power( x, r_ ) ;

    #ifdef DEBUG_PP_POLYNOMIAL
    cout << "r = " << r_ << endl ;
    cout << "x^r = " << x_to_r << endl ;
    cout << "is integer = " << x_to_r.isInteger() << endl ;
    #endif

    if (x_to_r.isInteger())
        //  Return the value a = constant term of g(x).
        return x_to_r[ 0 ] ;
    else
        return 0 ;

} /* ====================== end of function order_r ========================= */



/*==============================================================================
|                                  maximal_order                               |
================================================================================

DESCRIPTION
              k                                  n
    Check if x  = 1 (mod f(x), p) only when k = p  - 1 and not for any smaller
	power of k, i.e. that f(x) is a primitive polynomial.

INPUT

     f (int *)                Monic polynomial f(x).
     n      (int, n >= 1)     Degree of f(x).
     p      (int)             Modulo p coefficient arithmetic.

RETURNS

     true    if f( x ) is primitive.
     false   if it isn't.

EXAMPLE

               4
     f( x ) = x  + x  +  1 is a primitive polynomial modulo p = 2,
                                         4
     because it generates the group GF( 2  ) with the exception of
                            2         15
     zero.  The powers {x, x , ... , x  } modulo f(x), mod 2 are
                                        16       4
     distinct and not equal to 1 until x   (mod x  + x + 1, 2) = 1.

METHOD

    Confirm f(x) is primitive using the definition of primitive
    polynomial as a generator of the Galois group
         n                   n
    GF( p ) by testing that p - 1 is the smallest power for which
     k
    x = 1 (mod f(x), p).

BUGS

    None.

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

bool PolyOrder::maximal_order()
                 throw( PolynomialRangeError, bad_exception )
{
    /*  Highest possible order for x. */
    BigInt maxOrder = power( f_.modulus(), f_.deg()) - 1u ; // TODO allow int type.
    
    BigInt k = 1u ;
	Polynomial x1( "x", f_.modulus() ) ;
    PolyMod x( x1, f_ ) ;    // g(x) = x (mod f(x), p)
            
    while ( k <= maxOrder )
    {
        PolyMod x_to_power = power( x, k ) ; // x^k

        if (x_to_power.isInteger() &&
            x_to_power[0] == 1u &&
            k < maxOrder)
        {
            return false ;
        }

		++k ;

    } /* end for k */

    return true ;

} /* ================= end of function maximal_order ======================== */


/*==============================================================================
|                              has_multi_irred_factors                         |
================================================================================

DESCRIPTION

   Find out if the monic polynomial f( x ) has two or more distinct
   irreducible factors.

INPUT

    power_table (int **)   x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.
    n (int, n >= 1)        Degree of monic polynomial f(x).
    p (int, p >= 2)        Modulo p coefficient arithmetic.

RETURNS

   1 if the monic polynomial f( x ) has two or more distinct rreducible factors
   0 otherwise.

EXAMPLE

   Let n = 4, p = 5

             4    2
   f( x ) = x  + x  + 2 x + 3 is irreducible, so has one distinct factor.

             4    3   2                  4
   f( x ) = x + 4x + x + 4x + 1 = (x + 1)  has one distinct factor.

             3         2
   f( x ) = x  + 3 = (x + 3x + 4)(x + 2) has two distinct irreducible factors.

             4    3    2                          2
   f( x ) = x + 3x + 3x + 3x + 2 = (x + 1) (x + 2) (x + 3) has 3 distinct
   irreducible factors.

           4
   f(x) = x + 4 = (x+1)(x+2)(x+3)(x+4) has 4 distinct irreducible factors.

METHOD

      Berlekamp's method for factoring polynomials over GF( p ), modified to test
      for irreducibility only.

      See my notes;  I skip the polynomial GCD step which ensures polynomials
      are square-free due to time constraints, but this requires a proof that
      the method is still valid.

BUGS

    None.

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

bool PolyOrder::hasMultipleDistinctFactors() 
    throw( PolynomialRangeError, bad_exception )

{  
    // Generate the Q-I matrix.
    generate_Q_matrix() ;


    // Find nullity of Q-I
    findNullity() ;
    
    
    // If nullity_ >= 2, f( x ) is a reducible polynomial modulo p since it has
    // two or more distinct irreducible factors.
    //                                     e
    // Nullity of 1 implies f( x ) = p( x )  for some power e >= 1 so we cannot
    // determine reducibility.
    if (nullity_ >= 2)
        return true ;

    return false ;

} /* =============== end of function has_multi_irred_factors ================ */


// Find a "random" primitive polynomial.
bool findPrimitivePolynomial( uint p, uint n ) throw( PolynomialRangeError, bad_exception )
{
//
//   Generate and test all possible n th degree, monic, modulo p polynomials
//   f(x).  A polynomial is primitive if passes all the tests successfully.
//

/*                       n
     Initialize f(x) to x  + (-1).  Then, when f(x) passes through function 
                                                                          n
     next_trial_poly for the first time, it will have the correct value, x
*/
Polynomial f ;
f.initial_trial_poly( n, p ) ;

bool is_primitive_poly = false ;  // true as soon as a primitive 
                                  // polynomial is found.  
bool tried_all_poly = false ;

BigInt num_poly( 0u ) ;
BigInt max_num_poly = power( f.modulus(), f.deg() ) ;
PolyOrder order( f ) ;

do {
    f.next_trial_poly() ;      // Try next polynomal in sequence.

#ifdef DEBUG_PP_POLYNOMIAL
    cout << "Trial polynomial " << f << endl ;
#endif

    ++num_poly ;
    tried_all_poly = (num_poly > max_num_poly) ;
	
	order.newPolynomial( f ) ;
	is_primitive_poly = order.isPrimitive() ;

} while( !is_primitive_poly && !tried_all_poly ) ;

cout << endl << endl ;




/*
     Report on success or failure.
*/

if (is_primitive_poly)
{
    cout << "\n\nPrimitive polynomial modulo " << f.modulus() << " of degree " << f.deg() << "\n\n" ;
    cout << f ;
    cout << endl << endl ;
}
else 
{
    cerr << "Internal error:\n" 
            "Tested all " << max_num_poly << " possible polynomials, but "
            "failed to find a primitive polynomial.\n"
            "Please write to the author at artifex@seanerikoconnor.freeservers.com\n\n" ;

    return false ;
}


    return true ;
}

// Check if a given polynomial f(x) of degree n modulo p is primitive.
bool PolyOrder::isPrimitive() throw( PolynomialRangeError, bad_exception )
{

uint num_free_of_linear_factors = 0,// The number of polynomials which have  
                                   //   no linear factors.                    

    num_const_coeff_prim_root = 0, // Number of polynomials whose constant  
                                   //   is a primitive root of p.            

    num_passing_const_coeff_test = 0, // Number of polynomials whose constant
                                     //    term passes a consistency check.  
    num_irred_to_power = 0,        // Number of polynomials which are of the
                                   //   form irreducible poly to a power >= 1.
    num_order_m = 0,
    num_order_r = 0 ;              // The number of polynomials which pass
                                   //   the order_m and order_r tests.   

BigInt max_num_poly = power( p_, n_ ) ;
ArithModP modp( p_ ) ; 

bool printStatistics = false ;  // TODO link this.
bool isPrimitive = false ;

    // Constant coefficient of f(x) * (-1)^n must be a primitive root of p. 
    if (modp.const_coeff_is_primitive_root( f_[0], f_.deg() ))
    {
        ++num_const_coeff_prim_root ;

        if (printStatistics) 
            cout << "    (-1)^n const coeff " << f_[ 0 ] << " is primitive root of " << p_ << endl ;
		
        // f(x) can't have any linear factors. 
        if (!f_.hasLinearFactor())
        {
            ++num_free_of_linear_factors ;

            if (printStatistics) 
                cout << "    No linear factors" << endl ;

            // Do more in-depth checking.
            
            // f(x) can't have two or more distinct irreducible factors. 
            if (!hasMultipleDistinctFactors())
            {
                ++num_irred_to_power ;

                if (printStatistics)
                    cout << "    One distinct irreducible factor (possibly repeated)" << endl ;
							
                //  r 
                // x  (mod f(x), p) = a_ must be an integer.
                uint a = order_r() ;
                if (a != 0)
                {
                    ++num_order_r ;

                     if (printStatistics) 
                        cout << "    x^r = a (integer)" << endl ;

                     //              n
                     // Check if (-1)  (constant coefficient of f(x)) = a_ (mod p)
                     //       
                     if (modp.const_coeff_test( f_[ 0 ], a, f_.deg() ))
                     {
                         ++num_passing_const_coeff_test ;

                         if (printStatistics) 
                            cout << "    a (integer) = (-1)^n f[0]" << endl ;

                         //  x^m != integer for all m = r / q, q a prime divisor of r. 
                         if (order_m())
                         {
                             ++num_order_m ;
							 isPrimitive = true ;
                             goto Exit ;

                             if (printStatistics) 
                                cout << "    x^m != integer for m = r / prime divisor of r" << endl ;
                         }
                     } // end const coeff test 
                } // end order r 
            } // end can't determine if reducible 
        } // end no linear factors 
    } // end constant coefficient primitive. 
	
//  Print the statistics of the primitivity tests. 
Exit: 
if (printStatistics)
{
    cout << "+--------- Statistics ----------------------------\n" ;
    cout << "|\n" ;
    cout << "| Total num. degree " << f_.deg() << " poly mod " << f_.modulus() << " :     " << max_num_poly << endl ;
    cout << "|\n" ;
    cout << "| Actually tested :                     " << 1 << endl ; //num_poly << endl ;
	cout << "|\n" ;
    cout << "| a. Const. coeff. was primitive root : " << num_const_coeff_prim_root << endl ;
    cout << "| b. Free of linear factors :           " << num_free_of_linear_factors << endl ;
    cout << "| c. Irreducible or irred. to power :   " << num_irred_to_power << endl ;
    cout << "| d. Had order r (x^r = integer) :      " << num_order_r << endl ;
    cout << "| e. Passed const. coeff. test :        " << num_passing_const_coeff_test << endl ;
    cout << "| f. Had order m (x^m != integer) :     " << num_order_m << endl ;
    cout << "|\n" ;
    cout << "+-------------------------------------------------\n" ;
}
	
return isPrimitive ;
}


/*==============================================================================
|                              generate_Q_matrix                               |
================================================================================

DESCRIPTION


    Generate n x n matrix Q - I, where rows of Q are the powers,

        p                2p                      (n-1) p
    1, x  (mod f(x),p), x  (mod f(x), p), ... , x       (mod f(x), p)

    for monic polynomial f(x).


INPUT

    Q (int **)             Memory is allocated for this matrix already and
                           all entries are 0.

                            k
    power_table (int **)   x (mod f(x), p) for n <= k <= 2n-2, f monic.

    n (int, n >= 1)        Degree of monic polynomial f(x).

    p (int, p >= 2)        Modulo p coefficient arithmetic.

RETURNS


EXAMPLE

            4   2
    f(x) = x + x + 2 x + 3, n = 4, p = 5

        (    1     )    (    1              )       ( 1    0    0    0 )
		 
			  5                   2       3
    Q = (    x     )    (  2x + 3x   + 4 x  )       ( 0    2    3    4 )
	
	                              2     3
        (     10   ) =  (  3 + 4 x   + x    )   =  
        (    x     )                              ( 3    0    4    1 )
		
        (     15   )               2    3
        (    x     )    (   4x + 4x + 3x  )       ( 0    4    4    3 )

                                                              2    3
                                                    1   x    x    x

             ( 0 0 0 0 )
             ( 0 1 3 4 )
    Q - I =  ( 3 0 3 1 ) 
             ( 0 4 4 2 )
			 
			 
	The left nullspace has dimension = 1 with basis { (1 0 0 0) }.

METHOD

BUGS

    None.

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

void PolyOrder::generate_Q_matrix()
                  throw( PolynomialRangeError, bad_exception )
{
   // Check for invalid inputs.
if (n_ < 2 || p_ < 2)
    throw PolynomialRangeError( "generate Q matrix has n < 2 or p < 2" ) ;


// Row 0 of Q = (1 0 ... 0).
Q_[ 0 ][ 0 ] = 1 ;
for (int i = 1 ;  i < n_ ;  ++i)
    Q_[ 0 ][ i ] = 0 ;


/*             p
   Let q(x) = x (mod f(x),p)
   and Q[ 1 ] = coefficients of q(x).
*/
Polynomial x1( "x", p_ ) ;
PolyMod x( x1, f_ ) ;
PolyMod xp = power( x, static_cast< BigInt >( p_ ) ) ;

#ifdef DEBUG_PP_POLYNOMIAL
	cout << "x ^ p (mod f(x),p) = " << xp << endl ;
cout << "initial Q matrix " << printQMatrix() ;
#endif

PolyMod q = xp ;

for (int i = 0 ;  i < n_ ;  ++i)
    Q_[ 1 ][ i ] = xp[ i ] ;


/*               pk
   Row k of Q = x   (mod f(x), p) 2 <= k <= n-1, computed by
                                     p
   multiplying each previous row by x (mod f(x),p).
*/
for (uint k = 2 ;  k <= n_- 1 ;  ++k)
{
    q *= xp ;
	
	#ifdef DEBUG_PP_POLYNOMIAL
	cout << "x ^ pk (mod f(x),p) = " << q << " for k = " << k << endl ;
    #endif
	
    for (int i = 0 ;  i < n_ ;  ++i)
         Q_[ k ][ i ] = q[ i ] ;
}

#ifdef DEBUG_PP_POLYNOMIAL
cout << "computed Q matrix " << printQMatrix() ;
#endif

/*  Subtract Q - I */
for (uint row = 0 ;  row < n_ ;  ++row)
{
    Q_[ row ][ row ] = mod( Q_[ row ][ row ] - 1 ) ;
}

#ifdef DEBUG_PP_POLYNOMIAL
cout << "computed Q-I matrix " << printQMatrix() ;
#endif

return ;

} /* ================== end of function generate_Q_matrix =================== */


string PolyOrder::printQMatrix() const
{
    // Print the matrix as a string.
    ostringstream os ;
    
    os << endl ;
    for (int row = 0 ;  row < n_ ;  ++row)
	{
	    os << "( " ;
        for (int col = 0 ;  col < n_ ;  ++col)
        {
            os << setw( 4 ) << setfill( ' ' ) << Q_[ row ][ col ] ;
        }
		os << " )" << endl ;
	}

    return os.str() ;
}


/*==============================================================================
|                                  findNullity                                |
================================================================================

DESCRIPTION


INPUT

    Q (int **)          n x n matrix of integers mod p.
    n (int, n >= 1)     Degree of monic polynomial f(x).
    p (int, p >= 2)     Modulo p coefficient arithmetic.

RETURNS

    Nullity of Q.  However if the nullity is 2 or more, we just return 2.

EXAMPLE

    Let p = 5 and n = 3.  We use the facts that -1 = 4 (mod 5), 1/2 = 3, -1/2 = 2,
    1/3 = 2, -1/3 = 3, 1/4 = 4, -1/4 = 1.

    Consider the matrix

        ( 2 3 4 )
    Q = ( 0 2 1 )
        ( 3 3 3 )

    Begin with row 1.  No pivotal columns have been selected yet.  Scan row 1 and
    pick column 1 as the pivotal column because it contains a nonzero element.

    Normalizing column 1 by multiplying by -1/pivot = -1/2 = 2 gives

        ( 4 3 4 )
        ( 0 2 1 )
        ( 1 3 3 )

     Now perform column reduction on column 2 by multiplying the pivotal column 1
     by 3 (the column 2 element in the current row) and adding back to row 2.

        ( 4 0 4 )
        ( 0 2 1 )
        ( 1 1 3 )

     Column reduce column 3 by multiplying the pivotal column by 4 and adding back to row 3,

        ( 4 0 0 )
        ( 0 2 1 )
        ( 1 1 2 )

     For row 2, pick column 2 as the pivotal column, normalize it and reduce columns 1, then 3,

        ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )
        ( 0 2 1 ) => ( 0 4 1 ) => ( 0 4 1 ) => ( 0 4 0 )
        ( 1 1 2 )    ( 1 2 2 )    ( 1 2 2 )    ( 1 2 4 )
                 norm.         c.r. 1      c.r. 3

     For row 3, we must pick column 3 as pivotal column because we've used up columns 1 and 2,

        ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )
        ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 )
        ( 1 2 4 )    ( 1 2 4 )    ( 1 2 4 )    ( 0 0 4 )
                 norm.        c.r. 1        c.r. 2

     The nullity is zero, since we were always able to find a pivot in each row.

METHOD

      Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth, Addison-Wesley.

      We combine operations of normalization of columns,

      (       c       )                         (        C       )
      (       c       )                         (        C       )
      (       .       )                         (        C       )
      ( . . . q . . . ) row  ================>  ( . . . -1 . . . ) row
      (       c       )                         (        C       )
      (       c       )      column times       (        C       )
      (       c       )      -1/a modulo p      (        C       )
           pivotCol                                   pivotCol

      and column reduction,

      (        C      b       )                         (       C        B       )
      (        C      b       )                         (       C        B       )
      (        C      b       )                         (       C        B       )
      ( . . . -1 . . .e . . . ) row  ================>  ( . . . -1 . . . 0 . . . )
      (        C      b       )                         (       C        B       )
      (        C      b       )    pivotCol times       (       C        B       )
      (        C      b       )    e added back to col  (       C        B       )
           pivotCol  col                                       col

      to reduce the matrix to a form in which columns having pivots are zero until
      the pivotal row.

      The column operations don't change the left nullspace of the
      matrix.

      The matrix rank is the number of pivotal rows since they are linearly
      independent.  The nullity is then the number of non-pivotal rows.

BUGS


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

void PolyOrder::findNullity()
                  throw( PolynomialRangeError, bad_exception )
{
    InverseModP invmod( p_ ) ;
	
#ifdef DEBUG_PP_POLYNOMIAL
cout << "Q-I matrix " << printQMatrix() ;
#endif

    vector< bool >pivotInCol( n_, false ) ; // Is false if the column has no pivotal element. 

    nullity_ = 0 ;
    int pivotCol = -1 ; // No pivots yet.

    // Sweep through each row. 
    for (int row = 0 ;  row < n_ ;  ++row)
    {
        // Search for a pivot in this row:  a non-zero element
        // in a column which had no previous pivot.
        bool found = false ;
        for (int col = 0 ;  col < n_ ;  ++col)
        {
            if (Q_[ row ][ col ] > 0 && !pivotInCol[ col ])
            {
                found = true ;
                pivotCol = col ;
                break ;
            }
        }

        // No pivot;  increase nullity by 1. 
        if (found == false)
        {
            ++nullity_ ;

             // Early out
//            if (nullity_ >= 2)
//                return nullity ;
        }
        // Found a pivot, q. 
        else
        {
            int q = Q_[ row ][ pivotCol ] ;

            // Compute -1/q (mod p) 
            int t = mod( -invmod( q )) ;

            // Normalize the pivotal column. 
            for (int r = 0 ;  r < n_ ;  ++r)
            {
                Q_[ r ][ pivotCol ] = mod( t * Q_[ r ][ pivotCol ]) ;
            }

            // Do column reduction:  Add C times the pivotal column to the other
            // columns where C = element in the other column at current row.
            for (int col = 0 ;  col < n_ ;  ++col)
            {
                if (col != pivotCol)
                {
                    q = Q_[ row ][ col ] ;

                    for (int r = 0 ;  r < n_ ;  ++r)
                    {
                        t = mod( q * Q_[ r ][ pivotCol ]) ;
                        Q_[ r ][ col ] = mod( t + Q_[ r ][ col ] ) ;
                    }
                }
            }

            // Record the presence of a pivot in this column. 
            pivotInCol[ pivotCol ] = true ;

#ifdef DEBUG_PP_POLYNOMIAL
cout << "row = " << row << " pivot = " << q << " (-1/q) = " << t << " nullity = " << nullity_
     << " Row reduced Q-I matrix " << printQMatrix() ;
#endif

        } // found a pivot 

    } // end for row 


#ifdef DEBUG_PP_POLYNOMIAL
cout << "Row reduced Q-I matrix " << printQMatrix() ;
#endif

// Automagically free pivotInCol and mod objects.

} // ===================== end of function findNullity =====================