「ウォリス積」、「マダヴァの公式」による円周率

「やさしくわかる数学のはなし７７」という本に「ウォリス積」、「マダヴァの公式」による円周率の計算式が載っていましたのでやってみました。

• ウォリス積

import Ratio

-- ウォリス積
--         2・2・4・4・6・6・8・8...
-- 2 pi = -----------------------
--         1・3・3・5・5・7・7・9...

-- 無限数列を作る
-- x : 初期値
-- c : True と False の二進カウンタ
numSeque :: Integer -> Bool -> [Integer]
numSeque x c = x: case c of
                    False -> numSeque x     True
                    True  -> numSeque (x+2) False

-- [2,2,4,4,6,6,8,8,10,10,12,12,14,14...]
nume :: [Integer]
nume   = numSeque 2 False

-- [1,3,3,5,5,7,7,9,9,11,11,13,13,15...]
denomi  :: [Integer]
denomi = tail $ numSeque 1 False

-- n は数列の長さを指定します。 
p :: Int -> Double
p n =  fromRational $ ((product (take n nume)) % (product (take n denomi))) * 2

{-
*Main> mapM_ putStrLn $ map (show.p) $ take 10 [100,200..]
-- =>
3.1260789002154112
3.1337874906281624
3.1363784027059802
3.137677900950936
3.1384588976716157
3.138980103882128
3.139352659682266
3.1396322219293964
3.139849745446556
3.1400238186005973
-}

• マダヴァの公式

import Ratio

-- マダヴァの公式
-- pi / 4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 ...

p2 :: Integer -> Int -> Rational -> Double
p2 n count acc | count == 0  = fromRational $ acc * 4
               | countIsEven = p2 (n + 2) (count-1) $ acc + 1 % n
               | otherwise   = p2 (n + 2) (count-1) $ acc - 1 % n
     where countIsEven = count `mod` 2 == 0

pai :: Int -> Double
pai n = p2 1 n 0

{-
*Main> mapM_ putStrLn $ map (show.pai) $ take 10 [100,200..]
3.131592903558553
3.1365926848388166
3.1382593295155905
3.1390926574960125
3.1395926555897833
3.13992598808053
3.1401640828900828
3.1403426540780734
3.1404815428216173
3.140592653839793
-}

• ゼータ関数

import Ratio

-- ゼータ関数
-- pi**2 / 6 = 1/1**2 + 1/2**2 + 1/3**2 + 1/4**2 + 1/5**2 ....

p2 :: Integer -> Int -> Rational -> Double
p2 n 0     acc = sqrt $ fromRational $ acc * 6
p2 n count acc =  p2 (n + 1) (count-1) $ acc + 1 % (n * n)

pai :: Int -> Double
pai n = p2 1 n 0

{-
ghci>  mapM_ putStrLn $ map (show.pai) $ take 10 [100,200..]
3.132076531809106
3.136826306330968
3.1384132451584095
3.1392074056144135
3.139684123138722
3.140002027036458
3.1402291464241046
3.140399510673162
3.1405320308436044
3.1406380562059932
-}

import Ratio

-- ゼータ関数
-- pi**4 / 90 = 1/1**4 + 1/2**4 + 1/3**4 + 1/4**4 + 1/5**4

p2 :: Integer -> Int -> Rational -> Double
p2 n 0     acc = sqrt $ sqrt $ fromRational $ acc * 90
p2 n count acc =  p2 (n + 1) (count-1) $ acc + 1 % (n * n * n *n)

pai :: Int -> Double
pai n = p2 1 n 0

{-
ghci>  mapM_ putStrLn $ map (show.pai) $ take 10 [100,200..]
3.1415924153073678
3.141592623579992
3.141592644675728
3.141592649824466
3.1415926516604986
3.141592652472745
3.1415926528860942
3.1415926531182436
3.14159265325854
3.1415926533482694
-}