The Sieve of Eratosthenes / エラトステネスの篩

eratosthenes.scm

;;;;
;;;; Sieve of Eratosthenes
;;;;

(require srfi/2)      ; and-let*
(require srfi/4)      ; homogenous numeric vector
(require srfi/42)     ; list comprehension

;;;;
;;;;
;;;; index    0   1   2   3   4   5   6   7
;;;;        +---+---+---+---+---+---+---+---+
;;;; bitmap | 1 | 1 | 1 | 0 | 1 | 1 | 0 | 1 |
;;;;        +---+---+---+---+---+---+---+---+
;;;; means    3   5   7   9   11  13  15  17
;;;;
;;;;
;;;; The index k represents the odd 2 * k + 3.
;;;;
;;;; The table is the vector of u32 (unsigned 32bit integer).
;;;; The k-th bit of i-th slot means (32*i + k)-th bit of the table.
;;;; The odd indexed by (i, k) is 2*(32*i + k)+1.
;;;;
;;;; To solve the equation: n = 2*(32*i + k) + 1.
;;;;     i = (n - 1) div 64
;;;;     k = ((n - 1) mod 64) / 2
;;;; 


(define prime-table (make-u32vector 1000000 #xffffffff))

(define max-computed 2)  ; even number, the boundary.

(define (naive-prime? p)
  (and-let* ([n (and (odd? p)
                     (/ (- p 1) 2))])
     (let-values ([(i k) (quotient/remainder n 32)])
       (positive? (bitwise-and (u32vector-ref prime-table i)
                               (arithmetic-shift 1 k))))))

(define (clear! n)
  (and-let* ([k (and (odd? n)
                     (/ (- n 1) 2))])
     (let-values ([(i j) (quotient/remainder k 32)])
       (u32vector-set! prime-table
                       i
                       (bitwise-and (u32vector-ref prime-table i)
                                    (bitwise-not (arithmetic-shift 1 j)))))))

(define (sieve p . opt-limit)
  (let ([limit (if (null? opt-limit)
                   (* 64 (u32vector-length prime-table)) ; default
                   (car opt-limit))])
    (when (<= p
              limit
              (* 64 (u32vector-length prime-table)))
      (do-ec [: n (* 3 p) limit (* 2 p)]
             (clear! n)))))

(define (sieve-all)
  (do-ec [: p 3 (integer-sqrt (* 64 (u32vector-length prime-table))) 2]
         [if (naive-prime? p)]
         (sieve p)))

(define (primes n)
  (cons 2
        (list-ec [: p 3 (+ n 1) 2]
                 [if (naive-prime? p)]
                 p)))