無限ストリームと交代級数によるπの近似

以下はSICP 3.5.3を元にしている。

$\pi / 4 = 1 - (1/3) + (1/5) - (1/7) + (1/9) - \cdots$ を使って円周率を近似する。 その際に、Eulerによる「並び加速」を使って収束の速度を上げる。 最後にストリームのストリーム(タブロー)を使う。

\begin{code}
module Main where

pi_sums :: Double -> Double -> [Double]
pi_sums sgn n = (sgn / n) : (pi_sums (-sgn) (n + 2))

partial_sums :: [Double] -> [Double]
partial_sums xss = psum 0 xss
    where
        psum s [] = []
        psum s (x:xs) = (s + x) : psum (s + x) xs

pi_stream :: [Double]
pi_stream = map (*4.0) (partial_sums (pi_sums 1.0 1.0))

display_stream [] = putStr ""
display_stream (x:xs) = do
    putStrLn $ show x
    display_stream xs

euler_transform1 :: [Double] -> [Double]
euler_transform1 (s0:s1:s2:xs) = (s2 - ((s2-s1)^2 / (s0-2*s1+s2))) : (euler_transform1 xs)
euler_transform2 = euler_transform1 . euler_transform1
euler_transform3 = euler_transform2 . euler_transform1

make_tablau :: ([Double] -> [Double]) -> [Double] -> [[Double]]
make_tablau trans [] = []
make_tablau trans s@(x:xs) = (trans s) : (make_tablau trans (trans s))

accel_seq :: ([Double] -> [Double]) -> [Double] -> [Double]
accel_seq trans [] = []
accel_seq trans s = map head (make_tablau trans s)

euler_transformX = accel_seq euler_transform1

t0 = display_stream $ take 5 $ pi_stream
t1 = display_stream $ take 5 $ euler_transform1 pi_stream
t2 = display_stream $ take 5 $ euler_transform2 pi_stream
t3 = display_stream $ take 5 $ euler_transform3 pi_stream
t4 = display_stream $ take 5 $ euler_transformX pi_stream

\end{code}

Main> t0
4.0
2.66666666666667
3.46666666666667
2.8952380952381
3.33968253968254

Main> t1
3.16666666666667
3.13968253968254
3.14207181707182
3.1414067184965
3.14168318920776

Main> t2
3.14187746961483
3.14158972263034
3.14159290941745
3.14159260392652
3.14159266804457

Main> t3
3.14159287451049
3.1415926529004
3.14159265361311
3.14159265358745
3.14159265359021

Main> t4
3.16666666666667
3.14187746961483
3.14159287451049
3.14159265361083
3.14159265358979

[sicp] 2004年12月25日 00:33

■記事一覧

2006 [ 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 ]
2005 [ 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 ]
2004 [ 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 ]

• このページのトップへ
• 最新の記事へ

Bloglines

book(3) / char(2) / compiler(3) / craft(5) / data(7) / enum(1) / geb(2) / hawiki(1) / hugs(1) / info(1) / io(3) / list(2) / map(3) / monad(16) / nobsun(12) / report(4) / sicp(2) / soe(8) / suchthat(2) / yaht(8)

豊かな人生のための四つの法則