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.3 - A Program for Computing Primitive Polynomials.
17| Copyright (C) 1999-2024 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
55using 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.hpp" // Global functions.
70#include "ppArith.hpp" // Basic arithmetic functions.
71#include "ppBigInt.hpp" // Arbitrary precision integer arithmetic.
72#include "ppOperationCount.hpp" // OperationCount collection for factoring and poly finding.
73#include "ppFactor.hpp" // Prime factorization and Euler Phi.
74#include "ppPolynomial.hpp" // Polynomial operations and mod polynomial operations.
75#include "ppParser.hpp" // Parsing of polynomials and I/O services.
76#include "ppUnitTest.hpp" // 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
96template <typename IntType>
97bool 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
347template <typename IntType>
348Factorization<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
598template <typename IntType>
599void 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
736template <typename IntType>
737Factorization<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
755template <typename IntType>
756Factorization<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.
778template <typename IntType>
779Factorization<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
808template <typename IntType>
809PrimeFactor<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
835template <typename IntType>
836ppuint 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
880template <typename IntType>
881int 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
907template <typename IntType>
908vector<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
1027template <typename IntType>
1028Primality 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
1175template <typename IntType>
1176bool 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
1248template <typename IntType>
1249bool Factorization<IntType>::skipTest( ppuint p, int i ) const
1250{
1251
1252// p cannot divide the smaller number (p-1)
1253// i
1254if ( static_cast<IntType>( p-1 ) < static_cast<IntType>( factor_[ i ].prime_ ) )
1255 return false ;
1256else
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
1285template <typename IntType>
1286bool 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
1396Exit:
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
1410template Factorization<ppuint>::Factorization( const ppuint, FactoringAlgorithm, ppuint, ppuint ) ;
1411template Factorization<BigInt>::Factorization( const BigInt, FactoringAlgorithm, ppuint, ppuint ) ;
1412
1413template Factorization<ppuint>::~Factorization() ;
1414template Factorization<BigInt>::~Factorization() ;
1415
1416template Factorization<ppuint>::Factorization( const Factorization<ppuint> & ) ;
1417template Factorization<BigInt>::Factorization( const Factorization<BigInt> & ) ;
1418
1419template Factorization<ppuint> & Factorization<ppuint>::operator=( const Factorization<ppuint> & ) ;
1420template Factorization<BigInt> & Factorization<BigInt>::operator=( const Factorization<BigInt> & ) ;
1421
1422// TODO factorTable only uses BigInt type:
1423template bool Factorization<ppuint>::factorTable( ppuint, ppuint ) ;
1424template bool Factorization<BigInt>::factorTable( ppuint, ppuint ) ;
1425
1426template void Factorization<ppuint>::trialDivision() ;
1427template void Factorization<BigInt>::trialDivision() ;
1428
1429template bool Factorization<ppuint>::pollardRho( const ppuint & ) ;
1430template bool Factorization<BigInt>::pollardRho( const BigInt & ) ;
1431
1432template ppuint Factorization<ppuint>::numDistinctFactors() const ;
1433template ppuint Factorization<BigInt>::numDistinctFactors() const ;
1434
1435template PrimeFactor<ppuint> & Factorization<ppuint>::operator[]( int ) ;
1436template PrimeFactor<BigInt> & Factorization<BigInt>::operator[]( int ) ;
1437
1438template ppuint Factorization<ppuint>::primeFactor( int ) const ;
1439template BigInt Factorization<BigInt>::primeFactor( int ) const ;
1440
1441template int Factorization<ppuint>::multiplicity( int ) const ;
1442template int Factorization<BigInt>::multiplicity( int ) const ;
1443
1444template Primality isProbablyPrime( const ppuint &, const ppuint & ) ;
1445template Primality isProbablyPrime( const BigInt &, const BigInt & ) ;
1446
1447template bool isAlmostSurelyPrime( const ppuint & ) ;
1448template bool isAlmostSurelyPrime( const BigInt & ) ;
1449
1450template bool Factorization<int>::skipTest( ppuint, int ) const ;
1451template bool Factorization<BigInt>::skipTest( ppuint, int ) const ;
1452
1453template vector<ppuint>Factorization<ppuint>::getDistinctPrimeFactors() ;
1454template vector<BigInt>Factorization<BigInt>::getDistinctPrimeFactors() ;