π
~計算法の変遷~
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)  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
実は, 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
 12 arctan
 8 arctan
 5 arctan
4
18
57
239
(おそらく3項公式で最高効率)
(1985年にフェルトンによって1万21桁を得た)

1
1
1
1
 12 arctan
 20 arctan
 7 arctan
 24 arctan
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
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の2倍が傾きになると考えられる.
算術幾何平均とは・・・
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) で表す.
背景1/3
<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

背景2/3
<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 |
背景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 
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
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)
御静聴ありがとうございました。
ダウンロード

PPT - 明治大学数学科ホームページへ