codepad

Language:

; 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 (ceiling (/ n 6)))
         (b1 (make-vector len #t))
         (b2 (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! b1 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! b2 j #f))))
    (let loop ((t 0) (ts (list 4)))
      (cond ((= len t) (reverse ts))
            ((and (vector-ref b1 t) (vector-ref b2 t))
              (loop (+ t 1) (cons (* 6 (+ t 1)) ts)))
            (else (loop (+ t 1) ts))))))

(display (twins 1000)) (newline)
(display (length (twins 1000000))) (newline)

Private [?]

Run code