[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
good, bad = n/2, n
while bad-good > 1:
test = (good+bad)//2
if okay(test):
good = test
else:
bad = test
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.
More information about the Python-modules-team
mailing list