MAINT: core: Restore float128 power on Mac OSX.

WarrenWeckesser · WarrenWeckesser · commit 9875e47c3e14 · 2020-08-27T14:52:20.000-04:00
The use of the math library powl was blacklisted on Mac OSX (gh-8318) because powl(0, y) would trigger a divide-by-zero warning even when y > 0 (issue gh-8307). This change creates npy_powl on OSX that avoids the spurious warning by not calling powl when x = 0 and y > 0. The change in npy_math_internal.h.src moves the handling of pow out of a template repeat-loop, so powl and powf are handled individually. This allows the creation of npy_powl to be specialized when NPY_OS_DARWIN is defined. Closes gh-8608.
diff --git a/numpy/core/src/npymath/npy_math_internal.h.src b/numpy/core/src/npymath/npy_math_internal.h.src
@@ -55,6 +55,7 @@
  */
 #include "npy_math_private.h"
 
+
 /*
  *****************************************************************************
  **                     BASIC MATH FUNCTIONS                                **
@@ -398,8 +399,8 @@ NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
 /**end repeat1**/
 
 /**begin repeat1
- * #kind = atan2,hypot,pow,fmod,copysign#
- * #KIND = ATAN2,HYPOT,POW,FMOD,COPYSIGN#
+ * #kind = atan2,hypot,fmod,copysign#
+ * #KIND = ATAN2,HYPOT,FMOD,COPYSIGN#
  */
 #ifdef @kind@@c@
 #undef @kind@@c@
@@ -448,6 +449,37 @@ NPY_INPLACE @type@ npy_frexp@c@(@type@ x, int* exp)
 /**end repeat**/
 
 
+#ifdef powf
+#undef powf
+#endif
+#ifndef HAVE_POWF
+NPY_INPLACE npy_float npy_powf(npy_float x, npy_float y)
+{
+    return (npy_float) npy_pow((double) x, (double) y);
+}
+#endif
+
+#ifdef powl
+#undef powl
+#endif
+#ifndef HAVE_POWL
+#if defined(NPY_OS_DARWIN)
+NPY_INPLACE npy_longdouble npy_powl(npy_longdouble x, npy_longdouble y)
+{
+    if (x == 0 && y > 0) {
+        return 0.0L;
+    }
+    return powl(x, y);
+}
+#else
+// Fall back to the double precision version.
+NPY_INPLACE npy_longdouble npy_powl(npy_longdouble x, npy_longdouble y)
+{
+    return (npy_longdouble) npy_pow((npy_double) x, (npy_double) y);
+}
+#endif
+#endif
+
 /*
  * Decorate all the math functions which are available on the current platform
  */
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
@@ -546,6 +546,21 @@ def test_integer_to_negative_power(self):
             assert_raises(ValueError, np.power, one, b)
             assert_raises(ValueError, np.power, one, minusone)
 
+    @pytest.mark.skipif(not hasattr(np, 'float128'),
+                        reason='test requires float128')
+    def test_float128(self):
+        x = np.array([0, 1, 2], dtype=np.float128)
+        # On Mac OSX, this would trigger a divide-by-zero warning.
+        y = np.power(x, 4)
+        assert_array_equal(y, [0.0, 1.0, 16.0])
+
+    @pytest.mark.skipif(not sys.platform == 'darwin',
+                        reason='only test on Mac OSX (i.e. darwin)')
+    def test_float128_precision(self):
+        # Check that we are using more than double precision.
+        v = np.float128(2) ** 16383
+        assert_allclose(v, np.float128('5.9486574767861588254e+4931'))
+
 
 class TestFloat_power:
     def test_type_conversion(self):
