Chudnovsky algorithm

The Chudnovsky algorithm is a fast method for calculating the digits of π, based on Ramanujan's π formulae. Published by the Chudnovsky brothers in 1988,^[1] it was used to calculate π to a billion decimal places.^[2]

It was used in the world record calculations of 2.7 trillion digits of π in December 2009,^[3] 10 trillion digits in October 2011,^[4][5] 22.4 trillion digits in November 2016,^[6] 31.4 trillion digits in September 2018–January 2019,^[7] 50 trillion digits on January 29, 2020,^[8] 62.8 trillion digits on August 14, 2021,^[9] 100 trillion digits on March 21, 2022,^[10] 105 trillion digits on March 14, 2024,^[11] and 202 trillion digits on June 28, 2024.^[12]

Algorithm

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:^[13]$${\displaystyle {\frac {1}{\pi }}=12\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320)^{3k+3/2}}}}$$A detailed proof of this formula can be found here: ^[14]

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

The time complexity of the algorithm is ${\displaystyle O\left(n(\log n)^{3}\right)}$.^[15]

Optimizations

The optimization technique used for the world record computations is called binary splitting.^[16]

Binary splitting

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

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

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

${\displaystyle f(0)=1}$ from the original definition of ${\displaystyle f}$, so$${\displaystyle f(n)=\prod _{j=1}^{n}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}}$$

This definition of ${\displaystyle f}$ is not defined for ${\displaystyle n=0}$, so compute the first term of the sum and use the new definition of ${\displaystyle f}$$${\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)}{10939058860032000j^{3}}}{\Bigg )}\cdot (545140134k+13591409)}{\Bigg )}}$$

Let ${\displaystyle P(a,b)=\prod _{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)}}$ and ${\displaystyle Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}}$, so$${\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 )}}$$

Let ${\displaystyle S(a,b)=\sum _{k=a}^{b-1}{{\frac {P(a,k+1)}{Q(a,k+1)}}\cdot (545140134k+13591409)}}$ and ${\displaystyle R(a,b)=Q(a,b)\cdot S(a,b)}$$${\displaystyle \pi ={\frac {426880{\sqrt {10005}}}{13591409+S(1,\infty )}}}$$

${\displaystyle S(1,\infty )}$ can never be computed, so instead compute ${\displaystyle S(1,n)}$ and as ${\displaystyle n}$ approaches ${\displaystyle \infty }$, the ${\displaystyle \pi }$ approximation will get better.$${\displaystyle \pi =\lim _{n\to \infty }{\frac {426880{\sqrt {10005}}}{13591409+S(1,n)}}}$$

From the original definition of ${\displaystyle R}$, ${\displaystyle S(a,b)={\frac {R(a,b)}{Q(a,b)}}}$$${\displaystyle \pi =\lim _{n\to \infty }{\frac {426880{\sqrt {10005}}\cdot Q(1,n)}{13591409Q(1,n)+R(1,n)}}}$$

Recursively computing the functions

Consider a value ${\displaystyle m}$ such that ${\displaystyle a<m<b}$

• ${\displaystyle P(a,b)=P(a,m)\cdot P(m,b)}$

• ${\displaystyle Q(a,b)=Q(a,m)\cdot Q(m,b)}$
• ${\displaystyle S(a,b)=S(a,m)+{\frac {P(a,m)}{Q(a,m)}}S(m,b)}$
• ${\displaystyle R(a,b)=Q(m,b)R(a,m)+P(a,m)R(m,b)}$

Base case for recursion

Consider ${\displaystyle b=a+1}$

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

Python code

#Note: For extreme calculations, other code can be used to run on a GPU, which is much faster than this.
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):
    """Chudnovsky algorithm."""
    P1n, Q1n, R1n = binary_split(1, n)
    return (426880 * decimal.Decimal(10005).sqrt() * Q1n) / (13591409*Q1n + R1n)


print(f"1 = {chudnovsky(2)}")  # 3.141592653589793238462643384

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

Notes

${\displaystyle e^{\pi {\sqrt {163}}}\approx 640320^{3}+743.99999999999925\dots }$

${\displaystyle 640320^{3}/24=10939058860032000}$

${\displaystyle 545140134=163\cdot 127\cdot 19\cdot 11\cdot 7\cdot 3^{2}\cdot 2}$

${\displaystyle 13591409=13\cdot 1045493}$

See also

• Mathematics portal

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

External links

• How is π calculated to trillions of digits?

References

1. Chudnovsky, David; Chudnovsky, Gregory (1988), Approximation and complex multiplication according to Ramanujan, Ramanujan revisited: proceedings of the centenary conference
2. Warsi, Karl; Dangerfield, Jan; Farndon, John; Griffiths, Johny; Jackson, Tom; Patel, Mukul; Pope, Sue; Parker, Matt (2019). The Math Book: Big Ideas Simply Explained. New York: Dorling Kindersley Limited. p. 65. ISBN 978-1-4654-8024-8.
3. Baruah, Nayandeep Deka; Berndt, Bruce C.; Chan, Heng Huat (2009-08-01). "Ramanujan's Series for 1/π: A Survey". American Mathematical Monthly. 116 (7): 567–587. doi:10.4169/193009709X458555.
4. Yee, Alexander; Kondo, Shigeru (2011), 10 Trillion Digits of Pi: A Case Study of summing Hypergeometric Series to high precision on Multicore Systems, Technical Report, Computer Science Department, University of Illinois, hdl:2142/28348
5. Aron, Jacob (March 14, 2012), "Constants clash on pi day", New Scientist
6. "22.4 Trillion Digits of Pi". www.numberworld.org.
7. "Google Cloud Topples the Pi Record". www.numberworld.org/.
8. "The Pi Record Returns to the Personal Computer". www.numberworld.org/.
9. "Pi-Challenge - Weltrekordversuch der FH Graubünden - FH Graubünden". www.fhgr.ch. Retrieved 2021-08-17.
10. "Calculating 100 trillion digits of pi on Google Cloud". cloud.google.com. Retrieved 2022-06-10.
11. Yee, Alexander J. (2024-03-14). "Limping to a new Pi Record of 105 Trillion Digits". NumberWorld.org. Retrieved 2024-03-16.
12. Ranous, Jordan (2024-06-28). "StorageReview Lab Breaks Pi Calculation World Record with Over 202 Trillion Digits". StorageReview.com. Retrieved 2024-07-20.
13. Baruah, Nayandeep Deka; Berndt, Bruce C.; Chan, Heng Huat (2009), "Ramanujan's series for 1/π: a survey", American Mathematical Monthly, 116 (7): 567–587, doi:10.4169/193009709X458555, JSTOR 40391165, MR 2549375
14. Milla, Lorenz (2018), A detailed proof of the Chudnovsky formula with means of basic complex analysis, arXiv:1809.00533
15. "y-cruncher - Formulas". www.numberworld.org. Retrieved 2018-02-25.
16. Brent, Richard P.; Zimmermann, Paul (2010). Modern Computer Arithmetic. Vol. 18. Cambridge University Press. doi:10.1017/CBO9780511921698. ISBN 978-0-511-92169-8.