4

我遇到了一种叫做Chudnovsky 算法的 π 算法。Wikipedia上显示了 Python 实现，它使用decimalPython 附带的包。但是最近在测试Gauss Legendre 算法时，我发现mpmath包的运行效率比decimal处理高精度计算时要高得多，所以我希望该算法可以与mpmath. 这是我的代码：

#!/usr/bin/env python

from mpmath import *
import pi_compare # A module aim to compare result with standard pi

mp.dps = 1000
def pi():
    K, M, L, X, S = 6, mpf('1'), 13591409, 1, mpf('13591409')
    for i in xrange(0,100):
        M = (K**3 - K*16) * M / K**3
        L += 545140134
        X *= -262537412640768000
        S += (M * L) / X
        K += 12
    return mpf('426880') * mpf('10005').sqrt() / S

P = pi()
print P
print pi_compare.compare(str(P))

我确实相信 Python 本身可以处理大整数，所以我不会mpmath在 var 上使用，K, L, X因为在迭代中不会出现小数部分。我认为这件事发生在S += (M * L) / X一个X相当大的数字上。处理这么大的数字让我很困惑，希望您的建议，谢谢。

4

1 回答 1

3

除了 DSM 提到的错字之外，该代码还有另一个问题。除法k**3需要是地板除法，而不是浮动除法。这是使用decimal和mpmath模块的修复版本。此代码适用于 Python 2 和 Python 3。

from decimal import Decimal as Dec, getcontext as gc
from mpmath import mp

def pi_dec(maxK=70, prec=1008):
    gc().prec = prec
    K, M, L, X, S = 6, 1, 13591409, 1, 13591409 
    for k in range(1, maxK+1):
        M = (K**3 - (K<<4)) * M // k**3 
        L += 545140134
        X *= -262537412640768000
        S += Dec(M * L) / X
        K += 12
    pi = 426880 * Dec(10005).sqrt() / S
    return pi

def pi_mp(maxK=70, prec=1008):
    mp.dps = prec
    K, M, L, X, S = 6, 1, 13591409, 1, 13591409 
    for k in range(1, maxK+1):
        M = (K**3 - (K<<4)) * M // k**3 
        L += 545140134
        X *= -262537412640768000
        S += mp.mpf(M * L) / X
        K += 12
    pi = 426880 * mp.sqrt(10005) / S
    return pi


Pi = pi_dec()
print(Pi)

Pi = pi_mp()
print(Pi)

输出

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809533 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001902787649042899268

出于比较目的，以下是 的值mp.pi：

mp.dps = 1008
print(mp.pi)

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300198018762592625928

于 2017-01-26T08:54:19.957 回答