円周率でも求めてみるか

エルディンガー

なんだか、java-jaチャットで円周率の話になってた。
で、id:cactusman
4\int_0^1 {(1-x^2)^{\frac{1}{2}}}
でいいんじゃねぇの?って言ってたので、試してみた。
で、こんなコードを書いてみる

public class Pai {
    public static void main(String[] args) {
        double sum = 0;
        double d = 0.01;
        for(double x = 0; x <= 1; x += d){
            sum += Math.sqrt(1-x * x) * d;
        }
        System.out.println(sum * 4);
    }
}

3.1604170317790454になった。精度悪すぎ・・・
案の定、cactusmanに刻みが粗いとdisられた。
dを0.0001にしてみたら、3.1417914777848557になった。改善したとはいえ、せめて自分が覚えてる3.141592くらいは出したいなぁ。


ということで、シンプソン法で二次補正してみる。

public class Pai {
    public static void main(String[] args) {
        double sum = 0;
        int n = 10000;
        double d = 1. / n;
        for(int i = 0; i <= n; ++i){
            double x = i * d;
            double f = Math.sqrt(1-x * x);
            if(i == 0 || i == n){
                sum += f;
            }else if(i % 2 == 0){
                sum += f * 2;
            }else{
                sum += f * 4;
            }
        }
        System.out.println(sum*d/3*4);
    }
}


3.1415921943382403になった!
ソースを見てのとおり、シンプソン法は、端っこはそのまま、偶数番目は2倍、奇数番目は4倍にして足したあと、3で割ります。