BBP-Algorithm in Python - Working with Arbitrary-precision arithmetic

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 ?

Solution

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)