Haskellの「fib = 1:1:zipWith (+) fib (tail fib)」はとても遅い

ラボの昼休みに光成さん、中谷さんとご飯を食べながら話した内容を一応ざっくりとまとめておく。
発端はたしか最近Haskellを勉強の光成さんが、Haskellのかっこいいsieveは実はとても遅い(俺は Haskell の　sieve についてとんでもない思い違いをしていたようだ...)という話を見て、同様にかっこいいけど遅い下記のフィボナッチ数列の定義の速度を調べてみたら2.5乗くらいのオーダーになっていたという話だったかと思う。

fib = 1:1:zipWith (+) fib (tail fib)

僕も確認するために、コマンドライン引数でNを与えられるフィボナッチ数列のN番目を求めるコードを書いた。

import System

fib = 1:1:zipWith (+) fib (tail fib)

main = do
  args <- getArgs
  print $ (0 *) $ (fib !!) $ read $ args !! 0

(0 *)は大量の数字の列をコンソールに表示するところで時間がかかると測りたいものが測れないからつけてある。人間には常に0になるのがわかるけどもHaskell処理系にはわからないので、きちんとfibを計算してから0を掛けて結果の0を表示してくれる。これをghc -O3でコンパイルする。そのまま走らせると「Stack space overflow」というエラーで死んでしまうので+RTS -K200Mをつけて走らせる。

$ time ./a.out 100000 +RTS -K200M
0

real	0m1.080s
user	0m1.045s
sys	0m0.029s

$ time ./a.out 200000 +RTS -K200M
0

real	0m6.420s
user	0m6.313s
sys	0m0.072s

$ time ./a.out 300000 +RTS -K200M
0

real	0m18.477s
user	0m18.264s
sys	0m0.145s

これはおおよそO(N ** 2.6)のオーダーである。ちなみにfib(1000000)は5分経っても終わらなかったので中断した。

>> log(6.42 / 1.08) / log(2)
2.5715419849588343

>> log(18.477 / 6.42) / log(1.5)
2.6071505945442417

この時点で「次の定義は線形時間でフィボナッチ数列のリストを生成する」っていうWikipediaの記述は事実に反するということがわかった。それで終りにしてもよかったのだけどもまだお弁当が残っていたので「いつからこの間違った記述があるんだろう」という話題に。

それは履歴を見ればわかりますよ、と見つけてきたのがこれ。2006年7月5日の更新。そしてこの更新のコメントで英語版からの翻訳であることがわかる。当時の英語版を見ると確かにlinearって書いてある。そして2010年の2月に修正されている。"linear time"から"efficient implementation"へ。efficientって何に比べてだよ…。

ちなみにPythonで素朴にループするコードを書いてみる。

import sys

def fib(n):
    if n == 0 or n == 1: return 1
    x = y = 1
    while n > 1:
        x, y = y, x + y
        n -= 1
    return y
        
        
print 0 * fib(int(sys.argv[1]))

これの実行結果は下のようになる。驚いたことにN = 100000でHaskellよりわずかに速く、N = 300000だと2倍以上速い。むしろコンパイルしてあるくせにインタプリタのPythonに負けるくらいHaskellのコードが遅い。

$ time python fib.py 100000
0

real	0m0.968s
user	0m0.937s
sys	0m0.013s

$ time python fib.py 200000
0

real	0m3.687s
user	0m3.658s
sys	0m0.016s

$ time python fib.py 300000
0

real	0m8.205s
user	0m8.169s
sys	0m0.020s

$ time python fib.py 1000000
0

real	1m31.549s
user	1m29.910s
sys	0m0.148s

これはおよそO(N ** 2)のオーダーだ

>> log(3.687 / 0.968) / log(2)
1.9293684638192232

>> log(8.205 / 3.687) / log(1.5)
1.9728716214877173

>> log(91.549 / 8.205) / log(10 / 3.0)
2.0034760044903011

なぜO(N)ではなく2乗のオーダーになるかと言うと、フィボナッチ数がO(1.618 ** N)のオーダーで増えるからその数の桁数はO(N)のオーダーで増え、1回の足し算にかかる時間がO(N)のオーダーで増えるから。だから順番に足していく方式では2乗のオーダーまではしかたない。だけどHaskellの2.6乗の0.6はどこから来たんだ。

Haskellでefficientなコードを書いてみた。泥臭く自前で末尾再帰している。

import System

fib 0 = 1
fib 1 = 1
fib n = fib' n 1 1

fib' 1 _ y = y
fib' n x y = fib' (n - 1) y (x + y)

main = do
  args <- getArgs
  print $ (0 *) $ fib $ read $ args !! 0

これならPythonの素朴なコードに比べて10倍速いし、ちゃんと2乗のオーダーだ。

$ time ./a.out 100000
0

real	0m0.088s
user	0m0.075s
sys	0m0.006s
etude$ time ./a.out 200000
0

real	0m0.346s
user	0m0.329s
sys	0m0.011s
etude$ time ./a.out 300000
0

real	0m0.814s
user	0m0.789s
sys	0m0.019s
>> log(0.346 / 0.088) / log(2)
1.9751966089994275

>> log(0.814 / 0.346) / log(1.5)
2.1099758618849926

まとめ

  • Haskellの有名なフィボナッチ数列の定義「fib = 1:1:zipWith (+) fib (tail fib)」はWikipediaに「線形時間」と書かれているが事実ではない
  • 英語版には「An efficient implementation」と書かれているが、N = 100000でインタプリタ言語であるPythonより遅く、O(N ** 2.6)である。
  • この定義より10倍以上速く、O(N ** 2)であるような実装も容易

というわけで「fib = 1:1:zipWith (+) fib (tail fib)」は「Haskellらしさをアピールできるトリッキーなコードであるが、速度は遅い」と認識するのが正しい。

まとめ追記