What algorithm does Python employ in fractions.gcd()?

Question

I'm using the fractions module in Python v3.1 to compute the greatest common divisor. I would like to know what algorithm is used. I'm guessing the Euclidean method, but would like to be sure. The docs (http://docs.python.org/py3k/library/fractions.html?highlight=fractions.gcd#fractions.gcd) don't help. Can anybody clue me in?

Stedman

Answer 1

According to the 3.1.2 source code online, here's gcd as defined in Python-3.1.2/Lib/fractions.py:

def gcd(a, b):
    """Calculate the Greatest Common Divisor of a and b.

    Unless b==0, the result will have the same sign as b (so that when
    b is divided by it, the result comes out positive).
    """
    while b:
        a, b = b, a%b
    return a

So yes, it's the Euclidean algorithm, written in pure Python.

Brindled

If you're using IPython, you can see the source code immediately by typing gcd?? – Delozier

It's actually: import fractions, then: fractions.gcd?? in IPython. – Semivitreous

Answer 2

From fractions python

"Deprecated since version 3.5: Use math.gcd() instead."

I was looking for the algorithm as well. I hope it helped.

Recurve

Answer 3

Since Python 3.5, the GCD code has been moved to math.gcd. Since Python 3.9, math.gcd takes an arbitrary number of arguments.

The actual GCD code is now implemented in C (for CPython), making it significantly faster than the original pure Python implementation.

Boilerplate:

static PyObject *
math_gcd(PyObject *module, PyObject * const *args, Py_ssize_t nargs)
{
    PyObject *res, *x;
    Py_ssize_t i;

    if (nargs == 0) {
        return PyLong_FromLong(0);
    }
    res = PyNumber_Index(args[0]);
    if (res == NULL) {
        return NULL;
    }
    if (nargs == 1) {
        Py_SETREF(res, PyNumber_Absolute(res));
        return res;
    }

    PyObject *one = _PyLong_GetOne();  // borrowed ref
    for (i = 1; i < nargs; i++) {
        x = _PyNumber_Index(args[i]);
        if (x == NULL) {
            Py_DECREF(res);
            return NULL;
        }
        if (res == one) {
            /* Fast path: just check arguments.
               It is okay to use identity comparison here. */
            Py_DECREF(x);
            continue;
        }
        Py_SETREF(res, _PyLong_GCD(res, x));
        Py_DECREF(x);
        if (res == NULL) {
            return NULL;
        }
    }
    return res;
}

The actual computation uses Lehmer's GCD algorithm:

PyObject *
_PyLong_GCD(PyObject *aarg, PyObject *barg)
{
    PyLongObject *a, *b, *c = NULL, *d = NULL, *r;
    stwodigits x, y, q, s, t, c_carry, d_carry;
    stwodigits A, B, C, D, T;
    int nbits, k;
    Py_ssize_t size_a, size_b, alloc_a, alloc_b;
    digit *a_digit, *b_digit, *c_digit, *d_digit, *a_end, *b_end;

    a = (PyLongObject *)aarg;
    b = (PyLongObject *)barg;
    size_a = Py_SIZE(a);
    size_b = Py_SIZE(b);
    if (-2 <= size_a && size_a <= 2 && -2 <= size_b && size_b <= 2) {
        Py_INCREF(a);
        Py_INCREF(b);
        goto simple;
    }

    /* Initial reduction: make sure that 0 <= b <= a. */
    a = (PyLongObject *)long_abs(a);
    if (a == NULL)
        return NULL;
    b = (PyLongObject *)long_abs(b);
    if (b == NULL) {
        Py_DECREF(a);
        return NULL;
    }
    if (long_compare(a, b) < 0) {
        r = a;
        a = b;
        b = r;
    }
    /* We now own references to a and b */

    alloc_a = Py_SIZE(a);
    alloc_b = Py_SIZE(b);
    /* reduce until a fits into 2 digits */
    while ((size_a = Py_SIZE(a)) > 2) {
        nbits = bit_length_digit(a->ob_digit[size_a-1]);
        /* extract top 2*PyLong_SHIFT bits of a into x, along with
           corresponding bits of b into y */
        size_b = Py_SIZE(b);
        assert(size_b <= size_a);
        if (size_b == 0) {
            if (size_a < alloc_a) {
                r = (PyLongObject *)_PyLong_Copy(a);
                Py_DECREF(a);
            }
            else
                r = a;
            Py_DECREF(b);
            Py_XDECREF(c);
            Py_XDECREF(d);
            return (PyObject *)r;
        }
        x = (((twodigits)a->ob_digit[size_a-1] << (2*PyLong_SHIFT-nbits)) |
             ((twodigits)a->ob_digit[size_a-2] << (PyLong_SHIFT-nbits)) |
             (a->ob_digit[size_a-3] >> nbits));

        y = ((size_b >= size_a - 2 ? b->ob_digit[size_a-3] >> nbits : 0) |
             (size_b >= size_a - 1 ? (twodigits)b->ob_digit[size_a-2] << (PyLong_SHIFT-nbits) : 0) |
             (size_b >= size_a ? (twodigits)b->ob_digit[size_a-1] << (2*PyLong_SHIFT-nbits) : 0));

        /* inner loop of Lehmer's algorithm; A, B, C, D never grow
           larger than PyLong_MASK during the algorithm. */
        A = 1; B = 0; C = 0; D = 1;
        for (k=0;; k++) {
            if (y-C == 0)
                break;
            q = (x+(A-1))/(y-C);
            s = B+q*D;
            t = x-q*y;
            if (s > t)
                break;
            x = y; y = t;
            t = A+q*C; A = D; B = C; C = s; D = t;
        }

        if (k == 0) {
            /* no progress; do a Euclidean step */
            if (l_mod(a, b, &r) < 0)
                goto error;
            Py_SETREF(a, b);
            b = r;
            alloc_a = alloc_b;
            alloc_b = Py_SIZE(b);
            continue;
        }

        /*
          a, b = A*b-B*a, D*a-C*b if k is odd
          a, b = A*a-B*b, D*b-C*a if k is even
        */
        if (k&1) {
            T = -A; A = -B; B = T;
            T = -C; C = -D; D = T;
        }
        if (c != NULL) {
            Py_SET_SIZE(c, size_a);
        }
        else if (Py_REFCNT(a) == 1) {
            c = (PyLongObject*)Py_NewRef(a);
        }
        else {
            alloc_a = size_a;
            c = _PyLong_New(size_a);
            if (c == NULL)
                goto error;
        }

        if (d != NULL) {
            Py_SET_SIZE(d, size_a);
        }
        else if (Py_REFCNT(b) == 1 && size_a <= alloc_b) {
            d = (PyLongObject*)Py_NewRef(b);
            Py_SET_SIZE(d, size_a);
        }
        else {
            alloc_b = size_a;
            d = _PyLong_New(size_a);
            if (d == NULL)
                goto error;
        }
        a_end = a->ob_digit + size_a;
        b_end = b->ob_digit + size_b;

        /* compute new a and new b in parallel */
        a_digit = a->ob_digit;
        b_digit = b->ob_digit;
        c_digit = c->ob_digit;
        d_digit = d->ob_digit;
        c_carry = 0;
        d_carry = 0;
        while (b_digit < b_end) {
            c_carry += (A * *a_digit) - (B * *b_digit);
            d_carry += (D * *b_digit++) - (C * *a_digit++);
            *c_digit++ = (digit)(c_carry & PyLong_MASK);
            *d_digit++ = (digit)(d_carry & PyLong_MASK);
            c_carry >>= PyLong_SHIFT;
            d_carry >>= PyLong_SHIFT;
        }
        while (a_digit < a_end) {
            c_carry += A * *a_digit;
            d_carry -= C * *a_digit++;
            *c_digit++ = (digit)(c_carry & PyLong_MASK);
            *d_digit++ = (digit)(d_carry & PyLong_MASK);
            c_carry >>= PyLong_SHIFT;
            d_carry >>= PyLong_SHIFT;
        }
        assert(c_carry == 0);
        assert(d_carry == 0);

        Py_INCREF(c);
        Py_INCREF(d);
        Py_DECREF(a);
        Py_DECREF(b);
        a = long_normalize(c);
        b = long_normalize(d);
    }
    Py_XDECREF(c);
    Py_XDECREF(d);

simple:
    assert(Py_REFCNT(a) > 0);
    assert(Py_REFCNT(b) > 0);
/* Issue #24999: use two shifts instead of ">> 2*PyLong_SHIFT" to avoid
   undefined behaviour when LONG_MAX type is smaller than 60 bits */
#if LONG_MAX >> PyLong_SHIFT >> PyLong_SHIFT
    /* a fits into a long, so b must too */
    x = PyLong_AsLong((PyObject *)a);
    y = PyLong_AsLong((PyObject *)b);
#elif LLONG_MAX >> PyLong_SHIFT >> PyLong_SHIFT
    x = PyLong_AsLongLong((PyObject *)a);
    y = PyLong_AsLongLong((PyObject *)b);
#else
# error "_PyLong_GCD"
#endif
    x = Py_ABS(x);
    y = Py_ABS(y);
    Py_DECREF(a);
    Py_DECREF(b);

    /* usual Euclidean algorithm for longs */
    while (y != 0) {
        t = y;
        y = x % y;
        x = t;
    }
#if LONG_MAX >> PyLong_SHIFT >> PyLong_SHIFT
    return PyLong_FromLong(x);
#elif LLONG_MAX >> PyLong_SHIFT >> PyLong_SHIFT
    return PyLong_FromLongLong(x);
#else
# error "_PyLong_GCD"
#endif

error:
    Py_DECREF(a);
    Py_DECREF(b);
    Py_XDECREF(c);
    Py_XDECREF(d);
    return NULL;
}

Methylnaphthalene