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 ================ */
