数値計算をやってみた

はじめに

この記事は NSEG Advent Calendar 2015 の21日目の記事です。

なんとなくプログラムを書きたくなりました。Advent Calendarの時期でもあるし、何かネタがないかなあと考えてみたら、昔よくやっていた数値計算を久々にやりたくなりました。
というわけで、円周率の計算をさらっと書いてみました。昔はC言語でしたが、今回はHaskellで。

ガウスルジャンドルの公式

円周率の計算といえばこれでしょう。おなじみの公式ですね。

-- Gauss-Legendre algorithm

a :: Double -> Double
a 0 = 1
a n = ((a(n-1))+(b(n-1))) / 2

b :: Double -> Double
b 0 = 1 / sqrt 2
b n = sqrt((a(n-1))*(b(n-1)))

t :: Double -> Double
t 0 = 1 / 4
t n = (t(n-1)) - (p(n-1)) * (a(n-1) - a n) ** 2

p :: Double -> Double
p 0 = 1
p n = 2 * p(n-1)

pi' :: Double -> Double
pi' n = (a n + b n)**2 / (4*t n)

main :: IO ()
main = print $ pi' 3

実行結果

3.141592653589794

数式がそのまま書ける感じだしそのままでも読みやすいし、とても素敵な公式。

ボールウェインの4次の収束公式

収束効率が高くて良い感じです。

-- Borwein's quartically convergent algorithm

a :: Double -> Double
a 0 = 6 - 4 * sqrt 2
a n = a(n-1) * (1+y n)**4
      - 2**(2*(n-1) + 3) * y n * (1+y n+(y n)**2)

y :: Double -> Double
y 0 = (sqrt 2) - 1
y n = (1 - (1 - (y(n-1)) ** 4) ** (1 / 4))
      / (1 + (1 - (y(n-1)) ** 4) ** (1 / 4))

pi' :: Double -> Double
pi' n = 1 / a n

main :: IO ()
main = print $ pi' 2

実行結果

3.1415926535897927

1つあたりの計算はごちゃごちゃやってるけど、それでもだいぶシンプル。

あとがき

プログラムについてはツッコミどころは大量にあるとは思いますが、そこはまあ即席なんで勘弁してください。