How does sympy calculate pi?

What is the numerical background of sympy to calculate pi?

I know that SymPy uses mpmath in the background, which makes it possible to perform computations using arbitrary-precision arithmetic. That way, some special constants, like e, pi, oo, are treated as symbols and can be evaluated with arbitrary precision.

But how does Sympy determine the any number of decimal places? How does Sympy do it numerically?

Solution

mpmath computes pi using the Chudnovsky formula (https://en.wikipedia.org/wiki/Chudnovsky_algorithm).

Pi is approximated by an infinite series whose terms decrease like (1/151931373056000)^n, so each term adds roughly 14.18 digits. This makes it easy to select a number of terms N so that a desired accuracy is achieved.

The actual computation is done with integer arithmetic: that is, for a precision of prec bits, an approximation of pi * 2^(prec) is computed.

Here is the code, extracted from 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

This is just ordinary Python code, except that it depends on an extra function isqrt_fast() which computes the square root of a big integer. MPZ is the big integer type used: gmpy.fmpz if this is available, and Python's builtin long type otherwise. The @constant_memo decorator caches the computed value (pi is often needed repeatedly in a numerical calculation, so it makes sense to store it).

You can see that it computes pi by doing a radix conversion as follows:

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

The crucial trick to make the Chudnovsky formula fast is called binary splitting. The terms in the infinite series satisfy a recurrence relation with small coefficients (the recurrence equation can be seen in the b-a == 1 case in the bs_chudnovsky function). Instead of computing the terms sequentially, the sum is repeatedly split in two halves; the two halves are evaluated recursively, and the results are combined. In the end, one has two large integers p and q such that the sum of the first N terms of the series is exactly equal to p / q. Note that there is no rounding error in the binary splitting process, which is a very nice feature of the algorithm; the only roundings occur when computing the square root and doing the division at the very end.

Most fast programs that compute pi to high precision use a very similar strategy, though there are some complicated tricks that can speed up the process a bit further.

(Note: I'm the author of the code.)