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]

It 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] 62.8 trillion digits on August 14, 2021,[8] and 100 trillion digits on March 21, 2022.[9]

Algorithm

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

{\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: [10]

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

Optimizations

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

Binary splitting

A factor of ${\textstyle 1/{640320^{3/2}}}$ can be taken out of the sum and simplified to

{\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320)^{3k}}}

Let $f(n)={\frac {(-1)^{n}(6n)!}{(3n)!(n!)^{3}(640320)^{3n}}}$ , and substitute that into the sum.

{\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}\sum _{k=0}^{\infty }{f(k)\cdot (545140134k+13591409)}

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

f(n)=f(n-1)\cdot {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}

$f(0)=1$ from the original definition of $f$ , so

f(n)=\prod _{j=1}^{n}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}

This definition of $f$ isn't defined for $n=0$ , so compute the first term of the sum and use the new definition of $f$

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

{\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 $S(a,b)=\sum _{k=a}^{b-1}{{\frac {P(a,k+1)}{Q(a,k+1)}}\cdot (545140134k+13591409)}$ and $R(a,b)=Q(a,b)\cdot S(a,b)$

\pi ={\frac {426880{\sqrt {10005}}}{13591409+S(1,\infty )}}

$S(1,\infty )$ can never be computed, so instead compute $S(1,n)$ and as $n$ approaches $\infty$ , the $\pi$ approximation will get better.

\pi \approx {\frac {426880{\sqrt {10005}}}{13591409+S(1,n)}}

From the original definition of $R$ , $S(a,b)={\frac {R(a,b)}{Q(a,b)}}$

\pi \approx {\frac {426880{\sqrt {10005}}\cdot Q(1,n)}{13591409Q(1,n)+R(1,n)}}

Recursively computing the functions

Consider a value $m$ such that $a<m<b$

$P(a,b)=P(a,m)\cdot P(m,b)$

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

Base case for recursion

Consider $b=a+1$

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

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

Notes

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

640320^{3}/24=10939058860032000

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

13591409=13\cdot 1045493

References

Chudnovsky, David; Chudnovsky, Gregory (1988), Approximation and complex multiplication according to Ramanujan, Ramanujan revisited: proceedings of the centenary conference
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
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
Aron, Jacob (March 14, 2012), "Constants clash on pi day", New Scientist
"22.4 Trillion Digits of Pi". www.numberworld.org.
"Google Cloud Topples the Pi Record". www.numberworld.org/.
"The Pi Record Returns to the Personal Computer". www.numberworld.org/.
"Pi-Challenge - Weltrekordversuch der FH Graubünden - FH Graubünden". www.fhgr.ch. Retrieved 2021-08-17.
"Calculating 100 trillion digits of pi on Google Cloud". cloud.google.com. Retrieved 2022-06-10.
Milla, Lorenz (2018), A detailed proof of the Chudnovsky formula with means of basic complex analysis, arXiv:1809.00533
"y-cruncher - Formulas". www.numberworld.org. Retrieved 2018-02-25.
Rayton, Joshua (Sep 2023), How is π calculated to trillions of digits?, YouTube

This article is issued from Wikipedia. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.

[1] Chudnovsky, David; Chudnovsky, Gregory (1988), Approximation and complex multiplication according to Ramanujan, Ramanujan revisited: proceedings of the centenary conference

[baruah-2] 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

[3] 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

[4] Aron, Jacob (March 14, 2012), "Constants clash on pi day", New Scientist

[5] "22.4 Trillion Digits of Pi". www.numberworld.org.

[6] "Google Cloud Topples the Pi Record". www.numberworld.org/.

[7] "The Pi Record Returns to the Personal Computer". www.numberworld.org/.

[8] "Pi-Challenge - Weltrekordversuch der FH Graubünden - FH Graubünden". www.fhgr.ch. Retrieved 2021-08-17.

[9] "Calculating 100 trillion digits of pi on Google Cloud". cloud.google.com. Retrieved 2022-06-10.

[10] Milla, Lorenz (2018), A detailed proof of the Chudnovsky formula with means of basic complex analysis, arXiv:1809.00533

[11] "y-cruncher - Formulas". www.numberworld.org. Retrieved 2018-02-25.

[12] Rayton, Joshua (Sep 2023), How is π calculated to trillions of digits?, YouTube