The Chudnovsky algorithm is a fast method for calculating the digits of π, based on Ramanujan’s π formulae. It was published by the Chudnovsky brothers in 1988,^[1] and was used in the world record calculations of 2.7 trillion digits of π in December 2009,^[2] 10 trillion digits in October 2011,^[3][4] 22.4 trillion digits in November 2016,^[5] 31.4 trillion digits in September 2018–January 2019,^[6] and 50 trillion digits on January 29, 2020.^[7]

The algorithm is based on the negated Heegner number , the j-function , and on the following rapidly convergent generalized hypergeometric series:^[2]

A detailed proof of this formula can be found here:^[8]

For a high performance iterative implementation, this can be simplified to

There are 3 big integer terms (the multinomial term M_k, the linear term L_k, and the exponential term X_k) that make up the series and π equals the constant C divided by the sum of the series, as below:

, where:

The terms M_k, L_k, and X_k satisfy the following recurrences and can be computed as such:

The computation of M_k can be further optimized by introducing an additional term K_k as follows:

Note that

and

This identity is similar to some of Ramanujan's formulas involving π,^[2] and is an example of a Ramanujan–Sato series.

The time complexity of the algorithm is .^[9]

Example: Python implementation

^{[original research?]}

π can be computed to any precision using the above algorithm in any environment which supports arbitrary-precision arithmetic. As an example, here is a Python implementation:

from decimal import Decimal as Dec, getcontext as gc

def PI(maxK: int = 70, prec: int = 1008, disp: int = 1007):  # Parameter defaults chosen to gain 1000+ digits within a few seconds
    gc().prec = prec
    K, M, L, X, S = 6, 1, 13591409, 1, 13591409
    for k in range(1, maxK + 1):
        M *= Dec(K**3 - 16*K) / k**3 
        L += 545140134
        X *= -262537412640768000
        S += Dec(M * L) / X
        K += 12
    pi = 426880 * Dec(10005).sqrt() / S
    pi = Dec(str(pi)[:disp])  # Drop few digits of precision for accuracy
    print("PI(maxK={} iterations, gc().prec={}, disp={} digits) =\n{}".format(maxK, prec, disp, pi))
    return pi

Pi = PI()
print("\nFor greater precision and more digits (takes a few extra seconds) - Try")
print("Pi = PI(317, 4501, 4500)") 
print("Pi = PI(353, 5022, 5020)")

See also

• Borwein's algorithm
• Numerical approximations of π
• Ramanujan–Sato series

References