Implement isfinite, isinf, isnan ufuncs

Author: juntyr

quaddtype/numpy_quaddtype/src/ops.hpp

@@ -1,6 +1,12 @@
 #include <sleef.h>
 #include <sleefquad.h>
 #include <cmath>
+#include <limits>
+
+// Quad Constants, generated with qutil
+#define QUAD_ZERO sleef_q(+0x0000000000000LL, 0x0000000000000000ULL, -16383)
+#define QUAD_ONE sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 0)
+#define QUAD_POS_INF sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 16384)
 
 // Unary Quad Operations
 typedef Sleef_quad (*unary_op_quad_def)(const Sleef_quad *);
@@ -294,7 +300,27 @@ quad_signbit(const Sleef_quad *op)
     return Sleef_icmpltq1(one_signed, zero);
 }
 
-// Unary Quad properties
+static inline npy_bool
+quad_isfinite(const Sleef_quad *op)
+{
+    // isfinite(x) = abs(x) < inf
+    return Sleef_icmpltq1(Sleef_fabsq1(*op), QUAD_POS_INF);
+}
+
+static inline npy_bool
+quad_isinf(const Sleef_quad *op)
+{
+    // isinf(x) = abs(x) == inf
+    return Sleef_icmpeqq1(Sleef_fabsq1(*op), QUAD_POS_INF);
+}
+
+static inline npy_bool
+quad_isnan(const Sleef_quad *op)
+{
+    return Sleef_iunordq1(*op, *op);
+}
+
+// Unary long double properties
 typedef npy_bool (*unary_prop_longdouble_def)(const long double *);
 
 static inline npy_bool
@@ -303,6 +329,24 @@ ld_signbit(const long double *op)
     return signbit(*op);
 }
 
+static inline npy_bool
+ld_isfinite(const long double *op)
+{
+    return isfinite(*op);
+}
+
+static inline npy_bool
+ld_isinf(const long double *op)
+{
+    return isinf(*op);
+}
+
+static inline npy_bool
+ld_isnan(const long double *op)
+{
+    return isnan(*op);
+}
+
 // Binary Quad operations
 typedef Sleef_quad (*binary_op_quad_def)(const Sleef_quad *, const Sleef_quad *);
 

quaddtype/numpy_quaddtype/src/umath/unary_props.cpp

@@ -136,6 +136,15 @@ create_quad_unary_prop_ufunc(PyObject *numpy, const char *ufunc_name)
 int
 init_quad_unary_props(PyObject *numpy)
 {
+    if (create_quad_unary_prop_ufunc<quad_isfinite, ld_isfinite>(numpy, "isfinite") < 0) {
+        return -1;
+    }
+    if (create_quad_unary_prop_ufunc<quad_isinf, ld_isinf>(numpy, "isinf") < 0) {
+        return -1;
+    }
+    if (create_quad_unary_prop_ufunc<quad_isnan, ld_isnan>(numpy, "isnan") < 0) {
+        return -1;
+    }
     if (create_quad_unary_prop_ufunc<quad_signbit, ld_signbit>(numpy, "signbit") < 0) {
         return -1;
     }

quaddtype/tests/test_quaddtype.py

@@ -44,13 +44,7 @@ def test_binary_ops(op, other):
     quad_result = op_func(quad_a, quad_b)
     float_result = op_func(float_a, float_b)
 
-    # FIXME: @juntyr: replace with array_equal once isnan is supported
-    with np.errstate(invalid="ignore"):
-        assert (
-            (np.float64(quad_result) == float_result) or
-            (np.abs(np.float64(quad_result) - float_result) < 1e-10) or
-            ((float_result != float_result) and (quad_result != quad_result))
-        )
+    np.testing.assert_allclose(np.float64(quad_result), float_result, atol=1e-10, rtol=0, equal_nan=True)
 
 
 @pytest.mark.parametrize("op", ["eq", "ne", "le", "lt", "ge", "gt"])
@@ -95,9 +89,7 @@ def test_array_minmax(op, a, b):
     quad_res = op_func(quad_a, quad_b)
     float_res = op_func(float_a, float_b)
 
-    # FIXME: @juntyr: replace with array_equal once isnan is supported
-    with np.errstate(invalid="ignore"):
-        assert np.all((quad_res == float_res) | ((quad_res != quad_res) & (float_res != float_res)))
+    np.testing.assert_array_equal(quad_res.astype(float), float_res)
 
 
 @pytest.mark.parametrize("op", ["amin", "amax", "nanmin", "nanmax"])
@@ -114,12 +106,10 @@ def test_array_aminmax(op, a, b):
     quad_res = op_func(quad_ab)
     float_res = op_func(float_ab)
 
-    # FIXME: @juntyr: replace with array_equal once isnan is supported
-    with np.errstate(invalid="ignore"):
-        assert np.all((quad_res == float_res) | ((quad_res != quad_res) & (float_res != float_res)))
+    np.testing.assert_array_equal(np.array(quad_res).astype(float), float_res)
 
 
-@pytest.mark.parametrize("op", ["negative", "positive", "absolute", "sign", "signbit"])
+@pytest.mark.parametrize("op", ["negative", "positive", "absolute", "sign", "signbit", "isfinite", "isinf", "isnan"])
 @pytest.mark.parametrize("val", ["3.0", "-3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
 def test_unary_ops(op, val):
     op_func = dict(negative=operator.neg, positive=operator.pos, absolute=operator.abs).get(op, None)
@@ -135,14 +125,8 @@ def test_unary_ops(op, val):
         quad_result = of(quad_val)
         float_result = of(float_val)
 
-        # FIXME: @juntyr: replace with array_equal once isnan is supported
-        with np.errstate(invalid="ignore"):
-            assert (
-                (np.float64(quad_result) == float_result) or
-                ((float_result != float_result) and (quad_result != quad_result))
-            ) and (
-                np.signbit(float_result) == np.signbit(quad_result)
-            ), f"{op}({val}) should be {float_result}, but got {quad_result}"
+        np.testing.assert_array_equal(np.array(quad_result).astype(float), float_result)
+        assert np.signbit(float_result) == np.signbit(quad_result)
 
 
 def test_inf():