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.
Programming Praxis 
@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)
3+!997 [3 4 5 6 7 8 9 ... 999] {&/x!:\2+!_x^1%2}'3+!997 :" Which of the numbers are prime?" [1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 ... 0] a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Only keep the primes" [3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997] ,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Get each pair" [[3 5] [5 7] [7 11] [11 13] [13 17] [17 19] [19 23] [23 29] ... [991 997]] {2=(x@1)-x@0}'b::,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Which pairs have a difference of 2" [1 1 0 1 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] &{2=(x@1)-x@0}'b::,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Only keep the ones with the proper difference - these are indices" [0 1 3 5 8 11 15 18 24 26 31 33 39 41 43 47 50 55 58 62 67 79 81 87 96 102 107 111 114 118 138 140 142 146 150] b::b@&{2=(x@1)-x@0}'b::,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" These are the pairs" [[3 5] [5 7] [11 13] [17 19] [29 31] [41 43] [59 61] [71 73] [101 103] [107 109] [137 139] [149 151] [179 181] [191 193] [197 199] [227 229] [239 241] [269 271] [281 283] [311 313] [347 349] [419 421] [431 433] [461 463] [521 523] [569 571] [599 601] [617 619] [641 643] [659 661] [809 811] [821 823] [827 829] [857 859] [881 883]] mod::{:[x>0; x!y; (x+y)!y]} :dyad mod(12; 15) 12 mod(-3; 15) 12 {[m1 m2]; m1::x@0; m2::x@1; .p($x," = ",$mod(4*1+*/1+!m1-1; m1*m2)=mod(-m1; m1*m2))}'b;"" :" All pairs are twin primes" [3 5 = 1] [5 7 = 1] [11 13 = 1] [17 19 = 1] [29 31 = 1] [41 43 = 1] [59 61 = 1] [71 73 = 1] [101 103 = 1] [107 109 = 1] [137 139 = 1] [149 151 = 1] [179 181 = 1] [191 193 = 1] [197 199 = 1] [227 229 = 1] [239 241 = 1] [269 271 = 1] [281 283 = 1] [311 313 = 1] [347 349 = 1] [419 421 = 1] [431 433 = 1] [461 463 = 1] [521 523 = 1] [569 571 = 1] [599 601 = 1] [617 619 = 1] [641 643 = 1] [659 661 = 1] [809 811 = 1] [821 823 = 1] [827 829 = 1] [857 859 = 1] [881 883 = 1] ""In Python, using a “twin-wheel”. A twin can only happen at 11, 17 and 29 (repeating every 30).
from itertools import accumulate, cycle, chain def fac_mod(n, m): """calculates n! (mod m)""" if n >= m: return 0 fac = 1 for i in range(2, n+1): fac = fac * i % m if fac == 0: break return fac def is_twin(m): M = m * (m + 2) return (4 * (fac_mod(m-1, M) + 1) + m) % M == 0 def gtwin(limit=None): for cand in accumulate(chain((3, 2, 6), cycle((6, 12, 12)))): if limit is not None and cand > limit: break if is_twin(cand): yield cand print(list(gtwin(1000)))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.
def gtwin2(limit=None): for cand in accumulate(chain((3, 2, 6, 6), cycle((12, 12, 18, 12, 30, 6, 30, 12, 18, 12, 12, 6, 12, 12, 6)))): if limit is not None and cand > limit: break if is_twin(cand): yield cand