Test and fix comparisons with NaNs

quaddtype/numpy_quaddtype/src/ops.hpp

@@ -302,13 +302,17 @@ quad_mod(Sleef_quad *a, Sleef_quad *b)
 static inline Sleef_quad
 quad_minimum(Sleef_quad *in1, Sleef_quad *in2)
 {
-    return Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2;
+    return Sleef_iunordq1(*in1, *in2) ? (
+        Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
+    ) : Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2;
 }
 
 static inline Sleef_quad
 quad_maximum(Sleef_quad *in1, Sleef_quad *in2)
 {
-    return Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
+    return Sleef_iunordq1(*in1, *in2) ? (
+        Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
+    ) : Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
 }
 
 static inline Sleef_quad
@@ -386,7 +390,7 @@ quad_equal(const Sleef_quad *a, const Sleef_quad *b)
 static inline npy_bool
 quad_notequal(const Sleef_quad *a, const Sleef_quad *b)
 {
-    return Sleef_icmpneq1(*a, *b);
+    return Sleef_icmpneq1(*a, *b) || Sleef_iunordq1(*a, *b);
 }
 
 static inline npy_bool

quaddtype/numpy_quaddtype/src/scalar_ops.cpp

@@ -163,7 +163,7 @@ quad_richcompare(QuadPrecisionObject *self, PyObject *other, int cmp_op)
                 cmp = Sleef_icmpeqq1(self->value.sleef_value, other_quad->value.sleef_value);
                 break;
             case Py_NE:
-                cmp = Sleef_icmpneq1(self->value.sleef_value, other_quad->value.sleef_value);
+                cmp = Sleef_icmpneq1(self->value.sleef_value, other_quad->value.sleef_value) || Sleef_iunordq1(self->value.sleef_value, other_quad->value.sleef_value);
                 break;
             case Py_GT:
                 cmp = Sleef_icmpgtq1(self->value.sleef_value, other_quad->value.sleef_value);

quaddtype/tests/test_quaddtype.py

@@ -18,8 +18,11 @@ def test_basic_equality():
 
 
 @pytest.mark.parametrize("op", ["add", "sub", "mul", "truediv", "pow"])
-@pytest.mark.parametrize("other", ["3.0", "12.5", "100.0"])
+@pytest.mark.parametrize("other", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
 def test_binary_ops(op, other):
+    if op == "truediv" and float(other) == 0:
+        pytest.xfail("float division by zero")
+
     op_func = getattr(operator, op)
     quad_a = QuadPrecision("12.5")
     quad_b = QuadPrecision(other)
@@ -29,21 +32,80 @@ def test_binary_ops(op, other):
     quad_result = op_func(quad_a, quad_b)
     float_result = op_func(float_a, float_b)
 
-    assert np.abs(np.float64(quad_result) - float_result) < 1e-10
+    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))
+        )
 
 
 @pytest.mark.parametrize("op", ["eq", "ne", "le", "lt", "ge", "gt"])
-@pytest.mark.parametrize("other", ["3.0", "12.5", "100.0"])
-def test_comparisons(op, other):
+@pytest.mark.parametrize("a", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+@pytest.mark.parametrize("b", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+def test_comparisons(op, a, b):
     op_func = getattr(operator, op)
-    quad_a = QuadPrecision("12.5")
-    quad_b = QuadPrecision(other)
-    float_a = 12.5
-    float_b = float(other)
+    quad_a = QuadPrecision(a)
+    quad_b = QuadPrecision(b)
+    float_a = float(a)
+    float_b = float(b)
 
     assert op_func(quad_a, quad_b) == op_func(float_a, float_b)
 
 
+@pytest.mark.parametrize("op", ["eq", "ne", "le", "lt", "ge", "gt"])
+@pytest.mark.parametrize("a", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+@pytest.mark.parametrize("b", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+def test_array_comparisons(op, a, b):
+    op_func = getattr(operator, op)
+    quad_a = np.array(QuadPrecision(a))
+    quad_b = np.array(QuadPrecision(b))
+    float_a = np.array(float(a))
+    float_b = np.array(float(b))
+
+    assert np.array_equal(op_func(quad_a, quad_b), op_func(float_a, float_b))
+
+
+@pytest.mark.parametrize("op", ["minimum", "maximum", "fmin", "fmax"])
+@pytest.mark.parametrize("a", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+@pytest.mark.parametrize("b", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+def test_array_minmax(op, a, b):
+    if op in ["fmin", "fmax"]:
+        pytest.skip("fmin and fmax ufuncs are not yet supported")
+
+    op_func = getattr(np, op)
+    quad_a = np.array([QuadPrecision(a)])
+    quad_b = np.array([QuadPrecision(b)])
+    float_a = np.array([float(a)])
+    float_b = np.array([float(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)))
+
+
+@pytest.mark.parametrize("op", ["amin", "amax", "nanmin", "nanmax"])
+@pytest.mark.parametrize("a", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+@pytest.mark.parametrize("b", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+def test_array_aminmax(op, a, b):
+    if op in ["nanmin", "nanmax"]:
+        pytest.skip("fmin and fmax ufuncs are not yet supported")
+
+    op_func = getattr(np, op)
+    quad_ab = np.array([QuadPrecision(a), QuadPrecision(b)])
+    float_ab = np.array([float(a), float(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)))
+
+
 @pytest.mark.parametrize("op, val, expected", [
     ("neg", "3.0", "-3.0"),
     ("neg", "-3.0", "3.0"),