4 files changed, 153 insertions, 29 deletions

diff --git a/numpy/core/src/npymath/halffloat.c b/numpy/core/src/npymath/halffloat.c
index da92d3e12..79f77197b 100644
--- a/numpy/core/src/npymath/halffloat.c
+++ b/numpy/core/src/npymath/halffloat.c
@@ -77,11 +77,16 @@ npy_half npy_half_spacing(npy_half h)
     npy_half ret;
     npy_uint16 h_exp = h&0x7c00u;
     npy_uint16 h_sig = h&0x03ffu;
-    if (h_exp == 0x7c00u || h == 0x7bffu) {
+    if (h_exp == 0x7c00u) {
 #if NPY_HALF_GENERATE_INVALID
-        generate_invalid_error();
+        npy_set_floatstatus_invalid();
 #endif
         ret = NPY_HALF_NAN;
+    } else if (h == 0x7bffu) {
+#if NPY_HALF_GENERATE_OVERFLOW
+        npy_set_floatstatus_overflow();
+#endif
+        ret = NPY_HALF_PINF;
     } else if ((h&0x8000u) && h_sig == 0) { /* Negative boundary case */
         if (h_exp > 0x2c00u) { /* If result is normalized */
             ret = h_exp - 0x2c00u;
@@ -112,7 +117,7 @@ npy_half npy_half_nextafter(npy_half x, npy_half y)
 
     if (!npy_half_isfinite(x) || npy_half_isnan(y)) {
 #if NPY_HALF_GENERATE_INVALID
-        generate_invalid_error();
+        npy_set_floatstatus_invalid();
 #endif
         ret = NPY_HALF_NAN;
     } else if (npy_half_eq_nonan(x, y)) {
@@ -134,7 +139,7 @@ npy_half npy_half_nextafter(npy_half x, npy_half y)
     }
 #ifdef NPY_HALF_GENERATE_OVERFLOW
     if (npy_half_isinf(ret)) {
-        generate_overflow_error();
+        npy_set_floatstatus_overflow();
     }
 #endif
 
@@ -255,7 +260,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
         } else {
             /* overflow to signed inf */
 #if NPY_HALF_GENERATE_OVERFLOW
-            generate_overflow_error();
+            npy_set_floatstatus_overflow();
 #endif
             return (npy_uint16) (h_sgn + 0x7c00u);
         }
@@ -271,7 +276,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
 #if NPY_HALF_GENERATE_UNDERFLOW 
             /* If f != 0, it underflowed to 0 */
             if ((f&0x7fffffff) != 0) {
-                generate_underflow_error();
+                npy_set_floatstatus_underflow();
             }
 #endif
             return h_sgn;
@@ -282,7 +287,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
 #if NPY_HALF_GENERATE_UNDERFLOW 
         /* If it's not exactly represented, it underflowed */
         if ((f_sig&(((npy_uint32)1 << (126 - f_exp)) - 1)) != 0) {
-            generate_underflow_error();
+            npy_set_floatstatus_underflow();
         }
 #endif
         f_sig >>= (113 - f_exp);
@@ -334,7 +339,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
 #if NPY_HALF_GENERATE_OVERFLOW
     h_sig += h_exp;
     if (h_sig == 0x7c00u) {
-        generate_overflow_error();
+        npy_set_floatstatus_overflow();
     }
     return h_sgn + h_sig;
 #else
@@ -370,7 +375,7 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
         } else {
             /* overflow to signed inf */
 #if NPY_HALF_GENERATE_OVERFLOW
-            generate_overflow_error();
+            npy_set_floatstatus_overflow();
 #endif
             return h_sgn + 0x7c00u;
         }
@@ -385,8 +390,8 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
         if (d_exp < 0x3e60000000000000u) {
 #if NPY_HALF_GENERATE_UNDERFLOW 
             /* If d != 0, it underflowed to 0 */
-            if ((d&0x7fffffffffffffff) != 0) {
-                generate_underflow_error();
+            if ((d&0x7fffffffffffffffu) != 0) {
+                npy_set_floatstatus_underflow();
             }
 #endif
             return h_sgn;
@@ -397,7 +402,7 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
 #if NPY_HALF_GENERATE_UNDERFLOW 
         /* If it's not exactly represented, it underflowed */
         if ((d_sig&(((npy_uint64)1 << (1051 - d_exp)) - 1)) != 0) {
-            generate_underflow_error();
+            npy_set_floatstatus_underflow();
         }
 #endif
         d_sig >>= (1009 - d_exp);
@@ -450,7 +455,7 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
 #if NPY_HALF_GENERATE_OVERFLOW
     h_sig += h_exp;
     if (h_sig == 0x7c00u) {
-        generate_overflow_error();
+        npy_set_floatstatus_overflow();
     }
     return h_sgn + h_sig;
 #else
diff --git a/numpy/core/src/npymath/ieee754.c.src b/numpy/core/src/npymath/ieee754.c.src
index 8df903b2b..14d3f4a90 100644
--- a/numpy/core/src/npymath/ieee754.c.src
+++ b/numpy/core/src/npymath/ieee754.c.src
@@ -548,3 +548,122 @@ npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
     return nextafterl(x, y);
 }
 #endif
+
+/*
+ * Functions to set the floating point status word.
+ */
+
+#if defined(sun) || defined(__BSD__) || defined(__OpenBSD__) || (defined(__FreeBSD__) && (__FreeBSD_version < 502114)) || defined(__NetBSD__)
+#include <ieeefp.h>
+
+void npy_set_floatstatus_divbyzero(void)
+{
+    fpsetsticky(FP_X_DZ);
+} 
+
+void npy_set_floatstatus_overflow(void)
+{
+    fpsetsticky(FP_X_OFL);
+}
+
+void npy_set_floatstatus_underflow(void)
+{
+    fpsetsticky(FP_X_UFL);
+}
+
+void npy_set_floatstatus_invalid(void)
+{
+    fpsetsticky(FP_X_INV);
+}
+
+
+#elif defined(__GLIBC__) || defined(__APPLE__) || defined(__CYGWIN__) || defined(__MINGW32__) || (defined(__FreeBSD__) && (__FreeBSD_version >= 502114))
+
+# if defined(__GLIBC__) || defined(__APPLE__) || defined(__MINGW32__) || defined(__FreeBSD__)
+#  include <fenv.h>
+# elif defined(__CYGWIN__)
+#  include "fenv/fenv.c"
+# endif
+
+void npy_set_floatstatus_divbyzero(void)
+{
+    feraiseexcept(FE_DIVBYZERO);
+} 
+
+void npy_set_floatstatus_overflow(void)
+{
+    feraiseexcept(FE_OVERFLOW);
+}
+
+void npy_set_floatstatus_underflow(void)
+{
+    feraiseexcept(FE_UNDERFLOW);
+}
+
+void npy_set_floatstatus_invalid(void)
+{
+    feraiseexcept(FE_INVALID);
+}
+
+#elif defined(_AIX)
+#include <float.h>
+#include <fpxcp.h>
+
+void npy_set_floatstatus_divbyzero(void)
+{
+    fp_raise_xcp(FP_DIV_BY_ZERO);
+} 
+
+void npy_set_floatstatus_overflow(void)
+{
+    fp_raise_xcp(FP_OVERFLOW);
+}
+
+void npy_set_floatstatus_underflow(void)
+{
+    fp_raise_xcp(FP_UNDERFLOW);
+}
+
+void npy_set_floatstatus_invalid(void)
+{
+    fp_raise_xcp(FP_INVALID);
+}
+
+#else
+
+/*
+ * By using a volatile floating point value,
+ * the compiler is forced to actually do the requested
+ * operations because of potential concurrency.
+ *
+ * We shouldn't write multiple values to a single
+ * global here, because that would cause
+ * a race condition.
+ */
+static volatile double _npy_floatstatus_x,
+    _npy_floatstatus_zero = 0.0, _npy_floatstatus_big = 1e300,
+    _npy_floatstatus_small = 1e-300, _npy_floatstatus_inf;
+
+void npy_set_floatstatus_divbyzero(void)
+{
+    _npy_floatstatus_x = 1.0 / _npy_floatstatus_zero;
+}
+
+void npy_set_floatstatus_overflow(void)
+{
+    _npy_floatstatus_x = _npy_floatstatus_big * 1e300;
+}
+
+void npy_set_floatstatus_underflow(void)
+{
+    _npy_floatstatus_x = _npy_floatstatus_small * 1e-300;
+}
+
+void npy_set_floatstatus_invalid(void)
+{
+    _npy_floatstatus_inf = NPY_INFINITY;
+    _npy_floatstatus_x = _npy_floatstatus_inf - NPY_INFINITY;
+}
+
+#endif
+
diff --git a/numpy/core/src/scalarmathmodule.c.src b/numpy/core/src/scalarmathmodule.c.src
index 712932958..d1a66e101 100644
--- a/numpy/core/src/scalarmathmodule.c.src
+++ b/numpy/core/src/scalarmathmodule.c.src
@@ -145,7 +145,7 @@ static void
     if ((*out^a) >= 0 || (*out^b) >= 0) {
         return;
     }
-    generate_overflow_error();
+    npy_set_floatstatus_overflow();
     return;
 }
 static void
@@ -154,7 +154,7 @@ static void
     if ((*out^a) >= 0 || (*out^~b) >= 0) {
         return;
     }
-    generate_overflow_error();
+    npy_set_floatstatus_overflow();
     return;
 }
 /**end repeat**/
@@ -168,7 +168,7 @@ static void
     if (*out >= a && *out >= b) {
         return;
     }
-    generate_overflow_error();
+    npy_set_floatstatus_overflow();
     return;
 }
 static void
@@ -177,7 +177,7 @@ static void
     if (a >= b) {
         return;
     }
-    generate_overflow_error();
+    npy_set_floatstatus_overflow();
     return;
 }
 /**end repeat**/
@@ -206,7 +206,7 @@ static void
 #else
         if (temp > MAX_@NAME@)
 #endif
-            generate_overflow_error();
+            npy_set_floatstatus_overflow();
     return;
 }
 #endif
@@ -223,7 +223,7 @@ static void
 @name@_ctype_multiply(@name@ a, @name@ b, @name@ *out) {
     *out = a * b;
     if (@char@longlong_overflow(a, b)) {
-        generate_overflow_error();
+        npy_set_floatstatus_overflow();
     }
     return;
 }
@@ -239,12 +239,12 @@ static void
 static void
 @name@_ctype_divide(@name@ a, @name@ b, @name@ *out) {
     if (b == 0) {
-        generate_divbyzero_error();
+        npy_set_floatstatus_divbyzero();
         *out = 0;
     }
 #if @neg@
     else if (b == -1 && a < 0 && a == -a) {
-        generate_overflow_error();
+        npy_set_floatstatus_overflow();
         *out = a / b;
     }
 #endif
@@ -266,7 +266,7 @@ static void
 static void
 @name@_ctype_remainder(@name@ a, @name@ b, @name@ *out) {
     if (a == 0 || b == 0) {
-        if (b == 0) generate_divbyzero_error();
+        if (b == 0) npy_set_floatstatus_divbyzero();
         *out = 0;
         return;
     }
@@ -450,7 +450,7 @@ static void
 @name@_ctype_negative(npy_@name@ a, npy_@name@ *out)
 {
 #if @uns@
-    generate_overflow_error();
+    npy_set_floatstatus_overflow();
 #endif
     *out = -a;
 }
diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 0c777e464..8c98db6b7 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -795,7 +795,7 @@ NPY_NO_EXPORT void
         const @s@@type@ in1 = *(@s@@type@ *)ip1;
         const @s@@type@ in2 = *(@s@@type@ *)ip2;
         if (in2 == 0) {
-            generate_divbyzero_error();
+            npy_set_floatstatus_divbyzero();
             *((@s@@type@ *)op1) = 0;
         }
         else {
@@ -854,10 +854,10 @@ NPY_NO_EXPORT void
          * by -1 causes a SIFGPE (division overflow). We treat this case here
          * (to avoid a SIGFPE crash at python level), but a good solution would
          * be to treat integer division problems separately from FPU exceptions
-         * (i.e. fixing generate_divbyzero_error()).
+         * (i.e. a different approach than npy_set_floatstatus_divbyzero()).
          */
         if (in2 == 0 || (in1 == NPY_MIN_@TYPE@ && in2 == -1)) {
-            generate_divbyzero_error();
+            npy_set_floatstatus_divbyzero();
             *((@type@ *)op1) = 0;
         }
         else if (((in1 > 0) != (in2 > 0)) && (in1 % in2 != 0)) {
@@ -876,7 +876,7 @@ U@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func
         const u@type@ in1 = *(u@type@ *)ip1;
         const u@type@ in2 = *(u@type@ *)ip2;
         if (in2 == 0) {
-            generate_divbyzero_error();
+            npy_set_floatstatus_divbyzero();
             *((u@type@ *)op1) = 0;
         }
         else {
@@ -892,7 +892,7 @@ NPY_NO_EXPORT void
         const @type@ in1 = *(@type@ *)ip1;
         const @type@ in2 = *(@type@ *)ip2;
         if (in2 == 0) {
-            generate_divbyzero_error();
+            npy_set_floatstatus_divbyzero();
             *((@type@ *)op1) = 0;
         }
         else {
@@ -915,7 +915,7 @@ U@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(f
         const u@type@ in1 = *(u@type@ *)ip1;
         const u@type@ in2 = *(u@type@ *)ip2;
         if (in2 == 0) {
-            generate_divbyzero_error();
+            npy_set_floatstatus_divbyzero();
             *((@type@ *)op1) = 0;
         }
         else {