（これは Lisp Advent Calendar 2012 の 21 日目の記事です）
遅延リスト等を用いた無限素数列の実装には、素数を試し割りによって列挙する事が多いようです。 試し割りは簡潔に書けるので単に無限素数列を実装したいだけなら良いのですが、Project Euler などで実用するには速度が足りない場合があります。 ある程度の範囲以上に存在する素数を列挙するならば、試し割りよりもエラトステネスの篩を使った方がずっと速いので、このエントリではそれを紹介します。
僕は普段 Scheme でコードを書く時には R6RS しか使わないので Gauche の事はあまり知らないのですが、ここでは R6RS よりも人口の多い Gauche を使う事にしました。 Gauche には遅延シーケンス (Gauche ユーザリファレンス: 6.19 遅延評価) というものがあり、 無限素数列を書いた時にそれを普通のリストと同じように扱う事が出来、また遅延シーケンスを書く為の環境も整っており、何かと都合が良いという理由もあります。
まず、無限素数列の実装に当たって使用するアルゴリズムについて、それぞれの名前と参考となるリンクを以下に挙げます。
無限素数列は遅延シーケンスとして最終的に表現しますが、中身の実装にはジェネレータを使います。 ジェネレータとは引数を持たない手続きで、呼び出し毎に順次結果を返すものをいい、ここでは素数を次々と返すジェネレータを実装します。 なお、ジェネレータと遅延シーケンス、そして試し割りを使った無限素数列は Island Life - 無限素数列、続き の一番下に書かれており、このエントリを書くにあたって大いに参考にしました。 そこに書かれている試し割りによるジェネレータベースの無限素数列を、区間篩によるジェネレータベースの無限素数列に書き換えたものを今から作っていきます。
素数を次々と返すジェレータを prime-generator と名付けます。 prime-generator は今時点で探索し終えている数を何処かに保存しておき、呼び出された時それ以降に存在する素数を順次返すものです。 試し割りを使った prime-generator は、中でループを回して 1 つずつ素数判定を行い、素数があればそれを返す…という処理を書けば良いのですが、これをエラトステネスの篩を使ったものに変更します。 ただし、エラトステネスの篩は N までの素数を全て列挙するアルゴリズムであり、N が決まらない限り使う事は出来ません。 よってここではエラトステネスの篩の変形版である区間篩というアルゴリズムを使います。
区間篩について簡単に。これはある区間だけを篩に掛け素数を列挙するアルゴリズムで、合成数は必ずその数の平方根以下の素因数が含まれるという事を利用したものです。 詳しくは先ほど挙げた Spaghetti Source - 区間ふるい や 素数列挙について - MugiCha を読んで下さい。 区間篩を使えば N という上限が決まらなくとも、区間毎に篩を掛け素数をずらっと列挙し、それを順次返す事で prime-generator を篩によって実装する事が出来ます。
例えば small-primes という区間篩に掛ける為の素数列があるならば、(expt (apply max small-primes) 2) 以下であれば何処の区間であっても、そこを篩に掛けて素数の列挙が可能となります。 Island Life - 無限素数列、続き の prime-generator は、 遅延シーケンスで表現された無限素数列 *primes* を用いて試し割りを行なっており、*primes* と prime-generator は相互参照の形を取っていますが、 区間篩を使ったバージョンを書く場合も同様に、*primes* は prime-generator を参照し、prime-generator は small-primes として *primes* を使って区間篩を行い、素数を列挙しそれを返してやれば良いです。
区間篩を使った prime-generator の処理は次のような流れになります。
大まかな実装はこのようになります。
(use gauche.generator) ;; 篩に掛ける区間のサイズ (define segment-size ...) ;; 区間の開始位置の初期値 (define segment-start ...) ;; N までの素数列を得る (define (get-primes N) ...) ;; start から end までの閉区間を small-primes によって篩に掛け、素数を列挙しそれを返す (define (get-primes/segment start end small-primes) ...) ;; 以下の prime-generator を使った遅延シーケンスによる無限素数列 (define primes (generator->lseq (gappend (list->generator (get-primes (- segment-start 1))) (lambda () (prime-generator))))) ;; 遅延シーケンスで表現された primes を small-primes として使って区間篩を行ない、見つけた素数を順次返す (define prime-generator (generate (lambda (yield) (let recur ([start segment-start] [ps (list)]) (dolist (p ps) (yield p)) (recur (+ start segment-size) (get-primes/segment start (+ start segment-size -1) primes))))))
なお、get-primes/segment は small-primes によって区間篩を行い素数列を返しますが、その区間は small-primes で篩える範囲以下である必要があるので、 prime-generator の中で primes を small-primes として呼び出す時には注意しなければなりません。 その為 segment-size には好きな値を入れられますが、segment-start の値に関しては注意深く計算する必要があります。 僕は数学が苦手で考えるのが面倒になってしまったので、segment-start には (expt segment-size 0.6) の値を入れました。
このコードを動かすには後は segment-size, segment-start に値を入れて、get-primes, get-primes/segment の中身を書くだけで済みます。 ただこれは、単なる篩のコードであまり面白味がないのでここには直接書きません。その代わり gist へのリンクを貼っておきます。 また、普通の区間篩を使って書いたもの以外にも幾つかのバージョンを書いたので、それも合わせて貼っておきます。
最後に、それぞれのバージョンを用いて、(list-ref primes n) とした時に掛かった時間を計測しベンチマークを取ってみました。 区間篩で使う区間のサイズは全て #e1e5 で統一しており、実行環境は C2D E8500(3.16GHz) + Windows 7 x64 + Gauche 0.9.3.3 です。 naive という項目は、上述の shiro さんの書かれた試し割りによるジェネレータベースの無限素数列で測ったものですが、 naive(10^7) については、時間が掛かり過ぎたので途中で打ち切った為に空白となっています。
| 10^4 | 10^5 | 10^6 | 10^7
-----------------------------------+--------+--------+---------+----------
naive | 0.316s | 2.602s | 62.018s |
-----------------------------------+--------+--------+---------+----------
segment-sieve | 0.333s | 1.160s | 11.608s | 153.362s
-----------------------------------+--------+--------+---------+----------
segment-sieve + odd-sieve | 0.272s | 0.884s | 7.268s | 69.354s
-----------------------------------+--------+--------+---------+----------
segment-sieve + wheel-sieve(2,3,5) | 0.252s | 0.830s | 7.168s | 123.917s
奇数限定の篩を含む Wheel factorization には、篩の対象となる数列を省く事でメモリアクセスの回数を削り高速化出来るというメリットがありますが、 その反対に、剰余類の計算や index の計算等によって数値演算の回数が増えて、処理が遅くなるというデメリットもあります。 C や C++ といった高速な計算が可能でメモリアクセスが重荷となっている言語では、 エラトステネスの篩 - peria.jp の eratos003.c, eratos0031.c, eratos0032.c のようにわりと素直な効き目が出るのですが、 Gauche ではメモリアクセスよりも数値演算のコストが重いのかあまり有効ではないようです。
ちなみに遅延リスト + 区間篩 + Wheel factorization を使った無限素数列は、 以前に R6RS で Project Euler の為に書いた事があって (Project-Euler/pelib/primes.sls) こちらのベンチマークも貼っておきます。 実行環境は同じく C2D E8500(3.16GHz) + Windows 7 x64 で Racket v5.3 と Ypsilon 0.9.6-trunk/r503 です。Ypsilon には --heap-limit=1024 のオプションを渡しています。
| 10^4 | 10^5 | 10^6 | 10^7 | 10^8
------------------------------+--------+--------+--------+---------+----------
Racket + wheel-sieve (2) | 0.845s | 0.874s | 1.477s | 8.027s | 124.202s
------------------------------+--------+--------+--------+---------+----------
Racket + wheel-sieve (2,3) | 0.843s | 0.870s | 1.318s | 6.461s | 101.385s
------------------------------+--------+--------+--------+---------+----------
Racket + wheel-sieve (2,3,5) | 0.839s | 0.858s | 1.332s | 6.612s | 117.136s
------------------------------+--------+--------+--------+---------+----------
Ypsilon + wheel-sieve (2) | 0.460s | 0.518s | 3.749s | 44.223s |
------------------------------+--------+--------+--------+---------+----------
Ypsilon + wheel-sieve (2,3) | 0.209s | 0.363s | 2.483s | 29.667s |
------------------------------+--------+--------+--------+---------+----------
Ypsilon + wheel-sieve (2,3,5) | 0.181s | 0.302s | 2.041s | 26.086s |
pelib を使ったベンチマークは以下のコードで行なっています。
#!r6rs (import (rnrs) (pelib)) (define (ref xs n) (if (zero? n) (car xs) (ref (lcdr! xs) (- n 1)))) (print (ref lazy-primes #e1e7))
篩を使って素数を列挙するのは楽しいですね。 まあこれは篩に限った事ではないのでしょうが特に高速化を考えた時はとても面白く、 メモリアクセスの重さに気付かされたり、処理系の特性に気付かされたり、あちらを立てればこちらが立たずで悩んでみたり。 Gauche のジェネレータや遅延シーケンスには殆ど初めて触れたのですがとても使いやすく、ああこんな風にしてライブラリを書けば良いのかと参考になりました。
そんなこんなで以上です。時間が無くなってしまったのでここまでにします。 単にベンチマークを取りたかっただけじゃないの…って感じが自分で 7, 8 割を占めている感じがしますが、このエントリが誰かの役に立てれば幸いです。