3月14日は円周率の日。pythonで円周率を計算した話

概要

3月14日になると、コンビニやデパートがお菓子のセールスでにぎわっており、春も近づいてきたんだと感慨にふける人たちも多いと思います。
世間の皆さんもやはり円周率の日に注目を集めているわけで、経済効果もかなりあるのではないかと思ったりするわけです。
ですので、今日は円周率を計算するプログラムを書いてみます。

素朴なネタです。

ガウス・ルジャンドルのアルゴリズム

wikipediaのガウス＝ルジャンドルのアルゴリズムによると、割と簡単に実装できます。

こんなの考えるなんて、すごい天才なんだと思いました。

python実装

といっても、自分は円周率の計算方法はよくわからないので、いろいろな人のソースを見たりして出来上がりました。
特徴としては、多数の桁数を表すため、Decimalを利用したことくらい。

# -*- coding:utf-8 -*-
import math
import time
from decimal import *

import numpy as np


#文字数1002
PI_1000 = "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"


def get_pi(prec=1000, verbose=False):
    '''
    小数点以下の桁数がprecの円周率を文字列で返す
    '''
    prec = prec+1+1 #精度 "3."の部分があるので、1つ精度を増やして、さらにもう1けたないと丸目の関係でずれることがあるので、さらに1追加
    N = 2*math.ceil(math.log2(prec)) #繰り返し回数 wikipediaによると、log2(prec)程度でいいらしいので、念の為、その２倍程度にする

    getcontext().prec = prec 

    a = Decimal(1.0)
    b = Decimal(1.0) / Decimal(2.0).sqrt()
    t = Decimal(0.25)
    p = Decimal(1.0)


    start_time = time.clock() 

    # N回の試行を開始
    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)
    #実行時間を計算
    elapsed = time.clock()  - start_time

    #　実行結果を表示する
    if verbose:
        print("精度：",prec)
        print("円周率:", pi)
        print("実行時間:%f" % (elapsed))

    return str(pi)[0:prec]

下のように実行して1000桁くらいは正しいことを確かめました。（2億桁くらい計算している人もいるそうですけど、これくらいでいいでしょ。）

pi = get_pi(prec=1000, verbose=True)
print(pi)
#print(PI_1000)
#print(pi == PI_1000)

精度： 1002
円周率: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410
270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925
409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244
065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079
227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313
783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019890
実行時間:0.014581

統計情報

統計学入門に、円周率に出てくる数字の出現率を調べる問題があったので、解いてみました。

グラフを見たあと、カイ2乗検定も実施して、数字が一定に出ているといえるか調べてみます。

123桁目から100個の数字をとってきたとき

123桁目というのは適当です。

start = 123
pi_123 = pi[2+start:2+start+100]

freq = np.zeros(10, dtype=np.int)
for a in pi_123:
    freq[int(a)] += 1
print("freq(100)=",freq)

#chisquareを使って求める
print(chisquare(freq, 10*np.ones(10,dtype=np.int8)))

plt.xlim(-0.5, 9.5)
plt.bar(np.array([0,1,2,3,4,5,6,7,8,9]), freq, align='center')

グラフを見ると直感的には、気持ちが悪いけど、p値は0.37くらいなので、だいたい出現率は一定といっても矛盾はない。

freq(100)= [ 9 12 11  7 15 13  7  4 12 10]
Power_divergenceResult(statistic=9.8000000000000007, pvalue=0.3669177991127523)

1000桁の時

freq = np.zeros(10, dtype=np.int)
for a in pi:
    if a != ".":
        freq[int(a)] += 1
print("freq(1001)=",freq)
print(chisquare(freq, 100.1*np.ones(10,dtype=np.int8)))

plt.xlim(-0.5, 9.5)
plt.bar(np.array([0,1,2,3,4,5,6,7,8,9]), freq, align='center')

1000個くらいデータが集まると、出現率は数字に依らず一定と考えても矛盾はない。p値も0.85でいい感じ。

freq(1001)= [ 93 116 103 103  93  97  94  95 101 106]
Power_divergenceResult(statistic=4.7842157842157853, pvalue=0.85269838594638103)

他気になること

本当は、数字の並びの統計情報とかも知りたかったりしたけど、それは次ということで。
超越数とかいつか勉強しないとな。