; sieve of eratosthenes
(define (primes . args) ; (primes [lo] hi) inclusive at both ends
(let* ((lo (if (null? (cdr args)) 0 (car args)))
(hi (if (null? (cdr args)) (car args) (cadr args))))
(cond ((and (<= lo 100000) (<= hi 1000000)) ; simple sieve
(let* ((len (quotient (- hi 1) 2))
(bits (make-vector len #t)))
(let sieve ((i 0) (p 3) (ps (list 2)))
(cond ((< hi (* p p))
(let sweep ((i i) (p p) (ps ps))
(cond ((= i len)
(let loop ((ps (reverse ps)))
(if (<= lo (car ps)) ps
(loop (cdr ps)))))
((vector-ref bits i)
(sweep (+ i 1) (+ p 2) (cons p ps)))
(else (sweep (+ i 1) (+ p 2) ps)))))
((vector-ref bits i)
(let sift ((j (+ (* 2 i i) (* 6 i) 3)))
(when (< j len)
(vector-set! bits j #f)
(sift (+ j p))))
(sieve (+ i 1) (+ p 2) (cons p ps)))
(else (sieve (+ i 1) (+ p 2) ps))))))
((< lo (sqrt hi))
(let* ((r (inexact->exact (ceiling (sqrt hi))))
(r (if (even? r) r (+ r 1))))
(append (primes lo r) (primes r hi))))
(else ; segmented sieve
(let* ((lo (if (even? lo) lo (- lo 1)))
(b 50000) (bs (make-vector b #t))
(r (inexact->exact (ceiling (sqrt hi))))
(ps (cdr (primes r)))
(qs (map (lambda (p)
(modulo (* -1/2 (+ lo 1 p)) p)) ps))
(zs (list))
(z (lambda (p) (set! zs (cons p zs)))))
(do ((t lo (+ t b b))
(qs qs (map (lambda (p q) (modulo (- q b) p))
ps qs)))
((<= hi t)
(let loop ((zs zs))
(if (<= (car zs) hi) (reverse zs)
(loop (cdr zs)))))
(do ((i 0 (+ i 1))) ((= i b)) (vector-set! bs i #t))
(do ((ps ps (cdr ps)) (qs qs (cdr qs))) ((null? qs))
(do ((j (car qs) (+ j (car ps)))) ((<= b j))
(vector-set! bs j #f)))
(do ((j 0 (+ j 1))) ((= j b))
(if (vector-ref bs j) (z (+ t j j 1))))))))))
(time (display (length (primes #e1e9 (+ #e1e9 2043)))) (newline)) (newline)
(time (display (primes #e1e11 (+ #e1e11 100))) (newline)) (newline)
(time (display (primes #e1e12 (+ #e1e12 100))) (newline)) (newline)