In python, I wrote a program to work out the value of pi using the Chudnovsky Algorithm. It works on numbers under 800. However, if I use a number above 800, it returns a decimal.InvalidOperation error. This is my code:
from math import factorial
from decimal import Decimal, getcontext
getcontext().prec = 1000
pi_input = input("How Many Digits Of Pi Are To Be Represented: ")
num = int(pi_input)
def cal(n):
t = Decimal(0)
pi = Decimal(0)
deno = Decimal(0)
for k in range(n):
t = ((-1) ** k) * (factorial(6 * k)) * (13591409 + 545140134 * k)
deno = factorial(3 * k) * (factorial(k) ** 3) * (640320 ** (3 * k))
pi += Decimal(t) / Decimal(deno)
pi = pi * Decimal(12) / Decimal(640320 ** Decimal(1.5))
pi = 1 / pi
return round(pi, n)
print(cal(num))
Could Anyone Please Help Me?
I wouldn't recommend your method. Try using the gmpy module like this:
import sys
import math
import os.path
from gmpy2 import mpz, sqrt
from time import time, sleep
print_pi = input("Should The Value Of Pi Be Printed: ").lower()
print_pi_bool = 0
if "y" in print_pi:
print_pi_bool = True
elif "n" in print_pi:
print_pi = False
else:
print("Incorrect Input. Please Try Again.")
sys.exit()
def sqrt(n, one):
"""
Return the square root of n as a fixed point number with the one
passed in. It uses a second order Newton-Raphson convgence. This
doubles the number of significant figures on each iteration.
"""
# Use floating point arithmetic to make an initial guess
floating_point_precision = 10**16
n_float = float((n * floating_point_precision) // one) / floating_point_precision
x = (int(floating_point_precision * math.sqrt(n_float)) * one) // floating_point_precision
n_one = n * one
count = 0
while 1:
count += 1
print(count)
x_old = x
x = (x + n_one // x) // 2
if x == x_old:
break
return x
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 = sqrt(10005*one_squared, 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 = 100
pi = pi_chudnovsky_bs(digits)
#raise SystemExit
for log10_digits in range(1,11):
digits = 10**log10_digits
start =time()
pi = pi_chudnovsky_bs(digits)
pi = str(pi)
pi = pi[:(len(str(pi)) // 2) + 1]
length = int(len(str(pi))) - 1
print("Chudnovsky Binary Splitting Using GMPY: Digits",f"{digits:,}","\n--------------",time()-start,"--------------")
print("Length Of " + f"{length:,}" + " Digits")
It's my first answer, so I hope it helps!