Project Euler --- 素数を扱う問題

Project Euler では素数を扱う問題もいくつかみられますが、私が使っている素数関連のスクリプトを紹介してみます。

Petite Chez Scheme を使用していますので、他の処理系を使っている人は、適宜修正をしてください。ちなみに '[...]' は '(...)'に置き換えてください。

まずは、素数の判定法から……

これはオーソドックスに「試し割法」です。（実際に自分が使用しているものは大きな整数にも対応できるように、「Miler-Rabin 法」を組み合わせています）

(define prime?
  (lambda (n)
    (cond
      ([< n 2] #f)
      ([= n 2] #t)
      ([even? n] #f)
      (else
        (let ([max (sqrt n)])
          (let loop ([i 3])
            (cond
              ([> i max] #t)
              ([zero? (remainder n i)] #f)
              (else (loop (+ i 2))))))))))

もし、「n 番目までの素数のリストを求める」のであれば……

(define p-list
  (lambda (n)
    (let loop ([i 2] [a 3] [result '(2)])
      (cond
        ([> i n] (reverse result))
        ([prime? a] (loop (+ i 1) (+ a 2) (cons a result)))
        (else (loop i (+ a 2) result))))))

次は「n 以下の素数のリストを求める」場合です。

"prime?" を使ってもいいのですが、「エラトステネスの篩」を使うとかなり高速になります。

私が使っているものは Ruby の sample に含まれている "sieve.rb" を改良したものです。

(define sequence
  (lambda (start end . step)
    (let* ([step (if [null? step] 1 (car step))]
           [size (+ 1 (quotient (- end start) step))])
      (map (lambda (x) (+ (* x step) start)) (iota size)))))

(define prime-list
  (lambda (n)
    (cond
      ([= n 1] '())
      ([= n 2] '(2))
      (else
        ;; エラトステネスの篩
        (let* ([p-vct (list->vector (cons 0 (sequence 3 n 2)))]
               ;; p-vct = #(0 3 5 7 ...)
               [size (vector-length p-vct)]
               [max (sqrt n)])
          (let loop1 ([i 0])
            (let ([j (vector-ref p-vct i)])
              (cond
                ([zero? j] (loop1 (+ i 1)))
                ([> j max] (cons 2 (remv 0 (vector->list p-vct))))
                (else
                  (let loop2 ([k (* 2 i (+ i 1))])
                    (if (>= k size)
                        (loop1 (+ i 1))
                        (begin (vector-set! p-vct k 0)
                               (loop2 (+ j k))))))))))))))

ちなみに、私のマシン (T1300, 1.66GHz) で実行すると

> (time (car (prime-list 10000)))
(time (car (prime-list 10000)))
    no collections
    3 ms elapsed cpu time
    3 ms elapsed real time
    154832 bytes allocated

という結果が出ました。結構速いと思いません?