Reputation: 113
I am writing on a term paper about the calculation of pi. While i have finished the theoretical site, i`m now struggeling implementing the BBP-Algorithm in Python.
You can find the BBP-Algorithm here: http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
And this is my Implementation in Python:
from sympy.mpmath import *
pi = mpf(0)
mp.dps = 30;
for n in range(0,500):
pi = pi + mpf((1/16**n)*(4/(8*n+1)- 2/(8*n+4)- 1/(8*n+5)- 1/(8*n+6)) )
print(pi)
My Problem is, no matter how high i set k or how high i set the decimal places for pi i can`t get a higher precision than 16 digits.
I used mpmath for higher precision, because i encountered some problems before.
How can i improve my code, so i can get more digits ?
Upvotes: 5
Views: 1730
Reputation: 81
Used Willem Van Onsem's response
from mpmath import mp, mpf, nstr
# Bailey Borwein Plouffe algorithm for approximate pi computation
def BBP_pi(n: int, digits: int) -> str:
assert isinstance(n, int)
assert isinstance(digits, int)
_pi = mpf(0)
mp.dps = digits
for k in range(0,n):
u = mpf(4)/mpf(8*k+1)-mpf(2)/mpf(8*k+4)-mpf(1)/mpf(8*k+5)-mpf(1)/mpf(8*k+6)
d = mpf(16.0)**k
_pi += u/d
_pi_str = nstr(_pi, 1_000_000)
return _pi_str
def write_to_file(pi_str: str) -> None:
assert isinstance(pi_str, str)
with open(file='pi.txt', mode='w') as file:
file.write(pi_str)
def display_pi(pi_str:str, characters_per_row: int) -> None:
assert isinstance(pi_str, str)
assert isinstance(characters_per_row, int)
pi_str_size = len(pi_str)
for i,j in zip(range(0, pi_str_size-characters_per_row, characters_per_row), range(characters_per_row, pi_str_size, characters_per_row)):
print(pi_str[i:j])
digits_precision = 1_000_000
n_iterations = 100_000
approximate_pi = BBP_pi(n=n_iterations, digits=digits_precision)
print(f"Vaule of Pi is : {approximate_pi}")
Upvotes: 0
Reputation: 476594
By default, python will use the standard floating points as defined by IEEE-754. This has a precision of some 12 digits and can represent numbers as lows as 2-1022, now you can resolve this problem by calling the mpf
operator earlier in the process, a more readable and more precise version is thus:
from sympy.mpmath import *
pi = mpf(0)
mp.dps = 30;
for n in range(0,500):
u = 4.0/(8*n+1)-2.0/(8*n+4)-1.0/(8*n+5)-1.0/(8*n+6)
u = mpf(u)
d = mpf(16.0/1.0)**n
pi += u/d
print(pi)
This however still has problems when you calculate the u
part. To do it exactly, you can use:
from sympy.mpmath import *
pi = mpf(0)
mp.dps = 50;
for n in range(0,50000):
u = mpf(4)/mpf(8*n+1)-mpf(2)/mpf(8*n+4)-mpf(1)/mpf(8*n+5)-mpf(1)/mpf(8*n+6)
d = mpf(16.0)**n
pi += u/d
print(pi)
Which is correct for the first 50 digits:
3.1415926535 8979323846 2643383279 502884197 1693993751
(spaces added)
Upvotes: 4