From 896d2876aeae94e689f24fcb857eb0827d574612 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Oct 2026 23:02:12 +0000 Subject: [PATCH 1/6] =?UTF-8?q?test(vml):=20IEEE=20edge=20cases=20through?= =?UTF-8?q?=20the=2016-lane=20body=20(neg/ceil/trunc=20=C2=B10,=20exp=20?= =?UTF-8?q?=C2=B1inf,=20exp=20range)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every existing vsneg/vsceil/vstrunc/vsexp test uses fewer than 16 elements, so it only ever exercises the scalar tail. These run 0..=40 lengths so each value lands in the SIMD body too. All five fail on master. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg --- src/hpc/vml.rs | 90 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 89 insertions(+), 1 deletion(-) diff --git a/src/hpc/vml.rs b/src/hpc/vml.rs index 1d9513b0..132506cf 100644 --- a/src/hpc/vml.rs +++ b/src/hpc/vml.rs @@ -688,7 +688,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 +1215,92 @@ 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); + } + + #[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 are checked against the nearest subnormal spacing. + #[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 { + assert!((got - want).abs() <= 1.5e-45, "vsexp({v}): got {got:e}, want {want:e}"); + } + } + } } From 2b1f6e720c389abdb0741ad98de526da369994e8 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Oct 2026 23:10:39 +0000 Subject: [PATCH 2/6] fix(vml, simd_exp_f32): sign of zero in neg/ceil/trunc; exp range and reduction - vsneg, vsceil, vstrunc: the 16-lane bodies computed 0.0 - v, which turns -0.0 into +0.0 (vsneg(+0), vsceil(-0.5), vsceil(-0.0), vstrunc(-0.5), vstrunc(-0.0)). They now use Neg, a sign-bit flip on every backend, and match their own scalar tails bit for bit. - simd_exp_f32 (vsexp, and every caller of the polynomial): * +inf returned 1.7e38; anything above ln(f32::MAX) now returns +inf. * pow2n_from_int clamps n to 127, so every input above ~88.38 came back halved. Those lanes now scale by 2^(n-1) and then 2. * The lower clamp at -87.336 made -inf and every smaller input return ~1.18e-38. The clamp is now -110 and lanes with n < -126 scale by 2^(n+64) then 2^-64 (the first product is exact), so -inf gives +0.0 and the subnormal range is correct. * Range reduction uses Cody-Waite ln2 hi/lo with mul_add: x - n*ln2 in one f32 product had an error larger than the polynomial's. The polynomial is Horner with mul_add. - vstanh: -two_x instead of 0.0 - two_x (no observable change; exp of either zero is 1). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg --- src/hpc/vml.rs | 32 ++++++++++++--------------- src/simd.rs | 59 +++++++++++++++++++++++++++++++++++++------------- 2 files changed, 58 insertions(+), 33 deletions(-) diff --git a/src/hpc/vml.rs b/src/hpc/vml.rs index 132506cf..924b8e9f 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,10 @@ 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(); @@ -1286,7 +1280,7 @@ mod tests { /// 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 are checked against the nearest subnormal spacing. + /// 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(); @@ -1299,7 +1293,9 @@ mod tests { if want >= f64::from(f32::MIN_POSITIVE) { assert!(((got - want) / want).abs() < 4e-6, "vsexp({v}): got {got:e}, want {want:e}"); } else { - assert!((got - want).abs() <= 1.5e-45, "vsexp({v}): got {got:e}, want {want:e}"); + // 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..2678bb3a 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 above ln(f32::MAX) 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,21 @@ 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 + // Lanes above ln(f32::MAX) overflow; they are set to +Inf at the end. + // Below -110 the result rounds to +0.0 anyway, so clamping there only + // keeps `n` in range. + let overflow = x.simd_gt(F32x16::splat(OVERFLOW_X)); + 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, even + // where `mul_add` is not fused. 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. let n = (x * inv_ln2).round(); - let r = x - n * ln2; + let r = (-n).mul_add(F32x16::splat(LN2_HI), x); + let r = (-n).mul_add(F32x16::splat(LN2_LO), r); // Polynomial: exp(r) ≈ 1 + r + r²/2 + r³/6 + r⁴/24 + r⁵/120 let c2 = F32x16::splat(0.5); @@ -571,15 +577,38 @@ pub fn simd_exp_f32(x: F32x16) -> F32x16 { let c4 = F32x16::splat(1.0 / 24.0); let c5 = F32x16::splat(1.0 / 120.0); - let poly = one + r * (one + r * (c2 + r * (c3 + r * (c4 + r * c5)))); - - // Reconstruct: exp(x) = 2^n * poly + let poly = r.mul_add(c5, c4); + let poly = r.mul_add(poly, c3); + let poly = r.mul_add(poly, c2); + let poly = r.mul_add(poly, one); + let poly = r.mul_add(poly, one); + + // 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); + let result = overflow.select(F32x16::splat(f32::INFINITY), result); // Restore NaN in lanes where the input was NaN (clamp had destroyed them). nan_mask.select(F32x16::splat(f32::NAN), result) } +/// ln(f32::MAX): above it `exp` 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. From 43ce2c6f7c1bd97db995d0ff2e59836c31281493 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Oct 2026 23:33:25 +0000 Subject: [PATCH 3/6] fix(simd_exp_f32): drop the unreachable overflow mask; GELU survives exp = +inf - The overflow select never fired: clamping at OVERFLOW_X already makes poly*2^127*2 overflow to +inf (disable run stayed green), so it is gone. - simd::tests pinned the old saturating contract (exp(+inf) finite, exp(200) finite, exp(-inf) ~1.4e-38). Re-pinned to IEEE: +inf, +inf, +0.0. - The two GELU kernels (hpc::models::layers::gelu, gpt2 gelu_simd) computed tanh as (e-1)/(e+1), which is inf/inf = NaN once exp overflows (x >= ~10). They now use the equal form 1 - 2/(e+1). New tests fail with NaN on the quotient form. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg --- src/hpc/gpt2/inference.rs | 18 +++++++++++++++++- src/hpc/models/layers.rs | 18 +++++++++++++++++- src/hpc/vml.rs | 20 ++++++++++++++++---- src/simd.rs | 31 +++++++++++++------------------ 4 files changed, 63 insertions(+), 24 deletions(-) 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 924b8e9f..c6df2d05 100644 --- a/src/hpc/vml.rs +++ b/src/hpc/vml.rs @@ -318,7 +318,9 @@ fn vstrunc_slice(x: &[f32], out: &mut [f32]) { // 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]); + 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(); @@ -1233,7 +1235,9 @@ mod tests { /// 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) { + 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); @@ -1243,7 +1247,12 @@ mod tests { 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]); + assert_eq!( + got.to_bits(), + want.to_bits(), + "{name}({}) n={n} i={i}: got {got:?}, want {want:?}", + x[i] + ); } } } @@ -1283,7 +1292,10 @@ mod tests { /// 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 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()); diff --git a/src/simd.rs b/src/simd.rs index 2678bb3a..cd5c3a3a 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -531,7 +531,7 @@ 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: inputs above ln(f32::MAX) return +Inf, and inputs are clamped +/// 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`. @@ -555,10 +555,10 @@ 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); - // Lanes above ln(f32::MAX) overflow; they are set to +Inf at the end. - // Below -110 the result rounds to +0.0 anyway, so clamping there only - // keeps `n` in range. - let overflow = x.simd_gt(F32x16::splat(OVERFLOW_X)); + // 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 @@ -594,13 +594,12 @@ pub fn simd_exp_f32(x: F32x16) -> F32x16 { 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); - let result = overflow.select(F32x16::splat(f32::INFINITY), result); // Restore NaN in lanes where the input was NaN (clamp had destroyed them). nan_mask.select(F32x16::splat(f32::NAN), result) } -/// ln(f32::MAX): above it `exp` overflows to +inf. +/// 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; @@ -1678,29 +1677,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); } } @@ -1744,11 +1738,12 @@ 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); } } From c0b1882303bb3541f0c2284dfc784fdf2b00b40f Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Oct 2026 23:54:06 +0000 Subject: [PATCH 4/6] fix(simd_exp_f32): unfused ops and a tie-free n, so exp is bit-identical on every backend mul_add is fused on AVX-512/AVX2+FMA/NEON/scalar but unfused on wasm (and implementation-defined with relaxed-simd), and round() breaks ties to even on AVX-512/NEON/wasm but away from zero on AVX2/scalar/nightly. Either made exp differ by an ulp between backends. The Cody-Waite split already makes the fused form unnecessary for accuracy, so the reduction and Horner use plain * and +, and n = floor(x/ln2 + 0.5). Same rule as gemm_f64_tiled. New test: every lane equals the same op sequence in scalar f32, bit for bit. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg --- src/simd.rs | 73 ++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 61 insertions(+), 12 deletions(-) diff --git a/src/simd.rs b/src/simd.rs index cd5c3a3a..3e5c86d0 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -563,13 +563,22 @@ pub fn simd_exp_f32(x: F32x16) -> F32x16 { // 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, even - // where `mul_add` is not fused. 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. - let n = (x * inv_ln2).round(); - let r = (-n).mul_add(F32x16::splat(LN2_HI), x); - let r = (-n).mul_add(F32x16::splat(LN2_LO), r); + // 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); @@ -577,11 +586,7 @@ pub fn simd_exp_f32(x: F32x16) -> F32x16 { let c4 = F32x16::splat(1.0 / 24.0); let c5 = F32x16::splat(1.0 / 120.0); - let poly = r.mul_add(c5, c4); - let poly = r.mul_add(poly, c3); - let poly = r.mul_add(poly, c2); - let poly = r.mul_add(poly, one); - let poly = r.mul_add(poly, one); + let poly = one + r * (one + r * (c2 + r * (c3 + r * (c4 + r * c5)))); // 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 @@ -1747,6 +1752,50 @@ mod tests { } } + /// 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}"); + } + } + } + // ── W1a-#9 parity: U64x8/U32x16 andnot + ternlog ──────────────────────── // // The backends are compile-time exclusive, so "all three agree" is proven From 4519689f97bdf2e50307bcaaf0bffea0f52187f7 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 24 Oct 2026 00:01:08 +0000 Subject: [PATCH 5/6] test(vml): vsround's SIMD body must break ties like its tail Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg --- src/hpc/vml.rs | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/hpc/vml.rs b/src/hpc/vml.rs index c6df2d05..00c77e7a 100644 --- a/src/hpc/vml.rs +++ b/src/hpc/vml.rs @@ -1273,6 +1273,14 @@ mod tests { 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] { From 5f4636e760b21c436e01ca79444861afc17a7b5f Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 24 Oct 2026 00:17:32 +0000 Subject: [PATCH 6/6] fix(simd): round() breaks ties to even on every backend; wasm parity checks exp and round F32x16/F64x8::round used f32::round (ties away from zero) on the AVX2, scalar and nightly backends but ties-to-even on AVX-512, NEON and wasm, so vsround's SIMD body disagreed with its own round_ties_even tail on AVX2 (vsround(0.5) = 1.0). All backends now use ties-to-even. wasm-simd-parity selfcheck: simd_exp_f32 equals the unfused scalar sequence bit for bit under node (code 70), and round breaks ties to even (code 71). With the native unit test on the same sequence, exp is bit-identical between wasm and native. The nightly change is untested here (no nightly toolchain); CI checks it. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg --- crates/wasm-simd-parity/src/lib.rs | 76 ++++++++++++++++++++++++++++++ src/simd_avx2.rs | 8 +++- src/simd_nightly/f32_types.rs | 6 ++- src/simd_nightly/f64_types.rs | 6 ++- src/simd_scalar.rs | 4 +- src/simd_wasm.rs | 8 ++-- 6 files changed, 97 insertions(+), 11 deletions(-) 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/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