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

[Sanjoy Mahajan]

> 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.