diff --git a/.claude/blackboard.md b/.claude/blackboard.md index e8cc400f..5e44199b 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -3920,3 +3920,20 @@ not added, it is a u8/i8 crossing); recall on a real TurboVec-U9 index (owner's 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. + +## D-U9-ABACUS: schoolbook U9 x U9 on bit planes (2026-10-29, branch claude/u9-abacus, stacked on #388) + +Probe-only: `examples/u9_abacus/{abacus,main}.rs`, `tests/u9_abacus.rs`, report +`.claude/knowledge/u9-abacus-probe.md`. `carry_save3` (exact full adder on 64 bit columns) ADOPT-LAB; byte-abacus +as a B replacement REFUSE: best fold (per-plane popcount) 12.4 ns/64 arithmetic with free planes vs B 2.9 end-to-end +on VNNI, carry-save folds ~290 ns (LLVM does not unroll them; scalar, no vpternlog). No break-even for prepared +operands. On v4/AVX2 LLVM lowers `count_ones` to a `vpshufb` nibble table. Exact on all 7 tiers; 9 disable runs red. +D (A' products + 16-bit place-value column, carry by bit logic, no f32): react beats A' ~10 % on AVX-512 (5.2-5.7 vs +6.3-6.5 ns arith L1; LLVM folds the carry into 2 vpternlogq), anticipate wins on AVX2; B still 2.3x faster on VNNI because D +widens u8->u16 before mullo. ADOPT-LAB as A's f32 replacement, DEFER as B replacement. Exact on 7 tiers; 5 disables red. +E (B's split at 128, anchor term as masked sum << 7): E3 (3 byte-dots, +1/-1 sign bytes) ties B on VNNI (2.84-2.92 vs 2.81-2.89 e2e, +2.28 vs 2.34-2.41 arith), E4/Esum slower, AVX2 lost to mask_blend. No adoption. Facade gaps recorded: vertical masked vpsadbw +accumulate; sign plane -> +1/-1 bytes without mask_blend. Exact on 7 tiers; 5 disables red. +Re-measured after review (Codex, PR #389): compact timed operands (as #388) and the plane-extraction row consumes all 16 +planes (VNNI L1 2.52 -> 3.5 ns). Verdicts unchanged; streamed D react arith now level with A' (both memory-bound). + diff --git a/.claude/knowledge/u9-abacus-probe.md b/.claude/knowledge/u9-abacus-probe.md new file mode 100644 index 00000000..a588487c --- /dev/null +++ b/.claude/knowledge/u9-abacus-probe.md @@ -0,0 +1,362 @@ +# D-U9-ABACUS: schoolbook U9 x U9 with carry escalation + +Two different experiments, with two different verdicts. Do not read one +for the other: + +- **C, bit-plane multiplication** (`abacus.rs`): the product itself rebuilt + from binary partial products, compressed with carry-save logic. + **Refused** as a replacement for B. +- **D, the 64K carry boundary** (`carry16.rs`): the ordinary U8 x U8 -> U16 + product kept, added to a wrapping 16-bit word, the carry escalated into a + high tally. **Adopted (lab) as the exact replacement for A's f32 path**; + it beats A' on AVX-512 and loses to B only because of its multiply. +- **E, the anchor at 128** (`anchor.rs`): B's own split `|b| = low7 + 128 * + bit7`, with the anchor term a masked sum of `a` shifted by 7 instead of a + byte-dot against -128. **Ties B on AVX-512; no gain worth adopting.** + +READ BY: savant-architect, vector-synthesis, simd-savant, the TurboVec-U9 +owner. Follows `turbovec9-u9-head-to-head.md` (#388), whose B kernel and +`i64` oracle this probe reuses unchanged. + +Probe: `examples/u9_abacus/` (`abacus.rs` holds the arithmetic, `main.rs` +times it against #388's B). Correctness: `tests/u9_abacus.rs`. Plain Rust bit +operations and facade words only; no backend, facade or storage change. + +```sh +cargo run --release --example u9_abacus --features std +cargo test --release --test u9_abacus +``` + +## Candidate E: anchored at the power of two inside the operand + +No U9 operand reaches 64K, but every one has a known power-of-two anchor: +bit 7. B already splits there; E keeps the multiply only for the residual: + +`a * |b| = a * low7 + (sum of a over lanes with bit 7 of b set) << 7` + +The anchor term is a masked sum, weighted by 128 once, at flush. Three ways +to form it from existing facade words: + +| | per block | anchor term | +|---|---|---| +| B (#388) | 4 masked byte-dots | `a` x `top` (-128), a full byte-dot | +| E4 | 4 masked byte-dots | `a` x all-ones bytes, masked by sign AND bit 7, one per sign | +| E3 | 3 masked byte-dots | `a` x a `+1/-1` byte vector (from the sign plane via `mask_blend`), masked by bit 7 | +| Esum | 2 masked byte-dots + 2 `mask_blend` + 2 `sum_bytes_u64` | masked byte sums, reduced horizontally every block | + +Measured (ns per 64 values, 3 runs; L1 / streamed): + +| | AVX-512 VNNI | AVX-512, no VNNI | AVX2 | +|---|---|---|---| +| B end-to-end | 2.81 - 2.89 / 5.46 - 6.96 | 5.62 / 8.85 - 13.0 | 10.0 - 10.1 / 24.0 - 25.2 | +| **E3 end-to-end** | **2.84 - 2.92** / 5.35 - 6.91 | **5.52** / 7.28 - 12.9 | **37.6 - 43.6** / 217 - 222 | +| E4 end-to-end | 3.47 - 3.56 / 5.67 - 7.25 | 6.47 - 6.48 / 8.63 - 15.5 | 14.5 / 38.6 - 43.1 | +| Esum end-to-end | 4.33 - 4.44 / 6.51 - 8.95 | 5.95 / 8.08 - 13.5 | 61.2 - 66.8 / 278 - 279 | +| B arithmetic only | 2.34 - 2.41 / 16.4 - 18.1 | 5.02 / 22.4 - 23.8 | 9.17 - 9.23 / 27.3 - 29.3 | +| **E3 arithmetic only** | **2.28** / 8.29 - 10.3 | **4.94 - 4.95** / 14.6 - 17.7 | **37.5 - 45.3** / 217 - 230 | +| E4 arithmetic only | 2.67 - 2.73 / 11.0 - 13.1 | 5.81 / 16.0 - 18.9 | 12.9 / 25.0 - 25.4 | +| B scan | 2.66 / 2.85 - 3.10 | 4.73 - 4.77 / 5.31 - 8.01 | 9.16 - 9.17 / 26.1 - 29.1 | +| E3 scan | 2.64 / 2.83 - 3.74 | 4.95 - 4.97 / 4.89 - 5.66 | 35.5 - 43.3 / 213 - 214 | + +Reading it: + +- **E3 saves one byte-dot and pays it back building the sign vector.** On + VNNI it ties B end to end (2.84 - 2.92 vs 2.81 - 2.89) and in the scan, and is + about 4 % faster in arithmetic only (2.28 vs 2.34 - 2.41); without VNNI it is + about 2 % faster end to end and slower in the scan. Disassembly: 3 + `vpdpbusd` against B's 4, plus the `kmovq` and blend that make `+1/-1`. +- **E4 is slower everywhere.** A byte-dot against constant ones costs the + same instruction as one against -128; the extra mask ANDs are pure cost. +- **Esum is slower everywhere.** `sum_bytes_u64` reduces horizontally on + every call. +- **On AVX2, E3 and Esum lose to `mask_blend`'s per-lane realization** + (38 - 67 ns); E4, which needs no blend, is 1.4x B. + +Verdict for E: **no adoption.** The anchor idea is arithmetically right and +exact, but B already turns the anchor into one hardware instruction; the +saving is one byte-dot out of four, and every facade route to the masked sum +costs about that much. Two facade gaps would change it, both recorded and not +added: + +1. **A vertical masked byte sum** (`vpsadbw` into a `U64x8` accumulator, + reduced once at flush): the anchor term per block would be one SAD and + one add instead of a byte-dot. +2. **A sign plane to `+1/-1` bytes** that does not go through `mask_blend` + (which is a per-lane loop on AVX2). + +Correctness: the split `|b| = low7 + 128 * bit7` for every magnitude and the +`+1/-1` sign bytes for a mixed plane; all 511 x 511 pairs through E4, E3, +Esum and all three scans; full blocks; all-255 runs past the flush bound. +Passed on all seven tiers. Disable runs, each red: anchor weighted `<< 6`; +E4's anchor ignoring bit 7; sign bytes all +1; Esum's blend inverted; `low7` +keeping bit 7. + +Also found while checking E: on AVX2 the facade's masked byte-dot itself +contains 4 `vpshufb` (LLVM spreads the keep-mask bytes with them), so #388's +"no shuffle in B" holds on AVX-512 only. Corrected in #388's report. + +## Verdict for C (bit-plane multiplication only) + +- **`carry_save3`: ADOPT-LAB.** An exact, five-operation full adder over 64 + independent bit columns; LLVM lowers it to `xor`/`and` in general registers. +- **Bit-plane multiplication as a replacement for B: REFUSE.** No C variant + comes close. + The best one (one popcount per partial-product plane, no carry-save) is + 5.2x B's arithmetic on AVX-512 VNNI with both operands' planes prepared for + free, and the carry-save variants are about 120x. There is no break-even + for prepared-operand reuse: with planes given at no cost, C is still slower + than B paying all of its own preparation. +- **Deferring C's carry-save normalization: not a gain.** Normalizing once + at the end is no faster than normalizing every block, and slower when + streamed. This says nothing about D, whose carry is resolved on every add. + +The arithmetic is exact in every variant; the loss is purely in instruction +count and in what LLVM makes of the carry-save state. + +## Candidate D: the 16-bit place-value column (the narrower experiment) + +The question after the bit-plane result: keep #388's U8 multiplication as it +is (A': `U16x32::mullo`, sign routed by the polarity XOR) and replace only its +`f32` detour with exact positional accumulation. A product is below 2^16, so +it fits a `u16` lane. Each tally lane keeps `low` (wrapping `u16`) and `high` +(units of 2^16); the value is `low + 2^16 * high`. Crossing 2^16 is one carry +into `high`. 67,025 is `1,489 + 1 * 65,536`: the low word wraps, one unit +moves up. Nothing is widened before it must be. (`carry16.rs`.) + +Two ways to find the carry, both plain bit logic on `U16x32` lanes: + +- **react** (`carry_after`): add first, then read the carry out of the sum, + `msb((low & p) | ((low | p) & !sum))`; +- **anticipate** (`carry_before`): before adding, the carry is set exactly + when `p` exceeds the headroom `!low = 65535 - low`, the borrow out of + `!low - p`: `msb((low & p) | (!(!low ^ p) & (!low - p)))`. + +`high` gains at most one unit per block per lane, so it is flushed into the +`i64` total every 65,535 blocks. A flush moves `low + 2^16 * high` per lane. +Carries are counted, never XORed together: two carries are two units. + +### What LLVM made of it + +| | AVX-512 (VNNI and v4 alike) | AVX2 | +|---|---|---| +| react, per block (4 tallies x 32 lanes) | 53 insns: 8 `vpternlogq`, 8 `vpaddw`, 4 `vpsrlw`, 4 `vpmullw`, 4 `vpmovzxbw`, 2 `vextracti64x4` | 179 insns, incl. 16 `vmovq`, 6 `vpblendd` | +| anticipate, per block | 61 insns: 12 `vpternlogq`, 4 `vpsubw`, otherwise as react | 107 insns: `vpand`/`vpandn`/`vpor`/`vpxor` | + +On AVX-512 LLVM folds each carry formula into ternary logic: two +`vpternlogq` per carry (react) or three (anticipate). That is the XOR/AND +carry written in Rust and lowered by the compiler, with no facade change. +There is no shuffle, permute, pack or float instruction. The `vpmovzxbw` and +`vextracti64x4` are the u8 -> u16 widening that `mullo` needs, the same as +in A'. + +### Measured (ns per 64 values, 3 runs, L1 / streamed) + +| | AVX-512 VNNI | AVX-512, no VNNI | AVX2 | +|---|---|---|---| +| A' arithmetic only (f32 accumulate) | 6.34 - 6.50 / 18.4 - 18.8 | 6.28 - 6.69 / 19.0 - 19.6 | 10.7 - 11.0 / 21.7 - 22.0 | +| **D react arithmetic only** | **5.20 - 5.70** / 18.2 - 20.1 | **5.36 - 5.64** / 18.4 - 19.8 | **21.2 - 24.2** / 31.9 - 32.0 | +| D anticipate arithmetic only | 5.94 - 6.45 / 18.7 - 20.9 | 5.98 - 6.31 / 19.4 - 20.5 | 9.89 - 10.1 / 22.0 - 26.9 | +| B arithmetic only (#388) | 2.39 - 2.41 / 17.3 - 17.9 | 5.02 - 5.03 / 22.3 - 22.8 | 8.95 - 9.17 / 27.2 - 29.4 | +| A' end-to-end | 7.62 - 7.83 / 18.0 - 18.1 | 7.83 / 19.5 - 19.7 | 81.1 - 88.7 / 416 - 421 | +| **D react end-to-end** | **6.62 - 7.25** / 11.0 - 15.1 | **6.77 - 6.80** / 15.3 - 16.0 | **79.6 - 90.1** / 434 - 440 | +| D anticipate end-to-end | 7.78 - 7.97 / 15.3 - 19.5 | 7.78 - 7.80 / 19.1 - 19.7 | 77.6 - 79.7 / 422 - 426 | +| B end-to-end (#388) | 2.89 - 2.91 / 6.06 - 7.71 | 5.62 / 11.0 - 12.9 | 10.0 - 10.1 / 26.5 - 28.2 | +| D react scan, query widened once | 6.22 - 6.56 / 11.9 - 16.2 | 6.23 - 7.20 / 11.7 - 19.1 | 88.6 - 91.1 / 433 - 436 | +| D anticipate scan | 7.00 - 7.47 / 11.0 - 15.9 | 7.27 - 8.41 / 10.4 - 22.7 | 76.9 - 77.2 / 417 - 420 | +| B scan (#388) | 2.66 / 2.88 - 3.54 | 4.67 - 4.78 / 6.13 - 9.15 | 9.17 - 9.83 / 26.4 - 27.1 | + +### Reading it + +- **The carry column beats the f32 detour in L1.** On AVX-512, react is + about 10 % faster than A' in arithmetic (5.2 - 5.7 vs 6.3 - 6.5 ns) and in + the full path (6.6 - 7.3 vs 7.6 - 7.8). Streamed, on compact operands, + their arithmetic-only rows are level (18 - 20 vs 18 - 19 ns: both read + precomputed operands and are memory-bound), and react end to end is + faster (11 - 15 vs 18 ns). Anticipate costs one more ternary op per carry + and lands level with A'. +- **React vs anticipate depends on the tier.** On AVX-512, react wins (the + carry needs the sum anyway, and two ternary ops suffice). On AVX2, + anticipate wins (9.9 - 10.1 vs 21 - 24 ns); the react form gets 16 `vmovq` + and 6 `vpblendd` from LLVM there. +- **B still wins, and the reason is not the accumulation.** D removes the + float conversion, but it still widens u8 -> u16 (`vpmovzxbw`, + `vextracti64x4`) and runs four `vpmullw` per block before any add. B + multiplies straight from bytes and accumulates in `i32` in the same + instruction. On VNNI, B's arithmetic is 2.4 ns against D's 5.2 - 5.7. +- **Sign-routing overhead** (end-to-end minus arithmetic-only medians, so it + also contains `check` and widening `|b|`): D react on VNNI 1.45 ns, on AVX-512 + without VNNI 1.33 ns, on AVX2 about 64 ns. A' carries the same routing + (1.47 ns on VNNI, 77 ns on AVX2). B routes with the byte-dot's own mask + instead (0.48 ns on VNNI, 0.88 on AVX2). +- **On AVX2 the end-to-end and scan rows are dominated by routing.** D's + routing is A's `U8x64::mask_blend`, whose AVX2 realization is a per-lane + loop (77 - 91 ns). That cost has nothing to do with the carry. + +Verdict for D: **ADOPT-LAB as the exact replacement for A's f32 detour**, and +the right way to accumulate any U16 product stream (the compiler gives +ternary logic for free on AVX-512). **DEFER as a B replacement**: it would +need the multiply to stay in bytes, which is what B already does. + +Correctness: the carry rules against `u32` arithmetic for every `low` in +0..=65535, with `p` at, just above and just below the headroom, plus 0, +65,025 and random products; all 511 x 511 pairs; full blocks; all-255 runs of +131,075 blocks past the flush bound, with `high` near its limit. Passed on all +seven tiers. Disable runs, each red: react without the `(low | p) & !sum` +term; anticipate without the borrow term; carry not added to `high`; flush +at 70,000 blocks (`high` wraps); `high` weighted 2^15. + +## The arithmetic + +For magnitudes `a = sum a_i 2^i` and `b = sum b_j 2^j`, +`a * b = sum over i, j of (a_i AND b_j) * 2^(i + j)`. A block of 64 lanes is +read as eight bit planes (`planes`): bit `k` of `plane[i]` is bit `i` of +lane `k`. One AND of two planes gives one partial-product bit for all 64 +lanes at significance `i + j`. The sign is `negative = plane_a XOR plane_b`; +it only selects the positive or the negative tally. No magnitude is negated. + +`carry_save3(a, b, c)` returns `(sum, carry)` with, for every bit position +`k`, `a_k + b_k + c_k = sum_k + 2 * carry_k`. `sum` stays in its column, +`carry` belongs to the next; nothing is shifted across bit positions. + +Three folds of the 128 partial-product planes per block (64 per sign): + +| | fold | normalization | +|---|---|---| +| C popcnt | `popcount(plane) << (i + j)` into an `i64`, at once | every block (it is the fold) | +| C abacus | `Abacus`: at most two planes per column; a third compresses via `carry_save3`, the carry rippling up | once, at `finish` | +| C csa | one running plane per column and sign, folded on a fixed, data-independent schedule; carries leaving column 15 counted into `spill` | once, at `finish` | + +## Correctness + +`tests/u9_abacus.rs`, against #388's `i64` oracle and its B kernel: + +- `carry_save3` bit by bit for all eight input combinations at every bit + position, plus 10,000 random word triples (popcount identity); +- `planes` recovers every magnitude 0..=255 in every lane; +- carries through many columns: 2^20 all-ones planes in column 0 + (64 * 2^20 exactly), one all-ones plane per column 0..=39 three times, and + 5,000 mixed-pattern planes; +- all 511 x 511 pairs through C popcnt, C abacus, C csa, the C scan and B; +- full mixed blocks (1 to 300 blocks); +- all-255 runs of 49,157 blocks with all sign mixes, totals beyond `i32`; +- never-touched refused by `check` in either operand; touched zero scores 0. + +All five folds equal the oracle in every case. 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: `carry_save3` without `t & c`; the +abacus dropping its carry; carry kept in the same column; a plane extracted +from the wrong bit; the abacus sign routing swapped; the popcount weight +`i * j` instead of `i + j`; `spill` weighted `2^15`; the static schedule +dropping an odd plane; carries not forwarded to the next column. + +Exactness domain: the 511 numeric values, any number of blocks while the +`i64` total stays in range (about 2.2 * 10^12 worst-case blocks). The +abacus's column bound is by value: column `w` can only hold a set bit for a +lane whose total is at least `2^w`, so 64 columns hold any per-lane total +below 2^64; `Abacus::add` panics rather than drop a carry past them. + +## What LLVM made of it (`objdump` of `#[inline(never)]` copies) + +| function | AVX-512 VNNI | AVX-512, no VNNI | AVX2 | +|---|---|---|---| +| `carry_save3` | 3 `xor`, 2 `and`, 2 `mov`: scalar, no ternary logic | same | same | +| `planes` | 6 `vptestmb` + 2 `vpmovb2m` (+ `vpsllw` for bit 0): one instruction per plane | same | out of line (a `jmp`); extraction 11.1 ns vs 3.5 | +| C popcnt `mac` (140 / 196 / 216 insns) | vectorized: 32 `vpopcntq`, `vpand`, `vpsllvq` | **32 `vpshufb`** nibble-table popcount + `vpsadbw` | **32 `vpshufb`** + `vpsadbw` | +| C abacus `mac` (110 insns) | branchy loop over the column counters, columns in memory; scalar | same | same | +| C csa `mac` (437-488 insns) | not unrolled: loops over stack arrays, 122 `mov`, 37 `cmp`; scalar logic, no `vpternlog` | same | same | + +Three findings: + +1. LLVM lowers the plane extraction as well as possible and vectorizes the + per-plane popcounts. It does not turn either carry-save fold into + straight-line or vector code, even though the static schedule depends + only on loop indices. Writing it out by hand would be manual scheduling, + which this probe was told not to do. +2. **Without VPOPCNTDQ the popcount fold is a byte shuffle.** On v4 and AVX2, + `u64::count_ones` lowers to the `vpshufb` nibble table, the shape the + repository's no-shuffle rule forbids in hand-written code. Only the + AVX-512 VNNI host (which has VPOPCNTDQ) avoids it. B uses no shuffle on + any tier (#388). +3. No variant needed a facade word the facade lacks: `planes` is `BitAnd` + + `cmpeq_mask`, and LLVM fused them into `vptestmb`. + +## Measurements (ns per 64 values, medians of 3 runs, this host) + +One AVX-512 VNNI Xeon; each tier a separate native build (`config-v4`, +`config-v3`). L1 = 64 blocks per operand; streamed = 2^20 blocks (72 MiB). +Every kernel equals the oracle on the benchmark data before it is timed. + +Timed operands use #388's compact layout (64 magnitude bytes and one plane +per block in two arrays, 72 bytes); each `Block` is built in registers. The +first version streamed a `Vec`, which pads every block to 128 bytes, +and its extraction row consumed four of the sixteen planes, so the optimizer +could drop the rest (it read 2.52 ns on VNNI L1). Both were fixed after +review (Codex, PR #389) and every table here was re-measured; the plane +extraction row now passes both complete plane arrays through `black_box`. + +| L1-resident | AVX-512 VNNI | AVX-512, no VNNI | AVX2 | +|---|---|---|---| +| B end-to-end (#388) | 2.89 - 2.91 | 5.62 | 10.0 - 10.1 | +| C popcnt end-to-end | 28.1 - 28.2 | 47.1 - 48.2 | 42.4 | +| C abacus end-to-end | 286 - 288 | 290 - 292 | 317 - 319 | +| C csa end-to-end | 298 | 295 - 296 | 304 - 305 | +| B arithmetic only | 2.39 - 2.41 | 5.02 - 5.03 | 8.95 - 9.17 | +| C popcnt arithmetic only, planes given | 12.4 - 12.8 | 30.1 | 30.7 | +| C abacus arithmetic only, deferred | 266 - 276 | 289 - 292 | 292 - 298 | +| C abacus arithmetic only, normalized every block | 271 - 278 | 276 - 282 | 281 - 284 | +| C csa arithmetic only | 293 | 291 | 290 - 291 | +| C plane extraction, both operands (cold) | 3.45 - 3.53 | 3.45 - 3.58 | 11.1 | +| B scan (#388) | 2.66 | 4.67 - 4.78 | 9.17 - 9.83 | +| C popcnt scan, query planes once | 23.4 - 23.9 | 46.7 - 47.8 | 41.8 - 42.9 | +| C abacus scan | 286 - 291 | 296 - 298 | 313 - 329 | +| C csa scan | 301 - 307 | 300 - 301 | 302 - 303 | + +| Streamed | AVX-512 VNNI | AVX-512, no VNNI | AVX2 | +|---|---|---|---| +| B end-to-end | 6.06 - 7.71 | 11.0 - 12.9 | 26.5 - 28.2 | +| C popcnt end-to-end | 51.1 - 55.9 | 50.0 - 50.7 | 54.5 - 57.2 | +| C popcnt arithmetic only | 47.0 - 50.0 | 60.1 - 60.7 | 57.4 - 60.7 | +| C abacus, deferred / normalized every block | 438 - 446 / 283 - 285 | 404 - 405 / 293 - 309 | 398 - 408 / 292 - 300 | +| C plane extraction, both operands | 7.81 - 11.4 | 11.2 - 11.6 | 56.1 - 57.6 | +| B scan | 2.88 - 3.54 | 6.13 - 9.15 | 26.4 - 27.1 | +| C popcnt scan | 54.3 - 57.2 | 61.8 - 64.3 | 60.1 - 62.6 | + +Streamed "arithmetic only" rows read precomputed planes (another 64 bytes +per block per operand), so they measure more memory, not less arithmetic. + +### Break-even for prepared operands + +None, on any tier. Prepared reuse can only remove plane extraction +(3.5 ns per block for both operands on VNNI), and C popcnt's arithmetic +alone, with both operands' planes free, is 12.4 ns against B's 2.9 ns with +everything included. In the scan, where the query's planes are prepared once, +C popcnt costs 23.4 - 23.9 ns against B's 2.66 ns. + +## Why the gap is structural, not tuning + +B's byte-dot multiplies 64 byte pairs in four instructions per block. The +bit-serial form needs 64 partial-product planes per sign, one AND each, +before any reduction: 128 ANDs per block. Even packing eight blocks into +one 512-bit register, with the reduction perfectly vectorized (ternary logic +for `carry_save3`), it would take about 16 ANDs plus at least as many +reduction operations per block. That is an estimate, not a measurement, but +it is already several times B's instruction count. Carry-save trades +popcounts for logic operations; it cannot remove the 64 partial products +that the hardware multiplier forms in one step. + +## Falsification criteria (from the brief) + +| criterion | result | +|---|---| +| bit-plane preparation dominates the saved arithmetic | no saving to dominate: C's arithmetic alone exceeds B's total | +| carry compression generates more instructions than the byte-dot | yes: 437-488 instructions per block for C csa, against 4 masked byte-dots for B | +| register pressure or spills negate the gain | yes: the carry-save folds keep their columns in memory (`mov` dominated) | +| deferred normalization costs more than immediate | yes when streamed (438 - 446 vs 283 - 285 ns), equal in L1 | +| exactness or provenance violated | no: every fold is exact on every tier | +| competitiveness needs backend-specific code | yes, it would: a hand-unrolled ternary-logic schedule; not built | + +## Not done + +- No hand-unrolled or vectorized-across-blocks carry-save schedule (that is + manual instruction scheduling). +- No NEON or wasm timings (qemu timings mean nothing). diff --git a/Cargo.toml b/Cargo.toml index e7e552f6..aec13c9c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -71,6 +71,11 @@ name = "u9_head_to_head" path = "examples/u9_head_to_head/main.rs" required-features = ["std"] +[[example]] +name = "u9_abacus" +path = "examples/u9_abacus/main.rs" +required-features = ["std"] + [[example]] name = "cpu_guard_probe" required-features = ["std"] diff --git a/examples/u9_abacus/abacus.rs b/examples/u9_abacus/abacus.rs new file mode 100644 index 00000000..f10dbb2e --- /dev/null +++ b/examples/u9_abacus/abacus.rs @@ -0,0 +1,276 @@ +//! D-U9-ABACUS: exact U9 x U9 multiply-accumulate as schoolbook arithmetic on +//! bit planes, written in plain Rust bit operations. Probe-private; nothing +//! here is a facade word or a storage format. +//! +//! A block of 64 magnitudes is read as eight bit planes: bit `k` of +//! `plane[i]` is bit `i` of lane `k`. The product of two lanes is +//! `sum over i, j of (a_i AND b_j) * 2^(i + j)`, so one AND of two planes +//! yields one partial-product bit for all 64 lanes at significance `i + j`. +//! The sign is the XOR of the two polarity planes and only selects whether a +//! partial product joins the positive or the negative tally; no magnitude is +//! negated. +//! +//! Two ways to fold the partial products: +//! +//! * [`AccPopcnt`] counts each partial-product plane at once +//! (`popcount * 2^(i + j)`) into an `i64`; +//! * [`AccAbacus`] keeps them uncounted in carry-save columns +//! ([`Abacus`]): three planes in column `w` compress to one in `w` and one +//! carry in `w + 1` ([`carry_save3`]); counting happens once, at `finish`. +//! +//! Extracting the planes from the existing `U8x64` magnitudes uses facade +//! words only (`BitAnd`, `cmpeq_mask`) and is timed as part of the cost. + +#![allow(dead_code)] + +use ndarray::simd::U8x64; + +/// Full adder on 64 independent bit columns. For every bit position `k`, +/// `a_k + b_k + c_k = sum_k + 2 * carry_k`: `sum` stays in the current +/// significance, `carry` belongs to the next one. Nothing is shifted across +/// bit positions; each bit is one lane. +#[inline] +pub fn carry_save3(a: u64, b: u64, c: u64) -> (u64, u64) { + let t = a ^ b; + let sum = t ^ c; + let carry = (a & b) ^ (t & c); + (sum, carry) +} + +/// The eight bit planes of 64 magnitudes: bit `k` of `plane[i]` is bit `i` of +/// lane `k`. Facade words only: one AND and one byte compare per plane. +#[inline(always)] +pub fn planes(mag: U8x64) -> [u64; 8] { + core::array::from_fn(|i| { + let bit = U8x64::splat(1 << i); + (mag & bit).cmpeq_mask(bit) + }) +} + +/// Significance columns an [`Abacus`] holds. Column `w` can only hold a set +/// bit for a lane whose running total is at least `2^w`, so 64 columns cover +/// any per-lane total below `2^64`; [`Abacus::add`] panics rather than drop a +/// carry past the last column. +pub const COLUMNS: usize = 64; + +/// Uncounted carry-save tally of 64 independent lanes: column `w` holds at +/// most two planes of significance `2^w`. The represented per-lane value is +/// `sum over w of 2^w * (bits of that lane in column w)`. +#[derive(Clone)] +pub struct Abacus { + col: [[u64; 2]; COLUMNS], + n: [u8; COLUMNS], +} + +impl Abacus { + pub fn new() -> Self { + Self { + col: [[0; 2]; COLUMNS], + n: [0; COLUMNS], + } + } + + /// Adds plane `p` at significance `2^w`. A third plane in a column is + /// compressed with the two already there: the sum stays, the carry moves + /// up one column and is added there the same way. + #[inline(always)] + pub fn add(&mut self, mut w: usize, mut p: u64) { + loop { + let k = usize::from(self.n[w]); + if k < 2 { + self.col[w][k] = p; + self.n[w] += 1; + return; + } + let (s, c) = carry_save3(self.col[w][0], self.col[w][1], p); + self.col[w][0] = s; + self.n[w] = 1; + p = c; + w += 1; + assert!(w < COLUMNS, "abacus carry past column {COLUMNS}"); + } + } + + /// Sum over all 64 lanes, counted once: `popcount(plane) * 2^w`, exact in + /// `i128` for every abacus (at most 2 * 64 planes, 64 bits each, below 2^64). + pub fn value(&self) -> i128 { + let mut v = 0i128; + for w in 0..COLUMNS { + for k in 0..usize::from(self.n[w]) { + v += i128::from(self.col[w][k].count_ones()) << w; + } + } + v + } +} + +impl Default for Abacus { + fn default() -> Self { + Self::new() + } +} + +/// C, deferred: positive and negative partial products tallied in two +/// abaci, counted only at [`Self::finish`]. +pub struct AccAbacus { + pos: Abacus, + neg: Abacus, +} + +impl AccAbacus { + pub fn new() -> Self { + Self { + pos: Abacus::new(), + neg: Abacus::new(), + } + } + + /// One block from prepared planes. `negative` is `plane_a ^ plane_b`. + #[inline(always)] + pub fn mac(&mut self, a: &[u64; 8], negative: u64, b: &[u64; 8]) { + for (i, &ai) in a.iter().enumerate() { + let (ap, an) = (ai & !negative, ai & negative); + for (j, &bj) in b.iter().enumerate() { + self.pos.add(i + j, ap & bj); + self.neg.add(i + j, an & bj); + } + } + } + + pub fn finish(&self) -> i64 { + i64::try_from(self.pos.value() - self.neg.value()).expect("total outside i64") + } +} + +impl Default for AccAbacus { + fn default() -> Self { + Self::new() + } +} + +/// C, static carry-save: one running plane per column and sign, folded +/// with [`carry_save3`] on a fixed schedule that depends only on the column, +/// never on the data. Per block and sign, column `w` takes its partial +/// products and the carries out of column `w - 1`, two at a time, into its +/// running plane; carries leaving column 15 are counted at once into `spill`. +/// Represented value per sign: `sum over w of popcount(ones[w]) * 2^w + spill`. +pub struct AccCsa { + ones: [[u64; 16]; 2], + spill: [i64; 2], +} + +impl AccCsa { + pub fn new() -> Self { + Self { + ones: [[0; 16]; 2], + spill: [0; 2], + } + } + + #[inline(always)] + pub fn mac(&mut self, a: &[u64; 8], negative: u64, b: &[u64; 8]) { + for (s, mask) in [!negative, negative].into_iter().enumerate() { + let am: [u64; 8] = core::array::from_fn(|i| a[i] & mask); + let mut cin = [0u64; 8]; + let mut nc = 0usize; + for w in 0usize..16 { + // At most 8 partial products and 8 carries enter a column. + let mut items = [0u64; 16]; + let mut k = 0usize; + for i in w.saturating_sub(7)..=w.min(7) { + items[k] = am[i] & b[w - i]; + k += 1; + } + for &c in &cin[..nc] { + items[k] = c; + k += 1; + } + let mut o = self.ones[s][w]; + let mut cout = [0u64; 8]; + let mut no = 0usize; + let mut t = 0usize; + while t < k { + let q = if t + 1 < k { items[t + 1] } else { 0 }; + let (sum, carry) = carry_save3(o, items[t], q); + o = sum; + cout[no] = carry; + no += 1; + t += 2; + } + self.ones[s][w] = o; + cin = cout; + nc = no; + } + for &c in &cin[..nc] { + self.spill[s] += i64::from(c.count_ones()) << 16; + } + } + } + + pub fn finish(&self) -> i64 { + let side = |s: usize| -> i64 { + let cols: i64 = (0..16) + .map(|w| i64::from(self.ones[s][w].count_ones()) << w) + .sum(); + cols + self.spill[s] + }; + side(0) - side(1) + } +} + +impl Default for AccCsa { + fn default() -> Self { + Self::new() + } +} + +/// C, immediate: each partial-product plane counted at once into `i64`. +/// Per block the change is at most 64 * 65_025; the `i64` total is exact +/// while it stays in range (about 2.2 * 10^12 worst-case blocks). +pub struct AccPopcnt { + total: i64, +} + +impl AccPopcnt { + pub fn new() -> Self { + Self { total: 0 } + } + + #[inline(always)] + pub fn mac(&mut self, a: &[u64; 8], negative: u64, b: &[u64; 8]) { + let mut block = 0i64; + for (i, &ai) in a.iter().enumerate() { + let (ap, an) = (ai & !negative, ai & negative); + for (j, &bj) in b.iter().enumerate() { + let d = i64::from((ap & bj).count_ones()) - i64::from((an & bj).count_ones()); + block += d << (i + j); + } + } + self.total += block; + } + + pub fn finish(&self) -> i64 { + self.total + } +} + +impl Default for AccPopcnt { + fn default() -> Self { + Self::new() + } +} + +/// A query block prepared once for a scan: its planes and polarity. +pub struct QueryC { + pub planes: [u64; 8], + pub positive: u64, +} + +impl QueryC { + pub fn new(mag: U8x64, positive: u64) -> Self { + Self { + planes: planes(mag), + positive, + } + } +} diff --git a/examples/u9_abacus/anchor.rs b/examples/u9_abacus/anchor.rs new file mode 100644 index 00000000..895eef45 --- /dev/null +++ b/examples/u9_abacus/anchor.rs @@ -0,0 +1,230 @@ +//! D-U9-ANCHOR: candidate E. B's split of `|b|` at the power of two 128, +//! `|b| = low7 + 128 * bit7`, with the anchor term taken as a shift instead +//! of a multiplication: +//! +//! `a * |b| = a * low7 + (sum of a over the lanes where bit 7 of b is set) << 7` +//! +//! `low7` still needs a real byte multiply; the anchor term needs only a masked +//! sum of `a`, weighted by 128 once, at flush. B spends a full byte-dot against +//! `-128` on it. Three ways to form the masked sum from existing facade words: +//! +//! * E4 ([`AccE4`]): byte-dot of `a` against all-ones bytes, masked by +//! `sign mask AND bit-7 plane`, one per sign: still four byte-dots, but the +//! second operand is a constant; +//! * E3 ([`AccE3`]): one byte-dot of `a` against a `+1 / -1` byte vector +//! (built from the sign plane), masked by the bit-7 plane: three byte-dots; +//! * Esum ([`AccESum`]): `U8x64::mask_blend` + `sum_bytes_u64` per sign, which +//! reduces horizontally on every call. +//! +//! The polarity rule is unchanged: `negative = plane_a XOR plane_b`. + +#![allow(dead_code)] + +use crate::kernels::{bits_as_i8, lanes_i64, Block, B_FLUSH_BLOCKS}; +use ndarray::simd::{I32x16, I8x64, U8x64}; + +/// Lanes whose magnitude has bit 7 set, as a plane. +#[inline(always)] +pub fn bit7(mag: U8x64) -> u64 { + let top = U8x64::splat(0x80); + (mag & top).cmpeq_mask(top) +} + +/// The prepared `b` side: `low7` as signed bytes and the bit-7 plane. +#[inline(always)] +pub fn split_anchor(b: U8x64) -> (I8x64, u64) { + (bits_as_i8(b & U8x64::splat(0x7F)), bit7(b)) +} + +/// E4: four byte-dots, the two anchor dots against constant all-ones bytes. +/// Lanes: `[pos x low7, neg x low7, pos anchor, neg anchor]`. Per block a lane +/// moves by at most `4 * 255 * 127` (low7) or `4 * 255` (anchor), so +/// `B_FLUSH_BLOCKS` keeps every lane inside `i32`. +pub struct AccE4 { + acc: [I32x16; 4], + blocks: usize, + total: i64, +} + +impl AccE4 { + pub fn new() -> Self { + Self { + acc: [I32x16::splat(0); 4], + blocks: 0, + total: 0, + } + } + + pub fn flush(&mut self) { + let r = self.acc.map(lanes_i64); + self.total += (r[0] - r[1]) + ((r[2] - r[3]) << 7); + self.acc = [I32x16::splat(0); 4]; + self.blocks = 0; + } + + #[inline(always)] + pub fn mac(&mut self, a: U8x64, negative: u64, low7: I8x64, top: u64) { + let ones = I8x64::splat(1); + let (p, n) = (!negative, negative); + self.acc[0] = self.acc[0].dot4_u8_i8_add_masked(a, p, low7); + self.acc[1] = self.acc[1].dot4_u8_i8_add_masked(a, n, low7); + self.acc[2] = self.acc[2].dot4_u8_i8_add_masked(a, p & top, ones); + self.acc[3] = self.acc[3].dot4_u8_i8_add_masked(a, n & top, ones); + self.blocks += 1; + if self.blocks == B_FLUSH_BLOCKS { + self.flush(); + } + } + + #[inline(always)] + pub fn step(&mut self, a: Block, b: Block) { + let (low7, top) = split_anchor(b.mag); + self.mac(a.mag, a.positive ^ b.positive, low7, top); + } + + #[inline(always)] + pub fn scan_step(&mut self, x: Block, q: &QueryE) { + self.mac(x.mag, x.positive ^ q.positive, q.low7, q.top); + } + + pub fn finish(mut self) -> i64 { + self.flush(); + self.total + } +} + +/// `+1` where the product is positive, `-1` where it is negative, as signed +/// bytes. Built by `U8x64::mask_blend` from the sign plane; a sign, not a +/// negated magnitude. +#[inline(always)] +pub fn sign_bytes(negative: u64) -> I8x64 { + bits_as_i8(U8x64::mask_blend(negative, U8x64::splat(1), U8x64::splat(0xFF))) +} + +/// E3: three byte-dots, the anchor in one signed dot. +/// Lanes: `[pos x low7, neg x low7, signed anchor]`. +pub struct AccE3 { + acc: [I32x16; 3], + blocks: usize, + total: i64, +} + +impl AccE3 { + pub fn new() -> Self { + Self { + acc: [I32x16::splat(0); 3], + blocks: 0, + total: 0, + } + } + + pub fn flush(&mut self) { + let r = self.acc.map(lanes_i64); + self.total += (r[0] - r[1]) + (r[2] << 7); + self.acc = [I32x16::splat(0); 3]; + self.blocks = 0; + } + + #[inline(always)] + pub fn mac(&mut self, a: U8x64, negative: u64, low7: I8x64, top: u64) { + self.acc[0] = self.acc[0].dot4_u8_i8_add_masked(a, !negative, low7); + self.acc[1] = self.acc[1].dot4_u8_i8_add_masked(a, negative, low7); + self.acc[2] = self.acc[2].dot4_u8_i8_add_masked(a, top, sign_bytes(negative)); + self.blocks += 1; + if self.blocks == B_FLUSH_BLOCKS { + self.flush(); + } + } + + #[inline(always)] + pub fn step(&mut self, a: Block, b: Block) { + let (low7, top) = split_anchor(b.mag); + self.mac(a.mag, a.positive ^ b.positive, low7, top); + } + + #[inline(always)] + pub fn scan_step(&mut self, x: Block, q: &QueryE) { + self.mac(x.mag, x.positive ^ q.positive, q.low7, q.top); + } + + pub fn finish(mut self) -> i64 { + self.flush(); + self.total + } +} + +/// Esum: the two `low7` byte-dots, and the anchor sums as masked byte sums +/// reduced to scalars on every block (`sum_bytes_u64`). +pub struct AccESum { + acc: [I32x16; 2], + anchor: i64, + blocks: usize, + total: i64, +} + +impl AccESum { + pub fn new() -> Self { + Self { + acc: [I32x16::splat(0); 2], + anchor: 0, + blocks: 0, + total: 0, + } + } + + pub fn flush(&mut self) { + let r = self.acc.map(lanes_i64); + self.total += (r[0] - r[1]) + (self.anchor << 7); + self.acc = [I32x16::splat(0); 2]; + self.anchor = 0; + self.blocks = 0; + } + + #[inline(always)] + pub fn mac(&mut self, a: U8x64, negative: u64, low7: I8x64, top: u64) { + let zero = U8x64::splat(0); + self.acc[0] = self.acc[0].dot4_u8_i8_add_masked(a, !negative, low7); + self.acc[1] = self.acc[1].dot4_u8_i8_add_masked(a, negative, low7); + let pos = U8x64::mask_blend(!negative & top, zero, a).sum_bytes_u64(); + let neg = U8x64::mask_blend(negative & top, zero, a).sum_bytes_u64(); + self.anchor += pos as i64 - neg as i64; + self.blocks += 1; + if self.blocks == B_FLUSH_BLOCKS { + self.flush(); + } + } + + #[inline(always)] + pub fn step(&mut self, a: Block, b: Block) { + let (low7, top) = split_anchor(b.mag); + self.mac(a.mag, a.positive ^ b.positive, low7, top); + } + + #[inline(always)] + pub fn scan_step(&mut self, x: Block, q: &QueryE) { + self.mac(x.mag, x.positive ^ q.positive, q.low7, q.top); + } + + pub fn finish(mut self) -> i64 { + self.flush(); + self.total + } +} + +/// A query prepared once: `low7` bytes, bit-7 plane, polarity. +pub struct QueryE { + pub low7: I8x64, + pub top: u64, + pub positive: u64, +} + +impl QueryE { + pub fn new(q: Block) -> Self { + let (low7, top) = split_anchor(q.mag); + Self { + low7, + top, + positive: q.positive, + } + } +} diff --git a/examples/u9_abacus/carry16.rs b/examples/u9_abacus/carry16.rs new file mode 100644 index 00000000..32c2d1e8 --- /dev/null +++ b/examples/u9_abacus/carry16.rs @@ -0,0 +1,169 @@ +//! D-U9-CARRY16: candidate A' (`U16x32::mullo` products, sign routed by the +//! polarity XOR) with its `f32` detour replaced by a 16-bit place-value +//! column. Each tally lane keeps `low` (wrapping `u16`) and `high` (count of +//! 2^16 units): the represented value is `low + 2^16 * high`. A product +//! `p < 2^16` is added to `low` with a wrapping add; when the add crosses +//! 2^16, one unit moves into `high`. Nothing is widened before it must be. +//! +//! Two ways to find that carry, both plain bit logic on `U16x32` lanes +//! (no masks, no compares into mask registers, no shuffles): +//! +//! * react ([`carry_after`]): add first, then read the carry out of the +//! finished sum, `msb((low & p) | ((low | p) & !sum))`; +//! * anticipate ([`carry_before`]): before adding, the carry is set exactly +//! when `p` exceeds the headroom `!low = 65535 - low`, which is the borrow +//! out of `!low - p`: `msb((low & p) | (!(!low ^ p) & (!low - p)))`. + +#![allow(dead_code)] + +use crate::kernels::route_a; +use crate::kernels::Block; +use ndarray::simd::{U16x32, U8x64}; + +/// Carry out of the wrapping add `low + p = sum`, as 0 or 1 per lane. +#[inline(always)] +pub fn carry_after(low: U16x32, p: U16x32, sum: U16x32) -> U16x32 { + ((low & p) | ((low | p) & !sum)).shr(15) +} + +/// Carry out of `low + p`, as 0 or 1 per lane, computed from the operands +/// before the add: 1 exactly when `p > !low`, the borrow out of `!low - p`. +#[inline(always)] +pub fn carry_before(low: U16x32, p: U16x32) -> U16x32 { + let head = !low; + let d = head - p; + ((low & p) | (!(head ^ p) & d)).shr(15) +} + +/// One place-value column pair: `low + 2^16 * high` per lane. +#[derive(Clone, Copy)] +pub struct Col16 { + pub low: U16x32, + pub high: U16x32, +} + +impl Col16 { + pub fn zero() -> Self { + Self { + low: U16x32::splat(0), + high: U16x32::splat(0), + } + } + + /// Adds products `p` (each below 2^16). `ANTICIPATE` picks the carry rule. + #[inline(always)] + pub fn add(&mut self, p: U16x32) { + if ANTICIPATE { + let c = carry_before(self.low, p); + self.low += p; + self.high += c; + } else { + let sum = self.low + p; + self.high += carry_after(self.low, p, sum); + self.low = sum; + } + } + + /// Sum of all 32 lanes in `i64`, read lane by lane from a stored array. + pub fn total(&self) -> i64 { + let (lo, hi) = (core::hint::black_box(self.low.to_array()), core::hint::black_box(self.high.to_array())); + let mut t = 0i64; + for (&l, &h) in lo.iter().zip(&hi) { + t += i64::from(l) + (i64::from(h) << 16); + } + t + } +} + +/// Blocks before `flush`: each block adds one product per lane and tally, so +/// `high` grows by at most 1 per block and stays below 2^16 for 65_535 blocks. +pub const D_FLUSH_BLOCKS: usize = 65_535; + +/// Candidate D: positive and negative tallies, each for the low and the high +/// 32 lanes of a block, folded into `i64` at most every `D_FLUSH_BLOCKS`. +pub struct AccD { + col: [[Col16; 2]; 2], + blocks: usize, + total: i64, +} + +impl AccD { + pub fn new() -> Self { + Self { + col: [[Col16::zero(); 2]; 2], + blocks: 0, + total: 0, + } + } + + pub fn flush(&mut self) { + let side = |c: &[Col16; 2]| c[0].total() + c[1].total(); + self.total += side(&self.col[0]) - side(&self.col[1]); + self.col = [[Col16::zero(); 2]; 2]; + self.blocks = 0; + } + + /// The arithmetic alone, on magnitudes already routed by sign. + #[inline(always)] + pub fn mac(&mut self, a_pos: U8x64, a_neg: U8x64, b_lo: U16x32, b_hi: U16x32) { + let products = [ + [U16x32::from_u8x64_lo(a_pos).mullo(b_lo), U16x32::from_u8x64_hi(a_pos).mullo(b_hi)], + [U16x32::from_u8x64_lo(a_neg).mullo(b_lo), U16x32::from_u8x64_hi(a_neg).mullo(b_hi)], + ]; + for (cols, ps) in self.col.iter_mut().zip(products) { + for (c, p) in cols.iter_mut().zip(ps) { + c.add::(p); + } + } + self.blocks += 1; + if self.blocks == D_FLUSH_BLOCKS { + self.flush(); + } + } + + /// One block from the shared representation: route `|a|` as A does, + /// widen `|b|`, multiply, add with carry. + #[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, U16x32::from_u8x64_lo(b.mag), U16x32::from_u8x64_hi(b.mag)); + } + + /// One index block against a query whose magnitudes were widened once. + #[inline(always)] + pub fn scan_step(&mut self, x: Block, q: &QueryD) { + let negative = x.positive ^ q.positive; + let zero = U8x64::splat(0); + let x_pos = U8x64::mask_blend(negative, x.mag, zero); + let x_neg = U8x64::mask_blend(!negative, x.mag, zero); + self.mac(x_pos, x_neg, q.lo, q.hi); + } + + pub fn finish(mut self) -> i64 { + self.flush(); + self.total + } +} + +impl Default for AccD { + fn default() -> Self { + Self::new() + } +} + +/// A query block prepared once for a scan: magnitudes widened to `u16`. +pub struct QueryD { + pub lo: U16x32, + pub hi: U16x32, + pub positive: u64, +} + +impl QueryD { + pub fn new(q: Block) -> Self { + Self { + lo: U16x32::from_u8x64_lo(q.mag), + hi: U16x32::from_u8x64_hi(q.mag), + positive: q.positive, + } + } +} diff --git a/examples/u9_abacus/main.rs b/examples/u9_abacus/main.rs new file mode 100644 index 00000000..2050d129 --- /dev/null +++ b/examples/u9_abacus/main.rs @@ -0,0 +1,528 @@ +//! D-U9-ABACUS probe: schoolbook U9 x U9 on bit planes (`abacus.rs`) against +//! the U9 head-to-head's B kernel (`../u9_head_to_head/kernels.rs`) on the +//! same data. Every kernel is checked against the `i64` oracle before it is +//! timed. `cargo run --release --example u9_abacus --features std`; the tier +//! is whatever this build compiled (printed in the header). + +#[path = "../u9_head_to_head/kernels.rs"] +mod kernels; +mod abacus; +mod anchor; +mod carry16; + +use abacus::*; +use anchor::*; +use carry16::*; +use kernels::*; +use ndarray::simd::{I32x16, I8x64, U16x32, 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_arch = "aarch64") { + "aarch64 NEON" + } else { + "other" + } +} + +type Raw = ([u8; 64], u64); + +/// Uniform over the 511 numeric values; never-touched never appears. +fn data(blocks: usize, seed: u64) -> (Vec, Vec) { + let mut s = seed; + let mut mk = || { + let mut m = [0u8; 64]; + let mut p = 0u64; + for (i, mi) in m.iter_mut().enumerate() { + s = s + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + let v = ((s >> 33) % 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(); + (a, b) +} + +fn time(name: &str, blocks: usize, runs: usize, base: 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 med = t[runs / 2]; + let rel = base + .map(|b| format!(" {:6.2}x B", med / b)) + .unwrap_or_default(); + println!("{name:<60} {med:8.3} ns/64{rel}"); + med +} + +fn bench(label: &str, n: usize, runs: usize, seed: u64) { + let (ra, rb) = data(n, seed); + let (ab, bb) = (Operand::new(&ra), Operand::new(&rb)); + let want = oracle(&ra, &rb).unwrap(); + let qraw = vec![rb[0]; n]; + let want_scan = oracle(&ra, &qraw).unwrap(); + let q = bb.block(0); + + // Execution-form operands for the arithmetic-only rows, built once. + let split: Vec<(u64, I8x64, I8x64)> = ab + .blocks() + .zip(bb.blocks()) + .map(|(a, b)| prepare_b(a, b)) + .collect(); + let pa: Vec<[u64; 8]> = ab.blocks().map(|x| planes(x.mag)).collect(); + let pb: Vec<[u64; 8]> = bb.blocks().map(|x| planes(x.mag)).collect(); + let neg: Vec = ab + .blocks() + .zip(bb.blocks()) + .map(|(a, b)| a.positive ^ b.positive) + .collect(); + + let b_e2e = || { + 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 b_ari = || { + let mut acc = AccB::new(); + for (a, (ng, l, t)) in ab.mag.iter().zip(&split) { + acc.mac(*a, *ng, *l, *t); + } + acc.finish() + }; + let b_scan = || { + 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 base = || { + let w0 = I8x64::from_array(rb[0].0.map(|x| i8::from_ne_bytes([x]))); + 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 c_e2e_abacus = || { + let mut acc = AccAbacus::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.mac(&planes(a.mag), a.positive ^ b.positive, &planes(b.mag)); + } + acc.finish() + }; + let c_e2e_pop = || { + let mut acc = AccPopcnt::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.mac(&planes(a.mag), a.positive ^ b.positive, &planes(b.mag)); + } + acc.finish() + }; + let c_ari_abacus = || { + let mut acc = AccAbacus::new(); + for ((x, y), ng) in pa.iter().zip(&pb).zip(&neg) { + acc.mac(x, *ng, y); + } + acc.finish() + }; + let c_ari_abacus_norm_each = || { + // Deferred carries normalized after every block instead of once. + let mut total = 0i64; + for ((x, y), ng) in pa.iter().zip(&pb).zip(&neg) { + let mut acc = AccAbacus::new(); + acc.mac(x, *ng, y); + total += acc.finish(); + } + total + }; + let c_ari_csa = || { + let mut acc = AccCsa::new(); + for ((x, y), ng) in pa.iter().zip(&pb).zip(&neg) { + acc.mac(x, *ng, y); + } + acc.finish() + }; + let c_e2e_csa = || { + let mut acc = AccCsa::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.mac(&planes(a.mag), a.positive ^ b.positive, &planes(b.mag)); + } + acc.finish() + }; + let c_scan_csa = || { + let qc = QueryC::new(q.mag, q.positive); + let mut acc = AccCsa::new(); + for x in ab.blocks() { + check(x, q).unwrap(); + acc.mac(&planes(x.mag), x.positive ^ qc.positive, &qc.planes); + } + acc.finish() + }; + let c_ari_pop = || { + let mut acc = AccPopcnt::new(); + for ((x, y), ng) in pa.iter().zip(&pb).zip(&neg) { + acc.mac(x, *ng, y); + } + acc.finish() + }; + let c_prep = || { + // Both complete plane arrays go through the barrier, so none of the + // sixteen extractions can be dropped as unused. + for (a, b) in ab.mag.iter().zip(&bb.mag) { + black_box((planes(*a), planes(*b))); + } + 0 + }; + let c_scan_abacus = || { + let qc = QueryC::new(q.mag, q.positive); + let mut acc = AccAbacus::new(); + for x in ab.blocks() { + check(x, q).unwrap(); + acc.mac(&planes(x.mag), x.positive ^ qc.positive, &qc.planes); + } + acc.finish() + }; + let c_scan_pop = || { + let qc = QueryC::new(q.mag, q.positive); + let mut acc = AccPopcnt::new(); + for x in ab.blocks() { + check(x, q).unwrap(); + acc.mac(&planes(x.mag), x.positive ^ qc.positive, &qc.planes); + } + acc.finish() + }; + + // D: A' products with the f32 detour replaced by a 16-bit carry column. + let routed: Vec<(U8x64, U8x64)> = ab + .blocks() + .zip(bb.blocks()) + .map(|(a, b)| route_a(a, b)) + .collect(); + let bw: Vec<(U16x32, U16x32)> = bb + .blocks() + .map(|b| (U16x32::from_u8x64_lo(b.mag), U16x32::from_u8x64_hi(b.mag))) + .collect(); + let a2_e2e = || { + 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 a2_ari = || { + let mut acc = AccA::new(); + for ((p, q), b) in routed.iter().zip(&bb.mag) { + acc.mac_mullo(*p, *q, *b); + } + acc.finish() + }; + fn d_e2e(ab: &Operand, bb: &Operand) -> i64 { + let mut acc = AccD::::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.step(a, b); + } + acc.finish() + } + fn d_ari(routed: &[(U8x64, U8x64)], bw: &[(U16x32, U16x32)]) -> i64 { + let mut acc = AccD::::new(); + for ((p, q), (lo, hi)) in routed.iter().zip(bw) { + acc.mac(*p, *q, *lo, *hi); + } + acc.finish() + } + fn d_scan(ab: &Operand, q: Block) -> i64 { + let qd = QueryD::new(q); + let mut acc = AccD::::new(); + for x in ab.blocks() { + check(x, q).unwrap(); + acc.scan_step(x, &qd); + } + acc.finish() + } + let dr_e2e = || d_e2e::(&ab, &bb); + let da_e2e = || d_e2e::(&ab, &bb); + let dr_ari = || d_ari::(&routed, &bw); + let da_ari = || d_ari::(&routed, &bw); + let dr_scan = || d_scan::(&ab, q); + let da_scan = || d_scan::(&ab, q); + + for (name, f, w) in [ + ("A' mullo", &a2_e2e as &dyn Fn() -> i64, want), + ("A' mullo arith", &a2_ari, want), + ("D react", &dr_e2e, want), + ("D anticipate", &da_e2e, want), + ("D react arith", &dr_ari, want), + ("D anticipate arith", &da_ari, want), + ("D react scan", &dr_scan, want_scan), + ("D anticipate scan", &da_scan, want_scan), + ] { + assert_eq!(f(), w, "{name} disagrees with the oracle"); + } + for (name, f, w) in [ + ("B", &b_e2e as &dyn Fn() -> i64, want), + ("B arith", &b_ari, want), + ("C abacus", &c_e2e_abacus, want), + ("C popcnt", &c_e2e_pop, want), + ("C abacus arith", &c_ari_abacus, want), + ("C abacus norm each", &c_ari_abacus_norm_each, want), + ("C popcnt arith", &c_ari_pop, want), + ("C csa arith", &c_ari_csa, want), + ("C csa", &c_e2e_csa, want), + ("C csa scan", &c_scan_csa, want_scan), + ("B scan", &b_scan, want_scan), + ("C abacus scan", &c_scan_abacus, want_scan), + ("C popcnt scan", &c_scan_pop, want_scan), + ] { + assert_eq!(f(), w, "{name} disagrees with the oracle"); + } + + println!("-- {label}: {n} blocks of 64 values per operand"); + time("plain u8 x i8 dot4 against one fixed query (baseline)", n, runs, None, &base); + let be = time("B end-to-end: check + AND/XOR split + 4 masked dot4", n, runs, None, &b_e2e); + time("C abacus end-to-end: check + planes(a,b) + CSA columns", n, runs, Some(be), &c_e2e_abacus); + time("C popcnt end-to-end: check + planes(a,b) + 128 popcounts", n, runs, Some(be), &c_e2e_pop); + time("C csa end-to-end: check + planes(a,b) + static CSA schedule", n, runs, Some(be), &c_e2e_csa); + let ba = time("B arithmetic only (split operands given)", n, runs, None, &b_ari); + let ca = time("C abacus arithmetic only (planes given), deferred", n, runs, Some(ba), &c_ari_abacus); + time("C abacus arithmetic only, normalized every block", n, runs, Some(ba), &c_ari_abacus_norm_each); + let cp = time("C popcnt arithmetic only (planes given)", n, runs, Some(ba), &c_ari_pop); + let cc = time("C csa arithmetic only (planes given), static schedule", n, runs, Some(ba), &c_ari_csa); + let prep = time("C plane extraction alone, both operands (cold)", n, runs, None, &c_prep); + let bs = time("B scan: query split once", n, runs, None, &b_scan); + time("C abacus scan: query planes once, index planes per block", n, runs, Some(bs), &c_scan_abacus); + time("C popcnt scan: query planes once, index planes per block", n, runs, Some(bs), &c_scan_pop); + time("C csa scan: query planes once, index planes per block", n, runs, Some(bs), &c_scan_csa); + println!(" D: A' products, 16-bit place-value column with carry (no f32)"); + time("A' end-to-end: route + mullo + f32 add (#388)", n, runs, Some(be), &a2_e2e); + time("D react end-to-end: route + mullo + add, carry after", n, runs, Some(be), &dr_e2e); + time("D anticipate end-to-end: route + mullo + add, carry before", n, runs, Some(be), &da_e2e); + time("A' arithmetic only (routed operands given)", n, runs, Some(ba), &a2_ari); + time("D react arithmetic only (routed, b widened given)", n, runs, Some(ba), &dr_ari); + time("D anticipate arithmetic only (routed, b widened given)", n, runs, Some(ba), &da_ari); + time("D react scan: query widened once", n, runs, Some(bs), &dr_scan); + time("D anticipate scan: query widened once", n, runs, Some(bs), &da_scan); + let best_c = ca.min(cp).min(cc); + if best_c >= be { + println!( + "break-even: none. C with both operands' planes prepared for free ({best_c:.3}) is slower than B paying all its preparation ({be:.3}); extraction alone costs {prep:.3}." + ); + } else { + println!("break-even: C arithmetic {best_c:.3} < B end-to-end {be:.3}; extraction {prep:.3} per use"); + } +} + +/// Candidate E (anchored at 128) against B, on the same data. Run alone with +/// the argument `e`. +fn bench_e(label: &str, n: usize, runs: usize, seed: u64) { + let (ra, rb) = data(n, seed); + let (ab, bb) = (Operand::new(&ra), Operand::new(&rb)); + let want = oracle(&ra, &rb).unwrap(); + let want_scan = oracle(&ra, &vec![rb[0]; n]).unwrap(); + let q = bb.block(0); + let split: Vec<(u64, I8x64, I8x64)> = ab + .blocks() + .zip(bb.blocks()) + .map(|(a, b)| prepare_b(a, b)) + .collect(); + let anc: Vec<(u64, I8x64, u64)> = ab + .blocks() + .zip(bb.blocks()) + .map(|(a, b)| { + let (low7, top) = split_anchor(b.mag); + (a.positive ^ b.positive, low7, top) + }) + .collect(); + + let b_e2e = || { + 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 b_ari = || { + let mut acc = AccB::new(); + for (a, (ng, l, t)) in ab.mag.iter().zip(&split) { + acc.mac(*a, *ng, *l, *t); + } + acc.finish() + }; + let b_scan = || { + 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() + }; + macro_rules! e_rows { + ($acc:ident) => {{ + let e2e = || { + let mut acc = $acc::new(); + for (a, b) in ab.blocks().zip(bb.blocks()) { + check(a, b).unwrap(); + acc.step(a, b); + } + acc.finish() + }; + let ari = || { + let mut acc = $acc::new(); + for (a, (ng, l, t)) in ab.mag.iter().zip(&anc) { + acc.mac(*a, *ng, *l, *t); + } + acc.finish() + }; + let scan = || { + let qe = QueryE::new(q); + let mut acc = $acc::new(); + for x in ab.blocks() { + check(x, q).unwrap(); + acc.scan_step(x, &qe); + } + acc.finish() + }; + assert_eq!(e2e(), want, stringify!($acc)); + assert_eq!(ari(), want, stringify!($acc)); + assert_eq!(scan(), want_scan, stringify!($acc)); + (e2e, ari, scan) + }}; + } + let (e4_e2e, e4_ari, e4_scan) = e_rows!(AccE4); + let (e3_e2e, e3_ari, e3_scan) = e_rows!(AccE3); + let (es_e2e, es_ari, es_scan) = e_rows!(AccESum); + assert_eq!(b_e2e(), want); + assert_eq!(b_ari(), want); + assert_eq!(b_scan(), want_scan); + + println!("-- E, {label}: {n} blocks of 64 values per operand"); + let be = time("B end-to-end: 4 masked dot4 (top limb x -128)", n, runs, None, &b_e2e); + time("E4 end-to-end: 2 dot4 low7 + 2 dot4 x ones, << 7", n, runs, Some(be), &e4_e2e); + time("E3 end-to-end: 2 dot4 low7 + 1 dot4 x signed ones", n, runs, Some(be), &e3_e2e); + time("Esum end-to-end: 2 dot4 low7 + 2 masked sum_bytes", n, runs, Some(be), &es_e2e); + let ba = time("B arithmetic only", n, runs, None, &b_ari); + time("E4 arithmetic only", n, runs, Some(ba), &e4_ari); + time("E3 arithmetic only", n, runs, Some(ba), &e3_ari); + time("Esum arithmetic only", n, runs, Some(ba), &es_ari); + let bs = time("B scan", n, runs, None, &b_scan); + time("E4 scan", n, runs, Some(bs), &e4_scan); + time("E3 scan", n, runs, Some(bs), &e3_scan); + time("Esum scan", n, runs, Some(bs), &es_scan); +} + +#[inline(never)] +pub fn asm_e4_mac(acc: &mut AccE4, a: U8x64, negative: u64, low7: I8x64, top: u64) { + acc.mac(a, negative, low7, top) +} +#[inline(never)] +pub fn asm_e3_mac(acc: &mut AccE3, a: U8x64, negative: u64, low7: I8x64, top: u64) { + acc.mac(a, negative, low7, top) +} +#[inline(never)] +pub fn asm_b_mac(acc: &mut AccB, a: U8x64, negative: u64, low7: I8x64, top: I8x64) { + acc.mac(a, negative, low7, top) +} + +// Out-of-line copies for the disassembly (`objdump -d --no-show-raw-insn`). +#[inline(never)] +pub fn asm_carry_save3(a: u64, b: u64, c: u64) -> (u64, u64) { + carry_save3(a, b, c) +} +#[inline(never)] +pub fn asm_planes(mag: U8x64) -> [u64; 8] { + planes(mag) +} +#[inline(never)] +pub fn asm_abacus_mac(acc: &mut AccAbacus, a: &[u64; 8], negative: u64, b: &[u64; 8]) { + acc.mac(a, negative, b) +} +#[inline(never)] +pub fn asm_csa_mac(acc: &mut AccCsa, a: &[u64; 8], negative: u64, b: &[u64; 8]) { + acc.mac(a, negative, b) +} +#[inline(never)] +pub fn asm_d_react(acc: &mut AccD, a_pos: U8x64, a_neg: U8x64, lo: U16x32, hi: U16x32) { + acc.mac(a_pos, a_neg, lo, hi) +} +#[inline(never)] +pub fn asm_d_anticipate(acc: &mut AccD, a_pos: U8x64, a_neg: U8x64, lo: U16x32, hi: U16x32) { + acc.mac(a_pos, a_neg, lo, hi) +} +#[inline(never)] +pub fn asm_popcnt_mac(acc: &mut AccPopcnt, a: &[u64; 8], negative: u64, b: &[u64; 8]) { + acc.mac(a, negative, b) +} + +fn main() { + println!("tier: {} (target-cpu as compiled)", tier()); + { + let (a, l) = (black_box(U8x64::splat(200)), black_box(I8x64::splat(5))); + let mut e4 = AccE4::new(); + asm_e4_mac(&mut e4, a, black_box(0x55), l, black_box(0xF0)); + let mut e3 = AccE3::new(); + asm_e3_mac(&mut e3, a, black_box(0x55), l, black_box(0xF0)); + let mut b = AccB::new(); + asm_b_mac(&mut b, a, black_box(0x55), l, l); + black_box((e4.finish(), e3.finish(), b.finish())); + } + if std::env::args().any(|x| x == "e") { + bench_e("L1-resident, repeated 4096-value dot", 64, 2001, 1); + bench_e("streamed", 1 << 20, 5, 2); + return; + } + { + // Keep the out-of-line copies in the binary. + let p = black_box(asm_planes(black_box(U8x64::splat(0xA5)))); + let z = black_box(U16x32::splat(7)); + let mut dr = AccD::::new(); + asm_d_react(&mut dr, black_box(U8x64::splat(3)), black_box(U8x64::splat(4)), z, z); + let mut da = AccD::::new(); + asm_d_anticipate(&mut da, black_box(U8x64::splat(3)), black_box(U8x64::splat(4)), z, z); + black_box((dr.finish(), da.finish())); + let mut s = AccCsa::new(); + asm_csa_mac(&mut s, &p, black_box(0x55), &p); + black_box(s.finish()); + let mut a = AccAbacus::new(); + asm_abacus_mac(&mut a, &p, black_box(0x55), &p); + let mut c = AccPopcnt::new(); + asm_popcnt_mac(&mut c, &p, black_box(0x55), &p); + black_box((asm_carry_save3(black_box(1), black_box(2), black_box(3)), a.finish(), c.finish())); + } + bench("L1-resident, repeated 4096-value dot", 64, 2001, 1); + bench("streamed", 1 << 20, 5, 2); +} diff --git a/tests/u9_abacus.rs b/tests/u9_abacus.rs new file mode 100644 index 00000000..81ce2f99 --- /dev/null +++ b/tests/u9_abacus.rs @@ -0,0 +1,403 @@ +#![cfg(feature = "std")] +//! Exactness of the D-U9-ABACUS probe (`examples/u9_abacus/abacus.rs`) +//! against the `i64` oracle of the U9 head-to-head and against its B kernel. + +#[path = "../examples/u9_head_to_head/kernels.rs"] +mod kernels; +#[path = "../examples/u9_abacus/abacus.rs"] +mod abacus; +#[path = "../examples/u9_abacus/carry16.rs"] +mod carry16; +#[path = "../examples/u9_abacus/anchor.rs"] +mod anchor; + +use abacus::{carry_save3, planes, Abacus, AccAbacus, AccCsa, AccPopcnt, QueryC}; +use anchor::{bit7, sign_bytes, split_anchor, AccE3, AccE4, AccESum, QueryE}; +use carry16::{carry_after, carry_before, AccD, QueryD, D_FLUSH_BLOCKS}; +use kernels::{check, oracle, AccB, Block, NeverTouched, B_FLUSH_BLOCKS}; +use ndarray::simd::{U16x32, U8x64}; + +type Raw = ([u8; 64], u64); + +fn block(r: &Raw) -> Block { + Block { + mag: U8x64::from_array(r.0), + positive: r.1, + } +} + +fn encode(v: i32) -> (u8, bool) { + (v.unsigned_abs() as u8, v >= 0) +} + +/// The three C folds, the C scan (b as the prepared query) and B, on the same data. +fn run_all(a: &[Raw], b: &[Raw]) -> [i64; 5] { + let (mut abacus, mut pop, mut bk) = (AccAbacus::new(), AccPopcnt::new(), AccB::new()); + let mut csa = AccCsa::new(); + let mut scan = AccAbacus::new(); + for (ra, rb) in a.iter().zip(b) { + let (x, y) = (block(ra), block(rb)); + check(x, y).unwrap(); + let (pa, pb) = (planes(x.mag), planes(y.mag)); + let negative = x.positive ^ y.positive; + abacus.mac(&pa, negative, &pb); + pop.mac(&pa, negative, &pb); + csa.mac(&pa, negative, &pb); + let q = QueryC::new(y.mag, y.positive); + scan.mac(&planes(x.mag), x.positive ^ q.positive, &q.planes); + bk.step(x, y); + } + [abacus.finish(), pop.finish(), csa.finish(), scan.finish(), bk.finish()] +} + +/// The full-adder identity bit by bit, for all eight input combinations at +/// every bit position, and on random words (summed over all bits). +#[test] +fn carry_save3_is_a_full_adder_per_bit() { + for combo in 0u8..8 { + let bit = |n: u8| if combo >> n & 1 == 1 { u64::MAX } else { 0 }; + let (s, c) = carry_save3(bit(0), bit(1), bit(2)); + let want = u32::from(combo.count_ones() as u8); + for k in 0..64 { + let got = (s >> k & 1) + 2 * (c >> k & 1); + assert_eq!(got as u32, want, "combo {combo:03b} bit {k}"); + } + } + let mut x = 0x243F_6A88_85A3_08D3u64; + for _ in 0..10_000 { + let mut r = || { + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + x + }; + let (a, b, c) = (r(), r(), r()); + let (s, cy) = carry_save3(a, b, c); + assert_eq!(a.count_ones() + b.count_ones() + c.count_ones(), s.count_ones() + 2 * cy.count_ones()); + } +} + +/// Planes recover every magnitude 0..=255 in every lane. +#[test] +fn planes_recover_every_magnitude() { + for base in [0u8, 64, 128, 192] { + let m: [u8; 64] = core::array::from_fn(|i| (i as u8).wrapping_add(base)); + let p = planes(U8x64::from_array(m)); + for (k, &mk) in m.iter().enumerate() { + let back: u8 = (0..8).map(|i| ((p[i] >> k & 1) as u8) << i).sum(); + assert_eq!(back, mk); + } + } +} + +/// Carries forced through many consecutive columns: 2^20 all-ones planes at +/// column 0 must count exactly 64 * 2^20, and one all-ones plane at each of +/// columns 0..=62 must not lose the chain of carries it sets off. +#[test] +fn abacus_carries_through_many_columns() { + let mut ab = Abacus::new(); + for _ in 0..(1u32 << 20) { + ab.add(0, u64::MAX); + } + assert_eq!(ab.value(), 64i128 << 20); + let mut ab = Abacus::new(); + let mut want = 0i128; + for round in 0..3 { + for w in 0..40 { + ab.add(w, u64::MAX); + want += 64i128 << w; + } + assert_eq!(ab.value(), want, "round {round}"); + } + // A single-lane pattern: each lane gets a different count. + let mut ab = Abacus::new(); + let mut want = 0i128; + for t in 0..5000u64 { + let p = t.wrapping_mul(0x9E37_79B9_7F4A_7C15); + ab.add((t % 7) as usize, p); + want += i128::from(p.count_ones()) << (t % 7); + } + assert_eq!(ab.value(), want); +} + +#[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, mut bm) = ([0u8; 64], [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 want = i64::from(va) * i64::from(vb); + assert_eq!(run_all(&[(am, ap)], &[(bm, bp)]), [want; 5], "{va} x {vb}"); + n += 1; + } + } + assert_eq!(n, 511 * 511); +} + +#[test] +fn full_blocks_match_the_oracle_and_b() { + let mut s = 0x9E37_79B9_7F4A_7C15u64; + let mut val = move || { + s = s + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + ((s >> 33) % 511) as i32 - 255 + }; + for blocks in [1usize, 2, 7, 64, 300] { + let mut mk = || { + 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()).collect(); + let b: Vec<_> = (0..blocks).map(|_| mk()).collect(); + let want = oracle(&a, &b).unwrap(); + assert_eq!(run_all(&a, &b), [want; 5], "{blocks} blocks"); + } +} + +/// +-255 x +-255 in every lane, past B's flush bound, totals beyond `i32`. +#[test] +fn worst_case_runs_exceed_i32_exactly() { + let n = 3 * B_FLUSH_BLOCKS + 5; + 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); + assert_eq!(run_all(&a, &b), [want; 5]); + } +} + +/// Never-touched is refused before any arithmetic; touched zero is a number. +#[test] +fn never_touched_is_refused_and_zero_is_not() { + let zero = ([0u8; 64], u64::MAX); + assert_eq!(run_all(&[zero], &[zero]), [0; 5]); + let mut m = [9u8; 64]; + m[5] = 0; + let sentinel = block(&(m, !(1u64 << 5))); + assert_eq!(check(sentinel, block(&zero)), Err(NeverTouched(1 << 5))); + assert_eq!(check(block(&zero), sentinel), Err(NeverTouched(1 << 5))); +} + +/// D (react, anticipate) and their scans (b as the prepared query). +fn run_d(a: &[Raw], b: &[Raw]) -> [i64; 4] { + let (mut r, mut an) = (AccD::::new(), AccD::::new()); + let (mut sr, mut sa) = (AccD::::new(), AccD::::new()); + for (ra, rb) in a.iter().zip(b) { + let (x, y) = (block(ra), block(rb)); + check(x, y).unwrap(); + r.step(x, y); + an.step(x, y); + let q = QueryD::new(y); + sr.scan_step(x, &q); + sa.scan_step(x, &q); + } + [r.finish(), an.finish(), sr.finish(), sa.finish()] +} + +/// Both carry rules against `u32` arithmetic, for every `low` in 0..=65535, +/// with `p` on both sides of the headroom (`!low`, `!low + 1`) and spread +/// over the product range. +#[test] +fn carry_rules_match_u32_for_every_low() { + let mut x = 0x2545_F491_4F6C_DD1Du64; + for base in (0u32..65_536).step_by(32) { + let lows: [u16; 32] = core::array::from_fn(|k| (base + k as u32) as u16); + let low = U16x32::from_array(lows); + let mut cases: Vec<[u16; 32]> = vec![ + core::array::from_fn(|k| !lows[k]), + core::array::from_fn(|k| (!lows[k]).wrapping_add(1)), + core::array::from_fn(|k| (!lows[k]).wrapping_sub(1)), + [0; 32], + [65_025; 32], + ]; + for _ in 0..4 { + cases.push(core::array::from_fn(|_| { + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + (x % 65_026) as u16 + })); + } + for ps in cases { + let p = U16x32::from_array(ps); + let sum = low + p; + let (after, before) = (carry_after(low, p, sum).to_array(), carry_before(low, p).to_array()); + for k in 0..32 { + let want = u16::from(u32::from(lows[k]) + u32::from(ps[k]) > 0xFFFF); + assert_eq!(after[k], want, "react low {} p {}", lows[k], ps[k]); + assert_eq!(before[k], want, "anticipate low {} p {}", lows[k], ps[k]); + } + } + } +} + +#[test] +fn d_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, mut bm) = ([0u8; 64], [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 want = i64::from(va) * i64::from(vb); + assert_eq!(run_d(&[(am, ap)], &[(bm, bp)]), [want; 4], "{va} x {vb}"); + n += 1; + } + } + assert_eq!(n, 511 * 511); +} + +#[test] +fn d_full_blocks_match_the_oracle() { + let mut s = 0x0123_4567_89AB_CDEFu64; + let mut val = move || { + s = s + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + ((s >> 33) % 511) as i32 - 255 + }; + for blocks in [1usize, 2, 7, 64, 300] { + let mut mk = || { + 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()).collect(); + let b: Vec<_> = (0..blocks).map(|_| mk()).collect(); + assert_eq!(run_d(&a, &b), [oracle(&a, &b).unwrap(); 4], "{blocks} blocks"); + } +} + +/// All-255 runs past D's flush bound: a carry on almost every add, `high` +/// close to its limit, totals far beyond `i32`. +#[test] +fn d_worst_case_runs_cross_the_flush_bound() { + let n = 2 * D_FLUSH_BLOCKS + 5; + for (pa, pb) in [(u64::MAX, u64::MAX), (u64::MAX, 0), (0xF0F0_F0F0_0F0F_0F0F, 0x3333_CCCC_5555_AAAA)] { + let a = vec![([255u8; 64], pa); n]; + let b = vec![([255u8; 64], pb); n]; + assert_eq!(run_d(&a, &b), [oracle(&a, &b).unwrap(); 4]); + } +} + +/// E4, E3, Esum and their scans (b as the prepared query). +fn run_e(a: &[Raw], b: &[Raw]) -> [i64; 6] { + let (mut e4, mut e3, mut es) = (AccE4::new(), AccE3::new(), AccESum::new()); + let (mut s4, mut s3, mut ss) = (AccE4::new(), AccE3::new(), AccESum::new()); + for (ra, rb) in a.iter().zip(b) { + let (x, y) = (block(ra), block(rb)); + check(x, y).unwrap(); + e4.step(x, y); + e3.step(x, y); + es.step(x, y); + let q = QueryE::new(y); + s4.scan_step(x, &q); + s3.scan_step(x, &q); + ss.scan_step(x, &q); + } + [e4.finish(), e3.finish(), es.finish(), s4.finish(), s3.finish(), ss.finish()] +} + +/// The anchor split for every magnitude: `|b| = low7 + 128 * bit7`, and the +/// sign bytes are `+1` / `-1` exactly by the plane. +#[test] +fn anchor_split_and_sign_bytes() { + for base in [0u8, 64, 128, 192] { + let m: [u8; 64] = core::array::from_fn(|i| (i as u8).wrapping_add(base)); + let (low7, top) = split_anchor(U8x64::from_array(m)); + assert_eq!(top, bit7(U8x64::from_array(m))); + let low7 = low7.to_array(); + for k in 0..64 { + let t = i32::from((top >> k & 1) as u8); + assert!((0..=127).contains(&low7[k])); + assert_eq!(i32::from(low7[k]) + 128 * t, i32::from(m[k]), "magnitude {}", m[k]); + } + } + let plane = 0xF0F0_1234_8000_0001u64; + let s = sign_bytes(plane).to_array(); + for (k, &v) in s.iter().enumerate() { + assert_eq!(v, if plane >> k & 1 == 1 { -1 } else { 1 }); + } +} + +#[test] +fn e_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, mut bm) = ([0u8; 64], [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 want = i64::from(va) * i64::from(vb); + assert_eq!(run_e(&[(am, ap)], &[(bm, bp)]), [want; 6], "{va} x {vb}"); + n += 1; + } + } + assert_eq!(n, 511 * 511); +} + +#[test] +fn e_full_blocks_match_the_oracle() { + let mut s = 0x7F4A_7C15_9E37_79B9u64; + let mut val = move || { + s = s + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + ((s >> 33) % 511) as i32 - 255 + }; + for blocks in [1usize, 2, 7, 64, 300] { + let mut mk = || { + 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()).collect(); + let b: Vec<_> = (0..blocks).map(|_| mk()).collect(); + assert_eq!(run_e(&a, &b), [oracle(&a, &b).unwrap(); 6], "{blocks} blocks"); + } +} + +#[test] +fn e_worst_case_runs_cross_the_flush_bound() { + let n = 3 * B_FLUSH_BLOCKS + 5; + 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]; + assert_eq!(run_e(&a, &b), [oracle(&a, &b).unwrap(); 6]); + } +}