Programming Praxis

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 x^p + y^p = z^p.

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 f_p(c) = c * 3003 * 10^b – 1, f_q(c) = 2 f_p(c) + 1, and f_r(c) = (f_p(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 f_p(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 f_p(c) that survived all three sievings. If it is prime, check both f_q(c) and f_r(c) for primality. If f_q(c) is prime, then f_p(c) is a Sophie Germain prime. If f_r(c), then f_r(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.

Pages: 1 2