円周率を求めてみる

SML 30分プログラム

30分プログラム、その743。今日3/14は円周率の日らしいです。というわけで、円周率を求めてみました。
計算式は円周率の計算 - みずぴー日記で使ったやつと同じです。
直接計算するのは芸がないので、無限ストリームを定義して、その上でだんだんと精度が上っていくπのストリームを作りました。

使い方

- Stream.nth pi 100;
val it = 3.09616152646 : real
- Stream.nth pi 1000;
val it = 3.13659268484 : real
- Stream.nth pi 1000;
val it = 3.14109265362 : real

ソースコード

structure Stream =
struct
  type 'a thunk = unit -> 'a;
  datatype 'a stream = Cons of 'a * ('a stream) thunk;

  fun force f = f ();
  fun cons x y = Cons (x,y);
  fun hd (Cons(x,_)) = x;
  fun tl (Cons(_,xs)) = force xs;

  fun nth xs 0 = hd xs
    | nth xs n = nth (tl xs) (n-1);

  fun scanl f init (Cons (x,xs)) =
      cons init (fn () => scanl f (f (x,init)) (force xs));

  fun map f (Cons (x,xs)) =
      cons (f x) (fn () => (map f (force xs)));
end;

fun term i =
    let
	val a = (Real.fromInt i) * 4.0 + 1.0
	val b = a + 2.0
    in
	1.0 / (a * b)
    end;

fun pi_8 n =
    Stream.cons (term n) (fn () => pi_8 (n+1));

val pi =
    Stream.map (fn x => 8.0 * x) (Stream.scanl (fn (x,y) => x + y) 0.0 (pi_8 0));

参考