Friday, 28th April, 2000.
[1] 逆正接函数系
at(k) = arctan(1/k) と置く。但し tan(arctan x) = x である。すべて π /4 の値。
(1) John Machin (1680--1751) (1706): 4at(5)−at(239).
(2) C. Hutton (, L. Euler, G. Vega) (1776): 2at(3) + at(7).
(3) J. Hermann (, C. Hutton) (1706): 2at(2)−at(7).
(4) L. Euler (, C. Hutton) (1738): at(2) + at(3).
(5) Carl Friedlich Gauß (30th April 1777--23rd February 1855) (1863): 12at(18) + 8at(57)−5at(239).
(6) S. Klingenstierna (1730): 8at(10)−at(239)−4at(515).
(7) F. C. M. Störmer (1896): 4at(5)−2at(478) + at(54608393).
(8) J. W. Wrench. Jr. (1938): 5at(7) + 4at(53) + 2at(4443).
(9) F. C. M. Störmer (1896以下同じ): 4at(5)−at(240)−at(57361).
(10) 4at(5)−at(238) + at(56883).
(11) 4at(5)−at(241)−at(28800).
(12) 4at(5)−at(237) + at(28322).
(13) 6at(8) + 2at(57) + at(239). (普通これを Störmer の公式と呼ぶようである)
(14) 4at(5)−2at(577)−at(1393).
(15) 4at(5)−at(252)−at(4633).
これらのうち, 実際に計算機による桁数の記録樹立で用いられたものは (1), (5), (6), (7) だけのようである。 (1), (5) などによる計算は, 今のパソコン程度なら容易に 10,000 桁程度は計算できる。私の JAVA の page にもある。 勿論, arctan の計算には通常の Taylor-Maclaurin の展開公式を用いる (super computer の計算の時には, 掛算に FFT 即ち fast Fourier transformation と呼ばれる計算技法が用いられるが, 我々が計算するときにはそこまではやらない)。
因みにこれらの系列は Euler の 1737 年の arctan [1/ (p+q)] + arctan [q / (p2 + pq + 1)] = arctan [(1/ (p+q) + q / (p2 + pq +1)) / (1 - q / ((p+q) (p2 + pq + 1)))] = at(p) や, arctan a + arctan b = arctan ((a + b) / (1 - ab)) がもとになっている。
[2] Salamin-Brent の algorithm
金田等が 10 億桁求めた方法を書く。(一般にはもっと別のやり方もあるのである)
A = 1, B = 1/√2, T = 1/4, X = 1を初期値とする。A,Bの差が必要とする精度より大きい間, 次の { } 間の計算を繰り返す (勿論 A, B, T, Y の各々の値は, 必要とする精度以上に計算しておく必要がある。)
{
Y = A;
A = (A + B) / 2;
B = √(B*Y);
T = T−X*(Y−A)*(Y−A);
X = 2*X;
}
その結果, π = (A+B)*(A+B)/(4*T) となる。又, これで分かることは, 実は案外知られていないが π 以上に √2 が充分な精度で計算されているということである。 (実際に √2 も 1/√2 も 16 億桁以上計算されている)
不思議なことに, 次の方法よりも計算量は多いのに, 一億桁もの計算であると計算時間の方は短いという事実がある。
[3] J. M. Borwein, P. B. Borwein 兄弟の公式 (Canada, ダルハウジー大学)
二次, 三次, 四次, 七次の収束の公式があるが, D. Bailey (USA) は四次の収束公式を主計算で使用, 二次の収束公式を検証計算で用いた。金田, 田村 (1988。1987の時は他に久保, 小林, 花村) は四次の収束公式を検証の計算に用いている。
(1) 四次の収束を示す公式。
初期値としてa0 = 6−4×21/2, y0 = 21/2−1として, 以下
yk+1 = (1−(1−yk4)1/4) / (1 + (1−yk4)1/4), ak+1 = ak(1 + yk+1) 4−22k+3yk+1(1 + yk+1 + yk+12)
として, 求める精度以上の精度で計算しておくと, π = 1 / akとなる。実際, ak = 1/πは, 10 億桁の精度で保存されている。
(2) 二次の収束を示す公式。
今度は a0 = √2, b0 = 0, p0 = 2 + √2 として以下,
ak+1 = (ak1/2 + ak−1/2) / 2, bk+1 = ak1/2(1 + bk) / (ak + bk), pk+1 = pkbk+1(1 + ak+1) / (1 + bk+1)
として, 求める精度以上の精度で計算しておくと π = pk となる。
さて, 見てみると分かるのは, [2], [3] 共に, 平方根, 四乗根, 逆数の計算が必要になって来る。それは次の Newton 法で行われる。以下, 初期値としては通常の (C 言語で言うところの) float 或いは double で計算した値を X0 として反復計算する。
(1) 平方根
通常のものではなく f(x) = x−2−Aに対する Newton 法, Xi+1 = Xi(3−AXi4) / 4 を 用いると A−1/2が求まるので, √A = A×A−1/2として求める。
(2) 四乗根
先ず f(x) = x−4−Aに対する Newton 法 Xi+1 = Xi(5−AXi4) / 4 を用いると A−1/4 が求まるので, A1/4= A×(A−1/4)3 として求める。
(3) 逆数の計算
これは通常の f(x) = x−1−Aに関する Newton 法 Xi+1 = Xi(2−AXi) を用いて求める。
[4] Chudonovsky の公式
D.V. Chudonovsky (チュドノフスキー, Чудновский, USA), Grigoriĭ V. Chudonovsky の兄弟が 1989 年以降用いている公式
1/π = (6541681608/6403203/2)Σn=0∞[((13591409/545140134) + n)(-1)n(6n)!/((3n)! (n!)36423203n)].
以下の素因数分解が知られている。
6541681608 = 23･33･7･11･19･127･163,
640320 = 26･3･5･23･29,
13591409 = 13･1045493,
545140134 = 2･32･7･11･19･127･163.
この公式は項を一つ増やすと約 14 桁ずつ精度があがる。
[5] Ramanujan の公式
1985 年に W. Gosper (USA) が 17,526,000 桁の計算で使用した公式。
1/π = (2√2/9801) Σn=0∞[((4n)!(1103 + 26390n))/((4･99)n(n!))4].
前へ
次へ
円周率の目次へ
Miscellaneous の index へ
HOME