Chudnovsky algorithm

Short description: Fast method for calculating the digits of π

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.

It 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, 50 trillion digits on January 29, 2020, 62.8 trillion digits on August 14, 2021, and 100 trillion digits on March 21, 2022.

Algorithm

The algorithm is based on the negated Heegner number [math]\displaystyle{ d = -163 }[/math], the j-function [math]\displaystyle{ j \left(\tfrac{1 + i\sqrt{163}}{2}\right) = -640320^3 }[/math], and on the following rapidly convergent generalized hypergeometric series:[math]\displaystyle{ \frac{1}{\pi} = 12 \sum_{k=0}^{\infty} {\frac{(-1)^k (6k)! (545140134k + 13591409)}{(3k)! (k!)^3(640320)^{3k + 3/2}}} }[/math]A detailed proof of this formula can be found here:

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

Optimizations

The optimization technique used for the world record computations is called binary splitting.

Binary splitting

A factor of [math]\displaystyle{ 1/{640320^{3/2}} }[/math] can be taken out of the sum and simplified to[math]\displaystyle{ \frac{1}{\pi} = \frac{1}{426880 \sqrt{10005}} \sum_{k=0}^{\infty}{\frac{(-1)^k (6k)! (545140134k + 13591409)}{(3k)! (k!)^3 (640320)^{3k}}} }[/math]

Let [math]\displaystyle{ f(n) = \frac{(-1)^n (6n)!}{(3n)! (n!)^3 (640320)^{3n}} }[/math], and substitute that into the sum.[math]\displaystyle{ \frac{1}{\pi} = \frac{1}{426880 \sqrt{10005}} \sum_{k=0}^{\infty}{f(k) \cdot (545140134k + 13591409)} }[/math]

[math]\displaystyle{ \frac{f(n)}{f(n-1)} }[/math] can be simplified to [math]\displaystyle{ \frac{-(6n-1)(2n-1)(6n-5)}{10939058860032000 n^3} }[/math], so[math]\displaystyle{ f(n) = f(n-1) \cdot \frac{-(6n-1)(2n-1)(6n-5)}{10939058860032000 n^3} }[/math]

[math]\displaystyle{ f(0) = 1 }[/math] from the original definition of [math]\displaystyle{ f }[/math], so[math]\displaystyle{ f(n) = \prod_{j=1}^{n}{\frac{-(6j-1)(2j-1)(6j-5)}{10939058860032000 j^3}} }[/math]

This definition of [math]\displaystyle{ f }[/math] isn't defined for [math]\displaystyle{ n=0 }[/math], so compute the first term of the sum and use the new definition of [math]\displaystyle{ f }[/math][math]\displaystyle{ \frac{1}{\pi} = \frac{1}{426880 \sqrt{10005}} \Bigg( 13591409 + \sum_{k=1}^{\infty}{\Bigg( \prod_{j=1}^{k}{\frac{-(6j-1)(2j-1)(6j-5)}{10939058860032000 j^3}} \Bigg) \cdot (545140134k + 13591409)} \Bigg) }[/math]

Let [math]\displaystyle{ P(a,b) = \prod_{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)} }[/math] and [math]\displaystyle{ Q(a,b) = \prod_{j=a}^{b-1}{10939058860032000 j^3} }[/math], so[math]\displaystyle{ \frac{1}{\pi} = \frac{1}{426880 \sqrt{10005}} \Bigg( 13591409 + \sum_{k=1}^{\infty}{\frac{P(1, k+1)}{Q(1, k+1)} \cdot (545140134k + 13591409) } \Bigg) }[/math]

Let [math]\displaystyle{ S(a,b) = \sum_{k=a}^{b-1}{\frac{P(a, k+1)}{Q(a, k+1)} \cdot (545140134k + 13591409)} }[/math] and [math]\displaystyle{ R(a,b) = Q(a,b) \cdot S(a,b) }[/math][math]\displaystyle{ \pi = \frac{426880 \sqrt{10005}}{13591409 + S(1, \infty)} }[/math]

[math]\displaystyle{ S(1, \infty) }[/math] can never be computed, so instead compute [math]\displaystyle{ S(1, n) }[/math] and as [math]\displaystyle{ n }[/math] approaches [math]\displaystyle{ \infty }[/math], the [math]\displaystyle{ \pi }[/math] approximation will get better.[math]\displaystyle{ \pi \approx \frac{426880 \sqrt{10005}}{13591409 + S(1, n)} }[/math]

From the original definition of [math]\displaystyle{ R }[/math], [math]\displaystyle{ S(a,b) = \frac{R(a,b)}{Q(a,b)} }[/math][math]\displaystyle{ \pi \approx \frac{426880 \sqrt{10005} \cdot Q(1,n)}{13591409 Q(1,n) + R(1,n)} }[/math]

Recursively computing the functions

Consider a value [math]\displaystyle{ m }[/math] such that [math]\displaystyle{ a \lt m \lt b }[/math]

• [math]\displaystyle{ P(a,b) = P(a,m) \cdot P(m,b) }[/math]

• [math]\displaystyle{ Q(a,b) = Q(a,m) \cdot Q(m,b) }[/math]
• [math]\displaystyle{ S(a,b) = S(a,m) + \frac{P(a,m)}{Q(a,m)} S(m,b) }[/math]
• [math]\displaystyle{ R(a,b) = Q(m,b) R(a,m) + P(a,m) R(m,b) }[/math]

Base case for recursion

Consider [math]\displaystyle{ b = a + 1 }[/math]

• [math]\displaystyle{ P(a, a+1) = -(6a-1)(2a-1)(6a-5) }[/math]
• [math]\displaystyle{ Q(a, a+1) = 10939058860032000 a^3 }[/math]
• [math]\displaystyle{ S(a, a+1) = \frac{P(a, a+1)}{Q(a, a+1)} \cdot (545140134a + 13591409) }[/math]
• [math]\displaystyle{ R(a, a+1) = P(a, a+1) \cdot (545140134a + 13591409) }[/math]

Python code

import decimal


def binary_split(a, b):
    if b == a + 1:
        Pab = -(6*a - 5)*(2*a - 1)*(6*a - 1)
        Qab = 10939058860032000 * a**3
        Rab = Pab * (545140134*a + 13591409)
    else:
        m = (a + b) // 2
        Pam, Qam, Ram = binary_split(a, m)
        Pmb, Qmb, Rmb = binary_split(m, b)
        
        Pab = Pam * Pmb
        Qab = Qam * Qmb
        Rab = Qmb * Ram + Pam * Rmb
    return Pab, Qab, Rab


def chudnovsky(n):
    P1n, Q1n, R1n = binary_split(1, n)
    return (426880 * decimal.Decimal(10005).sqrt() * Q1n) / (13591409*Q1n + R1n)


print(chudnovsky(2))  # 3.141592653589793238462643384

decimal.getcontext().prec = 100
for n in range(2,10):
    print(f"{n=} {chudnovsky(n)}")  # 3.14159265358979323846264338...

Notes

[math]\displaystyle{ e^{\pi \sqrt{163}} \approx 640320^3+743.99999999999925\dots }[/math]

[math]\displaystyle{ 640320^3 / 24 = 10939058860032000 }[/math]

[math]\displaystyle{ 545140134 = 163 \cdot 127 \cdot 19 \cdot 11 \cdot 7 \cdot 3^2 \cdot 2 }[/math]

[math]\displaystyle{ 13591409 = 13 \cdot 1045493 }[/math]

See also

• Ramanujan–Sato series
• Bailey–Borwein–Plouffe formula
• Borwein's algorithm
• Approximations of π