diff --git a/.claude/blackboard.md b/.claude/blackboard.md index 4c67fc2f..e8cc400f 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -3898,3 +3898,25 @@ negative zero; both fixed, every exact kernel now checked against an i64 referen 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. + +## U9 x U9 head-to-head: U8 x U8 x sign x mul_add vs U8 x I8 byte-dot (2026-10-29, branch claude/u9-head-to-head) + +Probe-only, no production change: `examples/u9_head_to_head/{kernels,main}.rs`, `tests/u9_head_to_head.rs`, report +`.claude/knowledge/turbovec9-u9-head-to-head.md`. Both candidates use only facade words and share one i64 oracle: +511 x 511 exhaustive pairs, mixed full blocks, worst-case runs past both flush intervals, never-touched refusal. +Passed on native, v4, v3, avx (qemu), v2 (qemu), aarch64 (qemu), nightly; 12 disable runs all red. +Contract (report § arithmetic contract): unsigned product + polarity XOR as routing + exact fold; A and B are two +facade compositions of it. B is long multiplication, `m = low7 - top`, tested for all 256 magnitudes via `prepare_b`. + +Verdict: B (U8 x I8 byte-dot) is the cheaper exact I9 substrate in every row, symmetric included. B encodes with AND, +XOR and the presence bit only: `low7 = |b| & 0x7F`, `top = |b| ^ low7`, `negative = plane_a ^ plane_b`, four masked +byte-dots, `(pos_low7 - pos_top) - (neg_low7 - neg_top)`; no magnitude is negated. Disassembly (AVX-512): no shuffle, permute, +pack, extract or insert in B, the baselines or the flushes; on AVX2 LLVM puts 4 vpshufb in the facade's masked byte-dot. VNNI L1, ns/64: B end-to-end 3.05-3.06 vs A' 7.64-7.83; +B scan 2.75-2.81 vs A 7.02-7.19; plain byte-dot 1.94. AVX2: B end-to-end 10.1-10.2 vs A' 75-81 (A dominated by +`U8x64::mask_blend`'s per-lane AVX2 realization, not changed). + +Loose ends: the U8x64 -> I8x64 bit view (`bits_as_i8` is a store/load stand-in, 1.16-1.27 ns/block with the split; +not added, it is a u8/i8 crossing); recall on a real TurboVec-U9 index (owner's experiment, defined in the report). +Re-measured after review (Codex, PR #388): timed operands are now compact (72 B/block, not a padded 128-B `Vec`) +and the preparation row consumes every lane. Streamed VNNI B end-to-end 12.9-16.3 -> 6.0-6.5, B scan 5.9-6.9 -> 2.75-2.80; +verdict unchanged. diff --git a/.claude/knowledge/turbovec9-u9-head-to-head.md b/.claude/knowledge/turbovec9-u9-head-to-head.md new file mode 100644 index 00000000..b7e444a6 --- /dev/null +++ b/.claude/knowledge/turbovec9-u9-head-to-head.md @@ -0,0 +1,298 @@ +# TurboVec-U9: U8 x U8 x sign x mul_add against U8 x I8 byte-dot + +READ BY: savant-architect, vector-synthesis, simd-savant, the TurboVec-U9 +owner. Companion to `turbovec9-u9-substrate.md` (#387). + +Probe: `examples/u9_head_to_head/` (`kernels.rs` holds both candidates, +`main.rs` times them). Correctness: `tests/u9_head_to_head.rs`. Both use only +`ndarray::simd` facade words; no backend file changed. + +```sh +cargo run --release --example u9_head_to_head --features std +cargo test --release --test u9_head_to_head +``` + +## The question and the answer + +For full-range signed 9-bit arithmetic (-255..=255) on the shared 8+1 block +(64 `u8` magnitudes + one `u64` positive plane per operand), which composes +the cheaper bit-exact multiply-accumulate through the existing facade? + +**B, U8 x I8 byte-dot, in every row measured.** On AVX-512 VNNI: 3.05 ns per +64 values end to end against 7.64 - 7.83 for the best A, and 2.75 - 2.81 +against 7.02 - 7.19 in the scan form. A, as composed here, widens u8 -> u16 -> f32, +because the facade offers no u16 -> i32 widening to add the products in +integers. That is a cost of this composition, not of U8 x U8 multiplication. + +## The arithmetic contract both candidates implement + +> Sign is a routing decision. Multiplication is an exact decomposition. +> Accumulation is a fold. + +A value is an unsigned magnitude `m` in 0..=255 plus a polarity bit (Codec9: +`P=1, M=0` is numeric zero, `P=0, M=0` is never-touched and is refused). For +one lane pair: + +```text +p = m_a * m_b exact, always fits u16 (255 * 255 = 65_025) +negative = polarity_a XOR polarity_b the sign of the product +acc = acc - p if negative + acc + p otherwise +``` + +The magnitudes are never negated or rewritten. The polarity XOR only picks +which side of the fold the product joins; one subtraction per flush applies +it. + +A and B are two facade compositions of this one contract; neither has a +different numerical meaning, and neither is "more exact" than the other. +They differ only in how the product is broken into pieces the facade can +execute: + +- **A** forms `m_a * m_b` whole (in `f32`, or in `u16` with `mullo`) and adds + it to the positive or the negative sum. +- **B** does it like long multiplication (schriftliche Multiplikation): it + splits `m_b` into two weighted digits the byte-dot can consume, and adds + their partial products. With `low7 = m_b & 0x7F` and `top = m_b ^ low7` + read as a signed byte (0 or -128), `m_b = low7 - top`, so + `m_a * m_b = m_a * low7 - m_a * top`. For `255 x 255`: + `255 * 127 - 255 * (-128) = 32_385 + 32_640 = 65_025`. + +In a scan the query side of the decomposition (`low7`, `top`, its plane) is +prepared once and reused for every index block; only the routing XOR and the +partial products run per block. + +## Encoding rule: AND, XOR and the presence bit, never a shuffle + +No magnitude is ever negated. The sign lives in the presence plane and is +flipped with XOR; kernels apply it as a mask. + +- `negative = plane_a ^ plane_b` (XOR flips the presence bit). +- B splits `|b|` with AND and XOR: `low7 = |b| & 0x7F`, `top = |b| ^ low7` + (= `|b| & 0x80`). Read as signed bytes, `low7` is `|b| & 127` and `top` is + -128 where bit 7 is set, so `|a| * |b| = dot(|a|, low7) - dot(|a|, top)`. +- The masked byte-dot (`dot4_u8_i8_add_masked`, #387) ANDs `|a|` with + `!negative` or `negative`: four byte-dots per block, + `(pos_low7 - pos_top) - (neg_low7 - neg_top)`. +- A routes `|a|` with `U8x64::mask_blend` by the same XOR plane. +- Lane totals leave through a stored array read lane by lane (`black_box` + keeps the compiler from rebuilding a permute reduction). + +Disassembly of the timed closures (AVX-512 VNNI build): no shuffle, permute, +pack, unpack, align, extract or insert instruction in any B row, in either +baseline, or in any flush. A's rows keep `vpmovzxbw`/`vpmovzxwd` and the +upper-half `vextracti64x4`, which the facade's `U16x32::from_u8x64_hi` and +`to_f32x16_*` use to widen. That is how this composition of A lowers through +the current facade; it says nothing about U8 x U8 multiplication itself, nor +about other possible lowerings of A, which were not explored. + +**This holds on AVX-512 only (checked later, in #389).** On an AVX2 build +(`config-v3`) the masked byte-dot that B calls contains 4 `vpshufb`: LLVM +uses them to spread the bytes of the 64-bit keep mask across the register +inside `dot4_u8_i8_add_masked`'s AVX2 realization (the source uses a +variable shift, a multiply and a shift, not a shuffle). That is the facade's +realization, not this probe's code, and it was not changed here. + +## Candidates as built + +| | product | sign | accumulate | flush | +|---|---|---|---|---| +| A | `\|a\| * \|b\|` in f32 after `u8 -> u16 -> f32` (or `U16x32::mullo` then `-> f32`, A') | `\|a\|` routed into a positive or a negative sum by `plane_a ^ plane_b` | `F32x16::mul_add` per quarter, two sums | every 256 blocks (lane sums below 2^24, exact in f32), then `i64` | +| B | four `dot4_u8_i8_add_masked`: `\|a\|` against `low7` and `top` of `\|b\|` | mask `!negative` / `negative` from `plane_a ^ plane_b` | four `I32x16` accumulators | every 16_384 blocks (a lane moves at most 130_560 per block), then `i64` | +| B scan | as B, with `low7` / `top` of the query built once | mask from `plane_x ^ plane_q` | as B | as B | +| A scan | query widened to f32 once | `\|x\|` routed by `plane_x ^ plane_q` | as A | as A | + +Every row also runs `check`, which refuses a block holding the DR512 +never-touched codeword (plane 0, magnitude 0) in either operand, inside the +timed region, for both candidates alike. + +## Correctness + +One independent `i64` oracle (checked arithmetic) for both candidates: + +- all 511 x 511 numeric pairs, each alone in a block at a rotating lane, + through A, A', B and both scan forms; +- full blocks of 64 mixed pairs (catches any lane reordering between the + operands); +- worst-case runs (all magnitudes 255, same and mixed signs) of 49_157 blocks, + past both flush intervals, with totals beyond the `i32` range; +- never-touched refused in either operand, touched +0 and negative nonzero + magnitudes accepted, the oracle refusing the same lanes; +- the split itself, through the kernel's own `prepare_b`, for every magnitude + 0..=255: `low7 XOR top` gives back the original byte, `low7 - top` gives + back its value with `top` in {0, -128}, and changing either polarity bit + changes only the routing plane (`plane_a ^ plane_b`), never `low7` or `top`. + +Passed on native AVX-512 VNNI, config-v4, config-v3, config-avx (qemu), +config-v2 (qemu), aarch64 (qemu) and nightly-simd. Disable runs, each red: +A's routing swapped; A's flush at 4096 blocks (f32 loses exactness past +2^24); B's top limb added instead of subtracted; B's AND mask keeping bit 7; +B's presence plane not flipped by b; B scan ignoring the query plane; B's +flush at 2^20 blocks (i32 wraps); `check` ignoring operand b; the oracle +accepting never-touched; `top` keeping all bits; the plane formed with AND +instead of XOR; polarity leaking into the split bytes. + +CI runs this test file only in the test jobs on its own x86 runners; the +other tiers above were run locally, not in CI. + +### Exactness domain + +Both candidates are bit-exact, with no rounding, saturation or wrap, for: + +- operands from the 511 numeric values -255..=255 (never-touched refused by + `check`, which callers run before the kernel; the kernels do not check); +- any number of blocks whose running total fits `i64`, which at the worst case + of 64 * 65_025 per block is about 2.2 * 10^12 blocks. The final `i64` adds are + unchecked. + +Inside that domain each stage is bounded: + +- **A:** a lane gains one product, at most 65_025, per block. After 256 blocks + that is 16_646_400, below 2^24 = 16_777_216, where every `f32` integer + is exact; then A flushes. +- **B:** an `i32` lane gains at most 4 * 255 * 128 = 130_560 per block. After + 16_384 blocks that is 2_139_095_040, below `i32::MAX` = 2_147_483_647; then + B flushes. The byte-dot does not saturate: the worst-case run (every byte + 255, so `top = -128` in all four bytes of a lane) passes on every tier. +- The worst-case test runs past both limits. + +### Not proven here + +- Exactness beyond the `i64` total above. +- Any lowering of A other than the two measured (f32 `mul_add`, `u16 mullo`). +- The CI matrix running this test on every tier. +- Anything about Codebook512 / Lloyd-Max512 decoding or retrieval parity (see + the handoff). + +## Measurements (ns per 64 values, ranges over 3 runs, this host) + +One Xeon with AVX-512 VNNI. Each tier is a separate build run natively on it +(`--config .cargo/config-v4.toml`, `config-v3.toml`); nothing is projected. +L1 = 64 blocks per operand (4.5 KiB) dotted repeatedly; streamed = 2^20 blocks +per operand (72 MiB). Baseline = plain `dot4_u8_i8_add`; the scan rows compare +with the plain byte-dot against one fixed query. + +Timed input layout: each operand is two compact arrays, 64 magnitude bytes and +one 8-byte plane per block (72 bytes), and each `Block` is built in registers +inside the loop. A `Vec` would pad every block to 128 bytes, because +`U8x64` is 64-byte aligned; the first version of this probe streamed that +padded layout and its streamed rows were inflated by it (B end-to-end 12.9 - +16.3, B scan 5.9 - 6.9 on VNNI). The baseline streams 64 + 64 bytes per block +(no planes), B and A 72 + 72, and a scan only the index side. + +| L1-resident | AVX-512 VNNI | AVX-512, no VNNI | AVX2 | +|---|---|---|---| +| baseline: plain u8 x i8 dot4 | 1.94 | 2.02 - 2.03 | 2.11 - 2.56 | +| A end-to-end (check, route, widen, mul_add) | 8.78 - 8.88 | 8.70 - 8.89 | 74.3 - 77.2 | +| A' end-to-end (check, route, mullo, add) | 7.64 - 7.83 | 7.64 - 8.98 | 75.1 - 81.4 | +| **B end-to-end (check, AND/XOR split, 4 masked dot4)** | **3.05 - 3.06** | **5.61 - 7.23** | **10.1 - 10.2** | +| A arithmetic only | 7.30 - 7.31 | 7.33 - 8.14 | 13.7 - 13.8 | +| A' arithmetic only | 6.16 - 6.22 | 6.16 - 7.69 | 10.7 | +| **B arithmetic only** | **2.59** | **5.05 - 6.17** | **9.31** | +| B operand preparation alone (AND/XOR split) | 1.16 | 1.16 - 1.80 | 1.25 - 1.27 | +| baseline scan: plain dot4, fixed query | 1.88 - 1.92 | 2.09 - 2.27 | 1.59 - 1.61 | +| A scan | 7.02 - 7.19 | 7.05 - 7.39 | 76.5 - 83.5 | +| **B scan** | **2.75 - 2.81** | **4.91 - 4.95** | **9.27 - 9.28** | + +| Streamed | AVX-512 VNNI | AVX-512, no VNNI | AVX2 | +|---|---|---|---| +| baseline | 4.4 - 4.8 | 5.5 - 6.7 | 6.1 - 8.9 | +| A end-to-end | 22.4 - 25.7 | 22.2 - 25.0 | 421 - 426 | +| A' end-to-end | 16.6 - 19.6 | 18.2 - 19.8 | 420 - 427 | +| **B end-to-end** | **6.0 - 6.5** | **9.0 - 12.5** | **26.9 - 30.7** | +| B operand preparation alone | 2.77 - 2.80 | 2.77 - 2.83 | 2.94 - 3.11 | +| baseline scan | 2.23 - 2.30 | 2.46 - 2.52 | 2.60 - 2.95 | +| A scan | 10.6 - 12.0 | 10.7 - 11.2 | 422 - 424 | +| **B scan** | **2.75 - 2.80** | **4.68 - 4.75** | **27.6 - 32.6** | + +The streamed "arithmetic only" rows are left out: they read precomputed +operands (B: one plane plus two 64-byte vectors per block), which more than +doubles the bytes streamed, so they do not isolate arithmetic. + +Notes on the numbers: + +- "Arithmetic only" receives each candidate's execution-form operands + prepared before the timer (A: routed `|a|`; B: the XOR plane and the + `low7` / `top` bytes). The end-to-end and scan rows pay all preparation + inside the timer. +- "B operand preparation alone" passes the whole prepared tuple (plane, + `low7`, `top`) through `black_box` on every block, so no lane can be + dropped as unused; the row includes storing that tuple. The first version + consumed one lane of each vector and read 0.86 - 1.16. +- Keeping the sign out of the bytes costs B four byte-dots per block instead + of two: 2.59 ns against 1.94 for a single byte-dot on VNNI. +- Streamed on VNNI the B scan (2.75 - 2.80) equals its L1 figure: one + index-side block per 64 values does not saturate memory bandwidth, so the + scan stays compute-bound. +- On AVX2 every A row is dominated by `U8x64::mask_blend`, whose AVX2 + realization is a per-lane loop; B does not use it (the masked byte-dot + builds its mask in registers). Streamed AVX2 B rows (27 - 33) are slower + than L1 (9 - 10) by more than the extra bytes explain; not investigated. +- 4096 values per L1 dot: a full dot product costs 64 x the per-64 figure. + +## The one facade gap + +B needs `low7` and `top`, which are `U8x64`, as the `I8x64` operand of the +byte-dot. The facade has no `U8x64 -> I8x64` view, so `kernels::bits_as_i8` +stores the 64 bytes and loads them back as `i8` (every bit kept; no shuffle, +no arithmetic). Its cost is inside "B operand preparation" (1.16 - 1.27 ns per +block with the AND/XOR split, L1). Exact consumer need: a bit view +`U8x64 -> I8x64` (or a byte-dot taking the signed operand as raw bytes). +Recorded, not added: it is a u8/i8 crossing, which the repository rules keep +out of the facade unless decided deliberately. + +## What this means for TurboVec-U9 + +Two contracts share the 8+1 bits and are not the same thing: + +- **Arithmetic511**: signed magnitude, 511 numeric values; `P=0, M=0` is DR512 + never-touched and is refused. This probe proves both candidates exact over + exactly that domain. +- **Codebook512**: 512 centroid IDs, all valid, all distinguishable. 512 IDs + cannot map injectively onto 511 numeric values, so a numeric U9 index has + at most 511 distinct scores per dimension (or 512 patterns if the encoder + is allowed to write `P=0, M=0` as a second zero, which is then a TurboVec + contract, not DR512's, and must not go through `check`). + +Numeric U9 scoring (B scan) computes `sum q_d * x_d` over integer grid values +exactly. It reproduces Lloyd-Max512 scoring only if every centroid is an exact +integer multiple of one scale; arbitrary Lloyd-Max centroids are not. The #387 +probe measured snapping 512 Lloyd-Max centroids onto the 9-bit grid: 329 +distinct values remain, 44.7 dB against 49.9 dB for the exact codebook, and a +uniform 9-bit grid with a tuned clip reaches 46.0 dB. Arithmetic parity here +says nothing about retrieval parity. + +Exact Lloyd-Max512 scoring needs the centroid value of each ID: a 512-entry +decode. The facade offers that only as per-position table reads (`lut_u8` +shape, 38 - 40 ns per 64 IDs in #387); the `vpermi2b` form (about 6 ns) is not +a facade word. Neither cost is hidden in the B scan figure; it is extra. + +## Handoff to the TurboVec-U9 owner + +1. **Winning composition:** B scan. Per query block, once: + `low7 = |q| & 0x7F`, `top = |q| ^ low7`, viewed as `I8x64` + (`kernels::QueryB::new`); the query sign stays in its plane. Per index + block of 64 dims: `negative = index_plane ^ query_plane`, then four + `I32x16::dot4_u8_i8_add_masked(index_mag, !negative or negative, low7 or + top)`; flush into `i64` at most every 16_384 blocks (1,048,576 dims): + `(pos_low7 - pos_top) - (neg_low7 - neg_top)`. +2. **Inputs:** index stored as 64 `u8` magnitudes plus one `u64` positive + plane per 64 dims, passed as `U8x64` and `u64` with no conversion and no + negation. A query that is only 7-bit (|q| <= 127) has `top = 0`; its two + `top` byte-dots can be dropped, leaving the #387 masked pair. +3. **Native numeric U9 index:** yes, the kernel consumes it directly, as long + as the index holds numeric grid values (Arithmetic511), not centroid IDs. +4. **Lloyd-Max512:** needs a decode from ID to centroid, which is either + lossy (snap to the grid; IDs collide) or a lookup per dimension per query + (not amortized by this kernel, and not a fast facade word today). +5. **Unverified:** recall of any 9-bit index on real data; the clip and scale + of the grid; end-to-end search cost with TurboVec's layout, blocking and + rescoring; NEON, wasm and nightly timings (only their correctness ran, + under qemu or without timing). +6. **Smallest falsifying experiment:** in the TurboVec fork, encode one real + benchmark set twice: (a) numeric U9 on a uniform clipped grid, scored with + B scan; (b) Lloyd-Max512 IDs scored exactly in f32 through a decode. Report + recall@10 and ns per query for both against the current i4 path. The + recommendation fails if (a) loses more recall to (b) than the operator + accepts, or if (a)'s scan costs more per dimension than the i4 path by more + than the agreed margin. diff --git a/Cargo.toml b/Cargo.toml index 31bc45af..e7e552f6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -66,6 +66,11 @@ required-features = ["std"] name = "u9_cost_probe" required-features = ["std"] +[[example]] +name = "u9_head_to_head" +path = "examples/u9_head_to_head/main.rs" +required-features = ["std"] + [[example]] name = "cpu_guard_probe" required-features = ["std"] diff --git a/examples/u9_head_to_head/kernels.rs b/examples/u9_head_to_head/kernels.rs new file mode 100644 index 00000000..b0b5d440 --- /dev/null +++ b/examples/u9_head_to_head/kernels.rs @@ -0,0 +1,418 @@ +//! The two U9 x U9 multiply-accumulate candidates, written only with +//! `ndarray::simd` facade words. Shared by `examples/u9_head_to_head` (timing) +//! and `tests/u9_head_to_head.rs` (correctness). +//! +//! Both candidates read the same physical block per operand: 64 `u8` +//! magnitudes and one `u64` plane in which bit `i` set means lane `i` is +//! positive. A lane with plane bit 0 and magnitude 0 is the DR512 +//! never-touched sentinel; both candidates refuse a block that holds one +//! (`check`), so only the 511 numeric values -255..=255 enter the arithmetic +//! (touched zero is plane 1, magnitude 0). +//! +//! The product of lane `i` is negative when exactly one operand is negative: +//! `negative = plane_a ^ plane_b` (both planes mark the positive lanes). +//! +//! No magnitude is ever negated: the sign is the presence bit, flipped with +//! XOR, and every kernel applies it as a mask. Nothing here uses a byte +//! shuffle, permute or pack; lane sums leave through a store, not a +//! horizontal reduction. + +#![allow(dead_code)] + +use ndarray::simd::{F32x16, I32x16, I8x64, U16x32, U8x64}; + +/// One operand block: 64 magnitudes and the positive plane. +#[derive(Clone, Copy)] +pub struct Block { + pub mag: U8x64, + pub positive: u64, +} + +/// One operand in the compact layout the report describes: 64 magnitude bytes +/// and one 8-byte plane per block, in two separate arrays (72 bytes per block). +/// A `Block` holding a 64-byte-aligned `U8x64` next to its plane would pad to +/// 128 bytes, so the timed loops never stream `Block` arrays; they build each +/// `Block` in registers from these two arrays. +pub struct Operand { + pub mag: Vec, + pub plane: Vec, +} + +impl Operand { + pub fn new(v: &[([u8; 64], u64)]) -> Self { + Self { + mag: v.iter().map(|(m, _)| U8x64::from_array(*m)).collect(), + plane: v.iter().map(|(_, p)| *p).collect(), + } + } + + pub fn block(&self, i: usize) -> Block { + Block { + mag: self.mag[i], + positive: self.plane[i], + } + } + + pub fn blocks(&self) -> impl Iterator + '_ { + self.mag + .iter() + .zip(&self.plane) + .map(|(m, p)| Block { mag: *m, positive: *p }) + } +} + +/// A block held a never-touched lane (plane bit 0, magnitude 0); the mask +/// names the offending lanes. +#[derive(Debug, PartialEq, Eq)] +pub struct NeverTouched(pub u64); + +/// Never-touched lanes of one block. +#[inline(always)] +pub fn never_touched(b: Block) -> u64 { + b.mag.cmpeq_mask(U8x64::splat(0)) & !b.positive +} + +/// Refuses a pair of blocks if either holds a never-touched lane. +#[inline(always)] +pub fn check(a: Block, b: Block) -> Result<(), NeverTouched> { + let bad = never_touched(a) | never_touched(b); + if bad == 0 { + Ok(()) + } else { + Err(NeverTouched(bad)) + } +} + +/// The 64 bytes of `v` as four `F32x16` quarters, lanes 0..16, 16..32, +/// 32..48, 48..64. Every byte is exact in `f32`. +#[inline(always)] +fn f32_quarters(v: U8x64) -> [F32x16; 4] { + let (lo, hi) = (U16x32::from_u8x64_lo(v), U16x32::from_u8x64_hi(v)); + [lo.to_f32x16_lo(), lo.to_f32x16_hi(), hi.to_f32x16_lo(), hi.to_f32x16_hi()] +} + +// ── Candidate A: U8 x U8 x sign x mul_add ────────────────────────────────── + +/// A's sign routing: `|a|` where the product is positive, and where it is +/// negative, with `negative = plane_a ^ plane_b`; the other lanes are 0. +#[inline(always)] +pub fn route_a(a: Block, b: Block) -> (U8x64, U8x64) { + let negative = a.positive ^ b.positive; + let zero = U8x64::splat(0); + (U8x64::mask_blend(negative, a.mag, zero), U8x64::mask_blend(!negative, a.mag, zero)) +} + +/// Blocks an A accumulator may absorb before `flush`: each lane gains at most +/// 255 * 255 = 65_025 per block, so 256 blocks stay below 2^24, the range in +/// which every integer, and therefore every partial sum, is exact in `f32`. +pub const A_FLUSH_BLOCKS: usize = 256; + +/// Candidate A state: the unsigned products of positive and of negative +/// lanes, summed separately per lane in `f32` (exact below 2^24). +pub struct AccA { + pos: [F32x16; 4], + neg: [F32x16; 4], + blocks: usize, + total: i64, +} + +impl AccA { + pub fn new() -> Self { + let z = F32x16::splat(0.0); + Self { + pos: [z; 4], + neg: [z; 4], + blocks: 0, + total: 0, + } + } + + /// Moves the exact per-lane sums into the `i64` total and clears them. + #[inline(always)] + pub fn flush(&mut self) { + let mut lanes = I32x16::splat(0); + for k in 0..4 { + // Integers below 2^24: the difference and the conversion are exact. + lanes += (self.pos[k] - self.neg[k]).cast_i32(); + } + self.total += lanes_i64(lanes); + let z = F32x16::splat(0.0); + self.pos = [z; 4]; + self.neg = [z; 4]; + self.blocks = 0; + } + + /// One block, magnitudes kept unsigned: `|a| * |b|` is routed to the + /// positive or the negative sum by `negative = plane_a ^ plane_b`, with one + /// `mul_add` per quarter. Products are taken in `f32` after widening. + #[inline(always)] + pub fn step(&mut self, a: Block, b: Block) { + let (a_pos, a_neg) = route_a(a, b); + self.mac(a_pos, a_neg, b.mag); + } + + /// The arithmetic alone, on magnitudes already routed by sign. + #[inline(always)] + pub fn mac(&mut self, a_pos: U8x64, a_neg: U8x64, b_mag: U8x64) { + let (fp, fneg, fb) = (f32_quarters(a_pos), f32_quarters(a_neg), f32_quarters(b_mag)); + for k in 0..4 { + self.pos[k] = fp[k].mul_add(fb[k], self.pos[k]); + self.neg[k] = fneg[k].mul_add(fb[k], self.neg[k]); + } + self.blocks += 1; + if self.blocks == A_FLUSH_BLOCKS { + self.flush(); + } + } + + /// Same as [`Self::step`], but the products are formed in 16 bits with + /// `U16x32::mullo` (255 * 255 fits `u16`) and only then widened to `f32` + /// for the signed add. + #[inline(always)] + pub fn step_mullo(&mut self, a: Block, b: Block) { + let (a_pos, a_neg) = route_a(a, b); + self.mac_mullo(a_pos, a_neg, b.mag); + } + + /// The arithmetic of [`Self::step_mullo`] alone, on routed magnitudes. + #[inline(always)] + pub fn mac_mullo(&mut self, a_pos: U8x64, a_neg: U8x64, b_mag: U8x64) { + let (b_lo, b_hi) = (U16x32::from_u8x64_lo(b_mag), U16x32::from_u8x64_hi(b_mag)); + let pl = U16x32::from_u8x64_lo(a_pos).mullo(b_lo); + let ph = U16x32::from_u8x64_hi(a_pos).mullo(b_hi); + let nl = U16x32::from_u8x64_lo(a_neg).mullo(b_lo); + let nh = U16x32::from_u8x64_hi(a_neg).mullo(b_hi); + let p = [pl.to_f32x16_lo(), pl.to_f32x16_hi(), ph.to_f32x16_lo(), ph.to_f32x16_hi()]; + let n = [nl.to_f32x16_lo(), nl.to_f32x16_hi(), nh.to_f32x16_lo(), nh.to_f32x16_hi()]; + for k in 0..4 { + self.pos[k] += p[k]; + self.neg[k] += n[k]; + } + self.blocks += 1; + if self.blocks == A_FLUSH_BLOCKS { + self.flush(); + } + } + + pub fn finish(mut self) -> i64 { + self.flush(); + self.total + } +} + +impl Default for AccA { + fn default() -> Self { + Self::new() + } +} + +// ── Candidate B: U8 x I8 byte-dot ────────────────────────────────────────── +// +// No byte is ever negated and nothing is shuffled. The magnitude of b is split +// with AND and XOR, `low7 = |b| & 0x7F` and `top = |b| ^ low7` (= `|b| & 0x80`); +// the sign stays in the presence plane and is flipped with XOR, +// `negative = plane_a ^ plane_b`; the masked byte-dot applies it as an AND +// mask on `|a|`. Read as signed bytes, `low7` is `|b| & 127` and `top` is +// `-128` where bit 7 is set, so `|a| * |b| = dot(|a|, low7) - dot(|a|, top)`. + +/// Blocks a B accumulator may absorb before `flush`: per block an `i32` lane +/// moves by at most `4 * 255 * 128 = 130_560`, so 16_384 blocks stay below +/// `i32::MAX`. The flush sums the 16 lanes in `i64`. +pub const B_FLUSH_BLOCKS: usize = 16_384; + +/// Facade gap: the bytes of `v` as `I8x64`, every bit kept. `ndarray::simd` +/// has no `U8x64 -> I8x64` reinterpretation; this stand-in stores the 64 +/// bytes and loads them back (no shuffle, no arithmetic). The probe needs it +/// only to type `low7` and `top` for the byte-dot's signed operand. +#[inline(always)] +pub fn bits_as_i8(v: U8x64) -> I8x64 { + I8x64::from_array(v.to_array().map(|b| i8::from_ne_bytes([b]))) +} + +/// B's operands for one block pair: the negative-product plane (XOR of the +/// presence planes) and `|b|` split by AND / XOR into `low7` and `top`. +#[inline(always)] +pub fn prepare_b(a: Block, b: Block) -> (u64, I8x64, I8x64) { + let low7 = b.mag & U8x64::splat(0x7F); + let top = b.mag ^ low7; + (a.positive ^ b.positive, bits_as_i8(low7), bits_as_i8(top)) +} + +/// Sum of the 16 lanes in `i64`: the lanes are stored and each one is read +/// where it sits. `black_box` keeps the array in memory and the 16 adds +/// scalar; without it the compiler rebuilds a horizontal reduction from +/// permutes and lane extracts. +#[inline(always)] +pub fn lanes_i64(v: I32x16) -> i64 { + let lanes = core::hint::black_box(v.to_array()); + lanes + .iter() + .map(|&x| i64::from(core::hint::black_box(x))) + .sum() +} + +/// Candidate B state: four masked byte-dot accumulators, +/// `[pos x low7, pos x top, neg x low7, neg x top]`; the result is +/// `(pos_low7 - pos_top) - (neg_low7 - neg_top)`. +pub struct AccB { + acc: [I32x16; 4], + blocks: usize, + total: i64, +} + +impl AccB { + pub fn new() -> Self { + Self { + acc: [I32x16::splat(0); 4], + blocks: 0, + total: 0, + } + } + + #[inline(always)] + pub fn flush(&mut self) { + let r = self.acc.map(lanes_i64); + self.total += (r[0] - r[1]) - (r[2] - r[3]); + self.acc = [I32x16::splat(0); 4]; + self.blocks = 0; + } + + /// The arithmetic alone: `|a|` against the split of `|b|`, positive and + /// negative products separated by the plane `negative`. + #[inline(always)] + pub fn mac(&mut self, a_mag: U8x64, negative: u64, low7: I8x64, top: I8x64) { + let (p, n) = (!negative, negative); + self.acc[0] = self.acc[0].dot4_u8_i8_add_masked(a_mag, p, low7); + self.acc[1] = self.acc[1].dot4_u8_i8_add_masked(a_mag, p, top); + self.acc[2] = self.acc[2].dot4_u8_i8_add_masked(a_mag, n, low7); + self.acc[3] = self.acc[3].dot4_u8_i8_add_masked(a_mag, n, top); + self.blocks += 1; + if self.blocks == B_FLUSH_BLOCKS { + self.flush(); + } + } + + /// One block from the shared representation: [`prepare_b`], then [`Self::mac`]. + #[inline(always)] + pub fn step(&mut self, a: Block, b: Block) { + let (negative, low7, top) = prepare_b(a, b); + self.mac(a.mag, negative, low7, top); + } + + /// One index block against a prepared query: the sign is the XOR of the + /// two presence planes, applied as the byte-dot mask. + #[inline(always)] + pub fn scan_step(&mut self, x: Block, q: &QueryB) { + self.mac(x.mag, x.positive ^ q.positive, q.low7, q.top); + } + + pub fn finish(mut self) -> i64 { + self.flush(); + self.total + } +} + +impl Default for AccB { + fn default() -> Self { + Self::new() + } +} + +/// B's query: `|q|` split by AND / XOR into `low7` and `top`, prepared once; +/// the query's sign stays in its plane. +pub struct QueryB { + low7: I8x64, + top: I8x64, + positive: u64, +} + +impl QueryB { + pub fn new(q: Block) -> Self { + let low7 = q.mag & U8x64::splat(0x7F); + let top = q.mag ^ low7; + Self { + low7: bits_as_i8(low7), + top: bits_as_i8(top), + positive: q.positive, + } + } +} + +// ── Scan forms: the query operand prepared once, the index operand per block ─ +// +// In a search the query block is fixed while index blocks stream past, so each +// candidate may prepare the query side once and amortize it. Both scan forms +// get exactly that and nothing more; the index side is read from the shared +// U8 + plane representation every block. + +/// A's query: `|q|` as four `f32` quarters, plus its plane. +pub struct QueryA { + f: [F32x16; 4], + positive: u64, +} + +impl QueryA { + pub fn new(q: Block) -> Self { + Self { + f: f32_quarters(q.mag), + positive: q.positive, + } + } +} + +impl AccA { + /// One index block against a prepared query: route `|x|` by + /// `negative = plane_x ^ plane_q`, widen it, and `mul_add` with `|q|`. + #[inline(always)] + pub fn scan_step(&mut self, x: Block, q: &QueryA) { + let negative = x.positive ^ q.positive; + let zero = U8x64::splat(0); + let x_pos = f32_quarters(U8x64::mask_blend(negative, x.mag, zero)); + let x_neg = f32_quarters(U8x64::mask_blend(!negative, x.mag, zero)); + for k in 0..4 { + self.pos[k] = x_pos[k].mul_add(q.f[k], self.pos[k]); + self.neg[k] = x_neg[k].mul_add(q.f[k], self.neg[k]); + } + self.blocks += 1; + if self.blocks == A_FLUSH_BLOCKS { + self.flush(); + } + } +} + +// ── Baseline and oracle ──────────────────────────────────────────────────── + +/// Plain u8 x i8 byte-dot, the substrate both candidates are measured +/// against: one `dot4_u8_i8_add` per block. +#[inline(always)] +pub fn baseline_step(acc: I32x16, a: U8x64, w: I8x64) -> I32x16 { + acc.dot4_u8_i8_add(a, w) +} + +/// Signed value of lane `i` of a block, or `None` for never-touched. +pub fn value(mag: &[u8; 64], positive: u64, i: usize) -> Option { + let m = i64::from(mag[i]); + match (positive >> i & 1 == 1, m) { + (false, 0) => None, + (true, _) => Some(m), + (false, _) => Some(-m), + } +} + +/// Independent `i64` oracle over whole blocks; refuses never-touched lanes. +pub fn oracle(a: &[([u8; 64], u64)], b: &[([u8; 64], u64)]) -> Result { + let mut sum = 0i64; + for ((am, ap), (bm, bp)) in a.iter().zip(b) { + let mut bad = 0u64; + for i in 0..64 { + match (value(am, *ap, i), value(bm, *bp, i)) { + (Some(x), Some(y)) => sum = sum.checked_add(x.checked_mul(y).unwrap()).unwrap(), + _ => bad |= 1 << i, + } + } + if bad != 0 { + return Err(NeverTouched(bad)); + } + } + Ok(sum) +} diff --git a/examples/u9_head_to_head/main.rs b/examples/u9_head_to_head/main.rs new file mode 100644 index 00000000..78f037eb --- /dev/null +++ b/examples/u9_head_to_head/main.rs @@ -0,0 +1,240 @@ +//! U9 x U9 head-to-head through the `ndarray::simd` facade: candidate A +//! (U8 x U8, sign routed, `mul_add`) against candidate B (U8 x I8 byte-dot on +//! a signed 7+1 split), plus the plain u8 x i8 byte-dot baseline. +//! +//! Every kernel is checked against the `i64` oracle on the same data before +//! it is timed. `cargo run --release --example u9_head_to_head --features std`; +//! the tier is whatever this build compiled (printed in the header). + +mod kernels; + +use kernels::*; +use ndarray::simd::{I32x16, I8x64, U8x64}; +use std::hint::black_box; +use std::time::Instant; + +fn tier() -> &'static str { + if cfg!(target_feature = "avx512vnni") { + "x86_64 AVX-512 + VNNI" + } else if cfg!(target_feature = "avx512bw") { + "x86_64 AVX-512 (no VNNI)" + } else if cfg!(target_feature = "avx2") { + "x86_64 AVX2" + } else if cfg!(target_feature = "avx") { + "x86_64 AVX" + } else if cfg!(target_arch = "aarch64") { + "aarch64 NEON" + } else { + "scalar" + } +} + +struct Data { + a: Vec<([u8; 64], u64)>, + b: Vec<([u8; 64], u64)>, +} + +/// Uniform over the 511 numeric values; never-touched never appears. +fn data(blocks: usize, seed: u64) -> Data { + let mut s = seed; + 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 mut mk = || { + let mut m = [0u8; 64]; + let mut p = 0u64; + for (i, mi) in m.iter_mut().enumerate() { + let v = (next() % 511) as i32 - 255; + *mi = v.unsigned_abs() as u8; + p |= u64::from(v >= 0) << i; + } + (m, p) + }; + let a = (0..blocks).map(|_| mk()).collect(); + let b = (0..blocks).map(|_| mk()).collect(); + Data { a, b } +} + +/// Median of `runs` timings of `f` over `blocks` blocks, in ns per 64 values. +fn time(name: &str, blocks: usize, runs: usize, baseline: Option, f: &dyn Fn() -> i64) -> f64 { + black_box(f()); + let mut t: Vec = (0..runs) + .map(|_| { + let t0 = Instant::now(); + black_box(f()); + t0.elapsed().as_nanos() as f64 / blocks as f64 + }) + .collect(); + t.sort_by(f64::total_cmp); + let (min, med) = (t[0], t[runs / 2]); + let rel = baseline + .map(|b| format!(" {:5.2}x baseline", med / b)) + .unwrap_or_default(); + println!("{name:<52} {med:7.3} ns/64 (min {min:7.3}){rel}"); + med +} + +fn bench(label: &str, d: &Data, runs: usize) { + let n = d.a.len(); + let (ab, bb) = (Operand::new(&d.a), Operand::new(&d.b)); + let want = oracle(&d.a, &d.b).unwrap(); + + // Baseline operand: the b magnitudes reinterpreted as signed bytes are not + // the same numbers, but the byte-dot costs the same for any bytes. Built + // once, outside the timed region, because the baseline has no preparation. + let w: Vec = + d.b.iter() + .map(|(m, _)| I8x64::from_array(m.map(|x| i8::from_ne_bytes([x])))) + .collect(); + // Execution-form operands for the arithmetic-only rows, also built once. + let routed: Vec<(U8x64, U8x64)> = ab + .blocks() + .zip(bb.blocks()) + .map(|(a, b)| route_a(a, b)) + .collect(); + let split: Vec<(u64, I8x64, I8x64)> = ab + .blocks() + .zip(bb.blocks()) + .map(|(a, b)| prepare_b(a, b)) + .collect(); + + let e2e_a = || { + let mut acc = AccA::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.step(a, b); + } + acc.finish() + }; + let e2e_a2 = || { + let mut acc = AccA::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.step_mullo(a, b); + } + acc.finish() + }; + let e2e_b = || { + let mut acc = AccB::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.step(a, b); + } + acc.finish() + }; + let ari_a = || { + let mut acc = AccA::new(); + for ((p, q), b) in routed.iter().zip(&bb.mag) { + acc.mac(*p, *q, *b); + } + acc.finish() + }; + let ari_a2 = || { + let mut acc = AccA::new(); + for ((p, q), b) in routed.iter().zip(&bb.mag) { + acc.mac_mullo(*p, *q, *b); + } + acc.finish() + }; + let ari_b = || { + let mut acc = AccB::new(); + for (a, (neg, low7, top)) in ab.mag.iter().zip(&split) { + acc.mac(*a, *neg, *low7, *top); + } + acc.finish() + }; + let base = || { + let mut acc = I32x16::splat(0); + let mut total = 0i64; + for (k, (a, w)) in ab.mag.iter().zip(&w).enumerate() { + acc = baseline_step(acc, *a, *w); + if k % B_FLUSH_BLOCKS == B_FLUSH_BLOCKS - 1 { + total += lanes_i64(acc); + acc = I32x16::splat(0); + } + } + total + lanes_i64(acc) + }; + let prep_b = || { + // The whole prepared tuple goes through the barrier, so no lane of + // `low7` or `top` can be dropped as unused. + for (a, b) in ab.blocks().zip(bb.blocks()) { + black_box(prepare_b(a, b)); + } + 0 + }; + // Scan forms: operand b plays the query. One query block is reused for + // every index block, as in a search; its preparation runs once per call + // and is inside the timed region. + let q = bb.block(0); + let scan_a = || { + let qa = QueryA::new(q); + let mut acc = AccA::new(); + for x in ab.blocks() { + check(x, q).unwrap(); + acc.scan_step(x, &qa); + } + acc.finish() + }; + let scan_b = || { + let qb = QueryB::new(q); + let mut acc = AccB::new(); + for x in ab.blocks() { + check(x, q).unwrap(); + acc.scan_step(x, &qb); + } + acc.finish() + }; + let qs: Vec<_> = vec![d.b[0]; n]; + let want_scan = oracle(&d.a, &qs).unwrap(); + assert_eq!(scan_a(), want_scan, "A scan disagrees with the oracle"); + assert_eq!(scan_b(), want_scan, "B scan disagrees with the oracle"); + for (name, f) in [ + ("A", &e2e_a as &dyn Fn() -> i64), + ("A mullo", &e2e_a2), + ("B", &e2e_b), + ("A arith", &ari_a), + ("A mullo arith", &ari_a2), + ("B arith", &ari_b), + ] { + assert_eq!(f(), want, "{name} disagrees with the oracle on the benchmark data"); + } + + println!("-- {label}: {n} blocks of 64 values per operand ({} KiB per operand incl. plane)", n * 72 / 1024); + let b0 = time("baseline: plain u8 x i8 dot4 (one byte-dot)", n, runs, None, &base); + time("A end-to-end: check + route + widen + mul_add", n, runs, Some(b0), &e2e_a); + time("A' end-to-end: check + route + mullo + add", n, runs, Some(b0), &e2e_a2); + time("B end-to-end: check + AND/XOR split + 4 masked dot4", n, runs, Some(b0), &e2e_b); + time("A arithmetic only (routed operands given)", n, runs, Some(b0), &ari_a); + time("A' arithmetic only (routed operands given)", n, runs, Some(b0), &ari_a2); + time("B arithmetic only (split operands given)", n, runs, Some(b0), &ari_b); + time("B operand preparation alone (AND/XOR split)", n, runs, Some(b0), &prep_b); + let base_scan = || { + let w0 = w[0]; + let mut acc = I32x16::splat(0); + let mut total = 0i64; + for (k, a) in ab.mag.iter().enumerate() { + acc = baseline_step(acc, *a, w0); + if k % B_FLUSH_BLOCKS == B_FLUSH_BLOCKS - 1 { + total += lanes_i64(acc); + acc = I32x16::splat(0); + } + } + total + lanes_i64(acc) + }; + let bs = time("baseline scan: plain dot4, one fixed query", n, runs, None, &base_scan); + time("A scan: query widened once, index routed + mul_add", n, runs, Some(bs), &scan_a); + time("B scan: query split once, 4 masked dot4 (#387)", n, runs, Some(bs), &scan_b); +} + +fn main() { + println!("tier: {} (target-cpu as compiled)", tier()); + // 64 blocks = 4.5 KiB per operand: L1-resident. One dot product of 4096 + // values, repeated. + bench("L1-resident, repeated 4096-value dot", &data(64, 1), 4001); + // 2^20 blocks = 72 MiB per operand: streamed from DRAM. + bench("streamed", &data(1 << 20, 2), 9); +} diff --git a/tests/u9_head_to_head.rs b/tests/u9_head_to_head.rs new file mode 100644 index 00000000..9eeed3ee --- /dev/null +++ b/tests/u9_head_to_head.rs @@ -0,0 +1,186 @@ +#![cfg(feature = "std")] +//! Correctness of the two U9 x U9 candidates in +//! `examples/u9_head_to_head/kernels.rs` against one independent `i64` oracle. + +#[path = "../examples/u9_head_to_head/kernels.rs"] +mod kernels; + +use kernels::{check, oracle, AccA, AccB, Block, NeverTouched, QueryA, QueryB, A_FLUSH_BLOCKS, B_FLUSH_BLOCKS}; +use ndarray::simd::U8x64; + +fn block(mag: [u8; 64], positive: u64) -> Block { + Block { + mag: U8x64::from_array(mag), + positive, + } +} + +/// (magnitude, positive bit) of a numeric value; 0 is touched +0. +fn encode(v: i32) -> (u8, bool) { + (v.unsigned_abs() as u8, v >= 0) +} + +fn run_a(a: &[([u8; 64], u64)], b: &[([u8; 64], u64)], mullo: bool) -> i64 { + let mut acc = AccA::new(); + for ((am, ap), (bm, bp)) in a.iter().zip(b) { + let (x, y) = (block(*am, *ap), block(*bm, *bp)); + check(x, y).unwrap(); + if mullo { + acc.step_mullo(x, y) + } else { + acc.step(x, y) + } + } + acc.finish() +} + +fn run_b(a: &[([u8; 64], u64)], b: &[([u8; 64], u64)]) -> i64 { + let mut acc = AccB::new(); + for ((am, ap), (bm, bp)) in a.iter().zip(b) { + let (x, y) = (block(*am, *ap), block(*bm, *bp)); + check(x, y).unwrap(); + acc.step(x, y); + } + acc.finish() +} + +/// Scan forms: `b` is the query, prepared once per query block; `a` streams. +fn run_scan(a: &[([u8; 64], u64)], b: &[([u8; 64], u64)]) -> (i64, i64) { + let (mut sa, mut sb) = (AccA::new(), AccB::new()); + for ((am, ap), (bm, bp)) in a.iter().zip(b) { + let (x, q) = (block(*am, *ap), block(*bm, *bp)); + check(x, q).unwrap(); + sa.scan_step(x, &QueryA::new(q)); + sb.scan_step(x, &QueryB::new(q)); + } + (sa.finish(), sb.finish()) +} + +/// Every one of the 511 x 511 numeric pairs, each alone in its block (the +/// other 63 lanes are touched +0), at a lane that rotates with the pair so +/// every lane position is exercised. +#[test] +fn all_511_by_511_pairs_are_exact() { + let mut n = 0usize; + for va in -255..=255i32 { + for vb in -255..=255i32 { + let lane = n % 64; + let (ma, pa) = encode(va); + let (mb, pb) = encode(vb); + let mut am = [0u8; 64]; + let mut bm = [0u8; 64]; + am[lane] = ma; + bm[lane] = mb; + let ap = if pa { u64::MAX } else { !(1u64 << lane) }; + let bp = if pb { u64::MAX } else { !(1u64 << lane) }; + let (a, b) = ([(am, ap)], [(bm, bp)]); + let want = i64::from(va) * i64::from(vb); + assert_eq!(oracle(&a, &b), Ok(want), "oracle {va} x {vb}"); + assert_eq!(run_a(&a, &b, false), want, "A {va} x {vb}"); + assert_eq!(run_a(&a, &b, true), want, "A mullo {va} x {vb}"); + assert_eq!(run_b(&a, &b), want, "B {va} x {vb}"); + assert_eq!(run_scan(&a, &b), (want, want)); + n += 1; + } + } + assert_eq!(n, 511 * 511); +} + +/// Full blocks: 64 different pairs per block, so a lane mix-up between the +/// operands (an order-changing widening or narrowing) cannot cancel out. +#[test] +fn full_blocks_match_the_oracle() { + 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 mut val = || (next() % 511) as i32 - 255; + for blocks in [1usize, 2, 7, 64] { + let mk = |val: &mut dyn FnMut() -> i32| { + let mut m = [0u8; 64]; + let mut p = 0u64; + for (i, slot) in m.iter_mut().enumerate() { + let (mi, pi) = encode(val()); + *slot = mi; + p |= u64::from(pi) << i; + } + (m, p) + }; + let a: Vec<_> = (0..blocks).map(|_| mk(&mut val)).collect(); + let b: Vec<_> = (0..blocks).map(|_| mk(&mut val)).collect(); + let want = oracle(&a, &b).unwrap(); + assert_eq!(run_a(&a, &b, false), want, "A, {blocks} blocks"); + assert_eq!(run_a(&a, &b, true), want, "A mullo, {blocks} blocks"); + assert_eq!(run_b(&a, &b), want, "B, {blocks} blocks"); + assert_eq!(run_scan(&a, &b), (want, want)); + } +} + +/// Worst-case magnitudes past both flush intervals, all one sign and mixed: +/// nothing wraps, rounds or saturates. +#[test] +fn worst_case_runs_cross_the_flush_bounds_exactly() { + let n = 3 * B_FLUSH_BLOCKS + 5; + assert!(n > 4 * A_FLUSH_BLOCKS); + for (pa, pb) in [(u64::MAX, u64::MAX), (u64::MAX, 0), (0, 0), (0xF0F0_F0F0_0F0F_0F0F, 0x3333_CCCC_5555_AAAA)] { + let a = vec![([255u8; 64], pa); n]; + let b = vec![([255u8; 64], pb); n]; + let want = oracle(&a, &b).unwrap(); + assert!(want.unsigned_abs() > i32::MAX as u64 || pa != pb, "the run must leave the i32 range"); + assert_eq!(run_a(&a, &b, false), want); + assert_eq!(run_a(&a, &b, true), want); + assert_eq!(run_b(&a, &b), want); + assert_eq!(run_scan(&a, &b), (want, want)); + } +} + +/// The never-touched lane (plane 0, magnitude 0) is refused in either +/// operand; touched +0 (plane 1, magnitude 0) is accepted. +#[test] +fn never_touched_is_refused_and_touched_zero_is_not() { + let ok = block([0u8; 64], u64::MAX); + assert_eq!(check(ok, ok), Ok(())); + for lane in [0usize, 17, 63] { + let mut m = [7u8; 64]; + m[lane] = 0; + let sentinel = block(m, !(1u64 << lane)); + assert_eq!(check(sentinel, ok), Err(NeverTouched(1 << lane))); + assert_eq!(check(ok, sentinel), Err(NeverTouched(1 << lane))); + // Negative nonzero magnitudes are numbers, not the sentinel. + let negative = block([7u8; 64], 0); + assert_eq!(check(negative, ok), Ok(())); + assert_eq!(oracle(&[(m, !(1u64 << lane))], &[([0u8; 64], u64::MAX)]), Err(NeverTouched(1 << lane))); + } +} + +/// The schoolbook split behind B, for every magnitude 0..=255, through the +/// kernel's own `prepare_b`: the two partial operands recombine to the +/// original byte by XOR (no bit lost, nothing negated), the signed value is +/// `|b| = low7 - top` with `top` in {0, -128}, and flipping either polarity +/// bit changes only the routing plane, never the byte operands. +#[test] +fn split_recombines_and_polarity_only_routes() { + let mag: [u8; 64] = core::array::from_fn(|i| i as u8); + for base in [0u8, 64, 128, 192] { + let m: [u8; 64] = core::array::from_fn(|i| mag[i].wrapping_add(base)); + let a = block([1u8; 64], u64::MAX); + let (plane, low7, top) = kernels::prepare_b(a, block(m, u64::MAX)); + assert_eq!(plane, 0, "equal polarities route every lane positive"); + let (low7, top) = (low7.to_array(), top.to_array()); + for i in 0..64 { + let (l, t) = (low7[i], top[i]); + assert!(t == 0 || t == -128, "top of {} is {t}", m[i]); + assert_eq!(l.to_ne_bytes()[0] ^ t.to_ne_bytes()[0], m[i], "XOR recombination of {}", m[i]); + assert_eq!(i32::from(l) - i32::from(t), i32::from(m[i]), "low7 - top of {}", m[i]); + } + for (pa, pb) in [(0u64, u64::MAX), (u64::MAX, 0), (0x5555_5555_5555_5555, 0x0F0F_0F0F_0F0F_0F0F)] { + let (flipped, l2, t2) = kernels::prepare_b(block([1u8; 64], pa), block(m, pb)); + assert_eq!(flipped, pa ^ pb, "the plane is the XOR of the two polarities"); + assert_eq!(l2.to_array(), low7, "polarity must not touch low7"); + assert_eq!(t2.to_array(), top, "polarity must not touch top"); + } + } +}