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

[Sanjoy Mahajan]

Package: python-mpmath
Version: 0.11-2
Severity: normal

The following Python code uses mpmath to implement Salamin's algorithm
for computing pi (using arithmetic and geometric means).  The number
of digits of accuracy should roughly double with each iteration, and
it does so if mp.dps=1000 (the number of digits is printed at the end
of each swing through the loop).  

However, if mp.dps=4000, then the number of digits doubles for a few
iterations, then goes to infinity:

    3.0
    8.1
   18.7
   40.3
   83.6
    inf
    inf
    inf

Somehow the estimate became exactly equal to pi.  It seems strange to
me and I wonder if it's a bug in one of the math functions like sqrt
or whether it's a feature.  Or perhaps mpmath computes pi using the
same algorithm?


from mpmath import *

mp.dps = 4000

# 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)
    print "%7.1f" % digits


-- System Information:
Debian Release: squeeze/sid
  APT prefers unstable
  APT policy: (500, 'unstable'), (1, 'experimental')
Architecture: i386 (i686)

Kernel: Linux 2.6.29-2-686 (SMP w/2 CPU cores)
Locale: LANG=en_US.utf8, LC_CTYPE=en_US.utf8 (charmap=UTF-8)
Shell: /bin/sh linked to /bin/bash

Versions of packages python-mpmath depends on:
ii  python-support                1.0.3      automated rebuilding support for P

python-mpmath recommends no packages.

Versions of packages python-mpmath suggests:
ii  python-gmpy                   1.03.ds-1  Interfaces GMP to Python for fast,
ii  python-matplotlib             0.98.3-5   Python based plotting system in a 
ii  python-mpmath-doc             0.11-2     library for arbitrary-precision fl

-- no debconf information