diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index 171fc30..252fb2d 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -398,6 +398,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, *in2) ? (
+        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, *in2) ? (
+        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)
 {
@@ -452,13 +468,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
diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp
index e6fc9fa..8adfe4d 100644
--- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp
+++ b/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;
     }
diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md
index fa26176..5c998e9 100644
--- a/quaddtype/release_tracker.md
+++ b/quaddtype/release_tracker.md
@@ -75,8 +75,8 @@
 | logical_not   |       |                                                                         |
 | maximum       | ✅    | ✅                                                                      |
 | minimum       | ✅    | ✅                                                                      |
-| fmax          |       |                                                                         |
-| fmin          |       |                                                                         |
+| fmax          | #123  | ✅                                                                      |
+| fmin          | #123  | ✅                                                                      |
 | isfinite      | #121  | ✅                                                                      |
 | isinf         | #121  | ✅                                                                      |
 | isnan         | #221  | ✅                                                                      |
diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py
index d9583e9..ef7c6d0 100644
--- a/quaddtype/tests/test_quaddtype.py
+++ b/quaddtype/tests/test_quaddtype.py
@@ -77,9 +77,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)])
@@ -96,9 +93,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)])