diff --git a/numpy/core/setup.py b/numpy/core/setup.py index b87070d77ed6..22eb63f45092 100644 --- a/numpy/core/setup.py +++ b/numpy/core/setup.py @@ -753,6 +753,7 @@ def get_mathlib_info(*args): join('src', 'private', 'templ_common.h.src'), join('src', 'private', 'lowlevel_strided_loops.h'), join('src', 'private', 'mem_overlap.h'), + join('src', 'private', 'npy_longdouble.h'), join('src', 'private', 'ufunc_override.h'), join('src', 'private', 'binop_override.h'), join('src', 'private', 'npy_extint128.h'), @@ -827,6 +828,7 @@ def get_mathlib_info(*args): join('src', 'multiarray', 'vdot.c'), join('src', 'private', 'templ_common.h.src'), join('src', 'private', 'mem_overlap.c'), + join('src', 'private', 'npy_longdouble.c'), join('src', 'private', 'ufunc_override.c'), ] @@ -884,6 +886,7 @@ def generate_umath_c(ext, build_dir): join('src', 'umath', 'ufunc_type_resolution.c'), join('src', 'umath', 'override.c'), join('src', 'private', 'mem_overlap.c'), + join('src', 'private', 'npy_longdouble.c'), join('src', 'private', 'ufunc_override.c')] umath_deps = [ @@ -897,6 +900,7 @@ def generate_umath_c(ext, build_dir): join(codegen_dir, 'generate_ufunc_api.py'), join('src', 'private', 'lowlevel_strided_loops.h'), join('src', 'private', 'mem_overlap.h'), + join('src', 'private', 'npy_longdouble.h'), join('src', 'private', 'ufunc_override.h'), join('src', 'private', 'binop_override.h')] + npymath_sources diff --git a/numpy/core/src/multiarray/scalartypes.c.src b/numpy/core/src/multiarray/scalartypes.c.src index 3fa8f14dc3d4..306944d11231 100644 --- a/numpy/core/src/multiarray/scalartypes.c.src +++ b/numpy/core/src/multiarray/scalartypes.c.src @@ -27,6 +27,7 @@ #include "alloc.h" #include "npy_import.h" #include "dragon4.h" +#include "npy_longdouble.h" #include @@ -832,14 +833,7 @@ static PyObject * @char@longdoubletype_long(PyObject *self) { npy_longdouble val = PyArrayScalar_VAL(self, @CHAR@LongDouble)@POST@; - - /* - * This raises OverflowError when the cast overflows to infinity (gh-9964) - * - * Could fix this with a PyLong_FromLongDouble(longdouble ldval) - * but this would need some more work... - */ - return PyLong_FromDouble((double) val); + return npy_longdouble_to_PyLong(val); } #if !defined(NPY_PY3K) diff --git a/numpy/core/src/private/npy_longdouble.c b/numpy/core/src/private/npy_longdouble.c new file mode 100644 index 000000000000..d2f58c86e0bc --- /dev/null +++ b/numpy/core/src/private/npy_longdouble.c @@ -0,0 +1,99 @@ +#include + +#define NPY_NO_DEPRECATED_API NPY_API_VERSION +#include "numpy/ndarraytypes.h" +#include "numpy/npy_math.h" + +/* This is a backport of Py_SETREF */ +#define NPY_SETREF(op, op2) \ + do { \ + PyObject *_py_tmp = (PyObject *)(op); \ + (op) = (op2); \ + Py_DECREF(_py_tmp); \ + } while (0) + + +/* Heavily derived from PyLong_FromDouble + * Notably, we can't set the digits directly, so have to shift and or instead. + */ +PyObject * +npy_longdouble_to_PyLong(npy_longdouble ldval) +{ + PyObject *v; + PyObject *l_chunk_size; + // number of bits to extract at a time. CPython uses 30, but that's because + // it's tied to the internal long representation + const int chunk_size = NPY_BITSOF_LONGLONG; + npy_longdouble frac; + int i, ndig, expo, neg; + neg = 0; + + if (npy_isinf(ldval)) { + PyErr_SetString(PyExc_OverflowError, + "cannot convert longdouble infinity to integer"); + return NULL; + } + if (npy_isnan(ldval)) { + PyErr_SetString(PyExc_ValueError, + "cannot convert longdouble NaN to integer"); + return NULL; + } + if (ldval < 0.0) { + neg = 1; + ldval = -ldval; + } + frac = npy_frexpl(ldval, &expo); /* ldval = frac*2**expo; 0.0 <= frac < 1.0 */ + v = PyLong_FromLong(0L); + if (v == NULL) + return NULL; + if (expo <= 0) + return v; + + ndig = (expo-1) / chunk_size + 1; + + l_chunk_size = PyLong_FromLong(chunk_size); + if (l_chunk_size == NULL) { + Py_DECREF(v); + return NULL; + } + + /* Get the MSBs of the integral part of the float */ + frac = npy_ldexpl(frac, (expo-1) % chunk_size + 1); + for (i = ndig; --i >= 0; ) { + npy_ulonglong chunk = (npy_ulonglong)frac; + PyObject *l_chunk; + /* v = v << chunk_size */ + NPY_SETREF(v, PyNumber_Lshift(v, l_chunk_size)); + if (v == NULL) { + goto done; + } + l_chunk = PyLong_FromUnsignedLongLong(chunk); + if (l_chunk == NULL) { + Py_DECREF(v); + v = NULL; + goto done; + } + /* v = v | chunk */ + NPY_SETREF(v, PyNumber_Or(v, l_chunk)); + Py_DECREF(l_chunk); + if (v == NULL) { + goto done; + } + + /* Remove the msbs, and repeat */ + frac = frac - (npy_longdouble) chunk; + frac = npy_ldexpl(frac, chunk_size); + } + + /* v = -v */ + if (neg) { + NPY_SETREF(v, PyNumber_Negative(v)); + if (v == NULL) { + goto done; + } + } + +done: + Py_DECREF(l_chunk_size); + return v; +} diff --git a/numpy/core/src/private/npy_longdouble.h b/numpy/core/src/private/npy_longdouble.h new file mode 100644 index 000000000000..c0887eec8d33 --- /dev/null +++ b/numpy/core/src/private/npy_longdouble.h @@ -0,0 +1,17 @@ +#ifndef __NPY_LONGDOUBLE_H +#define __NPY_LONGDOUBLE_H + +#include "npy_config.h" +#include "numpy/ndarraytypes.h" + +/* Convert a npy_longdouble to a python `long` integer. + * + * Results are rounded towards zero. + * + * This performs the same task as PyLong_FromDouble, but for long doubles + * which have a greater range. + */ +NPY_NO_EXPORT PyObject * +npy_longdouble_to_PyLong(npy_longdouble ldval); + +#endif diff --git a/numpy/core/src/umath/scalarmath.c.src b/numpy/core/src/umath/scalarmath.c.src index 259adef0a0a5..3b23151f1a1c 100644 --- a/numpy/core/src/umath/scalarmath.c.src +++ b/numpy/core/src/umath/scalarmath.c.src @@ -24,6 +24,7 @@ #include "templ_common.h" #include "binop_override.h" +#include "npy_longdouble.h" /* Basic operations: * @@ -1391,45 +1392,42 @@ emit_complexwarning(void) * * #cmplx = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1# * #sign = (signed, unsigned)*5, , , , , , , # - * #unsigntyp = 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0*7# - * #ctype = long*8, PY_LONG_LONG*2, double*7# + * #ctype = long*8, PY_LONG_LONG*2, + * double*3, npy_longdouble, double*2, npy_longdouble# * #to_ctype = , , , , , , , , , , npy_half_to_double, , , , , , # - * #realtyp = 0*10, 1*7# * #func = (PyLong_FromLong, PyLong_FromUnsignedLong)*4, * PyLong_FromLongLong, PyLong_FromUnsignedLongLong, - * PyLong_FromDouble*7# + * PyLong_FromDouble*3, npy_longdouble_to_PyLong, + * PyLong_FromDouble*2, npy_longdouble_to_PyLong# */ static PyObject * @name@_int(PyObject *obj) { + PyObject *long_result; + #if @cmplx@ - @sign@ @ctype@ x= @to_ctype@(PyArrayScalar_VAL(obj, @Name@).real); - int ret; + @sign@ @ctype@ x = @to_ctype@(PyArrayScalar_VAL(obj, @Name@).real); #else - @sign@ @ctype@ x= @to_ctype@(PyArrayScalar_VAL(obj, @Name@)); -#endif - -#if @realtyp@ - double ix; - modf(x, &ix); - x = ix; + @sign@ @ctype@ x = @to_ctype@(PyArrayScalar_VAL(obj, @Name@)); #endif #if @cmplx@ - ret = emit_complexwarning(); - if (ret < 0) { + if (emit_complexwarning() < 0) { return NULL; } #endif -#if @unsigntyp@ - if(x < LONG_MAX) - return PyInt_FromLong(x); -#else - if(LONG_MIN < x && x < LONG_MAX) - return PyInt_FromLong(x); + long_result = @func@(x); + if (long_result == NULL){ + return NULL; + } + +#ifndef NPY_PY3K + /* Invoke long.__int__ to try to downcast */ + long_result = Py_TYPE(long_result)->tp_as_number->nb_int(long_result); #endif - return @func@(x); + + return long_result; } /**end repeat**/ @@ -1447,18 +1445,18 @@ static PyObject * * #to_ctype = (, , , , , , , , , , npy_half_to_double, , , , , , )*2# * #which = long*17, float*17# * #func = (PyLong_FromLongLong, PyLong_FromUnsignedLongLong)*5, - * PyLong_FromDouble*7, PyFloat_FromDouble*17# + * PyLong_FromDouble*3, npy_longdouble_to_PyLong, + * PyLong_FromDouble*2, npy_longdouble_to_PyLong, + * PyFloat_FromDouble*17# */ static NPY_INLINE PyObject * @name@_@which@(PyObject *obj) { #if @cmplx@ - int ret; - ret = emit_complexwarning(); - if (ret < 0) { + if (emit_complexwarning() < 0) { return NULL; } - return @func@(@to_ctype@((PyArrayScalar_VAL(obj, @Name@)).real)); + return @func@(@to_ctype@(PyArrayScalar_VAL(obj, @Name@).real)); #else return @func@(@to_ctype@(PyArrayScalar_VAL(obj, @Name@))); #endif diff --git a/numpy/core/tests/test_scalarmath.py b/numpy/core/tests/test_scalarmath.py index 1784c56949b9..d3cdd69dcb51 100644 --- a/numpy/core/tests/test_scalarmath.py +++ b/numpy/core/tests/test_scalarmath.py @@ -7,9 +7,10 @@ import numpy as np from numpy.testing import ( - run_module_suite, assert_, assert_equal, assert_raises, - assert_almost_equal, assert_allclose, assert_array_equal, IS_PYPY, - suppress_warnings, dec, _gen_alignment_data, + run_module_suite, + assert_, assert_equal, assert_raises, + assert_almost_equal, assert_allclose, assert_array_equal, + IS_PYPY, suppress_warnings, dec, _gen_alignment_data, ) types = [np.bool_, np.byte, np.ubyte, np.short, np.ushort, np.intc, np.uintc, @@ -398,7 +399,7 @@ def overflow_error_func(dtype): for code in 'lLqQ': assert_raises(OverflowError, overflow_error_func, code) - def test_longdouble_int(self): + def test_int_from_infinite_longdouble(self): # gh-627 x = np.longdouble(np.inf) assert_raises(OverflowError, int, x) @@ -410,7 +411,7 @@ def test_longdouble_int(self): @dec.knownfailureif(not IS_PYPY, "__int__ is not the same as int in cpython (gh-9972)") - def test_clongdouble___int__(self): + def test_int_from_infinite_longdouble___int__(self): x = np.longdouble(np.inf) assert_raises(OverflowError, x.__int__) with suppress_warnings() as sup: @@ -419,6 +420,21 @@ def test_clongdouble___int__(self): assert_raises(OverflowError, x.__int__) assert_equal(len(sup.log), 1) + @dec.skipif(np.finfo(np.double) == np.finfo(np.longdouble)) + def test_int_from_huge_longdouble(self): + # produce a longdouble that would overflow a double + exp = np.finfo(np.double).maxexp + huge_ld = 1234 * np.longdouble(2) ** exp + huge_i = 1234 * 2 ** exp + assert_(huge_ld != np.inf) + assert_equal(int(huge_ld), huge_i) + + def test_int_from_longdouble(self): + x = np.longdouble(1.5) + assert_equal(int(x), 1) + x = np.longdouble(-10.5) + assert_equal(int(x), -10) + def test_numpy_scalar_relational_operators(self): # All integer for dt1 in np.typecodes['AllInteger']: