prime.c
上传用户:zbbssh
上传日期:2007-01-08
资源大小:196k
文件大小:20k
- /*
- * Prime generation using the bignum library and sieving.
- *
- * Copyright (c) 1995 Colin Plumb. All rights reserved.
- * For licensing and other legal details, see the file legal.c.
- */
- #if HAVE_CONFIG_H
- #include "config.h"
- #endif
- #if !NO_ASSERT_H
- #include <assert.h>
- #else
- #define assert(x) (void)0
- #endif
- #include <stdarg.h> /* We just can't live without this... */
- #if BNDEBUG
- #include <stdio.h>
- #endif
- #include "bn.h"
- #include "lbnmem.h"
- #include "prime.h"
- #include "sieve.h"
- #include "kludge.h"
- /* Size of the shuffle table */
- #define SHUFFLE 256
- /* Size of the sieve area */
- #define SIEVE 32768u/16
- /*
- * Helper function that does the slow primality test.
- * bn is the input bignum; a and e are temporary buffers that are
- * allocated by the caller to save overhead.
- *
- * Returns 0 if prime, >0 if not prime, and -1 on error (out of memory).
- * If not prime, returns the number of modular exponentiations performed.
- * Calls the fiven progress function with a '*' for each primality test
- * that is passed.
- *
- * The testing consists of strong pseudoprimality tests, to the bases
- * 2, 3, 5, 7, 11, 13 and 17. (Also called Miller-Rabin, although that's
- * not technically correct if we're using fixed bases.) Some people worry
- * that this might not be enough. Number theorists may wish to generate
- * primality proofs, but for random inputs, this returns non-primes with a
- * probability which is quite negligible, which is good enough.
- *
- * It has been proved (see Carl Pomerance, "On the Distribution of
- * Pseudoprimes", Math. Comp. v.37 (1981) pp. 587-593) that the number
- * of pseudoprimes (composite numbers that pass a Fermat test to the
- * base 2) less than x is bounded by:
- * exp(ln(x)^(5/14)) <= P_2(x) ### CHECK THIS FORMULA - it looks wrong! ###
- * P_2(x) <= x * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))).
- * Thus, the local density of Pseudoprimes near x is at most
- * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))), and at least
- * exp(ln(x)^(5/14) - ln(x)). Here are some values of this function
- * for various k-bit numbers x = 2^k:
- * Bits Density <= Bit equivalent Density >= Bit equivalent
- * 128 3.577869e-07 21.414396 4.202213e-37 120.840190
- * 192 4.175629e-10 31.157288 4.936250e-56 183.724558
- * 256 5.804314e-13 40.647940 4.977813e-75 246.829095
- * 384 1.578039e-18 59.136573 3.938861e-113 373.400096
- * 512 5.858255e-24 77.175803 2.563353e-151 500.253110
- * 768 1.489276e-34 112.370944 7.872825e-228 754.422724
- * 1024 6.633188e-45 146.757062 1.882404e-304 1008.953565
- *
- * As you can see, there's quite a bit of slop between these estimates.
- * In fact, the density of pseudoprimes is conjectured to be closer
- * to the square of that upper bound. E.g. the density of pseudoprimes
- * of size 256 is around 3 * 10^-27. The density of primes is very
- * high, from 0.005636 at 256 bits to 0.001409 at 1024 bits, i.e.
- * more than 10^-3.
- *
- * For those people used to cryptographic levels of security where the
- * 56 bits of DES key space is too small because it's exhaustible with
- * custom hardware searching engines, note that you are not generating
- * 50,000,000 primes per second on each of 56,000 custom hardware chips
- * for several hours. The chances that another Dinosaur Killer asteroid
- * will land today is about 10^-11 or 2^-36, so it would be better to
- * spend your time worrying about *that*. Well, okay, there should be
- * some derating for the chance that astronomers haven't seen it yet, but
- * I think you get the idea. For a good feel about the probability of
- * various events, I have heard that a good book is by E'mile Borel,
- * "Les Probabilite's et la vie". (The 's are accents, not apostrophes.)
- *
- * For more on the subject, try "Finding Four Million Large Random Primes",
- * by Ronald Rivest, in Advancess in Cryptology: Proceedings of Crypto '90.
- * He used a small-divisor test, then a Fermat test to the base 2, and then
- * 8 iterations of a Miller-Rabin test. About 718 million random 256-bit
- * integers were generated, 43,741,404 passed the small divisor test,
- * 4,058,000 passed the Fermat test, and all 4,058,000 passed all 8
- * iterations of the Miller-Rabin test, proving their primality beyond most
- * reasonable doubts.
- *
- * If the probability of getting a pseudoprime is some small p, then
- * the probability of not getting it in t trials is (1-p)^t. Remember
- * that, for small p, (1-p)^(1/p) ~ 1/e, the base of natural logarithms.
- * (This is more commonly expressed as e = lim_{xtoinfty} (1+1/x)^x.)
- * Thus, (1-p)^t ~ e^(-p*t) = exp(-p*t). So the odds of being able to
- * do this many tests without seeing a pseudoprime if you assume that
- * p = 10^-6 (one in a million) is one in 57.86. If you assume that
- * p = 2*10^-6, it's one in 3347.6. So it's implausible that the
- * density of pseudoprimes is much more than one millionth the density
- * of primes.
- *
- * He also gives a theoretical argument that the chance of finding a
- * 256-bit non-prime which satisfies one Fermat test to the base 2 is less
- * than 10^-22. The small divisor test improves this number, and if the
- * numbers are 512 bits (as needed for a 1024-bit key) the odds of failure
- * shrink to about 10^-44. Thus, he concludes, for practical purposes
- * *one* Fermat test to the base 2 is sufficient.
- */
- static int
- primeTest(struct BigNum const *bn, struct BigNum *e, struct BigNum *a,
- int (*f)(void *arg, int c), void *arg)
- {
- unsigned i, j;
- unsigned k, l;
- int err;
- static unsigned const primes[] = {2, 3, 5, 7, 11, 13, 17};
- #if BNDEBUG /* Debugging */
- /*
- * This is debugging code to test the sieving stage.
- * If the sieving is wrong, it will let past numbers with
- * small divisors. The prime test here will still work, and
- * weed them out, but you'll be doing a lot more slow tests,
- * and presumably excluding from consideration some other numbers
- * which might be prime. This check just verifies that none
- * of the candidates have any small divisors. If this
- * code is enabled and never triggers, you can feel quite
- * confident that the sieving is doing its job.
- */
- i = bnModQ(bn, 30030); /* 30030 = 2 * 3 * 5 * 7 * 11 * 13 */
- if (!(i % 2)) printf("bn div by 2!");
- if (!(i % 3)) printf("bn div by 3!");
- if (!(i % 5)) printf("bn div by 5!");
- if (!(i % 7)) printf("bn div by 7!");
- if (!(i % 11)) printf("bn div by 11!");
- if (!(i % 13)) printf("bn div by 13!");
- i = bnModQ(bn, 7429); /* 7429 = 17 * 19 * 23 */
- if (!(i % 17)) printf("bn div by 17!");
- if (!(i % 19)) printf("bn div by 19!");
- if (!(i % 23)) printf("bn div by 23!");
- i = bnModQ(bn, 33263); /* 33263 = 29 * 31 * 37 */
- if (!(i % 29)) printf("bn div by 29!");
- if (!(i % 31)) printf("bn div by 31!");
- if (!(i % 37)) printf("bn div by 37!");
- #endif
- /*
- * Now, check that bn is prime. If it passes to the base 2,
- * it's prime beyond all reasonable doubt, and everything else
- * is just gravy, but it gives people warm fuzzies to do it.
- *
- * This starts with verifying Euler's criterion for a base of 2.
- * This is the fastest pseudoprimality test that I know of,
- * saving a modular squaring over a Fermat test, as well as
- * being stronger. 7/8 of the time, it's as strong as a strong
- * pseudoprimality test, too. (The exception being when bn ==
- * 1 mod 8 and 2 is a quartic residue, i.e. bn is of the form
- * a^2 + (8*b)^2.) The precise series of tricks used here is
- * not documented anywhere, so here's an explanation.
- * Euler's criterion states that if p is prime then a^((p-1)/2)
- * is congruent to Jacobi(a,p), modulo p. Jacobi(a,p) is
- * a function which is +1 if a is a square modulo p, and -1 if
- * it is not. For a = 2, this is particularly simple. It's
- * +1 if p == +/-1 (mod 8), and -1 if m == +/-3 (mod 8).
- * If p == 3 mod 4, then all a strong test does is compute
- * 2^((p-1)/2). and see if it's +1 or -1. (Euler's criterion
- * says *which* it should be.) If p == 5 (mod 8), then
- * 2^((p-1)/2) is -1, so the initial step in a strong test,
- * looking at 2^((p-1)/4), is wasted - you're not going to
- * find a +/-1 before then if it *is* prime, and it shouldn't
- * have either of those values if it isn't. So don't bother.
- *
- * The remaining case is p == 1 (mod 8). In this case, we
- * expect 2^((p-1)/2) == 1 (mod p), so we expect that the
- * square root of this, 2^((p-1)/4), will be +/-1 (mod p).
- * Evaluating this saves us a modular squaring 1/4 of the time.
- * If it's -1, a strong pseudoprimality test would call p
- * prime as well. Only if the result is +1, indicating that
- * 2 is not only a quadratic residue, but a quartic one as well,
- * does a strong pseudoprimality test verify more things than
- * this test does. Good enough.
- *
- * We could back that down another step, looking at 2^((p-1)/8)
- * if there was a cheap way to determine if 2 were expected to
- * be a quartic residue or not. Dirichlet proved that 2 is
- * a quartic residue iff p is of the form a^2 + (8*b^2).
- * All primes == 1 (mod 4) can be expressed as a^2 + (2*b)^2,
- * but I see no cheap way to evaluate this condition.
- */
- if (bnCopy(e, bn) < 0)
- return -1;
- (void)bnSubQ(e, 1);
- l = bnLSWord(e);
- j = 1; /* Where to start in prime array for strong prime tests */
- if (l & 7) {
- bnRShift(e, 1);
- if (bnTwoExpMod(a, e, bn) < 0)
- return -1;
- if ((l & 7) == 6) {
- /* bn == 7 mod 8, expect +1 */
- if (bnBits(a) != 1)
- return 1; /* Not prime */
- k = 1;
- } else {
- /* bn == 3 or 5 mod 8, expect -1 == bn-1 */
- if (bnAddQ(a, 1) < 0)
- return -1;
- if (bnCmp(a, bn) != 0)
- return 1; /* Not prime */
- k = 1;
- if (l & 4) {
- /* bn == 5 mod 8, make odd for strong tests */
- bnRShift(e, 1);
- k = 2;
- }
- }
- } else {
- /* bn == 1 mod 8, expect 2^((bn-1)/4) == +/-1 mod bn */
- bnRShift(e, 2);
- if (bnTwoExpMod(a, e, bn) < 0)
- return -1;
- if (bnBits(a) == 1) {
- j = 0; /* Re-do strong prime test to base 2 */
- } else {
- if (bnAddQ(a, 1) < 0)
- return -1;
- if (bnCmp(a, bn) != 0)
- return 1; /* Not prime */
- }
- k = 2 + bnMakeOdd(e);
- }
- /* It's prime! Now go on to confirmation tests */
- /*
- * Now, e = (bn-1)/2^k is odd. k >= 1, and has a given value
- * with probability 2^-k, so its expected value is 2.
- * j = 1 in the usual case when the previous test was as good as
- * a strong prime test, but 1/8 of the time, j = 0 because
- * the strong prime test to the base 2 needs to be re-done.
- */
- for (i = j; i < sizeof(primes)/sizeof(*primes); i++) {
- if (f && (err = f(arg, '*')) < 0)
- return err;
- (void)bnSetQ(a, primes[i]);
- if (bnExpMod(a, a, e, bn) < 0)
- return -1;
- if (bnBits(a) == 1)
- continue; /* Passed this test */
- l = k;
- for (;;) {
- if (bnAddQ(a, 1) < 0)
- return -1;
- if (bnCmp(a, bn) == 0) /* Was result bn-1? */
- break; /* Prime */
- if (!--l) /* Reached end, not -1? luck? */
- return i+2-j; /* Failed, not prime */
- /* This portion is executed, on average, once. */
- (void)bnSubQ(a, 1); /* Put a back where it was. */
- if (bnSquare(a, a) < 0 || bnMod(a, a, bn) < 0)
- return -1;
- if (bnBits(a) == 1)
- return i+2-j; /* Failed, not prime */
- }
- /* It worked (to the base primes[i]) */
- }
-
- /* Yes, we've decided that it's prime. */
- if (f && (err = f(arg, '*')) < 0)
- return err;
- return 0; /* Prime! */
- }
- /*
- * Modifies the bignum to return a nearby (slightly larger) number which
- * is a probable prime. Returns >=0 on success or -1 on failure (out of
- * memory). The return value is the number of unsuccessful modular
- * exponentiations performed. This never gives up searching.
- *
- * All other arguments are optional. They may be NULL. They are:
- *
- * unsigned (*rand)(unsigned limit)
- * For better distributed numbers, supply a non-null pointer to a
- * function which returns a random x, 0 <= x < limit. (It may make it
- * simpler to know that 0 < limit <= SHUFFLE, so you need at most a byte.)
- * The program generates a large window of sieve data and then does
- * pseudoprimality tests on the data. If a rand function is supplied,
- * the candidates which survive sieving are shuffled with a window of
- * size SHUFFLE before testing to increase the uniformity of the prime
- * selection. This isn't perfect, but it reduces the correlation between
- * the size of the prime-free gap before a prime and the probability
- * that that prime will be found by a sequential search.
- *
- * If rand is NULL, sequential search is used. If you want sequential
- * search, note that the search begins with the given number; if you're
- * trying to generate consecutive primes, you must increment the previous
- * one by two before calling this again.
- *
- * int (*f)(void *arg, int c), void *arg
- * The function f argument, if non-NULL, is called with progress indicator
- * characters for printing. A dot (.) is written every time a primality test
- * is failed, a star (*) every time one is passed, and a slash (/) in the
- * (very rare) case that the sieve was emptied without finding a prime
- * and is being refilled. f is also passed the void *arg argument for
- * private context storage. If f returns < 0, the test aborts and returns
- * that value immediately.
- *
- * The "exponent" argument, and following unsigned numbers, are exponents
- * for which an inverse is desired, modulo p. For a d to exist such that
- * (x^e)^d == x (mod p), then d*e == 1 (mod p-1), so gcd(e,p-1) must be 1.
- * The prime returned is constrained to not be congruent to 1 modulo
- * any of the zero-terminated list of 16-bit numbers. Note that this list
- * should contain all the small prime factors of e. (You'll have to test
- * for large prime factors of e elsewhere, but the chances of needing to
- * generate another prime are low.)
- *
- * The list is terminated by a 0, and may be empty.
- *
- */
- int
- primeGen(struct BigNum *bn, unsigned (*rand)(unsigned),
- int (*f)(void *arg, int c), void *arg, unsigned exponent, ...)
- {
- int retval;
- int modexps = 0;
- unsigned short offsets[SHUFFLE];
- unsigned i, j;
- unsigned p, q, prev;
- struct BigNum a, e;
- #ifdef MSDOS
- unsigned char *sieve;
- #else
- unsigned char sieve[SIEVE];
- #endif
- #ifdef MSDOS
- sieve = lbnMemAlloc(SIEVE);
- if (!sieve)
- return -1;
- #endif
- bnBegin(&a);
- bnBegin(&e);
- #if 0 /* Self-test (not used for production) */
- {
- struct BigNum t;
- static unsigned char const prime1[] = {5};
- static unsigned char const prime2[] = {7};
- static unsigned char const prime3[] = {11};
- static unsigned char const prime4[] = {1, 1}; /* 257 */
- static unsigned char const prime5[] = {0xFF, 0xF1}; /* 65521 */
- static unsigned char const prime6[] = {1, 0, 1}; /* 65537 */
- static unsigned char const prime7[] = {1, 0, 3}; /* 65539 */
- /* A small prime: 1234567891 */
- static unsigned char const prime8[] = {0x49, 0x96, 0x02, 0xD3};
- /* A slightly larger prime: 12345678901234567891 */
- static unsigned char const prime9[] = {
- 0xAB, 0x54, 0xA9, 0x8C, 0xEB, 0x1F, 0x0A, 0xD3 };
- /*
- * No, 123456789012345678901234567891 isn't prime; it's just a
- * lucky, easy-to-remember conicidence. (You have to go to
- * ...4567907 for a prime.)
- */
- static struct {
- unsigned char const *prime;
- unsigned size;
- } const primelist[] = {
- { prime1, sizeof(prime1) },
- { prime2, sizeof(prime2) },
- { prime3, sizeof(prime3) },
- { prime4, sizeof(prime4) },
- { prime5, sizeof(prime5) },
- { prime6, sizeof(prime6) },
- { prime7, sizeof(prime7) },
- { prime8, sizeof(prime8) },
- { prime9, sizeof(prime9) } };
- bnBegin(&t);
- for (i = 0; i < sizeof(primelist)/sizeof(primelist[0]); i++) {
- bnInsertBytes(&t, primelist[i].prime, 0,
- primelist[i].size);
- bnCopy(&e, &t);
- (void)bnSubQ(&e, 1);
- bnTwoExpMod(&a, &e, &t);
- p = bnBits(&a);
- if (p != 1) {
- printf(
- "Bug: Fermat(2) %u-bit output (1 expected)n", p);
- fputs("Prime = 0x", stdout);
- for (j = 0; j < primelist[i].size; j++)
- printf("%02X", primelist[i].prime[j]);
- putchar('n');
- }
- bnSetQ(&a, 3);
- bnExpMod(&a, &a, &e, &t);
- if (p != 1) {
- printf(
- "Bug: Fermat(3) %u-bit output (1 expected)n", p);
- fputs("Prime = 0x", stdout);
- for (j = 0; j < primelist[i].size; j++)
- printf("%02X", primelist[i].prime[j]);
- putchar('n');
- }
- }
- bnEnd(&t);
- }
- #endif
- /* First, make sure that bn is odd. */
- if ((bnLSWord(bn) & 1) == 0)
- (void)bnAddQ(bn, 1);
- retry:
- /* Then build a sieve starting at bn. */
- sieveBuild(sieve, SIEVE, bn, 2, 0);
- /* Do the extra exponent sieving */
- if (exponent) {
- va_list ap;
- unsigned t = exponent;
- va_start(ap, exponent);
- do {
- /* The exponent had better be odd! */
- assert(t & 1);
- i = bnModQ(bn, t);
- /* Find 1-i */
- if (i == 0)
- i = 1;
- else if (--i)
- i = t - i;
- /* Divide by 2, modulo the exponent */
- i = (i & 1) ? i/2 + t/2 + 1 : i/2;
- /* Remove all following multiples from the sieve. */
- sieveSingle(sieve, SIEVE, i, t);
- /* Get the next exponent value */
- t = va_arg(ap, unsigned);
- } while (t);
- va_end(ap);
- }
- /* Fill up the offsets array with the first SHUFFLE candidates */
- i = p = 0;
- /* Get first prime */
- if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0)
- offsets[i++] = p;
- /*
- * If we have a prime and we have a shuffle function so it makes
- * to fill up the table, so do.
- */
- if (rand && i) {
- do {
- p = sieveSearch(sieve, SIEVE, p);
- if (p == 0)
- break;
- offsets[i++] = p;
- } while (i < SHUFFLE);
- }
- /* Choose a random candidate for experimentation */
- prev = 0;
- while (i) {
- /* Pick a random entry from the shuffle table */
- j = rand ? rand(i) : 0;
- q = offsets[j];
- /* Replace the entry with some more data, if possible */
- if (p && (p = sieveSearch(sieve, SIEVE, p)) != 0)
- offsets[j] = p;
- else
- offsets[j] = offsets[--i];
- /* Adjust bn to have the right value */
- if (q > prev) {
- if (bnAddQ(bn, q-prev) < 0)
- goto failed;
- if (bnAddQ(bn, q-prev) < 0)
- goto failed;
- } else {
- if (bnSubQ(bn, prev-q))
- goto failed;
- if (bnSubQ(bn, prev-q))
- goto failed;
- }
- prev = q;
- /* Now do the Fermat tests */
- retval = primeTest(bn, &e, &a, f, arg);
- if (retval <= 0)
- goto done; /* Success or error */
- modexps += retval;
- if (f && (retval = f(arg, '.')) < 0)
- goto done;
- }
- /* Ran out of sieve space - increase bn and keep trying. */
- if (bnSubQ(bn, (SIEVE*8)-prev))
- goto failed;
- if (bnSubQ(bn, (SIEVE*8)-prev))
- goto failed;
- if (f && (retval = f(arg, '/')) < 0)
- goto done;
- goto retry;
- failed:
- retval = -1;
- done:
- bnEnd(&e);
- bnEnd(&a);
- lbnMemWipe(offsets, sizeof(offsets));
- #ifdef MSDOS
- lbnMemFree(sieve, SIEVE);
- #else
- lbnMemWipe(sieve, sizeof(sieve));
- #endif
- return retval < 0 ? retval : modexps;
- }
- /*
- * Similar, but searches forward from the given starting value in steps of
- * "step" rather than 1. The step size must be even, and bn must be odd.
- * Among other possibilities, this can be used to generate "strong"
- * primes, where p-1 has a large prime factor.
- */
- int
- primeGenStrong(struct BigNum *bn, struct BigNum const *step,
- int (*f)(void *arg, int c), void *arg)
- {
- int retval;
- unsigned p, prev;
- struct BigNum a, e;
- int modexps = 0;
- #ifdef MSDOS
- unsigned char *sieve;
- #else
- unsigned char sieve[SIEVE];
- #endif
- #ifdef MSDOS
- sieve = lbnMemAlloc(SIEVE);
- if (!sieve)
- return -1;
- #endif
- /* Step must be even and bn must be odd */
- assert((bnLSWord(step) & 1) == 0);
- assert((bnLSWord(bn) & 1) == 1);
- bnBegin(&a);
- bnBegin(&e);
- for (;;) {
- if (sieveBuildBig(sieve, SIEVE, bn, step, 0) < 0)
- goto failed;
- p = prev = 0;
- if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) {
- do {
- /*
- * Adjust bn to have the right value,
- * adding (p-prev) * 2*step.
- */
- assert(p >= prev);
- /* Compute delta into a */
- if (bnMulQ(&a, step, p-prev) < 0)
- goto failed;
- if (bnAdd(bn, &a) < 0)
- goto failed;
- prev = p;
- retval = primeTest(bn, &e, &a, f, arg);
- if (retval <= 0)
- goto done; /* Success! */
- modexps += retval;
- if (f && (retval = f(arg, '.')) < 0)
- goto done;
- /* And try again */
- p = sieveSearch(sieve, SIEVE, p);
- } while (p);
- }
- /* Ran out of sieve space - increase bn and keep trying. */
- #if SIEVE*8 == 65536
- /* Corner case that will never actually happen */
- if (!prev) {
- if (bnAdd(bn, step) < 0)
- goto failed;
- p = 65535;
- } else {
- p = (unsigned)(SIEVE*8 - prev);
- }
- #else
- p = SIEVE*8 - prev;
- #endif
- if (bnMulQ(&a, step, p) < 0 || bnAdd(bn, &a) < 0)
- goto failed;
- if (f && (retval = f(arg, '/')) < 0)
- goto done;
- } /* for (;;) */
- failed:
- retval = -1;
- done:
- bnEnd(&e);
- bnEnd(&a);
- #ifdef MSDOS
- lbnMemFree(sieve, SIEVE);
- #else
- lbnMemWipe(sieve, sizeof(sieve));
- #endif
- return retval < 0 ? retval : modexps;
- }