diff --git a/.claude/board/STATUS_BOARD.md b/.claude/board/STATUS_BOARD.md index e2e7bad14..afea71083 100644 --- a/.claude/board/STATUS_BOARD.md +++ b/.claude/board/STATUS_BOARD.md @@ -1,3 +1,11 @@ +## D-ART — Algebraic recipe table over concepts and terminals (2026-10-08) + +Board entry: `entries/2026-10-08-algebraic-recipe-probe.md`. A linear canonicaliser and a static guarded recipe table reduce (concept, transformation chain, terminal) to the minimal work; affine transforms of `CrossPowerSums`; projected order statistics of `m`-fold families. + +| D-id | scope | status | gate / falsifier | +|---|---|---|---| +| **D-ART-1** | Recipe table: constants, centre phasor, isotropic invariant, origin norm, sector quantile, vertex phasor, summary affine, materialise; PowerSums affine chains; order statistics by generate/analytic/profile LUT/sector permutation; lookup cost | Shipped (probe; `CrossPowerSums` affine transform proposed for ndarray) | 36,000 cases equal the materialising oracle; summary path bitwise equal in `i128`; dropping the isotropy guard, the half-turn fold, the rank reversal, the half-sector split, the population refusal or the general `‖L‖_F²` each fails | + ## D-PHT — Phasors without transcendental calls (2026-10-08) Board entry: `entries/2026-10-08-phasor-trig-probe.md`. Native `sin_cos`, a `u32`-turns phase LUT, CORDIC and complex recurrence on Wankel apex coordinates and coherent interference, plus an amplitude-bounded exact early exit. @@ -13,6 +21,7 @@ Board entry: `entries/2026-10-08-mexhat-bucket-cascade-probe.md`. Popcount stack | D-id | scope | status | gate / falsifier | |---|---|---|---| | **D-MHB-1** | Mexican-hat response: direct geometry (D), ring buckets (E), bit-sliced weight planes (F), per-row choice (G), mask-risc lane (B), exact-bound early exit (H) | Shipped (probe; R-MHB-1 proposed) | D = F = G = B = full scan asserted on every centre; H equals the full decision on every query; one zero crossing, one minimum, symmetry; dropping the unvisited-row bound decides 52 / 512 wrongly | +| **D-MHB-2** | Per-query D/F choice from the exact window count (Q); early exit over D rows (HD) and after a per-query choice (HQ) | Shipped (probe; Q proposed with R-MHB-1) | Q = D for cuts 0, 100, ∞ on every centre; HD/HQ decide as H after the same rows; a D row dropping one cell fails Q. Entry: `entries/2026-10-08-mexhat-df-choice.md` | ## D-RPF — Mask × fold over projections of resident bytes (2026-10-08) diff --git a/.claude/board/entries/2026-10-08-algebraic-recipe-probe.md b/.claude/board/entries/2026-10-08-algebraic-recipe-probe.md new file mode 100644 index 000000000..9eb68d45b --- /dev/null +++ b/.claude/board/entries/2026-10-08-algebraic-recipe-probe.md @@ -0,0 +1,183 @@ +# 2026-10-08 — Algebraic recipe table × percentile profiles × Wankel phasors (D-ART-1) + +**Status:** MEASURED. Probe only; no production code changes. +- **Probe:** `crates/lance-graph-mask-risc/examples/algebraic_recipe_probe.rs`. +- **Revision:** lance-graph `origin/main` 7f63ae9d. +- **Host:** Intel Xeon @ 2.8 GHz, 4 vCPU, `avx2` build, release, `debug = 0`, rustc 1.98.1. Medians of 5–31 runs. + +```text +CARGO_PROFILE_RELEASE_DEBUG=0 cargo run --release -p lance-graph-mask-risc --example algebraic_recipe_probe +``` + +## Question + +Given a concept (a shape family and its parameters), a chain of transformations and a requested terminal, how much geometric and trigonometric work can a small, guarded recipe table remove before anything is evaluated? And does it know when it cannot? + +## Inventory (read, not assumed) + +| item | state | +|---|---| +| `PowerSums`, `CrossPowerSums` (ndarray `simd_masking_ops.rs`, re-exported in `simd`) | **SHIPPED** as integer folds (`n`, `Σ`, `Σ²`, `Σxy`) with `checked_merge`. **MISSING:** any shift, scale or affine method; any mean/variance terminal on these types. | +| `MomentsU32` (ndarray `hpc/statistics.rs`) | **SHIPPED**, unsigned, with `mean`/`variance`. **SEMANTICALLY DIFFERENT** from `PowerSums`; not interchangeable. | +| Float terminals over power sums | **SHIPPED** consumer-side in lance-graph `crates/jc/src/stats.rs` (ANOVA, Pearson, covariance, regression). | +| `GroupPowerSumsI32`, `GroupCrossPowerSumsI32` | **SHIPPED** as mask-risc terminals (tiled path only). **PARTIALLY AVAILABLE** end to end: no Quack `Agg`/`GroupAgg` variant reaches them. | +| `RollingFloor`, `EmpiricalShape` (ndarray `hpc/rolling_floor.rs`) | **SHIPPED**, integer. The location-scale quantile transform `Q = μ + (σ/σ₀)(Q₀ − μ₀)`, Euclidean floor division, clamped to `u32`. Nearest-rank floor `rank_per_10000`; this probe uses that convention. | +| Arithmetic expression IR / algebraic simplifier | **MISSING.** Quack's typed IR is Boolean predicates only. The rewrites that exist (SQL three-valued logic, survivor-gate hoisting, ternlog fusion) are all Boolean. The planner's ten `OptimizerRule`s all return `Unchanged`. DataFusion is in its grace period. | +| Rotation / complex / phasor type | **MISSING** in ndarray `hpc` (only `GivensRotation`, no apply). D-PHT-1's LUT is a probe; ndarray PR #343 (`hpc::phase`) is open. | +| NARS inference recipes | **SEMANTICALLY INCOMPATIBLE.** Not touched; mathematical identities are not inference rules. | + +## What the probe holds + +- **Concepts:** a regular `m`-fold family `c + r·e^{i(θ + 2πk/m)}`, with either a Wankel centre `e·e^{i3θ}` or a fixed centre. No points are stored. +- **Canonicaliser:** one left fold of the chain. + - Rotations stay exact `u32` turns; `R(a)R(b) = R(a + b)` is a wrapping add. + - A negative similarity becomes a half turn plus a positive scale. + - Translations stay lazy, each tagged with the rotation and scale in force; trigonometry is spent only if a terminal reads the translation. + - The class drops back down `Affine ⊃ Similarity ⊃ Rotation ⊃ Identity` when scales multiply to 1 and rotations cancel. + - It is linear in the chain length, so it terminates; there is no search and no cycle. +- **Recipe table:** `static [[[Rid; 9]; 4]; 2]` over (shape, linear class, terminal). The key only proposes. Guards (`m ≥ 3` for isotropy, class, translation present) run at evaluation, and a failed guard falls through to `Materialize`. +- **Oracle:** materialises the `m` vertices and applies the raw chain to each point, independently of the canonical form. + +## Results + +### 1. The Wankel concept under `R(+1/8) · T(5, −3) · R(−1/8) · S(2) · S(0.5)` + +Canonical form: the identity plus one lazy translation (3 rewrites). The oracle spends 10 transcendental calls and materialises 3 vertices for every query. + +| query | recipe | trig | vertices | remaining work | +|---|---|---|---|---| +| A count | Const | 0 | 0 | `3` | +| B Σ positions, centroid | CenterPhasor | 2 | 0 | the centre phasor `3θ`, plus one rotation of the lazy translation | +| C centred Σ\|z − c\|², centred Σ(x−x̄)(y−ȳ) | IsoInvariant | 0 | 0 | `3R²`, `0` | +| D Σ\|z\|², translated | OriginNorm | 2 | 0 | centred invariant + `3\|c'\|²` | +| D′ Σ\|z\|² under `R · S(−1.5)`, untranslated | OriginNorm | 0 | 0 | `3s²(R² + e²)` | +| E apex 1 | VertexPhasor | 3 | 0 | one vertex phasor | +| E′ x-quantile, p = 0.5 | SectorQuantile | 3 | 0 | centre, translation, one cosine | +| F collision with (110, 0), d = 5 | Materialize | 10 | 3 | refused: no summary decides it | + +All equal the oracle. + +### 2. Guards, 36,000 random (concept, chain, terminal) cases + +- m ∈ {2, 3, 4, 5, 8, 16, 256}; chains of 0–4 rotations, signed scales, exact zero scales, translations and general linear maps. +- All equal the oracle within 1e-9 of the compared values plus an absolute floor of 1e-13 of the case's Σ|p|² + m (for results that cancel to zero). +- A translation applied after a zero scale survives it: `S(0) · T(5, −3)` has centroid (5, −3). + +| recipe | cases | trig/case | oracle trig/case | exactness | +|---|---|---|---|---| +| Const | 4,000 | 0.17 | 51.4 | real identity | +| CenterPhasor | 8,000 | 0.99 | 51.4 | real identity | +| IsoInvariant | 7,000 | 0.16 | 58.2 | real identity | +| OriginNorm | 3,500 | 0.64 | 58.2 | real identity | +| SectorQuantile | 4,000 | 1.99 | 51.4 | real identity | +| VertexPhasor | 4,000 | 2.19 | 51.4 | real identity | +| Materialize | 5,500 | 38.5 | 38.4 | reference | + +The nonzero trig on Const and IsoInvariant is the canonicaliser leaving the similarity form (a general linear map materialises `L` once). + +**Refusals, measured:** +- `m = 2` under `diag(2, 1)`: the isotropic formula gives 500.00, the truth is 777.16. The `m ≥ 3` guard refuses, and the recipe falls back to `Materialize`. +- A shear `[1 1.5; 0 1]` gives Σ|z − c|² = 637.50, not the rotation invariant 300.00. The class is `Affine`, and `‖L‖_F²` carries it. + +### 3. Populations: affine chains over `CrossPowerSums` + +T3∘T2∘T1 integer affine maps (composed `A = [0 −2; −1 −2]`, `t = (20, 4)`). Coordinates are in [−100, 100]. + +| rows | A: materialise each step + ndarray fold | A2: composed map per row + scalar fold | B: ndarray fold + 3 summary steps | C: ndarray fold + 1 composed step | +|---|---|---|---|---| +| 1,000 | 12.0 | 3.6 | 4.8 | 4.8 | +| 65,536 | 18.8 | 3.7 | 4.9 | 4.8 | +| 1,000,000 | 20.5 | 3.9 | 5.4 | 5.7 | + +ns per row. All four are **bitwise equal** (`i128`). A warm composed summary step costs **13.7 ns**, independent of the row count. + +- **f64:** the same identity is not bitwise. Max relative difference 1.5e-14; 1 of 6 fields bitwise equal. +- **Moments are not sufficient:** {0, 3, 3} and {1, 1, 4} on `y = 0` share `n, Σx, Σy, Σx², Σy², Σxy` = (3, 6, 0, 18, 0, 0). Their maxima are 3 vs 4, and a collision with (4, 0) is false vs true. The table maps Population × {quantile, vertex, collide} to `Materialize` for every class. +- **Group keys:** with the key `x ≥ 0` and the map `x + 50`, pushing the fold below the transform is wrong (group sizes 2,052 vs 3,059). With a key on an untouched column (id parity), it is exact. + +### 4. One projected order statistic per instance + +**Symmetry, measured.** +- `Q(θ + 1/m) = Q(θ) = Q(−θ)` to 2.6e-15 on 20,000 phases. +- Ties between `cos(θ + 2πk/m)` occur only at `θ ≡ 0 mod 1/(2m)` turn. So inside each half-sector the rank → vertex map is **fixed**, and a phase-conditioned percentile profile collapses to two permutations of `m` entries plus one cosine (arm E). +- 1/3 turn is not representable in `u32` turns (2^32 mod 3 = 1). Threefold symmetry holds in reals, and to one ulp-turn in `u32`. + +**Arms.** Error is the max |error| on R = 100 over 200,003 phases × every one of the `m` ranks (the oracle sorts all `m` projections once per phase). Timings are ns per instance at 1M instances; 65,536 instances agree within 10 %. + +| m | A generate + select | B analytic (m = 3) | C profile LUT, 64 buckets | E sector permutation + 1 `cos` | C table | E table | C error, 64 buckets | C error, 63 buckets | E error | +|---|---|---|---|---|---|---|---|---|---| +| 3 (Wankel) | 108 | 39 | 34 | 50 | 780 B | 24 B | 1.3e-2 | **1.4** | 2.0e-13 | +| 16 | 383 | — | 10.6 | 42 | 4.2 KB | 128 B | 4.7e-4 | 3.1e-1 | 2.1e-13 | +| 256 | 5,037 | — | 10.9 | 42 | 66.6 KB | 2 KB | 7.7e-6 | 1.9e-2 | 2.0e-13 | + +- B is one `sin_cos`, ±120° by constants, `cos 3θ = 4c³ − 3c`, and a three-element sorting network. +- E's time is two `f64` cosines: the value and the Wankel or fixed centre. +- A query that needs no profile (Σ|z − c|²) costs 1.2–3 ns. + +### 5. Lookup and rewrite cost + +- **Recipe lookup over 4,096 random keys:** static array 1.2 ns, `match` 2.0 ns, `HashMap` (SipHash) 27 ns. +- A `HashMap` whose hasher sends every key to one bucket still returns the right recipe (106 ns): structural key equality, not the hash, decides. +- **8-step chain:** canonicalising costs 49–72 ns; canonicalise + lookup + evaluate Σ|z − c|² costs 62–98 ns, against 301–496 ns for the oracle (3 vertices × 8 raw steps). Ranges are two runs. + +### 6. 1-D affine quantile + +`Q(aX + b)` is `b + a·X[rank]` for `a ≥ 0` and `b + a·X[n − 1 − rank]` for `a < 0`. It equals the sorted oracle on 2,000 random cases. + +## Disable runs + +Each disable ran after the commit, turned red, and was restored afterwards. + +| disable | effect | +|---|---| +| isotropy guard `m ≥ 3` → `true` | case 0, `m = 2` under a linear map: recipe ≠ oracle | +| negative scale not folded into a half turn | case 6, `m = 16`, Σ positions sign-flipped | +| 1-D quantile rank reversal for `a < 0` dropped | `a −2, b 3, p 2270` ≠ sorted oracle | +| sector permutation without the half-sector split | case 3, `m = 4` quantile ≠ oracle | +| Population × quantile mapped to a summary | the table assertion fires | +| `‖L‖_F²` of a general map treated as a similarity's | case 2, `m = 3` Σ\|z − c\|² 43,271 vs 136,940 | +| translation multiplier not updated by later scales | case 11, `T · S(0) · R`: Σ positions ≠ oracle | +| the first version's `scale / scale_at` translation factor | case 14, `S(0) · T`: Σ positions (0, 0) vs the oracle's translated sum | + +## Findings + +1. **The terminal decides the representation, and most terminals need no vertices.** + - Count, centred second moments and an untranslated Σ|z|² need no phase at all. + - Sum and centroid need one centre phasor; a vertex needs one phasor; an order statistic needs one cosine. + - Only collision falls back to materialisation, and the table refuses it rather than approximating. +2. **The centred second moment of any regular family with `m ≥ 3` is `(m r²/2) · ‖L‖_F²` for every linear map `L`, not only rotations.** + - It is phase-free even under shear. + - `m = 2` is a rank-one dyad, and the guard is load-bearing. +3. **The percentile profile is two permutations, not a table.** + - Ranks only change at `1/(2m)` turn, so the phase-conditioned profile reduces to a fixed vertex per (half-sector, rank). + - A bucketed LUT is fast (10 ns) but must align its bucket edges with the tie points: 63 buckets put a kink inside a bucket and the error jumps from 1e-2 to 1.4 (m = 3). + - For m = 3 no table pays: the analytic arm B is within 1.15× of the LUT, and exact. +4. **Summary transforms are exact and O(1), but the fold of the base population is the floor.** + - Folding once and transforming the summary removes every materialisation (2.5–4× against A). + - A fused per-row composed transform with a scalar fold (A2, 3.7 ns/row) is faster than the shipped grouped `masked_group_cross_power_sums_i32` without a transform (4.8 ns/row). + - The recipe's real value is the warm path: one fold, any number of affine queries at 14 ns each. +5. **Hashing costs more than many of the recipes it selects.** Σ|z − c|² is one multiply-add, and a SipHash lookup is 27 ns. A static array indexed by enum discriminants is 1.2 ns and cannot collide. + +## Recommendations + +| item | verdict | correctness contract | measurement | +|---|---|---|---| +| Affine transform of `CrossPowerSums` (`S' = AS + nt`, `M' = AMAᵀ + AStᵀ + tSᵀAᵀ + nttᵀ`), checked `i128` | **ADOPT** (a small ndarray PR, as a method on `CrossPowerSums`) | bitwise equal to materialise-and-fold within the documented widths; refuses when the group key reads a transformed column | 14 ns per summary step against 5–20 ns per row | +| Isotropic invariant `(m r²/2)·‖L‖_F²` for regular families | **ADOPT as a recipe** where such a concept exists | real identity; `m ≥ 3` guard | 0 trig against 56–64 for the oracle | +| Sector-permutation order statistic | **ADOPT as the method** for projected order statistics of `m`-fold families | exact up to `f64` rounding (2e-13) | 42 ns, flat in `m`; A is 5 µs at m = 256 | +| Sector permutation with the phase LUT from ndarray #343 instead of `f64` `cos` | **PROBE** | the LUT's stated bound times `r·ρ` | expected ≈ C's 10 ns; not measured | +| Phase-conditioned profile LUT | **PROBE** for m ≥ 16 only, bucket count a multiple of 2 per sector | bounded; edges must sit on the tie points | 10 ns, 4–67 KB | +| Profile LUT for the Wankel triangle (m = 3) | **REJECT** | — | the analytic arm is within 1.15× and exact | +| Static recipe table + linear canonicaliser as a production layer | **PARK** | — | no arithmetic expression IR exists to host it; building one is not justified until a consumer asks concept + terminal queries | +| `HashMap`-keyed recipe lookup | **REJECT** | — | 27 ns against 1.2 ns | +| e-graph or search-based rewriting | **REJECT** for this rule set | — | a single linear fold terminates and is stable | +| `RollingFloor` / `EmpiricalShape` for deterministic geometry | **REJECT the role** | — | static geometry needs no statistical frame | +| `RollingFloor` / `EmpiricalShape` for learned, drifting profiles | **PARK** | — | the location-scale transform matches, but no learned-shape workload exists here | +| NARS inference recipes | untouched | — | — | + +## OPEN + +- One host class only. +- The fused arm A2 beating the shipped grouped fold is one measurement; whether the group addressing or the mask walk costs the difference is not isolated. +- Multiple quantiles per instance (all `m` order statistics) were not timed. E then needs `m` cosines and the permutations, but no sort. +- Reflection symmetry of a projection along a general direction (the affine `β`) costs one `atan2`. It is counted, not removed. diff --git a/.claude/board/entries/2026-10-08-mexhat-df-choice.md b/.claude/board/entries/2026-10-08-mexhat-df-choice.md new file mode 100644 index 000000000..00aa7ce2f --- /dev/null +++ b/.claude/board/entries/2026-10-08-mexhat-df-choice.md @@ -0,0 +1,67 @@ +# 2026-10-08 — Mexican hat: per-query D/F choice and early exit on the D path (D-MHB-2) + +**Status:** MEASURED. Probe only. Closes the two OPEN items of D-MHB-1. +- **Probe:** `crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs`, new arms Q, HD, HQ and a cut sweep. +- **Revision:** lance-graph `origin/main` 7f63ae9d. +- **Host:** Intel Xeon @ 2.8 GHz, 4 vCPU, `avx2` build, release, `debug = 0`, rustc 1.98.1. Medians of 5 runs; the VM is noisy, about ±15 % between sections. + +## Arms + +- **Q:** fetch the window's 37 words once, count present disk cells exactly (37 popcounts), then run D if the count is ≤ the cut and F otherwise. Cut pinned at 100 cells (D-MHB-1 finding 3). +- **HD / HQ:** arm H's exact suffix-bound early exit with D's per-point row evaluator (HD), or with a per-query D/F choice (HQ). The bounds are the template's, so they hold for any exact row evaluator. + +## Correctness + +- Q equals D on every checked centre for cuts 0, 100 and `u32::MAX`, on both grids and every density. +- HD and HQ return the same decision as H, after the same number of rows, on every query. +- Disable run: a D row that drops one window cell makes `Q (cut 100) differs from D` fire on the first grid. + +## Results + +### Q against min(D, F): 1M grid, 2,048 queries, ns/query + +| ρ | cut 25 | cut 50 | cut 100 | cut 200 | cut 400 | D | F | +|---|---|---|---|---|---|---|---| +| 0.01 | 1.19 | 1.16 | 1.16 | 1.22 | 1.21 | 234 | 593 | +| 0.05 | 1.30 | 1.10 | 0.93 | 0.93 | 0.93 | 530 | 579 | +| 0.1 | 1.18 | 1.18 | 1.10 | 1.04 | 1.03 | 602 | 578 | +| 0.2 | 1.40 | 1.42 | 1.33 | 0.97 | 1.03 | 727 | 850 | +| 0.5 | 1.13 | 1.14 | 1.09 | 1.10 | 1.11 | 1,057 | 614 | +| 0.9 | 1.16 | 1.16 | 1.16 | 1.15 | 1.18 | 1,483 | 584 | + +Ratios are Q / min(D, F). Below 1 is noise or the shared word fetch. + +### Early exit at the median threshold, 1M grid, ns/query + +| fixture | H (F rows) | HD (D rows) | HQ | full F | full D | full Q | +|---|---|---|---|---|---|---| +| ρ 0.01 | 515 | 267 | 364 | 596 | 242 | 248 | +| ρ 0.1 | 520 | 559 | 606 | 602 | 642 | 618 | +| ρ 0.5 | 548 | 977 | 1,454 | 601 | 1,773 | 1,497 | +| ρ 0.9 | 481 | 1,258 | 608 | 597 | 1,479 | 687 | +| wavefront + 1 % noise | 487 | 272 | 357 | 589 | 262 | 247 | + +## Findings + +1. **A per-query choice reaches min(D, F) to within 0.93–1.22× at cut 200, and 0.93–1.33× at cut 100.** + - It removes the worst case of either fixed path: D is 2.5× slower than F at ρ 0.9; F is 2.5× slower than D at ρ 0.01. + - Its whole overhead is the exact 37-popcount count, about 40–90 ns. + - On this host cut 200 is slightly better than the pinned 100 (the ρ 0.2 row). The cut stays a host pin, not a constant. +2. **Early exit does not pay on the D path.** + - At sparse density a D row costs as many lookups as it has present cells, which is already few. The early exit's bound bookkeeping costs more than the rows it skips: HD 267 against full D 242. + - HQ pays the full window count up front, which the early exit cannot then save. It is worse than H or HD everywhere. +3. **Early exit stays a modest win on F only:** −10 to −20 % against full F, as in D-MHB-1. + +## Recommendations + +| item | verdict | why | +|---|---|---| +| Per-query D/F choice from the exact window count (Q) | **ADOPT** together with D and F (R-MHB-1's production PR) | within 1.0–1.3× of the better path at every density; removes a 2.5× worst case | +| Early exit over D rows (HD) | **REJECT** | no gain over full D where D is chosen | +| Early exit after a per-query choice (HQ) | **REJECT** | the count it needs cannot be skipped | +| Early exit over F rows (H) | stays **PROBE** | −10 to −20 % | + +## OPEN + +- One host. The cut must be re-measured per host class. +- A density estimate from the grid's global popcount would cost nothing per query, but it fails on local structure (the wavefront fixture). Not measured. diff --git a/.claude/board/entries/README.md b/.claude/board/entries/README.md index 71edbcb8c..5390f9dd0 100644 --- a/.claude/board/entries/README.md +++ b/.claude/board/entries/README.md @@ -25,7 +25,7 @@ index row, (3) no duplicate entry id. Checks 1 and 2 are deliberately opposite directions; the stranding this convention prevents shows up in exactly one of them, never both. -275 entries, 2026-08-06 .. 2026-10-08. +277 entries, 2026-08-06 .. 2026-10-08. | date | entry id | finding | file | |---|---|---|---| @@ -36,12 +36,14 @@ exactly one of them, never both. | 2026-10-08 | `D-PHT-1` | | [2026-10-08-phasor-trig-probe.md](2026-10-08-phasor-trig-probe.md) | | 2026-10-08 | `moore-nars16-isa-visible-representation` | | [2026-10-08-moore-nars16-isa-visible-representation.md](2026-10-08-moore-nars16-isa-visible-representation.md) | | 2026-10-08 | `D-MOORE-NARS-0` | | [2026-10-08-moore-nars-0-recipe-learning-gomoku.md](2026-10-08-moore-nars-0-recipe-learning-gomoku.md) | +| 2026-10-08 | `D-MHB-2` | | [2026-10-08-mexhat-df-choice.md](2026-10-08-mexhat-df-choice.md) | | 2026-10-08 | `D-MHB-1` | | [2026-10-08-mexhat-bucket-cascade-probe.md](2026-10-08-mexhat-bucket-cascade-probe.md) | | 2026-10-08 | `hhtl-nars-moore-value-tenants` | | [2026-10-08-hhtl-nars-moore-value-tenants.md](2026-10-08-hhtl-nars-moore-value-tenants.md) | | 2026-10-08 | `D-RPF-9` | | [2026-10-08-fold-join-deforestation-probe.md](2026-10-08-fold-join-deforestation-probe.md) | | 2026-10-08 | `D-RPF-0` | | [2026-10-08-d-rpf-0-resident-ce64-predicate.md](2026-10-08-d-rpf-0-resident-ce64-predicate.md) | | 2026-10-08 | `coresearch-ce64-moore-masking-wiring` | | [2026-10-08-coresearch-ce64-moore-masking-wiring.md](2026-10-08-coresearch-ce64-moore-masking-wiring.md) | | 2026-10-08 | `ce64-isa-register-contract` | | [2026-10-08-ce64-isa-register-contract.md](2026-10-08-ce64-isa-register-contract.md) | +| 2026-10-08 | `D-ART-1` | | [2026-10-08-algebraic-recipe-probe.md](2026-10-08-algebraic-recipe-probe.md) | | 2026-10-07 | `tinker-janus-fold-harvest` | TinkerPop bulk = K (GroupReduce), ONE_BULK = S, GValue pinning = bundle invalidation; JanusGraph slice = OrderedLaneWitness→Range (30–38×); no new V4 op | [2026-10-07-tinker-janus-fold-harvest.md](2026-10-07-tinker-janus-fold-harvest.md) | | 2026-10-07 | `selector-seq-aware-coverage` | | [2026-10-07-selector-seq-aware-coverage.md](2026-10-07-selector-seq-aware-coverage.md) | | 2026-10-07 | `D-PUZZLE-ATTN-0` | | [2026-10-07-puzzle-attn-0-entropy-attention-crosswords.md](2026-10-07-puzzle-attn-0-entropy-attention-crosswords.md) | diff --git a/crates/lance-graph-mask-risc/examples/algebraic_recipe_probe.rs b/crates/lance-graph-mask-risc/examples/algebraic_recipe_probe.rs new file mode 100644 index 000000000..c46a655bb --- /dev/null +++ b/crates/lance-graph-mask-risc/examples/algebraic_recipe_probe.rs @@ -0,0 +1,1621 @@ +//! D-ART-1 — algebraic recipe table × percentile profiles × Wankel phasors. +//! +//! Question: given a concept (a shape family with its parameters), a chain of +//! transformations and a requested terminal, how much of the geometric and +//! trigonometric work can a small, guarded recipe table eliminate before +//! anything is evaluated — and does it know when it cannot? +//! +//! Nothing here is a production primitive. The probe holds: +//! +//! - a canonicaliser that folds a transformation chain into one symbolic +//! affine (rotations stay exact `u32` turns; translations stay lazy); +//! - a static recipe table `[shape][linear class][terminal] → recipe`, with +//! every guard checked at evaluation, never trusted from the key; +//! - an independent oracle that materialises the vertices and applies the raw +//! chain point by point; +//! - affine transformation of `CrossPowerSums` summaries (exact `i128`); +//! - four implementations of a projected order statistic of an `m`-fold +//! symmetric point family, the Wankel apex triangle being `m = 3`. +//! +//! Run with: +//! +//! ```text +//! CARGO_PROFILE_RELEASE_DEBUG=0 cargo run --release -p lance-graph-mask-risc --example algebraic_recipe_probe +//! ``` + +use ndarray::hpc::rolling_floor::rank_per_10000; +use ndarray::simd::{masked_group_cross_power_sums_i32, CrossPowerSums}; +use std::collections::HashMap; +use std::f64::consts::TAU; +use std::hash::{BuildHasher, Hasher}; +use std::sync::atomic::{AtomicU64, Ordering::Relaxed}; +use std::time::Instant; + +const TURN: f64 = 4_294_967_296.0; + +// ───────────────────────── counters and rng ───────────────────────── + +/// Transcendental calls made by the planner and the oracle (not by the +/// benchmark arms, which use uncounted calls so the counter costs nothing). +static TRIG: AtomicU64 = AtomicU64::new(0); + +fn sc(rad: f64) -> (f64, f64) { + TRIG.fetch_add(1, Relaxed); + rad.sin_cos() +} + +fn trig() -> u64 { + TRIG.load(Relaxed) +} + +fn rad(p: u32) -> f64 { + p as f64 / TURN * TAU +} + +struct Rng(u64); +impl Rng { + fn next(&mut self) -> u64 { + self.0 ^= self.0 << 13; + self.0 ^= self.0 >> 7; + self.0 ^= self.0 << 17; + self.0 + } + fn f(&mut self) -> f64 { + (self.next() >> 11) as f64 / (1u64 << 53) as f64 + } + fn range(&mut self, lo: f64, hi: f64) -> f64 { + lo + (hi - lo) * self.f() + } + fn below(&mut self, n: u64) -> u64 { + self.next() % n + } +} + +fn median_ns(reps: usize, per: usize, mut f: impl FnMut() -> f64) -> f64 { + let mut v: Vec = (0..reps) + .map(|_| { + let t = Instant::now(); + std::hint::black_box(f()); + t.elapsed().as_nanos() as f64 / per as f64 + }) + .collect(); + v.sort_by(f64::total_cmp); + v[reps / 2] +} + +// ───────────────────────── concepts ───────────────────────── + +/// Where an `m`-fold family is centred. +#[derive(Clone, Copy, Debug)] +enum Center { + /// The Wankel apex triangle: centre `e · e^{i3θ}`, coupled to the phase. + Wankel { e: f64 }, + /// A fixed centre, independent of the phase. + Fixed(f64, f64), +} + +/// `m` points `c + r · e^{i(θ + 2πk/m)}`, `k = 0..m`. Holds no points. +#[derive(Clone, Copy, Debug)] +struct Regular { + m: u32, + r: f64, + phase: u32, + center: Center, +} + +/// Transformations, applied left to right. +#[derive(Clone, Copy, Debug)] +enum Xf { + Rotate(u32), + Scale(f64), + Translate(f64, f64), + /// A general linear map `[a, b, c, d]`: `(x, y) ↦ (ax + by, cx + dy)`. + Linear([f64; 4]), +} + +#[derive(Clone, Copy, Debug)] +enum Term { + Count, + Sum, + Centroid, + /// `Σ |p − centroid|²`. + SumSqCentered, + /// `Σ |p|²`. + SumSqOrigin, + /// `Σ (x − x̄)(y − ȳ)`. + CrossCentered, + /// Order statistic of the x projections, `rank_per_10000` convention. + QuantileX(u32), + Vertex(u32), + /// Does any point lie within `d` of `(px, py)`? + Collide(f64, f64, f64), +} + +#[derive(Clone, Copy, Debug, PartialEq)] +enum Answer { + Num(f64), + Pt(f64, f64), + Bool(bool), +} + +// ───────────────────────── recipe table ───────────────────────── + +#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)] +enum ShapeK { + Regular = 0, + Population = 1, +} +#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug, PartialOrd, Ord)] +enum XfK { + Identity = 0, + Rotation = 1, + Similarity = 2, + Affine = 3, +} +#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)] +enum TermK { + Count = 0, + Sum, + Centroid, + SumSqCentered, + SumSqOrigin, + CrossCentered, + Quantile, + Vertex, + Collide, +} +const NT: usize = 9; + +/// What remains after the rewrite. +#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)] +enum Rid { + /// A constant of the concept. + Const, + /// One centre phasor, no vertex. + CenterPhasor, + /// An invariant of the isotropic centred second moment: no phase. + IsoInvariant, + /// Centred invariant plus `m|c'|²`: phase-free unless translated. + OriginNorm, + /// One cosine through the symmetry-reduced sector permutation. + SectorQuantile, + /// One vertex phasor. + VertexPhasor, + /// The summary transforms; the population is never re-read. + SummaryAffine, + /// No sufficient summary: materialise and evaluate. + Materialize, +} + +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +enum Exactness { + /// An identity of real arithmetic; `f64` rounding order changes. + RealIdentity, + /// Bitwise exact in integers within the stated widths. + IntExact, + /// The oracle itself. + Reference, +} + +fn exactness(r: Rid) -> Exactness { + match r { + Rid::SummaryAffine => Exactness::IntExact, + Rid::Materialize => Exactness::Reference, + _ => Exactness::RealIdentity, + } +} + +fn term_k(t: &Term) -> TermK { + match t { + Term::Count => TermK::Count, + Term::Sum => TermK::Sum, + Term::Centroid => TermK::Centroid, + Term::SumSqCentered => TermK::SumSqCentered, + Term::SumSqOrigin => TermK::SumSqOrigin, + Term::CrossCentered => TermK::CrossCentered, + Term::QuantileX(_) => TermK::Quantile, + Term::Vertex(_) => TermK::Vertex, + Term::Collide(..) => TermK::Collide, + } +} + +/// The table for a regular family does not depend on the linear class: the +/// class changes the arithmetic inside a recipe, never whether it applies. +/// Guards that do depend on it run at evaluation. +const REGULAR_ROW: [Rid; NT] = [ + Rid::Const, + Rid::CenterPhasor, + Rid::CenterPhasor, + Rid::IsoInvariant, + Rid::OriginNorm, + Rid::IsoInvariant, + Rid::SectorQuantile, + Rid::VertexPhasor, + Rid::Materialize, +]; +/// A population carries only `n, Σx, Σy, Σx², Σy², Σxy`: order statistics, +/// vertices and collisions are not functions of those. +const POPULATION_ROW: [Rid; NT] = [ + Rid::SummaryAffine, + Rid::SummaryAffine, + Rid::SummaryAffine, + Rid::SummaryAffine, + Rid::SummaryAffine, + Rid::SummaryAffine, + Rid::Materialize, + Rid::Materialize, + Rid::Materialize, +]; +static TABLE: [[[Rid; NT]; 4]; 2] = [[REGULAR_ROW; 4], [POPULATION_ROW; 4]]; + +fn lookup_array(s: ShapeK, x: XfK, t: TermK) -> Rid { + TABLE[s as usize][x as usize][t as usize] +} + +fn lookup_match(s: ShapeK, _x: XfK, t: TermK) -> Rid { + match (s, t) { + (ShapeK::Regular, TermK::Count) => Rid::Const, + (ShapeK::Regular, TermK::Sum | TermK::Centroid) => Rid::CenterPhasor, + (ShapeK::Regular, TermK::SumSqCentered | TermK::CrossCentered) => Rid::IsoInvariant, + (ShapeK::Regular, TermK::SumSqOrigin) => Rid::OriginNorm, + (ShapeK::Regular, TermK::Quantile) => Rid::SectorQuantile, + (ShapeK::Regular, TermK::Vertex) => Rid::VertexPhasor, + (ShapeK::Population, TermK::Quantile | TermK::Vertex | TermK::Collide) => Rid::Materialize, + (ShapeK::Population, _) => Rid::SummaryAffine, + (ShapeK::Regular, TermK::Collide) => Rid::Materialize, + } +} + +const ALL_T: [TermK; NT] = [ + TermK::Count, + TermK::Sum, + TermK::Centroid, + TermK::SumSqCentered, + TermK::SumSqOrigin, + TermK::CrossCentered, + TermK::Quantile, + TermK::Vertex, + TermK::Collide, +]; +const ALL_X: [XfK; 4] = [XfK::Identity, XfK::Rotation, XfK::Similarity, XfK::Affine]; +const ALL_S: [ShapeK; 2] = [ShapeK::Regular, ShapeK::Population]; + +/// A hasher that sends every key to one bucket: the map must still answer +/// correctly, by comparing keys. +#[derive(Default, Clone)] +struct Collide; +impl Hasher for Collide { + fn finish(&self) -> u64 { + 0 + } + fn write(&mut self, _: &[u8]) {} +} +impl BuildHasher for Collide { + type Hasher = Collide; + fn build_hasher(&self) -> Collide { + Collide + } +} + +// ───────────────────────── canonical form ───────────────────────── + +/// `p ↦ L p + t`. While the linear part is a similarity it is held as exact +/// `(rot, scale)` and translations stay lazy, each tagged with the rotation +/// and scale in force when it was applied; nothing trigonometric is computed +/// until a terminal asks for it. +#[derive(Clone, Debug)] +struct Canon { + class: XfK, + rot: u32, + scale: f64, + /// Valid when `class == Affine`. + lin: [f64; 4], + /// `(tx, ty, rot_at, k)` while a similarity: `k` is the product of the + /// scales applied after the translation, kept per entry so a zero scale + /// zeroes only the translations it follows. + lazy_t: Vec<(f64, f64, u32, f64)>, + /// Valid when `class == Affine`. + t: (f64, f64), + rewrites: u32, +} + +fn mat_mul(a: [f64; 4], b: [f64; 4]) -> [f64; 4] { + [ + a[0] * b[0] + a[1] * b[2], + a[0] * b[1] + a[1] * b[3], + a[2] * b[0] + a[3] * b[2], + a[2] * b[1] + a[3] * b[3], + ] +} + +fn mat_vec(a: [f64; 4], v: (f64, f64)) -> (f64, f64) { + (a[0] * v.0 + a[1] * v.1, a[2] * v.0 + a[3] * v.1) +} + +fn rot_mat(p: u32, s: f64) -> [f64; 4] { + if p == 0 { + return [s, 0.0, 0.0, s]; + } + let (sn, cs) = sc(rad(p)); + [s * cs, -s * sn, s * sn, s * cs] +} + +impl Canon { + fn new() -> Self { + Canon { + class: XfK::Identity, + rot: 0, + scale: 1.0, + lin: [1.0, 0.0, 0.0, 1.0], + lazy_t: Vec::new(), + t: (0.0, 0.0), + rewrites: 0, + } + } + + fn of(chain: &[Xf]) -> Self { + let mut c = Canon::new(); + for x in chain { + c.push(*x); + } + c.settle(); + c + } + + fn push(&mut self, x: Xf) { + if self.class == XfK::Affine { + match x { + Xf::Rotate(p) => self.affine(rot_mat(p, 1.0)), + Xf::Scale(s) => self.affine([s, 0.0, 0.0, s]), + Xf::Linear(m) => self.affine(m), + Xf::Translate(a, b) => { + self.t = (self.t.0 + a, self.t.1 + b); + } + } + return; + } + match x { + // R(a)R(b) = R(a + b), exact in turns. + Xf::Rotate(p) => { + if self.class != XfK::Identity { + self.rewrites += 1; + } + self.rot = self.rot.wrapping_add(p); + self.class = self.class.max(XfK::Rotation); + } + // A negative similarity is a half turn and a positive scale: the + // rank order of a projection is preserved by construction. + Xf::Scale(s) => { + if s < 0.0 { + self.rot = self.rot.wrapping_add(1 << 31); + self.rewrites += 1; + } + self.scale *= s.abs(); + for e in &mut self.lazy_t { + e.3 *= s.abs(); + } + self.class = self.class.max(XfK::Similarity); + } + Xf::Translate(a, b) => self.lazy_t.push((a, b, self.rot, 1.0)), + Xf::Linear(m) => { + // Leave the similarity form: materialise L and t once. + self.lin = rot_mat(self.rot, self.scale); + self.t = self.lazy_translation(); + self.lazy_t.clear(); + self.class = XfK::Affine; + self.affine(m); + } + } + } + + fn affine(&mut self, m: [f64; 4]) { + self.lin = mat_mul(m, self.lin); + self.t = mat_vec(m, self.t); + } + + /// `R(α)R(−α) = I` and `s = 1`: drop back down the class lattice. + fn settle(&mut self) { + if self.class == XfK::Similarity && self.scale == 1.0 { + self.class = XfK::Rotation; + self.rewrites += 1; + } + if self.class == XfK::Rotation && self.rot == 0 { + self.class = XfK::Identity; + self.rewrites += 1; + } + } + + /// `Σ k_i R(rot − rot_i) t_i`; a term with no relative rotation costs + /// no trigonometry. + fn lazy_translation(&self) -> (f64, f64) { + let mut t = (0.0, 0.0); + for &(a, b, r, k) in &self.lazy_t { + let v = mat_vec(rot_mat(self.rot.wrapping_sub(r), k), (a, b)); + t = (t.0 + v.0, t.1 + v.1); + } + t + } + + fn translation(&self) -> (f64, f64) { + if self.class == XfK::Affine { + self.t + } else { + self.lazy_translation() + } + } + + fn has_translation(&self) -> bool { + if self.class == XfK::Affine { + self.t != (0.0, 0.0) + } else { + self.lazy_t.iter().any(|&(a, b, ..)| a != 0.0 || b != 0.0) + } + } + + /// `‖L‖_F²`, exact for a similarity without trigonometry. + fn frob2(&self) -> f64 { + match self.class { + XfK::Affine => self.lin.iter().map(|v| v * v).sum(), + _ => 2.0 * self.scale * self.scale, + } + } + + /// `(L Lᵀ)₀₁`: zero for every similarity. + fn llt01(&self) -> f64 { + match self.class { + XfK::Affine => self.lin[0] * self.lin[2] + self.lin[1] * self.lin[3], + _ => 0.0, + } + } + + fn apply(&self, p: (f64, f64)) -> (f64, f64) { + let l = match self.class { + XfK::Affine => self.lin, + _ => rot_mat(self.rot, self.scale), + }; + let v = mat_vec(l, p); + let t = self.translation(); + (v.0 + t.0, v.1 + t.1) + } +} + +// ───────────────────────── symmetry-reduced order statistics ───────────────────────── + +/// For `cos(φ + 2πk/m)`, ties occur only at `φ ≡ 0 mod 1/(2m)` turn, so inside +/// each half-sector the ascending order of `k` is fixed. Two permutations +/// replace every sort. +struct SectorPerm { + m: u32, + perm: [Vec; 2], +} + +impl SectorPerm { + fn new(m: u32) -> Self { + let mk = |half: f64| { + let psi = (0.25 + 0.5 * half) / m as f64; + let mut idx: Vec = (0..m).collect(); + idx.sort_by(|&a, &b| { + let va = (TAU * (psi + a as f64 / m as f64)).cos(); + let vb = (TAU * (psi + b as f64 / m as f64)).cos(); + va.total_cmp(&vb) + }); + idx + }; + SectorPerm { + m, + perm: [mk(0.0), mk(1.0)], + } + } + + /// The rank-th smallest of `cos(2π(phi_t + k/m))`, `phi_t` in turns. + #[inline] + fn value(&self, rank: usize, phi_t: f64) -> f64 { + let m = self.m as f64; + let f = (phi_t * m).rem_euclid(1.0); + let j = self.perm[usize::from(f >= 0.5)][rank]; + (TAU * ((f + j as f64) / m)).cos() + } +} + +// ───────────────────────── oracle ───────────────────────── + +fn center_of(s: &Regular) -> (f64, f64) { + match s.center { + Center::Wankel { e } => { + let (sn, cs) = sc(rad(s.phase.wrapping_mul(3))); + (e * cs, e * sn) + } + Center::Fixed(x, y) => (x, y), + } +} + +fn apply_raw(chain: &[Xf], mut p: (f64, f64)) -> (f64, f64) { + for x in chain { + p = match *x { + Xf::Rotate(a) => { + let (sn, cs) = sc(rad(a)); + (cs * p.0 - sn * p.1, sn * p.0 + cs * p.1) + } + Xf::Scale(s) => (s * p.0, s * p.1), + Xf::Translate(a, b) => (p.0 + a, p.1 + b), + Xf::Linear(m) => mat_vec(m, p), + }; + } + p +} + +/// Materialise every vertex, apply the raw chain to each, evaluate. +fn oracle(s: &Regular, chain: &[Xf], t: &Term) -> Answer { + let c = center_of(s); + let pts: Vec<(f64, f64)> = (0..s.m) + .map(|k| { + let (sn, cs) = sc(rad(s.phase) + TAU * k as f64 / s.m as f64); + apply_raw(chain, (c.0 + s.r * cs, c.1 + s.r * sn)) + }) + .collect(); + eval_points(&pts, t) +} + +/// The case's magnitude for [`close`]. +fn magnitude(s: &Regular, chain: &[Xf]) -> f64 { + match oracle(s, chain, &Term::SumSqOrigin) { + Answer::Num(v) => v + s.m as f64, + _ => unreachable!(), + } +} + +fn eval_points(pts: &[(f64, f64)], t: &Term) -> Answer { + let n = pts.len() as f64; + let (sx, sy) = pts.iter().fold((0.0, 0.0), |a, p| (a.0 + p.0, a.1 + p.1)); + let (mx, my) = (sx / n, sy / n); + match *t { + Term::Count => Answer::Num(n), + Term::Sum => Answer::Pt(sx, sy), + Term::Centroid => Answer::Pt(mx, my), + Term::SumSqCentered => Answer::Num( + pts.iter() + .map(|p| (p.0 - mx).powi(2) + (p.1 - my).powi(2)) + .sum(), + ), + Term::SumSqOrigin => Answer::Num(pts.iter().map(|p| p.0 * p.0 + p.1 * p.1).sum()), + Term::CrossCentered => Answer::Num(pts.iter().map(|p| (p.0 - mx) * (p.1 - my)).sum()), + Term::QuantileX(q) => { + let mut xs: Vec = pts.iter().map(|p| p.0).collect(); + xs.sort_by(f64::total_cmp); + Answer::Num(xs[rank_per_10000(xs.len(), q)]) + } + Term::Vertex(k) => { + let p = pts[k as usize]; + Answer::Pt(p.0, p.1) + } + Term::Collide(px, py, d) => Answer::Bool( + pts.iter() + .any(|p| (p.0 - px).powi(2) + (p.1 - py).powi(2) <= d * d), + ), + } +} + +// ───────────────────────── planner + evaluator ───────────────────────── + +#[derive(Debug)] +struct Plan { + rid: Rid, + /// The recipe the table proposed, before guards. + proposed: Rid, + trig: u64, + vertices: u32, +} + +/// The centre's image `L c + t`; for a Wankel centre under a similarity, one +/// phasor of `3θ + rot` and no matrix. +fn center_image(s: &Regular, c: &Canon) -> (f64, f64) { + let lc = match (s.center, c.class) { + (Center::Wankel { e }, k) if k != XfK::Affine => { + if e == 0.0 { + (0.0, 0.0) + } else { + let (sn, cs) = sc(rad(s.phase.wrapping_mul(3).wrapping_add(c.rot))); + (c.scale * e * cs, c.scale * e * sn) + } + } + _ => { + let c0 = center_of(s); + if c.class == XfK::Affine { + mat_vec(c.lin, c0) + } else { + mat_vec(rot_mat(c.rot, c.scale), c0) + } + } + }; + let t = c.translation(); + (lc.0 + t.0, lc.1 + t.1) +} + +/// `|c|` without a phase when the concept carries it. +fn center_norm(s: &Regular) -> f64 { + match s.center { + Center::Wankel { e } => e.abs(), + Center::Fixed(x, y) => x.hypot(y), + } +} + +/// Guarded evaluation. A guard that fails falls through to `Materialize`. +fn plan_eval( + s: &Regular, + chain: &[Xf], + t: &Term, + perms: &mut HashMap, +) -> (Answer, Plan) { + let t0 = trig(); + let c = Canon::of(chain); + let proposed = lookup_array(ShapeK::Regular, c.class, term_k(t)); + let m = s.m as f64; + let r2 = s.r * s.r; + // The centred second moment Σ (z−c)(z−c)ᵀ = (m r²/2)·I holds for m ≥ 3 + // only; m = 2 is a rank-one dyad. + let isotropic = s.m >= 3; + let mut rid = proposed; + let mut vertices = 0; + let ans = match (proposed, t) { + (Rid::Const, Term::Count) => Answer::Num(m), + (Rid::CenterPhasor, Term::Sum) => { + let p = center_image(s, &c); + Answer::Pt(m * p.0, m * p.1) + } + (Rid::CenterPhasor, Term::Centroid) => { + let p = center_image(s, &c); + Answer::Pt(p.0, p.1) + } + (Rid::IsoInvariant, Term::SumSqCentered) if isotropic => { + Answer::Num(0.5 * m * r2 * c.frob2()) + } + (Rid::IsoInvariant, Term::CrossCentered) if isotropic => { + Answer::Num(0.5 * m * r2 * c.llt01()) + } + (Rid::OriginNorm, Term::SumSqOrigin) if isotropic => { + let centred = 0.5 * m * r2 * c.frob2(); + let cn2 = if c.class != XfK::Affine && !c.has_translation() { + (c.scale * center_norm(s)).powi(2) + } else { + let p = center_image(s, &c); + p.0 * p.0 + p.1 * p.1 + }; + Answer::Num(centred + m * cn2) + } + (Rid::SectorQuantile, Term::QuantileX(q)) => { + // Σ-free: x'_k = c'_x + r·(L₀₀ cos φ_k + L₀₁ sin φ_k) + // = c'_x + r·ρ·cos(φ_k − β). + let (rho, shift_t) = if c.class == XfK::Affine { + let (a, b) = (c.lin[0], c.lin[1]); + (a.hypot(b), -b.atan2(a) / TAU) + } else { + (c.scale, c.rot as f64 / TURN) + }; + let cx = center_image(s, &c).0; + let rank = rank_per_10000(s.m as usize, *q); + let sp = perms.entry(s.m).or_insert_with(|| SectorPerm::new(s.m)); + TRIG.fetch_add(1, Relaxed); + let v = sp.value(rank, s.phase as f64 / TURN + shift_t); + Answer::Num(cx + s.r * rho * v) + } + (Rid::VertexPhasor, Term::Vertex(k)) => { + let c0 = center_of(s); + let (sn, cs) = sc(rad(s.phase) + TAU * *k as f64 / m); + let v = c.apply((c0.0 + s.r * cs, c0.1 + s.r * sn)); + Answer::Pt(v.0, v.1) + } + _ => { + rid = Rid::Materialize; + vertices = s.m; + oracle(s, chain, t) + } + }; + ( + ans, + Plan { + rid, + proposed, + trig: trig() - t0, + vertices, + }, + ) +} + +/// Relative `1e-9` of the compared values, plus an absolute floor of `1e-13` +/// of the case's magnitude `Σ|p|² + m` so results that cancel to zero are +/// judged against the size of the terms that cancelled. +fn close(a: Answer, b: Answer, mag: f64) -> bool { + let ok = |x: f64, y: f64| (x - y).abs() <= 1e-9 * x.abs().max(y.abs()) + 1e-13 * mag; + match (a, b) { + (Answer::Num(x), Answer::Num(y)) => ok(x, y), + (Answer::Pt(x0, y0), Answer::Pt(x1, y1)) => ok(x0, x1) && ok(y0, y1), + (Answer::Bool(x), Answer::Bool(y)) => x == y, + _ => false, + } +} + +// ───────────────────────── section 1: the Wankel table ───────────────────────── + +fn wankel_table() { + println!("== 1. Wankel concept + transformations: what each terminal needs =="); + let w = Regular { + m: 3, + r: 100.0, + phase: 0x1234_5678, + center: Center::Wankel { e: 14.0 }, + }; + let chain = [ + Xf::Rotate(0x2000_0000), + Xf::Translate(5.0, -3.0), + Xf::Rotate(0xE000_0000), // inverse of the first: cancels + Xf::Scale(2.0), + Xf::Scale(0.5), + ]; + let mut perms = HashMap::new(); + let queries: [(&str, Term); 9] = [ + ("A count", Term::Count), + ("B sum of positions", Term::Sum), + (" centroid", Term::Centroid), + ("C centred Σ|z−c|²", Term::SumSqCentered), + (" centred Σ(x−x̄)(y−ȳ)", Term::CrossCentered), + ("D Σ|z|² (translated)", Term::SumSqOrigin), + ("E apex 1", Term::Vertex(1)), + (" x-quantile p=0.5", Term::QuantileX(5000)), + ("F collide (110, 0) d 5", Term::Collide(110.0, 0.0, 5.0)), + ]; + println!( + " chain: R(+1/8) · T(5, −3) · R(−1/8) · S(2) · S(0.5) → canonical {:?}, {} rewrites", + Canon::of(&chain).class, + Canon::of(&chain).rewrites + ); + println!( + " {:<26} {:<15} {:<15} {:>5} {:>7} {:>6} oracle trig/vertices", + "query", "proposed", "executed", "trig", "verts", "equal" + ); + for (name, t) in queries { + let (a, p) = plan_eval(&w, &chain, &t, &mut perms); + let o0 = trig(); + let o = oracle(&w, &chain, &t); + let ot = trig() - o0; + let eq = close(a, o, magnitude(&w, &chain)); + assert!(eq, "{name}: recipe {a:?} vs oracle {o:?}"); + println!( + " {:<26} {:<15} {:<15} {:>5} {:>7} {:>6} {ot} / 3", + name, + format!("{:?}", p.proposed), + format!("{:?}", p.rid), + p.trig, + p.vertices, + eq + ); + } + // Untranslated: Σ|z|² = 3 s²(R² + e²) needs no phase at all. + let chain2 = [Xf::Rotate(0x1111_1111), Xf::Scale(-1.5)]; + let (a, p) = plan_eval(&w, &chain2, &Term::SumSqOrigin, &mut perms); + assert!(close( + a, + oracle(&w, &chain2, &Term::SumSqOrigin), + magnitude(&w, &chain2) + )); + assert_eq!(p.trig, 0, "untranslated Σ|z|² must be phase-free"); + println!( + " D' Σ|z|² under R·S(−1.5), no translation: {:?} with {} trig = 3·s²(R²+e²) = {:.3}", + p.rid, + p.trig, + 3.0 * 2.25 * (100.0f64.powi(2) + 14.0f64.powi(2)) + ); +} + +// ───────────────────────── section 2: guards and falsifiers ───────────────────────── + +fn guards() { + println!("\n== 2. guards: every recipe against the oracle, and where a recipe must refuse =="); + let mut rng = Rng(0x5EED); + let mut perms = HashMap::new(); + let mut checked = 0usize; + let mut by_rid: HashMap = HashMap::new(); + for case in 0..4000 { + let m = [2u32, 3, 3, 4, 5, 8, 16, 256][case % 8]; + let center = if case % 3 == 0 { + Center::Fixed(rng.range(-50.0, 50.0), rng.range(-50.0, 50.0)) + } else { + Center::Wankel { + e: rng.range(0.0, 20.0), + } + }; + let s = Regular { + m, + r: rng.range(1.0, 150.0), + phase: rng.next() as u32, + center, + }; + let n = rng.below(5) as usize; + let chain: Vec = (0..n) + .map(|_| match rng.below(5) { + 0 => Xf::Rotate(rng.next() as u32), + 1 => Xf::Scale(rng.range(-2.0, 2.0)), + 4 => Xf::Scale(0.0), + 2 => Xf::Translate(rng.range(-30.0, 30.0), rng.range(-30.0, 30.0)), + _ => Xf::Linear([ + rng.range(-2.0, 2.0), + rng.range(-2.0, 2.0), + rng.range(-2.0, 2.0), + rng.range(-2.0, 2.0), + ]), + }) + .collect(); + let terms = [ + Term::Count, + Term::Sum, + Term::Centroid, + Term::SumSqCentered, + Term::SumSqOrigin, + Term::CrossCentered, + Term::QuantileX(rng.below(10_001) as u32), + Term::Vertex(rng.below(m as u64) as u32), + Term::Collide( + rng.range(-200.0, 200.0), + rng.range(-200.0, 200.0), + rng.range(1.0, 100.0), + ), + ]; + let mag = magnitude(&s, &chain); + for t in terms { + let (a, p) = plan_eval(&s, &chain, &t, &mut perms); + let o0 = trig(); + let o = oracle(&s, &chain, &t); + let ot = trig() - o0; + assert!( + close(a, o, mag), + "case {case}: m {m} {chain:?} {t:?}: {:?} gave {a:?}, oracle {o:?}", + p.rid + ); + let e = by_rid.entry(p.rid).or_default(); + e.0 += 1; + e.1 += p.trig; + e.2 += ot; + checked += 1; + } + } + println!( + " {checked} (concept, chain, terminal) cases equal the oracle (rel. 1e-9 of the values + 1e-13 of Σ|p|² + m)" + ); + // A translation applied after a zero scale survives it. + let z = [Xf::Scale(0.0), Xf::Translate(5.0, -3.0)]; + let w0 = Regular { + m: 3, + r: 100.0, + phase: 0x1234_5678, + center: Center::Wankel { e: 14.0 }, + }; + let (a, _) = plan_eval(&w0, &z, &Term::Centroid, &mut perms); + assert_eq!(a, Answer::Pt(5.0, -3.0)); + assert_eq!(oracle(&w0, &z, &Term::Centroid), Answer::Pt(5.0, -3.0)); + println!(" S(0) · T(5, −3): centroid (5, −3) by recipe and oracle"); + let mut rows: Vec<_> = by_rid.into_iter().collect(); + rows.sort_by_key(|r| format!("{:?}", r.0)); + for (rid, (n, tr, ot)) in rows { + println!( + " {:<15} {:>6} cases trig/case {:>6.2} oracle trig/case {:>7.2} {:?}", + format!("{rid:?}"), + n, + tr as f64 / n as f64, + ot as f64 / n as f64, + exactness(rid) + ); + } + + // m = 2 is a dyad, not isotropic: the invariant must refuse. + let s2 = Regular { + m: 2, + r: 10.0, + phase: 0x0800_0000, + center: Center::Fixed(0.0, 0.0), + }; + let lin = [Xf::Linear([2.0, 0.0, 0.0, 1.0])]; + let iso = 0.5 * 2.0 * 100.0 * Canon::of(&lin).frob2(); + let truth = match oracle(&s2, &lin, &Term::SumSqCentered) { + Answer::Num(v) => v, + _ => unreachable!(), + }; + let (_, p) = plan_eval(&s2, &lin, &Term::SumSqCentered, &mut perms); + assert!((iso - truth).abs() > 1.0, "m = 2 must break isotropy"); + assert_eq!(p.rid, Rid::Materialize); + println!( + " m = 2 under a non-orthogonal map: isotropic formula {iso:.2} vs truth {truth:.2}; recipe refused → {:?}", + p.rid + ); + + // A non-orthogonal map is not a rotation: the invariant changes with it. + let w = Regular { + m: 3, + r: 10.0, + phase: 0, + center: Center::Wankel { e: 1.0 }, + }; + let shear = [Xf::Linear([1.0, 1.5, 0.0, 1.0])]; + let rot_only = 3.0 * 100.0; + let (a, _) = plan_eval(&w, &shear, &Term::SumSqCentered, &mut perms); + let Answer::Num(v) = a else { unreachable!() }; + assert!((v - rot_only).abs() > 1.0); + println!( + " shear [1 1.5; 0 1]: Σ|z−c|² = {v:.2}, not the rotation invariant {rot_only:.2} (class {:?})", + Canon::of(&shear).class + ); + + // Symmetry of the profile: period 1/m turn and reflection, measured. + let mut worst = 0.0f64; + for i in 0..20_000 { + let m = [3u32, 16][i % 2]; + let th = (i as f64 * 0.618_033_988_7).fract(); + let q = |phi: f64| { + let mut v: Vec = (0..m) + .map(|k| (TAU * (phi + k as f64 / m as f64)).cos()) + .collect(); + v.sort_by(f64::total_cmp); + v + }; + let (a, b, c) = (q(th), q(th + 1.0 / m as f64), q(-th)); + for k in 0..m as usize { + worst = worst.max((a[k] - b[k]).abs()).max((a[k] - c[k]).abs()); + } + } + assert!(worst < 1e-12); + println!( + " profile symmetry: Q(θ + 1/m) = Q(θ) = Q(−θ) to {worst:.1e} on 20,000 phases; one half-sector (1/(2m) turn) is the whole domain" + ); + println!( + " note: 1/3 turn is not representable in u32 turns (2^32 mod 3 = {}); threefold symmetry holds in reals, to 1 ulp-turn in u32", + (1u64 << 32) % 3 + ); +} + +// ───────────────────────── section 3: populations and PowerSums ───────────────────────── + +#[derive(Clone, Copy, PartialEq, Eq, Debug, Default)] +struct S2 { + n: i128, + x: i128, + y: i128, + xx: i128, + yy: i128, + xy: i128, +} + +impl From for S2 { + fn from(c: CrossPowerSums) -> Self { + S2 { + n: c.n as i128, + x: c.sum_x as i128, + y: c.sum_y as i128, + xx: c.sum_x_sq as i128, + yy: c.sum_y_sq as i128, + xy: c.sum_xy, + } + } +} + +/// Integer affine `(a b; c d) p + (tx, ty)`. +#[derive(Clone, Copy, Debug)] +struct IAff { + m: [i64; 4], + t: (i64, i64), +} + +impl IAff { + fn then(self, next: IAff) -> IAff { + let a = next.m; + let b = self.m; + IAff { + m: [ + a[0] * b[0] + a[1] * b[2], + a[0] * b[1] + a[1] * b[3], + a[2] * b[0] + a[3] * b[2], + a[2] * b[1] + a[3] * b[3], + ], + t: ( + a[0] * self.t.0 + a[1] * self.t.1 + next.t.0, + a[2] * self.t.0 + a[3] * self.t.1 + next.t.1, + ), + } + } + + #[inline] + fn point(&self, x: i64, y: i64) -> (i64, i64) { + ( + self.m[0] * x + self.m[1] * y + self.t.0, + self.m[2] * x + self.m[3] * y + self.t.1, + ) + } + + /// `S' = A S + n t`, `M' = A M Aᵀ + A S tᵀ + t Sᵀ Aᵀ + n t tᵀ`. + fn summary(&self, s: S2) -> S2 { + let [a, b, c, d] = self.m.map(i128::from); + let (tx, ty) = (i128::from(self.t.0), i128::from(self.t.1)); + let lx = a * s.x + b * s.y; + let ly = c * s.x + d * s.y; + S2 { + n: s.n, + x: lx + s.n * tx, + y: ly + s.n * ty, + xx: a * a * s.xx + 2 * a * b * s.xy + b * b * s.yy + 2 * tx * lx + s.n * tx * tx, + yy: c * c * s.xx + 2 * c * d * s.xy + d * d * s.yy + 2 * ty * ly + s.n * ty * ty, + xy: a * c * s.xx + + (a * d + b * c) * s.xy + + b * d * s.yy + + tx * ly + + ty * lx + + s.n * tx * ty, + } + } +} + +fn fold_ndarray(xs: &[i32], ys: &[i32], mask: &[u64], keys: &[u32]) -> S2 { + let mut out = [CrossPowerSums::default(); 1]; + masked_group_cross_power_sums_i32(mask, keys, xs, ys, &mut out); + out[0].into() +} + +fn full_mask(n: usize) -> Vec { + let mut m = vec![!0u64; n.div_ceil(64)]; + if !n.is_multiple_of(64) { + *m.last_mut().unwrap() = (1u64 << (n % 64)) - 1; + } + m +} + +fn population() { + println!("\n== 3. populations: affine chains over CrossPowerSums =="); + let ts = [ + IAff { + m: [1, 2, 0, 1], + t: (3, -5), + }, + IAff { + m: [0, -1, 1, 0], + t: (7, 2), + }, + IAff { + m: [2, 0, 0, -1], + t: (-4, 9), + }, + ]; + let comp = ts[0].then(ts[1]).then(ts[2]); + println!( + " T3∘T2∘T1 composed once: A = {:?}, t = {:?}", + comp.m, comp.t + ); + println!( + " {:>9} {:>10} {:>10} {:>10} {:>10} {:>10} ns/row, then ns per summary step", + "rows", "A mat+fold", "A2 fused", "B fold+3T", "C fold+1T", "warm 1T" + ); + for n in [1_000usize, 65_536, 1_000_000] { + let mut rng = Rng(0xAFF + n as u64); + let xs: Vec = (0..n).map(|_| rng.below(201) as i32 - 100).collect(); + let ys: Vec = (0..n).map(|_| rng.below(201) as i32 - 100).collect(); + let keys = vec![0u32; n]; + let mask = full_mask(n); + let mut bx = vec![0i32; n]; + let mut by = vec![0i32; n]; + let reps = if n > 100_000 { 7 } else { 31 }; + // A: materialise each step into lanes, then the shipped fold. + let mut a_res = S2::default(); + let ns_a = median_ns(reps, n, || { + bx.copy_from_slice(&xs); + by.copy_from_slice(&ys); + for t in &ts { + for i in 0..n { + let (x, y) = t.point(i64::from(bx[i]), i64::from(by[i])); + bx[i] = x as i32; + by[i] = y as i32; + } + } + a_res = fold_ndarray(&bx, &by, &mask, &keys); + a_res.n as f64 + }); + // A2: no lane, the composed map per row and a scalar fold. + let mut a2 = S2::default(); + let ns_a2 = median_ns(reps, n, || { + let mut s = S2::default(); + for i in 0..n { + let (x, y) = comp.point(i64::from(xs[i]), i64::from(ys[i])); + let (x, y) = (i128::from(x), i128::from(y)); + s.n += 1; + s.x += x; + s.y += y; + s.xx += x * x; + s.yy += y * y; + s.xy += x * y; + } + a2 = s; + s.n as f64 + }); + // B / C: fold X once, transform the summary. + let mut b_res = S2::default(); + let ns_b = median_ns(reps, n, || { + let mut s = fold_ndarray(&xs, &ys, &mask, &keys); + for t in &ts { + s = t.summary(s); + } + b_res = s; + s.n as f64 + }); + let mut c_res = S2::default(); + let ns_c = median_ns(reps, n, || { + c_res = comp.summary(fold_ndarray(&xs, &ys, &mask, &keys)); + c_res.n as f64 + }); + let base = fold_ndarray(&xs, &ys, &mask, &keys); + let ns_w = median_ns(101, 1000, || { + let mut acc = 0i128; + for _ in 0..1000 { + acc = acc.wrapping_add(comp.summary(std::hint::black_box(base)).xy); + } + acc as f64 + }); + assert_eq!(a_res, a2); + assert_eq!(a_res, b_res); + assert_eq!(a_res, c_res); + println!( + " {n:>9} {ns_a:>10.2} {ns_a2:>10.2} {ns_b:>10.2} {ns_c:>10.2} {ns_w:>10.2} all four equal, bitwise (i128)" + ); + } + + // f64: the same identity is not bitwise. + let mut rng = Rng(0xF64); + let n = 65_536; + let pts: Vec<(f64, f64)> = (0..n) + .map(|_| (rng.range(-100.0, 100.0), rng.range(-100.0, 100.0))) + .collect(); + let (a, b, c, d) = ( + 0.3f64.cos() * 1.7, + -(0.3f64.sin()) * 1.7, + 0.3f64.sin() * 1.7, + 0.3f64.cos() * 1.7, + ); + let (tx, ty) = (12.25, -7.5); + let fold = |pts: &mut dyn Iterator| { + pts.fold([0.0f64; 6], |s, (x, y)| { + [ + s[0] + 1.0, + s[1] + x, + s[2] + y, + s[3] + x * x, + s[4] + y * y, + s[5] + x * y, + ] + }) + }; + let mat = fold( + &mut pts + .iter() + .map(|&(x, y)| (a * x + b * y + tx, c * x + d * y + ty)), + ); + let s = fold(&mut pts.iter().copied()); + let (lx, ly) = (a * s[1] + b * s[2], c * s[1] + d * s[2]); + let sum = [ + s[0], + lx + s[0] * tx, + ly + s[0] * ty, + a * a * s[3] + 2.0 * a * b * s[5] + b * b * s[4] + 2.0 * tx * lx + s[0] * tx * tx, + c * c * s[3] + 2.0 * c * d * s[5] + d * d * s[4] + 2.0 * ty * ly + s[0] * ty * ty, + a * c * s[3] + (a * d + b * c) * s[5] + b * d * s[4] + tx * ly + ty * lx + s[0] * tx * ty, + ]; + let rel = (0..6) + .map(|i| (mat[i] - sum[i]).abs() / mat[i].abs().max(1.0)) + .fold(0.0f64, f64::max); + let bitwise = (0..6) + .filter(|&i| mat[i].to_bits() == sum[i].to_bits()) + .count(); + println!( + " f64, 65,536 rows, rotation·1.7 + t: summary path vs materialised fold max rel diff {rel:.1e}; {bitwise}/6 fields bitwise equal (a RealIdentity, not IntExact)" + ); + + // Moments are not sufficient for extrema or collision. + let p1 = [(0i32, 0i32), (3, 0), (3, 0)]; + let p2 = [(1i32, 0i32), (1, 0), (4, 0)]; + let fold3 = |p: &[(i32, i32)]| { + let xs: Vec = p.iter().map(|q| q.0).collect(); + let ys: Vec = p.iter().map(|q| q.1).collect(); + fold_ndarray(&xs, &ys, &full_mask(3), &[0, 0, 0]) + }; + assert_eq!(fold3(&p1), fold3(&p2)); + let max = |p: &[(i32, i32)]| p.iter().map(|q| q.0).max().unwrap(); + let hit = |p: &[(i32, i32)]| p.iter().any(|q| (q.0 - 4).abs() + q.1.abs() == 0); + assert_ne!(max(&p1), max(&p2)); + assert_ne!(hit(&p1), hit(&p2)); + for t in [TermK::Quantile, TermK::Vertex, TermK::Collide] { + for x in ALL_X { + assert_eq!(lookup_array(ShapeK::Population, x, t), Rid::Materialize); + } + } + println!( + " {{0,3,3}} vs {{1,1,4}} on y = 0: identical n, Σx, Σy, Σx², Σy², Σxy ({:?}); max {} vs {}, collide (4,0) {} vs {}; table: Population × {{quantile, vertex, collide}} → Materialize", + fold3(&p1), + max(&p1), + max(&p2), + hit(&p1), + hit(&p2) + ); + + // A fold may be pushed through a transform only if the group key does + // not read a transformed column. + let n = 4096; + let mut rng = Rng(0x6E7); + let xs: Vec = (0..n).map(|_| rng.below(201) as i32 - 100).collect(); + let ys: Vec = (0..n).map(|_| rng.below(201) as i32 - 100).collect(); + let ids: Vec = (0..n as u32).map(|i| i & 1).collect(); + let shift = IAff { + m: [1, 0, 0, 1], + t: (50, 0), + }; + let mask = full_mask(n); + let group = |xs: &[i32], ys: &[i32], keys: &[u32]| { + let mut out = [CrossPowerSums::default(); 2]; + masked_group_cross_power_sums_i32(&mask, keys, xs, ys, &mut out); + [S2::from(out[0]), S2::from(out[1])] + }; + let tx: Vec = xs.iter().map(|x| x + 50).collect(); + let key_x = |xs: &[i32]| xs.iter().map(|&x| u32::from(x >= 0)).collect::>(); + // key reads x, the transform moves x: pushing the fold down is wrong. + let pushed = group(&xs, &ys, &key_x(&xs)).map(|s| shift.summary(s)); + let truth = group(&tx, &ys, &key_x(&tx)); + assert_ne!(pushed, truth); + // key reads an untouched column: pushing down is exact. + let pushed_ok = group(&xs, &ys, &ids).map(|s| shift.summary(s)); + assert_eq!(pushed_ok, group(&tx, &ys, &ids)); + println!( + " group key x ≥ 0 under x + 50: pushed-down summaries ≠ truth (group n {} vs {}); key = id parity: pushed-down = truth", + pushed[1].n, truth[1].n + ); +} + +// ───────────────────────── section 4: projected order statistics ───────────────────────── + +#[derive(Clone, Copy)] +struct Inst { + /// Wankel `e` for m = 3, otherwise the fixed centre's x. + c: f64, + r: f64, + phase: u32, + rank: u32, +} + +/// Arm A: generate every vertex, project, select. +#[inline] +fn arm_a(m: u32, wankel: bool, it: &Inst, buf: &mut [f64; 256]) -> f64 { + let th = rad(it.phase); + let cx = if wankel { + it.c * rad(it.phase.wrapping_mul(3)).cos() + } else { + it.c + }; + let b = &mut buf[..m as usize]; + for (k, v) in b.iter_mut().enumerate() { + *v = cx + it.r * (th + TAU * k as f64 / m as f64).cos(); + } + let (_, v, _) = b.select_nth_unstable_by(it.rank as usize, f64::total_cmp); + *v +} + +/// Arm B, m = 3 only: one `sin_cos(θ)`, the ±120° rotations by constants, +/// `cos 3θ = 4c³ − 3c`, a three-element sorting network. +#[inline] +fn arm_b3(it: &Inst) -> f64 { + const H: f64 = 0.866_025_403_784_438_6; // √3/2 + let (s, c) = rad(it.phase).sin_cos(); + let mut u = [c, -0.5 * c - H * s, -0.5 * c + H * s]; + if u[0] > u[1] { + u.swap(0, 1); + } + if u[1] > u[2] { + u.swap(1, 2); + } + if u[0] > u[1] { + u.swap(0, 1); + } + it.c * (4.0 * c * c * c - 3.0 * c) + it.r * u[it.rank as usize] +} + +/// Arm C: a phase-conditioned profile `table[bucket][rank]` over one sector +/// `ψ ∈ [0, 1/m)`, linearly interpolated. +struct Profile { + m: u32, + b: usize, + t: Vec, +} + +impl Profile { + fn new(m: u32, b: usize) -> Self { + let mut t = vec![0f32; (b + 1) * m as usize]; + for i in 0..=b { + let psi = i as f64 / (b as f64 * m as f64); + let mut v: Vec = (0..m) + .map(|k| (TAU * (psi + k as f64 / m as f64)).cos()) + .collect(); + v.sort_by(f64::total_cmp); + for (k, x) in v.iter().enumerate() { + t[i * m as usize + k] = *x as f32; + } + } + Profile { m, b, t } + } + + #[inline] + fn value(&self, rank: usize, phi_t: f64) -> f64 { + let f = (phi_t * self.m as f64).rem_euclid(1.0) * self.b as f64; + let i = (f as usize).min(self.b - 1); + let w = (f - i as f64) as f32; + let m = self.m as usize; + let (a, b) = (self.t[i * m + rank], self.t[(i + 1) * m + rank]); + f64::from(a + w * (b - a)) + } + + fn bytes(&self) -> usize { + self.t.len() * 4 + } +} + +fn instances(n: usize, m: u32, seed: u64) -> Vec { + let mut rng = Rng(seed); + (0..n) + .map(|_| Inst { + c: if m == 3 { + rng.range(0.0, 20.0) + } else { + rng.range(-50.0, 50.0) + }, + r: rng.range(50.0, 150.0), + phase: rng.next() as u32, + rank: rng.below(u64::from(m)) as u32, + }) + .collect() +} + +fn quantiles() { + println!("\n== 4. one projected order statistic per instance: generate+select vs analytic vs profile LUT vs sector permutation =="); + for m in [3u32, 16, 256] { + let wankel = m == 3; + let sp = SectorPerm::new(m); + let prof = Profile::new(m, 64); + let prof_odd = Profile::new(m, 63); + let center = |it: &Inst| { + if wankel { + it.c * rad(it.phase.wrapping_mul(3)).cos() + } else { + it.c + } + }; + let arm_e = + |it: &Inst| center(it) + it.r * sp.value(it.rank as usize, it.phase as f64 / TURN); + let arm_c = |p: &Profile, it: &Inst| { + center(it) + it.r * p.value(it.rank as usize, it.phase as f64 / TURN) + }; + // Error over a dense phase sweep, every rank, including every tie. + let mut buf = [0f64; 256]; + let (mut eb, mut ec, mut eco, mut ee) = (0f64, 0f64, 0f64, 0f64); + let sweep = 200_003u32; + for i in 0..sweep { + let phase = (u64::from(i) * (1u64 << 32) / u64::from(sweep)) as u32; + // The oracle sorts all m projections once; every rank is checked. + let base = Inst { + c: 7.0, + r: 100.0, + phase, + rank: 0, + }; + let th = rad(phase); + let cx = center(&base); + let sorted = &mut buf[..m as usize]; + for (k, v) in sorted.iter_mut().enumerate() { + *v = cx + base.r * (th + TAU * k as f64 / m as f64).cos(); + } + sorted.sort_by(f64::total_cmp); + for rank in 0..m { + let it = Inst { rank, ..base }; + let o = buf[rank as usize]; + if wankel { + eb = eb.max((arm_b3(&it) - o).abs()); + } + ec = ec.max((arm_c(&prof, &it) - o).abs()); + eco = eco.max((arm_c(&prof_odd, &it) - o).abs()); + ee = ee.max((arm_e(&it) - o).abs()); + } + } + assert!(ee < 1e-9, "sector permutation is exact up to rounding"); + if wankel { + assert!(eb < 1e-9); + } + println!( + " m = {m:<3} max |error| on R = 100 ({sweep} phases × all {m} ranks): B {} C(64 even) {ec:.1e} C(63 odd) {eco:.1e} E {ee:.1e} C table {} B, E table {} B", + if wankel { format!("{eb:.1e}") } else { "—".into() }, + prof.bytes(), + 2 * m as usize * 4 + ); + for n in [65_536usize, 1_000_000] { + let inst = instances(n, m, 0x1A + u64::from(m) + n as u64); + let reps = 5; + let na = if m == 256 { n.min(65_536) } else { n }; + let ns_a = median_ns(reps, na, || { + let mut buf = [0f64; 256]; + inst[..na] + .iter() + .map(|it| arm_a(m, wankel, it, &mut buf)) + .sum() + }); + let ns_b = if wankel { + format!( + "{:.1}", + median_ns(reps, n, || inst.iter().map(arm_b3).sum()) + ) + } else { + "—".into() + }; + let ns_c = median_ns(reps, n, || inst.iter().map(|it| arm_c(&prof, it)).sum()); + let ns_e = median_ns(reps, n, || inst.iter().map(arm_e).sum()); + // D: terminals that need no profile at all. + let ns_d = median_ns(reps, n, || { + inst.iter() + .map(|it| 0.5 * m as f64 * it.r * it.r * 2.0) + .sum() + }); + println!( + " {n:>9} instances ns/instance: A gen+select {ns_a:>7.1} B analytic {ns_b:>6} C LUT {ns_c:>6.1} E sector {ns_e:>6.1} (D Σ|z−c|², no profile: {ns_d:.2})" + ); + } + } +} + +// ───────────────────────── section 5: lookup and rewrite cost ───────────────────────── + +fn lookup_cost() { + println!("\n== 5. recipe lookup and canonicalisation cost =="); + let mut rng = Rng(0x10C); + let keys: Vec<(ShapeK, XfK, TermK)> = (0..4096) + .map(|_| { + ( + ALL_S[rng.below(2) as usize], + ALL_X[rng.below(4) as usize], + ALL_T[rng.below(NT as u64) as usize], + ) + }) + .collect(); + let map: HashMap<(ShapeK, XfK, TermK), Rid> = ALL_S + .iter() + .flat_map(|&s| { + ALL_X.iter().flat_map(move |&x| { + ALL_T + .iter() + .map(move |&t| ((s, x, t), lookup_array(s, x, t))) + }) + }) + .collect(); + let mut bad: HashMap<(ShapeK, XfK, TermK), Rid, Collide> = HashMap::with_hasher(Collide); + for (k, v) in &map { + bad.insert(*k, *v); + } + for &(s, x, t) in &keys { + let a = lookup_array(s, x, t); + assert_eq!(a, lookup_match(s, x, t)); + assert_eq!(a, map[&(s, x, t)]); + assert_eq!( + a, + bad[&(s, x, t)], + "a colliding hasher must still return the right recipe" + ); + } + let n = keys.len(); + let sum = |f: &dyn Fn(ShapeK, XfK, TermK) -> Rid| { + keys.iter() + .map(|&(s, x, t)| f(s, x, t) as usize as f64) + .sum::() + }; + let ns_arr = median_ns(31, n, || { + sum(&|s, x, t| lookup_array(s, std::hint::black_box(x), t)) + }); + let ns_match = median_ns(31, n, || { + sum(&|s, x, t| lookup_match(s, std::hint::black_box(x), t)) + }); + let ns_map = median_ns(31, n, || sum(&|s, x, t| map[&(s, x, t)])); + let ns_bad = median_ns(31, n, || sum(&|s, x, t| bad[&(s, x, t)])); + println!( + " ns/lookup over 4,096 random keys: static array {ns_arr:.2} match {ns_match:.2} HashMap(SipHash) {ns_map:.2} HashMap(all keys collide) {ns_bad:.2}; all agree" + ); + let chain = [ + Xf::Rotate(0x2000_0000), + Xf::Translate(5.0, -3.0), + Xf::Rotate(0xE000_0000), + Xf::Scale(2.0), + Xf::Scale(0.5), + Xf::Rotate(7), + Xf::Translate(1.0, 1.0), + Xf::Scale(-1.0), + ]; + let ns_canon = median_ns(31, 1000, || { + (0..1000) + .map(|_| Canon::of(std::hint::black_box(&chain)).rewrites as f64) + .sum() + }); + let w = Regular { + m: 3, + r: 100.0, + phase: 0x1234_5678, + center: Center::Wankel { e: 14.0 }, + }; + let mut perms = HashMap::new(); + let ns_plan = median_ns(31, 1000, || { + (0..1000) + .map(|_| { + match plan_eval( + &w, + std::hint::black_box(&chain), + &Term::SumSqCentered, + &mut perms, + ) + .0 + { + Answer::Num(v) => v, + _ => 0.0, + } + }) + .sum() + }); + let ns_oracle = median_ns(31, 1000, || { + (0..1000) + .map( + |_| match oracle(&w, std::hint::black_box(&chain), &Term::SumSqCentered) { + Answer::Num(v) => v, + _ => 0.0, + }, + ) + .sum() + }); + println!( + " 8-step chain: canonicalise {ns_canon:.0} ns; canonicalise + lookup + evaluate Σ|z−c|² {ns_plan:.0} ns vs oracle (3 vertices × 8 raw steps) {ns_oracle:.0} ns" + ); +} + +// ───────────────────────── 1-D affine quantile guard ───────────────────────── + +/// `Quantile_p(a·X + b)`: for `a < 0` the order reverses, so the rank-th +/// smallest of `aX + b` is `a · X[n − 1 − rank] + b`. +fn affine_quantile(sorted: &[i64], a: i64, b: i64, per_10000: u32) -> i64 { + let r = rank_per_10000(sorted.len(), per_10000); + let i = if a < 0 { sorted.len() - 1 - r } else { r }; + a * sorted[i] + b +} + +fn quantile_guard() { + let mut rng = Rng(0x0D); + for _ in 0..2000 { + let n = 1 + rng.below(40) as usize; + let mut xs: Vec = (0..n).map(|_| rng.below(1000) as i64 - 500).collect(); + xs.sort_unstable(); + let a = rng.below(9) as i64 - 4; + let b = rng.below(100) as i64 - 50; + let p = rng.below(10_001) as u32; + let mut ys: Vec = xs.iter().map(|x| a * x + b).collect(); + ys.sort_unstable(); + assert_eq!( + affine_quantile(&xs, a, b, p), + ys[rank_per_10000(n, p)], + "a {a} b {b} p {p}" + ); + } + println!("\n== 6. 1-D affine quantile: t + s·Q(rank) for s ≥ 0, t + s·Q(n−1−rank) for s < 0 — 2,000 random cases equal the sorted oracle =="); +} + +fn main() { + println!( + "D-ART-1 algebraic recipe probe avx512f={} avx2={}", + cfg!(target_feature = "avx512f"), + cfg!(target_feature = "avx2") + ); + wankel_table(); + guards(); + population(); + quantiles(); + lookup_cost(); + quantile_guard(); + println!("\nall equalities and falsifiers hold"); +} diff --git a/crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs b/crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs index 7391b9788..98ae3b2e6 100644 --- a/crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs +++ b/crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs @@ -457,6 +457,58 @@ fn arm_g(g: &Grid, k: &Kernel, cx: usize, cy: usize, wk: &mut Work) -> i64 { s } +/// Present-cell crossover between D and F, in cells per window. +/// +/// A pin for this host (D-MHB-1 finding 3: ρ ≈ 0.1 of a 1,009-cell disk), +/// not a derived constant. `crossover_sweep` measures how sensitive arm Q is +/// to it. +const DF_CROSSOVER_CELLS: u32 = 100; + +/// One window row of arm D: the row's present disk cells through the LUT. +#[inline] +fn d_row(k: &Kernel, win: u64, yi: usize, wk: &mut Work) -> i64 { + let dy = yi as i64 - k.r as i64; + let mut m = win & k.disk[yi]; + let mut s = 0i64; + while m != 0 { + let dx = m.trailing_zeros() as i64 - k.r as i64; + m &= m - 1; + wk.geometry += 1; + s += i64::from(k.lut[(dx * dx + dy * dy) as usize]); + } + s +} + +/// The window's words, fetched once, and its exact present-cell count. +#[inline] +fn window_words(g: &Grid, k: &Kernel, cx: usize, cy: usize, wk: &mut Work) -> ([u64; 64], u32) { + let mut w = [0u64; 64]; + let mut n = 0u32; + for (yi, slot) in w.iter_mut().enumerate().take(k.wd()) { + *slot = g.window(cx - k.r, cy + yi - k.r, k.wd()); + wk.popcounts += 1; + n += (*slot & k.disk[yi]).count_ones(); + } + (w, n) +} + +/// Arm Q: one choice per query. The window's exact present-cell count picks +/// D below the crossover and F above it; the words are fetched once and +/// shared by the count and the chosen arm. +fn arm_q(g: &Grid, k: &Kernel, cx: usize, cy: usize, cut: u32, wk: &mut Work) -> i64 { + let (w, n) = window_words(g, k, cx, cy, wk); + if n <= cut { + (0..k.wd()).map(|yi| d_row(k, w[yi], yi, wk)).sum() + } else { + (0..k.wd()) + .map(|yi| { + wk.rows += 1; + f_row(k, w[yi], yi, wk) + }) + .sum() + } +} + /// Rows in the order arm H visits them: largest `|contribution|` bound first. fn h_order(k: &Kernel) -> Vec { let mut o: Vec = (0..k.wd()).collect(); @@ -515,6 +567,62 @@ fn arm_h( unreachable!("with no rows left the bounds are [0, 0], so one of the two tests fired"); } +/// The row evaluator an early exit runs. +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +enum RowEval { + /// F's popcounts on every row (arm H). + F, + /// D's per-point path on every row (arm HD). + D, + /// One choice per query from the window count (arm HQ). + Query(u32), +} + +/// Arms HD and HQ: arm H's exact suffix-bound early exit over any exact row +/// evaluator. The bounds are the template's, so they hold for every +/// evaluator; only the cost of a visited row changes. HQ pays the full +/// window count up front, which the early exit cannot then save. +#[allow(clippy::too_many_arguments)] +fn arm_hx( + g: &Grid, + k: &Kernel, + order: &[usize], + up: &[i64], + lo: &[i64], + cx: usize, + cy: usize, + t: i64, + ev: RowEval, + wk: &mut Work, +) -> (Decision, usize) { + let use_d = match ev { + RowEval::F => false, + RowEval::D => true, + RowEval::Query(cut) => window_words(g, k, cx, cy, wk).1 <= cut, + }; + let mut s = 0i64; + for i in 0..=order.len() { + if s + up[i] < t { + return (Decision::Reject, i); + } + if s + lo[i] >= t { + return (Decision::Accept, i); + } + if i == order.len() { + break; + } + let yi = order[i]; + let win = g.window(cx - k.r, cy + yi - k.r, k.wd()); + wk.rows += 1; + s += if use_d { + d_row(k, win, yi, wk) + } else { + f_row(k, win, yi, wk) + }; + } + unreachable!("with no rows left the bounds are [0, 0], so one of the two tests fired"); +} + // ───────────────────────── harness ───────────────────────── fn centres(g: &Grid, k: &Kernel, n: usize, seed: u64) -> Vec<(usize, usize)> { @@ -629,6 +737,13 @@ fn agreement(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)]) { d, "G differs from D at ({cx},{cy})" ); + for cut in [0, DF_CROSSOVER_CELLS, u32::MAX] { + assert_eq!( + arm_q(g, k, cx, cy, cut, &mut wk), + d, + "Q (cut {cut}) differs from D at ({cx},{cy})" + ); + } let a = arm_a(g, k, cx, cy, &mut wk); err_q += (a - d as f64).abs(); norm += a.abs().max(1.0); @@ -645,7 +760,7 @@ fn agreement(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)]) { } let n = qs.len() as f64; println!( - " {name}: Ascan = D = F = G on every checked centre; |quantised − continuous| mean {:.2} ({:.3} % of |S|) |E − D| mean: {}", + " {name}: Ascan = D = F = G = Q on every checked centre; |quantised − continuous| mean {:.2} ({:.3} % of |S|) |E − D| mean: {}", err_q / n, 100.0 * err_q / norm, rgs.iter() @@ -734,6 +849,11 @@ fn bench(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)], scan_queries: &mut |w, x, y| arm_g(g, k, x, y, w), nq, ); + run( + "Q per-query D/F choice", + &mut |w, x, y| arm_q(g, k, x, y, DF_CROSSOVER_CELLS, w), + nq, + ); // B: the shipped executor. The weight lane is filled per query over the // row band `cy − R ..= cy + R` (full width; zero outside the disk), then @@ -831,6 +951,14 @@ fn early_exit(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)]) { Decision::Reject }; assert_eq!(d, want, "H decided {d:?} at ({cx},{cy}), S {s}, T {t}"); + for ev in [RowEval::F, RowEval::D, RowEval::Query(DF_CROSSOVER_CELLS)] { + let (dx, rx) = arm_hx(g, k, &order, &up, &lo, cx, cy, t, ev, &mut Work::default()); + assert_eq!( + dx, want, + "{ev:?} decided {dx:?} at ({cx},{cy}), S {s}, T {t}" + ); + assert_eq!(rx, r, "{ev:?} visited a different number of rows"); + } accepted += usize::from(d == Decision::Accept); rows.push(r); } @@ -855,13 +983,78 @@ fn early_exit(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)]) { .map(|&(x, y)| arm_f(g, k, x, y, &mut Work::default())) .sum::() }); + let mut ns_x = Vec::new(); + for ev in [RowEval::D, RowEval::Query(DF_CROSSOVER_CELLS)] { + let (_, ns) = time(5, || { + qs.iter() + .map(|&(x, y)| arm_hx(g, k, &order, &up, &lo, x, y, t, ev, &mut Work::default()).1) + .sum::() + }); + ns_x.push(ns / qs.len() as f64); + } + let (_, ns_d) = time(5, || { + qs.iter() + .map(|&(x, y)| arm_d(g, k, x, y, &mut Work::default())) + .sum::() + }); + let (_, ns_q) = time(5, || { + qs.iter() + .map(|&(x, y)| arm_q(g, k, x, y, DF_CROSSOVER_CELLS, &mut Work::default())) + .sum::() + }); + let n = qs.len() as f64; println!( - " median T: H {:.0} ns/query against F {:.0} ns/query", - ns_h / qs.len() as f64, - ns_f / qs.len() as f64 + " median T: H {:.0} HD {:.0} HQ {:.0} ns/query against full F {:.0} D {:.0} Q {:.0}", + ns_h / n, + ns_x[0], + ns_x[1], + ns_f / n, + ns_d / n, + ns_q / n ); } +/// Arm Q's sensitivity to its one pin: the cut swept around the crossover on +/// every density, against the better of D and F measured on the same queries. +fn crossover_sweep(k: &Kernel) { + println!(" cut sweep, 1M grid, 2,048 queries; ns/query (Q / min(D, F))"); + let cuts = [25u32, 50, 100, 200, 400]; + print!(" {:<8}", "ρ"); + for c in cuts { + print!(" {:>12}", format!("cut {c}")); + } + println!(" {:>8} {:>8}", "D", "F"); + for rho in [0.01, 0.05, 0.1, 0.2, 0.5, 0.9] { + let g = random_grid(1024, 1024, rho, 0xC0 + (rho * 100.0) as u64); + let qs = centres(&g, k, 2048, 0xC1); + let med = |f: &dyn Fn() -> i64| { + let mut v: Vec = (0..5).map(|_| time(1, f).1).collect(); + v.sort_by(f64::total_cmp); + v[2] / qs.len() as f64 + }; + let d = med(&|| { + qs.iter() + .map(|&(x, y)| arm_d(&g, k, x, y, &mut Work::default())) + .sum() + }); + let f = med(&|| { + qs.iter() + .map(|&(x, y)| arm_f(&g, k, x, y, &mut Work::default())) + .sum() + }); + print!(" {rho:<8}"); + for c in cuts { + let q = med(&|| { + qs.iter() + .map(|&(x, y)| arm_q(&g, k, x, y, c, &mut Work::default())) + .sum() + }); + print!(" {:>12}", format!("{q:.0} / {:.2}", q / d.min(f))); + } + println!(" {d:>8.0} {f:>8.0}"); + } +} + /// Counterexamples the early exit must survive. fn counterexamples(k: &Kernel) { println!("== counterexamples =="); @@ -1118,6 +1311,9 @@ fn main() { let qs = centres(&g, &k, 4096, 0xB4); early_exit("1M wavefront ring + 1 % noise", &g, &k, &qs); + println!("\n== per-query D/F choice =="); + crossover_sweep(&k); + println!(); counterexamples(&k); println!();