diff --git a/crates/wasm-simd-parity/src/lib.rs b/crates/wasm-simd-parity/src/lib.rs index 80ec5c90..c4dc637d 100644 --- a/crates/wasm-simd-parity/src/lib.rs +++ b/crates/wasm-simd-parity/src/lib.rs @@ -35,9 +35,85 @@ pub extern "C" fn selfcheck() -> u32 { if let Err(code) = check_i32x16_compare() { return code; } + if let Err(code) = check_exp_and_round() { + return code; + } 0 } +/// `simd_exp_f32` and `F32x16::round` under wasm. exp must equal, bit for +/// bit, the same unfused op sequence in scalar f32; the native unit test +/// `simd_exp_f32_matches_the_unfused_scalar_sequence_bit_for_bit` pins the +/// native backends to that sequence, so together they give wasm == native. +/// `round` must break ties to even on every backend. +fn check_exp_and_round() -> Result<(), u32> { + // Constants of `simd_exp_f32` (private in ndarray::simd), restated. + const OVERFLOW_X: f32 = 88.722_84; + const LN2_HI: f32 = 0.693_359_4; + const LN2_LO: f32 = -2.121_944_4e-4; + const TWO_POW_NEG_64: f32 = 5.421_011e-20; + fn reference(x: f32) -> f32 { + if x.is_nan() { + return f32::NAN; + } + let xc = x.max(-110.0).min(OVERFLOW_X); + let n = (xc * (1.0 / core::f32::consts::LN_2) + 0.5).floor(); + let r = xc - n * LN2_HI; + let r = r - n * LN2_LO; + let poly = 1.0 + r * (1.0 + r * (0.5 + r * (1.0 / 6.0 + r * (1.0 / 24.0 + r * (1.0 / 120.0))))); + let (n, post) = if n > 127.0 { + (n - 1.0, 2.0) + } else if n < -126.0 { + (n + 64.0, TWO_POW_NEG_64) + } else { + (n, 1.0) + }; + let scale = f32::from_bits((((n as i32).clamp(-126, 127) + 127) as u32) << 23); + poly * scale * post + } + // -110 .. 88.72 in steps of 0.01, then ties of x/ln2 + 0.5 and specials. + let mut xs: Vec = (0..=19872).map(|j| -110.0 + j as f32 * 0.01).collect(); + xs.extend((-150..=128).map(|k| (k as f32 + 0.5) * core::f32::consts::LN_2)); + xs.extend([f32::INFINITY, f32::NEG_INFINITY, f32::NAN, OVERFLOW_X, -0.0, 0.0, f32::MIN_POSITIVE]); + while xs.len() % 16 != 0 { + xs.push(1.0); + } + for chunk in xs.chunks_exact(16) { + let got = ndarray::simd::simd_exp_f32(F32x16::from_slice(chunk)).to_array(); + for (&x, g) in chunk.iter().zip(got) { + let want = reference(x); + if g.to_bits() != want.to_bits() && !(g.is_nan() && want.is_nan()) { + return Err(70); + } + } + } + let ties: [f32; 16] = [ + 0.5, + 1.5, + 2.5, + -0.5, + -1.5, + -2.5, + 3.5, + 4.5, + -3.5, + 0.25, + 0.75, + -0.0, + 1e7 + 0.5, + 8388607.5, + -8388607.5, + 0.0, + ]; + let r = F32x16::from_array(ties).round().to_array(); + for k in 0..16 { + if r[k].to_bits() != ties[k].round_ties_even().to_bits() { + return Err(71); + } + } + Ok(()) +} + /// `U32x16` ARX triple (Add / BitXor / rotate_left) — the ChaCha20/BLAKE lane. fn check_u32x16() -> Result<(), u32> { let a_arr: [u32; 16] = [ diff --git a/src/hpc/gpt2/inference.rs b/src/hpc/gpt2/inference.rs index 1b44329e..2b5ac2cf 100644 --- a/src/hpc/gpt2/inference.rs +++ b/src/hpc/gpt2/inference.rs @@ -453,7 +453,9 @@ fn gelu_simd(x: &mut [f32]) { // tanh approximation via exp: tanh(x) = (exp(2x) - 1) / (exp(2x) + 1) let two_inner = inner + inner; let exp_2x = crate::simd::simd_exp_f32(two_inner); - let tanh_v = (exp_2x - one) / (exp_2x + one); + // 1 - 2/(e+1) equals (e-1)/(e+1) but stays 1 when exp overflows to + // +inf, where the quotient form gives inf/inf = NaN. + let tanh_v = one - (one + one) / (exp_2x + one); let result = v * half * (one + tanh_v); result.copy_to_slice(&mut x[off..off + 16]); } @@ -522,6 +524,20 @@ fn matmul_vec_simd(input: &[f32], weight: &[f32], bias: &[f32], output: &mut [f3 mod tests { use super::*; + /// FAILS IF: large inputs, where exp(2·inner) overflows to +inf, come + /// out NaN from the SIMD body instead of x (positive) or ~0 (negative). + #[test] + fn gelu_stays_finite_where_exp_overflows() { + for v in [12.0f32, 20.0, 1e3, -12.0, -20.0, -1e3] { + let mut x = [v; 17]; + gelu_simd(&mut x); + for (i, &g) in x.iter().enumerate() { + let want = if v > 0.0 { v } else { 0.0 }; + assert!((g - want).abs() <= 1e-4 * v.abs(), "gelu({v})[{i}] = {g}"); + } + } + } + #[test] fn test_layer_norm_identity() { let mut x = vec![1.0, 2.0, 3.0, 4.0]; diff --git a/src/hpc/models/layers.rs b/src/hpc/models/layers.rs index c7f28e1a..54504283 100644 --- a/src/hpc/models/layers.rs +++ b/src/hpc/models/layers.rs @@ -111,7 +111,9 @@ pub fn gelu(x: &mut [f32]) { let inner = sqrt_2_over_pi * (v + coeff * v3); let two_inner = inner + inner; let exp_2x = crate::simd::simd_exp_f32(two_inner); - let tanh_v = (exp_2x - one) / (exp_2x + one); + // 1 - 2/(e+1) equals (e-1)/(e+1) but stays 1 when exp overflows to + // +inf, where the quotient form gives inf/inf = NaN. + let tanh_v = one - (one + one) / (exp_2x + one); let result = v * half * (one + tanh_v); result.copy_to_slice(&mut x[off..off + 16]); } @@ -283,6 +285,20 @@ pub fn rope_apply(q: &mut [f32], k: &mut [f32], head_dim: usize, position: usize mod tests { use super::*; + /// FAILS IF: large inputs, where exp(2·inner) overflows to +inf, come + /// out NaN from the SIMD body instead of x (positive) or ~0 (negative). + #[test] + fn gelu_stays_finite_where_exp_overflows() { + for v in [12.0f32, 20.0, 1e3, -12.0, -20.0, -1e3] { + let mut x = [v; 17]; + gelu(&mut x); + for (i, &g) in x.iter().enumerate() { + let want = if v > 0.0 { v } else { 0.0 }; + assert!((g - want).abs() <= 1e-4 * v.abs(), "gelu({v})[{i}] = {g}"); + } + } + } + #[test] fn test_layer_norm_zero_mean() { let mut x = vec![1.0, 2.0, 3.0, 4.0]; diff --git a/src/hpc/vml.rs b/src/hpc/vml.rs index 1d9513b0..00c77e7a 100644 --- a/src/hpc/vml.rs +++ b/src/hpc/vml.rs @@ -243,7 +243,7 @@ fn vstanh_slice(x: &[f32], out: &mut [f32]) { // tanh(x) = 2·sigmoid(2x) - 1 let two_x = v * two; // sigmoid(z) = 1/(1+exp(-z)) - let neg_two_x = F32x16::splat(0.0) - two_x; + let neg_two_x = -two_x; let exp_neg = simd_exp_f32(neg_two_x); let sigmoid = one / (exp_neg + one); let tanh_v = sigmoid.mul_add(two, F32x16::splat(-1.0)); @@ -273,11 +273,9 @@ fn vsceil_slice(x: &[f32], out: &mut [f32]) { for i in 0..chunks { let off = i * 16; let v = F32x16::from_slice(&x[off..off + 16]); - // ceil = -floor(-x) - let neg = F32x16::splat(0.0) - v; - let floored = neg.floor(); - let ceiled = F32x16::splat(0.0) - floored; - ceiled.copy_to_slice(&mut out[off..off + 16]); + // ceil = -floor(-x). `Neg` flips the sign bit; `0.0 - v` would turn + // ceil(-0.5) and ceil(-0.0) into +0.0. + (-(-v).floor()).copy_to_slice(&mut out[off..off + 16]); } for i in (chunks * 16)..n { out[i] = x[i].ceil(); @@ -300,11 +298,11 @@ fn vsround_slice(x: &[f32], out: &mut [f32]) { fn vsneg_slice(x: &[f32], out: &mut [f32]) { let n = x.len().min(out.len()); let chunks = n / 16; - let zero = F32x16::splat(0.0); for i in 0..chunks { let off = i * 16; let v = F32x16::from_slice(&x[off..off + 16]); - (zero - v).copy_to_slice(&mut out[off..off + 16]); + // `Neg` flips the sign bit, so -(+0.0) is -0.0; `0.0 - v` gives +0.0. + (-v).copy_to_slice(&mut out[off..off + 16]); } for i in (chunks * 16)..n { out[i] = -x[i]; @@ -317,14 +315,12 @@ fn vstrunc_slice(x: &[f32], out: &mut [f32]) { for i in 0..chunks { let off = i * 16; let v = F32x16::from_slice(&x[off..off + 16]); - let abs_v = v.abs(); - let floored = abs_v.floor(); - let zero = F32x16::splat(0.0); - let mask = v.simd_lt(zero); // true where x < 0 - let pos_result = floored; - let neg_result = zero - floored; - let result = mask.select(neg_result, pos_result); - result.copy_to_slice(&mut out[off..off + 16]); + // Negative lanes take -floor(-x), the rest floor(x). -0.0 is not + // < 0, so it takes floor(-0.0) = -0.0 and keeps its sign. + let negative = v.simd_lt(F32x16::splat(0.0)); + negative + .select(-(-v).floor(), v.floor()) + .copy_to_slice(&mut out[off..off + 16]); } for i in (chunks * 16)..n { out[i] = x[i].trunc(); @@ -688,7 +684,7 @@ pub fn vstrunc(x: ArrayView<'_, f32, D>, mut out: ArrayViewMut<'_, #[cfg(test)] mod tests { use super::*; - use crate::{arr1, arr2, s, Array1, Array2}; + use crate::{arr1, arr2, s, Array1, Array2, ArrayView1, ArrayViewMut1}; // ----------------------------------------------------------------------- // vsexp @@ -1215,4 +1211,112 @@ mod tests { vstrunc(strided, out.view_mut()); assert_eq!(out.as_slice().unwrap(), &[1.0, -1.0, 0.0, 0.0]); } + + // ----------------------------------------------------------------------- + // IEEE edge cases through the 16-lane SIMD body, not just the scalar tail. + // Every length from 0 to 40 puts each value both in a SIMD chunk and in + // the tail, so a body that disagrees with its own tail cannot pass. + // ----------------------------------------------------------------------- + + const EDGES: [f32; 12] = [ + 0.0, + -0.0, + 0.5, + -0.5, + 1.5, + -1.5, + -0.25, + f32::INFINITY, + f32::NEG_INFINITY, + f32::MIN_POSITIVE, + -f32::MIN_POSITIVE, + f32::NAN, + ]; + + /// Checks `kernel` against `scalar` bit for bit (NaN only as NaN) for + /// every length 0..=40, cycling through `EDGES`. + fn assert_bitwise_like_scalar( + name: &str, kernel: fn(ArrayView1<'_, f32>, ArrayViewMut1<'_, f32>), scalar: fn(f32) -> f32, + ) { + for n in 0..=40usize { + let x = Array1::from_iter((0..n).map(|i| EDGES[i % EDGES.len()])); + let mut out = Array1::::zeros(n); + kernel(x.view(), out.view_mut()); + for i in 0..n { + let (got, want) = (out[i], scalar(x[i])); + if want.is_nan() { + assert!(got.is_nan(), "{name}({}) n={n} i={i}: got {got}", x[i]); + } else { + assert_eq!( + got.to_bits(), + want.to_bits(), + "{name}({}) n={n} i={i}: got {got:?}, want {want:?}", + x[i] + ); + } + } + } + } + + #[test] + fn vsneg_flips_the_sign_of_zero_in_the_simd_body() { + assert_bitwise_like_scalar("vsneg", vsneg, |v| -v); + } + + #[test] + fn vsceil_keeps_negative_zero_in_the_simd_body() { + assert_bitwise_like_scalar("vsceil", vsceil, f32::ceil); + } + + #[test] + fn vstrunc_keeps_negative_zero_in_the_simd_body() { + assert_bitwise_like_scalar("vstrunc", vstrunc, f32::trunc); + } + + /// FAILS IF: the SIMD body breaks ties differently from the scalar tail + /// (`round_ties_even`). On the AVX2, scalar and nightly arms + /// `F32x16::round` used `f32::round`, which rounds ties away from zero. + #[test] + fn vsround_breaks_ties_to_even_in_the_simd_body() { + assert_bitwise_like_scalar("vsround", vsround, f32::round_ties_even); + } + + #[test] + fn vsexp_maps_infinities_like_the_scalar_tail() { + for n in [1usize, 16, 17, 32] { + for (v, want) in [(f32::INFINITY, f32::INFINITY), (f32::NEG_INFINITY, 0.0)] { + let x = Array1::from_elem(n, v); + let mut out = Array1::::zeros(n); + vsexp(x.view(), out.view_mut()); + for i in 0..n { + assert_eq!(out[i].to_bits(), want.to_bits(), "vsexp({v}) n={n} i={i}: got {:?}", out[i]); + } + } + } + } + + /// FAILS IF: any finite input in [-110, 88.72] is off by more than the + /// polynomial's own error. Measured against f64 `exp`. Results below the + /// normal range may also be off by one subnormal step. + #[test] + fn vsexp_is_accurate_up_to_the_overflow_and_down_into_subnormals() { + let xs: Vec = (0..=19872) + .map(|i| -110.0 + i as f32 * 0.01) + .filter(|&x| x <= 88.72) + .collect(); + let x = Array1::from(xs); + let mut out = Array1::::zeros(x.len()); + vsexp(x.view(), out.view_mut()); + for (i, &v) in x.iter().enumerate() { + let want = f64::from(v).exp(); + let got = f64::from(out[i]); + if want >= f64::from(f32::MIN_POSITIVE) { + assert!(((got - want) / want).abs() < 4e-6, "vsexp({v}): got {got:e}, want {want:e}"); + } else { + // the polynomial's relative error plus one subnormal step + let step = f64::from(f32::from_bits(1)); + assert!((got - want).abs() <= 4e-6 * want + step, "vsexp({v}): got {got:e}, want {want:e}"); + } + } + } } diff --git a/src/simd.rs b/src/simd.rs index 3c618191..3e5c86d0 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -531,10 +531,10 @@ pub fn f32_to_bf16_batch(input: &[f32], output: &mut [u16]) { /// Max error ~2 ULP in [-10, 10]. Uses the standard range-reduction /// approach: exp(x) = 2^n * exp(r) where r = x - n*ln(2). /// -/// Domain: clamps input to [-87.336, 88.722] before reduction so that the -/// integer exponent `n` stays within the IEEE 754 f32 representable range. -/// Beyond the upper bound we'd hit `i32` overflow in `pow2n_from_int` and -/// silently return ~0.5 instead of +Inf (release) or panic (debug). +/// Domain: inputs from ln(f32::MAX) up return +Inf, and inputs are clamped +/// to [-110, ln(f32::MAX)] before reduction, so -Inf and anything below +/// about -103.9 return +0.0 and the range in between returns subnormals. +/// The clamp also keeps `n` far from `i32` overflow in `pow2n_from_int`. /// /// NaN handling: `simd_clamp` is `max(lo).min(hi)`, and `_mm512_max_ps` / /// `_mm512_min_ps` return the SECOND operand when the first is NaN (per @@ -548,7 +548,6 @@ pub fn f32_to_bf16_batch(input: &[f32], output: &mut [u16]) { #[inline(always)] #[allow(dead_code)] pub fn simd_exp_f32(x: F32x16) -> F32x16 { - let ln2 = F32x16::splat(core::f32::consts::LN_2); let inv_ln2 = F32x16::splat(1.0 / core::f32::consts::LN_2); let one = F32x16::splat(1.0); @@ -556,14 +555,30 @@ pub fn simd_exp_f32(x: F32x16) -> F32x16 { // Captured BEFORE the clamp because simd_clamp destroys NaN lanes. let nan_mask = x.simd_ne(x); - // Pre-clamp to the safe domain. Outside this band exp() is non-representable - // anyway (overflow → +Inf at ~88.7, underflow → +0 at ~-87.3) so the clamp - // is observable only at the saturation boundary. - let x = x.simd_clamp(F32x16::splat(-87.336_f32), F32x16::splat(88.722_f32)); - - // Range reduction: n = round(x / ln2), r = x - n * ln2 - let n = (x * inv_ln2).round(); - let r = x - n * ln2; + // At the upper clamp the true result already exceeds f32::MAX, so the + // final product overflows to +Inf for every input from there up. Below + // -110 the result rounds to +0.0 anyway, so clamping there only keeps + // `n` in range. + let x = x.simd_clamp(F32x16::splat(-110.0_f32), F32x16::splat(OVERFLOW_X)); + + // Range reduction: n = round(x / ln2), r = x - n * ln2, with ln2 split + // Cody-Waite style. LN2_HI has 9 significant bits, so n * LN2_HI is exact + // for |n| < 2^15 and x - n * LN2_HI is exact by Sterbenz's lemma. A + // 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. + // 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. + let n = (x * inv_ln2 + F32x16::splat(0.5)).floor(); + let r = x - n * F32x16::splat(LN2_HI); + let r = r - n * F32x16::splat(LN2_LO); // Polynomial: exp(r) ≈ 1 + r + r²/2 + r³/6 + r⁴/24 + r⁵/120 let c2 = F32x16::splat(0.5); @@ -573,13 +588,31 @@ pub fn simd_exp_f32(x: F32x16) -> F32x16 { let poly = one + r * (one + r * (c2 + r * (c3 + r * (c4 + r * c5)))); - // Reconstruct: exp(x) = 2^n * poly + // Reconstruct: exp(x) = 2^n * poly. 2^n is only a normal f32 for n in + // [-126, 127], but n reaches 128 near the overflow edge and about -159 + // at the lower clamp. Those lanes scale by 2^(n-1) then 2, or by + // 2^(n+64) then 2^-64, so the top results are not halved and the bottom + // ones land on the right subnormal instead of the 2^-126 floor. + let high = n.simd_gt(F32x16::splat(127.0)); + let low = n.simd_lt(F32x16::splat(-126.0)); + let n = high.select(n - one, low.select(n + F32x16::splat(64.0), n)); let result = poly * pow2n_from_int(n); + let result = high.select(result * F32x16::splat(2.0), result); + let result = low.select(result * F32x16::splat(TWO_POW_NEG_64), result); // Restore NaN in lanes where the input was NaN (clamp had destroyed them). nan_mask.select(F32x16::splat(f32::NAN), result) } +/// f32 just above ln(f32::MAX): `exp` of it overflows to +inf. +const OVERFLOW_X: f32 = 88.722_84; +/// High part of ln 2 (Cephes): 355/512, 9 significant bits. +const LN2_HI: f32 = 0.693_359_4; +/// ln 2 - LN2_HI. +const LN2_LO: f32 = -2.121_944_4e-4; +/// 2^-64, exact in f32. +const TWO_POW_NEG_64: f32 = 5.421_011e-20; + /// Compute 2^n where n is an integer stored as f32. /// /// Uses the IEEE 754 trick: set the exponent field directly. @@ -1649,29 +1682,24 @@ mod tests { #[test] fn simd_exp_f32_handles_positive_infinity() { - // Pre-fix: pow2n_from_int saturated f32::INFINITY to i32::MAX, - // (i32::MAX + 127) panicked in debug / wrapped in release to a - // garbage exponent, and simd_exp_f32(+Inf) silently returned ~0.5. - // Post-fix: input is clamped to 88.722 → exp(88.722) ≈ 3.4e38, - // representable but near f32::MAX. Saturated, not garbage. + // Once returned ~0.5 (i32 overflow in pow2n_from_int), then a + // finite 1.7e38 from the clamp. IEEE exp(+Inf) is +Inf. let inf = F32x16::splat(f32::INFINITY); let result = simd_exp_f32(inf); let arr = result.to_array(); for &v in &arr { - assert!(v.is_finite(), "exp(+Inf) must saturate to finite, got {}", v); - assert!(v > 1e30, "exp(+Inf) must saturate to a large value, got {}", v); + assert_eq!(v, f32::INFINITY, "exp(+Inf) must be +Inf, got {}", v); } } #[test] fn simd_exp_f32_handles_negative_infinity() { - // -Inf → clamped to -87.336 → exp ≈ 1.4e-38, near zero but representable. + // IEEE exp(-Inf) is +0.0 (it used to be ~1.2e-38 from the clamp). let neg_inf = F32x16::splat(f32::NEG_INFINITY); let result = simd_exp_f32(neg_inf); let arr = result.to_array(); for &v in &arr { - assert!(v.is_finite(), "exp(-Inf) must saturate to finite, got {}", v); - assert!(v >= 0.0 && v < 1e-30, "exp(-Inf) must saturate near 0, got {}", v); + assert_eq!(v.to_bits(), 0.0f32.to_bits(), "exp(-Inf) must be +0.0, got {}", v); } } @@ -1715,11 +1743,56 @@ mod tests { // is still in u32 range so didn't panic, but the resulting bits were // outside valid f32 exponent range, producing garbage that masqueraded // as a "valid" answer. + // IEEE: exp(200) overflows to +Inf. let big = F32x16::splat(200.0); let result = simd_exp_f32(big); let arr = result.to_array(); for &v in &arr { - assert!(v.is_finite(), "exp(200) must saturate, got {}", v); + assert_eq!(v, f32::INFINITY, "exp(200) must overflow to +Inf, got {}", v); + } + } + + /// The exact op sequence of `simd_exp_f32`, one lane, in plain scalar + /// f32. Scalar `*`, `+` and `floor` round identically on every target, + /// so a vector kernel that matches this bit for bit gives the same bits + /// on every backend, wasm included. + fn exp_reference_lane(x: f32) -> f32 { + if x.is_nan() { + return f32::NAN; + } + let xc = x.max(-110.0).min(super::OVERFLOW_X); + let n = (xc * (1.0 / core::f32::consts::LN_2) + 0.5).floor(); + let r = xc - n * super::LN2_HI; + let r = r - n * super::LN2_LO; + let poly = 1.0 + r * (1.0 + r * (0.5 + r * (1.0 / 6.0 + r * (1.0 / 24.0 + r * (1.0 / 120.0))))); + let (n, post) = if n > 127.0 { + (n - 1.0, 2.0) + } else if n < -126.0 { + (n + 64.0, super::TWO_POW_NEG_64) + } else { + (n, 1.0) + }; + let scale = f32::from_bits((((n as i32).clamp(-126, 127) + 127) as u32) << 23); + poly * scale * post + } + + /// FAILS IF: `simd_exp_f32` uses a fused or tie-dependent operation, so + /// its bits could differ between backends. + #[test] + fn simd_exp_f32_matches_the_unfused_scalar_sequence_bit_for_bit() { + let mut xs: Vec = (0..=20_000).map(|i| -110.0 + i as f32 * 0.01).collect(); + // ties of x/ln2 + 0.5, the overflow edge, specials + xs.extend((-150..=128).map(|k| (k as f32 + 0.5) * core::f32::consts::LN_2)); + xs.extend([f32::INFINITY, f32::NEG_INFINITY, super::OVERFLOW_X, 88.72, -103.9, 0.0, -0.0, f32::MIN_POSITIVE]); + while xs.len() % 16 != 0 { + xs.push(1.0); + } + for chunk in xs.chunks_exact(16) { + let got = simd_exp_f32(F32x16::from_slice(chunk)).to_array(); + for (&x, g) in chunk.iter().zip(got) { + let want = exp_reference_lane(x); + assert_eq!(g.to_bits(), want.to_bits(), "exp({x}): got {g:e}, want {want:e}"); + } } } diff --git a/src/simd_avx2.rs b/src/simd_avx2.rs index 2387f4a6..ea158dd0 100644 --- a/src/simd_avx2.rs +++ b/src/simd_avx2.rs @@ -583,7 +583,9 @@ impl F32x16 { let a = self.to_array(); let mut o = [0.0f32; 16]; for i in 0..16 { - o[i] = a[i].round(); + // Ties to even, like AVX-512 roundscale and NEON/wasm nearest; + // f32::round / f64::round round ties away from zero. + o[i] = a[i].round_ties_even(); } Self::from_array(o) } @@ -924,7 +926,9 @@ impl F64x8 { let a = self.to_array(); let mut o = [0.0f64; 8]; for i in 0..8 { - o[i] = a[i].round(); + // Ties to even, like AVX-512 roundscale and NEON/wasm nearest; + // f32::round / f64::round round ties away from zero. + o[i] = a[i].round_ties_even(); } Self::from_array(o) } diff --git a/src/simd_nightly/f32_types.rs b/src/simd_nightly/f32_types.rs index 3bbf1c7c..b78368cc 100644 --- a/src/simd_nightly/f32_types.rs +++ b/src/simd_nightly/f32_types.rs @@ -169,7 +169,8 @@ impl F32x16 { /// Per-lane round-to-nearest (ties to even). #[inline(always)] pub fn round(self) -> Self { - Self(self.0.round()) + // `StdFloat::round` rounds ties away from zero; per lane instead. + Self(self.0.to_array().map(f32::round_ties_even).into()) } /// Per-lane floor (toward negative infinity). @@ -359,7 +360,8 @@ impl F32x8 { /// Per-lane round-to-nearest (ties to even). #[inline(always)] pub fn round(self) -> Self { - Self(self.0.round()) + // `StdFloat::round` rounds ties away from zero; per lane instead. + Self(self.0.to_array().map(f32::round_ties_even).into()) } /// Per-lane floor (toward negative infinity). diff --git a/src/simd_nightly/f64_types.rs b/src/simd_nightly/f64_types.rs index 3d85881c..6d58a8e0 100644 --- a/src/simd_nightly/f64_types.rs +++ b/src/simd_nightly/f64_types.rs @@ -114,7 +114,8 @@ impl F64x8 { /// Per-lane round-to-nearest-even. #[inline(always)] pub fn round(self) -> Self { - Self(self.0.round()) + // `StdFloat::round` rounds ties away from zero; per lane instead. + Self(self.0.to_array().map(f64::round_ties_even).into()) } /// Per-lane floor. @@ -283,7 +284,8 @@ impl F64x4 { /// Per-lane round-to-nearest-even. #[inline(always)] pub fn round(self) -> Self { - Self(self.0.round()) + // `StdFloat::round` rounds ties away from zero; per lane instead. + Self(self.0.to_array().map(f64::round_ties_even).into()) } /// Per-lane floor. diff --git a/src/simd_scalar.rs b/src/simd_scalar.rs index 3308ee3e..a1fb5e62 100644 --- a/src/simd_scalar.rs +++ b/src/simd_scalar.rs @@ -132,7 +132,9 @@ macro_rules! impl_float_type { pub fn round(self) -> Self { let mut out = [0.0 as $elem; $lanes]; for i in 0..$lanes { - out[i] = self.0[i].round(); + // Ties to even, like the SIMD backends; `round` breaks + // ties away from zero. + out[i] = self.0[i].round_ties_even(); } Self(out) } diff --git a/src/simd_wasm.rs b/src/simd_wasm.rs index b0a2e0ca..0e2f5175 100644 --- a/src/simd_wasm.rs +++ b/src/simd_wasm.rs @@ -46,10 +46,10 @@ //! can differ in the last ULP — the same reduction-order tolerance the //! AVX/NEON backends already carry. //! * **round().** Uses `f32x4_nearest` / `f64x2_nearest` (round-half-to-even), -//! identical to NEON's `vrndnq_*`. The pure-scalar fallback uses -//! `f32::round` (round-half-away-from-zero); the two differ only on exact -//! `.5` ties, which the parity corpus avoids — the same situation that -//! already exists between the scalar and NEON backends. +//! identical to NEON's `vrndnq_*` and AVX-512 `roundscale`. The scalar, +//! AVX2 and nightly backends use `round_ties_even` for the same result +//! (they used `f32::round`, which breaks ties away from zero, until +//! `vsround_breaks_ties_to_even_in_the_simd_body` caught it). //! * **min/max NaN.** `f32x4_min` / `f32x4_max` follow IEEE minimum/maximum //! (NaN-propagating). The scalar fallback's `f32::min` returns the non-NaN //! operand. `simd_exp_f32` in `simd.rs` saves/restores NaN lanes around the