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[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] 50 trillion digits on January 29, 2020,[7] and 62.8 trillion digits on August 14, 2021.[8]

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

${\displaystyle {\frac {1}{\pi }}=12\sum _{q=0}^{\infty }{\frac {(-1)^{q}(6q)!(545140134q+13591409)}{(3q)!(q!)^{3}\left(640320\right)^{3q+{\frac {3}{2}}}}}}$

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

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

${\displaystyle {\frac {(640320)^{\frac {3}{2}}}{12\pi }}={\frac {426880{\sqrt {10005}}}{\pi }}=\sum _{q=0}^{\infty }{\frac {(6q)!(545140134q+13591409)}{(3q)!(q!)^{3}\left(-262537412640768000\right)^{q}}}}$

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

${\displaystyle \pi =C\left(\sum _{q=0}^{\infty }{\frac {M_{q}\cdot L_{q}}{X_{q}}}\right)^{-1}}$, where:
${\displaystyle C=426880{\sqrt {10005}}}$,
${\displaystyle M_{q}={\frac {(6q)!}{(3q)!(q!)^{3}}}}$,
${\displaystyle L_{q}=545140134q+13591409}$,
${\displaystyle X_{q}=(-262537412640768000)^{q}}$.

The terms Mq, Lq, and Xq satisfy the following recurrences and can be computed as such:

{\displaystyle {\begin{alignedat}{4}L_{q+1}&=L_{q}+545140134\,\,&&{\textrm {where}}\,\,L_{0}&&=13591409\\[4pt]X_{q+1}&=X_{q}\cdot (-262537412640768000)&&{\textrm {where}}\,\,X_{0}&&=1\\[4pt]M_{q+1}&=M_{q}\cdot \left({\frac {(12q-2)(12q-6)(12q-10)}{q^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[4pt]\end{alignedat}}}

The computation of Mq can be further optimized by introducing an additional term Kq as follows:

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

Note that

${\displaystyle e^{\pi {\sqrt {163}}}\approx 640320^{3}+743.99999999999925\dots }$ and
${\displaystyle 640320^{3}=262537412640768000}$
${\displaystyle 545140134=163\cdot 127\cdot 19\cdot 11\cdot 7\cdot 3^{2}\cdot 2}$
${\displaystyle 13591409=13\cdot 1045493}$

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 ${\displaystyle O\left(n\log(n)^{3}\right)}$.[10]