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 Sat Sep 27 11:09:15 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 Sat Sep 27 11:09:15 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 Sat Sep 27 11:09:15 2014 +0300 @@ -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 Sat Sep 27 11:09:15 2014 +0300 @@ -595,6 +595,45 @@ 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) + 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 Sat Sep 27 11:09:15 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 Sat Sep 27 11:09:15 2014 +0300 @@ -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) {