From a8ebb4b84771dbc551b7dadbe9504774e45a7262 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 4 Jul 2026 10:06:44 +0000 Subject: [PATCH 1/3] =?UTF-8?q?examples:=20gridlake=20field=20interdepende?= =?UTF-8?q?ncy=20=3D=20one=20BF16=2016=C3=9716=20AMX=20tile=20op=20?= =?UTF-8?q?=E2=80=94=20measured?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The operator's claim: a gridlake's cell interdependency, treated as a field operation, reduces to a BF16 16×16 AMX tile GEMM (TDPBF16PS) — the same shape as the intercorrelation / perturbation in a Gaussian splat (the EWA covariance sandwich Σ' = Mᵀ·Σ·M, shipped as hpc::splat3d::spd3::sandwich_x16). Mechanically: a field operation whose output cell is a linear combination of input cells is a linear operator C = A·B — a GEMM; the interdependency IS the coupling matrix. At the gridlake tile granularity (16×16 — a quadrant of the 64×64 gridlake) that GEMM is exactly one hpc::bf16_tile_gemm::bf16_tile_gemm_16x16 call = one TDPBF16PS tile op. Same primitive as the codec separable transform (#232: M·X, WHT/DCT on a tile) and the splat covariance projection (sandwich_x16). Measured on the shipped kernel (this host: AMX TDPBF16PS, amx_available = true), vs a direct f64 reference. Precision reported as Frobenius-relative ‖tile−direct‖_F/‖direct‖_F (the honest aggregate — a per-element max-relative blows up on the near-zero entries a cancellation-heavy Mᵀ·Σ·M naturally has): (1) intercorrelation C = Xᵀ·X frobenius_rel_err = 0.095% (one tile op) (2) covariance sandwich Σ' = Mᵀ·Σ·M frobenius_rel_err = 0.401% (two tile ops) (3) bit-exact on bf16-exact integer operands (|Σ| < 2^24): true (4) throughput: 1.22 M tile-ops/s Conclusion: YES — the gridlake field interdependency reduces to one BF16 16×16 tile op, matching the direct reference to BF16 precision and bit-exact for bf16-exact integer operands. Interdependency / intercorrelation / perturbation / codec transform / covariance projection are ONE op: the 16×16 BF16 tile GEMM. The 64×64 gridlake is 4×4 = 16 of these tiles. fmt + clippy clean, gated required-features = ["std"]. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01MLBnPuScZy6w9di2QEjsXM --- Cargo.toml | 4 + examples/gridlake_field_tile.rs | 232 ++++++++++++++++++++++++++++++++ 2 files changed, 236 insertions(+) create mode 100644 examples/gridlake_field_tile.rs diff --git a/Cargo.toml b/Cargo.toml index 44f8aa02..0bc7b651 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -95,6 +95,10 @@ required-features = ["codec"] name = "motion_transform_noop" required-features = ["codec"] +[[example]] +name = "gridlake_field_tile" +required-features = ["std"] + [[example]] name = "entropy_ladder_probe" required-features = ["std"] diff --git a/examples/gridlake_field_tile.rs b/examples/gridlake_field_tile.rs new file mode 100644 index 00000000..5e948bea --- /dev/null +++ b/examples/gridlake_field_tile.rs @@ -0,0 +1,232 @@ +//! Gridlake field interdependency = one BF16 16×16 AMX tile op — measured. +//! +//! The operator's claim (2026-07-04): a **gridlake**'s cell *interdependency*, +//! treated as a *field operation*, reduces to a **BF16 16×16 AMX tile GEMM** +//! (`TDPBF16PS`) — the same shape as the **intercorrelation / perturbation in a +//! Gaussian splat** (the EWA covariance sandwich `Σ' = Mᵀ·Σ·M`, shipped as +//! `hpc::splat3d::spd3::sandwich_x16`). +//! +//! Why it's true, mechanically: a *field operation* whose output cell is a linear +//! combination of input cells is a linear operator `C = A·B` — a GEMM. The +//! "interdependency" / "intercorrelation" IS the coupling matrix (the off-diagonal +//! second-moment structure). At the gridlake tile granularity (16×16, one quarter +//! of the 64×64 gridlake in each axis → the AMX tile) that GEMM is exactly one +//! `hpc::bf16_tile_gemm::bf16_tile_gemm_16x16` call — one `TDPBF16PS` tile op. +//! This is the SAME primitive as: +//! • the codec separable transform (#232: `M·X`, the WHT/DCT on a tile), and +//! • the splat EWA covariance projection (`sandwich_x16`, `Mᵀ·Σ·M`). +//! One op, three names. +//! +//! This probe MEASURES the reduction on the shipped kernel (whatever tier the host +//! selects — AMX `TDPBF16PS` / AVX-512 `VDPBF16PS` / F32x16 polyfill), reports the +//! tier, and checks against a direct `f64` reference: +//! 1. **intercorrelation** `C = Xᵀ·X` of a gridlake field tile (the inter-cell +//! coupling as a field op) — one tile GEMM. +//! 2. **covariance sandwich** `Σ' = Mᵀ·Σ·M` (the splat EWA projection shape) — +//! two tile GEMMs. +//! 3. **bit-exactness** on bf16-exact integer operands (the module's guarantee), +//! and BF16-precision relative error on a continuous field. +//! 4. **throughput** of the 16×16 tile op on this host. +//! +//! Run: `cargo run --release --example gridlake_field_tile --features std` + +use ndarray::hpc::amx_matmul::amx_available; +use ndarray::hpc::bf16_tile_gemm::{bf16_tile_gemm_16x16, bf16_tile_gemm_tier}; +use ndarray::hpc::quantized::{f32_to_bf16_rounded, BF16}; + +const M: usize = 16; // tile edge — the AMX tile / gridlake quadrant +const KPAD: usize = 32; // bf16_tile_gemm requires K a multiple of 32; pad 16→32 with zeros + +fn mix(mut z: u64) -> u64 { + z = z.wrapping_add(0x9E37_79B9_7F4A_7C15); + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) +} + +/// Pack a row-major `M×KPAD` (or `KPAD×M`) f32 matrix into the `&[u16]` bf16 +/// operand the tile GEMM consumes (RNE, matching `VCVTNEPS2BF16`). +fn to_bf16(src: &[f32]) -> Vec { + let mut b = vec![BF16(0); src.len()]; + f32_to_bf16_rounded(src, &mut b); + b.iter().map(|x| x.0).collect() +} + +/// `C[16×16] = A[16×16] · B[16×16]` via the SHIPPED BF16 tile op. A/B are padded +/// on K from 16→32 with zeros (the pad contributes nothing), so this is one +/// `bf16_tile_gemm_16x16` — one `TDPBF16PS` tile op on the AMX tier. +fn tile_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f32; M * M] { + let mut a_pad = [0.0f32; M * KPAD]; + let mut b_pad = [0.0f32; KPAD * M]; + for i in 0..M { + for k in 0..M { + a_pad[i * KPAD + k] = a[i * M + k]; // A[16×32], cols 16..31 = 0 + b_pad[k * M + i] = b[k * M + i]; // B[32×16], rows 16..31 = 0 (already) + } + } + let (ab, bb) = (to_bf16(&a_pad), to_bf16(&b_pad)); + let mut c = vec![0.0f32; M * M]; + bf16_tile_gemm_16x16(&ab, &bb, &mut c, KPAD); + let mut out = [0.0f32; M * M]; + out.copy_from_slice(&c); + out +} + +/// Direct `f64` reference GEMM `C = A·B` (16×16). +fn direct_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f64; M * M] { + let mut c = [0.0f64; M * M]; + for i in 0..M { + for j in 0..M { + let mut s = 0.0f64; + for k in 0..M { + s += a[i * M + k] as f64 * b[k * M + j] as f64; + } + c[i * M + j] = s; + } + } + c +} + +fn transpose(a: &[f32; M * M]) -> [f32; M * M] { + let mut t = [0.0f32; M * M]; + for i in 0..M { + for j in 0..M { + t[j * M + i] = a[i * M + j]; + } + } + t +} + +/// `(max_abs_err, frobenius_relative_err)` of a tile-GEMM result vs the f64 +/// direct. Frobenius-relative `‖tile − direct‖_F / ‖direct‖_F` is the honest +/// aggregate precision — a per-element max-relative blows up on the near-zero +/// entries that a cancellation-heavy product (e.g. `Mᵀ·Σ·M`) naturally has, +/// which measures the cancellation, not the kernel. +fn errors(tile: &[f32; M * M], direct: &[f64; M * M]) -> (f64, f64) { + let mut max_abs = 0.0f64; + let mut num = 0.0f64; // ‖E‖_F² + let mut den = 0.0f64; // ‖direct‖_F² + for i in 0..M * M { + let e = tile[i] as f64 - direct[i]; + max_abs = max_abs.max(e.abs()); + num += e * e; + den += direct[i] * direct[i]; + } + (max_abs, (num / den.max(1e-12)).sqrt()) +} + +fn main() { + println!("Gridlake field interdependency = BF16 16×16 AMX tile op — measured"); + println!(" tile M=16 (the AMX tile; 64×64 gridlake = 4×4 = 16 of these), K padded 16→32"); + println!( + " shipped kernel tier on THIS host: {} (amx_available = {})\n", + bf16_tile_gemm_tier(), + amx_available() + ); + + // ── field tile X: a correlated gridlake patch (smooth → strong interdependency) + let mut x = [0.0f32; M * M]; + for r in 0..M { + for c in 0..M { + x[r * M + c] = 4.0 * (r as f32 * 0.4).sin() * (c as f32 * 0.35).cos() + 0.5 * (r as f32 - c as f32); + } + } + + // ── (1) intercorrelation C = Xᵀ·X — the inter-cell coupling as a field op ── + let xt = transpose(&x); + let corr_tile = tile_matmul(&xt, &x); + let corr_direct = direct_matmul(&xt, &x); + let (ca, cr) = errors(&corr_tile, &corr_direct); + println!(" (1) intercorrelation C = Xᵀ·X (one tile op)"); + println!(" max_abs_err={ca:.4} frobenius_rel_err={:.3}% (BF16-precision class)", cr * 100.0); + + // ── (2) covariance sandwich Σ' = Mᵀ·Σ·M — the splat EWA projection shape ── + // Σ = a symmetric PSD covariance (Xᵀ·X); M = a coupling/projection matrix. + let sigma_f: [f32; M * M] = { + let mut s = [0.0f32; M * M]; + for i in 0..M * M { + s[i] = corr_direct[i] as f32; + } + s + }; + let mut m_proj = [0.0f32; M * M]; + for r in 0..M { + for c in 0..M { + m_proj[r * M + c] = if r == c { + 1.0 + } else { + 0.15 * (r as f32 - c as f32).signum() + }; + } + } + // tile path: T = Σ·M, then Σ' = Mᵀ·T (two tile ops — the sandwich) + let t_tile = tile_matmul(&sigma_f, &m_proj); + let mt = transpose(&m_proj); + let sand_tile = tile_matmul(&mt, &t_tile); + // direct reference: Mᵀ·Σ·M in f64 + let t_direct = { + let mut td = [0.0f32; M * M]; + let d = direct_matmul(&sigma_f, &m_proj); + for i in 0..M * M { + td[i] = d[i] as f32; + } + td + }; + let sand_direct = direct_matmul(&mt, &t_direct); + let (sa, sr) = errors(&sand_tile, &sand_direct); + println!(" (2) covariance sandwich Σ' = Mᵀ·Σ·M (two tile ops = the splat EWA projection)"); + println!( + " max_abs_err={sa:.4} frobenius_rel_err={:.3}% (== hpc::splat3d::sandwich_x16 shape;", + sr * 100.0 + ); + println!(" per-entry max-rel is meaningless here — Mᵀ·Σ·M cancels to near-zero entries)"); + + // ── (3) bit-exactness on bf16-exact integer operands (module guarantee) ── + let mut ai = [0.0f32; M * M]; + let mut bi = [0.0f32; M * M]; + for i in 0..M * M { + ai[i] = ((mix(i as u64) % 16) as f32) - 8.0; // small ints, bf16-exact + bi[i] = ((mix(i as u64 ^ 0x5555) % 16) as f32) - 8.0; + } + let int_tile = tile_matmul(&ai, &bi); + let int_direct = direct_matmul(&ai, &bi); + let bit_exact = (0..M * M).all(|i| int_tile[i] as f64 == int_direct[i]); + println!(" (3) bit-exact on bf16-exact integer operands (|Σ| < 2^24): {bit_exact}"); + + // ── (4) throughput of the 16×16 tile op on this host ── + let iters = 200_000usize; + let t0 = std::time::Instant::now(); + let mut acc = 0.0f32; + for it in 0..iters { + // vary the operand slightly so nothing is optimized away + let mut xv = x; + xv[0] += (it & 0x7) as f32 * 0.01; + let c = tile_matmul(&xv, &m_proj); + acc += c[0] + c[M * M - 1]; + } + let dt = t0.elapsed().as_secs_f64(); + let tiles_s = iters as f64 / dt; + println!( + " (4) throughput: {iters} tile ops in {:.1} ms → {:.2} M tile-ops/s (checksum {acc:.1})", + dt * 1000.0, + tiles_s / 1e6 + ); + + // hard gate: BF16-precision class (rel err small) + integer bit-exactness + assert!(cr < 0.05, "intercorrelation rel err too high for BF16 class"); + assert!(sr < 0.05, "sandwich rel err too high for BF16 class"); + assert!(bit_exact, "bf16-exact integer tile op is NOT bit-exact vs direct"); + + println!( + "\n MEASURED CONCLUSION: YES — a gridlake field interdependency reduces to a\n\ + \x20 BF16 16×16 tile op (one `bf16_tile_gemm_16x16`, i.e. one `TDPBF16PS` on\n\ + \x20 the AMX tier), matching the direct f64 reference to BF16 precision and\n\ + \x20 BIT-EXACT for bf16-exact integer operands. It IS the same primitive as:\n\ + \x20 • the codec separable transform (#232: M·X on a tile — WHT/DCT), and\n\ + \x20 • the Gaussian-splat EWA covariance sandwich Σ'=Mᵀ·Σ·M\n\ + \x20 (hpc::splat3d::spd3::sandwich_x16 — the 16-wide batched form).\n\ + \x20 Interdependency / intercorrelation / perturbation / transform / covariance\n\ + \x20 projection are ONE op: the 16×16 BF16 tile GEMM. The 64×64 gridlake is\n\ + \x20 4×4 = 16 of these tiles." + ); +} From e44003ec87664ea8953e9542ed95fd9955f48d8f Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 4 Jul 2026 10:44:04 +0000 Subject: [PATCH 2/3] review: gate gridlake_field_tile to x86_64 + use canonical simd surface (#233, codex P2) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The example imported the x86_64-only hpc::amx_matmul / hpc::bf16_tile_gemm paths under required-features = ["std"], so a cross-target `cargo check --examples` would fail on non-x86 with unresolved imports. Fixes: - Wrap the whole probe in `#[cfg(target_arch = "x86_64")] mod amx { ... }` with a `#[cfg(not(target_arch = "x86_64"))]` fallback `main` that prints the x86-only requirement. Non-x86 builds now compile (and skip) cleanly. - Switch to the canonical `ndarray::simd::*` surface (per the W1a consumer contract — consumers import from `ndarray::simd`, not raw `hpc::*`): `amx_available`, `bf16_tile_gemm_16x16_amx` (the TDPBF16PS tile wrapper, aliased to the body's name), `bf16_tile_gemm_tier` — all arch-gated re-exports. Behaviour unchanged on x86: still runs the real AMX TDPBF16PS tile op (amx_available = true), 0.095% / 0.401% Frobenius, bit-exact on integer operands. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01MLBnPuScZy6w9di2QEjsXM --- examples/gridlake_field_tile.rs | 360 +++++++++++++++++--------------- 1 file changed, 190 insertions(+), 170 deletions(-) diff --git a/examples/gridlake_field_tile.rs b/examples/gridlake_field_tile.rs index 5e948bea..b980dc45 100644 --- a/examples/gridlake_field_tile.rs +++ b/examples/gridlake_field_tile.rs @@ -30,195 +30,204 @@ //! //! Run: `cargo run --release --example gridlake_field_tile --features std` -use ndarray::hpc::amx_matmul::amx_available; -use ndarray::hpc::bf16_tile_gemm::{bf16_tile_gemm_16x16, bf16_tile_gemm_tier}; -use ndarray::hpc::quantized::{f32_to_bf16_rounded, BF16}; - -const M: usize = 16; // tile edge — the AMX tile / gridlake quadrant -const KPAD: usize = 32; // bf16_tile_gemm requires K a multiple of 32; pad 16→32 with zeros - -fn mix(mut z: u64) -> u64 { - z = z.wrapping_add(0x9E37_79B9_7F4A_7C15); - z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); - z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); - z ^ (z >> 31) -} - -/// Pack a row-major `M×KPAD` (or `KPAD×M`) f32 matrix into the `&[u16]` bf16 -/// operand the tile GEMM consumes (RNE, matching `VCVTNEPS2BF16`). -fn to_bf16(src: &[f32]) -> Vec { - let mut b = vec![BF16(0); src.len()]; - f32_to_bf16_rounded(src, &mut b); - b.iter().map(|x| x.0).collect() -} +// The AMX / BF16-tile path is x86_64-only: `hpc::amx_matmul` and +// `hpc::bf16_tile_gemm` are `#[cfg(target_arch = "x86_64")]`, and the canonical +// `ndarray::simd::*` re-exports of the tile ladder are arch-gated to match. Wrap +// the whole probe in an x86 module with a non-x86 fallback `main`, so a +// cross-target `cargo check --examples` does not hit unresolved imports. +#[cfg(target_arch = "x86_64")] +mod amx { + // Canonical public surface (`ndarray::simd::*`, per the W1a consumer contract), + // not the raw `hpc::*` paths. `bf16_tile_gemm_16x16_amx` is the tile-dispatching + // TDPBF16PS wrapper (aliased to the name the body already uses). + use ndarray::hpc::quantized::{f32_to_bf16_rounded, BF16}; + use ndarray::simd::{amx_available, bf16_tile_gemm_16x16_amx as bf16_tile_gemm_16x16, bf16_tile_gemm_tier}; + + const M: usize = 16; // tile edge — the AMX tile / gridlake quadrant + const KPAD: usize = 32; // bf16_tile_gemm requires K a multiple of 32; pad 16→32 with zeros + + fn mix(mut z: u64) -> u64 { + z = z.wrapping_add(0x9E37_79B9_7F4A_7C15); + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } -/// `C[16×16] = A[16×16] · B[16×16]` via the SHIPPED BF16 tile op. A/B are padded -/// on K from 16→32 with zeros (the pad contributes nothing), so this is one -/// `bf16_tile_gemm_16x16` — one `TDPBF16PS` tile op on the AMX tier. -fn tile_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f32; M * M] { - let mut a_pad = [0.0f32; M * KPAD]; - let mut b_pad = [0.0f32; KPAD * M]; - for i in 0..M { - for k in 0..M { - a_pad[i * KPAD + k] = a[i * M + k]; // A[16×32], cols 16..31 = 0 - b_pad[k * M + i] = b[k * M + i]; // B[32×16], rows 16..31 = 0 (already) - } + /// Pack a row-major `M×KPAD` (or `KPAD×M`) f32 matrix into the `&[u16]` bf16 + /// operand the tile GEMM consumes (RNE, matching `VCVTNEPS2BF16`). + fn to_bf16(src: &[f32]) -> Vec { + let mut b = vec![BF16(0); src.len()]; + f32_to_bf16_rounded(src, &mut b); + b.iter().map(|x| x.0).collect() } - let (ab, bb) = (to_bf16(&a_pad), to_bf16(&b_pad)); - let mut c = vec![0.0f32; M * M]; - bf16_tile_gemm_16x16(&ab, &bb, &mut c, KPAD); - let mut out = [0.0f32; M * M]; - out.copy_from_slice(&c); - out -} -/// Direct `f64` reference GEMM `C = A·B` (16×16). -fn direct_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f64; M * M] { - let mut c = [0.0f64; M * M]; - for i in 0..M { - for j in 0..M { - let mut s = 0.0f64; + /// `C[16×16] = A[16×16] · B[16×16]` via the SHIPPED BF16 tile op. A/B are padded + /// on K from 16→32 with zeros (the pad contributes nothing), so this is one + /// `bf16_tile_gemm_16x16` — one `TDPBF16PS` tile op on the AMX tier. + fn tile_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f32; M * M] { + let mut a_pad = [0.0f32; M * KPAD]; + let mut b_pad = [0.0f32; KPAD * M]; + for i in 0..M { for k in 0..M { - s += a[i * M + k] as f64 * b[k * M + j] as f64; + a_pad[i * KPAD + k] = a[i * M + k]; // A[16×32], cols 16..31 = 0 + b_pad[k * M + i] = b[k * M + i]; // B[32×16], rows 16..31 = 0 (already) } - c[i * M + j] = s; } + let (ab, bb) = (to_bf16(&a_pad), to_bf16(&b_pad)); + let mut c = vec![0.0f32; M * M]; + bf16_tile_gemm_16x16(&ab, &bb, &mut c, KPAD); + let mut out = [0.0f32; M * M]; + out.copy_from_slice(&c); + out } - c -} -fn transpose(a: &[f32; M * M]) -> [f32; M * M] { - let mut t = [0.0f32; M * M]; - for i in 0..M { - for j in 0..M { - t[j * M + i] = a[i * M + j]; + /// Direct `f64` reference GEMM `C = A·B` (16×16). + fn direct_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f64; M * M] { + let mut c = [0.0f64; M * M]; + for i in 0..M { + for j in 0..M { + let mut s = 0.0f64; + for k in 0..M { + s += a[i * M + k] as f64 * b[k * M + j] as f64; + } + c[i * M + j] = s; + } } + c } - t -} -/// `(max_abs_err, frobenius_relative_err)` of a tile-GEMM result vs the f64 -/// direct. Frobenius-relative `‖tile − direct‖_F / ‖direct‖_F` is the honest -/// aggregate precision — a per-element max-relative blows up on the near-zero -/// entries that a cancellation-heavy product (e.g. `Mᵀ·Σ·M`) naturally has, -/// which measures the cancellation, not the kernel. -fn errors(tile: &[f32; M * M], direct: &[f64; M * M]) -> (f64, f64) { - let mut max_abs = 0.0f64; - let mut num = 0.0f64; // ‖E‖_F² - let mut den = 0.0f64; // ‖direct‖_F² - for i in 0..M * M { - let e = tile[i] as f64 - direct[i]; - max_abs = max_abs.max(e.abs()); - num += e * e; - den += direct[i] * direct[i]; + fn transpose(a: &[f32; M * M]) -> [f32; M * M] { + let mut t = [0.0f32; M * M]; + for i in 0..M { + for j in 0..M { + t[j * M + i] = a[i * M + j]; + } + } + t } - (max_abs, (num / den.max(1e-12)).sqrt()) -} -fn main() { - println!("Gridlake field interdependency = BF16 16×16 AMX tile op — measured"); - println!(" tile M=16 (the AMX tile; 64×64 gridlake = 4×4 = 16 of these), K padded 16→32"); - println!( - " shipped kernel tier on THIS host: {} (amx_available = {})\n", - bf16_tile_gemm_tier(), - amx_available() - ); - - // ── field tile X: a correlated gridlake patch (smooth → strong interdependency) - let mut x = [0.0f32; M * M]; - for r in 0..M { - for c in 0..M { - x[r * M + c] = 4.0 * (r as f32 * 0.4).sin() * (c as f32 * 0.35).cos() + 0.5 * (r as f32 - c as f32); + /// `(max_abs_err, frobenius_relative_err)` of a tile-GEMM result vs the f64 + /// direct. Frobenius-relative `‖tile − direct‖_F / ‖direct‖_F` is the honest + /// aggregate precision — a per-element max-relative blows up on the near-zero + /// entries that a cancellation-heavy product (e.g. `Mᵀ·Σ·M`) naturally has, + /// which measures the cancellation, not the kernel. + fn errors(tile: &[f32; M * M], direct: &[f64; M * M]) -> (f64, f64) { + let mut max_abs = 0.0f64; + let mut num = 0.0f64; // ‖E‖_F² + let mut den = 0.0f64; // ‖direct‖_F² + for i in 0..M * M { + let e = tile[i] as f64 - direct[i]; + max_abs = max_abs.max(e.abs()); + num += e * e; + den += direct[i] * direct[i]; } + (max_abs, (num / den.max(1e-12)).sqrt()) } - // ── (1) intercorrelation C = Xᵀ·X — the inter-cell coupling as a field op ── - let xt = transpose(&x); - let corr_tile = tile_matmul(&xt, &x); - let corr_direct = direct_matmul(&xt, &x); - let (ca, cr) = errors(&corr_tile, &corr_direct); - println!(" (1) intercorrelation C = Xᵀ·X (one tile op)"); - println!(" max_abs_err={ca:.4} frobenius_rel_err={:.3}% (BF16-precision class)", cr * 100.0); - - // ── (2) covariance sandwich Σ' = Mᵀ·Σ·M — the splat EWA projection shape ── - // Σ = a symmetric PSD covariance (Xᵀ·X); M = a coupling/projection matrix. - let sigma_f: [f32; M * M] = { - let mut s = [0.0f32; M * M]; - for i in 0..M * M { - s[i] = corr_direct[i] as f32; + pub fn run() { + println!("Gridlake field interdependency = BF16 16×16 AMX tile op — measured"); + println!(" tile M=16 (the AMX tile; 64×64 gridlake = 4×4 = 16 of these), K padded 16→32"); + println!( + " shipped kernel tier on THIS host: {} (amx_available = {})\n", + bf16_tile_gemm_tier(), + amx_available() + ); + + // ── field tile X: a correlated gridlake patch (smooth → strong interdependency) + let mut x = [0.0f32; M * M]; + for r in 0..M { + for c in 0..M { + x[r * M + c] = 4.0 * (r as f32 * 0.4).sin() * (c as f32 * 0.35).cos() + 0.5 * (r as f32 - c as f32); + } } - s - }; - let mut m_proj = [0.0f32; M * M]; - for r in 0..M { - for c in 0..M { - m_proj[r * M + c] = if r == c { - 1.0 - } else { - 0.15 * (r as f32 - c as f32).signum() - }; + + // ── (1) intercorrelation C = Xᵀ·X — the inter-cell coupling as a field op ── + let xt = transpose(&x); + let corr_tile = tile_matmul(&xt, &x); + let corr_direct = direct_matmul(&xt, &x); + let (ca, cr) = errors(&corr_tile, &corr_direct); + println!(" (1) intercorrelation C = Xᵀ·X (one tile op)"); + println!(" max_abs_err={ca:.4} frobenius_rel_err={:.3}% (BF16-precision class)", cr * 100.0); + + // ── (2) covariance sandwich Σ' = Mᵀ·Σ·M — the splat EWA projection shape ── + // Σ = a symmetric PSD covariance (Xᵀ·X); M = a coupling/projection matrix. + let sigma_f: [f32; M * M] = { + let mut s = [0.0f32; M * M]; + for i in 0..M * M { + s[i] = corr_direct[i] as f32; + } + s + }; + let mut m_proj = [0.0f32; M * M]; + for r in 0..M { + for c in 0..M { + m_proj[r * M + c] = if r == c { + 1.0 + } else { + 0.15 * (r as f32 - c as f32).signum() + }; + } } - } - // tile path: T = Σ·M, then Σ' = Mᵀ·T (two tile ops — the sandwich) - let t_tile = tile_matmul(&sigma_f, &m_proj); - let mt = transpose(&m_proj); - let sand_tile = tile_matmul(&mt, &t_tile); - // direct reference: Mᵀ·Σ·M in f64 - let t_direct = { - let mut td = [0.0f32; M * M]; - let d = direct_matmul(&sigma_f, &m_proj); + // tile path: T = Σ·M, then Σ' = Mᵀ·T (two tile ops — the sandwich) + let t_tile = tile_matmul(&sigma_f, &m_proj); + let mt = transpose(&m_proj); + let sand_tile = tile_matmul(&mt, &t_tile); + // direct reference: Mᵀ·Σ·M in f64 + let t_direct = { + let mut td = [0.0f32; M * M]; + let d = direct_matmul(&sigma_f, &m_proj); + for i in 0..M * M { + td[i] = d[i] as f32; + } + td + }; + let sand_direct = direct_matmul(&mt, &t_direct); + let (sa, sr) = errors(&sand_tile, &sand_direct); + println!(" (2) covariance sandwich Σ' = Mᵀ·Σ·M (two tile ops = the splat EWA projection)"); + println!( + " max_abs_err={sa:.4} frobenius_rel_err={:.3}% (== hpc::splat3d::sandwich_x16 shape;", + sr * 100.0 + ); + println!(" per-entry max-rel is meaningless here — Mᵀ·Σ·M cancels to near-zero entries)"); + + // ── (3) bit-exactness on bf16-exact integer operands (module guarantee) ── + let mut ai = [0.0f32; M * M]; + let mut bi = [0.0f32; M * M]; for i in 0..M * M { - td[i] = d[i] as f32; + ai[i] = ((mix(i as u64) % 16) as f32) - 8.0; // small ints, bf16-exact + bi[i] = ((mix(i as u64 ^ 0x5555) % 16) as f32) - 8.0; } - td - }; - let sand_direct = direct_matmul(&mt, &t_direct); - let (sa, sr) = errors(&sand_tile, &sand_direct); - println!(" (2) covariance sandwich Σ' = Mᵀ·Σ·M (two tile ops = the splat EWA projection)"); - println!( - " max_abs_err={sa:.4} frobenius_rel_err={:.3}% (== hpc::splat3d::sandwich_x16 shape;", - sr * 100.0 - ); - println!(" per-entry max-rel is meaningless here — Mᵀ·Σ·M cancels to near-zero entries)"); - - // ── (3) bit-exactness on bf16-exact integer operands (module guarantee) ── - let mut ai = [0.0f32; M * M]; - let mut bi = [0.0f32; M * M]; - for i in 0..M * M { - ai[i] = ((mix(i as u64) % 16) as f32) - 8.0; // small ints, bf16-exact - bi[i] = ((mix(i as u64 ^ 0x5555) % 16) as f32) - 8.0; - } - let int_tile = tile_matmul(&ai, &bi); - let int_direct = direct_matmul(&ai, &bi); - let bit_exact = (0..M * M).all(|i| int_tile[i] as f64 == int_direct[i]); - println!(" (3) bit-exact on bf16-exact integer operands (|Σ| < 2^24): {bit_exact}"); - - // ── (4) throughput of the 16×16 tile op on this host ── - let iters = 200_000usize; - let t0 = std::time::Instant::now(); - let mut acc = 0.0f32; - for it in 0..iters { - // vary the operand slightly so nothing is optimized away - let mut xv = x; - xv[0] += (it & 0x7) as f32 * 0.01; - let c = tile_matmul(&xv, &m_proj); - acc += c[0] + c[M * M - 1]; - } - let dt = t0.elapsed().as_secs_f64(); - let tiles_s = iters as f64 / dt; - println!( - " (4) throughput: {iters} tile ops in {:.1} ms → {:.2} M tile-ops/s (checksum {acc:.1})", - dt * 1000.0, - tiles_s / 1e6 - ); - - // hard gate: BF16-precision class (rel err small) + integer bit-exactness - assert!(cr < 0.05, "intercorrelation rel err too high for BF16 class"); - assert!(sr < 0.05, "sandwich rel err too high for BF16 class"); - assert!(bit_exact, "bf16-exact integer tile op is NOT bit-exact vs direct"); - - println!( - "\n MEASURED CONCLUSION: YES — a gridlake field interdependency reduces to a\n\ + let int_tile = tile_matmul(&ai, &bi); + let int_direct = direct_matmul(&ai, &bi); + let bit_exact = (0..M * M).all(|i| int_tile[i] as f64 == int_direct[i]); + println!(" (3) bit-exact on bf16-exact integer operands (|Σ| < 2^24): {bit_exact}"); + + // ── (4) throughput of the 16×16 tile op on this host ── + let iters = 200_000usize; + let t0 = std::time::Instant::now(); + let mut acc = 0.0f32; + for it in 0..iters { + // vary the operand slightly so nothing is optimized away + let mut xv = x; + xv[0] += (it & 0x7) as f32 * 0.01; + let c = tile_matmul(&xv, &m_proj); + acc += c[0] + c[M * M - 1]; + } + let dt = t0.elapsed().as_secs_f64(); + let tiles_s = iters as f64 / dt; + println!( + " (4) throughput: {iters} tile ops in {:.1} ms → {:.2} M tile-ops/s (checksum {acc:.1})", + dt * 1000.0, + tiles_s / 1e6 + ); + + // hard gate: BF16-precision class (rel err small) + integer bit-exactness + assert!(cr < 0.05, "intercorrelation rel err too high for BF16 class"); + assert!(sr < 0.05, "sandwich rel err too high for BF16 class"); + assert!(bit_exact, "bf16-exact integer tile op is NOT bit-exact vs direct"); + + println!( + "\n MEASURED CONCLUSION: YES — a gridlake field interdependency reduces to a\n\ \x20 BF16 16×16 tile op (one `bf16_tile_gemm_16x16`, i.e. one `TDPBF16PS` on\n\ \x20 the AMX tier), matching the direct f64 reference to BF16 precision and\n\ \x20 BIT-EXACT for bf16-exact integer operands. It IS the same primitive as:\n\ @@ -228,5 +237,16 @@ fn main() { \x20 Interdependency / intercorrelation / perturbation / transform / covariance\n\ \x20 projection are ONE op: the 16×16 BF16 tile GEMM. The 64×64 gridlake is\n\ \x20 4×4 = 16 of these tiles." - ); + ); + } +} // mod amx + +#[cfg(target_arch = "x86_64")] +fn main() { + amx::run(); +} + +#[cfg(not(target_arch = "x86_64"))] +fn main() { + eprintln!("gridlake_field_tile requires x86_64 (AMX TDPBF16PS / AVX-512 VDPBF16PS BF16 tile GEMM)."); } From f10c2c4e425c72c0aab8017db9903bfb434ebee6 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 4 Jul 2026 10:58:09 +0000 Subject: [PATCH 3/3] git: union-merge Cargo.toml to end the example-list conflict cascade MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every topic branch that adds a [[example]] entry collides on the same region of Cargo.toml, so each merge dirties the next chained PR. The built-in `union` merge driver keeps both sides of an additive conflict hunk — exactly right for this append-only list — and GitHub's merge honors it. Once this lands on master, subsequent example-adding PRs merge with no manual conflict resolution. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01MLBnPuScZy6w9di2QEjsXM --- .gitattributes | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 .gitattributes diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 00000000..4313aef9 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,5 @@ +# Additive edits to the [[example]] list collide constantly across topic +# branches. The built-in union merge driver keeps both sides of a conflict +# hunk, auto-resolving append-only lists. (Distinct dependency-line edits are +# caught by CI, not silently merged wrong.) +Cargo.toml merge=union