Gauss-Legendreの2次の収束公式でπを計算してみました

算術幾何平均による円周率πの計算公式を参考にπを計算してみました。多倍長演算をごにょごにょするのが大変そうだったので、とりあえずDoubleの範囲で。

module Main (main) where

import Text.Printf (printf)

calcPi :: Int -> Double
calcPi n = calc n 1.0 (1.0 / (sqrt 2.0)) 1.0 4.0
  where
    calc 0 a b t _ = (a + b) ^ 2 / t
    calc n a b t x =
      let a' = (a + b) / 2.0
          b' = sqrt (a * b)
          t' = t - x * ((a - a') ^ 2)
          x' = x * 2
      in calc (n - 1) a' b' t' x'

main :: IO ()
main = mapM_ (starling (printf "n = %2d ... %f\n") calcPi) [1 .. 10]
  where
    starling f g x = f x $ g x

-- => n =  1 ... 3.1405792505221686
--    n =  2 ... 3.141592646213543
--    n =  3 ... 3.141592653589794
--    n =  4 ... 3.141592653589794
--    n =  5 ... 3.141592653589794
--    n =  6 ... 3.141592653589794
--    n =  7 ... 3.141592653589794
--    n =  8 ... 3.141592653589794
--    n =  9 ... 3.141592653589794
--    n = 10 ... 3.141592653589794

うーん、3回も繰り返すとDoubleの範囲では収束するのか…
参照: ガウス＝ルジャンドルのアルゴリズム - Wikipedia