意外と深い円周率を求めるプログラミング～数学を学ぼう～

突然ですが、初めて円周率を自分で計算したときのことを覚えていますか？私は鮮明に覚えています。そう、あれば小学校５年の算数の授業でした。

先生に急にダンボールと糸を渡され、「円周率３．１４を導け」と言われたのです。当時はうまく理解できていなかった気がしますがなんとか手伝ってもらいながら無事３．１４に近い値を算出したことを覚えています。そのときの手順は次の通り。ものすごくアナログです（笑）

1. コンパスでダンボールに円を描く
2. 円を切り抜く
3. 円の周りに沿って糸を一周巻いて切る
4. 糸の長さを測る
5. 糸の長さをダンボールの円の直径で割る

ちなみに手順３が一番難しいです。適当に切ると長すぎたり短くなったりするんです。コツはちょっと余分に巻いて重なっている部分をカッターで一気に切るんです！こうするとぴったり切れるんです。
もし暇すぎで死んでしまいそうだったらやってみてください。円周は直径に円周率をかけることで求められるので、円周率をもとめるには円周（糸の長さ）をダンボールの円の直径で割れば円周率っぽいものが得られるというものです。

\(\large r \): 半径
\(\large \pi\)：円周率
\(\large l \)：円周

上のように設定したとき、円の円周は次の式で求めることが可能です

\(\large 2\pi r = l\)

これを\(\large\pi\)で整理します

\(\large\pi = \frac{l}{2r}\)

大数の法則で近似値を求める方法

私もあれから大きくなりました。いまではもうおっさんです。おっさんがダンボールと糸でちまちまやっててもかっこよくないのでプログラムで円周率を求めてみましょう。

その前に問題です。
1辺が1の正方形があります。

この四角の半分に色を塗ります。

この四角形の上にランダムに点を打つと青い部分に入る確率は何%でしょうか？正解は、５０% です。
もちろん次のように５０%にならない場合もありますが点を増やすに従って結果は\(\large\frac{1}{2}\)に近づいていきます。いわゆる大数の法則というものです。

では、なぜ\(\large\frac{1}{2}\)に近づくか。答えは感覚的にも分かる通り、青い部分の面積が全体のちょうど５０%だからです。このことからわかるように面積の比率と穴が打たれる確率は一致するものと考えることが可能です。この考え方で円周率を求めていきます。

円周率を近似する方法

円周率の近似は次の図形を元に考えます。

この上にランダムな点をたくさん打ち込みます。

わかりやすいように円の内側に打たれた点と外側に打たれた点を
区別しておきましょう。

正方形に含まれる扇型の円の面積は半径1の円の\(\large\frac{1}{4}\)の大きさなので次のようになります。

\(\large \frac{1\times1\times\pi}{4} = \frac{\pi}{4}\)

つまり点をたくさん打ったとき、扇型の内側に点が入る確率は徐々に\(\large\frac{\pi}{4}\)になっていくということです。式にしてみます。全部の点を\(N\)、扇型に入った点を\(P\)とすると次のように書くことができます。

\(\large\frac{P}{N} = \frac{\pi}{4}\)

これを\(\large \pi\)で整理すると

\(\large\pi = \frac{4P}{N}\)

となります。
\(\large \pi\)を求めるには、扇型に入った点の数\((P\))を全ての点(\(N\))で割って４倍すれば求めることができます。

ということで、この内容をプログラムで実装していきます。ようやくプログラムの話ができますね。

プログラムでπを求めるのに必要なこと。

1辺が1の正方形に無数に点を打つにはどうすれば良いか？０から1までのランダムな数値を２つ用意した
それぞれを\(\large x\)座標、\(\large y\)座標を表すことで正方形の中に収まる点を打つことができます。

扇型の中に打たれた点を判断するにはどうすれば良いか？座標\(\large x\)と\(\large y\)がわかれば、三平方の定理で斜辺の長さが求められます。その長さが円の半径1以下であれば扇型の中に打たれた点と判断することができます。

上のことを踏まえて作成したのが次のプログラムです。

import random
import time

point = 0
N = 100
repeat = 5

start = time.time()
for count in range(repeat):
    N = N * 10
    point = 0
    for i in range(N):
        x = random.random()
        y = random.random()

        if x * x + y * y < 1.0:
            point += 1

        # 円周率の近似解
        pi = 4.0 * point / N
    process_time = time.time() - start
    print('N=' + str(N) + "  \t" + 'pi=' + str(pi) + "  \t" + 'time=' +
          str(process_time))

import random

import time

point = 0

N = 100

repeat = 5

start = time.time()

for count in range(repeat):

N = N * 10

point = 0

for i in range(N):

x = random.random()

y = random.random()

if x * x + y * y < 1.0:

point += 1

# 円周率の近似解

pi = 4.0 * point / N

process_time = time.time() - start

print('N=' + str(N) + " \t" + 'pi=' + str(pi) + " \t" + 'time=' +

str(process_time))

回数によって円周率の精度と処理時間を確認したいため、1000, 10000, 100000, 1000000, 10000000と点の数を変えて円周率を求める計算をしています。

実行結果はこちら

N=1000          pi=3.188        time=0.0013129711151123047
N=10000         pi=3.1612       time=0.015565156936645508
N=100000        pi=3.14284      time=0.14728093147277832
N=1000000       pi=3.138936     time=1.3955268859863281
N=10000000      pi=3.14151      time=13.65818190574646