/*==============================================================================
|
|  File Name:
|
|     ppFactor.cpp
|
|  Description:
|
|     Methods for factoring.
|
|  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 <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.


/*------------------------------------------------------------------------------
|                                Local Functions                               |
------------------------------------------------------------------------------*/

// Generic function.
template <typename IntType>
static
bool factorTable( const IntType &     n,
                  vector< PrimeFactor< IntType > > & factor,
                  uint &              numFactors )
{
    return false ;
}



// Specialization for unsigned integer.
template<>
static
bool factorTable< uint >( const uint &      n,
                          vector< PrimeFactor< uint > > &  factor,
                          uint &            numFactors )
{
    if (n == 1u)
    {
	    PrimeFactor<uint> f( 1u, 1u ) ;
        factor.push_back( f ) ;

        numFactors = 1u ;
        return true ;
    }

    return false ;
}



// Specialization for BigInt.
template<>
static
bool factorTable< BigInt >( const BigInt &      n,
                            vector< PrimeFactor<BigInt> > &  factor,
                            uint &              numFactors )
{
    if (n == static_cast<BigInt>(1u))
    {
	    PrimeFactor<BigInt> f( 1u, 1u ) ;
        factor.push_back( f ) ;

        numFactors = 1u ;
        return true ;
    }
    else if (n == static_cast<BigInt>("4611686018427387903")) /*  (2^62 - 1)/(2 - 1)  */
    {
	    PrimeFactor<BigInt> f1( 3u, 1u ) ;
		factor.push_back( f1 ) ;

	    PrimeFactor<BigInt> f2( static_cast<BigInt>("715827883"), 1u ) ;
		factor.push_back( f2 ) ;

	    PrimeFactor<BigInt> f3( static_cast<BigInt>("2147483647"), 1u ) ;
		factor.push_back( f3 ) ;

        numFactors = 3u ;
        return true ;
    }

    return false ;
}

template <typename IntType>
Factor<IntType>::Factor( const IntType n_arg )
{
    //  Handle special cases for speed.
    //  These are cases where r = (p^n - 1)/(p-1) has large prime factors.
    if (n_arg <= 0u || factorTable( n_arg, factor_, numFactors_ ))
        return ;
#if 1
    IntType crossoverPoint = 13693952u ;  // TODO this is largest uint type. Need 2^60 instead.
    if ( n_arg < crossoverPoint || !PollardRho( n_arg ))
    {
        // Tried and true, but slow.
        trialDivision( n_arg ) ;
    } 
#else
        trialDivision( n_arg ) ;
#endif
}


/*=============================================================================
|
| NAME
|
|    template <typename IntType> Factor<IntType>::Factor( const IntType n_arg )
|
| DESCRIPTION
|
     Factor n into primes.  Record all the distinct prime factors and how
     many times each occurs.  Return the number of primes - 1.
|    where n_arg >= 1
|
| EXAMPLE CALLING SEQUENCE
|
|    Factor<BigInt> f ;


     primes (bigint *)   List of distinct prime factors.
     count  (uint *)     List of how many times each factor occurred.

     When n = 1, t = 0, and primes[ 0 ] = count[ 0 ] = 1.

RETURNS

     t (int)             Number of prime factors - 1.
	                     Prime factors and their multiplicites are in locations 0 to t.

EXAMPLE
                   2
    For n = 156 = 2  * 3 * 13 we have

    k   primes[ k ]   count[ k ]
    ----------------------------
    0        2            2
    1        3            1
    2       13            1

METHOD

    Method described by D.E. Knuth, ART OF COMPUTER PROGRAMMING, vol. 2, 2nd ed.,
	                                                            -----
	in Algorithm A, pgs. 364-365.  The running time is O( max \/ p   , p  )
                                                                  t-1   t
    where p  is the largest prime divisor of n and p    is the next largest.
	       t                                        t-1

    (1)  First divide out all multiples of 2 and 3 and count them.

	(2)  Next, divide n by all integers d >= 5 except multiples of 2 and 3.

    (3)  Halt either when all prime factors have been divided out (leaving n = 1)
	     or when the current value of n is prime.  The stopping test

          (d > | n/d | AND r != 0)
               --   --
         detects when n is prime.

    In the example above, we divided out 2's and 3's first, leaving
    n = 13.  We continued with trial divisor d = 3.  But now the stopping
    test was activated because 5 > | 13/5 | = 2, and remainder = 3 (non-zero).
                                   --    --
    n = 13 itself is the last prime factor of 156.  We return t = 2.  There
    are 2 + 1 = 3 distinct primes.


    If we start with d = 5, and add 2 and 4 in succession, we will run through all
	the integers except multiples of 2 and 3.  To see this, observe the pattern:

    Integers:     1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17
	Mult. of 2:      x     x     x     x     x     x     x     x
	Mult. of 3:         x        x        x        x        x
	d:                        x     x           x     x          x
	Separation:                  2       4         2       4

    The sequence of divisors d, includes all the primes, so is safe to use for
	factorization testing.

    Theorem.  Suppose no divisor d in the sequence divides n.  Suppose at some point,
	q < d with r != 0.  Then n must be prime.

    Proof.  We begin with an integer n which has all powers of 2 and 3 removed.
    Assume, to the contrary, that n is composite:  n = p ... p
	      t                                             1     t
    n >= p    since p  is the smallest prime factor.
	      1          1

	      2
    n >= p    since n is composite, so has at least 2 prime factors.
	      1
		        2
	  >= (d + 1)   since p  > d implies p >= (d + 1).  p > d because we assumed
                          1              1              1
      that no d in the sequence divided n, so d couldn't be any of the prime divisors
	  p
	   i
         2            2               2
	  = d + 2d + 1 = d + d + d + 1 > d + d
         2
      > d + n mod d   since n mod d < d
                                                                              2
      > | n / d | d + n mod d  because our test said q < d, so | n / d | d < d
	    --     --                                              --     --

      = n.  So we get the contradiction n > n.  Thus n must be prime.
	                                               ---
    Note that this is sharper than the bound d > \/ n

    Theorem.  Conversely, if n is a prime, no d will divide it, so at some point,
	d will be large enough that q = | n / d | < d.  r != 0 since n is prime.
	                                --     --

    Theorem.  The factoring algorithm works.

	Proof.  If n == 1 we exit immediately.  If n is composite, we divide out all powers
	of 2 and 3 (if any).  Since d runs through all possible prime divisors, we divide
	these out also.  Composite trial d causes no problem;  prime factors of d are divided
	out of n before we try to divide by d, so such a d cannot divide n.

    We end in one of two ways (1) All divisors of n have been divided out in which case n = 1.
	(2) n is a prime, so the stopping test is activiated and n != 1 is recorded as a prime
	divisor.


BUGS

    Can be slow when n is a prime.  We could do a probabilistic test for
    the primality of n at the statement, "n_is_prime = (r != 0 && q < d) ;"
    which would speed up the test.

+============================================================================*/

template <typename IntType>
void Factor<IntType>::trialDivision( const IntType & n_arg )
{
    IntType n = n_arg ;

    // TODO clean this up.
    factor_.clear() ;


// Remove factors of 2.
uint cnt = 0u ;
while( n % 2u == 0u )
{
    n /= 2u ;
    ++cnt ;
}

if (cnt != 0u)
{
    PrimeFactor<IntType> f( 2u, cnt ) ;
    factor_.push_back( f ) ;
}

#ifdef DEBUG_PP_FACTOR
    cout << "after removing powers of 2, n =" <<  n << endl ;
#endif

// Remove factors of 3.
cnt = 0 ;
while( n % 3u == 0u )
{
    n /= 3 ;
    ++cnt ;
}

if (cnt != 0)
{
    PrimeFactor<IntType> f( 3u, cnt ) ;
	factor_.push_back( f ) ;
}

#ifdef DEBUG_PP_FACTOR
cout << "after removing powers of 3, n =" <<  n << endl ;
#endif


//  Factor the rest of n.  Continue until n = 1 (all factors divided out)
//  or until n is prime (so n itself is the last prime factor).

bool new_d = true ;       //  true if current divisor is different from
                          //  the previous one.
bool n_is_prime = false ; //  true if current value of n is prime.
bool flag = true ;        //  Alternately true and false.

IntType d = 5u ;  // First trial divisor.
do {
	/*  Integer divide to get quotient and remainder:  n = qd + r */
    IntType q = n / d ;
    IntType r = n % d ;

    #ifdef DEBUG_PP_FACTOR
    cout << "n = " ; printNumber( n ) ;
    cout << "d = " ; printNumber( d ) ;
    cout << "q = " ; printNumber( q ) ;
    cout << "r = " ; printNumber( r ) ;
    #endif

    // Stopping test.
    n_is_prime = (r != 0u && q < d) ;


    // d | n.
    if (r == 0u)
    {
        n = q ;     /*  Divide d out of n. */

        if (new_d)  /*  New prime factor. */
        {
		    PrimeFactor<IntType> f( d, 1u ) ;
			factor_.push_back( f ) ;

            new_d = false ;
            #ifdef DEBUG_PP_FACTOR
                cout << "after dividing out new prime divisor d = " << d << " n =" << n << endl ;
            #endif
         }
         else
         {
             ++factor_[ factor_.size() - 1 ].count_ ; //  Same divisor again.  Increment its count.

            #ifdef DEBUG_PP_FACTOR
                cout << "after dividing out repeated factor d = " << d << " n =" << n << endl ;
            #endif
         }
    }
    else {
        d += (flag ? 2u : 4u) ;  /* d did not divide n.  Try a new divisor. */
        flag = !flag ;
        new_d = true ;

        #ifdef DEBUG_PP_FACTOR
            cout << "next trial divisor d  = " << d << endl ;
        #endif
    }

} while ( ! n_is_prime && (n != 1u)) ;

if (n == 1u)       /*  All factors were divided out. */
{
    numFactors_ = factor_.size() ;
    return ;
}
else
{  //  Current value of n is prime.  It is the last prime factor.
    PrimeFactor<IntType> f( n, 1u ) ;
    factor_.push_back( f ) ;
    numFactors_ = factor_.size() ;
    return ;
}

} /* ========================= end of function factor ======================= */




template <typename IntType>
Factor<IntType>::~Factor()
{
   // Do nothing.
}


template <typename IntType>
Factor<IntType>::Factor( const Factor<IntType> & f )
		: factor_( f.factor_ )
		, numFactors_( f.numFactors_ )
{
}

// Assignment operator.
template <typename IntType>
Factor<IntType> & Factor<IntType>::operator=( const Factor<IntType> & g )
{
	// If assigned to itself, return the object g.
    if (this == &g)
        return *this ;
		
	// TODO make a safe copy.
	factor_     = g.factor_ ;
	numFactors_ = g.numFactors_ ;
	
    // Return a reference to assigned object to make chaining possible.
    return *this ;
}


template <typename IntType>
IntType Factor<IntType>::operator[]( uint i ) const
{
    return factor_[ i ].prime_ ;
}


// Return the number of distinct factors.
template <typename IntType>
uint Factor<IntType>::num() const
{
    return numFactors_ ;
}


// Return the count for the ith prime factor.
// TODO range check.
template <typename IntType>
IntType Factor<IntType>::count( uint i ) const
{
    return factor_[ i ].count_ ;
}



/*
      4
Phi[ 5  - 1 ]

   4              4
  5  - 1 = 624 = 2  3  13

  Phi = 624 (1 - 1/2) (1 - 1/3)(1 - 1/13) = 192



  Handle special cases.

  Warning: can take a long time due to factorization.
  */

template <typename IntType>
IntType Factor<IntType>::EulerPhi()
{
	// Compute Euler phi[ n ] =
	//
	//      num prime factors
	//        -------
	//     n   |  |   (1 - 1/p )
	//         i = 1          i
    //
    IntType phi = 1u ;
	for (uint i = 0u ;  i < numFactors_ ; ++i)
        for (uint j = 0u ;  j < factor_[ i ].count_ ; ++j)
		    phi *= factor_[ i ].prime_ ;

	for (uint i = 0u ;  i < numFactors_ ; ++i)
	{
		phi *= (factor_[ i ].prime_ - static_cast<IntType>( 1u )) ;
        phi /=  factor_[ i ].prime_ ;
	}

	return phi ;
}


/*==============================================================================
|                             is_probably_prime                                |
================================================================================

DESCRIPTION

     Test whether the number n is probably prime.  If n is composite
     we are correct always.  If n is prime, we will be wrong no more
     than about 1/4 of the time on average, probably less for
     any x and n.

INPUT

     n      (int, n >= 0)    Number to test for primality.
     x      (int, 1 < x < n for n > 6) A random integer.

RETURNS

     true    If n is probably prime.

     false   If n is definitely not prime or n < 0,
             or x is out of range.


METHOD

     Miller-Rabin probabilistic primality test, described by Knuth.

     If
                   k
         n = 1 + 2  q

     is prime and
           q
          x  mod n != 1

     the sequence
                                         k
            q          2q               2 q
         { x  mod n, x    mod n, ..., x     mod n}

     will end with 1 and the element in the sequence just before
     1 first appears = n-1, since

          2
         y = 1 (mod p)

     satisfies y = +1 or y = -1 only. The probability the algorithm fails is
     bounded above by about 1/4 for all n.

EXAMPLE

     If k = 5, q = 3 then n = 1 + 2^5 3 = 97, which is prime.  
     
     For x = 10, we get                    For x = 9, we get

      q                                   q  
     x      =   30                       x     =  50

         1                                   1
      q 2                                 q 2
     x      =   27                       x     = 75

         2                                   2
      q 2                                 q 2
     x      =   50                       x     =  96

         3                                   3
      q 2                                 q 2
     x      =   75                       x     = 1

         4                                   4
      q 2                                 q 2
     x      =   96                       x     = 1

         5                                   5
      q 2                                 q 2
     x      =    1                       x     =  1


     If k = 4, q = 3 then n = 1 + 2^4 3 = 49, which is composite.  
     
     For x = 10, we get

      q
     x      =   20

         1
      q 2
     x      =   8

         2
      q 2
     x      =   15

         3       
      q 2       
     x      =  29

         4      
      q 2      
     x      =  8
     

BUGS

    None.

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

template <typename IntType>
bool is_probably_prime( const IntType & n, const IntType & x )
{
/*  Not prime if input is out of range. */
if (n < static_cast<IntType>( 0u ))
    return false ;

/*  Handle special cases. */
if (n == static_cast<IntType>( 0u ) || n == static_cast<IntType>( 1u ) || n == static_cast<IntType>( 4u ))
    return false ;

if (n == static_cast<IntType>( 2u ) || n == static_cast<IntType>( 3u ) || n == static_cast<IntType>( 5u ))
    return true ;

/* Even numbers aren't prime. */
if (n % static_cast<IntType>( 2u ) == static_cast<IntType>( 0u ))
    return false ;

/* Return not prime if x is out of range. */
if (x <= static_cast<IntType>( 1u ) || x >= n)
    return false ;

// Factor out powers of 2 to get n = 1 + 2^k q, q odd.
IntType nm1 = n - static_cast<IntType>( 1u ) ;
IntType k   = static_cast<IntType>( 0u ) ;

while( nm1 % static_cast<IntType>( 2u ) == static_cast<IntType>( 0u ))
{
	nm1 /= static_cast<IntType>( 2u ) ;
	++k ;
}

IntType q = nm1 ;

#ifdef PP_PRIMALITY_DEBUG
cout << "Miller-Rabin primality test" << endl ;
cout << "    q = " << q << endl ;
#endif

//  y = x^q (mod n)
PowerMod<IntType> power_mod( n ) ;
IntType y = power_mod( x, q ) ;

#ifdef PP_PRIMALITY_DEBUG
cout << "    x^q = " << y << endl ;
#endif

for (IntType j = 0u ;  j < k ;  ++j)
{
    if ( (j == static_cast<IntType>( 0u ) && y == static_cast<IntType>( 1u )) || (y == n - static_cast<IntType>( 1u )))
		return true ;
	else if (j > static_cast<IntType>( 0u ) && y == static_cast<IntType>( 1 ))
		return false ;

	// Compute y^2 (mod n)
    y = power_mod( y, static_cast<IntType>( 2u ) ) ;

    #ifdef PP_PRIMALITY_DEBUG
    cout << "    x ^ (2^ " << j << " j q) = " << y << endl ;
    #endif
}

/*  Shouldn't get here, but do it anyway for safety. */
return false ;

} /* =============== end of function is_probably_prime ===================== */



/*==============================================================================
|                           is_almost_surely_prime                             |
================================================================================

DESCRIPTION

     Test whether the number n is almost surely prime.  If n is
     composite, we always return NO.  If n is prime, the
     probability of returning YES successfully is

                 NUM_PRIME_TEST_TRIALS
        1 - (1/4)

INPUT

     n      (int, n >= 0)

RETURNS


     YES    If n is probably prime with a very high degree of
            probability.
     NO     If n is definitely not prime.
	 -1     If inputs are out of range.

EXAMPLE

METHOD

    Use the Miller-Rabin probabilitic primality test for lots of
    random integers x.  If the test shows n is composite for any
    given x, is is non-prime for sure.

    If it passes the the primality test for several random values
    of x, the probability is it prime is approximately

         NUM_PRIME_TEST_TRIALS
    1 - 2

    assuming the random number generator is really random.

BUGS

    The system pseudorandom number generator needs to be a good one
    (i.e. passes the spectral test and other statistical tests).

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

template <typename IntType>
bool is_almost_surely_prime( const IntType & n )
{

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

IntType
    trial = 0u,
    x     = 3u ;

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

/* Seed the random-number generator. */
srand( 314159u );

for (trial = static_cast<IntType>( 1u ) ;  trial <= static_cast<IntType>( NUM_PRIME_TEST_TRIALS ) ;  ++trial)
{
    /*  Generate a new random integer such that 1 < x < n. */
    x = static_cast<IntType>( rand() ) % n ;
	
	// Clip away from 1.
    if (x <= static_cast<IntType>( 1u )) 
	    x = static_cast<IntType>( 3u ) ;

	/* Definitely not prime. */
	if (is_probably_prime( n, x ) == false)
		return false ;
}

/* Almost surely prime because it passed the probably prime test
 * above many times.
 */
return true ;

} // ============== end of function is_almost_surely_prime ================





/*==============================================================================
|                                    skip_test                                 |
================================================================================

DESCRIPTION

     Make the test p  | (p - 1).
                    i

INPUT

    i      (int)       ith prime divisor of r.

    primes (bigint *)  Array of distict prime factors of r.  primes[ i ] = p
                                                                            i
    p      (int)        p >= 2.

RETURNS

    true    if the test succeeds for some i.
    false   otherwise

EXAMPLE

    Suppose i = 0, primes[ 0 ] = 2 and p = 5.  Return true since 2 | 5-1.

METHOD

    Test if (p - 1) mod p = 0 for all i.
                         i
BUGS

    None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/
template <typename IntType>
bool Factor<IntType>::skip_test( int p, int i ) const
{

// p  cannot divide the smaller number (p-1)
//  i
if ( static_cast<IntType>( p-1 ) < static_cast<IntType>( factor_[ i ].prime_ ) )
    return false ;
else
    return( ( static_cast<IntType>(p - 1) % static_cast<IntType>( factor_[ i ].prime_ ) == static_cast<IntType>( 0u )) ? true : false ) ;

} // ====================== end of function skip_test =======================


template <typename IntType>
bool Factor<IntType>::PollardRho( const IntType & n, const IntType & c )
{
    IntType x  = static_cast<IntType>( 5u ) ;
    IntType xp = static_cast<IntType>( 2u ) ;
    IntType k  = static_cast<IntType>( 1u ) ;
    IntType l  = static_cast<IntType>( 1u ) ;
    IntType n1 = n ;
    IntType g = static_cast<IntType>( 0u ) ;

    bool status = true ;

    factor_.clear() ;

    // Seed 1 as a factor and early out if n == 1.
    PrimeFactor<IntType>f( 1u, 1u ) ;
    factor_.push_back( f ) ;
    numFactors_ = 1u ;

    if (n == 1u)
        return status ;
	
#ifdef DEBUG_PP_FACTOR
        cout << "Pollard rho n = " << n << " c = " << c << endl ;
#endif

    while( true )
    {
        if (is_almost_surely_prime( n1 ))
        {
#ifdef DEBUG_PP_FACTOR
            cout << "    >>>>>prime factor n1 = " << n1 << endl ;
#endif
            PrimeFactor<IntType>f( n1, 1u ) ;
            factor_.push_back( f  ) ;
            ++numFactors_ ;
            status = true ;
            goto Exit ;
        }

        while (true)
        {
            IntType absDiff = (xp > x) ? (xp - x) : (x - xp) ;
            g = gcd( absDiff, n1 ) ;

#ifdef DEBUG_PP_FACTOR
            cout << "    inner while loop, gcd = g = " << g << " |x -xp| = " << absDiff << " n1 = " << n1
                 << " k = " << k << " l = " << l << endl ;
#endif
            if  (g == 1u)
            {
                --k ;
                if (k == 0u)
                {
                    xp = x ;
                    l = l * 2u ;
                    k = l ;
                }
                x = (x * x + c) % n1 ;
            } 
            else if (g == n1)
            {
#ifdef DEBUG_PP_FACTOR
                cout << "    failure:  g equals non-prime n1 = " << n1 << endl ;
#endif
                status = false ;
                goto Exit ;
            }
            else if (is_almost_surely_prime( g ))
            {
#ifdef DEBUG_PP_FACTOR
                cout << "    >>>>>prime factor g = " << g << endl ;
#endif
                PrimeFactor<IntType> f( g, 1u ) ;
                factor_.push_back( f ) ;
                ++numFactors_ ;
            }
            else
            {
#ifdef DEBUG_PP_FACTOR
                cout << "    g = " << g << " isn't prime " << endl ;
#endif
                status = false ;
                goto Exit ;
            }

            n1 /= g ;
            x = x % n1 ;
            xp = xp % n1 ;

#ifdef DEBUG_PP_FACTOR
            cout << "    after division, n1 = " << n1 << " x = " << x << " xp = " << xp << endl ;
#endif
            break ;
        } // end gcd loop
    }


Exit:
#ifdef DEBUG_PP_FACTOR
    cout << "     Factors from Pollard, possibly duplicated." << endl ;
    for( typename vector< PrimeFactor<IntType> >::const_iterator i = 
         factor_.begin() ;  i < factor_.end() ;  ++i)
    {
        PrimeFactor<IntType> e = *i ;
        cout << e.prime_ << " " << e.count_ << endl ;
    }
#endif
    // Sort and remove duplicates and units.
    sort( factor_.begin(), factor_.end(), CompareFactor<IntType>() ) ;

    for (int i = 1 ;  i < factor_.size() ;  ++i)
    {
        // This term is a duplicate of the last.
        if (factor_[ i ].prime_ == factor_[ i-1 ].prime_)
        {
            // Consolidate.
            factor_[ i ].count_ += factor_[ i-1 ].count_ ;
            factor_[ i-1 ].count_ = static_cast<IntType>( 0u ) ;
        }
    }

    typename vector< PrimeFactor<IntType> >::iterator last =
        remove_if( factor_.begin(), factor_.end(), Unit<IntType>() ) ;
    factor_.erase( last, factor_.end() ) ;
    numFactors_ = factor_.size() ;

#ifdef DEBUG_PP_FACTOR
    numFactors_ = factor_.size() ;
    cout << numFactors_ << " sorted unique factors" << endl ;
    for (int i = 0 ;  i < numFactors_ ;  ++i)
        cout << factor_[ i ].prime_ << " ^ " << factor_[ i ].count_ << endl ;
#endif
	
    return status ;
}

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

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

template Factor<uint>::Factor( const uint ) ;
template Factor<BigInt>::Factor( const BigInt ) ;

template Factor<uint>::~Factor() ;
template Factor<BigInt>::~Factor() ;

template Factor<uint>::Factor( const Factor<uint> & f ) ;
template Factor<BigInt>::Factor( const Factor<BigInt> & f ) ;
		
template Factor<uint> & Factor<uint>::operator=( const Factor<uint> & g ) ;
template Factor<BigInt> & Factor<BigInt>::operator=( const Factor<BigInt> & g ) ;

template void Factor<uint>::trialDivision( const uint & n ) ;
template void Factor<BigInt>::trialDivision( const BigInt & n ) ;

template bool Factor<uint>::PollardRho( const uint & n, const uint & c ) ;
template bool Factor<BigInt>::PollardRho( const BigInt & n, const BigInt & c ) ;

template uint Factor<uint>::num() const ;
template uint Factor<BigInt>::num() const ;

template uint Factor<uint>::operator[]( uint ) const ;
template BigInt Factor<BigInt>::operator[]( uint ) const ;

template uint Factor<uint>::count( uint ) const ;
template BigInt Factor<BigInt>::count( uint ) const ;

template uint Factor<uint>::EulerPhi() ;
template BigInt Factor<BigInt>::EulerPhi() ;

template bool is_probably_prime( const uint &, const uint & ) ;
template bool is_probably_prime( const BigInt &, const BigInt & ) ;

template bool is_almost_surely_prime( const uint & ) ;
template bool is_almost_surely_prime( const BigInt & ) ;

template bool Factor<int>::skip_test( int p, int i ) const ;
template bool Factor<BigInt>::skip_test( int p, int i ) const ;