1/*==============================================================================
2|
3| NAME
4|
5| ppFactor.hpp
6|
7| DESCRIPTION
8|
9| Header for integer factoring classes.
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 the !DOT! replaced by . and the !AT! replaced by @
34|
35==============================================================================*/
36
37// Wrap this header file to prevent duplication if it is included
38// accidentally more than once.
39#ifndef __PPFACTOR_H__
40#define __PPFACTOR_H__
41
42
43/*=============================================================================
44|
45| NAME
46|
47| FactorError General factor error, including internal memory errors.
48| FactorRangeError Input range error.
49|
50| DESCRIPTION
51|
52| Exception classes for the Factor class
53| derived from the STL exception class runtime_error.
54|
55+============================================================================*/
56
57class FactorError : public runtime_error
58{
59 public:
60 // Throw with error message, file name and line number.
61 FactorError( const string & description, const string & file, const int & line )
62 : runtime_error( description + " in file " + file + " at line " + to_string(line) )
63 {
64 } ;
65
66 // Throw with an error message.
67 FactorError( const string & description )
68 : runtime_error( description )
69 {
70 } ;
71
72 // Default throw with no error message.
73 FactorError()
74 : runtime_error( "Factor error: " )
75 {
76 } ;
77
78} ; // end class BigIntMathError
79
80
81class FactorRangeError : public FactorError
82{
83 public:
84 // Throw with error message, file name and line number.
85 FactorRangeError( const string & description, const string & file, const int & line )
86 : FactorError( description + " in file " + file + " at line " + to_string(line) )
87 {
88 };
89
90 // Throw with an error message.
91 FactorRangeError( const string & description )
92 : FactorError( description )
93 {
94 } ;
95
96 // Default throw with no error message.
97 FactorRangeError()
98 : FactorError( "Factor range error: " )
99 {
100 } ;
101
102} ; // end class FactorRangeError
103
104
105/*=============================================================================
106 |
107 | NAME
108 |
109 | PrimeFactor
110 |
111 | DESCRIPTION
112 |
113 | Class to neatly package up unique prime factors to powers.
114 |
115 +============================================================================*/
116
117template <typename IntType>
118class PrimeFactor
119{
120 public:
121 explicit PrimeFactor( IntType prime = static_cast<IntType>( 0u ), int count = 0 )
122 : prime_( prime )
123 , count_( count )
124 {
125 } ;
126
127 ~PrimeFactor()
128 {
129 } ;
130
131 PrimeFactor( const PrimeFactor & factor )
132 : prime_( factor.prime_ )
133 , count_( factor.count_ )
134 {
135 }
136
137 PrimeFactor & operator=( const PrimeFactor & factor )
138 {
139 // Check for assigning to oneself: just pass back a reference to the unchanged object.
140 if (this == &factor)
141 return *this ;
142
143 prime_ = factor.prime_ ;
144 count_ = factor.count_ ;
145
146 return *this ;
147 } ;
148
149 // Print function for a factor.
150 friend ostream & operator<<( ostream & out, const PrimeFactor & factor )
151 {
152 out << factor.prime_ << " ^ " << factor.count_ << " " ;
153
154 return out ;
155 }
156
157 // Allow direct access to this simple data type for convenience.
158 public:
159 IntType prime_ ;
160 int count_ ;
161} ;
162
163
164/*=============================================================================
165 |
166 | NAME
167 |
168 | CompareFactor
169 |
170 | DESCRIPTION
171 |
172 | Class to sort unique prime factors to powers into order by prime size.
173 |
174 +============================================================================*/
175
176template <typename IntType>
177class CompareFactor
178{
179 public:
180 bool operator()( const PrimeFactor<IntType> & s1, const PrimeFactor<IntType> & s2 )
181 {
182 return s1.prime_ < s2.prime_ ;
183 }
184} ;
185
186
187/*=============================================================================
188 |
189 | NAME
190 |
191 | Unit
192 |
193 | DESCRIPTION
194 | 0 e
195 | Class to check for unit factors of the form p or 1 = 1
196 |
197 +============================================================================*/
198
199template <typename IntType>
200class Unit
201{
202 public:
203 bool operator()( const PrimeFactor<IntType> & s )
204 {
205 return s.count_ == 0u || s.prime_ == static_cast<IntType>( 1u ) ;
206 }
207} ;
208
209
210/*=============================================================================
211|
212| NAME
213|
214| Factorization
215|
216| DESCRIPTION
217|
218| Class for single and multiprecision factoring.
219|
220| EXAMPLE
221|
222| ppuint p{ 2 } ;
223| ppuint n{ 4 } ;
224| BigInt p_to_n_minus_1 { power( p, 4 ) - 1 } ;
225|
226| n n
227| This will give us the factors of p - 1. We can pass in either p and n or (p - 1), but p and n might be
228| the better choice for really large numbers because then we can look up the factorization in a table instead
229| of computing it.
230|
231| Factorization<BigInt> factors_of_p_to_n_minus_1( p_to_n_minus_1, FactoringAlgorithm::Automatic, p, n ) ;
232|
233| And so will this for a smaller sized integer.
234|
235| Factorization<ppuint> fact( power( p, 4 ) - 1 ) ;
236|
237| int numDistinctFactors = fact.num() ;
238| for (int i = 0 ; i < numDistinctFactors ; ++i)
239| {
240| BigInt primeFactor = fact[ i ] ;
241| int multiplicity = fact.multiplicity( i ) ;
242| }
243|
244| NOTES
245|
246| The member functions and friends are documented in detail ppBigInt.cpp
247|
248+============================================================================*/
249
250// Different flavors of factoring algorithm.
251enum class FactoringAlgorithm
252{
253 Automatic,
254 TrialDivisionAlgorithm,
255 pollardRhoAlgorithm,
256 FactorTable
257} ;
258
259
260// We need both single precision and multi-precision factoring, so use
261// a template with an integer type.
262template <typename IntType>
263class Factorization
264{
265 public:
266 // Factor n into distinct primes. Default constructor factors nothing.
267 Factorization( const IntType num = static_cast<IntType>( 1u ),
268 FactoringAlgorithm type = FactoringAlgorithm::Automatic, ppuint p = 0, ppuint m = 0 ) ;
269
270 ~Factorization() ;
271
272 // Copy constructor.
273 Factorization( const Factorization<IntType> & f ) ;
274
275 // Assignment operator.
276 Factorization<IntType> & operator=( const Factorization<IntType> & g ) ;
277
278 // Return the number of distinct factors.
279 ppuint numDistinctFactors() const ;
280
281 // Return the ith prime factor along with its multiplicity as an lvalue so we can change either.
282 PrimeFactor<IntType> & operator[]( int i ) ;
283
284 // Return the ith prime factor.
285 IntType primeFactor( int i ) const ;
286
287 // Return the multiplicity of the ith prime factor.
288 int multiplicity( int i ) const ;
289
290 // True if p | (p - 1).
291 // i
292 bool skipTest( ppuint p, int i ) const ;
293
294 // Factoring with tables.
295 bool factorTable( ppuint p, ppuint n ) ;
296
297 // ---
298 // Factoring by trial division up to \/ n
299 void trialDivision() ;
300
301 // Fast probabilistic factoring method.
302 bool pollardRho( const IntType & c = static_cast<IntType>( 2u )) ;
303
304 // Return a vector of only the distinct prime factors.
305 vector<IntType> getDistinctPrimeFactors() ;
306
307 // Allow direct access to this simple data type for convenience.
308 public:
309 OperationCount statistics_ ;
310
311 private:
312 // The unfactored remainder.
313 IntType n_ ;
314
315 // Total number of distinct factors.
316 ppuint numFactors_ ;
317
318 // Array of distinct prime factors of n with their multiplicities.
319 vector< PrimeFactor<IntType> > factor_ ;
320
321 // Distinct prime factors.
322 vector<IntType> distinctPrimeFactors_ ;
323} ;
324
325
326/*=============================================================================
327|
328| NAME
329|
330| Primality
331|
332| DESCRIPTION
333|
334| What confidence a number is prime?
335|
336+============================================================================*/
337
338enum class Primality
339{
340 Prime = 0,
341 Composite,
342 ProbablyPrime,
343 Undefined
344
345} ;
346
347
348/*=============================================================================
349 |
350 | NAME
351 |
352 | isProbablyPrime
353 |
354 | DESCRIPTION
355 |
356 | True if n is likely to be prime. Tests on a random integer x.
357 |
358 +============================================================================*/
359
360template <typename IntType>
361Primality isProbablyPrime( const IntType & n, const IntType & x ) ;
362
363
364/*=============================================================================
365 |
366 | NAME
367 |
368 | UniformRandomNumber
369 |
370 | DESCRIPTION
371 |
372 | Return a uniform random integer in [0, range).
373 |
374 | NOTES
375 |
376 | Use the JKISS random number generator from the article,
377 | Good Practice in (Pseudo) Random Number Generation for Bioinformatics Applications
378 | by David Jones, UCL Bioinformatics Group (E-mail: d.jones@cs.ucl.ac.uk)
379 | (Last revised May 7th 2010)
380 |
381 +============================================================================*/
382
383template <typename IntType>
384class UniformRandomIntegers
385{
386 public:
387 // Seed the JKISS generator parameters x, y, z, c, and set up the range of the generator [0, range)
388 explicit UniformRandomIntegers( const IntType range )
389 : range( range )
390 {
391 ++num_of_initializations ;
392
393 // Reseed at the first call to this class and only occasionally thereafter because it's slow.
394 if (num_of_initializations % how_often_to_reseed == 0)
395 {
396 x = trueRandomNumberFromDevice() ;
397
398 // Seed y must not be zero!
399 while ( !(y = trueRandomNumberFromDevice()) ) ;
400
401 z = trueRandomNumberFromDevice() ;
402
403 // We don’t really need to set c as well but let's anyway...
404 // NOTE: offset c by 1 to avoid z = c = 0
405 c = trueRandomNumberFromDevice() % JKISS_MAGIC_NUM + 1 ; // Should be less than 698769069
406 }
407 }
408
409 // Destructor has nothing to do.
410 UniformRandomIntegers()
411 {
412 } ;
413
414 IntType rand()
415 {
416 ppuint32 kiss = jKiss() ;
417
418 // [0, range) falls within our generator's range [0, JKISSMAX).
419 // To preserve a uniform distribution, map numbers in [0, JKISSMAX) down to [0, range) in the following way:
420 // Map [0, range), => [0, range)
421 // Map [range+1, 2 range) => [0, range)
422 // ...
423 // Map [(n-1) range + 1, n range) => [0, range)
424 // Discard numbers in [n range + 1, JKISSMAX) because mapping numbers in this interval < range would inject a non-uniform bias.
425 if (range < IntType( JKISSMAX ))
426 {
427 // It's safe to down convert the integer type.
428 ppuint32 range2{ static_cast<ppuint32>( range ) } ;
429
430 // Discard the random number unless it falls in a multiple of the range.
431 // Retry with a new random number. I hope we don't get an infinite loop!
432 ppuint32 within_multiple_of_range{ JKISSMAX - (JKISSMAX % static_cast<ppuint32>(range)) } ;
433 while (kiss > within_multiple_of_range)
434 kiss = jKiss() ;
435 kiss = jKiss() % range2 ;
436 }
437 // If the range is larger than the maximum of our generator, do nothing. We don't want to scale up
438 // and have a non-uniform distribution.
439 return static_cast<IntType>( kiss ) ;
440 } ;
441
442 private:
443 static int num_of_initializations ;
444 static const int how_often_to_reseed = 10000 ;
445 static const int JKISS_MAGIC_NUM = 698769068 ;
446
447 // Default seeds for JKISS.
448 ppuint32 x{ 123456789 } ;
449 ppuint32 y{ 987654321 } ;
450 ppuint32 z{ 43219876 } ;
451 ppuint32 c{ 6543217 } ;
452
453 // Top end of the jKiss() generator range.
454 const ppuint32 JKISSMAX { numeric_limits<ppuint32>::max() } ;
455
456 // Return uniform random numbers in the range [0, range).
457 IntType range ;
458
459 // Read from the pseudo device /dev/urandom which return true random integers.
460 unsigned int trueRandomNumberFromDevice()
461 {
462 // The device /dev/urandom only returns bytes, so grab four at a time to make on 32-bit unsigned integer.
463 union RandomDeviceByteToInteger
464 {
465 ppuint32 integer32 ;
466 char bytes4[ 4 ] ;
467 } ;
468
469 RandomDeviceByteToInteger numFromRandomDev ;
470
471 ifstream fin( "/dev/urandom", ifstream::binary ) ;
472 if (fin.good())
473 {
474 fin.read( numFromRandomDev.bytes4, 4 ) ;
475
476 if (!fin)
477 {
478 ostringstream os ;
479 os << "Cannot read from random device /dev/urandom" ;
480 throw FactorError( os.str(), __FILE__, __LINE__ ) ;
481 }
482 }
483 else
484 {
485 ostringstream os ;
486 os << "Cannot open random device /dev/urandom" ;
487 throw FactorError( os.str(), __FILE__, __LINE__ ) ;
488 }
489
490 return numFromRandomDev.integer32 ;
491 }
492
493 // 32-bit integer random number generator.
494 ppuint32 jKiss()
495 {
496 x = 314527869 * x + 1234567 ;
497
498 y ^= (y << 5) ;
499 y ^= (y >> 7) ;
500 y ^= (y << 22) ;
501
502 ppuint t{ 4294584393ULL * z + c } ; // This must be a 64-bit integer type.
503 c = t >> 32 ;
504 z = static_cast<ppuint32>( t ) ;
505
506 return x + y + z ;
507 }
508} ;
509
510
511/*==============================================================================
512| Static Class Variables Initialization |
513==============================================================================*/
514
515// Static class variables must be initialized outside the class.
516
517template <typename IntType> int UniformRandomIntegers<IntType>::num_of_initializations = 0 ;
518
519
520/*=============================================================================
521|
522| NAME
523|
524| isAlmostSurelyPrime
525|
526| DESCRIPTION
527|
528| True if n is probabilistically prime.
529|
530+============================================================================*/
531
532template <typename IntType>
533bool isAlmostSurelyPrime( const IntType & n ) ;
534
535#endif // __PPFACTOR_H__
