Document

π
～計算法の変遷～
2006年２月１７日
明治大学理工学部数学科
鎌田伊織
吉本清夏
円周率の計算法の歴史
• 円に内接・外接する正多角形の周長の利用
(アルキメデス（BC3C）・関孝和（１７C）・建部賢弘（１８C)）
• 逆正接関数(Arctan)のTaylor展開の利用
（Leibniz（１５C)・Sharp（（１７C)・Machin（１８C))
• ＡＧＭ(算術幾何平均：Arithmetic Geometric Mean)
(Salamin-BrentによるGauss-Legendre法（１９７６）・Borwein法）
• ＤＲＭ（分割有理数化:Divide and Rationalize Method）
（後保範（１９９８）＆金田康生（２００２））
Arctan級数とは？
f ( x)  arctan x の両辺を微分すると,

1
2 n
f ( x) 

(

x
)
2
1 x
n 0

 x  1
項別積分することで, ArctanのTaylor展開が得られる.
arctan  


n 1
(1) n 1 2 n 1

 x  1
2n  1
実は, ｘ＝1でも成立し, Leibniz級数を得る.

4
x＝
1
3
 1
1
1
1
1




3
5
7
9
（１）
とすると, Sharpの公式が得られる.

6
 arctan
1
3
（２）
・（１）は収束速度が非常に遅く, 効率が悪い
・（２）は（１）よりはずいぶん速くなる
・長い年月をかけて多くの人々が競って計算するようになった
Arctan級数の公式1/2
• Machin
(1706)

1
1
 4 arctan  arctan
4
5
239
（当時の最高記録１００桁を求め, その後
多くの人に利用された. ）
• Euler
(1737）
(1755)

1
1
 arctan  arctan
4
2
3

1
3
 5 arctan  2 arctan
4
7
79
Arctan級数の公式2/2
• Gauss
(1863)

1
1
1
 12 arctan
 8 arctan
 5 arctan
4
18
57
239
(おそらく３項公式で最高効率）
（１９８５年にフェルトンによって１万２１桁を得た）

1
1
1
1
 12 arctan
 20 arctan
 7 arctan
 24 arctan
4
38
57
239
268
• 高野喜久雄 (1982)
（２００２年に世界一の１兆２４００億桁を達成した公式！）

1
1
1
1
 12 arctan
 32 arctan
 5 arctan
 12 arctan
4
49
57
239
110443
実験結果1/3
部分和の項数ｊと誤差log１０（π－ S（ｊ））との関係
＊S（ｊ）＝各公式におけるTaylor
展開の第ｊ部分和
＊横軸＝項数ｊ
＊縦軸＝πとS( j )の誤差の
log 10 をとったもの
・ ■ マチンの公式
・ ■ オイラーの公式①
・ ■ オイラーの公式②
・ ■ ガウスの公式①
・ ■ ガウスの公式②
・ ■ 高野喜久雄の公式
実験結果2/3
＜グラフからわかったこと＞
• 傾きは, それぞれの公式において最も収束速度が遅い項による
（収束速度はArctan（ｘ）の｜ｘ｜が大きいほど遅くなる）
• Arctan（ｘ）の｜ｘ｜が大きい項をもつ公式ほど, |傾き|が小さい
• 傾きと｜ｘ｜のlog10 値を比べてみる
Arctanxの｜ｘ｜
傾き
log10 ( x)
１/2
－0.6
－0.301
1/5
－1.4
－0.698
1/7
－1.7
－0.845
1 / 18
－2.5
－1.255
1 / 38
－3.16
－1.579
1 / 49
－3.24
－1.690
＊傾きは log10（ｘ）の２倍
になっている
＊次に, この関係を調べる
実験結果3/3

1
1
＜Machinの公式（  4 arctan  arctan
）を例にとって考える＞
4
5
239
x
1
1
,y
5
239
とおき, Taylor展開をすると,

  (1) k 1 2 k 1
(1) k 1 2 k 1 
  44
x

y

2
k

1
2
k

1
k 1
 k 1

となる. また, j での部分和は,
j
 j (1) k 1 2 k 1
(1) k 1 2 k 1 
S ( j )  44
x

y

2
k

1
2
k

1
k 1
 k 1

よって, πと j での部分和との誤差は,
 ( 1) j 2 j 1 ( 1) j 2 j 1 
x 2 j 1
  S ( j )  4 4
x

y
 ≒ 16
2 j 1
2 j 1
 2 j 1

対数をとると
log10   S  j  ≒log10 16  log10 (2 j 1)  (2 j 1) log10 x
よって,log10xの２倍が傾きになると考えられる.
算術幾何平均とは・・・
ab
・ 2 , ab
のことを昔は算術平均・幾何平均と呼んでいた．
・ a0  a, b0  b とし、
an  bn
an 1 
, bn 1  anbn
2
(n  0,1,2)
・数列｛an｝と｛ bn｝は共通の極限に収束する．
・この値を a と b の算術幾何平均(arithmetic-geometric
mean)と呼び、 M (a, b) で表す．
背景１/３
＜Legendreの関係式＞
（ⅰ）第一種完全楕円積分
dx
1
K (k ) : 
0
(1  x )(1  k x )
2
2 2

2
0

d
1  k sin 
2
2
(0  k  1)
第二種完全楕円積分
1
1 k x
0
1 x2
E (k ) : 
2 2

2
0
dx   1  k 2 sin 2  d (0  k  1)
の間に成り立つLegendreの関係式

K (k ) E (k )  K (k ) E (k )  K (k ) K (k ) 
2
k 
1 k 2

背景２/３
＜Gauss,“the fundamental limit theorem”＞
第一種完全楕円積分、第二種完全楕円積分の二変数版を

I (a, b)   2
0

d
a 2 cos2   b 2 sin 2 
, J (a, b)   2 a 2 cos2   b2 sin 2  d
0
で定めると、


1
I ( a, b)  
, J (a, b)  (a   2 n 1 cn2 ) I (a, b)
2 M ( a, b)
n 0
が成り立つ．ただし、
cn  | an2  bn2 |
背景３/３
1
• Legendreの関係式で k 
の時、
2
 1   1 
 1  
2K 
 E
  K
  ①
2
 2  2
 2
2
)
K (k )  I (1, k ), E (k )  J (1, k ) (k   1  k 2 であるので、

1
 1  
 1  
 1 
n 1 2
K
, E
 
  (1   2 cn ) M 1,
 ②
2
 2  2 M 1, 1   2  2

n 0


2

①に②を代入して、

 1 
2 M 1,

2



1   2 n cn2
n 0
2
ガウス・ルジャンドルの公式
1
1
a0  1, b0 
, t0 
4
2
として、
以下の反復式を an と bn の差が所要桁以上になるま
で計算する．
an 1  bn 1
an 
, bn  an 1bn 1 , t n  t n 1  2 n 1 (an  an 1 ) 2
2
所要桁になったら円周率は、

と求められる．
≒
(an  bn ) 2
4tn
 ＝
・
・
 1 
2 M 1,

2


2

1   2 n 1 cn2
n 0
と
≒
(an  bn )  an  bn 
 1 
2


  an 1  M 1,
4
2
 2 

2
2
1 n k 1
t n    2 (ak  ak 1 ) 2
4 k 1
また、
(a n  bn ) 2
4t n
関係
2
ガウス・ルジャンドルの公式より

tn 
1   2 n cn2
n 0
2
1 1 n k 2
   2 ck
4 2 k 1
1 

c0 

2 

1 1 n k 2
    2 ck 
2 2  k 0

n
1
   2 k 1 ( ak2 bk2 )
4 k 1
n
1 1 2

   c0   2 k ck2 
2 2
k 1

n
1
   2 k 1 ( ak ak 1 ) 2
4 k 1
誤差の減少の速さ
n
πとの誤差
１
1.0134E-0003
２
7.3763E-0009
３
1.8313E-0019
４
5.4721E-0041
５
2.4061E-0084
６
2.3086E-0171
７
1.0586E-0345
８
1.1110E-0694
９
-1.5022E-1000
10
-1.5022E-1000
＜縦軸：log10 | log10（ との誤差
f i g.πとの誤差
）|
横軸：回数ｎ＞
まとめ
• 今回は取り上げられなかったが、ボールウェインの
４次式では計算精度が４倍である
• 現在の世界記録は高野喜久雄の（Ａｒｃｔａｎ）公式と
ＤＲＭ法を使って，２００２年１１月に後保範氏＆金
田康生氏によって計算された約１兆２４００億桁！
• 今後もコンピュータの発達によりπの計算記録の樹
立は変わってくると考えられる
参考文献
• （１）Ｅ．ハイラー, Ｇ．ワナー著, 蟹江幸博訳, 解析教程上, シュ
プリンガー・フェアラーク東京（１９９７）
• （２）桂田祐史, πノート, 明治大学数学科助教授（２００４）
• （３）清水康生, πの数値解析, 明治大学数学科２００３年度卒業
研究レポート（２００４）
• （４）ペートル・ベックマン著, 田尾陽一, 清水韶光訳, πの歴史,
蒼樹書房（１９７３）
• （５）数学文化, Ｖｏｌ.１, 日本数学協会（２００３）
• （６）金田康正, πのはなし, 東京書籍（１９９１）
• （７）梅村浩, 楕円関数論, 東京大学出版会（１９９９）
• （８）ドゥラエ・ジャン＝ポール（Ｊｅａｎ－ＰａｕｌＤｅｌａｈａｙｅ）著, 畑政義
訳, π-魅惑の数, 朝倉書店（２００１）
御静聴ありがとうございました。

Download Report