Jump to content

Chudnovsky algorithm

From Wikipedia, the free encyclopedia
This is an old revision of this page, as edited by Inputter1 (talk | contribs) at 02:02, 23 March 2025 (I added more code that can be used to calculate Pi with the Chudnovsky algorithm). The present address (URL) is a permanent link to this revision, which may differ significantly from the current revision.

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 , the j-function , and on the following rapidly convergent generalized hypergeometric series:[13]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 .[15]

Optimizations

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

Binary splitting

A factor of can be taken out of the sum and simplified to


Let , and substitute that into the sum.


can be simplified to , so

from the original definition of , so

This definition of is not defined for , so compute the first term of the sum and use the new definition of

Let and , so

Let and

can never be computed, so instead compute and as approaches , the approximation will get better.

From the original definition of ,

Recursively computing the functions

Consider a value such that

Base case for recursion

Consider

Python code

#Note: For extreme calculations, other code can be used to run on a GPU, which is much faster than this. Efficiency improvements lead to the possibility of the calculation of trillions of digits.
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...

Very Efficient Java Code

This can calculate to millions of digits in less than an hour, depending on the device. Incredibly powerful systems with dozens of GPUs and even more efficient code lead to the trillions of digits that have been calculated to this date.

import java.math.*;
public class PiCalculation {
    public static BigInteger[] compute(int start, int end) {
        if (end == start + 1) {
            BigInteger aValue = BigInteger.valueOf(start);
            BigInteger productP = BigInteger.valueOf(-(6 * start - 5))
                    .multiply(BigInteger.valueOf(2 * start - 1))
                    .multiply(BigInteger.valueOf(6 * start - 1));
            BigInteger productQ = new BigInteger("10939058860032000").multiply(aValue.pow(3));
            BigInteger sumR = productP.multiply(BigInteger.valueOf(545140134).multiply(aValue)
                    .add(BigInteger.valueOf(13591409)));
            return new BigInteger[]{productP, productQ, sumR};
        } else {
            int mid = (start + end) / 2;
            BigInteger[] leftPQR = computePQR(start, mid);
            BigInteger[] rightPQR = computePQR(mid, end);
            BigInteger productP = leftPQR[0].multiply(rightPQR[0]);
            BigInteger productQ = leftPQR[1].multiply(rightPQR[1]);
            BigInteger sumR = rightPQR[1].multiply(leftPQR[2]).add(leftPQR[0].multiply(rightPQR[2]));
            return new BigInteger[]{productP, productQ, sumR};
        }
    }

    public static BigDecimal computePi(int terms, MathContext mc) {
        BigInteger[] PQR = computePQR(1, terms);
        BigInteger productQ = PQR[1];
        BigInteger sumR = PQR[2];

        BigDecimal productQBD = new BigDecimal(productQ);
        BigDecimal sumRBD = new BigDecimal(sumR);
        BigDecimal denominator = BigDecimal.valueOf(13591409).multiply(productQBD).add(sumRBD);

        return new BigDecimal("426880")
                .multiply(sqrt(new BigDecimal("10005"), 10))
                .multiply(productQBD)
                .divide(denominator, mc);
    }

    public static BigDecimal sqrt(BigDecimal value, int precision) {
        BigDecimal x = new BigDecimal(Math.sqrt(value.doubleValue()));
        BigDecimal two = BigDecimal.valueOf(2);
        MathContext mc = new MathContext(precision + 9992);
        for (int i = 0; i < precision; i++) {
            x = x.add(value.divide(x, mc)).divide(two, mc);
        }
        return x;
    }

    public static void main(String[] args){
        int digits = 100000;
        String pi = computePi(digits, new MathContext(digits + 10);
        // print the 9000th digit of pi, for example
        System.out.print(pi.charAt(9001);

}

Notes

See also

  • [1] 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. ^ a b 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.