4回目のHaskell
ちょっとずつ慣れてきたがしかしまだ4回しか書いていないのでまだまだ初心者の領域。まああせらずに一歩ずつ進んでいきます。
今日のお題は"pi2.c"円周率の計算ですね。元のアルゴリズムを調べてみたところ、どうやらAGMというアルゴリズムみたいですね。
http://www.pluto.ai.kyutech.ac.jp/plt/matumoto/pi_small/node11.html
この式どっかで見たことあるので多分以前同じプログラムを書いたことがあると思われる。
とりあえず、この漸化式をそのまんま書いてみた。
main = print( p 1 ) >> print( p 2 ) >> print( p 3 ) a 0 = 1.0 a n = ( a (n-1) + b (n-1) ) / 2 b 0 = 1.0 / sqrt 2.0 b n = sqrt ( a (n-1) * b (n-1) ) t 0 = 1.0 / 4.0 t n = t (n-1) - 2**(n-1) * ( a (n-1) - a n ) ** 2 p n = ( a n + b n ) ** 2 / ( 4 * t n )
本当にそのまんまですね(^^; C言語で書かれたプログラムよりずっと可読性がよいです。数学者に好まれるという理由がよくわかる気がします。Haskellの威力が発揮できた例でしょうか?
さて、発展問題(?)として数列を無限リストとして作成する方法とかも考えたんだけど・・・挫折しましたorz
ab = (1.0, 1/sqrt 2.0) : map (\(a,b)->((a+b)/2,sqrt(a*b))) ab
aとbのリストを作ることはできたんですがtは単純な漸化式でなくnが入っているのでどうしても作れない。きっとやり方があるんだと思うんですが私の実力では限界でした。そのうち書ける様になるかなぁ〜(遠