Python：大きな数を分けるとmpmathの精度が失われる

Question

Chudnovskyアルゴリズムというというπアルゴリズムが出回った。 Pythonの実装はWikipediaに示されています。decimalパッケージにはPythonが付属しています。しかし最近私がGauss Legendreアルゴリズムをテストしたところ、mpmathパッケージは、高精度計算を扱うときにはdecimalよりも効率的に実行されることがわかったので、アルゴリズムがmpmathと動作することを願っています。ここに私のcodeです：

#!/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)) 

私は小数部は上げません反復のように、私はVAR K, L, Xにmpmathを使用していないので、Pythonの自体が大きな整数を扱う可能性が信じています。私はこの問題がS += (M * L)/Xにあると思うので、Xはかなり大きい数字です。 このような多数の問題に対処するには、私は多くのことを混乱させました。

Answer 1

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.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749上で動作します56735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809533 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309比較のために218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809523

、ここでmp.piの値です：

mp.dps = 1008 
print(mp.pi) 

3。1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989 3809525