Problems with modulo operation in Python for computing gcd

Question

I discovered some strange behavior of the python modulo operation. The command

a = 1.0 % 0.1

yields

a == 0.09999999999999995

which is quite a big error that causes me problems when calculating the greatest common divisor of two float numbers. I suppose the error is related to the non-representatilty of 0.1 and 1.0 in binary (http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems). Does a workaround for this problem exist or can someone point me at a reliably function to compute the gcd of two numbers? The code that I have been using so far for computing the gcd is

def gcd(a, b):
    a, b = np.broadcast_arrays(a, b)
    a = a.copy()
    b = b.copy()
    pos = np.nonzero(b)[0]
    while len(pos) > 0:
        b2 = b[pos]
        a[pos], b[pos] = b2, a[pos] % b2
        pos = pos[b[pos]!=0]
    return a

I have the same problem when I use the fractions module.

Answer 1

You can use the decimal module.

from fractions import gcd
from decimal import Decimal
a,b  = 1.0,0.1
print(gcd(Decimal(str(a)),Decimal(str(b))))
0.1
a,b  = 22.8,9.3
print(gcd(Decimal(str(a)),Decimal(str(b))))
0.3

Answer 2

You can try to use Decimal in python to not be dependent on float "errors".

>>> from decimal import Decimal as D
>>> a = D(1.0) % D(0.1)
>>> a
Decimal('0.09999999999999995003996389187')

Answer 3

A very simple hack that will work with normal numpy float arrays would be to multiply a and b by some large factor before taking the modulo:

def gcd(a, b, fac=1E12):
    a, b = np.broadcast_arrays(a, b)
    a = a.copy() * fac
    b = b.copy() * fac
    pos = np.nonzero(b)[0]
    while len(pos) > 0:
        b2 = b[pos]
        a[pos], b[pos] = b2, a[pos] % b2
        pos = pos[b[pos]!=0]
    return a / fac

print(gcd(np.r_[1.0, 22.8], np.r_[0.1, 9.3]))
# [ 0.1  0.3]

One advantage of this approach over using the decimal module is that you can still exploit the speed of numpy's vectorized array operations. Since decimal.Decimal can only represent scalars, you would otherwise have to loop over the elements in your input arrays and process them all individually.