bngermain.c

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

			* 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,
			* respectively.
			*/
			if (f && (err = f(arg, "-+"[(l >> 1) & 1])) < 0)
				 return err;
			/* Exponent is previous bn2 */
			if (bnCopy(e, bn2) < 0 || bnLShift(bn2, 1) < 0)
				 return -1;
			(void)bnAddQ(bn2, 1);	/* Can't overflow */
			if (bnTwoExpMod(a, e, bn2) < 0)
				 return -1;
			if (n | l) {	 /* Expect + */
				 if (bnBits(a) != 1)
				 	return 2+n;			/* Not prime */
			} else {
					if (bnAddQ(a, 1) < 0)
						return -1;
					if (bnCmp(a, bn2) != 0)
						return 2+n;		/* Not prime */
			}
			l = bnLSWord(bn2);
	}

/* Final success indicator - it's in the bag. */
if (f && (err = f(arg, '*')) < 0)
	return err;

/*
* Success! We have found a prime! Now go on to confirmation
* tests... k is an amount by which we know it's safe to shift
				* down e. j = 1 unless the test to the base 2 could stand to be
				* re-done (it wasn't *quite* a strong test), in which case it's 0.
				*
				* Here, we do the full strong pseudoprimality test. This proves
				* that a number is composite, or says that it's probably prime.
				*
				* For the given base a, find bn-1 = 2^k * e, then find
				* x == a^e (mod bn).
				* If x == +1 -> strong pseudoprime to base a
				* Otherwise, repeat k times:
				* If x == -1, -> strong pseudoprime
				* x = x^2 (mod bn)
				* If x = +1 -> composite
				* If we reach the end of the iteration and x is *not* +1, at the
				* end, it is composite. Which means that the squaring actually
				* only has to proceed k-1 times. If x is not -1 by then,
				* it's composite no matter what the result of the squaring is.
				*
				* For the multiples 2*bn+1, 4*bn+3, etc. then k = 1 (and e is
				* the previous multiple of bn) so there is no need to do any
				* looping at all.
				*/
				for (i = j; i < CONFIRMTESTS; i++) {
					if (bnCopy(e, bn) < 0)
							 return -1;
					bnRShift(e, k);
					k += bnMakeOdd(e);
					(void)bnSetQ(a, confirm[i]);
					if (bnExpMod(a, a, e, bn) < 0)
						return -1;

					if (bnBits(a) != 1) {
							l = k;
							for (;;) {
									if (bnAddQ(a, 1) < 0)
										return -1;
									if (bnCmp(a, bn) == 0)	 /* Was result bn-1? */
										break;	/* Prime */
									if (!--l)
										return (1+order)*i+2;	/* Fail */
									/* This part is executed once, on average. */
									(void)bnSubQ(a, 1);	/* Restore a */
									if (bnSquare(a, a) < 0 || bnMod(a, a, bn) < 0)
										return -1;
									if (bnBits(a) == 1)
										return (1+order)*i+1;	/* Fail */
							}
					}

					if (bnCopy(bn2, bn) < 0)
						return -1;
		
					/* Only do the following if we're not re-doing base 2 */
					if (i) for (n = 0; n < order; n++) {
						if (bnCopy(e, bn2) < 0 || bnLShift(bn2, 1) < 0)
							 return -1;
						(void)bnAddQ(bn2, 1);

							/* Print success indicator for previous test */
							j = bnJacobiQ(confirm[i], bn2);
							if (f && (err = f(arg, j < 0 ? '-' : '+')) < 0)
								 return err;

							/* Check that p^e == Jacobi(p,bn2) (mod bn2) */
							(void)bnSetQ(a, confirm[i]);
							if (bnExpMod(a, a, e, bn2) < 0)
								 return -1;
							/*
							* FIXME: Actually, we don't need to compute the
							* Jacobi symbol externally... it never happens that
							* a = +/-1 but it's the wrong one. So we can just
							* look at a and use its sign. Find a proof somewhere.
							*/
							if (j < 0) {
									/* Not a Q.R., should have a = bn2-1 */
									if (bnAddQ(a, 1) < 0)
										return -1;
									if (bnCmp(a, bn2) != 0)	/* Was result bn2-1? */
										return (1+order)*i+n+2;	/* Fail */
							} else {
									/* Quadratic residue, should have a = 1 */
									if (bnBits(a) != 1)
										return (1+order)*i+n+2;	/* Fail */
							}
					}
					/* Final success indicator for the base confirm[i]. */
					if (f && (err = f(arg, '*')) < 0)
						return err;
			}

		return 0;	 /* Prime! */
}

/*
* Add x*y to bn, which is usually (but not always) < 65536.
* Do it in a simple linear manner.
	*/
static int
bnAddMult(struct BigNum *bn, unsigned long x, unsigned y)
	{
			unsigned long z = (unsigned long)x * y;

			while (z > 65535) {
				 if (bnAddQ(bn, 65535) < 0)
				 	return -1;
				 z -= 65535;
			}
			return bnAddQ(bn, (unsigned)z);
	}

	/*
* Modifies the bignum to return the next Sophie Germain prime >= the
* input value. Sohpie Germain primes are number such that p is
* prime and 2*p+1 is also prime.
*
* This is actually parameterized: it generates primes p such that "order"
* multiples-plus-two are also prime, 2*p+1, 2*(2*p+1)+1 = 4*p+3, etc.
*
* Returns >=0 on success or -1 on failure (out of memory). On success,
* the return value is the number of modular exponentiations performed
* (excluding the final confirmations). This never gives up searching.
*
* The FILE *f argument, if non-NULL, has progress indicators written
* to it. A dot (.) is written every time a primeality test is failed,
* a plus (+) or minus (-) when the smaller prime of the pair passes a
* test, and a star (*) when the larger one does. Finally, a slash (/)
* is printed when the sieve was emptied without finding a prime and is
* being refilled.
*
* Apologies to structured programmers for all the GOTOs.
	*/
int
bnGermainPrimeGen(struct BigNum *bn, unsigned order,
			int (*f)(void *arg, int c), void *arg)
	{
			int retval;
			unsigned p, prev;
			unsigned inc;
			struct BigNum a, e, bn2;
			int modexps = 0;
#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);
			bnBegin(&bn2);

			/*
			* Obviously, the prime we find must be odd. Further, if 2*p+1
			* is also to be prime (order > 0) then p != 1 (mod 3), lest
			* 2*p+1 == 3 (mod 3). Added to p != 3 (mod 3), p == 2 (mod 3)
			* and p == 5 (mod 6).
			* If order > 2 and we care about 4*p+3 and 8*p+7, then similarly
			* p == 4 (mod 5), so p == 29 (mod 30).
			* So pick the step size for searching based on the order
			* and increse bn until it's == -1 (mod inc).
			*
			* mod 7 doesn't have a unique value for p because 2 -> 5 -> 4 -> 2,
			* nor does mod 11, and I don't want to think about things past
			* that. The required order would be impractically high, in any case.
			*/
			inc = order ? ((order > 2) ? 30 : 6) : 2;
			if (bnAddQ(bn, inc-1 - bnModQ(bn, inc)) < 0)
				 goto failed;

			for (;;) {
					if (sieveBuild(sieve, SIEVE, bn, inc, order) < 0)
						goto failed;

					p = prev = 0;
					if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) {
							do {
									/* Adjust bn to have the right value. */
									pgpAssert(p >= prev);
									if (bnAddMult(bn, p-prev, inc) < 0)
										goto failed;
									prev = p;

									/* Okay, do the strong tests. */
									retval = germainPrimeTest(bn, &bn2, &e, &a,
									order, f, arg);
									if (retval <= 0)
										goto done;
									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 (bnAddMult(bn, (unsigned long)SIEVE*8-prev, inc) < 0)
						goto failed;
					if (f && (retval = f(arg, '/')) < 0)
							goto done;
			} /* for (;;) */

	failed:
			retval = -1;
	done:
			bnEnd(&bn2);
			bnEnd(&e);
			bnEnd(&a);
#ifdef MSDOS
			bniMemFree(sieve, SIEVE);
#else
			bniMemWipe(sieve, sizeof(sieve));
#endif
			return retval < 0 ? retval : modexps+(order+1)*CONFIRMTESTS;
	}