python：浮点数的最大公约数（gcd），最好在 numpy 中答案

【问题标题】：python: Greatest common divisor (gcd) for floats, preferably in numpypython：浮点数的最大公约数（gcd），最好在 numpy 中
【发布时间】：2017-07-26 09:59:59
【问题描述】：

我正在寻找一种有效的方法来使用 python 确定两个浮点数的最大公约数。该例程应具有以下布局

gcd(a, b, rtol=1e-05, atol=1e-08)
"""
Returns the greatest common divisor of a and b

Parameters
----------
a,b : float
    two floats for gcd
rtol, atol : float, optional
    relative and absolute tolerance

Returns
-------
gcd : float
    Greatest common divisor such that for x in [a,b]:
    np.mod(x,gcd) < rtol*x + atol 

.. _PEP 484:
    https://www.python.org/dev/peps/pep-0484/

"""

例子：有理数和无理数的gcd

gcd(1., np.pi, rtol=0, atol=1e-5) 应该返回（大致）1e-5，因为

In [1]: np.mod(np.pi,1e-5)
Out[1]: 2.6535897928590063e-06

In [2]: np.mod(1.,1e-5)
Out[2]: 9.9999999999181978e-06

我更喜欢使用库实现，而不是自己编写。 fractions.gcd 函数在这里似乎不适合我，因为我不想使用分数，而且它（显然）没有公差参数。

【问题讨论】：

你能准确定义浮点数的最大公约数是什么意思吗？一对整数的 GCD 是众所周知且易于理解的事情。该定义很容易扩展到有理数，因此将浮点数视为有理数给出了对浮点数有效的 GCD 定义。但是鉴于您在这里包含宽容，我怀疑这就是您所追求的。一些示例可能会有所帮助。
比如1和pi的gcd是多少？
可能你必须自己实现它。你可以在 scipy/numpy 邮件列表中询问这个 SO 问题，我猜你可以在那里取得更大的成功。
我不确定这是否可以以一种漂亮的方式完成。多分辨率网格搜索可能是您的朋友。提供的示例：根据您的首选定义，gcd(1., np.pi, rtol=0, atol=1e-5) 大约不是 1e-5，因为大约 1e-4 已经击败了您。例如尝试 0.00010006。而且我很确定我可以找到一个适合您需要的更大的，如果我可以自己写整个网格搜索
请注意，fractions.gcd 无论如何都是纯 python 实现。如果您对atol 有一个很好的定义，那么将其合并到您自己的版本中应该很简单

标签： python numpy

【解决方案1】：

似乎您可以修改fractions.gcd 的代码以包含公差：

def float_gcd(a, b, rtol = 1e-05, atol = 1e-08):
    t = min(abs(a), abs(b))
    while abs(b) > rtol * t + atol:
        a, b = b, a % b
    return a

【讨论】：

我认为这不符合rtol 的要求，因为它不再与原版相关
也不要更改您的比较运算符 - OP 定义不正确，应该是 <=（使您的原始 > 正确）
我不清楚np.maximum 在这里是否正确。你能证明它的合理性吗？
另外，使用np.maximum over max, np.mod over %, np.absolute over abs 在这里是没有意义的，因为这在数组上无论如何都不起作用，因为循环
大部分情况下我只在numpy 工作，所以我使用我知道的功能（并且不必查找文档）。你的话我记住了。 max 也错了。

【解决方案2】：

以下代码可能有用。要调用的函数是 float_gdc(a, b)。

from math import gcd, log10, pow

def float_scale(x):
    max_digits = 14
    int_part = int(abs(x))
    magnitude = 1 if int_part == 0 else int(log10(int_part)) + 1
    if magnitude >= max_digits:
        return 0
    frac_part = abs(x) - int_part
    multiplier = 10 ** (max_digits - magnitude)
    frac_digits = multiplier + int(multiplier * frac_part + 0.5)
    while frac_digits % 10 == 0:
        frac_digits /= 10
    return int(log10(frac_digits))


def float_gcd(a, b):
    sc = float_scale(a)
    sc_b = float_scale(b)
    sc = sc_b if sc_b > sc else sc
    fac = pow(10, sc)

    a = int(round(a*fac))
    b = int(round(b*fac))

    return round(gcd(a, b)/fac, sc)

部分代码取自这里：Determine precision and scale of particular number in Python

【讨论】：