Is that a bug of GMP(gmpy)?

[Xiao-Qin Xia]

Programs:

import gmpy
gmpy.set_minprec(128)
def Agm_gmpy(a, b):
        """
        The arithmetic-geometric mean (often abbreviated AGM) M(a,b) of a and b.
        see http://mathworld.wolfram.com/Arithmetic-GeometricMean.html
        """
        if a < b: a, b = gmpy.mpf(b), gmpy.mpf(a)
        else: a, b = gmpy.mpf(a), gmpy.mpf(b)
        f2 = gmpy.mpf(2.0)
        d, d_old = a-b, gmpy.mpf(0.0)
        while d: 
                a, b = (a+b)/f2, gmpy.fsqrt(a*b)
                d = a - b
                print "a=",a
                print "b=",b
                print "d=",d
                if d == d_old:
                        break
                else: 
                        d_old = d
        return a

print Agm_gmpy(4.022e-5, 1.25)


Result:

d= 3.421138828918010427059886677953896804883e-49
a= 0.1673825781925968686521706972176572335454
b= 0.1673825781925968686521706972176572335454
d= 6.842277657836020854119773355907793609767e-49
a= 0.1673825781925968686521706972176572335454
b= 0.1673825781925968686521706972176572335454
d= 3.421138828918010427059886677953896804883e-49
a= 0.1673825781925968686521706972176572335454
b= 0.1673825781925968686521706972176572335454
d= 6.842277657836020854119773355907793609767e-49


Question:
arithmetic-geometric mean has a significant character:
if 
a_n >= b_n 
and a_(n+1) = (a_n + b_n)/2, 
    b_(n+a) = sqrt(a_n * b_n)
then
a_(n+1) - b_(n+1) <= (a_n - b_n)/2

that mean the difference between a_(n+1) and b_(n+1) reduce very quickly, 
i.e. a_(n+1) approach b_(n+1) very quickly.
however the result of the program show that d, the difference 
between a_(n+1) anb b_(n+1), can increase and reduce, falling in a dead 
cycle.

Is it a bug of gmpy (the add, minus, divide, multiple and sqrt function),
or the code above need to be improved?

Many thanks,
Xiao-Qin