5
二つ目の例題として円周率の計算方法を紹介します。
6
円周率とは、任意の円の
円周
の長さをその円の直径
の長さで割った数です。7
円周率の計算
緑の旗をクリックすると猫が円周率の値を教えてくれるプログラムを作りなさい。
8
平方根同様に様々な計算方法が知られています。ここでは代表的な方法を概説します。
9
6.2.1 正多角形に関する漸化式
10
まず、古典的な計算方法として円に内接する正n角形を考えましょう。正n角形の外周の長さ
11
を円の直径で割った数は、nが大きくなるに従い、円周率に近づいて行きます。
12
正n角形の外周の長さLn が分かっているならば、正2n角形の外周の長さL2n は Ln から比
13
較的容易に計算できることが知られています。まず、この関係式を導いてみましょう。
14
1 B
C
O D
図6.10: 正n角形と正2n角形の関係
今、半径1の円に内接する正n角形のひとつの辺が図6.10の線分ADであると仮定します。そ うすると、ADの中間点Cを通る直線は辺を正確に二分しますから、線分AB、BDは正2n角形 の二つの辺です。ここで̸ ACOは直角ですから、三平方の定理から以下の式が成り立ちます。
|AC|2+|CO|2 = 1
|AC|2+|BC|2 = |AB|2 この2式を以下のように変形します。
|CO| =
√
1− |AC|2 (6.1)
|BC| =
√
|AB|2− |AC|2 (6.2)
以下は明らかに成り立ちます。
|BC|+|CO| = 1
この|BC|、|CO|に式(6.1)、式(6.2)を代入すると以下の通りです。
√
|AB|2− |AC|2+
√
1− |AC|2= 1 この式を|AB|を左辺に置くように変形すると以下の式を得ます。
|AB|=
√ 2−2
√
1− |AC|2
(a)プログラム
(b)実行途中結果
図6.11: 正多角形による円周率の計算プログラム
この式に、|AC|=|AD|/2を代入すると以下の式を得ます。
|AB|=
√ 2−
√
4− |AD|2
正n角形の外周の長さLn はn|AD|、正2n角形の外周の長さL2n は2n|AB|です。これを上の式 へ代入すると、最終的に以下の式を得ます。
L2n = 2n
√ 2−√
4−(Ln/n)2 (6.3)
さて、半径1の円に内接する正6角形についてL6= 6が成り立ちます。そこで、上の式を用い
1
て L6= 6から出発し、L12、L24、...を求めて行けば、次第に2πに近づいて行くはずです。そし
2
て最後にその値を2で割れば、それが円周率の近似値です。
3
上の方法で円周率を計算するプログラムが図6.11です。このプログラムでは、途中結果も表示
4
するようにしています。とは言え、猫は3桁しか教えてくれませんから、より精密な値を知るため
5
にステージ左上の変数の値（図6.11の左上）を順に見ていくことにします。
6
正 角形
正96角形 3.141032
4
正192角形 3.141452
5
正384角形 3.141558
6
正768角形 3.141584
7
正1536角形 3.141590
8
正3072角形 3.141592
9
正6144角形 3.141593
10
正12288角形 3.141593
11
正24576角形 3.141593
12
正49152角形 3.141593
13
正98304角形 3.141593
14
正196608角形 3.141593
15
正393216角形 3.141594
16
正786432角形 3.141592
17
正1572864角形 3.141609
18
正3145728角形 3.141587
19
正6291456角形 3.141674
20
正12582912角形 3.141674
21
正25165824角形 3.143073
22
正50331648角形 3.159806
23
正100663296角形 3.181981
24
正201326592角形 3.354102
25
正402653184角形 4.242641
26
正805306368角形 6.000000
27
正1610612736角形 0.000000
28
正3221225472角形 0.000000
29
円周率の真値 3.14159265... と比較してみてください。正6144角形から正196608角形の範囲
30
では真値に最も近い値が求められています。しかし、さらに辺の数を増やすと、逆に真値から離れ
31
ていきます。これは計算誤差によるものです。結果として、この古典的な計算方法では精度の高い
32
計算は期待できません。なお、計算誤差については8章で学びます。
33
6.2.2 モンテカルロ法
34
平方根の場合と同様に円周率の計算にもモンテカルロ法が適用できます。基本的なアイディア
35
は、図6.12の正方形の面積と4分円（4分割した円）の面積の比を確率的に計算する方法です。
36
今、図6.12の正方形の1辺の長さを1とします。このとき、その正方形の面積は1です。次に 4分円の面積は半径1の円の面積の1/4ですから、π/4です。今、0≤x≤1、0≤y≤1の矩形領 域の中の点(x, y)をランダムに求めたと仮定します。そうすると、その点が4分円の中に含まれる 確率（すなわちx2+y2≤1が成り立つ確率）は
4分円の面積 正方形の面積 =π/4
1 = π 4
となるはずです。よって、その確率を4倍すれば円周率が求められます。
37
(x,y)
0
0 1
x
図6.12: 正方形の面積と4分円の面積
(a)プログラム
(b)実行途中結果
図 6.13: モンテカルロ法による円周率の計算プログラム
(a)プログラム
(b)実行途中結果
図6.14: 逆三角関数による円周率の計算プログラム
この方法によるプログラムは図6.13の通りです。以下は実際に計算した結果です8。
1
3.124800 ... 1万回の場合
2
3.149440 ... 10万回の場合
3
3.141964 ... 100万回の場合
4
3.141280 ... 1000万回の場合
5
3.141208 ... 1億回の場合
6
モンテカルロ法が実用に堪えないことは既に平方根の計算で見た通りです。しかし、素朴に円周率
7
を計算できる点で興味深い手法です。
8
6.2.3 逆三角関数の公式
9
円周率を多数桁計算する方法として逆三角関数の性質を用いる方法が知られており、近年の多く
10
の円周率計算はこの方法を基礎にしています。
11
まず、以下の式が成り立つことを思い出しましょう。
arctan(1) = π 4
よってarctan(1)を何らかの方法で計算し、その値を4倍すれば円周率が求めることができます。
12
8確率計算ですから、実行の度に求められる値は異なります。
2.6 2.8 3 3.2
1 2 3 4 5 6 7 8 9 10
図 6.15: モンテカルロ法による円周率の計算プログラム
テイラー展開
（Taylor expansion）を用いると、逆正弦関数arctan(x)は 以下のように多項式近似できることが知られています。arctan(x) =
∑∞ i=0
(−1)i
2i+ 1x2i+1=x−x3 3 +x5
5 −x7 7 +x9
9 − · · · 特にx= 1のときには、
arctan(1) = 1 1−1
3+1 5 −1
7+1
9− · · ·= π 4
となります。この式はライプニッツの級数またはグレゴリーの級数と呼ばれており、式自体は15
1
世紀から知られていたようです。この式を愚直に計算すれば円周率に辿り着けるはずです。上の
2
式は i = 0, ...,∞ の無限級数和ですが、もちろん実際には有限の計算しかできません。そこで、
3
i= 0, ..., N までで計算を打ち切ることにし、和∑N
i=0(−1)i/(2i+ 1)を求めてみます。
4
図6.14はN = 10 までの和を求めるプログラムです。このプログラムでは計算の途中結果
5
∑k
i=0(−1)i/(2i+ 1)（k= 0,1, ...,9）も表示します。以下がその途中結果です。
6
4.0
7
2.666667
8
3.466667
9
2.895238
10
3.339683
11
2.976046
12
3.283738
13
3.017072
14
3.252366
15
3.041840
16
この10個の数値の系列をグラフ化したものが図6.15です。数値が円周率を挟んで上下に振れながら
17
収束していきます。しかし収束は速いとは言えません。というのも、級数和∑N
i=0(−1)i/(2i+ 1)の
18
各項の絶対値(−1)i/(2i+ 1)= 1/(2i+ 1)はiに反比例する値であるため、ゆっくりとした速度で
19
上の問題点を改良すべく、現在ではより計算に適した公式を用いています。たとえば以下は18 世紀に発見されたマチンの公式です。
4 arctan (1
5 )
−arctan ( 1
239 )
= π 4 この公式の左辺各項をテイラー展開で置き換えると以下の式を得ます。
4
∑N
i=0
(−1)i 2i+ 1
(1 5
)2i+1
−
∑N
i=0
(−1)i 2i+ 1
( 1 239
)2i+1
≈π 4
第1項は(1/5)2i+1に比例しますから、iの増加と共に急速に小さくなります。第2項(1/239)2i+1
3
は第1項よりもさらに急速に小さくなる数です。よってライプニッツの級数よりも少ない繰り返し
4
回数で高い精度の値を求めることができます。
5
また、以下は高野9の公式です。マチンの公式よりもさらに高速な計算が可能です。
12 arctan(1
49) + 32 arctan( 1
57)−5 arctan( 1
239) + 12 arctan( 1
110443) = π 4 その他にも類似の公式が知られていますが、詳細は省略します。
6
6.2.4 ガウス＝ルジャンドルの反復計算
7
4番目に述べる円周率の計算方法は、特に最近の円周率の計算に用いられているものです。
8
計算方法は、四つの数列an、bn、tn、pnについて以下の初期値：
9
a1 = 1, b1 = 1
√2, t1 = 1
4, p1 = 1 と漸化式（n≥1）：
an+1 = an+bn 2 , bn+1 = √
an·bn,
tn+1 = tn−pn·(an−an+1)2, pn+1 = 2·pn
9高野喜久雄。日本の詩人、数学者。
(a)プログラム
(b)実行途中結果
図6.16: ガウス＝ルジャンドルの反復計算法による円周率の計算プログラム
で反復計算します。そして円周率の値を以下の式で求めます。
1
π≈ (an+bn)2 4·tn
この計算方法の理論的背景は複雑であるため、ここでは述べません。余力のある人は自ら調べてく
2
ださい。この計算式は計算手順が単純であるにも関わらず、計算効率がきわめて良い（真値への収
3
束が速い）ことが知られています。
4
図6.16にScratchのプログラムを載せます。以下はn= 1,2,3の場合の計算結果です。
5
3.140579
6
3.141593
7
3.141593
8
急速に真値へ近づくことが分かります。
9
このテキストのここまでの説明では、1匹の猫のスプライトのみを取り扱い、その猫の動作を
2
様々にプログラミングしてきました。しかし、Scratchでは同時に複数のスプライトをステージ上
3
で配置し、それぞれをそれぞれのプログラムで動かすことができます。また「ステージ」にもプロ
4
グラムを与えることができます。この説ではそのような二つの例を紹介します。前節までのプログ
5
ラミングとは様子が異なることに気づくでしょうか。
6
複数のプログラムが同時に動くことを一般に
並行実行
（concurrent execution）7
と呼びます10。しかし複数のプログラムがそれぞれ独立に動くだけでは並行実行は意味を持ちま
8
せん。プログラム間でデータのやり取りが行われることで、単独のプログラムとは異なる様相の複
9
雑なコンピュータ処理が可能になります。
10
データのやり取りを行う方法にはいくつかの方法があります。最も代表的な方法は
11
メッセージパッシング
（message passing）および12
共有メモリ
（shared memory）です。この節ではそれぞれの事例を紹介します。13
6.3.1 猫と犬の掛け合い
14
次のようなプログラムを考えましょう。
15
猫と犬の掛け合い
緑の旗をクリックすると、猫が「ニャー」と鳴き、次に犬が「ワン」と鳴き、次に猫が「ニャー」
と鳴き、...という動作を1秒間隔で繰り返すプログラムを、猫と犬をそれぞれスプライトと して作りなさい。
16
猫と犬のそれぞれの動作は簡単です。問題は、鳴くタイミングをうまくコントロールすることで
17
す。ここではそれを、「イベント」のブロックパレット（図5.6の右列）の下の方の「[ ]を受け取っ
18
たとき」ブロックと「[ ]を送る」ブロックを用いて実現します。[ ]の中にはメッセージの種類を指
19
定します。Scratchのオンラインエディタではあらかじめ「メッセージ1」が登録されていますが、
20
メッセージの種類はユーザが自由に登録することができます。ブロックの[ ]の中の「メッセージ
21
１▼」の▼をマウスでクリックして新規登録してください。
22
さて、猫と犬のスプライトに作成したプログラムの例を図6.17に示します。猫のプログラムでは、
23
1. 「犬が鳴いた」というメッセージを受信したときに猫のプログラムが起動し、
24
2. 猫のコスチュームを変更（スプライトの画像を少し変更）し、
25
3. 「ニャー」と１秒間言い、
26
10類似の術語に「並列実行」（parallel execution）があります。これら術語は区別されて使用されますが、その違いに ついてはこの授業の範囲を超えるので省略します。