numpy · ahaldane · Nov 8, 2017 · Nov 5, 2017 · eric-wieser · Nov 7, 2017
diff --git 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
@@ -27,6 +27,7 @@
 #include "alloc.h"
 #include "npy_import.h"
 #include "dragon4.h"
+#include "npy_longdouble.h"
 
 #include <stdlib.h>
 
@@ -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
@@ -0,0 +1,99 @@
+#include <Python.h>
+
+#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
@@ -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
@@ -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
@@ -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']: