bngermain.c

vc环境下的pgp源码

		 * 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(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(BigNum *bn, unsigned order,
	int (*f)(void *arg, int c), void *arg)
{
	int retval;
	unsigned p, prev;
	unsigned inc;
	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;
}