エラトステネスの篩

loop*1を使って、エラトステネスの篩を実装してみたメモ。
以下、処理系にはSBCLのver1.0.54(x86-64bit)を使用。

;; 引数nまでの範囲の素数のシーケンス(ジェネレータ)を作成する
(declaim (inline make-prime-sequence))
(defun make-prime-sequence (n)
  (let ((arr (make-array (1+ n) :element-type 'bit :initial-element 1)))
    (flet ((prime? (i) (= (bit arr i) 1))       
           (not-prime! (i) (setf (bit arr i) 0))) 
      (declare (inline prime? not-prime!))

      (loop:each (lambda (i)
                   (when (prime? i)
                     (loop:each #'not-prime! (loop:from (* i 2) :to n :by i))))
                 (loop:from 2 :to (floor (sqrt n))))
    
      (loop:filter #'prime? (loop:from 2 :to n)))))

;;; 実行例
;; 100以下の素数
(loop:collect (make-prime-sequence 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)

;; 1001から1010番目の素数
(loop:collect (loop:take 10 (loop:drop 1000 (make-prime-sequence 10000000))))
=> (7927 7933 7937 7949 7951 7963 7993 8009 8011 8017)

通常のループ(loopマクロ)を使った場合との速度比較。

;; 比較用に素数の合計値を求める関数を用意
(defun prime-sum1 (n)
  (declare (fixnum n)
           (optimize (speed 3) (safety 0) (debug 0)))
  (loop:sum #'identity (make-prime-sequence n)))

;; 一億以下の素数の合計値
(time (prime-sum1 100000000))
Evaluation took:
  1.675 seconds of real time  ; 1.675秒
  1.676105 seconds of total run time (1.676105 user, 0.000000 system)
  100.06% CPU
  3,342,591,038 processor cycles
  12,500,032 bytes consed
=> 279209790387276

;; loopマクロ版
(defun prime-sum2 (n)
  (declare (fixnum n)
           (optimize (speed 3) (safety 0) (debug 0)))
  (let ((arr (make-array (1+ n) :element-type 'bit :initial-element 1)))
    (flet ((prime? (i) (= (bit arr i) 1))
           (not-prime! (i) (setf (bit arr i) 0)))
      (declare (inline prime? not-prime!))

      (loop FOR i fixnum FROM 2 TO (floor (sqrt n))
            WHEN (prime? i)
        DO
        (loop FOR j fixnum FROM (* i 2) TO n BY i
          DO
          (not-prime! j)))

      (loop WITH sum OF-TYPE (unsigned-byte 64)
            FOR i fixnum FROM 2 TO n
            WHEN (prime? i)
        DO (incf sum i)
        FINALLY (return sum)))))

;; 一億以下の素数の合計値
(time (prime-sum2 100000000))
Evaluation took:
  1.476 seconds of real time  ; 1.476秒
  1.472092 seconds of total run time (1.468092 user, 0.004000 system)
  99.73% CPU
  2,944,592,020 processor cycles
  12,500,032 bytes consed
=> 279209790387276

*1:参照: ループ処理を関数型っぽく書いてみる(2) - sileの日記