## Twin Primes

### July 26, 2019

The primary purpose of this exercise is to demonstrate the `with-modulus`

environment that we wrote ten years ago, augmented here with a modular

factorial:

(define prime? ; strong pseudoprime to prime bases less than 100 (let* ((ps (list 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)) (p100 (apply * ps))) (lambda (n) (define (expm b e m) (let loop ((b b) (e e) (x 1)) (if (zero? e) x (loop (modulo (* b b) m) (quotient e 2) (if (odd? e) (modulo (* b x) m) x))))) (define (spsp? n a) ; #t if n is a strong pseudoprime base a (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1))) ((odd? d) (if (= (expm a d n) 1) #t (do ((r 0 (+ r 1))) ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s)) (< r s))))))) (if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f) (do ((ps ps (cdr ps))) ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))

(define (euclid x y) (let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w y)) (if (zero? w) (values a b g) (let ((q (quotient g w))) (loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w)))))))

(define (inverse x m) (if (not (= (gcd x m) 1)) (error 'inverse "divisor must be coprime to modulus") (call-with-values (lambda () (euclid x m)) (lambda (a b g) (modulo a m)))))

(define (expm b e m) (define (m* x y) (modulo (* x y) m)) (cond ((zero? e) 1) ((even? e) (expm (m* b b) (/ e 2) m)) (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(define (jacobi a n) (if (not (and (integer? a) (integer? n) (positive? n) (odd? n))) (error 'jacobi "modulus must be positive odd integer") (let jacobi ((a a) (n n)) (cond ((= a 0) 0) ((= a 1) 1) ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1))) ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n))) ((< n a) (jacobi (modulo a n) n)) ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a))) (else (jacobi n a))))))

(define (mod-fact n m) (if (<= m n) 0 (let loop ((k 2) (p 1)) (if (zero? p) 0 (if (< n k) p (loop (+ k 1) (modulo (* p k) m)))))))

(define (mod-sqrt a p) (define (both n) (list n (- p n))) (cond ((not (and (odd? p) (prime? p))) (error 'mod-sqrt "modulus must be an odd prime")) ((not (= (jacobi a p) 1)) (error 'mod-sqrt "must be a quadratic residual")) (else (let ((a (modulo a p))) (case (modulo p 8) ((3 7) (both (expm a (/ (+ p 1) 4) p))) ((5) (let* ((x (expm a (/ (+ p 3) 8) p)) (c (expm x 2 p))) (if (= a c) (both x) (both (modulo (* x (expm 2 (/ (- p 1) 4) p)) p))))) (else (let* ((d (let loop ((d 2)) (if (= (jacobi d p) -1) d (loop (+ d 1))))) (s (let loop ((p (- p 1)) (s 0)) (if (odd? p) s (loop (quotient p 2) (+ s 1))))) (t (quotient (- p 1) (expt 2 s))) (big-a (expm a t p)) (big-d (expm d t p)) (m (let loop ((i 0) (m 0)) (cond ((= i s) m) ((= (- p 1) (expm (* big-a (expm big-d m p)) (expt 2 (- s 1 i)) p)) (loop (+ i 1) (+ m (expt 2 i)))) (else (loop (+ i 1) m)))))) (both (modulo (* (expm a (/ (+ t 1) 2) p) (expm big-d (/ m 2) p)) p)))))))))

(define-syntax (with-modulus stx) (syntax-case stx () ((with-modulus e expr ...) (with-syntax ((modulus (datum->syntax-object (syntax with-modulus) 'modulus)) (== (datum->syntax-object (syntax with-modulus) '== )) (+ (datum->syntax-object (syntax with-modulus) '+ )) (- (datum->syntax-object (syntax with-modulus) '- )) (* (datum->syntax-object (syntax with-modulus) '* )) (/ (datum->syntax-object (syntax with-modulus) '/ )) (^ (datum->syntax-object (syntax with-modulus) '^ )) (! (datum->syntax-object (syntax with-modulus) '! )) (sqrt (datum->syntax-object (syntax with-modulus) 'sqrt ))) (syntax (letrec ((fold (lambda (op base xs) (if (null? xs) base (fold op (op base (car xs)) (cdr xs)))))) (let* ((modulus e) (mod (lambda (x) (if (not (integer? x)) (error 'with-modulus "all arguments must be integers") (modulo x modulus)))) (== (lambda (x y) (= (mod x) (mod y)))) (+ (lambda xs (fold (lambda (x y) (mod (+ x (mod y)))) 0 xs))) (- (lambda (x . xs) (if (null? xs) (mod (- 0 x)) (fold (lambda (x y) (mod (- x (mod y)))) x xs)))) (* (lambda xs (fold (lambda (x y) (mod (* x (mod y)))) 1 xs))) (/ (lambda (x . xs) (if (null? xs) (inverse x e) (fold (lambda (x y) (* x (inverse y e))) x xs)))) (^ (lambda (base exp) (expm base exp e))) (! (lambda (n) (mod-fact n modulus))) (sqrt (lambda (x) (mod-sqrt x e)))) expr ...)))))))

That’s a lot of code, but it’s a great illustration of the power of Scheme macros; you should look for opportunities to do things like that in your code, because it makes the code easily re-usable and simplifies the uses of the code. Here is the twin-prime predicate:

(define (twin? m) (with-modulus (* m (+ m 2)) (== (* 4 (+ (! (- m 1)) 1)) (- m))))

> (filter twin? (range 3 1000 2)) ; A001359 (3 5 11 17 29 41 59 71 101 107 137 149 179 191 197 227 239 269 281 311 347 419 431 461 521 569 599 617 641 659 809 821 827 857 881)

You can see the program at https://ideone.com/SZyqeb, but for some reason I can’t make it run with the `syntax-case`

macro; I assure you it works fine with Chez.

Advertisements

@programmingpraxis: Not understanding the filter. According to your calculations 3 is a solution. So substituting it into the filter — 4((m−1)! + 1) == –m (mod m (m+2)) — I get 4((3−1)! + 1) == –3 (mod 3 (3+2)) —> 4(2! + 1) == -3 (mod 3 5) —> 12 == -3 (mod 3 5) —> 12 == -3 3 —> 12 == -9, which fails.

What am I missing?

@bookofstevegraham: The left-hand side of the congruence is 12, which you computed correctly. The right-hand side of the congruence is -3, which you also computed correctly. The modulus of the congruence is m * (m+2) = 3 * 5 = 15. Since 12 and -3 differ by 15, they are congruent modulo 15 and 3 is the base of the twin prime pair (3, 5).

Cannot seem to upload solution (Steve)

Klong version (Try #5)

In Python, using a “twin-wheel”. A twin can only happen at 11, 17 and 29 (repeating every 30).

And just for fun, here a “twin-wheel” for the primes 2,3,5 and 7. Using this wheel only one in seven odd numbers have to be tested.