The Chudnovsky algorithm is a fast method for calculating the digits of π. It was published by the Chudnovsky brothers in 1989,^[1] and was used in the world record calculations of 2.7 trillion digits of π in December 2009,^[2] 5 trillion digits of π in August 2010,^[3] 10 trillion digits of π in October 2011,^[4][5] 12.1 trillion digits in December 2013^[6] and 22.4 trillion digits of π in November 2016. ^[7]

The algorithm is based on the negated Heegner number d = −163, the j-function j(1+√−163/2) = −7005640320000000000♠640320³, and on the following rapidly convergent generalized hypergeometric series:^[2]

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:

C = 426880√10005
Mk+1 = Mk * (12k+2) * (12k+6) * (12k+10) / (k+1)^3 and M0 = 1  [Mk = (6k)! / ((3k)! * (k!)^3)]
Lk+1 = Lk + 545,140,134 and L0 = 13,591,409  [Lk = 13591409 + 545140134*k]
Xk+1 = Xk * -262,537,412,640,768,000 and X0 = 1  [Xk = (-640320)^3k = (-262537412640768000)^k]

Mk can be optimized further:
Kk+1 = Kk + 12 and K0 = 6
Mk+1 = Mk * (Kk^3 - 16Kk) / (k+1)^3 and M0 = 1

Note that

and

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

Example: Python Implementation

π 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=70, prec=1008, disp=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 = (K**3 - (K<<4)) * M / 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=%d iterations, gc().prec=%d, disp=%d digits) =\n%s" % (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

• Numerical approximations of π
• Ramanujan–Sato series

References