pythonで円周率求めてみた！
はじめまして、にわこまです。今回は円周率を求めてみました。使うプログラミング言語はpythonです。
文章表現や内容に不備があるかもしれませんが、最後まで読んでいただけるとありがたいです。
スポンサードサーチ
円周率とは？
そもそも円周率とは何かわかりますか？
小学生くらいのときに教わるそうなのですが、私は忘れていました（笑）。
円周率は「円の円周の直径に対する比率」です。求めた方としては、円の円周を直径で割るだけです。
しかし、円の円周はどのように求めるのでしょうか？円の円周を正確に求めるのは難しいので、円周率を求めるアルゴリズムを使って円周率を求めます。
私はガウス＝ルジャンドルのアルゴリズムを使って求めてみました。この方法をpythonを使って解説していきたいと思います。
ちなみに、円周率の桁数はGoogleが2019年3月14日に発表した、小数点以下約31兆4000億桁が最長記録です。ギネスにも認定されています。
（https://www.itmedia.co.jp/news/articles/1903/14/news148.html）
ノートPCでも100桁くらいなら簡単に計算できます。
ガウス＝ルジャンドル
ガウス＝ルジャンドル
ガウス＝ルジャンドルのアルゴリズムとはなにか？
ガウス＝ルジャンドルのアルゴリズムは、初期値を反復式に代入して計算し、算出された値を円周率を求める式に代入して円周率を求める方法です。
くわしく見ていきましょう。
初期値は以下の通りです。
a0 = 1
b0 = 1/√2
t0 = 1/4
p0 = 1
反復式は以下の通りです。
an+1 = (an + bn)/2
bn+1 = √(anbn)
tn+1 = tn – pn(an – an+1)2
pn+1 = 2pn
円周率（π）を求める式は以下の通りです。
π = (a + b)2/4t
求めたい桁数になるように反復式を繰り返し、πを求める式に代入します。反復式の反復回数は求めたい小数点以下の桁数（ｎ）のとき、log2nです。
では実際にコードを見ていきましょう。
import time
import math
import numpy as np
from decimal import *
PI_1000 = "3.1415926535897932384626433832795028841971693993751058209749445923
07816406286208998628034825342117067982148086513282306647093844609550582231725
35940812848111745028410270193852110555964462294895493038196442881097566593344
61284756482337867831652712019091456485669234603486104543266482133936072602491
41273724587006606315588174881520920962829254091715364367892590360011330530548
82046652138414695194151160943305727036575959195309218611738193261179310511854
80744623799627495673518857527248912279381830119491298336733624406566430860213
94946395224737190702179860943702770539217176293176752384674818467669405132000
56812714526356082778577134275778960917363717872146844090122495343014654958537
10507922796892589235420199561121290219608640344181598136297747713099605187072
11349999998372978049951059731732816096318595024459455346908302642522308253344
68503526193118817101000313783875288658753320838142061717766914730359825349042
87554687311595628638823537875937519577818577805321712268066130019278766111959
092164201989"
def gauss(prec=1000):
# 反復回数
N = 2*math.ceil(math.log2(prec))
# 初期値の設定
a = Decimal(1)
b = Decimal(1) / Decimal(2).sqrt()
t = Decimal(0.25)
p = Decimal(1)
# 処理時間 計測開始
start_time = time.clock()
# 反復式
for _ in range(1, N):
a_next = (a + b) / Decimal(2)
b_next = (a*b).sqrt()
t_next = t - p * (a - a_next)**2
p_next = Decimal(2)*p
a = a_next
b = b_next
t = t_next
p = p_next
pi = (a + b)**2 / (Decimal(4)*t)
# 処理時間 計測終了
processing_time = time.clock() - start_time
print("小数点以下： ", prec)
print("繰り返し回数", N)
print("実行時間: %f" % processing_time)
return str(pi)[0:prec]
# 精度（桁数）
precision = 100
pi = gauss(prec=precision)
count = -2
for i in range(0, 1000):
if(count < precision):
if(PI_1000[i] == pi[i]):
count = count + 1
else:
break
else:
break
print("小数点以下", count, "桁まで正確")
上記のプログラムで実行すると以下のような結果になります。今回は小数点以下100桁を求めます。また上記のコードではPI_1000を改行していますが実際に使うときは改行を外してお使いください
小数点以下： 100
繰り返し回数： 14
実行時間：0.000088
小数点以下 26 桁まで正確
上記のプログラムでは小数点以下26桁までが正確に計算出来ているみたいです。
式や正誤確認のプログラムも間違っていないのに、なぜか小数点以下100桁まで正確に計算できません。
色々調べた結果、pythonでは整数を代入すると極微小な誤差があることがわかりました。普通の計算では気にならない誤差ですが、小数点以下100桁も計算するとかなりの誤差になります。
例
i = 1
i = 1.00000 … 000001111
iに1を代入しても、実際には小数点以下数10桁以降に誤差が生じます。
ではどのように小数点以下を定義するのか？
getcontext().precを使います。
（https://docs.python.org/ja/3/library/decimal.html）
getcontext().prec = 正確に計算したい桁数
で計算することが出来ます。
スポンサードサーチ
手直ししたコード
手直ししたコードは以下のとおりです。
import time
import math
import numpy as np
from decimal import *
PI_1000 = "3.14159265358979323846264338327950288419716939937510582097494459230
781640628620899862803482534211706798214808651328230664709384460955058223172535
940812848111745028410270193852110555964462294895493038196442881097566593344612
847564823378678316527120190914564856692346034861045432664821339360726024914127
372458700660631558817488152092096282925409171536436789259036001133053054882046
652138414695194151160943305727036575959195309218611738193261179310511854807446
237996274956735188575272489122793818301194912983367336244065664308602139494639
522473719070217986094370277053921717629317675238467481846766940513200056812714
526356082778577134275778960917363717872146844090122495343014654958537105079227
968925892354201995611212902196086403441815981362977477130996051870721134999999
837297804995105973173281609631859502445945534690830264252230825334468503526193
118817101000313783875288658753320838142061717766914730359825349042875546873115
95628638823537875937519577818577805321712268066130019278766111959092164201989"
def gauss(prec=1000):
# 反復回数
N = 2*math.ceil(math.log2(prec))
# 計算精度（桁数）
prec = prec + 2
getcontext().prec = prec
# 初期値の設定
a = Decimal(1)
b = Decimal(1) / Decimal(2).sqrt()
t = Decimal(0.25)
p = Decimal(1)
# 処理時間 計測開始
start_time = time.clock()
# 反復式
for _ in range(0, N):
a_next = (a + b) / Decimal(2)
b_next = (a*b).sqrt()
t_next = t - p * (a - a_next)**2
p_next = Decimal(2)*p
a = a_next
b = b_next
t = t_next
p = p_next
pi = (a + b)**2 / (Decimal(4)*t)
# 処理速度を求めるためのもの
# 終了時間から開始時間を引く
processing_time = time.clock() - start_time
print("小数点以下：", prec)
print("繰り返し回数： ", N)
print("処理時間: %f" % processing_time)
return str(pi)[0:prec]
# 精度（桁数）
precision = 100
pi = gauss(prec=precision)
count = -2
for i in range(0, 1000):
if(count < precision):
if(PI_1000[i] == pi[i]):
count = count + 1
else:
break
else:
break
print("小数点以下", count, "桁まで正確")
実行結果は以下の通りです。また、上記のコードではPI_1000を改行していますが実際に使う際は改行を外してお使いください。
小数点以下： 102
繰り返し回数： 14
処理時間: 0.000297
小数点以下 100 桁まで正確
処理時間はすこし遅くなっていますが、小数点以下100桁まで正確に計算できています。四捨五入などで値が変わらないように2桁ほど多く設定しています。
pythonでも円周率を求めることが出来る
ノートPCで使用言語pythonでも小数点以下100桁くらいなら余裕で計算できました。
10万桁くらいになるとかなり時間がかかります。pythonはプログラム言語のなかでは計算処理は早い方ですがそれでもかなりの処理時間がかかってしまいます。
処理時間を短くするためには他の方法が必要です。 pythonで円周率や小数点以下の計算をするときはgetcontext()を使う必要があると分かりました。
スポンサードサーチ
