Alexandre Lision | 7fd5d3d | 2013-12-04 13:06:40 -0500 | [diff] [blame] | 1 | /* |
| 2 | * Prime generation using the bignum library and sieving. |
| 3 | * |
| 4 | * Copyright (c) 1995 Colin Plumb. All rights reserved. |
| 5 | * For licensing and other legal details, see the file legal.c. |
| 6 | */ |
| 7 | #ifndef HAVE_CONFIG_H |
| 8 | #define HAVE_CONFIG_H 0 |
| 9 | #endif |
| 10 | #if HAVE_CONFIG_H |
| 11 | #include <bnconfig.h> |
| 12 | #endif |
| 13 | |
| 14 | /* |
| 15 | * Some compilers complain about #if FOO if FOO isn't defined, |
| 16 | * so do the ANSI-mandated thing explicitly... |
| 17 | */ |
| 18 | #ifndef NO_ASSERT_H |
| 19 | #define NO_ASSERT_H 0 |
| 20 | #endif |
| 21 | #if !NO_ASSERT_H |
| 22 | #include <assert.h> |
| 23 | #else |
| 24 | #define assert(x) (void)0 |
| 25 | #endif |
| 26 | |
| 27 | #include <stdarg.h> /* We just can't live without this... */ |
| 28 | |
| 29 | #ifndef BNDEBUG |
| 30 | #define BNDEBUG 1 |
| 31 | #endif |
| 32 | #if BNDEBUG |
| 33 | #include <stdio.h> |
| 34 | #endif |
| 35 | |
| 36 | #include "bn.h" |
| 37 | #include "lbnmem.h" |
| 38 | #include "prime.h" |
| 39 | #include "sieve.h" |
| 40 | |
| 41 | #include "kludge.h" |
| 42 | |
| 43 | /* Size of the shuffle table */ |
| 44 | #define SHUFFLE 256 |
| 45 | /* Size of the sieve area */ |
| 46 | #define SIEVE 32768u/16 |
| 47 | |
| 48 | /* Confirmation tests. The first one *must* be 2 */ |
| 49 | static unsigned const confirm[] = {2, 3, 5, 7, 11, 13, 17}; |
| 50 | #define CONFIRMTESTS (sizeof(confirm)/sizeof(*confirm)) |
| 51 | |
| 52 | /* |
| 53 | * Helper function that does the slow primality test. |
| 54 | * bn is the input bignum; a and e are temporary buffers that are |
| 55 | * allocated by the caller to save overhead. |
| 56 | * |
| 57 | * Returns 0 if prime, >0 if not prime, and -1 on error (out of memory). |
| 58 | * If not prime, returns the number of modular exponentiations performed. |
| 59 | * Calls the given progress function with a '*' for each primality test |
| 60 | * that is passed. |
| 61 | * |
| 62 | * The testing consists of strong pseudoprimality tests, to the bases given |
| 63 | * in the confirm[] array above. (Also called Miller-Rabin, although that's |
| 64 | * not technically correct if we're using fixed bases.) Some people worry |
| 65 | * that this might not be enough. Number theorists may wish to generate |
| 66 | * primality proofs, but for random inputs, this returns non-primes with |
| 67 | * a probability which is quite negligible, which is good enough. |
| 68 | * |
| 69 | * It has been proved (see Carl Pomerance, "On the Distribution of |
| 70 | * Pseudoprimes", Math. Comp. v.37 (1981) pp. 587-593) that the number of |
| 71 | * pseudoprimes (composite numbers that pass a Fermat test to the base 2) |
| 72 | * less than x is bounded by: |
| 73 | * exp(ln(x)^(5/14)) <= P_2(x) ### CHECK THIS FORMULA - it looks wrong! ### |
| 74 | * P_2(x) <= x * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))). |
| 75 | * Thus, the local density of Pseudoprimes near x is at most |
| 76 | * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))), and at least |
| 77 | * exp(ln(x)^(5/14) - ln(x)). Here are some values of this function |
| 78 | * for various k-bit numbers x = 2^k: |
| 79 | * Bits Density <= Bit equivalent Density >= Bit equivalent |
| 80 | * 128 3.577869e-07 21.414396 4.202213e-37 120.840190 |
| 81 | * 192 4.175629e-10 31.157288 4.936250e-56 183.724558 |
| 82 | * 256 5.804314e-13 40.647940 4.977813e-75 246.829095 |
| 83 | * 384 1.578039e-18 59.136573 3.938861e-113 373.400096 |
| 84 | * 512 5.858255e-24 77.175803 2.563353e-151 500.253110 |
| 85 | * 768 1.489276e-34 112.370944 7.872825e-228 754.422724 |
| 86 | * 1024 6.633188e-45 146.757062 1.882404e-304 1008.953565 |
| 87 | * |
| 88 | * As you can see, there's quite a bit of slop between these estimates. |
| 89 | * In fact, the density of pseudoprimes is conjectured to be closer to the |
| 90 | * square of that upper bound. E.g. the density of pseudoprimes of size |
| 91 | * 256 is around 3 * 10^-27. The density of primes is very high, from |
| 92 | * 0.005636 at 256 bits to 0.001409 at 1024 bits, i.e. more than 10^-3. |
| 93 | * |
| 94 | * For those people used to cryptographic levels of security where the |
| 95 | * 56 bits of DES key space is too small because it's exhaustible with |
| 96 | * custom hardware searching engines, note that you are not generating |
| 97 | * 50,000,000 primes per second on each of 56,000 custom hardware chips |
| 98 | * for several hours. The chances that another Dinosaur Killer asteroid |
| 99 | * will land today is about 10^-11 or 2^-36, so it would be better to |
| 100 | * spend your time worrying about *that*. Well, okay, there should be |
| 101 | * some derating for the chance that astronomers haven't seen it yet, |
| 102 | * but I think you get the idea. For a good feel about the probability |
| 103 | * of various events, I have heard that a good book is by E'mile Borel, |
| 104 | * "Les Probabilite's et la vie". (The 's are accents, not apostrophes.) |
| 105 | * |
| 106 | * For more on the subject, try "Finding Four Million Large Random Primes", |
| 107 | * by Ronald Rivest, in Advancess in Cryptology: Proceedings of Crypto |
| 108 | * '90. He used a small-divisor test, then a Fermat test to the base 2, |
| 109 | * and then 8 iterations of a Miller-Rabin test. About 718 million random |
| 110 | * 256-bit integers were generated, 43,741,404 passed the small divisor |
| 111 | * test, 4,058,000 passed the Fermat test, and all 4,058,000 passed all |
| 112 | * 8 iterations of the Miller-Rabin test, proving their primality beyond |
| 113 | * most reasonable doubts. |
| 114 | * |
| 115 | * If the probability of getting a pseudoprime is some small p, then the |
| 116 | * probability of not getting it in t trials is (1-p)^t. Remember that, |
| 117 | * for small p, (1-p)^(1/p) ~ 1/e, the base of natural logarithms. |
| 118 | * (This is more commonly expressed as e = lim_{x\to\infty} (1+1/x)^x.) |
| 119 | * Thus, (1-p)^t ~ e^(-p*t) = exp(-p*t). So the odds of being able to |
| 120 | * do this many tests without seeing a pseudoprime if you assume that |
| 121 | * p = 10^-6 (one in a million) is one in 57.86. If you assume that |
| 122 | * p = 2*10^-6, it's one in 3347.6. So it's implausible that the density |
| 123 | * of pseudoprimes is much more than one millionth the density of primes. |
| 124 | * |
| 125 | * He also gives a theoretical argument that the chance of finding a |
| 126 | * 256-bit non-prime which satisfies one Fermat test to the base 2 is |
| 127 | * less than 10^-22. The small divisor test improves this number, and |
| 128 | * if the numbers are 512 bits (as needed for a 1024-bit key) the odds |
| 129 | * of failure shrink to about 10^-44. Thus, he concludes, for practical |
| 130 | * purposes *one* Fermat test to the base 2 is sufficient. |
| 131 | */ |
| 132 | static int |
| 133 | primeTest(struct BigNum const *bn, struct BigNum *e, struct BigNum *a, |
| 134 | int (*f)(void *arg, int c), void *arg) |
| 135 | { |
| 136 | unsigned i, j; |
| 137 | unsigned k, l; |
| 138 | int err; |
| 139 | |
| 140 | #if BNDEBUG /* Debugging */ |
| 141 | /* |
| 142 | * This is debugging code to test the sieving stage. |
| 143 | * If the sieving is wrong, it will let past numbers with |
| 144 | * small divisors. The prime test here will still work, and |
| 145 | * weed them out, but you'll be doing a lot more slow tests, |
| 146 | * and presumably excluding from consideration some other numbers |
| 147 | * which might be prime. This check just verifies that none |
| 148 | * of the candidates have any small divisors. If this |
| 149 | * code is enabled and never triggers, you can feel quite |
| 150 | * confident that the sieving is doing its job. |
| 151 | */ |
| 152 | i = bnLSWord(bn); |
| 153 | if (!(i % 2)) printf("bn div by 2!"); |
| 154 | i = bnModQ(bn, 51051); /* 51051 = 3 * 7 * 11 * 13 * 17 */ |
| 155 | if (!(i % 3)) printf("bn div by 3!"); |
| 156 | if (!(i % 7)) printf("bn div by 7!"); |
| 157 | if (!(i % 11)) printf("bn div by 11!"); |
| 158 | if (!(i % 13)) printf("bn div by 13!"); |
| 159 | if (!(i % 17)) printf("bn div by 17!"); |
| 160 | i = bnModQ(bn, 63365); /* 63365 = 5 * 19 * 23 * 29 */ |
| 161 | if (!(i % 5)) printf("bn div by 5!"); |
| 162 | if (!(i % 19)) printf("bn div by 19!"); |
| 163 | if (!(i % 23)) printf("bn div by 23!"); |
| 164 | if (!(i % 29)) printf("bn div by 29!"); |
| 165 | i = bnModQ(bn, 47027); /* 47027 = 31 * 37 * 41 */ |
| 166 | if (!(i % 31)) printf("bn div by 31!"); |
| 167 | if (!(i % 37)) printf("bn div by 37!"); |
| 168 | if (!(i % 41)) printf("bn div by 41!"); |
| 169 | #endif |
| 170 | |
| 171 | /* |
| 172 | * Now, check that bn is prime. If it passes to the base 2, |
| 173 | * it's prime beyond all reasonable doubt, and everything else |
| 174 | * is just gravy, but it gives people warm fuzzies to do it. |
| 175 | * |
| 176 | * This starts with verifying Euler's criterion for a base of 2. |
| 177 | * This is the fastest pseudoprimality test that I know of, |
| 178 | * saving a modular squaring over a Fermat test, as well as |
| 179 | * being stronger. 7/8 of the time, it's as strong as a strong |
| 180 | * pseudoprimality test, too. (The exception being when bn == |
| 181 | * 1 mod 8 and 2 is a quartic residue, i.e. bn is of the form |
| 182 | * a^2 + (8*b)^2.) The precise series of tricks used here is |
| 183 | * not documented anywhere, so here's an explanation. |
| 184 | * Euler's criterion states that if p is prime then a^((p-1)/2) |
| 185 | * is congruent to Jacobi(a,p), modulo p. Jacobi(a,p) is |
| 186 | * a function which is +1 if a is a square modulo p, and -1 if |
| 187 | * it is not. For a = 2, this is particularly simple. It's |
| 188 | * +1 if p == +/-1 (mod 8), and -1 if m == +/-3 (mod 8). |
| 189 | * If p == 3 mod 4, then all a strong test does is compute |
| 190 | * 2^((p-1)/2). and see if it's +1 or -1. (Euler's criterion |
| 191 | * says *which* it should be.) If p == 5 (mod 8), then |
| 192 | * 2^((p-1)/2) is -1, so the initial step in a strong test, |
| 193 | * looking at 2^((p-1)/4), is wasted - you're not going to |
| 194 | * find a +/-1 before then if it *is* prime, and it shouldn't |
| 195 | * have either of those values if it isn't. So don't bother. |
| 196 | * |
| 197 | * The remaining case is p == 1 (mod 8). In this case, we |
| 198 | * expect 2^((p-1)/2) == 1 (mod p), so we expect that the |
| 199 | * square root of this, 2^((p-1)/4), will be +/-1 (mod p). |
| 200 | * Evaluating this saves us a modular squaring 1/4 of the time. |
| 201 | * If it's -1, a strong pseudoprimality test would call p |
| 202 | * prime as well. Only if the result is +1, indicating that |
| 203 | * 2 is not only a quadratic residue, but a quartic one as well, |
| 204 | * does a strong pseudoprimality test verify more things than |
| 205 | * this test does. Good enough. |
| 206 | * |
| 207 | * We could back that down another step, looking at 2^((p-1)/8) |
| 208 | * if there was a cheap way to determine if 2 were expected to |
| 209 | * be a quartic residue or not. Dirichlet proved that 2 is |
| 210 | * a quartic residue iff p is of the form a^2 + (8*b^2). |
| 211 | * All primes == 1 (mod 4) can be expressed as a^2 + (2*b)^2, |
| 212 | * but I see no cheap way to evaluate this condition. |
| 213 | */ |
| 214 | if (bnCopy(e, bn) < 0) |
| 215 | return -1; |
| 216 | (void)bnSubQ(e, 1); |
| 217 | l = bnLSWord(e); |
| 218 | |
| 219 | j = 1; /* Where to start in prime array for strong prime tests */ |
| 220 | |
| 221 | if (l & 7) { |
| 222 | bnRShift(e, 1); |
| 223 | if (bnTwoExpMod(a, e, bn) < 0) |
| 224 | return -1; |
| 225 | if ((l & 7) == 6) { |
| 226 | /* bn == 7 mod 8, expect +1 */ |
| 227 | if (bnBits(a) != 1) |
| 228 | return 1; /* Not prime */ |
| 229 | k = 1; |
| 230 | } else { |
| 231 | /* bn == 3 or 5 mod 8, expect -1 == bn-1 */ |
| 232 | if (bnAddQ(a, 1) < 0) |
| 233 | return -1; |
| 234 | if (bnCmp(a, bn) != 0) |
| 235 | return 1; /* Not prime */ |
| 236 | k = 1; |
| 237 | if (l & 4) { |
| 238 | /* bn == 5 mod 8, make odd for strong tests */ |
| 239 | bnRShift(e, 1); |
| 240 | k = 2; |
| 241 | } |
| 242 | } |
| 243 | } else { |
| 244 | /* bn == 1 mod 8, expect 2^((bn-1)/4) == +/-1 mod bn */ |
| 245 | bnRShift(e, 2); |
| 246 | if (bnTwoExpMod(a, e, bn) < 0) |
| 247 | return -1; |
| 248 | if (bnBits(a) == 1) { |
| 249 | j = 0; /* Re-do strong prime test to base 2 */ |
| 250 | } else { |
| 251 | if (bnAddQ(a, 1) < 0) |
| 252 | return -1; |
| 253 | if (bnCmp(a, bn) != 0) |
| 254 | return 1; /* Not prime */ |
| 255 | } |
| 256 | k = 2 + bnMakeOdd(e); |
| 257 | } |
| 258 | /* It's prime! Now go on to confirmation tests */ |
| 259 | |
| 260 | /* |
| 261 | * Now, e = (bn-1)/2^k is odd. k >= 1, and has a given value |
| 262 | * with probability 2^-k, so its expected value is 2. |
| 263 | * j = 1 in the usual case when the previous test was as good as |
| 264 | * a strong prime test, but 1/8 of the time, j = 0 because |
| 265 | * the strong prime test to the base 2 needs to be re-done. |
| 266 | */ |
| 267 | for (i = j; i < CONFIRMTESTS; i++) { |
| 268 | if (f && (err = f(arg, '*')) < 0) |
| 269 | return err; |
| 270 | (void)bnSetQ(a, confirm[i]); |
| 271 | if (bnExpMod(a, a, e, bn) < 0) |
| 272 | return -1; |
| 273 | if (bnBits(a) == 1) |
| 274 | continue; /* Passed this test */ |
| 275 | |
| 276 | l = k; |
| 277 | for (;;) { |
| 278 | if (bnAddQ(a, 1) < 0) |
| 279 | return -1; |
| 280 | if (bnCmp(a, bn) == 0) /* Was result bn-1? */ |
| 281 | break; /* Prime */ |
| 282 | if (!--l) /* Reached end, not -1? luck? */ |
| 283 | return i+2-j; /* Failed, not prime */ |
| 284 | /* This portion is executed, on average, once. */ |
| 285 | (void)bnSubQ(a, 1); /* Put a back where it was. */ |
| 286 | if (bnSquare(a, a) < 0 || bnMod(a, a, bn) < 0) |
| 287 | return -1; |
| 288 | if (bnBits(a) == 1) |
| 289 | return i+2-j; /* Failed, not prime */ |
| 290 | } |
| 291 | /* It worked (to the base confirm[i]) */ |
| 292 | } |
| 293 | |
| 294 | /* Yes, we've decided that it's prime. */ |
| 295 | if (f && (err = f(arg, '*')) < 0) |
| 296 | return err; |
| 297 | return 0; /* Prime! */ |
| 298 | } |
| 299 | |
| 300 | /* |
| 301 | * Add x*y to bn, which is usually (but not always) < 65536. |
| 302 | * Do it in a simple linear manner. |
| 303 | */ |
| 304 | static int |
| 305 | bnAddMult(struct BigNum *bn, unsigned x, unsigned y) |
| 306 | { |
| 307 | unsigned long z = (unsigned long)x * y; |
| 308 | |
| 309 | while (z > 65535) { |
| 310 | if (bnAddQ(bn, 65535) < 0) |
| 311 | return -1; |
| 312 | z -= 65535; |
| 313 | } |
| 314 | return bnAddQ(bn, (unsigned)z); |
| 315 | } |
| 316 | |
| 317 | static int |
| 318 | bnSubMult(struct BigNum *bn, unsigned x, unsigned y) |
| 319 | { |
| 320 | unsigned long z = (unsigned long)x * y; |
| 321 | |
| 322 | while (z > 65535) { |
| 323 | if (bnSubQ(bn, 65535) < 0) |
| 324 | return -1; |
| 325 | z -= 65535; |
| 326 | } |
| 327 | return bnSubQ(bn, (unsigned)z); |
| 328 | } |
| 329 | |
| 330 | /* |
| 331 | * Modifies the bignum to return a nearby (slightly larger) number which |
| 332 | * is a probable prime. Returns >=0 on success or -1 on failure (out of |
| 333 | * memory). The return value is the number of unsuccessful modular |
| 334 | * exponentiations performed. This never gives up searching. |
| 335 | * |
| 336 | * All other arguments are optional. They may be NULL. They are: |
| 337 | * |
| 338 | * unsigned (*rand)(unsigned limit) |
| 339 | * For better distributed numbers, supply a non-null pointer to a |
| 340 | * function which returns a random x, 0 <= x < limit. (It may make it |
| 341 | * simpler to know that 0 < limit <= SHUFFLE, so you need at most a byte.) |
| 342 | * The program generates a large window of sieve data and then does |
| 343 | * pseudoprimality tests on the data. If a rand function is supplied, |
| 344 | * the candidates which survive sieving are shuffled with a window of |
| 345 | * size SHUFFLE before testing to increase the uniformity of the prime |
| 346 | * selection. This isn't perfect, but it reduces the correlation between |
| 347 | * the size of the prime-free gap before a prime and the probability |
| 348 | * that that prime will be found by a sequential search. |
| 349 | * |
| 350 | * If rand is NULL, sequential search is used. If you want sequential |
| 351 | * search, note that the search begins with the given number; if you're |
| 352 | * trying to generate consecutive primes, you must increment the previous |
| 353 | * one by two before calling this again. |
| 354 | * |
| 355 | * int (*f)(void *arg, int c), void *arg |
| 356 | * The function f argument, if non-NULL, is called with progress indicator |
| 357 | * characters for printing. A dot (.) is written every time a primality test |
| 358 | * is failed, a star (*) every time one is passed, and a slash (/) in the |
| 359 | * (very rare) case that the sieve was emptied without finding a prime |
| 360 | * and is being refilled. f is also passed the void *arg argument for |
| 361 | * private context storage. If f returns < 0, the test aborts and returns |
| 362 | * that value immediately. (bn is set to the last value tested, so you |
| 363 | * can increment bn and continue.) |
| 364 | * |
| 365 | * The "exponent" argument, and following unsigned numbers, are exponents |
| 366 | * for which an inverse is desired, modulo p. For a d to exist such that |
| 367 | * (x^e)^d == x (mod p), then d*e == 1 (mod p-1), so gcd(e,p-1) must be 1. |
| 368 | * The prime returned is constrained to not be congruent to 1 modulo |
| 369 | * any of the zero-terminated list of 16-bit numbers. Note that this list |
| 370 | * should contain all the small prime factors of e. (You'll have to test |
| 371 | * for large prime factors of e elsewhere, but the chances of needing to |
| 372 | * generate another prime are low.) |
| 373 | * |
| 374 | * The list is terminated by a 0, and may be empty. |
| 375 | */ |
| 376 | int |
| 377 | primeGen(struct BigNum *bn, unsigned (*rand)(unsigned), |
| 378 | int (*f)(void *arg, int c), void *arg, unsigned exponent, ...) |
| 379 | { |
| 380 | int retval; |
| 381 | int modexps = 0; |
| 382 | unsigned short offsets[SHUFFLE]; |
| 383 | unsigned i, j; |
| 384 | unsigned p, q, prev; |
| 385 | struct BigNum a, e; |
| 386 | #ifdef MSDOS |
| 387 | unsigned char *sieve; |
| 388 | #else |
| 389 | unsigned char sieve[SIEVE]; |
| 390 | #endif |
| 391 | |
| 392 | #ifdef MSDOS |
| 393 | sieve = lbnMemAlloc(SIEVE); |
| 394 | if (!sieve) |
| 395 | return -1; |
| 396 | #endif |
| 397 | |
| 398 | bnBegin(&a); |
| 399 | bnBegin(&e); |
| 400 | |
| 401 | #if 0 /* Self-test (not used for production) */ |
| 402 | { |
| 403 | struct BigNum t; |
| 404 | static unsigned char const prime1[] = {5}; |
| 405 | static unsigned char const prime2[] = {7}; |
| 406 | static unsigned char const prime3[] = {11}; |
| 407 | static unsigned char const prime4[] = {1, 1}; /* 257 */ |
| 408 | static unsigned char const prime5[] = {0xFF, 0xF1}; /* 65521 */ |
| 409 | static unsigned char const prime6[] = {1, 0, 1}; /* 65537 */ |
| 410 | static unsigned char const prime7[] = {1, 0, 3}; /* 65539 */ |
| 411 | /* A small prime: 1234567891 */ |
| 412 | static unsigned char const prime8[] = {0x49, 0x96, 0x02, 0xD3}; |
| 413 | /* A slightly larger prime: 12345678901234567891 */ |
| 414 | static unsigned char const prime9[] = { |
| 415 | 0xAB, 0x54, 0xA9, 0x8C, 0xEB, 0x1F, 0x0A, 0xD3 }; |
| 416 | /* |
| 417 | * No, 123456789012345678901234567891 isn't prime; it's just a |
| 418 | * lucky, easy-to-remember conicidence. (You have to go to |
| 419 | * ...4567907 for a prime.) |
| 420 | */ |
| 421 | static struct { |
| 422 | unsigned char const *prime; |
| 423 | unsigned size; |
| 424 | } const primelist[] = { |
| 425 | { prime1, sizeof(prime1) }, |
| 426 | { prime2, sizeof(prime2) }, |
| 427 | { prime3, sizeof(prime3) }, |
| 428 | { prime4, sizeof(prime4) }, |
| 429 | { prime5, sizeof(prime5) }, |
| 430 | { prime6, sizeof(prime6) }, |
| 431 | { prime7, sizeof(prime7) }, |
| 432 | { prime8, sizeof(prime8) }, |
| 433 | { prime9, sizeof(prime9) } }; |
| 434 | |
| 435 | bnBegin(&t); |
| 436 | |
| 437 | for (i = 0; i < sizeof(primelist)/sizeof(primelist[0]); i++) { |
| 438 | bnInsertBytes(&t, primelist[i].prime, 0, |
| 439 | primelist[i].size); |
| 440 | bnCopy(&e, &t); |
| 441 | (void)bnSubQ(&e, 1); |
| 442 | bnTwoExpMod(&a, &e, &t); |
| 443 | p = bnBits(&a); |
| 444 | if (p != 1) { |
| 445 | printf( |
| 446 | "Bug: Fermat(2) %u-bit output (1 expected)\n", p); |
| 447 | fputs("Prime = 0x", stdout); |
| 448 | for (j = 0; j < primelist[i].size; j++) |
| 449 | printf("%02X", primelist[i].prime[j]); |
| 450 | putchar('\n'); |
| 451 | } |
| 452 | bnSetQ(&a, 3); |
| 453 | bnExpMod(&a, &a, &e, &t); |
| 454 | p = bnBits(&a); |
| 455 | if (p != 1) { |
| 456 | printf( |
| 457 | "Bug: Fermat(3) %u-bit output (1 expected)\n", p); |
| 458 | fputs("Prime = 0x", stdout); |
| 459 | for (j = 0; j < primelist[i].size; j++) |
| 460 | printf("%02X", primelist[i].prime[j]); |
| 461 | putchar('\n'); |
| 462 | } |
| 463 | } |
| 464 | |
| 465 | bnEnd(&t); |
| 466 | } |
| 467 | #endif |
| 468 | |
| 469 | /* First, make sure that bn is odd. */ |
| 470 | if ((bnLSWord(bn) & 1) == 0) |
| 471 | (void)bnAddQ(bn, 1); |
| 472 | |
| 473 | retry: |
| 474 | /* Then build a sieve starting at bn. */ |
| 475 | sieveBuild(sieve, SIEVE, bn, 2, 0); |
| 476 | |
| 477 | /* Do the extra exponent sieving */ |
| 478 | if (exponent) { |
| 479 | va_list ap; |
| 480 | unsigned t = exponent; |
| 481 | |
| 482 | va_start(ap, exponent); |
| 483 | |
| 484 | do { |
| 485 | /* The exponent had better be odd! */ |
| 486 | assert(t & 1); |
| 487 | |
| 488 | i = bnModQ(bn, t); |
| 489 | /* Find 1-i */ |
| 490 | if (i == 0) |
| 491 | i = 1; |
| 492 | else if (--i) |
| 493 | i = t - i; |
| 494 | |
| 495 | /* Divide by 2, modulo the exponent */ |
| 496 | i = (i & 1) ? i/2 + t/2 + 1 : i/2; |
| 497 | |
| 498 | /* Remove all following multiples from the sieve. */ |
| 499 | sieveSingle(sieve, SIEVE, i, t); |
| 500 | |
| 501 | /* Get the next exponent value */ |
| 502 | t = va_arg(ap, unsigned); |
| 503 | } while (t); |
| 504 | |
| 505 | va_end(ap); |
| 506 | } |
| 507 | |
| 508 | /* Fill up the offsets array with the first SHUFFLE candidates */ |
| 509 | i = p = 0; |
| 510 | /* Get first prime */ |
| 511 | if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) { |
| 512 | offsets[i++] = p; |
| 513 | p = sieveSearch(sieve, SIEVE, p); |
| 514 | } |
| 515 | /* |
| 516 | * Okay, from this point onwards, p is always the next entry |
| 517 | * from the sieve, that has not been added to the shuffle table, |
| 518 | * and is 0 iff the sieve has been exhausted. |
| 519 | * |
| 520 | * If we want to shuffle, then fill the shuffle table until the |
| 521 | * sieve is exhausted or the table is full. |
| 522 | */ |
| 523 | if (rand && p) { |
| 524 | do { |
| 525 | offsets[i++] = p; |
| 526 | p = sieveSearch(sieve, SIEVE, p); |
| 527 | } while (p && i < SHUFFLE); |
| 528 | } |
| 529 | |
| 530 | /* Choose a random candidate for experimentation */ |
| 531 | prev = 0; |
| 532 | while (i) { |
| 533 | /* Pick a random entry from the shuffle table */ |
| 534 | j = rand ? rand(i) : 0; |
| 535 | q = offsets[j]; /* The entry to use */ |
| 536 | |
| 537 | /* Replace the entry with some more data, if possible */ |
| 538 | if (p) { |
| 539 | offsets[j] = p; |
| 540 | p = sieveSearch(sieve, SIEVE, p); |
| 541 | } else { |
| 542 | offsets[j] = offsets[--i]; |
| 543 | offsets[i] = 0; |
| 544 | } |
| 545 | |
| 546 | /* Adjust bn to have the right value */ |
| 547 | if ((q > prev ? bnAddMult(bn, q-prev, 2) |
| 548 | : bnSubMult(bn, prev-q, 2)) < 0) |
| 549 | goto failed; |
| 550 | prev = q; |
| 551 | |
| 552 | /* Now do the Fermat tests */ |
| 553 | retval = primeTest(bn, &e, &a, f, arg); |
| 554 | if (retval <= 0) |
| 555 | goto done; /* Success or error */ |
| 556 | modexps += retval; |
| 557 | if (f && (retval = f(arg, '.')) < 0) |
| 558 | goto done; |
| 559 | } |
| 560 | |
| 561 | /* Ran out of sieve space - increase bn and keep trying. */ |
| 562 | if (bnAddMult(bn, SIEVE*8-prev, 2) < 0) |
| 563 | goto failed; |
| 564 | if (f && (retval = f(arg, '/')) < 0) |
| 565 | goto done; |
| 566 | goto retry; |
| 567 | |
| 568 | failed: |
| 569 | retval = -1; |
| 570 | done: |
| 571 | bnEnd(&e); |
| 572 | bnEnd(&a); |
| 573 | lbnMemWipe(offsets, sizeof(offsets)); |
| 574 | #ifdef MSDOS |
| 575 | lbnMemFree(sieve, SIEVE); |
| 576 | #else |
| 577 | lbnMemWipe(sieve, sizeof(sieve)); |
| 578 | #endif |
| 579 | |
| 580 | return retval < 0 ? retval : modexps + CONFIRMTESTS; |
| 581 | } |
| 582 | |
| 583 | /* |
| 584 | * Similar, but searches forward from the given starting value in steps of |
| 585 | * "step" rather than 1. The step size must be even, and bn must be odd. |
| 586 | * Among other possibilities, this can be used to generate "strong" |
| 587 | * primes, where p-1 has a large prime factor. |
| 588 | */ |
| 589 | int |
| 590 | primeGenStrong(struct BigNum *bn, struct BigNum const *step, |
| 591 | int (*f)(void *arg, int c), void *arg) |
| 592 | { |
| 593 | int retval; |
| 594 | unsigned p, prev; |
| 595 | struct BigNum a, e; |
| 596 | int modexps = 0; |
| 597 | #ifdef MSDOS |
| 598 | unsigned char *sieve; |
| 599 | #else |
| 600 | unsigned char sieve[SIEVE]; |
| 601 | #endif |
| 602 | |
| 603 | #ifdef MSDOS |
| 604 | sieve = lbnMemAlloc(SIEVE); |
| 605 | if (!sieve) |
| 606 | return -1; |
| 607 | #endif |
| 608 | |
| 609 | /* Step must be even and bn must be odd */ |
| 610 | assert((bnLSWord(step) & 1) == 0); |
| 611 | assert((bnLSWord(bn) & 1) == 1); |
| 612 | |
| 613 | bnBegin(&a); |
| 614 | bnBegin(&e); |
| 615 | |
| 616 | for (;;) { |
| 617 | if (sieveBuildBig(sieve, SIEVE, bn, step, 0) < 0) |
| 618 | goto failed; |
| 619 | |
| 620 | p = prev = 0; |
| 621 | if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) { |
| 622 | do { |
| 623 | /* |
| 624 | * Adjust bn to have the right value, |
| 625 | * adding (p-prev) * 2*step. |
| 626 | */ |
| 627 | assert(p >= prev); |
| 628 | /* Compute delta into a */ |
| 629 | if (bnMulQ(&a, step, p-prev) < 0) |
| 630 | goto failed; |
| 631 | if (bnAdd(bn, &a) < 0) |
| 632 | goto failed; |
| 633 | prev = p; |
| 634 | |
| 635 | retval = primeTest(bn, &e, &a, f, arg); |
| 636 | if (retval <= 0) |
| 637 | goto done; /* Success! */ |
| 638 | modexps += retval; |
| 639 | if (f && (retval = f(arg, '.')) < 0) |
| 640 | goto done; |
| 641 | |
| 642 | /* And try again */ |
| 643 | p = sieveSearch(sieve, SIEVE, p); |
| 644 | } while (p); |
| 645 | } |
| 646 | |
| 647 | /* Ran out of sieve space - increase bn and keep trying. */ |
| 648 | #if SIEVE*8 == 65536 |
| 649 | /* Corner case that will never actually happen */ |
| 650 | if (!prev) { |
| 651 | if (bnAdd(bn, step) < 0) |
| 652 | goto failed; |
| 653 | p = 65535; |
| 654 | } else { |
| 655 | p = (unsigned)(SIEVE*8 - prev); |
| 656 | } |
| 657 | #else |
| 658 | p = SIEVE*8 - prev; |
| 659 | #endif |
| 660 | if (bnMulQ(&a, step, p) < 0 || bnAdd(bn, &a) < 0) |
| 661 | goto failed; |
| 662 | if (f && (retval = f(arg, '/')) < 0) |
| 663 | goto done; |
| 664 | } /* for (;;) */ |
| 665 | |
| 666 | failed: |
| 667 | retval = -1; |
| 668 | |
| 669 | done: |
| 670 | |
| 671 | bnEnd(&e); |
| 672 | bnEnd(&a); |
| 673 | #ifdef MSDOS |
| 674 | lbnMemFree(sieve, SIEVE); |
| 675 | #else |
| 676 | lbnMemWipe(sieve, sizeof(sieve)); |
| 677 | #endif |
| 678 | return retval < 0 ? retval : modexps + CONFIRMTESTS; |
| 679 | } |