Alexandre Lision | 7fd5d3d | 2013-12-04 13:06:40 -0500 | [diff] [blame] | 1 | /* |
| 2 | * Sophie Germain 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 | #define BNDEBUG 1 |
| 28 | #ifndef BNDEBUG |
| 29 | #define BNDEBUG 0 |
| 30 | #endif |
| 31 | #if BNDEBUG |
| 32 | #include <stdio.h> |
| 33 | #endif |
| 34 | |
| 35 | #include "bn.h" |
| 36 | #include "germain.h" |
| 37 | #include "jacobi.h" |
| 38 | #include "lbnmem.h" /* For lbnMemWipe */ |
| 39 | #include "sieve.h" |
| 40 | |
| 41 | #include "kludge.h" |
| 42 | |
| 43 | /* Size of the sieve area (can be up to 65536/8 = 8192) */ |
| 44 | #define SIEVE 8192 |
| 45 | |
| 46 | static unsigned const confirm[] = {2, 3, 5, 7, 11, 13, 17}; |
| 47 | #define CONFIRMTESTS (sizeof(confirm)/sizeof(*confirm)) |
| 48 | |
| 49 | #if BNDEBUG |
| 50 | /* |
| 51 | * For sanity checking the sieve, we check for small divisors of the numbers |
| 52 | * we get back. This takes "rem", a partially reduced form of the prime, |
| 53 | * "div" a divisor to check for, and "order", a parameter of the "order" |
| 54 | * of Sophie Germain primes (0 = normal primes, 1 = Sophie Germain primes, |
| 55 | * 2 = 4*p+3 is also prime, etc.) and does the check. It just complains |
| 56 | * to stdout if the check fails. |
| 57 | */ |
| 58 | static void |
| 59 | germainSanity(unsigned rem, unsigned div, unsigned order) |
| 60 | { |
| 61 | unsigned mul = 1; |
| 62 | |
| 63 | rem %= div; |
| 64 | if (!rem) |
| 65 | printf("bn div by %u!\n", div); |
| 66 | while (order--) { |
| 67 | rem += rem+1; |
| 68 | if (rem >= div) |
| 69 | rem -= div; |
| 70 | mul += mul; |
| 71 | if (!rem) |
| 72 | printf("%u*bn+%u div by %u!\n", mul, mul-1, div); |
| 73 | } |
| 74 | } |
| 75 | #endif /* BNDEBUG */ |
| 76 | |
| 77 | /* |
| 78 | * Helper function that does the slow primality test. |
| 79 | * bn is the input bignum; a, e and bn2 are temporary buffers that are |
| 80 | * allocated by the caller to save overhead. bn2 is filled with |
| 81 | * a copy of 2^order*bn+2^order-1 if bn is found to be prime. |
| 82 | * |
| 83 | * Returns 0 if both bn and bn2 are prime, >0 if not prime, and -1 on |
| 84 | * error (out of memory). If not prime, the return value is the number |
| 85 | * of modular exponentiations performed. Prints a '+' or '-' on the |
| 86 | * given FILE (if any) for each test that is passed by bn, and a '*' |
| 87 | * for each test that is passed by bn2. |
| 88 | * |
| 89 | * The testing consists of strong pseudoprimality tests, to the bases given |
| 90 | * in the confirm[] array above. (Also called Miller-Rabin, although that's |
| 91 | * not technically correct if we're using fixed bases.) Some people worry |
| 92 | * that this might not be enough. Number theorists may wish to generate |
| 93 | * primality proofs, but for random inputs, this returns non-primes with |
| 94 | * a probability which is quite negligible, which is good enough. |
| 95 | * |
| 96 | * It has been proved (see Carl Pomerance, "On the Distribution of |
| 97 | * Pseudoprimes", Math. Comp. v.37 (1981) pp. 587-593) that the number of |
| 98 | * pseudoprimes (composite numbers that pass a Fermat test to the base 2) |
| 99 | * less than x is bounded by: |
| 100 | * exp(ln(x)^(5/14)) <= P_2(x) ### CHECK THIS FORMULA - it looks wrong! ### |
| 101 | * P_2(x) <= x * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))). |
| 102 | * Thus, the local density of Pseudoprimes near x is at most |
| 103 | * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))), and at least |
| 104 | * exp(ln(x)^(5/14) - ln(x)). Here are some values of this function |
| 105 | * for various k-bit numbers x = 2^k: |
| 106 | * Bits Density <= Bit equivalent Density >= Bit equivalent |
| 107 | * 128 3.577869e-07 21.414396 4.202213e-37 120.840190 |
| 108 | * 192 4.175629e-10 31.157288 4.936250e-56 183.724558 |
| 109 | * 256 5.804314e-13 40.647940 4.977813e-75 246.829095 |
| 110 | * 384 1.578039e-18 59.136573 3.938861e-113 373.400096 |
| 111 | * 512 5.858255e-24 77.175803 2.563353e-151 500.253110 |
| 112 | * 768 1.489276e-34 112.370944 7.872825e-228 754.422724 |
| 113 | * 1024 6.633188e-45 146.757062 1.882404e-304 1008.953565 |
| 114 | * |
| 115 | * As you can see, there's quite a bit of slop between these estimates. |
| 116 | * In fact, the density of pseudoprimes is conjectured to be closer to the |
| 117 | * square of that upper bound. E.g. the density of pseudoprimes of size |
| 118 | * 256 is around 3 * 10^-27. The density of primes is very high, from |
| 119 | * 0.005636 at 256 bits to 0.001409 at 1024 bits, i.e. more than 10^-3. |
| 120 | * |
| 121 | * For those people used to cryptographic levels of security where the |
| 122 | * 56 bits of DES key space is too small because it's exhaustible with |
| 123 | * custom hardware searching engines, note that you are not generating |
| 124 | * 50,000,000 primes per second on each of 56,000 custom hardware chips |
| 125 | * for several hours. The chances that another Dinosaur Killer asteroid |
| 126 | * will land today is about 10^-11 or 2^-36, so it would be better to |
| 127 | * spend your time worrying about *that*. Well, okay, there should be |
| 128 | * some derating for the chance that astronomers haven't seen it yet, |
| 129 | * but I think you get the idea. For a good feel about the probability |
| 130 | * of various events, I have heard that a good book is by E'mile Borel, |
| 131 | * "Les Probabilite's et la vie". (The 's are accents, not apostrophes.) |
| 132 | * |
| 133 | * For more on the subject, try "Finding Four Million Large Random Primes", |
| 134 | * by Ronald Rivest, in Advancess in Cryptology: Proceedings of Crypto |
| 135 | * '90. He used a small-divisor test, then a Fermat test to the base 2, |
| 136 | * and then 8 iterations of a Miller-Rabin test. About 718 million random |
| 137 | * 256-bit integers were generated, 43,741,404 passed the small divisor |
| 138 | * test, 4,058,000 passed the Fermat test, and all 4,058,000 passed all |
| 139 | * 8 iterations of the Miller-Rabin test, proving their primality beyond |
| 140 | * most reasonable doubts. |
| 141 | * |
| 142 | * If the probability of getting a pseudoprime is some small p, then the |
| 143 | * probability of not getting it in t trials is (1-p)^t. Remember that, |
| 144 | * for small p, (1-p)^(1/p) ~ 1/e, the base of natural logarithms. |
| 145 | * (This is more commonly expressed as e = lim_{x\to\infty} (1+1/x)^x.) |
| 146 | * Thus, (1-p)^t ~ e^(-p*t) = exp(-p*t). So the odds of being able to |
| 147 | * do this many tests without seeing a pseudoprime if you assume that |
| 148 | * p = 10^-6 (one in a million) is one in 57.86. If you assume that |
| 149 | * p = 2*10^-6, it's one in 3347.6. So it's implausible that the density |
| 150 | * of pseudoprimes is much more than one millionth the density of primes. |
| 151 | * |
| 152 | * He also gives a theoretical argument that the chance of finding a |
| 153 | * 256-bit non-prime which satisfies one Fermat test to the base 2 is |
| 154 | * less than 10^-22. The small divisor test improves this number, and |
| 155 | * if the numbers are 512 bits (as needed for a 1024-bit key) the odds |
| 156 | * of failure shrink to about 10^-44. Thus, he concludes, for practical |
| 157 | * purposes *one* Fermat test to the base 2 is sufficient. |
| 158 | */ |
| 159 | static int |
| 160 | germainPrimeTest(struct BigNum const *bn, struct BigNum *bn2, struct BigNum *e, |
| 161 | struct BigNum *a, unsigned order, int (*f)(void *arg, int c), void *arg) |
| 162 | { |
| 163 | int err; |
| 164 | unsigned i; |
| 165 | int j; |
| 166 | unsigned k, l, n; |
| 167 | |
| 168 | #if BNDEBUG /* Debugging */ |
| 169 | /* |
| 170 | * This is debugging code to test the sieving stage. |
| 171 | * If the sieving is wrong, it will let past numbers with |
| 172 | * small divisors. The prime test here will still work, and |
| 173 | * weed them out, but you'll be doing a lot more slow tests, |
| 174 | * and presumably excluding from consideration some other numbers |
| 175 | * which might be prime. This check just verifies that none |
| 176 | * of the candidates have any small divisors. If this |
| 177 | * code is enabled and never triggers, you can feel quite |
| 178 | * confident that the sieving is doing its job. |
| 179 | */ |
| 180 | i = bnLSWord(bn); |
| 181 | if (!(i % 2)) printf("bn div by 2!"); |
| 182 | i = bnModQ(bn, 51051); /* 51051 = 3 * 7 * 11 * 13 * 17 */ |
| 183 | germainSanity(i, 3, order); |
| 184 | germainSanity(i, 7, order); |
| 185 | germainSanity(i, 11, order); |
| 186 | germainSanity(i, 13, order); |
| 187 | germainSanity(i, 17, order); |
| 188 | i = bnModQ(bn, 63365); /* 63365 = 5 * 19 * 23 * 29 */ |
| 189 | germainSanity(i, 5, order); |
| 190 | germainSanity(i, 19, order); |
| 191 | germainSanity(i, 23, order); |
| 192 | germainSanity(i, 29, order); |
| 193 | i = bnModQ(bn, 47027); /* 47027 = 31 * 37 * 41 */ |
| 194 | germainSanity(i, 31, order); |
| 195 | germainSanity(i, 37, order); |
| 196 | germainSanity(i, 41, order); |
| 197 | #endif |
| 198 | /* |
| 199 | * First, check whether bn is prime. This uses a fast primality |
| 200 | * test which usually obviates the need to do one of the |
| 201 | * confirmation tests later. See prime.c for a full explanation. |
| 202 | * We check bn first because it's one bit smaller, saving one |
| 203 | * modular squaring, and because we might be able to save another |
| 204 | * when testing it. (1/4 of the time.) A small speed hack, |
| 205 | * but finding big Sophie Germain primes is *slow*. |
| 206 | */ |
| 207 | if (bnCopy(e, bn) < 0) |
| 208 | return -1; |
| 209 | (void)bnSubQ(e, 1); |
| 210 | l = bnLSWord(e); |
| 211 | |
| 212 | j = 1; /* Where to start in prime array for strong prime tests */ |
| 213 | |
| 214 | if (l & 7) { |
| 215 | bnRShift(e, 1); |
| 216 | if (bnTwoExpMod(a, e, bn) < 0) |
| 217 | return -1; |
| 218 | if ((l & 7) == 6) { |
| 219 | /* bn == 7 mod 8, expect +1 */ |
| 220 | if (bnBits(a) != 1) |
| 221 | return 1; /* Not prime */ |
| 222 | k = 1; |
| 223 | } else { |
| 224 | /* bn == 3 or 5 mod 8, expect -1 == bn-1 */ |
| 225 | if (bnAddQ(a, 1) < 0) |
| 226 | return -1; |
| 227 | if (bnCmp(a, bn) != 0) |
| 228 | return 1; /* Not prime */ |
| 229 | k = 1; |
| 230 | if (l & 4) { |
| 231 | /* bn == 5 mod 8, make odd for strong tests */ |
| 232 | bnRShift(e, 1); |
| 233 | k = 2; |
| 234 | } |
| 235 | } |
| 236 | } else { |
| 237 | /* bn == 1 mod 8, expect 2^((bn-1)/4) == +/-1 mod bn */ |
| 238 | bnRShift(e, 2); |
| 239 | if (bnTwoExpMod(a, e, bn) < 0) |
| 240 | return -1; |
| 241 | if (bnBits(a) == 1) { |
| 242 | j = 0; /* Re-do strong prime test to base 2 */ |
| 243 | } else { |
| 244 | if (bnAddQ(a, 1) < 0) |
| 245 | return -1; |
| 246 | if (bnCmp(a, bn) != 0) |
| 247 | return 1; /* Not prime */ |
| 248 | } |
| 249 | k = 2 + bnMakeOdd(e); |
| 250 | } |
| 251 | |
| 252 | |
| 253 | /* |
| 254 | * It's prime! Now check higher-order forms bn2 = 2*bn+1, 4*bn+3, |
| 255 | * etc. Since bn2 == 3 mod 4, a strong pseudoprimality test boils |
| 256 | * down to looking at a^((bn2-1)/2) mod bn and seeing if it's +/-1. |
| 257 | * (+1 if bn2 is == 7 mod 8, -1 if it's == 3) |
| 258 | * Of course, that exponent is just the previous bn2 or bn... |
| 259 | */ |
| 260 | if (bnCopy(bn2, bn) < 0) |
| 261 | return -1; |
| 262 | for (n = 0; n < order; n++) { |
| 263 | /* |
| 264 | * Print a success indicator: the sign of Jacobi(2,bn2), |
| 265 | * which is available to us in l. bn2 = 2*bn + 1. Since bn |
| 266 | * is odd, bn2 must be == 3 mod 4, so the options modulo 8 |
| 267 | * are 3 and 7. 3 if l == 1 mod 4, 7 if l == 3 mod 4. |
| 268 | * The sign of the Jacobi symbol is - and + for these cases, |
| 269 | * respectively. |
| 270 | */ |
| 271 | if (f && (err = f(arg, "-+"[(l >> 1) & 1])) < 0) |
| 272 | return err; |
| 273 | /* Exponent is previous bn2 */ |
| 274 | if (bnCopy(e, bn2) < 0 || bnLShift(bn2, 1) < 0) |
| 275 | return -1; |
| 276 | (void)bnAddQ(bn2, 1); /* Can't overflow */ |
| 277 | if (bnTwoExpMod(a, e, bn2) < 0) |
| 278 | return -1; |
| 279 | if (n | l) { /* Expect + */ |
| 280 | if (bnBits(a) != 1) |
| 281 | return 2+n; /* Not prime */ |
| 282 | } else { |
| 283 | if (bnAddQ(a, 1) < 0) |
| 284 | return -1; |
| 285 | if (bnCmp(a, bn2) != 0) |
| 286 | return 2+n; /* Not prime */ |
| 287 | } |
| 288 | l = bnLSWord(bn2); |
| 289 | } |
| 290 | |
| 291 | /* Final success indicator - it's in the bag. */ |
| 292 | if (f && (err = f(arg, '*')) < 0) |
| 293 | return err; |
| 294 | |
| 295 | /* |
| 296 | * Success! We have found a prime! Now go on to confirmation |
| 297 | * tests... k is an amount by which we know it's safe to shift |
| 298 | * down e. j = 1 unless the test to the base 2 could stand to be |
| 299 | * re-done (it wasn't *quite* a strong test), in which case it's 0. |
| 300 | * |
| 301 | * Here, we do the full strong pseudoprimality test. This proves |
| 302 | * that a number is composite, or says that it's probably prime. |
| 303 | * |
| 304 | * For the given base a, find bn-1 = 2^k * e, then find |
| 305 | * x == a^e (mod bn). |
| 306 | * If x == +1 -> strong pseudoprime to base a |
| 307 | * Otherwise, repeat k times: |
| 308 | * If x == -1, -> strong pseudoprime |
| 309 | * x = x^2 (mod bn) |
| 310 | * If x = +1 -> composite |
| 311 | * If we reach the end of the iteration and x is *not* +1, at the |
| 312 | * end, it is composite. But it's also composite if the result |
| 313 | * *is* +1. Which means that the squaring actually only has to |
| 314 | * proceed k-1 times. If x is not -1 by then, it's composite |
| 315 | * no matter what the result of the squaring is. |
| 316 | * |
| 317 | * For the multiples 2*bn+1, 4*bn+3, etc. then k = 1 (and e is |
| 318 | * the previous multiple of bn) so the squaring loop is never |
| 319 | * actually executed at all. |
| 320 | */ |
| 321 | for (i = j; i < CONFIRMTESTS; i++) { |
| 322 | if (bnCopy(e, bn) < 0) |
| 323 | return -1; |
| 324 | bnRShift(e, k); |
| 325 | k += bnMakeOdd(e); |
| 326 | (void)bnSetQ(a, confirm[i]); |
| 327 | if (bnExpMod(a, a, e, bn) < 0) |
| 328 | return -1; |
| 329 | |
| 330 | if (bnBits(a) != 1) { |
| 331 | l = k; |
| 332 | for (;;) { |
| 333 | if (bnAddQ(a, 1) < 0) |
| 334 | return -1; |
| 335 | if (bnCmp(a, bn) == 0) /* Was result bn-1? */ |
| 336 | break; /* Prime */ |
| 337 | if (!--l) |
| 338 | return (1+order)*i+2; /* Fail */ |
| 339 | /* This part is executed once, on average. */ |
| 340 | (void)bnSubQ(a, 1); /* Restore a */ |
| 341 | if (bnSquare(a, a) < 0 || bnMod(a, a, bn) < 0) |
| 342 | return -1; |
| 343 | if (bnBits(a) == 1) |
| 344 | return (1+order)*i+1; /* Fail */ |
| 345 | } |
| 346 | } |
| 347 | |
| 348 | if (bnCopy(bn2, bn) < 0) |
| 349 | return -1; |
| 350 | |
| 351 | /* Only do the following if we're not re-doing base 2 */ |
| 352 | if (i) for (n = 0; n < order; n++) { |
| 353 | if (bnCopy(e, bn2) < 0 || bnLShift(bn2, 1) < 0) |
| 354 | return -1; |
| 355 | (void)bnAddQ(bn2, 1); |
| 356 | |
| 357 | /* Print success indicator for previous test */ |
| 358 | j = bnJacobiQ(confirm[i], bn2); |
| 359 | if (f && (err = f(arg, j < 0 ? '-' : '+')) < 0) |
| 360 | return err; |
| 361 | |
| 362 | /* Check that p^e == Jacobi(p,bn2) (mod bn2) */ |
| 363 | (void)bnSetQ(a, confirm[i]); |
| 364 | if (bnExpMod(a, a, e, bn2) < 0) |
| 365 | return -1; |
| 366 | /* |
| 367 | * FIXME: Actually, we don't need to compute the |
| 368 | * Jacobi symbol externally... it never happens that |
| 369 | * a = +/-1 but it's the wrong one. So we can just |
| 370 | * look at a and use its sign. Find a proof somewhere. |
| 371 | */ |
| 372 | if (j < 0) { |
| 373 | /* Not a Q.R., should have a = bn2-1 */ |
| 374 | if (bnAddQ(a, 1) < 0) |
| 375 | return -1; |
| 376 | if (bnCmp(a, bn2) != 0) /* Was result bn2-1? */ |
| 377 | return (1+order)*i+n+2; /* Fail */ |
| 378 | } else { |
| 379 | /* Quadratic residue, should have a = 1 */ |
| 380 | if (bnBits(a) != 1) |
| 381 | return (1+order)*i+n+2; /* Fail */ |
| 382 | } |
| 383 | } |
| 384 | /* Final success indicator for the base confirm[i]. */ |
| 385 | if (f && (err = f(arg, '*')) < 0) |
| 386 | return err; |
| 387 | } |
| 388 | |
| 389 | return 0; /* Prime! */ |
| 390 | } |
| 391 | |
| 392 | /* |
| 393 | * Add x*y to bn, which is usually (but not always) < 65536. |
| 394 | * Do it in a simple linear manner. |
| 395 | */ |
| 396 | static int |
| 397 | bnAddMult(struct BigNum *bn, unsigned long x, unsigned y) |
| 398 | { |
| 399 | unsigned long z = (unsigned long)x * y; |
| 400 | |
| 401 | while (z > 65535) { |
| 402 | if (bnAddQ(bn, 65535) < 0) |
| 403 | return -1; |
| 404 | z -= 65535; |
| 405 | } |
| 406 | return bnAddQ(bn, (unsigned)z); |
| 407 | } |
| 408 | |
| 409 | /* |
| 410 | * Modifies the bignum to return the next Sophie Germain prime >= the |
| 411 | * input value. Sohpie Germain primes are number such that p is |
| 412 | * prime and 2*p+1 is also prime. |
| 413 | * |
| 414 | * This is actually parameterized: it generates primes p such that "order" |
| 415 | * multiples-plus-two are also prime, 2*p+1, 2*(2*p+1)+1 = 4*p+3, etc. |
| 416 | * |
| 417 | * Returns >=0 on success or -1 on failure (out of memory). On success, |
| 418 | * the return value is the number of modular exponentiations performed |
| 419 | * (excluding the final confirmations). This never gives up searching. |
| 420 | * |
| 421 | * The FILE *f argument, if non-NULL, has progress indicators written |
| 422 | * to it. A dot (.) is written every time a primeality test is failed, |
| 423 | * a plus (+) or minus (-) when the smaller prime of the pair passes a |
| 424 | * test, and a star (*) when the larger one does. Finally, a slash (/) |
| 425 | * is printed when the sieve was emptied without finding a prime and is |
| 426 | * being refilled. |
| 427 | * |
| 428 | * Apologies to structured programmers for all the GOTOs. |
| 429 | */ |
| 430 | int |
| 431 | germainPrimeGen(struct BigNum *bn, unsigned order, |
| 432 | int (*f)(void *arg, int c), void *arg) |
| 433 | { |
| 434 | int retval; |
| 435 | unsigned p, prev; |
| 436 | unsigned inc; |
| 437 | struct BigNum a, e, bn2; |
| 438 | int modexps = 0; |
| 439 | #ifdef MSDOS |
| 440 | unsigned char *sieve; |
| 441 | #else |
| 442 | unsigned char sieve[SIEVE]; |
| 443 | #endif |
| 444 | |
| 445 | #ifdef MSDOS |
| 446 | sieve = lbnMemAlloc(SIEVE); |
| 447 | if (!sieve) |
| 448 | return -1; |
| 449 | #endif |
| 450 | |
| 451 | bnBegin(&a); |
| 452 | bnBegin(&e); |
| 453 | bnBegin(&bn2); |
| 454 | |
| 455 | /* |
| 456 | * Obviously, the prime we find must be odd. Further, if 2*p+1 |
| 457 | * is also to be prime (order > 0) then p != 1 (mod 3), lest |
| 458 | * 2*p+1 == 3 (mod 3). Added to p != 3 (mod 3), p == 2 (mod 3) |
| 459 | * and p == 5 (mod 6). |
| 460 | * If order > 2 and we care about 4*p+3 and 8*p+7, then similarly |
| 461 | * p == 4 (mod 5), so p == 29 (mod 30). |
| 462 | * So pick the step size for searching based on the order |
| 463 | * and increse bn until it's == -1 (mod inc). |
| 464 | * |
| 465 | * mod 7 doesn't have a unique value for p because 2 -> 5 -> 4 -> 2, |
| 466 | * nor does mod 11, and I don't want to think about things past |
| 467 | * that. The required order would be impractically high, in any case. |
| 468 | */ |
| 469 | inc = order ? ((order > 2) ? 30 : 6) : 2; |
| 470 | if (bnAddQ(bn, inc-1 - bnModQ(bn, inc)) < 0) |
| 471 | goto failed; |
| 472 | |
| 473 | for (;;) { |
| 474 | if (sieveBuild(sieve, SIEVE, bn, inc, order) < 0) |
| 475 | goto failed; |
| 476 | |
| 477 | p = prev = 0; |
| 478 | if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) { |
| 479 | do { |
| 480 | /* Adjust bn to have the right value. */ |
| 481 | assert(p >= prev); |
| 482 | if (bnAddMult(bn, p-prev, inc) < 0) |
| 483 | goto failed; |
| 484 | prev = p; |
| 485 | |
| 486 | /* Okay, do the strong tests. */ |
| 487 | retval = germainPrimeTest(bn, &bn2, &e, &a, |
| 488 | order, f, arg); |
| 489 | if (retval <= 0) |
| 490 | goto done; |
| 491 | modexps += retval; |
| 492 | if (f && (retval = f(arg, '.')) < 0) |
| 493 | goto done; |
| 494 | |
| 495 | /* And try again */ |
| 496 | p = sieveSearch(sieve, SIEVE, p); |
| 497 | } while (p); |
| 498 | } |
| 499 | |
| 500 | /* Ran out of sieve space - increase bn and keep trying. */ |
| 501 | if (bnAddMult(bn, (unsigned long)SIEVE*8-prev, inc) < 0) |
| 502 | goto failed; |
| 503 | if (f && (retval = f(arg, '/')) < 0) |
| 504 | goto done; |
| 505 | } /* for (;;) */ |
| 506 | |
| 507 | failed: |
| 508 | retval = -1; |
| 509 | done: |
| 510 | bnEnd(&bn2); |
| 511 | bnEnd(&e); |
| 512 | bnEnd(&a); |
| 513 | #ifdef MSDOS |
| 514 | lbnMemFree(sieve, SIEVE); |
| 515 | #else |
| 516 | lbnMemWipe(sieve, sizeof(sieve)); |
| 517 | #endif |
| 518 | return retval < 0 ? retval : modexps+(order+1)*CONFIRMTESTS; |
| 519 | } |
| 520 | |
| 521 | int |
| 522 | germainPrimeGenStrong(struct BigNum *bn, struct BigNum const *step, |
| 523 | unsigned order, int (*f)(void *arg, int c), void *arg) |
| 524 | { |
| 525 | int retval; |
| 526 | unsigned p, prev; |
| 527 | struct BigNum a, e, bn2; |
| 528 | int modexps = 0; |
| 529 | #ifdef MSDOS |
| 530 | unsigned char *sieve; |
| 531 | #else |
| 532 | unsigned char sieve[SIEVE]; |
| 533 | #endif |
| 534 | |
| 535 | #ifdef MSDOS |
| 536 | sieve = lbnMemAlloc(SIEVE); |
| 537 | if (!sieve) |
| 538 | return -1; |
| 539 | #endif |
| 540 | bnBegin(&a); |
| 541 | bnBegin(&e); |
| 542 | bnBegin(&bn2); |
| 543 | |
| 544 | for (;;) { |
| 545 | if (sieveBuildBig(sieve, SIEVE, bn, step, order) < 0) |
| 546 | goto failed; |
| 547 | |
| 548 | p = prev = 0; |
| 549 | if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) { |
| 550 | do { |
| 551 | /* |
| 552 | * Adjust bn to have the right value, |
| 553 | * adding (p-prev) * 2*step. |
| 554 | */ |
| 555 | assert(p >= prev); |
| 556 | /* Compute delta into a */ |
| 557 | if (bnMulQ(&a, step, p-prev) < 0) |
| 558 | goto failed; |
| 559 | if (bnAdd(bn, &a) < 0) |
| 560 | goto failed; |
| 561 | prev = p; |
| 562 | |
| 563 | /* Okay, do the strong tests. */ |
| 564 | retval = germainPrimeTest(bn, &bn2, &e, &a, |
| 565 | order, f, arg); |
| 566 | if (retval <= 0) |
| 567 | goto done; |
| 568 | modexps += retval; |
| 569 | if (f && (retval = f(arg, '.')) < 0) |
| 570 | goto done; |
| 571 | |
| 572 | /* And try again */ |
| 573 | p = sieveSearch(sieve, SIEVE, p); |
| 574 | } while (p); |
| 575 | } |
| 576 | |
| 577 | /* Ran out of sieve space - increase bn and keep trying. */ |
| 578 | #if SIEVE*8 == 65536 |
| 579 | /* Corner case that will never actually happen */ |
| 580 | if (!prev) { |
| 581 | if (bnAdd(bn, step) < 0) |
| 582 | goto failed; |
| 583 | p = 65535; |
| 584 | } else { |
| 585 | p = (unsigned)(SIEVE*8 - prev); |
| 586 | } |
| 587 | #else |
| 588 | p = SIEVE*8 - prev; |
| 589 | #endif |
| 590 | if (bnMulQ(&a, step, p) < 0 || bnAdd(bn, &a) < 0) |
| 591 | goto failed; |
| 592 | if (f && (retval = f(arg, '/')) < 0) |
| 593 | goto done; |
| 594 | } /* for (;;) */ |
| 595 | |
| 596 | failed: |
| 597 | retval = -1; |
| 598 | done: |
| 599 | bnEnd(&bn2); |
| 600 | bnEnd(&e); |
| 601 | bnEnd(&a); |
| 602 | #ifdef MSDOS |
| 603 | lbnMemFree(sieve, SIEVE); |
| 604 | #else |
| 605 | lbnMemWipe(sieve, sizeof(sieve)); |
| 606 | #endif |
| 607 | return retval < 0 ? retval : modexps+(order+1)*CONFIRMTESTS; |
| 608 | } |