6

I am looking for an efficient way to determine the greatest common divisor of two floats with python. The routine should have the following layout

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/

"""

Example: gcd of rational and irrational number

The gcd(1., np.pi, rtol=0, atol=1e-5) should return (roughly) 1e-5, as

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

I would prefer to use a library implementation and not to write it myself. The fractions.gcd function does not seem appropriate to me here, as I do not want to work with fractions and it (obviously) does not have the tolerance parameters.

Markus Dutschke
  • 9,341
  • 4
  • 63
  • 58
  • 6
    Can you define exactly what you mean by the greatest common divisor of floats? The GCD of a pair of integers is a well-known and well-understood thing. The definition extends easily to rational numbers, so regarding floats as rationals gives a definition of GCD that's valid for floats. But given your inclusion of tolerance here, I doubt that's what you're after. Some examples might help. – Mark Dickinson Jul 26 '17 at 10:04
  • 1
    For example, what is the gcd of 1 and pi ? – Gribouillis Jul 26 '17 at 10:05
  • Probably you'll have to implement it yourself. You can ask in the scipy/numpy mailing list making reference to this SO question, my guess is that you can have more success there. – Ignacio Vergara Kausel Jul 26 '17 at 10:31
  • 1
    I am not sure this can be done in a beautiful way. Probably multi-resolution grid search is your friend here. The example provided: according to your preferred definition, the gcd(1., np.pi, rtol=0, atol=1e-5) is not roughly 1e-5, as roughly 1e-4 already beats you to the punch. Try 0.00010006 for instance. And I am quite sure I can find a bigger one suiting your needs if I could be bothered to write up the whole grid search myself – Uvar Jul 26 '17 at 10:44
  • Note that [`fractions.gcd`](https://github.com/python/cpython/blob/2.7/Lib/fractions.py#L18-L26) is a pure-python implementation anyway. If you have a good definition of `atol`, then it should be straightforward to incorporate it into your own version – Eric Jul 26 '17 at 11:05
  • Your definition should be `<=`, not `<`, else `rtol=atol=0` (the default for integers) is ill-posed – Eric Jul 26 '17 at 11:32

2 Answers2

8

Seems like you could just modify the code of fractions.gcd to include the tolerances:

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
Daniel F
  • 13,620
  • 2
  • 29
  • 55
  • I don't think this meets the requirements for `rtol`, since it's no longer relative to the original – Eric Jul 26 '17 at 11:27
  • Don't change your comparison operator either - the OPs definition is incorrect, and should be `<=` (making your original `>` correct) – Eric Jul 26 '17 at 11:34
  • It's not clear to me that `np.maximum` is correct here. Can you justify it? – Eric Jul 26 '17 at 11:41
  • Also, using `np.maximum` over `max`, `np.mod` over `%`, `np.absolute` over `abs` is pointless here, because this won't work on arrays anyway due to the loop – Eric Jul 26 '17 at 11:42
  • Mostly I only work in `numpy` so I use the functions I know (and don't have to look up the documentation for). I'll take your word for it. Also `max` was wrong. – Daniel F Jul 26 '17 at 11:44
  • `min` is intuitively better, but I'm unconvinced it's correct - might be overly strict – Eric Jul 26 '17 at 11:55
1

The following code may be useful. The function to call is 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)

A part of the code is taken from here: Determine precision and scale of particular number in Python

Lorenzo Celli
  • 11
  • 1