ちょっと変わった円周率の計算法(ついでに平方根常用対数の値も)

〜(有限のメモリで)計算できる桁数が最大となる方法〜5/23/09

 よく使われる Machin の公式ではなく

   π = 4 Arctan 1

で計算しますが,Arctan の級数展開は Gregory (, Leibniz) の級数ではなく

   Arctan x = {x / (x 2 + 1)} Σ {x 2 / (x 2 + 1)} n (2 n)!! / (2 n + 1)!!

を使います.ただし,Σ は n = 0 から ∞ までの総和で,

   (2 n)!! = 2 n n! ,
   (2 n - 1)!! = (2 n)! / (2 n)!! .

すなわち,

   π = 2 Σ (1/2) n (2 n)!! / (2 n + 1)!!

ですから,2 n の order で収束する Gregory の級数に対して,2 n の order で収束します.
 実際は,例えば n = 4 までの和の場合,

   π ≈ (((2 · 4/9 + 2) 3/7 + 2) 2/5 + 2) 1/3 + 2

として計算すると,計算が簡単になり,桁落ちも少なくなります.(1/2) n < 10 -15 となる自然数 n の最小値は 50 ですから(実際は 49 まででよい),JavaScript では
   p = 2;
   for (n = 49; n > 0; n--) p = p * n / (2 * n + 1) + 2
とするだけで,π の近似値が 16 桁得られます.実際に,

n = (から 0)まで  すると,π =

 この方法は 1979 年頃にプログラム電卓で可能な限り多数の桁を計算するために考案したもので,他の方法に比べて時間はかかりますが,2倍以上の桁数が計算できます(現代のコンピュータではメモリよりも計算時間の方が問題になるのでしょうが).当時,SHARP PC-1200 で小数点以下 77 桁(7 桁 × 11 個のメモリ)を 29'30" かかって計算しました.というわけで,

1 桁 (クリックする前に下の注意を)すると,

π = 3.

ここでは 10 桁単位に分けて計算していますので(徐算では余りに 10 10 をかけて下位の 10 桁に加える),JavaScript の整数(IEEE double precision 型の significand)の最大値は 2 53 - 1 ですから,正しく計算できる桁数 d は概算で,d < n log10 2 と 2 n · 10 10 < 2 53 より,135001 桁となります(JavaScript では計算時間の点で非現実的ですが).

 なお,

   π = 6 Arcsin (1/2) = 6 Σ {(1/2) 2 n + 1 / (2 n + 1)} (2 n - 1)!! / (2 n)!!

で計算すると,2 n · 4 n の order で収束しますので上記の方法より収束は速いのですが,各項の演算数はほぼ倍になります.例えば n = 4 までの和の場合,

   π ≈ (((3/9 · 7/8 · 1/4 + 3/7) 5/6 · 1/4 + 3/5) 3/4 · 1/4 + 3/3) 1/2 · 1/4 + 3

と計算し,π の近似値を 16 桁得るには
   p = 3 / (2 * 22 + 1);
   for (n = 22; n > 0; n--) p = p * (2 * n - 1) / (2 * n * 4) + 3 / (2 * n - 1)
となります.実際に,

n = (から 0)まで  すると,π =

 以下の平方根と常用対数の値も当時 PC-1200 で計算した方法です.


 平方根の計算には

   (a 2 - b) 1/2 = a - a Σ (b / a 2) n (2 n - 3)!! / (2 n)!!

を使います(Σ は n = 1 から ∞ までの総和,a > 0,|b| < a 2,(-1)!! = 1).例えば,

   2 √2 = 3 - 3 Σ (1/9) n (2 n - 3)!! / (2 n)!!

ですから,JavaScript では

   r = 3;
   for (n = 14; n > 0; n--) r = r * (2 * n - 3) / (18 * n) + 3;
   r /= 2

とすると,√2 の近似値が 16 桁得られます.√3 は a = 2, b = 1 とするより,a = 7, b = 1(4 √3)とすると計算が速くなります(16 桁のとき,前者の和は n = 24 まで,後者は n = 9 まで).

 a b /
2312
3714
5914
6512
7813
1010103

として,

1 桁 すると,

正しく計算できる桁数 d が最も少ないのは √5 で,d < n log10 81 と 162 n · 10 10 < 2 53 より 10001 桁,最も多いのは √2 で,d < n log10 9 と 18 n · 10 10 < 2 53 より 47001 桁となります.√10 は 3 √10 より 6 √10(a = 19, b = 1)の方が収束は速いのですが,正しく計算できる桁数を優先しました(6 √10 の場合 3200 桁程度であるのに対して,3 √10 の場合は 45001 桁以上).


 自然対数の計算には

   log((1 + x) / (1 - x)) = 2 x Σ x 2 (n - 1) (2 n - 3)!! / (2 n - 1)!!

を使います(Σ は n = 1 から ∞ までの総和,|x| < 1).log 2 は x = 1/3 として(log 8 より),log 3 は log 9 = log (9/8) + log 8 より,log 10 は log (10/9) + log 9 として,JavaScript では

   l8 = 2;
   for (n = 15; n > 1; n--) l8 = l8 * (2 * n - 3) / (18 * n - 9) + 2;
   l9 = 2;
   for (n = 6; n > 1; n--) l9 = l9 * (2 * n - 3) / (578 * n - 289) + 2;
   l9 = l9 / 17 + l8;
   l10 = 2;
   for (n = 6; n > 1; n--) l10 = l10 * (2 * n - 3) / (722 * n - 361) + 2;
   l10 = l10 / 19 + l9;
   l2 = l8 / l10 / 3;
   l3 = l9 / l10 / 2

とすると,log10 2 が 15 桁,log10 3 が 16 桁得られます.ですが,√10 の場合と同じ理由で,log 3 は x = 1/2 として,log 10 は log (5/4) + log 8 として,

0 桁  すると( log10 2 だけ),

log10 2 = 0.

log10 3 = 0.

正しく計算できる桁数は √5 の場合と同じ計算(log 10 が最も少ない)で 10000 桁となります.なお,最後の徐算では,除数 log 10 も多桁ですから 5 桁単位に分け直して計算していますので,そこでの桁数制限はありません.


Home へ戻る