From b44e85ff3b82a41ce7675b36ab0f0bd82e9e9a32 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 8 Oct 2026 18:35:46 +0000 Subject: [PATCH] D-PHT-1: phasors without transcendental calls (probe + board) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Measures native sin_cos against a u32-turns phase LUT, CORDIC and a complex recurrence on two phasor-only observables: Wankel apex coordinates and coherent interference over 1,000 sources and 65,536 detectors. Probe only; no production primitive. - A 2^12 LUT with linear interpolation matches f32 sin_cos accuracy at about 4.5x its speed and 12x f64; 3θ is one wrapping_mul. - CORDIC is 2.6-4.8x slower than native f64 on this CPU. - ndarray vml sin/cos call scalar cos/sin per lane. - Two waves reduce to one cos(Δφ): 2.8x, 6x with the LUT. - The amplitude-bounded early exit equals the full fold on every detector; it saves 8-36x of the terms with heavy-tailed amplitudes and almost nothing with comparable ones. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01EFw2WdKr1oxvaKCJC2ua2R --- .claude/board/STATUS_BOARD.md | 8 + .../entries/2026-10-08-phasor-trig-probe.md | 161 ++++ .claude/board/entries/README.md | 3 +- .../examples/phasor_trig_probe.rs | 756 ++++++++++++++++++ 4 files changed, 927 insertions(+), 1 deletion(-) create mode 100644 .claude/board/entries/2026-10-08-phasor-trig-probe.md create mode 100644 crates/lance-graph-mask-risc/examples/phasor_trig_probe.rs diff --git a/.claude/board/STATUS_BOARD.md b/.claude/board/STATUS_BOARD.md index a103a47a7..969d8dd76 100644 --- a/.claude/board/STATUS_BOARD.md +++ b/.claude/board/STATUS_BOARD.md @@ -1,3 +1,11 @@ +## 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. + +| D-id | scope | status | gate / falsifier | +|---|---|---|---| +| **D-PHT-1** | Phasor evaluation: f64/f32 `sin_cos`, LUT nearest/interpolated, CORDIC, recurrence, ndarray `vml`; two-wave closed form; amplitude-bounded early exit with LUT error and exact fallback | Shipped (probe; LUT + integer phase proposed for ndarray) | Wankel invariants to 1e-11; early exit equals the full fold on every detector; dropping the remaining amplitude, the LUT error or the signed cast each fails | + ## D-RPF — Mask × fold over projections of resident bytes (2026-10-08) Plan: `.claude/plans/2026-10-08-resident-projection-fold-mask-v1.md`. Mask = admissibility, fold = one terminal, CE64 register = one instruction; no materialised hop between them. diff --git a/.claude/board/entries/2026-10-08-phasor-trig-probe.md b/.claude/board/entries/2026-10-08-phasor-trig-probe.md new file mode 100644 index 000000000..a0bd29fb6 --- /dev/null +++ b/.claude/board/entries/2026-10-08-phasor-trig-probe.md @@ -0,0 +1,161 @@ +# 2026-10-08 — Phasors without transcendental calls: sin_cos vs phase LUT vs CORDIC vs recurrence (D-PHT-1) + +**Status:** MEASURED. Probe only, no production primitive. +- **Probe:** `crates/lance-graph-mask-risc/examples/phasor_trig_probe.rs`. +- **Revision:** lance-graph `origin/main` d127f7dc. +- **Host:** Intel Xeon @ 2.8 GHz, 4 vCPU. +- **Build:** `x86-64-v3`, release, `debug = 0`, rustc 1.98.1. +- **Timings:** medians of 5 runs. + +```text +CARGO_PROFILE_RELEASE_DEBUG=0 cargo run --release -p lance-graph-mask-risc --example phasor_trig_probe +``` + +## Question + +Can two observables that need only phasors avoid every transcendental call in the hot path? + +1. **Wankel apex coordinates:** `z_k(θ) = e·e^{i3θ} + R·e^{i(θ + 2πk/3)}`. This is the ideal apex triangle, not the rotor flanks and not sealing. +2. **Coherent interference:** `I(P) = |Σ A_i e^{iφ_i(P)}|²`, with an amplitude-bounded exact early exit. + +The phase is a `u32` in turns, so `3θ` is `phase.wrapping_mul(3)` and wraparound is free. + +## Inventory + +- **MISSING:** CORDIC, a phase LUT, or a SIMD `sin`/`cos` anywhere in ndarray or lance-graph. +- **SHIPPED, but not what it looks like:** `ndarray::hpc::vml::vscos`/`vssin` exist, but they call scalar `cos`/`sin` per lane inside an `F32x16` wrapper (`hpc/vml.rs`). That is not a SIMD transcendental. + +## Wankel invariants (f64 reference) + +All hold, checked on 100,000 angles in ±1000 rad, for (R, e) ∈ {(100, 14), (100, 1e-9), (100, 0), (1, 0.15)}: +- every apex lies on the housing curve at `θ + 2πk/3`; +- the triangle side stays `R√3`; +- the rotor centre stays at distance `e`; +- threefold symmetry holds. + +Maximum error is 5.8e-12 for R = 100. The 3:1 law as `wrapping_mul(3)` equals `3θ mod 2π` on 100,000 phases. + +## Results + +### Wankel apex: 1M random phases + +ns per apex; error in housing units, with R = 100. + +| arm | ns | max error | table | +|---|---|---|---| +| native f64 `sin_cos` ×2 | 66.8 | reference | — | +| native f32 `sin_cos` ×2 | 25.4 | 5.1e-5 | — | +| **L12 LUT, nearest** | **3.1** | 8.7e-2 | 32 KB | +| L16 LUT, nearest | 4.0 | 5.5e-3 | 512 KB | +| **L12i LUT, linear interpolation** | **5.6** | 4.6e-5 | 32 KB | +| L10i LUT, linear interpolation | 5.5 | 5.5e-4 | 8 KB | +| CORDIC Q30, 24 iterations | 258 | 1.4e-5 | 192 B | +| ndarray `vml` (4 calls) | 50.9 | 3.4e-5 | materialises 24 MB of arrays | + +The 64K run gives the same ordering: f64 61.9, f32 25.2, L12 2.5, L12i 5.8, `vml` 45.6. + +### Trajectory, 1M equal steps + +| arm | max error | +|---|---| +| f64 complex recurrence, never renormalised | 4.8e-9 | +| f32 recurrence, never renormalised | 1.95 | +| f32 recurrence, renormalised every 1,024 steps | 1.2e-2 | +| integer `u32` phase accumulator + L12i | no drift by construction; 3.8e-5, the table's own error | + +The integer step differs from the requested Δθ by 5.7e-10 rad. That is a frequency choice, not accumulation. + +### Two waves: 65,536 detectors + +| arm | ns | error / peak | +|---|---|---| +| two `sin_cos` + `\|Σ\|²` | 80.7 | reference | +| one `cos(Δφ)`: `I = A1² + A2² + 2A1A2 cos Δφ` | 28.9 | 1.7e-14 (an identity) | +| one L12i lookup of Δφ | 13.4 | 1.7e-7 | + +### 1,000 sources × 65,536 detectors + +ns per (source, detector) term: + +| arm | ns | max \|ΔI\| / peak I | +|---|---|---| +| native f64 `sin_cos` | 45.9 | reference | +| native f32 `sin_cos` | 36.4 | 3.5e-6 | +| **L12 LUT, nearest** | **9.95** | 4.5e-4 | +| L10i LUT, linear interpolation | 12.9 | 6.0e-6 | +| CORDIC 24 | 156.7 | 6.1e-8 | +| path length only (one `√`) | 2.36 | the floor every arm pays | + +### Amplitude-bounded early exit for `I ≥ T` + +Sources are visited in descending amplitude. The bound is `max(0, |S| − R)² ≤ I ≤ (|S| + R)²`. + +**Correctness, two arms:** +- **f64 arm:** asserted equal to the full fold on every detector. +- **L12 arm:** carries its own per-term error `ε = π/2^12 + 2·f32::EPSILON` in the bound. Any detector still undecided after all terms falls back to the f64 reference. + +**Terms evaluated** (of 1,000; mean, p95 in brackets): + +| T quantile | uniform A ∈ [0.5, 1] | `A_i = i^-1.5` | L12 fallbacks (uniform / heavy) | +|---|---|---|---| +| 50 % | 983 (999) | 122 (550) | 2,679 / 278 | +| 90 % | 968 (996) | 76 (390) | 961 / 160 | +| 99 % | 941 (980) | 28 (100) | 124 / 35 | + +## Falsifiers and disable runs + +Each disable run was made after the commit, turned red, and was restored afterwards. + +| disable | effect | +|---|---| +| early exit ignores the unevaluated amplitude | `exact early exit disagrees at detector 4` | +| LUT bound ignores the table's own error | `decided wrongly at detector 7088` | +| `as u64` instead of `as i64` for a phase difference | two-wave LUT error 2.77 of peak, assertion fires | + +**This probe's own first run had two bugs, both caught by its assertions:** +- `as u64` saturates a negative phase difference to 0, giving a 99 % error. +- Deciding at `R = 0` through `√` then a square disagreed with the full fold at the threshold detector. The fix decides on `re² + im²` exactly as the fold computes it, and absorbs the `√` rounding in a relative margin of 1e-12 before that. + +## Findings + +1. **A phase LUT removes the transcendental cost.** + - L12i matches f32 `sin_cos` accuracy (4.6e-5 vs 5.1e-5 on R = 100) at about a quarter of the time (5.6 vs 25.4 ns), and is about 12× faster than f64. + - Nearest L12 is 20× faster than f64 at 8.7e-4 relative error. + - The `u32` turns representation makes `3θ` one `wrapping_mul` and wraparound free. +2. **CORDIC loses on this CPU.** It is 2.6–4.8× slower than native f64 here (16–30 iterations): a serial, branchy shift-add loop against hardware multiply. It is the right tool only without a multiplier or without table space. +3. **ndarray `vml` is not a SIMD `sin`/`cos`.** It runs as fast as native calls and, used through arrays, materialises its inputs and outputs. +4. **Algebra beats lookup where it applies.** For two waves, the phase difference turns two complex phasors into one cosine (2.8×), and a LUT on Δφ makes it 6×. +5. **The bound early exit is a random-walk problem.** + - With N comparable amplitudes, `|S| ~ √N·A` while the remaining bound is `~N·A`, so the bound only closes at the very end: 94–98 % of the terms are still needed. + - With heavy-tailed amplitudes it saves 8–36× (28–122 of 1,000 terms). + - Whether it pays is a property of the amplitude distribution, not of the bound. +6. **After the LUT, the path length is the next floor.** One `√` per term is 2.36 ns against ~10 ns for the LUT arm. The rest is table access and the multiply-add. +7. **Incremental trajectories:** + - an integer phase accumulator plus a LUT has zero drift by construction; + - an f64 recurrence drifts 4.8e-9 in 1M steps, negligible; + - an f32 recurrence is unusable without renormalisation, and still 1e-2 with it. + +## Recommendations + +| item | verdict | benefit | falsifier | +|---|---|---|---| +| `u32`-turns phase + 2^12 `(cos, sin)` LUT with linear interpolation | **ADOPT** (separate PR, ndarray, exposed through `ndarray::simd`) | ~4.5× vs f32 `sin_cos`, ~12× vs f64, same accuracy as f32 | max error ≤ 5e-5 relative on 1M phases; `3θ = wrapping_mul(3)` exact | +| Closed form for two waves, `cos Δφ` | **ADOPT** where exactly two coherent terms meet | 2.8× (6× with LUT) | identity to 1e-14 | +| Amplitude-bounded early exit | **PROBE**, only for heavy-tailed amplitude sets | 8–36× fewer terms when heavy-tailed; none when uniform | equals the full fold on every detector; dropping R fails | +| LUT + error-inflated bound with exact fallback | **PROBE** | exact decisions on top of the LUT speed | dropping ε fails; fallback count reported | +| Integer phase accumulator for trajectories | **ADOPT**, together with the LUT | no drift | error equals the table bound after 1M steps | +| CORDIC | **REJECT** on CPU | — | 2.6–4.8× slower than native f64 | +| f32 complex recurrence | **REJECT** without renormalisation | — | 1.95 units after 1M steps | +| ndarray `vml` sin/cos as a SIMD path | **REJECT the equivalence** | — | per-lane scalar; no speed-up | + +## Not done here + +- **SIMD gather of the LUT** (AVX2 `vpgatherdd`) was not measured. +- **The Wankel rotor as a moving mask** was not built: containment, `point_inside_rotor`, rotating aperture × wavefront, visibility. +- **VSA bind/bundle against complex multiply/add** was not compared. +- **Log₂(r²) overlap buckets** were not built. D-MHB-1 counterexample 5 already shows that distance buckets cannot see phase, and log buckets are coarser still. + +## OPEN + +- One host class only. +- The LUT's place in ndarray, and whether its interpolation should be `f32` or fixed point, are open. diff --git a/.claude/board/entries/README.md b/.claude/board/entries/README.md index 3ffc42e98..b2b0f1c39 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. -273 entries, 2026-08-06 .. 2026-10-08. +274 entries, 2026-08-06 .. 2026-10-08. | date | entry id | finding | file | |---|---|---|---| @@ -33,6 +33,7 @@ exactly one of them, never both. | 2026-10-08 | `rbac-nested-scope-path` | | [2026-10-08-rbac-nested-scope-path.md](2026-10-08-rbac-nested-scope-path.md) | | 2026-10-08 | `rbac-membership-scope` | | [2026-10-08-rbac-membership-scope.md](2026-10-08-rbac-membership-scope.md) | | 2026-10-08 | `rbac-hotplug-socket` | | [2026-10-08-rbac-hotplug-socket.md](2026-10-08-rbac-hotplug-socket.md) | +| 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 | `hhtl-nars-moore-value-tenants` | | [2026-10-08-hhtl-nars-moore-value-tenants.md](2026-10-08-hhtl-nars-moore-value-tenants.md) | diff --git a/crates/lance-graph-mask-risc/examples/phasor_trig_probe.rs b/crates/lance-graph-mask-risc/examples/phasor_trig_probe.rs new file mode 100644 index 000000000..1488ab2ee --- /dev/null +++ b/crates/lance-graph-mask-risc/examples/phasor_trig_probe.rs @@ -0,0 +1,756 @@ +//! D-PHT-1 — Phasor evaluation without transcendental calls: native +//! `sin_cos` against a phase LUT, CORDIC and a complex recurrence, measured on +//! two observables that need nothing but phasors: +//! +//! 1. **Wankel apex coordinates** `z_k(θ) = e·e^{i3θ} + R·e^{i(θ + 2πk/3)}` +//! (the ideal apex triangle; not the rotor flanks, not sealing). +//! 2. **Coherent interference** `I(P) = |Σ_i A_i e^{iφ_i(P)}|²` over many +//! sources and detectors, with an amplitude-bounded exact early exit: +//! `max(0, |S| − R)² ≤ I ≤ (|S| + R)²`, `R = Σ_{unevaluated} A_i`. +//! +//! The phase is a `u32` in TURNS: one full rotation is `2^32`, so `3θ` is +//! `phase.wrapping_mul(3)` and wraparound is free. +//! +//! Nothing here is a production primitive. Every approximate arm is measured +//! against an `f64` `sin_cos` reference that shares no code with it. Timings +//! are printed, never asserted; invariants and decisions are asserted. +//! +//! ```text +//! CARGO_PROFILE_RELEASE_DEBUG=0 cargo run --release -p lance-graph-mask-risc --example phasor_trig_probe +//! ``` + +// A source index addresses several parallel tables at once (position, +// amplitude, phase offset, suffix bound); an iterator per table would hide +// that the loops walk ONE source list. +#![allow(clippy::needless_range_loop)] + +use std::f64::consts::TAU; +use std::hint::black_box; +use std::time::Instant; + +// ───────────────────────── rng ───────────────────────── + +struct Rng(u64); +impl Rng { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = self.0; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } + fn unit(&mut self) -> f64 { + (self.next() >> 11) as f64 / (1u64 << 53) as f64 + } +} + +const TWO32: f64 = 4_294_967_296.0; + +fn turns_to_rad(p: u32) -> f64 { + f64::from(p) / TWO32 * TAU +} + +// ───────────────────────── the four phasor evaluators ───────────────────────── + +/// A `(cos, sin)` table over `2^bits` phases, `f32`. +struct Lut { + bits: u32, + t: Vec<(f32, f32)>, +} + +impl Lut { + fn new(bits: u32) -> Self { + let n = 1usize << bits; + let t = (0..n) + .map(|i| { + let a = TAU * i as f64 / n as f64; + (a.cos() as f32, a.sin() as f32) + }) + .collect(); + Lut { bits, t } + } + /// Nearest entry. + #[inline] + fn nearest(&self, p: u32) -> (f32, f32) { + let sh = 32 - self.bits; + let i = (p.wrapping_add(1 << (sh - 1)) >> sh) as usize; + self.t[i & (self.t.len() - 1)] + } + /// Linear interpolation between the two neighbouring entries. + #[inline] + fn lerp(&self, p: u32) -> (f32, f32) { + let sh = 32 - self.bits; + let i = (p >> sh) as usize; + let f = (p & ((1 << sh) - 1)) as f32 / (1u32 << sh) as f32; + let m = self.t.len() - 1; + let (a, b) = (self.t[i & m], self.t[(i + 1) & m]); + (a.0 + (b.0 - a.0) * f, a.1 + (b.1 - a.1) * f) + } + fn bytes(&self) -> usize { + self.t.len() * 8 + } + /// A conservative bound on `|lut(p) − e^{ip}|` for `nearest`: half a step + /// of arc plus the `f32` rounding of both components. + fn nearest_bound(&self) -> f64 { + std::f64::consts::PI / (1u64 << self.bits) as f64 + 2.0 * f64::from(f32::EPSILON) + } +} + +/// Circular CORDIC in rotation mode, `i64` fixed point Q30, angle in turns. +struct Cordic { + iters: usize, + /// `atan(2^-i)` in turns · 2^32. + atan: Vec, + /// `K = Π 1/√(1 + 2^-2i)` in Q30: the start vector absorbs the gain. + k_q30: i64, +} + +impl Cordic { + fn new(iters: usize) -> Self { + let atan = (0..iters) + .map(|i| ((2f64.powi(-(i as i32))).atan() / TAU * TWO32).round() as i64) + .collect(); + let k: f64 = (0..iters) + .map(|i| 1.0 / (1.0 + 4f64.powi(-(i as i32))).sqrt()) + .product(); + Cordic { + iters, + atan, + k_q30: (k * (1u64 << 30) as f64).round() as i64, + } + } + /// `(cos, sin)` of a phase in turns, as `f64` from Q30. + #[inline] + fn eval(&self, p: u32) -> (f64, f64) { + // Quadrant from the top two bits; the residual is in [0, π/2), inside + // CORDIC's ~1.74 rad convergence range. + let q = p >> 30; + let mut z = i64::from(p & 0x3FFF_FFFF); + let (mut x, mut y) = (self.k_q30, 0i64); + for i in 0..self.iters { + let (dx, dy) = (y >> i, x >> i); + if z >= 0 { + x -= dx; + y += dy; + z -= self.atan[i]; + } else { + x += dx; + y -= dy; + z += self.atan[i]; + } + } + let (c, s) = match q { + 0 => (x, y), + 1 => (-y, x), + 2 => (-x, -y), + _ => (y, -x), + }; + let s30 = (1u64 << 30) as f64; + (c as f64 / s30, s as f64 / s30) + } +} + +// ───────────────────────── Wankel ───────────────────────── + +struct Wankel { + r: f64, + e: f64, +} + +impl Wankel { + /// Reference apex `k` at rotor angle `θ` (radians), `f64` `sin_cos`. + fn apex_ref(&self, theta: f64, k: usize) -> (f64, f64) { + let (s3, c3) = (3.0 * theta).sin_cos(); + let (s, c) = (theta + TAU * k as f64 / 3.0).sin_cos(); + (self.e * c3 + self.r * c, self.e * s3 + self.r * s) + } + /// The housing curve `R e^{iθ} + e e^{i3θ}`. + fn housing(&self, theta: f64) -> (f64, f64) { + self.apex_ref(theta, 0) + } +} + +/// Invariants of the ideal geometry, on the reference implementation. +fn wankel_invariants() { + println!("== Wankel invariants (f64 reference) =="); + for &(r, e) in &[(100.0, 14.0), (100.0, 1e-9), (100.0, 0.0), (1.0, 0.15)] { + let w = Wankel { r, e }; + let mut rng = Rng(0x3A11); + let (mut on_curve, mut side, mut centre, mut sym) = (0f64, 0f64, 0f64, 0f64); + for _ in 0..100_000 { + // Includes angles far outside [0, 2π) to exercise wrapping. + let th = (rng.unit() - 0.5) * 2000.0; + let a: Vec<(f64, f64)> = (0..3).map(|k| w.apex_ref(th, k)).collect(); + for (k, ak) in a.iter().enumerate() { + // 1. Every apex is on the housing curve at θ + 2πk/3. + let h = w.housing(th + TAU * k as f64 / 3.0); + on_curve = on_curve.max(((ak.0 - h.0).powi(2) + (ak.1 - h.1).powi(2)).sqrt()); + // 2. The apex triangle keeps its side R√3. + let b = a[(k + 1) % 3]; + let d = ((ak.0 - b.0).powi(2) + (ak.1 - b.1).powi(2)).sqrt(); + side = side.max((d - r * 3f64.sqrt()).abs()); + } + // 3. The rotor centre (mean of the apexes) is at distance e. + let c = ( + (a[0].0 + a[1].0 + a[2].0) / 3.0, + (a[0].1 + a[1].1 + a[2].1) / 3.0, + ); + centre = centre.max(((c.0 * c.0 + c.1 * c.1).sqrt() - e).abs()); + // 4. Threefold symmetry: advancing θ by 2π/3 moves apex 0 to apex 1. + let n = w.apex_ref(th + TAU / 3.0, 0); + sym = sym.max(((n.0 - a[1].0).powi(2) + (n.1 - a[1].1).powi(2)).sqrt()); + } + let tol = 1e-9 * r; + assert!( + on_curve < tol && side < tol && centre < tol && sym < tol, + "R {r} e {e}" + ); + println!( + " R {r:>5} e {e:>7}: max off-curve {on_curve:.1e} side error {side:.1e} centre error {centre:.1e} symmetry {sym:.1e}" + ); + } + // The 3:1 law in integer turns is exact: 3·p wraps exactly as 3θ mod 2π. + let mut rng = Rng(0x31); + for _ in 0..100_000 { + let p = rng.next() as u32; + let lhs = turns_to_rad(p.wrapping_mul(3)); + let rhs = (3.0 * turns_to_rad(p)).rem_euclid(TAU); + assert!((lhs - rhs).abs() < 1e-9 || (lhs - rhs).abs() > TAU - 1e-9); + } + println!( + " all hold; the 3:1 law as `phase.wrapping_mul(3)` matches 3θ mod 2π on 100,000 phases" + ); +} + +fn bench_ns(n: usize, reps: usize, mut f: impl FnMut() -> T) -> f64 { + black_box(f()); + let t = Instant::now(); + for _ in 0..reps { + black_box(f()); + } + t.elapsed().as_nanos() as f64 / (reps * n) as f64 +} + +/// Apex coordinates for `n` phases: speed and max error against the reference. +fn wankel_eval(n: usize) { + println!( + "\n== Wankel apex coordinates, {n} random phases (R 100, e 14; error in housing units) ==" + ); + let w = Wankel { r: 100.0, e: 14.0 }; + let mut rng = Rng(0xA9); + let phases: Vec = (0..n).map(|_| rng.next() as u32).collect(); + let reference: Vec<(f64, f64)> = phases + .iter() + .map(|p| w.apex_ref(turns_to_rad(*p), 0)) + .collect(); + let reps = 5; + let mut out = vec![(0f64, 0f64); n]; + + let report = |name: &str, ns: f64, out: &[(f64, f64)], trig: &str, table: usize| { + let err = out + .iter() + .zip(&reference) + .map(|(a, b)| ((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt()) + .fold(0f64, f64::max); + println!(" {name:<30} {ns:>6.2} ns/apex max error {err:>9.2e} trig calls {trig:<8} table {table:>7} B"); + }; + + let ns = bench_ns(n, reps, || { + for (o, p) in out.iter_mut().zip(&phases) { + let th = turns_to_rad(*p); + let (s, c) = th.sin_cos(); + let (s3, c3) = (3.0 * th).sin_cos(); + *o = (w.r * c + w.e * c3, w.r * s + w.e * s3); + } + }); + report("S64 native f64 sin_cos", ns, &out, "2/apex", 0); + + let (r32, e32) = (w.r as f32, w.e as f32); + let ns = bench_ns(n, reps, || { + for (o, p) in out.iter_mut().zip(&phases) { + let th = turns_to_rad(*p) as f32; + let (s, c) = th.sin_cos(); + let (s3, c3) = (3.0 * th).sin_cos(); + *o = (f64::from(r32 * c + e32 * c3), f64::from(r32 * s + e32 * s3)); + } + }); + report("S32 native f32 sin_cos", ns, &out, "2/apex", 0); + + for bits in [10u32, 12, 16] { + let lut = Lut::new(bits); + let ns = bench_ns(n, reps, || { + for (o, p) in out.iter_mut().zip(&phases) { + let (c, s) = lut.nearest(*p); + let (c3, s3) = lut.nearest(p.wrapping_mul(3)); + *o = (f64::from(r32 * c + e32 * c3), f64::from(r32 * s + e32 * s3)); + } + }); + report(&format!("L{bits} LUT nearest"), ns, &out, "0", lut.bytes()); + } + for bits in [8u32, 10, 12] { + let lut = Lut::new(bits); + let ns = bench_ns(n, reps, || { + for (o, p) in out.iter_mut().zip(&phases) { + let (c, s) = lut.lerp(*p); + let (c3, s3) = lut.lerp(p.wrapping_mul(3)); + *o = (f64::from(r32 * c + e32 * c3), f64::from(r32 * s + e32 * s3)); + } + }); + report( + &format!("L{bits}i LUT linear interpolation"), + ns, + &out, + "0", + lut.bytes(), + ); + } + for iters in [16usize, 24, 30] { + let cd = Cordic::new(iters); + let ns = bench_ns(n, reps, || { + for (o, p) in out.iter_mut().zip(&phases) { + let (c, s) = cd.eval(*p); + let (c3, s3) = cd.eval(p.wrapping_mul(3)); + *o = (w.r * c + w.e * c3, w.r * s + w.e * s3); + } + }); + report(&format!("C{iters} CORDIC Q30"), ns, &out, "0", iters * 8); + } + + // ndarray's vector math over arrays: θ and 3θ must be materialised as + // arrays first, and the outputs read back. Its cos/sin are per-lane + // scalar calls inside an F32x16 wrapper (hpc/vml.rs). + { + use ndarray::Array1; + let th: Array1 = phases.iter().map(|p| turns_to_rad(*p) as f32).collect(); + let th3: Array1 = phases + .iter() + .map(|p| turns_to_rad(p.wrapping_mul(3)) as f32) + .collect(); + let (mut c, mut s, mut c3, mut s3) = ( + Array1::::zeros(n), + Array1::::zeros(n), + Array1::::zeros(n), + Array1::::zeros(n), + ); + let ns = bench_ns(n, reps, || { + ndarray::hpc::vml::vscos(th.view(), c.view_mut()); + ndarray::hpc::vml::vssin(th.view(), s.view_mut()); + ndarray::hpc::vml::vscos(th3.view(), c3.view_mut()); + ndarray::hpc::vml::vssin(th3.view(), s3.view_mut()); + for (i, o) in out.iter_mut().enumerate() { + *o = ( + f64::from(r32 * c[i] + e32 * c3[i]), + f64::from(r32 * s[i] + e32 * s3[i]), + ); + } + }); + report("V ndarray vml vscos/vssin", ns, &out, "4/apex", 0); + println!( + " (V also materialises θ, 3θ and four output arrays: {} B)", + 6 * n * 4 + ); + } +} + +/// A trajectory of `steps` equal increments: the complex recurrence against +/// the reference, with and without renormalisation. +fn recurrence(steps: usize) { + println!("\n== Wankel trajectory, {steps} steps of Δθ = 2π/4096·φ (recurrence; error at the last step and max) =="); + let w = Wankel { r: 100.0, e: 14.0 }; + // An irrational step, so the trajectory never revisits a phase exactly. + let dth = TAU / 4096.0 * 1.618_033_988_749_895; + let reference = |n: usize| w.apex_ref(n as f64 * dth, 0); + + // f64 and f32 recurrences, renormalising every `every` steps (0 = never). + for (label, every) in [("never", 0usize), ("every 1024", 1024), ("every 64", 64)] { + let (w1, w3) = ( + (dth.cos(), dth.sin()), + ((3.0 * dth).cos(), (3.0 * dth).sin()), + ); + let (mut u, mut v) = ((1.0f64, 0.0f64), (1.0f64, 0.0f64)); + let (mut u32_, mut v32) = ((1.0f32, 0.0f32), (1.0f32, 0.0f32)); + let (w1f, w3f) = ((w1.0 as f32, w1.1 as f32), (w3.0 as f32, w3.1 as f32)); + let (mut e64, mut e32) = (0f64, 0f64); + let mul = |a: (f64, f64), b: (f64, f64)| (a.0 * b.0 - a.1 * b.1, a.0 * b.1 + a.1 * b.0); + let mulf = |a: (f32, f32), b: (f32, f32)| (a.0 * b.0 - a.1 * b.1, a.0 * b.1 + a.1 * b.0); + let t = Instant::now(); + for n in 1..=steps { + u = mul(u, w1); + v = mul(v, w3); + u32_ = mulf(u32_, w1f); + v32 = mulf(v32, w3f); + if every != 0 && n % every == 0 { + let m = (u.0 * u.0 + u.1 * u.1).sqrt(); + u = (u.0 / m, u.1 / m); + let m = (v.0 * v.0 + v.1 * v.1).sqrt(); + v = (v.0 / m, v.1 / m); + let m = (u32_.0 * u32_.0 + u32_.1 * u32_.1).sqrt(); + u32_ = (u32_.0 / m, u32_.1 / m); + let m = (v32.0 * v32.0 + v32.1 * v32.1).sqrt(); + v32 = (v32.0 / m, v32.1 / m); + } + if n % 997 == 0 || n == steps { + let r = reference(n); + let z = (w.r * u.0 + w.e * v.0, w.r * u.1 + w.e * v.1); + e64 = e64.max(((z.0 - r.0).powi(2) + (z.1 - r.1).powi(2)).sqrt()); + let zf = ( + w.r * f64::from(u32_.0) + w.e * f64::from(v32.0), + w.r * f64::from(u32_.1) + w.e * f64::from(v32.1), + ); + e32 = e32.max(((zf.0 - r.0).powi(2) + (zf.1 - r.1).powi(2)).sqrt()); + } + } + let ns = t.elapsed().as_nanos() as f64 / steps as f64; + println!( + " renormalise {label:<10}: max error f64 {e64:>9.2e} f32 {e32:>9.2e} ({ns:.2} ns/step for both, incl. the sampled reference)" + ); + } + // The same trajectory as an integer phase: an exact u32 step accumulates + // no drift at all; only the LUT's own error remains. + let step = (dth / TAU * TWO32).round() as u32; + let lut = Lut::new(12); + let (mut p, mut emax) = (0u32, 0f64); + let step_err = (f64::from(step) / TWO32 * TAU - dth).abs(); + for n in 1..=steps { + p = p.wrapping_add(step); + if n % 997 == 0 || n == steps { + // Compare against the trajectory the integer step actually takes. + let th = turns_to_rad(p); + let (c, s) = lut.lerp(p); + let (c3, s3) = lut.lerp(p.wrapping_mul(3)); + let z = ( + w.r * f64::from(c) + w.e * f64::from(c3), + w.r * f64::from(s) + w.e * f64::from(s3), + ); + let r = w.apex_ref(th, 0); + emax = emax.max(((z.0 - r.0).powi(2) + (z.1 - r.1).powi(2)).sqrt()); + } + } + println!( + " integer phase + L12i : max error {emax:.2e} against its own trajectory; the u32 step differs from Δθ by {step_err:.1e} rad per step (a frequency choice, not drift)" + ); +} + +// ───────────────────────── interference ───────────────────────── + +struct Field { + /// Source positions, amplitudes (sorted descending), phase offsets in turns. + sx: Vec, + sy: Vec, + a: Vec, + p0: Vec, + /// 1 / λ. + inv_lambda: f64, +} + +fn field(n: usize, heavy_tail: bool, seed: u64) -> Field { + let mut r = Rng(seed); + let mut a: Vec = (0..n) + .map(|i| { + if heavy_tail { + ((i + 1) as f64).powf(-1.5) + } else { + 0.5 + 0.5 * r.unit() + } + }) + .collect(); + a.sort_by(|x, y| y.partial_cmp(x).unwrap()); + Field { + sx: (0..n).map(|_| -64.0 + 384.0 * r.unit()).collect(), + sy: (0..n).map(|_| -64.0 + 384.0 * r.unit()).collect(), + a, + p0: (0..n).map(|_| r.unit()).collect(), + inv_lambda: 1.0 / 7.3, + } +} + +/// Phase of source `i` at detector `(x, y)`, in turns (`f64`, unwrapped). +#[inline] +fn turns(f: &Field, i: usize, x: f64, y: f64) -> f64 { + let (dx, dy) = (x - f.sx[i], y - f.sy[i]); + (dx * dx + dy * dy).sqrt() * f.inv_lambda + f.p0[i] +} + +/// Turns to a `u32` phase, modulo one turn. Through `i64`, not `u64`: a +/// negative turn count (a phase DIFFERENCE) saturates to 0 under `as u64`, +/// which the first run of this probe did, giving a 99 % error on two waves. +#[inline] +fn to_u32(t: f64) -> u32 { + (t * TWO32) as i64 as u32 +} + +/// The bound tests for `I ≥ t`, given the partial sum `(re, im)` and an upper +/// bound `r` on what the unevaluated terms (plus any evaluation error) can +/// still add. `None` while undecided. +/// +/// With nothing left (`r == 0`) the decision is `re² + im² ≥ t` computed +/// exactly as the full fold computes it, never through `√` then a square. +/// Before that, the `√` rounding is absorbed by a relative margin. +#[inline] +fn bound_decide(re: f64, im: f64, r: f64, t: f64) -> Option { + let i2 = re * re + im * im; + if r == 0.0 { + return Some(i2 >= t); + } + let m = i2.sqrt(); + const MARGIN: f64 = 1e-12; + if (m + r) * (m + r) < t * (1.0 - MARGIN) { + Some(false) + } else if (m - r).max(0.0).powi(2) >= t * (1.0 + MARGIN) { + Some(true) + } else { + None + } +} + +/// Full fold at one detector with the given phasor evaluator. +#[inline] +fn fold(f: &Field, x: f64, y: f64, ev: &impl Fn(f64) -> (f64, f64)) -> f64 { + let (mut re, mut im) = (0.0, 0.0); + for i in 0..f.a.len() { + let (c, s) = ev(turns(f, i, x, y)); + re += f.a[i] * c; + im += f.a[i] * s; + } + re * re + im * im +} + +fn interference(n_src: usize, side: usize, heavy_tail: bool) { + let f = field(n_src, heavy_tail, if heavy_tail { 0x7A } else { 0x7B }); + let det: Vec<(f64, f64)> = (0..side * side) + .map(|i| ((i % side) as f64, (i / side) as f64)) + .collect(); + let law = if heavy_tail { + "A_i = i^-1.5" + } else { + "A_i uniform in [0.5, 1]" + }; + println!( + "\n== interference: {n_src} coherent sources, {} detectors, {law}, λ 7.3 ==", + det.len() + ); + + let reference = |t: f64| { + let (s, c) = (TAU * t).sin_cos(); + (c, s) + }; + let t0 = Instant::now(); + let truth: Vec = det + .iter() + .map(|&(x, y)| fold(&f, x, y, &reference)) + .collect(); + let ns_ref = t0.elapsed().as_nanos() as f64 / (det.len() * n_src) as f64; + let peak = truth.iter().cloned().fold(0f64, f64::max); + println!(" S64 native f64 sin_cos {ns_ref:>6.2} ns/term (reference)"); + + let err = |v: &[f64]| { + v.iter() + .zip(&truth) + .map(|(a, b)| (a - b).abs()) + .fold(0f64, f64::max) + / peak + }; + let time_fold = |name: &str, ev: &dyn Fn(f64) -> (f64, f64)| { + let t = Instant::now(); + let v: Vec = det.iter().map(|&(x, y)| fold(&f, x, y, &ev)).collect(); + let ns = t.elapsed().as_nanos() as f64 / (det.len() * n_src) as f64; + println!( + " {name:<28} {ns:>6.2} ns/term max |ΔI| / peak I {:.2e}", + err(&v) + ); + }; + time_fold("S32 native f32 sin_cos", &|t: f64| { + let (s, c) = ((TAU * t) as f32).sin_cos(); + (f64::from(c), f64::from(s)) + }); + let l12 = Lut::new(12); + time_fold("L12 LUT nearest", &|t: f64| { + let (c, s) = l12.nearest(to_u32(t)); + (f64::from(c), f64::from(s)) + }); + let l10 = Lut::new(10); + time_fold("L10i LUT linear interpolation", &|t: f64| { + let (c, s) = l10.lerp(to_u32(t)); + (f64::from(c), f64::from(s)) + }); + let cd = Cordic::new(24); + time_fold("C24 CORDIC Q30", &|t: f64| cd.eval(to_u32(t))); + // The path length alone: the square root every arm shares. + let t = Instant::now(); + let mut acc = 0.0; + for &(x, y) in &det { + for i in 0..n_src { + acc += turns(&f, i, x, y); + } + } + black_box(acc); + let ns = t.elapsed().as_nanos() as f64 / (det.len() * n_src) as f64; + println!(" geometry only (√ per term) {ns:>6.2} ns/term (the floor every arm pays)"); + + // Early exit for I ≥ T, sources in descending amplitude. + let suffix: Vec = { + let mut s = vec![0.0; n_src + 1]; + for i in (0..n_src).rev() { + s[i] = s[i + 1] + f.a[i]; + } + s + }; + let mut sorted = truth.clone(); + sorted.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let eps = l12.nearest_bound(); + for q in [0.5, 0.9, 0.99] { + let t_thr = sorted[((sorted.len() - 1) as f64 * q) as usize]; + // Exact path: f64 reference phasors. + let mut used = Vec::with_capacity(det.len()); + // LUT path: per-term error ε inflates the bound; undecided at the end + // falls back to the reference for that detector. + let mut used_lut = Vec::with_capacity(det.len()); + let mut fallback = 0usize; + for (k, &(x, y)) in det.iter().enumerate() { + let want = truth[k] >= t_thr; + // f64. + let (mut re, mut im) = (0.0f64, 0.0f64); + let mut decided = None; + for i in 0..=n_src { + if let Some(d) = bound_decide(re, im, suffix[i], t_thr) { + decided = Some((d, i)); + break; + } + if i == n_src { + break; + } + let (c, s) = reference(turns(&f, i, x, y)); + re += f.a[i] * c; + im += f.a[i] * s; + } + let (d, i) = decided.expect("with nothing left the bounds are the value itself"); + assert_eq!(d, want, "exact early exit disagrees at detector {k}"); + used.push(i); + // LUT with conservative error. + let (mut re, mut im, mut aerr) = (0.0f64, 0.0f64, 0.0f64); + let mut decided = None; + for i in 0..=n_src { + // The LUT sum is NOT the reference value even with nothing + // left: its own error `aerr` keeps the bound open. The floor + // keeps `bound_decide` off its exact `r == 0` branch. + let r = (suffix[i] + aerr).max(f64::MIN_POSITIVE); + if let Some(d) = bound_decide(re, im, r, t_thr) { + decided = Some((d, i)); + break; + } + if i == n_src { + break; + } + let (c, s) = l12.nearest(to_u32(turns(&f, i, x, y))); + re += f.a[i] * f64::from(c); + im += f.a[i] * f64::from(s); + aerr += f.a[i] * eps; + } + match decided { + Some((d, i)) => { + assert_eq!( + d, want, + "LUT early exit with error bound decided wrongly at detector {k}" + ); + used_lut.push(i); + } + None => { + fallback += 1; + used_lut.push(n_src); + } + } + } + let mean = |v: &[usize]| v.iter().sum::() as f64 / v.len() as f64; + let p95 = |v: &mut Vec| { + v.sort_unstable(); + v[(v.len() - 1) * 95 / 100] + }; + println!( + " early exit, T at the {:>2.0} % quantile: f64 terms mean {:>6.1} p95 {:>4} | L12 + error bound terms mean {:>6.1} p95 {:>4}, fallbacks {} of {}", + q * 100.0, + mean(&used), + p95(&mut used), + mean(&used_lut), + p95(&mut used_lut), + fallback, + det.len() + ); + } +} + +/// Two waves: the closed form needs one cosine of the phase DIFFERENCE, +/// `I = A1² + A2² + 2·A1·A2·cos(Δφ)`, instead of two complex phasors. +fn two_waves(side: usize) { + println!("\n== two waves: two phasors against one cos(Δφ) =="); + let f = field(2, false, 0x22); + let det: Vec<(f64, f64)> = (0..side * side) + .map(|i| ((i % side) as f64, (i / side) as f64)) + .collect(); + let (a1, a2) = (f.a[0], f.a[1]); + let lut = Lut::new(12); + let n = det.len(); + let mut out = vec![0f64; n]; + let ns_full = bench_ns(n, 5, || { + for (o, &(x, y)) in out.iter_mut().zip(&det) { + let (s1, c1) = (TAU * turns(&f, 0, x, y)).sin_cos(); + let (s2, c2) = (TAU * turns(&f, 1, x, y)).sin_cos(); + let (re, im) = (a1 * c1 + a2 * c2, a1 * s1 + a2 * s2); + *o = re * re + im * im; + } + }); + let truth = out.clone(); + let ns_diff = bench_ns(n, 5, || { + for (o, &(x, y)) in out.iter_mut().zip(&det) { + let d = turns(&f, 0, x, y) - turns(&f, 1, x, y); + *o = a1 * a1 + a2 * a2 + 2.0 * a1 * a2 * (TAU * d).cos(); + } + }); + let e_diff = out + .iter() + .zip(&truth) + .map(|(a, b)| (a - b).abs()) + .fold(0f64, f64::max); + let ns_lut = bench_ns(n, 5, || { + for (o, &(x, y)) in out.iter_mut().zip(&det) { + let d = turns(&f, 0, x, y) - turns(&f, 1, x, y); + let (c, _) = lut.lerp(to_u32(d)); + *o = a1 * a1 + a2 * a2 + 2.0 * a1 * a2 * f64::from(c); + } + }); + let e_lut = out + .iter() + .zip(&truth) + .map(|(a, b)| (a - b).abs()) + .fold(0f64, f64::max); + let peak = (a1 + a2) * (a1 + a2); + // The closed form is an identity; the lookup only adds the table's error. + assert!( + e_diff / peak < 1e-12 && e_lut / peak < 1e-5, + "two waves: {e_diff} / {e_lut}" + ); + println!(" two sin_cos + |Σ|² {ns_full:>6.2} ns/detector (reference)"); + println!( + " one cos(Δφ), closed form {ns_diff:>6.2} ns/detector max error {:.1e} of peak", + e_diff / peak + ); + println!( + " one L12i lookup of Δφ {ns_lut:>6.2} ns/detector max error {:.1e} of peak", + e_lut / peak + ); +} + +fn main() { + println!( + "D-PHT-1 phasor trig probe avx512f={} avx2={}", + cfg!(target_feature = "avx512f"), + cfg!(target_feature = "avx2") + ); + wankel_invariants(); + wankel_eval(65_536); + wankel_eval(1_000_000); + recurrence(1_000_000); + two_waves(256); + interference(1000, 256, false); + interference(1000, 256, true); + println!("\nall invariants and decisions hold"); +}