diff --git a/.claude/knowledge/synergy-map-2026-10-29.md b/.claude/knowledge/synergy-map-2026-10-29.md index 8f7b1647..ec55f199 100644 --- a/.claude/knowledge/synergy-map-2026-10-29.md +++ b/.claude/knowledge/synergy-map-2026-10-29.md @@ -323,6 +323,31 @@ and `lance_graph_onnx::` name in `registry.rs` against `simd.rs` / `cargo update`; `quote-harvester` for ORT's reduction kernels when V13 starts. Spot-check: re-read every hit used, re-run one count. +## R6. H11 and L18 implemented (2026-10-30) + +The two ndarray-side gaps R5 named for the ONNX registry probe: + +- **H11 → CLOSED.** `hpc::vml::vserf` (ONNX `Erf`): the Cephes `erff` + polynomial for `|x| < 1`, `1 − exp(−x²)/x · P(1/x²)` up to 4, and `1.0` + from 4. At most 2 ULP from the correctly rounded value over `[1e-38, 100]`, + both signs, chunk and tail (tested against an `f64` series). It is exactly + odd, keeps `erf(−0) = −0` and maps `±inf` to `±1`; NaN stays NaN. The Cephes + `erfcf` polynomial for `x ≥ 2` turned out unnecessary for `erf` (a disable + run that dropped it stayed green, and a prototype measured the same 2 ULP + with one polynomial), so it was removed. +- **L18 → CLOSED.** `simd::blend_i32_scalar(mask, a, b: i32, dst)`: ONNX + `Where(C, X, y)` with a scalar `y`, without the caller filling a `b` + vector. Pinned against `blend_i32` with a filled vector over 131 lanes + (two full words, a partial word, stray bits past the end). + +Disable runs, each restored from the commit: sign restore removed (3 tests +red), NaN clamped by a `< 4` select (special-values test red), small branch +on `x` instead of `|x|` (2 red), scalar ignored in `blend_i32_scalar` (red). + +What is left for the probe is lance-graph work: register `Erf` over `vserf`, +switch its `Where` DO to `blend_i32_scalar` when `Y` is a scalar (this also +closes half of L17), plus L11–L15. V13 stays a probe. + ## 0. What the six ndarray PRs did | PR | closed | opened or left | diff --git a/src/hpc/vml.rs b/src/hpc/vml.rs index d5499616..da2d158c 100644 --- a/src/hpc/vml.rs +++ b/src/hpc/vml.rs @@ -357,6 +357,77 @@ fn vstrunc_slice(x: &[f32], out: &mut [f32]) { } // =========================================================================== +/// 16-lane `erf`, computed on `a = |x|` and negated where `x < 0`, so the +/// function is exactly odd. Lanes where `x == 0` return `x` itself, which keeps +/// `erf(-0.0) = -0.0`. Float compares and `select` only, like [`tanh_f32x16`]. +/// +/// Three ranges, chosen per lane (the Cephes `erff`/`erfcf` minimax +/// polynomials, the `erfcf` one for `x < 2`): +/// - `a < 1`: `a · T(a²)`. No cancellation; for subnormal or tiny `a` the +/// result is `a · 1.1283792` (`2/√π · a`), the correct leading term. +/// - `1 ≤ a < 4`: `1 − exp(−a²)/a · P(1/a²)`. `erfc ≤ 0.16` here, so the +/// subtraction loses no bits. Cephes uses a second polynomial for `erfc` +/// from 2; `erf` does not need it, because `erfc < 0.005` there and the +/// error of `1 − erfc` stays within the same bound (measured: 2 ULP worst +/// case over `[1, 4)` with one polynomial or two). +/// - `a ≥ 4` (including `+inf`): exactly `1.0`; `erf(4) = 1 − 1.5e-8` rounds +/// to `1.0` in `f32`. +/// +/// A NaN input takes the middle branch and stays NaN. +#[inline(always)] +fn erf_f32x16(v: F32x16) -> F32x16 { + let zero = F32x16::splat(0.0); + let one = F32x16::splat(1.0); + let a = v.abs(); + let z = a * a; + let small = F32x16::splat(7.853_861_353_153_693e-5) + .mul_add(z, F32x16::splat(-8.010_193_625_184_903e-4)) + .mul_add(z, F32x16::splat(5.188_327_685_732_524e-3)) + .mul_add(z, F32x16::splat(-2.685_381_193_529_856e-2)) + .mul_add(z, F32x16::splat(1.128_358_514_861_418e-1)) + .mul_add(z, F32x16::splat(-3.761_262_582_423_300e-1)) + .mul_add(z, F32x16::splat(1.128_379_165_726_710)) + * a; + let q = one / a; + let y = q * q; + let p = F32x16::splat(2.326_819_970_068_386e-2) + .mul_add(y, F32x16::splat(-1.387_039_388_740_657e-1)) + .mul_add(y, F32x16::splat(3.687_424_674_597_105e-1)) + .mul_add(y, F32x16::splat(-5.824_733_027_278_666e-1)) + .mul_add(y, F32x16::splat(6.210_004_621_745_983e-1)) + .mul_add(y, F32x16::splat(-4.944_515_323_274_145e-1)) + .mul_add(y, F32x16::splat(3.404_879_937_665_872e-1)) + .mul_add(y, F32x16::splat(-2.741_127_028_184_656e-1)) + .mul_add(y, F32x16::splat(5.638_259_427_386_472e-1)); + let mid = one - simd_exp_f32(-z) * q * p; + // `a < 4` is false for NaN, so test `a >= 4` the other way round: NaN + // keeps `mid` (NaN) instead of being clamped to 1. + let big = a.simd_ge(F32x16::splat(4.0)); + let r = big.select(one, mid); + let r = a.simd_lt(one).select(small, r); + let r = v.simd_lt(zero).select(-r, r); + v.simd_eq(zero).select(v, r) +} + +fn vserf_slice(x: &[f32], out: &mut [f32]) { + let n = x.len().min(out.len()); + let chunks = n / 16; + for i in 0..chunks { + let off = i * 16; + erf_f32x16(F32x16::from_slice(&x[off..off + 16])).copy_to_slice(&mut out[off..off + 16]); + } + // The tail goes through the same 16-lane kernel, so a value gives the same + // bits wherever it sits in the array. + let rest = n - chunks * 16; + if rest > 0 { + let mut buf = [0.0f32; 16]; + buf[..rest].copy_from_slice(&x[chunks * 16..n]); + let r = erf_f32x16(F32x16::from_slice(&buf)); + r.copy_to_slice(&mut buf); + out[chunks * 16..n].copy_from_slice(&buf[..rest]); + } +} + // Internal helpers for the contiguous fast path. // =========================================================================== @@ -653,6 +724,36 @@ pub fn vstanh(x: ArrayView<'_, f32, D>, mut out: ArrayViewMut<'_, .for_each(|o, &v| *o = tanh_f32x16(F32x16::splat(v)).to_array()[0]); } +/// Element-wise error function `erf(x)` (ONNX `Erf`). +/// +/// At most 2 ULP from the correctly rounded value over `[-5, 5]` (measured +/// against an `f64` series in the tests). Exactly odd, `erf(±0) = ±0`, +/// `erf(±inf) = ±1`, NaN stays NaN. Every element, contiguous or strided, +/// goes through the same 16-lane kernel. +/// +/// # Example +/// ``` +/// use ndarray::{array, Array1}; +/// use ndarray::hpc::vml::vserf; +/// let x = array![0.0f32, 0.5, -1.0, 10.0]; +/// let mut out = Array1::::zeros(4); +/// vserf(x.view(), out.view_mut()); +/// assert!((out[1] - 0.520_499_9).abs() < 1e-6); +/// assert_eq!(out[3], 1.0); +/// ``` +/// +/// # Panics +/// Panics if `x` and `out` shapes differ. +pub fn vserf(x: ArrayView<'_, f32, D>, mut out: ArrayViewMut<'_, f32, D>) { + assert_eq!(x.shape(), out.shape(), "vserf: x/out shape mismatch"); + if dispatch_unary_contig(x.view(), &mut out, vserf_slice) { + return; + } + Zip::from(&mut out) + .and(&x) + .for_each(|o, &v| *o = erf_f32x16(F32x16::splat(v)).to_array()[0]); +} + /// Element-wise `floor`. Hot path uses `F32x16::floor()` (VRNDSCALEPS on /// AVX-512). /// @@ -1540,4 +1641,101 @@ mod tests { } } } + + /// `erf` in `f64` by its Maclaurin series, rounded once to `f32`. The + /// largest term at `|x| = 4.5` is about `1e8`, so the series keeps ~8 + /// digits more than `f32` needs. Beyond 4.5, `erf` rounds to `±1.0`. + fn erf_ref(x: f32) -> f32 { + let xf = x as f64; + if xf.is_nan() { + return f32::NAN; + } + if xf.abs() >= 4.5 { + return xf.signum() as f32; + } + let x2 = xf * xf; + let mut term = xf; // x^(2n+1) (-1)^n / n! + let mut sum = 0.0f64; + for n in 0..400u32 { + let t = term / (2 * n + 1) as f64; + sum += t; + if t.abs() < 1e-30 * sum.abs().max(1e-300) { + break; + } + term *= -x2 / (n + 1) as f64; + } + (sum * core::f64::consts::FRAC_2_SQRT_PI) as f32 + } + + fn vserf_vec(xs: &[f32]) -> Vec { + let x = Array1::from(xs.to_vec()); + let mut out = Array1::::zeros(xs.len()); + vserf(x.view(), out.view_mut()); + out.to_vec() + } + + #[test] + fn test_vserf_within_2_ulp_across_magnitudes() { + // The bands meet at the kernel's branch points 1, 2 and 4. The length + // is not a multiple of 16, so the last values take the tail path. + let bands = [(1.0e-38f32, 1.0e-10f32), (1.0e-10, 1.0e-3), (1.0e-3, 1.0), (1.0, 2.0), (2.0, 4.0), (4.0, 100.0)]; + for (lo, hi) in bands { + let n = 4099; + let mut xs = Vec::with_capacity(2 * n + 1); + for i in 0..n { + let t = i as f64 / (n - 1) as f64; + let v = ((lo as f64).ln() + t * ((hi as f64).ln() - (lo as f64).ln())).exp() as f32; + xs.push(v); + xs.push(-v); + } + xs.push(hi); + assert_ne!(xs.len() % 16, 0, "fixture must reach the tail path"); + let ys = vserf_vec(&xs); + let (worst, at) = xs + .iter() + .zip(&ys) + .map(|(&x, &y)| (ulp_distance(y, erf_ref(x)), x)) + .max_by_key(|&(u, _)| u) + .unwrap(); + assert!(worst <= 2, "band [{lo:e}, {hi:e}]: {worst} ulp at x = {at:e}"); + } + } + + #[test] + fn test_vserf_special_values_in_chunk_and_tail() { + let cases: [(f32, f32); 6] = + [(0.0, 0.0), (-0.0, -0.0), (f32::INFINITY, 1.0), (f32::NEG_INFINITY, -1.0), (4.0, 1.0), (-7.5, -1.0)]; + // 17 lanes: the first 16 fill one chunk, the 17th is the tail. + for (x, want) in cases { + let ys = vserf_vec(&[x; 17]); + for (i, &y) in ys.iter().enumerate() { + assert_eq!(y.to_bits(), want.to_bits(), "lane {i}: erf({x:e}) = {y:e}, want {want:e}"); + } + } + for (i, y) in vserf_vec(&[f32::NAN; 17]).iter().enumerate() { + assert!(y.is_nan(), "lane {i}: erf(NaN) = {y}"); + } + } + + #[test] + fn test_vserf_is_exactly_odd() { + let xs: Vec = (0..2000).map(|i| (i as f32) * 0.003_1 + 1.0e-7).collect(); + let neg: Vec = xs.iter().map(|&v| -v).collect(); + let (p, n) = (vserf_vec(&xs), vserf_vec(&neg)); + for ((&x, &a), &b) in xs.iter().zip(&p).zip(&n) { + assert_eq!((-a).to_bits(), b.to_bits(), "erf(-{x}) != -erf({x})"); + } + } + + #[test] + fn test_vserf_strided_matches_contiguous() { + let x = Array2::from_shape_fn((37, 5), |(i, j)| (i as f32 - 18.0) * 0.21 + j as f32 * 0.05); + let xt = x.t(); + let mut strided = Array2::::zeros((5, 37)); + vserf(xt, strided.view_mut()); + for ((i, j), &y) in strided.indexed_iter() { + let want = vserf_vec(&[xt[[i, j]]])[0]; + assert_eq!(y.to_bits(), want.to_bits(), "({i}, {j})"); + } + } } diff --git a/src/simd.rs b/src/simd.rs index caf56d91..348961f2 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -881,6 +881,7 @@ pub use crate::simd_masking_ops::{ // is NOT free — `U64x8` is 8 lanes, so eight groups share a word and each // lands at its own byte. blend_i32, + blend_i32_scalar, eq_f32_to_mask, eq_f64_to_mask, eq_i32_to_mask, diff --git a/src/simd_masking_ops.rs b/src/simd_masking_ops.rs index a880a06d..1cb80850 100644 --- a/src/simd_masking_ops.rs +++ b/src/simd_masking_ops.rs @@ -5840,6 +5840,43 @@ pub fn blend_i32(mask_words: &[u64], a: &[i32], b: &[i32], dst: &mut [i32]) { } } +/// `dst[i] = if mask bit i { a[i] } else { b }`: [`blend_i32`] with a scalar +/// `else` value, the ONNX `Where(C, X, y)` shape with a scalar `y`. It saves the +/// caller filling a `b` vector. Bits at or past `a.len()` are ignored. Plain +/// index loop, for the same reason as [`blend_i32`]. +/// +/// # Panics +/// +/// Panics unless `a.len() == dst.len()` and +/// `mask_words.len() >= a.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::blend_i32_scalar; +/// +/// let a = [1i32, 2, 3, 4]; +/// let mut dst = [0i32; 4]; +/// // bits 0 and 2 (0b0101) pick from `a`; bits 1 and 3 take the scalar +/// blend_i32_scalar(&[0b0101], &a, -1, &mut dst); +/// assert_eq!(dst, [1, -1, 3, -1]); +/// ``` +#[inline] +pub fn blend_i32_scalar(mask_words: &[u64], a: &[i32], b: i32, dst: &mut [i32]) { + assert_eq!(a.len(), dst.len(), "blend_i32_scalar: a/dst length mismatch"); + let n = a.len(); + assert!( + mask_words.len() >= mask_words_for(n), + "blend_i32_scalar: mask_words.len()={} < required {}", + mask_words.len(), + mask_words_for(n) + ); + for i in 0..n { + let bit = (mask_words[i / 64] >> (i % 64)) & 1; + dst[i] = if bit == 1 { a[i] } else { b }; + } +} + // ──────────────────────────────────────────────────────────────────────── // Morton lattice axis shift — the D-GTM-1m word-level op (§15); the hex // probe composes its six rails from these four Cartesian axis moves @@ -8777,6 +8814,31 @@ mod tests { assert_eq!(masked_max_i32(&v, &[u64::MAX]), Some(9)); } + #[test] + fn blend_i32_scalar_equals_blend_i32_with_a_filled_vector() { + // 131 lanes: two full words and a partial third, with stray bits set + // past the end that must be ignored. + let n = 131; + let a: Vec = (0..n as i32).map(|i| i * 7 - 300).collect(); + let mut m = [0u64; 3]; + for i in 0..n { + if (i * 37) % 5 < 2 { + m[i / 64] |= 1 << (i % 64); + } + } + m[2] |= !0u64 << (n % 64); + let (mut want, mut got) = (vec![0i32; n], vec![0i32; n]); + blend_i32(&m, &a, &vec![17; n], &mut want); + blend_i32_scalar(&m, &a, 17, &mut got); + assert_eq!(got, want); + let taken = got + .iter() + .zip(&a) + .filter(|(g, x)| g == x && **x != 17) + .count(); + assert!(taken > 0 && taken < n, "fixture must take both sides: {taken}"); + } + #[test] fn blend_i32_selects_by_bit_and_ignores_bits_past_len() { for &n in &LENS {