Add I32x16::dot4_u8_i8_add_masked for exact 9-bit sign-magnitude dot products

.claude/blackboard.md

@@ -3876,3 +3876,25 @@ Loose ends: AVX-VNNI (`avxvnni`, 256-bit `vpdpbusd`) is not used on the AVX2 arm
 is compiled in. Neither is measured. The signed form (`i8 x i8`, turbovec's `xor 0x80` + `zero` correction) is not
 a facade word yet.
 
+
+## dot4_u8_i8_add_masked + 9-bit TurboVec9 cost probe (2026-10-29, branch claude/u9-turbovec9-substrate)
+
+Follow-up to #386. Added on every backend: `I32x16::dot4_u8_i8_add_masked(a, keep, b)`; two calls with the positive
+plane and its complement give an exact 9-bit sign-magnitude x i8 dot product (all i8 incl. -128, i32 accumulators,
+no widening). Crate-private `simd_masking_ops::{BIT_IN_BYTE, mask_byte_spread}` expand a u64 mask to bytes with one
+AND + one compare per register. Measurements, precision tables and the TurboVec handoff:
+`.claude/knowledge/turbovec9-u9-substrate.md`; probe `examples/u9_cost_probe.rs` (R1-excluded: it times shuffle
+lookups on purpose).
+
+Witnessed: 8/8 dot4 tests on native (VNNI), v4, v3, avx (qemu), v2 (qemu), aarch64 (qemu), nightly; neon/wasm parity
+programs pass (0x504, 0x505). clippy -D warnings clean on native, v4, v3. Disable runs all red: unmasked AVX-512 move,
+reversed BIT_IN_BYTE (AVX2), mask dropped in the fused AVX2 loop.
+
+Key numbers (AVX-512 VNNI, L1): masked pair 1.27-1.33 ns/64 lanes vs dot4 1.10-1.19; 512-entry lookup 6.1 ns
+(vpermi2b) / 38-40 ns (per-position reads). Stored uniform SM9 46.0 dB vs i4 20.2 dB at the same per-dim L1 cost and
+1.4x per dim streamed. AVX2: masked pair 6.35-6.44 ns vs dot4 1.59-1.79 ns (corrected after review: the first
+facade timing returned lane 0 only and let the optimizer drop half the AVX2 work; probe kernel B also mishandled a
+negative zero; both fixed, every exact kernel now checked against an i64 reference before it is timed).
+
+Loose ends: fused single-widen SM9 op for non-VNNI tiers; NEON SDOT/USDOT and AVX-VNNI unused; no NEON/wasm timings;
+recall on a real TurboVec index is the TurboVec owner's measurement.

.claude/hooks/facet-stride-exclude.txt

@@ -4,3 +4,4 @@
 src/hpc/perm.rs | R1 | Perm64 is a 64-lane permutation map; composing and materializing it is a genuine byte move, the operation itself (module doc: "composable maps").
 examples/ternlog_popcnt_gap_probe.rs | R1 | Benchmark that measures the nibble-LUT popcount as its baseline against ternlog + popcnt; the shuffle is the thing being measured.
 vendor/chacha20/* | R1 | ChaCha20 rotate-left by 8 and 16 bits is a fixed byte rotation (pshufb / vqtbl1q with a constant mask): the cipher's defined operation, in vendored upstream code.
+examples/u9_cost_probe.rs | R1 | Cost probe for 9-bit TurboVec scoring: the vpshufb / vpermb / vpermi2b table lookups are the candidates being timed against the lookup-free layout.

.claude/knowledge/turbovec9-u9-substrate.md

@@ -0,0 +1,159 @@
+# TurboVec9 on ndarray: 9-bit sign-magnitude scoring, measured
+
+READ BY: savant-architect, vector-synthesis, simd-savant, anyone wiring a
+9-bit index (TurboVec9) or a DR512 dot product onto `ndarray::simd`.
+
+Status: probe + one facade word (`I32x16::dot4_u8_i8_add_masked`). Numbers
+below come from `examples/u9_cost_probe.rs` on one AVX-512 VNNI + VBMI Xeon
+(target-cpu=native), three runs, reported as ranges. Every exact SM9 kernel
+must match an independent i64 reference before its time is printed, and
+every facade kernel folds all 16 accumulator lanes into its result. They are not TurboVec
+recall numbers; the codebooks in the probe are probe-local Lloyd-Max fits.
+
+## Two contracts, never mixed
+
+- **DR512 / Codec9** (`hpc::codec9`, `simd_stride::lut9`): u8 magnitude plus
+  one P bit per value, -255..=255, and `P=0, M=0` is the never-touched
+  sentinel. `lut9` belongs here and nowhere else.
+- **TurboVec9**: an ID 0..=511 selects one of 512 centroids. Every ID is
+  valid. An ID is not a Codec9 value and is never read as sign-magnitude.
+
+The facade word below is pure arithmetic over a u8 magnitude and a u64
+plane. It does not know which contract produced them. For DR512 data both
+`P=0, M=0` and `P=1, M=0` contribute 0 to a dot product; a caller that must
+tell them apart reads the plane itself.
+
+## The facade word
+
+```rust
+// Lane j gains the sum of a[i] * b[i] over bytes i in 4j..4j+4 with bit i of keep set.
+I32x16::dot4_u8_i8_add_masked(self, a: U8x64, keep: u64, b: I8x64) -> I32x16
+```
+
+9-bit sign-magnitude times i8, exact for every magnitude and every i8
+(-128 included), i32 accumulators, operands never widened:
+
+```rust
+let (mut pos, mut neg) = (I32x16::splat(0), I32x16::splat(0));
+for blk in blocks {                       // 64 dims per block
+    pos = pos.dot4_u8_i8_add_masked(blk.magnitude, blk.positive, w);
+    neg = neg.dot4_u8_i8_add_masked(blk.magnitude, !blk.positive, w);
+}
+let score_lanes = pos - neg;              // 16 lanes, each a 4-dim partial sum
+```
+
+Overflow: each accumulator moves by at most `4 * 255 * 128 = 130_560` per
+call, so 16_448 calls (65_792 dims per lane) cannot wrap; the test
+`dot4_masked_worst_case_holds_16448_calls` pins both sides of that bound.
+No saturating instruction is involved (no `vpmaddubsw`, no `vpdpbusds`).
+
+Realizations: AVX-512 one `vmovdqu8{z}` before `vpdpbusd` (or the widen +
+`vpmaddwd` path without VNNI); AVX2 builds the byte mask in registers
+(variable shift, multiply by `0x0101_0101`, shift by 32) and applies it
+inside the dot loop; AVX, NEON and wasm expand a stored spread with one AND
+against `BIT_IN_BYTE` and one byte compare per register (no shuffle, no
+per-lane loop); nightly selects through `core::simd::Mask`;
+scalar is the reference loop.
+
+Witnessed: `simd::dot4_nibble_facade_parity` 8/8 on native (AVX-512 VNNI),
+config-v4 (AVX-512 without VNNI), config-v3 (AVX2), config-avx (qemu
+SandyBridge), config-v2 (scalar, qemu Nehalem), aarch64 (qemu) and
+nightly-simd; `neon-parity.sh` and `wasm-parity.sh` pass (codes 0x504,
+0x505). Disable runs, all red: unmasked AVX-512 move (3 tests), reversed
+`BIT_IN_BYTE` on AVX2 (2 tests), mask dropped in the fused AVX2 loop (3),
+wrong byte shifts for the upper AVX2 half (2).
+
+## Arithmetic: three exact forms compared (AVX-512 VNNI, L1, ns per 64 lanes)
+
+| kernel | ns / 64 lanes |
+|---|---|
+| u8 x i8 `vpdpbusd` (#386 `dot4_u8_i8_add`) | 1.12 - 1.14 |
+| i8 x i8 via `xor 0x80` bias | 1.20 - 1.25 |
+| A: sign folded into w (exact only for w in -127..=127) | 1.23 - 1.26 |
+| B: two's-complement limbs, `dot(lo) - 256 * dot(sign)`, sign limb skips magnitude 0 | 1.89 - 1.92 |
+| C: `dot(m) - 2 * dot(m & neg)`, one accumulator | 1.79 - 1.88 |
+| **C': `dot(m & pos) - dot(m & neg)`, two accumulators** | **1.28 - 1.30** |
+
+The loop is one latency-bound chain. B and C put two dependent ops on the
+accumulator; C' keeps two independent chains and lands within ~10% of the
+plain u8 dot. A is not exact for w = -128 and was not added.
+
+Through the facade, same host: `dot4_u8_i8_add` 1.08 - 1.15, two masked
+accumulators 1.34 - 1.38. On config-v3 (AVX2, no VNNI): `dot4_u8_i8_add`
+1.59 - 1.79, two plain accumulators 2.75 - 2.91, two masked accumulators
+6.35 - 6.44. Without VNNI each masked call adds the mask expansion (about
+1.4 ns, measured as one masked accumulator against one plain one) on top of
+its own 16-bit widening, so the pair costs about 4x one dot4. The loop is
+instruction-bound, not memory-bound: building the mask in registers instead
+of from a stored spread did not change the time. A fused SM9 op that widens
+once and negates in 16 bits would avoid both the second widening and the
+mask, and is not built.
+
+## 512-entry lookup (WP2): keep it out of the scan
+
+ID = low byte + high-bit plane, decoded to magnitude + positive plane:
+
+| lookup | ns / 64 lanes |
+|---|---|
+| per-position reads (`lut_u8` shape, two 512-entry tables) | 38 - 39 |
+| 4x `vpermi2b` + blends for magnitude, same for the sign | 6.03 - 6.09 |
+| 2x `vpermi2b` on an antisymmetric codebook (`255 - lo` for negatives) | 6.72 - 6.77 |
+
+`lut9` and `lut_u8` keep their contracts; no 512-entry facade word was
+added. The lookup alone costs ~5x the arithmetic it feeds. TurboVec's
+levels do not depend on the query, so the decode can happen once at encode
+time; then the scan stores values, not IDs (next section).
+
+## Precision vs cost: the layouts (N(0,1) coordinates)
+
+| layout | SNR | scan cost (AVX-512 VNNI) |
+|---|---|---|
+| i4, 16 Lloyd-Max levels (today) | 20.22 dB | 2.65 - 2.68 ns / 128 dims (L1) |
+| i8 uniform grid, best clip | 40.54 dB | 1.12 - 1.14 ns / 64 dims |
+| 512 Lloyd-Max IDs, executed by `vpermi2b` lookup | 49.86 dB | 5.37 - 5.70 ns / 64 dims |
+| 512 Lloyd-Max levels snapped to the SM9 grid | 44.74 dB (329 distinct) | 1.29 - 1.32 ns / 64 dims |
+| **SM9 uniform grid, best clip (4.21 sigma), stored** | **46.02 dB** | **1.29 - 1.32 ns / 64 dims** |
+| 4+5 additive C16 + R32, two i8 grids | 43.25 dB | 2.61 - 2.65 ns / 64 dims |
+
+Streamed from DRAM (64 MiB of codes, weights in L1), per dimension:
+i4 0.035 - 0.037 ns, stored SM9 0.048 - 0.052 ns (72 bytes per 64 dims).
+
+Reading the table:
+
+- Stored values make the 512 Lloyd-Max IDs pointless: snapping them to
+  the 9-bit grid keeps 329 distinct values and loses to a plain uniform
+  9-bit grid with a tuned clip. The Lloyd-Max levels only pay off when
+  executed as IDs, and the lookup makes that path ~4x slower than the
+  stored one.
+- The 4+5 additive competitor is worse than stored SM9 on both axes in this
+  probe. Its residual codebook is one shared Lloyd fit after a 16-level
+  split; a joint fit could do better, but it cannot beat the 512-point
+  optimum and it pays two lookups.
+
+## Answer to the acceptance question
+
+On AVX-512 VNNI: yes, for the stored-value layout. 64 dims cost one
+masked-pair of `vpdpbusd`, ~1.3 ns in L1 (about the same per dimension as
+today's i4 permute-dot) and 1.4x i4 per dimension when streamed from DRAM,
+for +25.8 dB over i4. With IDs plus an in-scan 512 lookup: no, the lookup
+dominates. On AVX2 the pair costs ~4x one plain dot4; a fused SM9 op is the
+next step there if AVX2 matters.
+
+## Handoff to the TurboVec owner (owns ABI, encode, pack, search, mutation, serialization, Python)
+
+Proposed scan layout, per 64 dimensions of one vector, vector-major as
+today: 64 magnitude bytes plus one u64 positive plane (72 bytes). Scoring:
+the two-accumulator loop above, then the existing `scale` / `zero`
+correction (`zero` is 0 here: magnitudes are unbiased u8, weights stay i8).
+Open decisions on the TurboVec side: grid clip per dimension or global, how
+the encoder writes 0 (either plane bit scores 0), and recall measured on
+real data with TurboVec's own codebook solver. None of this was changed
+from this session.
+
+## Not done
+
+- No NEON or wasm timings (qemu timings mean nothing); NEON `SDOT` /
+  `USDOT` and AVX-VNNI (256-bit) are not used by the masked word.
+- No fused single-widen SM9 op for non-VNNI tiers.
+- No signed i8 x i8 facade word.
+- Recall on a real index.

Cargo.toml

@@ -62,6 +62,10 @@ required-features = ["std"]
 name = "sm9_presence_probe"
 required-features = ["std"]
 
+[[example]]
+name = "u9_cost_probe"
+required-features = ["std"]
+
 [[example]]
 name = "cpu_guard_probe"
 required-features = ["std"]

crates/neon-simd-parity/src/main.rs

@@ -18,7 +18,7 @@ fn main() {
             eprintln!("neon-simd-parity FAILED: lane/op code = {rc}");
             std::process::exit(rc as i32);
         }
-        println!("neon-simd-parity OK: U32x16 / F32x16 / I8x16 lanes, ternlog (256 tables, U32x16 + U64x8), U64x8 algebra, I32x16 compares, dot4 + nibbles — bit-identical to scalar");
+        println!("neon-simd-parity OK: U32x16 / F32x16 / I8x16 lanes, ternlog (256 tables, U32x16 + U64x8), U64x8 algebra, I32x16 compares, dot4 (+ masked) + nibbles — bit-identical to scalar");
     }
     #[cfg(not(target_arch = "aarch64"))]
     {
@@ -530,6 +530,39 @@ mod checks {
                 return Err(0x503);
             }
         }
+        // dot4_u8_i8_add_masked: each keep bit gates its own byte, and the
+        // pos - neg form is exact for every SM9 magnitude in both polarities.
+        for i in 0..64 {
+            let got = I32x16::splat(0)
+                .dot4_u8_i8_add_masked(U8x64::splat(7), !(1u64 << i), I8x64::splat(-3))
+                .to_array();
+            if got != core::array::from_fn(|j| if j == i / 4 { -63 } else { -84 }) {
+                return Err(0x504);
+            }
+        }
+        for w in [-128i8, 127, -1] {
+            for chunk in 0..4usize {
+                let mags: [u8; 64] = core::array::from_fn(|i| (64 * chunk + i) as u8);
+                let plane = 0x5A5A_C3C3_0F0F_9696u64.rotate_left(chunk as u32);
+                let (m, wv) = (U8x64::from_array(mags), I8x64::splat(w));
+                let pos = I32x16::splat(0)
+                    .dot4_u8_i8_add_masked(m, plane, wv)
+                    .to_array();
+                let neg = I32x16::splat(0)
+                    .dot4_u8_i8_add_masked(m, !plane, wv)
+                    .to_array();
+                for j in 0..16 {
+                    let mut want = 0i32;
+                    for i in 4 * j..4 * j + 4 {
+                        let v = i32::from(mags[i]);
+                        want += if plane >> i & 1 == 1 { v } else { -v } * i32::from(w);
+                    }
+                    if pos[j] - neg[j] != want {
+                        return Err(0x505);
+                    }
+                }
+            }
+        }
         Ok(())
     }
 }

crates/wasm-simd-parity/src/lib.rs

@@ -596,5 +596,38 @@ fn check_dot4_and_nibbles() -> Result<(), u32> {
             return Err(0x503);
         }
     }
+    // dot4_u8_i8_add_masked: each keep bit gates its own byte, and the
+    // pos - neg form is exact for every SM9 magnitude in both polarities.
+    for i in 0..64 {
+        let got = I32x16::splat(0)
+            .dot4_u8_i8_add_masked(U8x64::splat(7), !(1u64 << i), I8x64::splat(-3))
+            .to_array();
+        if got != core::array::from_fn(|j| if j == i / 4 { -63 } else { -84 }) {
+            return Err(0x504);
+        }
+    }
+    for w in [-128i8, 127, -1] {
+        for chunk in 0..4usize {
+            let mags: [u8; 64] = core::array::from_fn(|i| (64 * chunk + i) as u8);
+            let plane = 0x5A5A_C3C3_0F0F_9696u64.rotate_left(chunk as u32);
+            let (m, wv) = (U8x64::from_array(mags), I8x64::splat(w));
+            let pos = I32x16::splat(0)
+                .dot4_u8_i8_add_masked(m, plane, wv)
+                .to_array();
+            let neg = I32x16::splat(0)
+                .dot4_u8_i8_add_masked(m, !plane, wv)
+                .to_array();
+            for j in 0..16 {
+                let mut want = 0i32;
+                for i in 4 * j..4 * j + 4 {
+                    let v = i32::from(mags[i]);
+                    want += if plane >> i & 1 == 1 { v } else { -v } * i32::from(w);
+                }
+                if pos[j] - neg[j] != want {
+                    return Err(0x505);
+                }
+            }
+        }
+    }
     Ok(())
 }