[issue22486] Add math.gcd()
Serhiy Storchaka
report at bugs.python.org
Thu Sep 25 14:53:53 CEST 2014
Serhiy Storchaka added the comment:
Well, here is a patch which keeps the same weird behavior of fractions.gcd().
----------
Added file: http://bugs.python.org/file36723/lehmer_gcd_5.patch
_______________________________________
Python tracker <report at bugs.python.org>
<http://bugs.python.org/issue22486>
_______________________________________
-------------- next part --------------
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
<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
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)
{
More information about the Python-bugs-list
mailing list