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

もう一個だけ、円周率を求めておく
算術幾何平均というのを使う。
これは、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

ソースはこれ

import java.math.BigDecimal;
import java.math.MathContext;

public class PaiGaussBig {
    public static void main(String[] args){
        MathContext mc = new MathContext(100);
        BigDecimal b2 = new BigDecimal(2);
        BigDecimal b4 = new BigDecimal(4);
        BigDecimal a = BigDecimal.ONE;
        BigDecimal b = sqrt(b2, mc.getPrecision()).divide(b2, mc); //Math.sqrt(2) / 2;
        BigDecimal x = BigDecimal.ONE;
        BigDecimal t = BigDecimal.ONE.divide(b4, mc);
        for (int i = 0; i < 6; ++i) {
            BigDecimal y = a;
            a = (a.add(b)).divide(b2, mc);
            b = sqrt(b.multiply(y), mc.getPrecision());
            t = t.subtract(x.multiply((y.subtract(a)).pow(2)));
            x = x.multiply(b2);
        }
        BigDecimal pi = a.add(b).pow(2).divide(b4.multiply(t), mc);
        System.out.println(pi);
    }

    private static BigDecimal sqrt(BigDecimal a, int scale){
        //とりあえずdoubleのsqrtを求める
        BigDecimal x = new BigDecimal(
            Math.sqrt(a.doubleValue()), MathContext.DECIMAL64);
        if(scale < 17) return x;

        BigDecimal b2 = new BigDecimal(2);
        for(int tempScale = 16; tempScale < scale; tempScale *= 2){
            //x = x - (x * x - a) / (2 * x);
            x = x.subtract(
                    x.multiply(x).subtract(a).divide(
                    x.multiply(b2), Math.min(scale, tempScale * 2),
                    BigDecimal.ROUND_HALF_EVEN));
        }
        return x;
    }
}