; twin primes
; primes n -- list of primes not greater than n in ascending order
(define (primes n) ; assumes n is an integer greater than one
(let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
(let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
(cond ((< n (* p p))
(do ((i i (+ i 1)) (p p (+ p 2))
(ps ps (if (vector-ref bits i) (cons p ps) ps)))
((= i len) (reverse ps))))
((vector-ref bits i)
(do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
(vector-set! bits j #f)))
(else (loop (+ i 1) (+ p 2) ps))))))
; inverse x n -- modular inverse: y such that x*y == 1 (mod n)
(define (inverse x n)
(let loop ((x (modulo x n)) (a 1))
(cond ((zero? x) (error 'inverse "division by zero"))
((= x 1) a)
(else (let ((q (- (quotient n x))))
(loop (+ n (* q x)) (modulo (* q a) n)))))))
(define (twins n)
(let* ((len (quotient n 6))
(bits (make-vector len #t))
(ps (cddr (primes (inexact->exact (ceiling (sqrt n)))))))
(do ((ps ps (cdr ps))) ((null? ps))
(let* ((x (inverse 6 (car ps)))
(x (+ x (if (= (car ps) (- (* 6 x) 1)) (car ps) 0))))
(do ((i (- x 1) (+ i (car ps)))) ((<= len i))
(vector-set! bits i #f)))
(let* ((x (inverse -6 (car ps)))
(x (+ x (if (= (car ps) (+ (* 6 x) 1)) (car ps) 0))))
(do ((j (- x 1) (+ j (car ps)))) ((<= len j))
(vector-set! bits j #f))))
(let loop ((t 0) (ts (list 4)))
(cond ((= len t) (reverse ts))
((vector-ref bits t)
(loop (+ t 1) (cons (* 6 (+ t 1)) ts)))
(else (loop (+ t 1) ts))))))
(display (twins 1000)) (newline)
(display (length (twins 1000000))) (newline)