bngermain.c

著名的加密软件的应用于电子邮件中

/*
* bngermain.c - Sophie Germain prime generation using the bignum library
* and sieving.
*
* Copyright (C) 1995-1997 Pretty Good Privacy, Inc. All rights reserved.
*
* Written by Colin Plumb.
*
* $Id: bngermain.c,v 1.2.2.1 1997/06/07 09:49:28 mhw Exp $
	*/
#ifndef HAVE_CONFIG_H
#define HAVE_CONFIG_H 0
#endif
#if HAVE_CONFIG_H
#include "config.h"
#endif

#define BNDEBUG 1
#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 "rem", a partially reduced form of the prime,
* "div" 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 rem, unsigned div, unsigned order)
	{
			unsigned mul = 1;

rem %= div;
if (!rem)
	printf("bn div by %u!\n", div);
while (order--) {
rem += rem+1;
if (rem >= div)
	rem -= div;
mul += mul;
if (!rem)
			printf("%u*bn+%u div by %u!\n", mul, mul-1, div);
	}
	}
	#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(struct BigNum const *bn, struct BigNum *bn2, struct BigNum *e,
			struct 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++) {
			/*