generating provable primes

[Wei Dai <weidai@eskimo.com>]

Several days ago someone (I forgot who he was) asked about code for 
primality tests.  I just implemented an algorithm to generate random 
provable primes that is only about 50% slower than generating probable 
primes.  It will be in the next version of Crypto++, but I've attached 
code for the main function in case anyone is interested in the 
algorithm.  Full description can be found in "Fast Generation of Prime 
Numbers and Secure Public-Key Cryptographic Parameters" by U.M. Maurer in 
Journal of Cryptology, Volume 8 Number 3, 1995.  The paper also describes 
a more complicated algorithm that produces primes with a more uniform 
distribution.

There was discussion some days ago about generating strong primes 
for DH exchange moduli.  Eric Young reported that he spent tens of hours 
of CPU time to generate a 4096 bit prime p such that (p-1)/2 is also 
prime.  However, there is really no reason why DH exchange moduli must be of 
the form 2q+1 where q is a prime.  It should be sufficient that they are 
of the form rq+1, where q is a large enough prime (say more than 256 
bits).  The following algorithm generates a provable prime p=2rq+1, where q 
is a prime with at least half the length of p.

bignum ProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
{
	const unsigned smallPrimeBound = 29, c_opt=10;
	bignum p;

	BuildPrimeTable();
	if (bits < smallPrimeBound)
	{
		do
			p.Randomize(rng, bignum::Power2(bits-1), 
bignum::Power2(bits)-1, ODD);
		while (TrialDivision(p, 1 << ((bits+1)/2)));
	}
	else
	{
		const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
		double relativeSize;
		do
			relativeSize = pow(2.0, 
double(rng.GetLong())/0xffffffff - 1);
		while (bits * relativeSize >= bits - margin);
		
		bignum a,b;
		bignum q = ProvablePrime(rng, unsigned(bits*relativeSize));
		bignum I = bignum::Power2(bits-2)/q;
		bignum I2 = I << 1;
		unsigned int trialDivisorBound = (unsigned 
int)min((unsigned long)primeTable[primeTableSize-1], (unsigned 
long)bits*bits/c_opt);
		boolean success = FALSE;
		do
		{
			p.Randomize(rng, I, I2, ANY);
			p *= q; p <<= 1; ++p;
			if (!TrialDivision(p, trialDivisorBound))
			{
				a.Randomize(rng, 2, p-1, ANY);
				b = a.ExponentiateMod((p-1)/q, p);
				success = (Gcd(b-1, p) == 1) && 
(b.ExponentiateMod(q, p) == 1);
			}
		}
		while (!success);
	}
	return p;
}