Transcript PPT - 明治大学

π
～計算法の変遷～
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)  arctanx の両辺を微分すると,

1
2 n
f ( x) 

(

x
)
2
1 x
n0

 x  1
項別積分することで, ArctanのTaylor展開が得られる.

(1) n1 2n1
arctan  

 x  1
2
n

1
n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
 12arctan  8 arctan  5 arctan
4
18
57
239
(おそらく３項公式で最高効率）
（１９８５年にフェルトンによって１万２１桁を得た）

1
1
1
1
 12arctan  20arctan  7 arctan
 24arctan
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
1
1
x , y
5
239
とおき, Taylor展開をすると,
  (1) k 1 2k 1  (1) k 1 2k 1 
  44
x

y

2
k

1
2
k

1
k 1
 k 1

となる. また, j での部分和は,
j
 j (1) k 1 2k 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
・ , ab のことを昔は算術平均・幾何平均と呼んでいた．
2
・ 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の関係式＞
（ⅰ）第一種完全楕円積分
K (k ) : 
dx
1
0
(1 x2 )(1 k 2 x2 )

2
0

d
1 k 2 sin 2 
(0  k  1)
第二種完全楕円積分
1
1  k 2 x2
0
1  x2
E(k ) : 

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
a2 cos2   b2 sin 2 
, J (a, b)   2 a2 cos2   b2 sin 2  d
0
で定めると、


1
I ( a, b)  
, J (a, b)  (a   2n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 M 1, 1   2  2
 2
n 0
 2
①に②を代入して、

 1 
2M 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 , tn  tn1  2n1 (an  an1 )2
2
所要桁になったら円周率は、
≒
と求められる．
(an  bn )2
4tn

・
・
 1 
2M 1,

2


＝
2

1   2 n 1 cn2
n 0
と
≒
(an  bn )2  an  bn 
 1 
2


  an1  M 1,
4
2


 2
2
1 n k 1
tn    2 (ak  ak 1 )2
4 k 1
また、
(a n  bn ) 2
4t n
関係
2
ガウス・ルジャンドルの公式より

tn 
1   2n 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

1 n k 1 2 2
   2 (ak bk )
4 k 1
1 1 2 n k 2
   c0   2 ck 
2 2
k 1

1 n k 1
   2 (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.πとの誤差
横軸：回数ｎ＞
まとめ
• 今回は取り上げられなかったが、ボールウェインの
４次式では計算精度が４倍である
• 現在の世界記録は高野喜久雄の（Ａｒｃｔａｎ）公式と
ＤＲＭ法を使って，２００２年１１月に後 保範氏＆金
田 康生氏によって計算された約１兆２４００億桁！
• 今後もコンピュータの発達によりπの計算記録の樹
立は変わってくると考えられる
参考文献
• （１） Ｅ．ハイラー, Ｇ．ワナー 著, 蟹江幸博 訳, 解析教程 上, シュ
プリンガー・フェアラーク東京 （１９９７）
• （２） 桂田祐史, πノート, 明治大学数学科助教授 （２００４）
• （３） 清水康生, πの数値解析, 明治大学数学科 ２００３年度卒業
研究レポート （２００４）
• （４） ペートル・ベックマン 著, 田尾陽一, 清水韶光 訳, πの歴史,
蒼樹書房 （１９７３）
• （５） 数学文化, Ｖｏｌ.１, 日本数学協会 （２００３）
• （６）金田康正, πのはなし, 東京書籍 （１９９１）
• （７） 梅村浩, 楕円関数論, 東京大学出版会 （１９９９）
• （８）ドゥラエ・ジャン＝ポール（Ｊｅａｎ－Ｐａｕｌ Ｄｅｌａｈａｙｅ）著, 畑政義
訳, π-魅惑の数, 朝倉書店 （２００１）
御静聴ありがとうございました。