算術幾何平均(Arithmetic Geometric Mean)による円周率πの計算公式


πの計算には π/4 = arctan(1) に基づく各種の方法が多く使われていたが、現在ではガウス(Johann Gauss 1777-1855)の算術幾何平均による完全楕円積分のアルゴリズムと、楕円積分に関するルジャンドル(Adrien-Marie Legendre 1752-1833)の関係式に基づく算術幾何平均アルゴリズムが主流*Note1になっている。

【ガウス・ルジャンドル(Gauss-Legendre)の2次の収束公式によるプログラム例】
very_long_float
pi2(iter)
    int             iter;
{
    int             i;
    very_long_float a, b, t, x, y;

    a = 1;
    b = 1 / sqrt(2);
    t = 1;
    x = 4;
    for (i = 0; i < iter; i++) {
        y = a;
        a = (a + b) / 2;
        b = sqrt(y * b);
        t -= x * sqr(y - a);
        x *= 2;
    }
    return (sqr(a + b) / t);
}
(但し、sqr(x) は x2, sqrt(x) は x1/2)
この公式は2次の収束を示し、n桁の答えが log2(n) 回程度の反復で求まる。
繰り返しごとの誤差は次のようになる。
err( 1) = -3.2257622e-4
err( 2) = -2.3479336e-9
err( 3) = -5.8292283e-20
err( 4) = -1.7418264e-41
err( 5) = -7.6589246e-85
err( 6) = -7.3484406e-172
err( 7) = -3.3697129e-346
err( 8) = -3.5362794e-695
err( 9) = -1.9454428e-1393
err(10) = -2.9425892e-2790

【ボールウェイン(P.Borwein)の4次の収束公式によるプログラム例】
very_long_float
pi(iter)
    int             iter;
{
    int             i;
    very_long_float a, b, y, t;

    y = sqrt(2) - 1;
    a = 6 - sqrt(32);
    b = 8;
    for (i = 0; i < iter; i++) {
        t = qroot(1 - pow4(y));
        y = (1 - t) / (1 + t);
        t = sqr(1 + y);
        a = sqr(t) * a - y * (t - y) * b;
        b *= 4;
    }
    return (1 / a);
}
(但し、sqr(x) は x2, pow4(x) は x4, sqrt(x) は x1/2, qroot(x) は x1/4)
この公式は4次の収束を示し、繰り返しごとの誤差は次のようになる。
err(1) = -2.3479336e-9
err(2) = -1.7418264e-41
err(3) = -7.3484406e-172
err(4) = -3.5362794e-695
err(5) = -2.9425892e-2790

おまけ。ニュートン法的な Beeler (1972) の公式
very_long_float h;
very_long_float a0 = 355.0/113.0;

/* ai+1 = ai + sin(ai) は π に収束する。
 * ai+1 = ai - tan(ai) でも同じ。
 */ 
while(abs(h = sin(a)) > epsilon) {
    a += h;
}
この公式は3次の収束を示し、繰り返しごとの誤差は次のようになる。
err(1) = 1.007e-21
err(2) = 1.680e-63
err(3) = 7.804e-189

Note1
1999年9月時点において。2002年12月6日 東京大学の金田らによるπの世界新記録の更新(1兆2400億桁)には分割有理数化法(Divide and Rationalize methode)が使用された。 2009年8月の筑波大学の高橋大介らによるπの世界記録の更新にもガウス・ルジャンドルのアルゴリズムとボールウェインのアルゴリズムが使用された。

しかし 2009年12月のたった1台のデスクトップPCを使ったFabrice Bellard によるπの世界新記録の更新(2兆6999億桁)には Ramanujan のπの無限級数の変種 Chudnovsky 兄弟の式が使用された。
Chudnovsky's formula
この公式の繰り返し毎の誤差は次のようになる。
err(0)= 1.8790e-14
err(1)= -9.7991e-29
err(2)= 5.4766e-43
err(3)= -3.1840e-57
err(4)= 1.8969e-71
err(5)= -1.1488e-85
err(6)= 7.0403e-100

SEE ALSO:
逆数と平方根を求める高次収束アルゴリズム
参考文献
  1. J.Borwein and P.Borwein, Ramanujan and Pi, IEEE, Science and Applications, Supercomputing 88, Volume II, 117-128, 1988
  2. D.Takahashi, Y.Kanada, Improvement of the Algorithms for π Calculation:The Gauss-Legendre Algorithm and the Borwein's Quartically Convergent Algorithm, IPSJ JOURNAL Abstract Vol.38 No.11 - 033
  3. 後 保範, 金田康正, 高橋大介, 級数に基づく多数計算の演算量削減を実現する分割有理数化法, 情報処理学会論文誌 第41巻 第6号(2002年6月) P.1,811〜1,819
  4. David H. Bailey, Peter Borwein, and Simon Plouffe, On the rapid computaion of various polylogarithmic constants, Mathematics of Computation 66 (1997) 903-913.

External LINKS
5 Trillion Digits of Pi - New World Record(Aug. 7, 2010)
Numerical approximations of π - Wikipedia

www.finetune.co.jp [Mail]