概要
世の中には円周率を確率的に求める方法があります。有名所だとビュフォンの針でしょうか。これをコンピュータでシミュレーションするには乱数を使います。通常はコンピュータで発生させた乱数を使うのですが、今回はその乱数に円周率を使ってみました。
方法
今回はモンテカルロ法(という言い方で良いのかよくわからないのですが)を使いました。その詳細は以下のサイトに譲りますが、簡単に言うと正方形に内接する円を描いておいて、正方形内部にランダムに点を打っていき、円の内部に入った点の個数を数えることで円周率を求めていきます。
http://hp.vector.co.jp/authors/VA014765/pi/montecalro.html
今回は100x100ピクセルの架空の正方形を描き、その中にx,y座標それぞれ0-99までの乱数を座標とした点を置いていきます。ここで、その乱数に円周率を3141592...という数列にして、2桁ごとに区切った数を使いました。
円周率の入手
乱数として円周率を使うので、どこかから円周率を入手しなくてはなりません。今回は100万桁まではWebから入手し(下のURL)、それを超える場合は自前で計算したもの(結局25億桁計算しました)をつかうことにしました。
http://www.tstcl.jp/ja/randd/pi.php
↑100万桁の入手先
http://www.numberworld.org/y-cruncher/
↑円周率計算ソフト
言語
私が比較的得意という理由で(遅いですが)Pythonを使いました。
乱数としての性能の比較
本当は乱数に関する検定をして円周率の乱数としての性能を評価すべきですが、せっかく円周率を求めたので、その円周率の正確性をコンピュータで生成した乱数によるものと比較して評価しようと思います。
円周率を乱数とした場合の実験結果
使用した円周率の桁数は試行回数の4倍です(x,y座標それぞれに2桁ずつ使っているので)。
求まった円周率は求められた最高の有効数字で記述し、相対誤差は測定値と理論値共にその有効数字で計算した上で有効数字を2桁に統一しました。
ちなみに正確な円周率は3.14159265...です(そこそこ覚えてる)。
| 試行回数 | 求まった円周率 | 相対誤差 |
|---|---|---|
| 1.0*10^4 | 3.0548 | -2.8% |
| 1.0*10^6 | 3.068044 | -2.3% |
| 1.0*10^8 | 3.06708068 | -2.4% |
| 6.2*10^8 | 3.0671944387096772 | -2.4% |
コンピュータ生成の乱数による実験結果
データの処理方法は上と同じです。時間的な問題で表のそれぞれの実験は1回ずつしか行っていません。
| 試行回数 | 求まった円周率 | 相対誤差 |
|---|---|---|
| 1.0*10^4 | 3.1512 | 0.31% |
| 1.0*10^6 | 3.143292 | 0.054% |
| 1.0*10^8 | 3.14128316 | -0.0099% |
| 6.2*10^8 | 3.1415816256 | -0.00035% |
結論
円周率を乱数として使ってもそれなりに円周率は求まる。
ただ、コンピュータで生成した乱数のほうが比較的良い精度で求まりそうです。悲しいですね…なんででしょうか…しかも円周率使う方は特に円周率の読み込みでめっちゃ時間かかって辛い…
それと、今回の実験は単に、円周率の無駄遣い！6億桁使って3.14が求まらないとか、しょぼすぎですね。
今度もし気が向いたら乱数としてネイピア数を使ってみようかとも思います。Twitterでいただいたアイデアでネイピア数と円周率の組み合わせでまた別の乱数(?)を作るのも面白そうなので、また気が向いたらやってみます。
コード
GitHubで他のシミュレータと共に公開しています。
https://github.com/Nyanyan/Small-Projects/tree/master/CalculatingPI
コンピュータの乱数使用ver
全てのコードの基本的な構造です。
import sys
import random
num = 10000
r = 99/2
count=0
b = 200
j=1
for i in range(num):
x = random.uniform(0,2*r)
y = random.uniform(0,2*r)
if (x-r)**2 + (y-r)**2 <= r**2:
count+=1
if i % (num // b) == 0:
sys.stdout.write("\r")
for k in range(j):
sys.stdout.write("=")
for k in range(b - j):
sys.stdout.write(" ")
sys.stdout.write("|")
sys.stdout.flush()
j+=1
pi = 4 * count / num
print('')
print(count,num,pi)
Webから円周率を取ってくるver
あまり真面目に高速化していません。
import sys
import mechanize
num = 10000
r = 99/2
count=0
browser = mechanize.Browser()
response = browser.open('http://www.tstcl.jp/ja/randd/pi.php')
a = str(response.read())
start = a.find('3.')
pirandom = []
i=start
j=1
b = 50
while len(pirandom)<=num*4:
if ord(a[i]) >= 48 and ord(a[i]) <= 57:
pirandom.append(int(a[i]))
if len(pirandom) % (num*4 // b) == 0:
sys.stdout.write("\r")
for k in range(j):
sys.stdout.write("=")
for k in range(b - j):
sys.stdout.write(" ")
sys.stdout.write("|")
sys.stdout.flush()
j+=1
i+=1
print('')
j=1
for i in range(num):
x = pirandom[4*i]*10 + pirandom[4*i+1]
y = pirandom[4*i+2]*10 + pirandom[4*i+3]
if (x-r)**2 + (y-r)**2 <= r**2:
count+=1
if i % (num // b) == 0:
sys.stdout.write("\r")
for k in range(j):
sys.stdout.write("=")
for k in range(b - j):
sys.stdout.write(" ")
sys.stdout.write("|")
sys.stdout.flush()
j+=1
pi = 4 * count / num
print('')
print(count,num,pi)
pi.txtに円周率を格納し置いておくのを読み込むver
ちょっとだけ真面目に高速化しました。
import sys
num = 10000
r = 99/2
count=0
path = 'pi.txt'
with open(path) as f:
a = f.read()
pirandom = [0]*4
i=0
j=1
b = 200
c = 0
countnum=0
while countnum<=num:
if ord(a[i]) >= 48 and ord(a[i]) <= 57:
pirandom[c] = int(a[i])
c+=1
if c == 4:
x = pirandom[0]*10 + pirandom[1]
y = pirandom[2]*10 + pirandom[3]
if (x-r)**2 + (y-r)**2 <= r**2:
count+=1
countnum+=1
c=0
if countnum % (num // b) == 0:
sys.stdout.write("\r")
for k in range(j):
sys.stdout.write("=")
for k in range(b - j):
sys.stdout.write(" ")
sys.stdout.write("|")
sys.stdout.flush()
j+=1
i+=1
pi = 4 * count / num
print('')
print(count,num,pi)