patternpythonMinor
Pi-calculating program
Viewed 0 times
calculatingprogramstackoverflow
Problem
I saw this question and answer about calculating pi on Stack Overflow, and I decided to write my own program for calculating pi. I used Python and only integers (I didn't want to use floating point numbers), and used the Gauss–Legendre algorithm because it was the simplest to implement (I considered using the Borwein's algorithm, but I didn't want to calculate third roots of numbers, and the Chudnovsky algorithm seemed a little complicated, although maybe I'll give it a try). I would like to know how efficient are the algorithms above, in \$\mathcal{O}(n)\$ (n is the number of digits to calculate), and how efficient is my program?
```
# Calculate next square root approximation of the number.
def calculate_next_square_root(number, square_root):
next_square_root = ((number / square_root) + square_root) / 2
return next_square_root
# Calculate the square root of the number.
def calculate_square_root(number, digits, add):
# Start with 1, followed by (digits + add) zeros.
square_root = 1 * (10 ** (digits + add))
# Calculate next square root approximation of the number.
next_square_root = calculate_next_square_root(number=number, square_root=square_root)
while (next_square_root != square_root):
# Replace square root with next square root.
square_root = next_square_root
# Calculate next square root approximation of the number.
next_square_root = calculate_next_square_root(number=number, square_root=square_root)
return square_root
# Calculate next pi approximation.
def calculate_next_pi(a, b, t, digits, add):
next_pi = ((10 (digits + add)) ((a + b) 2)) / (4 t)
return next_pi
# Calculate pi to 5,000 digits.
digits = 5000
add = 500
a = 10 ** (digits + add)
b = calculate_square_root(number=(10 ** ((digits + add) * 2)) / 2, digits=digits, add=add)
t = (10 ** ((digits + add) * 2)) / 4
p = 1
pi = -1 # pi must be different than next_pi
next_pi = calculate_next_pi(a=a, b=b, t=t, digits=digits, add=
```
# Calculate next square root approximation of the number.
def calculate_next_square_root(number, square_root):
next_square_root = ((number / square_root) + square_root) / 2
return next_square_root
# Calculate the square root of the number.
def calculate_square_root(number, digits, add):
# Start with 1, followed by (digits + add) zeros.
square_root = 1 * (10 ** (digits + add))
# Calculate next square root approximation of the number.
next_square_root = calculate_next_square_root(number=number, square_root=square_root)
while (next_square_root != square_root):
# Replace square root with next square root.
square_root = next_square_root
# Calculate next square root approximation of the number.
next_square_root = calculate_next_square_root(number=number, square_root=square_root)
return square_root
# Calculate next pi approximation.
def calculate_next_pi(a, b, t, digits, add):
next_pi = ((10 (digits + add)) ((a + b) 2)) / (4 t)
return next_pi
# Calculate pi to 5,000 digits.
digits = 5000
add = 500
a = 10 ** (digits + add)
b = calculate_square_root(number=(10 ** ((digits + add) * 2)) / 2, digits=digits, add=add)
t = (10 ** ((digits + add) * 2)) / 4
p = 1
pi = -1 # pi must be different than next_pi
next_pi = calculate_next_pi(a=a, b=b, t=t, digits=digits, add=
Solution
I think the most efficient method is the Chudnovsky algoritm (100 million digits of Pi, in under 10 minutes!)
To know the math behind it:
For further details:http://www.craig-wood.com/nick/articles/pi-chudnovsky/
To know the math behind it:
"""
Python3 program to calculate Pi using python long integers, BINARY
splitting and the Chudnovsky algorithm
"""
import math
from gmpy2 import mpz
from time import time
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()
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 = 100
pi = pi_chudnovsky_bs(digits)
print(pi)
#raise SystemExit
for log10_digits in range(1,9):
digits = 10**log10_digits
start =time()
pi = pi_chudnovsky_bs(digits)
print("chudnovsky_gmpy_mpz_bs: digits",digits,"time",time()-start)
if digits in check_digits:
last_five_digits = pi % 100000
if check_digits[digits] == last_five_digits:
print("Last 5 digits %05d OK" % last_five_digits)
else:
print("Last 5 digits %05d wrong should be %05d" % (last_five_digits, check_digits[digits]))For further details:http://www.craig-wood.com/nick/articles/pi-chudnovsky/
Code Snippets
"""
Python3 program to calculate Pi using python long integers, BINARY
splitting and the Chudnovsky algorithm
"""
import math
from gmpy2 import mpz
from time import time
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()
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 = 100
pi = pi_chudnovsky_bs(digits)
print(pi)
#raise SystemExit
for log10_digits in range(1,9):
digits = 10**log10_digits
start =time()
pi = pi_chudnovsky_bs(digits)
print("chudnovsky_gmpy_mpz_bs: digits",digits,"time",time()-start)
if digits in check_digits:
last_five_digits = pi % 100000
if check_digits[digits] == last_five_digits:
print("Last 5 digits %05d OK" % last_five_digits)
else:
print("Last 5 digits %05d wrong should be %05d" % (last_five_digits, check_digits[digits]))Context
StackExchange Code Review Q#79974, answer score: 4
Revisions (0)
No revisions yet.