Fix incorrect rounding for subnormals in widening FMA

libm/src/math/fma.rs

@@ -168,6 +168,73 @@ mod tests {
         assert_eq!(fmaf(a, b, c), expected);
     }
 
+    #[test]
+    fn fmaf_subnormal_double_rounding() {
+        // https://github.com/rust-lang/compiler-builtins/issues/1262
+        let a = f32::from_bits(0x9700_0800);
+        let b = f32::from_bits(0x1cff_f001);
+        let c = f32::from_bits(0x0001_0002);
+
+        // The exact result is (65537.5 - 2^-37) * 2^-149, just below the
+        // midpoint between these two subnormal values.
+        let expected = f32::from_bits(0x0001_0001);
+        let result = generic::fma_wide_round::<f32, f64>(a, b, c, Round::Nearest).val;
+        assert_biteq!(result, expected);
+    }
+
+    #[test]
+    fn fmaf_subnormal_round_to_odd_parity() {
+        let a = f32::from_bits(0x15ef_b8d7);
+        let b = f32::from_bits(0x9e08_b110);
+        let c = f32::from_bits(0x8004_c5ce);
+
+        // The widened sum is inexact but already odd. Moving it another ULP
+        // toward the residual would incorrectly round back to c.
+        let product = f64::from(a) * f64::from(b);
+        let widened_sum = product + f64::from(c);
+        assert_eq!(widened_sum.to_bits() & 1, 1);
+
+        let expected = f32::from_bits(0x8004_c5cf);
+        let result = generic::fma_wide_round::<f32, f64>(a, b, c, Round::Nearest).val;
+        assert_biteq!(result, expected);
+    }
+
+    #[test]
+    fn fmaf_round_to_odd_before_overflow() {
+        let a = f32::from_bits(0x5f78_0000);
+        let b = f32::from_bits(0x5f84_2108);
+        let c = f32::from_bits(0x8000_0001);
+
+        // The exact product is 2^128 - 2^103, the round-to-nearest overflow midpoint. Subtracting
+        // the minimum subnormal puts the fused result just below that midpoint, so it rounds to the
+        // maximum finite value rather than infinity.
+        let result = generic::fma_wide_round::<f32, f64>(a, b, c, Round::Nearest).val;
+        assert_biteq!(result, f32::MAX);
+    }
+
+    #[test]
+    fn fmaf_wide_infinity() {
+        let result =
+            generic::fma_wide_round::<f32, f64>(f32::INFINITY, 1.0, 1.0, Round::Nearest).val;
+        assert_biteq!(result, f32::INFINITY);
+    }
+
+    // 32-bit ARM has ABI bugs around f128 as of rustc 1.99.0-nightly, so it's disabled for this test
+    #[test]
+    #[cfg(all(f128_enabled, not(target_arch = "arm")))]
+    fn fma_wide_f64_subnormal_double_rounding() {
+        let a = f64::from_bits(0x9e55_2156_d547_5b4a);
+        let b = f64::from_bits(0x1e58_3b0e_3ea1_8955);
+        let c = f64::from_bits(0x0000_0040_0000_0002);
+
+        // If q = 2^-1074, the exact product is
+        // -(1/2 + 0x5615e992 * 2^-106) * q. The fused result is therefore
+        // just below the midpoint between the expected value and c.
+        let expected = f64::from_bits(0x0000_0040_0000_0001);
+        let result = generic::fma_wide_round::<f64, f128>(a, b, c, Round::Nearest).val;
+        assert_biteq!(result, expected);
+    }
+
     #[test]
     fn fma_segfault() {
         // These two inputs cause fma to segfault on release due to overflow:

libm/src/math/generic/fma_wide.rs

@@ -1,12 +1,16 @@
 /* SPDX-License-Identifier: MIT */
 /* origin: musl src/math/fmaf.c Ported to generic Rust algorithm in 2025, TG. */
+/* The musl subnormal rounding bug is fixed using the formally proven algorithm from */
+/* "Emulation of FMA and correctly-rounded sums: proved algorithms using rounding to odd" */
+/* by Sylvie Boldo and Guillaume Melquiond, https://guillaume.melquiond.fr/doc/08-tc.pdf */
 
 use crate::support::{
     CastFrom, CastInto, Float, FpResult, IntTy, MinInt, NarrowFloat, Round, Status, WideFloat,
 };
 
-/// Fma implementation when a hardware-backed larger float type is available. For `f32` and `f64`,
-/// `f64` has enough precision to represent the `f32` in its entirety, except for double rounding.
+/// Fma implementation when a hardware-backed larger float type is available.
+/// The larger type has enough precision and exponent range to represent the exact product,
+/// leaving only the addition and the final narrowing susceptible to double rounding.
 #[inline]
 pub fn fma_wide_round<F, B>(x: F, y: F, z: F, round: Round) -> FpResult<F>
 where
@@ -26,10 +30,11 @@ where
     let prec_diff = B::SIG_BITS - F::SIG_BITS;
     let excess_prec = ui & ((one << prec_diff) - one);
     let halfway = one << (prec_diff - 1);
+    let min_normal_exp = (B::EXP_BIAS as i32 + F::EXP_MIN) as u32;
 
     // Common case: the larger precision is fine if...
-    // This is not a halfway case
-    if excess_prec != halfway
+    // This is a normal result and not a halfway case
+    if (re >= min_normal_exp && excess_prec != halfway)
         // Or the result is NaN
         || re == B::EXP_SAT
         // Or the result is exact
@@ -38,11 +43,10 @@ where
         || round != Round::Nearest
     {
         let min_inexact_exp = (B::EXP_BIAS as i32 + F::EXP_MIN_SUBNORM) as u32;
-        let max_inexact_exp = (B::EXP_BIAS as i32 + F::EXP_MIN) as u32;
 
         let mut status = Status::OK;
 
-        if (min_inexact_exp..max_inexact_exp).contains(&re) && status.inexact() {
+        if (min_inexact_exp..min_normal_exp).contains(&re) && status.inexact() {
             // This branch is never hit; requires previous operations to set a status
             status.set_inexact(false);
 
@@ -60,12 +64,20 @@ where
         };
     }
 
+    // FastTwoSum recovers the exact residual of the widened addition. If the addition was inexact
+    // and its rounded significand is even, move it one ULP toward the residual to produce a
+    // round-to-odd intermediate. Theorem 3 proves that rounding this intermediate to nearest in `F`
+    // gives the correctly rounded result, including for subnormals and underflow.
     let neg = ui >> (B::BITS - 1) != IntTy::<B>::ZERO;
     let err = if neg == (zb > xy) {
         xy - result + zb
     } else {
         zb - result + xy
     };
+    // Exact sums need no correction, and odd inexact sums are already round-to-odd.
+    if err == B::ZERO || (ui & one) != IntTy::<B>::ZERO {
+        return FpResult::ok(result.narrow());
+    }
     if neg == (err < B::ZERO) {
         ui += one;
     } else {