From 48e47aaa28d6dfdae128142ffcbc4b0642422e90 Mon Sep 17 00:00:00 2001 From: Serhiy Storchaka Date: Wed, 13 May 2015 00:19:51 +0300 Subject: Issue #22486: Added the math.gcd() function. The fractions.gcd() function now is deprecated. Based on patch by Mark Dickinson. --- Doc/library/fractions.rst | 3 + Doc/library/math.rst | 8 ++ Include/longobject.h | 3 + Lib/fractions.py | 21 ++++- Lib/test/test_fractions.py | 33 +++++--- Lib/test/test_math.py | 51 +++++++++++ Misc/NEWS | 3 + Modules/mathmodule.c | 28 +++++++ Objects/longobject.c | 205 +++++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 342 insertions(+), 13 deletions(-) diff --git a/Doc/library/fractions.rst b/Doc/library/fractions.rst index c2c7401..e0f0682 100644 --- a/Doc/library/fractions.rst +++ b/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:: diff --git a/Doc/library/math.rst b/Doc/library/math.rst index eda0056..11389e6 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -100,6 +100,14 @@ Number-theoretic and representation functions `_\. +.. 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 --git a/Include/longobject.h b/Include/longobject.h index ff43309..aed59ce 100644 --- a/Include/longobject.h +++ b/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 diff --git a/Lib/fractions.py b/Lib/fractions.py index 5ddc84c..60b0728 100644 --- a/Lib/fractions.py +++ b/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 diff --git a/Lib/test/test_fractions.py b/Lib/test/test_fractions.py index 2dd528f..1699852 100644 --- a/Lib/test/test_fractions.py +++ b/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): diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 023dea9..fcd78d5 100644 --- a/Lib/test/test_math.py +++ b/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) diff --git a/Misc/NEWS b/Misc/NEWS index d714016..58a1daa 100644 --- a/Misc/NEWS +++ b/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. diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 153d152..a65de47 100644 --- a/Modules/mathmodule.c +++ b/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}, diff --git a/Objects/longobject.c b/Objects/longobject.c index 27bee50..40da0b1 100644 --- a/Objects/longobject.c +++ b/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) { -- cgit v0.12