The Algorithm They Used to Compute the Value of π up to a Trillion Decimal Places

The Chudnovsky algorithm and how it can be implemented in Python

One of the oldest and most challenging problems in mathematics is to compute the value of π with high accuracy and efficiency. Many methods have been devised over the centuries, ranging from simple geometric approximations, such as using polygons or circles, to sophisticated analytical formulas, such as using infinite series or integrals.

However, most of these methods are either slow, complex, or impractical, and they require a lot of computational resources and time to produce a large number of digits of π.

I’ve talked more about the computation in the following article,

In 1988, two brothers, David and Gregory Chudnovsky, published a remarkable algorithm that revolutionized the computation of π. The Chudnovsky algorithm is a fast and elegant method that can generate 14 or more digits of π for every summation step, based on some of the formulas of the Indian mathematical genius Srinivasa Ramanujan.

The Chudnovsky algorithm has been used to achieve numerous world record calculations of π, such as 2.7 trillion digits in 2009, 10 trillion digits in 2011, 22.4 trillion digits in 2016, 31.4 trillion digits in 2019, 50 trillion digits in 2020, 62.8 trillion digits in 2021, and 100 trillion digits in 2022.

In this story, I’ll try to explain the main idea and the steps of the Chudnovsky algorithm, and I will show how it can be implemented in Python code.

The algorithm is based on the following identity,

This identity is similar to some of the formulas that Ramanujan discovered for π, and it is an example of a Ramanujan-Sato series.

The main idea of this algorithm is to invert this identity and use it to compute π as follows:

The advantage of this method is that the series converges very fast, meaning that each term adds a lot of digits to the approximation of π. In fact, each term adds about 14.18 digits, which is the logarithm of the ratio of two consecutive terms. This means that, for example, to compute 100 digits of π, we only need to sum about 8 terms of the series.

The Steps

• Define the constants and the variables that are used in the series, such as,

• Define a function that computes the square root of a number using Newton’s method, which iteratively improves an initial guess until it reaches the desired accuracy.
• Define a function that computes the arctangent of a number using Taylor series, which approximates the function by adding infinitely many terms of increasing powers and decreasing coefficients.
• Then define a function that computes the terms of the series using binary splitting, which is a technique that divides the series into smaller subseries and multiplies their partial sums, avoiding the need to compute large factorials and powers.
• Compute π using the series and the functions defined above, and print the result with the desired number of digits.

The Code

# Import modules
import decimal

# Set the precision
decimal.getcontext().prec = 100 # adjust as needed

# Define constants
K = 6 # the constant term in the series
C = 640320 # the coefficient in the series
L = 13591409 # the linear term in the series
X = 1 # the exponential term in the series
M = 1 # the factorial term in the series
S = L # the sum of the series
N = 0 # the index of the series

# Define functions
def sqrt (x):
    # Square root function using Newton's method
    return x ** decimal.Decimal ('0.5')

def arctan (x):
    # Arctangent function using Taylor series
    x = decimal.Decimal (x)
    y = x # the first term
    z = x # the sum of the series
    n = 1 # the index of the series
    sign = 1 # the sign of the terms
    while y != 0:
        y *= x * x # the next power of x
        n += 2 # the next odd number
        sign *= -1 # the next sign
        z += sign * y / n # the next term
    return z

def binary_split(a, b):
    # Binary splitting function for the series
    if b == a + 1:
        # Base case
        Pab = -(6*a - 5)*(2*a - 1)*(6*a - 1)
        Qab = C**3 * a**3
        Rab = Pab * (L*a + K)
    else:
        # Recursive case
        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

# Calculate the digits of pi
P1n, Q1n, R1n = binary_split(1, N)
pi = decimal.Decimal(C) * sqrt(C) / 12 * (Q1n * L + R1n) / P1n

# Print pi
print(pi)

The Optimizations and Challenges

The Chudnovsky algorithm is a very efficient and elegant method for computing π, but it also has some optimizations and challenges that need to be considered, such as:

The algorithm requires high-precision arithmetic, which means that the numbers and operations involved have more digits than the standard ones. This can be achieved by using special libraries or modules, such as the decimal module in Python, that allow arbitrary-precision arithmetic. However, this also increases the computational complexity and the memory usage of the algorithm, and it may introduce some errors or limitations due to rounding or truncation.

It uses binary splitting, which is a technique that reduces the number of multiplications and divisions needed to compute the series. Binary splitting works by dividing the series into smaller subseries and multiplying their partial sums, avoiding the need to compute large factorials and powers. However, binary splitting also requires recursion, which is a process that calls a function within itself until a base case is reached. Recursion can be very elegant and simple, but it can also cause some problems, such as stack overflow or infinite loops, if not implemented correctly or efficiently.

The algorithm is based on some of the formulas of Ramanujan, which are very beautiful and mysterious, but also very difficult and obscure. Ramanujan was a self-taught mathematical genius who discovered many amazing results and identities, often without providing any proofs or explanations. His formulas for π are examples of Ramanujan-Sato series, which are rapidly convergent series that involve complex numbers, modular forms, and elliptic functions. These concepts and connections are very deep and fascinating, but they are also very advanced and challenging, and they require a lot of background knowledge and understanding. Read more about the mathematical genius in the following stories,

Final thoughts,

I hope that you have enjoyed this article and that you have learned something new and interesting about π and the Chudnovsky algorithm. If you want to try the algorithm yourself, you can copy and paste the code above and run it in your Python interpreter. You can also modify the code and experiment with different values and parameters. Have fun and happy mathematical coding!

References:

Borwein, Jonathan M., and Peter B. Borwein. “Ramanujan and pi.” Scientific American 258.2 (1988): 112–117.

Borwein, Peter. “The amazing number Pi.” Nieuw Archief Voor Wiskunde 1 (2000): 254–258.

Thank you so much for reading. If you liked this story don’t forget to press that clap icon as many times as you want. If you like my work and want to support me then you can Buy me a coffee ☕️. Keep following for more such stories!