codepad

Language:

; 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)

Private [?]

Run code