Sophie Germain Primes
August 6, 2013
Let’s begin by reviewing how to sieve the polynomial f(x) over the the sequence f(1), f(2), … f(n). If p is prime and f(a) ≡ 0 (mod p), then f(a + k p) ≡ 0 (mod p) for all integers k. In our case, f(x) = c t − 1 for some t in each of the three polynomials, so we compute a as c t ≡ 1 (mod p) by calculating the inverse of t mod p. Thus we can sieve by setting each element of an array of length m to #t, calculating a and setting the ath element of the sieve to #f, then iterating through the sieve in steps of p, setting each a + k pth item to #f, repeating for each of the sieving primes. The three polynomials share a single sieve, since we are looking for cases where it is likely that at least two of them are prime:
(define (germain b f m)
(define (p c) (- (* c 3003 (expt 10 b)) 1))
(define (q c) (- (* c 6006 (expt 10 b)) 1))
(define (r c) (- (* c 15015 (expt 10 (- b 1))) 1))
(define (p-start p) (inverse (* 3003 (expt 10 b)) p))
(define (q-start p) (inverse (* 6006 (expt 10 b)) p))
(define (r-start p) (inverse (* 15015 (expt 10 (- b 1))) p))
(let ((sieve (make-vector (+ m 1) #t)) ; ignore sieve[0]
(ps (drop 6 (primes f)))) ; sieving primes
; sieve on p
(do ((ps ps (cdr ps))) ((null? ps))
(do ((i (p-start (car ps)) (+ i (car ps)))) ((< m i))
(vector-set! sieve i #f)))
; sieve on q
(do ((ps ps (cdr ps))) ((null? ps))
(do ((i (q-start (car ps)) (+ i (car ps)))) ((< m i))
(vector-set! sieve i #f)))
; sieve on r
(do ((ps ps (cdr ps))) ((null? ps))
(do ((i (r-start (car ps)) (+ i (car ps)))) ((< m i))
(vector-set! sieve i #f)))
; collect and test primes
(let ((gs (list))) ; sophie-germain primes
(do ((c 1 (+ c 1))) ((< m c) (sort < gs))
(when (and (vector-ref sieve c) (prime? (p c)))
(when (prime? (q c)) (set! gs (cons (p c) gs)))
(when (prime? (r c)) (set! gs (cons (r c) gs))))))))
Here b is the exponent of 10 in the formula, f is the largest sieving prime we want to use, and m is the size of the sieve. After all three polynomials are sieved, any numbers that survive the sieve are tested for primality, and Sophie Germain primes are retained to be output. Here’s an example:
> (for-each
(lambda (p) (display p) (newline))
(germain 25 100000 100000))
1831829999999999999999999999999
48873824999999999999999999999999
82912829999999999999999999999999
… 204 lines elided …
2900807909999999999999999999999999
2916273359999999999999999999999999
2940717779999999999999999999999999
You can use any convenient primality test. In fact, it’s not hard to prove the results prime, since adding 1 to each prime makes it very easy to factor, so a Lucas test will be both trivial and definitive. Dubner talks about using f = 30000000 and m = 4000000, and finds Sophie Germain primes for various b from 1000 to 5000; at the time it was written, Dubner’s program held the record for finding the largest Sophie Germain prime. Today, the largest known Sophie Germain prime is p = 18543637900515 · 2666668 − 1 which has 200701 decimal digits.
We used several functions from number theory that we have described in previous exercises. You can see and run the entire program at http://programmingpraxis.codepad.org/Ryqvrx2j.
Programming Praxis 