From 97d8cdd5d04e22b327fe4c233e58f971fd558ac2 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 12:38:55 +0000 Subject: [PATCH] D-SIMD-MICRO-MACRO-1 step 1: packed microtile oracle (lab fixture) One 128-bit tile carrying 5xU24/I24 (+8 metadata bits) or 14xU9/I9 (+2 spare) isolated fields: wrapping add/sub, per-field carry/borrow/signed overflow, an intra-tile stencil, and an exact unsigned+signed fold. Four realizations (u128, SoA U64x8, widened U32x16, bit planes over full_add_mask) are checked against one unpacked oracle on six SIMD realizations, with disable runs D1-D13. No realization wins every case; report in .claude/knowledge/microtile-probe.md. NEON/wasm U32x16 facade gaps recorded, not fixed. No public type, no backend change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY --- .claude/blackboard.md | 15 + .claude/knowledge/microtile-probe.md | 114 +++++++ Cargo.toml | 5 + examples/microtile/main.rs | 135 ++++++++ examples/microtile/tile.rs | 479 +++++++++++++++++++++++++++ tests/microtile.rs | 220 ++++++++++++ 6 files changed, 968 insertions(+) create mode 100644 .claude/knowledge/microtile-probe.md create mode 100644 examples/microtile/main.rs create mode 100644 examples/microtile/tile.rs create mode 100644 tests/microtile.rs diff --git a/.claude/blackboard.md b/.claude/blackboard.md index 5e091424..257873f1 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -4010,3 +4010,18 @@ laws, ExtentUnsupported, row bounds) and WITNESS in OrderedLaneWitness; no new t - PR A #402 (draft): `half_add_mask` / `full_add_mask` in `simd_masking_ops`, re-exported; vectorized `pow2n_from_int` (vsexp 4x, vserf 2.4x on AVX-512); `U32x16` `Shl`/`Shr` added to AVX2, AVX, NEON, wasm, nightly. - Lab: `Abacus` with transactional refusal in `examples/u9_domino`; cost equal to eager F. PR B and PR C deferred: no consumer, planes already cheap with existing words. Report: `.claude/knowledge/binary-abacus-polyfill-probe.md`. - Loose end: NEON `simd_exp_f32` has ~47 scalar fcmp/csel (pre-existing, select/compare path); AVX2 mask packing for `select` still lane by lane. + +## 2026-10-30 — D-SIMD-MICRO-MACRO-1 step 1: packed microtile oracle (lab only) + +- `examples/microtile/` + `tests/microtile.rs`: 128-bit tiles of `5×U24/I24` (+8 metadata + bits) and `14×U9/I9` (+2 spare), field-isolated add/sub/carry/borrow/signed overflow, + intra-tile stencil, exact unsigned+signed fold; u128, SoA `U64x8`, widened `U32x16`, + bit-plane (`full_add_mask`) realizations, all against one unpacked oracle. +- 9 tests green on AVX-512, AVX2, AVX, v2, NEON (qemu), nightly. Disable runs D1–D13 red + except D6 (redundant mask, removed). +- Cost gate: no single winner. U24x5 → widened for add and fold; U9x14 → packed SoA for + L2-resident add, bit planes for the fold. u128 alone never wins (not vectorized). +- Found: NEON/wasm `U32x16` lacks `&`, `|`, `-`, `!`, `+=`, `from_slice`, `copy_to_slice`. + Recorded, not fixed. +- Open: steps 2–4 (KPEA method binding, Thought-0 delta reuse, interlacing). H13 stays + deferred. Report: `.claude/knowledge/microtile-probe.md`. diff --git a/.claude/knowledge/microtile-probe.md b/.claude/knowledge/microtile-probe.md new file mode 100644 index 00000000..ad765c2f --- /dev/null +++ b/.claude/knowledge/microtile-probe.md @@ -0,0 +1,114 @@ +# D-SIMD-MICRO-MACRO-1, step 1: the packed microtile oracle + +Lab fixture only: `examples/microtile/` (realizations + cost gate) and +`tests/microtile.rs` (exactness). No public type, no backend change, no +production caller. Steps 2+ (binding to a KPEA method, Thought-0 delta reuse, +interlacing) are not started. + +## What the microtile is + +One 128-bit tile carries several isolated arithmetic fields: + +| layout | fields | payload | other bits | +|---|---|---|---| +| `U24x5` | 5 × 24 bit (`U24`/`I24`) | 120 bits | 8 metadata bits, 120..128 | +| `U9x14` | 14 × 9 bit (`U9`/`I9`) | 126 bits | 2 spare bits, 126..128 | + +Operations: field-isolated wrapping add and sub, per-field unsigned carry and +borrow out, per-field signed overflow, an intra-tile stencil +`f[j-1] + 2 f[j] + f[j+1]`, and an exact fold (unsigned and signed per-field +sums over many tiles). Every result keeps the left operand's metadata. + +Isolation is the SWAR law: clear each field's top bit in the addends so a +carry stops at its own top bit, then fix that bit with +`(a ^ b) & HIGH`; for sub, set each top bit of `a` so a borrow stops there. +Signed fold: `Σ sext(f) = Σ f - 2^W · #(top bit set)`, so one pass yields +both sums. + +## Four realizations, one oracle + +The oracle unpacks every field and does the arithmetic in `u64`/`i64`. + +1. **u128**: the tile as one 128-bit integer. +2. **SoA `U64x8`**: 8 tiles per op, the tile split into a low and a high + 64-bit plane. The fields that straddle bit 64 (`U24` field 2, `U9` field + 7) take an explicit carry/borrow out of bit 63, computed with lane-wise + logic only (no shuffle). +3. **Widened `U32x16`**: field-major `u32` lanes. The fold flushes to `u64` + every `floor((2^32-1)/(2^W-1))` steps (256 for `U24`). +4. **Bit planes**: one bit per tile, 64 tiles per word, add by ripple of + `simd::full_add_mask`, fold by `Σ_k 2^k popcount(plane_k)`. + +Exactness: 9 tests over random tiles plus every pairing of edge field values +(0, 1, sign bit - 1, sign bit, max - 1, max), the two straddling fields' +cross-word carry and borrow, 8192 maximal `U24` values in one `u32` lane, and +metadata against an all-ones right operand. Green on AVX-512, AVX2, +AVX-without-AVX2 (qemu SandyBridge), x86-64-v2 (qemu Nehalem), NEON (qemu +aarch64) and nightly `core::simd`. + +Disable runs, each red: top bit not cleared in add (D1) or set in sub (D8), +no carry (D2) or borrow (D7) across bit 64, metadata dropped (D3), signed +fold without the top-bit correction (D4) and its plane form (D9), widened +fold never flushing (D5), carry-out as `a & b` only (D10), straddling field +read from the low word only (D11), fold tail skipped (D12), widened add +unmasked (D13). D6 (an extra mask on the stencil's left shift) stayed green: +`add` already masks its right operand, so the mask was removed rather than +kept as a guard that does not bind. + +## Cost gate + +`cargo run --release --example microtile --features std`, median ns per tile +over 31 runs, three runs (two AVX-512 v4, one AVX2 v3) on a 4-core host. +4096 tiles = 64 KB per operand (L2-resident); 65536 tiles = 1 MB. + +| op | layout | u128 | SoA U64x8 | widened | planes | +|---|---|---|---|---|---| +| add, 4K | U24x5 | 2.00–2.10 | 0.95–1.73 | **0.69–0.71** | 1.87–1.93 | +| add, 4K | U9x14 | 2.03–2.10 | **1.08–1.66** | 1.92–1.97 | 1.83–2.03 | +| add, 64K | U24x5 | 2.29–2.31* | 1.94–3.16 | 2.18–2.35 | **1.88–2.07** | +| add, 64K | U9x14 | 2.04–2.41 | 1.54–2.81 | 6.19–6.39 | **1.87–2.06** | +| fold, 4K | U24x5 | — | 0.60–1.01 | **0.29–0.30** | 0.85–0.90 | +| fold, 4K | U9x14 | — | 2.71–3.60 | 0.76–1.20 | **0.73–1.42** | +| fold, 64K | U24x5 | — | 0.53–0.93 | **0.28–0.31** | 0.75–0.82 | +| fold, 64K | U9x14 | — | 2.81–3.22 | 1.85–2.16 | **0.60–0.76** | + +\* one outlier at 4.19 excluded. + +Storage: packed and planes 16 B/tile; widened 20 B/tile (`U24x5`) and +56 B/tile (`U9x14`). + +Read: + +- **No realization wins everything, so packed is not assumed.** Wide fields + (`U24x5`) favour widened lanes for both add and fold: only 25% more bytes, + and one native add per field. Narrow fields (`U9x14`) favour packed SoA for + cache-resident add and bit planes for the fold; widened costs 3.5× the + bytes and is 3× slower once the data leaves L2. +- **The single u128 is never the fast form.** It runs at about 2 ns/tile + whatever the layout: rustc does not vectorize the two-word add. The packed + win comes from the SoA split, which puts 8 tiles in each register. +- **At 1 MB every add is about 2 ns/tile** except widened `U9x14`: the + operation is bandwidth-bound and the bytes per tile decide. +- **Bit planes earn their place on folds, not adds.** A narrow-field fold is + W popcounts per 64 tiles; their add needs W full adders and is no faster + than packed. +- Unexplained: SoA add is faster on AVX2 (0.95–1.08 at 4K) than on AVX-512 + (1.61–1.73). Not diagnosed. + +## Findings outside the fixture + +- **NEON and wasm `U32x16` lack most of the facade:** no `&`, `|`, `-`, + `!`, `+=`, `from_slice` or `copy_to_slice`; they implement `+`, `^`, the + two shifts and `from_array`/`to_array`. AVX-512, AVX2, scalar and nightly + have the full set. The fixture works around it (mask by `<< k >> k`, load + through `from_array`), so this is recorded, not fixed here. It is a real + parity gap for any consumer that builds on `U32x16`. + +## What this does not show yet + +- No consumer: the gates for binding to a KPEA mask/count method (step 2), + Thought-0 invariant reuse across thought deltas (step 3) and interlacing + (step 4) are open. The fold here is the candidate primitive for step 2. +- Pack/unpack between layouts is untimed; a method that converts would pay it. +- Only x86 was timed; NEON and the other tiers ran under emulation for + exactness only. diff --git a/Cargo.toml b/Cargo.toml index 20a1cb07..15c4ad21 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -90,6 +90,11 @@ name = "u9_domino" path = "examples/u9_domino/main.rs" required-features = ["std"] +[[example]] +name = "microtile" +path = "examples/microtile/main.rs" +required-features = ["std"] + [[example]] name = "cpu_guard_probe" required-features = ["std"] diff --git a/examples/microtile/main.rs b/examples/microtile/main.rs new file mode 100644 index 00000000..975201ab --- /dev/null +++ b/examples/microtile/main.rs @@ -0,0 +1,135 @@ +//! D-SIMD-MICRO-MACRO-1 step 1 cost gate: field-isolated add and exact fold +//! over 64K microtiles, in four realizations. Every realization is checked +//! against the unpacked oracle before it is timed. +//! `cargo run --release --example microtile --features std`. + +// The tests exercise sub and the stencil; the bench times add and fold only. +#[allow(dead_code)] +mod tile; + +use std::hint::black_box; +use std::time::Instant; +use tile::*; + +fn tier() -> &'static str { + if cfg!(target_feature = "avx512bw") { + "x86_64 AVX-512" + } else if cfg!(target_feature = "avx2") { + "x86_64 AVX2" + } else if cfg!(target_arch = "aarch64") { + "aarch64 NEON" + } else { + "other" + } +} + +fn xorshift(s: &mut u64) -> u64 { + *s ^= *s << 13; + *s ^= *s >> 7; + *s ^= *s << 17; + *s +} + +/// Median ns per tile over 31 runs. +fn time(name: &str, tiles: usize, base: Option, f: &mut dyn FnMut() -> u64) -> f64 { + black_box(f()); + let mut t: Vec = (0..31) + .map(|_| { + let t0 = Instant::now(); + black_box(f()); + t0.elapsed().as_nanos() as f64 / tiles as f64 + }) + .collect(); + t.sort_by(f64::total_cmp); + let med = t[15]; + let rel = base + .map(|b| format!(" {:6.2}x", med / b)) + .unwrap_or_default(); + println!(" {name:<44} {med:8.3} ns/tile{rel}"); + med +} + +fn run(label: &str, tiles: usize) { + let mut s = 0x9E37_79B9_7F4A_7C15u64; + let a: Vec = (0..tiles) + .map(|_| (xorshift(&mut s) as u128) << 64 | xorshift(&mut s) as u128) + .collect(); + let b: Vec = (0..tiles) + .map(|_| (xorshift(&mut s) as u128) << 64 | xorshift(&mut s) as u128) + .collect(); + let ((alo, ahi), (blo, bhi)) = (split(&a), split(&b)); + let (sa, sb) = (Soa { lo: &alo, hi: &ahi }, Soa { lo: &blo, hi: &bhi }); + let (wa, wb) = (widen::

(&a), widen::

(&b)); + let (pa, pb) = (to_planes::

(&a), to_planes::

(&b)); + let (n, w) = (P::N as usize, P::W as usize); + + // Exactness before timing. + let want: Vec = a + .iter() + .zip(&b) + .map(|(&x, &y)| oracle_add::

(x, y).0) + .collect(); + let mut out128 = vec![0u128; tiles]; + for t in 0..tiles { + out128[t] = add::

(a[t], b[t]); + } + assert_eq!(out128, want, "{label} u128 add"); + let (mut ol, mut oh) = (vec![0u64; tiles], vec![0u64; tiles]); + soa_add_slices::

(&sa, &sb, &mut ol, &mut oh); + assert!((0..tiles).all(|t| (ol[t] as u128 | (oh[t] as u128) << 64) == want[t]), "{label} soa add"); + let mut ow = vec![vec![0u32; tiles]; n]; + for j in 0..n { + widened_add::

(&wa[j], &wb[j], &mut ow[j]); + assert!((0..tiles).all(|t| ow[j][t] == get::

(want[t], j as u32)), "{label} widened add"); + } + let (fu, fs) = oracle_fold::

(&a); + assert_eq!(soa_fold::

(&sa), (fu.clone(), fs.clone()), "{label} soa fold"); + for j in 0..n { + assert_eq!(widened_fold::

(&wa[j]), fu[j], "{label} widened fold"); + assert_eq!(planes_fold(&pa[j], w), (fu[j], fs[j]), "{label} planes fold"); + } + + println!("{label}: {n} fields x {w} bits, {tiles} tiles, all realizations exact"); + println!(" storage: packed 16 B/tile, widened {} B/tile, planes 16 B/tile", 4 * n); + let base = time("add u128 (one 128-bit word)", tiles, None, &mut || { + for t in 0..tiles { + out128[t] = add::

(a[t], b[t]); + } + out128[tiles - 1] as u64 + }); + time("add SoA U64x8 (8 tiles per op)", tiles, Some(base), &mut || { + soa_add_slices::

(&sa, &sb, &mut ol, &mut oh); + ol[tiles - 1] + }); + time("add widened U32x16 (one field per lane)", tiles, Some(base), &mut || { + for j in 0..n { + widened_add::

(&wa[j], &wb[j], &mut ow[j]); + } + ow[0][tiles - 1] as u64 + }); + let mut op = vec![0u64; w]; + time("add bit planes, full_add_mask ripple", tiles, Some(base), &mut || { + let mut acc = 0; + for j in 0..n { + for (ak, bk) in pa[j].chunks_exact(w).zip(pb[j].chunks_exact(w)) { + acc ^= planes_add(ak, bk, &mut op); + } + } + acc + }); + let base = time("fold SoA U64x8 (unsigned + signed)", tiles, None, &mut || soa_fold::

(&sa).0[0]); + time("fold widened U32x16 (unsigned)", tiles, Some(base), &mut || { + (0..n).map(|j| widened_fold::

(&wa[j])).sum() + }); + time("fold bit planes, popcount (unsigned + signed)", tiles, Some(base), &mut || { + (0..n).map(|j| planes_fold(&pa[j], w).0).sum() + }); +} + +fn main() { + println!("D-SIMD-MICRO-MACRO-1 step 1, {}", tier()); + for tiles in [1 << 12, 1 << 16] { + run::("U24x5", tiles); + run::("U9x14", tiles); + } +} diff --git a/examples/microtile/tile.rs b/examples/microtile/tile.rs new file mode 100644 index 00000000..e9ee895a --- /dev/null +++ b/examples/microtile/tile.rs @@ -0,0 +1,479 @@ +//! D-SIMD-MICRO-MACRO-1, step 1: one 128-bit microtile carrying several +//! isolated arithmetic fields, in four realizations that are checked against +//! one unpacked oracle. +//! +//! Layouts (`Packing`): +//! - `U24x5`: five 24-bit fields in bits 0..120, 8 metadata bits 120..128. +//! - `U9x14`: fourteen 9-bit fields in bits 0..126, 2 spare bits 126..128. +//! +//! Field `j` occupies bits `j*W .. j*W+W`. Its top bit is the sign bit when +//! the field is read as `I24`/`I9`. Bits outside the fields (`META`) are never +//! touched by arithmetic: every result carries the LEFT operand's metadata. +//! +//! Realizations: +//! - `u128` (`add`, `sub`, ...): the whole tile in one 128-bit integer. This is +//! the "one register" form; rustc lowers it to a two-word add/adc pair. +//! - `soa` (`Soa*`): 8 tiles at a time over `ndarray::simd::U64x8`, with the +//! tile split into a low and a high 64-bit plane. Fields that straddle bit 64 +//! (`U24` field 2, `U9` field 7) get their carry/borrow from an explicit +//! 64-bit carry-out, computed with lane-wise logic only. +//! - `widened`: field-major `u32` lanes over `ndarray::simd::U32x16`. +//! - `planes`: bit planes over 64 tiles, added with `ndarray::simd::full_add_mask`. +//! +//! Lab fixture only: no public type, no production caller. + +use ndarray::simd::{full_add_mask, U32x16, U64x8}; + +/// A packed layout of `N` fields of width `W` in 128 bits. +pub trait Packing { + const W: u32; + const N: u32; + /// One field's value mask, `2^W - 1`. + const VMASK: u128 = (1u128 << Self::W) - 1; + /// All field bits. + const FIELDS: u128 = { + let mut m = 0u128; + let mut j = 0; + while j < Self::N { + m |= Self::VMASK << (j * Self::W); + j += 1; + } + m + }; + /// The top (sign) bit of every field. + const HIGH: u128 = { + let mut m = 0u128; + let mut j = 0; + while j < Self::N { + m |= 1u128 << (j * Self::W + Self::W - 1); + j += 1; + } + m + }; + /// Field bits below each field's top bit. + const LOW: u128 = Self::FIELDS & !Self::HIGH; + /// Bits outside the fields: metadata (`U24x5`) or spare (`U9x14`). + const META: u128 = !Self::FIELDS; +} + +/// Five 24-bit fields + 8 metadata bits. +pub struct U24x5; +impl Packing for U24x5 { + const W: u32 = 24; + const N: u32 = 5; +} + +/// Fourteen 9-bit fields + 2 spare bits. +pub struct U9x14; +impl Packing for U9x14 { + const W: u32 = 9; + const N: u32 = 14; +} + +// ── Oracle: unpacked per-field arithmetic ──────────────────────────────── + +/// Field `j` as an unsigned value. +pub fn get(t: u128, j: u32) -> u32 { + ((t >> (j * P::W)) & P::VMASK) as u32 +} + +/// Field `j` sign-extended (`I24` / `I9`). +pub fn get_signed(t: u128, j: u32) -> i32 { + let s = 32 - P::W; + ((get::

(t, j) << s) as i32) >> s +} + +/// Build a tile from field values (each masked to `W` bits) and metadata. +pub fn pack(fields: &[u32], meta: u128) -> u128 { + assert_eq!(fields.len(), P::N as usize); + let mut t = meta & P::META; + for (j, &f) in fields.iter().enumerate() { + t |= (f as u128 & P::VMASK) << (j as u32 * P::W); + } + t +} + +/// Oracle add: per-field wrapping sum, per-field unsigned carry-out and +/// signed overflow, metadata from `a`. +pub fn oracle_add(a: u128, b: u128) -> (u128, u32, u32) { + let m = (1u64 << P::W) - 1; + let (mut fields, mut carry, mut ovf) = (Vec::new(), 0u32, 0u32); + for j in 0..P::N { + let s = get::

(a, j) as u64 + get::

(b, j) as u64; + fields.push((s & m) as u32); + carry |= ((s > m) as u32) << j; + let ss = get_signed::

(a, j) as i64 + get_signed::

(b, j) as i64; + let lim = 1i64 << (P::W - 1); + ovf |= ((ss >= lim || ss < -lim) as u32) << j; + } + (pack::

(&fields, a), carry, ovf) +} + +/// Oracle sub: per-field wrapping difference and unsigned borrow-out. +pub fn oracle_sub(a: u128, b: u128) -> (u128, u32) { + let m = (1u64 << P::W) - 1; + let (mut fields, mut borrow) = (Vec::new(), 0u32); + for j in 0..P::N { + let (x, y) = (get::

(a, j) as u64, get::

(b, j) as u64); + fields.push((x.wrapping_sub(y) & m) as u32); + borrow |= ((x < y) as u32) << j; + } + (pack::

(&fields, a), borrow) +} + +/// Oracle stencil inside one tile: `f[j-1] + 2 f[j] + f[j+1]` mod `2^W`, +/// zero beyond the first and last field, metadata kept. +pub fn oracle_stencil(t: u128) -> u128 { + let n = P::N as i64; + let f = |j: i64| { + if j < 0 || j >= n { + 0u64 + } else { + get::

(t, j as u32) as u64 + } + }; + let m = (1u64 << P::W) - 1; + let fields: Vec = (0..n) + .map(|j| ((f(j - 1) + 2 * f(j) + f(j + 1)) & m) as u32) + .collect(); + pack::

(&fields, t) +} + +/// Oracle fold: exact per-field sums, unsigned and signed. +pub fn oracle_fold(tiles: &[u128]) -> (Vec, Vec) { + let mut u = vec![0u64; P::N as usize]; + let mut s = vec![0i64; P::N as usize]; + for &t in tiles { + for j in 0..P::N { + u[j as usize] += get::

(t, j) as u64; + s[j as usize] += get_signed::

(t, j) as i64; + } + } + (u, s) +} + +/// Collect a per-field flag word (bit `j` = field `j`) from top-bit positions. +pub fn flags(at_high: u128) -> u32 { + (0..P::N).fold(0, |acc, j| acc | ((at_high >> (j * P::W + P::W - 1) & 1) as u32) << j) +} + +// ── Realization 1: one u128 ────────────────────────────────────────────── + +/// Field-isolated wrapping add, metadata from `a`. Carries never leave a +/// field: the top bit of every addend is cleared, so a field's carry lands +/// in its own top bit, which is then fixed by XOR. +#[inline] +pub fn add(a: u128, b: u128) -> u128 { + let s = (a & P::LOW) + (b & P::LOW); + ((s ^ ((a ^ b) & P::HIGH)) & P::FIELDS) | (a & P::META) +} + +/// Per-field unsigned carry-out of `a + b` given its result `r`, at top-bit +/// positions: the majority of `aH`, `bH` and the carry into `H`. +#[inline] +pub fn carry_out(a: u128, b: u128, r: u128) -> u128 { + ((a & b) | ((a | b) & !r)) & P::HIGH +} + +/// Per-field signed overflow of `a + b` (I24/I9), at top-bit positions. +#[inline] +pub fn signed_overflow(a: u128, b: u128, r: u128) -> u128 { + (!(a ^ b) & (a ^ r)) & P::HIGH +} + +/// Field-isolated wrapping sub, metadata from `a`. Setting each field's top +/// bit in `a` stops every borrow at its own field. +#[inline] +pub fn sub(a: u128, b: u128) -> u128 { + let d = (a | P::HIGH) - (b & P::LOW); + ((d ^ ((a ^ !b) & P::HIGH)) & P::FIELDS) | (a & P::META) +} + +/// Per-field unsigned borrow-out of `a - b` given its result `r`. +#[inline] +pub fn borrow_out(a: u128, b: u128, r: u128) -> u128 { + ((!a & b) | (!(a ^ b) & r)) & P::HIGH +} + +/// Intra-tile stencil `f[j-1] + 2 f[j] + f[j+1]`: neighbours are moved by a +/// whole-tile shift of `W` bits and masked back to the field bits. +#[inline] +pub fn stencil(t: u128) -> u128 { + let f = t & P::FIELDS; + // f[j-1] into field j. The last field lands outside the fields; `add` + // masks its right operand to the field bits, so no mask is needed here. + let left = f << P::W; + let right = f >> P::W; // f[j+1] into field j + let two = add::

(f, f); + let r = add::

(add::

(two, left), right); + (r & P::FIELDS) | (t & P::META) +} + +// ── Realization 2: SoA low/high planes over U64x8 ──────────────────────── + +/// The `u128` constants split into low and high 64-bit words. +const fn lo(m: u128) -> u64 { + m as u64 +} +const fn hi(m: u128) -> u64 { + (m >> 64) as u64 +} + +/// Carry out of bit 63 of `x + y = s`, as 0 or 1 per lane. +#[inline] +fn carry63(x: U64x8, y: U64x8, s: U64x8) -> U64x8 { + ((x & y) | ((x | y) & !s)) >> U64x8::splat(63) +} + +/// Borrow out of bit 63 of `x - y = d`, as 0 or 1 per lane. +#[inline] +fn borrow63(x: U64x8, y: U64x8, d: U64x8) -> U64x8 { + ((!x & y) | ((!x | y) & d)) >> U64x8::splat(63) +} + +/// Field-isolated add of 8 tiles. Returns the result planes and the carry-out +/// planes at top-bit positions. +#[inline] +pub fn soa_add(alo: U64x8, ahi: U64x8, blo: U64x8, bhi: U64x8) -> [U64x8; 4] { + let (ll, lh) = (U64x8::splat(lo(P::LOW)), U64x8::splat(hi(P::LOW))); + let (hl, hh) = (U64x8::splat(lo(P::HIGH)), U64x8::splat(hi(P::HIGH))); + let (fl, fh) = (U64x8::splat(lo(P::FIELDS)), U64x8::splat(hi(P::FIELDS))); + let (ml, mh) = (U64x8::splat(lo(P::META)), U64x8::splat(hi(P::META))); + let (xl, yl) = (alo & ll, blo & ll); + let sl = xl + yl; + let sh = (ahi & lh) + (bhi & lh) + carry63(xl, yl, sl); + let rl = ((sl ^ ((alo ^ blo) & hl)) & fl) | (alo & ml); + let rh = ((sh ^ ((ahi ^ bhi) & hh)) & fh) | (ahi & mh); + let cl = ((alo & blo) | ((alo | blo) & !rl)) & hl; + let ch = ((ahi & bhi) | ((ahi | bhi) & !rh)) & hh; + [rl, rh, cl, ch] +} + +/// Field-isolated sub of 8 tiles. Returns the result planes. +#[inline] +pub fn soa_sub(alo: U64x8, ahi: U64x8, blo: U64x8, bhi: U64x8) -> [U64x8; 2] { + let (ll, lh) = (U64x8::splat(lo(P::LOW)), U64x8::splat(hi(P::LOW))); + let (hl, hh) = (U64x8::splat(lo(P::HIGH)), U64x8::splat(hi(P::HIGH))); + let (fl, fh) = (U64x8::splat(lo(P::FIELDS)), U64x8::splat(hi(P::FIELDS))); + let (ml, mh) = (U64x8::splat(lo(P::META)), U64x8::splat(hi(P::META))); + let (xl, yl) = (alo | hl, blo & ll); + let dl = xl - yl; + let dh = (ahi | hh) - (bhi & lh) - borrow63(xl, yl, dl); + let rl = ((dl ^ ((alo ^ !blo) & hl)) & fl) | (alo & ml); + let rh = ((dh ^ ((ahi ^ !bhi) & hh)) & fh) | (ahi & mh); + [rl, rh] +} + +/// Borrowed SoA tiles: tile `t` is `lo[t] | hi[t] << 64`. +pub struct Soa<'a> { + pub lo: &'a [u64], + pub hi: &'a [u64], +} + +impl Soa<'_> { + pub fn len(&self) -> usize { + self.lo.len() + } + pub fn tile(&self, t: usize) -> u128 { + self.lo[t] as u128 | (self.hi[t] as u128) << 64 + } +} + +/// Split tiles into owned SoA planes (setup only, not timed). +pub fn split(tiles: &[u128]) -> (Vec, Vec) { + (tiles.iter().map(|&t| t as u64).collect(), tiles.iter().map(|&t| (t >> 64) as u64).collect()) +} + +/// `out = a + b` field-wise over whole slices; the tail uses the u128 form. +pub fn soa_add_slices(a: &Soa, b: &Soa, out_lo: &mut [u64], out_hi: &mut [u64]) { + let n = a.len(); + let body = n - n % 8; + for i in (0..body).step_by(8) { + let [rl, rh, _, _] = soa_add::

( + U64x8::from_slice(&a.lo[i..]), + U64x8::from_slice(&a.hi[i..]), + U64x8::from_slice(&b.lo[i..]), + U64x8::from_slice(&b.hi[i..]), + ); + rl.copy_to_slice(&mut out_lo[i..]); + rh.copy_to_slice(&mut out_hi[i..]); + } + for t in body..n { + let r = add::

(a.tile(t), b.tile(t)); + out_lo[t] = r as u64; + out_hi[t] = (r >> 64) as u64; + } +} + +/// `out = a - b` field-wise over whole slices. +pub fn soa_sub_slices(a: &Soa, b: &Soa, out_lo: &mut [u64], out_hi: &mut [u64]) { + let n = a.len(); + let body = n - n % 8; + for i in (0..body).step_by(8) { + let [rl, rh] = soa_sub::

( + U64x8::from_slice(&a.lo[i..]), + U64x8::from_slice(&a.hi[i..]), + U64x8::from_slice(&b.lo[i..]), + U64x8::from_slice(&b.hi[i..]), + ); + rl.copy_to_slice(&mut out_lo[i..]); + rh.copy_to_slice(&mut out_hi[i..]); + } + for t in body..n { + let r = sub::

(a.tile(t), b.tile(t)); + out_lo[t] = r as u64; + out_hi[t] = (r >> 64) as u64; + } +} + +/// Exact per-field fold over SoA tiles: unsigned sums and signed sums. +/// +/// The signed sum uses `sum(sext f) = sum(f) - 2^W * #(top bit set)`, so both +/// come from one pass: one masked field read and one top-bit read per field. +pub fn soa_fold(a: &Soa) -> (Vec, Vec) { + let (n, nf) = (a.len(), P::N as usize); + let body = n - n % 8; + let vm = U64x8::splat(P::VMASK as u64); + let one = U64x8::splat(1); + let top = U64x8::splat((P::W - 1) as u64); + // One pass over the tiles, one accumulator pair per field. + let mut acc = [U64x8::splat(0); 16]; + let mut cnt = [U64x8::splat(0); 16]; + for i in (0..body).step_by(8) { + let l = U64x8::from_slice(&a.lo[i..]); + let h = U64x8::from_slice(&a.hi[i..]); + for j in 0..nf { + let off = j as u32 * P::W; + let v = if off + P::W <= 64 { + l >> U64x8::splat(off as u64) + } else if off >= 64 { + h >> U64x8::splat((off - 64) as u64) + } else { + (l >> U64x8::splat(off as u64)) | (h << U64x8::splat((64 - off) as u64)) + } & vm; + acc[j] += v; + cnt[j] += (v >> top) & one; + } + } + let mut u: Vec = acc[..nf].iter().map(|x| x.reduce_sum()).collect(); + let mut neg: Vec = cnt[..nf].iter().map(|x| x.reduce_sum()).collect(); + for t in body..n { + for j in 0..P::N { + let f = get::

(a.tile(t), j); + u[j as usize] += f as u64; + neg[j as usize] += (f >> (P::W - 1)) as u64; + } + } + let s = u + .iter() + .zip(&neg) + .map(|(&s, &c)| s as i64 - ((c as i64) << P::W)) + .collect(); + (u, s) +} + +// ── Realization 3: widened u32 lanes over U32x16 ───────────────────────── + +/// Load 16 lanes. `U32x16::from_slice` is not on every realization (NEON +/// and wasm lack it), `from_array` is. +#[inline] +fn load16(s: &[u32]) -> U32x16 { + U32x16::from_array(s[..16].try_into().expect("16 lanes")) +} + +/// Field-major widened layout: `fields[j][t]` = field `j` of tile `t`. +pub fn widen(tiles: &[u128]) -> Vec> { + (0..P::N) + .map(|j| tiles.iter().map(|&t| get::

(t, j)).collect()) + .collect() +} + +/// `out = (a + b) mod 2^W` per lane; lengths must be a multiple of 16. +/// +/// The `mod 2^W` is two shifts, not a mask: on NEON and wasm `U32x16` has +/// only `+`, `^` and the shifts, so `&` would not compile there. +pub fn widened_add(a: &[u32], b: &[u32], out: &mut [u32]) { + let k = U32x16::splat(32 - P::W); + for i in (0..a.len()).step_by(16) { + let s = (load16(&a[i..]) + load16(&b[i..])) << k >> k; + s.to_array() + .iter() + .zip(&mut out[i..i + 16]) + .for_each(|(&v, o)| *o = v); + } +} + +/// Exact unsigned fold of one widened field. A `u32` lane can absorb +/// `floor((2^32 - 1) / (2^W - 1))` values before it overflows, so the lanes are +/// flushed to `u64` after that many steps. +pub fn widened_fold(a: &[u32]) -> u64 { + let block = (u32::MAX / (P::VMASK as u32)) as usize; + let (mut total, mut acc, mut steps) = (0u64, U32x16::splat(0), 0usize); + for i in (0..a.len()).step_by(16) { + // Not `+=`: `U32x16` has no `AddAssign` on NEON or wasm. + #[allow(clippy::assign_op_pattern)] + { + acc = acc + load16(&a[i..]); + } + steps += 1; + if steps == block { + total += acc.to_array().iter().map(|&x| x as u64).sum::(); + acc = U32x16::splat(0); + steps = 0; + } + } + total + acc.to_array().iter().map(|&x| x as u64).sum::() +} + +// ── Realization 4: bit planes over 64 tiles ────────────────────────────── + +/// `planes[j][blk * W + k]`: bit `k` of field `j` for the 64 tiles of block +/// `blk`, one bit per tile. A block's `W` planes are contiguous. +pub fn to_planes(tiles: &[u128]) -> Vec> { + let (w, blocks) = (P::W as usize, tiles.len().div_ceil(64)); + let mut p = vec![vec![0u64; blocks * w]; P::N as usize]; + for (t, &tile) in tiles.iter().enumerate() { + for j in 0..P::N { + let f = get::

(tile, j); + for k in 0..w { + p[j as usize][t / 64 * w + k] |= ((f >> k & 1) as u64) << (t % 64); + } + } + } + p +} + +/// Ripple add of two `W`-plane fields over one 64-tile block; returns the +/// carry-out mask. +#[inline] +pub fn planes_add(a: &[u64], b: &[u64], out: &mut [u64]) -> u64 { + let mut c = 0; + for k in 0..a.len() { + let (s, c2) = full_add_mask(a[k], b[k], c); + out[k] = s; + c = c2; + } + c +} + +/// Exact fold of one plane-stored field: `sum = Σ_k 2^k popcount(plane_k)`; +/// the signed sum gives the top plane the weight `-2^(W-1)`. +pub fn planes_fold(planes: &[u64], w: usize) -> (u64, i64) { + let mut pop = [0u64; 32]; + for blk in planes.chunks_exact(w) { + for (k, &x) in blk.iter().enumerate() { + pop[k] += x.count_ones() as u64; + } + } + let mut u = 0u64; + let mut s = 0i64; + for (k, &c) in pop[..w].iter().enumerate() { + u += c << k; + s += if k + 1 == w { + -((c as i64) << k) + } else { + (c as i64) << k + }; + } + (u, s) +} diff --git a/tests/microtile.rs b/tests/microtile.rs new file mode 100644 index 00000000..7d3e1e42 --- /dev/null +++ b/tests/microtile.rs @@ -0,0 +1,220 @@ +#![cfg(feature = "std")] +//! D-SIMD-MICRO-MACRO-1 step 1 exactness (`examples/microtile/tile.rs`): +//! every realization against the unpacked per-field oracle. + +#[path = "../examples/microtile/tile.rs"] +mod tile; + +use ndarray::simd::U64x8; +use tile::*; + +fn xorshift(s: &mut u64) -> u64 { + *s ^= *s << 13; + *s ^= *s >> 7; + *s ^= *s << 17; + *s +} + +fn rand128(s: &mut u64) -> u128 { + (xorshift(s) as u128) << 64 | xorshift(s) as u128 +} + +/// Field values that hit carry and sign boundaries: 0, 1, the sign bit, the +/// sign bit - 1, max - 1, max. +fn edges() -> Vec { + let m = P::VMASK as u32; + let h = 1u32 << (P::W - 1); + vec![0, 1, h - 1, h, m - 1, m] +} + +/// Random tiles plus every tile whose fields all share one edge value, plus +/// tiles that alternate two edge values (neighbour fields differ). +fn corpus(seed: u64, random: usize) -> Vec { + let mut s = seed; + let mut v: Vec = (0..random).map(|_| rand128(&mut s)).collect(); + let e = edges::

(); + for &x in &e { + for &y in &e { + let fields: Vec = (0..P::N).map(|j| if j % 2 == 0 { x } else { y }).collect(); + v.push(pack::

(&fields, rand128(&mut s))); + } + } + v +} + +fn check_scalar(name: &str) { + let c = corpus::

(0x9E37_79B9_7F4A_7C15, 4000); + for (i, &a) in c.iter().enumerate() { + let b = c[(i * 7 + 3) % c.len()]; + let r = add::

(a, b); + let (want, carry, ovf) = oracle_add::

(a, b); + assert_eq!(r, want, "{name} add tile {i}: {a:#x} + {b:#x}"); + assert_eq!(flags::

(carry_out::

(a, b, r)), carry, "{name} carry tile {i}"); + assert_eq!(flags::

(signed_overflow::

(a, b, r)), ovf, "{name} signed overflow tile {i}"); + let d = sub::

(a, b); + let (want, borrow) = oracle_sub::

(a, b); + assert_eq!(d, want, "{name} sub tile {i}: {a:#x} - {b:#x}"); + assert_eq!(flags::

(borrow_out::

(a, b, d)), borrow, "{name} borrow tile {i}"); + assert_eq!(stencil::

(a), oracle_stencil::

(a), "{name} stencil tile {i}: {a:#x}"); + } +} + +#[test] +fn u24x5_u128_matches_the_oracle() { + check_scalar::("U24x5"); +} + +#[test] +fn u9x14_u128_matches_the_oracle() { + check_scalar::("U9x14"); +} + +/// The fields that straddle bit 64 must take a carry and a borrow across the +/// two words: low part all ones plus one, and zero minus one. +fn check_straddle(name: &str, j: u32) { + assert!(j * P::W < 64 && j * P::W + P::W > 64, "{name}: field {j} does not straddle bit 64"); + let low_bits = 64 - j * P::W; + let mut f = vec![0u32; P::N as usize]; + f[j as usize] = (1 << low_bits) - 1; + let a = pack::

(&f, 0); + let mut g = vec![0u32; P::N as usize]; + g[j as usize] = 1; + let b = pack::

(&g, 0); + let (alo, ahi) = (U64x8::splat(a as u64), U64x8::splat((a >> 64) as u64)); + let (blo, bhi) = (U64x8::splat(b as u64), U64x8::splat((b >> 64) as u64)); + let [rl, rh, _, _] = soa_add::

(alo, ahi, blo, bhi); + let r = rl.to_array()[0] as u128 | (rh.to_array()[0] as u128) << 64; + assert_eq!(get::

(r, j), 1 << low_bits, "{name}: carry did not cross bit 64"); + let [dl, dh] = soa_sub::

(U64x8::splat(0), U64x8::splat(0), blo, bhi); + let d = dl.to_array()[0] as u128 | (dh.to_array()[0] as u128) << 64; + assert_eq!(get::

(d, j), P::VMASK as u32, "{name}: borrow did not cross bit 64"); +} + +#[test] +fn straddling_fields_carry_and_borrow_across_the_word_boundary() { + check_straddle::("U24x5", 2); + check_straddle::("U9x14", 7); +} + +fn check_soa(name: &str) { + // 1001 tiles: not a multiple of 8, so the u128 tail runs too. + let a = corpus::

(11, 1001 - 36); + let b = corpus::

(12, 1001 - 36); + assert_eq!(a.len(), 1001); + let ((alo, ahi), (blo, bhi)) = (split(&a), split(&b)); + let (sa, sb) = (Soa { lo: &alo, hi: &ahi }, Soa { lo: &blo, hi: &bhi }); + let (mut ol, mut oh) = (vec![0u64; a.len()], vec![0u64; a.len()]); + soa_add_slices::

(&sa, &sb, &mut ol, &mut oh); + for t in 0..a.len() { + let r = ol[t] as u128 | (oh[t] as u128) << 64; + assert_eq!(r, oracle_add::

(a[t], b[t]).0, "{name} soa add tile {t}"); + } + soa_sub_slices::

(&sa, &sb, &mut ol, &mut oh); + for t in 0..a.len() { + let r = ol[t] as u128 | (oh[t] as u128) << 64; + assert_eq!(r, oracle_sub::

(a[t], b[t]).0, "{name} soa sub tile {t}"); + } + // Carry-out planes from the vector form against the oracle. + for i in (0..a.len() - a.len() % 8).step_by(8) { + let [rl, rh, cl, ch] = soa_add::

( + U64x8::from_slice(&alo[i..]), + U64x8::from_slice(&ahi[i..]), + U64x8::from_slice(&blo[i..]), + U64x8::from_slice(&bhi[i..]), + ); + let (rl, rh, cl, ch) = (rl.to_array(), rh.to_array(), cl.to_array(), ch.to_array()); + for k in 0..8 { + let r = rl[k] as u128 | (rh[k] as u128) << 64; + let c = cl[k] as u128 | (ch[k] as u128) << 64; + assert_eq!( + flags::

(c), + flags::

(carry_out::

(a[i + k], b[i + k], r)), + "{name} soa carry tile {}", + i + k + ); + assert_eq!(flags::

(c), oracle_add::

(a[i + k], b[i + k]).1, "{name} soa carry tile {}", i + k); + } + } + let (u, s) = soa_fold::

(&sa); + let (wu, ws) = oracle_fold::

(&a); + assert_eq!(u, wu, "{name} soa unsigned fold"); + assert_eq!(s, ws, "{name} soa signed fold"); +} + +#[test] +fn u24x5_soa_matches_the_oracle() { + check_soa::("U24x5"); +} + +#[test] +fn u9x14_soa_matches_the_oracle() { + check_soa::("U9x14"); +} + +fn check_widened_and_planes(name: &str, n: usize) { + let a = corpus::

(21, n - 36); + let b = corpus::

(22, n - 36); + let (wa, wb) = (widen::

(&a), widen::

(&b)); + let (wu, ws) = oracle_fold::

(&a); + let (pa, pb) = (to_planes::

(&a), to_planes::

(&b)); + for j in 0..P::N as usize { + let mut out = vec![0u32; n]; + widened_add::

(&wa[j], &wb[j], &mut out); + for t in 0..n { + assert_eq!( + out[t], + get::

(oracle_add::

(a[t], b[t]).0, j as u32), + "{name} widened add field {j} tile {t}" + ); + } + assert_eq!(widened_fold::

(&wa[j]), wu[j], "{name} widened fold field {j}"); + + let w = P::W as usize; + for blk in 0..n / 64 { + let ak = &pa[j][blk * w..(blk + 1) * w]; + let bk = &pb[j][blk * w..(blk + 1) * w]; + let mut out = vec![0u64; w]; + let c = planes_add(ak, bk, &mut out); + for lane in 0..64 { + let t = blk * 64 + lane; + let v = (0..w).fold(0u32, |acc, k| acc | ((out[k] >> lane & 1) as u32) << k); + let (r, carry, _) = oracle_add::

(a[t], b[t]); + assert_eq!(v, get::

(r, j as u32), "{name} planes add field {j} tile {t}"); + assert_eq!(c >> lane & 1, (carry >> j & 1) as u64, "{name} planes carry field {j} tile {t}"); + } + } + assert_eq!(planes_fold(&pa[j], w), (wu[j], ws[j]), "{name} planes fold field {j}"); + } +} + +#[test] +fn u24x5_widened_and_planes_match_the_oracle() { + check_widened_and_planes::("U24x5", 1024); +} + +#[test] +fn u9x14_widened_and_planes_match_the_oracle() { + check_widened_and_planes::("U9x14", 1024); +} + +/// A u32 lane holds only 256 maximal U24 values; the widened fold must flush +/// before then. 8192 maximal tiles put 512 values in every lane. +#[test] +fn widened_fold_does_not_wrap_its_u32_lanes() { + let m = U24x5::VMASK as u32; + let field = vec![m; 8192]; + assert_eq!(widened_fold::(&field), m as u64 * 8192); +} + +/// Metadata of the left operand survives every operation, including when the +/// right operand's metadata is all ones and every field overflows. +#[test] +fn metadata_is_the_left_operands() { + let meta_a = 0xA5u128 << 120; + let a = pack::(&[U24x5::VMASK as u32; 5], meta_a); + let b = pack::(&[1; 5], !0); + assert_eq!(add::(a, b) & U24x5::META, meta_a); + assert_eq!(sub::(a, b) & U24x5::META, meta_a); + assert_eq!(stencil::(a) & U24x5::META, meta_a); + assert_eq!(add::(a, b) & U24x5::FIELDS, 0, "every field wraps to zero"); +}