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