sympy如何计算pi？

Question

sympy计算圆周率的数值背景是什么？

我知道 SymPy 在后台使用 mpmath，这使得使用任意精度算法执行计算成为可能。这样，一些特殊常量，如 e、pi、oo，将被视为符号，并且可以以任意精度进行计算。

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

Answer 1

mpmath 使用 Chudnovsky 公式 (https://en.wikipedia.org/wiki/Chudnovsky_algorithm) 计算圆周率。

Pi 近似为一个无限级数，其项递减为 (1/151931373056000)^n，因此每项相加大约 14.18 位。这使得 select 多个术语 N 变得容易，从而达到所需的精度。

实际计算是用整数运算完成的：也就是说，对于 prec 位的精度，pi * 2^(prec) 被计算出来。

这是从 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 类型。 @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 的快速程序都使用非常相似的策略，尽管有一些复杂的技巧可以进一步加快该过程。

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

sympy如何计算pi？

How does sympy calculate pi?

python

pi

sympy