Issue #22486: Added the math.gcd() function. The fractions.gcd() function now is

deprecated.  Based on patch by Mark Dickinson.

Doc/library/fractions.rst

@@ -172,6 +172,9 @@ another rational number, or from a string.
    sign as *b* if *b* is nonzero; otherwise it takes the sign of *a*.  ``gcd(0,
    0)`` returns ``0``.
 
+   .. deprecated:: 3.5
+      Use :func:`math.gcd` instead.
+
 
 .. seealso::
 

Doc/library/math.rst

@@ -100,6 +100,14 @@ Number-theoretic and representation functions
    <http://code.activestate.com/recipes/393090/>`_\.
 
 
+.. function:: gcd(a, b)
+
+   Return the greatest common divisor of the integers *a* and *b*.  If either
+   *a* or *b* is nonzero, then the value of ``gcd(a, b)`` is the largest
+   positive integer that divides both *a* and *b*.  ``gcd(0, 0)`` returns
+   ``0``.
+
+
 .. function:: isfinite(x)
 
    Return ``True`` if *x* is neither an infinity nor a NaN, and

Include/longobject.h

@@ -198,6 +198,9 @@ PyAPI_FUNC(int) _PyLong_FormatAdvancedWriter(
 PyAPI_FUNC(unsigned long) PyOS_strtoul(const char *, char **, int);
 PyAPI_FUNC(long) PyOS_strtol(const char *, char **, int);
 
+/* For use by the gcd function in mathmodule.c */
+PyAPI_FUNC(PyObject *) _PyLong_GCD(PyObject *, PyObject *);
+
 #ifdef __cplusplus
 }
 #endif

Lib/fractions.py

@@ -20,6 +20,17 @@ def gcd(a, 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).
     """
+    import warnings
+    warnings.warn('fractions.gcd() is deprecated. Use math.gcd() instead.',
+                  DeprecationWarning, 2)
+    if type(a) is int is type(b):
+        if (b or a) < 0:
+            return -math.gcd(a, b)
+        return math.gcd(a, b)
+    return _gcd(a, b)
+
+def _gcd(a, b):
+    # Supports non-integers for backward compatibility.
     while b:
         a, b = b, a%b
     return a
@@ -159,7 +170,7 @@ class Fraction(numbers.Rational):
                                 "or a Rational instance")
 
         elif type(numerator) is int is type(denominator):
-            pass  # *very* normal case
+            pass # *very* normal case
 
         elif (isinstance(numerator, numbers.Rational) and
             isinstance(denominator, numbers.Rational)):
@@ -174,7 +185,13 @@ class Fraction(numbers.Rational):
         if denominator == 0:
             raise ZeroDivisionError('Fraction(%s, 0)' % numerator)
         if _normalize:
-            g = gcd(numerator, denominator)
+            if type(numerator) is int is type(denominator):
+                # *very* normal case
+                g = math.gcd(numerator, denominator)
+                if denominator < 0:
+                    g = -g
+            else:
+                g = _gcd(numerator, denominator)
             numerator //= g
             denominator //= g
         self._numerator = numerator

Lib/test/test_fractions.py

@@ -8,6 +8,7 @@ import operator
 import fractions
 import sys
 import unittest
+import warnings
 from copy import copy, deepcopy
 from pickle import dumps, loads
 F = fractions.Fraction
@@ -49,7 +50,7 @@ class DummyRational(object):
     """Test comparison of Fraction with a naive rational implementation."""
 
     def __init__(self, num, den):
-        g = gcd(num, den)
+        g = math.gcd(num, den)
         self.num = num // g
         self.den = den // g
 
@@ -83,16 +84,26 @@ class DummyFraction(fractions.Fraction):
 class GcdTest(unittest.TestCase):
 
     def testMisc(self):
-        self.assertEqual(0, gcd(0, 0))
-        self.assertEqual(1, gcd(1, 0))
-        self.assertEqual(-1, gcd(-1, 0))
-        self.assertEqual(1, gcd(0, 1))
-        self.assertEqual(-1, gcd(0, -1))
-        self.assertEqual(1, gcd(7, 1))
-        self.assertEqual(-1, gcd(7, -1))
-        self.assertEqual(1, gcd(-23, 15))
-        self.assertEqual(12, gcd(120, 84))
-        self.assertEqual(-12, gcd(84, -120))
+        # fractions.gcd() is deprecated
+        with self.assertWarnsRegex(DeprecationWarning, r'fractions\.gcd'):
+            gcd(1, 1)
+        with warnings.catch_warnings():
+            warnings.filterwarnings('ignore', r'fractions\.gcd',
+                                    DeprecationWarning)
+            self.assertEqual(0, gcd(0, 0))
+            self.assertEqual(1, gcd(1, 0))
+            self.assertEqual(-1, gcd(-1, 0))
+            self.assertEqual(1, gcd(0, 1))
+            self.assertEqual(-1, gcd(0, -1))
+            self.assertEqual(1, gcd(7, 1))
+            self.assertEqual(-1, gcd(7, -1))
+            self.assertEqual(1, gcd(-23, 15))
+            self.assertEqual(12, gcd(120, 84))
+            self.assertEqual(-12, gcd(84, -120))
+            self.assertEqual(gcd(120.0, 84), 12.0)
+            self.assertEqual(gcd(120, 84.0), 12.0)
+            self.assertEqual(gcd(F(120), F(84)), F(12))
+            self.assertEqual(gcd(F(120, 77), F(84, 55)), F(12, 385))
 
 
 def _components(r):

Lib/test/test_math.py

@@ -175,6 +175,14 @@ def parse_testfile(fname):
                    flags
                   )
 
+# Class providing an __index__ method.
+class MyIndexable(object):
+    def __init__(self, value):
+        self.value = value
+
+    def __index__(self):
+        return self.value
+
 class MathTests(unittest.TestCase):
 
     def ftest(self, name, value, expected):
@@ -595,6 +603,49 @@ class MathTests(unittest.TestCase):
             s = msum(vals)
             self.assertEqual(msum(vals), math.fsum(vals))
 
+    def testGcd(self):
+        gcd = math.gcd
+        self.assertEqual(gcd(0, 0), 0)
+        self.assertEqual(gcd(1, 0), 1)
+        self.assertEqual(gcd(-1, 0), 1)
+        self.assertEqual(gcd(0, 1), 1)
+        self.assertEqual(gcd(0, -1), 1)
+        self.assertEqual(gcd(7, 1), 1)
+        self.assertEqual(gcd(7, -1), 1)
+        self.assertEqual(gcd(-23, 15), 1)
+        self.assertEqual(gcd(120, 84), 12)
+        self.assertEqual(gcd(84, -120), 12)
+        self.assertEqual(gcd(1216342683557601535506311712,
+                             436522681849110124616458784), 32)
+        c = 652560
+        x = 434610456570399902378880679233098819019853229470286994367836600566
+        y = 1064502245825115327754847244914921553977
+        a = x * c
+        b = y * c
+        self.assertEqual(gcd(a, b), c)
+        self.assertEqual(gcd(b, a), c)
+        self.assertEqual(gcd(-a, b), c)
+        self.assertEqual(gcd(b, -a), c)
+        self.assertEqual(gcd(a, -b), c)
+        self.assertEqual(gcd(-b, a), c)
+        self.assertEqual(gcd(-a, -b), c)
+        self.assertEqual(gcd(-b, -a), c)
+        c = 576559230871654959816130551884856912003141446781646602790216406874
+        a = x * c
+        b = y * c
+        self.assertEqual(gcd(a, b), c)
+        self.assertEqual(gcd(b, a), c)
+        self.assertEqual(gcd(-a, b), c)
+        self.assertEqual(gcd(b, -a), c)
+        self.assertEqual(gcd(a, -b), c)
+        self.assertEqual(gcd(-b, a), c)
+        self.assertEqual(gcd(-a, -b), c)
+        self.assertEqual(gcd(-b, -a), c)
+
+        self.assertRaises(TypeError, gcd, 120.0, 84)
+        self.assertRaises(TypeError, gcd, 120, 84.0)
+        self.assertEqual(gcd(MyIndexable(120), MyIndexable(84)), 12)
+
     def testHypot(self):
         self.assertRaises(TypeError, math.hypot)
         self.ftest('hypot(0,0)', math.hypot(0,0), 0)

Misc/NEWS

@@ -42,6 +42,9 @@ Core and Builtins
 Library
 -------
 
+- Issue #22486: Added the math.gcd() function.  The fractions.gcd() function now is
+  deprecated.  Based on patch by Mark Dickinson.
+
 - Issue #22681: Added support for the koi8_t encoding.
 
 - Issue #22682: Added support for the kz1048 encoding.

Modules/mathmodule.c

@@ -685,6 +685,33 @@ m_log10(double x)
 }
 
 
+static PyObject *
+math_gcd(PyObject *self, PyObject *args)
+{
+    PyObject *a, *b, *g;
+
+    if (!PyArg_ParseTuple(args, "OO:gcd", &a, &b))
+        return NULL;
+
+    a = PyNumber_Index(a);
+    if (a == NULL)
+        return NULL;
+    b = PyNumber_Index(b);
+    if (b == NULL) {
+        Py_DECREF(a);
+        return NULL;
+    }
+    g = _PyLong_GCD(a, b);
+    Py_DECREF(a);
+    Py_DECREF(b);
+    return g;
+}
+
+PyDoc_STRVAR(math_gcd_doc,
+"gcd(x, y) -> int\n\
+greatest common divisor of x and y");
+
+
 /* Call is_error when errno != 0, and where x is the result libm
  * returned.  is_error will usually set up an exception and return
  * true (1), but may return false (0) without setting up an exception.
@@ -1987,6 +2014,7 @@ static PyMethodDef math_methods[] = {
     {"frexp",           math_frexp,     METH_O,         math_frexp_doc},
     {"fsum",            math_fsum,      METH_O,         math_fsum_doc},
     {"gamma",           math_gamma,     METH_O,         math_gamma_doc},
+    {"gcd",             math_gcd,       METH_VARARGS,   math_gcd_doc},
     {"hypot",           math_hypot,     METH_VARARGS,   math_hypot_doc},
     {"isfinite",        math_isfinite,  METH_O,         math_isfinite_doc},
     {"isinf",           math_isinf,     METH_O,         math_isinf_doc},

Objects/longobject.c

@@ -4327,6 +4327,211 @@ long_long(PyObject *v)
     return v;
 }
 
+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 = bits_in_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_divmod(a, b, NULL, &r) < 0)
+                goto error;
+            Py_DECREF(a);
+            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_SIZE(c) = size_a;
+        else if (Py_REFCNT(a) == 1) {
+            Py_INCREF(a);
+            c = a;
+        }
+        else {
+            alloc_a = size_a;
+            c = _PyLong_New(size_a);
+            if (c == NULL)
+                goto error;
+        }
+
+        if (d != NULL)
+            Py_SIZE(d) = size_a;
+        else if (Py_REFCNT(b) == 1 && size_a <= alloc_b) {
+            Py_INCREF(b);
+            d = b;
+            Py_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);
+#if LONG_MAX >> 2*PyLong_SHIFT
+    /* a fits into a long, so b must too */
+    x = PyLong_AsLong((PyObject *)a);
+    y = PyLong_AsLong((PyObject *)b);
+#elif defined(PY_LONG_LONG) && PY_LLONG_MAX >> 2*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 >> 2*PyLong_SHIFT
+    return PyLong_FromLong(x);
+#elif defined(PY_LONG_LONG) && PY_LLONG_MAX >> 2*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;
+}
+
 static PyObject *
 long_float(PyObject *v)
 {