Implement the fmin and fmax ufuncs

Author: juntyr

quaddtype/numpy_quaddtype/src/ops.hpp

@@ -358,6 +358,22 @@ quad_maximum(const Sleef_quad *in1, const Sleef_quad *in2)
     ) : Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
 }
 
+static inline Sleef_quad
+quad_fmin(const Sleef_quad *in1, const Sleef_quad *in2)
+{
+    return Sleef_iunordq1(*in1, *in1) ? (
+        Sleef_iunordq1(*in2, *in2) ? *in1 : *in2
+    ) : Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2;
+}
+
+static inline Sleef_quad
+quad_fmax(const Sleef_quad *in1, const Sleef_quad *in2)
+{
+    return Sleef_iunordq1(*in1, *in1) ? (
+        Sleef_iunordq1(*in2, *in2) ? *in1 : *in2
+    ) : Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
+}
+
 static inline Sleef_quad
 quad_atan2(const Sleef_quad *in1, const Sleef_quad *in2)
 {
@@ -412,13 +428,25 @@ ld_mod(const long double *a, const long double *b)
 static inline long double
 ld_minimum(const long double *in1, const long double *in2)
 {
-    return (*in1 < *in2) ? *in1 : *in2;
+    return isnan(*in1) ? *in1 : (*in1 < *in2) ? *in1 : *in2;
 }
 
 static inline long double
 ld_maximum(const long double *in1, const long double *in2)
 {
-    return (*in1 > *in2) ? *in1 : *in2;
+    return isnan(*in1) ? *in1 : (*in1 > *in2) ? *in1 : *in2;
+}
+
+static inline long double
+ld_fmin(const long double *in1, const long double *in2)
+{
+    return fmin(*in1, *in2);
+}
+
+static inline long double
+ld_fmax(const long double *in1, const long double *in2)
+{
+    return fmax(*in1, *in2);
 }
 
 static inline long double

quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp

@@ -228,6 +228,12 @@ init_quad_binary_ops(PyObject *numpy)
     if (create_quad_binary_ufunc<quad_maximum, ld_maximum>(numpy, "maximum") < 0) {
         return -1;
     }
+    if (create_quad_binary_ufunc<quad_fmin, ld_fmin>(numpy, "fmin") < 0) {
+        return -1;
+    }
+    if (create_quad_binary_ufunc<quad_fmax, ld_fmax>(numpy, "fmax") < 0) {
+        return -1;
+    }
     if (create_quad_binary_ufunc<quad_atan2, ld_atan2>(numpy, "arctan2") < 0) {
         return -1;
     }

quaddtype/tests/test_quaddtype.py

@@ -83,9 +83,6 @@ def test_array_comparisons(op, a, b):
 @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)])
@@ -104,9 +101,6 @@ def test_array_minmax(op, a, b):
 @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)])