Chudnovsky algorithm

The Chudnovsky algorithm is a fast method for calculating the digits of $π$ . It was published by the Chudnovsky brothers in 1989,^[1] and was used in the world record calculations of 2.7 trillion digits of $π$ in December 2009,^[2] 5 trillion digits in August 2010,^[3] 10 trillion digits in October 2011,^[4]^[5] 12.1 trillion digits in December 2013^[6] and 22.4 trillion digits of $π$ in November 2016.^[7]

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:^[2]

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

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

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

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

\pi =C\left(\sum _{k=0}^{\infty }{\frac {M_{k}\cdot L_{k}}{X_{k}}}\right)^{-1}

, where:

C=426880{\sqrt {10005}},\quad \quad M_{k}={\frac {(6k)!}{(3k)!(k!)^{3}}},\quad \quad L_{k}=545140134k+13591409,\quad \quad X_{k}=(-262537412640768000)^{k}

The terms M_k, L_k, and X_k satisfy the following recurrences and can be computed as such:

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

The computation of M_k can be further optimized by introducing an additional term K_k as follows:

{\begin{alignedat}{4}K_{k+1}&=K_{k}+12\,\,&&{\textrm {where}}\,\,K_{0}&&=6\\[4pt]M_{k+1}&=M_{k}\cdot \left({\frac {K_{k}^{3}-16K_{k}}{(k+1)^{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 $π$ ,^[2] and is an example of a Ramanujan–Sato series.

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

Example: Python Implementation

$π$ can be computed to any precision using the above algorithm in any environment which supports arbitrary-precision arithmetic. As an example, here is a Python implementation:

from decimal import Decimal as Dec, getcontext as gc

def PI(maxK=70, prec=1008, disp=1007): # parameter defaults chosen to gain 1000+ digits within a few seconds
    gc().prec = prec
    K, M, L, X, S = 6, 1, 13591409, 1, 13591409
    for k in range(1, maxK+1):
        M = (K**3 - 16*K) * M // k**3 
        L += 545140134
        X *= -262537412640768000
        S += Dec(M * L) / X
        K += 12
    pi = 426880 * Dec(10005).sqrt() / S
    pi = Dec(str(pi)[:disp]) # drop few digits of precision for accuracy
    print("PI(maxK=%d iterations, gc().prec=%d, disp=%d digits) =\n%s" % (maxK, prec, disp, pi))
    return pi

Pi = PI()
print("\nFor greater precision and more digits (takes a few extra seconds) - Try")
print("Pi = PI(317,4501,4500)") 
print("Pi = PI(353,5022,5020)")

References

↑ Chudnovsky, David V.; Chudnovsky, Gregory V. (1989), "The Computation of Classical Constants", Proceedings of the National Academy of Sciences of the United States of America, 86 (21): 8178–8182, doi:10.1073/pnas.86.21.8178, ISSN 0027-8424, JSTOR 34831, PMC 298242, PMID 16594075
1 2 3 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
↑ Geeks slice pi to 5 trillion decimal places, Australian Broadcasting Corporation, August 6, 2010
↑ 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
↑ Aron, Jacob (March 14, 2012), "Constants clash on pi day", New Scientist
↑ Yee, Alexander J.; Kondo, Shigeru (28 December 2013). "12.1 Trillion Digits of Pi". www.numberworld.org.
↑ "22.4 Trillion Digits of Pi". www.numberworld.org.
↑ "y-cruncher - Formulas". www.numberworld.org. Retrieved 2018-02-25.

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 V.; Chudnovsky, Gregory V. (1989), "The Computation of Classical Constants", Proceedings of the National Academy of Sciences of the United States of America, 86 (21): 8178–8182, doi:10.1073/pnas.86.21.8178, ISSN 0027-8424, JSTOR 34831, PMC 298242, PMID 16594075

[baruah-2] 1 2 3 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] Geeks slice pi to 5 trillion decimal places, Australian Broadcasting Corporation, August 6, 2010

[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

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

[6] Yee, Alexander J.; Kondo, Shigeru (28 December 2013). "12.1 Trillion Digits of Pi". www.numberworld.org.

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

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

Chudnovsky algorithm

Example: Python Implementation

See also

References