PI (円周率) を求める

ガウス・ルジャンドルの方法での円周率というページからですが、高速な円周率の計算をやってみます。 円周率の計算には様々な方法がありますが、ここではガウス＝ルジャンドルのアルゴリズムを使ってみます。

このアルゴリズムは非常に収束性が高いことが特徴ですが、awk で書くと以下のような感じになります。

#! /bin/awk -f
# gauss_legendre.awk
#   ガウス＝ルジャンドルのアルゴリズムで PI を求める
#   usage: awk -f gauss_legendre.awk [num]
#   Ref: http://ja.wikipedia.org/wiki/ガウス＝ルジャンドルのアルゴリズム

BEGIN {
    # 引数から繰り返し回数の設定
    if (ARGV[1] == "") {
        num = 3;
    } else {
        num = ARGV[1];
    }

    printf("%1.16f\n", gauss_legendre_pi(num));
}

# gauss_legendre_pi():  ガウス＝ルジャンドルのアルゴリズムで num 回の収
#                       束で得られる PI を返す
#   in:     整数 num
#   out:    ガウス＝ルジャンドルのアルゴリズムで num 回の収束で得られる PI
function gauss_legendre_pi(num,     a, b, x, t, y, i) {
    a = 1;
    b = sqrt(2) / 2;
    x = 1;
    t = 1 / 4;

    for (i = 0; i < num; i++) {
        y = a;
        a = (a + b) / 2;
        b = sqrt(b * y);
        t = t - x * (y - a) * (y - a);
        x = x * 2;
    }

    return (a + b) * (a + b) / 4 / t;
}

式そのものも非常に簡単であることも大きな特徴と言えるでしょう。

$ nawk -f gauss_legendre.awk 1
3.1405792505221686

$ nawk -f gauss_legendre.awk 2
3.1415926462135428

$ nawk -f gauss_legendre.awk 3
3.1415926535897940

$ nawk -f gauss_legendre.awk 4
3.1415926535897940

計算結果から分かるように 3 回のループで awk の倍精度の限界まで計算し切れていることが分かります。

一方、良く子供の教科書に載っているモンテカルロ法での計算を合わせて載せておきます。

#! /bin/awk -f
# pi.awk
#   モンテカルロ法で PI を求める
#   usage: awk -f pi.awk [num]

BEGIN {
    # 引数から繰り返し回数の設定
    if (ARGV[1] == "") {
        num = 3;
    } else {
        num = ARGV[1];
    }

    printf("%1.16f\n", monte_carlo_pi(num));

}

# monte_carlo_pi():     モンテカルロ法で PI を求めて返す
#   in:     整数 num
#   out:    モンテカルロ法で num 回の収束で得られる PI
function monte_carlo_pi(num,    n, i, x, y) {
    n = 0;
    srand();
    for (i = 0; i < num; i++) {
        x = rand();
        y = rand();
        if (x ^ 2 + y ^ 2 < 1) {
            n++;
        }
    }

    return  n / num * 4;
}

モンテカルロ法には収束性を向上させるため、様々な工夫や条件を課すことがありますが、ここでは最も収束が遅そうなものを挙げています。

$ nawk -f pi.awk 1
4.0000000000000000

$ nawk -f pi.awk 10
3.6000000000000001

$ nawk -f pi.awk 100
3.3599999999999999

$ nawk -f pi.awk 1000
3.1320000000000001

$ nawk -f pi.awk 10000
3.1156000000000001

$ nawk -f pi.awk 100000
3.1351200000000001

$ nawk -f pi.awk 1000000
3.1411799999999999

$ nawk -f pi.awk 10000000
3.1412635999999998

まだまだ収束していません。 その上、ループ回数が増えるに従い時間も多く必要としています。

コンピューターは 0, 1 の世界と言われることもありますが、PI のように割り切れない計算 (終わりのない計算) を探求してみることも面白いと思います。