-
Notifications
You must be signed in to change notification settings - Fork 0
/
atkin-sieve.scm
70 lines (43 loc) · 1.56 KB
/
atkin-sieve.scm
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
;; atkin-sieve.scm
;;
;; determine primes up to given limit by Atkin/Bernstein
;; algorithm
;;
;; FIXME: It's not any faster than Sieve of Eratosthenes (sieve-primes.scm).
;; FIXME: How come?
(define primes-bitvector #f)
(define (prime? n)
(and (> n 1)
(< n (bitvector-length primes-bitvector))
(bitvector-ref primes-bitvector n)))
(define (atkin-sieve limit)
(define (flip-prime n)
(bitvector-set! primes-bitvector n (not (bitvector-ref primes-bitvector n))))
(define (not-a-prime n)
(bitvector-set! primes-bitvector n #f))
(set! primes-bitvector (make-bitvector (+ limit 1) #f))
(for-each (lambda (n) (bitvector-set! primes-bitvector n #t)) '(1 2 3 4))
(let ((sqr_lim (inexact->exact (floor (expt limit 0.5))))
(x2 0)
(y2 0)
(n 0))
(do ((x 1 (+ x 1))) ((> x sqr_lim))
(set! x2 (* x x))
(do ((y 1 (+ y 1))) ((> y sqr_lim))
(set! y2 (* y y))
(set! n (+ (* 4 x2) y2))
(if (and (<= n limit)
(or (= 1 (remainder n 12))
(= 5 (remainder n 12))))
(flip-prime n))
(set! n (+ (* 3 x2) y2))
(if (and (<= n limit) (= 7 (remainder n 12)))
(flip-prime n))
(set! n (- (* 3 x2) y2))
(if (and (> x y) (<= n limit) (= 11 (remainder n 12)))
(flip-prime n))))
(do ((n 1 (+ n 1))) ((> n sqr_lim))
(if (prime? n)
(do ((nn (* n n) (+ nn (* n n)))) ((> nn limit))
(not-a-prime nn))))))
;; end of file