--- a/mlx/ops.cpp
+++ b/mlx/ops.cpp
@@ -2410,10 +2410,20 @@
   if (flat_size % 2 == 0) {
     start.back() = mp - 1;
     stop.back() = mp;
-    median_a = multiply(
-        add(median_a, astype(slice(sorted_a, start, stop, s), dtype, s), s),
-        array(0.5, dtype),
-        s);
+    auto lower = astype(slice(sorted_a, start, stop, s), dtype, s);
+    auto half = array(0.5, dtype);
+    auto average = multiply(add(median_a, lower, s), half, s);
+    if (dtype == float32 || dtype == float16 || dtype == bfloat16) {
+      float limit = dtype == float16 ? float(numeric_limits<float16_t>::max())
+          : dtype == bfloat16 ? float(numeric_limits<bfloat16_t>::max())
+          : numeric_limits<float>::max();
+      auto large = logical_or(
+          greater(abs(median_a, s), array(limit * 0.5f, dtype), s),
+          greater(abs(lower, s), array(limit * 0.5f, dtype), s), s);
+      auto scaled = add(multiply(median_a, half, s), multiply(lower, half, s), s);
+      average = where(large, scaled, average, s);
+    }
+    median_a = average;
   }
   // Sorting moves NaN to the end, so the midpoint slice never selects it.
   // Propagate it explicitly to stay consistent with max, min and mean.