Chudnovsky algorithm

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

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:

1 π = 12 ∑ k = 0 ∞ ( − 1 ) k ( 6 k ) ! ( 545140134 k + 13591409 ) ( 3 k ) ! ( k ! ) 3 ( 640320 ) 3 k + 3 / 2

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

( 640320 ) 3 / 2 12 π = 426880 10005 π = ∑ k = 0 ∞ ( 6 k ) ! ( 545140134 k + 13591409 ) ( 3 k ) ! ( k ! ) 3 ( − 262537412640768000 ) k

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:

π = C ( ∑ k = 0 ∞ M k ⋅ L k X k ) − 1 , where:C = 426880√10005
M_k+1 = M_k * (12k+2) * (12k+6) * (12k+10) / (k+1)^3 and M₀ = 1 [M_k = (6k)! / ((3k)! * (k!)^3)]
L_k+1 = L_k + 545,140,134 and L₀ = 13,591,409 [L_k = 13591409 + 545140134*k]
X_k+1 = X_k * -262,537,412,640,768,000 and X₀ = 1 [X_k = (-640320)^3k = (-262537412640768000)^k]M_k can be optimized further:
K_k+1 = K_k + 12 and K₀ = 6
M_k+1 = M_k * (K_k^3 - 16K_k) / (k+1)^3 and M₀ = 1

Note that

e π 163 ≈ 640320 3 + 743.99999999999925 … and 640320 3 = 262537412640768000 545140134 = 163 ⋅ 127 ⋅ 19 ⋅ 11 ⋅ 7 ⋅ 3 2 ⋅ 2 13591409 = 13 ⋅ 1045493

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