August 6, 2013
A prime number p is a Sophie Germain prime if 2p+1 is also prime. The list of Sophie Germain primes begins 2, 3, 5, 11, 23, 29, 41, 53, 83, 89, 113, 131, (A005384); it is not known if the list is infinite. Sophie Germain primes are named after the French mathematician Marie-Sophie Germain, who in 1825 proved Fermat’s last theorem for such primes; thus, if p is a Sophie Germain prime, then there do not exist integers x, y and z not equal to zero and not a multiple of p such that xp + yp = zp.
One way to find Sophie Germain primes is to enumerate the primes p and check the primality of each 2p+1. That works fine for relatively small p. For larger p, amateur mathematician Harvey Dubner developed a sieve that speeds finding Sophie Germain primes:
Define the three polynomials fp(c) = c * 3003 * 10b – 1, fq(c) = 2 fp(c) + 1, and fr(c) = (fp(c) – 1) / 2. Choose b which is the approximate number of digits for the range of Sophie Germain primes you are searching.
Sieve each of the polynomials over the range 1 ≤ c ≤ m for some suitable m. Use as sieving primes the primes greater than 13 but less than some suitable f which is conveniently large but much less than the square root of fp(1). All three polynomials are sieved over the same sieving array. Sieving is done in a similar manner to sieving the polynomials in the multiple polynomial quadratic sieve.
After sieving is complete, apply a primality test to each fp(c) that survived all three sievings. If it is prime, check both fq(c) and fr(c) for primality. If fq(c) is prime, then fp(c) is a Sophie Germain prime. If fr(c), then fr(c) is a Sophie Germain prime.
Your task is to write a sieve that finds Sophie Germain primes as described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.