円周率を求める

Python

この項目は次のページを参考にしました。円周率以外にもPythonコードでさまざまな計算を紹介していますので、教育関係者に限らず是非読んでみてください。

高知大学 共通教育 情報処理 II 第6回の教材 (3)

以下は私が少し変更を加えたソースです。

## arctan(1/n)をTaylor展開で求める
def arctan(p, n):
    x = p/n
    nn = n*n
    c = 1
    s = x
    k = 1
    while x > 0:
        x /= nn
        k += 2
        c = -c
        s += c*(x/k)
    return s

## 円周率を求める
## 整数部分「3」も含む仮数表現の長整数で返します
def pi_mantissa(digit, redund=10):
    p = 10**digit
    q = 10**redund
    p *= q
    pi = 4*(12*arctan(p,18)+8*arctan(p,57)-5*arctan(p,239))
    pi /= q
    return pi

print pi_mantissa(1000)

3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450
2841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881
7488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272
4891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609
1736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185
9502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053
21712268066130019278766111959092164201989

引数が1000だと1001桁の長整数を返します。最後の「9」が小数点以下1000桁目となります。

decimalを使って円周率を求める

実はPythonライブラリリファレンスのdecimalモジュールの解説には、これとは別の円周率計算方法が載っています。

decimal ― Decimal fixed point and floating point arithmetic

以下は私が少し変更を加えたソースです。

from __future__ import with_statement
import decimal

def pi(prec):
    with decimal.localcontext() as ctx:
        ctx.prec = prec+2
        lasts = 0
        t = decimal.Decimal(3)
        s = 3
        n = 1
        na = 0
        d = 0
        da = 24
        while s != lasts:
            lasts = s
            n += na
            na += 8
            d += da
            da += 32
            t = (t * n) / d
            s += t
    return s

print pi(1000)

これは連分数表現のひとつを使って計算する方法みたいです。上のarctanと長整数演算を使ったアルゴリズムが爆速なのに比べると、収束するまで時間がかかります。

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745
0284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588
1748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527
2489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960
9173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318
5950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805
3217122680661300192787661119590921642019895

こちらは引数に1000を指定すると、decimalの精度を1002桁として計算します。decimalのオリジナルの方では計算した後、最後の1桁を切り捨ててちょうど小数点以下1000桁になるようにしていますが、ここでは省略しています。最後の桁の「5」は間違っているだろうと思って無視してください（桁数を増やして計算してみれば分かりますが、正しくは1001桁目は「3」です）。