From cdbdd1a01a9e6db792a8d249dcd7ec70c0fbac82 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 15:53:18 +0000
Subject: [PATCH 01/11] simd(wasm): mul_add rounds once on F32x16 and F64x8

simd128 has no FMA, and relaxed-simd's relaxed_madd may or may not fuse
depending on the engine. wasm mul_add was mul-then-add (two roundings), or
relaxed_madd under relaxed-simd, so it differed from every other backend.

- F32x16: each lane in f64, where a * b is exact; the f64 sum is rounded to
  odd (two-sum error, then one ulp toward it when inexact and even), so the
  single rounding to f32 is correct, subnormals included (Boldo & Melquiond,
  IEEE TC 2008).
- F64x8: Dekker exact product, two-sum with c, rounded-to-odd sum of the error
  terms, one final rounding (the same paper's emulated FMA). A v128 pair with
  a lane outside 2^-460 <= |a|, |b| <= 2^460, |c| <= 2^960 (or NaN, inf, zero)
  is computed lane by lane with f64::mul_add.

The relaxed-simd branches are gone. Parity group 0x11xx compares both against
the element mul_add bit for bit (NaN matches NaN): random bit patterns over
every exponent, near-cancelling c, ties that double rounding breaks the wrong
way, and lanes that take the f64 fallback; anti-vacuity requires the unfused
result to differ on at least 256 lanes per type.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 crates/simd-masking-parity/src/lib.rs | 119 ++++++++++++++++++++++++-
 src/simd.rs                           |  10 +--
 src/simd_wasm.rs                      | 120 +++++++++++++++++++++-----
 3 files changed, 219 insertions(+), 30 deletions(-)

diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs
index efbce508..c495e02f 100644
--- a/crates/simd-masking-parity/src/lib.rs
+++ b/crates/simd-masking-parity/src/lib.rs
@@ -33,7 +33,8 @@
 //! bytes 12..16); `0xFxx` the 16-lane byte vectors `I8x16` / `U8x16` (wrapping
 //! `add`/`sub`, signed/unsigned `min`/`max`, the load/store round trips); `0x10xx`
 //! the saturating narrows down the 512-bit ladder (`I32x16 → I16x32 → I8x64 /
-//! U8x64`) and the `F32x16` ↔ `I32x16` casts. `main.rs` (native / qemu) and
+//! U8x64`) and the `F32x16` ↔ `I32x16` casts; `0x11xx` `F32x16` / `F64x8`
+//! `mul_add` rounding once, bit for bit against the element `mul_add`. `main.rs` (native / qemu) and
 //! `selfcheck()` (the wasm cdylib export, driven by `run.mjs`) both call
 //! [`run`].
 
@@ -54,7 +55,7 @@ use ndarray::simd::{
 };
 
 /// Number of check groups [`run`] executes (for the log line only).
-pub const CHECKS: usize = 17;
+pub const CHECKS: usize = 18;
 
 /// The wasm export: identical to [`run`], `extern "C"` so `run.mjs` can call it.
 #[no_mangle]
@@ -69,6 +70,7 @@ pub fn run() -> u32 {
         check_mask_algebra, check_care_match, check_masked_reductions, check_blend, check_morton_shift,
         check_predicates_under, check_set_range, check_unsigned_compare_to_mask, check_gather_scatter_group,
         check_f64x8_compare, check_i8x16_u8x16_lanes, check_care_match16, check_narrow_ladder,
+        check_mul_add_single_rounding,
     ];
     for g in groups {
         if let Err(code) = g() {
@@ -1937,3 +1939,116 @@ fn check_narrow_ladder() -> Result<(), u32> {
     }
     Ok(())
 }
+
+/// `F32x16::mul_add` and `F64x8::mul_add` round once: bit-identical to the
+/// element `mul_add` (NaN matches any NaN). The fixtures include ties that
+/// double rounding resolves the wrong way, and products the f64 emulation must
+/// hand to its lane fallback.
+fn check_mul_add_single_rounding() -> Result<(), u32> {
+    let mut rng = SplitMix64(0x1100_0000_0018);
+    let same32 = |x: f32, y: f32| x.to_bits() == y.to_bits() || (x.is_nan() && y.is_nan());
+    let same64 = |x: f64, y: f64| x.to_bits() == y.to_bits() || (x.is_nan() && y.is_nan());
+    let mut unfused_differs = 0usize;
+
+    // f32: random bit patterns over every exponent (subnormals, inf, NaN), plus
+    // a near-cancelling c and a tie broken by a tiny c.
+    let x = 1.0f32 + f32::EPSILON * 2048.0; // 1 + 2^-12: x * x = 1 + 2^-11 + 2^-24
+    for round in 0..4096u32 {
+        let r32 = |r: &mut SplitMix64| f32::from_bits(r.next() as u32);
+        let (a, b, c): ([f32; 16], [f32; 16], [f32; 16]) = match round % 4 {
+            0 => (
+                core::array::from_fn(|_| r32(&mut rng)),
+                core::array::from_fn(|_| r32(&mut rng)),
+                core::array::from_fn(|_| r32(&mut rng)),
+            ),
+            1 => {
+                let a: [f32; 16] =
+                    core::array::from_fn(|_| f32::from_bits((rng.next() as u32 & 0x807F_FFFF) | 0x3F00_0000));
+                let b: [f32; 16] =
+                    core::array::from_fn(|_| f32::from_bits((rng.next() as u32 & 0x807F_FFFF) | 0x3F00_0000));
+                (a, b, core::array::from_fn(|i| -(a[i] * b[i])))
+            }
+            _ => {
+                // a * b = (1 + 2^-11 + 2^-24) * 4^k: a tie in f32 that c (far below
+                // half an ulp, alternating sign) must break.
+                let k = (rng.next() % 100) as i32 - 50;
+                let scale = 2f32.powi(k);
+                let sign = |i: usize| if i % 2 == 0 { 1.0f32 } else { -1.0 };
+                (
+                    [x * scale; 16],
+                    [x * scale; 16],
+                    core::array::from_fn(|i| sign(i) * scale * scale * 2f32.powi(-40 - i as i32)),
+                )
+            }
+        };
+        let got = F32x16::from_array(a)
+            .mul_add(F32x16::from_array(b), F32x16::from_array(c))
+            .to_array();
+        for i in 0..16 {
+            let want = a[i].mul_add(b[i], c[i]);
+            if !same32(got[i], want) {
+                return Err(0x1100 + (round % 4) * 16 + i as u32);
+            }
+            if !same32(a[i] * b[i] + c[i], want) {
+                unfused_differs += 1;
+            }
+        }
+    }
+    if unfused_differs < 256 {
+        return Err(0x11F0);
+    }
+
+    // f64: in-range operands (the vector path), a cancelling c, the same tie
+    // trick, and out-of-range lanes (the lane fallback). In the tie rounds c sits
+    // below the 53 bits after the product's error term pl, so `tl + pl` is
+    // inexact and only its round-to-odd keeps the tie broken.
+    let mut unfused_differs = 0usize;
+    // (1 + 2^-26)(1 + 2^-27) = 1 + 2^-26 + 2^-27 + 2^-53: 2^-53 is half an ulp, a tie.
+    let (ya, yb) = (1.0f64 + f64::EPSILON * 2f64.powi(26), 1.0f64 + f64::EPSILON * 2f64.powi(25));
+    for round in 0..4096u32 {
+        let rin = |r: &mut SplitMix64| {
+            f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | ((1023 - 400 + r.next() % 800) << 52))
+        };
+        let (a, b, c): ([f64; 8], [f64; 8], [f64; 8]) = match round % 4 {
+            0 => (
+                core::array::from_fn(|_| rin(&mut rng)),
+                core::array::from_fn(|_| rin(&mut rng)),
+                core::array::from_fn(|_| rin(&mut rng)),
+            ),
+            1 => {
+                let a: [f64; 8] = core::array::from_fn(|_| rin(&mut rng));
+                let b: [f64; 8] = core::array::from_fn(|_| rin(&mut rng));
+                (a, b, core::array::from_fn(|i| -(a[i] * b[i])))
+            }
+            2 => {
+                let k = (rng.next() % 400) as i32 - 200;
+                (
+                    core::array::from_fn(|_| ya * 2f64.powi(k)),
+                    [yb; 8],
+                    core::array::from_fn(|i| if i % 2 == 0 { 1.0 } else { -1.0 } * 2f64.powi(k - 114 - i as i32)),
+                )
+            }
+            _ => (
+                core::array::from_fn(|_| f64::from_bits(rng.next())),
+                core::array::from_fn(|_| f64::from_bits(rng.next())),
+                core::array::from_fn(|_| f64::from_bits(rng.next())),
+            ),
+        };
+        let got = F64x8::from_array(a)
+            .mul_add(F64x8::from_array(b), F64x8::from_array(c))
+            .to_array();
+        for i in 0..8 {
+            let want = a[i].mul_add(b[i], c[i]);
+            if !same64(got[i], want) {
+                return Err(0x1140 + (round % 4) * 8 + i as u32);
+            }
+            if !same64(a[i] * b[i] + c[i], want) {
+                unfused_differs += 1;
+            }
+        }
+    }
+    if unfused_differs < 256 {
+        return Err(0x11F1);
+    }
+    Ok(())
+}
diff --git a/src/simd.rs b/src/simd.rs
index d7ffdc7b..24d5d1f4 100644
--- a/src/simd.rs
+++ b/src/simd.rs
@@ -572,12 +572,10 @@ pub fn simd_exp_f32(x: F32x16) -> F32x16 {
     // single-constant `x - n * ln2` rounds n * ln2 at |x| ~ 90, where one ulp
     // is ~7.6e-6: that error was larger than the polynomial's.
     //
-    // Plain `*` and `+`, never `mul_add`: wasm's `mul_add` is unfused (or,
-    // with relaxed-simd, implementation-defined) while every other backend
-    // fuses, so `mul_add` here would make the result differ by an ulp
-    // between wasm and native. Unfused ops round identically on every
-    // backend (the same rule as `gemm_f64_tiled` in `simd_ops`), and the
-    // split above makes the fused form unnecessary for accuracy.
+    // Plain `*` and `+`: unfused ops round identically on every backend (the
+    // same rule as `gemm_f64_tiled` in `simd_ops`), and the split above makes
+    // the fused form unnecessary for accuracy. (`mul_add` is fused on every
+    // backend too, wasm included, but costs an emulation there.)
     // floor(t + 0.5), not round(): `round` breaks ties to even on AVX-512,
     // NEON and wasm but away from zero on AVX2, scalar and nightly, so a tie
     // would pick a different n per backend. floor is the same everywhere.
diff --git a/src/simd_wasm.rs b/src/simd_wasm.rs
index 6d535b09..7c1882e2 100644
--- a/src/simd_wasm.rs
+++ b/src/simd_wasm.rs
@@ -34,13 +34,20 @@
 //!
 //! # Semantic notes (cross-backend parity)
 //!
-//! * **FMA.** Base simd128 has no fused multiply-add. `mul_add` uses
-//!   `mul`-then-`add` (two roundings) unless the `relaxed-simd` target
-//!   feature is enabled, in which case `f32x4_relaxed_madd` is used. The
-//!   unfused path differs from the fused scalar/NEON/AVX `mul_add` by ~1 ULP
-//!   in the common case, but can diverge by more under catastrophic
-//!   cancellation. Tests use a tolerance, matching the existing
-//!   cross-backend tests.
+//! * **FMA.** Base simd128 has no fused multiply-add, and `relaxed-simd`'s
+//!   `relaxed_madd` may or may not fuse depending on the engine. `mul_add` is
+//!   therefore emulated and rounds once, like the other backends:
+//!   - `F32x16`: each lane is computed in f64, where `a * b` is exact, and the
+//!     sum `a * b + c` is rounded to odd (two-sum error term, then one ulp
+//!     toward it when the f64 sum is inexact and even). Rounding that to f32
+//!     gives the correctly rounded f32 result (Boldo & Melquiond, "Emulation of
+//!     FMA and correctly rounded sums: proved algorithms using rounding to odd",
+//!     IEEE TC 2008), including for subnormal results.
+//!   - `F64x8`: Dekker's exact product, two-sum with `c`, a rounded-to-odd sum
+//!     of the two error terms, then one rounding (the same paper's emulated
+//!     FMA). It needs no overflow or underflow in the split, so a v128 pair whose
+//!     lanes leave `2^-460 <= |a|, |b| <= 2^460`, `|c| <= 2^960` (or are NaN,
+//!     infinite or zero) is computed lane by lane with `f64::mul_add`.
 //! * **reduce_sum.** Uses a balanced tree reduction (pairwise v128 adds)
 //!   rather than the scalar fallback's sequential left-fold, so the result
 //!   can differ in the last ULP — the same reduction-order tolerance the
@@ -186,18 +193,21 @@ pub mod wasm32_simd {
             Self([f32x4_floor(self.0[0]), f32x4_floor(self.0[1]), f32x4_floor(self.0[2]), f32x4_floor(self.0[3])])
         }
 
+        /// Fused multiply-add, rounded once (see the module notes on FMA).
         #[inline(always)]
         pub fn mul_add(self, b: Self, c: Self) -> Self {
             #[inline(always)]
             fn madd(a: v128, b: v128, c: v128) -> v128 {
-                #[cfg(target_feature = "relaxed-simd")]
-                {
-                    f32x4_relaxed_madd(a, b, c)
-                }
-                #[cfg(not(target_feature = "relaxed-simd"))]
-                {
-                    f32x4_add(f32x4_mul(a, b), c)
-                }
+                // Lanes 2..4 moved down so `promote_low` reaches them.
+                let up = |v: v128| i64x2_shuffle::<1, 1>(v, v);
+                let lo =
+                    fma_f32_in_f64(f64x2_promote_low_f32x4(a), f64x2_promote_low_f32x4(b), f64x2_promote_low_f32x4(c));
+                let hi = fma_f32_in_f64(
+                    f64x2_promote_low_f32x4(up(a)),
+                    f64x2_promote_low_f32x4(up(b)),
+                    f64x2_promote_low_f32x4(up(c)),
+                );
+                i64x2_shuffle::<0, 2>(f32x4_demote_f64x2_zero(lo), f32x4_demote_f64x2_zero(hi))
             }
             Self([
                 madd(self.0[0], b.0[0], c.0[0]),
@@ -535,18 +545,27 @@ pub mod wasm32_simd {
             Self([f64x2_floor(self.0[0]), f64x2_floor(self.0[1]), f64x2_floor(self.0[2]), f64x2_floor(self.0[3])])
         }
 
+        /// Fused multiply-add, rounded once (see the module notes on FMA).
         #[inline(always)]
         pub fn mul_add(self, b: Self, c: Self) -> Self {
             #[inline(always)]
             fn madd(a: v128, b: v128, c: v128) -> v128 {
-                #[cfg(target_feature = "relaxed-simd")]
-                {
-                    f64x2_relaxed_madd(a, b, c)
-                }
-                #[cfg(not(target_feature = "relaxed-simd"))]
-                {
-                    f64x2_add(f64x2_mul(a, b), c)
+                let in_range = |x: v128, lo: f64, hi: f64| {
+                    let m = f64x2_abs(x);
+                    v128_and(f64x2_ge(m, f64x2_splat(lo)), f64x2_le(m, f64x2_splat(hi)))
+                };
+                let ok = v128_and(
+                    v128_and(in_range(a, TWO_M460, TWO_460), in_range(b, TWO_M460, TWO_460)),
+                    f64x2_le(f64x2_abs(c), f64x2_splat(TWO_960)),
+                );
+                if i64x2_all_true(ok) {
+                    return fma_f64_emulated(a, b, c);
                 }
+                let lane = |i: usize| match i {
+                    0 => f64x2_extract_lane::<0>(a).mul_add(f64x2_extract_lane::<0>(b), f64x2_extract_lane::<0>(c)),
+                    _ => f64x2_extract_lane::<1>(a).mul_add(f64x2_extract_lane::<1>(b), f64x2_extract_lane::<1>(c)),
+                };
+                f64x2(lane(0), lane(1))
             }
             Self([
                 madd(self.0[0], b.0[0], c.0[0]),
@@ -2398,6 +2417,63 @@ pub mod wasm32_simd {
         }
     }
 
+    /// `2^-460`, `2^460`, `2^960`: the range in which [`fma_f64_emulated`] is exact.
+    const TWO_M460: f64 = f64::from_bits((1023 - 460) << 52);
+    const TWO_460: f64 = f64::from_bits((1023 + 460) << 52);
+    const TWO_960: f64 = f64::from_bits((1023 + 960) << 52);
+
+    /// Two-sum: `s + e == x + y` exactly, `s = RN(x + y)`, for finite operands.
+    #[inline(always)]
+    fn two_sum(x: v128, y: v128) -> (v128, v128) {
+        let s = f64x2_add(x, y);
+        let yy = f64x2_sub(s, x);
+        let e = f64x2_add(f64x2_sub(x, f64x2_sub(s, yy)), f64x2_sub(y, yy));
+        (s, e)
+    }
+
+    /// Round-to-odd of `x + y`: `RN(x + y)`, moved one ulp toward the exact sum
+    /// when it is inexact and its last mantissa bit is even. Lanes whose sum is
+    /// not finite pass through unchanged.
+    #[inline(always)]
+    fn add_round_to_odd(x: v128, y: v128) -> v128 {
+        let (s, e) = two_sum(x, y);
+        let inexact = f64x2_ne(e, f64x2_splat(0.0));
+        let finite = f64x2_lt(f64x2_abs(s), f64x2_splat(f64::INFINITY));
+        let even = i64x2_eq(v128_and(s, i64x2_splat(1)), i64x2_splat(0));
+        // -1 when `e` and `s` differ in sign (the exact sum is smaller in
+        // magnitude), +1 otherwise; added to the bits it steps the magnitude.
+        let step = v128_or(i64x2_shr(v128_xor(s, e), 63), i64x2_splat(1));
+        i64x2_add(s, v128_and(step, v128_and(v128_and(inexact, finite), even)))
+    }
+
+    /// `a * b + c` for f32 lanes promoted to f64: the product is exact in f64,
+    /// and the sum rounded to odd leaves the single rounding to f32 correct.
+    #[inline(always)]
+    fn fma_f32_in_f64(a: v128, b: v128, c: v128) -> v128 {
+        add_round_to_odd(f64x2_mul(a, b), c)
+    }
+
+    /// Emulated f64 FMA (Boldo & Melquiond 2008): `ph + pl == a * b` (Dekker),
+    /// `th + tl == c + ph` (two-sum), then `RN(th + RO(tl + pl))`. Exact for the
+    /// operand range checked by `F64x8::mul_add`.
+    #[inline(always)]
+    fn fma_f64_emulated(a: v128, b: v128, c: v128) -> v128 {
+        let split = |x: v128| {
+            let t = f64x2_mul(f64x2_splat(134_217_729.0), x); // 2^27 + 1
+            let hi = f64x2_sub(t, f64x2_sub(t, x));
+            (hi, f64x2_sub(x, hi))
+        };
+        let (ah, al) = split(a);
+        let (bh, bl) = split(b);
+        let ph = f64x2_mul(a, b);
+        let pl = f64x2_add(
+            f64x2_add(f64x2_add(f64x2_sub(f64x2_mul(ah, bh), ph), f64x2_mul(ah, bl)), f64x2_mul(al, bh)),
+            f64x2_mul(al, bl),
+        );
+        let (th, tl) = two_sum(c, ph);
+        f64x2_add(th, add_round_to_odd(tl, pl))
+    }
+
     impl I32x16 {
         pub const LANES: usize = 16;
 

From f3a877bc7f1dac6f0213a97dbd3827660121c7b2 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:01:24 +0000
Subject: [PATCH 02/11] parity: is_multiple_of in the mul_add group (clippy)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 crates/simd-masking-parity/src/lib.rs | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs
index c495e02f..bfff4996 100644
--- a/crates/simd-masking-parity/src/lib.rs
+++ b/crates/simd-masking-parity/src/lib.rs
@@ -1973,7 +1973,7 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
                 // half an ulp, alternating sign) must break.
                 let k = (rng.next() % 100) as i32 - 50;
                 let scale = 2f32.powi(k);
-                let sign = |i: usize| if i % 2 == 0 { 1.0f32 } else { -1.0 };
+                let sign = |i: usize| if i.is_multiple_of(2) { 1.0f32 } else { -1.0 };
                 (
                     [x * scale; 16],
                     [x * scale; 16],
@@ -2025,7 +2025,9 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
                 (
                     core::array::from_fn(|_| ya * 2f64.powi(k)),
                     [yb; 8],
-                    core::array::from_fn(|i| if i % 2 == 0 { 1.0 } else { -1.0 } * 2f64.powi(k - 114 - i as i32)),
+                    core::array::from_fn(
+                        |i| if i.is_multiple_of(2) { 1.0 } else { -1.0 } * 2f64.powi(k - 114 - i as i32),
+                    ),
                 )
             }
             _ => (

From a03450e468122ef6504b2a504252f263a27b2373 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:05:38 +0000
Subject: [PATCH 03/11] ci(simd-matrix): skip libm's CPUID-dispatched fma in
 the v2 no-AVX scan

The mul_add parity group calls scalar f32/f64::mul_add as its reference,
which links compiler-builtins' libm fma/fmaf. Those carry fma_with_fma /
fma_with_fma4 variants picked at runtime by CPUID (the Nehalem qemu run
takes fma_fallback and passes), and their vfmadd* lines tripped the
"no AVX instruction" scan. The scan now skips exactly those four
compiler_builtins functions; any other AVX op still fails it (checked
against libc.so.6: 9620 hits remain).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 .github/workflows/simd-matrix.yaml | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/simd-matrix.yaml b/.github/workflows/simd-matrix.yaml
index 7acea754..bc590aab 100644
--- a/.github/workflows/simd-matrix.yaml
+++ b/.github/workflows/simd-matrix.yaml
@@ -313,9 +313,13 @@ jobs:
         run: |
           bin=crates/simd-masking-parity/target/release/simd-masking-parity
           # Nehalem has no AVX: no ymm/zmm or mask register, and no VEX-encoded
-          # (v-prefixed) op on xmm either. The parity binary has no
-          # runtime-gated paths, so any hit came from an ungated intrinsic.
-          bad=$(objdump -d --no-show-raw-insn "$bin" | grep -E '%ymm|%zmm|%k[0-7]|\sv[a-z0-9]+\s.*%xmm' || true)
+          # (v-prefixed) op on xmm either. The one runtime-gated path is
+          # compiler-builtins' libm fma/fmaf (what scalar `mul_add` calls),
+          # whose `*_with_fma` / `*_with_fma4` variants are picked by CPUID and
+          # skipped here; any other hit came from an ungated intrinsic.
+          bad=$(objdump -d --no-show-raw-insn "$bin" \
+            | awk '/>:$/ { skip = ($0 ~ /compiler_builtins.*fmaf?_with_fma4?>:$/) } !skip' \
+            | grep -E '%ymm|%zmm|%k[0-7]|\sv[a-z0-9]+\s.*%xmm' || true)
           if [ -n "$bad" ]; then echo "$bad" | head -20; echo "::error::AVX instruction in the v2 parity binary"; exit 1; fi
           echo "no AVX instruction found"
       - name: SIMD unit tests under qemu -cpu Nehalem

From ebefc2ffa75b1d47862780f15b8797a47f18ec58 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:07:35 +0000
Subject: [PATCH 04/11] ci(simd-matrix): same libm fma skip in the AVX-arm
 no-AVX2/FMA scan

The Sandy Bridge scan hit the same four compiler_builtins fma/fmaf
CPUID-dispatched variants as the v2 scan. Locally: 4 hits before, 0 after
on the avx-qemu parity binary; the filtered scan still reports 4481
AVX2/FMA/F16C lines in libc.so.6.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 .github/workflows/simd-matrix.yaml | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/simd-matrix.yaml b/.github/workflows/simd-matrix.yaml
index bc590aab..add78fdf 100644
--- a/.github/workflows/simd-matrix.yaml
+++ b/.github/workflows/simd-matrix.yaml
@@ -268,9 +268,12 @@ jobs:
         run: |
           bin=crates/simd-masking-parity/target/release/simd-masking-parity
           # Sandy Bridge has AVX but neither AVX2 (integer ops on ymm), FMA3
-          # nor F16C. The parity binary has no runtime-gated paths, so any of
-          # them here came from an intrinsic the AVX arm failed to gate.
-          dis=$(objdump -d --no-show-raw-insn "$bin")
+          # nor F16C. The one runtime-gated path is compiler-builtins' libm
+          # fma/fmaf (what scalar `mul_add` calls), whose `*_with_fma` /
+          # `*_with_fma4` variants are picked by CPUID and skipped here; any
+          # other hit came from an intrinsic the AVX arm failed to gate.
+          dis=$(objdump -d --no-show-raw-insn "$bin" \
+            | awk '/>:$/ { skip = ($0 ~ /compiler_builtins.*fmaf?_with_fma4?>:$/) } !skip')
           bad=$( { echo "$dis" | grep -E '\svp[a-z0-9]+\s.*%ymm' | grep -vE '\svptest\s';
                    echo "$dis" | grep -E '\sv(fn?m(add|sub)|fmaddsub|fmsubadd)[0-9a-z]*\s|\svcvt(ph2ps|ps2ph)\s'; } || true)
           if [ -n "$bad" ]; then echo "$bad" | head -20; echo "::error::AVX2/FMA/F16C instruction in the AVX-arm parity binary"; exit 1; fi

From b21226b419c7ae1158bf45a3c60cc6cc48b561a4 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:09:00 +0000
Subject: [PATCH 05/11] Revert "ci(simd-matrix): same libm fma skip in the
 AVX-arm no-AVX2/FMA scan"

This reverts commit d3dbea2fb0f2d3bc98138f9760b6e95987425ee1.
---
 .github/workflows/simd-matrix.yaml | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/.github/workflows/simd-matrix.yaml b/.github/workflows/simd-matrix.yaml
index add78fdf..bc590aab 100644
--- a/.github/workflows/simd-matrix.yaml
+++ b/.github/workflows/simd-matrix.yaml
@@ -268,12 +268,9 @@ jobs:
         run: |
           bin=crates/simd-masking-parity/target/release/simd-masking-parity
           # Sandy Bridge has AVX but neither AVX2 (integer ops on ymm), FMA3
-          # nor F16C. The one runtime-gated path is compiler-builtins' libm
-          # fma/fmaf (what scalar `mul_add` calls), whose `*_with_fma` /
-          # `*_with_fma4` variants are picked by CPUID and skipped here; any
-          # other hit came from an intrinsic the AVX arm failed to gate.
-          dis=$(objdump -d --no-show-raw-insn "$bin" \
-            | awk '/>:$/ { skip = ($0 ~ /compiler_builtins.*fmaf?_with_fma4?>:$/) } !skip')
+          # nor F16C. The parity binary has no runtime-gated paths, so any of
+          # them here came from an intrinsic the AVX arm failed to gate.
+          dis=$(objdump -d --no-show-raw-insn "$bin")
           bad=$( { echo "$dis" | grep -E '\svp[a-z0-9]+\s.*%ymm' | grep -vE '\svptest\s';
                    echo "$dis" | grep -E '\sv(fn?m(add|sub)|fmaddsub|fmsubadd)[0-9a-z]*\s|\svcvt(ph2ps|ps2ph)\s'; } || true)
           if [ -n "$bad" ]; then echo "$bad" | head -20; echo "::error::AVX2/FMA/F16C instruction in the AVX-arm parity binary"; exit 1; fi

From c6b285d1ecd1a0efb52e22d662141f578ae21220 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:09:00 +0000
Subject: [PATCH 06/11] Revert "ci(simd-matrix): skip libm's CPUID-dispatched
 fma in the v2 no-AVX scan"

This reverts commit b0fd6a658eaaab5098dd2b5b0c58b0398e413ffd.
---
 .github/workflows/simd-matrix.yaml | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/.github/workflows/simd-matrix.yaml b/.github/workflows/simd-matrix.yaml
index bc590aab..7acea754 100644
--- a/.github/workflows/simd-matrix.yaml
+++ b/.github/workflows/simd-matrix.yaml
@@ -313,13 +313,9 @@ jobs:
         run: |
           bin=crates/simd-masking-parity/target/release/simd-masking-parity
           # Nehalem has no AVX: no ymm/zmm or mask register, and no VEX-encoded
-          # (v-prefixed) op on xmm either. The one runtime-gated path is
-          # compiler-builtins' libm fma/fmaf (what scalar `mul_add` calls),
-          # whose `*_with_fma` / `*_with_fma4` variants are picked by CPUID and
-          # skipped here; any other hit came from an ungated intrinsic.
-          bad=$(objdump -d --no-show-raw-insn "$bin" \
-            | awk '/>:$/ { skip = ($0 ~ /compiler_builtins.*fmaf?_with_fma4?>:$/) } !skip' \
-            | grep -E '%ymm|%zmm|%k[0-7]|\sv[a-z0-9]+\s.*%xmm' || true)
+          # (v-prefixed) op on xmm either. The parity binary has no
+          # runtime-gated paths, so any hit came from an ungated intrinsic.
+          bad=$(objdump -d --no-show-raw-insn "$bin" | grep -E '%ymm|%zmm|%k[0-7]|\sv[a-z0-9]+\s.*%xmm' || true)
           if [ -n "$bad" ]; then echo "$bad" | head -20; echo "::error::AVX instruction in the v2 parity binary"; exit 1; fi
           echo "no AVX instruction found"
       - name: SIMD unit tests under qemu -cpu Nehalem

From b27d40159eb1d42aa90e82cf0d1e0ec062a1191e Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:17:37 +0000
Subject: [PATCH 07/11] simd: exact scalar FMA instead of libm on backends
 without an FMA instruction

On targets without a compile-time FMA instruction, f32/f64::mul_add is a
libm call that picks an FMA3, FMA4 or software path by CPUID at run time.
The polyfill selects its realization at compile time, so the per-lane
mul_add of the scalar, AVX-without-FMA and x86-64-v2 backends, and the wasm
F64x8 lanes outside the vector range, now use simd_fma::ExactFma:

- with target_feature = "fma" (x86) or on aarch64: mul_add, which lowers to
  the instruction;
- otherwise f32 through f64 with the sum rounded to odd, and f64 in integer
  arithmetic (exact sum in a 128-bit window with a sticky bit, one rounding).

The 0x11xx parity group now compares against its own exact reference (a
4608-bit integer sum rounded to nearest-even) instead of the element
mul_add, so the parity binary links no libm fma at all. Both the emulation
and the reference are checked against hardware FMA in unit tests
(400k + 300k cases plus all special-value triples).

This replaces the CI-scan exclusions of libm's CPUID-dispatched fma, which
are reverted.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 crates/simd-masking-parity/src/lib.rs | 271 ++++++++++++++++++-
 src/lib.rs                            |   3 +
 src/simd_avx2.rs                      |   4 +-
 src/simd_fma.rs                       | 358 ++++++++++++++++++++++++++
 src/simd_scalar.rs                    |   2 +-
 src/simd_wasm.rs                      |  15 +-
 6 files changed, 643 insertions(+), 10 deletions(-)
 create mode 100644 src/simd_fma.rs

diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs
index bfff4996..e2031cbf 100644
--- a/crates/simd-masking-parity/src/lib.rs
+++ b/crates/simd-masking-parity/src/lib.rs
@@ -34,7 +34,8 @@
 //! `add`/`sub`, signed/unsigned `min`/`max`, the load/store round trips); `0x10xx`
 //! the saturating narrows down the 512-bit ladder (`I32x16 → I16x32 → I8x64 /
 //! U8x64`) and the `F32x16` ↔ `I32x16` casts; `0x11xx` `F32x16` / `F64x8`
-//! `mul_add` rounding once, bit for bit against the element `mul_add`. `main.rs` (native / qemu) and
+//! `mul_add` rounding once, bit for bit against an exact integer reference (no
+//! `mul_add`, so no libm in the binary). `main.rs` (native / qemu) and
 //! `selfcheck()` (the wasm cdylib export, driven by `run.mjs`) both call
 //! [`run`].
 
@@ -1941,7 +1942,7 @@ fn check_narrow_ladder() -> Result<(), u32> {
 }
 
 /// `F32x16::mul_add` and `F64x8::mul_add` round once: bit-identical to the
-/// element `mul_add` (NaN matches any NaN). The fixtures include ties that
+/// exact reference [`oracle_fma`] (NaN matches any NaN). The fixtures include ties that
 /// double rounding resolves the wrong way, and products the f64 emulation must
 /// hand to its lane fallback.
 fn check_mul_add_single_rounding() -> Result<(), u32> {
@@ -1985,7 +1986,7 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
             .mul_add(F32x16::from_array(b), F32x16::from_array(c))
             .to_array();
         for i in 0..16 {
-            let want = a[i].mul_add(b[i], c[i]);
+            let want = oracle_fma_f32(a[i], b[i], c[i]);
             if !same32(got[i], want) {
                 return Err(0x1100 + (round % 4) * 16 + i as u32);
             }
@@ -2040,7 +2041,7 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
             .mul_add(F64x8::from_array(b), F64x8::from_array(c))
             .to_array();
         for i in 0..8 {
-            let want = a[i].mul_add(b[i], c[i]);
+            let want = oracle_fma_f64(a[i], b[i], c[i]);
             if !same64(got[i], want) {
                 return Err(0x1140 + (round % 4) * 8 + i as u32);
             }
@@ -2054,3 +2055,265 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
     }
     Ok(())
 }
+
+/// Exact `a * b + c`, rounded once to nearest-even: the reference for the
+/// `0x11xx` group. It holds the exact sum in a 4608-bit integer (bit 0 is
+/// `2^ORACLE_BASE`, below every product of two subnormals), so it shares no
+/// arithmetic with the realization under test and calls no `mul_add`.
+const ORACLE_WORDS: usize = 72;
+const ORACLE_BASE: i32 = -2300;
+
+/// An IEEE binary format: `p` significand bits, smallest step `2^gmin`,
+/// largest exponent `emax`.
+#[derive(Clone, Copy)]
+struct Format {
+    p: u32,
+    gmin: i32,
+    emax: i32,
+}
+const BINARY32: Format = Format {
+    p: 24,
+    gmin: -149,
+    emax: 127,
+};
+const BINARY64: Format = Format {
+    p: 53,
+    gmin: -1074,
+    emax: 1023,
+};
+
+/// A finite value as `(negative, m, e)` with value `m * 2^e`, from its bits.
+fn oracle_split(bits: u64, f: Format) -> (bool, u128, i32) {
+    let frac_bits = f.p - 1;
+    let total = if f.p == 24 { 32 } else { 64 };
+    let neg = bits >> (total - 1) & 1 == 1;
+    let exp = (bits >> frac_bits) & ((1u64 << (total - 1 - frac_bits)) - 1);
+    let frac = bits & ((1u64 << frac_bits) - 1);
+    if exp == 0 {
+        (neg, frac as u128, f.gmin)
+    } else {
+        (neg, (frac | 1u64 << frac_bits) as u128, exp as i32 - 1 + f.gmin)
+    }
+}
+
+fn oracle_place(m: u128, e: i32) -> [u64; ORACLE_WORDS] {
+    let mut w = [0u64; ORACLE_WORDS];
+    let off = (e - ORACLE_BASE) as usize;
+    let (word, bit) = (off / 64, off % 64);
+    let lo = m << bit;
+    let hi = if bit == 0 { 0 } else { m >> (128 - bit) };
+    w[word] = lo as u64;
+    w[word + 1] = (lo >> 64) as u64;
+    w[word + 2] = hi as u64;
+    w
+}
+
+fn oracle_cmp(x: &[u64; ORACLE_WORDS], y: &[u64; ORACLE_WORDS]) -> core::cmp::Ordering {
+    x.iter().rev().cmp(y.iter().rev())
+}
+
+fn oracle_add(x: &[u64; ORACLE_WORDS], y: &[u64; ORACLE_WORDS]) -> [u64; ORACLE_WORDS] {
+    let mut out = [0u64; ORACLE_WORDS];
+    let mut carry = 0u64;
+    for i in 0..ORACLE_WORDS {
+        let (s1, c1) = x[i].overflowing_add(y[i]);
+        let (s2, c2) = s1.overflowing_add(carry);
+        out[i] = s2;
+        carry = (c1 | c2) as u64;
+    }
+    out
+}
+
+/// `x - y` for `x >= y`.
+fn oracle_sub(x: &[u64; ORACLE_WORDS], y: &[u64; ORACLE_WORDS]) -> [u64; ORACLE_WORDS] {
+    let mut out = [0u64; ORACLE_WORDS];
+    let mut borrow = 0u64;
+    for i in 0..ORACLE_WORDS {
+        let (d1, b1) = x[i].overflowing_sub(y[i]);
+        let (d2, b2) = d1.overflowing_sub(borrow);
+        out[i] = d2;
+        borrow = (b1 | b2) as u64;
+    }
+    out
+}
+
+fn oracle_bit(w: &[u64; ORACLE_WORDS], i: usize) -> bool {
+    w[i / 64] >> (i % 64) & 1 == 1
+}
+
+/// Round the exact magnitude `w * 2^ORACLE_BASE` to the format; returns the
+/// result's bit pattern (sign applied).
+fn oracle_round(neg: bool, w: &[u64; ORACLE_WORDS], f: Format) -> u64 {
+    let total = if f.p == 24 { 32 } else { 64 };
+    let sign = (neg as u64) << (total - 1);
+    let Some(h) = (0..ORACLE_WORDS * 64).rev().find(|&i| oracle_bit(w, i)) else {
+        return sign;
+    };
+    let top = h as i32 + ORACLE_BASE;
+    let g = (top - (f.p as i32 - 1)).max(f.gmin);
+    let shift = (g - ORACLE_BASE) as usize;
+    let mut q: u64 = 0;
+    for i in (shift..=h).rev() {
+        q = q << 1 | oracle_bit(w, i) as u64;
+    }
+    let round = shift > 0 && oracle_bit(w, shift - 1);
+    let sticky = (0..shift.saturating_sub(1)).any(|i| oracle_bit(w, i));
+    if round && (sticky || q & 1 == 1) {
+        q += 1;
+    }
+    let (mut q, mut g) = (q, g);
+    if q == 1u64 << f.p {
+        q >>= 1;
+        g += 1;
+    }
+    let frac_bits = f.p - 1;
+    if q >> frac_bits == 0 {
+        // Subnormal (or zero): only reachable at the smallest step.
+        return sign | q;
+    }
+    let t = g + frac_bits as i32;
+    if t > f.emax {
+        return sign | (((1u64 << (total - 1 - frac_bits)) - 1) << frac_bits);
+    }
+    let biased = (t - (f.gmin + frac_bits as i32) + 1) as u64;
+    sign | biased << frac_bits | (q & ((1u64 << frac_bits) - 1))
+}
+
+/// The oracle on bit patterns; `None` means the exact result is NaN.
+fn oracle_fma(a: u64, b: u64, c: u64, f: Format, value: impl Fn(u64) -> f64) -> Option<u64> {
+    let total = if f.p == 24 { 32 } else { 64 };
+    let (va, vb, vc) = (value(a), value(b), value(c));
+    if va.is_nan() || vb.is_nan() || vc.is_nan() {
+        return None;
+    }
+    let sign_of = |x: u64| x >> (total - 1) & 1 == 1;
+    let sp = sign_of(a) != sign_of(b);
+    let inf = |neg: bool| ((neg as u64) << (total - 1)) | (((1u64 << (total - f.p)) - 1) << (f.p - 1));
+    if va.is_infinite() || vb.is_infinite() {
+        if va == 0.0 || vb == 0.0 || (vc.is_infinite() && sign_of(c) != sp) {
+            return None;
+        }
+        return Some(inf(sp));
+    }
+    if vc.is_infinite() {
+        return Some(c);
+    }
+    if va == 0.0 || vb == 0.0 {
+        if vc == 0.0 {
+            // Round to nearest: -0 only when both terms are -0.
+            return Some(((sp && sign_of(c)) as u64) << (total - 1));
+        }
+        return Some(c);
+    }
+    let (_, ma, ea) = oracle_split(a, f);
+    let (_, mb, eb) = oracle_split(b, f);
+    let prod = oracle_place(ma * mb, ea + eb);
+    if vc == 0.0 {
+        return Some(oracle_round(sp, &prod, f));
+    }
+    let (sc, mc, ec) = oracle_split(c, f);
+    let cw = oracle_place(mc, ec);
+    let (neg, mag) = if sp == sc {
+        (sp, oracle_add(&prod, &cw))
+    } else if oracle_cmp(&prod, &cw) != core::cmp::Ordering::Less {
+        (sp, oracle_sub(&prod, &cw))
+    } else {
+        (sc, oracle_sub(&cw, &prod))
+    };
+    if mag.iter().all(|&x| x == 0) {
+        return Some(0);
+    }
+    Some(oracle_round(neg, &mag, f))
+}
+
+fn oracle_fma_f32(a: f32, b: f32, c: f32) -> f32 {
+    let v = |x: u64| f32::from_bits(x as u32) as f64;
+    match oracle_fma(a.to_bits() as u64, b.to_bits() as u64, c.to_bits() as u64, BINARY32, v) {
+        Some(bits) => f32::from_bits(bits as u32),
+        None => f32::NAN,
+    }
+}
+
+fn oracle_fma_f64(a: f64, b: f64, c: f64) -> f64 {
+    match oracle_fma(a.to_bits(), b.to_bits(), c.to_bits(), BINARY64, f64::from_bits) {
+        Some(bits) => f64::from_bits(bits),
+        None => f64::NAN,
+    }
+}
+
+#[cfg(test)]
+mod oracle_tests {
+    use super::*;
+
+    /// The oracle against the host's `mul_add` (a hardware FMA on the hosts
+    /// these tests run on): random patterns, cancellation, ties, the
+    /// subnormal and overflow edges, and every pair of special values.
+    #[test]
+    fn oracle_matches_hardware_fma() {
+        let mut rng = SplitMix64(0x0AC1_E000_0001);
+        let specials64 = [
+            0.0,
+            -0.0,
+            1.0,
+            -1.0,
+            f64::MIN_POSITIVE,
+            5e-324,
+            -5e-324,
+            f64::MAX,
+            -f64::MAX,
+            f64::INFINITY,
+            f64::NEG_INFINITY,
+            f64::NAN,
+            1.5,
+        ];
+        for &a in &specials64 {
+            for &b in &specials64 {
+                for &c in &specials64 {
+                    let (w, g) = (a.mul_add(b, c), oracle_fma_f64(a, b, c));
+                    assert!(
+                        w.to_bits() == g.to_bits() || (w.is_nan() && g.is_nan()),
+                        "f64 {a:e} {b:e} {c:e}: {g:e} vs {w:e}"
+                    );
+                    let (x, y, z) = (a as f32, b as f32, c as f32);
+                    let (w, g) = (x.mul_add(y, z), oracle_fma_f32(x, y, z));
+                    assert!(
+                        w.to_bits() == g.to_bits() || (w.is_nan() && g.is_nan()),
+                        "f32 {x:e} {y:e} {z:e}: {g:e} vs {w:e}"
+                    );
+                }
+            }
+        }
+        for round in 0..300_000u32 {
+            let (a, b, c) = match round % 3 {
+                0 => (f64::from_bits(rng.next()), f64::from_bits(rng.next()), f64::from_bits(rng.next())),
+                1 => {
+                    let ea = 1 + rng.next() % 2046;
+                    let t = if rng.next() & 1 == 0 {
+                        rng.next() % 70
+                    } else {
+                        1990 + rng.next() % 70
+                    };
+                    let eb = (t + 1023).saturating_sub(ea).clamp(1, 2046);
+                    let a = f64::from_bits((rng.next() & 0x800F_FFFF_FFFF_FFFF) | ea << 52);
+                    let b = f64::from_bits((rng.next() & 0x800F_FFFF_FFFF_FFFF) | eb << 52);
+                    (a, b, f64::from_bits(rng.next() & 0x801F_FFFF_FFFF_FFFF))
+                }
+                _ => {
+                    let a = f64::from_bits(rng.next());
+                    let b = f64::from_bits(rng.next());
+                    (a, b, -(a * b))
+                }
+            };
+            let (w, g) = (a.mul_add(b, c), oracle_fma_f64(a, b, c));
+            assert!(w.to_bits() == g.to_bits() || (w.is_nan() && g.is_nan()), "f64 {a:e} {b:e} {c:e}: {g:e} vs {w:e}");
+            let (x, y) = (f32::from_bits(rng.next() as u32), f32::from_bits(rng.next() as u32));
+            let z = if round % 2 == 0 {
+                f32::from_bits(rng.next() as u32)
+            } else {
+                -(x * y)
+            };
+            let (w, g) = (x.mul_add(y, z), oracle_fma_f32(x, y, z));
+            assert!(w.to_bits() == g.to_bits() || (w.is_nan() && g.is_nan()), "f32 {x:e} {y:e} {z:e}: {g:e} vs {w:e}");
+        }
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
index 482cb9af..60cd0a1b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -486,6 +486,9 @@ pub mod simd_ops;
 #[allow(clippy::all, missing_docs, dead_code, unused_variables, unused_imports)]
 pub mod simd_half;
 
+/// Correctly rounded scalar FMA for backends without an FMA instruction.
+pub(crate) mod simd_fma;
+
 /// Runtime SIMD dispatch — release-binary distribution path.
 ///
 /// Gated under `--features runtime-dispatch`. Mutually exclusive with
diff --git a/src/simd_avx2.rs b/src/simd_avx2.rs
index d76a2674..102e0c18 100644
--- a/src/simd_avx2.rs
+++ b/src/simd_avx2.rs
@@ -605,7 +605,7 @@ impl F32x16 {
         let ca = c.to_array();
         let mut o = [0.0f32; 16];
         for i in 0..16 {
-            o[i] = a[i].mul_add(ba[i], ca[i]);
+            o[i] = crate::simd_fma::ExactFma::exact_fma(a[i], ba[i], ca[i]);
         }
         Self::from_array(o)
     }
@@ -948,7 +948,7 @@ impl F64x8 {
         let ca = c.to_array();
         let mut o = [0.0f64; 8];
         for i in 0..8 {
-            o[i] = a[i].mul_add(ba[i], ca[i]);
+            o[i] = crate::simd_fma::ExactFma::exact_fma(a[i], ba[i], ca[i]);
         }
         Self::from_array(o)
     }
diff --git a/src/simd_fma.rs b/src/simd_fma.rs
new file mode 100644
index 00000000..9fa26ef2
--- /dev/null
+++ b/src/simd_fma.rs
@@ -0,0 +1,358 @@
+//! Correctly rounded scalar fused multiply-add for the backends that have no
+//! FMA instruction at compile time.
+//!
+//! `f32::mul_add` / `f64::mul_add` on such a target is a libm call that picks
+//! an FMA3, FMA4 or software path by CPUID at run time. The polyfill selects
+//! its realization at compile time, so the per-lane `mul_add` of the scalar,
+//! AVX-without-FMA and x86-64-v2 backends (and the wasm `F64x8` lanes outside
+//! the vector path's range) use [`ExactFma`] instead:
+//!
+//! * with an FMA instruction (`target_feature = "fma"` on x86, always on
+//!   aarch64) it is `mul_add`, which lowers to that instruction;
+//! * otherwise `f32` is computed in f64, where the product is exact, with the
+//!   sum rounded to odd before the final rounding (Boldo & Melquiond, IEEE TC
+//!   2008), and `f64` in integer arithmetic.
+//!
+//! Both round once, to nearest-even, for every input including subnormals,
+//! zeros, infinities and NaN (any NaN result is a NaN, payload unspecified).
+
+// With an FMA instruction only the tests reach the emulation.
+#![cfg_attr(
+    any(
+        all(any(target_arch = "x86", target_arch = "x86_64"), target_feature = "fma"),
+        target_arch = "aarch64"
+    ),
+    allow(dead_code)
+)]
+
+/// `a * b + c` rounded once.
+pub(crate) trait ExactFma: Copy {
+    fn exact_fma(self, b: Self, c: Self) -> Self;
+}
+
+#[cfg(any(
+    all(any(target_arch = "x86", target_arch = "x86_64"), target_feature = "fma"),
+    target_arch = "aarch64"
+))]
+impl ExactFma for f32 {
+    #[inline(always)]
+    fn exact_fma(self, b: Self, c: Self) -> Self {
+        self.mul_add(b, c)
+    }
+}
+
+#[cfg(any(
+    all(any(target_arch = "x86", target_arch = "x86_64"), target_feature = "fma"),
+    target_arch = "aarch64"
+))]
+impl ExactFma for f64 {
+    #[inline(always)]
+    fn exact_fma(self, b: Self, c: Self) -> Self {
+        self.mul_add(b, c)
+    }
+}
+
+#[cfg(not(any(
+    all(any(target_arch = "x86", target_arch = "x86_64"), target_feature = "fma"),
+    target_arch = "aarch64"
+)))]
+impl ExactFma for f32 {
+    #[inline(always)]
+    fn exact_fma(self, b: Self, c: Self) -> Self {
+        fma_f32(self, b, c)
+    }
+}
+
+#[cfg(not(any(
+    all(any(target_arch = "x86", target_arch = "x86_64"), target_feature = "fma"),
+    target_arch = "aarch64"
+)))]
+impl ExactFma for f64 {
+    #[inline(always)]
+    fn exact_fma(self, b: Self, c: Self) -> Self {
+        fma_f64(self, b, c)
+    }
+}
+
+/// f32 FMA through f64: the product of two f32 is exact in f64, and the sum
+/// rounded to odd makes the single rounding to f32 correct.
+#[inline]
+pub(crate) fn fma_f32(a: f32, b: f32, c: f32) -> f32 {
+    let p = a as f64 * b as f64;
+    let c = c as f64;
+    let s = p + c;
+    // Two-sum error of p + c (exact for finite operands).
+    let cc = s - p;
+    let e = (p - (s - cc)) + (c - cc);
+    let mut bits = s.to_bits();
+    if e != 0.0 && s.is_finite() && bits & 1 == 0 {
+        // One ulp toward the exact sum: up in magnitude when e has the sign of s.
+        if (e < 0.0) == (s < 0.0) {
+            bits += 1;
+        } else {
+            bits -= 1;
+        }
+    }
+    f64::from_bits(bits) as f32
+}
+
+/// `x = (-1)^neg * m * 2^e` for a finite nonzero f64.
+#[inline]
+fn split_f64(x: f64) -> (bool, u64, i32) {
+    let bits = x.to_bits();
+    let exp = ((bits >> 52) & 0x7ff) as i32;
+    let frac = bits & ((1u64 << 52) - 1);
+    if exp == 0 {
+        (bits >> 63 != 0, frac, -1074)
+    } else {
+        (bits >> 63 != 0, frac | (1u64 << 52), exp - 1075)
+    }
+}
+
+/// `m * 2^shift` as a u128 scaled to `2^base`, with every bit shifted out
+/// OR-ed into bit 0 (sticky).
+#[inline]
+fn align(m: u128, exp: i32, base: i32) -> u128 {
+    let d = exp - base;
+    if d >= 0 {
+        m << d
+    } else {
+        let k = (-d) as u32;
+        if k >= 128 {
+            (m != 0) as u128
+        } else {
+            (m >> k) | ((m & ((1u128 << k) - 1) != 0) as u128)
+        }
+    }
+}
+
+/// `2^e` for `-1074 <= e <= 1023`.
+#[inline]
+fn pow2(e: i32) -> f64 {
+    if e >= -1022 {
+        f64::from_bits(((e + 1023) as u64) << 52)
+    } else {
+        f64::from_bits(1u64 << (e + 1074))
+    }
+}
+
+/// f64 FMA in integer arithmetic: the exact sum held in a 128-bit window
+/// with a sticky bit, then one round to nearest-even.
+#[inline]
+pub(crate) fn fma_f64(a: f64, b: f64, c: f64) -> f64 {
+    if !a.is_finite() || !b.is_finite() || a == 0.0 || b == 0.0 {
+        // The product is an exact zero, an infinity or NaN; one rounding.
+        return a * b + c;
+    }
+    if !c.is_finite() {
+        // A finite product never turns an infinite or NaN c into anything else.
+        return c + c;
+    }
+    if c == 0.0 {
+        // a*b is RN(a*b); adding ±0 to it would lose the sign of an underflow.
+        return a * b;
+    }
+    let (sa, ma, ea) = split_f64(a);
+    let (sb, mb, eb) = split_f64(b);
+    let (sc, mc, ec) = split_f64(c);
+    let sp = sa != sb;
+    let p = ma as u128 * mb as u128;
+    let ep = ea + eb;
+    let lp = 128 - p.leading_zeros() as i32;
+    let lc = 64 - mc.leading_zeros() as i32;
+    // The larger term ends at bit 126 of the window, so it is never shifted
+    // right and a carry still fits; only a term far below it gets a sticky bit.
+    let top = (ep + lp).max(ec + lc);
+    let base = top - 126;
+    let pa = align(p, ep, base);
+    let ca = align(mc as u128, ec, base);
+    let (neg, s) = if sp == sc {
+        (sp, pa + ca)
+    } else if pa >= ca {
+        (sp, pa - ca)
+    } else {
+        (sc, ca - pa)
+    };
+    if s == 0 {
+        return 0.0;
+    }
+    let len = 128 - s.leading_zeros() as i32;
+    // Keep 53 bits, or fewer where the result is subnormal (granularity 2^-1074).
+    let shift = (len - 53).max(-1074 - base);
+    let (q, e2) = if shift <= 0 {
+        (s << (-shift) as u32, base + shift)
+    } else {
+        let k = shift as u32;
+        let q = if k > 128 {
+            // s < 2^128 <= 2^(k-1): below half the smallest step.
+            0
+        } else {
+            let (q, rem, half) = if k == 128 {
+                (0u128, s, 1u128 << 127)
+            } else {
+                (s >> k, s & ((1u128 << k) - 1), 1u128 << (k - 1))
+            };
+            if rem > half || (rem == half && q & 1 == 1) {
+                q + 1
+            } else {
+                q
+            }
+        };
+        (q, base + shift)
+    };
+    // q <= 2^53 is exact in f64 and q * 2^e2 is representable or overflows.
+    let mag = if e2 > 1023 {
+        q as f64 * pow2(1023) * pow2((e2 - 1023).min(1023))
+    } else {
+        q as f64 * pow2(e2)
+    };
+    if neg {
+        -mag
+    } else {
+        mag
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    struct Rng(u64);
+    impl Rng {
+        fn next(&mut self) -> u64 {
+            self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
+            let mut z = self.0;
+            z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+            z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+            z ^ (z >> 31)
+        }
+    }
+
+    fn same32(x: f32, y: f32) -> bool {
+        x.to_bits() == y.to_bits() || (x.is_nan() && y.is_nan())
+    }
+    fn same64(x: f64, y: f64) -> bool {
+        x.to_bits() == y.to_bits() || (x.is_nan() && y.is_nan())
+    }
+
+    /// Against the platform `mul_add` (hardware FMA on the CI hosts), over
+    /// random bit patterns, near-cancelling c, ties and subnormal results.
+    #[test]
+    fn emulated_fma_matches_mul_add() {
+        let mut r = Rng(0xF3A_0001);
+        let mut differs_unfused = 0usize;
+        for round in 0..200_000u32 {
+            let (a, b, c) = match round % 4 {
+                0 => {
+                    (f32::from_bits(r.next() as u32), f32::from_bits(r.next() as u32), f32::from_bits(r.next() as u32))
+                }
+                1 => {
+                    let a = f32::from_bits((r.next() as u32 & 0x807F_FFFF) | 0x3F00_0000);
+                    let b = f32::from_bits((r.next() as u32 & 0x807F_FFFF) | 0x3F00_0000);
+                    (a, b, -(a * b))
+                }
+                2 => {
+                    let x = 1.0f32 + f32::EPSILON * 2048.0;
+                    let s = 2f32.powi((r.next() % 100) as i32 - 50);
+                    (
+                        x * s,
+                        x * s,
+                        s * s * 2f32.powi(-40 - (r.next() % 16) as i32) * if r.next() & 1 == 0 { 1.0 } else { -1.0 },
+                    )
+                }
+                _ => {
+                    // Results in and around the subnormal range.
+                    let a = f32::from_bits((r.next() as u32 & 0x807F_FFFF) | (((r.next() % 40) as u32 + 40) << 23));
+                    let b = f32::from_bits((r.next() as u32 & 0x807F_FFFF) | (((r.next() % 40) as u32 + 20) << 23));
+                    (a, b, f32::from_bits(r.next() as u32 & 0x80FF_FFFF))
+                }
+            };
+            let want = a.mul_add(b, c);
+            assert!(
+                same32(fma_f32(a, b, c), want),
+                "f32 fma({a:e}, {b:e}, {c:e}) = {:e}, want {want:e}",
+                fma_f32(a, b, c)
+            );
+            if !same32(a * b + c, want) {
+                differs_unfused += 1;
+            }
+        }
+        assert!(differs_unfused > 1000, "fixtures must separate fused from unfused ({differs_unfused})");
+
+        let mut differs_unfused = 0usize;
+        let (ya, yb) = (1.0f64 + f64::EPSILON * 2f64.powi(26), 1.0f64 + f64::EPSILON * 2f64.powi(25));
+        for round in 0..200_000u32 {
+            let (a, b, c) = match round % 5 {
+                0 => (f64::from_bits(r.next()), f64::from_bits(r.next()), f64::from_bits(r.next())),
+                1 => {
+                    let a = f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | 0x3FE0_0000_0000_0000);
+                    let b = f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | 0x3FE0_0000_0000_0000);
+                    (a, b, -(a * b))
+                }
+                2 => {
+                    let k = (r.next() % 1200) as i32 - 600;
+                    let sign = if r.next() & 1 == 0 { 1.0 } else { -1.0 };
+                    (ya * 2f64.powi(k), yb, sign * 2f64.powi(k - 60 - (r.next() % 80) as i32))
+                }
+                3 => {
+                    // Products whose exponent lands near the subnormal range or overflow.
+                    let ea = 1 + r.next() % 2046;
+                    let target = if r.next() & 1 == 0 {
+                        r.next() % 70
+                    } else {
+                        1990 + r.next() % 70
+                    };
+                    let eb = (target + 1023).saturating_sub(ea).clamp(1, 2046);
+                    let a = f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | (ea << 52));
+                    let b = f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | (eb << 52));
+                    (a, b, f64::from_bits(r.next() & 0x801F_FFFF_FFFF_FFFF))
+                }
+                _ => {
+                    let e = |r: &mut Rng| (r.next() % 2047) << 52;
+                    let a = f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | e(&mut r));
+                    let b = f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | e(&mut r));
+                    (a, b, -(a * b) * (1.0 + f64::EPSILON))
+                }
+            };
+            let want = a.mul_add(b, c);
+            assert!(
+                same64(fma_f64(a, b, c), want),
+                "f64 fma({a:e}, {b:e}, {c:e}) = {:e}, want {want:e}",
+                fma_f64(a, b, c)
+            );
+            if !same64(a * b + c, want) {
+                differs_unfused += 1;
+            }
+        }
+        assert!(differs_unfused > 1000, "fixtures must separate fused from unfused ({differs_unfused})");
+    }
+
+    #[test]
+    fn emulated_fma_special_values() {
+        let v64 = [
+            0.0,
+            -0.0,
+            1.0,
+            -1.0,
+            f64::MIN_POSITIVE,
+            5e-324,
+            -5e-324,
+            f64::MAX,
+            -f64::MAX,
+            f64::INFINITY,
+            f64::NEG_INFINITY,
+            f64::NAN,
+            1.5,
+            3.0,
+        ];
+        for &a in &v64 {
+            for &b in &v64 {
+                for &c in &v64 {
+                    assert!(same64(fma_f64(a, b, c), a.mul_add(b, c)), "f64 fma({a:e}, {b:e}, {c:e})");
+                    let (x, y, z) = (a as f32, b as f32, c as f32);
+                    assert!(same32(fma_f32(x, y, z), x.mul_add(y, z)), "f32 fma({x:e}, {y:e}, {z:e})");
+                }
+            }
+        }
+    }
+}
diff --git a/src/simd_scalar.rs b/src/simd_scalar.rs
index c34569d9..a2461acb 100644
--- a/src/simd_scalar.rs
+++ b/src/simd_scalar.rs
@@ -114,7 +114,7 @@ macro_rules! impl_float_type {
             pub fn mul_add(self, b: Self, c: Self) -> Self {
                 let mut out = [0.0 as $elem; $lanes];
                 for i in 0..$lanes {
-                    out[i] = self.0[i].mul_add(b.0[i], c.0[i]);
+                    out[i] = crate::simd_fma::ExactFma::exact_fma(self.0[i], b.0[i], c.0[i]);
                 }
                 Self(out)
             }
diff --git a/src/simd_wasm.rs b/src/simd_wasm.rs
index 7c1882e2..925cdfac 100644
--- a/src/simd_wasm.rs
+++ b/src/simd_wasm.rs
@@ -47,7 +47,8 @@
 //!     of the two error terms, then one rounding (the same paper's emulated
 //!     FMA). It needs no overflow or underflow in the split, so a v128 pair whose
 //!     lanes leave `2^-460 <= |a|, |b| <= 2^460`, `|c| <= 2^960` (or are NaN,
-//!     infinite or zero) is computed lane by lane with `f64::mul_add`.
+//!     infinite or zero) is computed lane by lane with the integer FMA in
+//!     `simd_fma.rs` (no libm call, so no run-time CPU dispatch).
 //! * **reduce_sum.** Uses a balanced tree reduction (pairwise v128 adds)
 //!   rather than the scalar fallback's sequential left-fold, so the result
 //!   can differ in the last ULP — the same reduction-order tolerance the
@@ -562,8 +563,16 @@ pub mod wasm32_simd {
                     return fma_f64_emulated(a, b, c);
                 }
                 let lane = |i: usize| match i {
-                    0 => f64x2_extract_lane::<0>(a).mul_add(f64x2_extract_lane::<0>(b), f64x2_extract_lane::<0>(c)),
-                    _ => f64x2_extract_lane::<1>(a).mul_add(f64x2_extract_lane::<1>(b), f64x2_extract_lane::<1>(c)),
+                    0 => crate::simd_fma::fma_f64(
+                        f64x2_extract_lane::<0>(a),
+                        f64x2_extract_lane::<0>(b),
+                        f64x2_extract_lane::<0>(c),
+                    ),
+                    _ => crate::simd_fma::fma_f64(
+                        f64x2_extract_lane::<1>(a),
+                        f64x2_extract_lane::<1>(b),
+                        f64x2_extract_lane::<1>(c),
+                    ),
                 };
                 f64x2(lane(0), lane(1))
             }

From 1723eddde717250275ef6e747b591e3b4f80ca2e Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:20:52 +0000
Subject: [PATCH 08/11] simd: FMA instruction on AVX2+FMA halves; scalar
 mul_add_f32/f64 in simd_ops

- simd_avx2 F32x16/F64x8::mul_add were a per-lane scalar loop even on
  AVX2 + FMA. They now call the halves' F32x8/F64x4::mul_add:
  _mm256_fmadd_ps/_pd with target_feature = "fma", the exact scalar FMA per
  lane without it. F64x4 gains mul_add; F32x8's non-FMA branch uses the
  exact scalar FMA instead of f32::mul_add.
- simd_ops gains mul_add_f32 / mul_add_f64 (re-exported from ndarray::simd):
  the polyfill's single-rounding scalar FMA, so tails and consumers (wasm
  included) do not hand-roll one or fall back to f32::mul_add. The add_mul
  and gemm tails use it; the gemm doc no longer claims wasm is unfused.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 src/simd.rs        |  4 +--
 src/simd_avx2.rs   | 22 ++++---------
 src/simd_avx512.rs | 25 ++++++++++++--
 src/simd_ops.rs    | 82 ++++++++++++++++++++++++++++++++--------------
 4 files changed, 87 insertions(+), 46 deletions(-)

diff --git a/src/simd.rs b/src/simd.rs
index 24d5d1f4..a84f5e6d 100644
--- a/src/simd.rs
+++ b/src/simd.rs
@@ -993,8 +993,8 @@ pub use crate::bitwise::popcount_batch_u64;
 #[cfg(feature = "std")]
 pub use crate::simd_ops::{
     add_f32, add_f32_inplace, add_f64, add_f64_inplace, add_mul_f32, add_mul_f64, add_scalar_f32, bf16_tile_gemm_16x16,
-    div_f32, div_f32_inplace, mul_f32, mul_f32_inplace, mul_f64, scale_f32, scale_f32_inplace, sub_f32,
-    sub_f32_inplace,
+    div_f32, div_f32_inplace, mul_add_f32, mul_add_f64, mul_f32, mul_f32_inplace, mul_f64, scale_f32,
+    scale_f32_inplace, sub_f32, sub_f32_inplace,
 };
 
 // ChaCha20 keystream is NO LONGER an `ndarray::simd` surface: the AdaWorldAPI
diff --git a/src/simd_avx2.rs b/src/simd_avx2.rs
index 102e0c18..ea1fb91d 100644
--- a/src/simd_avx2.rs
+++ b/src/simd_avx2.rs
@@ -599,15 +599,10 @@ impl F32x16 {
         Self::from_array(o)
     }
     #[inline(always)]
+    /// Fused multiply-add, single rounding: `_mm256_fmadd_ps` on each half
+    /// with FMA3, the exact scalar FMA per lane without it.
     pub fn mul_add(self, b: Self, c: Self) -> Self {
-        let a = self.to_array();
-        let ba = b.to_array();
-        let ca = c.to_array();
-        let mut o = [0.0f32; 16];
-        for i in 0..16 {
-            o[i] = crate::simd_fma::ExactFma::exact_fma(a[i], ba[i], ca[i]);
-        }
-        Self::from_array(o)
+        Self(self.0.mul_add(b.0, c.0), self.1.mul_add(b.1, c.1))
     }
     #[inline(always)]
     pub fn simd_min(self, other: Self) -> Self {
@@ -942,15 +937,10 @@ impl F64x8 {
         Self::from_array(o)
     }
     #[inline(always)]
+    /// Fused multiply-add, single rounding: `_mm256_fmadd_pd` on each half
+    /// with FMA3, the exact scalar FMA per lane without it.
     pub fn mul_add(self, b: Self, c: Self) -> Self {
-        let a = self.to_array();
-        let ba = b.to_array();
-        let ca = c.to_array();
-        let mut o = [0.0f64; 8];
-        for i in 0..8 {
-            o[i] = crate::simd_fma::ExactFma::exact_fma(a[i], ba[i], ca[i]);
-        }
-        Self::from_array(o)
+        Self(self.0.mul_add(b.0, c.0), self.1.mul_add(b.1, c.1))
     }
     #[inline(always)]
     pub fn simd_min(self, other: Self) -> Self {
diff --git a/src/simd_avx512.rs b/src/simd_avx512.rs
index fa389a2f..794d9014 100644
--- a/src/simd_avx512.rs
+++ b/src/simd_avx512.rs
@@ -1440,12 +1440,13 @@ impl F32x8 {
             // the binary may run on has it (`cpu_guard` checks at startup).
             Self(unsafe { _mm256_fmadd_ps(self.0, a.0, b.0) })
         }
-        // AVX without FMA (the `simd_avx.rs` arm, e.g. Sandy Bridge): a fused
-        // multiply-add per lane, same single rounding as the FMA3 form.
+        // AVX without FMA (the `simd_avx.rs` arm, e.g. Sandy Bridge): the exact
+        // scalar FMA per lane (`simd_fma`), same single rounding as the FMA3 form.
         #[cfg(not(target_feature = "fma"))]
         {
+            use crate::simd_fma::ExactFma;
             let (x, y, z) = (self.to_array(), a.to_array(), b.to_array());
-            Self::from_array(core::array::from_fn(|i| x[i].mul_add(y[i], z[i])))
+            Self::from_array(core::array::from_fn(|i| x[i].exact_fma(y[i], z[i])))
         }
     }
 
@@ -2682,6 +2683,24 @@ impl F64x4 {
             Self(_mm256_and_pd(self.0, mask))
         }
     }
+
+    /// Fused multiply-add: `self * b + c`, single rounding (`_mm256_fmadd_pd`).
+    #[inline(always)]
+    pub fn mul_add(self, b: Self, c: Self) -> Self {
+        #[cfg(target_feature = "fma")]
+        {
+            // SAFETY: `_mm256_fmadd_pd` needs FMA3, which `cfg(target_feature
+            // = "fma")` proves is enabled for this whole build.
+            Self(unsafe { _mm256_fmadd_pd(self.0, b.0, c.0) })
+        }
+        // AVX without FMA: the exact scalar FMA per lane, one rounding.
+        #[cfg(not(target_feature = "fma"))]
+        {
+            use crate::simd_fma::ExactFma;
+            let (x, y, z) = (self.to_array(), b.to_array(), c.to_array());
+            Self::from_array(core::array::from_fn(|i| x[i].exact_fma(y[i], z[i])))
+        }
+    }
 }
 
 impl Add for F64x4 {
diff --git a/src/simd_ops.rs b/src/simd_ops.rs
index fde4a5af..cfd990e2 100644
--- a/src/simd_ops.rs
+++ b/src/simd_ops.rs
@@ -130,19 +130,52 @@ pub fn scale_f32_inplace(a: &mut [f32], scalar: f32) {
 // `add_mul_*` is the slice-level FMA-into-accumulator: `acc[i] += a[i] * b[i]`.
 // Built directly on `F32x16::mul_add` / `F64x8::mul_add`, the polyfilled
 // FMA primitives on the SIMD types (AVX-512: vfmadd231ps; AVX2 + FMA:
-// _mm256_fmadd_ps; NEON: vfmaq_f32; scalar: f32::mul_add → also fused
-// on FMA-capable hosts). Single rounding step, no intermediate
-// `a[i] * b[i]` temporary in scalar emulation — same semantic as
-// BLAS-1 `axpy` with a vector multiplier, and the dominant inner-loop
-// shape in the BLAS-graph GEMM kernels.
+// _mm256_fmadd_ps; NEON: vfmaq_f32; wasm simd128: exact emulation; scalar
+// and targets without FMA: the exact scalar FMA), and on [`mul_add_f32`] /
+// [`mul_add_f64`] for the tails. Single rounding step on every backend —
+// same semantic as BLAS-1 `axpy` with a vector multiplier, and the dominant
+// inner-loop shape in the BLAS-graph GEMM kernels.
+
+/// Scalar fused multiply-add `a * b + c`, rounded once, on every backend.
+///
+/// The FMA instruction where the target has one at compile time
+/// (`target_feature = "fma"` on x86, aarch64); otherwise an exact emulation
+/// (no libm call, no run-time CPU dispatch). Use this instead of
+/// `f32::mul_add` in code that must behave the same on every polyfill arm.
+///
+/// # Examples
+///
+/// ```
+/// use ndarray::simd::mul_add_f32;
+/// let x = 1.0f32 + f32::EPSILON * 2048.0; // x * x = 1 + 2^-11 + 2^-24 exactly
+/// assert_eq!(mul_add_f32(x, x, -(1.0 + f32::EPSILON * 4096.0)), f32::EPSILON / 2.0);
+/// ```
+#[inline(always)]
+pub fn mul_add_f32(a: f32, b: f32, c: f32) -> f32 {
+    crate::simd_fma::ExactFma::exact_fma(a, b, c)
+}
+
+/// f64 sibling of [`mul_add_f32`].
+///
+/// # Examples
+///
+/// ```
+/// use ndarray::simd::mul_add_f64;
+/// let x = 1.0f64 + f64::EPSILON * 2f64.powi(25); // x * x = 1 + 2^-26 + 2^-54
+/// assert_eq!(mul_add_f64(x, x, -(1.0 + f64::EPSILON * 2f64.powi(26))), 2f64.powi(-54));
+/// ```
+#[inline(always)]
+pub fn mul_add_f64(a: f64, b: f64, c: f64) -> f64 {
+    crate::simd_fma::ExactFma::exact_fma(a, b, c)
+}
 
 /// Fused multiply-add into accumulator: `acc[i] += a[i] * b[i]`.
 ///
 /// Operates on the prefix `min(acc.len(), a.len(), b.len())` lanes.
 /// Each lane is computed with a single rounding step via the polyfilled
-/// `F32x16::mul_add` (hardware FMA where available). On scalar fallback
-/// targets, `f32::mul_add` is also fused if the host has FMA — never
-/// less precise than the manual `acc[i] += a[i] * b[i]` two-step.
+/// `F32x16::mul_add` (hardware FMA where available), and [`mul_add_f32`]
+/// for the tail — never less precise than the manual
+/// `acc[i] += a[i] * b[i]` two-step.
 ///
 /// # Examples
 ///
@@ -166,7 +199,7 @@ pub fn add_mul_f32(acc: &mut [f32], a: &[f32], b: &[f32]) {
         i += 16;
     }
     while i < n {
-        acc[i] = a[i].mul_add(b[i], acc[i]);
+        acc[i] = mul_add_f32(a[i], b[i], acc[i]);
         i += 1;
     }
 }
@@ -174,7 +207,7 @@ pub fn add_mul_f32(acc: &mut [f32], a: &[f32], b: &[f32]) {
 /// Fused multiply-add into accumulator (f64): `acc[i] += a[i] * b[i]`.
 ///
 /// f64 sibling of [`add_mul_f32`]. Uses `F64x8::mul_add` (8-wide on
-/// AVX-512 / 4-wide on AVX2-FMA / 2-wide on NEON / scalar `f64::mul_add`).
+/// AVX-512 / 4-wide on AVX2-FMA / 2-wide on NEON / [`mul_add_f64`] otherwise).
 #[inline]
 pub fn add_mul_f64(acc: &mut [f64], a: &[f64], b: &[f64]) {
     let n = acc.len().min(a.len()).min(b.len());
@@ -187,7 +220,7 @@ pub fn add_mul_f64(acc: &mut [f64], a: &[f64], b: &[f64]) {
         i += 8;
     }
     while i < n {
-        acc[i] = a[i].mul_add(b[i], acc[i]);
+        acc[i] = mul_add_f64(a[i], b[i], acc[i]);
         i += 1;
     }
 }
@@ -959,21 +992,16 @@ pub fn gemm_f64_tiled(
 /// Fast tier of [`gemm_f64_tiled`]: same tiling, same ascending-`p`
 /// per-element order, but each step is a **fused** multiply-add
 /// (`c = fma(alpha·A[i,p], B[p,j], c)` — one rounding per step via
-/// `F64x8::mul_add`; fused semantics on every backend, though the
-/// lowering varies: native `vfmadd` on the AVX-512 arm, `vfmaq_f64` on
-/// NEON, per-lane `f64::mul_add` on the AVX2 polyfill and the scalar
-/// tail — which can become a libm `fma()` call on targets built
-/// without the `fma` feature: correct, but slow. Prefer the reference
-/// tier, or an FMA-capable target pin, in that situation).
+/// `F64x8::mul_add`, and [`mul_add_f64`] for the scalar tail; fused
+/// semantics on every backend, though the lowering varies: native
+/// `vfmadd` on AVX-512 and AVX2 + FMA, `vfmaq_f64` on NEON, an exact
+/// emulation on wasm and on targets built without the `fma` feature:
+/// correct, but slower. Prefer the reference tier, or an FMA-capable
+/// target pin, in that situation).
 ///
 /// Semantics relative to the reference tier:
-/// - **Deterministic per (build, runtime)**, but NOT bit-stable across
-///   backends: fusion differs (WASM without `relaxed-simd` has no
-///   fused vector `mul_add`, so its vector lanes match the unfused
-///   reference while the scalar tail is still fused; WASM WITH
-///   `relaxed-simd` uses `f64x2_relaxed_madd`, whose fusion is
-///   implementation-defined — the same wasm binary may round
-///   differently on different runtimes).
+/// - **Bit-identical across backends**: every backend's `mul_add`
+///   rounds once (simd-masking-parity group `0x11xx`).
 /// - **Bit-identical to [`gemm_f64_tiled`]** whenever every
 ///   `a_val·b` product and running sum is exactly representable —
 ///   e.g. integer-valued data with products and accumulation inside
@@ -1092,7 +1120,11 @@ fn gemm_f64_tiled_impl<const FMA: bool>(
                         }
                         for (t, tv) in tail.iter_mut().enumerate().take(ntail) {
                             let bv = b[b_row + nvec * LANES + t];
-                            *tv = if FMA { a_val.mul_add(bv, *tv) } else { *tv + a_val * bv };
+                            *tv = if FMA {
+                                mul_add_f64(a_val, bv, *tv)
+                            } else {
+                                *tv + a_val * bv
+                            };
                         }
                     }
                     for (v, av) in acc.iter().enumerate().take(nvec) {

From a119298a850888db30206dcaddac7d0105c2f739 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:25:40 +0000
Subject: [PATCH 09/11] simd_fma: gate on std like every module that uses it

The no-std builds (tests/1.99.0 without std, nostd/thumbv6m, blas-msrv)
compiled simd_fma while every caller (simd, simd_ops, simd_wasm,
simd_avx512) is std-gated, so -D warnings failed on dead code.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 src/lib.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/lib.rs b/src/lib.rs
index 60cd0a1b..d8cc364f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -487,6 +487,7 @@ pub mod simd_ops;
 pub mod simd_half;
 
 /// Correctly rounded scalar FMA for backends without an FMA instruction.
+#[cfg(feature = "std")]
 pub(crate) mod simd_fma;
 
 /// Runtime SIMD dispatch — release-binary distribution path.

From efdc7d3d894fe657a9850d8d07d82a3c8a44f994 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:34:32 +0000
Subject: [PATCH 10/11] hpc, simd_half: scalar mul_add through
 simd::mul_add_f32/f64

Nine scalar f32/f64::mul_add calls in hpc (randomized_signature,
signature_pde, renderer) and the two lane loops of BF16x16/F16x16::fma
called the std method, which on a target without the fma feature is a
call into compiler_builtins' run-time-dispatched fma. They now go
through the polyfill's exact scalar FMA. Results are unchanged: both
paths are the correctly rounded a*b+c. The BF16/F16 double rounding
(f32 FMA, then narrowing) is untouched.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 src/hpc/randomized_signature.rs | 4 ++--
 src/hpc/renderer.rs             | 8 ++++----
 src/hpc/signature_pde.rs        | 8 ++++----
 src/simd_half.rs                | 4 ++--
 4 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/hpc/randomized_signature.rs b/src/hpc/randomized_signature.rs
index f7465a86..e9962abf 100644
--- a/src/hpc/randomized_signature.rs
+++ b/src/hpc/randomized_signature.rs
@@ -109,7 +109,7 @@ fn dot_simd(row: &[f64], z: &[f64]) -> f64 {
     }
     let mut sum = acc.reduce_sum();
     while c < row.len() {
-        sum = row[c].mul_add(z[c], sum);
+        sum = crate::simd::mul_add_f64(row[c], z[c], sum);
         c += 1;
     }
     sum
@@ -130,7 +130,7 @@ fn axpy_simd(scale: f64, src: &[f64], out: &mut [f64]) {
         r += LANES;
     }
     while r < src.len() {
-        out[r] = scale.mul_add(src[r], out[r]);
+        out[r] = crate::simd::mul_add_f64(scale, src[r], out[r]);
         r += 1;
     }
 }
diff --git a/src/hpc/renderer.rs b/src/hpc/renderer.rs
index 9122e975..d9ddd60c 100644
--- a/src/hpc/renderer.rs
+++ b/src/hpc/renderer.rs
@@ -29,7 +29,7 @@
 //! | x86 AVX2         | 8         | `_mm256_fmadd_ps`  |
 //! | x86 AMX          | 16+tile   | `_tile_dpbf16ps`   |
 //! | aarch64 NEON     | 4         | `vfmaq_f32`        |
-//! | scalar fallback  | 16 (loop) | `f32::mul_add`     |
+//! | no FMA (scalar)  | 16 (loop) | exact `simd_fma`   |
 //!
 //! Consumer writes `crate::simd::F32x16`. The polyfill picks the path.
 //!
@@ -338,12 +338,12 @@ pub fn apply_uniform_force(velocities: &mut [f32], force: [f32; 3], dt: f32) {
             let _ = (f_v, dt_v);
             for k in 0..16 {
                 let idx = (i + k) * POSITION_DIMS;
-                velocities[idx] = force[0].mul_add(dt, velocities[idx]);
+                velocities[idx] = crate::simd::mul_add_f32(force[0], dt, velocities[idx]);
             }
         } else {
             for k in 0..(n_nodes - i) {
                 let idx = (i + k) * POSITION_DIMS;
-                velocities[idx] = force[0].mul_add(dt, velocities[idx]);
+                velocities[idx] = crate::simd::mul_add_f32(force[0], dt, velocities[idx]);
             }
         }
     }
@@ -351,7 +351,7 @@ pub fn apply_uniform_force(velocities: &mut [f32], force: [f32; 3], dt: f32) {
     for axis in 1..POSITION_DIMS {
         for n in 0..n_nodes {
             let idx = n * POSITION_DIMS + axis;
-            velocities[idx] = force[axis].mul_add(dt, velocities[idx]);
+            velocities[idx] = crate::simd::mul_add_f32(force[axis], dt, velocities[idx]);
         }
     }
 }
diff --git a/src/hpc/signature_pde.rs b/src/hpc/signature_pde.rs
index 990746b0..c3d38bbf 100644
--- a/src/hpc/signature_pde.rs
+++ b/src/hpc/signature_pde.rs
@@ -150,11 +150,11 @@ pub fn signature_pde_sweep(x: &[Vec<f64>], y: &[Vec<f64>]) -> f64 {
             debug_assert!(r < m - 1, "signature_pde_sweep: dyr scalar index out of range");
             let mut c = 0.0f64;
             for a in 0..dim {
-                c = dx[a][i - 1].mul_add(dyr[a][r], c);
+                c = crate::simd::mul_add_f64(dx[a][i - 1], dyr[a][r], c);
             }
-            let t = 1.0f64.mul_add(prev1[i], prev1[i - 1]);
-            let u = (-1.0f64).mul_add(prev2[i - 1], t);
-            cur[i] = c.mul_add(prev2[i - 1], u);
+            let t = crate::simd::mul_add_f64(1.0, prev1[i], prev1[i - 1]);
+            let u = crate::simd::mul_add_f64(-1.0, prev2[i - 1], t);
+            cur[i] = crate::simd::mul_add_f64(c, prev2[i - 1], u);
             i += 1;
         }
         std::mem::swap(&mut prev2, &mut prev1);
diff --git a/src/simd_half.rs b/src/simd_half.rs
index a6121c46..ea723666 100644
--- a/src/simd_half.rs
+++ b/src/simd_half.rs
@@ -96,7 +96,7 @@ impl BF16x16 {
             let av = BF16(self.0[i]).to_f32();
             let bv = BF16(b.0[i]).to_f32();
             let cv = BF16(c.0[i]).to_f32();
-            out[i] = BF16::from_f32_rounded(av.mul_add(bv, cv)).0;
+            out[i] = BF16::from_f32_rounded(crate::simd::mul_add_f32(av, bv, cv)).0;
         }
         BF16x16(out)
     }
@@ -219,7 +219,7 @@ impl F16x16 {
             let av = F16(self.0[i]).to_f32();
             let bv = F16(b.0[i]).to_f32();
             let cv = F16(c.0[i]).to_f32();
-            out[i] = F16::from_f32_rounded(av.mul_add(bv, cv)).0;
+            out[i] = F16::from_f32_rounded(crate::simd::mul_add_f32(av, bv, cv)).0;
         }
         F16x16(out)
     }

From b41760d01420149280a4ef575dac015d2c7f396a Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sat, 24 Oct 2026 16:49:27 +0000
Subject: [PATCH 11/11] simd-masking-parity: exact f64 midpoints in the mul_add
 group

The f64 tie rounds always added a tiny c that broke the tie, so the
final rounding never saw an exact midpoint: a ties-away final rounding
in simd_fma::fma_f64 passed the whole group on the v2 arm (it was caught
only by simd_fma's own unit test). Round 4 now cancels the product's
middle terms so a*b + c is exactly halfway between two f64 neighbours,
with the even neighbour below on some lanes and above on others, in both
signs, over 2^-900..2^900.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
---
 crates/simd-masking-parity/src/lib.rs | 31 ++++++++++++++++++++++-----
 1 file changed, 26 insertions(+), 5 deletions(-)

diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs
index e2031cbf..52ba1bea 100644
--- a/crates/simd-masking-parity/src/lib.rs
+++ b/crates/simd-masking-parity/src/lib.rs
@@ -2002,15 +2002,17 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
     // f64: in-range operands (the vector path), a cancelling c, the same tie
     // trick, and out-of-range lanes (the lane fallback). In the tie rounds c sits
     // below the 53 bits after the product's error term pl, so `tl + pl` is
-    // inexact and only its round-to-odd keeps the tie broken.
+    // inexact and only its round-to-odd keeps the tie broken. Round 4 makes the
+    // exact sum itself a midpoint, so the final rounding must pick the even
+    // neighbour (ties-away or ties-toward-zero fail it).
     let mut unfused_differs = 0usize;
     // (1 + 2^-26)(1 + 2^-27) = 1 + 2^-26 + 2^-27 + 2^-53: 2^-53 is half an ulp, a tie.
     let (ya, yb) = (1.0f64 + f64::EPSILON * 2f64.powi(26), 1.0f64 + f64::EPSILON * 2f64.powi(25));
-    for round in 0..4096u32 {
+    for round in 0..5120u32 {
         let rin = |r: &mut SplitMix64| {
             f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | ((1023 - 400 + r.next() % 800) << 52))
         };
-        let (a, b, c): ([f64; 8], [f64; 8], [f64; 8]) = match round % 4 {
+        let (a, b, c): ([f64; 8], [f64; 8], [f64; 8]) = match round % 5 {
             0 => (
                 core::array::from_fn(|_| rin(&mut rng)),
                 core::array::from_fn(|_| rin(&mut rng)),
@@ -2031,11 +2033,30 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
                     ),
                 )
             }
-            _ => (
+            3 => (
                 core::array::from_fn(|_| f64::from_bits(rng.next())),
                 core::array::from_fn(|_| f64::from_bits(rng.next())),
                 core::array::from_fn(|_| f64::from_bits(rng.next())),
             ),
+            _ => {
+                // a * b = (1 + 2^-26 + 2^-27 + 2^-53) * 2^k; each c cancels the middle
+                // terms so the exact sum is a midpoint between two f64 neighbours:
+                //   lane 0: 1 + 2^-53                  -> 1            (even below)
+                //   lane 1: 1 + 2^-52 + 2^-53          -> 1 + 2^-51    (even above)
+                //   lane 2: 1 + 2^-26 + 2^-53          -> 1 + 2^-26    (even below)
+                //   lane 3: 1 + 2^-26 + 2^-52 + 2^-53  -> 1 + 2^-26 + 2^-51 (even above)
+                // Lanes 4..7 are the same sums negated.
+                let k = (rng.next() % 1800) as i32 - 900;
+                let s = 2f64.powi(k);
+                let (e26, e27, e52) = (2f64.powi(-26), 2f64.powi(-27), 2f64.powi(-52));
+                let cut = [-(e26 + e27), -(e26 + e27 - e52), -e27, -(e27 - e52)];
+                let neg = |i: usize| if i < 4 { 1.0 } else { -1.0 };
+                (
+                    core::array::from_fn(|i| neg(i) * ya * s),
+                    [yb; 8],
+                    core::array::from_fn(|i| neg(i) * cut[i % 4] * s),
+                )
+            }
         };
         let got = F64x8::from_array(a)
             .mul_add(F64x8::from_array(b), F64x8::from_array(c))
@@ -2043,7 +2064,7 @@ fn check_mul_add_single_rounding() -> Result<(), u32> {
         for i in 0..8 {
             let want = oracle_fma_f64(a[i], b[i], c[i]);
             if !same64(got[i], want) {
-                return Err(0x1140 + (round % 4) * 8 + i as u32);
+                return Err(0x1140 + (round % 5) * 8 + i as u32);
             }
             if !same64(a[i] * b[i] + c[i], want) {
                 unfused_differs += 1;