ベルヌーイ数を求める関数を少し高速化する

エイダの誕生日らしいのでベルヌーイ数を求めるプログラムを書いてみた

前回のエントリでベルヌーイ数を馬鹿正直に求めるプログラムを書いた．正しく動くけど，20番目あたりからかなり計算に時間がかかるようになってしまった．B(n+1) を求めるのに B(n) .. B(0) を毎回計算しているので，遅いのは当然である． そこで，とりあえず小手先の小細工で少し高速化を図る． B(1) 以外の奇数番目の数が実行した範囲では全部 0 だったので少し調べてみると，どうやら B(3) 以降の奇数番目のベルヌーイ数はすべて 0 になるらしい．

http://homepage3.nifty.com/y_sugi/gf/gf11.htm

なので，次のようにしてしまって差し支えない．

import Data.Ratio

bernoulliNum :: Integer -> Ratio Integer
bernoulliNum 0 = toRational 1
bernoulliNum 1 = -1 % 2
bernoulliNum n
    | even n = (-1) % (n+1) * (sum $ map (\i -> bernoulliNum i * combi (n+1) i) [0,1..n-1])
    | otherwise = toRational 0
    where
        combi n k = let fact n = product [n, n-1 .. 1] in
                    fact n % (fact (n-k) * fact k)

これで少し早くなったが，計算量はせいぜい半分にしかならないので，やはり 40 番目あたりで計算がきつくなってきてしまう．

これに対しては，やはりメモ化などのもう少し真っ当なやり方が必要だけど，Haskell だと普通は変数を変更できないので古いメモ化用テーブルを取って更新したメモ化用テーブルを返す関数などを書かなければならず，実装が思いつかないので今度調べてみる．