bngermain.c

vc环境下的pgp源码

/*
 * Sophie Germain prime generation using the bignum library and sieving.
 *
 * Written by Colin Plumb.
 *
 * $Id: bngermain.c,v 1.6 1997/12/11 23:12:01 lloyd Exp $
 */
#include "pgpConfig.h"

#ifndef BNDEBUG
#define BNDEBUG 0
#endif
#if BNDEBUG
#include <stdio.h>
#endif

#include "bn.h"
#include "bngermain.h"
#include "bnjacobi.h"
#include "bnimem.h"	/* For bniMemWipe */
#include "bnsieve.h"

#include "bnkludge.h"
#include "pgpDebug.h"

/* Size of the sieve area (can be up to 65536/8 = 8192) */
#define SIEVE 8192

static unsigned const confirm[] = {2, 3, 5, 7, 11, 13, 17};
#define CONFIRMTESTS (sizeof(confirm)/sizeof(*confirm))

#if BNDEBUG
/*
 * For sanity checking the sieve, we check for small divisors of the numbers
 * we get back.  This takes "remainder", a partially reduced form of the
 * prime (the prime modulo some multiple of the divisor) a divisor to check
 * for, and "order", a parameter of the "order" of Sophie Germain primes
 * (0 = normal primes, 1 = Sophie Germain primes, 2 = 4*p+3 is also prime,
 * etc.) and does the check.  It just complains to stdout if the check fails.
 */
static void
germainSanity(unsigned remainder, unsigned divisor, unsigned order)
{
	unsigned o;

	remainder %= divisor;
	if (!remainder)
		printf("bn is divisible by %u!\n", divisor);
	for (o = 1; o < order; o++) {
		remainder += remainder+1;
		if (remainder >= divisor)
			remainder -= divisor;
		mul += mul;
		if (!remainder)
			printf("%u*bn+%u is divisible by %u!\n",
			       1u<<o, (1u<<o)-1, divisor);
	}
}
#endif /* BNDEBUG */

/*
 * Helper function that does the slow primality test.
 * bn is the input bignum; a, e and bn2 are temporary buffers that are
 * allocated by the caller to save overhead.  bn2 is filled with
 * a copy of 2^order*bn+2^order-1 if bn is found to be prime.
 *
 * Returns 0 if both bn and bn2 are prime, >0 if not prime, and -1 on
 * error (out of memory).  If not prime, the return value is the number
 * of modular exponentiations performed.   Prints a '+' or '-' on the
 * given FILE (if any) for each test that is passed by bn, and a '*'
 * for each test that is passed by bn2.
 *
 * The testing consists of strong pseudoprimality tests, to the bases given
 * in the confirm[] array above.  (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_{x\to\infty} (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
germainPrimeTest(BigNum const *bn, BigNum *bn2, BigNum *e,
	BigNum *a, unsigned order, int (*f)(void *arg, int c),
	void *arg)
{
	int err;
	unsigned i;
	int j;
	unsigned k, l, n;

#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 = bnLSWord(bn);
	if (!(i % 2)) printf("bn div by 2!");
	i = bnModQ(bn, 51051);	/* 51051 = 3 * 7 * 11 * 13 * 17 */
	germainSanity(i, 3, order);
	germainSanity(i, 7, order);
	germainSanity(i, 11, order);
	germainSanity(i, 13, order);
	germainSanity(i, 17, order);
	i = bnModQ(bn, 63365);	/* 63365 = 5 * 19 * 23 * 29 */
	germainSanity(i, 5, order);
	germainSanity(i, 19, order);
	germainSanity(i, 23, order);
	germainSanity(i, 29, order);
	i = bnModQ(bn, 47027);	/* 47027 = 31 * 37 * 41 */
	germainSanity(i, 31, order);
	germainSanity(i, 37, order);
	germainSanity(i, 41, order);
#endif
	/*
	 * First, check whether bn is prime.  This uses a fast primality
	 * test which usually obviates the need to do one of the
	 * confirmation tests later.  See bnprime.c for a full explanation.
	 * We check bn first because it's one bit smaller, saving one
	 * modular squaring, and because we might be able to save another
	 * when testing it.  (1/4 of the time.)  A small speed hack,
	 * but finding big Sophie Germain primes is *slow*.
	 */
	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 check higher-order forms bn2 = 2*bn+1, 4*bn+3,
	 * etc.  Since bn2 == 3 mod 4, a strong pseudoprimality test boils
	 * down to looking at a^((bn2-1)/2) mod bn and seeing if it's +/-1.
	 * (+1 if bn2 is == 7 mod 8, -1 if it's == 3)
	 * Of course, that exponent is just the previous bn2 or bn...
	 */
	if (bnCopy(bn2, bn) < 0)
			return -1;
	for (n = 0; n < order; n++) {
		/*
		 * Print a success indicator: the sign of Jacobi(2,bn2),
		 * which is available to us in l.  bn2 = 2*bn + 1.  Since bn
		 * is odd, bn2 must be == 3 mod 4, so the options modulo 8
		 * are 3 and 7.  3 if l == 1 mod 4, 7 if l == 3 mod 4.
		 * The sign of the Jacobi symbol is - and + for these cases,