; the sum of the first billion primes

(define (last-pair xs)
  (if (null? (cdr xs)) xs
    (last-pair (cdr xs))))

(define (sum-primes gen n)
  (let loop ((p (gen)) (n n) (s 0))
    (if (zero? n) s
      (loop (gen) (- n 1) (+ s p)))))

(define (gen-brute)
  (define (prime? n)
    (if (even? n) (= n 2)
      (let loop ((d 3))
        (cond ((< n (* d d)) #t)
              ((zero? (modulo n d)) #f)
              (else (loop (+ d 2)))))))
  (let ((next 2))
    (lambda ()
      (let loop ((n (+ next 1)))
        (if (prime? n)
            (let ((p next))
              (set! next n)
              p)
            (loop (+ n 1)))))))

(time (display "brute: ")
      (display (sum-primes (gen-brute) 10000))
      (newline))
(newline)

(define (gen-wheel)
  (define (prime? n) ; n < 2^32
    (define (expm b e m)
      (define (times x y) (modulo (* x y) m))
      (let loop ((b b) (e e) (r 1))
        (if (zero? e) r
          (loop (times b b) (quotient e 2)
                (if (odd? e) (times b r) r)))))
    (define (witness? a n) ; composite if #t
      (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
          ((odd? d)
            (let ((t (expm a d n)))
              (if (or (= t 1) (= t (- n 1))) #f
                (do ((s (- s 1) (- s 1))
                     (t (expm t 2 n) (expm t 2 n)))
                    ((or (zero? s) (= t (- n 1)))
                      (zero? s))))))))
    (cond ((zero? (modulo n 2)) (= n 2))
          ((zero? (modulo n 7)) (= n 7))
          ((zero? (modulo n 61)) (= n 61))
          (else (not (or (witness? 2 n)
            (witness? 7 n) (witness? 61 n))))))
  (let ((wheel (cons 1 (cons 2 (cons 2 (cons 4
          (let ((xs (list 2 4 2 4 6 2 6 4 2 4 6 6 2 6
                 4 2 6 4 6 8 4 2 4 2 4 8 6 4 6 2 4 6
                 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10)))
            (set-cdr! (last-pair xs) xs) xs))))))
        (next 2))
    (lambda ()
      (let loop ((n (+ next (car wheel))))
        (set! wheel (cdr wheel))
        (if (not (prime? n)) (loop (+ n (car wheel)))
          (let ((result next)) (set! next n) result))))))

(time (display "wheel: ")
      (display (sum-primes (gen-wheel) 10000))
      (newline))
(newline)

(define (gen-oneill)
  (define (pq-rank pq) (vector-ref pq 0))
  (define (pq-item pq) (vector-ref pq 1))
  (define (pq-lkid pq) (vector-ref pq 2))
  (define (pq-rkid pq) (vector-ref pq 3))
  (define pq-empty (vector 0 'pq-empty 'pq-empty 'pq-empty))
  (define (pq-empty? pq) (eqv? pq pq-empty))
  (define (pq-merge lt? p1 p2)
    (define (pq-swap item lkid rkid)
      (if (< (pq-rank rkid) (pq-rank lkid))
          (vector (+ (pq-rank rkid) 1) item lkid rkid)
          (vector (+ (pq-rank lkid) 1) item rkid lkid)))
    (cond ((pq-empty? p1) p2)
          ((pq-empty? p2) p1)
          ((lt? (pq-item p2) (pq-item p1))
            (pq-swap (pq-item p2) (pq-lkid p2)
                     (pq-merge lt? p1 (pq-rkid p2))))
          (else (pq-swap (pq-item p1) (pq-lkid p1)
                         (pq-merge lt? (pq-rkid p1) p2)))))
  (define (pq-insert lt? x pq)
    (pq-merge lt? (vector 1 x pq-empty pq-empty) pq))
  (define (pq-first pq)
    (if (pq-empty? pq)
        (error 'pq-first "empty priority queue")
        (pq-item pq)))
  (define (pq-rest lt? pq)
    (if (pq-empty? pq)
        (error 'pq-rest "empty priority queue")
        (pq-merge lt? (pq-lkid pq) (pq-rkid pq))))
  (define (lt? a b)
    (or (< (car a) (car b))
        (and (= (car a) (car b))
             (< (cdr a) (cdr b)))))
  (let ((first? 2) (second? 3) (prev 3)
        (pq (pq-insert lt? (cons 9 6) pq-empty)))
    (lambda ()
      (if first? (begin (set! first? #f) 2)
        (if second? (begin (set! second? #f) 3)
          (let loop1 ((n (+ prev 2)))
            (if (< n (car (pq-first pq)))
                (let ((c (* n n)) (s (+ n n)))
                  (set! prev n)
                  (set! pq (pq-insert lt? (cons c s) pq))
                  n)
                (let loop2 ()
                  (if (< n (car (pq-first pq)))
                      (loop1 (+ n 2))
                      (let* ((c (car (pq-first pq)))
                             (s (cdr (pq-first pq)))
                             (cs (cons (+ c s) s)))
                        (set! pq (pq-insert lt? cs (pq-rest lt? pq)))
                        (loop2)))))))))))

(time (display "oneill: ")
      (display (sum-primes (gen-oneill) 10000))
      (newline))
(newline)

(define (gen-sieve)
  (define (prime? n)
    (if (even? n) (= n 2)
      (let loop ((d 3))
        (cond ((< n (* d d)) #t)
              ((zero? (modulo n d)) #f)
              (else (loop (+ d 2)))))))
  (define (init)
    (let ((ps (list)) (qs (list))
          (sv (make-vector 50000 #t)))
      (do ((p 3 (+ p 2))) ((< 100000 (* p p)))
        (when (prime? p) (set! ps (cons p ps))))
      (set! qs (map (lambda (p) (+ (* 1/2 (- p 1)) p)) ps))
      (do ((ps ps (cdr ps)) (qs qs (cdr qs)))
          ((null? ps))
        (let ((p (car ps)) (q (car qs)))
          (do ((i q (+ i p))) ((<= 50000 i))
            (vector-set! sv i #f))))
      (values ps qs sv)))
  (define (advance ps qs sv bot)
    (do ((i 0 (+ i 1))) ((= i 50000))
      (vector-set! sv i #t))
    (set! qs (map (lambda (p q) (modulo (- q 50000) p)) ps qs))
    (do ((p (+ (car ps) 2) (+ p 2)))
        ((< (+ bot 100000) (* p p)))
      (when (prime? p)
        (set! ps (cons p ps))
        (set! qs (cons (/ (- (* p p) bot 1) 2) qs))))
    (do ((ps ps (cdr ps)) (qs qs (cdr qs)))
        ((null? ps))
      (let ((p (car ps)) (q (car qs)))
        (do ((i q (+ i p))) ((<= 50000 i))
          (vector-set! sv i #f))))
    (values ps qs sv))
  (let-values (((ps qs sv) (init))
        ((bot next-i next-p) (values 0 1 2)))
    (lambda ()
      (let ((result next-p))
        (let loop ((i next-i))
          (cond ((= i 50000)
                  (set! bot (+ bot 100000))
                  (let-values (((p q s)
                      (advance ps qs sv bot)))
                    (set! ps p) (set! qs q)
                    (set! sv s) (set! next-i 0))
                  (loop 0))
                ((vector-ref sv i)
                  (set! next-i (+ i 1))
                  (set! next-p (+ bot i i 1))
                   result)
                (else (loop (+ i 1)))))))))

(time (display "sieve: ")
      (display (sum-primes (gen-sieve) 10000))
      (newline))

brute: 496165411
cpu time: 90 real time: 438 gc time: 0

wheel: 496165411
cpu time: 150 real time: 865 gc time: 20

oneill: 496165411
cpu time: 230 real time: 4648 gc time: 0

sieve: 496165411
cpu time: 0 real time: 84 gc time: 0