/*==============================================================================
|
|  File Name:
|
|     ppUnitTest.cpp
|
|  Description:
|
|     Unit test routines.
|
|  Legal:
|
|     Primpoly Version 8.0 - 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
|
==============================================================================*/

/*------------------------------------------------------------------------------
|                                Includes                                      |
------------------------------------------------------------------------------*/

#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 <limits>       // Numeric limits.
#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 ;   // So we don't need to say std::vector everywhere.

/*------------------------------------------------------------------------------
|                                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.


// Unit test program for checking classes and member functions.
bool unitTest()
{
    // Place results into a log file in the current directory.
    const char * unitTestLogFileName = "unitTest.log" ;
    bool status = true ;

    ofstream fout( unitTestLogFileName ) ;
    if (!fout)
    {
        cerr << "Unit test:  cannot open output log file" << unitTestLogFileName << (string)__FILE__ << endl ;
        cerr << "Skipping the self check." << endl ;
        return false ;
    }

    fout << legalNotice ;

    try
    {
        // Create a parser with parse tables.
        string s ;
        Value v ;

        Parser p ;

        ////////////////////////////////////////////////////////////////////////
        // Test the parser.
        ////////////////////////////////////////////////////////////////////////

		fout << "\n-------------------------> Testing the parser" ;


        // Parse sentences.
        fout << "\nTEST:  parsing constant 0" ;
        s = "0" ;
        v = p.parse( s ) ;
        if (!(v.scalar_ == 2 && (v.f_.size()-1) == 0 && v.f_[0] == 0))
        {
            fout << "\nERROR: Parser failed " << endl ;
            fout << "while parsing input " << s << endl << " result = " ;  v.print() ; fout << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  parsing polynomial, modulus:  2 x ^ 3 + 3 x + 4, 5" ;
        s = "2 x ^ 3 + 3 x + 4, 5" ;
        v = p.parse( s ) ;
        if (!(v.scalar_ == 5 && (v.f_.size()-1) == 3 &&
              v.f_[0] == 4 && v.f_[1] == 3 && v.f_[2] == 0 && v.f_[3] == 2))
        {
            fout << "\nERROR: Parser failed " << endl ;
            fout << "while parsing input " << s << endl << " result = " ;  v.print() ; fout << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  parsing polynomial, modulus:  2x, 7." ;
        s = "2x, 7" ;
        v = p.parse( s ) ;
        if (!(v.scalar_ == 7 && (v.f_.size()-1) == 1 &&
              v.f_[0] == 0 && v.f_[1] == 2 ))
        {
            fout << "\nERROR: Parser failed " << endl ;
            fout << "while parsing input " << s << endl << " result = " ;  v.print() ; fout << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        //fout << "Skip the parser test by throwing a bogus exception." << endl ;
        //throw "Skip the parser test" ;

        fout << "\nTEST:  parsing constant 2" ;
        s = "2" ;
        v = p.parse( s ) ;
        if (!(v.scalar_ == 2 && (v.f_.size()-1) == 0 &&
              v.f_[0] == 2))
        {
            fout << "\nERROR: Parser failed " << endl ;
            fout << "while parsing input " << s << endl << " result = " ;  v.print() ; fout << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  parsing bad syntax:  x 2 ^ + x" ;
        try
        {
            s = "x 2 ^ + x" ;
            v = p.parse( s ) ;
            if (!(v.scalar_ == 0 && (v.f_.size()-1) == 0 &&
                  v.f_[0] == 0 && v.f_[1] == 2 ))
            {
                fout << "\nERROR: Parser failed to throw parsing error " << endl ;
                fout << "while parsing input " << s << endl << " result = " ;  v.print() ; fout << endl ;
                status = false ;
            }
        }
        catch( ParserError & e )
        {
            if ((string)(e.what()) != "Expecting to see x^ or x or x ^ integer")
            {
                fout << "\nERROR: Parser failed " << endl ;
                fout << "while parsing input " << s << endl << " result = " ;  v.print() ; fout << endl ;
                fout << "incorrect parse error message:  " << endl ;
                fout << "Parser error:  |" << e.what() << "|" << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

        ////////////////////////////////////////////////////////////////////////
        // Test BigInt
        ////////////////////////////////////////////////////////////////////////

  		fout << "\n-------------------------> Testing mulitple precision arithmetic." ;

        uint oldBase = 0 ;

        // ------------------- Set to base 10 ----------------------

        // Set all numbers to base 10 for easier testing.
        try
        {
            BigInt w ;
            oldBase = getBase( w ) ;
            setBase( w, 10 ) ;
            BigInt v ;
            if (getBase(v) != 10)
            {
                fout << "\nERROR: Changing to base 10 for all BigInt objects "
                     << "       failed.  base = " << getBase( v ) << endl ;
                status = false ;
            }
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR: Changing to base 10 for all BigInt objects failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

  		fout << "\n-------------------------> Testing BigInt constructors." ;

        fout << "\nTEST:  BigInt() default constructor." ;
        try
        {
            BigInt u ;
            if (getNumDigits(u) != 0)
            {
                fout << "\nERROR: BigInt default constructor and destructor failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: BigInt default constructor and destructor failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt( 0 ) constructor from zero." ;
        try
        {
            BigInt u( 0 ) ;
            if (getNumDigits(u) != 1 || getDigit( u, 0 ) != 0u)
            {
                fout << "\nERROR:  BigInt constructor from special constant 0 failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt constructor from special constant 0 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt( 1234 ) constructor from uint" ;
        try
        {
            uint d = 1234 ;
            BigInt u( d ) ;
            if (getNumDigits( u ) != 4 ||
                getDigit( u, 3 )  != 1 || getDigit( u, 2 ) != 2 ||
                getDigit( u, 1 )  != 3 || getDigit( u, 0 ) != 4
               )
            {
                fout << "\nERROR:  BigInt constructor from uint = 1234 failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt constructor from uint = 1234 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt u( \"1234\" ) constructor from <decimal number string>" ;
        try
        {
            string s = "1234" ;
            BigInt u( s ) ;
            if (getNumDigits( u ) != 4 ||
                getDigit( u, 3 )  != 1 || getDigit( u, 2 ) != 2 ||
                getDigit( u, 1 )  != 3 || getDigit( u, 0 ) != 4
               )
            {
                fout << "\nERROR:  BigInt constructor from number string 1234 failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntRangeError & e )
        {
            fout << "\nERROR:  BigInt constructor from number string 1234 failed." << endl ;
            fout << "BigIntRangeError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt u( \"12x34\" ) constructor from invalid <decimal number string>" ;
        try
        {
            string s = "12x34" ;
            BigInt u( s ) ;
            fout << "\nERROR:  BigInt constructor from invalid number string 12x34 failed." << endl ;
            status = false ;
        }
        catch( BigIntRangeError & e )
        {
            fout << ".........PASS!" ;
        }

  		fout << "\n-------------------------> Testing BigInt copy and assign." ;

        fout << "\nTEST:  BigInt copy constructor: u( 123 ) ; v( u )." ;
        try
        {
            BigInt u( 123 ) ;
            BigInt v( u ) ;

            if (getNumDigits( u ) != getNumDigits( v ) ||
                getDigit( u, 0 ) != getDigit( v, 0 ) ||
                getDigit( u, 1 ) != getDigit( v, 1 ) ||
                getDigit( u, 2 ) != getDigit( v, 2 ))
            {
                fout << "\nERROR:  BigInt copy constructor failed." << endl ;
                printNumber( u ) ;
                printNumber( v ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR:  BigInt copy constructor failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt u = v assignment for BigInt" ;
        try
        {
            BigInt u( 123 ) ;
            BigInt v ; // Don't use BigInt v = u because that would call the copy constructor.
            v = u ;

            if (getNumDigits( u ) != getNumDigits( v ) ||
                getDigit( u, 0 ) != getDigit( v, 0 ) ||
                getDigit( u, 1 ) != getDigit( v, 1 ) ||
                getDigit( u, 2 ) != getDigit( v, 2 ))
            {
                fout << "\nERROR:  BigInt operator= failed." << endl ;
                printNumber( u ) ;
                printNumber( v ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt operator= failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

  		fout << "\n-------------------------> Testing BigInt conversions." ;

        try
        {
            fout << "\nTEST:  cast BigInt to uint implicitly Bigint u ; uint d = u." ;
            BigInt u( "01234" ) ;
            uint d = u ;
            if (d != 1234u)
            {
                fout << "\nERROR:  Implicit casting of BigInt to uint failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;

            fout << "\nTEST:  cast BigInt to uint explicitly: BigInt u ; static_cast<uint>(u)." ;
            if (static_cast<uint>(u) != 1234u)
            {
                fout << "\nERROR:  static casting BigInt to uint failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  Casting BigInt to uint failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  overflow during cast from BigInt to uint?" ;
        try
        {
            BigInt u( "3141592653589793238462643383279" ) ;

            uint ui = static_cast<uint>(u) ;

            fout << "\nERROR:  Casting BigInt didn't overflow." << endl ;
            printNumber( u ) ;
            fout << "casted number = " << ui << endl ;
            status = false ;
        }
        catch( BigIntOverflow & e )
        {
            // Should overflow!
            fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  Casting BigInt didn't throw overflow exception." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

  		fout << "\n-------------------------> Testing BigInt I/O." ;

        fout << "\nTEST:  stream output << BigInt" ;
        try
        {
            BigInt u( "1234567890" ) ;
            ostringstream os ;
            os << u ;
            if (os.str() != "1234567890")
            {
                fout << "\nERROR:  BigInt stream output failed." << endl ;
                fout << "BigInt = |" << u << "|" << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt stream output failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  stream input >> BigInt" ;
        try
        {
            BigInt u ;
            istringstream is ;
            is.clear() ;
            is.str( "314159265358979323846264" ) ;
            is >> u ;

            ostringstream os ;
            os << u ;
            if (os.str() != "314159265358979323846264" )
            {
                fout << "\nERROR:  BigInt stream input failed." << endl ;
                fout << "BigInt = |" << u << "|" << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt stream input failed." << endl ;
            status = false ;
        }

  		fout << "\n-------------------------> Testing BigInt comparisions." ;

        fout << "\nTEST:  BigInt == uint" ;
        {
            BigInt u( "9" ) ;
            uint d = 9 ;
            if (u != d)
            {
                fout << "\nERROR:  BigInt to uint equality test failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

        fout << "\nTEST:  BigInt == BigInt" ;
        {
            BigInt u( "1234" ) ;
            BigInt v( "1234" ) ;
            if (u != v)
            {
                fout << "\nERROR:  BigInt to BigInt equality test failed." << endl ;
                printNumber( u ) ;
                printNumber( v ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

        fout << "\nTEST:  BigInt > uint" ;
        {
            BigInt u( "933134" ) ;
            BigInt v( "3844035" ) ;
	
            if (u > v)
            {
                fout << "\nERROR:  BigInt to BigInt > test failed.  u = " << u  << "v = " << v << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

  		fout << "\n-------------------------> Testing BigInt operator=(uint)." ;

        fout << "\nTEST:  BigInt -= uint" ;
        try
        {
            BigInt u( "1234" ) ;
            u -= 5u ;
            if (getNumDigits( u ) != 4 ||
                getDigit( u, 3 )  != 1 || getDigit( u, 2 ) != 2 ||
                getDigit( u, 1 )  != 2 || getDigit( u, 0 ) != 9
               )
            {
                fout << "\nERROR:  BigInt - uint failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt - uint failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt -= uint underflow" ;
        try
        {
            try
            {
                BigInt u( "4" ) ;
                u -= 5u ;

                fout << "\nERROR:  BigInt - uint failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            catch( BigIntUnderflow & e )
            {
               // Caught underflow;  works correctly.
               // fout << "Underflow: " << e.what() << endl ;
                fout << ".........PASS!" ;
            }
            catch( ... )
            {
                fout << "\nERROR:  BigInt - uint failed." << endl ;
                status = false ;
            }
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt - uint failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt += uint" ;
        try {
            BigInt u( "9994" ) ;
            u += 6 ;
            if (getNumDigits( u ) != 5 ||
                getDigit( u, 4 )  != 1 ||
                getDigit( u, 3 )  != 0 || getDigit( u, 2 ) != 0 ||
                getDigit( u, 1 )  != 0 || getDigit( u, 0 ) != 0
               )
            {
                fout << "\nERROR:  BigInt + uint failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt + uint failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt *= uint" ;
        try
        {
            BigInt u( "123" ) ;
            uint d = 4 ;
            BigInt v = u * d ;

            if (getNumDigits( v ) != 3 ||
                getDigit( v, 2 ) != 4 ||
                getDigit( v, 1 ) != 9 || getDigit( v, 0 ) != 2
               )
            {
                fout << "\nERROR:  BigInt * uint failed." << endl ;
                printNumber( u ) ;
                printNumber( v ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt * uint failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt /= uint" ;
        try
        {
            BigInt u( "12" ) ;
            uint d = 4 ;
            u /= d ;

            if (getNumDigits( u ) != 1 || getDigit( u, 0 ) != 3
               )
            {
                fout << "\nERROR:  BigInt / uint failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt / uint failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt /= uint underflow to zero." ;
        try
        {
            BigInt u( "3" ) ;
            uint d = 4 ;
            u /= d ;

            if (getNumDigits( u ) != 1 || getDigit( u, 0 ) != 0)
            {
                fout << "\nERROR:  BigInt / uint integer truncate to zero failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt / uint integer truncate to zero failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

  		fout << "\n-------------------------> Testing BigInt ++ and --." ;

        fout << "\nTEST:  ++BigInt" ;
        try
        {
            BigInt u( "123" ) ;
            BigInt v = ++u ;

            if (getNumDigits( u ) != 3 ||
                getDigit( u, 2 ) != 1 || getDigit( u, 1 )  != 2 || getDigit( u, 0 ) != 4)
            {
                fout << "\nERROR:  ++BigInt failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            // u better have changed!
            else if (u != v)
            {
                fout << "\nERROR:  v = ++BigInt failed." << endl ;
                printNumber( v ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  ++BigInt failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        try
        {
            fout << "\nTEST:  --BigInt" ;
            BigInt u( "123" ) ;
            BigInt v = --u ;

            if (getNumDigits( u ) != 3 ||
                getDigit( u, 2 ) != 1 || getDigit( u, 1 )  != 2 || getDigit( u, 0 ) != 2)
            {
                fout << "\nERROR:  --BigInt failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            // u better have changed!
            else if (u != v)
            {
                fout << "\nERROR:  v = BigInt-- failed." << endl ;
                printNumber( v ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  v = --BigInt failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt++" ;
        try
        {
            BigInt u( "123" ) ;
            BigInt v = u++ ;

            if (getNumDigits( u ) != 3 ||
                getDigit( u, 2 ) != 1 || getDigit( u, 1 )  != 2 || getDigit( u, 0 ) != 4
               )
            {
                fout << "\nERROR:  ++BigInt failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else if (getNumDigits( v ) != 3 ||
                getDigit( v, 2 ) != 1 || getDigit( v, 1 )  != 2 || getDigit( v, 0 ) != 3
               )
            {
                fout << "\nERROR:  BigInt++ failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt++ failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt--" ;
        try
        {
            BigInt u( "123" ) ;
            BigInt v = u-- ;

            if (getNumDigits( u ) != 3 ||
                getDigit( u, 2 ) != 1 || getDigit( u, 1 )  != 2 || getDigit( u, 0 ) != 2
               )
            {
                fout << "\nERROR:  BigInt-- failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else if (getNumDigits( v ) != 3 ||
                getDigit( v, 2 ) != 1 || getDigit( v, 1 )  != 2 || getDigit( v, 0 ) != 3
               )
            {
                fout << "\nERROR:  BigInt-- failed." << endl ;
                printNumber( u ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  v = BigInt-- failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

  		fout << "\n-------------------------> Testing BigInt + and -." ;

        fout << "\nTEST:  one digit BigInt + uint" ;
        try
        {
            BigInt u( "3" ) ;
            uint d = 4 ;
            BigInt w = u + d ;

            if (getNumDigits( w ) != 1 || getDigit( w, 0 ) != 7)
            {
                fout << "\nERROR:  BigInt + BigInt 3 + 4 = 7 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt + BigInt 3 + 4 = 7 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  two digit BigInt + uint" ;
        try
        {
            BigInt u( "3" ) ;
            uint d = 9 ;
            BigInt w = u + d ;

            if (getNumDigits( w ) != 2 || getDigit( w, 1 ) != 1
                                       || getDigit( w, 0 ) != 2)
            {
                fout << "\nERROR:  BigInt + BigInt 3 + 9 = 12 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt + BigInt 3 + 9 = 12 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt + BigInt" ;
        try
        {
            BigInt u( "9999" ) ;
            BigInt v(  "999" ) ;
            BigInt w = u + v ;

            if (getNumDigits( w ) != 5 || getDigit( w, 4 ) != 1
                                       || getDigit( w, 3 ) != 0
                                       || getDigit( w, 2 ) != 9
                                       || getDigit( w, 1 ) != 9
                                       || getDigit( w, 0 ) != 8)
            {
                fout << "\nERROR:  BigInt + BigInt 9999 + 999 = 10998 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt + BigInt 9999 + 999 = 10998 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt + BigInt" ;
        try
        {
            BigInt u(  "999" ) ;
            BigInt v( "9999" ) ;
            BigInt w = u + v ;

            if (getNumDigits( w ) != 5 || getDigit( w, 4 ) != 1
                                       || getDigit( w, 3 ) != 0
                                       || getDigit( w, 2 ) != 9
                                       || getDigit( w, 1 ) != 9
                                       || getDigit( w, 0 ) != 8)
            {
                fout << "\nERROR:  BigInt + BigInt 999 + 9999 = 10998 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt + BigInt 999 + 9999 = 10998 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt - BigInt" ;
        try
        {
            BigInt u( "103" ) ;
            BigInt v(   "9" ) ;
            BigInt w = u - v ;

            if (getNumDigits( w ) != 2 || getDigit( w, 1 ) != 9
                                       || getDigit( w, 0 ) != 4)
            {
                fout << "\nERROR:  BigInt - BigInt 103 - 9 = 94 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt - BigInt 103 - 9 = 94 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt - BigInt < 0" ;
        try
        {
            BigInt u( "9" ) ;
            BigInt v( "103" ) ;
            BigInt w = u - v ;

            fout << "\nERROR:  BigInt - BigInt 9 - 103 failed didn't catch range exception." << endl ;
            printNumber( w ) ;
            status = false ;
        }
        catch( BigIntUnderflow & e )
        {
            // Caught underflow;  works correctly.
            // fout << "Underflow: " << e.what() << endl ;
            fout << ".........PASS!" ;
        }

        fout << "\nTEST:  BigInt - uint" ;
        try
        {
            BigInt u( "103" ) ;
            uint d = 9 ;
            BigInt w = u - d ;

            if (getNumDigits( w ) != 2 || getDigit( w, 1 ) != 9
                                       || getDigit( w, 0 ) != 4)
            {
                fout << "\nERROR:  BigInt - uint 103 - 9 = 94 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt - uint 103 - 9 = 94 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

  		fout << "\n-------------------------> Testing BigInt * and /." ;

        fout << "\nTEST:  one digit BigInt * BigInt" ;
        try
        {
            BigInt u( "3" ) ;
            BigInt v( "3" ) ;
            BigInt w = u * v ;

            if (getNumDigits( w ) != 1 || getDigit( w, 0 ) != 9)
            {
                fout << "\nERROR:  BigInt * BigInt 3 * 3 = 9 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt * BigInt failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  two digit BigInt * BigInt" ;
        try
        {
            BigInt u( "3" ) ;
            BigInt v( "4" ) ;
            BigInt w = u * v ;

            if (getNumDigits( w ) != 2 || getDigit( w, 1 ) != 1 || getDigit( w, 0 ) != 2)
            {
                fout << "\nERROR:  BigInt * BigInt = 3 * 4 = 12 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt * BigInt = 3 * 4 = 12 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt / BigInt one digit divisor." ;
        try
        {
            BigInt u( "12" ) ;
            BigInt v( "4" ) ;
            BigInt w = u / v ;

            if (getNumDigits( w ) != 1 || getDigit( w, 0 ) != 3)
            {
                fout << "\nERROR:  BigInt / BigInt = 12/4 = 3 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
               fout << "\nERROR:  BigInt / BigInt = 12/4 failed." << endl ;
               fout << "BigIntMathError:  " << e.what() << endl ;
                status = false ;
        }

        fout << "\nTEST:  BigInt / BigInt multidigit" ;
        try
        {
            BigInt u( "398765" ) ;
            BigInt v( "3457" ) ;
            BigInt w = u / v ;

            if (getNumDigits( w ) != 3 ||
                    getDigit( w, 2 ) != 1 || getDigit( w, 1 ) != 1 || getDigit( w, 0 ) != 5)
            {
                fout << "\nERROR:  BigInt / BigInt = 398765/3457 = 215 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
               fout << "\nERROR:  BigInt / BigInt = 398765/3457 = 215 failed!" << endl ;
               fout << "BigIntMathError:  " << e.what() << endl ;
                status = false ;
        }

        fout << "\nTEST:  BigInt / BigInt leading zero digit." ;
        try
        {
            BigInt u( "120" ) ;
            BigInt v( "40" ) ;
            BigInt w = u / v ;

            if (getNumDigits( w ) != 1 || getDigit( w, 0 ) != 3)
            {
                fout << "\nERROR:  BigInt / BigInt = 120/40 = 3 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt / BigInt = 120/40 = 3 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt / 0 " ;
        try
        {
            BigInt u( "120" ) ;
            BigInt v( "0" ) ;
            BigInt w = u / v ;

            fout << "\nERROR:  BigInt / 0 = 120/0 failed." << endl ;
            status = false ;
        }
        catch( BigIntZeroDivide & e )
        {
            // Should catch zero divide here.
            fout << ".........PASS!" ;
        }

  		fout << "\n-------------------------> Testing BigInt % (mod)." ;

        fout << "\nTEST:  BigInt % BigInt with u > v" ;
        {
            BigInt u( "398765" ) ;
            BigInt v( "3457" ) ;
            BigInt r = u % v ;

            if (getNumDigits( r ) != 4 ||
                    getDigit( r, 3 ) != 1 || getDigit( r, 2 ) != 2 ||
                    getDigit( r, 1 ) != 1 || getDigit( r, 0 ) != 0)
            {
                fout << "\nERROR:  BigInt % BigInt = 398765 / 3457 = 1210 failed." << endl ;
                printNumber( r ) ;
                status = false ;
            }
            else
			{
                fout << ".........PASS!" ;
			}
        }

        fout << "\nTEST:  BigInt multidigit mod with normalizing constant d = 1" ;
        {
            BigInt w ;
            BigInt u( "1369244731822264511994463394" ) ;
            BigInt v(  "954901783703457032047844259" ) ;

            w = u % v ;

            string s = w.to_string() ;

            if (s != "414342948118807479946619135" )
            {
                fout << "\nERROR:  BigInt multidigit mod failed." << endl ;
                fout << "u = " ; printNumber( u ) ;
                fout << "v = " ; printNumber( v ) ;
                fout << "w = " ; printNumber( w ) ;
                status = false ;
            }
            else
			{
                fout << ".........PASS!" ;
			}
        }

        fout << "\nTEST:  BigInt % BigInt with u < v" ;
        try
        {
            BigInt u( "12" ) ;
            BigInt v( "34567" ) ;
            BigInt r = u % v ;

            if (getNumDigits( r ) != 2 ||
                getDigit( r, 1 ) != 1 || getDigit( r, 0 ) != 2)
            {
                fout << "\nERROR:  BigInt % BigInt = 12 mod 345 = 12 failed." << endl ;
                printNumber( r ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt % BigInt = 12 mod 345 = 12 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  BigInt % uint = 314159 / 9 = 5 with uint < base " ;
        try
        {
            BigInt u( "314159" ) ;
            uint v = 9u  ;
            BigInt r = u % v ;

            if (getNumDigits( r ) != 1 || getDigit( r, 0 ) != 5)
            {
                fout << "\nERROR:  BigInt % uint = 314159 / 9 = 5 failed." << endl ;
                printNumber( r ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
               fout << "\nERROR:  BigInt % uint failed 314159 / 9 = 5 ." << endl ;
               fout << "BigIntMathError:  " << e.what() << endl ;
                status = false ;
        }

        fout << "\nTEST:  BigInt % uint = 398765 % 3457u with uint > base overflow?" ;
        try
        {
            BigInt u( "398765" ) ;
            uint v = 3457u  ;
            BigInt r = u % v ;

            fout << "\nERROR:  BigInt % uint = 398765 % 3457u with uint > base didn't overflow." << endl ;
        }
        catch( BigIntOverflow & e )
        {
	        fout << ".........PASS!" ;
        }

        fout << "\nTEST:  BigInt / BigInt low probability if branch." ;
        try
        {
            BigInt u( "4100" ) ;
            BigInt v( "588" ) ;
            BigInt w = u / v ;

            if (w != BigInt("6"))
                fout << "error" << endl;

            if (getNumDigits( w ) != 1 || getDigit( w, 0 ) != 6)
            {
                fout << "\nERROR:  BigInt / BigInt = 4100/588 = 6 failed." << endl ;
                printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt / BigInt = 4100/588 = 6 failed." << endl ;
            fout << "BigIntMathError:  " << e.what() << endl ;
            status = false ;
        }


  		fout << "\n-------------------------> Testing BigInt with default large base." ;

        BigInt dummy ;
        setBase( dummy, oldBase ) ;

        fout << "\nTEST:  Decimal string to BigInt to string with large base." ;
        {
            BigInt x( "3141592653589793238462643383279" ) ;
            string s = x.to_string() ;

            if (s != "3141592653589793238462643383279" )
            {
                fout << "\nERROR:  BigInt large base conversion failed." << endl ;
                fout << "x = " << x << " " ; printNumber( x ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

        fout << "\nTEST:  BigInt */ BigInt multidigit with large base." ;
        {
            BigInt w ;
            BigInt x( "3141592653589793238462643383279" ) ;
            BigInt y( "2718281828459045" ) ;
            BigInt z ;

            z = x * y ;
            w = z / y ;

            if (w != x)
            {
                fout << "\nERROR:  BigInt multiply then inverse divide = failed." << endl ;
                fout << "x = " ; printNumber( x ) ;
                fout << "y = " ; printNumber( y ) ;
                fout << "z = " ; printNumber( z ) ;
                fout << "w = " ; printNumber( w ) ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

  		fout << "\n-------------------------> Testing BigInt bit operations." ;

        fout << "\nTEST:  BigInt testBit" ;
        try
        {
            BigInt u( "31415926535897932" ) ;

            if (u.testBit( 0u ) == false  &&
                u.testBit( 1u ) == false  &&
                u.testBit( 2u ) == true   &&
                u.testBit( 3u ) == true   &&
                u.testBit( 4u ) == false  &&
                u.testBit( 5u ) == false  &&
                u.testBit( 6u ) == true   &&
                u.testBit( 7u ) == false)
            {
                fout << ".........PASS!" ;
            }
            else
            {
                fout << "\nERROR:  BigInt testBit failed." << endl ;
                printNumber( u ) ;
                fout << "testBit 0 = " << (u.testBit( 0u ) == true) << endl ;
                fout << "testBit 1 = " << (u.testBit( 1u ) == true) << endl ;
                fout << "testBit 2 = " << (u.testBit( 2u ) == true) << endl ;
                fout << "testBit 3 = " << (u.testBit( 3u ) == true) << endl ;
                fout << "testBit 4 = " << (u.testBit( 4u ) == true) << endl ;
                fout << "testBit 5 = " << (u.testBit( 5u ) == true) << endl ;
                fout << "testBit 6 = " << (u.testBit( 6u ) == true) << endl ;
                fout << "testBit 7 = " << (u.testBit( 7u ) == true) << endl ;
                status = false ;
            }
        }
        catch( BigIntMathError & e )
        {
               fout << "\nERROR:  BigInt testBit Failed." << endl ;
               fout << "BigIntMathError:  " << e.what() << endl ;
                status = false ;
        }

        fout << "\nTEST:  uint testBit" ;
        try
        {
            uint u = 0x1 | 0x4 | 0x20 ; // Set bits 0, 2, 5

            if (testBit( u, 0u ) == true   &&
                testBit( u, 1u ) == false  &&
                testBit( u, 2u ) == true   &&
                testBit( u, 3u ) == false  &&
                testBit( u, 4u ) == false  &&
                testBit( u, 5u ) == true   &&
                testBit( u, 6u ) == false  &&
                testBit( u, 7u ) == false)
            {
                fout << ".........PASS!" ;
            }
            else
            {
                fout << "\nERROR:  uint testBit failed for u = " << u << endl ;
                fout << "testBit 0 = " << (testBit( u, 0u ) == true) << endl ;
                fout << "testBit 1 = " << (testBit( u, 1u ) == true) << endl ;
                fout << "testBit 2 = " << (testBit( u, 2u ) == true) << endl ;
                fout << "testBit 3 = " << (testBit( u, 3u ) == true) << endl ;
                fout << "testBit 4 = " << (testBit( u, 4u ) == true) << endl ;
                fout << "testBit 5 = " << (testBit( u, 5u ) == true) << endl ;
                fout << "testBit 6 = " << (testBit( u, 6u ) == true) << endl ;
                fout << "testBit 7 = " << (testBit( u, 7u ) == true) << endl ;
                status = false ;
            }
        }
        catch( BigIntMathError & e )
        {
               fout << "\nERROR:  BigInt testBit Failed." << endl ;
               fout << "BigIntMathError:  " << e.what() << endl ;
                status = false ;
        }


  		fout << "\n-------------------------> Testing BigInt power." ;

        fout << "\nTEST:  BigInt power( uint 2, uint 100 )" ;
        try
        {
            uint p = 2u  ;
            uint n = 100u  ;
            BigInt u = power( p, n ) ;
            string s = u.to_string() ;

            BigInt v = 1u ;
            for (uint i = 1 ; i <= n ; ++i)
                v *= p ;
            string sv = v.to_string() ;

            if (s != sv)
            {
                fout << "\nERROR:  BigInt power( 2, 100 ) = " << u << endl ;
                fout << "correct answer = " << v << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt power( 2, 100 ) internal error." << endl ;
                status = false ;
        }

      fout << "\nTEST:  BigInt ceilLg( 6 )" ;
        try
        {
            BigInt u = 6 ;

			int ceilingOfLog2 = u.ceilLg() ;
			if (ceilingOfLog2 != 3)
            {
                fout << "\nERROR:  BigInt ceilingOfLog2( 6 ) = " << ceilingOfLog2 << endl ;
                fout << "correct answer = 3" << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( BigIntMathError & e )
        {
            fout << "\nERROR:  BigInt ceilingOfLog2 BigIntMathError." << endl ;
                status = false ;
        }


        ////////////////////////////////////////////////////////////////////////
        // Test mod p arithmetic
        ////////////////////////////////////////////////////////////////////////

		fout << "\n-------------------------> Testing mod p arithmetic." ;

        fout << "\nTEST:  ModP 10 = 3 (mod 7)" ;
        ModP modp( 7 ) ;
        if (modp( 10 ) != 3)
        {
            fout << "\nERROR:  ModP modp( 7 );  modp( 10 ) = " << modp( 10 ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  ModP -10 = 4 (mod 7)" ;
        if (modp( -10 ) != 4)
        {
            fout << "\nERROR:  ModP modp( 7 );  modp( -10 ) = " << modp( -10 ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  uint gcd( 85, 25 ) = 5" ;
        {
            uint u = 85u ;
            uint v = 25u ;
            uint g = gcd( u, v ) ;
            if (g != 5u)
            {
                fout << "\nERROR:  uint gcd( 85, 25 ) = " << g << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

        fout << "\nTEST:  BigInt gcd( 779953197883173551166308319545, 1282866356929526866866376009397 ) = 1" ;
        {
            BigInt u( "779953197883173551166308319545" ) ;
            BigInt v( "1282866356929526866866376009397" ) ;
            BigInt g = gcd( u, v ) ;
            if (g != 1u)
            {
                fout << "\nERROR:  BigInt gcd = " << g << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }


        fout << "\nTEST:  PowerMod uint 3^10 = 4 (mod 7)" ;
        PowerMod<uint> powermodu( 7u ) ;
        if (powermodu( 3u, 10u ) != 4u)
        {
            fout << "\nERROR:  PowerMod powermod( 7 );  powermod( 3, 10 ) = " << powermodu( 3, 10 ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  PowerMod BigInt 3^10 = 4 (mod 7)" ;
        PowerMod<BigInt> powermod( static_cast<BigInt>(7u) ) ;
        if (powermod( static_cast<BigInt>(3u), static_cast<BigInt>(10u) ) != static_cast<BigInt>(4u))
        {
            fout << "\nERROR:  PowerMod powermod( 7 );  powermod( 3, 10 ) = " << powermod( 3, 10 ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;


        fout << "\nTEST:  InverseModP 3 * 5 = 1 (mod 7)" ;
        InverseModP imodp( 7 ) ;
        if (imodp( 3 ) != 5)
        {
            fout << "\nERROR:  InverseModP imodp( 7 );  imodp( 3 ) = " << imodp( 3 ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  IsPrimitiveRoot.   3 is a primitive root of 7." ;
        IsPrimitiveRoot isroot( 7 ) ;
        if (isroot( 3 ))
            fout << ".........PASS!" ;
        else
        {
            fout << "\nERROR:  IsPrimitiveRoot( 7 ) isroot ;  isroot( 3 ) = " << isroot( 3 ) << " failed." << endl ;
            status = false ;
        }

        fout << "\nTEST:  IsPrimitiveRoot.   2 is a primitive root of 11." ;
        IsPrimitiveRoot isroot11( 11 ) ;
        if (isroot11( 2 ))
            fout << ".........PASS!" ;
        else
        {
            fout << "\nERROR:  IsPrimitiveRoot( 11 ) isroot11 ;  isroot11( 2 ) = " << isroot11( 2 ) << " failed." << endl ;
            status = false ;
        }

        fout << "\nTEST:  IsPrimitiveRoot.   3 is NOT a primitive root of 11." ;
        if (isroot11( 3 ))
        {
            fout << "\nERROR:  IsPrimitiveRoot( 11 ) isroot11 ;  isroot11( 3 ) = " << isroot11( 3 ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  constant coefficient test." ;
        ArithModP arith1( 5 ) ;
        if (arith1.const_coeff_test( 4, 1, 11 ))
            fout << ".........PASS!" ;
        else
        {
            fout << "\nERROR:  constant coefficient test failed = " << arith1.const_coeff_test( 4, 1, 11 ) << " failed." << endl ;
            status = false ;
        }

        fout << "\nTEST:  constant coefficient is primitive root." ;
        ArithModP arith2( 7 ) ;
        if (arith2.const_coeff_is_primitive_root( 4, 11 ))
            fout << ".........PASS!" ;
        else
        {
            fout << "\nERROR:  constant coefficient test failed = " << arith1.const_coeff_is_primitive_root( 4, 11 ) << " failed." << endl ;
            status = false ;
        }

        fout << "\nTEST:  is_probably_prime 97 with random x = 10 says yes." ;
        if ( is_probably_prime( 97u, 10u ))
            fout << ".........PASS!" ;
        else
        {
            fout << "\nERROR:  is_probably_prime( 97u, 10u ) = " << is_probably_prime( 97u, 10u ) << " failed." << endl ;
            status = false ;
        }

        fout << "\nTEST:  is_probably_prime 97 with random x = 9 says yes." ;
        if ( is_probably_prime( 97u, 9u ))
            fout << ".........PASS!" ;
        else
        {
            fout << "\nERROR:  is_probably_prime( 97u, 9u ) = " << is_probably_prime( 97u, 9u ) << " failed." << endl ;
            status = false ;
        }

        fout << "\nTEST:  is_almost_surely_prime 97 says yes." ;
        if ( is_almost_surely_prime( 97u ))
            fout << ".........PASS!" ;
        else
        {
            fout << "\nERROR:  is_almost_surely_prime( 97u ) = " << is_almost_surely_prime( 97u ) << " failed." << endl ;
            status = false ;
        }

        fout << "\nTEST:  is_probably_prime 49 with random x = 10 says no." ;
        if (is_probably_prime( 49u, 10u ))
        {
            fout << "\nERROR:  is_probably_prime( 49u, 10u ) = " << is_probably_prime( 49u, 10u ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  is_almost_surely_prime 49 says no." ;
        if (is_almost_surely_prime( 49u ))
        {
            fout << "\nERROR:  is_almost_surely_prime( 49u ) = " << is_almost_surely_prime( 49u ) << " failed." << endl ;
            status = false ;
        }
        else
            fout << ".........PASS!" ;



        ////////////////////////////////////////////////////////////////////////
        // Test Factoring
        ////////////////////////////////////////////////////////////////////////

		fout << "\n-------------------------> Testing factoring." ;

        fout << "\nTEST:  Factor unsigned int 337500 = 2^2 3^3 5^5" ;
        Factor<uint> f( 337500u ) ;
        if (!(f.count(0) == 2u  && f[0] == 2u &&
              f.count(1) == 3u  && f[1] == 3u &&
              f.count(2) == 5u  && f[2] == 5u ))
        {
            fout << "\nERROR:  Factor<uint> failed on 337500 = 2^2 3^3 5^5." << endl ;
            status = false ;

            fout << "Number of factors = " << f.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < f.num() ;  ++i)
                fout << f[ i ] << " ^ " << f.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  Factor BigInt 337500 = 2^2 3^3 5^5." ;
        Factor<BigInt> f1( static_cast<BigInt>("337500") ) ;
        if (!(f1.count(0) == 2u  && f1[0] == 2u &&
              f1.count(1) == 3u  && f1[1] == 3u &&
              f1.count(2) == 5u  && f1[2] == 5u ))
        {
            fout << "\nERROR:  Factor<BigInt> failed on 337500 = 2^2 3^3 5^5." << endl ;
            status = false ;

            fout << "Number of factors = " << f1.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < f1.num() ;  ++i)
                fout << f1[ i ] << " ^ " << f1.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  Pollard Rho factor unsigned int 25852 = 2^2 23 281" ;
        Factor<uint> fr( 1u ) ;
        bool pollardStatus = fr.PollardRho( 25852 ) ;
        if (pollardStatus == false ||
            !(fr.count(0) == 2u  && fr[0] ==  2u &&
              fr.count(1) == 1u  && fr[1] == 23u &&
              fr.count(2) == 1u  && fr[2] == 281u))
        {
            fout << "\nERROR:  Factor<uint> failed on 25852 = 2^2 23 281." << endl ;

            status = false ;
            fout << "Pollard rho status " << (pollardStatus ? "Success" : "Failure") << endl ;
            fout << "Number of factors = " << fr.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < fr.num() ;  ++i)
                fout << fr[ i ] << " ^ " << fr.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  Pollard Rho factor BigInt 25852 = 2^2 23 281" ;
        Factor<BigInt> frb( 1u ) ;
        pollardStatus = frb.PollardRho( 25852 ) ;
        if (pollardStatus == false ||
            !(frb.count(0) == 2u  && frb[0] ==  2u &&
              frb.count(1) == 1u  && frb[1] == 23u &&
              frb.count(2) == 1u  && frb[2] == 281u))
        {
            fout << "\nERROR:  Factor<BigInt> failed on 25852 = 2^2 23 281." << endl ;

            status = false ;
            fout << "Pollard rho status " << (pollardStatus ? "Success" : "Failure") << endl ;
            fout << "Number of factors = " << frb.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < frb.num() ;  ++i)
                fout << frb[ i ] << " ^ " << frb.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;

#if 0 // Takes too long for a self-check
        fout << "\nTEST:  Pollard Rho Factor BigInt on 10865143687 = 104233 104239" ;
        Factor<BigInt> f10( 1u ) ;
        BigInt x( "10865143687" ) ;
        pollardStatus = f10.PollardRho( static_cast<BigInt>( x ) ) ;
        if (pollardStatus == false ||
            !(f10.count(0) == 1u  && f10[0] == static_cast<BigInt>(104233u) &&
              f10.count(1) == 1u  && f10[1] == static_cast<BigInt>(104239u) ))
        {
            fout << "\nERROR:  Factor<BigInt> failed on x = 104233 104239" << endl ;
            status = false ;

            fout << "Pollard rho status " << (pollardStatus ? "Success" : "Failure") << endl ;
            fout << "Number of factors = " << f10.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < f10.num() ;  ++i)
                fout << f10[ i ] << " ^ " << f10.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;
#endif	

#if 0 // Compare trial division vs Pollard for really big n.
        {
            cout << "\nTEST:  Pollard Rho Factor BigInt on 2^62 - 1" ;
            BigInt n = power( 2u, 62u ) - 1u ;
            cout << " = " << n << endl ;
            Factor<BigInt> f( 1u ) ;
            f.PollardRho( n ) ;
            cout << "Pollard Rho:  Number of factors = " << f.num() << endl ;
            cout << "Factors = " << endl ;
            for (int i = 0 ;  i < f.num() ;  ++i)
                cout << f[ i ] << " ^ " << f.count( i ) << " " ;
            cout << endl ;

            f.trialDivision( n ) ;
            cout << "Trial division factoring:  Number of factors = " << f.num() << endl ;
            cout << "Factors = " << endl ;
            for (int i = 0 ;  i < f.num() ;  ++i)
                cout << f[ i ] << " ^ " << f.count( i ) << " " ;
            cout << endl ;
        }
#endif

        fout << "\nTEST:  Factor Copy constructor" ;
        Factor<BigInt> f2( f1 ) ;
        if (!(f2.count(0) == 2u  && f2[0] == 2u &&
              f2.count(1) == 3u  && f2[1] == 3u &&
              f2.count(2) == 5u  && f2[2] == 5u ))
        {
            fout << "\nERROR:  Factor copy constructor failed on 337500 = 2^2 3^3 5^5." << endl ;
            status = false ;

            fout << "Number of factors = " << f2.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < f2.num() ;  ++i)
                fout << f2[ i ] << " ^ " << f2.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;

		fout << "\nTEST:  Factor assignment operator" ;
        Factor<BigInt> f3 ;
		f3 = f2 ;
        if (!(f3.count(0) == 2u  && f3[0] == 2u &&
              f3.count(1) == 3u  && f3[1] == 3u &&
              f3.count(2) == 5u  && f3[2] == 5u ))
        {
            fout << "\nERROR:  Factor assignment operator failed on 337500 = 2^2 3^3 5^5." << endl ;
            status = false ;

            fout << "Number of factors = " << f3.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < f3.num() ;  ++i)
                fout << f3[ i ] << " ^ " << f3.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  Factor BigInt 24018350267611933650627567399079537500 = 2^2 3^3 5^5 7^7 11^11 13^13." ;
        Factor<BigInt> f6( static_cast<BigInt>("24018350267611933650627567399079537500") ) ;
        if (!(f6.count(0) == 2u  && f6[0] == 2u &&
              f6.count(1) == 3u  && f6[1] == 3u &&
              f6.count(2) == 5u  && f6[2] == 5u &&
              f6.count(3) == 7u  && f6[3] == 7u &&
              f6.count(4) ==11u  && f6[4] ==11u &&
              f6.count(5) ==13u  && f6[5] ==13u ))
        {
            fout << "\nERROR:  Factor<BigInt> failed on 24018350267611933650627567399079537500 = 2^2 3^3 5^5 7^7 11^11 13^13." << endl ;
            status = false ;

            fout << "Number of factors = " << f6.num() << endl ;
            fout << "Factors = " << endl ;
            for (int i = 0 ;  i < f6.num() ;  ++i)
                fout << f6[ i ] << " ^ " << f6.count( i ) << " " ;
            fout << endl ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  uint Euler Phi( 5^4 - 1 ) = 192" ;
        Factor<uint> f4( (uint)(5 * 5 * 5 * 5 - 1) ) ;
        uint phi = f4.EulerPhi() ;
        if ( phi != 192u)
        {
            fout << "\nERROR:  EulerPhi( 5^4 - 1) = " << phi << " failed." << endl ;
                status = false ;
        }
        else
            fout << ".........PASS!" ;

        fout << "\nTEST:  BigInt Euler Phi( 1234567890 ) = 329040288" ;
        Factor<BigInt> f5( (uint)(1234567890) ) ;
        BigInt phi2 = f5.EulerPhi() ;
        if (phi2 != static_cast<BigInt>(329040288u))
        {
            fout << "\nERROR:  EulerPhi( 1234567890 ) = " << phi2 << " failed." << endl ;
                status = false ;
        }
        else
            fout << ".........PASS!" ;


        ////////////////////////////////////////////////////////////////////////
        // Test polynomials
        ////////////////////////////////////////////////////////////////////////

		fout << "\n-------------------------> Testing polynomial arithmetic." ;

        fout << "\nTEST:  Polynomial() default constructor." ;
        try
        {
            Polynomial p ;
            if (p.deg() != 0)
            {
                fout << "\nERROR: Polynomial default constructor failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: Polynomial default constructor failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

        fout << "\nTEST:  Polynomial() from string." ;
        {
            Polynomial p( "2x^2 + 1, 3" ) ;
            if (p.deg() != 2 || p.modulus() != 3 || p[0] != 1 || p[1] != 0 || p[2] != 2)
            {
                fout << "\nERROR: Polynomial p( \"2x^2 + 1, 3\" ) failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

        fout << "\nTEST:  Polynomial() from string with negative constant." ;
        {
            Polynomial p( "x^4-1, 5" ) ;
            if (p.deg() != 4 || p.modulus() != 5 || p[0] != 4) // because -1 = 4 (mod 5)
            {
                fout << "\nERROR: Polynomial p( \"x^4-1, 5\" ) failed." << endl ;
                fout << " p = " << p << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

		fout << "\nTEST:  Polynomial() to string." ;
        try
        {
            Polynomial p ; // ( "2x^2 + 1, 3" )
			Polynomial q( p ) ;
			p[0] = 1 ; p[2] = 2 ;
			p.setModulus( 3 ) ;
			string s = p ;
            if (s != "2 x ^ 2  + 1, 3")
			{
                fout << "\nERROR: Polynomial p( \"2x^2 + 1, 3\" ) to string s = " << s << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  p( \"2x^2 + 1, 3\" ) failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: p( \"2x^2 + 1, 3\" ) failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  Polynomial() copy constructor." ;
        try
        {
            Polynomial p( "2x^2 + 1, 3" ) ;
			Polynomial q( p ) ;
            if (static_cast<string>(q) != "2 x ^ 2  + 1, 3")
			{
                fout << "\nERROR: Polynomial copy constructor p( q ) = " << q << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  copy constructor failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: copy constructor failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

	    fout << "\nTEST:  Polynomial assignment operator." ;
        try
        {
            Polynomial p( "2x^2 + 1, 3" ) ;
			Polynomial q ;
			q = p ; // If we did Polynomial q = p it would call the copy constructor.
            if (static_cast<string>(q) != "2 x ^ 2  + 1, 3")
			{
                fout << "\nERROR: Polynomial assignment operator p = q " << q << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  assignment operator p = q failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: assignment operator p = q failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

	    fout << "\nTEST:  Polynomial()[] lvalue operator." ;
        try
        {
            Polynomial p( "2x^2 + 1, 3" ) ;
			p[ 5 ] = 2 ;
			p[ 1 ] = 1 ;
            if (static_cast<string>(p) != "2 x ^ 5  + 2 x ^ 2  +  x  + 1, 3")
			{
                fout << "\nERROR: Polynomial [] operator " << p << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  [] failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: += failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

	    fout << "\nTEST:  Polynomial() + operator." ;
        try
        {
            Polynomial p( "2x^2 + 1, 3" ) ;
			Polynomial q( " x^2 + 1, 3" ) ;
			Polynomial r = p + q ;
            if (static_cast<string>(r) != "2, 3")
			{
                fout << "\nERROR: Polynomial += operator " << r << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  += failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: += failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  Polynomial * scalar" ;
        try
        {
            Polynomial p( "2x^2 + 1, 3" ) ;
			uint k = 2 ;

			Polynomial q = p * k ;
            if (static_cast<string>(q) != " x ^ 2  + 2, 3")
			{
                fout << "\nERROR: Polynomial *= operator " << q << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  += failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: += failed." << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }


	    fout << "\nTEST:  Polynomial evaluation x^4 + 3x + 3 (mod 5)" ;
        {
            Polynomial p( "x^4 + 3x + 3, 5" ) ;
			uint f2 = p( 2 ) ;
			uint f3 = p( 3 ) ;
			uint f0 = p( 0 ) ;
            if (f2 != 0 || f3 != 3 || f0 != 3)
			{
                fout << "\nERROR: Polynomial operator() = " << f2 << f3 << f0 << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

		fout << "\nTEST:  Polynomial evaluation x^4 + x + 1 (mod 2)" ;
        {
            Polynomial p( "x^4 + x + 1, 2" ) ;
			uint f0 = p( 0 ) ;
			uint f1 = p( 1 ) ;
            if (f0 != 1 || f1 != 1)
			{
                fout << "\nERROR: Polynomial operator() = " << f0 << f1 << " failed." << endl ;
                status = false ;
            }
            else
                fout << ".........PASS!" ;
        }

	    fout << "\nTEST:  Polynomial hasLinearFactor is true" ;
		{
            Polynomial p( "x^4 + 3x + 3, 5" ) ;
		    bool linFac = p.hasLinearFactor() ;
                if (linFac)
                    fout << ".........PASS!" ;
                else
                {
                    fout << "\nERROR: Polynomial hasLinearFactor = " << linFac << " failed." << endl ;
                    status = false ;
                }
			}

		fout << "\nTEST:  Polynomial hasLinearFactor is false" ;
        {
            Polynomial p( "x^4 + 3x^2 + x + 1, 5" ) ;
			bool linFac = p.hasLinearFactor() ;
            if (!linFac)
                fout << ".........PASS!" ;
            else
            {
                fout << "\nERROR: Polynomial hasLinearFactor = " << linFac << " failed." << endl ;
                status = false ;
            }
        }

	    fout << "\nTEST:  Polynomial isInteger" ;
        try
        {
            Polynomial p( "x^4 + 3x + 3, 5" ) ;
			bool isInt = p.isInteger() ;
            if (!isInt)
                fout << ".........PASS!" ;
            else
            {
                fout << "\nERROR: Polynomial " << p << " isInteger = " << isInt << " failed." << endl ;
                status = false ;
            }

            Polynomial q( "3, 5" ) ;
			isInt = q.isInteger() ;
            if (isInt)
                fout << ".........PASS!" ;
            else
            {
                fout << "\nERROR: Polynomial " << q << " isInteger = " << isInt << " failed." << endl ;
                status = false ;
            }

        }
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error: polynomial operator() failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: += failed." << endl ;
            status = false ;
            fout << "bad_exception:  " << e.what() << endl ;
        }

		fout << "\nTEST:  Polynomial initial and next trial polynomials" ;
		try {
		    Polynomial p ;
		    p.initial_trial_poly( 4, 5 ) ;

		    int i = 0 ;
		    for (i = 1 ;  i <= 20 ; ++i)
			{
		        p.next_trial_poly() ;
			}
			if (static_cast<string>(p) == " x ^ 4  + 3 x  + 4, 5")
                fout << ".........PASS!" ;
			else
            {
                fout << "\nERROR: Polynomial " << p << " (20th iteration from initial) failed." << endl ;
                status = false ;
            }

		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error: polynomial operator() failed." << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: += failed." << endl ;
            status = false ;
            fout << "bad_exception:  " << e.what() << endl ;
        }

		////////////////////////////////////////////////////////////////////////
        // Test polynomial mod
        ////////////////////////////////////////////////////////////////////////

		fout << "\n-------------------------> Testing polynomial modulo arithmetic." ;

		fout << "\nTEST:  PolyMod constructor from polynomials." ;
		try {
		    Polynomial g( "x^4 + x^2 + 1,2" ) ;
		    Polynomial f( "x^4 + x + 1,2" ) ;

			PolyMod p( g, f ) ;
			if (static_cast<string>(p) == " x ^ 2  +  x , 2" &&
			    static_cast<string>( p.getf() ) == " x ^ 4  +  x  + 1, 2" && p.getModulus() == 2)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod constructor from polynomials" << p << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyMod constructor from string and polynomial." ;
		try {
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
			PolyMod p( "x^6 + 2x^2 + 3x + 2, 5", f ) ;

			if (static_cast<string>(p) == "3 x ^ 3 , 5" &&
			    static_cast<string>( p.getf() ) == " x ^ 4  +  x ^ 2  + 2 x  + 3, 5" && p.getModulus() == 5)
                fout << ".........PASS!" ;
			else
			{
		        fout << "\nERROR: PolyMod constructor from string and polynomial failed." << endl ;
		        fout << "\ng(x) mod f(x), p = " << p << endl ;
		        fout << "\nf(x) = " << f << endl ;
                status = false ;
			}
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            status = false ;
            fout << "bad_exception:  " << e.what() << endl ;
        }

		fout << "\nTEST:  PolyMod timesX." ;
		try {
		    Polynomial g( "2x^3 + 4x^2 + 3x, 5" ) ;
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;

			PolyMod p( g, f ) ;

			p.timesX() ;
			if (static_cast<string>(p) == "4 x ^ 3  +  x ^ 2  +  x  + 4, 5" )
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod timesX " << p << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyMod autoconvolve." ;
		try {
		    Polynomial t( "4x^3 + x^2 + 3x + 3, 5" ) ;
			uint k = 3 ; uint lower = 1 ; uint upper = 3 ;
			uint c = autoConvolve( t, k, lower, upper ) ;

            if (c == 3)
			    fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod autoconvolve = " << c
                     << " for t = " << t << " deg = " << t.deg() << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError autoconvolve:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            status = false ;
            fout << "bad_exception:  " << e.what() << endl ;
        }

		fout << "\nTEST:  PolyMod convolve." ;
		try {
		    Polynomial s( "4x^3 + x^2 + 3x + 3, 5" ) ;
		    Polynomial t( "4x^3 + x^2 + 3x + 3, 5" ) ;
			uint k = 3 ; uint lower = 1 ; uint upper = 3 ;
			uint c = convolve( s, t, k, lower, upper ) ;

            if (c == 3)
			    fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod convolve = " << c
                     << " for t = " << t << " deg = " << t.deg() << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError convolve:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyMod coeffOfSquare." ;
		try {
		    Polynomial g( "4x^3 + x^2 + 3x + 3, 5" ) ;
            int n = 4 ;
			uint c0 = coeffOfSquare( g, 0, n ) ;
			uint c1 = coeffOfSquare( g, 1, n ) ;
			uint c2 = coeffOfSquare( g, 2, n ) ;
			uint c3 = coeffOfSquare( g, 3, n ) ;
			uint c4 = coeffOfSquare( g, 4, n ) ;
			uint c5 = coeffOfSquare( g, 5, n ) ;
			uint c6 = coeffOfSquare( g, 6, n ) ;

            if (c0 == 4 && c1 == 3 && c2 == 0 && c3 == 0 && c4 == 0 && c5 == 3 && c6 == 1)
			    fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod coeffOfSquare (c0 ... c6) = " << c0 << " " << c1 << " " << c2
                     << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError coeffOfSquare:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            status = false ;
            fout << "bad_exception:  " << e.what() << endl ;
        }

		fout << "\nTEST:  PolyMod coeffOfProduct." ;
		try {
		    Polynomial s( "4x^3 + x^2 + 4, 5" ) ;
		    Polynomial t( "3x^2 + x + 2, 5" ) ;
            int n = 4 ;
			uint c0 = coeffOfProduct( s, t, 0, n ) ;
			uint c1 = coeffOfProduct( s, t, 1, n ) ;
			uint c2 = coeffOfProduct( s, t, 2, n ) ;
			uint c3 = coeffOfProduct( s, t, 3, n ) ;
			uint c4 = coeffOfProduct( s, t, 4, n ) ;
			uint c5 = coeffOfProduct( s, t, 5, n ) ;
			uint c6 = coeffOfProduct( s, t, 6, n ) ;

            if (c0 == 3 && c1 == 4 && c2 == 4 && c3 == 4 && c4 == 2 && c5 == 2 && c6 == 0)
			    fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod coeffOfProduct (c0 ... c6) = " << c0 << " " << c1 << " " << c2
                     << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError coeffOfProduct:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            status = false ;
            fout << "bad_exception:  " << e.what() << endl ;
        }

		fout << "\nTEST:  PolyMod square." ;
		try {
		    Polynomial g( "4x^3 + x^2 + 4, 5" ) ;
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;

			PolyMod p( g, f ) ;

			p.square() ;
			if (static_cast<string>(p) == "2 x ^ 3  + 4 x ^ 2  +  x  + 1, 5" )
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod square " << p << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyMod operator* and implicitly, operator*=" ;
		try {
		    Polynomial s( "4x^3 + x^2 + 4, 5" ) ;
		    Polynomial t( "3x^2 + x + 2, 5" ) ;
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;

			PolyMod smodf( s, f ) ;
			PolyMod tmodf( t, f ) ;

			PolyMod p = smodf * tmodf ;
			if (static_cast<string>(p) == "2 x ^ 3  + 3 x ^ 2  + 4 x  + 2, 5" )
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod operator* " << p << " failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyMod x_to_power and isInteger()" ;
		try {
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
			PolyMod x( "x, 5", f ) ; // g(x) = 0, modulus = 5.

			PolyMod p = power( x, 156u ) ;
			if (static_cast<string>(p) == "3, 5" && p.isInteger())
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyMod x_to_power = |" << static_cast<string>(p) << "| failed." << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

        ////////////////////////////////////////////////////////////////////////
        // Test polynomial order
        ////////////////////////////////////////////////////////////////////////

		fout << "\nTEST:  PolyOrder reduced Q-I matrix" ;
		try {
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
			PolyOrder order( f ) ;

		    order.hasMultipleDistinctFactors() ;
			string s = order.printQMatrix() ;
			string t = "\n(    0   0   0   0 )\n(    0   4   0   0 )\n(    4   0   0   0 )\n(    0   0   4   0 )\n" ;
			if (s == t)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder reduced Q-I failed = " << s << endl ;
		        fout << "\n                   true reduced Q-I = " << t << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyOrder 3 distinct factors out of 4" ;
		try {
		    Polynomial f( "x^4 + 3 x^3 + 3 x^2 + 3 x + 2, 5" ) ;
			PolyOrder order( f ) ;

		    bool multipleFactors = order.hasMultipleDistinctFactors() ;
			if (multipleFactors && order.getNullity() == 3)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder 3 distinct factors out of 4 failed" << endl ;
				fout << " f( x ) = " << f << endl ;
				fout << " nullity = " << order.getNullity() << endl ;
                fout << "\n    reduced Q-I matrix " << order.printQMatrix() << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyOrder 2 distinct factors" ;
		try {
		    Polynomial f( "x^3 + 3, 5" ) ;
			PolyOrder order( f ) ;

		    bool multipleFactors = order.hasMultipleDistinctFactors() ;
			if (multipleFactors && order.getNullity() == 2)
                fout << ".........PASS!" ;
			else
            {
                status = false ;
		        fout << "\nERROR: PolyOrder 2 distinct factors failed" << endl ;
				fout << " f( x ) = " << f << endl ;
				fout << " nullity = " << order.getNullity() << endl ;
                fout << "\n    reduced Q-I matrix " << order.printQMatrix() << endl ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }

		fout << "\nTEST:  PolyOrder irreducible" ;
		try {
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
			PolyOrder order( f ) ;

		    bool multipleFactors = order.hasMultipleDistinctFactors() ;
			if (multipleFactors == false && order.getNullity() == 1)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder irreducible" << endl ;
				fout << " f( x ) = " << f << endl ;
				fout << " nullity = " << order.getNullity() << endl ;
                fout << "\n    reduced Q-I matrix " << order.printQMatrix() << endl ;
                status = false ;
            }
		}
	    catch( PolynomialRangeError & e )
        {
            fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
            status = false ;
        }
        catch( bad_exception & e )
        {
            fout << "\nERROR: " << endl ;
            fout << "bad_exception:  " << e.what() << endl ;
            status = false ;
        }


		fout << "\nTEST:  PolyOrder 1 distinct factor 4 times" ;
		{
		    Polynomial f( "x^4 + 4x^3 + x^2 + 4x + 1, 5" ) ;
			PolyOrder order( f ) ;

		    bool multipleFactors = order.hasMultipleDistinctFactors() ;
			if (multipleFactors == false && order.getNullity() == 1)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder 1 distinct factor 4 times" << endl ;
				fout << " f( x ) = " << f << endl ;
				fout << " nullity = " << order.getNullity() << endl ;
                fout << "\n    reduced Q-I matrix " << order.printQMatrix() << endl ;
                status = false ;
            }
		}

		fout << "\nTEST:  PolyOrder order_m()" ;
		{
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
			PolyOrder order( f ) ;

			if (order.order_m())
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder order_m failed." << endl ;
                status = false ;
            }
		}

		fout << "\nTEST:  PolyOrder order_r() is true" ;
		{
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
			PolyOrder order( f ) ;

			uint a = order.order_r() ;
			if (a == 3)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder order_r failed." << endl ;
                status = false ;
            }
		}

		fout << "\nTEST:  PolyOrder order_r() is false" ;
		{
		    Polynomial f( "x^4 + x + 3, 5" ) ;
			PolyOrder order( f ) ;

			uint a = order.order_r() ;
			if (a == 0)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder order_r failed." << endl ;
                status = false ;
            }
		}

		fout << "\nTEST:  PolyOrder isPrimitive on non-primitive poly" ;
		{
		    Polynomial f( "x^5 + x + 1, 2" ) ;
			PolyOrder order( f ) ;

		    bool isPrimitive = order.isPrimitive() ;
			if (!isPrimitive)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder isPrimitive on non-primitive poly" << endl ;
				fout << " f( x ) = " << f << endl ;
                status = false ;
            }
		}

		fout << "\nTEST:  PolyOrder isPrimitive on primitive poly" ;
		{
		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
			PolyOrder order( f ) ;

		    bool isPrimitive = order.isPrimitive() ;
			if (isPrimitive)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder isPrimitive on primitive poly" << endl ;
				fout << " f( x ) = " << f << endl ;
                status = false ;
            }
		}

		fout << "\nTEST:  PolyOrder isPrimitive on primitive poly, part II" ;
		{
		    Polynomial f0( "x^4-1, 5" ) ;
			PolyOrder order( f0 ) ;

		    Polynomial f( "x^4 + x^2 + 2x + 3, 5" ) ;
            order.newPolynomial( f ) ;

		    bool isPrimitive = order.isPrimitive() ;
			if (isPrimitive)
                fout << ".........PASS!" ;
			else
            {
		        fout << "\nERROR: PolyOrder isPrimitive on primitive poly, part II" << endl ;
				fout << " f( x ) = " << f << endl ;
                status = false ;
            }
		}

    // These catch blocks should be in order from more specific to less specific.
	// i.e. up the class hierarchy.
    }
    catch( PrimpolyError & e )
    {
        fout << "\nERROR:  Primpoly error:  " << e.what() << endl ;
        status = false ;
    }
    catch( ParserError & e )
    {
        fout << "\nERROR:  Parser error:  " << e.what() << endl ;
        status = false ;
    }
   catch( PolynomialRangeError & e )
    {
        fout << "\nERROR:  PolynomialRangeError error:  " << e.what() << endl ;
        status = false ;
    }
   catch( PolynomialError & e )
    {
        fout << "\nERROR:  PolynomialError error:  " << e.what() << endl ;
        status = false ;
    }
    catch( BigIntOverflow & e )
    {
        fout << "\nERROR:  BigIntOverflow:  " << e.what() << endl ;
        status = false ;
    }
    catch( BigIntMathError & e )
    {
        fout << "\nERROR:  BigIntMathError:  " << e.what() << endl ;
        status = false ;
    }
    catch( bad_alloc & e )
    {
        fout << "\nERROR:  Caught bad_alloc exception:  " << e.what() << endl ;
        status = false ;
    }
    // Catch all uncaught exceptions, which would otherwise call terminate()
    // which would call abort()
    // Can't catch those thrown by initializing or
    // destructing global variables.
    // Also terminate() will be called if exception mechanism finds a corrupt
    // stack or catches a destructor throwing an exception.
    catch( bad_exception & e )
    {
           fout << "\nERROR:  Internal error.  bad_exeception. " << e.what() << endl ;
           status = false ;
    }
    //  Standard library exceptions.
    catch( exception & e )
    {
        fout << "\nERROR:  Standard library error: " << e.what() << endl ;
        status = false ;
    }
    catch( ... )
    {
           fout << "\nERROR:  uncaught exception" << endl ;
           status = false ;
    }

    fout << "\nDone!" << endl ;
    fout.close() ;

    return status ;
}