mpprime.c (15056B)
1 /* 2 * mpprime.c 3 * 4 * Utilities for finding and working with prime and pseudo-prime 5 * integers 6 * 7 * This Source Code Form is subject to the terms of the Mozilla Public 8 * License, v. 2.0. If a copy of the MPL was not distributed with this 9 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ 10 11 #include "mpi-priv.h" 12 #include "mpprime.h" 13 #include "mplogic.h" 14 #include <stdlib.h> 15 #include <string.h> 16 17 #define SMALL_TABLE 0 /* determines size of hard-wired prime table */ 18 19 #define RANDOM() rand() 20 21 #include "primes.c" /* pull in the prime digit table */ 22 23 /* 24 Test if any of a given vector of digits divides a. If not, MP_NO 25 is returned; otherwise, MP_YES is returned and 'which' is set to 26 the index of the integer in the vector which divided a. 27 */ 28 mp_err s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which); 29 30 /* {{{ mpp_divis(a, b) */ 31 32 /* 33 mpp_divis(a, b) 34 35 Returns MP_YES if a is divisible by b, or MP_NO if it is not. 36 */ 37 38 mp_err 39 mpp_divis(mp_int *a, mp_int *b) 40 { 41 mp_err res; 42 mp_int rem; 43 44 if ((res = mp_init(&rem)) != MP_OKAY) 45 return res; 46 47 if ((res = mp_mod(a, b, &rem)) != MP_OKAY) 48 goto CLEANUP; 49 50 if (mp_cmp_z(&rem) == 0) 51 res = MP_YES; 52 else 53 res = MP_NO; 54 55 CLEANUP: 56 mp_clear(&rem); 57 return res; 58 59 } /* end mpp_divis() */ 60 61 /* }}} */ 62 63 /* {{{ mpp_divis_d(a, d) */ 64 65 /* 66 mpp_divis_d(a, d) 67 68 Return MP_YES if a is divisible by d, or MP_NO if it is not. 69 */ 70 71 mp_err 72 mpp_divis_d(mp_int *a, mp_digit d) 73 { 74 mp_err res; 75 mp_digit rem; 76 77 ARGCHK(a != NULL, MP_BADARG); 78 79 if (d == 0) 80 return MP_NO; 81 82 if ((res = mp_mod_d(a, d, &rem)) != MP_OKAY) 83 return res; 84 85 if (rem == 0) 86 return MP_YES; 87 else 88 return MP_NO; 89 90 } /* end mpp_divis_d() */ 91 92 /* }}} */ 93 94 /* {{{ mpp_random(a) */ 95 96 /* 97 mpp_random(a) 98 99 Assigns a random value to a. This value is generated using the 100 standard C library's rand() function, so it should not be used for 101 cryptographic purposes, but it should be fine for primality testing, 102 since all we really care about there is good statistical properties. 103 104 As many digits as a currently has are filled with random digits. 105 */ 106 107 mp_err 108 mpp_random(mp_int *a) 109 110 { 111 mp_digit next = 0; 112 unsigned int ix, jx; 113 114 ARGCHK(a != NULL, MP_BADARG); 115 116 for (ix = 0; ix < USED(a); ix++) { 117 for (jx = 0; jx < sizeof(mp_digit); jx++) { 118 next = (next << CHAR_BIT) | (RANDOM() & UCHAR_MAX); 119 } 120 DIGIT(a, ix) = next; 121 } 122 123 return MP_OKAY; 124 125 } /* end mpp_random() */ 126 127 /* }}} */ 128 129 static mpp_random_fn mpp_random_insecure = &mpp_random; 130 131 /* {{{ mpp_random_size(a, prec) */ 132 133 mp_err 134 mpp_random_size(mp_int *a, mp_size prec) 135 { 136 mp_err res; 137 138 ARGCHK(a != NULL && prec > 0, MP_BADARG); 139 140 if ((res = s_mp_pad(a, prec)) != MP_OKAY) 141 return res; 142 143 return (*mpp_random_insecure)(a); 144 145 } /* end mpp_random_size() */ 146 147 /* }}} */ 148 149 /* {{{ mpp_divis_vector(a, vec, size, which) */ 150 151 /* 152 mpp_divis_vector(a, vec, size, which) 153 154 Determines if a is divisible by any of the 'size' digits in vec. 155 Returns MP_YES and sets 'which' to the index of the offending digit, 156 if it is; returns MP_NO if it is not. 157 */ 158 159 mp_err 160 mpp_divis_vector(mp_int *a, const mp_digit *vec, int size, int *which) 161 { 162 ARGCHK(a != NULL && vec != NULL && size > 0, MP_BADARG); 163 164 return s_mpp_divp(a, vec, size, which); 165 166 } /* end mpp_divis_vector() */ 167 168 /* }}} */ 169 170 /* {{{ mpp_divis_primes(a, np) */ 171 172 /* 173 mpp_divis_primes(a, np) 174 175 Test whether a is divisible by any of the first 'np' primes. If it 176 is, returns MP_YES and sets *np to the value of the digit that did 177 it. If not, returns MP_NO. 178 */ 179 mp_err 180 mpp_divis_primes(mp_int *a, mp_digit *np) 181 { 182 int size, which; 183 mp_err res; 184 185 ARGCHK(a != NULL && np != NULL, MP_BADARG); 186 187 size = (int)*np; 188 if (size > prime_tab_size) 189 size = prime_tab_size; 190 191 res = mpp_divis_vector(a, prime_tab, size, &which); 192 if (res == MP_YES) 193 *np = prime_tab[which]; 194 195 return res; 196 197 } /* end mpp_divis_primes() */ 198 199 /* }}} */ 200 201 /* {{{ mpp_fermat(a, w) */ 202 203 /* 204 Using w as a witness, try pseudo-primality testing based on Fermat's 205 little theorem. If a is prime, and (w, a) = 1, then w^a == w (mod 206 a). So, we compute z = w^a (mod a) and compare z to w; if they are 207 equal, the test passes and we return MP_YES. Otherwise, we return 208 MP_NO. 209 */ 210 mp_err 211 mpp_fermat(mp_int *a, mp_digit w) 212 { 213 mp_int base, test; 214 mp_err res; 215 216 if ((res = mp_init(&base)) != MP_OKAY) 217 return res; 218 219 mp_set(&base, w); 220 221 if ((res = mp_init(&test)) != MP_OKAY) 222 goto TEST; 223 224 /* Compute test = base^a (mod a) */ 225 if ((res = mp_exptmod(&base, a, a, &test)) != MP_OKAY) 226 goto CLEANUP; 227 228 if (mp_cmp(&base, &test) == 0) 229 res = MP_YES; 230 else 231 res = MP_NO; 232 233 CLEANUP: 234 mp_clear(&test); 235 TEST: 236 mp_clear(&base); 237 238 return res; 239 240 } /* end mpp_fermat() */ 241 242 /* }}} */ 243 244 /* 245 Perform the fermat test on each of the primes in a list until 246 a) one of them shows a is not prime, or 247 b) the list is exhausted. 248 Returns: MP_YES if it passes tests. 249 MP_NO if fermat test reveals it is composite 250 Some MP error code if some other error occurs. 251 */ 252 mp_err 253 mpp_fermat_list(mp_int *a, const mp_digit *primes, mp_size nPrimes) 254 { 255 mp_err rv = MP_YES; 256 257 while (nPrimes-- > 0 && rv == MP_YES) { 258 rv = mpp_fermat(a, *primes++); 259 } 260 return rv; 261 } 262 263 /* {{{ mpp_pprime(a, nt) */ 264 265 /* 266 mpp_pprime(a, nt) 267 268 Performs nt iteration of the Miller-Rabin probabilistic primality 269 test on a. Returns MP_YES if the tests pass, MP_NO if one fails. 270 If MP_NO is returned, the number is definitely composite. If MP_YES 271 is returned, it is probably prime (but that is not guaranteed). 272 */ 273 274 mp_err 275 mpp_pprime(mp_int *a, int nt) 276 { 277 return mpp_pprime_ext_random(a, nt, mpp_random_insecure); 278 } 279 280 mp_err 281 mpp_pprime_ext_random(mp_int *a, int nt, mpp_random_fn random) 282 { 283 mp_err res; 284 mp_int x, amo, m, z; /* "amo" = "a minus one" */ 285 int iter; 286 unsigned int jx; 287 mp_size b; 288 289 ARGCHK(a != NULL, MP_BADARG); 290 291 MP_DIGITS(&x) = 0; 292 MP_DIGITS(&amo) = 0; 293 MP_DIGITS(&m) = 0; 294 MP_DIGITS(&z) = 0; 295 296 /* Initialize temporaries... */ 297 MP_CHECKOK(mp_init(&amo)); 298 /* Compute amo = a - 1 for what follows... */ 299 MP_CHECKOK(mp_sub_d(a, 1, &amo)); 300 301 b = mp_trailing_zeros(&amo); 302 if (!b) { /* a was even ? */ 303 res = MP_NO; 304 goto CLEANUP; 305 } 306 307 MP_CHECKOK(mp_init_size(&x, MP_USED(a))); 308 MP_CHECKOK(mp_init(&z)); 309 MP_CHECKOK(mp_init(&m)); 310 MP_CHECKOK(mp_div_2d(&amo, b, &m, 0)); 311 312 /* Do the test nt times... */ 313 for (iter = 0; iter < nt; iter++) { 314 315 /* Choose a random value for 1 < x < a */ 316 MP_CHECKOK(s_mp_pad(&x, USED(a))); 317 MP_CHECKOK((*random)(&x)); 318 MP_CHECKOK(mp_mod(&x, a, &x)); 319 if (mp_cmp_d(&x, 1) <= 0) { 320 iter--; /* don't count this iteration */ 321 continue; /* choose a new x */ 322 } 323 324 /* Compute z = (x ** m) mod a */ 325 MP_CHECKOK(mp_exptmod(&x, &m, a, &z)); 326 327 if (mp_cmp_d(&z, 1) == 0 || mp_cmp(&z, &amo) == 0) { 328 res = MP_YES; 329 continue; 330 } 331 332 res = MP_NO; /* just in case the following for loop never executes. */ 333 for (jx = 1; jx < b; jx++) { 334 /* z = z^2 (mod a) */ 335 MP_CHECKOK(mp_sqrmod(&z, a, &z)); 336 res = MP_NO; /* previous line set res to MP_YES */ 337 338 if (mp_cmp_d(&z, 1) == 0) { 339 break; 340 } 341 if (mp_cmp(&z, &amo) == 0) { 342 res = MP_YES; 343 break; 344 } 345 } /* end testing loop */ 346 347 /* If the test passes, we will continue iterating, but a failed 348 test means the candidate is definitely NOT prime, so we will 349 immediately break out of this loop 350 */ 351 if (res == MP_NO) 352 break; 353 354 } /* end iterations loop */ 355 356 CLEANUP: 357 mp_clear(&m); 358 mp_clear(&z); 359 mp_clear(&x); 360 mp_clear(&amo); 361 return res; 362 363 } /* end mpp_pprime() */ 364 365 /* }}} */ 366 367 /* Produce table of composites from list of primes and trial value. 368 ** trial must be odd. List of primes must not include 2. 369 ** sieve should have dimension >= MAXPRIME/2, where MAXPRIME is largest 370 ** prime in list of primes. After this function is finished, 371 ** if sieve[i] is non-zero, then (trial + 2*i) is composite. 372 ** Each prime used in the sieve costs one division of trial, and eliminates 373 ** one or more values from the search space. (3 eliminates 1/3 of the values 374 ** alone!) Each value left in the search space costs 1 or more modular 375 ** exponentations. So, these divisions are a bargain! 376 */ 377 mp_err 378 mpp_sieve(mp_int *trial, const mp_digit *primes, mp_size nPrimes, 379 unsigned char *sieve, mp_size nSieve) 380 { 381 mp_err res; 382 mp_digit rem; 383 mp_size ix; 384 unsigned long offset; 385 386 memset(sieve, 0, nSieve); 387 388 for (ix = 0; ix < nPrimes; ix++) { 389 mp_digit prime = primes[ix]; 390 mp_size i; 391 if ((res = mp_mod_d(trial, prime, &rem)) != MP_OKAY) 392 return res; 393 394 if (rem == 0) { 395 offset = 0; 396 } else { 397 offset = prime - rem; 398 } 399 400 for (i = offset; i < nSieve * 2; i += prime) { 401 if (i % 2 == 0) { 402 sieve[i / 2] = 1; 403 } 404 } 405 } 406 407 return MP_OKAY; 408 } 409 410 #define SIEVE_SIZE 32 * 1024 411 412 mp_err 413 mpp_make_prime(mp_int *start, mp_size nBits, mp_size strong) 414 { 415 return mpp_make_prime_ext_random(start, nBits, strong, mpp_random_insecure); 416 } 417 418 mp_err 419 mpp_make_prime_ext_random(mp_int *start, mp_size nBits, mp_size strong, mpp_random_fn random) 420 { 421 mp_digit np; 422 mp_err res; 423 unsigned int i = 0; 424 mp_int trial; 425 mp_int q; 426 mp_size num_tests; 427 unsigned char *sieve; 428 429 ARGCHK(start != 0, MP_BADARG); 430 ARGCHK(nBits > 16, MP_RANGE); 431 432 sieve = malloc(SIEVE_SIZE); 433 ARGCHK(sieve != NULL, MP_MEM); 434 435 MP_DIGITS(&trial) = 0; 436 MP_DIGITS(&q) = 0; 437 MP_CHECKOK(mp_init(&trial)); 438 MP_CHECKOK(mp_init(&q)); 439 /* values originally taken from table 4.4, 440 * HandBook of Applied Cryptography, augmented by FIPS-186 441 * requirements, Table C.2 and C.3 */ 442 if (nBits >= 2000) { 443 num_tests = 3; 444 } else if (nBits >= 1536) { 445 num_tests = 4; 446 } else if (nBits >= 1024) { 447 num_tests = 5; 448 } else if (nBits >= 550) { 449 num_tests = 6; 450 } else if (nBits >= 450) { 451 num_tests = 7; 452 } else if (nBits >= 400) { 453 num_tests = 8; 454 } else if (nBits >= 350) { 455 num_tests = 9; 456 } else if (nBits >= 300) { 457 num_tests = 10; 458 } else if (nBits >= 250) { 459 num_tests = 20; 460 } else if (nBits >= 200) { 461 num_tests = 41; 462 } else if (nBits >= 100) { 463 num_tests = 38; /* funny anomaly in the FIPS tables, for aux primes, the 464 * required more iterations for larger aux primes */ 465 } else 466 num_tests = 50; 467 468 if (strong) 469 --nBits; 470 MP_CHECKOK(mpl_set_bit(start, nBits - 1, 1)); 471 MP_CHECKOK(mpl_set_bit(start, 0, 1)); 472 for (i = mpl_significant_bits(start) - 1; i >= nBits; --i) { 473 MP_CHECKOK(mpl_set_bit(start, i, 0)); 474 } 475 /* start sieveing with prime value of 3. */ 476 MP_CHECKOK(mpp_sieve(start, prime_tab + 1, prime_tab_size - 1, 477 sieve, SIEVE_SIZE)); 478 479 #ifdef DEBUG_SIEVE 480 res = 0; 481 for (i = 0; i < SIEVE_SIZE; ++i) { 482 if (!sieve[i]) 483 ++res; 484 } 485 fprintf(stderr, "sieve found %d potential primes.\n", res); 486 #define FPUTC(x, y) fputc(x, y) 487 #else 488 #define FPUTC(x, y) 489 #endif 490 491 res = MP_NO; 492 for (i = 0; i < SIEVE_SIZE; ++i) { 493 if (sieve[i]) /* this number is composite */ 494 continue; 495 MP_CHECKOK(mp_add_d(start, 2 * i, &trial)); 496 FPUTC('.', stderr); 497 /* run a Fermat test */ 498 res = mpp_fermat(&trial, 2); 499 if (res != MP_OKAY) { 500 if (res == MP_NO) 501 continue; /* was composite */ 502 goto CLEANUP; 503 } 504 505 FPUTC('+', stderr); 506 /* If that passed, run some Miller-Rabin tests */ 507 res = mpp_pprime_ext_random(&trial, num_tests, random); 508 if (res != MP_OKAY) { 509 if (res == MP_NO) 510 continue; /* was composite */ 511 goto CLEANUP; 512 } 513 FPUTC('!', stderr); 514 515 if (!strong) 516 break; /* success !! */ 517 518 /* At this point, we have strong evidence that our candidate 519 is itself prime. If we want a strong prime, we need now 520 to test q = 2p + 1 for primality... 521 */ 522 MP_CHECKOK(mp_mul_2(&trial, &q)); 523 MP_CHECKOK(mp_add_d(&q, 1, &q)); 524 525 /* Test q for small prime divisors ... */ 526 np = prime_tab_size; 527 res = mpp_divis_primes(&q, &np); 528 if (res == MP_YES) { /* is composite */ 529 mp_clear(&q); 530 continue; 531 } 532 if (res != MP_NO) 533 goto CLEANUP; 534 535 /* And test with Fermat, as with its parent ... */ 536 res = mpp_fermat(&q, 2); 537 if (res != MP_YES) { 538 mp_clear(&q); 539 if (res == MP_NO) 540 continue; /* was composite */ 541 goto CLEANUP; 542 } 543 544 /* And test with Miller-Rabin, as with its parent ... */ 545 res = mpp_pprime_ext_random(&q, num_tests, random); 546 if (res != MP_YES) { 547 mp_clear(&q); 548 if (res == MP_NO) 549 continue; /* was composite */ 550 goto CLEANUP; 551 } 552 553 /* If it passed, we've got a winner */ 554 mp_exch(&q, &trial); 555 mp_clear(&q); 556 break; 557 558 } /* end of loop through sieved values */ 559 if (res == MP_YES) 560 mp_exch(&trial, start); 561 CLEANUP: 562 mp_clear(&trial); 563 mp_clear(&q); 564 if (sieve != NULL) { 565 memset(sieve, 0, SIEVE_SIZE); 566 free(sieve); 567 } 568 return res; 569 } 570 571 /*========================================================================*/ 572 /*------------------------------------------------------------------------*/ 573 /* Static functions visible only to the library internally */ 574 575 /* {{{ s_mpp_divp(a, vec, size, which) */ 576 577 /* 578 Test for divisibility by members of a vector of digits. Returns 579 MP_NO if a is not divisible by any of them; returns MP_YES and sets 580 'which' to the index of the offender, if it is. Will stop on the 581 first digit against which a is divisible. 582 */ 583 584 mp_err 585 s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which) 586 { 587 mp_err res; 588 mp_digit rem; 589 590 int ix; 591 592 for (ix = 0; ix < size; ix++) { 593 if ((res = mp_mod_d(a, vec[ix], &rem)) != MP_OKAY) 594 return res; 595 596 if (rem == 0) { 597 if (which) 598 *which = ix; 599 return MP_YES; 600 } 601 } 602 603 return MP_NO; 604 605 } /* end s_mpp_divp() */ 606 607 /* }}} */ 608 609 /*------------------------------------------------------------------------*/ 610 /* HERE THERE BE DRAGONS */