Gaucheでアトキンの篩で素数リスト

  • 内包表記、百花繚乱。
  • ナイーブな実装なので、がんばればもっと速くなる。参考文献をどうぞ。
  • うっかり list:ec とtypoして、気づくのに20分ほど無駄にしたのはナイショだ。
(use srfi-42)

(define-syntax vector-flip!
  (syntax-rules ()
    ((_ vec i)
     (vector-set! vec i (not (vector-ref vec i))))))

(define (atkin-sieve limit)
  (let ((is-prime (make-vector (+ limit 1) #f))
        (sqrt-limit (ceiling->exact (sqrt limit))))
    (do-ec (: x 1 sqrt-limit) (: y 1 sqrt-limit)
           (begin
             (let1 n (+ (* 4 (expt x 2)) (expt y 2))
               (when (and (<= n limit) 
                          (or (= 1 (modulo n 12))
                              (= 5 (modulo n 12))))
                 (vector-flip! is-prime n)))
             (let1 n (+ (* 3 (expt x 2)) (expt y 2))
               (when (and (<= n limit) 
                          (= 7 (modulo n 12)))
                 (vector-flip! is-prime n)))
             (let1 n (- (* 3 (expt x 2)) (expt y 2))
               (when (and (> x y)
                          (<= n limit)
                          (= 11 (modulo n 12)))
                 (vector-flip! is-prime n)))))
    (do-ec (: n 5 sqrt-limit)
           (if (vector-ref is-prime n))
           (:range x (expt n 2) limit (expt n 2))
           (vector-set! is-prime x #f))
    (cons 2 (cons 3 (list-ec (: x 5 limit)
                             (if (vector-ref is-prime x))
                             x)))))

(print (atkin-sieve 100))
=> (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)