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] and 105 trillion digits on March 14, 2024.[11]
The optimization technique used for the world record computations is called binary splitting.[14]
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 isn't 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
importdecimaldefbinary_split(a,b):ifb==a+1:Pab=-(6*a-5)*(2*a-1)*(6*a-1)Qab=10939058860032000*a**3Rab=Pab*(545140134*a+13591409)else:m=(a+b)//2Pam,Qam,Ram=binary_split(a,m)Pmb,Qmb,Rmb=binary_split(m,b)Pab=Pam*PmbQab=Qam*QmbRab=Qmb*Ram+Pam*RmbreturnPab,Qab,Rabdefchudnovsky(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.141592653589793238462643384decimal.getcontext().prec=100# number of digits of decimal precisionforninrange(2,10):print(f"{n} = {chudnovsky(n)}")# 3.14159265358979323846264338...