計算法の背景

(細かい話になるけれど、一応書いておく。)

「発掘」された方法の基礎となっているのは、 19世紀数学の華とも呼ばれる楕円関数論からの次の二つの事実である。

(i)

第一種完全楕円積分

$$K(k):=\int_0^{1}\frac{\Dx}{\sqrt{(1-x^2)(1-k^2x^2)}\;} =\int_0^{\pi/2}\frac{\D\theta}{\sqrt{1-k^2\sin^2\theta\;}\;}$$   $$\mbox{($0\le k<1$)}$$

と第二種完全楕円積分

$$E(k):=\int_0^{1}\frac{\sqrt{1-k^2x^2}}{\sqrt{1-x^2}} \Dx =\int_0^{\pi/2}\sqrt{1-k^2\sin^2\theta\;} \D\theta$$   $$\mbox{($0\le k\le $)}$$

の間に成り立つ ルジャンドルLegendre の関係式 (Legendre's relation) と呼ばれる

$$K(k)E(k')+K(k')E(k)-K(k)K(k')=\frac{\pi}{2} \mbox{(ただし $k':=\sqrt{1-k^2}$)}$$ (5)

という式の、特に $k=1/\sqrt{2}$ の場合の

$$2K\left(\frac{1}{\sqrt{2}}\right)^2 -E\left(\frac{1}{\sqrt{2}}\right) K\left(\frac{1}{\sqrt{2}}\right) =\frac{\pi}{2}.$$ (6)

(ii)

$K(k)$, $E(k)$ はいわゆる算術幾何平均 (arithmetic geometric mean, AGM) アルゴリズムで計算できる。 第一種完全楕円積分、 第二種完全楕円積分の二変数版 $I(a,b)$, $J(a,b)$ を

$$I(a,b):=\int_0^{\pi/2}\frac{\D\theta} {\sqrt{a^2\cos^2\theta+b^2\... ...quad J(a,b):=\int_0^{\pi/2} \sqrt{a^2\cos^2\theta+b^2\sin^2\theta\;} \D\theta$$

で定めるとき、

$$I(a,b)M(a,b)=\frac{\pi}{2},\quad J(a,b)=\left(a-\sum_{n=0}^\infty 2^{n-1}c_n^2\right)I(a,b)$$

が成り立つ。 ここで $M(a,b)$ は $a$, $b$ の算術幾何平均と呼ばれる量で、 (3) で定義される数列 $\{a_n\}$, $\{b_n\}$ の 共通の極限として定義される:

$$M(a,b):=\lim_{n\to\infty}a_n=\lim_{n\to\infty}b_n.$$

また数列 $\{c_n\}$ は

$$c_n:=\sqrt{\left\vert a_n^2-b_n^2\right\vert}$$

で定義される。 容易に

$$K(k)=I\left(1,k'\right),\quad E(k)=J\left(1,k'\right), \quad k'=\sqrt{1-k^2}$$

であることが分かるので、

$$K\left(\frac{1}{\sqrt{2}}\right)=\frac{\pi}{2}\frac{1}{M(1,1/\sqr... ...\pi}{2}\left(1-\sum_{n=0}^\infty 2^{n-1}c_n^2\right)\frac{1}{M(1,1/\sqrt{2})} .$$ (7)

(7) を (6) に代入して整理すると

$$\pi=\frac{2M\left(1,1/\sqrt{2}\right)^2}{1-\dsp\sum_{n=0}^\infty 2^n c_n^2}.$$

ここで述べたことはすべて19世紀の段階で分かっていたわけであるが、 $\pi$ の計算法として、 このアルゴリズムが 1976 年という時点で初めて注目 (「発掘」) されるようになった理由は、 この方法が最初から長い桁の数の掛け算、平方根を必要とするため、 $\arctan$ の級数展開を利用する方法と比べて不利だと考えられたせいであろう。 1971 年の Strassen と Schonhage による高速乗算法 (これは高速 Fourier 変換に基づいている) の発見により、 その立場が逆転してしまった。 二つの $n$ 桁の数の積の計算が $O(n^2)$ ではなく、 $O(n\log n)$ 程度の計算量で実行可能というのは発見当時は 驚くべきことだったと思われる。