Software FMA for f32 returns wrong result for subnormals on x86_64

[Shnatsel]

libm's software FMA disagrees with hardware FMA result for rounding subnormals.

Run this to reproduce it on x86_64, and with libm depended on with default-features=false to get the software FMA rather than the hardware one:

fn main() {
    let a = std::hint::black_box(f32::from_bits(0x9700_0800));
    let b = std::hint::black_box(f32::from_bits(0x1cff_f001));
    let c = std::hint::black_box(f32::from_bits(0x0001_0002));

    let expected_bits = 0x0001_0001;
    let software_bits = libm::fmaf(a, b, c).to_bits();

    println!("a                 = {:#010x}", a.to_bits());
    println!("b                 = {:#010x}", b.to_bits());
    println!("c                 = {:#010x}", c.to_bits());
    println!("libm::fmaf        = {software_bits:#010x}");
    println!("correctly rounded = {expected_bits:#010x}");

    #[cfg(target_arch = "x86_64")]
    if std::is_x86_feature_detected!("fma") {
        // SAFETY: guarded by runtime FMA feature detection.
        let hardware_bits = unsafe { hardware_fmaf(a, b, c) }.to_bits();
        println!("x86 hardware FMA  = {hardware_bits:#010x}");
        assert_eq!(hardware_bits, expected_bits);
    }

    assert_eq!(software_bits, expected_bits);
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "fma")]
unsafe fn hardware_fmaf(a: f32, b: f32, c: f32) -> f32 {
    use core::arch::x86_64::{_mm_cvtss_f32, _mm_fmadd_ss, _mm_set_ss};

    _mm_cvtss_f32(_mm_fmadd_ss(_mm_set_ss(a), _mm_set_ss(b), _mm_set_ss(c)))
}

The output is:

a                 = 0x97000800
b                 = 0x1cfff001
c                 = 0x00010002
libm::fmaf        = 0x00010002
correctly rounded = 0x00010001
x86 hardware FMA  = 0x00010001

thread 'main' (365592) panicked at src/main.rs:23:5:
assertion `left == right` failed
  left: 65538
 right: 65537

This is happening on x86_64, so it has SSE2 and proper f64 math.

As far as I can tell, the root cause of the issue is this:

compiler-builtins/libm/src/math/generic/fma_wide.rs

Lines 26 to 28 in 7a3101d

	let prec_diff = B::SIG_BITS - F::SIG_BITS;
	let excess_prec = ui & ((one << prec_diff) - one);
	let halfway = one << (prec_diff - 1);

This detection of points exactly halfway between representable ones only works for normal values, but not for subnormals.

This 2008 paper describes correct handling of subnormals, and their algorithm comes with a Coq proof: https://guillaume.melquiond.fr/doc/08-tc.pdf
A Rust implementation of this algorithm handles this case correctly for me in SIMD code without FMA (SSE4.2): linebender/fearless_simd#323

[Shnatsel]

This bug is also present in musl libc from which this implementation was ported: https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c?id=f21a96538f78fa8e2040831b4209b35f2fb581da

[Shnatsel]

Upstream bug report for musl: https://www.openwall.com/lists/musl/2026/08/10/1