算術幾何平均による円周率πの計算公式

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

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

【ガウス・ルジャンドル(Gauss-Legendre)の２次の収束公式によるプログラム例】

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) は x², sqrt(x) は x^1/2)

この公式は２次の収束を示し、ｎ桁の答えが log₂(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)の４次の収束公式によるプログラム例】

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) は x², pow4(x) は x⁴, sqrt(x) は x^1/2, qroot(x) は x^1/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 a₀ = 355.0/113.0;

/* a_i+1 = a_i + sin(a_i) は π に収束する。
 * a_i+1 = a_i - tan(a_i) でも同じ。
 */ 
while(abs(h = sin(a)) > epsilon) {
    a += h;
}

この公式は３次の収束を示し、繰り返しごとの誤差は次のようになる。

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月のたった１台のデスクトップＰＣを使った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