Add I32x16::dot4_u8_i8_add and the U8x64 nibble split on every backend

.claude/blackboard.md

@@ -3849,3 +3849,30 @@ clippy errors are all pre-existing (same set on master). Disable run: AVX2 cmpgt
 Trap recorded: `cargo --config <file> clippy` DROPS the config (clippy is an external subcommand); the v3 lint pass
 reported clean in 0.05s while linting the native build. Put `--config` after the subcommand:
 `cargo clippy --config .cargo/config-v3.toml`. `cargo --config <file> test` does apply it.
+
+## dot4_u8_i8_add + U8x64 nibble split on every backend (2026-10-29, branch claude/dot4-nibble-facade)
+
+Trigger: turbovec's VNNI "permute-dot" kernels are raw `vpdpbusd` on x86 and `SDOT` on aarch64 (the same
+operation written twice), plus a GFNI-only macro variant for the high nibble. Adding the operation to the facade
+lets those kernels move onto `ndarray::simd`.
+
+Added on AVX-512, AVX2, AVX, scalar, NEON, wasm and nightly: `I32x16::dot4_u8_i8_add(a: U8x64, b: I8x64)`
+(`vpdpbusd` when avx512vnni is compiled in, otherwise widen-then-`madd`; see the contract doc's P0 section on why
+`vpmaddubsw` is wrong) and `U8x64::low_nibbles` / `high_nibbles` (GFNI `gf2p8affine` when compiled in).
+
+Fixed on the way: on aarch64 and wasm, `U16x32` comes from the scalar file but `F32x16` is native, so the scalar
+`to_f32x16_lo/hi` (PR #385) returned a type that is not `ndarray::simd::F32x16`, and the aarch64 lib tests did not
+compile at all. They now return `crate::simd::F32x16`. CI never builds aarch64/wasm lib tests, only the two
+parity programs, which is why it slipped through; both programs now check dot4 and the nibbles (codes 0x500..).
+
+Witnessed: `simd::dot4_nibble_facade_parity` 4/4 on native (AVX-512 + VNNI + GFNI), config-v4 (AVX-512 emulation),
+config-v3 (AVX2), config-avx (qemu SandyBridge), config-v2 (scalar, qemu Nehalem), aarch64 (qemu, plus
+`u16x32_facade_parity`) and nightly-simd; `neon-parity.sh` and `wasm-parity.sh` pass. Lib suite 2657 passed
+(native). clippy -D warnings clean on native/v4/v3/v2 (`--config` AFTER the subcommand) and aarch64; the AVX arm's
+8 errors are pre-existing, none on added lines. Disable runs, all red as expected: `vpmaddubsw` AVX2 realization
+(saturation test + reference), wrong GFNI matrix (nibble test), `vpdpbusds` for `vpdpbusd` (reference, i32::MAX acc).
+
+Loose ends: AVX-VNNI (`avxvnni`, 256-bit `vpdpbusd`) is not used on the AVX2 arm; NEON could use `USDOT` when i8mm
+is compiled in. Neither is measured. The signed form (`i8 x i8`, turbovec's `xor 0x80` + `zero` correction) is not
+a facade word yet.
+

.claude/knowledge/vertical-simd-consumer-contract.md

@@ -80,6 +80,37 @@ Any saturating-abs primitive in ndarray that does NOT produce `i8::MAX` for `i8:
 
 ---
 
+## Four-byte dot product: `vpmaddubsw` is NOT an emulation of `vpdpbusd` (P0, 2026-10-29)
+
+`I32x16::dot4_u8_i8_add(acc, a: U8x64, b: I8x64)` adds `a[4j..4j+4] · b[4j..4j+4]`
+(u8 times signed i8) to lane `j`. AVX-512 VNNI does it in one `vpdpbusd`; every
+other realization must give the same answer, and the obvious AVX2 emulation does not:
+
+- `vpmaddubsw` sums each pair of `u8 * i8` products into an **i16 with signed
+  saturation**. Two products reach `2 * 255 * 127 = 64_770`, past `i16::MAX`, so the
+  pair clamps to 32_767 and the dot product is wrong. `vpmaddwd` with ones afterwards
+  cannot recover it.
+- The correct realization widens first: even and odd bytes into separate 16-bit lanes
+  (`a` zero-extended with `& 0x00FF` / logical `>> 8`, `b` sign-extended with
+  `<< 8` then arithmetic `>> 8` / arithmetic `>> 8`), then `vpmaddwd` on each, then
+  add. The even lanes give bytes 4j and 4j+2, the odd lanes 4j+1 and 4j+3; every
+  product and partial sum is exact in i32. The same shape is used on AVX-512 without
+  VNNI, AVX2, AVX (SSE2 quarters) and wasm (`i32x4_dot_i16x8`); NEON widens with
+  `vmovl` and sums with `vmull_s16` + `vpaddq_s32`; nightly takes byte k of each
+  32-bit word by shifts.
+- Only the final addition to the accumulator wraps, as `vpdpbusd` does. `vpdpbusds`
+  (the saturating form) is a different operation.
+
+**Mandatory tests:** `simd::dot4_nibble_facade_parity::dot4_does_not_saturate_pair_sums`
+(all-255 times all-127 must give `4 * 255 * 127` per lane, times all -128 `-4 * 255 * 128`) and the
+`i32::MAX` accumulator cases in `dot4_matches_reference`. Disable runs on 2026-10-29:
+a `vpmaddubsw` AVX2 realization fails both; `vpdpbusds` in place of `vpdpbusd` fails
+the second.
+
+`U8x64::low_nibbles` / `high_nibbles` are the nibble split every 4-bit LUT kernel
+starts with. `high_nibbles` is one `gf2p8affine` (matrix `0x1020_4080_0000_0000`) when
+GFNI is compiled in, otherwise a 16-bit shift by 4 and a mask.
+
 ## W1a queue — 5 primitives ndarray must ship
 
 Each is a tight-scope PR. Recommended: one branch per primitive, parallel review.

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 — bit-identical to scalar");
+        println!("neon-simd-parity OK: U32x16 / F32x16 / I8x16 lanes, ternlog (256 tables, U32x16 + U64x8), U64x8 algebra, I32x16 compares, dot4 + nibbles — bit-identical to scalar");
     }
     #[cfg(not(target_arch = "aarch64"))]
     {
@@ -50,6 +50,9 @@ mod checks {
         if let Err(code) = check_i32x16_compare() {
             return code;
         }
+        if let Err(code) = check_dot4_and_nibbles() {
+            return code;
+        }
         0
     }
 
@@ -460,4 +463,73 @@ mod checks {
         }
         Ok(())
     }
+
+    /// `I32x16::dot4_u8_i8_add` and the `U8x64` nibble split against inline
+    /// scalar arithmetic. Bytes cover `0`/`255` and `-128`/`127` (a pair of
+    /// `255 * 127` products overflows `i16`, so a saturating emulation fails
+    /// here), accumulators cover `i32::MIN`/`i32::MAX` (only the final addition
+    /// may wrap), and a single nonzero byte must move exactly its own lane.
+    /// Return codes `0x500 + op`.
+    fn check_dot4_and_nibbles() -> Result<(), u32> {
+        use ndarray::simd::{I32x16, I8x64, U8x64};
+        let mut s = 0x9E37_79B9_7F4A_7C15u64;
+        let mut next = move || {
+            s = s
+                .wrapping_mul(6_364_136_223_846_793_005)
+                .wrapping_add(1_442_695_040_888_963_407);
+            (s >> 33) as u32
+        };
+        let edges_u = [0u8, 255, 128, 127];
+        let edges_i = [-128i8, 127, -1, 0];
+        for round in 0..40 {
+            let a: [u8; 64] = core::array::from_fn(|i| {
+                if round < 4 {
+                    edges_u[(i + round) % 4]
+                } else {
+                    next() as u8
+                }
+            });
+            let b: [i8; 64] = core::array::from_fn(|i| {
+                if round < 4 {
+                    edges_i[(i + round) % 4]
+                } else {
+                    i8::from_ne_bytes([next() as u8])
+                }
+            });
+            let rand: [i32; 16] = core::array::from_fn(|_| next() as i32);
+            for acc in [[0i32; 16], [i32::MAX; 16], [i32::MIN; 16], rand] {
+                let got = I32x16::from_array(acc)
+                    .dot4_u8_i8_add(U8x64::from_array(a), I8x64::from_array(b))
+                    .to_array();
+                let want: [i32; 16] = core::array::from_fn(|j| {
+                    let mut sum = 0i32;
+                    for i in 4 * j..4 * j + 4 {
+                        sum += i32::from(a[i]) * i32::from(b[i]);
+                    }
+                    acc[j].wrapping_add(sum)
+                });
+                if got != want {
+                    return Err(0x500);
+                }
+            }
+            let x = U8x64::from_array(a);
+            if x.low_nibbles().to_array() != a.map(|v| v & 0x0F) {
+                return Err(0x501);
+            }
+            if x.high_nibbles().to_array() != a.map(|v| v >> 4) {
+                return Err(0x502);
+            }
+        }
+        for i in 0..64 {
+            let mut a = [0u8; 64];
+            a[i] = 3;
+            let got = I32x16::splat(0)
+                .dot4_u8_i8_add(U8x64::from_array(a), I8x64::splat(-5))
+                .to_array();
+            if got != core::array::from_fn(|j| if j == i / 4 { -15 } else { 0 }) {
+                return Err(0x503);
+            }
+        }
+        Ok(())
+    }
 }

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

@@ -35,6 +35,9 @@ pub extern "C" fn selfcheck() -> u32 {
     if let Err(code) = check_i32x16_compare() {
         return code;
     }
+    if let Err(code) = check_dot4_and_nibbles() {
+        return code;
+    }
     if let Err(code) = check_exp_and_round() {
         return code;
     }
@@ -526,3 +529,72 @@ fn check_i32x16_compare() -> Result<(), u32> {
     }
     Ok(())
 }
+
+/// `I32x16::dot4_u8_i8_add` and the `U8x64` nibble split against inline
+/// scalar arithmetic. Bytes cover `0`/`255` and `-128`/`127` (a pair of
+/// `255 * 127` products overflows `i16`, so a saturating emulation fails
+/// here), accumulators cover `i32::MIN`/`i32::MAX` (only the final addition
+/// may wrap), and a single nonzero byte must move exactly its own lane.
+/// Return codes `0x500 + op`.
+fn check_dot4_and_nibbles() -> Result<(), u32> {
+    use ndarray::simd::{I32x16, I8x64, U8x64};
+    let mut s = 0x9E37_79B9_7F4A_7C15u64;
+    let mut next = move || {
+        s = s
+            .wrapping_mul(6_364_136_223_846_793_005)
+            .wrapping_add(1_442_695_040_888_963_407);
+        (s >> 33) as u32
+    };
+    let edges_u = [0u8, 255, 128, 127];
+    let edges_i = [-128i8, 127, -1, 0];
+    for round in 0..40 {
+        let a: [u8; 64] = core::array::from_fn(|i| {
+            if round < 4 {
+                edges_u[(i + round) % 4]
+            } else {
+                next() as u8
+            }
+        });
+        let b: [i8; 64] = core::array::from_fn(|i| {
+            if round < 4 {
+                edges_i[(i + round) % 4]
+            } else {
+                i8::from_ne_bytes([next() as u8])
+            }
+        });
+        let rand: [i32; 16] = core::array::from_fn(|_| next() as i32);
+        for acc in [[0i32; 16], [i32::MAX; 16], [i32::MIN; 16], rand] {
+            let got = I32x16::from_array(acc)
+                .dot4_u8_i8_add(U8x64::from_array(a), I8x64::from_array(b))
+                .to_array();
+            let want: [i32; 16] = core::array::from_fn(|j| {
+                let mut sum = 0i32;
+                for i in 4 * j..4 * j + 4 {
+                    sum += i32::from(a[i]) * i32::from(b[i]);
+                }
+                acc[j].wrapping_add(sum)
+            });
+            if got != want {
+                return Err(0x500);
+            }
+        }
+        let x = U8x64::from_array(a);
+        if x.low_nibbles().to_array() != a.map(|v| v & 0x0F) {
+            return Err(0x501);
+        }
+        if x.high_nibbles().to_array() != a.map(|v| v >> 4) {
+            return Err(0x502);
+        }
+    }
+    for i in 0..64 {
+        let mut a = [0u8; 64];
+        a[i] = 3;
+        let got = I32x16::splat(0)
+            .dot4_u8_i8_add(U8x64::from_array(a), I8x64::splat(-5))
+            .to_array();
+        if got != core::array::from_fn(|j| if j == i / 4 { -15 } else { 0 }) {
+            return Err(0x503);
+        }
+    }
+    Ok(())
+}

src/simd.rs

@@ -2313,3 +2313,114 @@ mod u16x32_facade_parity {
         assert_eq!(U16x32::default().to_array(), [0u16; 32]);
     }
 }
+
+/// Parity of `I32x16::dot4_u8_i8_add` and the `U8x64` nibble split against a
+/// plain-array reference, on whichever backend this build selected. The
+/// corpus mixes a fixed-seed generator with the byte extremes (`0`, `255`,
+/// `-128`, `127`) and accumulators at `i32::MIN` / `i32::MAX`, where only the
+/// final addition may wrap.
+#[cfg(test)]
+mod dot4_nibble_facade_parity {
+    use super::{I32x16, I8x64, U8x64};
+
+    fn next(state: &mut u64) -> u64 {
+        *state = state
+            .wrapping_mul(6_364_136_223_846_793_005)
+            .wrapping_add(1_442_695_040_888_963_407);
+        *state >> 33
+    }
+
+    fn bytes() -> Vec<[u8; 64]> {
+        let mut s = 0x9E37_79B9_7F4A_7C15u64;
+        let mut out: Vec<[u8; 64]> = (0..48)
+            .map(|_| core::array::from_fn(|_| next(&mut s) as u8))
+            .collect();
+        for e in [0u8, 1, 0x0F, 0x10, 0x7F, 0x80, 0xF0, 0xFF] {
+            out.push([e; 64]);
+        }
+        out.push(core::array::from_fn(|i| [0u8, 255, 128, 127][i % 4]));
+        out
+    }
+
+    fn signed() -> Vec<[i8; 64]> {
+        let mut s = 0xD1B5_4A32_D192_ED03u64;
+        let mut out: Vec<[i8; 64]> = (0..48)
+            .map(|_| core::array::from_fn(|_| i8::from_ne_bytes([next(&mut s) as u8])))
+            .collect();
+        for e in [0i8, 1, -1, 127, -128] {
+            out.push([e; 64]);
+        }
+        out.push(core::array::from_fn(|i| [-128i8, 127, -1, 0][i % 4]));
+        out
+    }
+
+    fn reference(acc: [i32; 16], a: &[u8; 64], b: &[i8; 64]) -> [i32; 16] {
+        core::array::from_fn(|j| {
+            let sum: i32 = (4 * j..4 * j + 4)
+                .map(|i| i32::from(a[i]) * i32::from(b[i]))
+                .sum();
+            acc[j].wrapping_add(sum)
+        })
+    }
+
+    #[test]
+    fn dot4_matches_reference() {
+        let (us, is) = (bytes(), signed());
+        let mut s = 0x2545_F491_4F6C_DD1Du64;
+        let mut checked = 0;
+        for a in &us {
+            for b in &is {
+                let rand: [i32; 16] = core::array::from_fn(|_| next(&mut s) as i32);
+                for acc in [[0i32; 16], [i32::MAX; 16], [i32::MIN; 16], rand] {
+                    let got = I32x16::from_array(acc)
+                        .dot4_u8_i8_add(U8x64::from_array(*a), I8x64::from_array(*b))
+                        .to_array();
+                    assert_eq!(got, reference(acc, a, b), "acc {acc:?}\na {a:?}\nb {b:?}");
+                    checked += 1;
+                }
+            }
+        }
+        assert!(checked > 10_000, "corpus too small: {checked}");
+    }
+
+    /// A pair of `255 * 127` products is 64_770, which a saturating 16-bit
+    /// pair sum (`vpmaddubsw`) would clamp to 32_767. These lanes would expose
+    /// such an emulation in either sign.
+    #[test]
+    fn dot4_does_not_saturate_pair_sums() {
+        let acc = I32x16::splat(0);
+        let a = U8x64::splat(255);
+        assert_eq!(acc.dot4_u8_i8_add(a, I8x64::splat(127)).to_array(), [4 * 255 * 127; 16]);
+        assert_eq!(acc.dot4_u8_i8_add(a, I8x64::splat(-128)).to_array(), [-4 * 255 * 128; 16]);
+    }
+
+    /// Lane `j` must read bytes `4j..4j+4` and nothing else: a single nonzero
+    /// byte moves exactly one lane.
+    #[test]
+    fn dot4_keeps_bytes_in_their_lane() {
+        for i in 0..64 {
+            let mut a = [0u8; 64];
+            a[i] = 3;
+            let got = I32x16::splat(0)
+                .dot4_u8_i8_add(U8x64::from_array(a), I8x64::splat(-5))
+                .to_array();
+            let want: [i32; 16] = core::array::from_fn(|j| if j == i / 4 { -15 } else { 0 });
+            assert_eq!(got, want, "byte {i}");
+        }
+    }
+
+    #[test]
+    fn nibbles_match_reference() {
+        for a in bytes() {
+            let x = U8x64::from_array(a);
+            assert_eq!(x.low_nibbles().to_array(), a.map(|v| v & 0x0F), "low_nibbles {a:?}");
+            assert_eq!(x.high_nibbles().to_array(), a.map(|v| v >> 4), "high_nibbles {a:?}");
+        }
+        let every: [u8; 64] = core::array::from_fn(|i| (i as u8).wrapping_mul(4).wrapping_add(3));
+        for v in [every, every.map(|x| !x)] {
+            let x = U8x64::from_array(v);
+            assert_eq!(x.high_nibbles().to_array(), v.map(|b| b >> 4));
+            assert_eq!(x.low_nibbles().to_array(), v.map(|b| b & 0x0F));
+        }
+    }
+}