Python 中用于计算 gcd 的模运算问题

Question

我发现 python 模运算有一些奇怪的行为。命令

a = 1.0 % 0.1

产量

a == 0.09999999999999995

这是一个相当大的错误，导致我在计算两个浮点数的最大公约数时遇到问题。我想这个错误与二进制中 0.1 和 1.0 的非代表性有关 (http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems)。是否存在解决此问题的方法，或者有人可以指出一个可靠的函数来计算两个数字的 gcd？到目前为止，我一直使用的计算 gcd 的代码是

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

我在使用分数模块时遇到了同样的问题。

Answer 1

您可以尝试在 python 中使用 Decimal 来不依赖于 float "errors"。

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

Answer 2

您可以使用 decimal 模块。

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 3

适用于普通 numpy 浮点数组的一个非常简单的 hack 是在取模之前将 a 和 b 乘以某个大的因子：

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]

与使用 decimal 模块相比，这种方法的一个优点是您仍然可以利用 numpy 的向量化数组操作的速度。由于 decimal.Decimal 只能表示标量，否则您将不得不遍历输入数组中的元素并单独处理它们。

Python 中用于计算 gcd 的模运算问题

Problems with modulo operation in Python for computing gcd

python

numpy