bnprime.c

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

* 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. (bn is set to the last value tested, so you
* can increment bn and continue.)
*
* 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
bnPrimeGen(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 = bniMemAlloc(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);
						p = bnBits(&a);
						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! */
						pgpAssert(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;
			p = sieveSearch(sieve, SIEVE, p);
		}
		/*
		* Okay, from this point onwards, p is always the next entry
		* from the sieve, that has not been added to the shuffle table,
		* and is 0 iff the sieve has been exhausted.
		*
		* If we want to shuffle, then fill the shuffle table until the
		* sieve is exhausted or the table is full.
		*/
		if (rand && p) {
				do {
					offsets[i++] = p;
					p = sieveSearch(sieve, SIEVE, p);
				} while (p && 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];	/* The entry to use */

		/* Replace the entry with some more data, if possible */
		if (p) {
			offsets[j] = p;
			p = sieveSearch(sieve, SIEVE, p);
		} else {
			offsets[j] = offsets[--i];
			offsets[i] = 0;
		}

			/* Adjust bn to have the right value */
			if ((q > prev ? bnAddMult(bn, q-prev, 2)
					: bnSubMult(bn, prev-q, 2)) < 0)
			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 (bnAddMult(bn, SIEVE*8-prev, 2) < 0)
		goto failed;
	if (f && (retval = f(arg, '/')) < 0)
	goto done;
		goto retry;

failed:
		retval = -1;
done:
		bnEnd(&e);
		bnEnd(&a);
		bniMemWipe(offsets, sizeof(offsets));
#ifdef MSDOS
			bniMemFree(sieve, SIEVE);
#else
			bniMemWipe(sieve, sizeof(sieve));
#endif

		return retval < 0 ? retval : modexps + CONFIRMTESTS;
}

/*
* 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
bnPrimeGenStrong(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 = bniMemAlloc(SIEVE);
		if (!sieve)
			return -1;
#endif

		/* Step must be even and bn must be odd */
		pgpAssert((bnLSWord(step) & 1) == 0);
		pgpAssert((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.
								*/
								pgpAssert(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
		bniMemFree(sieve, SIEVE);
#else
		bniMemWipe(sieve, sizeof(sieve));
#endif
		return retval < 0 ? retval : modexps + CONFIRMTESTS;
}