PI (円周率) を求める (2)

PI (円周率) を求めるで既にガウス・ルジャンドルの方法とモンテカルロ法で円周率を求めていますが、円周率の計算 - みずぴー日記に記載されている方法で円周率を求めてみます。 これは SICP と呼ばれている計算機プログラムの構造と解釈に掲載されているそうですが、以下のようにすると π/8 が求まるそうです。

 1     1     1
--- + --- + ---- + ....
1*3   5*7   9*11

この級数が (i - 1) * 4 + 1 で表されることが分かれば簡単に awk で計算できるでしょう。

#! /usr/bin/gawk -f
# pi_3.awk

BEGIN {
    num = ARGV[1] ? ARGV[1] : 10;

    for (i = 1; i <= num; i++) {
        num1 = 4 * (i - 1) + 1;
        num2 = num1 + 2;
        pi = pi + 1 / (num1 * num2);
        printf("%f\n", pi * 8);
    }
}

実行してみます。

$ nawk -f pi_3.awk
2.666667
2.895238
2.976046
3.017072
3.041840
3.058403
3.070255
3.079153
3.086080
3.091624

実際に円周率は 3.141592 ということが知られていますが、何回くらいで収束するかを調べてみます。

$ nawk -f pi_3.awk 10 | tail -n 1
3.091624

$ nawk -f pi_3.awk 100 | tail -n 1
3.136593

$ nawk -f pi_3.awk 1000 | tail -n 1
3.141093

$ nawk -f pi_3.awk 10000 | tail -n 1
3.141543

$ nawk -f pi_3.awk 100000 | tail -n 1
3.141588

$ nawk -f pi_3.awk 1000000 | tail -n 1
3.141592

何と 100 万回もかかってしまいました。