python怎么进行计算_python – sympy如何计算pi？

mpmath使用Chudnovsky公式(

https://en.wikipedia.org/wiki/Chudnovsky_algorithm)计算pi.

Pi由无穷大系列近似,其项减少如(1/151931373056000)^ n,因此每个项大约增加14.18位.这使得容易选择多个项N,从而实现期望的精度.

实际计算是用整数运算完成的：即,对于预比特的精度,计算pi * 2 ^(prec)的近似值.

# Constants in Chudnovsky's series

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's series (see comments in

libelefun.py for details).

"""

# The Chudnovsky series gives 14.18 digits per term

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)

使Chudnovsky公式快速运行的关键技巧称为二进制分裂.无穷级数中的项满足小系数的递归关系(递归方程可以在bs_chudnovsky函数中的b-a == 1情况下看到).不是按顺序计算术语,而是将总和重复分成两半;递归地评估这两半,并将结果组合.最后,一个具有两个大整数p和q,使得该系列的前N个项的和恰好等于p / q.请注意,二进制分割过程中没有舍入误差,这是算法的一个非常好的特性;唯一的舍入发生在计算平方根并在最后进行除法时.

计算pi到高精度的大多数快速程序使用非常类似的策略,尽管有一些复杂的技巧可以进一步加快进程.

(注意：我是代码的作者.)