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