Implement isfinite, isinf, isnan ufuncs

Author: juntyr

quaddtype/numpy_quaddtype/src/ops.hpp

@@ -294,6 +294,28 @@ quad_signbit(const Sleef_quad *op)
     return Sleef_icmpltq1(one_signed, zero);
 }
 
+static inline npy_bool
+quad_isfinite(const Sleef_quad *op)
+{
+    // isfinite(x) = abs(x) < inf
+    Sleef_quad inf = Sleef_cast_from_doubleq1(std::numeric_limits<double>::infinity());
+    return Sleef_icmpltq1(Sleef_fabsq1(*op), inf);
+}
+
+static inline npy_bool
+quad_isinf(const Sleef_quad *op)
+{
+    // isinf(x) = abs(x) == inf
+    Sleef_quad inf = Sleef_cast_from_doubleq1(std::numeric_limits<double>::infinity());
+    return Sleef_icmpeqq1(Sleef_fabsq1(*op), inf);
+}
+
+static inline npy_bool
+quad_isnan(const Sleef_quad *op)
+{
+    return Sleef_iunordq1(*op, *op);
+}
+
 // Unary Quad properties
 typedef npy_bool (*unary_prop_longdouble_def)(const long double *);
 
@@ -303,6 +325,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.cpp

@@ -387,6 +387,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

@@ -32,13 +32,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"])
@@ -83,9 +77,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"])
@@ -102,12 +94,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,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign"), (None, "signbit")])
+@pytest.mark.parametrize("op,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign"), (None, "signbit"), (None, "isfinite"), (None, "isinf"), (None, "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, nop, val):
     op_func = None if op is None else getattr(operator, op)
@@ -123,14 +113,8 @@ def test_unary_ops(op, nop, val):
         quad_result = op_func(quad_val)
         float_result = op_func(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():