エイダの誕生日らしいのでベルヌーイ数を求めるプログラムを書いてみた
今日はエイダ・ラブレス の誕生日らしい.彼女は世界で初めてコンピュータプログラムを書いた人物で,プログラミング言語名などで有名だけど,その世界初のプログラムがベルヌーイ数を求めるプログラムだと初めて知ったのでとりあえず書いてみることにした.ベルヌーイ数を使うとマクローリン展開の難しい一部の関数が簡単に表せるらしい.
ベルヌーイ数は有理数であるので,一定の精度までしか計算出来ない Double では駄目だから,Haskell 標準の有理数ライブラリ Data.Ratio
を使うことにした.
import Data.Ratio bernoulliNum :: Integer -> Ratio Integer bernoulliNum 0 = toRational 1 bernoulliNum n = (-1) % (n+1) * (sum $ map (\i -> bernoulliNum i * combi (n+1) i) [0,1..n-1]) where combi n k = let fact n = product [n, n-1 .. 1] in fact n % (fact (n-k) * fact k)
試しに 20 項まで求めてみる.
import Control.Monad (forM_) main :: IO () main = forM_ [0..20] $ \i -> do print $ bernoulliNum i return i
計算結果:
1 % 1 (-1) % 2 1 % 6 0 % 1 (-1) % 30 0 % 1 1 % 42 0 % 1 (-1) % 30 0 % 1 5 % 66 0 % 1 (-691) % 2730 0 % 1 7 % 6 0 % 1 (-3617) % 510 0 % 1 43867 % 798 0 % 1 (-174611) % 330
どうやら正しく求められているらしい. ベルヌーイ数の定義から
らしいので,気が向いたら検算してみる.
追記
現実逃避検算する気になったのでやってみた.
main :: IO () main = forM_ [1..6] $ \i -> do putStrLn $ "左辺:" ++ (show $ lhs i) putStrLn $ "右辺:" ++ (show $ rhs i) putStrLn "" return i where e = 2.718281828459045235360287471352 fact n = product [n, n-1 .. 1] lhs :: Integer -> Double lhs x = realToFrac x / (e^x - 1) rhs x = fromRational $ sum $ map (\n -> bernoulliNum n * ((x^n) % fact n)) [0..18]
lhs
が上の式の真ん中の項,rhs
が上の式の右の項を表している.ただ,ネイピア数を使うので左辺は Double で計算し,右辺はベルヌーイ数の第 18 項までで近似する.
計算結果:
左辺:0.5819767068693265 右辺:0.5819767068693267 左辺:0.31303528549933135 右辺:0.31303528570640216 左辺:0.1571870894737679 右辺:0.15718770702004148 左辺:7.46294414550962e-2 右辺:7.479960583486796e-2 左辺:3.3918274531521166e-2 右辺:4.661754185811595e-2 左辺:1.4909469941067517e-2 右辺:0.43081551362013865
早くも 5 項目で大きく食い違ってしまった.やっぱりたった18項目までの計算 & Double では駄目だなぁと思いつつ,数値計算は大変だと思った.