2012/07/13

Unbounded spigot algorithm

ふとしたきっかけでπなどの値を「頭から順番に」求めるSpigotアルゴリズムを 見直していたら、「メモリの許す限り順番に値を計算し続ける」アルゴリズムを示した 論文を見つけた。

• Jeremy Gibbons, Unbounded Spigot Algorithms for the Digits of Pi

よくあるspigotアルゴリズムは計算したい桁数を最初に決めてデータ構造を初期化する必要がある。Gaucheのexamples/spigotに入ってるのもそれ。 一方、こちらのアルゴリズムはあらかじめ精度を決めておく必要がない。

論文のコードはHaskellで書かれている。

piとpiLは使っている級数の違い。後者の方が生成が速い。論文には他にもう一種類出ている。

実行例。結果は無限数列で返ってくるので欲しいところまで取れば良い。

*Main> take 100 Main.piL
[3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3,3,8,3,2,7,9,5,0,2,8,8,4,1,9,7,1,6,9,3,9,9,3,7,5,1,0,5,8,2,0,9,7,4,9,4,4,5,9,2,3,0,7,8,1,6,4,0,6,2,8,6,2,0,8,9,9,8,6,2,8,0,3,4,8,2,5,3,4,2,1,1,7,0,6,7]

遅延評価に大きく依存したコードなので通常のSchemeに直すのは面倒そうだが、 gauche.lazy を使えばわりと素直にGaucheに移植できる。 元コードのパターンマッチ部分を手で展開するのはさすがにしんどかったので ちょっとマクロの手を借りた (define*のところ)。

あと、上の素直なコードだと状態配列がすぐに巨大なbignumになって、多倍長演算の遅い Gaucheには負担なので、可能な限り約分するようにしている (lftのところ)。

pi、piL はジェネレータを返すので、こんな感じで結果を取り出せる。

gosh> (generator->list (piL) 100)
(3 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 8 3 2 7 9 5 0 2 8 8 4 1 9 7 1 6 9 3 9 9 3 7 5 1 0 5 8 2 0 9 7 4 9 4 4 5 9 2 3 0 7 8 1 6 4 0 6 2 8 6 2 0 8 9 9 8 6 2 8 0 3 4 8 2 5 3 4 2 1 1 7 0 6 7)

先頭から次々に、いくらでも桁を取り出せるのは魅力的なんだけれど、 実際に走らせてみると状態配列(lft)の値が(適宜約分しても)どんどん大きくなるので、 実は計算量はO(n^2)とかになってるっぽい。桁数が進むとどんどん遅くなる。 なのでおもしろいけれど実用的ではなさそうだ。

Tags: Gauche, Haskell, Programming