ガウス・ルジャンドルの方法での円周率

java 数学

もう一個だけ、円周率を求めておく
算術幾何平均というのを使う。
これは、a, bについて
a_1=\frac{a+b}{2}, b_1=\sqrt{ab}
とすると
a_{n+1}=\frac{a_n+b_n}{2}, b_{n+1}=\sqrt{a_nb_n}
がnが無限のときanとbnが等しくなるというもの。


で、なんかいろいろ調べてたらどこからひっぱってきたか忘れたのだけど、次のようなプログラムで円周率が出せる。

public class PaiGauss {
    public static void main(String[] args){
        double a = 1;
        double b = Math.sqrt(2) / 2;
        double x = 1;
        double t = 1./4;

        for(int i = 0; i < 3; ++i){
            double y = a;
            a = (a + b) / 2;
            b = Math.sqrt(b * y);
            t = t - x * (y - a) * (y - a);
            x = x * 2;
        }
        System.out.println((a + b) * (a + b) / 4 / t);
        System.out.println(Math.PI);


    }
}


結果はこう。

3.141592653589794
3.141592653589793


なんと、3回の繰り返しでdoubleの限界まで求めれた。ただ、これにはちょっと落とし穴があって、sqrtを使っているので、この中にも繰り返しがある。


まあ、とりあえず、BigDecimalを使って100桁まで求めてみる。
これで、次のように100桁もとめれた。ちょっと改行した。

3.1415926535 8979323846 2643383279 5028841971 939937510
5820974944 5923078164 0628620899 8628034825 342117068

ソースはこれ

続きを読む