続・エラトステネスの無限の篩

本題の前に、前回のコードを少し書き換えた*1ので、まずはそちらから。

♯エラトステネスの無限の篩 Ver.1

アルゴリズム上の変更点は倍数リスト*2を「n を除く n の倍数」ではなく、「n^2 を初期値とする n の倍数」にしただけです。

(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 [[s1a s1b] (split-with (partial > x2) s1)]
                 (lazy-cat s1a (diff-seq s1b s2))))))

(defn- arithmetic-seq
  [start step]
  (iterate (partial + step) start))

(defn- prime-numbers*
  ([nums prms-que]
     (let [p (first prms-que)
           pp (* p p)
           sieved-nums (diff-seq nums (arithmetic-seq pp p))
           [prms rest-nums] (split-with (partial > pp) sieved-nums)
           prms-que (lazy-cat (rest prms-que) prms)]
       (lazy-cat prms (prime-numbers* rest-nums prms-que)))))

(defn prime-numbers
  []
  (prime-numbers* (arithmetic-seq 2 1) '(2)))

※gistはこちら。
前回の muls は名称を変えて arithmetic-seq になりました。英語で等差数列を「arithmetic sequence」と言うそうです。
「n の倍数」は「公差 n の等差数列」ですから目的は一緒です。muls と異なるのは、数列の初期値を与えるようにした点です。
また前回あった sieve は関数そのものを無くしました。初期値を与える関係で arithmetic-seq (と diff-seq) をメインルーチンから直接呼びだしています。
前回のコードでメインのルーチンだった、primes-sieve は、prime-numbers と名を改め、外部公開用と内部処理用の２つに分けました。内部用は prime-numbers* です。*3
以上、基本的にやってることは変っていません。

♯エラトステネスの無限の篩 Ver.2

ここから本題。diff-seq と arithmetic-seq は Ver.1 と全く同じなので省略します。
prms-que に関連した部分が変わっています。

(defn- prime-numbers*
  [nums prms-que]
  (let [p (first prms-que)
        pp (* p p)
        sieved-nums (diff-seq nums (arithmetic-seq pp p))
        [prms rest-nums] (split-with (partial > pp) sieved-nums)]
    (lazy-cat prms (prime-numbers* rest-nums (rest prms-que)))))

(defn prime-numbers
  []
  (lazy-seq
   (prime-numbers* (arithmetic-seq 2 1)
                   (cons 2 (drop 1 (prime-numbers))))))

※gistはこちら
prms-que に prms を追加するコードが消えました。
では prms-que への素数の補充はどこでするか？
答えは「補充する必要が無くなった」です。
上のコードの一番最後の関数 prime-numbers をアイディアの本質だけ見えるようにすると次のようになります。

;; prime-numbers (Ver.2) の擬似コード
(defn prime-numbers
  []
  (prime-numbers* (arithmetic-seq 2 1)
                  (prime-numbers)))

prime-numbers* 呼び出しをしています。第１引数には「２以上の整数*4」、第２引数に無限素数リストを与えてます。初めから無限素数リストを与えているので、prime-numbers* 内では素数を追加する必要もないわけです。
このコードの面白いところは、prime-numbers の中で、素数列生成のために素数列が必要になり、その素数列生成に prime-numbers を使っている点です。
自分の生成物（素数列）を作るために、自分の生成物を使っているんです。目的と手段が再帰しています。


さて、上で書いた擬似コードでは実際には期待通りに動きません。prime-numbers の自己再帰呼出が何も値を返さないまま、次の自己再帰をするので無限の再帰に陥いってしまいます。そこで素数列の最初の1項(=2)だけは計算させずに規定値として書いてあげます。

(defn prime-numbers
  []
  (prime-numbers* (arithmetic-seq 2 1)
                  (cons 2 (drop 1 (prime-numbers)))))

第１項の 2 だけは prime-numbers に期待せず、cons しています。そのかわり prime-numbers の返すリストの最初の１項は drop します。これで辻褄があいます。
しかし、それだけでは遅延評価にならないので全体を lazy-seq にします。これで完成です。

(defn prime-numbers
  []
  (lazy-seq
   (prime-numbers* (arithmetic-seq 2 1)
                   (cons 2 (drop 1 (prime-numbers))))))

♯prime-numbers のバリエーション

lazy-seq ではなく lazy-cat を使った書き方もできます。

(defn prime-numbers
  []
  (prime-numbers* (arithmetic-seq 2 1)
                  (lazy-cat [2] (drop 1 (prime-numbers)))))

(defn prime-numbers
  []
  (prime-numbers* (arithmetic-seq 2 1)
                  (lazy-cat [2] (rest (prime-numbers)))))

ただし、次のような書き方はダメ。この場合 rest ではオーバーフローになります。

;; NG (Stack Overflow)
(defn prime-numbers
  []
  (lazy-seq
   (prime-numbers* (arithmetic-seq 2 1)
                   (cons 2 (rest (prime-numbers))))))

遅延評価になってないということでしょうか。何故これだけダメなのかは、まだ追っていません。いずれ気が向いたら....

♯おわりに

Ver.1 のコードで一番気になっていたのが、待ち行列リスト prms-que でした。このリストは要素を取り出すペースよりも、追加するペースの方が早く、どんどん大きくなってゆきます。使用している遅延リストはどれも余再帰、あるいはそれと同等のもので、不要になった要素はどんどん解放できる*5のですが、prms-que だけはデータを貯め込む形になっていました。*6 *7 *8
そんな気持ちの悪さを解消しようと今回の Ver.2 を作ったのですが、思わぬ不思議なプログラムになってしまいました。「無限素数リストを作るために、無限素数リストを使う」という冗談のようなプログラムが期待どおりに動いています。

♯補足

Ver.2 の方がメモリ使用効率がよい、と断定はできません。Ver.2 は prms-que に素数を pool することはありませんが、その代わり素数を得るたび*9に新たな素数リストを作ります。新たな素数リストを作るということはそれに伴って倍数リストも作り直しますので、トータルでは冗長な処理やインスタンスが積み重なって行きます。*10

♯追記

素数確定リスト prms を切り出す境界条件は実はもう少し大きく取れます。n 番目の素数を P(n) とし P(0)〜P(n) まで篩にかけたとき。これまではリストの中の「P(n)の二乗未満」を境界にしていましたが、実際は「P(n+1)の二乗未満」の値まで素数確定です。
この点を反映した Ver.2 のソースコードを修正しました。ここにあります。

♯さらに追記

ちょっとだけ無駄があったので修正。コードはこちら。

*1: @omasanori さんのコードを一部参考にさせていただきました。

*2:本当はClojureではリストではなくシーケンスというべきなのですが、この記事ではリストで統一します。

*3:この様に外部公開目的の関数と内部処理用の関数を分け、公開用の関数名の最後に * をつけたものを内部用関数名とするのは Clojure の習慣です。

*4:初期値２公差１の等差数列としてリストを作成。

*5:理論上の話です。実際に解放されるタイミングは処理系の都合なので分かりません。

*6:とは言え、全てのデータを貯め込むわけではなく、少しずつスローペースで消費はしています。

*7:取得する素数の値の１０％程の数の要素が prms-que に貯まります。１００万台の素数を得ると、１０万要素くらい貯まります。

*8:prms-que も、そこへ追加される prms も遅延リストなので全部の要素がメモリに置かれるのかどうかは分かりません。

*9:厳密には「prms-queで新たに要素を現実化するたび」にです。大きな素数リストを作れば作るほど全体における再帰の頻度（新たな素数を現実化する頻度）は希薄になって行きます。

*10:確認したところ、１００万以下の素数を列挙するのに行われた prime-numbers の再帰回数はたったの３回でした。作られる素数リストの数は想像していたよりはずっと少いようです。