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 on January 29, 2020.

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

${\frac {1}{\pi }}=12\sum _{q=0}^{\infty }{\frac {(-1)^{q}(6q)!(545140134q+13591409)}{(3q)!(q!)^{3}\left(640320\right)^{3q+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 _{q=0}^{\infty }{\frac {(6q)!(545140134q+13591409)}{(3q)!(q!)^{3}\left(-262537412640768000\right)^{q}}}$

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

$\pi =C\left(\sum _{q=0}^{\infty }{\frac {M_{q}\cdot L_{q}}{X_{q}}}\right)^{-1}$, where:

$C=426880{\sqrt {10005}}$,

$M_{q}={\frac {(6q)!}{(3q)!(q!)^{3}}}$,

$L_{q}=545140134q+13591409$,

$X_{q}=(-262537412640768000)^{q}$.

The terms M_q, L_q, and X_q satisfy the following recurrences and can be computed as such:

${\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 M_q can be further optimized by introducing an additional term K_q as follows:

${\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

$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})$.

Programming example

An example implementation of the Chudnovsky algorithm, written in Python 3, is provided below.

#!/usr/bin/env python3
import decimal
from decimal import Decimal as D  # For internal use only.


def pi(precision: int = decimal.getcontext().prec,
       rounding: decimal.Context.rounding = decimal.ROUND_FLOOR
       ) -> decimal.Decimal:
    """Return pi as a decimal.Decimal to the specified precision.

    This function uses the Chudnovsky algorithm
    (https://wikipedia.org/wiki/Chudnovsky_algorithm) and implements
    acceleration by using the Mq (and Kq), Lq, and Xq recurrence
    patterns.

    Note that the first argument is the precision, NOT the number of
    decimal places. The relationship is:
        decimal places = precision + 1

    Positional/keyword arguments:
        precision -- the precision to which pi is calculated
                (default decimal.getcontext().prec)
        rounding -- how to round the final value of pi
                (default decimal.ROUND_FLOOR)
    """
    # Add 2 to the precision to ensure the accuracy of the answer.
    precision += 2
    decimal.getcontext().prec = precision

    # Use incremental calculations.
    # Set initial values for q, C, L, X, K, and M.
    q = D(0)
    C = D(426880)*D(10005).sqrt()
    L = D(13591409)
    X = D(1)
    K = D(-6)
    M = D(1)
    running_sum = D(0)
    pi = D('NaN')

    # Iteration 0.
    lastpi = pi
    running_sum += M*L/X
    pi = C * running_sum**D(-1)

    # Loop until pi and lastpi converge to the current precision.
    while pi - lastpi:
        lastpi = pi
        q += D(1)
        L += D(545140134)
        X *= D(-262537412640768000)
        K += D(12)
        M *= (K**D(3) - D(16)*K) / q**D(3)
        running_sum += M*L/X
        pi = C * running_sum**D(-1)

    # Reduce precision back to the specified precision.
    precision -= 2
    decimal.getcontext().prec = precision

    # quantize pi to sepcified precision with specified rounding.
    return pi.quantize(D('1e-'+str(precision-1)), rounding)


if __name__ == '__main__':
    print(pi())  # Calculate pi to current precision.
    print(pi(51))  # Calculate pi to 50 decimal places.
    # Expected output (https://oeis.org/A000796):
    #
    # >>> print(pi())
    # 3.141592653589793238462643383
    # >>> print(pi(51))
    # 3.14159265358979323846264338327950288419716939937510

See also

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