Skip to content

Instantly share code, notes, and snippets.

@sasaki-shigeo
Created September 16, 2020 18:48
Show Gist options
  • Save sasaki-shigeo/54a57a7e30ebd90c6eeda1b2d3b3cb88 to your computer and use it in GitHub Desktop.
Save sasaki-shigeo/54a57a7e30ebd90c6eeda1b2d3b3cb88 to your computer and use it in GitHub Desktop.
The Sieve of Eratosthenes / エラトステネスの篩
;;;;
;;;; 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)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment