simd: mul_add rounds once on every backend (wasm emulation, exact scalar FMA, no libm)

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;
     }
 }

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,20 +338,20 @@ 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]);
             }
         }
     }
     // Axes 1 (Y), 2 (Z)
     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]);
         }
     }
 }

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);

src/lib.rs

@@ -486,6 +486,10 @@ 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.
+#[cfg(feature = "std")]
+pub(crate) mod simd_fma;
+
 /// Runtime SIMD dispatch — release-binary distribution path.
 ///
 /// Gated under `--features runtime-dispatch`. Mutually exclusive with

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.
@@ -995,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

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] = a[i].mul_add(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] = a[i].mul_add(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 {

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 {