エラトステネスの無限の篩
♯エラトステネスの篩とは
一般にエラトステネスの篩（ふるい）というと求める素数の上限を決めてそれ以下の整数（２以上）から素数の倍数を消すことで素数を篩い出すアルゴリズムです。詳細はwikipediaに分かり易く書かかれていますのでそちらをどうぞ。
このエラトステネスの篩に上限をあたえず無限素数リスト作ってみようという試みです。
Clojure は遅延シーケンスにより無限長のリストを扱うことができます。無限素数リストをつくる例はWEB上でいくつも見られるのですが、エラトステネスの篩で書いているのは少ないのではないかと思います。
エラトステネスの篩といいつつも実は試し割りで素数を求めている例をちらほらと見掛けます。
剰余を求めるのは試し割りです。N 以下の素数で割ってみて余りが非ゼロのものを残す、というアルゴリズムはエラトステネスの篩ではないと思います。
♯基本方針
基本方針としては、
- ２以上の整数リスト nums を作る
- nums から「２の倍数リスト」を除去する。
- nums から「３の倍数リスト」を除去する。
- nums から「５の倍数リスト」を除去する。
- 以下同様に nums から「p の倍数リスト」を除去する操作を繰り返す。（pは素数）。
これで nums には素数だけが残ります。
「p の倍数リスト」は p の倍数を要素とする無限リストです。p そのものは含まないものとします。
「p の倍数リストを除去する」とは、p の倍数リストに現れる要素を nums からすべて除去することを意味します。
nums も 倍数リストも、どちらも無限リストなので無限リスト同士の演算になります。
♯課題
注意しなければならない点は、無限の素数リストを作る場合、２の倍数、３の倍数、５の倍数....と、nums から除去しなければならない倍数リストも無限にあるという事です。
「無限リストから無限個の要素を除去する操作」を「無限個の倍数リストを対象に」行うことになります。
遅延評価と計算範囲の限定を行っていかないと、有限時間で処理できないプログラムになってしまいます。
♯パーツ作り
まずは部品作りから。
(defn- diff-seq [s1 s2] (let [x1 (first s1), x2 (first s2)] (cond (= x1 x2) (recur (rest s1) (rest s2)) (> x1 x2) (recur s1 (drop-while (partial > x1) s2)) (< x1 x2) (let [[s11 s12] (split-with (partial > x2) s1)] (lazy-cat s11 (diff-seq s12 s2)))))) (defn- muls [n] (rest (iterate (partial + n) n))) (defn- sieve [nums n] (diff-seq nums (muls n)))
diff-seq は 第１引数の無限リストから、第２引数の無限リストにある要素を除去する関数です。ざっと見たところ標準モジュールに該当関数が無い*1ようなので作りました*2。
muls は n を含まない n の倍数リストを生成する関数。
sieve は「無限リスト nums」から「n 以外の n の倍数」を除去する関数です。
;; sieve 使用例 (let [nums1 (iterate inc 2) ;=> ２以上の整数 nums2 (sieve nums1 2) ;=> ２以外の２の倍数除去 nums3 (sieve nums2 3) ;=> ３以外の３の倍数除去 nums4 (sieve nums3 5)] ;=> ５以外の５の倍数除去 (take 20 nums4) ;=> (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 49 53 59 61 67) )
7 の倍数を除去していないのでリストには合成数 49 が残っています。
♯素数リストを作り続ける仕組み
上限のあるエラトステネスの篩では、上限値の平方根以下の素数について篩いを行えばよい事になっています*3。
「上限値の平方根以下」などというといかにも有限の考え方に思えますが、裏を返せば「素数 p まで篩をかけると、p の２乗未満の数まで素数が確定している」ともいえます。
この事実を使うとリストの素数確定部分を取得しながら、nums を篩にかけ続けることができます。
具体的な手順は次の通り。
- 【1-a】２以上の整数リスト nums から ２以外の２の倍数を除去。
- 【1-b】 リスト nums を値 2^2(=4) を境界とし二つに分割する。分割した最初の部分（有限長）は素数確定リストであり、二つ目の部分（無限長）は素数未確定部分である。前者 を prms 後者を rest-nums とする。
- 【2-a】 [1-b] で得た rest-nums を次の整数リスト nums とし、nums から ３以外の３の倍数を除去する。
- 【2-b】 再び [1-b] の操作を行う。ただし境界値は 3^2(=9) とします。
- 【3-a】 [2-b] で得た rest-nums を次の整数リスト nums として nums から ５以外の５の倍数を除去する。
- 【3-b】 再び [1-b] の操作を行う。境界値は 5^2(=25) とする。
以降 a,b を繰返して素数確定リストを切り出して行きます。
♯素数リスト生成関数
Clojure のコードにするとこんな感じです。
(defn primes-sieve
([]
(primes-sieve (iterate inc 2) '(2)))
([nums prms-que]
(let [p (first prms-que)
nums (sieve nums p)
[prms rest-nums] (split-with #(> (* p p) %) nums)
prms-que (lazy-cat (rest prms-que) prms)]
(lazy-cat prms (primes-sieve rest-nums prms-que)))))利用者は primes-sieve を無引数で評価します。すると規定値による２引数評価がおこなわれます。
２引数でこの関数を評価すると、終了判定のない再帰呼出しに入ります。ただし再帰呼出は lazy-cat 式の中で行うので遅延評価になります。そのお陰でビジーループになることもなく、また非末尾再帰であるにもかかわらず再帰によるスタック溢れを起すこともありません。
prms-que は確定済み素数を追加してゆく待ち行列リストです。倍数リストはこの prms-que の先頭から取り出した素数 p を元に作成していきます。
関数全体をトップダウンで眺めると、この関数は lazy-cat によって連結された遅延リストを返します。
遅延リストは 最初に素数リスト prms から要素を一つずつ返し、prms の要素が無くなると、primes-sieve の再帰呼出しを行い、新たな prms を得ます。これを繰り返すことで無限に素数を返し続けます。
♯おわりに
以上が Clojure による「エラトステネスの無限の篩」です。剰余計算や除算は一切行っていません。
とはいえ、割り算以上に計算コストがかかりそうなリスト操作を行っていますので、パフォーマンスはあまり保証できないかもしれません。
実用的かと言われると甚だ疑問ですが、でも「無限リストに無限個の無限リストを適用する」なんてちょっと凄そうじゃないですか？
♯追記(2011/9/30)
muls の作る「p の倍数リスト」は「p を除く p の倍数」ではなく「p の二乗以上の p の倍数」でいいですね。p を小さい方から順に処理するなら「p の二乗未満の p の倍数」は既にみんな除かれてますから。
といっても、このプログラムだと処理コストへの影響はほとんど無いですけど。
♯修正(2011/10/1)
diff-seq の末尾再帰箇所を recur に書き換え。
♯追記(2011/10/1)
続編を書きました。
