円周率の数値計算

電気通信大学オープンキャンパス

緒方 秀教
電気通信大学

大学院情報理工学研究科 情報・ネットワーク工学専攻

2020 年 7 月 19 日（日）円周率の数値計算
緒方秀教

電気通信大学オープンキャン

はじめに数学においで，古代から現代まで多くの数学者が円周率 π の計算に挑戦してきた．
ここでは，その円周率計算の歴史を振り返り，同時に「数値計算」という学術分野の一端の紹介も兼ねたい．
円周率 π の計算法は，大まかに次のように大別される．
1

正多角形の周長により評価する（古代，和算）．

2

arctan x のグレゴリー級数による計算法．

3連分数．

4

算術幾何平均．

5

ラマヌジャンの公式．

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

古代の円周率紀元前 2000 年頃のリンドパピルス．

π=

1692= 3.1604 . . .

旧約聖書：π = 3 を示す記述：
彼はまた海を鋳（い）で造った．縁から縁まで 10 キュッピットであって，周囲は円形をなし，高さは 5 キュッピットで，その周囲は網をもって測ると，30 キュッピットあった．
（歴代志下，4.2）
アルキメデス (Archimedes)：
円に内接／外接する正 96 角形の周長より，3.1410 . . . = 3

110< π < 3 = 3.1427 . . .
717
円周率の数値計算

緒方秀教

電気通信大学オープンキャン

古代の円周率

古代における円周率 π の計算は，正多角形の周長により評価する方法が主流であった．

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

古代の円周率

古代における円周率 π の計算は，正多角形の周長により評価する方法が主流であった．
（余談）2003 年，東京大学入試問題
円周率が 3.05 より大きいことを証明せよ．
多くの受験参考書等に載っている解答は，円に内接する正 12 ／ 8 角形の周長により評価する方法である．

円周率の数値計算
緒方秀教電気通信大学オープンキャン

和算：関孝和の円周率計算
関孝和 (1640?–1708) の円周率計算法
単位円に内接する正 2n 角形の一辺の長さを an とすれば，

π = lim 2n−1 an .
n→∞

an に対する漸化式を求める．
p
OD = 1 − (an /2)2 ,AD = 1 − ODp
= 1 − 1 − (an /2)2 ,AC

B
D

CD = an /2,

an+1 2 = AC2 = AD2 + CD2
p= (1 − 1 − (an /2)2 )2 + (an /2)2p
= 2 − 4 − an 2 ,qp∴ an+1 = 2 − 4 − an 2 .

BCBA, ACO

円周率の数値計算

緒方秀教

O
正 2n 角形の一辺
正 2n+1 角形の辺単位円の中心

電気通信大学オープンキャン

和算：関孝和の円周率計算
円周率の近似数列 {πn }（関孝和）
qp√
a2 = 2, an+1 = 2 − 4 − an 2

πn = 2n−1

an → π( n = 2, 3, . . . ),( n → ∞ ).

これを使ってパソコンで π の近似値を計算した．

∗nπn

相対誤差

10203040
π の厳密値

3.14158 77252 77159 70062 88542 62702
3.14159 26535 85093 23106 89057 95336
3.14159 26535 89793 23398 03670 44950
3.14159 26535 89793 23846 26391 08648
3.14159 26535 89793 23846 26433 83280 . . .

2 × 10−6
2 × 10−121 × 10−181 × 10−24

関孝和は n = 17 まで計算したそうである．
円周率の数値計算
緒方秀教

電気通信大学オープンキャン

和算：関孝和の円周率計算とエイトケン加速
さらに，関孝和は次の事実に気づいた．

π15 = 3.14159 26487 76985 66948 51079 69277 . . . ,π16 = 3.14159 26523 86591 34580 35255 21058 . . . ,π17 = 3.14159 26532 88992 76527 19430 42174 . . . ,π17 − π16= 0.25000 00001 07716 23496 75673 20818 . . .
π16 − π15πn+1 − πn は公比≒ 0.25 の等比数列に近づいていくようだ！
関はこの事実から，π17 までの値を使ってさらに良い π の近似値を求めた．r = (π17 − π16 )/(π16 − π15 ) とおいて，

π = π16 + (π17 − π16 ) + (π18 − π17 ) + · · ·≃ π16 + (π16 − π15 )(r + r 2 + r 3 + · · · )r,= π16 + (π16 − π15 )1−r円周率の数値計算
緒方秀教

電気通信大学オープンキャン

和算：関孝和の円周率計算とエイトケン加速r1−r
π17 − π16
π16 − π15
= π16 + (π16 − π15 )
π17 − π161−
π16 − π15
(π16 − π15 )(π17 − π16 )π ≃ π16 −.
(π17 − π16 ) − (π16 − π15 )π ≃ π16 + (π16 − π15 )今の議論を一般化すると，πn よりもっと速く収束する
π の近似数列 πn′ が得られる：

πn′ = πn+1 −(πn+1 − πn )(πn+2 − πn+1 ),(πn+2 − πn+1 ) − (πn+1 − πn )円周率の数値計算
緒方秀教

電気通信大学オープンキャン

和算：関孝和の円周率計算とエイトケン加速
改良された近似数列 πn′ で π の近似値を求めた．

n10203038

πn （改良前）

相対誤差3.14158 77252 77159 70062 88542 62702
3.14159 26535 85093 23106 89057 95336
3.14159 26535 89793 23398 03670 44950
3.14159 26535 89793 23846 25749 89170

nπn′ （改良後）10203038厳密値

3.14159 26535 89938 19856 11615 72583
3.14159 26535 89793 23846 26435 15120
3.14159 26535 89793 23846 26433 83280
3.14159 26535 89793 23846 26433 83280
3.14159 26535 89793 23846 26433 83280 . . .

2 × 10−6
2 × 10−12
1 × 10−18
2 × 10−23
相対誤差

5 × 10−14
4 × 10−26
4 × 10−38
9 × 10−48

確かに，改良された数列 πn′ は πn より速く収束している．

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

和算：関孝和の円周率計算とエイトケン加速

エイトケン (Aitken) 加速
数列 {an } の極限 a = limn→∞ an を求めたいとき，{an } よりもっと速くa に収束する数列 {an′ } を下記で求め，極限 a を効率よく求める方法．

an′ = an −(an+1 − an )2.
an+2 − 2an+1 + an

日本では関孝和が 1680 年に発見した（当時は「増約術」と呼んだ）．
欧州では 1926 年にエイトケンが最初にエイトケン加速法を用いた．
1876 年に H. von Nägelsbach が最初に用いたという話も．

関がエイトケン加速を発見した事実は，Brezinski が再発見した．

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

和算で最初の円周率公式
建部賢弘 (1772)

1212 · 22
12 · 22 · 32
π2 = 9 1 ++++ ··· .
3·4 3·4·5·6 3·4·5·6·7·8

πn ：右辺を第 0 ∼ n 項で打ち切って計算された π の近似値．
nπn

相対誤差

4812162040100188

3.14148 50901 17183 96036 66075
3.14159 24555 12120 93228 34703
3.14159 26531 19234 49522 63975
3.14159 26535 88526 71739 89080
3.14159 26535 89789 56596 43702
3.14159 26535 89793 23846 26434
3.14159 26535 89793 23846 26434
3.14159 26535 89793 23846 26434

3 × 10−5
6 × 10−8
1 × 10−10
4 × 10−13
1 × 10−15
4 × 10−28
8 × 10−65
7 × 10−116円周率の数値計算
緒方秀教

電気通信大学オープンキャン

arctan x のグレゴリー級数による方法
17 世紀の微積分の発見により，様々な関数の無限級数展開が得られた．
arctan x ：tan の逆関数．y = arctan x ⇔ x = tan y .
arctan x に対するグレゴリー (Gregory) 級数
arctan x =∞X(−1)nn=02n + 1x 2n+1 = x −x3x5x7+−+ ··· .
357x = 1 を代入すると（tan(π/4) = 1 より arctan 1 = π/4 であるから），
1 1 1π= 1 − + − + ··· .
43 5 7
ただし，この級数は収束が遅く，円周率計算には使えない．
（式自体は美しく神秘的であるが…）円周率の数値計算
緒方秀教

電気通信大学オープンキャン

arctan x のグレゴリー級数による方法
マチン (Machine) の公式

11π= 4 arctan − arctan.
45239
Machine はこの等式において，Gregory 級数を有限桁で打ち切ってarctan を計算することにより，π を 100 桁計算した．
グレゴリー級数を第 0 ∼ n 項で打ち切って arctan を計算することにより得られる π の近似値を πn とする．

nπn

相対誤差

2468101270

3.14162 10293 25034 42504 68325
3.14159 26824 04399 51724 02598
3.14159 26536 23554 76199 55046
3.14159 26535 89835 84748 57007
3.14159 26535 89793 29474 73749
3.14159 26535 89793 23853 93246
3.14159 26535 89793 23846 26434

9 × 10−6
9 × 10−96 × 10−11
1 × 10−14
2 × 10−17
2 × 10−20
4 × 10−102

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

arctan x のグレゴリー級数による方法
同種の方法．
クリンゲンシュテルナ (Klingenstierna,1730)

π1
1
1= 8 arctan− arctan− 4 arctan.
410239515
Gauss(1863)
π1
1
1= 12 arctan+ 8 arctan− 5 arctan.
41857239
シュテルメル (Störmer,1896)π1
1
1= 6 arctan + 2 arctan+ arctan.
4857239
高野喜久雄 (1982)

1111
π= 12 arctan+ 32 arctan− 5 arctan+ 12 arctan
.
44957
239
110443円周率の数値計算
緒方秀教

電気通信大学オープンキャン

arctan x のグレゴリー級数による方法
グレゴリー級数による方法
有理数の四則演算だけでできる．
コンピュータが出現する以前の π の数値計算法．
高野の方法で π の近似値を計算した．πn の定義は前と同様．

nπn

相対誤差

2468

3.14159 26536 09251 91213 79438
3.14159 26535 89793 24014 20633
3.14159 26535 89793 23846 26436
3.14159 26535 89793 23846 26434···
3.14159 26535 89793 23846 26434

6 × 10−12
5 × 10−19
6 × 10−26
7 × 10−3330

8 × 10−108

意外と速く収束する．

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

連分数
arctan x の連分数展開
arctan x =x.
12 x 21+22 x 23+32 x 25+
42 x 27+.
9 + ..

√x = 1, 1/ 3 を代入することにより，π の連分数表示が得られる．
tan

π= 1,4π1より= √631πarctan √ = .
63tan

arctan 1 =π,4

緒方

秀教

円周率の数値計算

電気通信大学オープンキャン

連分数
π の連分数表示
4(1) π =,121+223+32

5+7+42.
9 + ..

√2 3(2) π =.
12 /3
1+22 /3
3+32 /3
5+42 /3
7+.
9 + ..
円周率の数値計算
緒方秀教

電気通信大学オープンキャン

連分数
(1)

πn ：π の連分数展開 (1) を第 0 ∼ n 項で打ち切った π の近似値(2)πn ：π の連分数展開 (2) を第 0 ∼ n 項で打ち切った π の近似値
(1)nπn相対誤差10203040

3.14159 26730 30334 59326 95704 02201
3.14159 26535 89793 67126 22736 56036
3.14159 26535 89793 23846 26529 75192
3.14159 26535 89793 23846 26433 83280

nπn

6 × 10−9
1 × 10−16
3 × 10−24
7 × 10−32

10203040

3.14159 26535 91088 45342 46022 38255
3.14159 26535 89793 23846 26481 38655
3.14159 26535 89793 23846 26433 83280
3.14159 26535 89793 23846 26433 83280

厳密値

3.14159 26535 89793 23846 26433 83280 . . .

(2)

相対誤差

4 × 10−13
2 × 10−24
6 × 10−36
2 × 10−47

連分数展開 (2) のほうが速く収束する（見た目は (1) の方が美しいが…）．
円周率の数値計算 電気通信大学オープンキャン
緒方秀教

算術幾何平均による方法
a, b > 0 の算術平均（相加平均），幾何平均（相乗平均）．
算術平均 :

1(a + b),2

幾何平均 :√ab.

算術平均・幾何平均をとる操作を繰り返して数列 {an }, {bn } を生成する：

a0 = a,1
an+1 = (an + bn ),2

b0 = b,pbn+1 = an bn ( n = 0, 1, 2, . . . ).

b0 ≤ b1 ≤ b2 ≤ · · · ≤ bn ≤ · · · ≤ an ≤ · · · ≤ a2 ≤ a1 ≤ a0 ,1an − bn ≤ n (a0 − b0 ) → 0 ( n → ∞ ).
2∴n → ∞ で an , bn は同じ値に収束する…算術幾何平均a, b の算術幾何平均 M(a, b) = lim an = lim bn .
n→∞円周率の数値計算
緒方秀教

n→∞

電気通信大学オープンキャン

算術幾何平均による方法

an , bn → M(a, b) の収束は極めて速い．
√（例）a = a0 = 1, b = b0 = 1/ 2 の場合：
n012
345

an
1.00000 00000 00000 00000 00000
0.85355 33905 93273 76220 04221
0.84722 49029 23494 15261 57738
0.84721 30848 35192 80650 97026
0.84721 30847 93979 08660 700030.84721 30847 93979 08660 64991

bn
0.70710 67811 86547 52440 08443
0.84089 64152 53714 54303 11254
0.84720 12667 46891 46040 36315
0.84721 30847 52765 36670 42981
0.84721 30847 93979 08660 59979
0.84721 30847 93979 08660 64991

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

算術幾何平均による方法
楕円積分の理論より，算術幾何平均による円周率の計算法が得られる．
√
a0 = 1, b0 = 1/ 2,pp
1
an+1 = (an + bn ), bn+1 = an bn , cn = an2 − bn2 (n = 0, 1, 2, . . .)2
算術幾何平均による円周率の表式

π=√2M(1, 1/ 2)2.
∞
X
n 22 cn
1−n=0

算術幾何平均による円周率の近似式π ≃ πN =2aN+1 2.
NXn 21−2 cnn=0

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

算術幾何平均による方法
πN → π は極めて収束が速い．
N

πN

相対誤差

0
1234562.91421 35623 73095 04880 16887 24209 69807 85696 71875
3.14057 92505 22168 24831 13312 68975 82331 17734 40237
3.14159 26462 13542 28214 93444 31982 69577 43144 37223
3.14159 26535 89793 23827 95127 74801 86397 43812 25504
3.14159 26535 89793 23846 26433 83279 50288 41971 14678
3.14159 26535 89793 23846 26433 83279 50288 41971 69399
3.14159 26535 89793 23846 26433 83279 50288 41971 69399···
3.14159 26535 89793 23846 26433 83279 50288 41971 69399

7 × 10−2
3 × 10−4
2 × 10−9
6 × 10−20
2 × 10−41
8 × 10−85
7 × 10−172

9

8 × 10−1001∗ π ≃ πN の誤差評価：
N+1

π − πN ≤π 2 2N+4 e−π2√,M(1, 1/ 2)2π − πN+1 ≤2−(N+1)(π − πN )2 .
π2各段で πN ≃ π の有効桁数が倍になっている．

円周率の数値計算

緒方秀教

電気通信大学オープンキャン

ラマヌジャン

ラマヌジャン (Srinivasa Aiyangar Ramanujan)

1887∼1920．インド出身．
港湾事務所の事務員の仕事をしながら，独学で数学の研究に勤しむ．
英国数学者ハーディ(Hardy) にその才能を見いだされて渡英，共同研究．
常人ではとても思いつかない摩訶不思議な多くの公式．
ラマヌジャンいわく，ナーマギリ女神が夢の中で教えてくれたと…

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

ラマヌジャン（郵便切手にもなっている）

画像は，MacTutor History of Mathematics(https://mathshistory.st-andrews.ac.uk/) より転載．
円周率の数値計算
緒方秀教

電気通信大学オープンキャン

ラマヌジャン (Ramanujan)

円周率の公式（ラマヌジャン, 1914）√ ∞18 X (4n)! 1103 + 26390n.
= 2π99(4n n!)4994nn=0∞X4(−1)n (4n)! 1123 + 21460n
1.
=π882 n=0 (4n n!)48822n
こんな公式，どうすれば思いつくのやら…？

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

ラマヌジャン (Ramanujan)

1 番目の公式で π を計算してみた．
第 0 ∼ n 項の部分和から求めた π の近似値を πn とする．
nπn

誤差

151050π の厳密値

3.14159 26535 89793 87799 89058 26306
3.14159 26535 89793 23846 26433 83280···
驚異的な収束の速さ！
3.14159 26535 89793 23846 26433 83280 . . .

2 × 10−16
2 × 10−48
1 × 10−88
3 × 10−408

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

Chudnovski の公式：ラマヌジャンタイプの公式
Chudnovski の公式 (1987)∞
X
12(−1)n (6n)!(13591409 + 545140134n)1.
=π(3n)!(n!)3 6403203n
6403203/2 n=0πn : 第 0 ∼ n 項の部分和から得られる π の近似値．
n

πn

相対誤差

12
31020
50
70

3.14159 26535 89793 23846 26433 83587 35068 84758 66345
3.14159 26535 89793 23846 26433 83279 50288 41971 67678
3.14159 26535 89793 23846 26433 83279 50288 41971 69399

1 × 10−286 × 10−43
3 × 10−57
1 × 10−156
1 × 10−298
3 × 10−724
7 × 10−1002

ラマヌジャンの公式よりもっと驚異的な速さで収束する！

円周率の数値計算
緒方秀教

電気通信大学オープンキャン

円周率計算の記録
年

1706
1844
1947
1949
1958
1961
1973
1982
2010
2019
2020

計算者

桁数

Archimedes
3
van Ceulen
34
Machin100
Dase205
Ferguson808（電子計算機時代）
Reitwiesner
2,037
Genuys
10,000Shanks, Wrench
100,000
Gouilloud, Bouyer
1,000,000
金田，吉野，田村
16,000,000···
Yee, 近藤
5,000,000,000,000
Yee, 岩尾
31,415,926,535,897
Yee, Mullican
50,000,000,000,000
円周率の数値計算
緒方

秀教

要した時間
一生

8 時間
70 時間
100 分9 時間
24 時間
30 時間
90 日
118 日
303 日

電気通信大学オープンキャン

終わりに
この資料を作るため，π のいろんな数値計算公式を計算した．中には初めて試した公式もあるが，驚くほど速い収束を示すものもあり，楽しい作業であった．
私は数値計算を専攻しているが，ここで紹介した諸公式で円周率を計算してみるのも，数値計算という分野に触れる方法のひとつであると思う．数値計算は実際にプログラミングして性能を実感するものである．
数値計算は実験数学である．

π を何万・何億・何兆桁したところで，実用には役に立たないと思うかもしれないが，円周率計算の歴史は数学・計算機科学の発展とともにある．その意味で，円周率の計算は重要な学術分野である．
なお，何百・何千桁もの数値計算は多倍長計算ライブラリが必要であるが，私は exflib というライブラリを愛用している．これは初心者にも使いやすく，パソコンレベルの計算機で手軽に長い桁数の数値計算が体験できる．使い方はネットにいくつか解説記事があるのでそれを参照されたい．
円周率の数値計算
緒方秀教

電気通信大学オープンキャン

参考文献・Web サイト

梅村浩：楕円関数論–楕円曲線の解析学，東京大学出版会，
2000 年．

P. Henrici: Applied and Computational Complex Analysis
Vol.2, John Wiley & Suns, 1977.
円周率.jp：http://円周率.jp/formula/ramanujan.html
建部賢弘の「綴術算経」(Imujii’s Page)：https://sites.google.com/site/seijiimura/home/16-he-suannoyan-jiu/03-jian-bu-xian-hongno-zhui-shu-suan-jing
コラム 円周率（江戸の数学）：
https://www.ndl.go.jp/math/s1/c4.html

円周率の数値計算
緒方秀教

電気通信大学オープンキャン