Transcript PPT - 明治大学
π
~計算法の変遷~
2006年2月17日
明治大学理工学部数学科
鎌田 伊織
吉本 清夏
円周率の計算法の歴史
• 円に内接・外接する正多角形の周長の利用
(アルキメデス(BC3C)・関 孝和(17C)・建部 賢弘(18C))
• 逆正接関数(Arctan)のTaylor展開の利用
(Leibniz(15C)・Sharp((17C)・Machin(18C))
• AGM(算術幾何平均:Arithmetic Geometric Mean)
(Salamin-BrentによるGauss-Legendre法(1976)・Borwein法)
• DRM(分割有理数化:Divide and Rationalize Method)
(後 保範(1998) & 金田 康生(2002))
Arctan級数とは?
f (x) arctanx の両辺を微分すると,
1
2 n
f ( x)
(
x
)
2
1 x
n0
x 1
項別積分することで, ArctanのTaylor展開が得られる.
(1) n1 2n1
arctan
x 1
2
n
1
n1
実は, x=1でも成立し, Leibniz級数を得る.
4
x=
1
3
1
1 1 1 1
3 5 7 9
(1)
とすると, Sharpの公式が得られる.
6
arctan
1
3
(2)
・ (1)は収束速度が非常に遅く, 効率が悪い
・ (2)は(1)よりはずいぶん速くなる
・ 長い年月をかけて多くの人々が競って計算するようになった
Arctan級数の公式1/2
• Machin
(1706)
1
1
4 arctan arctan
4
5
239
(当時の最高記録100桁を求め, その後
多くの人に利用された. )
• 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
(おそらく3項公式で最高効率)
(1985年にフェルトンによって1万21桁を得た)
1
1
1
1
12arctan 20arctan 7 arctan
24arctan
4
38
57
239
268
• 高野喜久雄 (1982)
(2002年に世界一の1兆2400億桁を達成した公式!)
1
1
1
1
12 arctan 32 arctan 5 arctan
12 arctan
4
49
57
239
110443
実験結果1/3
部分和の項数 j と 誤差log10(π- S( j ))との関係
*S( j )=各公式におけるTaylor
展開の第 j 部分和
*横軸 =項数 j
*縦軸 =πとS( j )の誤差の
log 10 をとったもの
・ ■ マチンの公式
・ ■ オイラーの公式①
・ ■ オイラーの公式②
・ ■ ガウスの公式①
・ ■ ガウスの公式②
・ ■ 高野喜久雄の公式
実験結果2/3
<グラフからわかったこと>
• 傾きは, それぞれの公式において最も収束速度が遅い項による
(収束速度はArctan(x)の|x|が大きいほど遅くなる)
• Arctan(x)の|x|が大きい項をもつ公式ほど, |傾き|が小さい
• 傾きと|x|のlog10 値を比べてみる
Arctanxの|x|
傾き
log10 ( x)
1/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(x )の2倍
になっている
* 次に, この関係を調べる
実験結果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
44
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 ) 44
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の2倍が傾きになると考えられる.
算術幾何平均とは・・・
a b
・ , ab のことを昔は算術平均・幾何平均と呼んでいた.
2
・ a0 a, b0 b とし、
an bn
an1
, bn1 anbn
2
(n 0,1,2)
・数列{an}と{ bn}は共通の極限に収束する.
・この値を a と b の算術幾何平均(arithmetic-geometric
mean)と呼び、 M (a, b) で表す.
背景1/3
<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
背景2/3
<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 2n1 cn2 ) I (a, b)
2 M ( a, b)
n 0
が成り立つ.ただし、
cn | an2 bn2 |
背景3/3
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
n1 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 の差が所要桁以上になるま
で計算する.
an1 bn1
an
, bn an1bn1 , tn tn1 2n1 (an an1 )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
an1 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
1.0134E-0003
2
7.3763E-0009
3
1.8313E-0019
4
5.4721E-0041
5
2.4061E-0084
6
2.3086E-0171
7
1.0586E-0345
8
1.1110E-0694
9
-1.5022E-1000
10
-1.5022E-1000
<縦軸:log10 | log10( との誤差 )|
f i g.πとの誤差
横軸:回数n>
まとめ
• 今回は取り上げられなかったが、ボールウェインの
4次式では計算精度が4倍である
• 現在の世界記録は高野喜久雄の(Arctan)公式と
DRM法を使って,2002年11月に後 保範氏&金
田 康生氏によって計算された約1兆2400億桁!
• 今後もコンピュータの発達によりπの計算記録の樹
立は変わってくると考えられる
参考文献
• (1) E.ハイラー, G.ワナー 著, 蟹江幸博 訳, 解析教程 上, シュ
プリンガー・フェアラーク東京 (1997)
• (2) 桂田祐史, πノート, 明治大学数学科助教授 (2004)
• (3) 清水康生, πの数値解析, 明治大学数学科 2003年度卒業
研究レポート (2004)
• (4) ペートル・ベックマン 著, 田尾陽一, 清水韶光 訳, πの歴史,
蒼樹書房 (1973)
• (5) 数学文化, Vol.1, 日本数学協会 (2003)
• (6)金田康正, πのはなし, 東京書籍 (1991)
• (7) 梅村浩, 楕円関数論, 東京大学出版会 (1999)
• (8)ドゥラエ・ジャン=ポール(Jean-Paul Delahaye)著, 畑政義
訳, π-魅惑の数, 朝倉書店 (2001)
御静聴ありがとうございました。