From dca345d154d2fee8eb4960861252dcf4b74cc23a Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 13:45:52 +0000 Subject: [PATCH 1/6] =?UTF-8?q?simd:=20NaN/presence=20bijection=20?= =?UTF-8?q?=E2=80=94=20not=5Fnan=5F*=5Fto=5Fmask,=20nan=5Fwhere=5Fclear=5F?= =?UTF-8?q?*,=20presence=5Fnan=5Fmismatches=5F*?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit On a float lane a clear presence bit and a NaN payload are one state: presence bit i == 0 <=> values[i] is NaN. Three primitives hold it: not_nan_{f64,f32}_to_mask (lane -> mask, ordered simd_eq against itself), nan_where_clear_{f64,f32} (mask -> lane, NaN stamped under clear bits) and presence_nan_mismatches_{f64,f32} (the witness: popcount(presence XOR derived) over live rows, no allocation). F64Mask8::to_bitmask added on the AVX-512, AVX2/AVX, NEON, wasm and scalar backends. Parity group 0x13xx checks all three against a bit-pattern NaN reference (quiet/signalling/negative/full-payload NaN, +-inf, +-0, subnormal, MAX) at every lane-group boundary, plus the round trip and anti-vacuity counts. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- .../plans/nan-exclusion-mindmap-scratch.md | 37 +- crates/simd-masking-parity/src/lib.rs | 161 ++++++- src/simd.rs | 6 + src/simd_avx2.rs | 5 + src/simd_avx512.rs | 5 + src/simd_masking_ops.rs | 426 ++++++++++++++++++ src/simd_neon.rs | 5 + src/simd_scalar.rs | 11 +- src/simd_wasm.rs | 5 + 9 files changed, 654 insertions(+), 7 deletions(-) diff --git a/.claude/plans/nan-exclusion-mindmap-scratch.md b/.claude/plans/nan-exclusion-mindmap-scratch.md index ee1368b9..df526783 100644 --- a/.claude/plans/nan-exclusion-mindmap-scratch.md +++ b/.claude/plans/nan-exclusion-mindmap-scratch.md @@ -1,7 +1,12 @@ # NaN exclusion for masked statistics: scratch mindmap -> **Status: SCRATCH, for later examination.** Nothing here is built, measured -> or decided. It maps purposes against the benefits and costs of each option, +> **Status: PARTLY BUILT (2026-10-27).** The NaN / presence bijection is +> shipped: `not_nan_{f64,f32}_to_mask` (lane → mask), `nan_where_clear_{f64,f32}` +> (mask → lane) and `presence_nan_mismatches_{f64,f32}` (the witness), with +> parity group `0x13xx` on every realization. See § "Decision: presence IS NaN" +> at the end. The masked float fold (S2) and float-output determinism are +> still open. The rest of this page is the original scratch: it maps purposes +> against the benefits and costs of each option. It maps purposes against the benefits and costs of each option, > so a later session can pick a design without re-deriving the trade-offs. > Context: the follow-up to #378 (wasm `mul_add` rounds once on every > realization) and lance-graph `D-WFL-DET` (replay determinism). @@ -212,3 +217,31 @@ are struck. **Council status:** halted before Phase 1; nothing was implemented. Q2 alone is probably too small for a council (two mask primitives and a convention). Q1 and Q3 are not ndarray's to decide. + +## Decision: presence IS NaN (operator, 2026-10-27) + +On a float lane the presence bit and NaN are one state with two spellings, +held as a bijection: + +``` +presence bit i == 0 ⇔ values[i] is NaN (for every live row i) +``` + +This replaces the premise audit's Q1/Q2 split for float lanes. "Missing" and +"arithmetic produced NaN" are deliberately merged: a fold can read either +spelling and gets the same population, and `n` is the popcount of the plane it +folds, so dilution cannot occur. "Is NaN" is the predicate `x != x`, never a +bit pattern, so backends that do not preserve NaN payloads (wasm) cannot break +it. Integer lanes have no NaN; there the presence plane stays the only carrier. + +Shipped primitives (`src/simd_masking_ops.rs`, re-exported from `ndarray::simd`): + +| direction | function | cost | +|---|---|---| +| lane → mask | `not_nan_f64_to_mask`, `not_nan_f32_to_mask` | one ordered `simd_eq` per 8 / 16 lanes | +| mask → lane | `nan_where_clear_f64`, `nan_where_clear_f32` | one word per 64 rows plus one store per clear bit | +| witness | `presence_nan_mismatches_f64`, `presence_nan_mismatches_f32` | `popcount(presence XOR derived)` over live rows, no allocation | + +`F64Mask8::to_bitmask` was added on the AVX-512, AVX2/AVX, NEON, wasm and +scalar backends (nightly already had it) so the f64 packer reads one +spelling on every arm. diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs index b1d51969..e625fbc4 100644 --- a/crates/simd-masking-parity/src/lib.rs +++ b/crates/simd-masking-parity/src/lib.rs @@ -36,7 +36,9 @@ //! U8x64`) and the `F32x16` ↔ `I32x16` casts; `0x11xx` `F32x16` / `F64x8` //! `mul_add` rounding once, bit for bit against an exact integer reference (no //! `mul_add`, so no libm in the binary); `0x12xx` the `I32x16` gemmlowp requantization -//! methods against an independent `i128` reference. `main.rs` (native / qemu) and +//! methods against an independent `i128` reference; `0x13xx` the NaN / presence +//! bijection (`not_nan_*_to_mask`, `nan_where_clear_*`, `presence_nan_mismatches_*`) +//! against a bit-pattern NaN reference. `main.rs` (native / qemu) and //! `selfcheck()` (the wasm cdylib export, driven by `run.mjs`) both call //! [`run`]. @@ -50,14 +52,15 @@ use ndarray::simd::{ mask_shift_morton, mask_ternlog, mask_ternlog_any, mask_ternlog_assign, mask_ternlog_popcount, mask_xor, mask_xor_assign, masked_group_sum_i32, masked_group_sum_i32_via, masked_key_run_count_u32, masked_max_i32, masked_min_i32, masked_strided_group_sum, masked_sum_i32, masked_sum_wrapping_add_i32, ne_i32_to_mask, - ne_i32_to_mask_under, ne_u32_to_mask, ne_u32_to_mask_under, ne_u64_to_mask, ne_u8_to_mask, + nan_where_clear_f32, nan_where_clear_f64, ne_i32_to_mask_under, ne_u32_to_mask, ne_u32_to_mask_under, ne_u64_to_mask, + ne_u8_to_mask, not_nan_f32_to_mask, not_nan_f64_to_mask, presence_nan_mismatches_f32, presence_nan_mismatches_f64, ternary_match_strided16_to_mask, ternary_match_strided_to_mask, ternary_match_u32_to_mask, ternary_match_u32_to_mask_under, ternary_match_u64_to_mask, ternary_match_u64_to_mask_under, ternlog, F32x16, F64x8, I16x32, I32x16, I8x16, KeyRunCarry, MortonDir, U32x16, U64x8, U8x16, }; /// Number of check groups [`run`] executes (for the log line only). -pub const CHECKS: usize = 19; +pub const CHECKS: usize = 20; /// The wasm export: identical to [`run`], `extern "C"` so `run.mjs` can call it. #[no_mangle] @@ -72,7 +75,7 @@ pub fn run() -> u32 { check_mask_algebra, check_care_match, check_masked_reductions, check_blend, check_morton_shift, check_predicates_under, check_set_range, check_unsigned_compare_to_mask, check_gather_scatter_group, check_f64x8_compare, check_i8x16_u8x16_lanes, check_care_match16, check_narrow_ladder, - check_mul_add_single_rounding, check_fixedpoint_requant, + check_mul_add_single_rounding, check_fixedpoint_requant, check_nan_presence_bijection, ]; for g in groups { if let Err(code) = g() { @@ -2713,3 +2716,153 @@ mod oracle_tests { } } } + +// ── 0x13xx: the NaN / presence bijection ───────────────────────────────────── +// +// `presence bit i == 0 ⇔ values[i] is NaN`, held by three primitives: the +// lane → mask derivation (`not_nan_*_to_mask`), the mask → lane stamp +// (`nan_where_clear_*`) and the witness (`presence_nan_mismatches_*`). The +// reference classifies NaN from the BIT PATTERN (exponent all ones, mantissa +// non-zero) and never compares a float, so it shares nothing with the +// realization's ordered `simd_eq`. Edge values: quiet and signalling-pattern +// NaN, negative NaN, full payload, ±inf (exponent all ones but NOT NaN), ±0, +// the smallest subnormal, MAX. Lengths cover every lane-group boundary of the +// 8-lane f64 and 16-lane f32 packers. + +const NAN_LENS: [usize; 12] = [0, 1, 7, 8, 15, 16, 17, 63, 64, 65, 130, 200]; + +fn nan_bits_f64(b: u64) -> bool { + (b >> 52) & 0x7FF == 0x7FF && b & ((1u64 << 52) - 1) != 0 +} + +fn nan_bits_f32(b: u32) -> bool { + (b >> 23) & 0xFF == 0xFF && b & ((1u32 << 23) - 1) != 0 +} + +fn nan_lane_f64(n: usize, rng: &mut SplitMix64) -> Vec { + const EDGES: [u64; 11] = [ + 0x7FF8_0000_0000_0000, 0xFFF8_0000_0000_0000, 0x7FF0_0000_0000_0001, 0x7FFF_FFFF_FFFF_FFFF, + 0x7FF0_0000_0000_0000, 0xFFF0_0000_0000_0000, 0x0000_0000_0000_0000, 0x8000_0000_0000_0000, + 0x0000_0000_0000_0001, 0x7FEF_FFFF_FFFF_FFFF, 0x3FF0_0000_0000_0000, + ]; + (0..n) + .map(|_| { + let r = rng.next(); + if r % 3 == 0 { f64::from_bits(EDGES[(r >> 8) as usize % EDGES.len()]) } else { f64::from_bits(rng.next()) } + }) + .collect() +} + +fn nan_lane_f32(n: usize, rng: &mut SplitMix64) -> Vec { + const EDGES: [u32; 11] = [ + 0x7FC0_0000, 0xFFC0_0000, 0x7F80_0001, 0x7FFF_FFFF, 0x7F80_0000, 0xFF80_0000, 0x0000_0000, 0x8000_0000, + 0x0000_0001, 0x7F7F_FFFF, 0x3F80_0000, + ]; + (0..n) + .map(|_| { + let r = rng.next(); + if r % 3 == 0 { f32::from_bits(EDGES[(r >> 8) as usize % EDGES.len()]) } else { f32::from_bits(rng.next() as u32) } + }) + .collect() +} + +fn check_nan_presence_bijection() -> Result<(), u32> { + let mut rng = SplitMix64(0x0BAD_C0DE_F00D_0013); + let mut nan_rows = 0usize; + let mut witnessed = 0u64; + for (li, &n) in NAN_LENS.iter().enumerate() { + let code = li as u32; + let words = words_for(n); + let v64 = nan_lane_f64(n, &mut rng); + let v32 = nan_lane_f32(n, &mut rng); + + // lane → mask, against the bit-pattern reference; surplus word and tail zero. + let mut got = vec![u64::MAX; words + 1]; + not_nan_f64_to_mask(&v64, &mut got); + if got != reference_mask(n, words + 1, |i| !nan_bits_f64(v64[i].to_bits())) { + return Err(0x1300 | code); + } + nan_rows += v64.iter().filter(|x| nan_bits_f64(x.to_bits())).count(); + let mut got32 = vec![u64::MAX; words + 1]; + not_nan_f32_to_mask(&v32, &mut got32); + if got32 != reference_mask(n, words + 1, |i| !nan_bits_f32(v32[i].to_bits())) { + return Err(0x1310 | code); + } + + // A random presence plane with set bits past `n` (dead bits). + let mut presence: Vec = (0..words + 1).map(|_| rng.next()).collect(); + if n % 64 != 0 { + presence[n / 64] |= !0u64 << (n % 64); + } + let bit = |p: &[u64], i: usize| (p[i / 64] >> (i % 64)) & 1 == 1; + + // Witness, against the reference count over live rows only. + let want = (0..n).filter(|&i| bit(&presence, i) == nan_bits_f64(v64[i].to_bits())).count() as u64; + if presence_nan_mismatches_f64(&v64, &presence) != want { + return Err(0x1340 | code); + } + witnessed += want; + let want32 = (0..n).filter(|&i| bit(&presence, i) == nan_bits_f32(v32[i].to_bits())).count() as u64; + if presence_nan_mismatches_f32(&v32, &presence) != want32 { + return Err(0x1350 | code); + } + + // mask → lane: NaN exactly under clear bits, set rows bit-identical. + let mut s64 = v64.clone(); + nan_where_clear_f64(&presence, &mut s64); + for i in 0..n { + let ok = if bit(&presence, i) { s64[i].to_bits() == v64[i].to_bits() } else { nan_bits_f64(s64[i].to_bits()) }; + if !ok { + return Err(0x1320 | code); + } + } + let mut s32 = v32.clone(); + nan_where_clear_f32(&presence, &mut s32); + for i in 0..n { + let ok = if bit(&presence, i) { s32[i].to_bits() == v32[i].to_bits() } else { nan_bits_f32(s32[i].to_bits()) }; + if !ok { + return Err(0x1330 | code); + } + } + + // The bijection: a derived plane agrees with its lane; stamping NaN under + // (random AND derived) makes the lane agree with that plane, and + // re-deriving returns the plane bit for bit. + let mut derived = vec![0u64; words]; + not_nan_f64_to_mask(&v64, &mut derived); + if presence_nan_mismatches_f64(&v64, &derived) != 0 { + return Err(0x1360 | code); + } + let mut plane = vec![0u64; words]; + mask_and(&presence[..words], &derived, &mut plane); + let mut stamped = v64.clone(); + nan_where_clear_f64(&plane, &mut stamped); + let mut again = vec![0u64; words]; + not_nan_f64_to_mask(&stamped, &mut again); + if presence_nan_mismatches_f64(&stamped, &plane) != 0 || again != plane { + return Err(0x1370 | code); + } + let mut derived32 = vec![0u64; words]; + not_nan_f32_to_mask(&v32, &mut derived32); + let mut plane32 = vec![0u64; words]; + mask_and(&presence[..words], &derived32, &mut plane32); + let mut stamped32 = v32.clone(); + nan_where_clear_f32(&plane32, &mut stamped32); + let mut again32 = vec![0u64; words]; + not_nan_f32_to_mask(&stamped32, &mut again32); + if presence_nan_mismatches_f32(&v32, &derived32) != 0 + || presence_nan_mismatches_f32(&stamped32, &plane32) != 0 + || again32 != plane32 + { + return Err(0x1380 | code); + } + } + // Anti-vacuity: NaN rows of the edge kinds exist, and the witness fires. + if nan_rows < 40 { + return Err(0x13F0); + } + if witnessed < 100 { + return Err(0x13F1); + } + Ok(()) +} diff --git a/src/simd.rs b/src/simd.rs index e028b232..9334ee67 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -962,12 +962,18 @@ pub use crate::simd_masking_ops::{ masked_strided_group_sum, masked_sum_i32, masked_sum_wrapping_add_i32, + nan_where_clear_f32, + nan_where_clear_f64, ne_i32_to_mask, ne_i32_to_mask_under, ne_u32_to_mask, ne_u32_to_mask_under, ne_u64_to_mask, ne_u8_to_mask, + not_nan_f32_to_mask, + not_nan_f64_to_mask, + presence_nan_mismatches_f32, + presence_nan_mismatches_f64, ternary_match_strided16_to_mask, ternary_match_strided_to_mask, ternary_match_u32_to_mask, diff --git a/src/simd_avx2.rs b/src/simd_avx2.rs index 20ed7997..2818e681 100644 --- a/src/simd_avx2.rs +++ b/src/simd_avx2.rs @@ -1144,6 +1144,11 @@ impl Default for F64x8 { #[derive(Copy, Clone, Debug)] pub struct F64Mask8(pub u8); impl F64Mask8 { + /// Lane `i` of the comparison is bit `i` of the result. + #[inline(always)] + pub fn to_bitmask(self) -> u8 { + self.0 + } #[inline(always)] pub fn select(self, true_val: F64x8, false_val: F64x8) -> F64x8 { let t = true_val.to_array(); diff --git a/src/simd_avx512.rs b/src/simd_avx512.rs index cc77d6e2..641b4bf1 100644 --- a/src/simd_avx512.rs +++ b/src/simd_avx512.rs @@ -507,6 +507,11 @@ impl PartialEq for F64x8 { pub struct F64Mask8(pub __mmask8); impl F64Mask8 { + /// Lane `i` of the comparison is bit `i` of the result. + #[inline(always)] + pub fn to_bitmask(self) -> u8 { + self.0 + } #[inline(always)] pub fn select(self, true_val: F64x8, false_val: F64x8) -> F64x8 { F64x8(unsafe { _mm512_mask_blend_pd(self.0, false_val.0, true_val.0) }) diff --git a/src/simd_masking_ops.rs b/src/simd_masking_ops.rs index 66f191f9..0791d4f5 100644 --- a/src/simd_masking_ops.rs +++ b/src/simd_masking_ops.rs @@ -3999,6 +3999,235 @@ pub fn ne_u64_to_mask(values: &[u64], needle: u64, out_words: &mut [u64]) { clear_mask_tail(out_words, values.len()); } +// ── The NaN / presence bijection ───────────────────────────────────────────── +// +// On a float lane, "this row is absent" and "this value is NaN" are ONE state +// with two spellings: a clear presence bit and a NaN payload. The invariant is +// +// presence bit i == 0 ⇔ values[i] is NaN (for every i < len) +// +// Three primitives keep it: one per direction and one witness. +// +// * lane → mask: [`not_nan_f64_to_mask`] / [`not_nan_f32_to_mask`] derive the +// presence plane from the lane (`x == x` is false only for NaN). +// * mask → lane: [`nan_where_clear_f64`] / [`nan_where_clear_f32`] stamp NaN +// into every row whose presence bit is clear. +// * witness: [`presence_nan_mismatches_f64`] / [`presence_nan_mismatches_f32`] +// count the rows where the two spellings disagree; `0` is the invariant. +// +// "Is NaN" is the IEEE predicate `x != x`, never a bit pattern: every NaN +// counts (quiet, signalling-pattern, either sign, any payload), so a backend +// that does not preserve NaN payloads (wasm) cannot break the bijection. The +// compare is the facade's ordered `simd_eq` on every realization, so `±0`, +// `±inf`, subnormals and `MAX` are all present. +// +// Holding the invariant merges "missing" and "arithmetic produced NaN" into +// one state. That is the point: a masked fold can read either spelling and +// get the same population, and `n` is the popcount of the plane it folds. + +/// Packs `values[i] == values[i]` — "is not NaN", the presence bit of a float +/// lane under the NaN/presence bijection — into `out_words`, one bit per +/// element, LSB-first within each `u64` word. +/// +/// `out_words` is **fully overwritten**. Trailing bits in the final word +/// beyond `values.len()`, and any surplus words past `ceil(len / 64)`, are +/// written as `0`. +/// +/// Every NaN (quiet or signalling pattern, either sign, any payload) packs +/// `0`; every other value (`±0`, `±inf`, subnormals, `MAX`) packs `1`. +/// +/// Runs 8 lanes at a time through [`crate::simd::F64x8::simd_eq`], the +/// ordered compare, against itself; the final partial group is zero-padded +/// into one register and its padding bits masked off. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::not_nan_f64_to_mask; +/// +/// let values = [1.0, f64::NAN, -0.0, f64::INFINITY]; +/// let mut words = [0u64; 1]; +/// not_nan_f64_to_mask(&values, &mut words); +/// // only element 1 is NaN → bits 0, 2, 3 → 0b1101 +/// assert_eq!(words[0], 0b1101); +/// ``` +#[inline] +pub fn not_nan_f64_to_mask(values: &[f64], out_words: &mut [u64]) { + pack::("not_nan_f64_to_mask", values, out_words, |lanes, live| { + let v = crate::simd::F64x8::from_array(lanes); + (v.simd_eq(v).to_bitmask() & live8(live)) as u64 + }); +} + +/// Packs `values[i] == values[i]` — "is not NaN" — for an `f32` lane. The +/// `f32` sibling of [`not_nan_f64_to_mask`]; same contract. +/// +/// Runs 16 lanes at a time through [`crate::simd::F32x16::simd_eq`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::not_nan_f32_to_mask; +/// +/// let values = [f32::NAN, 2.5, f32::NEG_INFINITY, f32::NAN]; +/// let mut words = [0u64; 1]; +/// not_nan_f32_to_mask(&values, &mut words); +/// assert_eq!(words[0], 0b0110); +/// ``` +#[inline] +pub fn not_nan_f32_to_mask(values: &[f32], out_words: &mut [u64]) { + pack::("not_nan_f32_to_mask", values, out_words, |lanes, live| { + let v = crate::simd::F32x16::from_array(lanes); + (v.simd_eq(v).to_bitmask() & live16(live)) as u64 + }); +} + +/// Writes `f64::NAN` into every row of `values` whose `presence` bit is +/// clear: the mask → lane direction of the NaN/presence bijection. Rows whose +/// bit is set are not touched, NaN or not. +/// +/// After this call every absent row holds NaN. Together with a presence plane +/// that has no set bit over a NaN row (for example one built by +/// [`not_nan_f64_to_mask`]), [`presence_nan_mismatches_f64`] is `0`. +/// +/// Bits of `presence` beyond `values.len()`, and surplus words, are ignored. +/// Cost is one word read per 64 rows plus one store per clear bit: it walks +/// the clear bits with `trailing_zeros`, so a fully present word costs one +/// compare. This is a mask-gated write of a constant, with nothing to +/// vectorize per lane and nothing a backend could realize differently. +/// +/// # Panics +/// +/// Panics if `presence.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::nan_where_clear_f64; +/// +/// let mut values = [1.0, 9.9, 3.0, 4.0]; +/// nan_where_clear_f64(&[0b1101], &mut values); // row 1 is absent +/// assert!(values[1].is_nan()); +/// assert_eq!([values[0], values[2], values[3]], [1.0, 3.0, 4.0]); +/// ``` +#[inline] +pub fn nan_where_clear_f64(presence: &[u64], values: &mut [f64]) { + nan_where_clear("nan_where_clear_f64", presence, values, f64::NAN); +} + +/// Writes `f32::NAN` into every row of `values` whose `presence` bit is +/// clear. The `f32` sibling of [`nan_where_clear_f64`]; same contract. +/// +/// # Panics +/// +/// Panics if `presence.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::nan_where_clear_f32; +/// +/// let mut values = [1.0f32, 2.0, 3.0]; +/// nan_where_clear_f32(&[0b011], &mut values); +/// assert!(values[2].is_nan()); +/// ``` +#[inline] +pub fn nan_where_clear_f32(presence: &[u64], values: &mut [f32]) { + nan_where_clear("nan_where_clear_f32", presence, values, f32::NAN); +} + +/// The shared mask → lane writer: `nan` into every row whose bit is clear. +#[inline(always)] +fn nan_where_clear(name: &str, presence: &[u64], values: &mut [T], nan: T) { + let words = mask_words_for(values.len()); + assert!(presence.len() >= words, "{name}: presence.len()={} < required {words}", presence.len()); + for (k, chunk) in values.chunks_mut(64).enumerate() { + let mut absent = !presence[k] & live64(chunk.len()); + while absent != 0 { + chunk[absent.trailing_zeros() as usize] = nan; + absent &= absent - 1; + } + } +} + +/// Counts the rows of a float lane where the NaN/presence bijection does not +/// hold: `presence` bit `i` set while `values[i]` is NaN, or bit `i` clear +/// while `values[i]` is not NaN. `0` means the two spellings agree on every +/// row; it is the invariant check and the parity witness. +/// +/// Computed as `popcount(presence XOR not_nan(values))` over the live rows, +/// 64 rows at a time through [`not_nan_f64_to_mask`], with no allocation. +/// Bits of `presence` beyond `values.len()`, and surplus words, are ignored. +/// +/// # Panics +/// +/// Panics if `presence.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::{nan_where_clear_f64, presence_nan_mismatches_f64}; +/// +/// let presence = [0b1011u64]; // row 2 absent +/// let mut values = [1.0, f64::NAN, 7.0, 4.0]; // row 1: NaN under a set bit +/// // Rows 1 (present but NaN) and 2 (absent but not NaN) disagree. +/// assert_eq!(presence_nan_mismatches_f64(&values, &presence), 2); +/// nan_where_clear_f64(&presence, &mut values); +/// // Row 2 now holds NaN; row 1 still contradicts its set bit. +/// assert_eq!(presence_nan_mismatches_f64(&values, &presence), 1); +/// ``` +#[inline] +pub fn presence_nan_mismatches_f64(values: &[f64], presence: &[u64]) -> u64 { + presence_nan_mismatches("presence_nan_mismatches_f64", values, presence, not_nan_f64_to_mask) +} + +/// Counts the rows of an `f32` lane where the NaN/presence bijection does +/// not hold. The `f32` sibling of [`presence_nan_mismatches_f64`]; same +/// contract. +/// +/// # Panics +/// +/// Panics if `presence.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::{not_nan_f32_to_mask, presence_nan_mismatches_f32}; +/// +/// let values = [f32::NAN, 1.0, 2.0]; +/// let mut presence = [0u64; 1]; +/// not_nan_f32_to_mask(&values, &mut presence); +/// // A plane derived from the lane always agrees with it. +/// assert_eq!(presence_nan_mismatches_f32(&values, &presence), 0); +/// ``` +#[inline] +pub fn presence_nan_mismatches_f32(values: &[f32], presence: &[u64]) -> u64 { + presence_nan_mismatches("presence_nan_mismatches_f32", values, presence, not_nan_f32_to_mask) +} + +/// The shared witness: `popcount(presence XOR derived)` over the live rows, +/// one 64-row word at a time so no scratch plane is allocated. +#[inline(always)] +fn presence_nan_mismatches(name: &str, values: &[T], presence: &[u64], derive: fn(&[T], &mut [u64])) -> u64 { + let words = mask_words_for(values.len()); + assert!(presence.len() >= words, "{name}: presence.len()={} < required {words}", presence.len()); + let mut derived = [0u64; 1]; + let mut total = 0u64; + for (k, chunk) in values.chunks(64).enumerate() { + derive(chunk, &mut derived); + total += ((derived[0] ^ presence[k]) & live64(chunk.len())).count_ones() as u64; + } + total +} + /// `dst = !src` over `n_rows` elements — the tail-aware complement. Bits at /// or past `n_rows` are written `0`, so a conforming input yields a /// conforming output (plain `!` on the words would set every tail bit). @@ -5529,6 +5758,203 @@ pub fn mask_shift_morton(src: &[u64], dir: MortonDir, dst: &mut [u64]) { mod tests { use super::*; + // ── NaN / presence bijection ──────────────────────────────────────────── + // + // The oracle classifies NaN from the bit pattern (exponent all ones, + // mantissa non-zero), never with a float compare, so it shares nothing + // with the `simd_eq` path under test. + + fn is_nan_bits_f64(x: f64) -> bool { + let b = x.to_bits(); + (b >> 52) & 0x7FF == 0x7FF && b & ((1u64 << 52) - 1) != 0 + } + + fn is_nan_bits_f32(x: f32) -> bool { + let b = x.to_bits(); + (b >> 23) & 0xFF == 0xFF && b & ((1u32 << 23) - 1) != 0 + } + + const F64_EDGES: [u64; 12] = [ + 0x7FF8_0000_0000_0000, // quiet NaN + 0xFFF8_0000_0000_0000, // negative quiet NaN + 0x7FF0_0000_0000_0001, // signalling-pattern NaN, smallest payload + 0x7FFF_FFFF_FFFF_FFFF, // NaN, all payload bits + 0x7FF0_0000_0000_0000, // +inf (exponent all ones, mantissa zero: NOT NaN) + 0xFFF0_0000_0000_0000, // -inf + 0x0000_0000_0000_0000, // +0 + 0x8000_0000_0000_0000, // -0 + 0x0000_0000_0000_0001, // smallest subnormal + 0x7FEF_FFFF_FFFF_FFFF, // MAX + 0x3FF0_0000_0000_0000, // 1.0 + 0xBFF0_0000_0000_0000, // -1.0 + ]; + + const F32_EDGES: [u32; 12] = [ + 0x7FC0_0000, 0xFFC0_0000, 0x7F80_0001, 0x7FFF_FFFF, 0x7F80_0000, 0xFF80_0000, 0x0000_0000, 0x8000_0000, + 0x0000_0001, 0x7F7F_FFFF, 0x3F80_0000, 0xBF80_0000, + ]; + + fn f64_lane(n: usize, seed: u64) -> Vec { + let mut s = seed; + (0..n) + .map(|i| { + s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + let pick = (s >> 33) as usize; + if pick % 3 == 0 { + f64::from_bits(F64_EDGES[(i + pick) % F64_EDGES.len()]) + } else { + // Arbitrary bit patterns: about 1 in 2048 is a NaN on its own. + f64::from_bits(s) + } + }) + .collect() + } + + fn f32_lane(n: usize, seed: u64) -> Vec { + let mut s = seed; + (0..n) + .map(|i| { + s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + let pick = (s >> 33) as usize; + if pick % 3 == 0 { + f32::from_bits(F32_EDGES[(i + pick) % F32_EDGES.len()]) + } else { + f32::from_bits((s >> 16) as u32) + } + }) + .collect() + } + + const NAN_LENS: [usize; 9] = [0, 1, 7, 8, 15, 63, 64, 65, 200]; + + #[test] + fn not_nan_to_mask_matches_bit_pattern_oracle_with_clean_tail() { + let mut nans = 0usize; + for (s, &n) in NAN_LENS.iter().enumerate() { + let words = n.div_ceil(64) + 1; // one surplus word + let v64 = f64_lane(n, 0xA5A5 + s as u64); + let mut got = vec![u64::MAX; words]; + not_nan_f64_to_mask(&v64, &mut got); + let mut want = vec![0u64; words]; + for (i, &x) in v64.iter().enumerate() { + if !is_nan_bits_f64(x) { + want[i / 64] |= 1 << (i % 64); + } else { + nans += 1; + } + } + assert_eq!(got, want, "f64 n={n}"); + + let v32 = f32_lane(n, 0x5A5A + s as u64); + let mut got = vec![u64::MAX; words]; + not_nan_f32_to_mask(&v32, &mut got); + let mut want = vec![0u64; words]; + for (i, &x) in v32.iter().enumerate() { + if !is_nan_bits_f32(x) { + want[i / 64] |= 1 << (i % 64); + } + } + assert_eq!(got, want, "f32 n={n}"); + } + // Anti-vacuity: the corpus really contains NaNs of the edge kinds. + assert!(nans >= 20, "only {nans} NaN rows in the f64 corpus"); + } + + #[test] + fn nan_where_clear_writes_nan_exactly_at_clear_bits() { + for (s, &n) in NAN_LENS.iter().enumerate() { + let original = f64_lane(n, 0x1234 + s as u64); + // Random presence with garbage in the bits past `n` and a surplus word. + let mut presence: Vec = (0..n.div_ceil(64) + 1) + .map(|w| 0x9E37_79B9_7F4A_7C15u64.rotate_left((w as u32 + s as u32) * 7)) + .collect(); + if n % 64 != 0 { + presence[n / 64] |= !0u64 << (n % 64); + } + let mut values = original.clone(); + nan_where_clear_f64(&presence, &mut values); + for i in 0..n { + let set = (presence[i / 64] >> (i % 64)) & 1 == 1; + if set { + assert_eq!(values[i].to_bits(), original[i].to_bits(), "set row {i} touched, n={n}"); + } else { + assert!(values[i].is_nan(), "clear row {i} not NaN, n={n}"); + } + } + + let original32 = f32_lane(n, 0x4321 + s as u64); + let mut values32 = original32.clone(); + nan_where_clear_f32(&presence, &mut values32); + for i in 0..n { + let set = (presence[i / 64] >> (i % 64)) & 1 == 1; + if set { + assert_eq!(values32[i].to_bits(), original32[i].to_bits()); + } else { + assert!(values32[i].is_nan()); + } + } + } + } + + #[test] + fn mismatches_count_every_disagreement_and_ignore_dead_bits() { + let mut fired = 0u64; + for (s, &n) in NAN_LENS.iter().enumerate() { + let values = f64_lane(n, 0xBEEF + s as u64); + let mut presence: Vec = (0..n.div_ceil(64) + 1) + .map(|w| 0xD1B5_4A32_D192_ED03u64.rotate_left((w as u32 * 13) + s as u32)) + .collect(); + if n % 64 != 0 { + presence[n / 64] |= !0u64 << (n % 64); // dead bits must not count + } + let want = (0..n) + .filter(|&i| ((presence[i / 64] >> (i % 64)) & 1 == 1) == is_nan_bits_f64(values[i])) + .count() as u64; + assert_eq!(presence_nan_mismatches_f64(&values, &presence), want, "n={n}"); + fired += want; + + let values32 = f32_lane(n, 0xCAFE + s as u64); + let want32 = (0..n) + .filter(|&i| ((presence[i / 64] >> (i % 64)) & 1 == 1) == is_nan_bits_f32(values32[i])) + .count() as u64; + assert_eq!(presence_nan_mismatches_f32(&values32, &presence), want32, "f32 n={n}"); + } + // Can-fire: random planes disagree with random lanes somewhere. + assert!(fired > 100, "witness fired only {fired} times"); + } + + #[test] + fn bijection_holds_after_either_direction_and_stays_silent() { + for (s, &n) in NAN_LENS.iter().enumerate() { + let values = f64_lane(n, 0x7777 + s as u64); + // lane → mask: a derived plane always agrees with its lane. + let mut derived = vec![0u64; n.div_ceil(64)]; + not_nan_f64_to_mask(&values, &mut derived); + assert_eq!(presence_nan_mismatches_f64(&values, &derived), 0, "derived n={n}"); + + // mask → lane: stamping NaN under clear bits of a plane that has no + // set bit over a NaN row makes the lane agree with that plane. + let random: Vec = (0..n.div_ceil(64)) + .map(|w| 0xF00D_F00D_1234_5678u64.rotate_left(w as u32 * 5 + s as u32)) + .collect(); + let mut plane = vec![0u64; n.div_ceil(64)]; + mask_and(&random, &derived, &mut plane); + let mut stamped = values.clone(); + nan_where_clear_f64(&plane, &mut stamped); + assert_eq!(presence_nan_mismatches_f64(&stamped, &plane), 0, "stamped n={n}"); + // Re-deriving from the stamped lane gives the plane back. + let mut again = vec![0u64; n.div_ceil(64)]; + not_nan_f64_to_mask(&stamped, &mut again); + assert_eq!(again, plane, "round trip n={n}"); + } + } + + #[test] + #[should_panic(expected = "presence.len()")] + fn mismatches_refuse_a_short_presence_plane() { + presence_nan_mismatches_f64(&[1.0; 65], &[0u64; 1]); + } + // ── Packed-bitmask predicates + mask algebra ──────────────────────────── // // Every test compares the shipped path against an INDEPENDENT scalar diff --git a/src/simd_neon.rs b/src/simd_neon.rs index bfc8a81f..62b25a0c 100644 --- a/src/simd_neon.rs +++ b/src/simd_neon.rs @@ -1200,6 +1200,11 @@ pub mod aarch64_simd { #[derive(Copy, Clone, Debug)] pub struct F64Mask8(pub u8); impl F64Mask8 { + /// Lane `i` of the comparison is bit `i` of the result. + #[inline(always)] + pub fn to_bitmask(self) -> u8 { + self.0 + } #[inline(always)] pub fn select(self, true_val: F64x8, false_val: F64x8) -> F64x8 { let t = true_val.to_array(); diff --git a/src/simd_scalar.rs b/src/simd_scalar.rs index 6a049d3c..c7fd501b 100644 --- a/src/simd_scalar.rs +++ b/src/simd_scalar.rs @@ -1181,7 +1181,8 @@ impl Neg for I32x16 { // x86-shared API on the mask: `to_bitmask` exists on every x86 arm (and NEON); // the x86-64-v2 arm takes this realization, so it needs it too. Only on -// `F32Mask16`, as on x86 — `F64Mask8` has no `to_bitmask` there. +// `F32Mask16`; `F64Mask8` below carries it on every arm since the NaN/presence +// bijection (`not_nan_f64_to_mask`) reads it. impl F32Mask16 { /// Lane `i` of the comparison is bit `i` of the result. /// @@ -1199,6 +1200,14 @@ impl F32Mask16 { } } +impl F64Mask8 { + /// Lane `i` of the comparison is bit `i` of the result. + #[inline(always)] + pub fn to_bitmask(self) -> u8 { + self.0 + } +} + // Extra for F32x16: to_bits/from_bits/cast_i32 impl F32x16 { /// Gather 16 `f32` values from `base_ptr` at the signed element offsets in diff --git a/src/simd_wasm.rs b/src/simd_wasm.rs index 46200966..3bd33f83 100644 --- a/src/simd_wasm.rs +++ b/src/simd_wasm.rs @@ -790,6 +790,11 @@ pub mod wasm32_simd { #[derive(Copy, Clone, Debug)] pub struct F64Mask8(pub u8); impl F64Mask8 { + /// Lane `i` of the comparison is bit `i` of the result. + #[inline(always)] + pub fn to_bitmask(self) -> u8 { + self.0 + } #[inline(always)] pub fn select(self, true_val: F64x8, false_val: F64x8) -> F64x8 { let t = true_val.to_array(); From ee8171fa85bc5d728515ff606b450de1deb2c787 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 14:07:09 +0000 Subject: [PATCH 2/6] simd: SQL comparisons over a NULL=NaN float lane + mask_not_under (NOT under presence) {lt,le,gt,ge,eq,ne}_{f64,f32}_to_mask pack the SQL TRUE plane: a NaN row is NULL and never matches (<> included, unlike IEEE !=), and a NaN threshold is the NULL literal and matches nothing. mask_not_under(pred, presence) = presence & !pred is SQL NOT under three-valued logic: the FALSE plane, with UNKNOWN rows in neither. presence is the PREDICATE's presence (AND of its operands'), empty for a NULL literal; the unit test caught the column's plane being wrong there (x > NULL reported every present row FALSE). Parity group 0x139xx-0x13Cxx checks every operator x 9 thresholds (NaN, +-0, +-inf, subnormal, +-1, MAX) against a bit-pattern ordering, plus the FALSE plane. Also: clippy is_multiple_of fixes in the bijection tests. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- crates/simd-masking-parity/src/lib.rs | 203 ++++++++++-- src/simd.rs | 13 + src/simd_masking_ops.rs | 427 +++++++++++++++++++++++++- 3 files changed, 614 insertions(+), 29 deletions(-) diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs index e625fbc4..a8b41fa8 100644 --- a/crates/simd-masking-parity/src/lib.rs +++ b/crates/simd-masking-parity/src/lib.rs @@ -38,29 +38,34 @@ //! `mul_add`, so no libm in the binary); `0x12xx` the `I32x16` gemmlowp requantization //! methods against an independent `i128` reference; `0x13xx` the NaN / presence //! bijection (`not_nan_*_to_mask`, `nan_where_clear_*`, `presence_nan_mismatches_*`) -//! against a bit-pattern NaN reference. `main.rs` (native / qemu) and +//! against a bit-pattern NaN reference; `0x139xx`–`0x13Cxx` the SQL comparisons over +//! a NULL=NaN lane and `mask_not_under` (SQL `NOT` under presence) against a +//! bit-pattern ordering. `main.rs` (native / qemu) and //! `selfcheck()` (the wasm cdylib export, driven by `run.mjs`) both call //! [`run`]. use ndarray::simd::{ - blend_i32, eq_i32_to_mask, eq_i32_to_mask_under, eq_u32_strided_to_mask, eq_u32_to_mask, eq_u32_to_mask_under, - eq_u32_via_to_mask, eq_u64_to_mask, eq_u8_to_mask, ge_i32_to_mask, ge_i32_to_mask_under, ge_u64_to_mask, - ge_u8_to_mask, gt_i32_to_mask, gt_i32_to_mask_under, gt_u64_to_mask, gt_u8_to_mask, le_i32_to_mask, - le_i32_to_mask_under, le_u64_to_mask, le_u8_to_mask, lt_i32_to_mask, lt_i32_to_mask_under, lt_u64_to_mask, - lt_u8_to_mask, mask_all, mask_and, mask_and_assign, mask_andnot, mask_andnot_assign, mask_any, mask_gather_u32, - mask_gather_u32_under, mask_not, mask_not_assign, mask_or, mask_or_assign, mask_scatter_or_u16, mask_scatter_or_u32, mask_set_range, + blend_i32, eq_f32_to_mask, eq_f64_to_mask, eq_i32_to_mask, eq_i32_to_mask_under, eq_u32_strided_to_mask, + eq_u32_to_mask, eq_u32_to_mask_under, eq_u32_via_to_mask, eq_u64_to_mask, eq_u8_to_mask, ge_f32_to_mask, + ge_f64_to_mask, ge_i32_to_mask, ge_i32_to_mask_under, ge_u64_to_mask, ge_u8_to_mask, gt_f32_to_mask, + gt_f64_to_mask, gt_i32_to_mask, gt_i32_to_mask_under, gt_u64_to_mask, gt_u8_to_mask, le_f32_to_mask, + le_f64_to_mask, le_i32_to_mask, le_i32_to_mask_under, le_u64_to_mask, le_u8_to_mask, lt_f32_to_mask, + lt_f64_to_mask, lt_i32_to_mask, lt_i32_to_mask_under, lt_u64_to_mask, lt_u8_to_mask, mask_all, mask_and, + mask_and_assign, mask_andnot, mask_andnot_assign, mask_any, mask_gather_u32, mask_gather_u32_under, mask_not, + mask_not_assign, mask_not_under, mask_or, mask_or_assign, mask_scatter_or_u16, mask_scatter_or_u32, mask_set_range, mask_shift_morton, mask_ternlog, mask_ternlog_any, mask_ternlog_assign, mask_ternlog_popcount, mask_xor, mask_xor_assign, masked_group_sum_i32, masked_group_sum_i32_via, masked_key_run_count_u32, masked_max_i32, - masked_min_i32, masked_strided_group_sum, masked_sum_i32, masked_sum_wrapping_add_i32, ne_i32_to_mask, - nan_where_clear_f32, nan_where_clear_f64, ne_i32_to_mask_under, ne_u32_to_mask, ne_u32_to_mask_under, ne_u64_to_mask, - ne_u8_to_mask, not_nan_f32_to_mask, not_nan_f64_to_mask, presence_nan_mismatches_f32, presence_nan_mismatches_f64, - ternary_match_strided16_to_mask, ternary_match_strided_to_mask, ternary_match_u32_to_mask, - ternary_match_u32_to_mask_under, ternary_match_u64_to_mask, ternary_match_u64_to_mask_under, ternlog, F32x16, - F64x8, I16x32, I32x16, I8x16, KeyRunCarry, MortonDir, U32x16, U64x8, U8x16, + masked_min_i32, masked_strided_group_sum, masked_sum_i32, masked_sum_wrapping_add_i32, nan_where_clear_f32, + nan_where_clear_f64, ne_f32_to_mask, ne_f64_to_mask, ne_i32_to_mask, ne_i32_to_mask_under, ne_u32_to_mask, + ne_u32_to_mask_under, ne_u64_to_mask, ne_u8_to_mask, not_nan_f32_to_mask, not_nan_f64_to_mask, + presence_nan_mismatches_f32, presence_nan_mismatches_f64, ternary_match_strided16_to_mask, + ternary_match_strided_to_mask, ternary_match_u32_to_mask, ternary_match_u32_to_mask_under, + ternary_match_u64_to_mask, ternary_match_u64_to_mask_under, ternlog, F32x16, F64x8, I16x32, I32x16, I8x16, + KeyRunCarry, MortonDir, U32x16, U64x8, U8x16, }; /// Number of check groups [`run`] executes (for the log line only). -pub const CHECKS: usize = 20; +pub const CHECKS: usize = 21; /// The wasm export: identical to [`run`], `extern "C"` so `run.mjs` can call it. #[no_mangle] @@ -76,6 +81,7 @@ pub fn run() -> u32 { check_predicates_under, check_set_range, check_unsigned_compare_to_mask, check_gather_scatter_group, check_f64x8_compare, check_i8x16_u8x16_lanes, check_care_match16, check_narrow_ladder, check_mul_add_single_rounding, check_fixedpoint_requant, check_nan_presence_bijection, + check_sql_null_predicates, ]; for g in groups { if let Err(code) = g() { @@ -954,7 +960,7 @@ fn check_mask_algebra() -> Result<(), u32> { if mask_all(&missing, n) || (n > 1 && !mask_any(&missing)) { return Err(0x624); } - if n % 64 != 0 { + if !n.is_multiple_of(64) { let mut tail_only = zero.clone(); tail_only[nw - 1] = 1u64 << (n % 64); if !mask_any(&tail_only) || mask_all(&tail_only, n) { @@ -2748,7 +2754,11 @@ fn nan_lane_f64(n: usize, rng: &mut SplitMix64) -> Vec { (0..n) .map(|_| { let r = rng.next(); - if r % 3 == 0 { f64::from_bits(EDGES[(r >> 8) as usize % EDGES.len()]) } else { f64::from_bits(rng.next()) } + if r.is_multiple_of(3) { + f64::from_bits(EDGES[(r >> 8) as usize % EDGES.len()]) + } else { + f64::from_bits(rng.next()) + } }) .collect() } @@ -2761,7 +2771,11 @@ fn nan_lane_f32(n: usize, rng: &mut SplitMix64) -> Vec { (0..n) .map(|_| { let r = rng.next(); - if r % 3 == 0 { f32::from_bits(EDGES[(r >> 8) as usize % EDGES.len()]) } else { f32::from_bits(rng.next() as u32) } + if r.is_multiple_of(3) { + f32::from_bits(EDGES[(r >> 8) as usize % EDGES.len()]) + } else { + f32::from_bits(rng.next() as u32) + } }) .collect() } @@ -2791,18 +2805,22 @@ fn check_nan_presence_bijection() -> Result<(), u32> { // A random presence plane with set bits past `n` (dead bits). let mut presence: Vec = (0..words + 1).map(|_| rng.next()).collect(); - if n % 64 != 0 { + if !n.is_multiple_of(64) { presence[n / 64] |= !0u64 << (n % 64); } let bit = |p: &[u64], i: usize| (p[i / 64] >> (i % 64)) & 1 == 1; // Witness, against the reference count over live rows only. - let want = (0..n).filter(|&i| bit(&presence, i) == nan_bits_f64(v64[i].to_bits())).count() as u64; + let want = (0..n) + .filter(|&i| bit(&presence, i) == nan_bits_f64(v64[i].to_bits())) + .count() as u64; if presence_nan_mismatches_f64(&v64, &presence) != want { return Err(0x1340 | code); } witnessed += want; - let want32 = (0..n).filter(|&i| bit(&presence, i) == nan_bits_f32(v32[i].to_bits())).count() as u64; + let want32 = (0..n) + .filter(|&i| bit(&presence, i) == nan_bits_f32(v32[i].to_bits())) + .count() as u64; if presence_nan_mismatches_f32(&v32, &presence) != want32 { return Err(0x1350 | code); } @@ -2811,7 +2829,11 @@ fn check_nan_presence_bijection() -> Result<(), u32> { let mut s64 = v64.clone(); nan_where_clear_f64(&presence, &mut s64); for i in 0..n { - let ok = if bit(&presence, i) { s64[i].to_bits() == v64[i].to_bits() } else { nan_bits_f64(s64[i].to_bits()) }; + let ok = if bit(&presence, i) { + s64[i].to_bits() == v64[i].to_bits() + } else { + nan_bits_f64(s64[i].to_bits()) + }; if !ok { return Err(0x1320 | code); } @@ -2819,7 +2841,11 @@ fn check_nan_presence_bijection() -> Result<(), u32> { let mut s32 = v32.clone(); nan_where_clear_f32(&presence, &mut s32); for i in 0..n { - let ok = if bit(&presence, i) { s32[i].to_bits() == v32[i].to_bits() } else { nan_bits_f32(s32[i].to_bits()) }; + let ok = if bit(&presence, i) { + s32[i].to_bits() == v32[i].to_bits() + } else { + nan_bits_f32(s32[i].to_bits()) + }; if !ok { return Err(0x1330 | code); } @@ -2866,3 +2892,136 @@ fn check_nan_presence_bijection() -> Result<(), u32> { } Ok(()) } + +// ── 0x139x–0x13Cx: SQL comparisons over a NULL=NaN lane, and NOT under presence ─ +// +// The reference orders floats by their BIT PATTERN (sign-magnitude to a signed +// key, `-0` and `+0` both key 0) and treats a NaN pattern as NULL, so it never +// performs a float compare. TRUE(x op t) = neither is NULL and key(x) op key(t); +// FALSE = neither is NULL and not key(x) op key(t). `mask_not_under(T, presence)` +// must equal FALSE exactly when `presence` is the PREDICATE's presence (the +// column's, or empty for a NULL threshold): a bare complement would wrongly set +// every NULL row, and the column's plane under a NULL threshold every present one. + +fn key64(b: u64) -> Option { + if nan_bits_f64(b) { + return None; + } + let mag = (b & 0x7FFF_FFFF_FFFF_FFFF) as i128; + Some(if b >> 63 == 1 { -mag } else { mag }) +} + +fn key32(b: u32) -> Option { + if nan_bits_f32(b) { + return None; + } + let mag = (b & 0x7FFF_FFFF) as i128; + Some(if b >> 31 == 1 { -mag } else { mag }) +} + +/// `op` 0..6 = lt, le, gt, ge, eq, ne over the reference keys. +fn ref_rel(op: usize, a: i128, b: i128) -> bool { + match op { + 0 => a < b, + 1 => a <= b, + 2 => a > b, + 3 => a >= b, + 4 => a == b, + _ => a != b, + } +} + +fn check_sql_null_predicates() -> Result<(), u32> { + type F64Op = fn(&[f64], f64, &mut [u64]); + type F32Op = fn(&[f32], f32, &mut [u64]); + let ops64: [F64Op; 6] = + [lt_f64_to_mask, le_f64_to_mask, gt_f64_to_mask, ge_f64_to_mask, eq_f64_to_mask, ne_f64_to_mask]; + let ops32: [F32Op; 6] = + [lt_f32_to_mask, le_f32_to_mask, gt_f32_to_mask, ge_f32_to_mask, eq_f32_to_mask, ne_f32_to_mask]; + let thr64: [u64; 9] = [ + 0x7FF8_0000_0000_0000, 0x0000_0000_0000_0000, 0x8000_0000_0000_0000, 0x7FF0_0000_0000_0000, + 0xFFF0_0000_0000_0000, 0x0000_0000_0000_0001, 0x3FF0_0000_0000_0000, 0xBFF0_0000_0000_0000, + 0x7FEF_FFFF_FFFF_FFFF, + ]; + let thr32: [u32; 9] = [ + 0x7FC0_0000, 0x0000_0000, 0x8000_0000, 0x7F80_0000, 0xFF80_0000, 0x0000_0001, 0x3F80_0000, 0xBF80_0000, + 0x7F7F_FFFF, + ]; + let mut rng = SplitMix64(0x5A1_0000_0000_0013); + let mut ieee_not_would_leak = 0usize; + let mut true_bits = 0usize; + for &n in NAN_LENS.iter() { + let words = words_for(n); + let mut v64 = nan_lane_f64(n, &mut rng); + let mut v32 = nan_lane_f32(n, &mut rng); + // Plant the thresholds into the lane so `=` and the boundaries of `<=`/`>=` fire. + for i in (0..n).step_by(5) { + v64[i] = f64::from_bits(thr64[i % thr64.len()]); + v32[i] = f32::from_bits(thr32[i % thr32.len()]); + } + let mut present = vec![0u64; words]; + not_nan_f64_to_mask(&v64, &mut present); + let mut present32 = vec![0u64; words]; + not_nan_f32_to_mask(&v32, &mut present32); + for (ti, (&t64, &t32)) in thr64.iter().zip(&thr32).enumerate() { + for op in 0..6 { + let code = ((ti as u32) << 4) | op as u32; + // f64 TRUE plane, with a surplus word that must come back zero. + let mut got = vec![u64::MAX; words + 1]; + ops64[op](&v64, f64::from_bits(t64), &mut got); + let truth = |i: usize| match (key64(v64[i].to_bits()), key64(t64)) { + (Some(a), Some(b)) => Some(ref_rel(op, a, b)), + _ => None, + }; + if got != reference_mask(n, words + 1, |i| truth(i) == Some(true)) { + return Err(0x1_3900 | code); + } + true_bits += (0..n).filter(|&i| truth(i) == Some(true)).count(); + // NOT under the predicate's presence (empty for a NULL + // threshold) is exactly the FALSE plane. + let pred_present = if key64(t64).is_some() { + present.clone() + } else { + vec![0u64; words] + }; + let mut neg = vec![0u64; words]; + mask_not_under(&got[..words], &pred_present, &mut neg); + if neg != reference_mask(n, words, |i| truth(i) == Some(false)) { + return Err(0x1_3B00 | code); + } + // Rows a bare complement would wrongly report as FALSE: NULL rows. + ieee_not_would_leak += (0..n) + .filter(|&i| truth(i).is_none() && key64(t64).is_some()) + .count(); + + let mut got32 = vec![u64::MAX; words + 1]; + ops32[op](&v32, f32::from_bits(t32), &mut got32); + let truth32 = |i: usize| match (key32(v32[i].to_bits()), key32(t32)) { + (Some(a), Some(b)) => Some(ref_rel(op, a, b)), + _ => None, + }; + if got32 != reference_mask(n, words + 1, |i| truth32(i) == Some(true)) { + return Err(0x1_3A00 | code); + } + let pred_present32 = if key32(t32).is_some() { + present32.clone() + } else { + vec![0u64; words] + }; + let mut neg32 = vec![0u64; words]; + mask_not_under(&got32[..words], &pred_present32, &mut neg32); + if neg32 != reference_mask(n, words, |i| truth32(i) == Some(false)) { + return Err(0x1_3C00 | code); + } + } + } + } + // Anti-vacuity: TRUE rows exist, and NULL rows exist where a bare NOT would leak. + if true_bits < 1000 { + return Err(0x13F2); + } + if ieee_not_would_leak < 500 { + return Err(0x13F3); + } + Ok(()) +} diff --git a/src/simd.rs b/src/simd.rs index 9334ee67..a64112d4 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -881,6 +881,8 @@ 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, + eq_f32_to_mask, + eq_f64_to_mask, eq_i32_to_mask, eq_i32_to_mask_under, eq_u32_strided_to_mask, @@ -891,18 +893,26 @@ pub use crate::simd_masking_ops::{ eq_u8_to_mask, fold_bounded_cross_power_sums_tiles, fold_bounded_power_sums_tiles, + ge_f32_to_mask, + ge_f64_to_mask, ge_i32_to_mask, ge_i32_to_mask_under, ge_u64_to_mask, ge_u8_to_mask, + gt_f32_to_mask, + gt_f64_to_mask, gt_i32_to_mask, gt_i32_to_mask_under, gt_u64_to_mask, gt_u8_to_mask, + le_f32_to_mask, + le_f64_to_mask, le_i32_to_mask, le_i32_to_mask_under, le_u64_to_mask, le_u8_to_mask, + lt_f32_to_mask, + lt_f64_to_mask, lt_i32_to_mask, lt_i32_to_mask_under, lt_u64_to_mask, @@ -917,6 +927,7 @@ pub use crate::simd_masking_ops::{ mask_gather_u32_under, mask_not, mask_not_assign, + mask_not_under, mask_or, mask_or_assign, mask_scatter_or_u16, @@ -964,6 +975,8 @@ pub use crate::simd_masking_ops::{ masked_sum_wrapping_add_i32, nan_where_clear_f32, nan_where_clear_f64, + ne_f32_to_mask, + ne_f64_to_mask, ne_i32_to_mask, ne_i32_to_mask_under, ne_u32_to_mask, diff --git a/src/simd_masking_ops.rs b/src/simd_masking_ops.rs index 0791d4f5..b5d6fa5c 100644 --- a/src/simd_masking_ops.rs +++ b/src/simd_masking_ops.rs @@ -4228,6 +4228,353 @@ fn presence_nan_mismatches(name: &str, values: &[T], presence: &[u64], derive total } +// ── SQL comparisons over a NULL=NaN float lane ────────────────────────────── +// +// Under the NaN/presence bijection a NaN row is a SQL NULL, and a comparison +// with NULL is UNKNOWN, never TRUE. Each `*_f64_to_mask` / `*_f32_to_mask` +// below packs the TRUE plane of `values[i] threshold`: +// +// * a NaN row packs `0` for every operator, `<>` included (IEEE says +// `NaN != t` is true; SQL says `NULL <> t` is UNKNOWN); +// * a NaN threshold packs `0` everywhere (`x NULL` is UNKNOWN); +// * `-0.0 = 0.0` is TRUE, as in both IEEE and SQL. +// +// NOT is where IEEE and SQL disagree: `!(NaN > 5)` is true, but +// `NOT (NULL > 5)` is UNKNOWN. So never complement a TRUE plane with a bare +// `mask_not`; negate it inside the presence plane with [`mask_not_under`] +// (`presence & !pred`), which yields SQL's FALSE plane and leaves UNKNOWN rows +// out of both. + +/// Shared packer for the f64 SQL comparisons: `rel(x, t)` is the ordered IEEE +/// relation, ANDed with `x == x` so a NaN (NULL) row never matches. A NaN +/// threshold makes every row UNKNOWN, so the result is all zeros. +#[inline(always)] +fn sql_cmp_f64( + name: &str, values: &[f64], threshold: f64, out_words: &mut [u64], + rel: impl Fn(crate::simd::F64x8, crate::simd::F64x8) -> u8, +) { + let t = crate::simd::F64x8::splat(threshold); + let unknown = threshold.is_nan(); + pack::(name, values, out_words, |lanes, live| { + if unknown { + return 0; + } + let v = crate::simd::F64x8::from_array(lanes); + (rel(v, t) & v.simd_eq(v).to_bitmask() & live8(live)) as u64 + }); +} + +/// Shared packer for the f32 SQL comparisons; see [`sql_cmp_f64`]. +#[inline(always)] +fn sql_cmp_f32( + name: &str, values: &[f32], threshold: f32, out_words: &mut [u64], + rel: impl Fn(crate::simd::F32x16, crate::simd::F32x16) -> u16, +) { + let t = crate::simd::F32x16::splat(threshold); + let unknown = threshold.is_nan(); + pack::(name, values, out_words, |lanes, live| { + if unknown { + return 0; + } + let v = crate::simd::F32x16::from_array(lanes); + (rel(v, t) & v.simd_eq(v).to_bitmask() & live16(live)) as u64 + }); +} + +/// Packs the SQL TRUE plane of `values[i] < threshold` over a NULL=NaN `f64` +/// lane: bit `i` is set iff `values[i]` is not NaN, `threshold` is not NaN and +/// `values[i] < threshold`. Full overwrite, trailing bits zero. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::lt_f64_to_mask; +/// +/// let values = [1.0, f64::NAN, 9.0, -0.0]; +/// let mut words = [0u64; 1]; +/// lt_f64_to_mask(&values, 5.0, &mut words); +/// assert_eq!(words[0], 0b1001); // the NaN row is UNKNOWN, not TRUE +/// ``` +#[inline] +pub fn lt_f64_to_mask(values: &[f64], threshold: f64, out_words: &mut [u64]) { + sql_cmp_f64("lt_f64_to_mask", values, threshold, out_words, |v, t| v.simd_lt(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] <= threshold` over a NULL=NaN `f64` +/// lane. Same contract as [`lt_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::le_f64_to_mask; +/// +/// let mut words = [0u64; 1]; +/// le_f64_to_mask(&[5.0, f64::NAN, 6.0], 5.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn le_f64_to_mask(values: &[f64], threshold: f64, out_words: &mut [u64]) { + sql_cmp_f64("le_f64_to_mask", values, threshold, out_words, |v, t| v.simd_le(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] > threshold` over a NULL=NaN `f64` +/// lane. Same contract as [`lt_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::gt_f64_to_mask; +/// +/// let mut words = [0u64; 1]; +/// gt_f64_to_mask(&[9.0, f64::NAN, f64::INFINITY], 5.0, &mut words); +/// assert_eq!(words[0], 0b101); +/// ``` +#[inline] +pub fn gt_f64_to_mask(values: &[f64], threshold: f64, out_words: &mut [u64]) { + sql_cmp_f64("gt_f64_to_mask", values, threshold, out_words, |v, t| v.simd_gt(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] >= threshold` over a NULL=NaN `f64` +/// lane. Same contract as [`lt_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::ge_f64_to_mask; +/// +/// let mut words = [0u64; 1]; +/// ge_f64_to_mask(&[5.0, f64::NAN, 4.0], 5.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn ge_f64_to_mask(values: &[f64], threshold: f64, out_words: &mut [u64]) { + sql_cmp_f64("ge_f64_to_mask", values, threshold, out_words, |v, t| v.simd_ge(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] = needle` over a NULL=NaN `f64` +/// lane. `-0.0 = 0.0` is TRUE; a NaN row or a NaN needle never matches. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::eq_f64_to_mask; +/// +/// let mut words = [0u64; 1]; +/// eq_f64_to_mask(&[-0.0, f64::NAN, 1.0], 0.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn eq_f64_to_mask(values: &[f64], needle: f64, out_words: &mut [u64]) { + sql_cmp_f64("eq_f64_to_mask", values, needle, out_words, |v, t| v.simd_eq(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] <> needle` over a NULL=NaN `f64` +/// lane. Unlike IEEE `!=`, a NaN row packs `0`: `NULL <> needle` is UNKNOWN. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::ne_f64_to_mask; +/// +/// let mut words = [0u64; 1]; +/// ne_f64_to_mask(&[1.0, f64::NAN, 2.0], 2.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn ne_f64_to_mask(values: &[f64], needle: f64, out_words: &mut [u64]) { + sql_cmp_f64("ne_f64_to_mask", values, needle, out_words, |v, t| v.simd_ne(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] < threshold` over a NULL=NaN `f32` +/// lane. The `f32` sibling of [`lt_f64_to_mask`]; same contract. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::lt_f32_to_mask; +/// +/// let mut words = [0u64; 1]; +/// lt_f32_to_mask(&[1.0, f32::NAN, 9.0], 5.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn lt_f32_to_mask(values: &[f32], threshold: f32, out_words: &mut [u64]) { + sql_cmp_f32("lt_f32_to_mask", values, threshold, out_words, |v, t| v.simd_lt(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] <= threshold` over a NULL=NaN `f32` +/// lane. Same contract as [`lt_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::le_f32_to_mask; +/// +/// let mut words = [0u64; 1]; +/// le_f32_to_mask(&[5.0, f32::NAN, 6.0], 5.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn le_f32_to_mask(values: &[f32], threshold: f32, out_words: &mut [u64]) { + sql_cmp_f32("le_f32_to_mask", values, threshold, out_words, |v, t| v.simd_le(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] > threshold` over a NULL=NaN `f32` +/// lane. Same contract as [`lt_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::gt_f32_to_mask; +/// +/// let mut words = [0u64; 1]; +/// gt_f32_to_mask(&[9.0, f32::NAN, 1.0], 5.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn gt_f32_to_mask(values: &[f32], threshold: f32, out_words: &mut [u64]) { + sql_cmp_f32("gt_f32_to_mask", values, threshold, out_words, |v, t| v.simd_gt(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] >= threshold` over a NULL=NaN `f32` +/// lane. Same contract as [`lt_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::ge_f32_to_mask; +/// +/// let mut words = [0u64; 1]; +/// ge_f32_to_mask(&[5.0, f32::NAN, 4.0], 5.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn ge_f32_to_mask(values: &[f32], threshold: f32, out_words: &mut [u64]) { + sql_cmp_f32("ge_f32_to_mask", values, threshold, out_words, |v, t| v.simd_ge(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] = needle` over a NULL=NaN `f32` +/// lane. Same contract as [`eq_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::eq_f32_to_mask; +/// +/// let mut words = [0u64; 1]; +/// eq_f32_to_mask(&[-0.0, f32::NAN, 1.0], 0.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn eq_f32_to_mask(values: &[f32], needle: f32, out_words: &mut [u64]) { + sql_cmp_f32("eq_f32_to_mask", values, needle, out_words, |v, t| v.simd_eq(t).to_bitmask()); +} + +/// Packs the SQL TRUE plane of `values[i] <> needle` over a NULL=NaN `f32` +/// lane. Same contract as [`ne_f64_to_mask`]. +/// +/// # Panics +/// +/// Panics if `out_words.len() < values.len().div_ceil(64)`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::ne_f32_to_mask; +/// +/// let mut words = [0u64; 1]; +/// ne_f32_to_mask(&[1.0, f32::NAN, 2.0], 2.0, &mut words); +/// assert_eq!(words[0], 0b001); +/// ``` +#[inline] +pub fn ne_f32_to_mask(values: &[f32], needle: f32, out_words: &mut [u64]) { + sql_cmp_f32("ne_f32_to_mask", values, needle, out_words, |v, t| v.simd_ne(t).to_bitmask()); +} + +/// SQL `NOT` under three-valued logic: `dst = presence & !pred`. +/// +/// `pred` is a TRUE plane (for example from [`gt_f64_to_mask`]) and +/// `presence` is **the predicate's** presence: the rows where every operand +/// is non-NULL, i.e. the AND of the operands' presence planes (for a +/// NULL=NaN float lane, [`not_nan_f64_to_mask`]; for any other column, its +/// validity plane). A NULL literal operand (a NaN threshold) makes it +/// all-zero, since `x > NULL` is UNKNOWN on every row; passing the column's +/// plane there would report every present row FALSE. The result is the FALSE +/// plane: rows where the predicate was evaluated and is false. Rows whose +/// presence bit is clear are UNKNOWN and stay out of both planes, which is +/// what SQL requires and what a bare [`mask_not`] gets wrong (`!(NaN > 5)` +/// is true in IEEE, `NOT (NULL > 5)` is UNKNOWN in SQL). +/// +/// It is [`mask_andnot`]`(presence, pred)` under the name of the rule it +/// carries, so a lowering that negates a nullable predicate cannot reach for +/// the bare complement by accident. Tail bits are zero whenever `presence`'s +/// tail is zero. +/// +/// # Panics +/// +/// Panics unless `pred.len() == presence.len() == dst.len()`. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::{gt_f64_to_mask, mask_not_under, not_nan_f64_to_mask}; +/// +/// let x = [9.0, f64::NAN, 1.0]; +/// let (mut gt, mut present, mut not_gt) = ([0u64; 1], [0u64; 1], [0u64; 1]); +/// gt_f64_to_mask(&x, 5.0, &mut gt); // WHERE x > 5 → row 0 +/// not_nan_f64_to_mask(&x, &mut present); +/// mask_not_under(>, &present, &mut not_gt); // WHERE NOT x > 5 → row 2 only +/// assert_eq!((gt[0], not_gt[0]), (0b001, 0b100)); // the NULL row is in neither +/// ``` +#[inline] +pub fn mask_not_under(pred: &[u64], presence: &[u64], dst: &mut [u64]) { + mask_andnot(presence, pred, dst); +} + /// `dst = !src` over `n_rows` elements — the tail-aware complement. Bits at /// or past `n_rows` are written `0`, so a conforming input yields a /// conforming output (plain `!` on the words would set every tail bit). @@ -5798,9 +6145,11 @@ mod tests { let mut s = seed; (0..n) .map(|i| { - s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + s = s + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); let pick = (s >> 33) as usize; - if pick % 3 == 0 { + if pick.is_multiple_of(3) { f64::from_bits(F64_EDGES[(i + pick) % F64_EDGES.len()]) } else { // Arbitrary bit patterns: about 1 in 2048 is a NaN on its own. @@ -5814,9 +6163,11 @@ mod tests { let mut s = seed; (0..n) .map(|i| { - s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + s = s + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); let pick = (s >> 33) as usize; - if pick % 3 == 0 { + if pick.is_multiple_of(3) { f32::from_bits(F32_EDGES[(i + pick) % F32_EDGES.len()]) } else { f32::from_bits((s >> 16) as u32) @@ -5868,7 +6219,7 @@ mod tests { let mut presence: Vec = (0..n.div_ceil(64) + 1) .map(|w| 0x9E37_79B9_7F4A_7C15u64.rotate_left((w as u32 + s as u32) * 7)) .collect(); - if n % 64 != 0 { + if !n.is_multiple_of(64) { presence[n / 64] |= !0u64 << (n % 64); } let mut values = original.clone(); @@ -5904,7 +6255,7 @@ mod tests { let mut presence: Vec = (0..n.div_ceil(64) + 1) .map(|w| 0xD1B5_4A32_D192_ED03u64.rotate_left((w as u32 * 13) + s as u32)) .collect(); - if n % 64 != 0 { + if !n.is_multiple_of(64) { presence[n / 64] |= !0u64 << (n % 64); // dead bits must not count } let want = (0..n) @@ -5949,6 +6300,68 @@ mod tests { } } + fn sql_key_f64(x: f64) -> Option { + let b = x.to_bits(); + if is_nan_bits_f64(x) { + return None; + } + let mag = (b & 0x7FFF_FFFF_FFFF_FFFF) as i128; + Some(if b >> 63 == 1 { -mag } else { mag }) + } + + #[test] + fn sql_predicates_are_true_planes_and_not_under_is_the_false_plane() { + type Op = fn(&[f64], f64, &mut [u64]); + type Rel = fn(i128, i128) -> bool; + let ops: [(Op, Rel); 6] = [ + (lt_f64_to_mask, |a, b| a < b), + (le_f64_to_mask, |a, b| a <= b), + (gt_f64_to_mask, |a, b| a > b), + (ge_f64_to_mask, |a, b| a >= b), + (eq_f64_to_mask, |a, b| a == b), + (ne_f64_to_mask, |a, b| a != b), + ]; + let mut leaks = 0usize; + for (s, &n) in NAN_LENS.iter().enumerate() { + let mut values = f64_lane(n, 0x5E1 + s as u64); + for i in (0..n).step_by(4) { + values[i] = f64::from_bits(F64_EDGES[i % F64_EDGES.len()]); + } + let words = n.div_ceil(64); + let mut column = vec![0u64; words]; + not_nan_f64_to_mask(&values, &mut column); + for &tb in F64_EDGES.iter() { + let t = f64::from_bits(tb); + // The predicate's presence: the column's, and nothing at all + // when the threshold is the NULL literal. + let present = if t.is_nan() { vec![0u64; words] } else { column.clone() }; + for (op, rel) in ops { + let truth = |i: usize| match (sql_key_f64(values[i]), sql_key_f64(t)) { + (Some(a), Some(b)) => Some(rel(a, b)), + _ => None, + }; + let mut got = vec![u64::MAX; words + 1]; + op(&values, t, &mut got); + let mut want = vec![0u64; words + 1]; + let mut not_want = vec![0u64; words]; + for i in 0..n { + match truth(i) { + Some(true) => want[i / 64] |= 1 << (i % 64), + Some(false) => not_want[i / 64] |= 1 << (i % 64), + None => leaks += usize::from(sql_key_f64(t).is_some()), + } + } + assert_eq!(got, want, "TRUE plane n={n} t={t:?}"); + let mut neg = vec![0u64; words]; + mask_not_under(&got[..words], &present, &mut neg); + assert_eq!(neg, not_want, "FALSE plane n={n} t={t:?}"); + } + } + } + // Anti-vacuity: NULL rows a bare complement would have reported as FALSE. + assert!(leaks > 100, "only {leaks} NULL rows exercised the NOT rule"); + } + #[test] #[should_panic(expected = "presence.len()")] fn mismatches_refuse_a_short_presence_plane() { @@ -7218,7 +7631,7 @@ mod tests { let src_bits: Vec = (0..n).map(|_| splitmix(&mut seed) & 1 == 1).collect(); let mut src = bits_to_words(&src_bits); // Dirty tail bits past n must be ignored by both. - if n % 64 != 0 { + if !n.is_multiple_of(64) { if let Some(last) = src.last_mut() { *last |= !0u64 << (n % 64); } From a6803bc84a8775f4477cdc917e3c02873df607a0 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 14:18:57 +0000 Subject: [PATCH 3/6] board: NaN/presence bijection + SQL NULL=NaN predicates (blackboard, plan status) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- .claude/blackboard.md | 9 +++++++++ .claude/plans/nan-exclusion-mindmap-scratch.md | 13 +++++++++++++ 2 files changed, 22 insertions(+) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index deff2874..316a522b 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -1,3 +1,12 @@ +## 2026-10-27 — NaN / presence bijection + SQL NULL=NaN predicates + +- **Decision (operator):** on a float lane, presence IS NaN, held as a bijection: `presence bit i == 0 ⇔ values[i] is NaN` (predicate `x != x`, never a bit pattern, so wasm NaN canonicalization cannot break it). Missing and arithmetic NaN are deliberately merged. Integer lanes keep the presence plane as the only carrier. Supersedes the premise audit's Q1/Q2 split for float lanes (`.claude/plans/nan-exclusion-mindmap-scratch.md`). +- **Shipped (`simd_masking_ops.rs`, re-exported from `ndarray::simd`):** `not_nan_{f64,f32}_to_mask` (lane → mask, ordered `simd_eq` against itself), `nan_where_clear_{f64,f32}` (mask → lane, NaN under clear bits, cost ∝ clear bits), `presence_nan_mismatches_{f64,f32}` (witness, `popcount(presence XOR derived)` over live rows, no allocation), `{lt,le,gt,ge,eq,ne}_{f64,f32}_to_mask` (SQL TRUE plane: a NaN row never matches, `<>` included; a NaN threshold matches nothing), `mask_not_under(pred, presence)` = `presence & !pred` (SQL NOT, the FALSE plane). `F64Mask8::to_bitmask` added on AVX-512, AVX2/AVX, NEON, wasm and scalar (nightly had it). +- **Finding from the unit test:** `presence` for NOT must be the PREDICATE's presence (AND of its operands'), empty for a NULL literal. With the column's plane, `NOT (x > NULL)` reported every present row FALSE. Documented on `mask_not_under`. +- **Parity:** groups `0x13xx` (bijection: both directions, witness, round trip) and `0x139xx–0x13Cxx` (6 ops × 9 thresholds × 12 lengths, f64 and f32, plus the FALSE plane) against a bit-pattern reference (NaN by exponent/mantissa, order by sign-magnitude key, `-0 == +0`); no float compare in the reference. Arms green: native AVX-512, v3 AVX2, wasm, wasm-scalar, wasm-relaxed, neon-qemu, avx-qemu, v2-qemu, nightly. +- **Disable runs (each restored), all red:** f64 packer tail mask dropped (0x1301); `x == x` replaced by `x < inf` (0x1304); `nan_where_clear` skips the partial chunk (0x1322); witness counts dead bits (0x1341); SQL NULL-row gate dropped (0x13915, `<>`); NULL-literal threshold not refused (0x13905); NOT as a bare complement (0x13b00). +- **Not done:** the masked float fold (`n, mean, M2`, fixed order) and bit-identical float output across backends; three-valued AND/OR of two nullable predicates (TRUE/FALSE planes combine as `Ta & Tb` / `Fa | Fb` and dually; no primitive needed beyond mask_and/or, but no test either); GROUP BY over float keys (NaN != NaN, key on presence). + ## 2026-10-27 — wasm `mul_add` certification audit (#378): relaxed-SIMD, proof obligations, X11 - **Relaxed SIMD, tested:** no source path selects on `relaxed-simd` (only the `simd_wasm.rs` module note mentions it). The parity wasm built with `+simd128` and with `+simd128,+relaxed-simd` has a byte-identical code section (621,383 bytes); only the `name` and `target_features` custom sections differ. Byte scan of the code section for relaxed-SIMD opcodes (`0xFD`, sub-opcode `0x100..0x113`): 0 hits; a control calling `f64x2_relaxed_madd` scans to 1 hit. The relaxed build passes under node 22. New arm `scripts/masking-parity.sh wasm-relaxed` and a CI step in `realization/wasm + scalar × wasm32`, as a guard against a relaxed path appearing later. What it establishes: the shipped arithmetic is the same under both flags for this compiler; it does not test an engine executing relaxed instructions, since none are emitted. diff --git a/.claude/plans/nan-exclusion-mindmap-scratch.md b/.claude/plans/nan-exclusion-mindmap-scratch.md index df526783..305a661f 100644 --- a/.claude/plans/nan-exclusion-mindmap-scratch.md +++ b/.claude/plans/nan-exclusion-mindmap-scratch.md @@ -245,3 +245,16 @@ Shipped primitives (`src/simd_masking_ops.rs`, re-exported from `ndarray::simd`) `F64Mask8::to_bitmask` was added on the AVX-512, AVX2/AVX, NEON, wasm and scalar backends (nightly already had it) so the f64 packer reads one spelling on every arm. + +### SQL over a NULL=NaN lane (2026-10-27) + +With NULL=NaN the SQL aggregates need no caller clean-up: `COUNT(*)` = +`popcount(filter)`, `COUNT(x)` = `popcount(filter & presence)`, and +`SUM`/`AVG`/`VAR` fold over the same plane, with `n == 0` giving NULL (NaN). +Comparisons are the `{lt,le,gt,ge,eq,ne}_{f64,f32}_to_mask` TRUE planes. NOT +is where IEEE and SQL differ (`!(NaN > 5)` is true, `NOT (NULL > 5)` is +UNKNOWN), so a lowering negates with `mask_not_under(pred, presence)`, where +`presence` is the predicate's presence: the AND of its operands' planes, empty +for a NULL literal. Known differences from PostgreSQL/DuckDB: a stored NaN +value reads back as NULL; GROUP BY and DISTINCT on float keys must key on +presence, since `NaN != NaN`. From 1a155e482312d943948df209a9c748a7dd2ec0ae Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 15:43:10 +0000 Subject: [PATCH 4/6] =?UTF-8?q?lab:=20sm9=5Fpresence=5Fprobe=20=E2=80=94?= =?UTF-8?q?=20sign-magnitude=209-bit=20lanes=20with=20NULL=20on=20the=20pr?= =?UTF-8?q?esence=20plane?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Follow-up to #380: NaN is a cleared presence bit, not the 'negative zero' codeword, so all of -255..=+255 are values and the writer stores +0 for zero and for missing rows. Masked dot splits lanes by sign_a ^ sign_b into two unsigned u8 x u8 accumulators, NULL excluded by the presence AND, n the popcount of the same mask. Decode applies nan_where_clear, so the float NaN/presence bijection holds on the decoded lane. tanh exactness vs i8 Sym (-127..=127), measured over 2,000,001 points: output quantization alone 2.01x lower max error, 1.94x lower RMS; a 256-entry table (sign-magnitude indexed by magnitude, sign passed through) 1.92-1.99x max, 1.97-2.02x RMS for input ranges 2..8; saturation onset 3.11 -> 3.46. This is the one extra bit. Scalar reference only; SIMD cost not measured. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- Cargo.toml | 4 + examples/sm9_presence_probe.rs | 385 +++++++++++++++++++++++++++++++++ 2 files changed, 389 insertions(+) create mode 100644 examples/sm9_presence_probe.rs diff --git a/Cargo.toml b/Cargo.toml index ae75923a..afdcfdff 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -58,6 +58,10 @@ required-features = ["std"] name = "codec8_asm_probe" required-features = ["std"] +[[example]] +name = "sm9_presence_probe" +required-features = ["std"] + [[example]] name = "cpu_guard_probe" required-features = ["std"] diff --git a/examples/sm9_presence_probe.rs b/examples/sm9_presence_probe.rs new file mode 100644 index 00000000..8305ec0b --- /dev/null +++ b/examples/sm9_presence_probe.rs @@ -0,0 +1,385 @@ +//! Sign-magnitude 9-bit lanes with NULL carried by the presence plane. +//! +//! A follow-up to the Palette512 / SignedMagnitude9 probe (#380). There, the +//! codeword "sign set, magnitude 0" was reserved as the NaN sentinel. Here NaN +//! is not a codeword: a row is NaN (NULL) because its presence bit is clear, +//! whatever its payload holds. That is the NaN/presence rule `ndarray::simd` +//! ships for float lanes (`not_nan_*_to_mask`, `nan_where_clear_*`, +//! `presence_nan_mismatches_*`), applied to an integer lane, where the presence +//! plane is the only carrier. +//! +//! Per 64-lane block: 64 `u8` magnitudes, one `u64` sign plane, one `u64` +//! presence plane. Every magnitude/sign pair is a value, -255..=+255; the +//! writer stores +0 for zero and for missing rows, so each row has one encoding. +//! +//! What it measures: +//! +//! 1. The masked dot product: lanes split by `sign_a ^ sign_b` into a positive +//! and a negative unsigned accumulator, NULL rows excluded by the presence +//! AND, `n` the popcount of the same mask. Checked exhaustively against a +//! decoded oracle. Scalar reference only: the SIMD cost of the u8 x u8 +//! accumulation is not measured here. +//! 2. tanh exactness against i8 symmetric (-127..=127), both as a 256-entry +//! lookup table and as output quantization alone. tanh is odd, so the +//! sign-magnitude table is indexed by magnitude and the sign plane passes +//! through untouched. +//! +//! Run `cargo run --release --example sm9_presence_probe` for the table and +//! `cargo test --release --example sm9_presence_probe` for the checks. + +use ndarray::simd::nan_where_clear_f64; + +/// One 64-lane block. `presence` bit i clear means lane i is NULL (NaN). +#[derive(Clone, Debug)] +struct Sm9Block { + magnitudes: [u8; 64], + signs: u64, + presence: u64, +} + +impl Sm9Block { + /// Store `values[i]` (`None` = NULL). A NULL or zero row is written as +0. + fn from_values(values: &[Option; 64]) -> Self { + let mut b = Self { + magnitudes: [0; 64], + signs: 0, + presence: 0, + }; + for (i, v) in values.iter().enumerate() { + if let Some(v) = *v { + assert!((-255..=255).contains(&v), "value {v} outside -255..=255"); + b.presence |= 1 << i; + b.magnitudes[i] = v.unsigned_abs() as u8; + if v < 0 { + b.signs |= 1 << i; + } + } + } + b + } + + /// The value of lane `i`, `None` when the presence bit is clear. + fn value(&self, i: usize) -> Option { + if (self.presence >> i) & 1 == 0 { + return None; + } + let m = i16::from(self.magnitudes[i]); + Some(if (self.signs >> i) & 1 == 1 { -m } else { m }) + } + + /// Masked dot product over the rows both blocks hold (`filter` narrows it + /// further). Returns `(dot, n)`. The sign is a mask, not arithmetic: lanes + /// whose signs differ go to the negative accumulator. + fn masked_dot(&self, rhs: &Self, filter: u64) -> (i64, u32) { + let sel = self.presence & rhs.presence & filter; + let negative = self.signs ^ rhs.signs; + let (mut pos, mut neg) = (0u64, 0u64); + let mut bits = sel; + while bits != 0 { + let i = bits.trailing_zeros() as usize; + let p = u64::from(self.magnitudes[i]) * u64::from(rhs.magnitudes[i]); + if (negative >> i) & 1 == 1 { + neg += p; + } else { + pos += p; + } + bits &= bits - 1; + } + (pos as i64 - neg as i64, sel.count_ones()) + } + + /// Decode to an f64 lane at `scale`, NaN where the presence bit is clear. + fn decode(&self, scale: f64) -> [f64; 64] { + let mut out = [0.0f64; 64]; + for (i, o) in out.iter_mut().enumerate() { + let m = f64::from(self.magnitudes[i]) * scale; + *o = if (self.signs >> i) & 1 == 1 { -m } else { m }; + } + nan_where_clear_f64(&[self.presence], &mut out); + out + } +} + +// ── tanh exactness ────────────────────────────────────────────────────────── + +/// i8 symmetric code for `x` at `scale`: round half to even, saturate at ±127. +fn q_i8(x: f64, scale: f64) -> i32 { + ((x / scale).round_ties_even() as i32).clamp(-127, 127) +} + +/// Sign-magnitude code for `x` at `scale`: (negative, magnitude 0..=255); +/// zero is always stored as +0. +fn q_sm9(x: f64, scale: f64) -> (bool, u8) { + let m = ((x.abs() / scale).round_ties_even() as i32).clamp(0, 255) as u8; + (x < 0.0 && m > 0, m) +} + +/// Error statistics of an approximation `f` to tanh over a grid on [-r, r]. +#[derive(Debug, Clone, Copy)] +struct Err { + max: f64, + rms: f64, +} + +fn tanh_error(r: f64, points: usize, f: impl Fn(f64) -> f64) -> Err { + let (mut max, mut sq) = (0.0f64, 0.0f64); + for k in 0..points { + let x = -r + 2.0 * r * (k as f64) / ((points - 1) as f64); + let e = (f(x) - x.tanh()).abs(); + max = max.max(e); + sq += e * e; + } + Err { + max, + rms: (sq / points as f64).sqrt(), + } +} + +/// The two lookup tables for inputs on [-r, r]. +struct Tables { + r: f64, + /// i8 table: input code -127..=127 at offset +127, output code. + lut_i8: [i8; 255], + /// Sign-magnitude table: input magnitude -> output magnitude. tanh is odd, + /// so the sign is not looked up. + lut_sm9: [u8; 256], +} + +impl Tables { + fn new(r: f64) -> Self { + let mut lut_i8 = [0i8; 255]; + for (k, o) in lut_i8.iter_mut().enumerate() { + let x = (k as f64 - 127.0) * (r / 127.0); + *o = q_i8(x.tanh(), 1.0 / 127.0) as i8; + } + let mut lut_sm9 = [0u8; 256]; + for (m, o) in lut_sm9.iter_mut().enumerate() { + let x = m as f64 * (r / 255.0); + *o = q_sm9(x.tanh(), 1.0 / 255.0).1; + } + Self { r, lut_i8, lut_sm9 } + } + + fn tanh_i8(&self, x: f64) -> f64 { + let c = q_i8(x, self.r / 127.0); + f64::from(self.lut_i8[(c + 127) as usize]) / 127.0 + } + + fn tanh_sm9(&self, x: f64) -> f64 { + let (neg, m) = q_sm9(x, self.r / 255.0); + let y = f64::from(self.lut_sm9[m as usize]) / 255.0; + if neg { + -y + } else { + y + } + } +} + +/// The largest |x| at which the quantized tanh output is still below the top +/// code, i.e. how far into saturation the format still tells values apart. +fn saturation_onset(top_code: f64) -> f64 { + // tanh(x) rounds to the top code once tanh(x) >= (top - 0.5) / top. + ((top_code - 0.5) / top_code).atanh() +} + +fn main() { + const POINTS: usize = 2_000_001; + println!("tanh over [-r, r], {POINTS} points; errors are |approx - tanh(x)| in absolute units"); + println!(); + println!("output quantization only (exact input):"); + let out_i8 = tanh_error(4.0, POINTS, |x| f64::from(q_i8(x.tanh(), 1.0 / 127.0)) / 127.0); + let out_sm9 = tanh_error(4.0, POINTS, |x| { + let (neg, m) = q_sm9(x.tanh(), 1.0 / 255.0); + let y = f64::from(m) / 255.0; + if neg { + -y + } else { + y + } + }); + println!( + " i8 sym max {:.6} rms {:.6}\n sm9 max {:.6} rms {:.6} (max ratio {:.2}x, rms ratio {:.2}x)", + out_i8.max, + out_i8.rms, + out_sm9.max, + out_sm9.rms, + out_i8.max / out_sm9.max, + out_i8.rms / out_sm9.rms + ); + println!(); + println!("256-entry lookup table (input and output quantized):"); + println!(" r i8 max i8 rms sm9 max sm9 rms max x rms x"); + for r in [2.0, 3.0, 4.0, 6.0, 8.0] { + let t = Tables::new(r); + let a = tanh_error(r, POINTS, |x| t.tanh_i8(x)); + let b = tanh_error(r, POINTS, |x| t.tanh_sm9(x)); + println!( + " {r:>4.1} {:.6} {:.6} {:.6} {:.6} {:>5.2}x {:>5.2}x", + a.max, + a.rms, + b.max, + b.rms, + a.max / b.max, + a.rms / b.rms + ); + } + println!(); + println!( + "saturation onset (largest |x| still below the top output code): i8 {:.3}, sm9 {:.3}", + saturation_onset(127.0), + saturation_onset(255.0) + ); + + // The masked dot on one block: every third row NULL, signs alternating. + let a: [Option; 64] = std::array::from_fn(|i| (i % 3 != 0).then_some(if i % 2 == 0 { 255 } else { -255 })); + let b: [Option; 64] = std::array::from_fn(|_| Some(-255)); + let (ba, bb) = (Sm9Block::from_values(&a), Sm9Block::from_values(&b)); + let (dot, n) = ba.masked_dot(&bb, u64::MAX); + let nulls = ba.decode(1.0 / 255.0).iter().filter(|x| x.is_nan()).count(); + let first = ba.value(1).map_or(0, i16::from); + println!(); + println!("masked dot: dot {dot}, n {n}, NULL rows {nulls} (decode NaN), lane 1 = {first}"); +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::simd::presence_nan_mismatches_f64; + + /// Every (value or NULL) x (value or NULL) pair: the masked product + /// equals the decoded product, and NULL contributes nothing to dot or n. + #[test] + fn exhaustive_masked_product_matches_decoded_oracle() { + let mut checked = 0u64; + let all: Vec> = std::iter::once(None) + .chain((-255..=255).map(Some)) + .collect(); + for &a in &all { + for &b in &all { + let mut va = [None; 64]; + let mut vb = [None; 64]; + va[37] = a; + vb[37] = b; + let (ba, bb) = (Sm9Block::from_values(&va), Sm9Block::from_values(&vb)); + let want = match (a, b) { + (Some(x), Some(y)) => (i64::from(x) * i64::from(y), 1), + _ => (0, 0), + }; + assert_eq!(ba.masked_dot(&bb, u64::MAX), want, "a={a:?} b={b:?}"); + checked += 1; + } + } + // 511 values plus NULL on each side. + assert_eq!(checked, 512 * 512); + } + + /// Zero and NULL rows are stored as +0, so no row has two encodings, and + /// the payload under a clear presence bit never reaches the result. + #[test] + fn null_and_zero_store_canonical_plus_zero_and_payload_is_ignored() { + let mut v = [Some(7i16); 64]; + v[3] = None; + v[4] = Some(0); + let mut b = Sm9Block::from_values(&v); + assert_eq!((b.magnitudes[3], (b.signs >> 3) & 1), (0, 0)); + assert_eq!((b.magnitudes[4], (b.signs >> 4) & 1), (0, 0)); + let ones = Sm9Block::from_values(&[Some(1i16); 64]); + let before = b.masked_dot(&ones, u64::MAX); + // Garbage under the clear bit changes nothing. + b.magnitudes[3] = 255; + b.signs |= 1 << 3; + assert_eq!(b.masked_dot(&ones, u64::MAX), before); + assert_eq!(before, (7 * 62, 63)); + } + + /// Random blocks with NULLs and a filter, against a wide decoded oracle. + #[test] + fn random_blocks_match_wide_oracle_and_n_is_the_selection() { + let mut s = 0x9E37_79B9_7F4A_7C15u64; + let mut next = || { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + s + }; + let mut nulls_skipped = 0u32; + for _ in 0..2000 { + let mut va = [None; 64]; + let mut vb = [None; 64]; + for i in 0..64 { + let r = next(); + va[i] = (r % 5 != 0).then(|| ((r >> 8) % 511) as i16 - 255); + vb[i] = ((r >> 32) % 7 != 0).then(|| ((r >> 40) % 511) as i16 - 255); + } + let filter = next(); + let (a, b) = (Sm9Block::from_values(&va), Sm9Block::from_values(&vb)); + let mut want = (0i64, 0u32); + for i in 0..64 { + if (filter >> i) & 1 == 1 { + match (a.value(i), b.value(i)) { + (Some(x), Some(y)) => { + want.0 += i64::from(x) * i64::from(y); + want.1 += 1; + } + _ => nulls_skipped += 1, + } + } + } + assert_eq!(a.masked_dot(&b, filter), want); + } + assert!(nulls_skipped > 10_000, "only {nulls_skipped} NULL rows skipped"); + } + + /// Decoding gives a float lane on which the NaN/presence bijection holds. + #[test] + fn decode_satisfies_the_float_bijection() { + let mut v = [Some(-200i16); 64]; + for i in (0..64).step_by(3) { + v[i] = None; + } + let b = Sm9Block::from_values(&v); + let lane = b.decode(1.0 / 255.0); + assert_eq!(presence_nan_mismatches_f64(&lane, &[b.presence]), 0); + assert_eq!(lane.iter().filter(|x| x.is_nan()).count(), 22); + } + + /// The sign-magnitude table is odd by construction; check it, and check + /// that it is more exact than the i8 table at every tested range. + #[test] + fn sm9_tanh_table_is_odd_and_more_exact_than_i8() { + for r in [2.0, 4.0, 8.0] { + let t = Tables::new(r); + for k in 0..=4000 { + let x = r * k as f64 / 4000.0; + assert_eq!(t.tanh_sm9(-x), -t.tanh_sm9(x), "r={r} x={x}"); + } + let a = tanh_error(r, 200_001, |x| t.tanh_i8(x)); + let b = tanh_error(r, 200_001, |x| t.tanh_sm9(x)); + assert!(b.max < a.max && b.rms < a.rms, "r={r}: sm9 {b:?} vs i8 {a:?}"); + } + } + + /// The inline i8 quantizer used as the baseline is the shipped codec. + #[test] + fn i8_baseline_matches_codec8_sym() { + use ndarray::hpc::codec8::{Calibration, Codec8, NoNaN, Sym, TiesEven, I8}; + struct Out; + impl Calibration for Out { + const SCALE: f32 = 1.0 / 127.0; + const ZERO_POINT: i32 = 0; + } + let xs: Vec = (0..=2000).map(|k| -1.0 + k as f32 / 1000.0).collect(); + let mut codes = vec![0i8; xs.len()]; + Codec8::::encode_checked::(&xs, &mut codes).unwrap(); + let mut differ = 0; + for (&x, &c) in xs.iter().zip(&codes) { + // The codec divides in f32; allow the inline f64 version to differ + // only at an exact f32/f64 half-tie. + if q_i8(f64::from(x), 1.0 / 127.0) != i32::from(c) { + differ += 1; + } + } + assert!(differ <= 2, "{differ} codes differ from Codec8"); + } +} From 2dfbcedcbf6a41c46895f042622c10bc8fdd3c29 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 15:51:47 +0000 Subject: [PATCH 5/6] probe: 3-sigma coverage gate for SM9 vs i8 tanh table At a tolerance of 0.010 fixed in advance (between the two formats' measured 99.73% error quantiles at r = 4), sign-magnitude covers at least 0.9973 of 200,000 independent uniform inputs and i8 does not; the paired mean error difference is far above 3 standard errors. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- examples/sm9_presence_probe.rs | 43 ++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/examples/sm9_presence_probe.rs b/examples/sm9_presence_probe.rs index 8305ec0b..16ae9605 100644 --- a/examples/sm9_presence_probe.rs +++ b/examples/sm9_presence_probe.rs @@ -344,6 +344,49 @@ mod tests { assert_eq!(lane.iter().filter(|x| x.is_nan()).count(), 22); } + /// The 3σ gate: at a tolerance fixed before measuring, the fraction of + /// inputs whose table error stays within it must reach 0.9973 (the + /// two-sided 3σ coverage). On 200,000 independent uniform inputs over + /// [-4, 4], sign-magnitude passes and i8 fails, and the paired mean + /// difference of the errors is far above 3 standard errors. + /// + /// The tolerance 0.010 sits between the two formats' measured 99.73% + /// error quantiles at r = 4 (about 0.0079 and 0.0157), so neither half of + /// the gate is true by construction. + #[test] + fn three_sigma_coverage_gate_passes_for_sm9_and_fails_for_i8() { + const TOL: f64 = 0.010; + const N: usize = 200_000; + const COVERAGE_3SIGMA: f64 = 0.9973; + let t = Tables::new(4.0); + let mut s = 0xD1B5_4A32_D192_ED03u64; + let (mut in8, mut in9) = (0usize, 0usize); + let (mut sum, mut sum_sq) = (0.0f64, 0.0f64); + for _ in 0..N { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + // 53 random bits to a uniform x in [-4, 4). + let x = -4.0 + 8.0 * ((s >> 11) as f64 / (1u64 << 53) as f64); + let e8 = (t.tanh_i8(x) - x.tanh()).abs(); + let e9 = (t.tanh_sm9(x) - x.tanh()).abs(); + in8 += usize::from(e8 <= TOL); + in9 += usize::from(e9 <= TOL); + let d = e8 - e9; + sum += d; + sum_sq += d * d; + } + let (cov8, cov9) = (in8 as f64 / N as f64, in9 as f64 / N as f64); + let mean = sum / N as f64; + let var = (sum_sq - N as f64 * mean * mean) / (N - 1) as f64; + let z = mean / (var / N as f64).sqrt(); + assert!(cov9 >= COVERAGE_3SIGMA, "sm9 coverage {cov9:.4} below 3σ at tol {TOL}"); + assert!(cov8 < COVERAGE_3SIGMA, "i8 coverage {cov8:.4} unexpectedly reaches 3σ at tol {TOL}"); + assert!(z > 3.0, "paired z {z:.1} not significant at 3σ"); + // Anti-vacuity: i8 is not merely a hair below the line. + assert!(cov8 < 0.99, "i8 coverage {cov8:.4}: the gate barely separates the formats"); + } + /// The sign-magnitude table is odd by construction; check it, and check /// that it is more exact than the i8 table at every tested range. #[test] From 7be0676737a8de958c78204e7cbe2eaf2057206b Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 15:52:10 +0000 Subject: [PATCH 6/6] probe: print the 3-sigma gate's measured coverage and z Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- examples/sm9_presence_probe.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/sm9_presence_probe.rs b/examples/sm9_presence_probe.rs index 16ae9605..895fb475 100644 --- a/examples/sm9_presence_probe.rs +++ b/examples/sm9_presence_probe.rs @@ -380,6 +380,7 @@ mod tests { let mean = sum / N as f64; let var = (sum_sq - N as f64 * mean * mean) / (N - 1) as f64; let z = mean / (var / N as f64).sqrt(); + eprintln!("3σ gate: tol {TOL}, n {N}: sm9 {cov9:.4}, i8 {cov8:.4}, paired z {z:.1}"); assert!(cov9 >= COVERAGE_3SIGMA, "sm9 coverage {cov9:.4} below 3σ at tol {TOL}"); assert!(cov8 < COVERAGE_3SIGMA, "i8 coverage {cov8:.4} unexpectedly reaches 3σ at tol {TOL}"); assert!(z > 3.0, "paired z {z:.1} not significant at 3σ");