Fast method for calculating the digits of π
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]
Algorithm [ edit ]
The algorithm is based on the negated Heegner number
d
=
−
163
{\displaystyle d=-163}
, the j -function
j
(
1
+
i
163
2
)
=
−
640320
3
{\displaystyle j\left({\tfrac {1+i{\sqrt {163}}}{2}}\right)=-640320^{3}}
, and on the following rapidly convergent generalized hypergeometric series :[12]
1
π
=
12
∑
k
=
0
∞
(
−
1
)
k
(
6
k
)
!
(
545140134
k
+
13591409
)
(
3
k
)
!
(
k
!
)
3
(
640320
)
3
k
+
3
/
2
{\displaystyle {\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:
[13]
This identity is similar to some of Ramanujan 's formulas involving π ,[12] and is an example of a Ramanujan–Sato series .
The time complexity of the algorithm is
O
(
n
(
log
n
)
3
)
{\displaystyle O\left(n(\log n)^{3}\right)}
.[14]
Optimizations [ edit ]
The optimization technique used for the world record computations is called binary splitting .[15]
Binary splitting [ edit ]
A factor of
1
/
640320
3
/
2
{\textstyle 1/{640320^{3/2}}}
can be taken out of the sum and simplified to
1
π
=
1
426880
10005
∑
k
=
0
∞
(
−
1
)
k
(
6
k
)
!
(
545140134
k
+
13591409
)
(
3
k
)
!
(
k
!
)
3
(
640320
)
3
k
{\displaystyle {\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
)
=
(
−
1
)
n
(
6
n
)
!
(
3
n
)
!
(
n
!
)
3
(
640320
)
3
n
{\displaystyle f(n)={\frac {(-1)^{n}(6n)!}{(3n)!(n!)^{3}(640320)^{3n}}}}
, and substitute that into the sum.
1
π
=
1
426880
10005
∑
k
=
0
∞
f
(
k
)
⋅
(
545140134
k
+
13591409
)
{\displaystyle {\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}\sum _{k=0}^{\infty }{f(k)\cdot (545140134k+13591409)}}
f
(
n
)
f
(
n
−
1
)
{\displaystyle {\frac {f(n)}{f(n-1)}}}
can be simplified to
−
(
6
n
−
1
)
(
2
n
−
1
)
(
6
n
−
5
)
10939058860032000
n
3
{\displaystyle {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}}
, so
f
(
n
)
=
f
(
n
−
1
)
⋅
−
(
6
n
−
1
)
(
2
n
−
1
)
(
6
n
−
5
)
10939058860032000
n
3
{\displaystyle f(n)=f(n-1)\cdot {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}}
f
(
0
)
=
1
{\displaystyle f(0)=1}
from the original definition of
f
{\displaystyle f}
, so
f
(
n
)
=
∏
j
=
1
n
−
(
6
j
−
1
)
(
2
j
−
1
)
(
6
j
−
5
)
10939058860032000
j
3
{\displaystyle f(n)=\prod _{j=1}^{n}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}}
This definition of
f
{\displaystyle f}
isn't defined for
n
=
0
{\displaystyle n=0}
, so compute the first term of the sum and use the new definition of
f
{\displaystyle f}
1
π
=
1
426880
10005
(
13591409
+
∑
k
=
1
∞
(
∏
j
=
1
k
−
(
6
j
−
1
)
(
2
j
−
1
)
(
6
j
−
5
)
10939058860032000
j
3
)
⋅
(
545140134
k
+
13591409
)
)
{\displaystyle {\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
)
=
∏
j
=
a
b
−
1
−
(
6
j
−
1
)
(
2
j
−
1
)
(
6
j
−
5
)
{\displaystyle P(a,b)=\prod _{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)}}
and
Q
(
a
,
b
)
=
∏
j
=
a
b
−
1
10939058860032000
j
3
{\displaystyle Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}}
, so
1
π
=
1
426880
10005
(
13591409
+
∑
k
=
1
∞
P
(
1
,
k
+
1
)
Q
(
1
,
k
+
1
)
⋅
(
545140134
k
+
13591409
)
)
{\displaystyle {\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
)
=
∑
k
=
a
b
−
1
P
(
a
,
k
+
1
)
Q
(
a
,
k
+
1
)
⋅
(
545140134
k
+
13591409
)
{\displaystyle 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
)
⋅
S
(
a
,
b
)
{\displaystyle R(a,b)=Q(a,b)\cdot S(a,b)}
π
=
426880
10005
13591409
+
S
(
1
,
∞
)
{\displaystyle \pi ={\frac {426880{\sqrt {10005}}}{13591409+S(1,\infty )}}}
S
(
1
,
∞
)
{\displaystyle S(1,\infty )}
can never be computed, so instead compute
S
(
1
,
n
)
{\displaystyle S(1,n)}
and as
n
{\displaystyle n}
approaches
∞
{\displaystyle \infty }
, the
π
{\displaystyle \pi }
approximation will get better.
π
≈
426880
10005
13591409
+
S
(
1
,
n
)
{\displaystyle \pi \approx {\frac {426880{\sqrt {10005}}}{13591409+S(1,n)}}}
From the original definition of
R
{\displaystyle R}
,
S
(
a
,
b
)
=
R
(
a
,
b
)
Q
(
a
,
b
)
{\displaystyle S(a,b)={\frac {R(a,b)}{Q(a,b)}}}
π
≈
426880
10005
⋅
Q
(
1
,
n
)
13591409
Q
(
1
,
n
)
+
R
(
1
,
n
)
{\displaystyle \pi \approx {\frac {426880{\sqrt {10005}}\cdot Q(1,n)}{13591409Q(1,n)+R(1,n)}}}
Recursively computing the functions [ edit ]
Consider a value
m
{\displaystyle m}
such that
a
<
m
<
b
{\displaystyle a<m<b}
P
(
a
,
b
)
=
P
(
a
,
m
)
⋅
P
(
m
,
b
)
{\displaystyle P(a,b)=P(a,m)\cdot P(m,b)}
Q
(
a
,
b
)
=
Q
(
a
,
m
)
⋅
Q
(
m
,
b
)
{\displaystyle Q(a,b)=Q(a,m)\cdot Q(m,b)}
S
(
a
,
b
)
=
S
(
a
,
m
)
+
P
(
a
,
m
)
Q
(
a
,
m
)
S
(
m
,
b
)
{\displaystyle 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
)
{\displaystyle R(a,b)=Q(m,b)R(a,m)+P(a,m)R(m,b)}
Base case for recursion [ edit ]
Consider
b
=
a
+
1
{\displaystyle b=a+1}
P
(
a
,
a
+
1
)
=
−
(
6
a
−
1
)
(
2
a
−
1
)
(
6
a
−
5
)
{\displaystyle P(a,a+1)=-(6a-1)(2a-1)(6a-5)}
Q
(
a
,
a
+
1
)
=
10939058860032000
a
3
{\displaystyle Q(a,a+1)=10939058860032000a^{3}}
S
(
a
,
a
+
1
)
=
P
(
a
,
a
+
1
)
Q
(
a
,
a
+
1
)
⋅
(
545140134
a
+
13591409
)
{\displaystyle S(a,a+1)={\frac {P(a,a+1)}{Q(a,a+1)}}\cdot (545140134a+13591409)}
R
(
a
,
a
+
1
)
=
P
(
a
,
a
+
1
)
⋅
(
545140134
a
+
13591409
)
{\displaystyle R(a,a+1)=P(a,a+1)\cdot (545140134a+13591409)}
Python code [ edit ]
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...
e
π
163
≈
640320
3
+
743.99999999999925
…
{\displaystyle e^{\pi {\sqrt {163}}}\approx 640320^{3}+743.99999999999925\dots }
640320
3
/
24
=
10939058860032000
{\displaystyle 640320^{3}/24=10939058860032000}
545140134
=
163
⋅
127
⋅
19
⋅
11
⋅
7
⋅
3
2
⋅
2
{\displaystyle 545140134=163\cdot 127\cdot 19\cdot 11\cdot 7\cdot 3^{2}\cdot 2}
13591409
=
13
⋅
1045493
{\displaystyle 13591409=13\cdot 1045493}
See also [ edit ]
References [ edit ]
^ Chudnovsky, David; Chudnovsky, Gregory (1988), Approximation and complex multiplication according to Ramanujan , Ramanujan revisited: proceedings of the centenary conference
^ 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 .
^
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 .
^ 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 .
^ Yee, Alexander J. (2024-03-14). "Limping to a new Pi Record of 105 Trillion Digits" . NumberWorld.org . Retrieved 2024-03-16 .
^ Jump up to: 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
^ 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