さらに円周率を求めてみる

こないだ、積分の練習に円周率を求めてみたんだけど、もうちょっとちゃんと円周率を求めてみる
とりあえず、モンテカルロでってことなので、やってみる。
円周率をモンテカルロで求めるときには、ランダムに選んだ点が円に入るかどうかで判定します。
つまりランダムにx、yを選んで

を満たす割合を求めます。

public class PaiMontecarlo {
    public static void main(String[] args){
        Random r = new MTRandom();
        int in = 0;
        int total = 100000;
        for(int i = 0; i < total; ++i){
            double x = r.nextDouble();
            double y = r.nextDouble();
            if(x * x + y * y < 1) ++in;
        }
        System.out.println((double)in / total * 4);
    }
}

そしたらこう。

3.14888

3桁精度出すのがやっと。計算結果として桁数を増やすにも実行回数の桁数を増やす必要がある。あと、乱数がだいたいちゃんとばらついていると円周率出るので、あんまり乱数の性能は関係ない気がする。

ということで、まじめに級数展開したもので計算。

(はてなのtexモジュールの動きが変なので、0足してます。数学的な意味はないです。)

public class PaiNormal {
    public static void main(String[] args){
        double d = 0;
        double sig = 1;
        for(int i = 1; i < 1000000; i += 2){
            d += sig / i;
            sig *= -1;
        }
        System.out.println(d * 4);
        System.out.println(Math.PI);
    }
}

そしたらこう。

3.141590653589692
3.141592653589793

50万回計算してやっと小数点以下5桁がもとまった。

これはこれでちょっと収束遅すぎなので、オイラーの加速を使って収束を速くしてみる。
π/2=1+2/3×1/2+2・4/3・5×1/2^2+2・4・6/3・5・7×1/2^3+・・・
(もうどうやってもtexモジュールが表示されない)

public class PaiEularAcc {
    public static void main(String[] args){
        double pi = 1;
        double a = 1;
        double b = 1. / 2;
        for(int i = 2; i < 99; i += 2){
            a *= i / (double)(i + 1);
            pi += a * b;
            b /= 2;
        }

        System.out.println(pi * 2);
        System.out.println(Math.PI);
    }
}

そしたらこう。

3.1415926535897922
3.141592653589793

50回くらいでdoubleの限界まで来た。速い。

まあ、でもここまでならSICPにも書いてあるので、あまり面白くない。
ということでマチンの公式を使う。

public class PaiMachin {
    public static void main(String[] args){
        double sig = 1;
        double a = 0;
        double b = 0;
        double s1 = 1. / 5;
        double s2 = 1. / 239;
        for(int i = 1; i < 22; i += 2){
            a += sig * s1 / i;
            b += sig * s2 / i;
            s1 /= 5 * 5;
            s2 /= 239 * 239;
            sig *= -1;
        }
        System.out.println(a * 16 - b * 4);
        System.out.println(Math.PI);
    }
}

そしたらこう。

3.141592653589794
3.141592653589793

10回程度の繰り返しでdoubleの限界まできた。

じゃあ、ということで、もちょっと桁数を増やす。doulbeの限界を超えるので、BigDecimalを使うことにする。

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

public class PaiMachinBig {
    public static void main(String[] args){
        MathContext mc = new MathContext(100);
        BigDecimal sig = BigDecimal.ONE; // 1
        BigDecimal a = BigDecimal.ZERO;  // 0
        BigDecimal b = BigDecimal.ZERO;  // 0
        BigDecimal d5 = new BigDecimal(5, mc); // 5
        BigDecimal d239 = new BigDecimal(239, mc); // 239
        BigDecimal s1 = BigDecimal.ONE.divide(d5, mc); // 1 / 5
        BigDecimal s2 = BigDecimal.ONE.divide(d239, mc); // 1 / 239
        BigDecimal d5_2 = d5.multiply(d5); // 5*5
        BigDecimal d239_2 = d239.multiply(d239); // 239 * 239
        System.out.println("開始");
        long start = System.nanoTime();
        for(int i = 1; i < 140; i += 2){
            BigDecimal bi = new BigDecimal(i, mc);
            a = a.add(s1.divide(bi, mc).multiply(sig)); // a + s1 / i * sig
            b = b.add(s2.divide(bi, mc).multiply(sig)); // b + s2 / i * sig
            s1 = s1.divide(d5_2, mc);   // s1 / 5^2
            s2 = s2.divide(d239_2, mc); // s2 / 239^2
            sig = sig.negate(); // -sig
        }
        System.out.println((System.nanoTime() - start) / 1000 / 1000 / 1000.);
        a = a.multiply(new BigDecimal(16, mc));
        b = b.multiply(new BigDecimal(4, mc));

        System.out.println(a.subtract(b, mc));
    }
}

こんな感じで出た(10桁で区切ってます)

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 342117068

75回程度でここまで来た。マチンすごい！
kisaragiweb.jpに載ってるのと同じ数字が出た。
ループの条件をi < 300くらいにして150回繰り返すようにすると、200桁までちゃんと求まる。
けど、BigDecimalのプログラムめんどすぎ泣ける。Scala使ったほうがいいね。
ところで、円周率100万桁のページを表示するとFireFox3のCPU使用率が100%になる。どうにかならんもんか

2021/9/17 追記 リンク切れをそのままにしてますが、円周率の確認はこちらあたりを。
https://www.tstcl.jp/ja/randd/pi.php