1 /*==============================================================================
  2 | 
  3 |  NAME     
  4 |
  5 |     ppFactor.h
  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.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 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 
 57 class 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 
 81 class 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 
117 template <typename IntType>
118 class 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 
176 template <typename IntType>
177 class 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 
199 template <typename IntType>
200 class 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.
251 enum 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.
262 template <typename IntType>
263 class 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 
338 enum 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 
360 template <typename IntType>
361 Primality 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 
383 template <typename IntType>
384 class 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 // Static class variables must be initialized outside the class.
511 template <typename IntType> int UniformRandomIntegers<IntType>::num_of_initializations = 0 ;
512 
513 
514 /*=============================================================================
515 |
516 | NAME
517 |
518 |     isAlmostSurelyPrime
519 |
520 | DESCRIPTION
521 |
522 |     True if n is probabilistically prime.
523 |
524 +============================================================================*/
525 
526 template <typename IntType>
527 bool isAlmostSurelyPrime( const IntType & n ) ;
528 
529 #endif // __PPFACTOR_H__