2011/09/29

無限素数列、続き

「「よりClojureらしい素数列」をGaucheで」のエントリでは、Gaucheのlazy sequenceの例を簡単に示すのに良いと思って単純な素数列生成のClojureコードをそのまま移植したのだけれど、何やら素数列生成が盛り上がっているようなのでもうちょい続けてみる。なおこのエントリのコードも、実行にはGaucheの開発版HEADが必要。

前エントリでは、無限の自然数列から、素数であるものだけを抽出して素数列を求めていたが、nが素数であるかの判定は極めてナイーブで、2から(+ 1 (/ n 2))までの自然数でnを割ってみる、というものだった。

(define (prime? n)
  (not (any zero? (lmap (pa$ mod n) (lrange 2 (+ 1 (/ n 2)))))))

でもこれはかなり無駄なことをしている。2と3で割れなければ6でも割れるわけないのだから、(ceiling √n) 以下の素数についてだけ試せばいい。

今、仮に、*primes* に無限素数列が既に束縛されているとしよう。 すると、prime? の定義は次のように書ける。

(define (prime? n)
  (not (any zero? (lmap (pa$ mod n) (take-while (^k (<= (* k k) n)) *primes*)))))

ある数k (kは3以上の奇数とする) が与えられた時、そこから先の素数列は次のように書ける。

(define (primes-from k)
  (if (prime? k)
    (lcons k (primes-from (+ k 2)))
    (primes-from (+ k 2))))

(+ k 2) してるのは、3以上の素数はすべて奇数なので偶数を試す必要はないからだ。

ここで、lcons はlazyなcons。Clojureなら (lazy-seq (cons n (recur (+ k 2)))) などと書くところだが、Gaucheの場合lconsというプリミティブでlazyなペアを作る。ちょっと注意が必要なのは、lconsはcdr部分についてだけ評価を遅延するということ。つまり (cons k (delay (loop (+ k 2)))) に近い。ただしdelayと違って、lconsで作られたlazyなペアのcdr部は必要になったら自動的に評価されるので、いちいちforce する必要はない。

primes-fromを使えば、*primes*を定義できる。primes-from の中でprime?を使ってて、その中で*primes*を参照しているから循環定義になりそうなものだが、 primes-from でk以上の素数を計算する時には、(ceiling √k) 以下の素数しか参照しない。つまりある程度kが大きければ、既に計算済みの部分しか使わないようにできる。

(define *primes* (list* 2 3 (lcons 5 (primes-from 7))))

(primes-from 7) を計算するには、*primes* の頭から2と3だけを試せば良いので、2と3をあらかじめリストの頭にくっつけて置けば良いというわけだ。全部まとめるとこうなる。

(define (prime? n)
  (not (any zero? (lmap (pa$ mod n) (take-while (^k (<= (* k k) n)) *primes*)))))

(define (primes-from k)
  (if (prime? k)
    (lcons k (primes-from (+ k 2)))
    (primes-from (+ k 2))))

(define *primes* (list* 2 3 (lcons 5 (primes-from 7))))

(lazyな計算に慣れてる人は、(cons 2 (lcons 3 (primes-from 5))) でも良いんじゃないかと思うかもしれない。 だが、前エントリ でちょっと触れたように、Gaucheのlazy sequenceは一つ先の要素を計算してしまう。(primes-from 5) を計算するには2と3を試す必要があり、lazy pairのcarである3を取り出そうとした瞬間、Gaucheはそのcdrの要素を一つ計算しようとしてしまうので、計算が終わらなくなる。「一つ先読み」が落とし穴になりそうだというのは、こういうケースがあるから。)

なお、Clojureのcontribにも無限素数列のコードがあるが(clojure.contrib.lazy-seqs/primes)、そちらはもうちょい賢いことをやっている。

(defvar primes
  (concat
   [2 3 5 7]
   (lazy-seq
    (let [primes-from
          (fn primes-from [n [f & r]]
              (if (some #(zero? (rem n %))
                        (take-while #(<= (* % %) n) primes))
                (recur (+ n f) r)
                (lazy-seq (cons n (primes-from (+ n f) r)))))
          wheel (cycle [2 4 2 4 6 2 6 4 2 4 6 6 2 6  4  2
                        6 4 6 8 4 2 4 2 4 8 6 4 6 2  4  6
                        2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10])]
      (primes-from 11 wheel))))
  "Lazy sequence of all the prime numbers.")

次の素数を試すのに、単純に2を足すのではなく、何やら妙な循環数列wheelを使っている。 これはざっくり言えば「ある素数から次の素数までの差は、この数の倍数になる」という列。 詳しくは Recursion on the cycle of gaps を参照。ここで使ってるのはG(7#)なので[2 3 5 7]をあらかじめ計算済みとしている。

さて最後に、「Gaucheらしい」無限素数列を。普通のlazyな言語で無限リストを作るには、余再帰と呼ばれる、(lazy-cons <要素> <再帰呼び出し>) というパターンが使われる。上に示したコードもこのパターンになっている。ただ、このパターンでは <再帰呼び出し> の部分を後で計算するために、その計算環境を保存しておかねばならない。計算環境を保存するのは、Gaucheではクロージャを作るのと同じコストがかかる (もともとlazyな言語なら処理系がもっと効率よく処理するだろうけど)。リストの1要素ごとに新しいクロージャを作るのは、性能上はあんまり嬉しくない。

そこで、Gaucheのlazy sequenceでは、最下層にジェネレータを置いた。ジェネレータは引数無しの手続きで、呼ばれる度に新しい値を返す。値が尽きたら#<eof>を返す。 例えばread-charは引数無しで呼べば次々と#<eof>に達するまで文字を返すので、ジェネレータとして使える。

プリミティブgenerator->lseqはジェネレータを取って、生成される値のlazy sequenceを返す。(generator->lseq read-char) とすれば、文字のlazy sequenceが返ってくるわけだ。ジェネレータの前にあらかじめわかっている要素を付け足すこともできる。(generator->lseq #\a read-char) とすると、最初に#\aが、その後に現在の入力ポートから読んだ文字がくるlazy sequenceとなる。

ジェネレータは作り方によってはnon-consing、つまり全くアロケーションを行わないように書けるので、一要素生成するごとにクロージャを作るのに比べるとかなり速い。

こちらがジェネレータベースの無限素数列。gen-primes-above がジェネレータを作って返す。

(define (gen-primes-above last)
  (^() (set! last (next-prime (+ last 2))) last))

(define *primes* (generator->lseq 2 3 5 (gen-primes-above 5)))

(define (next-prime last)
  (let loop ([ps *primes*])
    (let1 p (car ps)
      (cond [(> (* p p) last) last]
            [(zero? (modulo last p)) (next-prime (+ last 2))]
            [else (loop (cdr ps))]))))

前エントリのコード、lconsを使った余再帰によるコード、ジェネレータベースのコードの ベンチマークを示しておこう。無限素数列から最初の5000個の素数を取ってくるのに要した 時間 (Core2 quad 2.4GHz):

なお、エラトステネスの無限の篩 - 水底で思うことでは「無限リストに無限個の無限リストを適用する」という鮮やかな技を使ってる。時間ができたらこれもGaucheでやってみたい。

(追記2011/11/20 13:07:07 UTC: ジェネレータから遅延シーケンスを作るプリミティブは当初lseqという名前だったが、何を引数とするかわかりにくいのでgenerator->lseqと改名した。本文もそれに合わせて修正してある。lseqという名前は<sequence>から遅延シーケンスを作るコンストラクタにするかもしれん。)

Tags: Gauche, Clojure, Programming