# Chudnovsky algorithm

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, and was used in the world record calculations of 2.7 trillion digits of π in December 2009, 10 trillion digits in October 2011,, 22.4 trillion digits in November 2016, 31.4 trillion digits in September 2018–January 2019, and 50 trillion digits in January 29, 2020.

The algorithm is based on the negated Heegner number $d=-163$ , the j-function $j\left({\frac {1+{\sqrt {-163}}}{2}}\right)=-640320^{3}$ , and on the following rapidly convergent generalized hypergeometric series:

${\frac {1}{\pi }}=12\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}\left(640320\right)^{3k+3/2}}}$ A detailed proof of this formula can be found here:

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

${\frac {(640320)^{3/2}}{12\pi }}={\frac {426880{\sqrt {10005}}}{\pi }}=\sum _{k=0}^{\infty }{\frac {(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}\left(-262537412640768000\right)^{k}}}$ There are 3 big integer terms (the multinomial term Mk, the linear term Lk, and the exponential term Xk) that make up the series and π equals the constant C divided by the sum of the series, as below:

$\pi =C\left(\sum _{k=0}^{\infty }{\frac {M_{k}\cdot L_{k}}{X_{k}}}\right)_{k}^{-1}$ , where:
$C=426880{\sqrt {10005}},\quad \quad M_{k}={\frac {(6k)!}{(3k)!(k!)^{3}}},\quad \quad L_{k}=545140134k+13591409,\quad \quad X_{k}=(-262537412640768000)^{k}$ The terms Mk, Lk, and Xk satisfy the following recurrences and can be computed as such:

{\begin{alignedat}{4}L_{k+1}&=L_{k}+545140134\,\,&&{\textrm {where}}\,\,L_{0}&&=13591409\\[4pt]X_{k+1}&=X_{k}\cdot (-262537412640768000)&&{\textrm {where}}\,\,X_{0}&&=1\\[4pt]M_{k+1}&=M_{k}\cdot \left({\frac {(12k+2)(12k+6)(12k+10)}{(k+1)^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[4pt]\end{alignedat}} The computation of Mk can be further optimized by introducing an additional term Kk as follows:

{\begin{alignedat}{4}K_{k+1}&=K_{k}+12\,\,&&{\textrm {where}}\,\,K_{0}&&=6\\[4pt]M_{k+1}&=M_{k}\cdot \left({\frac {K_{k}^{3}-16K_{k}}{(k+1)^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[12pt]\end{alignedat}} Note that

$e^{\pi {\sqrt {163}}}\approx 640320^{3}+743.99999999999925\dots$ and
$640320^{3}=262537412640768000$ $545140134=163\cdot 127\cdot 19\cdot 11\cdot 7\cdot 3^{2}\cdot 2$ $13591409=13\cdot 1045493$ This identity is similar to some of Ramanujan's formulas involving π, and is an example of a Ramanujan–Sato series.

The time complexity of the algorithm is $O(n\log(n)^{3})$ .

## 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: 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 = (K**3 - 16*K) * 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={} 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)")