mpmath computes pi using the Chudnovsky formula ( https://en.wikipedia.org/wiki/Chudnovsky_algorithm ).
Pi , (1/151931373056000) ^ n, 14,18 . N, .
: pi * 2 ^ (prec).
, mpmath/libmp/libelefun.py(https://github.com/fredrik-johansson/mpmath/blob/master/mpmath/libmp/libelefun.py):
CHUD_A = MPZ(13591409)
CHUD_B = MPZ(545140134)
CHUD_C = MPZ(640320)
CHUD_D = MPZ(12)
def bs_chudnovsky(a, b, level, verbose):
"""
Computes the sum from a to b of the series in the Chudnovsky
formula. Returns g, p, q where p/q is the sum as an exact
fraction and g is a temporary value used to save work
for recursive calls.
"""
if b-a == 1:
g = MPZ((6*b-5)*(2*b-1)*(6*b-1))
p = b**3 * CHUD_C**3 // 24
q = (-1)**b * g * (CHUD_A+CHUD_B*b)
else:
if verbose and level < 4:
print(" binary splitting", a, b)
mid = (a+b)//2
g1, p1, q1 = bs_chudnovsky(a, mid, level+1, verbose)
g2, p2, q2 = bs_chudnovsky(mid, b, level+1, verbose)
p = p1*p2
g = g1*g2
q = q1*p2 + q2*g1
return g, p, q
@constant_memo
def pi_fixed(prec, verbose=False, verbose_base=None):
"""
Compute floor(pi * 2**prec) as a big integer.
This is done using Chudnovsky series (see comments in
libelefun.py for details).
"""
N = int(prec/3.3219280948/14.181647462 + 2)
if verbose:
print("binary splitting with N =", N)
g, p, q = bs_chudnovsky(0, N, 0, verbose)
sqrtC = isqrt_fast(CHUD_C<<(2*prec))
v = p*CHUD_C*sqrtC//((q+CHUD_A*p)*CHUD_D)
return v
Python, , isqrt_fast(), . MPZ - : gmpy.fmpz, , Python - . @constant_memo (pi , ).
, pi, :
>>> pi_fixed(53) * 10**16 // 2**53
mpz(31415926535897932)
, , . ( b-a == 1 bs_chudnovsky). , , ; , . , p q , N p/q. , , ; .
, pi , , , .
(: .)