diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs index efbce508..52ba1bea 100644 --- a/crates/simd-masking-parity/src/lib.rs +++ b/crates/simd-masking-parity/src/lib.rs @@ -33,7 +33,9 @@ //! 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 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`]. @@ -54,7 +56,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 +71,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 +1940,401 @@ fn check_narrow_ladder() -> Result<(), u32> { } Ok(()) } + +/// `F32x16::mul_add` and `F64x8::mul_add` round once: bit-identical to the +/// 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> { + 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.is_multiple_of(2) { 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 = oracle_fma_f32(a[i], 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. 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..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 % 5 { + 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.is_multiple_of(2) { 1.0 } else { -1.0 } * 2f64.powi(k - 114 - i as i32), + ), + ) + } + 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)) + .to_array(); + 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 % 5) * 8 + i as u32); + } + if !same64(a[i] * b[i] + c[i], want) { + unfused_differs += 1; + } + } + } + if unfused_differs < 256 { + return Err(0x11F1); + } + 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 { + 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/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], y: &[Vec]) -> 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/lib.rs b/src/lib.rs index 482cb9af..d8cc364f 100644 --- a/src/lib.rs +++ b/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 diff --git a/src/simd.rs b/src/simd.rs index d7ffdc7b..a84f5e6d 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. @@ -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 diff --git a/src/simd_avx2.rs b/src/simd_avx2.rs index d76a2674..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] = 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 { 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_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_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) } 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( } 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) { 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 6d535b09..925cdfac 100644 --- a/src/simd_wasm.rs +++ b/src/simd_wasm.rs @@ -34,13 +34,21 @@ //! //! # 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 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 @@ -186,18 +194,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 +546,35 @@ 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 => 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)) } Self([ madd(self.0[0], b.0[0], c.0[0]), @@ -2398,6 +2426,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;