diff -r e9d4288c32de Doc/library/math.rst --- a/Doc/library/math.rst Wed Sep 24 13:29:27 2014 +0300 +++ b/Doc/library/math.rst Thu Sep 25 15:51:26 2014 +0300 @@ -100,6 +100,14 @@ Number-theoretic and representation func `_\. +.. 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 diff -r e9d4288c32de Include/longobject.h --- a/Include/longobject.h Wed Sep 24 13:29:27 2014 +0300 +++ b/Include/longobject.h Thu Sep 25 15:51:26 2014 +0300 @@ -198,6 +198,9 @@ PyAPI_FUNC(int) _PyLong_FormatAdvancedWr 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 diff -r e9d4288c32de Lib/fractions.py --- a/Lib/fractions.py Wed Sep 24 13:29:27 2014 +0300 +++ b/Lib/fractions.py Thu Sep 25 15:51:26 2014 +0300 @@ -20,9 +20,9 @@ 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). """ - while b: - a, b = b, a%b - return a + if (b or a) < 0: + return -math.gcd(a, b) + return math.gcd(a, b) # Constants related to the hash implementation; hash(x) is based # on the reduction of x modulo the prime _PyHASH_MODULUS. @@ -174,9 +174,12 @@ class Fraction(numbers.Rational): if denominator == 0: raise ZeroDivisionError('Fraction(%s, 0)' % numerator) if _normalize: - g = gcd(numerator, denominator) + g = math.gcd(numerator, denominator) numerator //= g denominator //= g + if denominator < 0: + numerator = -numerator + denominator = -denominator self._numerator = numerator self._denominator = denominator return self diff -r e9d4288c32de Lib/test/test_math.py --- a/Lib/test/test_math.py Wed Sep 24 13:29:27 2014 +0300 +++ b/Lib/test/test_math.py Thu Sep 25 15:51:26 2014 +0300 @@ -595,6 +595,24 @@ class MathTests(unittest.TestCase): s = msum(vals) self.assertEqual(msum(vals), math.fsum(vals)) + def testGcd(self): + 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(190738355881570558882299312308821696901058000, + 76478560266291874249006856460326062498333440), + 652560) + self.assertEqual(gcd(83763289342793979220453055528167457860243376086879213707165435635135627040075, + 33585776402955145260404154387726204875807368546078094789530226423049489520976), + 286573572687563623189610484223662247799) + def testHypot(self): self.assertRaises(TypeError, math.hypot) self.ftest('hypot(0,0)', math.hypot(0,0), 0) diff -r e9d4288c32de Modules/mathmodule.c --- a/Modules/mathmodule.c Wed Sep 24 13:29:27 2014 +0300 +++ b/Modules/mathmodule.c Thu Sep 25 15:51:26 2014 +0300 @@ -656,6 +656,22 @@ m_log10(double x) } +static PyObject * +math_gcd(PyObject *self, PyObject *args) +{ + PyObject *a, *b; + + if (!PyArg_ParseTuple(args, "O!O!:gcd", &PyLong_Type, &a, &PyLong_Type, &b)) + return NULL; + + return _PyLong_GCD(a, b); +} + +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. @@ -1958,6 +1974,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}, diff -r e9d4288c32de Objects/longobject.c --- a/Objects/longobject.c Wed Sep 24 13:29:27 2014 +0300 +++ b/Objects/longobject.c Thu Sep 25 15:51:26 2014 +0300 @@ -4327,6 +4327,188 @@ long_long(PyObject *v) return v; } +static PyLongObject * +long_gcd(PyLongObject *a, PyLongObject *b) +{ + PyLongObject *c, *d; + stwodigits x, y, q, s, t, c_carry, d_carry; + digit A, B, C, D; + int nbits, k; + Py_ssize_t size_a, size_b; + digit *a_digit, *b_digit, *c_digit, *d_digit, *a_end, *b_end; + + /* Initial reduction: make sure that 0 <= b <= a. */ + a = (PyLongObject *)long_abs(a); + b = (PyLongObject *)long_abs(b); + if (long_compare(a, b) < 0) { + d = a; + a = b; + b = d; + } + /* We now own references to a and 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); + x = ((a->ob_digit[size_a-1] << (2*PyLong_SHIFT-nbits)) | + (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 ? b->ob_digit[size_a-2] << (PyLong_SHIFT-nbits) : 0) | + (size_b >= size_a ? 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 = (digit)s; D = (digit)t; + } + + if (k == 0) { + /* no progress; do a Euclidean step */ + if (Py_SIZE(b) == 0) { + Py_DECREF(b); + return a; + } + if (l_divmod(a, b, NULL, &d) < 0) { + Py_DECREF(a); + Py_DECREF(b); + return NULL; + } + Py_DECREF(a); + a = b; + b = d; + 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 + */ + c = _PyLong_New(size_a); + if (c == NULL) { + Py_DECREF(a); + Py_DECREF(b); + return NULL; + } + + d = _PyLong_New(size_a); + if (d == NULL) { + Py_DECREF(a); + Py_DECREF(b); + Py_DECREF(c); + return NULL; + } + 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; + if (k&1) { + while (b_digit < b_end) { + c_carry += (A * *b_digit) - (B * *a_digit); + d_carry += (D * *a_digit++) - (C * *b_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 -= B * *a_digit; + d_carry += D * *a_digit++; + *c_digit++ = (digit)(c_carry & PyLong_MASK); + *d_digit++ = (digit)(d_carry & PyLong_MASK); + c_carry >>= PyLong_SHIFT; + d_carry >>= PyLong_SHIFT; + } + } + else { + 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_DECREF(a); + Py_DECREF(b); + a = long_normalize(c); + b = long_normalize(d); + } + + /* a fits into a long, so b must too */ + x = PyLong_AsLong((PyObject *)a); + y = PyLong_AsLong((PyObject *)b); + Py_DECREF(a); + Py_DECREF(b); + + /* usual Euclidean algorithm for longs */ + while (y != 0) { + t = y; + y = x % y; + x = t; + } + return (PyLongObject *)PyLong_FromLong(x); +} + +PyObject * +_PyLong_GCD(PyObject *a, PyObject *b) +{ + long x, y, t; + int overflow; + + x = PyLong_AsLongAndOverflow(a, &overflow); + if (!overflow && !(x == -1 && PyErr_Occurred()) && x >= -LONG_MAX) { + y = PyLong_AsLongAndOverflow(b, &overflow); + if (!overflow && !(y == -1 && PyErr_Occurred()) && y >= -LONG_MAX) { + /* Both a and b are small integers; + use the usual gcd algorithm. */ + if (x < 0) + x = -x; + if (y < 0) + y = -y; + while (y != 0) { + t = x % y; + x = y; + y = t; + } + return PyLong_FromLong(x); + } + } + + return (PyObject *)long_gcd((PyLongObject *)a, (PyLongObject *)b); +} + static PyObject * long_float(PyObject *v) {