sympy如何计算pi？

Question

sympy计算pi的数值背景是什么？

我知道SymPy在后台使用mpmath，这使得使用任意精度算法来执行计算成为可能。这样，一些特殊的常量，如e，pi，oo就被当作符号来处理，并且可以用任意的精度进行评估。

但是，Sympy如何确定任意小数位数？ Sympy如何在数字上做到这一点？

Answer 1

mpmath使用Chudnovsky公式计算pi（https://en.wikipedia.org/wiki/Chudnovsky_algorithm）。

Pi近似为一个无限序列，其项减少为（1/151931373056000）^ n，所以每个项增加大约14.18位数。这使得易于选择多个术语N，从而达到期望的精度。

实际计算与整数运算来完成：即，为了PREC位精度，* 2 ^（PREC）被计算pi的近似值。

下面是代码，从mpmath/libmp/libelefun.py萃取（https://github.com/fredrik-johansson/mpmath/blob/master/mpmath/libmp/libelefun.py）：

# 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的内置long类型。装饰器缓存计算值（在数值计算中通常需要重复使用pi，所以存储它是有意义的）。

你可以看到它做一个基数转换如下计算圆周率：

>>> pi_fixed(53) * 10**16 // 2**53 
mpz(31415926535897932) 

的关键技巧，使Chudnovsky式快速被称为二元分裂。无穷级数中的项满足小系数的递推关系（递推方程可以在bs_chudnovsky函数的b-a == 1情况下看到）。不是按顺序计算条件，而是将总和分成两半;递归评估两半，并将结果合并。最后，一个具有两个大整数p和q，使得该系列的第一Ñ项之和恰好等于p/q。请注意，二进制拆分过程中不存在舍入误差，这是该算法的一个非常好的功能;当计算平方根并在最后进行分割时，唯一的舍入发生。

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

（注意：我是代码的作者。）