Submit your draft for review!
Fast method for calculating the digits of π
The Chudnovsky algorithm is a fast method for calculating the digits of π Ramanujan 's π formulaeChudnovsky brothers in 1988.[ 1]
It was used in the world record calculations of 2.7 trillion digits of π in December 2009,[ 2] [ 3] [ 4] [ 5] [ 6] [ 7] [ 8] [ 9]
Algorithm
The algorithm is based on the negated Heegner number
d
=
−
163
{\displaystyle d=-163}
j -function
j
(
1
+
i
163
2
)
=
−
640320
3
{\displaystyle j\left({\tfrac {1+i{\sqrt {163}}}{2}}\right)=-640320^{3}}
generalized hypergeometric series :[ 2]
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}}}}
[ 10]
This identity is similar to some of Ramanujan 's formulas involving π ,[ 2] Ramanujan–Sato series .
The time complexity of the algorithm is
O
(
n
(
log
n
)
3
)
{\displaystyle O\left(n(\log n)^{3}\right)}
[ 11]
Optimizations
The optimization technique used for the world record computations is called binary splitting . Here is a video that explains the binary splitting of the algorithm: [ 12]
Binary splitting
A factor of
1
/
640320
3
/
2
{\textstyle 1/{640320^{3/2}}}
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}}}}
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}}}}
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)}}}
−
(
6
n
−
1
)
(
2
n
−
1
)
(
6
n
−
5
)
10939058860032000
n
3
{\displaystyle {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}}
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}
f
{\displaystyle f}
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}
n
=
0
{\displaystyle n=0}
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)}}
Q
(
a
,
b
)
=
∏
j
=
a
b
−
1
10939058860032000
j
3
{\displaystyle Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}}
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)}}
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 )}
S
(
1
,
n
)
{\displaystyle S(1,n)}
n
{\displaystyle n}
∞
{\displaystyle \infty }
π
{\displaystyle \pi }
π
≈
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
Consider a value
m
{\displaystyle m}
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
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
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 ):
P1n , Q1n , R1n = binary_split ( 1 , n )
return ( 426880 * decimal . Decimal ( 10005 ) . sqrt () * Q1n ) / ( 13591409 * Q1n + R1n )
print ( chudnovsky ( 2 )) # 3.141592653589793238462643384
Notes
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
References
^ Chudnovsky, David; Chudnovsky, Gregory (1988), Approximation and complex multiplication according to Ramanujan , Ramanujan revisited: proceedings of the centenary conference ^ a b c 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 ^ 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 .^ 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, How is π calculated to trillions of digits? {{citation }}: CS1 maint: url-status (link ) 
