-
Notifications
You must be signed in to change notification settings - Fork 0
/
googlepicrunch.py
82 lines (73 loc) · 2.25 KB
/
googlepicrunch.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
"""
Python3 program to calculate Pi using python long integers, binary
splitting and the Chudnovsky algorithm
See: http://www.craig-wood.com/nick/articles/pi-chudnovsky/ for more
info
Nick Craig-Wood <nick@craig-wood.com>
"""
import math
from gmpy2 import mpz
from time import time
import gmpy2
def pi_chudnovsky_bs(digits):
"""
Compute int(pi * 10**digits)
This is done using Chudnovsky's series with binary splitting
"""
C = 640320
C3_OVER_24 = C**3 // 24
def bs(a, b):
"""
Computes the terms for binary splitting the Chudnovsky infinite series
a(a) = +/- (13591409 + 545140134*a)
p(a) = (6*a-5)*(2*a-1)*(6*a-1)
b(a) = 1
q(a) = a*a*a*C3_OVER_24
returns P(a,b), Q(a,b) and T(a,b)
"""
if b - a == 1:
# Directly compute P(a,a+1), Q(a,a+1) and T(a,a+1)
if a == 0:
Pab = Qab = mpz(1)
else:
Pab = mpz((6*a-5)*(2*a-1)*(6*a-1))
Qab = mpz(a*a*a*C3_OVER_24)
Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)
if a & 1:
Tab = -Tab
else:
# Recursively compute P(a,b), Q(a,b) and T(a,b)
# m is the midpoint of a and b
m = (a + b) // 2
# Recursively calculate P(a,m), Q(a,m) and T(a,m)
Pam, Qam, Tam = bs(a, m)
# Recursively calculate P(m,b), Q(m,b) and T(m,b)
Pmb, Qmb, Tmb = bs(m, b)
# Now combine
Pab = Pam * Pmb
Qab = Qam * Qmb
Tab = Qmb * Tam + Pam * Tmb
return Pab, Qab, Tab
# how many terms to compute
DIGITS_PER_TERM = math.log10(C3_OVER_24/6/2/6)
N = int(digits/DIGITS_PER_TERM + 1)
# Calclate P(0,N) and Q(0,N)
P, Q, T = bs(0, N)
one_squared = mpz(10)**(2*digits)
#sqrtC = (10005*one_squared).sqrt()
sqrtC = gmpy2.isqrt(10005*one_squared)
return (Q*426880*sqrtC) // T
# The last 5 digits or pi for various numbers of digits
check_digits = {
100 : 70679,
1000 : 1989,
10000 : 75678,
100000 : 24646,
1000000 : 58151,
10000000 : 55897,
}
if __name__ == "__main__":
digits = 1000000000
pi = pi_chudnovsky_bs(digits)
print(pi)
print('\nB)\n')