読者です 読者をやめる 読者になる 読者になる

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

Haskell

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

前回のエントリでベルヌーイ数を馬鹿正直に求めるプログラムを書いた.正しく動くけど,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 だと普通は変数を変更できないので古いメモ化用テーブルを取って更新したメモ化用テーブルを返す関数などを書かなければならず,実装が思いつかないので今度調べてみる.