1 /*==============================================================================
  2 | 
  3 |  File Name:     
  4 |
  5 |     ppFactor.c
  6 |
  7 |  Description:   
  8 |
  9 |     Routines for integer factorization and primality prime testing.
 10 | 
 11 |  Functions:
 12 |
 13 |     factor
 14 |     is_probably_prime
 15 |     is_almost_surely_prime
 16 |
 17 |  LEGAL
 18 |
 19 |     Primpoly Version 16.1 - A Program for Computing Primitive Polynomials.
 20 |     Copyright (C) 1999-2021 by Sean Erik O'Connor.  All Rights Reserved.
 21 |
 22 |     This program is free software: you can redistribute it and/or modify
 23 |     it under the terms of the GNU General Public License as published by
 24 |     the Free Software Foundation, either version 3 of the License, or
 25 |     (at your option) any later version.
 26 |
 27 |     This program is distributed in the hope that it will be useful,
 28 |     but WITHOUT ANY WARRANTY; without even the implied warranty of
 29 |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 30 |     GNU General Public License for more details.
 31 |
 32 |     You should have received a copy of the GNU General Public License
 33 |     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 34 |     
 35 |     The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com
 36 |     with !DOT! replaced by . and the !AT! replaced by @
 37 |
 38 ==============================================================================*/
 39 
 40 /*------------------------------------------------------------------------------
 41 |                                Include Files                                 |
 42 ------------------------------------------------------------------------------*/
 43 
 44 #include <stdio.h>
 45 #include <stdlib.h> /* for rand() function */
 46 
 47 #include "Primpoly.h"
 48 
 49 
 50 /*==============================================================================
 51 |                                    factor                                    |
 52 ================================================================================
 53 
 54 DESCRIPTION
 55 
 56      Factor n into primes.  Record all the distinct prime factors and how
 57      many times each occurs.  Return the number of primes - 1.
 58 
 59 INPUT
 60                    
 61      n      (int, n >= 1)     The integer to be factored.
 62 
 63 OUTPUT
 64 
 65      primes (bigint *)   List of distinct prime factors.
 66      count  (int *)      List of how many times each factor occurred.
 67 
 68      When n = 1, t = 0, and primes[ 0 ] = count[ 0 ] = 1.
 69 
 70 RETURNS
 71 
 72      t (int)             Number of prime factors - 1.  
 73                          Prime factors and their multiplicites are in locations 0 to t.
 74 
 75 EXAMPLE 
 76                    2
 77     For n = 156 = 2  * 3 * 13 we have
 78 
 79     k   primes[ k ]   count[ k ]
 80     ----------------------------
 81     0        2            2
 82     1        3            1
 83     2       13            1
 84 
 85 METHOD
 86 
 87     Method described by D.E. Knuth, ART OF COMPUTER PROGRAMMING, vol. 2, 2nd ed.,
 88                                                                 -----
 89     in Algorithm A, pgs. 364-365.  The running time is O( max \/ p   , p  )
 90                                                                   t-1   t
 91     where p  is the largest prime divisor of n and p    is the next largest.
 92            t                                        t-1
 93 
 94     (1)  First divide out all multiples of 2 and 3 and count them.  
 95     
 96     (2)  Next, divide n by all integers d >= 5 except multiples of 2 and 3. 
 97 
 98     (3)  Halt either when all prime factors have been divided out (leaving n = 1) 
 99          or when the current value of n is prime.  The stopping test 
100 
101           (d > | n/d | AND r != 0)
102                --   --
103          detects when n is prime.
104 
105     In the example above, we divided out 2's and 3's first, leaving
106     n = 13.  We continued with trial divisor d = 3.  But now the stopping
107     test was activated because 5 > | 13/5 | = 2, and remainder = 3 (non-zero).
108                                    --    --
109     n = 13 itself is the last prime factor of 156.  We return t = 2.  There
110     are 2 + 1 = 3 distinct primes.
111 
112 
113     If we start with d = 5, and add 2 and 4 in succession, we will run through all
114     the integers except multiples of 2 and 3.  To see this, observe the pattern:
115 
116     Integers:     1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17
117     Mult. of 2:      x     x     x     x     x     x     x     x
118     Mult. of 3:         x        x        x        x        x
119     d:                        x     x           x     x          x
120     Separation:                  2       4         2       4
121 
122     The sequence of divisors d, includes all the primes, so is safe to use for
123     factorization testing.
124 
125     Theorem.  Suppose no divisor d in the sequence divides n.  Suppose at some point,
126     q < d with r != 0.  Then n must be prime.
127 
128     Proof.  We begin with an integer n which has all powers of 2 and 3 removed.
129     Assume, to the contrary, that n is composite:  n = p ... p
130           t                                             1     t
131     n >= p    since p  is the smallest prime factor.
132           1          1
133 
134           2                                      
135     n >= p    since n is composite, so has at least 2 prime factors.
136           1
137                 2
138       >= (d + 1)   since p  > d implies p >= (d + 1).  p > d because we assumed
139                           1              1              1
140       that no d in the sequence divided n, so d couldn't be any of the prime divisors
141       p
142        i
143          2            2               2
144       = d + 2d + 1 = d + d + d + 1 > d + d
145          2
146       > d + n mod d   since n mod d < d
147                                                                               2
148       > | n / d | d + n mod d  because our test said q < d, so | n / d | d < d
149         --     --                                              --     --
150 
151       = n.  So we get the contradiction n > n.  Thus n must be prime. 
152                                                    ---
153     Note that this is sharper than the bound d > \/ n
154 
155     Theorem.  Conversely, if n is a prime, no d will divide it, so at some point, 
156     d will be large enough that q = | n / d | < d.  r != 0 since n is prime.
157                                     --     --
158 
159     Theorem.  The factoring algorithm works.
160     
161     Proof.  If n == 1 we exit immediately.  If n is composite, we divide out all powers
162     of 2 and 3 (if any).  Since d runs through all possible prime divisors, we divide
163     these out also.  Composite trial d causes no problem;  prime factors of d are divided
164     out of n before we try to divide by d, so such a d cannot divide n.
165 
166     We end in one of two ways (1) All divisors of n have been divided out in which case n = 1.
167     (2) n is a prime, so the stopping test is activiated and n != 1 is recorded as a prime
168     divisor.
169 
170 
171 BUGS
172 
173     Can be slow when n is a prime.  We could do a probabilistic test for
174     the primality of n at the statement, "n_is_prime = (r != 0 && q < d) ;"
175     which would speed up the test.  
176     
177     Or we could implement the Pollard rho algorithm instead.
178 
179 --------------------------------------------------------------------------------
180 |                                Function Call                                 |
181 ------------------------------------------------------------------------------*/
182 
183 int
184 factor( bigint n, bigint * primes, int * count )
185 {
186 
187 /*------------------------------------------------------------------------------
188 |                               Local Variables                                |
189 ------------------------------------------------------------------------------*/
190 
191 int
192     t = 0,              /*  Array index for primes and count.         */
193     cnt = 0,            /*  Counter for factors of 2 and 3.           */
194     new_d = YES,        /*  YES if current divisor is different from
195                             the previous one.                         */
196     n_is_prime = NO,    /*  YES if current value of n is prime.       */
197     flag = YES ;        /*  Alternately YES and NO.                   */
198 
199 bigint
200     q,                  /*  The quotient, floor( n / d )              */
201     r,                  /*  The remainder, n mod d                    */
202     d = 5 ;             /*  Trial divisor and its first value.        */
203 
204 /*------------------------------------------------------------------------------
205 |                                Function Body                                 |
206 ------------------------------------------------------------------------------*/
207 
208 /*  Handle special cases for speed. 
209     These are cases where r = (p^n - 1)/(p-1) is a large prime or has large prime 
210     factors and p^n is less than the maximum size allowed by this program.
211     (I don't have all these cases).
212 */
213 if (n == 1)
214 {
215     primes[ 0 ] = count[ 0 ] = 1 ;
216     return( 0 ) ;
217 }
218 else if (n == (bigint)4611686018427387903ULL) /*  (2^62 - 1)/(2 - 1)  */
219 {
220     primes[ 0 ] = 3 ;
221     count [ 0 ] = 1 ;
222 
223     primes[ 1 ] = 715827883ULL ;
224     count [ 1 ] = 1 ;
225 
226     primes[ 2 ] = 2147483647ULL ;
227     count [ 2 ] = 1 ;
228 
229     return( 2 ) ;
230 }
231 else if (n == (bigint)2305843009213693951ULL) /*  (2^61 - 1)/(2 - 1)  */
232 {
233     primes[ 0 ] = (bigint)2305843009213693951ULL ;
234     count [ 0 ] = 1 ;
235 
236     return( 0 ) ;
237 }
238 
239 
240 while( n % 2 == 0 )       /*  Remove factors of 2. */
241 {
242     n /= 2 ;
243     ++cnt ;
244 }
245 
246 if (cnt != 0)
247 {
248     primes[ t ] = 2 ;
249     count[ t++ ] = cnt ;
250 }
251 
252 cnt = 0 ;
253 
254 while( n % 3 == 0 )       /*  Remove factors of 3. */
255 {
256     n /= 3 ;
257     ++cnt ;
258 }
259 
260 if (cnt != 0)
261 {
262     primes[ t ] = 3 ;
263     count[ t++ ] = cnt ;
264 }
265 
266 
267 /*
268     Factor the rest of n.  Continue until n = 1 (all factors divided out)
269     or until n is prime (so n itself is the last prime factor).
270 */
271 
272 do {
273 
274     /*  Integer divide to get quotient and remainder:  n = qd + r */
275 
276     q = n / d ;
277 
278     r = n % d ;
279 
280     n_is_prime = (r != 0 && q < d) ;   /* n is prime. */
281 
282     if (r == 0)
283     {
284         n = q ;     /*  Divide d out of n. */
285 
286         if (new_d)  /*  New prime factor. */
287         {
288             primes[ t ] = d ;
289             count[ t++ ] = 1 ;
290             new_d = NO ;
291          }
292          else
293 
294              ++count[ t-1 ] ;   /*  Same divisor again.  Increment its count. */
295     }
296     else {
297 
298         d += (flag ? 2 : 4) ;  /* d did not divide n.  Try a new divisor. */
299 
300         flag = !flag ;
301         new_d = YES ;
302     }
303 
304 } while ( ! n_is_prime && (n != 1)) ;
305 
306 if (n == 1)       /*  All factors were divided out. */
307 
308     return( t - 1 ) ;
309 
310 else {  /*  Current value of n is prime.  It is the last prime factor. */
311 
312     primes[ t ] = n ;
313     count[ t ] = 1 ;
314     return( t ) ;
315 }
316 
317 } /* ========================= end of function factor ======================= */
318 
319 /*
320       4
321 Phi[ 5  - 1 ] / 4 
322 
323    4              4
324   5  - 1 = 624 = 2  3  13
325 
326   Phi = 624 (1 - 1/2) (1 - 1/3)(1 - 1/13) = 192
327 
328   192 / 4 = 48
329 
330 
331   Handle special cases.
332 
333   Warning: can take a long time due to factorization.
334   */
335 
336 bigint
337 EulerPhi( bigint n )
338 {
339     int prime_count = 0 ;
340     int i = 0 ;
341     bigint phi ;
342 
343     /* Factor n into primes. */
344     bigint primes[ MAXNUMPRIMEFACTORS ] ;
345     int    count [ MAXNUMPRIMEFACTORS ] ;
346 
347     /* Error return and special cases. */
348     if (n <= 0)
349         return 0 ;
350     else if (n == 1)
351         return 1 ;
352 
353     /* Factor n >= 2 into distinct primes. */
354     prime_count = factor( n, primes, count ) ;
355 
356     /* Compute Euler phi[ n ] =   */
357     /*                            */
358     /*      num prime factors     */
359     /*        -------             */
360     /*     n   |  |   (1 - 1/p )  */
361     /*         i = 1          i   */
362     phi = n ;
363     for (i = 0 ;  i <= prime_count ; ++i)
364     {
365         phi /=  primes[ i ] ;
366         phi *= (primes[ i ] - 1) ;
367     }
368 
369     return phi ;
370 }
371 
372 
373 /*==============================================================================
374 |                             is_probably_prime                                |
375 ================================================================================
376 
377 DESCRIPTION
378 
379      Test whether the number n is probably prime.  If n is composite
380      we are correct always.  If n is prime, we will be wrong no more
381      than about 1/4 of the time on average, probably less for
382      any x and n.
383 
384 INPUT
385                    
386      n      (int, n >= 0)    Number to test for primality.
387      x      (int, 1 < x < n for n > 6) A random integer.
388 
389 RETURNS
390 
391 
392      YES    If n is probably prime.
393      NO     If n is definitely not prime or n < 0,
394             or x is out of range.
395 
396 EXAMPLE 
397  
398 METHOD
399 
400      Miller-Rabin probabilistic primality test, described by Knuth.
401      If n = 1 + 2^k q is prime, and x^q mod n != 1 the sequence, 
402      { x^q mod n, x^(2q) mod n, ..., x^(2^k q) mod n} 
403      will end with 1 and the element in the sequence just before 
404      1 first appears = n-1, since y^2 = 1 (mod p) satisfies 
405      y = +1 or y = -1 only. The probability the algorithm fails is 
406      bounded above by about 1/4 for all n.
407 
408 BUGS
409 
410     None.
411 
412 --------------------------------------------------------------------------------
413 |                                Function Call                                 |
414 ------------------------------------------------------------------------------*/
415 
416 int is_probably_prime( int n, int x )
417 {
418 
419 /*------------------------------------------------------------------------------
420 |                               Local Variables                                |
421 ------------------------------------------------------------------------------*/
422 
423 int
424     j,
425     k = 0,
426     nm1,
427     q = 0,
428     y = 0 ;
429 
430 /*------------------------------------------------------------------------------
431 |                                Function Body                                 |
432 ------------------------------------------------------------------------------*/
433 
434 /*  Not prime if input is out of range. */
435 if (n < 0)
436     return NO ;
437 
438 /*  Handle special cases. */
439 if (n == 0 || n == 1 || n == 4)
440     return NO ;
441 
442 if (n == 2 || n == 3 || n == 5)
443     return YES ;
444 
445 /* Even numbers aren't prime. */
446 if (n % 2 == 0)
447     return NO ;
448 
449 /* Return not prime if x is out of range. */
450 if (x <= 1 || x >= n)
451     return NO ;
452 
453 /* Factor out powers of 2 to get n = 1 + 2^k q, q odd. */
454 nm1 = n - 1 ;
455 k = 0 ;
456 while( nm1 % 2 == 0)
457 {
458     nm1 /= 2 ;
459     ++k ;
460 }
461 
462 q = nm1 ;
463 
464 
465 
466 /*  y = x^q (mod n) */
467 y = power_mod( x, q, n ) ;
468 
469 for (j = 0 ;  j < k ;  ++j)
470 {
471     if ( (j == 0 && y == 1) || (y == n-1))
472         return YES ;
473     else if (j > 0 && y == 1)
474         return NO ;
475 
476     /* Compute y^2 (mod n) */
477     y = power_mod( y, 2, n ) ;
478 }
479 
480 /*  Shouldn't get here, but do it anyway for safety. */
481 return NO ;
482 
483 } /* =============== end of function is_probably_prime ===================== */
484 
485 
486 
487 /*==============================================================================
488 |                           is_almost_surely_prime                             |
489 ================================================================================
490 
491 DESCRIPTION
492 
493      Test whether the number n is almost surely prime.  If n is
494      composite, we always return NO.  If n is prime, the 
495      probability of returning YES successfully is
496 
497                  NUM_PRIME_TEST_TRIALS
498         1 - (1/4)
499 
500 INPUT
501                    
502      n      (int, n >= 0)
503 
504 RETURNS
505 
506 
507      YES    If n is probably prime with a very high degree of
508             probability.
509      NO     If n is definitely not prime.
510      -1     If inputs are out of range.
511 
512 EXAMPLE
513   
514 METHOD
515 
516     Use the Miller-Rabin probabilitic primality test for lots of
517     random integers x.  If the test shows n is composite for any 
518     given x, is is non-prime for sure.  
519     
520     If it passes the the primality test for several random values
521     of x, the probability is it prime is approximately
522     
523          NUM_PRIME_TEST_TRIALS
524     1 - 2
525 
526     assuming the random number generator is really random.
527 
528 BUGS
529 
530     The system pseudorandom number generator needs to be a good one
531     (i.e. passes the spectral test and other statistical tests).
532 
533 --------------------------------------------------------------------------------
534 |                                Function Call                                 |
535 ------------------------------------------------------------------------------*/
536 
537 int is_almost_surely_prime( int n )
538 {
539 
540 /*------------------------------------------------------------------------------
541 |                               Local Variables                                |
542 ------------------------------------------------------------------------------*/
543 
544 int
545     trial = 0,
546     x     = 3 ;
547 
548 /*------------------------------------------------------------------------------
549 |                                Function Body                                 |
550 ------------------------------------------------------------------------------*/
551 
552 /* Seed the random-number generator. */
553 srand( 314159 );
554 
555 for (trial = 1 ;  trial <= NUM_PRIME_TEST_TRIALS ;  ++trial)
556 {
557     /*  Generate a new random integer such that 1 < x < n. */
558     x = rand() % n ;
559     if (x <= 1) x = 3 ;
560 
561     /* Definitely not prime. */
562     if (is_probably_prime( n, x ) == NO)
563         return NO ;
564 }
565 
566 /* Almost surely prime because it passed the probably prime test
567  * above many times.
568  */
569 return YES ;
570 
571 } /* ============== end of function is_almost_surely_prime ================ */