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

1843 年のコンピューター先見者エイダ ラブレスを讃えて

今日はエイダ・ラブレス の誕生日らしい.彼女は世界で初めてコンピュータプログラムを書いた人物で,プログラミング言語名などで有名だけど,その世界初のプログラムがベルヌーイ数を求めるプログラムだと初めて知ったのでとりあえず書いてみることにした.ベルヌーイ数を使うとマクローリン展開の難しい一部の関数が簡単に表せるらしい.

ベルヌーイ数は有理数であるので,一定の精度までしか計算出来ない 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 では駄目だなぁと思いつつ,数値計算は大変だと思った.