From 02f901808d7484df56d5667c393f6847e04af7cf Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 11:39:53 +0000 Subject: [PATCH 1/4] Vectorize pow2n_from_int; add U32x16 variable shifts to every backend simd_exp_f32 (and with it vsexp, vstanh, vserf) built 2^n one lane at a time through to_array. It now stays in registers through the facade's typed words: clamp n, add 2^23 + 127 so the biased exponent sits in the mantissa field, and shift the bits left by 23. That needs U32x16 << U32x16, which only the AVX-512 and scalar backends had. Added Shl/Shr for U32x16 on AVX2 (VPSLLVD/VPSRLVD per 256-bit half), AVX without AVX2 (per lane; no AVX1 form exists), NEON (vshlq_u32), wasm simd128 (lane extracts; simd128 shifts are uniform) and the nightly core::simd backend, following each backend's U64x8 shifts. Tests: pow2n_from_int bit for bit against the scalar exponent rule over n in [-200, 200]; U32x16 shifts against per-lane shifts with distinct counts per lane. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- src/simd.rs | 51 ++++++++++++++++++++---- src/simd_avx.rs | 24 +++++++++++- src/simd_avx2.rs | 66 ++++++++++++++++++++++++++++++++ src/simd_neon.rs | 24 ++++++++++++ src/simd_nightly/u_word_types.rs | 19 +++++++++ src/simd_wasm.rs | 21 ++++++++++ 6 files changed, 196 insertions(+), 9 deletions(-) diff --git a/src/simd.rs b/src/simd.rs index 072dc41f..5b02fbab 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -639,14 +639,14 @@ const TWO_POW_NEG_64: f32 = 5.421_011e-20; #[inline(always)] #[allow(dead_code)] fn pow2n_from_int(n: F32x16) -> F32x16 { - let arr = n.to_array(); - let mut out = [0.0f32; 16]; - for i in 0..16 { - let ni = (arr[i] as i32).clamp(-126, 127); - let bits = ((ni + 127) as u32) << 23; - out[i] = f32::from_bits(bits); - } - F32x16::from_array(out) + // `n` is integral (it comes from `round`), so after the clamp + // `2^23 + 127 + n` is an exact f32 whose mantissa field holds the biased + // exponent `k = n + 127` in [1, 254] and whose exponent field is + // 0x4B00_0000. Shifting the bits left by 23 pushes 0x4B00_0000 out of the + // word and moves `k` into the exponent field: the bits of 2^n. No lane + // leaves its register. + let biased = n.simd_clamp(F32x16::splat(-126.0), F32x16::splat(127.0)) + F32x16::splat(8_388_608.0 + 127.0); + F32x16::from_bits(biased.to_bits() << U32x16::splat(23)) } /// Fast natural log for F32x16. @@ -1800,6 +1800,41 @@ mod tests { assert_eq!(&out[..16], &data[..16]); } + #[test] + fn u32x16_variable_shifts_match_per_lane_shifts() { + // Distinct value and count per lane, every count 0..32 used, so a + // backend that shifted by lane 0's count or a uniform count fails. + for base in 0..2u32 { + let a: [u32; 16] = core::array::from_fn(|i| 0x9E37_79B9u32.rotate_left(i as u32 * 5 + base)); + let n: [u32; 16] = core::array::from_fn(|i| (i as u32 * 7 + base * 16) % 32); + let (va, vn) = (U32x16::from_array(a), U32x16::from_array(n)); + let l = (va << vn).to_array(); + let r = (va >> vn).to_array(); + for i in 0..16 { + assert_eq!(l[i], a[i] << n[i], "shl lane {i}"); + assert_eq!(r[i], a[i] >> n[i], "shr lane {i}"); + } + } + } + + #[test] + fn pow2n_from_int_is_exactly_two_to_the_clamped_n() { + // Every integral n in [-200, 200] (so both clamp ends are crossed), + // checked bit for bit against the scalar exponent-field rule. + let ns: Vec = (-200..=200).map(|n| n as f32).collect(); + for chunk in ns.chunks(16) { + let mut lanes = [0.0f32; 16]; + lanes[..chunk.len()].copy_from_slice(chunk); + let got = pow2n_from_int(F32x16::from_array(lanes)).to_array(); + for (&n, &g) in lanes.iter().zip(&got) { + let k = ((n as i32).clamp(-126, 127) + 127) as u32; + assert_eq!(g.to_bits(), k << 23, "n={n}"); + } + } + assert_eq!(pow2n_from_int(F32x16::splat(0.0)).to_array(), [1.0; 16]); + assert_eq!(pow2n_from_int(F32x16::splat(-126.0)).to_array(), [f32::MIN_POSITIVE; 16]); + } + #[test] fn simd_exp_f32_basic() { let zero = F32x16::splat(0.0); diff --git a/src/simd_avx.rs b/src/simd_avx.rs index c4585f99..71c06291 100644 --- a/src/simd_avx.rs +++ b/src/simd_avx.rs @@ -28,7 +28,7 @@ use core::ops::{ Add, AddAssign, BitAnd, BitAndAssign, BitOr, BitOrAssign, BitXor, BitXorAssign, Mul, Not, Shl, Shr, Sub, SubAssign, }; -use crate::simd_avx2::{F32x16, I32x16, I8x64, U16x16, U16x32, U64x8, U8x32, U8x64}; +use crate::simd_avx2::{F32x16, I32x16, I8x64, U16x16, U16x32, U32x16, U64x8, U8x32, U8x64}; use crate::simd_avx512::{F32x8, I16x16, I8x32}; // --------------------------------------------------------------------------- @@ -308,6 +308,28 @@ impl Shr for U64x8 { } } +/// Lane-wise variable left shift. Counts of 32 or more zero the lane, as +/// `VPSLLVD` does on the AVX2 arm. No SSE/AVX1 form exists, so this is per +/// lane, like `U64x8`'s above. +impl Shl for U32x16 { + type Output = Self; + #[inline(always)] + fn shl(self, rhs: Self) -> Self { + debug_assert!(rhs.0.iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + Self(core::array::from_fn(|i| if rhs.0[i] < 32 { self.0[i] << rhs.0[i] } else { 0 })) + } +} + +/// Lane-wise variable right shift; see the `Shl` impl above. +impl Shr for U32x16 { + type Output = Self; + #[inline(always)] + fn shr(self, rhs: Self) -> Self { + debug_assert!(rhs.0.iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + Self(core::array::from_fn(|i| if rhs.0[i] < 32 { self.0[i] >> rhs.0[i] } else { 0 })) + } +} + // --------------------------------------------------------------------------- // U8x64 (array-backed) // --------------------------------------------------------------------------- diff --git a/src/simd_avx2.rs b/src/simd_avx2.rs index 52724f33..ecbeec2f 100644 --- a/src/simd_avx2.rs +++ b/src/simd_avx2.rs @@ -1921,6 +1921,72 @@ impl Shr for U64x8 { } } +#[cfg(not(all(target_feature = "avx", not(target_feature = "avx2"))))] +impl U32x16 { + /// The two 256-bit halves of the 64-byte array; see `U64x8::avx2_halves` + /// for the run-time AVX2 obligation and why both loads are in bounds. + #[inline(always)] + fn avx2_halves(self) -> (__m256i, __m256i) { + // SAFETY: as for `U64x8::avx2_halves`: 64-byte array, two 32-byte + // unaligned loads, AVX2 present by this arm's contract. + unsafe { + let p = self.0.as_ptr() as *const __m256i; + (_mm256_loadu_si256(p), _mm256_loadu_si256(p.add(1))) + } + } + + /// Store two 256-bit halves into a fresh `[u32; 16]`. + #[inline(always)] + fn from_avx2_halves(lo: __m256i, hi: __m256i) -> Self { + let mut o = [0u32; 16]; + // SAFETY: see `avx2_halves`; two 32-byte stores into a 64-byte array. + unsafe { + let p = o.as_mut_ptr() as *mut __m256i; + _mm256_storeu_si256(p, lo); + _mm256_storeu_si256(p.add(1), hi); + } + Self(o) + } +} + +#[cfg(not(all(target_feature = "avx", not(target_feature = "avx2"))))] +/// Lane-wise variable left shift (`VPSLLVD` per 256-bit half). Counts stay +/// below 32 by the callers' contract, as on every other arm. +/// +/// ``` +/// use ndarray::simd::U32x16; +/// let x = U32x16::splat(1); +/// let n = U32x16::from_array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 23, 31]); +/// let left = x << n; +/// assert_eq!(left.to_array()[14], 1 << 23); +/// assert_eq!((left >> n).to_array(), [1; 16]); +/// ``` +impl Shl for U32x16 { + type Output = Self; + #[inline(always)] + fn shl(self, rhs: Self) -> Self { + debug_assert!(rhs.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + let (lo, hi) = self.avx2_halves(); + let (clo, chi) = rhs.avx2_halves(); + // SAFETY: AVX2 present by the arm's contract; register-only ops. + unsafe { Self::from_avx2_halves(_mm256_sllv_epi32(lo, clo), _mm256_sllv_epi32(hi, chi)) } + } +} + +#[cfg(not(all(target_feature = "avx", not(target_feature = "avx2"))))] +/// Lane-wise variable right shift (`VPSRLVD`); see `Shl` above. +impl Shr for U32x16 { + type Output = Self; + #[inline(always)] + fn shr(self, rhs: Self) -> Self { + debug_assert!(rhs.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + let (lo, hi) = self.avx2_halves(); + let (clo, chi) = rhs.avx2_halves(); + // SAFETY: as for `shl`. + unsafe { Self::from_avx2_halves(_mm256_srlv_epi32(lo, clo), _mm256_srlv_epi32(hi, chi)) } + } +} + /// BLAKE3 `hash_many` shuffle surface — the transpose network's four unpacks /// and two half-concatenations, at DEGREE 16. /// diff --git a/src/simd_neon.rs b/src/simd_neon.rs index f89b58fa..ae9717e0 100644 --- a/src/simd_neon.rs +++ b/src/simd_neon.rs @@ -3245,6 +3245,30 @@ impl core::ops::Not for U64x8 { })) } } +/// Lane-wise `self << rhs` for `U32x16` (per-lane counts below 32). +#[cfg(target_arch = "aarch64")] +impl core::ops::Shl for U32x16 { + type Output = Self; + #[inline(always)] + fn shl(self, r: Self) -> Self { + debug_assert!(r.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + // SAFETY: NEON baseline; pure register ops. + Self(core::array::from_fn(|p| U32x4(unsafe { vshlq_u32(self.0[p].0, vreinterpretq_s32_u32(r.0[p].0)) }))) + } +} +/// Lane-wise `self >> rhs` for `U32x16` (`vshlq_u32` with negated counts). +#[cfg(target_arch = "aarch64")] +impl core::ops::Shr for U32x16 { + type Output = Self; + #[inline(always)] + fn shr(self, r: Self) -> Self { + debug_assert!(r.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + // SAFETY: NEON baseline; pure register ops. + Self(core::array::from_fn(|p| { + U32x4(unsafe { vshlq_u32(self.0[p].0, vnegq_s32(vreinterpretq_s32_u32(r.0[p].0))) }) + })) + } +} /// Lane-wise `self << rhs` (per-lane counts; a count of 64 or more yields 0, /// which is what the shift instruction does — the scalar backend's `<<` /// panics there in debug builds, so callers already stay inside `0..64`). diff --git a/src/simd_nightly/u_word_types.rs b/src/simd_nightly/u_word_types.rs index f03d26e7..33992a3b 100644 --- a/src/simd_nightly/u_word_types.rs +++ b/src/simd_nightly/u_word_types.rs @@ -562,6 +562,25 @@ impl Default for U32x8 { #[repr(transparent)] pub struct U32x16(pub u32x16); +/// Lane-wise variable shifts; counts stay below 32 by the callers' contract. +impl core::ops::Shl for U32x16 { + type Output = Self; + #[inline(always)] + fn shl(self, rhs: Self) -> Self { + debug_assert!(rhs.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + Self(self.0 << rhs.0) + } +} + +impl core::ops::Shr for U32x16 { + type Output = Self; + #[inline(always)] + fn shr(self, rhs: Self) -> Self { + debug_assert!(rhs.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + Self(self.0 >> rhs.0) + } +} + /// BLAKE3 `hash_many` shuffle surface — the transpose network's four unpacks /// and two half-concatenations, at DEGREE 16. /// diff --git a/src/simd_wasm.rs b/src/simd_wasm.rs index cd345a19..cb5b7635 100644 --- a/src/simd_wasm.rs +++ b/src/simd_wasm.rs @@ -2392,6 +2392,27 @@ pub mod wasm32_simd { Self(core::array::from_fn(|p| U64x2(v128_not(self.0[p].0)))) } } + /// Lane-wise `U32x16 << rhs` with per-lane counts below 32. simd128 has + /// only a uniform shift, so this goes through lane extracts like `U64x8`. + impl core::ops::Shl for U32x16 { + type Output = Self; + #[inline(always)] + fn shl(self, r: Self) -> Self { + debug_assert!(r.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + let (a, n) = (self.to_array(), r.to_array()); + Self::from_array(core::array::from_fn(|i| a[i].wrapping_shl(n[i]))) + } + } + /// Lane-wise `U32x16 >> rhs` with per-lane counts (see `Shl`). + impl core::ops::Shr for U32x16 { + type Output = Self; + #[inline(always)] + fn shr(self, r: Self) -> Self { + debug_assert!(r.to_array().iter().all(|&n| n < 32), "U32x16 shift counts are a caller contract: < 32"); + let (a, n) = (self.to_array(), r.to_array()); + Self::from_array(core::array::from_fn(|i| a[i].wrapping_shr(n[i]))) + } + } /// Lane-wise `self << rhs` with PER-LANE counts. simd128 has only a uniform /// shift, so this one operator goes through lane extracts; the scalar /// backend's `<<` panics on a count of 64 or more in debug builds, so From 9aea2fc1769e46523a995168e2b262c2205b6d54 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 11:47:21 +0000 Subject: [PATCH 2/4] Add half_add_mask and full_add_mask over 64 one-bit lanes; wasm parity for U32x16 shifts half_add_mask(a, carry) = (a ^ carry, a & carry) and full_add_mask(a, b, c) = (a ^ b ^ c, majority) are the flip and carry_save3 laws proven in #391 and #389, now in simd_masking_ops and re-exported through ndarray::simd. Each u64 is one significance column of 64 independent counters; the carry belongs one column up in the same lane. Tests: both truth tables in every lane position, lanes 0/31/32/63 alone, 2000 random words (with zero, all-ones and alternating masks) against per-lane addition, and a ripple counter that keeps lanes apart while moving the carry to a neighbouring bit corrupts it. The wasm parity harness now also checks U32x16 << and >> with a distinct count per lane. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- crates/wasm-simd-parity/src/lib.rs | 12 +++ src/simd.rs | 3 + src/simd_masking_ops.rs | 151 +++++++++++++++++++++++++++++ 3 files changed, 166 insertions(+) diff --git a/crates/wasm-simd-parity/src/lib.rs b/crates/wasm-simd-parity/src/lib.rs index 1efa6c79..a957b0ae 100644 --- a/crates/wasm-simd-parity/src/lib.rs +++ b/crates/wasm-simd-parity/src/lib.rs @@ -161,6 +161,18 @@ fn check_u32x16() -> Result<(), u32> { } } } + // Per-lane variable shifts: a distinct count in every lane. + let n_arr: [u32; 16] = core::array::from_fn(|i| (i as u32 * 7) % 32); + let n = U32x16::from_array(n_arr); + let (l, r) = ((a << n).to_array(), (a >> n).to_array()); + for i in 0..16 { + if l[i] != a_arr[i] << n_arr[i] { + return Err(14); + } + if r[i] != a_arr[i] >> n_arr[i] { + return Err(15); + } + } Ok(()) } diff --git a/src/simd.rs b/src/simd.rs index 5b02fbab..2b297fba 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -1042,6 +1042,9 @@ pub use crate::simd_masking_ops::{ BOUNDED_TILE_ROWS, SYM_EMPTY_I64, }; +// Binary positional adders over 64 one-bit lanes (the #389/#391 `carry_save3` +// and `flip` laws). +pub use crate::simd_masking_ops::{full_add_mask, half_add_mask}; #[cfg(feature = "std")] pub use crate::simd_stride::{lut9, lut_u8, FACET_BYTES}; // The popcount that closes the loop on the masks above: `mask_count` in ABI diff --git a/src/simd_masking_ops.rs b/src/simd_masking_ops.rs index ebfa444c..453cadf5 100644 --- a/src/simd_masking_ops.rs +++ b/src/simd_masking_ops.rs @@ -6494,6 +6494,59 @@ pub fn mask_shift_morton(src: &[u64], dir: MortonDir, dst: &mut [u64]) { } } +// ════════════════════════════════════════════════════════════════════════ +// Binary positional adders over 64 independent one-bit lanes +// ════════════════════════════════════════════════════════════════════════ +// +// A `u64` here is one significance column of 64 counters: bit `j` belongs to +// lane `j`, and every lane is added independently. These are the half adder +// (`flip`, examples/u9_domino) and full adder (`carry_save3`, +// examples/u9_abacus) laws proven in #389 and #391. The carry plane belongs +// one column up (significance x2) in the SAME lanes: callers move it to the +// next column, never shift it to a neighbouring bit. + +/// Half adder on 64 independent lanes: `(a ^ carry, a & carry)`. +/// +/// For every lane `j`, `a_j + carry_j = sum_j + 2·carry_out_j` exactly. The +/// carry out belongs to the next significance column, in lane `j` again. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::half_add_mask; +/// +/// let (sum, carry) = half_add_mask(0b0110, 0b0011); +/// assert_eq!((sum, carry), (0b0101, 0b0010)); +/// ``` +#[inline(always)] +#[must_use] +pub const fn half_add_mask(a: u64, carry: u64) -> (u64, u64) { + (a ^ carry, a & carry) +} + +/// Full adder (3:2 carry-save compressor) on 64 independent lanes. +/// +/// Returns `(sum, carry)` with `sum = a ^ b ^ c` and `carry` the per-lane +/// majority of `a`, `b`, `c`, so for every lane `j` +/// `a_j + b_j + c_j = sum_j + 2·carry_j` exactly. The carry belongs to the +/// next significance column, in lane `j` again. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::full_add_mask; +/// +/// // Lane 0: 1+1+1 = 1 + 2·1. Lane 1: 0+1+0 = 1. Lane 2: 1+1+0 = 0 + 2·1. +/// let (sum, carry) = full_add_mask(0b101, 0b111, 0b001); +/// assert_eq!((sum, carry), (0b011, 0b101)); +/// ``` +#[inline(always)] +#[must_use] +pub const fn full_add_mask(a: u64, b: u64, c: u64) -> (u64, u64) { + let t = a ^ b; + (t ^ c, (a & b) ^ (t & c)) +} + #[cfg(test)] mod tests { use super::*; @@ -11477,3 +11530,101 @@ mod bounded_register_tests { assert_eq!(out[0], PowerSums::default()); } } + +#[cfg(test)] +mod binary_adder_tests { + use super::{full_add_mask, half_add_mask}; + + fn bit(x: u64, j: usize) -> u64 { + (x >> j) & 1 + } + + /// Every row of both truth tables, in every lane position: each input bit + /// pattern is spread so lane j sees row (j % rows), so a lane mix-up or a + /// wrong row fails. + #[test] + fn truth_tables_hold_in_every_lane() { + let spread = |row_bit: usize, rows: usize| -> u64 { + (0..64).fold(0u64, |w, j| w | ((((j % rows) >> row_bit) & 1) as u64) << j) + }; + let (a, c) = (spread(0, 4), spread(1, 4)); + let (s, k) = half_add_mask(a, c); + for j in 0..64 { + let (x, y) = (bit(a, j), bit(c, j)); + assert_eq!((bit(s, j), bit(k, j)), ((x + y) & 1, (x + y) >> 1), "half lane {j}"); + } + let (a, b, c) = (spread(0, 8), spread(1, 8), spread(2, 8)); + let (s, k) = full_add_mask(a, b, c); + for j in 0..64 { + let t = bit(a, j) + bit(b, j) + bit(c, j); + assert_eq!((bit(s, j), bit(k, j)), (t & 1, t >> 1), "full lane {j}"); + } + } + + /// Lanes 0, 31, 32 and 63 alone: the result stays in the input's lane. + #[test] + fn a_single_lane_stays_in_its_lane() { + for j in [0, 31, 32, 63] { + let m = 1u64 << j; + assert_eq!(half_add_mask(m, m), (0, m), "half lane {j}"); + assert_eq!(half_add_mask(m, 0), (m, 0), "half lane {j}"); + assert_eq!(full_add_mask(m, m, m), (m, m), "full lane {j}"); + assert_eq!(full_add_mask(m, 0, m), (0, m), "full lane {j}"); + assert_eq!(full_add_mask(0, m, 0), (m, 0), "full lane {j}"); + } + } + + /// Random words against per-lane integer addition; the masks include zero, + /// all ones and alternating patterns. + #[test] + fn random_words_match_per_lane_addition() { + let mut s = 0x9E37_79B9_7F4A_7C15u64; + let mut next = || { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + s + }; + let fixed = [0, u64::MAX, 0xAAAA_AAAA_AAAA_AAAA, 0x5555_5555_5555_5555]; + for i in 0..2000 { + let pick = |i: usize, r: u64| if i % 5 < 4 { r } else { fixed[(r % 4) as usize] }; + let (a, b, c) = (pick(i, next()), pick(i + 1, next()), pick(i + 2, next())); + let (s2, k2) = half_add_mask(a, b); + let (s3, k3) = full_add_mask(a, b, c); + for j in 0..64 { + assert_eq!(bit(a, j) + bit(b, j), bit(s2, j) + 2 * bit(k2, j)); + assert_eq!(bit(a, j) + bit(b, j) + bit(c, j), bit(s3, j) + 2 * bit(k3, j)); + } + } + } + + /// A ripple counter built from the words keeps lanes apart: counting + /// lane j up `j + 1` times yields `j + 1` in lane j and nothing in any + /// other lane, across carry chains of up to 7 columns. Moving the carry + /// to a neighbouring bit instead (`carry << 1`) corrupts it (anti-test). + #[test] + fn rippling_carries_preserve_lane_identity() { + fn count(cols: &mut [u64; 8], m: u64, shift_carry: bool) { + let mut c = m; + for col in cols.iter_mut() { + let (s, k) = half_add_mask(*col, c); + *col = s; + c = if shift_carry { k << 1 } else { k }; + if c == 0 { + break; + } + } + } + let value = |cols: &[u64; 8], j: usize| (0..8).map(|k| bit(cols[k], j) << k).sum::(); + for shift in [false, true] { + let mut cols = [0u64; 8]; + for j in 0..64 { + for _ in 0..=j { + count(&mut cols, 1 << j, shift); + } + } + let ok = (0..64).all(|j| value(&cols, j) == j as u64 + 1); + assert_eq!(ok, !shift, "shift_carry={shift}"); + } + } +} From a70f336ec9a90b4ca43f85c519b97af6d4201026 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 12:06:39 +0000 Subject: [PATCH 3/4] u9_domino: candidate C, Abacus with transactional overflow refusal (lab) D-BINARY-ABACUS-POLYFILL-0 cost-gate candidate, kept in the research fixture. F's eager ripple built from ndarray::simd::half_add_mask; on overflow the touched columns are restored by recomputing the carry chain forward from the input mask, so the success path keeps no undo log. BITS <= 122 keeps lane values and the 64-lane total exact in u128. Tests: lanes 0/31/32/63 alone, carry chains of 1..64 columns, repeated carries, refusal leaves the counter unchanged (overflow and out-of-range significance), 24,000 random steps against 64 u128 counters with both outcomes common, and reversed traversal. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- examples/u9_domino/domino.rs | 106 +++++++++++++++++++++++++ tests/u9_domino.rs | 149 +++++++++++++++++++++++++++++++++++ 2 files changed, 255 insertions(+) diff --git a/examples/u9_domino/domino.rs b/examples/u9_domino/domino.rs index 72a17717..f7780d58 100644 --- a/examples/u9_domino/domino.rs +++ b/examples/u9_domino/domino.rs @@ -293,6 +293,112 @@ impl SubstLazy { } } +/// Why [`Abacus::try_add_at`] refused. The counter is unchanged either way. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum AbacusError { + /// `significance >= BITS`. + SignificanceOutOfRange { significance: usize, bits: usize }, + /// A carry left the highest column in these lanes. + Overflow { lanes: u64 }, +} + +/// C (D-BINARY-ABACUS-POLYFILL-0 candidate): F's eager ripple built from +/// `ndarray::simd::half_add_mask`, with transactional overflow refusal. +/// +/// Each of the 64 lanes is an independent unsigned counter of `BITS` bits. +/// `BITS <= 122`, so a lane value (< 2^122) and the total over 64 lanes +/// (< 2^128) are both exact in `u128`. +/// +/// The success path is F's. On overflow the touched columns are restored by +/// recomputing the carry chain forward from the input mask: +/// `old_k = new_k ^ m_k` and `m_{k+1} = old_k & m_k`. No undo log is kept. +#[derive(Clone)] +pub struct Abacus { + col: [u64; BITS], +} + +impl Abacus { + const BOUND: () = assert!(BITS >= 1 && BITS <= 122, "Abacus: 1 <= BITS <= 122"); + + pub fn new() -> Self { + let () = Self::BOUND; + Self { col: [0; BITS] } + } + + /// Adds `2^significance` to every lane set in `lanes`. + #[inline(always)] + pub fn try_add_at(&mut self, significance: usize, lanes: u64) -> Result<(), AbacusError> { + if significance >= BITS { + return Err(AbacusError::SignificanceOutOfRange { + significance, + bits: BITS, + }); + } + let (mut k, mut m) = (significance, lanes); + while m != 0 { + if k == BITS { + self.undo(significance, lanes); + return Err(AbacusError::Overflow { lanes: m }); + } + let (next, carry) = ndarray::simd::half_add_mask(self.col[k], m); + self.col[k] = next; + m = carry; + k += 1; + } + Ok(()) + } + + /// Restores every column from `start` up after an add of `lanes` whose + /// carry escaped the top. Cold: runs only on refusal. + #[cold] + fn undo(&mut self, start: usize, lanes: u64) { + let mut m = lanes; + for k in start..BITS { + let old = self.col[k] ^ m; + self.col[k] = old; + m &= old; + if m == 0 { + break; + } + } + } + + /// Value of lane `lane` (< 64). + pub fn lane_value(&self, lane: usize) -> u128 { + (0..BITS) + .map(|k| u128::from(self.col[k] >> lane & 1) << k) + .sum() + } + + /// Sum over all 64 lanes. + pub fn total(&self) -> u128 { + (0..BITS) + .map(|k| u128::from(self.col[k].count_ones()) << k) + .sum() + } + + pub fn columns(&self) -> [u64; BITS] { + self.col + } +} + +impl Counter for Abacus { + fn new() -> Self { + Abacus::new() + } + #[inline(always)] + fn add_at(&mut self, k: usize, m: u64) -> Result<(), Overflow> { + self.try_add_at(k, m).map_err(|e| match e { + AbacusError::Overflow { lanes } => Overflow(lanes), + AbacusError::SignificanceOutOfRange { .. } => Overflow(m), + }) + } + fn columns(&self) -> [u64; N] { + self.col + } + fn normalize(&mut self) {} +} + /// Native baseline: 64 per-lane integers. #[derive(Clone)] pub struct Native { diff --git a/tests/u9_domino.rs b/tests/u9_domino.rs index 8725e299..645fce07 100644 --- a/tests/u9_domino.rs +++ b/tests/u9_domino.rs @@ -275,3 +275,152 @@ fn xor_made_full_column_is_full_and_carries_by_state() { assert_eq!(pl.lane(j), 2, "lane {j}"); } } + +// ── D-BINARY-ABACUS-POLYFILL-0: candidate C (`Abacus`) ───────────────── + +/// 64 independent u128 counters: the reference for `Abacus`. +fn reference_add(want: &mut [u128; 64], k: usize, m: u64) { + for (j, w) in want.iter_mut().enumerate() { + *w += u128::from(m >> j & 1) << k; + } +} + +fn assert_matches(a: &Abacus, want: &[u128; 64], ctx: &str) { + for (j, &w) in want.iter().enumerate() { + assert_eq!(a.lane_value(j), w, "{ctx}: lane {j}"); + } + assert_eq!(a.total(), want.iter().sum::(), "{ctx}: total"); +} + +#[test] +fn abacus_single_lanes_stay_in_their_lane() { + for j in [0usize, 31, 32, 63] { + let mut a = Abacus::<16>::new(); + for _ in 0..1000 { + a.try_add_at(0, 1 << j).unwrap(); + } + for l in 0..64 { + assert_eq!(a.lane_value(l), if l == j { 1000 } else { 0 }, "lane {j} leaked into {l}"); + } + } +} + +#[test] +fn abacus_carry_chains_of_every_length() { + for len in [1usize, 2, 8, 16, 32, 64] { + for j in [0usize, 31, 32, 63] { + let mut a = Abacus::<96>::new(); + // Preset lane j to 2^len - 1 (columns 0..len all ones in lane j). + for k in 0..len { + a.try_add_at(k, 1 << j).unwrap(); + } + assert_eq!(a.lane_value(j), (1u128 << len) - 1); + a.try_add_at(0, 1 << j).unwrap(); + assert_eq!(a.lane_value(j), 1u128 << len, "chain {len} lane {j}"); + let cols = a.columns(); + assert_eq!(cols[len], 1 << j, "carry landed in column {len}"); + assert!(cols[..len].iter().all(|&c| c == 0), "chain {len}: lower columns cleared"); + } + } +} + +#[test] +fn abacus_repeated_carries_add_up() { + // Two carries into the same column are two units, never an XOR cancel. + let mut a = Abacus::<8>::new(); + a.try_add_at(3, u64::MAX).unwrap(); + a.try_add_at(3, u64::MAX).unwrap(); + for j in 0..64 { + assert_eq!(a.lane_value(j), 16); + } + a.try_add_at(4, u64::MAX).unwrap(); + for j in 0..64 { + assert_eq!(a.lane_value(j), 32); + } +} + +#[test] +fn abacus_overflow_refuses_and_leaves_the_counter_unchanged() { + let mut s = 0xC0FF_EE00_DEAD_BEEFu64; + let mut a = Abacus::<8>::new(); + // Lanes 0, 31, 32, 63 at 255 (top of 8 bits), the rest random below 255. + let full = (1u64 << 0) | (1 << 31) | (1 << 32) | (1 << 63); + for k in 0..8 { + a.try_add_at(k, full | (xorshift(&mut s) & !full & !(1 << 5))) + .unwrap(); + } + let before = a.columns(); + // Adding 1 to a mask containing a full lane must carry out of column 7. + for m in [full, 1 << 63, (1 << 0) | 0xF0, u64::MAX] { + let err = a.try_add_at(0, m); + assert!(matches!(err, Err(AbacusError::Overflow { lanes }) if lanes & full != 0), "{m:#x}: {err:?}"); + assert_eq!(a.columns(), before, "refused add {m:#x} changed the counter"); + } + // A lane with headroom still accepts. + a.try_add_at(0, 1 << 5).unwrap(); + assert_ne!(a.columns(), before); + // Significance past the top is refused without touching anything. + let before = a.columns(); + assert_eq!( + a.try_add_at(8, 1), + Err(AbacusError::SignificanceOutOfRange { + significance: 8, + bits: 8 + }) + ); + assert_eq!(a.columns(), before); +} + +#[test] +fn abacus_random_sequences_match_64_u128_counters() { + let mut s = 0x0123_4567_89AB_CDEFu64; + const B: usize = 12; + let (mut refused, mut accepted) = (0, 0); + for seq in 0..40 { + let mut a = Abacus::::new(); + let mut want = [0u128; 64]; + for step in 0..600 { + let r = xorshift(&mut s); + let m = match r % 6 { + 0 => 0, + 1 => u64::MAX, + 2 => 0xAAAA_AAAA_AAAA_AAAA, + 3 => 1 << (r >> 8 & 63), + _ => xorshift(&mut s), + }; + let k = (xorshift(&mut s) % B as u64) as usize; + let mut next = want; + reference_add(&mut next, k, m); + let fits = next.iter().all(|&w| w < 1 << B); + let got = a.try_add_at(k, m); + if fits { + assert_eq!(got, Ok(()), "seq {seq} step {step}"); + want = next; + accepted += 1; + } else { + assert!(matches!(got, Err(AbacusError::Overflow { .. })), "seq {seq} step {step}: {got:?}"); + refused += 1; + } + assert_matches(&a, &want, &format!("seq {seq} step {step}")); + } + } + // Anti-vacuity: both outcomes occur often. + assert!(accepted > 5000 && refused > 1000, "accepted={accepted} refused={refused}"); +} + +#[test] +fn abacus_reversed_order_gives_the_same_counters() { + // Addition is commutative, so the fold may run in either order. + let mut s = 77u64; + let adds: Vec<(usize, u64)> = (0..3000) + .map(|_| ((xorshift(&mut s) % 40) as usize, xorshift(&mut s))) + .collect(); + let (mut fwd, mut rev) = (Abacus::<72>::new(), Abacus::<72>::new()); + for &(k, m) in &adds { + fwd.try_add_at(k, m).unwrap(); + } + for &(k, m) in adds.iter().rev() { + rev.try_add_at(k, m).unwrap(); + } + assert_eq!(fwd.columns(), rev.columns()); +} From de0300103e419c298f02ab75e89f57e070aa9038 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 12:10:28 +0000 Subject: [PATCH 4/4] u9_domino: C in the counter benchmark; D-BINARY-ABACUS-POLYFILL-0 report Adds candidate C next to native and F in the counter rows, and records the cost gate: C matches F, both beat 64 per-lane adds ~4x on dense and single-lane masks. PR B (public BinaryAbacus64) and PR C (facet bit-plane helper) are deferred: no consumer yet, and the planes already cost one AND and one cmpeq_mask each on existing facade words. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- .claude/blackboard.md | 5 + .../knowledge/binary-abacus-polyfill-probe.md | 96 +++++++++++++++++++ examples/u9_domino/main.rs | 4 +- 3 files changed, 104 insertions(+), 1 deletion(-) create mode 100644 .claude/knowledge/binary-abacus-polyfill-probe.md diff --git a/.claude/blackboard.md b/.claude/blackboard.md index 3947e959..5e091424 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -4005,3 +4005,8 @@ laws, ExtentUnsupported, row bounds) and WITNESS in OrderedLaneWitness; no new t | P2 | E2 as the AVX-512-without-VNNI U9 kernel (joins #393's ADOPT-LAB row) | low | integration owner | | P3 | AVX2 P1 anomaly: setup+exec 8 ns apart, 42 ns together | low | asm of the combined path | | P4 | RPR stale cases X1-X9 / #1848 B0 fixture can reuse the A5 scalar oracle | high | lance-graph owners | + +## 2026-10-30 — D-BINARY-ABACUS-POLYFILL-0 +- PR A #402 (draft): `half_add_mask` / `full_add_mask` in `simd_masking_ops`, re-exported; vectorized `pow2n_from_int` (vsexp 4x, vserf 2.4x on AVX-512); `U32x16` `Shl`/`Shr` added to AVX2, AVX, NEON, wasm, nightly. +- Lab: `Abacus` with transactional refusal in `examples/u9_domino`; cost equal to eager F. PR B and PR C deferred: no consumer, planes already cheap with existing words. Report: `.claude/knowledge/binary-abacus-polyfill-probe.md`. +- Loose end: NEON `simd_exp_f32` has ~47 scalar fcmp/csel (pre-existing, select/compare path); AVX2 mask packing for `select` still lane by lane. diff --git a/.claude/knowledge/binary-abacus-polyfill-probe.md b/.claude/knowledge/binary-abacus-polyfill-probe.md new file mode 100644 index 00000000..755d6704 --- /dev/null +++ b/.claude/knowledge/binary-abacus-polyfill-probe.md @@ -0,0 +1,96 @@ +# D-BINARY-ABACUS-POLYFILL-0: binary positional arithmetic over the facade + +Status: PR A shipped as #402 (draft). PR B and PR C deferred, with reasons +below. Host: AVX-512 VNNI Xeon. Date: 2026-10-30. + +## What shipped (PR A, #402) + +- `ndarray::simd::half_add_mask(a, carry) = (a ^ carry, a & carry)` and + `full_add_mask(a, b, c) = (a ^ b ^ c, majority)`: the #391 `flip` and + #389 `carry_save3` laws, moved unchanged into `simd_masking_ops` and + re-exported. A `u64` is one significance column of 64 independent lanes; + the carry belongs one column up in the same lane. +- The in-flight exponent change: `pow2n_from_int` no longer goes through + `to_array`. It clamps `n` to [-126, 127], adds `2^23 + 127` so the biased + exponent sits in the mantissa field, and shifts the bits left by 23. That + needed `U32x16 << U32x16` on every backend; it existed only on AVX-512 + and scalar. Lowering, timings and disable runs are in #402. + +## Candidate C: `Abacus` (lab only, `examples/u9_domino/domino.rs`) + +- F's eager ripple, built from `half_add_mask`. +- Each of the 64 lanes is an independent unsigned counter of `BITS` bits. + `1 <= BITS <= 122`, so a lane value (< 2^122) and the 64-lane total + (< 2^128) are both exact in `u128`. +- `try_add_at(significance, lanes)` is transactional. On overflow the + touched columns are restored by recomputing the carry chain forward from + the input mask: `old_k = new_k ^ m_k`, `m_{k+1} = old_k & m_k`. The + success path keeps no undo log. A significance at or past `BITS` is + refused without touching anything. +- F (`Eager`) is NOT transactional: on overflow its lower columns stay + updated (its doc says so). C is the first counter in the fixtures that + refuses without side effects. + +Tests (`tests/u9_domino.rs`, 6 new): +- lanes 0/31/32/63 alone; +- carry chains of 1, 2, 8, 16, 32 and 64 columns in those lanes; +- repeated carries into one column add up; +- refusal on overflow and on an out-of-range significance leaves the + counter unchanged; +- 24,000 random steps against 64 `u128` counters, with zero, all-ones, + alternating, single-lane and random masks, refusals included + (accepted > 5000 and refused > 1000 asserted); +- reversed traversal gives identical columns. + +Disable runs, each red: C1 no rollback; C2 rollback stops after one +column; C3 no significance range check; C4 carry moved to the neighbouring +lane. + +## Cost gate (ns per add, 4096 adds, median of 2001 runs; range over 3 processes) + +| tier | workload | A native per-lane | B = F eager | C abacus | +|---|---|---|---|---| +| AVX-512 VNNI | dense all-ones | 8.79 – 9.08 | 2.46 – 2.66 | 2.27 – 2.35 | +| | random masks | 8.78 – 8.80 | 7.79 – 8.05 | 8.34 – 8.80 | +| | sparse one lane | 8.02 – 9.08 | 2.14 – 2.44 | 2.21 – 3.90 | +| | weighted (2^0/16/17/32) | 8.79 – 9.39 | 6.99 – 7.69 | 7.27 – 9.96 | +| AVX-512, no VNNI | dense | 8.01 – 8.78 | 2.39 – 2.45 | 2.21 – 2.22 | +| | random | 8.02 – 9.39 | 7.56 – 8.50 | 6.88 – 7.84 | +| | sparse | 8.78 – 8.79 | 2.30 – 2.35 | 3.72 – 5.26 | +| | weighted | 8.78 – 9.09 | 7.09 – 8.24 | 6.89 – 7.19 | +| AVX2 | dense | 9.04 – 9.92 | 2.25 – 2.52 | 2.33 – 2.43 | +| | random | 9.69 – 10.23 | 7.22 – 8.32 | 8.53 – 10.57 | +| | sparse | 9.91 – 10.33 | 2.09 – 2.36 | 2.29 – 3.42 | +| | weighted | 9.91 – 11.59 | 7.39 – 7.59 | 8.08 – 8.21 | + +Reading it: +- Both bit-plane counters beat 64 per-lane adds by about 4x on dense and + single-lane masks and are within 0.75 – 1.1x on random masks, the same + shape #391 found. +- C costs the same as F on dense masks. On sparse masks it is slower and + noisier (most visibly 3.7 – 5.3 vs 2.3 ns on AVX-512 without VNNI), + although the rollback never runs in these rows. Not diagnosed further. +- Not timed: AVX, SSE2 and NEON (qemu only here) and wasm (node). The + exactness tests run on every tier; timings from emulators mean nothing. +- D (U9 byte-dot) and E (vertical masked byte sum) are the existing + kernels for their consumers (#388, #390, #393). The counter is not a + replacement for either and was not imposed on them. + +## Verdicts + +- **PR B (public `BinaryAbacus64`): DEFER.** Semantics are settled and + tested, and the cost matches the eager domino. What is missing is a real + consumer: nothing in ndarray, and no lance-graph caller found, counts + 64-lane masks at chosen significances. A public type waits for one. C + stays in the fixture as the reference. +- **PR C (facet bit-plane helper): DEFER, not built.** The eight planes of + a `U8x64` are already one AND and one `cmpeq_mask` each on existing + facade words (`examples/u9_abacus/abacus.rs::planes`). #389 measured + plane preparation at 3.5 ns per block for both operands on VNNI against + the full B kernel at 2.9 ns, and reuse never broke even. There are not + two actual consumers of prepared planes to amortize across. A helper + would duplicate an existing two-word expression. +- Presence stays with the caller. The counters carry no presence; a lane + that was never touched and a lane that summed to zero both read 0, so a + consumer that needs the distinction keeps its own presence plane (as + `lut9` and Codec9 do). diff --git a/examples/u9_domino/main.rs b/examples/u9_domino/main.rs index 95b3903d..25a77089 100644 --- a/examples/u9_domino/main.rs +++ b/examples/u9_domino/main.rs @@ -132,6 +132,7 @@ fn counters(label: &str, adds: &[(usize, u64)], runs: usize) { } n.lanes.iter().map(|&x| u128::from(x)).sum() }; + let c = || run::>(adds); let (f, p, l, pl) = ( || run::>(adds), || run::>(adds), @@ -142,7 +143,7 @@ fn counters(label: &str, adds: &[(usize, u64)], runs: usize) { if fits { assert_eq!(native(), want_total, "native"); } - for (name, g) in [("F", &f as &dyn Fn() -> u128), ("P", &p), ("L", &l), ("PL", &pl)] { + for (name, g) in [("F", &f as &dyn Fn() -> u128), ("C", &c), ("P", &p), ("L", &l), ("PL", &pl)] { assert_eq!(g(), want_total, "{name}"); } println!("-- B, {label}: {} adds", adds.len()); @@ -152,6 +153,7 @@ fn counters(label: &str, adds: &[(usize, u64)], runs: usize) { None }; time("F eager", adds.len(), runs, "add", b, &f); + time("C abacus: F + transactional refusal", adds.len(), runs, "add", b, &c); time("P EMPTY substitution", adds.len(), runs, "add", b, &p); time("L lazy FULL", adds.len(), runs, "add", b, &l); time("PL both", adds.len(), runs, "add", b, &pl);