# [Python-modules-team] Bug#530381: Bug#530381: python-mpmath: strangely high accuracy in computing pi

Sanjoy Mahajan sanjoy at MIT.EDU
Fri Jun 5 12:54:33 UTC 2009

```> This problem only happens when gmpy is installed: could you please
> verify that removing it solves the issue?

I removed python-gmpy and the problem seemed to be fixed.  Here is the
number of digits at each iteration (with mp.dps=4000):

3.0
8.1
18.7
40.3
83.6
170.6
345.0
694.0
1392.2
2789.0
3997.3
3996.2
3995.8

But when I set mp.dps = 10000 and it failed as before.  So I wrapped the
computation in a short binary search to find the lowest failing value of
mp.dps.  The output of that test is:

Okay:   100  # with mp.dps=100, the computation of pi was fine
Okay:   200  # etc.
Okay:   400
Okay:   800
Okay:  1600
Okay:  3200
Fail:  6400  # but with mp.dps=6400, digits of accuracy went to infinity
Okay:  4800
Okay:  5600
Okay:  6000
Okay:  6200
Okay:  6300
Fail:  6350
Fail:  6325
Fail:  6312
Okay:  6306
Okay:  6309
Okay:  6310
Fail:  6311

So mps.dps=6310 is the highest good value, and 6311 is the lowest
failing value (assuming that the good and failing values are each
contiguous).

Here is the testing code:

from mpmath import *

# compute pi using n digits of accuracy, return True if okay, False if buggy
def okay(precision):
mp.dps = precision
# a sequence is successive arithmetic means;
# g sequence is successive geometric means
a, b = 1, sqrt(0.5)
n = 1
sum = 0
factor = 2
while n < mp.dps-10:
n = -log(a*a-b*b)/log(10)
a,b = (a+b)/2, sqrt(a*b)
factor *= 2
sum += factor * (a*a-b*b)
M = (a+b)/2
estimate = 4*M*M/(1-sum)
digits = -log(abs(estimate-pi))/log(10)
if digits > precision + 100: # bogusly high precision
print "Fail: %5d" % precision
return False
print "Okay: %5d" % precision
return True

n = 100
while okay(n):
n *= 2
if okay(test):
good = test
else:

If I reinstall python-gmpy and rerun the above script, I get:

Okay:   100
Okay:   200
Okay:   400
Okay:   800
Okay:  1600
Fail:  3200
Fail:  2400
Fail:  2000
Okay:  1800
Okay:  1900
Fail:  1950
Okay:  1925
Okay:  1937
Okay:  1943
Fail:  1946
Okay:  1944
Okay:  1945

So 1946 is now the lowest failing value.

```