円周率を求めてみる - gifnksmの雑多なメモ

なんとなくやりたくなったからやってみた。
意外と素直にあっさり書けた。やっぱり有理数型があるとこういうの楽だなぁ。
使ったアルゴリズムは、円周率 - Wikipediaで見つけたこの式。
$\pi = 6\sum_{n=0}^{\infty}\frac{(2n)!}{2^{4n+1}(n!)^2(2n+1)}$
arcsin(x)のx=1/2におけるテイラー展開で、収束が速いらしい。へぇー。
で、できたのがこれ。

open Num
let pi length = 
  let min = Int 1 // Int 10 **/ Int length in
  let zero = Int 0 and one = Int 1
  and two = Int 2 and four = Int 4 in
  let prod n m =
    let rec prod' n x =
      if n >/ m then x
      else prod' (n +/ one) (n */ x)
    in prod' n one in
  let factorial n = 
    let rec factorial' n x =
      if n =/ zero then x
      else factorial' (n -/ one) x */ n
    in factorial' n one in
  let a n = 
    (prod (n +/ one) (two */ n)) // 
    ((two **/ (four */ n +/ one)) */ 
     (factorial n) */ (two */ n +/ one)) in
  let rec sum i x =
    let ai = a i in
    if ai </ min then x +/ ai
    else sum (i +/ one) (x +/ ai) in
  Int 6 */ sum zero zero;;

ちょっとばかし解説。
まず、 $a_n=\frac{(2n)!}{2^{4n+1}(n!)^2(2n+1)}$ なので、 $\pi=6\sum_{n=0}^{\infty}a_n$ となります。
有限時間で処理を終わらせる為には、無限個の項を足すわけにはいかないので、 $a_n \lt min$ となるような項まで足したら終了します。それをやってるのがsumですね。
あとは、 $a_n=\frac{{}_{2n}P_{n}}{2^{4n+1}n!(2n+1)}$ として、少し計算量を減らしています。
prod n mはn以上m以下の整数を掛け合わせる関数、factorial nはnの階乗を求める関数です。
実行結果。

# approx_num_fix 500 (pi 500)
- : string =
"+3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194911"

インタプリタ上で計算したら6秒ほどかかりました。小数点以下499桁目まで合っているようです。アルゴリズムはどうやら正しいようですね。

今思ったんですが、これ、 $a_n$ を漸化式で表現したらもっと速くなりそうですね。
後でやってみよう。