1 /*==============================================================================
   2 |
   3 |  NAME
   4 |
   5 |     ppFactor.cpp
   6 |
   7 |  DESCRIPTION
   8 |
   9 |     Classes for factoring integers.
  10 |
  11 |     User manual and technical documentation are described in detail in my web page at
  12 |     http://seanerikoconnor.freeservers.com/Mathematics/AbstractAlgebra/PrimitivePolynomials/overview.html
  13 |
  14 |  LEGAL
  15 |
  16 |     Primpoly Version 16.1 - A Program for Computing Primitive Polynomials.
  17 |     Copyright (C) 1999-2021 by Sean Erik O'Connor.  All Rights Reserved.
  18 |
  19 |     This program is free software: you can redistribute it and/or modify
  20 |     it under the terms of the GNU General Public License as published by
  21 |     the Free Software Foundation, either version 3 of the License, or
  22 |     (at your option) any later version.
  23 |
  24 |     This program is distributed in the hope that it will be useful,
  25 |     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26 |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  27 |     GNU General Public License for more details.
  28 |
  29 |     You should have received a copy of the GNU General Public License
  30 |     along with this program.  If not, see http://www.gnu.org/licenses/.
  31 |     
  32 |     The author's address is seanerikoconnor!AT!gmail!DOT!com
  33 |     with !DOT! replaced by . and the !AT! replaced by @
  34 |
  35 ==============================================================================*/
  36 
  37 /*------------------------------------------------------------------------------
  38 |                                Include Files                                 |
  39 ------------------------------------------------------------------------------*/
  40 
  41 #include <cstdlib>      // abort()
  42 #include <iostream>     // Basic stream I/O.
  43 #include <new>          // set_new_handler()
  44 #include <cmath>        // Basic math functions e.g. sqrt()
  45 #include <complex>      // Complex data type and operations.
  46 #include <fstream>      // File stream I/O.
  47 #include <sstream>      // String stream I/O.
  48 #include <vector>       // STL vector class.
  49 #include <string>       // STL string class.
  50 #include <algorithm>    // Iterators, sorting, merging, union.
  51 #include <stdexcept>    // Exceptions.
  52 #include <cassert>      // assert()
  53 #include <regex>        // Regular expressions.
  54 
  55 using namespace std ;
  56 
  57 // Current working directory, recursive dir search.
  58 #ifdef USE_EXPERIMENTAL_FILESYSTEM_GNU_CPP
  59     #include <experimental/filesystem>
  60     using namespace experimental::filesystem ;
  61 #else
  62     #include <filesystem>
  63     using namespace filesystem ;
  64 #endif
  65 /*------------------------------------------------------------------------------
  66 |                                PP Include Files                              |
  67 ------------------------------------------------------------------------------*/
  68 
  69 #include "Primpoly.h"         // Global functions.
  70 #include "ppArith.h"          // Basic arithmetic functions.
  71 #include "ppBigInt.h"         // Arbitrary precision integer arithmetic.
  72 #include "ppOperationCount.h" // OperationCount collection for factoring and poly finding.
  73 #include "ppFactor.h"         // Prime factorization and Euler Phi.
  74 #include "ppPolynomial.h"     // Polynomial operations and mod polynomial operations.
  75 #include "ppParser.h"         // Parsing of polynomials and I/O services.
  76 #include "ppUnitTest.h"       // Complete unit test.
  77 
  78 /*=============================================================================
  79  |
  80  | NAME
  81  |
  82  |    Factorization<IntType>::factorTable( ppuint p, ppuint n )
  83  |
  84  | DESCRIPTION
  85  |                                                             n
  86  |    Table lookup for the prime factorization of the integer p
  87  |    False means the integer wasn't found in any of the factorization tables,
  88  |    i.e. p or n was too large for the tables.
  89  |
  90  |    Throw FactorError if the expected table files can't be
  91  |    located, or the factorization bad:  factors aren't primes or product of factors
  92  |    doesn't equal the integer.
  93  |
  94  +============================================================================*/
  95 
  96 template <typename IntType>
  97 bool Factorization<IntType>::factorTable( ppuint p, ppuint n )
  98 {
  99     // Clear out the factorization.
 100     factor_.clear() ;
 101 
 102     // List of factor table names for each prime p.
 103     vector<string> factorTableName
 104     {
 105         "",
 106         "",
 107         "c02minus.txt", // prime p = 2
 108         "c03minus.txt",
 109         "",             // p = 4 isn't a prime, so no table for it.
 110         "c05minus.txt",
 111         "c06minus.txt",
 112         "c07minus.txt",
 113         "",
 114         "",
 115         "c10minus.txt",
 116         "c11minus.txt",
 117         "c12minus.txt"
 118     } ;
 119 
 120     // Check if p is in the range of any of the above tables.  If not, return immediately.
 121     if (p > factorTableName.size() - 1 || factorTableName[ p ].length() == 0)
 122         return false ;
 123 
 124     // All the factor tables should be in the current working directory (the location of the executable) or in some subdirectory.
 125     const string factorTableFileExtension = ".txt" ;
 126     string factorTableFilePath = "" ; // If file isn't found, we'll error out later when we try to open the file.
 127 
 128     for (auto & fileOrDir : recursive_directory_iterator( PolyParser<PolySymbol, PolyValue>::current_working_dir_ ))
 129     {
 130         #ifdef DEBUG_PP_FACTOR
 131         cout << "Searching path = " << fileOrDir.path() << endl ;
 132         cout << "    File or dir name = " << fileOrDir.path().filename() << endl ;
 133         #ifdef USE_EXPERIMENTAL_FILESYSTEM_GNU_CPP
 134         {
 135         path pathObj( fileOrDir ) ;
 136         cout << "    Regular file = " << is_regular_file( pathObj ) << endl ;
 137         }
 138         #else
 139         cout << "    Regular file = " << fileOrDir.is_regular_file() << endl ;
 140         #endif
 141         cout << "    Extension = " << fileOrDir.path().extension() << endl ;
 142         cout << "    Looking for file = " << factorTableName[ p ] << endl ;
 143         cout << "    Found it = " << ((fileOrDir.path().filename() == factorTableName[ p ]) ? "yes" : "no") << endl ;
 144         #endif // DEBUG_PP_FACTOR
 145 
 146         #ifdef USE_EXPERIMENTAL_FILESYSTEM_GNU_CPP
 147         path pathObj( fileOrDir ) ;
 148         if (is_regular_file( pathObj ))
 149         #else
 150         if (fileOrDir.is_regular_file())
 151         #endif
 152         {
 153             if (fileOrDir.path().extension() == factorTableFileExtension)
 154             {
 155                 if (fileOrDir.path().filename() == factorTableName[ p ])
 156                 {
 157                     factorTableFilePath = fileOrDir.path().string() ;
 158                     #ifdef DEBUG_PP_FACTOR
 159                     cout << "    Early out from recursive directory loop" << endl ;
 160                     #endif
 161 
 162                     break ;
 163                 }
 164             }
 165         }
 166     }
 167 
 168     #ifdef DEBUG_PP_FACTOR
 169     if (factorTableFilePath == "")
 170         cout << "Could not find the factor table for p = " << p << endl ;
 171     #endif // DEBUG_PP_FACTOR
 172 
 173     // Open an input stream from the file.
 174     ifstream fin( factorTableFilePath ) ;
 175 
 176     // At this point, we ought to have a table for p as one of the named files above in the same directory
 177     // as the executable.
 178     if (!fin)
 179     {
 180         ostringstream os ;
 181         os << "Missing the factor table for p = " << p << " named " << factorTableName[ p ]
 182         << " Please copy it into the current directory " << PolyParser<PolySymbol, PolyValue>::current_working_dir_ << " which contains the executable!" ;
 183         throw FactorError( os.str(), __FILE__, __LINE__ ) ;
 184     }
 185 
 186     // The header pattern right before the factorizations, e.g.
 187     //    n  #Fac  Factorisation
 188     regex headerPattern                      { R"(^\s*n\s*#Fac\s+Factorisation)" } ;
 189 
 190     // The first factorization line.  It must have at least three numbers, whitespace separated, where third number and subsequent numbers are separated by dots and carets, e.g.
 191     //     84    14  3^2.5.7^2.13.29.43.113.127.337.1429.5419.14449
 192     regex beginFactorizationLinePattern      { R"(\s*\d+\s+\d+\s+(\d+|\^|\.)+)" } ;
 193 
 194     // A continuation line either ends in a backslash, e.g.
 195     //     306    19  3^3.7.19.73.103.307.919.2143.2857.6529.11119.43691.123931.131071.26159806891.27439122228481.755824884241793\
 196     //                47083438319
 197     // Or it ends with a period separating the factors, e.g.,
 198     //     300    28  3^2.5^3.7.11.13.31.41.61.101.151.251.331.601.1201.1321.1801.4051.8101.63901.100801.268501.10567201.13334701.
 199     //                1182468601.1133836730401
 200     regex continuationLinePattern { R"(.*(\\|\.)$)" } ;
 201 
 202     // Scan through the lines of the file.
 203     bool foundHeader           { false } ;
 204     bool continuationLineState { false } ;
 205     vector<string> lineOfTable ;
 206 
 207     for (string line ;  getline( fin, line ) ; )
 208     {
 209         // Skip initial lines in the comment section until we encounter the header line
 210         //   n  #Fac  Factorisation
 211         // Skip this line too, but tag as having found the factorization lines.
 212         if (regex_match( line, headerPattern ))
 213         {
 214             foundHeader = true ;
 215             continue ;
 216         }
 217         else if (!foundHeader)
 218             continue ;
 219 
 220         // Not in a continuation sequence.
 221         if (!continuationLineState)
 222         {
 223             // NOW we are in a continuation pattern.
 224             if (regex_match( line, continuationLinePattern ))
 225                 continuationLineState = true ;
 226 
 227             // Either not a continuation line (i.e. list of factors) or the first continuation line in the sequence.
 228             // Either way, it's the start of the factorization line.
 229             lineOfTable.push_back( line ) ;
 230         }
 231         // In a continuation sequence.
 232         else
 233         {
 234             // Could be a continuation line, or the non-continuation line which terminates the sequence.
 235             // Either way, concatenate to the end of the current line.
 236             lineOfTable.back() += line ;
 237 
 238             // This is the first non-continuation line, so ends the sequence.
 239             if (!regex_match( line, continuationLinePattern ))
 240                 continuationLineState = false ;
 241         }
 242     }
 243 
 244     // Set up the factorization parser.
 245     FactorizationParser< FactorizationSymbol, FactorizationValue<IntType> > parser ;
 246 
 247     // Parse the factorization lines until we see the one which matches p and n.
 248     for (auto & line : lineOfTable)
 249     {
 250         #ifdef DEBUG_PP_FACTOR
 251         // Look for + in the table (incomplete factorization).
 252         // We don't handle incomplete factorizations in the table.  The composite number is likely to be too large to factor
 253         // with Pollard's method.
 254         cout << "Line of factor table for " << p << " ^ n - 1 = " << line << ((line.find( "+" ) != string::npos) ? "WARNING:  incomplete factorization" : "") << endl ;
 255         #endif // DEBUG_PP_FACTOR
 256 
 257         // The factorization is complete.
 258         if (line.find( "+" ) == string::npos)
 259         {
 260             // Parse a factorization line.  For example p = 3 and n = 295 has the line
 261             //  295     9  2.5^2.1181.3221.106185841.70845409351.33083146850190391025301565142735000331370209599...68349655382191991676711\3033493902702...03521
 262             FactorizationValue<IntType> v = parser.parse( line ) ;
 263 
 264             // Did we find an entry for n?
 265             if (FactorizationValue<IntType>::numberStringToInteger( v.numberString_ ) == static_cast<IntType>( n ))
 266             {
 267                 // Copy the factors over.
 268                 factor_.clear() ;
 269                 for (auto & f : v.factor_)
 270                     factor_.push_back( f ) ;
 271 
 272                 // Multiply the factors together, whilst testing each distinct prime factor
 273                 // is probably prime.
 274                 IntType prod = static_cast<IntType>( 1u ) ;
 275 
 276                 for (unsigned int i = 0 ;  i < factor_.size() ;  ++i)
 277                 {
 278                     // Test if the prime factor is really prime.
 279                     if (!isAlmostSurelyPrime( factor_[ i ].prime_ ))
 280                     {
 281                         ostringstream os ;
 282                         os << "Distinct prime factor p = " << factor_[ i ].prime_ << " fails the primality test" ;
 283                         throw FactorError( os.str(), __FILE__, __LINE__ ) ;
 284                     }
 285                     else
 286                         for (int j = 1 ;  j <= factor_[ i ].count_ ;  ++j)
 287                             prod *= factor_[ i ].prime_ ;
 288                 }
 289 
 290                 IntType prod2 = static_cast<IntType>( 1u ) ;
 291 
 292                 //          n
 293                 // Compute p  - 1
 294                 for (int j = 1 ;  j <= n ;  ++j)
 295                     prod2 *= p ;
 296 
 297                 prod2 -= static_cast<IntType>( 1u ) ;
 298 
 299                 // Test whether the product of factors equals the original number again.
 300                 if (prod == prod2)
 301                     return true ;
 302                 else
 303                 {
 304                     ostringstream os ;
 305                     os << "Product of factors doesn't equal the number " << " p^n - 1 " ;
 306                     throw FactorError( os.str(), __FILE__, __LINE__ ) ;
 307                 }
 308             } // Found n
 309         } // no plus sign
 310     } // line of table
 311 
 312     // If we got here we didn't find n in the table.
 313     #ifdef DEBUG_PP_FACTOR
 314     cout << "Table was too small and had no entry for n = " << n << endl ;
 315     #endif // DEBUG_PP_FACTOR
 316 
 317     return false ;
 318 }
 319 
 320 
 321 /*=============================================================================
 322  |
 323  | NAME
 324  |
 325  |    Factorization
 326  |
 327  | DESCRIPTION
 328  |
 329  |    Factor a large integer into primes using tables of prime factors, trial
 330  |    division and Pollard's rho methods.  Factoring algorithm type defaults to 
 331  |    Automatic, which decides when to use each algorithm for best speed.
 332  |
 333  | NOTES
 334  |
 335  |    Tables of prime factors takes negligible time.
 336  |
 337  |                                       ---          1/2
 338  |    Trial division takes max( p    , \/ p  )  = O( N   ) operations, while
 339  |                               t-1       t
 340  |
 341  |                                       ---          1/4
 342  |    Pollard rho takes                \/ p       = O( N  ) operations.
 343  |                                         t-1
 344  | 
 345  +============================================================================*/
 346 
 347 template <typename IntType>
 348 Factorization<IntType>::Factorization( const IntType n, FactoringAlgorithm type, ppuint p, ppuint m )
 349     : n_()
 350     , numFactors_()
 351     , factor_()
 352     , statistics_()
 353     , distinctPrimeFactors_()
 354 {
 355     factor_.clear() ;
 356     distinctPrimeFactors_.clear() ;
 357 
 358     // Initialize the unfactored remainder to n to begin with.
 359     n_ = n ;
 360 
 361     // For unit testing only, specify the type of algorithm.
 362     if (type == FactoringAlgorithm::FactorTable)
 363     {
 364         if (!factorTable( p, m ))
 365         #ifdef DEBUG_PP_FACTOR
 366         cout << "Table lookup factoring failed for p = " << p << " ^ m =" << m << endl ;
 367         #endif
 368         ;
 369     }
 370     else if (type == FactoringAlgorithm::pollardRhoAlgorithm)
 371     {
 372         if (!pollardRho())
 373         #ifdef DEBUG_PP_FACTOR
 374         cout << "Pollard Rho factoring failed for p = " << p << " ^ m =" << m << endl ;
 375         #endif
 376         ;
 377     }
 378     else if (type == FactoringAlgorithm::TrialDivisionAlgorithm)
 379     {
 380         trialDivision() ;
 381         #ifdef DEBUG_PP_FACTOR
 382         cout << "Trial division factoring always succeeds for p = " << p << " ^ m =" << m << endl ;
 383         #endif
 384         ;
 385     }
 386     // type == Automatic
 387     else
 388     {
 389         #ifdef DEBUG_PP_FACTOR
 390         cout << "Try best factoring methods on n = " << n_ << endl ;
 391         #endif
 392 
 393         // Try table lookup first.
 394         if (!factorTable( p, m ))
 395         {
 396             #ifdef DEBUG_PP_FACTOR
 397             cout << "Table lookup factoring failed on n =" << n_ << endl ;
 398             cout << "Try again using Pollard Rho" << endl ;
 399             #endif
 400 
 401             // Try Pollard's method first.  If it fails, keep the
 402             // factors found so far, and retry on the unfactored
 403             // remainder.
 404             if (!pollardRho())
 405             {
 406                 #ifdef DEBUG_PP_FACTOR
 407                 cout << "Pollard Rho failed on n =" << n_ << endl ;
 408                 cout << "Try again with different random seed" << endl ;
 409                 #endif
 410 
 411                 // Try again with random c, but avoid using c in { 0, 1, -2 }
 412                 if (!pollardRho( static_cast<IntType>( 5u )))
 413                 {
 414                     #ifdef DEBUG_PP_FACTOR
 415                     cout << "Pollard Rho failed once again on n =" << n_ << endl ;
 416                     cout << "switching to trial division" << endl ;
 417                     #endif
 418 
 419                     // Still fails?  Fall back onto trial division.  This WILL succeed
 420                     // but can be very slow.
 421                     trialDivision() ;
 422                 }
 423             }
 424         }
 425      }
 426 
 427     #ifdef DEBUG_PP_FACTOR
 428     cout << "     Unsorted prime factors." << endl ;
 429     for (auto & i : factor_)
 430         cout << i.prime_ << " " << i.count_ << endl ;
 431     #endif
 432 
 433     // Sort primes into ascending order.
 434     sort( factor_.begin(), factor_.end(), CompareFactor<IntType>() ) ;
 435 
 436     // Merge duplicated prime factors into unique primes to powers. e.g.
 437     //  2     1           5      0     3     0     6
 438     // 2     2     3     3  =>  2     2     3     3
 439     for (unsigned int i = 1 ;  i < factor_.size() ;  ++i)
 440     {
 441         // This prime factor is a duplicate of the previous factor.
 442         if (factor_[ i ].prime_ == factor_[ i-1 ].prime_)
 443         {
 444             // Merge into current prime power and zero the previous power.
 445             factor_[ i ].count_ += factor_[ i-1 ].count_ ;
 446             factor_[ i-1 ].count_ = 0 ;
 447         }
 448     }
 449 
 450     // Remove factors of 1 (including powers of 0 introduced in the previous step).
 451     typename vector< PrimeFactor<IntType> >::iterator last =
 452         remove_if( factor_.begin(), factor_.end(), Unit<IntType>() ) ;
 453     factor_.erase( last, factor_.end() ) ;
 454     numFactors_ = factor_.size() ;
 455 
 456     // Record a vector of the distinct prime factors.
 457     for (unsigned int i = 0 ;  i < numFactors_ ;  ++i)
 458         distinctPrimeFactors_.push_back( factor_[ i ].prime_ ) ;
 459 
 460     #ifdef DEBUG_PP_FACTOR
 461     numFactors_ = factor_.size() ;
 462     cout << numFactors_ << " sorted unique factors" << endl ;
 463     for (int i = 0 ;  i < numFactors_ ;  ++i)
 464         cout << factor_[ i ].prime_ << " ^ " << factor_[ i ].count_ << endl ;
 465     #endif
 466  }
 467 
 468 
 469 
 470 /*=============================================================================
 471  |
 472  | NAME
 473  |
 474  |    Factorization< IntType >
 475  |
 476  | DESCRIPTION
 477  |
 478  |   Factor a generic integer type n into primes.  Record all the distinct
 479  |   prime factors and how many times each occurs.  Return the number of 
 480  |   primes - 1.
 481  |
 482  |    Factorization<BigInt> f ;
 483  |
 484  |     primes (BigInt *)     List of distinct prime factors.
 485  |     count  (ppuint *)     List of how many times each factor occurred.
 486  | 
 487  |      When n = 1, t = 0, and primes[ 0 ] = count[ 0 ] = 1.
 488  | 
 489  | RETURNS
 490  | 
 491  |      t (int)          Number of prime factors - 1.
 492  |                       Prime factors and their multiplicites are in locations 0 to t.
 493  | 
 494  | EXAMPLE
 495  |                    2
 496  |     For n = 156 = 2  * 3 * 13 we have
 497  | 
 498  |     k   primes[ k ]   count[ k ]
 499  |     ----------------------------
 500  |     0        2            2
 501  |     1        3            1
 502  |     2       13            1
 503  | 
 504  | METHOD
 505  | 
 506  |     Method described by D.E. Knuth, ART OF COMPUTER PROGRAMMING, vol. 2, 3rd ed.,
 507  |                                                              -----
 508  |     in Algorithm A, pgs. 364-365.  The running time is O( max \/ p   , p  )
 509  |                                                                   t-1   t
 510  |     where p  is the largest prime divisor of n and p    is the next largest.
 511  |            t                                        t-1
 512  | 
 513  |     (1)  First divide out all multiples of 2 and 3 and count them.
 514  | 
 515  |     (2)  Next, divide n by all integers d >= 5 except multiples of 2 and 3.
 516  | 
 517  |     (3)  Halt either when all prime factors have been divided out (leaving n = 1)
 518  |          or when the current value of n is prime.  The stopping test
 519  | 
 520  |           (d > | n/d | AND r != 0)
 521  |                --   --
 522  |          detects when n is prime.
 523  | 
 524  |     In the example above, we divided out 2's and 3's first, leaving
 525  |     n = 13.  We continued with trial divisor d = 3.  But now the stopping
 526  |     test was activated because 5 > | 13/5 | = 2, and remainder = 3 (non-zero).
 527  |                                    --    --
 528  |     n = 13 itself is the last prime factor of 156.  We return t = 2.  There
 529  |     are 2 + 1 = 3 distinct primes.
 530  | 
 531  | 
 532  |     If we start with d = 5, and add 2 and 4 in succession, we will run through all
 533  |     the integers except multiples of 2 and 3.  To see this, observe the pattern:
 534  | 
 535  |        Integers:     1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17
 536  |     Mult. of 2:      x     x     x     x     x     x     x     x
 537  |     Mult. of 3:         x        x        x        x        x
 538  |     d:                        x     x           x     x          x
 539  |     Separation:                  2       4         2       4
 540  | 
 541  |     The sequence of divisors d, includes all the primes, so is safe to use for
 542  |     factorization testing.
 543  | 
 544  |     Theorem.  Suppose no divisor d in the sequence divides n.  Suppose at some point,
 545  |     q < d with r != 0.  Then n must be prime.
 546  | 
 547  |     Proof.  We begin with an integer n which has all powers of 2 and 3 removed.
 548  |     Assume, to the contrary, that n is composite:  n = p ... p
 549  |        t                                             1     t
 550  |     n >= p    since p  is the smallest prime factor.
 551  |        1          1
 552  | 
 553  |        2
 554  |     n >= p    since n is composite, so has at least 2 prime factors.
 555  |        1
 556  |              2
 557  |     >= (d + 1)   since p  > d implies p >= (d + 1).  p > d because we assumed
 558  |                           1              1              1
 559  |    that no d in the sequence divided n, so d couldn't be any of the prime divisors p
 560  |                                                                                     i
 561  |            2           2               2
 562  |    = d + 2d + 1 = d + d + d + 1 > d + d
 563  |
 564  |       2
 565  |    > d + n mod d   since n mod d < d
 566  |                                                                               2
 567  |    > | n / d | d + n mod d  because our test said q < d, so | n / d | d < d
 568  |      --     --                                              --     --
 569  | 
 570  |    = n.  So we get the contradiction n > n.  Thus n must be prime.
 571  |                                                 ---
 572  |    Note that this is sharper than the bound d > \/ n
 573  | 
 574  |    Theorem.  Conversely, if n is a prime, no d will divide it, so at some point,
 575  |    d will be large enough that q = | n / d | < d.  r != 0 since n is prime.
 576  |                                  --     --
 577  | 
 578  |    Theorem.  The factoring algorithm works.
 579  | 
 580  |    Proof.  If n == 1 we exit immediately.  If n is composite, we divide out all powers
 581  |    of 2 and 3 (if any).  Since d runs through all possible prime divisors, we divide
 582  |    these out also.  Composite trial d causes no problem;  prime factors of d are divided
 583  |    out of n before we try to divide by d, so such a d cannot divide n.
 584  | 
 585  |     We end in one of two ways (1) All divisors of n have been divided out in which case n = 1.
 586  |    (2) n is a prime, so the stopping test is activiated and n != 1 is recorded as a prime
 587  |    divisor.
 588  | 
 589  | 
 590  | TODO:
 591  | 
 592  |     Can be slow when n is a prime.  We could do a probabilistic test for
 593  |     the primality of n at the statement, "n_is_prime = (r != 0 && q < d) ;"
 594  |     which might speed up the test.
 595  | 
 596  +============================================================================*/
 597 
 598 template <typename IntType>
 599 void Factorization<IntType>::trialDivision()
 600 {
 601     #ifdef DEBUG_PP_FACTOR
 602     cout << "trial division on n =" << n_ << endl ;
 603     #endif
 604 
 605     // Remove factors of 2.
 606     int cnt = 0 ;
 607     while( n_ % static_cast<IntType>( 2u ) == static_cast<IntType>( 0u ))
 608     {
 609         n_ /= static_cast<IntType>( 2u ) ;
 610         ++cnt ;
 611         ++statistics_.num_trial_divides ;
 612     }
 613 
 614     if (cnt != 0)
 615     {
 616         PrimeFactor<IntType> f( static_cast<IntType>( 2u ), cnt ) ;
 617         factor_.push_back( f ) ;
 618     }
 619 
 620     #ifdef DEBUG_PP_FACTOR
 621     cout << "after removing powers of 2, n =" << n_ << endl ;
 622     #endif
 623 
 624     // Remove factors of 3.
 625     cnt = 0 ;
 626     while( n_ % static_cast<IntType>( 3u ) == static_cast<IntType>( 0u ))
 627     {
 628         n_ /= static_cast<IntType>( 3u ) ;
 629         ++cnt ;
 630         ++statistics_.num_trial_divides ;
 631     }
 632 
 633     if (cnt != 0)
 634     {
 635         PrimeFactor<IntType> f( static_cast<IntType>( 3u ), cnt ) ;
 636         factor_.push_back( f ) ;
 637     }
 638 
 639     #ifdef DEBUG_PP_FACTOR
 640     cout << "after removing powers of 3, n =" << n_ << endl ;
 641     #endif
 642 
 643 
 644     //  Factor the rest of n.  Continue until n = 1 (all factors divided out)
 645     //  or until n is prime (so n itself is the last prime factor).
 646 
 647     bool new_d = true ;       //  true if current divisor is different from
 648                               //  the previous one.
 649     bool n_is_prime = false ; //  true if current value of n is prime.
 650     bool flag = true ;        //  Alternately true and false.
 651     IntType d = static_cast<IntType>( 5u ) ;  // First trial divisor.
 652 
 653     do {
 654         //  Integer divide to get quotient and remainder:  n = qd + r
 655         // TODO:  We can speed this up by 2x using a combined divMod call.
 656         IntType q = n_ / d ;
 657         IntType r = n_ % d ;
 658         ++statistics_.num_trial_divides ;
 659 
 660         #ifdef DEBUG_PP_FACTOR
 661         cout << "n = " << n_ << endl ;
 662         cout << "d = " << d  << endl ;
 663         cout << "q = " << q  << endl ;
 664         cout << "r = " << r  << endl ;
 665         #endif
 666 
 667         // Stopping test.
 668         n_is_prime = (r != static_cast<IntType>( 0u ) && q < d) ;
 669 
 670         // d | n.
 671         if (r == static_cast<IntType>( 0u ))
 672         {
 673             n_ = q ;     //  Divide d out of n.
 674 
 675             if (new_d)  //  New prime factor.
 676             {
 677                 PrimeFactor<IntType> f( d, 1u ) ;
 678                 factor_.push_back( f ) ;
 679 
 680                 new_d = false ;
 681                 #ifdef DEBUG_PP_FACTOR
 682                 cout << "after dividing out new prime divisor d = " << d << " n =" << n_ << endl ;
 683                 #endif
 684              }
 685              else
 686              {
 687                  ++factor_[ factor_.size() - 1 ].count_ ; //  Same divisor again.  Increment its count.
 688 
 689                 #ifdef DEBUG_PP_FACTOR
 690                 cout << "after dividing out repeated factor d = " << d << " n =" << n_ << endl ;
 691                 #endif
 692              }
 693         }
 694         else
 695         {
 696             d += (flag ? static_cast<IntType>( 2u ) : static_cast<IntType>( 4u )) ;  // d did not divide n.  Try a new divisor.
 697             flag = !flag ;
 698             new_d = true ;
 699 
 700             #ifdef DEBUG_PP_FACTOR
 701             cout << "next trial divisor d  = " << d << endl ;
 702             #endif
 703         }
 704 
 705     } while( !n_is_prime && (n_ != static_cast<IntType>( 1u )) ) ;
 706 
 707     if (n_ == static_cast<IntType>( 1u ))       //  All factors were divided out.
 708     {
 709         numFactors_ = factor_.size() ;
 710         return ;
 711     }
 712     else
 713     {  //  Current value of n is prime.  It is the last prime factor.
 714         PrimeFactor<IntType> f( n_, 1u ) ;
 715         factor_.push_back( f ) ;
 716         numFactors_ = factor_.size() ;
 717         return ;
 718     }
 719 }
 720 
 721 
 722 
 723 
 724 /*=============================================================================
 725  |
 726  | NAME
 727  |
 728  |    ~Factor
 729  |
 730  | DESCRIPTION
 731  |
 732  |    Destructor for factor.
 733  |
 734  +============================================================================*/
 735 
 736 template <typename IntType>
 737 Factorization<IntType>::~Factorization()
 738 {
 739    // Do nothing.
 740 }
 741 
 742 
 743 /*=============================================================================
 744  |
 745  | NAME
 746  |
 747  |    Factor copy constructor
 748  |
 749  | DESCRIPTION
 750  |
 751  |    Copy a factor object.
 752  |
 753  +============================================================================*/
 754 
 755 template <typename IntType>
 756 Factorization<IntType>::Factorization( const Factorization<IntType> & f )
 757         : n_( f.n_ )
 758         , factor_( f.factor_ )
 759         , distinctPrimeFactors_( f.distinctPrimeFactors_ )
 760         , numFactors_( f.numFactors_ )
 761         , statistics_( f.statistics_ )
 762 {
 763 }
 764 
 765 /*=============================================================================
 766  |
 767  | NAME
 768  |
 769  |    Factor assignment operator
 770  |
 771  | DESCRIPTION
 772  |
 773  |    Assignment constructor.
 774  |
 775  +============================================================================*/
 776 
 777 // Assignment operator.
 778 template <typename IntType>
 779 Factorization<IntType> & Factorization<IntType>::operator=( const Factorization<IntType> & g )
 780 {
 781     // Check for assigning to oneself:  just pass back a reference to the unchanged object.
 782     if (this == &g)
 783         return *this ;
 784 
 785     n_          = g.n_ ;
 786     factor_     = g.factor_ ;
 787     distinctPrimeFactors_ = g.distinctPrimeFactors_ ;
 788     numFactors_ = g.numFactors_ ;
 789     statistics_ = g.statistics_ ;
 790 
 791     // Return a reference to assigned object to make chaining possible.
 792     return *this ;
 793 }
 794 
 795 
 796 /*=============================================================================
 797  |
 798  | NAME
 799  |
 800  |    operator[]
 801  |
 802  | DESCRIPTION
 803  |
 804  |    Return a reference to the ith prime factor and its multiplicity.
 805  |
 806  +============================================================================*/
 807 
 808 template <typename IntType>
 809 PrimeFactor<IntType> & Factorization<IntType>::operator[]( int i )
 810 {
 811     // We throw our own exception here.
 812     if (i > numFactors_)
 813      {
 814         ostringstream os ;
 815         os << "Error accessing Factor at index i = " << i ;
 816         throw FactorRangeError( os.str(), __FILE__, __LINE__ ) ;
 817      }
 818 
 819     return factor_[ i ] ;
 820 }
 821 
 822 
 823 /*=============================================================================
 824  |
 825  | NAME
 826  |
 827  |    numDistinctFactors
 828  |
 829  | DESCRIPTION
 830  |
 831  |    Return the number of distinct prime factors.
 832  |
 833  +============================================================================*/
 834 
 835 template <typename IntType>
 836 ppuint Factorization<IntType>::numDistinctFactors() const
 837 {
 838     return numFactors_ ;
 839 }
 840 
 841  /*=============================================================================
 842   |
 843   | NAME
 844   |
 845   |    primeFactor
 846   |
 847   | DESCRIPTION
 848   |
 849   |    Return the ith distinct prime factor.
 850   |
 851   +============================================================================*/
 852 
 853  template <typename IntType>
 854  IntType Factorization<IntType>::primeFactor( int i ) const
 855  {
 856      // We throw our own exception here.
 857      if (i > numFactors_)
 858      {
 859          ostringstream os ;
 860          os << "Error accessing distinct prime factor at index i = " << i ;
 861          throw FactorRangeError( os.str(), __FILE__, __LINE__ ) ;
 862      }
 863 
 864      return factor_[ i ].prime_ ;
 865  }
 866 
 867 
 868 /*=============================================================================
 869  |
 870  | NAME
 871  |
 872  |    multiplicity
 873  |
 874  | DESCRIPTION
 875  |
 876  |    Return the multiplicity for the ith prime factor.
 877  |
 878  +============================================================================*/
 879 
 880 template <typename IntType>
 881 int Factorization<IntType>::multiplicity( int i ) const
 882 {
 883     // We throw our own exception here.
 884     if (i > numFactors_)
 885      {
 886         ostringstream os ;
 887         os << "Error accessing multiplicity at index i = " << i ;
 888         throw FactorRangeError( os.str(), __FILE__, __LINE__ ) ;
 889      }
 890 
 891     return factor_[ i ].count_ ;
 892 }
 893 
 894 
 895 /*=============================================================================
 896  |
 897  | NAME
 898  |
 899  |    getDistinctPrimeFactors
 900  |
 901  | DESCRIPTION
 902  |
 903  |    Return a vector of only the distinct prime factors.
 904  |
 905  +============================================================================*/
 906 
 907 template <typename IntType>
 908 vector<IntType>Factorization<IntType>::getDistinctPrimeFactors()
 909 {
 910     // TODO:  do I want really copy on the stack?
 911     return distinctPrimeFactors_ ;
 912 }
 913 
 914 
 915 /*=============================================================================
 916  |
 917  | NAME
 918  |
 919  |    isProbablyPrime
 920  |
 921  | DESCRIPTION
 922  |
 923  |      Test whether the number n is probably prime.  If n is composite
 924  |      we are correct always.  If n is prime, we will be wrong no more
 925  |      than about 1/4 of the time on average, probably less for
 926  |      any x and n.
 927  | 
 928  | INPUT
 929  | 
 930  |      n      (int, n >= 0)    Number to test for primality.
 931  |      x      (int, 1 < x < n for n > 6) A random integer.
 932  | 
 933  | RETURNS
 934  | 
 935  |      Primality::Prime          n is definitely prime (table of primes lookup)
 936  |      Primality::Composite      n is definitely composite.
 937  |      Primality::ProbablyPrime  n is probably prime with probability ~3/4.
 938  |
 939  | METHOD
 940  | 
 941  |      Miller-Rabin probabilistic primality test, described by 
 942  |
 943  |          D. E. Knuth, THE ART OF COMPUTER PROGRAMMING, VOL. 2, 3rd ed.,
 944  |          Addison-Wesley, 1981, pgs. 250-265.
 945  |
 946  |          Errata for Volume 2: http://www-cs-faculty.stanford.edu/~knuth/taocp.html
 947  | 
 948  |      Let
 949  |                   k
 950  |          n - 1 = 2  q
 951  | 
 952  |      for odd q.  Now suppose n is prime and
 953  |
 954  |            q
 955  |           x  mod n != 1
 956  | 
 957  |      Then the sequence
 958  |                                          k
 959  |                q          2q               2 q          n-1
 960  |         y = { x  mod n, x    mod n, ..., x     mod n = x    mod n }
 961  | 
 962  |      must end with 1 by Fermat's theorem, and the element in the sequence just before the 1 first appears 
 963  |      must be n-1, since in the field GF( n ), the polynomial equation
 964  | 
 965  |           2
 966  |          y = 1 (mod n)
 967  | 
 968  |      has only two solutions, y = +1 or y = -1.  So if these conditions fail, n is definitely composite.  
 969  |      If the conditions succeed, we cannot tell for sure n is prime, but the probability the algorithm was 
 970  |      fooled in this case is bounded above by about 1/4 for any n.
 971  | 
 972  | EXAMPLE
 973  | 
 974  |      If k = 5, q = 3 then n = 1 + 2^5 3 = 97, which is prime.  
 975  |      
 976  |      For x = 10, we get                    For x = 9, we get
 977  | 
 978  |       q                                   q  
 979  |      x      =   30                       x     =  50
 980  | 
 981  |          1                                   1
 982  |       q 2                                 q 2
 983  |      x      =   27                       x     = 75
 984  | 
 985  |          2                                   2
 986  |       q 2                                 q 2
 987  |      x      =   50                       x     =  96 = n-1 => prob prime
 988  | 
 989  |          3                                   3
 990  |       q 2                                 q 2
 991  |      x      =   75                       x     = 1
 992  | 
 993  |          4                                   4
 994  |       q 2                                 q 2
 995  |      x      =   96 = n-1 => prob prime   x     = 1
 996  | 
 997  |          5                                   5
 998  |       q 2                                 q 2
 999  |      x      =    1                       x     =  1
1000  | 
1001  | 
1002  |      If k = 4, q = 3 then n = 1 + 2^4 3 = 49, which is composite.  
1003  |      
1004  |      For x = 10, we get
1005  | 
1006  |       q
1007  |      x      =   20
1008  | 
1009  |          1
1010  |       q 2
1011  |      x      =   8
1012  | 
1013  |          2
1014  |       q 2
1015  |      x      =   15
1016  | 
1017  |          3       
1018  |       q 2       
1019  |      x      =  29
1020  | 
1021  |          4      
1022  |       q 2      
1023  |      x      =  8  didn't find n-1 => composite.
1024  | 
1025  +============================================================================*/
1026 
1027 template <typename IntType>
1028 Primality isProbablyPrime( const IntType & n, const IntType & x )
1029 {
1030     // Small composite numbers.
1031     if (n == static_cast<IntType>( 0u ) || n == static_cast<IntType>( 1u ) || n == static_cast<IntType>( 4u ))
1032         return Primality::Composite ;
1033 
1034     // Small primes.
1035     if (n == static_cast<IntType>( 2u ) || n == static_cast<IntType>( 3u ) || n == static_cast<IntType>( 5u ))
1036         return Primality::Prime ;
1037 
1038     // Composites aren't prime.
1039     if (n % static_cast<IntType>( 2u ) == static_cast<IntType>( 0u ) ||
1040         n % static_cast<IntType>( 3u ) == static_cast<IntType>( 0u ) ||
1041         n % static_cast<IntType>( 5u ) == static_cast<IntType>( 0u ))
1042         return Primality::Composite ;
1043 
1044     #ifdef DEBUG_PP_PRIMALITY_TESTING
1045     cout << "Begin the Miller-Rabin primality test on n = " << n << " for random x = " << x << endl ;
1046     #endif
1047 
1048     //                                                     k
1049     // Factor out powers of 2 to find odd q where n - 1 = 2  q.
1050     IntType reduced_n = n - static_cast<IntType>( 1u ) ;
1051     IntType k   = static_cast<IntType>( 0u ) ;
1052 
1053     while( reduced_n % static_cast<IntType>( 2u ) == static_cast<IntType>( 0u ))
1054     {
1055         reduced_n /= static_cast<IntType>( 2u ) ;
1056         ++k ;
1057     }
1058 
1059     IntType q = reduced_n ;
1060 
1061     #ifdef DEBUG_PP_PRIMALITY_TESTING
1062     cout << "After factoring out powers of 2, q = " << q << endl ;
1063     cout << " n = 1 + 2 ^ k q " << endl ;
1064     cout << n << " = 1 + 2 ^ " << k << " " << q << endl ;
1065     #endif
1066 
1067     //       q
1068     //  y = x  (mod n)
1069     PowerMod<IntType> power_mod( n ) ;
1070     IntType y = power_mod( x, q ) ;
1071 
1072     #ifdef DEBUG_PP_PRIMALITY_TESTING
1073     cout << "Compute y = x ^ q (mod n) " << endl ;
1074     cout << y << " = " << x << " ^ " << q << " (mod " << n << ")" << endl ;
1075     cout << "Iterate j = 0 ... " << k - static_cast<IntType>( 1u ) << endl ;
1076     #endif
1077 
1078     for (IntType j = static_cast<IntType>( 0u ) ;  j < k ;  ++j)
1079     {
1080         #ifdef DEBUG_PP_PRIMALITY_TESTING
1081         cout << "j = " << j << endl ;
1082         #endif
1083 
1084         //           q
1085         // Check if x  = 1 mod n immediately.
1086         if ( j == static_cast<IntType>( 0u ) && y == static_cast<IntType>( 1u ) )
1087         {
1088             #ifdef DEBUG_PP_PRIMALITY_TESTING
1089             cout << "y = 1 immediately:  probably prime!  j = " << j << " y = " << y << endl ;
1090             #endif
1091 
1092             return Primality::ProbablyPrime ;
1093 
1094         }
1095         //           q   2q
1096         // Check if x , x  ,... = n-1 mod n at any point.
1097         else if (y == n - static_cast<IntType>( 1u ))
1098         {
1099             #ifdef DEBUG_PP_PRIMALITY_TESTING
1100             cout << "y = n-1:  probably prime!  j = " << j << " y = " << y << endl ;
1101             #endif
1102 
1103             return Primality::ProbablyPrime ;
1104         }
1105 
1106         // Found a 1 but never found an n-1 term before it:  n can't be prime.
1107         else if (j > static_cast<IntType>( 0u ) && y == static_cast<IntType>( 1u ))
1108         {
1109             #ifdef DEBUG_PP_PRIMALITY_TESTING
1110             cout << "Found 1 but not n-1 term:  composite!  j = " << j << " y = " << y << endl ;
1111             #endif
1112 
1113             return Primality::Composite ;
1114         }
1115 
1116         //          2
1117         // Compute y  (mod n)
1118         y = power_mod( y, static_cast<IntType>( 2u ) ) ;
1119 
1120         #ifdef DEBUG_PP_PRIMALITY_TESTING
1121         cout << "Looping again:  y^2 (mod n) = " <<  y << endl ;
1122         #endif
1123     }
1124 
1125     #ifdef DEBUG_PP_PRIMALITY_TESTING
1126     cout << "After looping finished, we saw no 1 or (n-1) terms:  composite!" << y << endl ;
1127     #endif
1128 
1129     // If we got here j = k-1 and the sequence had no 1 or n-1 terms so n is composite.
1130     return Primality::Composite ;
1131 }
1132 
1133 
1134 
1135 
1136 /*=============================================================================
1137  |
1138  | NAME
1139  |
1140  |    isAlmostSurelyPrime
1141  |
1142  | DESCRIPTION
1143  |
1144  |      Test whether the number n >= 0 is almost surely prime.  If n is
1145  |      composite, we always return false.  If n is prime, the
1146  |      probability of returning true successfully is
1147  | 
1148  |                  NUM_PRIME_TEST_TRIALS
1149  |         1 - (1/4)
1150  | 
1151  | METHOD
1152  | 
1153  |     Use the Miller-Rabin probabilistic primality test for lots of
1154  |     random integers x.  If the test shows n is composite for any
1155  |     given x, is is non-prime for sure.
1156  | 
1157  |     But if n is prime, P( failure | n=prime ) <= 1/4, on each
1158  |     trial, so if the random numbers are independent,
1159  | 
1160  |                                    NUM_PRIME_TEST_TRIALS 
1161  |     P( failure | n=prime ) <= (1/4)
1162  | 
1163  |     for 25 trials, P <= 0.8817841970012523e-16
1164  |     and for 14 trials, P <= 3.7252902984619141e-09
1165  |
1166  |     The algorithm is from
1167  |
1168  |        D. E. Knuth, THE ART OF COMPUTER PROGRAMMING, VOL. 2, 3rd ed.,
1169  |        Addison-Wesley, 1981, pgs. 250-265.
1170  |
1171  |        Errata for Volume 2: http://www-cs-faculty.stanford.edu/~knuth/taocp.html
1172  |
1173  +============================================================================*/
1174 
1175 template <typename IntType>
1176 bool isAlmostSurelyPrime( const IntType & n )
1177 {
1178     IntType
1179         trial{ 0u },
1180         x{ 3u } ;
1181 
1182     constexpr ppuint NUM_PRIME_TEST_TRIALS { 14u }  ;
1183 
1184     // Generate uniform random numbers in the range [0, n).
1185     UniformRandomIntegers< IntType > randum( n ) ;
1186 
1187     for (trial = static_cast<IntType>( 1u ) ;  trial <= static_cast<IntType>( NUM_PRIME_TEST_TRIALS ) ;  ++trial)
1188     {
1189         //  Generate a new random integer such that 1 < x < n.
1190         x = randum.rand() ;
1191 
1192         // Random number has to be > 1
1193         if (x <= static_cast<IntType>( 1u ))
1194             x = static_cast<IntType>( 3u ) ;
1195 
1196         #ifdef DEBUG_PP_PRIMALITY_TESTING
1197         cout << "isAlmostSurelyPrime( n = " << n << " ) at trial = " << trial << " with random integer x = " << x << endl ;
1198         #endif
1199 
1200         if (isProbablyPrime( n, x ) == Primality::Prime)
1201             return true ;
1202         else if (isProbablyPrime( n, x ) == Primality::Composite)
1203             return false ;
1204         // else it's probably prime.
1205     }
1206 
1207     // Almost surely prime because it passed the probably prime test
1208     // above enough times in the loop.
1209     return true ;
1210 }
1211 
1212 
1213 
1214 /*=============================================================================
1215  | 
1216  | NAME
1217  | 
1218  |     skipTest
1219  | 
1220  | DESCRIPTION
1221  | 
1222  |      Make the test p  | (p - 1).
1223  |                     i
1224  | 
1225  | INPUT
1226  | 
1227  |     i      (int)       ith prime divisor of r.
1228  | 
1229  |     primes (bigint *)  Array of distict prime factors of r.  primes[ i ] = p
1230  |                                                                             i
1231  |     p      (int)        p >= 2.
1232  | 
1233  | RETURNS
1234  | 
1235  |     true    if the test succeeds for some i.
1236  |     false   otherwise
1237  | 
1238  | EXAMPLE
1239  | 
1240  |     Suppose i = 0, primes[ 0 ] = 2 and p = 5.  Return true since 2 | 5-1.
1241  | 
1242  | METHOD
1243  | 
1244  |     Test if (p - 1) mod p = 0 for all i.
1245  |                          i
1246  +============================================================================*/
1247 
1248 template <typename IntType>
1249 bool Factorization<IntType>::skipTest( ppuint p, int i ) const
1250 {
1251 
1252 // p  cannot divide the smaller number (p-1)
1253 //  i
1254 if ( static_cast<IntType>( p-1 ) < static_cast<IntType>( factor_[ i ].prime_ ) )
1255     return false ;
1256 else
1257     return(
1258         (static_cast<IntType>(p - 1) % static_cast<IntType>( factor_[ i ].prime_ ) ==
1259          static_cast<IntType>( 0u ))
1260         ? true : false ) ;
1261 
1262 } // ====================== end of function skipTest =======================
1263 
1264 
1265 
1266 /*=============================================================================
1267  |
1268  | NAME
1269  |
1270  |    pollardRho
1271  |
1272  | DESCRIPTION
1273  |
1274  |    Factor an integer using Pollard's rho method as modified by Brent and
1275  |    detailed in
1276  |
1277  |        D. E. Knuth, THE ART OF COMPUTER PROGRAMMING, VOL. 2, 3rd ed.,
1278  |        Addison-Wesley, 1981, pgs. 250-265.
1279  |
1280  |        Errata for Volume 2:
1281  |        http://www-cs-faculty.stanford.edu/~knuth/taocp.html
1282  | 
1283  +============================================================================*/
1284 
1285 template <typename IntType>
1286 bool Factorization<IntType>::pollardRho( const IntType & c )
1287 {
1288     IntType x  = static_cast<IntType>( 5u ) ;
1289     IntType xp = static_cast<IntType>( 2u ) ;
1290     IntType k  = static_cast<IntType>( 1u ) ;
1291     IntType l  = static_cast<IntType>( 1u ) ;
1292     IntType g  = static_cast<IntType>( 0u ) ;
1293 
1294     bool status = true ;
1295 
1296     #ifdef DEBUG_PP_FACTOR
1297     cout << "Pollard rho n = " << n_ << " c = " << c << endl ;
1298     #endif
1299 
1300     // Seed 1 as a factor and early out if n == 1.
1301     PrimeFactor<IntType>f( static_cast<IntType>( 1u ), 1u ) ;
1302     factor_.push_back( f ) ;
1303     numFactors_ = 1u ;
1304 
1305     if (n_ == static_cast<IntType>( 1u ))
1306         return status ;
1307 
1308     while( true )
1309     {
1310         if (isAlmostSurelyPrime( n_ ))
1311         {
1312             #ifdef DEBUG_PP_FACTOR
1313             cout << "    >>>>>prime factor n_ = " << n_ << endl ;
1314             #endif
1315 
1316             PrimeFactor<IntType>f( n_, 1u ) ;
1317             factor_.push_back( f ) ;
1318             ++numFactors_ ;
1319 
1320             status = true ;
1321             ++statistics_.num_primality_tests ;
1322             goto Exit ;
1323         }
1324 
1325         while (true)
1326         {
1327             // TODO:   We can speed up by not checking gcd when k > l/2
1328             // if (k > l/2u)
1329             //     continue ;
1330             // TODO:  We can speed up by accumulating gcd products.
1331             IntType absDiff = (xp > x) ? (xp - x) : (x - xp) ;
1332             g = gcd( absDiff, n_ ) ;
1333             ++statistics_.num_gcds ;
1334 
1335             #ifdef DEBUG_PP_FACTOR
1336             cout << "    inner while loop, gcd = g = " << g << " |x -xp| = " << absDiff << " n_ = " << n_
1337                  << " k = " << k << " l = " << l << endl ;
1338             #endif
1339 
1340             if  (g == static_cast<IntType>( 1u ))
1341             {
1342                 --k ;
1343                 if (k == static_cast<IntType>( 0u ))
1344                 {
1345                     xp = x ;
1346                     l = l * static_cast<IntType>( 2u ) ;
1347                     k = l ;
1348                 }
1349                 x = (x * x + c) % n_ ;
1350                 ++statistics_.num_squarings ;
1351             }
1352             else if (g == n_)
1353             {
1354                 #ifdef DEBUG_PP_FACTOR
1355                 cout << "    failure:  g equals non-prime n_ = " << n_ << endl ;
1356                 #endif
1357 
1358                 status = false ;
1359                 goto Exit ;
1360             }
1361             else if (isAlmostSurelyPrime( g ))
1362             {
1363                 #ifdef DEBUG_PP_FACTOR
1364                 cout << "    >>>>>prime factor g = " << g << endl ;
1365                 #endif
1366 
1367                 PrimeFactor<IntType> f( g, 1u ) ;
1368                 factor_.push_back( f ) ;
1369                 ++numFactors_ ;
1370 
1371                 ++statistics_.num_primality_tests ;
1372             }
1373             else
1374             {
1375                 #ifdef DEBUG_PP_FACTOR
1376                 cout << "    g = " << g << " isn't prime " << endl ;
1377                 #endif
1378 
1379                 status = false ;
1380                 goto Exit ;
1381             }
1382 
1383             n_ /= g ;
1384             x = x % n_ ;
1385             xp = xp % n_ ;
1386 
1387             #ifdef DEBUG_PP_FACTOR
1388             cout << "    after division, n_ = " << n_ << " x = " << x << " xp = " << xp << endl ;
1389             #endif
1390 
1391             break ;
1392         } // end gcd loop
1393     }
1394 
1395 
1396 Exit:
1397     return status ;
1398 }
1399 
1400 /*==============================================================================
1401 |                     Forced Template Instantiations                           |
1402 ==============================================================================*/
1403 
1404 // C++ doesn't automatically generate templated classes or functions until they
1405 // are first used.  So depending on the order of compilation, the linker may say
1406 // the templated functions are undefined.
1407 //
1408 // Therefore, explicitly instantiate ALL templates here:
1409 
1410 template Factorization<ppuint>::Factorization( const ppuint, FactoringAlgorithm, ppuint p, ppuint m ) ;
1411 template Factorization<BigInt>::Factorization( const BigInt, FactoringAlgorithm, ppuint p, ppuint m ) ;
1412 
1413 template Factorization<ppuint>::~Factorization() ;
1414 template Factorization<BigInt>::~Factorization() ;
1415 
1416 template Factorization<ppuint>::Factorization( const Factorization<ppuint> & ) ;
1417 template Factorization<BigInt>::Factorization( const Factorization<BigInt> & ) ;
1418 
1419 template Factorization<ppuint> & Factorization<ppuint>::operator=( const Factorization<ppuint> & ) ;
1420 template Factorization<BigInt> & Factorization<BigInt>::operator=( const Factorization<BigInt> & ) ;
1421 
1422 // TODO factorTable only uses BigInt type:
1423 template bool Factorization<ppuint>::factorTable( ppuint, ppuint ) ;
1424 template bool Factorization<BigInt>::factorTable( ppuint, ppuint ) ;
1425 
1426 template void Factorization<ppuint>::trialDivision() ;
1427 template void Factorization<BigInt>::trialDivision() ;
1428 
1429 template bool Factorization<ppuint>::pollardRho( const ppuint & ) ;
1430 template bool Factorization<BigInt>::pollardRho( const BigInt & ) ;
1431 
1432 template ppuint Factorization<ppuint>::numDistinctFactors() const ;
1433 template ppuint Factorization<BigInt>::numDistinctFactors() const ;
1434 
1435 template PrimeFactor<ppuint> & Factorization<ppuint>::operator[]( int ) ;
1436 template PrimeFactor<BigInt> & Factorization<BigInt>::operator[]( int ) ;
1437 
1438 template ppuint Factorization<ppuint>::primeFactor( int ) const ;
1439 template BigInt Factorization<BigInt>::primeFactor( int ) const ;
1440 
1441 template int Factorization<ppuint>::multiplicity( int ) const  ;
1442 template int Factorization<BigInt>::multiplicity( int ) const  ;
1443 
1444 template Primality isProbablyPrime( const ppuint &, const ppuint & ) ;
1445 template Primality isProbablyPrime( const BigInt &, const BigInt & ) ;
1446 
1447 template bool isAlmostSurelyPrime( const ppuint & ) ;
1448 template bool isAlmostSurelyPrime( const BigInt & ) ;
1449 
1450 template bool Factorization<int>::skipTest( ppuint, int ) const ;
1451 template bool Factorization<BigInt>::skipTest( ppuint, int ) const ;
1452 
1453 template vector<ppuint>Factorization<ppuint>::getDistinctPrimeFactors() ;
1454 template vector<BigInt>Factorization<BigInt>::getDistinctPrimeFactors()  ;