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.4 - A Program for Computing Primitive Polynomials.
20| Copyright (C) 1999-2025 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 seanerikoconnor!AT!gmail!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
54DESCRIPTION
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
59INPUT
60
61 n (int, n >= 1) The integer to be factored.
62
63OUTPUT
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
70RETURNS
71
72 t (int) Number of prime factors - 1.
73 Prime factors and their multiplicites are in locations 0 to t.
74
75EXAMPLE
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
85METHOD
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
171BUGS
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
183int
184factor( bigint n, bigint * primes, int * count )
185{
186
187/*------------------------------------------------------------------------------
188| Local Variables |
189------------------------------------------------------------------------------*/
190
191int
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
199bigint
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*/
213if (n == 1)
214{
215 primes[ 0 ] = count[ 0 ] = 1 ;
216 return( 0 ) ;
217}
218else 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}
231else 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
240while( n % 2 == 0 ) /* Remove factors of 2. */
241{
242 n /= 2 ;
243 ++cnt ;
244}
245
246if (cnt != 0)
247{
248 primes[ t ] = 2 ;
249 count[ t++ ] = cnt ;
250}
251
252cnt = 0 ;
253
254while( n % 3 == 0 ) /* Remove factors of 3. */
255{
256 n /= 3 ;
257 ++cnt ;
258}
259
260if (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
272do {
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
306if (n == 1) /* All factors were divided out. */
307
308 return( t - 1 ) ;
309
310else { /* 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
321Phi[ 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
336bigint
337EulerPhi( 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
377DESCRIPTION
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
384INPUT
385
386 n (int, n >= 0) Number to test for primality.
387 x (int, 1 < x < n for n > 6) A random integer.
388
389RETURNS
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
396EXAMPLE
397
398METHOD
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
408BUGS
409
410 None.
411
412--------------------------------------------------------------------------------
413| Function Call |
414------------------------------------------------------------------------------*/
415
416int is_probably_prime( int n, int x )
417{
418
419/*------------------------------------------------------------------------------
420| Local Variables |
421------------------------------------------------------------------------------*/
422
423int
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. */
435if (n < 0)
436 return NO ;
437
438/* Handle special cases. */
439if (n == 0 || n == 1 || n == 4)
440 return NO ;
441
442if (n == 2 || n == 3 || n == 5)
443 return YES ;
444
445/* Even numbers aren't prime. */
446if (n % 2 == 0)
447 return NO ;
448
449/* Return not prime if x is out of range. */
450if (x <= 1 || x >= n)
451 return NO ;
452
453/* Factor out powers of 2 to get n = 1 + 2^k q, q odd. */
454nm1 = n - 1 ;
455k = 0 ;
456while( nm1 % 2 == 0)
457{
458 nm1 /= 2 ;
459 ++k ;
460}
461
462q = nm1 ;
463
464
465
466/* y = x^q (mod n) */
467y = power_mod( x, q, n ) ;
468
469for (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. */
481return NO ;
482
483} /* =============== end of function is_probably_prime ===================== */
484
485
486
487/*==============================================================================
488| is_almost_surely_prime |
489================================================================================
490
491DESCRIPTION
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
500INPUT
501
502 n (int, n >= 0)
503
504RETURNS
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
512EXAMPLE
513
514METHOD
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
528BUGS
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
537int is_almost_surely_prime( int n )
538{
539
540/*------------------------------------------------------------------------------
541| Local Variables |
542------------------------------------------------------------------------------*/
543
544int
545 trial = 0,
546 x = 3 ;
547
548/*------------------------------------------------------------------------------
549| Function Body |
550------------------------------------------------------------------------------*/
551
552/* Seed the random-number generator. */
553srand( 314159 );
554
555for (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 */
569return YES ;
570
571} /* ============== end of function is_almost_surely_prime ================ */
