python：浮点的最大公约数（gcd），最好是numpy

Question

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

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，如下所示

^{pr2}$

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

Answer 1

似乎您可以修改^{}的代码以包含公差：

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

Answer 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