diff --git a/doc/release/1.15.0-notes.rst b/doc/release/1.15.0-notes.rst index 283c992ea489..be6d201a3926 100644 --- a/doc/release/1.15.0-notes.rst +++ b/doc/release/1.15.0-notes.rst @@ -10,6 +10,9 @@ Highlights New functions ============= +* `np.gcd` and `np.lcm`, to compute the greatest common divisor and least + common multiple. + Deprecations ============ @@ -40,6 +43,12 @@ C API changes New Features ============ +``np.gcd`` and ``np.lcm`` ufuncs added for integer and objects types +-------------------------------------------------------------------- +These compute the greatest common divisor, and lowest common multiple, +respectively. These work on all the numpy integer types, as well as the +builtin arbitrary-precision `Decimal` and `long` types. + Improvements ============ diff --git a/doc/source/reference/routines.math.rst b/doc/source/reference/routines.math.rst index 4c2f2800a780..821363987acc 100644 --- a/doc/source/reference/routines.math.rst +++ b/doc/source/reference/routines.math.rst @@ -101,6 +101,14 @@ Floating point routines nextafter spacing +Rational routines +----------------- +.. autosummary:: + :toctree: generated/ + + lcm + gcd + Arithmetic operations --------------------- .. autosummary:: diff --git a/doc/source/reference/ufuncs.rst b/doc/source/reference/ufuncs.rst index 38f2926f7ef8..3711f660f8e7 100644 --- a/doc/source/reference/ufuncs.rst +++ b/doc/source/reference/ufuncs.rst @@ -550,6 +550,8 @@ Math operations square cbrt reciprocal + gcd + lcm .. tip:: diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py index af058b4befef..9287be095d0b 100644 --- a/numpy/core/code_generators/generate_umath.py +++ b/numpy/core/code_generators/generate_umath.py @@ -875,6 +875,20 @@ def english_upper(s): TypeDescription('d', None, 'd', 'di'), TypeDescription('g', None, 'g', 'gi'), ], + ), +'gcd' : + Ufunc(2, 1, Zero, + docstrings.get('numpy.core.umath.gcd'), + "PyUFunc_SimpleBinaryOperationTypeResolver", + TD(ints), + TD('O', f='npy_ObjectGCD'), + ), +'lcm' : + Ufunc(2, 1, None, + docstrings.get('numpy.core.umath.lcm'), + "PyUFunc_SimpleBinaryOperationTypeResolver", + TD(ints), + TD('O', f='npy_ObjectLCM'), ) } diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py index 5626f50d8fde..504d7e6a977f 100644 --- a/numpy/core/code_generators/ufunc_docstrings.py +++ b/numpy/core/code_generators/ufunc_docstrings.py @@ -3679,3 +3679,63 @@ def add_newdoc(place, name, doc): array([ 0., 1., 2., 3., 4., 5.]) """) + +add_newdoc('numpy.core.umath', 'gcd', + """ + Returns the greatest common divisor of |x1| and |x2| + + Parameters + ---------- + x1, x2 : array_like, int + Arrays of values + + Returns + ------- + y : ndarray or scalar + The greatest common divisor of the absolute value of the inputs + + See Also + -------- + lcm : The lowest common multiple + + Examples + -------- + >>> np.gcd(12, 20) + 4 + >>> np.gcd.reduce([15, 25, 35]) + 5 + >>> np.gcd(np.arange(6), 20) + array([20, 1, 2, 1, 4, 5]) + + """) + +add_newdoc('numpy.core.umath', 'lcm', + """ + Returns the lowest common multiple of |x1| and |x2| + + Parameters + ---------- + x1, x2 : array_like, int + Arrays of values + + Returns + ------- + y : ndarray or scalar + The lowest common multiple of the absolute value of the inputs + + See Also + -------- + gcd : The greatest common divisor + + Examples + -------- + >>> np.lcm(12, 20) + 60 + >>> np.lcm.reduce([3, 12, 20]) + 60 + >>> np.lcm.reduce([40, 12, 20]) + 120 + >>> np.lcm(np.arange(6), 20) + array([ 0, 20, 20, 60, 20, 20]) + + """) diff --git a/numpy/core/src/npymath/npy_math_internal.h.src b/numpy/core/src/npymath/npy_math_internal.h.src index 093e51b2dcdf..f2e5229b0ae1 100644 --- a/numpy/core/src/npymath/npy_math_internal.h.src +++ b/numpy/core/src/npymath/npy_math_internal.h.src @@ -678,3 +678,41 @@ npy_divmod@c@(@type@ a, @type@ b, @type@ *modulus) #undef DEG2RAD /**end repeat**/ + +/**begin repeat + * + * #type = npy_uint, npy_ulong, npy_ulonglong# + * #c = u,ul,ull# + */ +NPY_INPLACE @type@ +npy_gcd@c@(@type@ a, @type@ b) +{ + @type@ c; + while (a != 0) { + c = a; + a = b%a; + b = c; + } + return b; +} + +NPY_INPLACE @type@ +npy_lcm@c@(@type@ a, @type@ b) +{ + @type@ gcd = npy_gcd@c@(a, b); + return gcd == 0 ? 0 : a / gcd * b; +} +/**end repeat**/ + +/**begin repeat + * + * #type = (npy_int, npy_long, npy_longlong)*2# + * #c = (,l,ll)*2# + * #func=gcd*3,lcm*3# + */ +NPY_INPLACE @type@ +npy_@func@@c@(@type@ a, @type@ b) +{ + return npy_@func@u@c@(a < 0 ? -a : a, b < 0 ? -b : b); +} +/**end repeat**/ diff --git a/numpy/core/src/umath/funcs.inc.src b/numpy/core/src/umath/funcs.inc.src index 5613c30ee60e..da2ab07f8b33 100644 --- a/numpy/core/src/umath/funcs.inc.src +++ b/numpy/core/src/umath/funcs.inc.src @@ -8,6 +8,7 @@ #define NPY_NO_DEPRECATED_API NPY_API_VERSION #include "npy_pycompat.h" +#include "npy_import.h" /* @@ -158,6 +159,73 @@ npy_ObjectLogicalNot(PyObject *i1) } } +static PyObject * +npy_ObjectGCD(PyObject *i1, PyObject *i2) +{ + PyObject *gcd = NULL; + + /* use math.gcd if available, and valid on the provided types */ +#if PY_VERSION_HEX >= 0x03050000 + { + static PyObject *math_gcd_func = NULL; + + npy_cache_import("math", "gcd", &math_gcd_func); + if (math_gcd_func == NULL) { + return NULL; + } + gcd = PyObject_CallFunction(math_gcd_func, "OO", i1, i2); + if (gcd != NULL) { + return gcd; + } + /* silence errors, and fall back on pure-python gcd */ + PyErr_Clear(); + } +#endif + + /* otherwise, use our internal one, written in python */ + { + static PyObject *internal_gcd_func = NULL; + + npy_cache_import("numpy.core._internal", "_gcd", &internal_gcd_func); + if (internal_gcd_func == NULL) { + return NULL; + } + gcd = PyObject_CallFunction(internal_gcd_func, "OO", i1, i2); + if (gcd == NULL) { + return NULL; + } + /* _gcd has some unusual behaviour regarding sign */ + return PyNumber_Absolute(gcd); + } +} + +static PyObject * +npy_ObjectLCM(PyObject *i1, PyObject *i2) +{ + /* lcm(a, b) = abs(a // gcd(a, b) * b) */ + + PyObject *gcd = npy_ObjectGCD(i1, i2); + PyObject *tmp; + if(gcd == NULL) { + return NULL; + } + /* Floor divide preserves integer types - we know the division will have + * no remainder + */ + tmp = PyNumber_FloorDivide(i1, gcd); + if(tmp == NULL) { + return NULL; + } + + tmp = PyNumber_Multiply(tmp, i2); + if(tmp == NULL) { + return NULL; + } + + /* even though we fix gcd to be positive, we need to do it again here */ + return PyNumber_Absolute(tmp); +} + /* ***************************************************************************** ** COMPLEX FUNCTIONS ** diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src index 789717555738..8791788d0862 100644 --- a/numpy/core/src/umath/loops.c.src +++ b/numpy/core/src/umath/loops.c.src @@ -1041,6 +1041,7 @@ NPY_NO_EXPORT void /**begin repeat * #TYPE = BYTE, SHORT, INT, LONG, LONGLONG# * #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong# + * #c = ,,,l,ll# */ NPY_NO_EXPORT NPY_GCC_OPT_3 void @@ -1132,11 +1133,26 @@ NPY_NO_EXPORT void } } +/**begin repeat1 + * #kind = gcd, lcm# + **/ +NPY_NO_EXPORT void +@TYPE@_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) +{ + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const @type@ in2 = *(@type@ *)ip2; + *((@type@ *)op1) = npy_@kind@@c@(in1, in2); + } +} +/**end repeat1**/ + /**end repeat**/ /**begin repeat * #TYPE = UBYTE, USHORT, UINT, ULONG, ULONGLONG# * #type = npy_ubyte, npy_ushort, npy_uint, npy_ulong, npy_ulonglong# + * #c = u,u,u,ul,ull# */ NPY_NO_EXPORT void @@ -1204,6 +1220,20 @@ NPY_NO_EXPORT void } } +/**begin repeat1 + * #kind = gcd, lcm# + **/ +NPY_NO_EXPORT void +@TYPE@_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) +{ + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const @type@ in2 = *(@type@ *)ip2; + *((@type@ *)op1) = npy_@kind@@c@(in1, in2); + } +} +/**end repeat1**/ + /**end repeat**/ /* diff --git a/numpy/core/src/umath/loops.h.src b/numpy/core/src/umath/loops.h.src index a978b03ee349..a01ef152953f 100644 --- a/numpy/core/src/umath/loops.h.src +++ b/numpy/core/src/umath/loops.h.src @@ -140,6 +140,12 @@ NPY_NO_EXPORT void NPY_NO_EXPORT void @S@@TYPE@_divmod(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +@S@@TYPE@_gcd(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)); + +NPY_NO_EXPORT void +@S@@TYPE@_lcm(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)); + /**end repeat1**/ /**end repeat**/ diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py index bebeddc92fd3..93764e7b7ced 100644 --- a/numpy/core/tests/test_umath.py +++ b/numpy/core/tests/test_umath.py @@ -2203,6 +2203,105 @@ def test_mixed(self): assert_equal(np.choose(c, (a, 1)), np.array([1, 1])) +class TestRationalFunctions(object): + def test_lcm(self): + self._test_lcm_inner(np.int16) + self._test_lcm_inner(np.uint16) + + def test_lcm_object(self): + self._test_lcm_inner(np.object_) + + def test_gcd(self): + self._test_gcd_inner(np.int16) + self._test_lcm_inner(np.uint16) + + def test_gcd_object(self): + self._test_gcd_inner(np.object_) + + def _test_lcm_inner(self, dtype): + # basic use + a = np.array([12, 120], dtype=dtype) + b = np.array([20, 200], dtype=dtype) + assert_equal(np.lcm(a, b), [60, 600]) + + if not issubclass(dtype, np.unsignedinteger): + # negatives are ignored + a = np.array([12, -12, 12, -12], dtype=dtype) + b = np.array([20, 20, -20, -20], dtype=dtype) + assert_equal(np.lcm(a, b), [60]*4) + + # reduce + a = np.array([3, 12, 20], dtype=dtype) + assert_equal(np.lcm.reduce([3, 12, 20]), 60) + + # broadcasting, and a test including 0 + a = np.arange(6).astype(dtype) + b = 20 + assert_equal(np.lcm(a, b), [0, 20, 20, 60, 20, 20]) + + def _test_gcd_inner(self, dtype): + # basic use + a = np.array([12, 120], dtype=dtype) + b = np.array([20, 200], dtype=dtype) + assert_equal(np.gcd(a, b), [4, 40]) + + if not issubclass(dtype, np.unsignedinteger): + # negatives are ignored + a = np.array([12, -12, 12, -12], dtype=dtype) + b = np.array([20, 20, -20, -20], dtype=dtype) + assert_equal(np.gcd(a, b), [4]*4) + + # reduce + a = np.array([15, 25, 35], dtype=dtype) + assert_equal(np.gcd.reduce(a), 5) + + # broadcasting, and a test including 0 + a = np.arange(6).astype(dtype) + b = 20 + assert_equal(np.gcd(a, b), [20, 1, 2, 1, 4, 5]) + + def test_lcm_overflow(self): + # verify that we don't overflow when a*b does overflow + big = np.int32(np.iinfo(np.int32).max // 11) + a = 2*big + b = 5*big + assert_equal(np.lcm(a, b), 10*big) + + def test_gcd_overflow(self): + for dtype in (np.int32, np.int64): + # verify that we don't overflow when taking abs(x) + # not relevant for lcm, where the result is unrepresentable anyway + a = dtype(np.iinfo(dtype).min) # negative power of two + q = -(a // 4) + assert_equal(np.gcd(a, q*3), q) + assert_equal(np.gcd(a, -q*3), q) + + def test_decimal(self): + from decimal import Decimal + a = np.array([1, 1, -1, -1]) * Decimal('0.20') + b = np.array([1, -1, 1, -1]) * Decimal('0.12') + + assert_equal(np.gcd(a, b), 4*[Decimal('0.04')]) + assert_equal(np.lcm(a, b), 4*[Decimal('0.60')]) + + def test_float(self): + # not well-defined on float due to rounding errors + assert_raises(TypeError, np.gcd, 0.3, 0.4) + assert_raises(TypeError, np.lcm, 0.3, 0.4) + + def test_builtin_long(self): + # sanity check that array coercion is alright for builtin longs + assert_equal(np.array(2**200).item(), 2**200) + + # expressed as prime factors + a = np.array(2**100 * 3**5) + b = np.array([2**100 * 5**7, 2**50 * 3**10]) + assert_equal(np.gcd(a, b), [2**100, 2**50 * 3**5]) + assert_equal(np.lcm(a, b), [2**100 * 3**5 * 5**7, 2**100 * 3**10]) + + assert_equal(np.gcd(2**100, 3**100), 1) + + def is_longdouble_finfo_bogus(): info = np.finfo(np.longcomplex) return not np.isfinite(np.log10(info.tiny/info.eps))