diff --git a/.claude/blackboard.md b/.claude/blackboard.md index 5e44199b..88a610ab 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -3937,3 +3937,15 @@ accumulate; sign plane -> +1/-1 bytes without mask_blend. Exact on 7 tiers; 5 di 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). +## D-DOMINO-MOV-0: bit flip, positional substitution, lazy FULL carry (2026-10-29, branch claude/u9-domino, stacked on #389) + +Probe-only: `examples/u9_domino/{domino,main}.rs`, `tests/u9_domino.rs`, report `.claude/knowledge/u9-domino-probe.md`. +F (eager flip ripple on 64-lane bit planes) ADOPT-LAB: beats native per-lane adds on dense/single-lane masks (1.94 vs +6.99 ns/add), at the 64K boundary beats D only on AVX2 (5.7 vs 9.4 ns/64, D's array-backed U16x32 shuffles) and still +loses to native u64 lanes there. P/L/PL REFUSE: flag tests cost more than the XOR/AND they skip (1.6-1.9 vs 0.5 ns per +FULL column). R1 (two u64) and R2 (u128) compile to identical code. Coverage audit of #389: C's tests are lane-blind +(totals only); EMPTY flags only matter over stale words (Y2 was green until stale-word runs were added). 7 tests on 5 +tiers; 6 disable runs red. +Review fix (Codex, PR #391): PL's ordinary ripple now marks an all-ones result FULL (as L); new test + per-add FULL +invariant in check_all, Y7 red; 8 tests on 5 tiers. Re-measured on 3 tiers: PL ~2.2 ns per FULL column; verdicts unchanged. + diff --git a/.claude/knowledge/u9-domino-probe.md b/.claude/knowledge/u9-domino-probe.md new file mode 100644 index 00000000..611404bc --- /dev/null +++ b/.claude/knowledge/u9-domino-probe.md @@ -0,0 +1,225 @@ +# D-DOMINO-MOV-0: bit flip, positional substitution and lazy FULL carry + +READ BY: savant-architect, vector-synthesis, simd-savant, the TurboVec-U9 owner. Follows `u9-abacus-probe.md` (#389). + +Scope: branch `claude/u9-domino`, stacked on #389 (base `claude/u9-abacus`). Files: `examples/u9_domino/{domino.rs, main.rs}` and `tests/u9_domino.rs`. Plain Rust, no facade, backend or storage change, no intrinsics. #388 and #389 are untouched. + +All numbers were measured on one AVX-512 VNNI Xeon. A range is the min - max of 3 runs of medians. + +The tables are a second full measurement, taken after the PL fix below (review, Codex, PR #391): PL's ordinary ripple stored an all-ones result without marking the column FULL, unlike L, so a later all-ones carry read and wrote the word instead of passing it as FULL -> EMPTY. The fix makes PL's L1 dense row 6.4 - 9.6 ns (was 6.3 - 6.4) and its per-FULL-column cost about 2.2 ns (was 1.85); verdicts are unchanged. The second run also differs where no code changed: F on random masks reads 8.1 - 8.2 ns on VNNI (was 6.4), and the streamed rows sit higher for every counter alike, so those differences are host and code-layout variance, not the fix. + +## Verdict + +- **F (Eager): ADOPT-LAB.** A bit-plane counter for dense and sparse masks. At the 64K boundary it pays off only on AVX2, where D's lowering is poor, and even there it does not beat native u64 lanes. +- **P (Subst): REFUSE.** Assignment into an EMPTY column compiles to a store anyway, and the `empty` bitmap check costs more than the XOR/AND it avoids. +- **L (LazyFull): REFUSE.** The FULL bookkeeping costs more per column than the plain ripple it replaces. +- **PL (SubstLazy): REFUSE.** It pays both bookkeeping costs and is slower than F in every measured case. + +Answer to the final research question: the bit flip crosses 64K exactly, and R1 (two u64 planes) and R2 (one u128) cost the same. But a two-plane abacus does not make carry escalation cheaper than native addition in general. LLVM already compiles native lane adds to vector code, and the lazy and substitution bookkeeping costs more than the XOR/AND it avoids. The measured exception is dense or single-lane masks fed one at a time, where F's short ripple beats 64 per-lane adds. + +## The question + +**H1, positional substitution.** When a carry lands in a column that is known to be empty, the result is just the carry mask. The column can be assigned without being read, and no further carry is produced. + +**H2, lazy binary carry.** A column that an XOR leaves all ones is FULL. Adding a mask m to a FULL column gives `!m` and carries m, and an all-ones carry through a FULL column only changes its state to EMPTY, with no data read or written. + +## Invariants + +- Bit flip: `b + c = (b ^ c) + 2(b & c)` per lane. +- Substitution: EMPTY column + m = m, no carry. +- Lazy FULL: FULL + m = `!m`, carry m. FULL + all-ones = EMPTY, carry all-ones. +- Finding: with zero-initialized words, the physical word under an EMPTY or FULL flag is always 0 in PL. The EMPTY flag therefore carries information only when words can hold stale data (reused buffers). The tests also run P and PL over stale words (`0xDEAD_BEEF_F00D_CAFE`). + +## Implementations + +All in `domino.rs`. + +- `flip(old, carry) = (old ^ carry, old & carry)`: a half adder on 64 lanes, bit j = lane j. +- `Counter`: N columns of u64. Column k has significance 2^k for all 64 lanes. `add_at(k, m)` adds 2^k to the lanes in m. The carry mask moves to column k+1 unshifted. +- **F `Eager`:** loops `(col[k], m) = flip(col[k], m); k += 1` while m != 0. Returns `Err(Overflow(m))` if a carry leaves column N-1. +- **P `Subst`:** F plus an `empty: u128` bitmap. Adding into an EMPTY column assigns `col[k] = m` without reading it and stops. `Subst::stale(x)` builds a zero counter over words holding garbage x. +- **L `LazyFull`:** F plus a `full: u128` bitmap. A column becomes FULL when an XOR leaves it all ones. Adding m to FULL writes `!m` without reading and carries m. +- **PL `SubstLazy`:** both bitmaps. An all-ones carry through a FULL column is a state change FULL -> EMPTY only. EMPTY + m assigns, or sets FULL if m is all ones. +- **Native:** `[u64; 64]` per-lane integers, `lane += (m >> j & 1) << k`. +- **64K boundary:** low16 as two U16x32 (wrapping add). The carry vector comes from #389's `carry_after`, becomes a 64-lane mask with `U16x32::cmpeq_mask`, and is added at column 0 (= bit 16) of a bit-sliced counter. Exact value = sum(low lanes) + 2^16 * counter total. +- **R1 `Pair128 {lo, hi}`:** two significance planes (64 independent 2-bit counters). **R2:** the same planes in one u128 (low half = 2^0 plane). + +## Coverage matrix of #388/#389 + +Audited from source: 15 tests in `tests/u9_abacus.rs`, 5 in `tests/u9_head_to_head.rs`. + +| case | status | evidence | +|---|---|---| +| crossing 2^16 | PROVEN | `carry_rules_match_u32_for_every_low`: every low 0..=65535, p at !low, !low+1, !low-1, 0, 65025, random; D 511x511, full blocks, worst case | +| flipping the bit-16 position | MISSING | D keeps high as a u16 count (`high += carry`); no bit-16 plane exists or is tested | +| carry into bit 17 when bit 16 is set | PARTIAL | implied by D's worst-case run (high counts near 65k) and C's 2^20 run; never asserted at bit level | +| repeated carries must not cancel | PARTIAL | numeric totals only (D worst case, C 2^20 all-ones) | +| eager vs deferred normalization | PARTIAL | C timed both; correctness only for deferred; D deferred (flush 65,535) only | +| consecutive FULL columns | PARTIAL | `abacus_carries_through_many_columns`: up to ~20 columns, all-lane masks only; no 32/64 chains, no single-lane chains | +| assignment into a known-empty destination | PARTIAL | `Abacus::add` stores into an empty slot without reading; numerically covered, never isolated or measured | +| carry mask forwarded without moving lanes | PARTIAL | every C test compares only lane-blind totals (`value()` = sum of popcount * 2^w); a carry rotated into another lane would pass every C test | +| two 64-bit planes in one 128-bit carrier | MISSING | | +| read-after-write / observation between adds | PARTIAL | abacus test reads `value()` between 3 rounds; D only at finish | + +D's `high += carry`: numeric correctness is PROVEN. The bit-flip formulation was never measured (confirmed). + +The harvest was verified: 5 quoted lines were re-read at their stated places, and the counts 15 and 5 were re-counted. + +## New tests + +`tests/u9_domino.rs`, 8 tests, all pass on native (AVX-512 VNNI), v4, v3, v2 (qemu) and aarch64 (qemu); clippy `-D warnings` and the shuffle guard are clean. Test counters use 96 columns (oracle u128). + +- `flip_is_a_half_adder_in_every_lane`: 4 input pairs, full and zero masks, single lanes 0/17/63, 1000 random words checked per lane. +- `random_adds_match_the_oracle_in_any_order`: 3000 adds, masks zero/full/one lane/random, significances 0, 1, 15, 16, 17, 31, 63, 64, 65. A per-lane u128 oracle with checked arithmetic is compared after every add on lanes 0, 1, 31, 32, 62, 63, for F, P, L, PL, P-stale and PL-stale. The same run in reverse order. All six give identical columns in both orders. +- `carries_run_through_long_full_chains`: chains of 1, 2, 8, 16, 32, 64 FULL columns, all-lane and single-lane (lane 37) masks, then two successive carries into column 0. Exact column contents are asserted (2^n + 1 in the mask lanes, 0 elsewhere, no other lane touched). +- `the_64k_boundary_flips_bit_16_then_carries_into_bit_17`: low16 in {0, 1, 510, 511, 512, 65000, 65534, 65535} x p in {0, 1, 511, 65024, 65025} x bit 16 preset or not. Value `low + 2^16*high` is exact in lanes 0, 31, 32, 63. Bit 16 flips exactly when crossing XOR preset. Bit 17 is set exactly when crossing AND preset. +- `product_sums_match_native_and_d`: 1, 300 random and 300, 70,000 dense (65,025) blocks. Native, D, F, P, L, PL, each with deferred and per-block normalization, all equal the u128 sum. +- `overflow_past_the_last_column_is_refused`: 4-column counters of every variant. A carry past the top returns `Err(Overflow(lanes))`. An add at k=N with m=0 is Ok. +- `xor_made_full_column_is_full_and_carries_by_state` (added after review): complementary masks at one significance leave a column all ones by XOR; L and PL both mark it FULL, and the next all-ones carry turns it EMPTY and lands in column 1 as FULL, lanes exact. `check_all` also asserts after every add that L, PL and PL-stale hold exactly the logically all-ones columns as FULL. +- `pair_and_packed_agree_and_roundtrip`: 10,000 random adds, R1 and R2 identical, roundtrip identity, per-lane 2-bit counter check (64 independent counters, not one integer). + +## Disable runs + +Each disable must turn tests red. + +| id | disable | result | +|---|---|---| +| Y1 | F carry not propagated | RED (5 tests) | +| Y2 | PL lost the FULL->EMPTY transition | GREEN at first, RED after adding stale-word runs (2 tests) | +| Y3 | P assigns into a non-empty column | RED (4 tests) | +| Y4 | L treats FULL as unchanged after a partial carry | RED (3 tests) | +| Y5 | carry mask rotated one lane (R1) | RED (1 test) | +| Y6 | P reads the stale word instead of assigning | RED (2 tests) | +| Y7 | PL ordinary ripple does not mark an all-ones result FULL (the pre-review code) | RED (2 tests) | + +**The Y2 story.** The first Y2 run was green. With zero-initialized words, the physical word under an EMPTY or FULL flag is always 0 in PL, so a lost FULL->EMPTY transition left the data identical to the correct result and nothing could tell the difference. Only after P and PL were also run over stale words (`0xDEAD_BEEF_F00D_CAFE`) did the same disable turn red, in 2 tests. + +This shows two things. A green disable run on a lazy structure may mean the tests cannot see the state, not that the guard is unneeded. And the EMPTY flag is the only thing that makes stale data harmless, so a test that never feeds stale data cannot cover it. + +## Measurements + +How to run: + +```sh +cargo run --release --example u9_domino --features std +cargo test --release --test u9_domino +``` + +### A. u16 products summed per lane across the 64K boundary (ns per 64 products) + +L1 = 64 blocks. Streamed = 2^20 blocks (128 MiB). Dense = products 65,024-65,025 (a carry on almost every add). Sparse = products 0..=255 (a carry about once per 260 blocks per lane). + +| L1 dense | AVX-512 VNNI | AVX-512 (v4) | AVX2 (v3) | +|---|---|---|---| +| native u64 lanes | 4.81 - 4.83 | 4.83 - 4.84 | 4.83 - 4.84 | +| D (#389) | 2.28 - 2.56 | 2.30 - 2.66 | 9.39 - 9.42 | +| F | 2.70 - 2.88 | 3.05 - 3.14 | 5.66 - 5.67 | +| P | 5.44 - 5.47 | 5.78 - 5.95 | 7.61 - 7.81 | +| L | 5.89 - 6.22 | 6.55 - 6.61 | 8.50 - 8.56 | +| PL | 6.36 - 9.58 | 6.67 - 6.91 | 8.80 - 8.98 | +| P normalized every block | 74.7 - 74.8 | 75.0 | 68.0 - 68.8 | +| PL normalized every block | 107 - 110 | 111 - 112 | 115 - 116 | + +| L1 sparse | VNNI | v4 | v3 | +|---|---|---|---| +| native | 4.83 - 5.58 | 4.81 - 4.86 | 4.83 - 4.84 | +| D | 2.31 - 3.19 | 2.22 - 2.67 | 9.42 | +| F | 2.59 - 2.94 | 2.94 - 3.03 | 5.61 - 5.77 | +| P | 3.64 - 4.91 | 3.95 - 3.98 | 6.53 - 6.70 | +| L | 3.62 - 3.83 | 3.95 - 3.97 | 6.61 - 6.75 | +| PL | 3.92 - 4.00 | 4.20 - 4.36 | 7.05 - 7.06 | + +| streamed dense | VNNI | v4 | v3 | +|---|---|---|---| +| native | 15.5 - 16.7 | 15.7 - 16.1 | 14.2 - 16.6 | +| D | 11.3 - 12.4 | 11.3 - 12.2 | 31.6 - 32.4 | +| F | 12.5 - 13.7 | 13.6 - 14.4 | 15.8 - 20.4 | +| P | 13.8 - 15.9 | 14.6 - 14.8 | 14.6 - 25.8 | +| L | 13.9 - 16.5 | 13.4 - 14.7 | 13.6 - 27.3 | +| PL | 13.3 - 15.2 | 13.8 - 15.1 | 12.2 - 26.2 | + +| streamed sparse | VNNI | v4 | v3 | +|---|---|---|---| +| native | 15.5 - 16.5 | 15.5 - 16.5 | 16.0 - 16.7 | +| D | 10.7 - 12.1 | 10.2 - 11.4 | 31.4 - 32.5 | +| F | 10.3 - 13.2 | 10.2 - 12.9 | 17.2 - 17.8 | +| P | 11.2 - 12.9 | 9.90 - 13.0 | 16.8 - 19.0 | +| L | 10.9 - 12.1 | 10.2 - 12.7 | 16.0 - 19.6 | +| PL | 11.1 - 12.3 | 10.9 - 12.5 | 16.2 - 19.3 | + +The streamed rows are noisy. The ranges overlap between F, P, L and PL. + +### B. Counters alone, ns per add (4,096 adds, 64 lanes) + +VNNI figures. v4 and v3 are similar unless noted. + +| workload | native | F | P | L | PL | +|---|---|---|---|---|---| +| dense all-ones at 2^0 (binary counting) | 6.99 - 7.04 | 1.93 - 1.96 (v4 1.94, v3 1.95 - 1.98) | 2.92 - 2.93 | 3.42 - 3.49 | 3.22 - 3.29 | +| random masks at 2^0 | 6.99 | 8.12 - 8.21 (v4 8.15 - 10.1, v3 8.05 - 8.21) | 10.5 - 10.8 | 12.4 | 15.5 - 15.6 | +| sparse one lane at 2^0 | 6.99 | 5.77 - 6.41 (v4 5.79 - 6.21, v3 6.10 - 6.29) | 6.17 - 6.33 | 6.14 - 7.05 | 7.14 - 7.65 | +| weighted: random masks entering at 2^0, 2^16, 2^17, 2^32 | 6.99 | 7.57 - 7.71 (v4 7.73 - 8.83, v3 7.64 - 8.21) | 10.2 - 10.8 | 11.8 - 12.8 | 14.5 - 14.8 | + +### B'. One add whose carry runs through n FULL columns + +ns per rep = clone + add. Each cell is "clone-only / clone+add", so subtract the clone-only figure. VNNI; v4 and v3 are within a few percent. + +| n | F | P | L | PL | +|---|---|---|---|---| +| 1 | 5.27 - 5.46 / 5.88 | 4.93 - 5.03 / 12.5 | 5.65 / 11.5 - 11.8 | 5.66 / 12.1 - 12.3 | +| 8 | 5.14 - 5.47 / 11.8 | 4.91 - 5.04 / 18.0 - 18.4 | 5.51 - 5.65 / 24.0 - 24.4 | 5.65 - 5.66 / 24.6 - 25.2 | +| 32 | 5.27 - 5.47 / 30.1 - 30.4 | 5.03 - 5.04 / 61.9 - 70.6 | 5.62 - 7.39 / 75.2 - 94.3 | 5.52 - 5.66 / 79.2 - 82.5 | +| 64 | 5.46 - 7.47 / 38.6 - 60.4 | 5.03 - 6.76 / 108 - 143 | 5.64 - 7.53 / 127 - 168 | 5.66 - 7.18 / 146 - 150 | + +Per FULL column (64-column row minus clone, medians, divided by 64): F about 0.52 ns, P about 1.6 ns, L about 1.9 ns, PL about 2.2 ns. + +"Constant state" rows (counter built from literals inside the timed loop, then add + total): LLVM did not fold the counter to a constant. F measured 101 ns and PL 347 - 349 ns at n=64 on VNNI. They include building the preset, so they are not comparable to the opaque rows above. + +### C. R1 vs R2 (ns per add, 4,096 adds) + +| | VNNI | v4 | v3 | +|---|---|---|---| +| R1 two u64 | 0.412 - 0.415 | 0.415 - 0.424 | 0.411 - 0.414 | +| R2 one u128 | 0.417 - 0.453 | 0.442 - 0.453 | 0.443 - 0.517 | + +## Assembly + +objdump of `#[inline(never)]` copies. Counts came from a Haiku census with the stated grep; 2 rows were re-run by hand and matched. + +| function | instructions | notes | +|---|---|---| +| flip | 5 | mov, mov, xor, and, ret; scalar GPR, no ternlog | +| R1 `pair_add` and R2 `packed_add` | 10 each, identical | native and v3: mov, xor, and, mov, xor/and, and/xor, 3 stores, ret. LLVM sees through the u128 packing: no shifts, no extracts, no extra moves | +| `eager_add` (F) | 47 - 59 | 10 - 12 branches, loop over the column array in memory; logic 14 - 18; no ternlog, no shuffles | +| `subst_add` (P) | 99 - 123 | 13 - 16 branches, 12 - 16 shifts (u128 bitmap tests), 6 - 9 push/pop | +| `lazy_add` (L) | 51 - 88 | 11 - 15 branches, 6 - 8 compares | +| `pl_add` (PL) | 86 | 15 branches, 8 compares, 8 push/pop | +| `native_add` (64 lanes += bit << k) | 68 - 69 | vectorized: 16 adds, 16 shifts, 16 logic, 1 branch, 2 vpbroadcastq | +| 64K block, AVX-512 (native and v4): D | 20 | 4 vpternlogq (the carry rule), no shuffle, 1 branch | +| 64K block, AVX-512 (native and v4): F | 68 - 79 | 4 vpternlogq, 10 - 12 branches (the ripple loop), 27 - 33 moves incl. mask transfers | +| 64K block, AVX2 (v3): D | 107 | 13 shuffle-class instructions (3 vextracti128, 1 vinserti128, 5 vpbroadcastq, 1 vperm2i128, 1 vpermq, 1 vpunpcklqdq) from the array-backed U16x32 | +| 64K block, AVX2 (v3): F | 142 | 18 shuffle-class, 10 branches | + +A first R1/R2 anchor was constant-propagated because it was called with a known-zero high half. It was fixed with an opaque input. + +## Reading it + +- **F vs D at the 64K boundary, AVX-512.** D is faster in L1 (2.2 - 2.6 vs 2.7 - 3.3 ns) because its whole high column is 4 vpternlogq plus adds, with no branch. F's ripple is a scalar loop with branches. Streamed, the two overlap (the run is memory-bound). +- **F vs D, AVX2.** F is 1.6 - 1.7x faster than D in L1 and about 2x streamed, because D's AVX2 U16x32 lowering moves halves around (13 shuffle-class ops). The native u64-per-lane sum (4.8 ns) still beats both F and D on AVX2 in L1. +- **The flip itself is cheap.** The cost is the data-dependent ripple loop and getting the carry out of the vector into a GPR mask. +- **P, L, PL are slower than F in every measured case.** Every metadata check is a branch plus a u128 shift/test. The gap grows with chain length: about 1.6 - 2.2 ns per FULL column against 0.5 ns for F. The deferred FULL->EMPTY transition skips a data read and write, but the flag read and write costs more. +- **Normalizing P or PL every block costs 70 - 120 ns per block** (normalize rewrites all 72 columns). Only the deferred form is usable, and even deferred it loses to F. +- **Positional substitution does not disappear in the compiler.** Assignment into an EMPTY column compiles to a store: the column array lives in memory indexed by a runtime k, so no register renaming happens. LLVM does not eliminate it. +- **R1 and R2 cost the same** (section C), and LLVM compiles them to identical code. +- **Limits.** One machine, one compiler build. The streamed rows are noisy. The "native beats F" statement on AVX2 is for the L1 rows (dense and sparse); streamed, the ranges overlap. + +## Break-even + +- No break-even for P, L or PL against F in any workload measured. +- F breaks even against native integer addition on dense and single-lane masks (3.6x faster on dense counting). It does not on random masks: about 16% slower on every tier (8.1 - 8.2 vs 7.0 ns). + +## Not done + +- NEON and wasm timings. +- Unrolled or vectorized ripple. +- Hand scheduling (out of scope by design). diff --git a/Cargo.toml b/Cargo.toml index aec13c9c..01fcd5af 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -76,6 +76,11 @@ name = "u9_abacus" path = "examples/u9_abacus/main.rs" required-features = ["std"] +[[example]] +name = "u9_domino" +path = "examples/u9_domino/main.rs" +required-features = ["std"] + [[example]] name = "cpu_guard_probe" required-features = ["std"] diff --git a/examples/u9_domino/domino.rs b/examples/u9_domino/domino.rs new file mode 100644 index 00000000..72a17717 --- /dev/null +++ b/examples/u9_domino/domino.rs @@ -0,0 +1,412 @@ +//! D-DOMINO-MOV-0: positional substitution and lazy binary carry on +//! bit-sliced columns. Probe-private. +//! +//! A column is one `u64`: bit `j` belongs to lane `j`, and every column of a +//! counter carries the same 64 lane identities at its own significance. A +//! carry mask moves to the next column as it is; no lane position changes. +//! +//! Invariants: +//! +//! * bit flip (half adder, per lane): adding carry `c` to bit `b` gives +//! `b' = b ^ c` and carries `b & c` to the next significance; +//! `b + c = b' + 2 * (b & c)` in every lane. +//! * positional substitution: if a column is logically zero (EMPTY), adding +//! mask `m` gives exactly `m` and no carry, so the column may be assigned +//! without reading or clearing it. +//! * lazy FULL: if a column is logically all ones (FULL), adding `m` gives +//! `!m` and carries `m`; when `m` is all ones the column becomes EMPTY and +//! the whole mask moves on, so only the state changes. +//! +//! Four counters over the same columns, all exact: +//! +//! * [`Eager`] (F): ripple `col ^= m; m = old & m` until the carry is empty; +//! * [`Subst`] (P): `Eager` plus an EMPTY bitmap, assignment into EMPTY columns; +//! * [`LazyFull`] (L): `Eager` plus a FULL bitmap, an all-ones carry passes a +//! FULL column by writing zero, without reading it; +//! * [`SubstLazy`] (PL): both bitmaps; an all-ones carry passes FULL columns +//! by changing state bits only, the data words are not touched. + +#![allow(dead_code)] + +/// Half adder on 64 lanes: `(old ^ carry, old & carry)`. +#[inline(always)] +pub fn flip(old: u64, carry: u64) -> (u64, u64) { + (old ^ carry, old & carry) +} + +/// A carry left the highest column: the lanes in the mask would need one more +/// column. The counter then holds its value modulo `2^N` in those lanes. +#[derive(Debug, PartialEq, Eq)] +pub struct Overflow(pub u64); + +/// The interface all four counters share. +pub trait Counter { + fn new() -> Self; + /// Adds 1 at significance `2^k` to every lane set in `m`. + fn add_at(&mut self, k: usize, m: u64) -> Result<(), Overflow>; + /// The logical columns (EMPTY / FULL resolved). + fn columns(&self) -> [u64; N]; + /// Writes the logical columns back and clears any state metadata. + fn normalize(&mut self); + + /// Value of lane `j`, exact up to 128 columns. + fn lane(&self, j: usize) -> u128 { + let c = self.columns(); + (0..N).map(|k| u128::from(c[k] >> j & 1) << k).sum() + } + /// Sum over all 64 lanes. + fn total(&self) -> u128 { + let c = self.columns(); + (0..N).map(|k| u128::from(c[k].count_ones()) << k).sum() + } +} + +/// F: plain domino ripple. +#[derive(Clone)] +pub struct Eager { + col: [u64; N], +} + +impl Counter for Eager { + fn new() -> Self { + Self { col: [0; N] } + } + #[inline(always)] + fn add_at(&mut self, mut k: usize, mut m: u64) -> Result<(), Overflow> { + while m != 0 { + if k == N { + return Err(Overflow(m)); + } + let (next, carry) = flip(self.col[k], m); + self.col[k] = next; + m = carry; + k += 1; + } + Ok(()) + } + fn columns(&self) -> [u64; N] { + self.col + } + fn normalize(&mut self) {} +} + +/// P: an EMPTY bitmap (bit `k` set: column `k` is logically zero whatever the +/// word holds). Adding into an EMPTY column assigns the mask. +#[derive(Clone)] +pub struct Subst { + col: [u64; N], + empty: u128, +} + +impl Counter for Subst { + fn new() -> Self { + // The words are never read while EMPTY, so they need no clearing; + // they are zero here only because Rust requires initialized memory. + Self { + col: [0; N], + empty: u128::MAX, + } + } + #[inline(always)] + fn add_at(&mut self, mut k: usize, mut m: u64) -> Result<(), Overflow> { + while m != 0 { + if k == N { + return Err(Overflow(m)); + } + if self.empty >> k & 1 == 1 { + self.col[k] = m; + self.empty &= !(1u128 << k); + return Ok(()); + } + let (next, carry) = flip(self.col[k], m); + self.col[k] = next; + m = carry; + k += 1; + } + Ok(()) + } + fn columns(&self) -> [u64; N] { + core::array::from_fn(|k| if self.empty >> k & 1 == 1 { 0 } else { self.col[k] }) + } + fn normalize(&mut self) { + self.col = self.columns(); + self.empty = 0; + } +} + +impl Subst { + /// A zero counter over words that still hold `stale` (a reused buffer): + /// every column is EMPTY, so the words are never read before written. + pub fn stale(stale: u64) -> Self { + Self { + col: [stale; N], + empty: u128::MAX, + } + } +} + +/// L: a FULL bitmap (bit `k` set: column `k` is logically all ones). An +/// all-ones carry into a FULL column writes zero without reading the word; any +/// other carry materializes `!m`. A column becomes FULL when an update leaves +/// it all ones. +#[derive(Clone)] +pub struct LazyFull { + col: [u64; N], + full: u128, +} + +impl Counter for LazyFull { + fn new() -> Self { + Self { col: [0; N], full: 0 } + } + #[inline(always)] + fn add_at(&mut self, mut k: usize, mut m: u64) -> Result<(), Overflow> { + while m != 0 { + if k == N { + return Err(Overflow(m)); + } + let bit = 1u128 << k; + if self.full & bit != 0 { + // Logical value is all ones: the result is !m, the carry is m. + self.col[k] = !m; + self.full &= !bit; + } else { + let (next, carry) = flip(self.col[k], m); + self.col[k] = next; + if next == u64::MAX { + self.full |= bit; + } + m = carry; + } + k += 1; + } + Ok(()) + } + fn columns(&self) -> [u64; N] { + core::array::from_fn(|k| if self.full >> k & 1 == 1 { u64::MAX } else { self.col[k] }) + } + fn normalize(&mut self) { + self.col = self.columns(); + self.full = 0; + } +} + +/// PL: EMPTY and FULL bitmaps (never both set for one column). An all-ones +/// carry passes a FULL column as a state change FULL -> EMPTY; an EMPTY column +/// takes the mask by assignment (and turns FULL if the mask is all ones). +#[derive(Clone)] +pub struct SubstLazy { + col: [u64; N], + empty: u128, + full: u128, +} + +impl Counter for SubstLazy { + fn new() -> Self { + Self { + col: [0; N], + empty: u128::MAX, + full: 0, + } + } + #[inline(always)] + fn add_at(&mut self, mut k: usize, mut m: u64) -> Result<(), Overflow> { + while m != 0 { + if k == N { + return Err(Overflow(m)); + } + let bit = 1u128 << k; + if self.empty & bit != 0 { + self.empty &= !bit; + if m == u64::MAX { + self.full |= bit; + } else { + self.col[k] = m; + } + return Ok(()); + } + if self.full & bit != 0 { + self.full &= !bit; + if m == u64::MAX { + // FULL + all ones = EMPTY, carry all ones: state only. + self.empty |= bit; + } else { + self.col[k] = !m; + } + } else { + let (next, carry) = flip(self.col[k], m); + self.col[k] = next; + if next == u64::MAX { + // Same FULL policy as L: an ordinary ripple that leaves the + // column all ones marks it, so a later all-ones carry passes + // it as a state change. + self.full |= bit; + } + m = carry; + } + k += 1; + } + Ok(()) + } + fn columns(&self) -> [u64; N] { + core::array::from_fn(|k| { + if self.empty >> k & 1 == 1 { + 0 + } else if self.full >> k & 1 == 1 { + u64::MAX + } else { + self.col[k] + } + }) + } + fn normalize(&mut self) { + self.col = self.columns(); + self.empty = 0; + self.full = 0; + } +} + +impl LazyFull { + /// Bit `k` set: column `k` is held as FULL (logically all ones). + pub fn full_bits(&self) -> u128 { + self.full + } +} + +impl SubstLazy { + /// Bit `k` set: column `k` is held as FULL (logically all ones). + pub fn full_bits(&self) -> u128 { + self.full + } + /// Bit `k` set: column `k` is held as EMPTY (logically zero). + pub fn empty_bits(&self) -> u128 { + self.empty + } + + /// As [`Subst::stale`]: logically zero, words holding `stale`. + pub fn stale(stale: u64) -> Self { + Self { + col: [stale; N], + empty: u128::MAX, + full: 0, + } + } +} + +/// Native baseline: 64 per-lane integers. +#[derive(Clone)] +pub struct Native { + pub lanes: [u64; 64], +} + +impl Native { + pub fn new() -> Self { + Self { lanes: [0; 64] } + } + /// Adds `2^k` to every lane set in `m` (`k < 64`). + #[inline(always)] + pub fn add_at(&mut self, k: usize, m: u64) { + for (j, l) in self.lanes.iter_mut().enumerate() { + *l = l.wrapping_add((m >> j & 1) << k); + } + } +} + +// ── Two significance planes in 128 bits ───────────────────────────────────── + +/// R1: two adjacent significance planes as two words. `lo` is the 2^0 plane, +/// `hi` the 2^1 plane, of the same 64 lanes (64 independent 2-bit counters). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct Pair128 { + pub lo: u64, + pub hi: u64, +} + +/// Adds 1 to every lane set in `m`; returns the carry out of the 2-bit lanes. +#[inline(always)] +pub fn pair_add(p: Pair128, m: u64) -> (Pair128, u64) { + let (lo, c) = flip(p.lo, m); + let (hi, out) = flip(p.hi, c); + (Pair128 { lo, hi }, out) +} + +/// R2: the same two planes packed into one `u128` (low half = 2^0 plane). +#[inline(always)] +pub fn packed_add(x: u128, m: u64) -> (u128, u64) { + let (lo, hi) = (x as u64, (x >> 64) as u64); + let (lo, c) = flip(lo, m); + let (hi, out) = flip(hi, c); + (u128::from(hi) << 64 | u128::from(lo), out) +} + +/// R1 <-> R2. +pub fn pack(p: Pair128) -> u128 { + u128::from(p.hi) << 64 | u128::from(p.lo) +} +pub fn unpack(x: u128) -> Pair128 { + Pair128 { + lo: x as u64, + hi: (x >> 64) as u64, + } +} + +// ── The 64K boundary: low 16 bits wrap, the carry enters the bit-16 column ── + +use crate::carry16::carry_after; +use ndarray::simd::U16x32; + +/// Carry-out lanes of two 32-lane halves as one 64-lane mask (lane `i` of +/// the low half is bit `i`, lane `i` of the high half bit `32 + i`). +#[inline(always)] +pub fn carry_mask(c_lo: U16x32, c_hi: U16x32) -> u64 { + let one = U16x32::splat(1); + u64::from(c_lo.cmpeq_mask(one)) | u64::from(c_hi.cmpeq_mask(one)) << 32 +} + +/// Accumulates `u16` products per lane as `low16 + 2^16 * high`, with `high` +/// held by a bit-sliced counter `C` whose column 0 is the bit-16 plane. +/// `normalize_each` normalizes the counter after every block (eager). +pub fn sum_bitsliced, const N: usize>(prods: &[[U16x32; 2]], normalize_each: bool) -> u128 { + let mut low = [U16x32::splat(0); 2]; + let mut high = C::new(); + for p in prods { + let s0 = low[0] + p[0]; + let s1 = low[1] + p[1]; + let m = carry_mask(carry_after(low[0], p[0], s0), carry_after(low[1], p[1], s1)); + low = [s0, s1]; + high.add_at(0, m).expect("bit-sliced high counter overflow"); + if normalize_each { + high.normalize(); + } + } + let lows: u128 = low.iter().flat_map(|v| v.to_array()).map(u128::from).sum(); + lows + (high.total() << 16) +} + +/// D (#389): the same low word, `high` a `u16` count per lane (`high += carry`), +/// flushed every 65,535 blocks. +pub fn sum_d(prods: &[[U16x32; 2]]) -> u128 { + let mut col = [crate::carry16::Col16::zero(); 2]; + let mut total = 0u128; + for (n, p) in prods.iter().enumerate() { + col[0].add::(p[0]); + col[1].add::(p[1]); + if n % crate::carry16::D_FLUSH_BLOCKS == crate::carry16::D_FLUSH_BLOCKS - 1 { + total += (col[0].total() + col[1].total()) as u128; + col = [crate::carry16::Col16::zero(); 2]; + } + } + total + (col[0].total() + col[1].total()) as u128 +} + +/// Native: one `u64` accumulator per lane, plain Rust (the compiler picks the +/// instructions). +pub fn sum_native(prods: &[[u16; 64]]) -> u128 { + let mut acc = [0u64; 64]; + for p in prods { + for (a, &x) in acc.iter_mut().zip(p) { + *a += u64::from(x); + } + } + acc.iter().map(|&a| u128::from(a)).sum() +} diff --git a/examples/u9_domino/main.rs b/examples/u9_domino/main.rs new file mode 100644 index 00000000..95b3903d --- /dev/null +++ b/examples/u9_domino/main.rs @@ -0,0 +1,357 @@ +//! D-DOMINO-MOV-0 probe: bit flip, positional substitution and lazy FULL +//! forwarding (`domino.rs`) against native integer addition and #389's D. +//! Every row is checked against an exact reference before it is timed. +//! `cargo run --release --example u9_domino --features std`. + +#[path = "../u9_head_to_head/kernels.rs"] +mod kernels; +#[path = "../u9_abacus/carry16.rs"] +mod carry16; +mod domino; + +use domino::*; +use ndarray::simd::U16x32; +use std::hint::black_box; +use std::time::Instant; + +const N: usize = 72; + +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 { + "other" + } +} + +fn xorshift(s: &mut u64) -> u64 { + *s ^= *s << 13; + *s ^= *s >> 7; + *s ^= *s << 17; + *s +} + +/// Median ns per unit over `runs`. +fn time(name: &str, units: usize, runs: usize, unit: &str, base: Option, f: &dyn Fn() -> u128) -> f64 { + black_box(f()); + let mut t: Vec = (0..runs) + .map(|_| { + let t0 = Instant::now(); + black_box(f()); + t0.elapsed().as_nanos() as f64 / units as f64 + }) + .collect(); + t.sort_by(f64::total_cmp); + let med = t[runs / 2]; + let rel = base + .map(|b| format!(" {:6.2}x", med / b)) + .unwrap_or_default(); + println!("{name:<58} {med:9.3} ns/{unit}{rel}"); + med +} + +/// Section A: u16 products summed exactly per lane across the 64K boundary. +fn boundary(label: &str, blocks: usize, runs: usize, dense: bool, seed: u64) { + let mut s = seed; + let raw: Vec<[u16; 64]> = (0..blocks) + .map(|_| { + core::array::from_fn(|_| { + let a = xorshift(&mut s); + if dense { + 65_025 - (a & 1) as u16 + } else { + // Small products: a carry about once per 260 blocks per lane. + (a & 0xFF) as u16 + } + }) + }) + .collect(); + let prods: Vec<[U16x32; 2]> = raw + .iter() + .map(|p| [U16x32::from_slice(&p[..32]), U16x32::from_slice(&p[32..])]) + .collect(); + let want: u128 = raw.iter().flatten().map(|&x| u128::from(x)).sum(); + let native = || sum_native(&raw); + let d = || sum_d(&prods); + let f = || sum_bitsliced::, N>(&prods, false); + let p = || sum_bitsliced::, N>(&prods, false); + let l = || sum_bitsliced::, N>(&prods, false); + let pl = || sum_bitsliced::, N>(&prods, false); + let p_each = || sum_bitsliced::, N>(&prods, true); + let pl_each = || sum_bitsliced::, N>(&prods, true); + for (name, g) in [ + ("native", &native as &dyn Fn() -> u128), + ("D", &d), + ("F", &f), + ("P", &p), + ("L", &l), + ("PL", &pl), + ("P each", &p_each), + ("PL each", &pl_each), + ] { + assert_eq!(g(), want, "{name}"); + } + println!( + "-- A, {label}: {blocks} blocks x 64 u16 products, {}", + if dense { "dense carries" } else { "sparse carries" } + ); + let b = time("native: u64 per lane, plain Rust adds", blocks, runs, "64", None, &native); + time("D (#389): wrapping low + carry rule + high += carry", blocks, runs, "64", Some(b), &d); + time("F: low + carry mask -> bit-16 plane, eager flip ripple", blocks, runs, "64", Some(b), &f); + time("P: F + EMPTY bitmap, assignment into empty planes", blocks, runs, "64", Some(b), &p); + time("L: F + FULL bitmap, all-ones carry writes zero unread", blocks, runs, "64", Some(b), &l); + time("PL: both bitmaps, all-ones carry through FULL = state", blocks, runs, "64", Some(b), &pl); + time("P, normalized every block", blocks, runs, "64", Some(b), &p_each); + time("PL, normalized every block", blocks, runs, "64", Some(b), &pl_each); +} + +/// Section B: the counters alone, one add per step, against per-lane integers. +fn counters(label: &str, adds: &[(usize, u64)], runs: usize) { + let mut want = [0u128; 64]; + for &(k, m) in adds { + for (j, w) in want.iter_mut().enumerate() { + *w += u128::from(m >> j & 1) << k; + } + } + let want_total: u128 = want.iter().sum(); + fn run>(adds: &[(usize, u64)]) -> u128 { + let mut c = C::new(); + for &(k, m) in adds { + c.add_at(k, black_box(m)).unwrap(); + } + c.total() + } + let native = || { + let mut n = Native::new(); + for &(k, m) in adds { + n.add_at(k, black_box(m)); + } + n.lanes.iter().map(|&x| u128::from(x)).sum() + }; + let (f, p, l, pl) = ( + || run::>(adds), + || run::>(adds), + || run::>(adds), + || run::>(adds), + ); + let fits = adds.iter().all(|&(k, _)| k < 48); + if fits { + assert_eq!(native(), want_total, "native"); + } + for (name, g) in [("F", &f as &dyn Fn() -> u128), ("P", &p), ("L", &l), ("PL", &pl)] { + assert_eq!(g(), want_total, "{name}"); + } + println!("-- B, {label}: {} adds", adds.len()); + let b = if fits { + Some(time("native: 64 lanes += bit << k", adds.len(), runs, "add", None, &native)) + } else { + None + }; + time("F eager", adds.len(), runs, "add", b, &f); + time("P EMPTY substitution", adds.len(), runs, "add", b, &p); + time("L lazy FULL", adds.len(), runs, "add", b, &l); + time("PL both", adds.len(), runs, "add", b, &pl); +} + +/// Section B': one add whose carry runs through `len` FULL columns. Each rep +/// clones a preset counter (opaque to the compiler) and adds once; a +/// clone-only row per variant gives the cost to subtract. Correctness is +/// checked outside the timed loop. +fn chains(len: usize, runs: usize) { + const REPS: usize = 1000; + fn preset>(len: usize) -> C { + let mut c = C::new(); + for k in 0..len { + c.add_at(k, u64::MAX).unwrap(); + } + c + } + fn check>(len: usize) { + let mut c: C = preset(len); + c.add_at(0, u64::MAX).unwrap(); + assert_eq!(c.total(), 64u128 << len, "chain {len}"); + } + fn add_rep + Clone>(base: &C) -> u128 { + for _ in 0..REPS { + let mut c = black_box(base.clone()); + let _ = c.add_at(0, black_box(u64::MAX)); + black_box(&c); + } + 0 + } + fn clone_rep + Clone>(base: &C) -> u128 { + for _ in 0..REPS { + let c = black_box(base.clone()); + black_box(&c); + } + 0 + } + fn constant_rep>(len: usize) -> u128 { + let mut acc = 0u128; + for _ in 0..REPS { + let mut c: C = preset(len); + let _ = c.add_at(0, u64::MAX); + acc += c.total() >> len; + } + acc + } + check::>(len); + check::>(len); + check::>(len); + check::>(len); + assert_eq!(constant_rep::>(len), 64 * REPS as u128); + assert_eq!(constant_rep::>(len), 64 * REPS as u128); + println!("-- B', one add carrying through {len} FULL columns (per rep: clone + add)"); + let (f, p, l, pl): (Eager, Subst, LazyFull, SubstLazy) = + (preset(len), preset(len), preset(len), preset(len)); + for (name, g) in [ + ("F clone only", &(|| clone_rep(&f)) as &dyn Fn() -> u128), + ("F clone + add", &|| add_rep(&f)), + ("P clone only", &|| clone_rep(&p)), + ("P clone + add", &|| add_rep(&p)), + ("L clone only", &|| clone_rep(&l)), + ("L clone + add", &|| add_rep(&l)), + ("PL clone only", &|| clone_rep(&pl)), + ("PL clone + add", &|| add_rep(&pl)), + ("F constant state: build + add + total", &|| constant_rep::>(len)), + ("PL constant state: build + add + total", &|| constant_rep::>(len)), + ] { + time(name, REPS, runs, "rep", None, g); + } +} + +/// Section C: R1 (two words) against R2 (one u128), same planes. +fn planes128(runs: usize) { + let mut s = 99u64; + let masks: Vec = (0..4096).map(|_| xorshift(&mut s)).collect(); + let r1 = || { + let (mut p, mut out) = (Pair128 { lo: 0, hi: 0 }, 0u64); + for &m in &masks { + let (np, o) = pair_add(p, m); + p = np; + out ^= o; + } + pack(p) ^ u128::from(out) + }; + let r2 = || { + let (mut x, mut out) = (0u128, 0u64); + for &m in &masks { + let (nx, o) = packed_add(x, m); + x = nx; + out ^= o; + } + x ^ u128::from(out) + }; + assert_eq!(r1(), r2()); + println!("-- C, two planes in 128 bits: {} adds", masks.len()); + let b = time("R1: two u64 planes", masks.len(), runs, "add", None, &r1); + time("R2: one u128 carrier", masks.len(), runs, "add", Some(b), &r2); +} + +// Out-of-line copies for the disassembly. +#[inline(never)] +pub fn asm_flip(a: u64, b: u64) -> (u64, u64) { + flip(a, b) +} +#[inline(never)] +pub fn asm_pair_add(p: Pair128, m: u64) -> (Pair128, u64) { + pair_add(p, m) +} +#[inline(never)] +pub fn asm_packed_add(x: u128, m: u64) -> (u128, u64) { + packed_add(x, m) +} +#[inline(never)] +pub fn asm_eager_add(c: &mut Eager, m: u64) -> Result<(), Overflow> { + c.add_at(0, m) +} +#[inline(never)] +pub fn asm_subst_add(c: &mut Subst, m: u64) -> Result<(), Overflow> { + c.add_at(0, m) +} +#[inline(never)] +pub fn asm_lazy_add(c: &mut LazyFull, m: u64) -> Result<(), Overflow> { + c.add_at(0, m) +} +#[inline(never)] +pub fn asm_pl_add(c: &mut SubstLazy, m: u64) -> Result<(), Overflow> { + c.add_at(0, m) +} +#[inline(never)] +pub fn asm_native_add(c: &mut Native, m: u64) { + c.add_at(0, m) +} +#[inline(never)] +pub fn asm_boundary_f(low: &mut [U16x32; 2], hi: &mut Eager, p: &[U16x32; 2]) { + let s0 = low[0] + p[0]; + let s1 = low[1] + p[1]; + let m = carry_mask(carry16::carry_after(low[0], p[0], s0), carry16::carry_after(low[1], p[1], s1)); + *low = [s0, s1]; + hi.add_at(0, m).unwrap(); +} +#[inline(never)] +pub fn asm_boundary_d(col: &mut [carry16::Col16; 2], p: &[U16x32; 2]) { + col[0].add::(p[0]); + col[1].add::(p[1]); +} + +fn main() { + println!("tier: {} (target-cpu as compiled)", tier()); + { + let m = black_box(0x5555u64); + let mut e = Eager::::new(); + let mut su = Subst::::new(); + let mut la = LazyFull::::new(); + let mut pl = SubstLazy::::new(); + let mut na = Native::new(); + black_box(( + asm_eager_add(&mut e, m).is_ok(), + asm_subst_add(&mut su, m).is_ok(), + asm_lazy_add(&mut la, m).is_ok(), + asm_pl_add(&mut pl, m).is_ok(), + )); + asm_native_add(&mut na, m); + let mut low = [U16x32::splat(black_box(9)); 2]; + asm_boundary_f(&mut low, &mut e, &[U16x32::splat(black_box(7)); 2]); + let mut col = [carry16::Col16::zero(); 2]; + asm_boundary_d(&mut col, &[U16x32::splat(black_box(7)); 2]); + black_box(( + asm_flip(m, m), + asm_pair_add(black_box(Pair128 { lo: !m, hi: m }), m), + asm_packed_add(black_box(u128::from(m) << 64 | u128::from(!m)), m), + na.lanes[0], + low[0].to_array()[0], + col[0].total(), + e.total(), + )); + } + if std::env::args().any(|a| a == "quick") { + return; + } + boundary("L1-resident", 64, 2001, true, 1); + boundary("L1-resident", 64, 2001, false, 2); + boundary("streamed", 1 << 20, 5, true, 3); + boundary("streamed", 1 << 20, 5, false, 4); + + let mut s = 5u64; + let dense: Vec<(usize, u64)> = (0..4096).map(|_| (0, u64::MAX)).collect(); + let random: Vec<(usize, u64)> = (0..4096).map(|_| (0, xorshift(&mut s))).collect(); + let sparse: Vec<(usize, u64)> = (0..4096) + .map(|_| (0, 1 << (xorshift(&mut s) % 64))) + .collect(); + let weighted: Vec<(usize, u64)> = (0..4096) + .map(|t| ([0, 16, 17, 32][t % 4], xorshift(&mut s))) + .collect(); + counters("dense: all-ones at 2^0 (binary counting, long FULL runs)", &dense, 2001); + counters("random masks at 2^0", &random, 2001); + counters("sparse: one lane at 2^0", &sparse, 2001); + counters("weighted: random masks entering at 2^0, 2^16, 2^17, 2^32", &weighted, 2001); + for len in [1, 8, 32, 64] { + chains(len, 201); + } + planes128(2001); +} diff --git a/tests/u9_domino.rs b/tests/u9_domino.rs new file mode 100644 index 00000000..8725e299 --- /dev/null +++ b/tests/u9_domino.rs @@ -0,0 +1,277 @@ +#![cfg(feature = "std")] +//! D-DOMINO-MOV-0 exactness (`examples/u9_domino/domino.rs`) against +//! independent per-lane `u128` arithmetic. + +#[path = "../examples/u9_head_to_head/kernels.rs"] +mod kernels; +#[path = "../examples/u9_abacus/carry16.rs"] +mod carry16; +#[path = "../examples/u9_domino/domino.rs"] +mod domino; + +use domino::*; +use ndarray::simd::U16x32; + +const N: usize = 96; + +fn xorshift(s: &mut u64) -> u64 { + *s ^= *s << 13; + *s ^= *s >> 7; + *s ^= *s << 17; + *s +} + +fn full_tracked(name: &str, step: usize, cols: &[u64; N], full: u128) { + for (k, &c) in cols.iter().enumerate() { + assert_eq!( + c == u64::MAX, + full >> k & 1 == 1, + "{name} step {step}: column {k} is {c:#x} but its FULL bit is {}", + full >> k & 1 + ); + } +} + +/// Applies `adds` to all four counters, checking every lane against a `u128` +/// oracle after each add (observations between additions). +/// P and PL also run over stale words (`STALE` under every EMPTY column), so +/// a lost EMPTY transition or a write skipped into a non-empty column shows. +fn check_all(adds: &[(usize, u64)]) -> [[u64; N]; 6] { + const STALE: u64 = 0xDEAD_BEEF_F00D_CAFE; + let (mut f, mut p) = (Eager::::new(), Subst::::new()); + let (mut l, mut pl) = (LazyFull::::new(), SubstLazy::::new()); + let (mut ps, mut pls) = (Subst::::stale(STALE), SubstLazy::::stale(STALE)); + let mut oracle = [0u128; 64]; + for (step, &(k, m)) in adds.iter().enumerate() { + f.add_at(k, m).unwrap(); + p.add_at(k, m).unwrap(); + l.add_at(k, m).unwrap(); + pl.add_at(k, m).unwrap(); + ps.add_at(k, m).unwrap(); + pls.add_at(k, m).unwrap(); + for (j, o) in oracle.iter_mut().enumerate() { + *o = o.checked_add(u128::from(m >> j & 1) << k).unwrap(); + } + // L and PL hold every logically all-ones column as FULL, whichever + // path produced it (the FULL -> EMPTY state-only carry depends on it). + full_tracked("L", step, &l.columns(), l.full_bits()); + full_tracked("PL", step, &pl.columns(), pl.full_bits()); + full_tracked("PL stale", step, &pls.columns(), pls.full_bits()); + for j in [0usize, 1, 31, 32, 62, 63] { + let want = oracle[j]; + assert_eq!(f.lane(j), want, "F step {step} lane {j}"); + assert_eq!(p.lane(j), want, "P step {step} lane {j}"); + assert_eq!(l.lane(j), want, "L step {step} lane {j}"); + assert_eq!(pl.lane(j), want, "PL step {step} lane {j}"); + assert_eq!(ps.lane(j), want, "P stale step {step} lane {j}"); + assert_eq!(pls.lane(j), want, "PL stale step {step} lane {j}"); + } + } + for (j, &o) in oracle.iter().enumerate() { + assert_eq!(pl.lane(j), o, "PL final lane {j}"); + assert_eq!(f.lane(j), o, "F final lane {j}"); + } + [f.columns(), p.columns(), l.columns(), pl.columns(), ps.columns(), pls.columns()] +} + +/// The bit flip on every input pair, in every lane and in single lanes. +#[test] +fn flip_is_a_half_adder_in_every_lane() { + for (b, c, nb, out) in [(0u64, 0u64, 0u64, 0u64), (0, 1, 1, 0), (1, 0, 1, 0), (1, 1, 0, 1)] { + let fill = |x: u64| if x == 1 { u64::MAX } else { 0 }; + assert_eq!(flip(fill(b), fill(c)), (fill(nb), fill(out))); + for j in [0, 17, 63] { + assert_eq!(flip(b << j, c << j), (nb << j, out << j)); + } + } + let mut s = 7u64; + for _ in 0..1000 { + let (b, c) = (xorshift(&mut s), xorshift(&mut s)); + let (nb, out) = flip(b, c); + for j in 0..64 { + assert_eq!((b >> j & 1) + (c >> j & 1), (nb >> j & 1) + 2 * (out >> j & 1)); + } + } +} + +/// Random masks (zero, full, one lane, random) at random significances, +/// including 16, 63, 64 and beyond, in two different orders. +#[test] +fn random_adds_match_the_oracle_in_any_order() { + let mut s = 0x9E37_79B9_7F4A_7C15u64; + let mut adds = Vec::new(); + for t in 0..3000u64 { + let r = xorshift(&mut s); + let m = match t % 5 { + 0 => 0, + 1 => u64::MAX, + 2 => 1 << (r % 64), + _ => r, + }; + let k = [0usize, 1, 15, 16, 17, 31, 63, 64, 65][(r >> 20) as usize % 9]; + adds.push((k, m)); + } + let fwd = check_all(&adds); + adds.reverse(); + let rev = check_all(&adds); + for v in 0..6 { + assert_eq!(fwd[v], rev[v], "variant {v}: order changed the result"); + assert_eq!(fwd[v], fwd[0], "variant {v} disagrees with F"); + } +} + +/// Carry through 1, 2, 8, 16, 32 and 64 consecutive FULL columns, for all +/// lanes and for one lane, and two successive carries into the same column. +#[test] +fn carries_run_through_long_full_chains() { + for n in [1usize, 2, 8, 16, 32, 64] { + for m in [u64::MAX, 1 << 37] { + let mut adds: Vec<(usize, u64)> = (0..n).map(|k| (k, m)).collect(); + adds.push((0, m)); + adds.push((0, m)); + let cols = check_all(&adds); + // Value 2^n + 1 in the lanes of m, 0 elsewhere. + for c in cols { + assert_eq!(c[0], m); + assert_eq!(c[n], m, "chain {n}"); + for (k, &w) in c.iter().enumerate() { + if k != 0 && k != n { + assert_eq!(w, 0, "chain {n} column {k}"); + } + } + } + } + } +} + +/// The 64K boundary: low16 + p crosses 2^16 exactly when the sum is at least +/// 65,536; the carry flips bit 16 and moves on to bit 17 when bit 16 is set. +#[test] +fn the_64k_boundary_flips_bit_16_then_carries_into_bit_17() { + for low0 in [0u16, 1, 510, 511, 512, 65_000, 65_534, 65_535] { + for p in [0u16, 1, 511, 65_024, 65_025] { + for bit16_set in [false, true] { + let mut lows = [U16x32::splat(low0); 2]; + let mut hi = Eager::::new(); + if bit16_set { + hi.add_at(0, u64::MAX).unwrap(); + } + let pv = [U16x32::splat(p); 2]; + let s = [lows[0] + pv[0], lows[1] + pv[1]]; + let m = + carry_mask(carry16::carry_after(lows[0], pv[0], s[0]), carry16::carry_after(lows[1], pv[1], s[1])); + lows = s; + hi.add_at(0, m).unwrap(); + let sum = u32::from(low0) + u32::from(p) + if bit16_set { 65_536 } else { 0 }; + for j in [0usize, 31, 32, 63] { + let low = lows[j / 32].to_array()[j % 32]; + assert_eq!(u32::from(low) + ((hi.lane(j) as u32) << 16), sum, "low {low0} p {p} bit16 {bit16_set}"); + } + let crossed = u32::from(low0) + u32::from(p) >= 65_536; + let c = hi.columns(); + assert_eq!(c[0] == u64::MAX, crossed != bit16_set, "bit 16"); + assert_eq!(c[1] == u64::MAX, crossed && bit16_set, "bit 17"); + } + } + } +} + +/// u16 products summed through every high representation against `u128`. +#[test] +fn product_sums_match_native_and_d() { + let mut s = 0x0123_4567_89AB_CDEFu64; + for (blocks, dense) in [(1usize, false), (300, false), (300, true), (70_000, true)] { + let raw: Vec<[u16; 64]> = (0..blocks) + .map(|_| { + core::array::from_fn(|_| { + if dense { + 65_025 + } else { + let a = xorshift(&mut s); + ((a & 0xFF) * (a >> 8 & 0xFF)) as u16 + } + }) + }) + .collect(); + let prods: Vec<[U16x32; 2]> = raw + .iter() + .map(|p| [U16x32::from_slice(&p[..32]), U16x32::from_slice(&p[32..])]) + .collect(); + let want: u128 = raw.iter().flatten().map(|&x| u128::from(x)).sum(); + assert_eq!(sum_native(&raw), want); + assert_eq!(sum_d(&prods), want, "D {blocks}"); + for each in [false, true] { + assert_eq!(sum_bitsliced::, N>(&prods, each), want, "F {blocks} {each}"); + assert_eq!(sum_bitsliced::, N>(&prods, each), want, "P {blocks} {each}"); + assert_eq!(sum_bitsliced::, N>(&prods, each), want, "L {blocks} {each}"); + assert_eq!(sum_bitsliced::, N>(&prods, each), want, "PL {blocks} {each}"); + } + } +} + +/// A carry past the last column is refused, in every variant. +#[test] +fn overflow_past_the_last_column_is_refused() { + fn run>() { + let mut c = C::new(); + for k in 0..4 { + c.add_at(k, 0b1010).unwrap(); + } + assert_eq!(c.add_at(0, 0b0011), Err(Overflow(0b0010))); + assert_eq!(c.add_at(4, 1), Err(Overflow(1))); + assert_eq!(c.add_at(4, 0), Ok(())); + } + run::>(); + run::>(); + run::>(); + run::>(); +} + +/// R1 (two words) and R2 (one u128) hold the same two planes and agree on +/// every add; the roundtrip is the identity. +#[test] +fn pair_and_packed_agree_and_roundtrip() { + let mut s = 42u64; + let mut p = Pair128 { lo: 0, hi: 0 }; + let mut x = 0u128; + for _ in 0..10_000 { + let m = xorshift(&mut s); + let (np, o1) = pair_add(p, m); + let (nx, o2) = packed_add(x, m); + assert_eq!((pack(np), o1), (nx, o2)); + assert_eq!(unpack(nx), np); + // 64 independent 2-bit counters, not one 128-bit integer. + for j in 0..64 { + let before = (p.lo >> j & 1) + 2 * (p.hi >> j & 1); + let after = (np.lo >> j & 1) + 2 * (np.hi >> j & 1) + 4 * (o1 >> j & 1); + assert_eq!(after, before + (m >> j & 1), "lane {j}"); + } + p = np; + x = nx; + } +} + +/// A column made all ones by an ordinary XOR (complementary masks at one +/// significance) is FULL in PL as in L, and the next all-ones carry passes it +/// as FULL -> EMPTY without touching the word. +#[test] +fn xor_made_full_column_is_full_and_carries_by_state() { + let m = 0x0F0F_0F0F_0F0F_0F0Fu64; + let mut l = LazyFull::::new(); + let mut pl = SubstLazy::::new(); + for (k, mask) in [(0, m), (0, !m)] { + l.add_at(k, mask).unwrap(); + pl.add_at(k, mask).unwrap(); + } + assert_eq!(l.full_bits() & 1, 1, "L marks the XOR-made column FULL"); + assert_eq!(pl.full_bits() & 1, 1, "PL marks the XOR-made column FULL"); + pl.add_at(0, u64::MAX).unwrap(); + assert_eq!(pl.full_bits() & 1, 0); + assert_eq!(pl.empty_bits() & 1, 1, "FULL + all ones leaves column 0 EMPTY"); + // Column 1 was EMPTY, so it takes the all-ones carry by assignment and + // is FULL in turn. + assert_eq!(pl.full_bits() >> 1 & 1, 1, "the all-ones carry lands in column 1 as FULL"); + for j in 0..64 { + assert_eq!(pl.lane(j), 2, "lane {j}"); + } +}