From 58bfbfe5e5c6218b29ca4b3d1fb2a39e93104d21 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 12 Jul 2026 06:09:24 +0000 Subject: [PATCH 01/15] =?UTF-8?q?motivic:=20A=5FC,=20the=20C-motivic=20Ste?= =?UTF-8?q?enrod=20algebra=20engine=20over=20F=E2=82=82[=CF=84]?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Foundation layer for computing the C-motivic Adams E₂ by deformation: the coefficient ring and the product engine, with no wiring into the resolution engine yet — the mod-τ reduction that the engine resolves comes in a follow-up. - `tau`: F₂[τ] as a small homogeneous scalar (`Tau`). Every structure constant in the motivic world is a single power of τ, so a coefficient is a one-integer valuation rather than a polynomial. - `milnor`: `MotivicMilnorAlgebra` = A_C as a free F₂[τ]-module on the Milnor basis Q(E)P(R). The product is computed two ways — a duality oracle (dualize the coproduct ψ) and the closed-form Kong–Lin Theorem 5.1 (arXiv:2411.12890, ρ = 0) — and the fast path is validated exhaustively against the oracle. It is intentionally not an `Algebra`: that trait is over F_p, and this is the F₂[τ] engine the deformation lift builds on. Built against the bit-packed classical Milnor basis (#280), so the two classical cross-check tests share one `classical_mul` helper plus the paper/`PPart` index conversions rather than each carrying its own, and `Monomial` is a named type rather than `(u32, Vec)` spelled out 19 times. Tests: rewrite_tau identities, product associativity, weight-homogeneity, and closed-form-vs-duality agreement over a range. Co-Authored-By: Claude Opus 4.8 Co-Authored-By: Claude Opus 5 --- ext/crates/algebra/src/algebra/mod.rs | 3 + .../algebra/src/algebra/motivic/milnor.rs | 1737 +++++++++++++++++ ext/crates/algebra/src/algebra/motivic/mod.rs | 28 + ext/crates/algebra/src/algebra/motivic/tau.rs | 142 ++ 4 files changed, 1910 insertions(+) create mode 100644 ext/crates/algebra/src/algebra/motivic/milnor.rs create mode 100644 ext/crates/algebra/src/algebra/motivic/mod.rs create mode 100644 ext/crates/algebra/src/algebra/motivic/tau.rs diff --git a/ext/crates/algebra/src/algebra/mod.rs b/ext/crates/algebra/src/algebra/mod.rs index 67baed2b04..44c0f4454a 100644 --- a/ext/crates/algebra/src/algebra/mod.rs +++ b/ext/crates/algebra/src/algebra/mod.rs @@ -18,6 +18,9 @@ pub use field::Field; pub mod milnor_algebra; pub use milnor_algebra::MilnorAlgebra; +pub mod motivic; +pub use motivic::{MotivicMilnorAlgebra, Tau}; + mod steenrod_algebra; pub use steenrod_algebra::{AlgebraType, SteenrodAlgebra}; diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs new file mode 100644 index 0000000000..c901ff5b17 --- /dev/null +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -0,0 +1,1737 @@ +//! The C-motivic (over $\mathbb{C}$, prime 2) dual Steenrod algebra $A_{**}$ and +//! the Steenrod algebra $A_C$ in the Milnor basis, over $\mathbb{F}_2[\tau]$. +//! +//! This follows Kong–Lin, *Product formulas for motivic Milnor basis* +//! (arXiv:2411.12890), specialized to the C-motivic point where $\rho = 0$ (the +//! motivic cohomology of a point over $\mathbb{C}$ is $\mathbb{F}_2[\tau]$). We +//! use the *conjugate* generators of that paper, so the coproduct and +//! Milnor-matrix product match Milnor's classical formulas on the polynomial +//! (`ξ`) part. +//! +//! The dual algebra is $A_{**} = \mathbb{F}_2[\tau][\xi_1, \xi_2, \dots] \otimes +//! E(\tau_0, \tau_1, \dots)$ with $\tau_i^2 = \tau \xi_{i+1}$ (the $\rho = 0$ +//! reduction of $\tau_i^2 = \tau\xi_{i+1} + \rho\tau_{i+1}$). Bidegrees: +//! $|\xi_i| = (2(2^i-1), 2^i-1)$, $|\tau_i| = (2^{i+1}-1, 2^i-1)$, +//! $|\tau| = (0, -1)$ in (stem, weight). +//! +//! Basis elements are $Q(E)P(R)$ with $E \in \mathrm{Seq}_1$ (entries in +//! $\{0,1\}$, encoded as a bitmask: bit $i$ is $\tau_i$) and $R \in \mathrm{Seq}$ +//! (the exponent vector of the $\xi_i$). +//! +//! The product in the Steenrod algebra $A_C$ is computed by **duality**: the +//! dual algebra $A_{**}$ is a commutative $\mathbb{F}_2[\tau]$-algebra whose +//! multiplication ([`dual_mul`]) reduces $\tau_i^2 = \tau\xi_{i+1}$ via +//! [`rewrite_tau`] (Kong–Lin Theorem 3.4), and its coproduct $\psi$ +//! ([`coproduct`], Kong–Lin §2.2) is a $\tau$-free algebra map. The product +//! $a \cdot b$ in $A_C$ is then read off from $\psi$: the coefficient of $z$ in +//! $a \cdot b$ is the coefficient of $\mathrm{mon}(a) \otimes \mathrm{mon}(b)$ in +//! $\psi(\mathrm{mon}(z))$ ([`multiply`]). This is a first, correctness-oriented +//! implementation; the closed-form product (Kong–Lin Theorem 5.1) can replace it +//! later for speed, validated against this one. +//! +//! Coefficients are homogeneous, so each one is a single power of $\tau$: a +//! [`Tau`] scalar (see [`super::tau`]). This is the $A_C$ product **engine**, +//! over $\mathbb{F}_2[\tau]$; the mod-$\tau$ reduction $A_C/\tau$ is presented to +//! the resolution engine as an ordinary $\mathbb{F}_2$-algebra in a follow-up +//! layer built on top of this engine. +//! +//! Weight convention: the motivic weight of an algebra basis element is the +//! *negative* of the weight of the dual monomial it pairs with, so that products +//! are weight-homogeneous with $\tau$ (weight $-1$) absorbing the difference. + +use std::{ + cell::RefCell, + collections::BTreeMap, + rc::Rc, + sync::{Arc, OnceLock, RwLock}, +}; + +use fp::prime::Binomial; +use itertools::Itertools; +use maybe_rayon::prelude::*; +use once::OnceVec; +use rustc_hash::FxHashMap; + +use super::Tau; + +/// Drop the trailing zeros from an exponent vector. +/// +/// Exponent vectors are compared and used as map keys throughout, so they are kept trimmed to +/// make that comparison canonical: `ξ(R)` does not see the slots past the last non-zero one. +fn trim(mut v: Vec) -> Vec { + while let Some(&0) = v.last() { + v.pop(); + } + v +} + +/// A basis monomial $\tau(E)\xi(R)$ of the dual algebra: the square-free exterior part `E` as a +/// bitmask (bit `i` is $\tau_i$), and the $\xi$ exponent vector `R` in the paper's indexing, +/// where `R[0]` belongs to $\xi_0 = 1$ and is therefore ignored. Trailing zeros are trimmed, +/// so equal monomials have equal representations and the tuple can be used as a map key. +pub type Monomial = (u32, Vec); + +/// $\Sigma(R) = \sum_i r_i$. +fn sigma(r: &[u32]) -> u32 { + r.iter().sum() +} + +/// $\Sigma_2(R) = \sum_i r_i 2^i$. +fn sigma2(r: &[u32]) -> u64 { + r.iter().enumerate().map(|(i, &v)| (v as u64) << i).sum() +} + +/// The coefficient $c(S, R)$ of Notation 3.3 / Theorem 3.4, reduced mod 2. +/// +/// $c(S, R) = \prod_{n \ge 1} \binom{\lfloor \sum_{i=0}^{n-1} 2^{i-n}(s_i - r_i)\rfloor}{r_n}$ +/// when $r_0 = 0$, and $0$ otherwise. The inner floor is computed exactly with integer +/// (Euclidean) division: $\lfloor \sum_{i=0}^{n-1} 2^{i-n}(s_i-r_i)\rfloor +/// = \lfloor (\sum_{i=0}^{n-1} (s_i-r_i)2^i) / 2^n \rfloor$. Each factor is +/// [`Binomial::binomial2`], which is $0$ when the top is negative or below the bottom. +fn c_coeff(s: &[u32], r: &[u32]) -> u32 { + if r.first().copied().unwrap_or(0) > 0 { + return 0; + } + let get = |seq: &[u32], i: usize| seq.get(i).copied().unwrap_or(0) as i32; + let len = s.len().max(r.len()); + let mut prod = 1u32; + for n in 1..len { + let mut num: i32 = 0; + for i in 0..n { + num += (get(s, i) - get(r, i)) << i; + } + let floor = num.div_euclid(1i32 << n); + prod *= i32::binomial2(floor, get(r, n)) as u32; + if prod == 0 { + return 0; + } + } + prod +} + +/// A term $\tau^{\text{tau\_pow}} \cdot Q(E) P(R)$ (mod 2, C-motivic so no $\rho$), the output +/// of rewriting a monomial into the Milnor basis. `e_mask` encodes $E \in \mathrm{Seq}_1$. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct MotivicTerm { + pub tau_pow: u32, + pub e_mask: u32, + pub r: Vec, +} + +/// Rewrite $\tau(S) = \tau_0^{s_0}\tau_1^{s_1}\cdots$ (with `s` an arbitrary exponent vector) +/// into the C-motivic Milnor basis, using $\tau_i^2 = \tau\xi_{i+1}$ (Theorem 3.4 at +/// $\rho = 0$): +/// $$\tau(S) = \sum_{\substack{E \in \mathrm{Seq}_1,\ R \in \mathrm{Seq}\\ \Sigma_2(E)+\Sigma_2(R)=\Sigma_2(S)\\ \Sigma(E)+2\Sigma(R)=\Sigma(S)}} c(S,R)\, \tau^{\Sigma(R)}\, \tau(E)\xi(R),$$ +/// where the second constraint is the vanishing of the $\rho$-exponent +/// $\Sigma(S) - \Sigma(E) - 2\Sigma(R)$. +pub fn rewrite_tau(s: &[u32]) -> Vec { + let target2 = sigma2(s); + let target_sum = sigma(s); + // Highest index we must consider: any nonzero e_i or r_i contributes at least 2^i, so + // 2^i <= target2. + let max_idx = if target2 == 0 { + 0 + } else { + (63 - target2.leading_zeros()) as usize + }; + + let mut out = Vec::new(); + let mut e = vec![0u32; max_idx + 1]; + let mut r = vec![0u32; max_idx + 1]; + rewrite_tau_dfs(s, 0, max_idx, target2, target_sum, &mut e, &mut r, &mut out); + out +} + +#[allow(clippy::too_many_arguments)] +fn rewrite_tau_dfs( + s: &[u32], + idx: usize, + max_idx: usize, + rem2: u64, // remaining Σ2(E) + Σ2(R) budget + target_sum: u32, // Σ(S), for the ρ = 0 constraint Σ(E) + 2Σ(R) = Σ(S) + e: &mut [u32], + r: &mut [u32], + out: &mut Vec, +) { + if idx > max_idx { + if rem2 != 0 { + return; + } + // ρ-exponent must vanish (C-motivic): Σ(E) + 2Σ(R) = Σ(S). + if sigma(e) + 2 * sigma(r) != target_sum { + return; + } + let coeff = c_coeff(s, r); + if coeff == 0 { + return; + } + let mut e_mask = 0u32; + for (i, &ei) in e.iter().enumerate() { + if ei != 0 { + e_mask |= 1 << i; + } + } + out.push(MotivicTerm { + tau_pow: sigma(r), + e_mask, + r: trim(r.to_vec()), + }); + return; + } + + let weight = 1u64 << idx; + // e_idx ∈ {0, 1} (E ∈ Seq_1), r_idx ≥ 0, with (e_idx + r_idx) * 2^idx ≤ rem2. + for e_idx in 0..=1u32 { + if (e_idx as u64) * weight > rem2 { + break; + } + let after_e = rem2 - (e_idx as u64) * weight; + let max_r = (after_e / weight) as u32; + for r_idx in 0..=max_r { + e[idx] = e_idx; + r[idx] = r_idx; + rewrite_tau_dfs( + s, + idx + 1, + max_idx, + after_e - (r_idx as u64) * weight, + target_sum, + e, + r, + out, + ); + } + } + e[idx] = 0; + r[idx] = 0; +} + +// --------------------------------------------------------------------------- +// The dual algebra A_** (used to compute the product in A_C by duality). +// --------------------------------------------------------------------------- + +/// An element of the dual algebra $A_{**}$: [`Monomial`]s mapped to their +/// $\mathbb{F}_2[\tau]$ coefficients. +/// +/// Everything we compute is homogeneous — $\psi$ and the products are graded — so every +/// coefficient is a single power of $\tau$: a [`Tau`] scalar. Zero coefficients are never stored, +/// so equality is canonical, and coefficient arithmetic is exactly [`Tau`]'s arithmetic — `mul` +/// adds valuations, `add` cancels equal powers mod 2 (unequal powers would be inhomogeneous and +/// cannot arise). +pub type DualElement = BTreeMap; + +/// Elementwise sum of two exponent vectors, trimmed of trailing zeros. +fn vec_add(a: &[u32], b: &[u32]) -> Vec { + let n = a.len().max(b.len()); + trim( + (0..n) + .map(|i| a.get(i).copied().unwrap_or(0) + b.get(i).copied().unwrap_or(0)) + .collect(), + ) +} + +/// Add `coeff * key` into a sparse $\mathbb{F}_2[\tau]$-linear combination `acc`, dropping the +/// entry if it cancels to zero. Used for both [`DualElement`] and [`TensorElement`]; the +/// coefficient bookkeeping is entirely [`Tau`] arithmetic. +fn add_term(acc: &mut BTreeMap, key: K, coeff: Tau) { + use std::collections::btree_map::Entry; + if coeff.is_zero() { + return; + } + match acc.entry(key) { + Entry::Occupied(mut o) => { + let sum = *o.get() + coeff; + if sum.is_zero() { + o.remove(); + } else { + *o.get_mut() = sum; + } + } + Entry::Vacant(v) => { + v.insert(coeff); + } + } +} + +/// Multiply two basis monomials of $A_{**}$ and accumulate `coeff * (m1 * m2)` into `acc`. +/// +/// $A_{**}$ is commutative, so $\tau(E_1)\xi(R_1)\cdot\tau(E_2)\xi(R_2)$ is obtained by adding +/// exponents: the exterior parts form $S = E_1 + E_2$ (entries in $\{0,1,2\}$), which is +/// rewritten into the square-free basis via [`rewrite_tau`] ($\tau_i^2 = \tau\xi_{i+1}$), and +/// the resulting $\xi$ exponents are added to $R_1 + R_2$. +fn mul_monomials(m1: &Monomial, m2: &Monomial, coeff: Tau, acc: &mut DualElement) { + let (e1, e2) = (m1.0, m2.0); + let bits = u32::BITS - (e1 | e2).leading_zeros(); + let s: Vec = (0..bits) + .map(|i| ((e1 >> i) & 1) + ((e2 >> i) & 1)) + .collect(); + let r12 = vec_add(&m1.1, &m2.1); + for term in rewrite_tau(&s) { + let r = vec_add(&term.r, &r12); + add_term(acc, (term.e_mask, r), coeff * Tau::power(term.tau_pow)); + } +} + +/// The product of two elements of $A_{**}$. +pub fn dual_mul(a: &DualElement, b: &DualElement) -> DualElement { + let mut out = DualElement::new(); + for (m1, &c1) in a { + for (m2, &c2) in b { + mul_monomials(m1, m2, c1 * c2, &mut out); + } + } + out +} + +/// The generator $\tau_i \in A_{**}$. +pub fn tau_gen(i: usize) -> DualElement { + DualElement::from([((1u32 << i, vec![]), Tau::one())]) +} + +/// The generator $\xi_i \in A_{**}$ (for $i \ge 1$). +pub fn xi_gen(i: usize) -> DualElement { + let mut r = vec![0u32; i + 1]; + r[i] = 1; + DualElement::from([((0, r), Tau::one())]) +} + +/// The unit $1 \in A_{**}$. +pub fn dual_one() -> DualElement { + DualElement::from([((0, vec![]), Tau::one())]) +} + +// --------------------------------------------------------------------------- +// The coproduct ψ: A_** → A_** ⊗ A_**, an algebra map (Kong–Lin §2.2). +// --------------------------------------------------------------------------- + +/// An element of $A_{**} \otimes A_{**}$: pairs of [`Monomial`]s → $\mathbb{F}_2[\tau]$ +/// coefficients. Zero coefficients are never stored. +pub type TensorElement = BTreeMap<(Monomial, Monomial), Tau>; + +/// The unit $1 \otimes 1$. +fn tensor_one() -> TensorElement { + TensorElement::from([(((0, vec![]), (0, vec![])), Tau::one())]) +} + +/// Multiply in $A_{**} \otimes A_{**}$: $(a_L \otimes a_R)(b_L \otimes b_R) = (a_L b_L) \otimes (a_R b_R)$. +fn tensor_mul(t1: &TensorElement, t2: &TensorElement) -> TensorElement { + let mut out = TensorElement::new(); + for ((al, ar), &c1) in t1 { + for ((bl, br), &c2) in t2 { + let mut left = DualElement::new(); + mul_monomials(al, bl, Tau::one(), &mut left); + let mut right = DualElement::new(); + mul_monomials(ar, br, Tau::one(), &mut right); + let c = c1 * c2; + for (ml, &cl) in &left { + for (mr, &cr) in &right { + add_term(&mut out, (ml.clone(), mr.clone()), c * (cl * cr)); + } + } + } + } + out +} + +/// The monomial $\xi_j^p$ (with $\xi_0 = 1$). +fn xi_pow_mon(j: usize, p: u32) -> Monomial { + if j == 0 || p == 0 { + (0, vec![]) + } else { + let mut r = vec![0u32; j + 1]; + r[j] = p; + (0, r) + } +} + +/// $\psi(\tau_k) = 1 \otimes \tau_k + \sum_{i=0}^{k} \tau_i \otimes \xi_{k-i}^{2^i}$. +fn coprod_tau(k: usize) -> TensorElement { + let mut out = TensorElement::new(); + add_term(&mut out, ((0, vec![]), (1 << k, vec![])), Tau::one()); + for i in 0..=k { + add_term( + &mut out, + ((1 << i, vec![]), xi_pow_mon(k - i, 1 << i)), + Tau::one(), + ); + } + out +} + +/// $\psi(\xi_k) = \sum_{i=0}^{k} \xi_i \otimes \xi_{k-i}^{2^i}$ (with $\xi_0 = 1$). +fn coprod_xi(k: usize) -> TensorElement { + let mut out = TensorElement::new(); + for i in 0..=k { + add_term( + &mut out, + (xi_pow_mon(i, 1), xi_pow_mon(k - i, 1 << i)), + Tau::one(), + ); + } + out +} + +/// The coproduct of a basis monomial $\tau(E)\xi(R)$, computed as an algebra map: the product +/// of the coproducts of its generators. +fn coproduct_monomial(e_mask: u32, r: &[u32]) -> TensorElement { + let mut acc = tensor_one(); + for i in 0..u32::BITS { + if (e_mask >> i) & 1 != 0 { + acc = tensor_mul(&acc, &coprod_tau(i as usize)); + } + } + for (j, &rj) in r.iter().enumerate() { + for _ in 0..rj { + acc = tensor_mul(&acc, &coprod_xi(j)); + } + } + acc +} + +/// The coproduct of an arbitrary element of $A_{**}$ (extended $\mathbb{F}_2[\tau]$-linearly). +pub fn coproduct(elt: &DualElement) -> TensorElement { + let mut out = TensorElement::new(); + for ((e, r), &c) in elt { + for (key, cc) in coproduct_monomial(*e, r) { + add_term(&mut out, key, c * cc); + } + } + out +} + +// --------------------------------------------------------------------------- +// The antipode χ: A_** → A_** (conjugate generators ⟷ standard Milnor generators). +// --------------------------------------------------------------------------- + +/// $\xi_j^p \in A_{**}$ (with $\xi_0 = 1$). +fn xi_pow_elt(j: usize, p: u32) -> DualElement { + DualElement::from([(xi_pow_mon(j, p), Tau::one())]) +} + +/// $\chi(\xi_k)$, from the antipode axiom with this module's coproduct: +/// $\chi(\xi_k) = \sum_{i=0}^{k-1} \chi(\xi_i)\,\xi_{k-i}^{2^i}$, $\chi(\xi_0) = 1$. +fn chi_xi(k: usize) -> DualElement { + if k == 0 { + return dual_one(); + } + let mut acc = DualElement::new(); + for i in 0..k { + let term = dual_mul(&chi_xi(i), &xi_pow_elt(k - i, 1 << i)); + for (mon, c) in term { + add_term(&mut acc, mon, c); + } + } + acc +} + +/// $\chi(\tau_k)$, from the antipode axiom: +/// $\chi(\tau_k) = \tau_k + \sum_{i=0}^{k-1} \chi(\tau_i)\,\xi_{k-i}^{2^i}$. +fn chi_tau(k: usize) -> DualElement { + let mut acc = tau_gen(k); + for i in 0..k { + let term = dual_mul(&chi_tau(i), &xi_pow_elt(k - i, 1 << i)); + for (mon, c) in term { + add_term(&mut acc, mon, c); + } + } + acc +} + +/// The antipode $\chi$ of $A_{**}$ — an algebra map (since $A_{**}$ is commutative) extended +/// from `chi_tau`/`chi_xi` on the generators. It converts the paper's conjugate generators to +/// the standard Milnor/Voevodsky generators (`χ(ξ_2) = ξ_2 + ξ_1^3`, etc.) and is an +/// involution. +pub fn antipode(elt: &DualElement) -> DualElement { + let mut out = DualElement::new(); + for ((e, r), &c) in elt { + let mut m = dual_one(); + for i in 0..u32::BITS { + if (e >> i) & 1 != 0 { + m = dual_mul(&m, &chi_tau(i as usize)); + } + } + for (j, &rj) in r.iter().enumerate() { + for _ in 0..rj { + m = dual_mul(&m, &chi_xi(j)); + } + } + for (mon, cc) in m { + add_term(&mut out, mon, c * cc); + } + } + out +} + +// --------------------------------------------------------------------------- +// The product in the Steenrod algebra A_C, computed by dualizing ψ. +// --------------------------------------------------------------------------- + +/// The C-motivic (prime 2) bidegree `(t, w)` of the basis monomial `τ(E)ξ(R)`, in the +/// *paper's* R-convention (`r[j]` is the exponent of `ξ_j`, with `ξ_0 = 1` so `r[0]` is +/// ignored). `|τ_i| = (2^{i+1}-1, 2^i-1)`, `|ξ_j| = (2^{j+1}-2, 2^j-1)`. +fn paper_bidegree(e_mask: u32, r: &[u32]) -> (i32, i32) { + let (mut t, mut w) = (0i32, 0i32); + for i in 0..u32::BITS { + if (e_mask >> i) & 1 != 0 { + t += (1 << (i + 1)) - 1; + w += (1 << i) - 1; + } + } + for (j, &rj) in r.iter().enumerate().skip(1) { + t += rj as i32 * ((1 << (j + 1)) - 2); + w += rj as i32 * ((1 << j) - 1); + } + (t, w) +} + +/// All basis monomials `τ(E)ξ(R)` of a given topological degree `target`. +pub fn enum_basis(target: i32) -> Vec { + // Generators available up to `target`: τ_i of degree 2^{i+1}-1, ξ_j (j≥1) of degree + // 2^{j+1}-2. Represented as (is_tau, index, degree). + let mut gens: Vec<(bool, usize, i32)> = Vec::new(); + let mut i = 0; + while (1i32 << (i + 1)) - 1 <= target { + gens.push((true, i, (1 << (i + 1)) - 1)); + i += 1; + } + let mut j = 1; + while (1i32 << (j + 1)) - 2 <= target { + gens.push((false, j, (1 << (j + 1)) - 2)); + j += 1; + } + let r_len = gens + .iter() + .filter(|g| !g.0) + .map(|g| g.1 + 1) + .max() + .unwrap_or(0); + + let mut out = Vec::new(); + let mut e = 0u32; + let mut r = vec![0u32; r_len]; + enum_basis_dfs(&gens, 0, target, &mut e, &mut r, &mut out); + out +} + +fn enum_basis_dfs( + gens: &[(bool, usize, i32)], + idx: usize, + rem: i32, + e: &mut u32, + r: &mut [u32], + out: &mut Vec, +) { + if idx == gens.len() { + if rem == 0 { + out.push((*e, trim(r.to_vec()))); + } + return; + } + let (is_tau, index, deg) = gens[idx]; + let max_mult = if is_tau { 1 } else { rem / deg }; + for m in 0..=max_mult { + if m * deg > rem { + break; + } + if is_tau && m == 1 { + *e |= 1 << index; + } else if !is_tau { + r[index] = m as u32; + } + enum_basis_dfs(gens, idx + 1, rem - m * deg, e, r, out); + if is_tau { + *e &= !(1 << index); + } else { + r[index] = 0; + } + } +} + +/// The product `a · b` in the C-motivic Steenrod algebra `A_C`, where `a`, `b` are Milnor basis +/// elements `Q(E)P(R)` given as the monomials `(E, R)` they are dual to. The result is a map +/// from basis monomials to their `𝔽₂[τ]` coefficients. +/// +/// Computed by duality: the coefficient of `z` in `a · b` is the coefficient of +/// `mon(a) ⊗ mon(b)` in `ψ(mon(z))`, summed over the basis `z` of the appropriate +/// topological degree. +pub fn multiply(a: &Monomial, b: &Monomial) -> DualElement { + let key = (a.clone(), b.clone()); + let t = paper_bidegree(a.0, &a.1).0 + paper_bidegree(b.0, &b.1).0; + let mut out = DualElement::new(); + for z in enum_basis(t) { + if let Some(&c) = coproduct_monomial(z.0, &z.1).get(&key) + && !c.is_zero() + { + out.insert(z, c); + } + } + out +} + +// --------------------------------------------------------------------------- +// The closed-form product (Kong–Lin Theorem 5.1, C-motivic ρ = 0). +// --------------------------------------------------------------------------- +// +// Theorem 5.1 computes Q(E₁)P(R₁)·Q(E₂)P(R₂) as a sum over two Milnor-style +// matrices X, Y with the output monomial Q(E₂ + T(Y)) P(T(X)). The subtlety that +// makes the motivic case "much more complex than classical" (and that trips a +// naive reading) is the ξ₀ = 1 / τ₀ ≠ 1 asymmetry (Remark 4.3): **ξ-type** +// sequences (R₁, R₂, S(X), R(X), T(X)) are compared on indices ≥ 1 — index 0 is +// absorbed — while **τ-type** sequences (E, S(Y), T(Y)) keep index 0. +// +// The τ-rewriting piece of the formula — the coefficient c(S(Y), R₁−S(X)), the +// Σ₂ constraint, and the ρ = 0 filter — is exactly [`rewrite_tau`], which we +// reuse rather than re-derive. +// +// Mod τ this collapses towards the classical shape: A_C/τ ≅ 𝔽₂[ξᵢ] ⊗ E(τᵢ) has the +// odd-primary dual's structure (with 2ⁱ powers rather than pⁱ), and the X constraint +// tightens from S(X) ≤ R₁ to S(X) = R₁ — the classical admissible matrix. What remains +// over the classical product is the second matrix Y, carrying the τᵢ/Q-part interaction. + +/// All columns `(v₀, v₁, …)` of non-negative integers with `Σᵢ 2ⁱ vᵢ == target`. +fn columns_eq(target: u32) -> Vec> { + if target == 0 { + return vec![vec![]]; + } + let mut out = Vec::new(); + let mut v0 = target % 2; + while v0 <= target { + for rest in columns_eq((target - v0) / 2) { + let mut col = vec![v0]; + col.extend(rest); + out.push(trim(col)); + } + v0 += 2; + } + out +} + +/// All columns with `Σᵢ 2ⁱ vᵢ ≤ bound` (distinct weighted sums, so no duplicates). +fn columns_le(bound: u32) -> Vec> { + (0..=bound).flat_map(columns_eq).collect() +} + +thread_local! { + /// Per-thread memo of [`columns_eq`] / [`columns_le`], keyed by target/bound. + /// The candidate lists are pure functions of one integer and are requested + /// over and over across products; caching them avoids re-enumerating and + /// re-allocating. (`Rc` is confined to a single product's call stack.) + static COLS_EQ: RefCell>>>> = RefCell::new(FxHashMap::default()); + static COLS_LE: RefCell>>>> = RefCell::new(FxHashMap::default()); +} + +fn columns_eq_cached(target: u32) -> Rc>> { + COLS_EQ.with(|c| { + Rc::clone( + c.borrow_mut() + .entry(target) + .or_insert_with(|| Rc::new(columns_eq(target))), + ) + }) +} + +fn columns_le_cached(bound: u32) -> Rc>> { + COLS_LE.with(|c| { + Rc::clone( + c.borrow_mut() + .entry(bound) + .or_insert_with(|| Rc::new(columns_le(bound))), + ) + }) +} + +/// The column-0 options for the `X` matrix: `x_{0,0} = 0`, and `x_{i,0} ≤ R₁[i]` +/// for `i ≥ 1` (bounded because it feeds the row sum `S(X)_i ≤ R₁[i]`). +fn col0_x_options(r1: &[u32]) -> Vec> { + let mut result = vec![vec![0u32]]; + for &bound in r1.iter().skip(1) { + let mut next = Vec::new(); + for base in &result { + for v in 0..=bound { + let mut c = base.clone(); + c.push(v); + next.push(c); + } + } + result = next; + } + result.into_iter().map(trim).collect() +} + +/// Buffer size for row / column / anti-diagonal indices. All are bounded by +/// ~2·log₂(degree) + (number of ξ generators), far below this for any feasible +/// stem (a ξ or τ at index `i` already has degree ≥ 2^i). +const NB: usize = 64; + +/// Immutable context for the closed-form product recursion (Theorem 5.1, ρ = 0). +struct Closed<'a> { + r1v: &'a [u32], + r2v: &'a [u32], + l: usize, + e1_mask: u32, + e2_mask: u32, + sigma2_e1: i64, + /// Per-column candidate columns for the `X` matrix. + x_cands: &'a [Rc>>], +} + +/// Incrementally maintained matrix state: row sums `S`, a running XOR per +/// anti-diagonal (`or`, which equals the OR because we only keep disjoint entries) +/// for the `b`-coefficient pruning, and the anti-diagonal sums `T`. +struct Acc { + rows: [u32; NB], + or: [u32; NB], + sum: [u32; NB], +} + +impl Acc { + fn zero() -> Self { + Self { + rows: [0; NB], + or: [0; NB], + sum: [0; NB], + } + } +} + +impl<'a> Closed<'a> { + /// Enumerate the `X` matrix column by column. A column is rejected — pruning + /// the whole subtree — if a row sum would exceed `R₁` (rows ≥ 1; row 0 is + /// absorbed by ξ₀ = 1) or if any entry collides on its anti-diagonal (which + /// would make `b(X) = 0`). So only `b = 1` matrices are visited, and `b`/`T`/`S` + /// never need recomputing. The accumulator `acc` and `xcols` are reused. + fn enum_x(&self, j: usize, xcols: &mut Vec<&'a [u32]>, acc: &mut Acc, out: &mut DualElement) { + if j == self.l { + self.on_x(xcols, acc, out); + return; + } + for cand in self.x_cands[j].iter() { + let fits = cand.iter().enumerate().all(|(i, &v)| { + v == 0 + || (!(i >= 1 && acc.rows[i] + v > self.r1v.get(i).copied().unwrap_or(0)) + && acc.or[i + j] & v == 0) + }); + if !fits { + continue; + } + for (i, &v) in cand.iter().enumerate() { + acc.rows[i] += v; + acc.or[i + j] ^= v; + acc.sum[i + j] += v; + } + xcols.push(cand); + self.enum_x(j + 1, xcols, acc, out); + xcols.pop(); + for (i, &v) in cand.iter().enumerate() { + acc.rows[i] -= v; + acc.or[i + j] ^= v; + acc.sum[i + j] -= v; + } + } + } + + /// A complete `X` (already `b = 1`): form `S′`, the `Y` column targets, and the + /// output P-part `T(X) = acc.sum`, then enumerate `Y`. + fn on_x(&self, xcols: &[&'a [u32]], acc: &Acc, out: &mut DualElement) { + // S′ = R₁ − S(X) on indices ≥ 1 (S′[0] = 0); τ-power = Σ(S′). + let sprime: Vec = (0..self.l) + .map(|j| if j == 0 { 0 } else { self.r1v[j] - acc.rows[j] }) + .collect(); + let sigma2_sprime: i64 = sprime + .iter() + .enumerate() + .map(|(j, &v)| (v as i64) << j) + .sum(); + let sigma_sprime: u32 = sprime.iter().sum(); + + // RY targets: column j ≥ 1 of Y has weighted sum R₂[j] − R(X)[j]. + let mut ry: Vec = vec![0; self.l]; + for (j, ry_j) in ry.iter_mut().enumerate().skip(1) { + let rxj: u32 = xcols[j].iter().enumerate().map(|(i, &v)| v << i).sum(); + *ry_j = self.r2v[j] as i64 - rxj as i64; + if *ry_j < 0 { + return; + } + } + // Σ₂(S(Y)) = Σ₂(E₁) + Σ₂(S′) is forced and equals Σ_{j≥0} R(Y)[j], so column + // 0 of Y has a determined weighted sum. + let ry0 = (self.sigma2_e1 + sigma2_sprime) - ry[1..].iter().sum::(); + if ry0 < 0 { + return; + } + + // The output P-part T(X) (index 0 dropped: ξ₀ = 1). + let out_r = trim( + std::iter::once(0) + .chain((1..NB).map(|d| acc.sum[d])) + .collect(), + ); + + let mut y_cands: Vec>>> = Vec::with_capacity(self.l); + y_cands.push(columns_eq_cached(ry0 as u32)); + for &t in ry.iter().skip(1) { + y_cands.push(columns_eq_cached(t as u32)); + } + + ClosedY { + y_cands: &y_cands, + e1_mask: self.e1_mask, + e2_mask: self.e2_mask, + sprime: &sprime, + sigma_sprime, + out_r: &out_r, + } + .enum_y(0, &mut Acc::zero(), out); + } +} + +/// Immutable context for the `Y` half of the recursion, the counterpart of [`Closed`]. `X` is +/// already fixed by the time this is built, so `sprime` and `out_r` are constants here. +struct ClosedY<'a> { + y_cands: &'a [Rc>>], + e1_mask: u32, + e2_mask: u32, + sprime: &'a [u32], + sigma_sprime: u32, + out_r: &'a [u32], +} + +impl ClosedY<'_> { + /// Enumerate the `Y` matrix column by column (each column an exact weighted sum), + /// pruning on anti-diagonal collisions (`b(Y) = 0`) and accumulating each matching + /// contribution. + fn enum_y(&self, j: usize, acc: &mut Acc, out: &mut DualElement) { + if j == self.y_cands.len() { + self.on_y(acc, out); + return; + } + for cand in self.y_cands[j].iter() { + if !cand + .iter() + .enumerate() + .all(|(i, &v)| acc.or[i + j] & v == 0) + { + continue; + } + for (i, &v) in cand.iter().enumerate() { + acc.rows[i] += v; + acc.or[i + j] ^= v; + acc.sum[i + j] += v; + } + self.enum_y(j + 1, acc, out); + for (i, &v) in cand.iter().enumerate() { + acc.rows[i] -= v; + acc.or[i + j] ^= v; + acc.sum[i + j] -= v; + } + } + } + + /// A complete `Y` (already `b = 1`): accept it if it contributes, and add its term. + fn on_y(&self, acc: &Acc, out: &mut DualElement) { + // The rewriting τ(S(Y)) contains the term with exterior part E₁ and ξ-part + // S′ iff the ρ = 0 degree equation Σ(E₁) + 2Σ(S′) = Σ(S(Y)) holds (with + // Σ(E₁) = popcount) and c(S(Y), S′) ≠ 0. (Σ₂(E₁) + Σ₂(S′) = Σ₂(S(Y)) is + // forced by the Y-column construction.) That term contributes at τ^{Σ(S′)}. + let sigma_sy: u32 = acc.rows.iter().sum(); + if self.e1_mask.count_ones() + 2 * self.sigma_sprime != sigma_sy { + return; + } + let sy_len = (0..NB) + .rev() + .find(|&i| acc.rows[i] > 0) + .map_or(0, |i| i + 1); + if c_coeff(&acc.rows[..sy_len], self.sprime) == 0 { + return; + } + // E_out = E₂ + T(Y) must be square-free (Seq₁). Q indices fit in a u32. + let mut out_e_mask = 0u32; + for i in 0..u32::BITS as usize { + let val = ((self.e2_mask >> i) & 1) + acc.sum[i]; + if val > 1 { + return; + } + if val == 1 { + out_e_mask |= 1 << i; + } + } + add_term( + out, + (out_e_mask, self.out_r.to_vec()), + Tau::power(self.sigma_sprime), + ); + } +} + +/// The product `a · b` in `A_C` via Kong–Lin Theorem 5.1 (ρ = 0). Same contract +/// as [`multiply`] (the duality oracle it is validated against). +pub fn multiply_closed(a: &Monomial, b: &Monomial) -> DualElement { + let (e1_mask, r1) = (a.0, &a.1); + let (e2_mask, r2) = (b.0, &b.1); + let l = r1.len().max(r2.len()).max(1); + let r1v: Vec = (0..l).map(|j| r1.get(j).copied().unwrap_or(0)).collect(); + let r2v: Vec = (0..l).map(|j| r2.get(j).copied().unwrap_or(0)).collect(); + + // Per-column candidate columns for X: column 0 bounded by R₁ rows, columns + // j ≥ 1 by R₂[j] weighted. + let mut x_cands: Vec>>> = Vec::with_capacity(l); + x_cands.push(Rc::new(col0_x_options(&r1v))); + for &bound in r2v.iter().skip(1) { + x_cands.push(columns_le_cached(bound)); + } + + let ctx = Closed { + r1v: &r1v, + r2v: &r2v, + l, + e1_mask, + e2_mask, + sigma2_e1: e1_mask as i64, // Σ (bit i)·2ⁱ = e1_mask + x_cands: &x_cands, + }; + let mut out = DualElement::new(); + let mut xcols: Vec<&[u32]> = Vec::with_capacity(l); + let mut acc = Acc::zero(); + ctx.enum_x(0, &mut xcols, &mut acc, &mut out); + out +} + +// --------------------------------------------------------------------------- +// The C-motivic Steenrod algebra A_C as a free 𝔽₂[τ]-module on the Milnor basis. +// --------------------------------------------------------------------------- + +/// The C-motivic (prime 2) Steenrod algebra $A_C$, presented as a free +/// $\mathbb{F}_2[\tau]$-module on the Milnor basis $\{Q(E)P(R)\}$ with lazy +/// per-(topological-)degree basis indexing. +/// +/// This is the $A_C$ product **engine** over $\mathbb{F}_2[\tau]$: the dual-based product is +/// exposed through [`product_indexed`](MotivicMilnorAlgebra::product_indexed), whose coefficients +/// are [`Tau`] scalars. It is deliberately not an [`Algebra`](crate::algebra::Algebra) +/// implementation, because that trait is over $\mathbb{F}_p$; the mod-$\tau$ reduction, which *is* +/// such an algebra, is a follow-up layer on top of this engine, and the honest +/// $\mathbb{F}_2[\tau]$ resolution is built by lifting against this engine (Phase 2). +/// +/// Weight convention: [`bidegree`](MotivicMilnorAlgebra::bidegree) returns `(t, w)` where `t` is +/// the topological degree and `w` is the motivic weight in the presentation where $\tau$ has +/// weight $-1$ and products are weight-homogeneous (i.e. `w = -(dual monomial weight)`; see the +/// module-level note). +#[derive(Default)] +pub struct MotivicMilnorAlgebra { + /// `basis[t]` is the $\mathbb{F}_2[\tau]$-basis in topological degree `t`, sorted for stable + /// indexing. + basis: OnceVec>, + /// Memoized basis-element products, one dense [`ProductBlock`] per pair of + /// topological degrees `(t1, t2)`. The duality product is expensive and a + /// resolution asks for the same structure constants repeatedly, so we cache + /// them (the role the classical Milnor algebra's `cache-multiplication` table + /// plays). Blocking by degree pair — rather than a flat `(t1, idx1, t2, idx2)` + /// map — makes the cache the natural *batch unit* (see [`ProductBlock`]): the + /// per-degree-pair registry is small and read-mostly, and each block's entries + /// fill through independent [`OnceLock`]s, so a concurrent resolution reads + /// hits lock-free instead of serializing on one global product lock. + blocks: RwLock>>, +} + +/// A dense block of basis-element products for one pair of topological degrees +/// `(t1, t2)`. Entry `(idx1, idx2)`, at flat position `idx1 * dim2 + idx2`, is the +/// product of the `idx1`-th basis element in degree `t1` with the `idx2`-th in +/// degree `t2` — the `(Tau, index-in-degree-(t1+t2))` list [`multiply_closed`] +/// returns. +/// +/// This is the natural **batch unit** for the product: every structure constant +/// the resolution needs for a given degree pair lives in one block. Entries fill +/// lazily (each an independent [`OnceLock`], so reads are lock-free and a +/// concurrent resolution never serializes on a global product lock), and +/// [`MotivicMilnorAlgebra::fill_block`] computes the whole block at once — the +/// shape a GPU kernel would fill in a single launch. +pub struct ProductBlock { + /// The number of basis elements in degree `t2` (the block's row stride). + dim2: usize, + entries: Vec>>, +} + +impl MotivicMilnorAlgebra { + pub fn new() -> Self { + Self { + basis: OnceVec::new(), + blocks: RwLock::new(FxHashMap::default()), + } + } + + /// Compute and cache the basis in every topological degree up to and including `degree`. + /// Idempotent and cheap to re-call. + pub fn compute_basis(&self, degree: i32) { + for t in self.basis.len() as i32..=degree { + let mut b = enum_basis(t); + b.sort(); + self.basis.push(b); + } + } + + /// The $\mathbb{F}_2[\tau]$-rank of `A_C` in topological degree `degree`. + pub fn dimension(&self, degree: i32) -> usize { + if degree < 0 { + return 0; + } + self.basis[degree as usize].len() + } + + /// The `idx`-th basis monomial `(E, R)` in degree `degree`. + pub fn basis_element(&self, degree: i32, idx: usize) -> &Monomial { + &self.basis[degree as usize][idx] + } + + /// The index of a basis monomial in its degree, if present. + pub fn index_of(&self, degree: i32, elt: &Monomial) -> Option { + self.basis[degree as usize].binary_search(elt).ok() + } + + /// The `(topological degree, motivic weight)` of a basis element. + pub fn bidegree(&self, degree: i32, idx: usize) -> (i32, i32) { + let (e, r) = &self.basis[degree as usize][idx]; + let (t, w) = paper_bidegree(*e, r); + (t, -w) + } + + /// Parse a basis element written as `basis_element_to_string` prints it — + /// `1`, `Q_i`, `P(r_0, r_1, …)`, or a space-separated product `Q_i … P(R)` — + /// into its `(topological degree, index)`. The inverse of + /// [`Self::basis_element_to_string`]; returns `None` on a malformed string or a + /// monomial absent from the basis. This is what lets `.json` module descriptors + /// be written over the motivic Steenrod algebra. + pub fn basis_element_from_string(&self, elt: &str) -> Option<(i32, usize)> { + let mut e_mask: u32 = 0; + let mut r: Vec = Vec::new(); + // The `P(…)` block (if any) is last and may contain spaces after commas, so + // split it off before whitespace-tokenizing the `Q_i` factors. + let (q_part, p_part) = match elt.split_once("P(") { + Some((q, rest)) => (q, Some(rest.strip_suffix(')')?)), + None => (elt, None), + }; + for token in q_part.split_whitespace() { + if token == "1" { + continue; // the unit factor + } + let i: u32 = token.strip_prefix("Q_")?.parse().ok()?; + e_mask |= 1 << i; + } + if let Some(p) = p_part { + r = p + .split(',') + .map(|x| x.trim().parse::()) + .collect::>() + .ok()?; + } + + let degree = paper_bidegree(e_mask, &r).0; + if degree < 0 { + return None; + } + self.compute_basis(degree); + // The basis is stored trimmed, so trimming the parsed vector makes this an exact match. + self.index_of(degree, &(e_mask, trim(r))) + .map(|idx| (degree, idx)) + } + + /// A display string for a basis element (`Q_i … P(R)`). + pub fn basis_element_to_string(&self, degree: i32, idx: usize) -> String { + let (e_mask, r) = self.basis_element(degree, idx); + let mut parts = Vec::new(); + for i in fp::prime::iter::BitflagIterator::set_bit_iterator(*e_mask as u64) { + parts.push(format!("Q_{i}")); + } + if !r.is_empty() { + parts.push(format!("P({})", r.iter().format(", "))); + } + if parts.is_empty() { + "1".to_string() + } else { + parts.join(" ") + } + } + + /// The dense [`ProductBlock`] for the degree pair `(t1, t2)`, created (with all + /// entries empty) on first request and shared thereafter. Creation is the only + /// step that touches the block registry's write lock; entry computation happens + /// lock-free through the block's [`OnceLock`]s. + fn block(&self, t1: i32, t2: i32) -> Arc { + if let Some(block) = self.blocks.read().unwrap().get(&(t1, t2)) { + return Arc::clone(block); + } + // Ensure the operand and output bases exist before sizing/indexing. + self.compute_basis(t1); + self.compute_basis(t2); + self.compute_basis(t1 + t2); + let dim1 = self.dimension(t1); + let dim2 = self.dimension(t2); + let mut w = self.blocks.write().unwrap(); + Arc::clone(w.entry((t1, t2)).or_insert_with(|| { + Arc::new(ProductBlock { + dim2, + entries: (0..dim1 * dim2).map(|_| OnceLock::new()).collect(), + }) + })) + } + + /// The closed-form product (Kong–Lin Theorem 5.1) of the two basis elements at + /// `(t1, idx1)` and `(t2, idx2)`, mapped into the degree-`(t1+t2)` basis. The + /// bases must already be computed (the caller through [`Self::block`] ensures + /// this). + fn compute_product(&self, t1: i32, idx1: usize, t2: i32, idx2: usize) -> Vec<(Tau, usize)> { + let a = &self.basis[t1 as usize][idx1]; + let b = &self.basis[t2 as usize][idx2]; + let t = t1 + t2; + multiply_closed(a, b) + .into_iter() + .map(|(z, c)| { + ( + c, + self.index_of(t, &z) + .expect("product landed outside the basis"), + ) + }) + .collect() + } + + /// The product of two basis elements, as an $\mathbb{F}_2[\tau]$-linear combination of basis + /// elements in degree `t1 + t2`: a list of `(coefficient, index)` pairs. The product of two + /// homogeneous basis elements is weight-homogeneous, so each coefficient is a single power of + /// $\tau$ — a [`Tau`] scalar. + pub fn product_indexed(&self, t1: i32, idx1: usize, t2: i32, idx2: usize) -> Vec<(Tau, usize)> { + self.cached_product(&self.block(t1, t2), t1, idx1, t2, idx2) + .to_vec() + } + + /// One entry of `block`, computed on first request. `block` must be the block for + /// `(t1, t2)`; it is passed in so [`Self::fill_block`] can look it up once for a whole + /// block instead of per entry. + fn cached_product<'a>( + &self, + block: &'a ProductBlock, + t1: i32, + idx1: usize, + t2: i32, + idx2: usize, + ) -> &'a [(Tau, usize)] { + block.entries[idx1 * block.dim2 + idx2] + .get_or_init(|| self.compute_product(t1, idx1, t2, idx2)) + } + + /// Compute *every* entry of the `(t1, t2)` product block at once, filling it in + /// parallel (under the `concurrent` feature) — the batch shape a GPU kernel + /// fills in a single launch. Idempotent: entries already present are kept, and + /// a second call is a cheap no-op. This is the batching boundary the GPU port + /// hooks into; on CPU it lets a caller warm a whole degree pair up front + /// instead of paying a lock-free-but-serial miss per structure constant. + pub fn fill_block(&self, t1: i32, t2: i32) { + let block = self.block(t1, t2); + let dim2 = block.dim2; + if dim2 == 0 || block.entries.is_empty() { + return; + } + let dim1 = block.entries.len() / dim2; + (0..dim1).into_maybe_par_iter().for_each(|idx1| { + for idx2 in 0..dim2 { + self.cached_product(&block, t1, idx1, t2, idx2); + } + }); + } +} + +impl std::fmt::Display for MotivicMilnorAlgebra { + fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { + write!(f, "MotivicMilnorAlgebra(C, p=2)") + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// The classical mod-2 Milnor product of `P(a) * P(b)`, as a list of exponent sequences. + /// + /// The motivic layer works with the paper's `R = [r_0, r_1, ...]` indexing, where the + /// leading slot belongs to ξ_0 = 1; the classical [`PPart`] starts at r_1. Callers convert + /// with [`pp_to_paper`]/[`paper_to_pp`]. + fn classical_mul(a: &[u32], b: &[u32]) -> Vec> { + use fp::{prime::TWO, vector::FpVector}; + + use crate::algebra::{ + Algebra, + milnor_algebra::{MilnorAlgebra, MilnorBasisElement, PPart}, + }; + + // `MilnorAlgebra` is cheap to build but caches its basis, so keep one per process + // rather than recomputing the tables for every product in a loop. + thread_local! { + static ALG: MilnorAlgebra = MilnorAlgebra::new(TWO, false); + } + ALG.with(|alg| { + let mk = |p: &[u32]| { + let mut m = MilnorBasisElement { + q_part: 0, + p_part: PPart::try_from_slice(p).unwrap(), + degree: 0, + }; + m.compute_degree(TWO); + m + }; + let (m1, m2) = (mk(a), mk(b)); + let degree = m1.degree + m2.degree; + alg.compute_basis(degree); + let mut res = FpVector::new(TWO, alg.dimension(degree)); + alg.multiply(res.as_slice_mut(), 1, m1, m2); + res.iter_nonzero() + // `PPart::iter` stops at the last non-zero entry, so this is already trimmed. + .map(|(idx, _)| { + alg.basis_element_from_index(degree, idx) + .p_part + .iter() + .collect() + }) + .collect() + }) + } + + /// A classical exponent sequence as the paper's `R`, with the ξ_0 slot prepended. + fn pp_to_paper(pp: &[u32]) -> Vec { + let mut r = vec![0]; + r.extend_from_slice(pp); + trim(r) + } + + /// The inverse of [`pp_to_paper`]: drop the ξ_0 slot. + fn paper_to_pp(r: &[u32]) -> Vec { + trim(r.get(1..).unwrap_or_default().to_vec()) + } + + #[test] + fn test_algebra_basis_and_multiply() { + let alg = MotivicMilnorAlgebra::new(); + alg.compute_basis(8); + + // Degree 0 is the unit; degrees 1 and 2 are 1-dimensional (Q_0 and P(ξ_1)). + assert_eq!(alg.dimension(0), 1); + assert_eq!(alg.basis_element(0, 0), &(0u32, vec![])); + assert_eq!(alg.dimension(1), 1); + assert_eq!(alg.basis_element(1, 0), &(0b1u32, vec![])); // Q_0 + assert_eq!(alg.dimension(2), 1); + assert_eq!(alg.basis_element(2, 0), &(0u32, vec![0, 1])); // P(ξ_1) + + // bidegree: Q_0 is (1, 0), P(ξ_1) is (2, -1) in this presentation. + assert_eq!(alg.bidegree(1, 0), (1, 0)); + assert_eq!(alg.bidegree(2, 0), (2, -1)); + + // Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1), reconstructed from indices. + let terms: DualElement = alg + .product_indexed(1, 0, 2, 0) + .into_iter() + .map(|(c, idx)| (alg.basis_element(3, idx).clone(), c)) + .collect(); + assert_eq!( + terms, + DualElement::from([ + ((0b10, vec![]), Tau::one()), + ((0b1, vec![0, 1]), Tau::one()) + ]) + ); + } + + #[test] + fn test_algebra_multiply_matches_raw_and_is_homogeneous() { + let alg = MotivicMilnorAlgebra::new(); + alg.compute_basis(12); + for t1 in 0..=6 { + for idx1 in 0..alg.dimension(t1) { + for t2 in 0..=6 { + for idx2 in 0..alg.dimension(t2) { + let a = alg.basis_element(t1, idx1).clone(); + let b = alg.basis_element(t2, idx2).clone(); + // Indexed product agrees with the raw monomial product. + let indexed: DualElement = alg + .product_indexed(t1, idx1, t2, idx2) + .into_iter() + .map(|(c, idx)| (alg.basis_element(t1 + t2, idx).clone(), c)) + .collect(); + assert_eq!(indexed, multiply(&a, &b)); + + // Weight-homogeneous: w_out - (τ-power) = w_a + w_b. + let (wa, wb) = (alg.bidegree(t1, idx1).1, alg.bidegree(t2, idx2).1); + for (c, idx) in alg.product_indexed(t1, idx1, t2, idx2) { + let w = alg.bidegree(t1 + t2, idx).1; + assert_eq!(w - c.valuation().unwrap() as i32, wa + wb); + } + } + } + } + } + } + + #[test] + fn test_dual_mul_relations() { + // τ_0^2 = τ ξ_1 (the defining relation at ρ = 0). + assert_eq!( + dual_mul(&tau_gen(0), &tau_gen(0)), + DualElement::from([((0, vec![0, 1]), Tau::power(1))]) // τ^1 · ξ_1 + ); + // ξ_1^2 is just the monomial ξ_1^2. + assert_eq!( + dual_mul(&xi_gen(1), &xi_gen(1)), + DualElement::from([((0, vec![0, 2]), Tau::one())]) + ); + // Distinct τ's commute and stay square-free: τ_0 τ_1. + assert_eq!( + dual_mul(&tau_gen(0), &tau_gen(1)), + DualElement::from([((0b11, vec![]), Tau::one())]) + ); + // τ_1^2 = τ ξ_2. + assert_eq!( + dual_mul(&tau_gen(1), &tau_gen(1)), + DualElement::from([((0, vec![0, 0, 1]), Tau::power(1))]) + ); + // Multiplication by the unit is the identity. + assert_eq!(dual_mul(&dual_one(), &tau_gen(2)), tau_gen(2)); + } + + #[test] + fn test_pure_xi_matches_classical_milnor() { + // Motivically the even part is NOT closed: P(R_1)·P(R_2) can leak into Q-terms + // carrying τ (e.g. Sq^2·Sq^2 ∋ τ Q_0 Q_1). But the pure-ξ *outputs* (E = 0) are τ-free + // and must agree with the ordinary mod-2 Milnor product. Compare those against the + // codebase's `MilnorAlgebra` (paper R = [0, p_part...]; drop the ξ_0 slot to convert). + let classical = |p1: &[u32], p2: &[u32]| -> Vec> { + let mut out = classical_mul(p1, p2); + out.sort(); + out + }; + + let motivic = |p1: &[u32], p2: &[u32]| -> Vec> { + let mut out: Vec> = multiply(&(0, pp_to_paper(p1)), &(0, pp_to_paper(p2))) + .into_iter() + .filter(|((e, _), _)| *e == 0) // keep the pure-ξ outputs + .map(|((_, r), c)| { + assert_eq!(c, Tau::one(), "pure-ξ output carried a τ power"); + paper_to_pp(&r) + }) + .collect(); + out.sort(); + out + }; + + // The paper uses *conjugate* generators, so P(paper R) differs from the codebase's + // standard Milnor basis by the antipode χ for ξ_2 and higher (e.g. χ(ξ_2)=ξ_2+ξ_1^3). + // They agree on the self-conjugate ξ_1-power sector (ξ_1 is primitive), so restrict the + // comparison to outputs that are pure powers of ξ_1 (a `p_part` of length ≤ 1). + let xi1_only = |mut v: Vec>| { + v.retain(|pp| pp.len() <= 1); + v + }; + for (p1, p2) in [ + (vec![1u32], vec![1u32]), + (vec![2], vec![1]), + (vec![1], vec![2]), + (vec![3], vec![1]), + (vec![2, 1], vec![1]), + (vec![0, 1], vec![0, 1]), + (vec![2], vec![2]), + (vec![3], vec![3]), + ] { + assert_eq!( + xi1_only(motivic(&p1, &p2)), + xi1_only(classical(&p1, &p2)), + "ξ_1-sector mismatch for {p1:?} * {p2:?}" + ); + } + } + + #[test] + fn test_full_xi_matches_classical_via_antipode() { + // Full (all-ξ) cross-check against the codebase's classical Milnor product, using the + // antipode χ to reconcile the paper's conjugate generators with the codebase's standard + // ones. Identity checked (pure-ξ, τ-free, so over 𝔽₂): expressing both inputs and + // outputs of the motivic conjugate product in the standard basis via χ reproduces the + // codebase product. `A[μ][W] = coeff of conj-mon μ in std ξ(W) = [χ(mon_W)]_μ`. + let mdeg = |pp: &[u32]| paper_bidegree(0, &pp_to_paper(pp)).0; + let xi_pps = |t: i32| -> Vec> { + enum_basis(t) + .into_iter() + .filter(|(e, _)| *e == 0) + .map(|(_, r)| paper_to_pp(&r)) + .collect() + }; + let a_coeff = |mu_pp: &[u32], w_pp: &[u32]| -> u32 { + let ap = antipode(&DualElement::from([( + (0u32, pp_to_paper(w_pp)), + Tau::one(), + )])); + // The antipode of a pure-ξ element is pure-ξ and τ-free, so a present monomial has + // coefficient τ^0. + ap.get(&(0u32, pp_to_paper(mu_pp))) + .map_or(0, |&c| u32::from(c == Tau::one())) + }; + let cb_mul = |a_pp: &[u32], b_pp: &[u32], acc: &mut BTreeMap, u32>| { + for pp in classical_mul(a_pp, b_pp) { + *acc.entry(pp).or_insert(0) ^= 1; + } + }; + + for (r1, r2) in [ + (vec![2u32], vec![1u32]), + (vec![1], vec![2]), + (vec![0, 1], vec![1]), + (vec![2], vec![2]), + (vec![3], vec![1]), + (vec![1, 1], vec![1]), + (vec![4], vec![1]), + (vec![2, 1], vec![2]), + ] { + // LHS: motivic conjugate product (E=0), converted to std via χ. + let conj: Vec> = multiply(&(0, pp_to_paper(&r1)), &(0, pp_to_paper(&r2))) + .into_iter() + .filter(|((e, _), _)| *e == 0) + .map(|((_, r), _)| paper_to_pp(&r)) + .collect(); + let mut lhs: BTreeMap, u32> = BTreeMap::new(); + for w in xi_pps(mdeg(&r1) + mdeg(&r2)) { + let s = conj.iter().fold(0u32, |acc, mu| acc ^ a_coeff(mu, &w)); + if s != 0 { + lhs.insert(w, s); + } + } + + // RHS: convert inputs to std via χ, multiply in the codebase. + let in1: Vec> = xi_pps(mdeg(&r1)) + .into_iter() + .filter(|w| a_coeff(&r1, w) == 1) + .collect(); + let in2: Vec> = xi_pps(mdeg(&r2)) + .into_iter() + .filter(|w| a_coeff(&r2, w) == 1) + .collect(); + let mut rhs: BTreeMap, u32> = BTreeMap::new(); + for w1 in &in1 { + for w2 in &in2 { + cb_mul(w1, w2, &mut rhs); + } + } + rhs.retain(|_, v| *v != 0); + + assert_eq!(lhs, rhs, "full ξ cross-check failed for {r1:?} * {r2:?}"); + } + } + + #[test] + fn test_product_associative() { + // Associativity (a·b)·c = a·(b·c) is convention-independent and a strong global check. + // Extend the basis×basis product to element×basis, F_2[τ]-linearly (add_term handles the + // coefficient sum + zero cancellation via Tau arithmetic). + fn mul_elt_basis(x: &DualElement, c: &Monomial) -> DualElement { + let mut out = DualElement::new(); + for (z, &cz) in x { + for (w, cw) in multiply(z, c) { + add_term(&mut out, w, cz * cw); + } + } + out + } + + let basis: Vec = (0..=5).flat_map(enum_basis).collect(); + for a in &basis { + for b in &basis { + let ab = multiply(a, b); + for c in &basis { + let lhs = mul_elt_basis(&ab, c); + let bc = multiply(b, c); + // a · (b·c): multiply basis `a` on the left of each term of bc. + let mut rhs = DualElement::new(); + for (z, &cz) in &bc { + for (w, cw) in multiply(a, z) { + add_term(&mut rhs, w, cz * cw); + } + } + assert_eq!(lhs, rhs, "associativity failed at {a:?},{b:?},{c:?}"); + } + } + } + } + + #[test] + fn test_product_q0_p_xi1() { + // The case my reading of Kong–Lin Theorem 5.1 got wrong. By duality: + // Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1) (both with coefficient τ^0). + let q0 = (0b1, vec![]); + let p_xi1 = (0, vec![0, 1]); + assert_eq!( + multiply(&q0, &p_xi1), + DualElement::from([ + ((0b10, vec![]), Tau::one()), // Q_1 + ((0b1, vec![0, 1]), Tau::one()), // Q_0 P(ξ_1) + ]) + ); + } + + #[test] + fn test_product_unit_and_squares() { + // 1 · x = x. + let x = (0b101, vec![0, 2]); + assert_eq!( + multiply(&(0, vec![]), &x), + DualElement::from([(x.clone(), Tau::one())]) + ); + // Q_i^2 = 0. (Q_0 = Sq^1 is the motivic Bockstein; P(ξ_1) = Sq^2 does NOT square to + // zero — its square is τ Q_0 Q_1 + …, a genuine motivic feature.) + for i in 0..3 { + assert!( + multiply(&(1 << i, vec![]), &(1 << i, vec![])).is_empty(), + "Q_{i}^2 ≠ 0" + ); + } + } + + #[test] + fn test_product_weight_homogeneous_and_tau_appears() { + // Every product is weight-homogeneous: term z occurs with a single power τ^k, and by + // weight-preservation of ψ, k = (w_a + w_b) - w_z ≥ 0 in *dual*-monomial weights (the + // algebra weight of a basis element is the negative of its dual monomial's weight, so + // τ, of weight -1, absorbs the difference). Also confirm τ genuinely enters at least + // one product (i.e. the τ_i^2 = τξ_{i+1} relation fires). + let basis: Vec = (0..=6).flat_map(enum_basis).collect(); + let mut saw_tau = false; + for a in &basis { + for b in &basis { + let (wa, wb) = (paper_bidegree(a.0, &a.1).1, paper_bidegree(b.0, &b.1).1); + for (z, c) in multiply(a, b) { + // Each coefficient is a single τ power by construction (a homogeneous Tau). + let k = c.valuation().expect("stored coefficients are never zero") as i32; + let wz = paper_bidegree(z.0, &z.1).1; + assert_eq!(wa + wb - wz, k, "weight mismatch: {a:?}*{b:?} → z={z:?}"); + assert!(k >= 0); + saw_tau |= k > 0; + } + } + } + assert!( + saw_tau, + "no product produced a τ coefficient — τ_i^2 relation never exercised" + ); + } + + #[test] + fn test_antipode() { + // χ(ξ_1) = ξ_1 (primitive); χ(ξ_2) = ξ_2 + ξ_1^3; χ(τ_0) = τ_0. + assert_eq!(antipode(&xi_gen(1)), xi_gen(1)); + assert_eq!( + antipode(&xi_gen(2)), + DualElement::from([ + ((0, vec![0, 0, 1]), Tau::one()), + ((0, vec![0, 3]), Tau::one()) + ]) + ); + assert_eq!(antipode(&tau_gen(0)), tau_gen(0)); + + // χ is an involution: χ(χ(x)) = x. + for x in [ + xi_gen(1), + xi_gen(2), + xi_gen(3), + tau_gen(0), + tau_gen(1), + tau_gen(2), + dual_mul(&tau_gen(0), &xi_gen(2)), + ] { + assert_eq!(antipode(&antipode(&x)), x, "χ not an involution on {x:?}"); + } + + // Antipode axiom: m(χ ⊗ id)ψ(x) = η ε(x); for positive-degree x this is 0. + let counit_via_antipode = |x: &DualElement| -> DualElement { + let mut out = DualElement::new(); + for ((l, r), c) in coproduct(x) { + let prod = dual_mul( + &antipode(&DualElement::from([(l, Tau::one())])), + &DualElement::from([(r, Tau::one())]), + ); + for (mon, cc) in prod { + add_term(&mut out, mon, c * cc); + } + } + out + }; + for x in [ + xi_gen(2), + xi_gen(3), + tau_gen(1), + tau_gen(2), + dual_mul(&tau_gen(0), &tau_gen(1)), + ] { + assert!( + counit_via_antipode(&x).is_empty(), + "antipode axiom failed on {x:?}" + ); + } + } + + #[test] + fn test_coproduct_generators() { + // ψ(τ_1) = 1⊗τ_1 + τ_0⊗ξ_1 + τ_1⊗1 (Kong–Lin §2.2). + assert_eq!( + coproduct(&tau_gen(1)), + TensorElement::from([ + (((0, vec![]), (0b10, vec![])), Tau::one()), // 1 ⊗ τ_1 + (((0b1, vec![]), (0, vec![0, 1])), Tau::one()), // τ_0 ⊗ ξ_1 + (((0b10, vec![]), (0, vec![])), Tau::one()), // τ_1 ⊗ 1 + ]) + ); + // ψ(ξ_2) = 1⊗ξ_2 + ξ_1⊗ξ_1^2 + ξ_2⊗1 (matches Milnor). + assert_eq!( + coproduct(&xi_gen(2)), + TensorElement::from([ + (((0, vec![]), (0, vec![0, 0, 1])), Tau::one()), // 1 ⊗ ξ_2 + (((0, vec![0, 1]), (0, vec![0, 2])), Tau::one()), // ξ_1 ⊗ ξ_1^2 + (((0, vec![0, 0, 1]), (0, vec![])), Tau::one()), // ξ_2 ⊗ 1 + ]) + ); + } + + #[test] + fn test_coproduct_is_algebra_map() { + // ψ(xy) = ψ(x)ψ(y); in particular this checks that the τ_i^2 = τξ_{i+1} reduction is + // compatible with the (τ-free) coproduct. + for (x, y) in [ + (tau_gen(0), tau_gen(0)), + (tau_gen(1), tau_gen(1)), + (tau_gen(0), tau_gen(1)), + (xi_gen(1), xi_gen(1)), + (tau_gen(0), xi_gen(1)), + (tau_gen(1), xi_gen(2)), + ] { + let lhs = coproduct(&dual_mul(&x, &y)); + let rhs = tensor_mul(&coproduct(&x), &coproduct(&y)); + assert_eq!(lhs, rhs, "ψ not an algebra map on {x:?} * {y:?}"); + } + } + + #[test] + fn test_rewrite_tau_square_free_is_identity() { + // A square-free τ(E) is already a basis element: τ_0 τ_2 → Q(E)P(0), τ^0. + let terms = rewrite_tau(&[1, 0, 1]); + assert_eq!( + terms, + vec![MotivicTerm { + tau_pow: 0, + e_mask: 0b101, + r: vec![], + }] + ); + } + + #[test] + fn test_rewrite_tau_example_3_1() { + // Kong–Lin Example 3.1, specialized to ρ = 0 (C-motivic): + // τ_0^2 τ_1 = τ · τ_1 ξ_1 + // Encoded: S = (2, 1); expected single term τ^1 · Q({1}) P((0,1)). + let terms = rewrite_tau(&[2, 1]); + assert_eq!( + terms, + vec![MotivicTerm { + tau_pow: 1, + e_mask: 0b10, // τ_1 + r: vec![0, 1], // ξ_1^1 + }], + "τ_0^2 τ_1 = τ τ_1 ξ_1 (ρ=0)" + ); + + // τ_0^4 = τ^2 ξ_1^2 + // Encoded: S = (4); expected single term τ^2 · Q(∅) P((0,2)). + let terms = rewrite_tau(&[4]); + assert_eq!( + terms, + vec![MotivicTerm { + tau_pow: 2, + e_mask: 0, + r: vec![0, 2], // ξ_1^2 + }], + "τ_0^4 = τ^2 ξ_1^2 (ρ=0)" + ); + } + + #[test] + fn test_closed_form_small_cases() { + // The case a naive reading of Theorem 5.1 gets wrong (the ξ₀ = 1 index + // absorption): Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1). + let q0 = (0b1, vec![]); + let p_xi1 = (0, vec![0, 1]); + assert_eq!(multiply_closed(&q0, &p_xi1), multiply(&q0, &p_xi1)); + + // The τ-generating case: P(ξ_1)² has a τ Q_0 Q_1 term. + assert_eq!(multiply_closed(&p_xi1, &p_xi1), multiply(&p_xi1, &p_xi1)); + + // Q_i² = 0. + for i in 0..3 { + assert_eq!( + multiply_closed(&(1 << i, vec![]), &(1 << i, vec![])), + multiply(&(1 << i, vec![]), &(1 << i, vec![])) + ); + } + } + + #[test] + fn test_closed_form_matches_duality_oracle() { + // Exhaustive fuzz: the closed form (Theorem 5.1) agrees with the duality + // oracle on every ordered pair of basis elements with deg(a)+deg(b) ≤ 18, + // including the τ-carrying and multi-Q cases. (The bound is kept modest + // because the *oracle* is slow; the closed form is not.) + let basis: Vec<(i32, Monomial)> = (0..=14) + .flat_map(|t| enum_basis(t).into_iter().map(move |m| (t, m))) + .collect(); + for (ta, a) in &basis { + for (tb, b) in &basis { + if ta + tb > 18 { + continue; + } + assert_eq!( + multiply_closed(a, b), + multiply(a, b), + "closed form ≠ oracle for {a:?} · {b:?}" + ); + } + } + } + + #[test] + fn test_fill_block_matches_lazy_products() { + // The batch path (`fill_block`, which a GPU kernel would drive) must agree + // entry-for-entry with the lazy per-product path. Fill whole degree-pair + // blocks up front on one algebra, compute the same products lazily on + // another, and compare every structure constant. + let batched = MotivicMilnorAlgebra::new(); + let lazy = MotivicMilnorAlgebra::new(); + batched.compute_basis(12); + lazy.compute_basis(12); + for t1 in 0..=6 { + for t2 in 0..=6 { + batched.fill_block(t1, t2); + for idx1 in 0..batched.dimension(t1) { + for idx2 in 0..batched.dimension(t2) { + assert_eq!( + batched.product_indexed(t1, idx1, t2, idx2), + lazy.product_indexed(t1, idx1, t2, idx2), + "batch ≠ lazy at ({t1},{idx1})·({t2},{idx2})" + ); + } + } + } + } + } + + #[test] + fn test_rewrite_tau_unit() { + assert_eq!( + rewrite_tau(&[]), + vec![MotivicTerm { + tau_pow: 0, + e_mask: 0, + r: vec![], + }] + ); + } +} diff --git a/ext/crates/algebra/src/algebra/motivic/mod.rs b/ext/crates/algebra/src/algebra/motivic/mod.rs new file mode 100644 index 0000000000..94f1fa1bd3 --- /dev/null +++ b/ext/crates/algebra/src/algebra/motivic/mod.rs @@ -0,0 +1,28 @@ +//! The C-motivic (over $\mathbb{C}$, prime 2) Steenrod algebra and its mod-$\tau$ +//! reduction. +//! +//! This layer implements the *deformation* view of the C-motivic Adams $E_2$: the +//! C-motivic dual Steenrod algebra $A_C$ over $\mathbb{F}_2[\tau]$, its mod-$\tau$ +//! reduction $A_C/\tau$ (a connected finite-type $\mathbb{F}_2$-algebra), and the +//! coefficient ring $\mathbb{F}_2[\tau]$ itself. +//! +//! This module provides the foundation layer for that computation: +//! +//! - [`Tau`] — the coefficient ring $\mathbb{F}_2[\tau]$ as a small homogeneous +//! scalar ([`tau`]); the whole $\tau$-tower is carried here rather than threaded +//! through the resolution engine. +//! - [`MotivicMilnorAlgebra`] — $A_C$, a free $\mathbb{F}_2[\tau]$-module on the +//! Milnor basis, with the Kong–Lin product ([`milnor`]). This is the product +//! **engine**; it is deliberately not an [`Algebra`](crate::algebra::Algebra) +//! (that trait is over $\mathbb{F}_p$). +//! +//! The mod-$\tau$ reduction $A_C/\tau$ — the connected finite-type +//! $\mathbb{F}_2$-algebra that the existing resolution engine resolves to yield the +//! algebraic Novikov $E_2$ — is presented as an [`Algebra`](crate::algebra::Algebra) +//! in a follow-up on top of this engine. + +pub mod tau; +pub use tau::Tau; + +pub mod milnor; +pub use milnor::MotivicMilnorAlgebra; diff --git a/ext/crates/algebra/src/algebra/motivic/tau.rs b/ext/crates/algebra/src/algebra/motivic/tau.rs new file mode 100644 index 0000000000..1a6e53c526 --- /dev/null +++ b/ext/crates/algebra/src/algebra/motivic/tau.rs @@ -0,0 +1,142 @@ +//! The coefficient ring $\mathbb{F}_2[\tau]$, represented by $\tau$-valuations. +//! +//! Every computation over the C-motivic Steenrod algebra that we perform is +//! **homogeneous** — differentials, products, and coproducts all preserve the +//! motivic weight — and a homogeneous element of $\mathbb{F}_2[\tau]$ in a fixed +//! bidegree is exactly $\tau^k$ for a single $k$ (or $0$). So we never need a +//! polynomial: a coefficient is one integer, its $\tau$-valuation, with a +//! sentinel for zero. That is [`Tau`]. +//! +//! This is a small, self-contained scalar type for the motivic layer only. It is +//! deliberately **not** an [`fp`] replacement threaded through the resolution +//! engine: the engine stays $\mathbb{F}_p$-only (see the `motivic` module docs), +//! and `Tau` is used solely by the motivic algebra and the Phase 2 lift. + +/// A homogeneous $\mathbb{F}_2[\tau]$ coefficient: either $0$ or $\tau^k$ for a +/// unique $k \ge 0$, stored as the $\tau$-valuation with a zero sentinel. +/// +/// Over $\mathbb{F}_2$ a nonzero homogeneous coefficient has unit $\mathbb{F}_2$ +/// part, so the valuation is the whole story: `None` is $0$, `Some(k)` is +/// $\tau^k$. Multiplication adds valuations; addition of two coefficients in the +/// *same* bidegree can only ever combine equal powers (homogeneity), so it +/// cancels mod 2. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] +pub struct Tau(Option); + +impl Tau { + /// The unit $\tau^0 = 1$. + pub const ONE: Self = Self(Some(0)); + /// The zero coefficient. + pub const ZERO: Self = Self(None); + + /// The zero coefficient (method form, for symmetry with [`Tau::one`]). + pub fn zero() -> Self { + Self::ZERO + } + + /// The unit $1 = \tau^0$. + pub fn one() -> Self { + Self::ONE + } + + /// The coefficient $\tau^k$. + pub fn power(k: u32) -> Self { + Self(Some(k)) + } + + /// Whether this coefficient is zero. + pub fn is_zero(self) -> bool { + self.0.is_none() + } + + /// The $\tau$-valuation: `Some(k)` if this is $\tau^k$, `None` if zero. + pub fn valuation(self) -> Option { + self.0 + } + + /// The $n$-th power $(\tau^k)^n = \tau^{kn}$, with the convention $x^0 = 1$. + pub fn pow(self, n: u32) -> Self { + if n == 0 { + Self::ONE + } else { + Self(self.0.map(|k| k * n)) + } + } + + /// Multiply by $\tau^k$, i.e. raise the valuation by `k` (zero stays zero). + pub fn shift(self, k: u32) -> Self { + Self(self.0.map(|v| v + k)) + } +} + +impl std::ops::Mul for Tau { + type Output = Self; + + /// Product: $\tau^i \cdot \tau^j = \tau^{i+j}$; zero absorbs. + // Multiplying τ-powers genuinely *adds* their valuations, so the `+` here is correct. + #[allow(clippy::suspicious_arithmetic_impl)] + fn mul(self, other: Self) -> Self { + match (self.0, other.0) { + (Some(i), Some(j)) => Self(Some(i + j)), + _ => Self::ZERO, + } + } +} + +impl std::ops::Add for Tau { + type Output = Self; + + /// Sum of two homogeneous coefficients in one bidegree. By homogeneity the + /// two nonzero powers must be equal, so $\tau^k + \tau^k = 0$ over + /// $\mathbb{F}_2$; a zero summand is absorbed. Adding *distinct* nonzero + /// powers is inhomogeneous and cannot arise (checked in debug builds). + fn add(self, other: Self) -> Self { + match (self.0, other.0) { + (None, x) | (x, None) => Self(x), + (Some(i), Some(j)) => { + debug_assert_eq!(i, j, "adding inhomogeneous tau powers τ^{i} + τ^{j}"); + Self::ZERO + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_zero_one() { + assert!(Tau::zero().is_zero()); + assert!(!Tau::one().is_zero()); + assert_eq!(Tau::one().valuation(), Some(0)); + assert_eq!(Tau::power(3).valuation(), Some(3)); + assert_eq!(Tau::zero().valuation(), None); + } + + #[test] + fn test_mul_adds_valuations() { + assert_eq!(Tau::power(2) * Tau::power(3), Tau::power(5)); + assert_eq!(Tau::one() * Tau::power(4), Tau::power(4)); + assert_eq!(Tau::zero() * Tau::power(4), Tau::zero()); + assert_eq!(Tau::power(4) * Tau::zero(), Tau::zero()); + } + + #[test] + fn test_add_cancels_equal_powers() { + assert_eq!(Tau::power(2) + Tau::power(2), Tau::zero()); + assert_eq!(Tau::power(2) + Tau::zero(), Tau::power(2)); + assert_eq!(Tau::zero() + Tau::power(5), Tau::power(5)); + assert_eq!(Tau::zero() + Tau::zero(), Tau::zero()); + } + + #[test] + fn test_pow_and_shift() { + assert_eq!(Tau::power(3).pow(2), Tau::power(6)); + assert_eq!(Tau::power(3).pow(0), Tau::one()); + assert_eq!(Tau::zero().pow(0), Tau::one()); + assert_eq!(Tau::zero().pow(3), Tau::zero()); + assert_eq!(Tau::power(2).shift(3), Tau::power(5)); + assert_eq!(Tau::zero().shift(3), Tau::zero()); + } +} From a6f269d6448af3afa2a96ea0620be490b0d16a97 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Thu, 27 Aug 2026 21:35:31 -0400 Subject: [PATCH 02/15] motivic: add a criterion bench for the C-motivic engine MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three groups, from the kernel outwards: a single `multiply_closed`, a whole `fill_block` (the batch unit a resolution asks for, with throughput in structure constants), and `enum_basis`. The GPU handoff note calls `multiply_closed` the arithmetic bottleneck of the deformation pipeline, and the review raises the cost of the `BTreeMap` behind `DualElement`; neither had a number attached. This is the measurement both need, and it replaces the ad-hoc `PRODUCT_NANOS` counter that used to stand in for it. Baseline on this machine (mean): motivic_product/xi_small 1.40 µs motivic_product/xi_medium 11.44 µs motivic_product/xi_large 417.44 µs motivic_product/q_small 0.68 µs motivic_product/q_medium 2.40 µs motivic_product/q_large 15.41 µs motivic_block/12 348.69 µs motivic_block/16 2.09 ms motivic_block/20 10.36 ms motivic_basis/20 2.04 µs motivic_basis/30 7.13 µs motivic_basis/40 14.87 µs Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01F25ZVbsP7ULg41iY3MP6FX --- ext/crates/algebra/Cargo.toml | 4 + ext/crates/algebra/benches/motivic.rs | 101 ++++++++++++++++++++++++++ 2 files changed, 105 insertions(+) create mode 100644 ext/crates/algebra/benches/motivic.rs diff --git a/ext/crates/algebra/Cargo.toml b/ext/crates/algebra/Cargo.toml index 3b98825cb5..4dbf1fc02a 100644 --- a/ext/crates/algebra/Cargo.toml +++ b/ext/crates/algebra/Cargo.toml @@ -55,3 +55,7 @@ harness = false [[bench]] name = "nassau_milnor" harness = false + +[[bench]] +name = "motivic" +harness = false diff --git a/ext/crates/algebra/benches/motivic.rs b/ext/crates/algebra/benches/motivic.rs new file mode 100644 index 0000000000..45d496e9dc --- /dev/null +++ b/ext/crates/algebra/benches/motivic.rs @@ -0,0 +1,101 @@ +//! Benchmarks for the C-motivic Steenrod algebra engine. +//! +//! Three levels, from the kernel outwards: +//! +//! - `motivic_product` — a single [`multiply_closed`], the Kong–Lin Theorem 5.1 product. This is +//! the arithmetic the whole layer is built on. +//! - `motivic_block` — [`MotivicMilnorAlgebra::fill_block`], the batch unit a resolution actually +//! asks for: every structure constant for one pair of topological degrees. Throughput is in +//! structure constants, so the numbers are comparable across degrees. +//! - `motivic_basis` — [`enum_basis`], the basis enumeration each new degree pays once. +//! +//! The `motivic_block` group is the one to watch when changing the coefficient representation: +//! it is the only group that exercises the `DualElement` map, the index lookup, and the product +//! together, in the proportion a resolution hits them. + +use algebra::{ + MotivicMilnorAlgebra, + motivic::milnor::{Monomial, enum_basis, multiply_closed}, +}; +use criterion::{ + BenchmarkGroup, BenchmarkId, Criterion, Throughput, criterion_group, criterion_main, + measurement::WallTime, +}; +use pprof::criterion::{Output, PProfProfiler}; + +/// `Q(E)P(R)` from the `Q` indices and the ξ exponents, in the paper's indexing where `R[0]` +/// belongs to ξ_0 = 1 and is skipped. +fn elt(q: &[u32], xi: &[u32]) -> Monomial { + let mut r = vec![0]; + r.extend_from_slice(xi); + (q.iter().map(|i| 1 << i).sum(), r) +} + +fn bench_product(g: &mut BenchmarkGroup, name: &str, a: Monomial, b: Monomial) { + g.bench_function(name, |bench| { + bench.iter(|| std::hint::black_box(multiply_closed(&a, &b))); + }); +} + +fn product(c: &mut Criterion) { + let mut g = c.benchmark_group("motivic_product"); + + // Pure ξ: the classical Milnor-matrix part of the formula, with no Y enumeration. + bench_product(&mut g, "xi/small", elt(&[], &[2]), elt(&[], &[2])); + bench_product(&mut g, "xi/medium", elt(&[], &[4, 1]), elt(&[], &[2, 1])); + bench_product(&mut g, "xi/large", elt(&[], &[6, 2, 1]), elt(&[], &[4, 1])); + + // With a Q-part, which is what turns on the second (`Y`) matrix and the τ-rewriting. + bench_product(&mut g, "q/small", elt(&[0], &[1]), elt(&[1], &[1])); + bench_product(&mut g, "q/medium", elt(&[0, 1], &[2]), elt(&[2], &[1, 1])); + bench_product(&mut g, "q/large", elt(&[0, 2], &[3, 1]), elt(&[1], &[2, 1])); + + g.finish(); +} + +fn block(c: &mut Criterion) { + let mut g = c.benchmark_group("motivic_block"); + + for t in [12, 16, 20] { + // Size the throughput by the number of structure constants in the block, so the + // per-product cost is comparable across degrees. + let dims = { + let alg = MotivicMilnorAlgebra::new(); + alg.compute_basis(t); + alg.dimension(t) + }; + g.throughput(Throughput::Elements((dims * dims) as u64)); + g.bench_with_input(BenchmarkId::from_parameter(t), &t, |bench, &t| { + // A fresh algebra per iteration: `fill_block` is memoized, so reusing one would + // measure the cache rather than the product. + bench.iter_batched( + MotivicMilnorAlgebra::new, + |alg| alg.fill_block(t, t), + criterion::BatchSize::SmallInput, + ); + }); + } + + g.finish(); +} + +fn basis(c: &mut Criterion) { + let mut g = c.benchmark_group("motivic_basis"); + + for t in [20, 30, 40] { + g.bench_with_input(BenchmarkId::from_parameter(t), &t, |bench, &t| { + bench.iter(|| std::hint::black_box(enum_basis(t))); + }); + } + + g.finish(); +} + +criterion_group! { + name = benches; + config = Criterion::default() + .measurement_time(std::time::Duration::from_secs(3)) + .with_profiler(PProfProfiler::new(100, Output::Flamegraph(None))); + targets = product, block, basis +} +criterion_main!(benches); From 4588992c68bffea3338bd04fd13e0cb2beb99116 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Thu, 27 Aug 2026 21:40:29 -0400 Subject: [PATCH 03/15] motivic: rework the monomial and coefficient representation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Four steps on one arc, two of them answering review comments. Named fields first — "Could we switch from `(u32, Vec)` to a struct with named fields? Would be clearer what they mean." `Monomial { q_part, p_part }`, field names matching `MilnorBasisElement`'s. The derived `Ord` is lexicographic in that field order, exactly the tuple ordering the per-degree bases were already sorted and binary-searched by, so indexing is unchanged. Then the xi exponents become the classical bit-packed `PPart`, making `Monomial` `Copy` and 12 bytes with no heap. Kong–Lin index from ξ₀ = 1, so their `R[0]` is identically zero for every monomial and carries no information; dropping it is what lets the motivic exponent sequence *be* a classical one, with `from_paper`/`paper_p_part` converting at the two boundaries where the paper's indexing is the natural one. Then linear combinations, twice reviewed — "this `BTreeMap` is probably really bad for performance" and "should be consistent with how we handle linear combinations of milnor basis elements in the normal steenrod algebra". `SparseSum` is a `Vec<(K, Tau)>` kept sorted by key, so equality stays canonical and lookup stays logarithmic. Worth recording that its measured effect was a wash — not what the review, or I, expected. Finally the allocation traffic, from a pprof profile of `motivic_block/20` rather than a guess: a quarter of it was walking `Vec>`. The candidate column lists become `Columns` — every column end to end in one allocation with an offset table, iterated as `&[u32]` — and `NB` drops from 64 to 32. Co-Authored-By: Claude Opus 4.8 Co-Authored-By: Claude Opus 5 --- ext/crates/algebra/benches/motivic.rs | 4 +- .../algebra/src/algebra/milnor_algebra.rs | 8 +- .../algebra/src/algebra/motivic/milnor.rs | 681 +++++++++++------- 3 files changed, 438 insertions(+), 255 deletions(-) diff --git a/ext/crates/algebra/benches/motivic.rs b/ext/crates/algebra/benches/motivic.rs index 45d496e9dc..50a39a8430 100644 --- a/ext/crates/algebra/benches/motivic.rs +++ b/ext/crates/algebra/benches/motivic.rs @@ -28,12 +28,12 @@ use pprof::criterion::{Output, PProfProfiler}; fn elt(q: &[u32], xi: &[u32]) -> Monomial { let mut r = vec![0]; r.extend_from_slice(xi); - (q.iter().map(|i| 1 << i).sum(), r) + Monomial::from_paper(q.iter().map(|i| 1 << i).sum(), &r).unwrap() } fn bench_product(g: &mut BenchmarkGroup, name: &str, a: Monomial, b: Monomial) { g.bench_function(name, |bench| { - bench.iter(|| std::hint::black_box(multiply_closed(&a, &b))); + bench.iter(|| std::hint::black_box(multiply_closed(a, b))); }); } diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index afbe59bf8c..d62ec0793b 100644 --- a/ext/crates/algebra/src/algebra/milnor_algebra.rs +++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs @@ -122,7 +122,13 @@ pub type PPartEntry = u32; /// [`Self::MAX_DEGREE`], which [`MilnorAlgebra::compute_basis`] enforces up front, so the packing /// can never silently truncate. [`Self::set`] asserts it anyway, and [`Self::try_from_slice`] /// reports failure instead of panicking for input that has not been through that gate. -#[derive(Clone, Copy, PartialEq, Eq, Hash, Default)] +/// +/// The derived ordering compares the packed words. That is a total order and consistent with +/// equality, which is all a sorted table needs, but it is *not* lexicographic in the exponents: +/// entry `i` sits at [`Self::shift`]`(i)`, so `r_1` is the least significant field and therefore +/// the last tie-breaker. Callers that need the exponents ordered lexicographically must compare +/// [`Self::iter`] instead. +#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] pub struct PPart(u64); impl PPart { diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index c901ff5b17..ce271df9cf 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -41,18 +41,18 @@ use std::{ cell::RefCell, - collections::BTreeMap, rc::Rc, sync::{Arc, OnceLock, RwLock}, }; -use fp::prime::Binomial; +use fp::prime::{Binomial, iter::BitflagIterator}; use itertools::Itertools; use maybe_rayon::prelude::*; use once::OnceVec; use rustc_hash::FxHashMap; use super::Tau; +use crate::algebra::milnor_algebra::PPart; /// Drop the trailing zeros from an exponent vector. /// @@ -65,11 +65,77 @@ fn trim(mut v: Vec) -> Vec { v } -/// A basis monomial $\tau(E)\xi(R)$ of the dual algebra: the square-free exterior part `E` as a -/// bitmask (bit `i` is $\tau_i$), and the $\xi$ exponent vector `R` in the paper's indexing, -/// where `R[0]` belongs to $\xi_0 = 1$ and is therefore ignored. Trailing zeros are trimmed, -/// so equal monomials have equal representations and the tuple can be used as a map key. -pub type Monomial = (u32, Vec); +/// A basis monomial $\tau(E)\xi(R)$ of the dual algebra. +/// +/// The derived ordering is lexicographic in `(q_part, p_part)`, which is what +/// [`MotivicMilnorAlgebra`] sorts and binary-searches its per-degree bases by. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub struct Monomial { + /// The square-free exterior part $E$ as a bitmask: bit `i` is $\tau_i$. + pub q_part: u32, + /// The $\xi$ exponents, as the classical bit-packed [`PPart`]: entry `i` is the exponent of + /// $\xi_{i+1}$. + /// + /// The Kong–Lin formulas index $R$ from $\xi_0 = 1$, so their `R[0]` is always zero for a + /// monomial and carries no information. Dropping it is what lets the exponents share the + /// classical representation; [`Self::from_paper`] and [`Self::paper_p_part`] convert at the + /// boundaries where the paper's indexing is the natural one. + pub p_part: PPart, +} + +impl Monomial { + /// The monomial $\tau(E)\xi(R)$. + pub fn new(q_part: u32, p_part: PPart) -> Self { + Self { q_part, p_part } + } + + /// The monomial whose $\xi$ exponents are given in the paper's indexing, `r[j]` being the + /// exponent of $\xi_j$. `r[0]` must be zero, since $\xi_0 = 1$. + /// + /// Returns `None` if an exponent is too large to pack, which for the degrees this engine + /// reaches means the caller has gone out of range. + pub fn from_paper(q_part: u32, r: &[u32]) -> Option { + debug_assert_eq!( + r.first().copied().unwrap_or(0), + 0, + "xi_0 = 1 slot must be 0" + ); + Some(Self::new( + q_part, + PPart::try_from_slice(r.get(1..).unwrap_or_default())?, + )) + } + + /// The $\xi$ exponents in the paper's indexing, with the $\xi_0$ slot restored. + pub fn paper_p_part(&self) -> Vec { + if self.p_part.is_empty() { + return Vec::new(); + } + std::iter::once(0).chain(self.p_part.iter()).collect() + } + + /// The unit monomial $1$. + pub fn one() -> Self { + Self::default() + } + + /// The C-motivic (prime 2) bidegree `(t, w)`: `|τ_i| = (2^{i+1}-1, 2^i-1)` and + /// `|ξ_{i+1}| = (2^{i+2}-2, 2^{i+1}-1)`. + /// + /// [`MotivicMilnorAlgebra::bidegree`] negates `w` for the algebra's own weight convention. + pub fn bidegree(&self) -> (i32, i32) { + let (mut t, mut w) = (0i32, 0i32); + for i in BitflagIterator::set_bit_iterator(self.q_part as u64) { + t += (1 << (i + 1)) - 1; + w += (1 << i) - 1; + } + for (i, r) in self.p_part.iter().enumerate() { + t += r as i32 * ((1 << (i + 2)) - 2); + w += r as i32 * ((1 << (i + 1)) - 1); + } + (t, w) + } +} /// $\Sigma(R) = \sum_i r_i$. fn sigma(r: &[u32]) -> u32 { @@ -210,73 +276,150 @@ fn rewrite_tau_dfs( // The dual algebra A_** (used to compute the product in A_C by duality). // --------------------------------------------------------------------------- -/// An element of the dual algebra $A_{**}$: [`Monomial`]s mapped to their -/// $\mathbb{F}_2[\tau]$ coefficients. +/// A sparse $\mathbb{F}_2[\tau]$-linear combination of `K`s, as a `Vec` sorted by key. /// /// Everything we compute is homogeneous — $\psi$ and the products are graded — so every -/// coefficient is a single power of $\tau$: a [`Tau`] scalar. Zero coefficients are never stored, -/// so equality is canonical, and coefficient arithmetic is exactly [`Tau`]'s arithmetic — `mul` -/// adds valuations, `add` cancels equal powers mod 2 (unequal powers would be inhomogeneous and -/// cannot arise). -pub type DualElement = BTreeMap; - -/// Elementwise sum of two exponent vectors, trimmed of trailing zeros. -fn vec_add(a: &[u32], b: &[u32]) -> Vec { - let n = a.len().max(b.len()); - trim( - (0..n) - .map(|i| a.get(i).copied().unwrap_or(0) + b.get(i).copied().unwrap_or(0)) - .collect(), - ) +/// coefficient is a single power of $\tau$: a [`Tau`] scalar. Zero coefficients are never stored +/// and the terms are kept sorted, so equality is canonical, and coefficient arithmetic is exactly +/// [`Tau`]'s arithmetic — `mul` adds valuations, `add` cancels equal powers mod 2 (unequal powers +/// would be inhomogeneous and cannot arise). +/// +/// A sorted `Vec` rather than a `BTreeMap`: these combinations are small — a product of two basis +/// elements has a handful of terms — and the keys are now `Copy` words, so one contiguous +/// allocation searched by [`slice::binary_search_by`] beats a tree that allocates per node. +#[derive(Debug, Clone, Default, PartialEq, Eq)] +pub struct SparseSum { + terms: Vec<(K, Tau)>, } -/// Add `coeff * key` into a sparse $\mathbb{F}_2[\tau]$-linear combination `acc`, dropping the -/// entry if it cancels to zero. Used for both [`DualElement`] and [`TensorElement`]; the -/// coefficient bookkeeping is entirely [`Tau`] arithmetic. -fn add_term(acc: &mut BTreeMap, key: K, coeff: Tau) { - use std::collections::btree_map::Entry; - if coeff.is_zero() { - return; +impl SparseSum { + /// The zero combination. + pub fn new() -> Self { + Self { terms: Vec::new() } } - match acc.entry(key) { - Entry::Occupied(mut o) => { - let sum = *o.get() + coeff; - if sum.is_zero() { - o.remove(); - } else { - *o.get_mut() = sum; + + /// The number of non-zero terms. + pub fn len(&self) -> usize { + self.terms.len() + } + + /// Whether this is zero. + pub fn is_empty(&self) -> bool { + self.terms.is_empty() + } + + /// The coefficient of `key`, if non-zero. + pub fn get(&self, key: &K) -> Option<&Tau> { + self.position(key).ok().map(|i| &self.terms[i].1) + } + + /// Add `coeff * key`, dropping the term if it cancels to zero. + pub fn add_term(&mut self, key: K, coeff: Tau) { + if coeff.is_zero() { + return; + } + match self.position(&key) { + Ok(i) => { + let sum = self.terms[i].1 + coeff; + if sum.is_zero() { + self.terms.remove(i); + } else { + self.terms[i].1 = sum; + } } + Err(i) => self.terms.insert(i, (key, coeff)), } - Entry::Vacant(v) => { - v.insert(coeff); + } + + /// The index of `key`, or where it would be inserted. + fn position(&self, key: &K) -> Result { + self.terms.binary_search_by(|(k, _)| k.cmp(key)) + } + + /// The terms, in key order. + pub fn iter(&self) -> std::slice::Iter<'_, (K, Tau)> { + self.terms.iter() + } +} + +impl FromIterator<(K, Tau)> for SparseSum { + fn from_iter>(iter: I) -> Self { + let mut out = Self::new(); + for (key, coeff) in iter { + out.add_term(key, coeff); } + out + } +} + +impl From<[(K, Tau); N]> for SparseSum { + fn from(terms: [(K, Tau); N]) -> Self { + terms.into_iter().collect() + } +} + +impl IntoIterator for SparseSum { + type IntoIter = std::vec::IntoIter<(K, Tau)>; + type Item = (K, Tau); + + fn into_iter(self) -> Self::IntoIter { + self.terms.into_iter() } } +impl<'a, K> IntoIterator for &'a SparseSum { + type IntoIter = std::slice::Iter<'a, (K, Tau)>; + type Item = &'a (K, Tau); + + fn into_iter(self) -> Self::IntoIter { + self.terms.iter() + } +} + +/// An element of the dual algebra $A_{**}$: [`Monomial`]s with $\mathbb{F}_2[\tau]$ +/// coefficients. +pub type DualElement = SparseSum; + +/// Elementwise sum of two exponent sequences, i.e. the product of the two $\xi$ monomials. +fn ppart_add(a: PPart, b: PPart) -> PPart { + let mut out = PPart::zero(); + for i in 0..a.len().max(b.len()) { + out.set(i, a.get(i) + b.get(i)); + } + out +} + /// Multiply two basis monomials of $A_{**}$ and accumulate `coeff * (m1 * m2)` into `acc`. /// /// $A_{**}$ is commutative, so $\tau(E_1)\xi(R_1)\cdot\tau(E_2)\xi(R_2)$ is obtained by adding /// exponents: the exterior parts form $S = E_1 + E_2$ (entries in $\{0,1,2\}$), which is /// rewritten into the square-free basis via [`rewrite_tau`] ($\tau_i^2 = \tau\xi_{i+1}$), and /// the resulting $\xi$ exponents are added to $R_1 + R_2$. -fn mul_monomials(m1: &Monomial, m2: &Monomial, coeff: Tau, acc: &mut DualElement) { - let (e1, e2) = (m1.0, m2.0); +fn mul_monomials(m1: Monomial, m2: Monomial, coeff: Tau, acc: &mut DualElement) { + let (e1, e2) = (m1.q_part, m2.q_part); let bits = u32::BITS - (e1 | e2).leading_zeros(); let s: Vec = (0..bits) .map(|i| ((e1 >> i) & 1) + ((e2 >> i) & 1)) .collect(); - let r12 = vec_add(&m1.1, &m2.1); + let r12 = ppart_add(m1.p_part, m2.p_part); for term in rewrite_tau(&s) { - let r = vec_add(&term.r, &r12); - add_term(acc, (term.e_mask, r), coeff * Tau::power(term.tau_pow)); + // `term.r` is paper-indexed and `term.r[0]` is zero (`c_coeff` vanishes otherwise); it + // is also trimmed, so it can be empty. + let term_r = PPart::try_from_slice(term.r.get(1..).unwrap_or_default()) + .expect("rewrite_tau exponent out of packing range"); + let r = ppart_add(term_r, r12); + acc.add_term( + Monomial::new(term.e_mask, r), + coeff * Tau::power(term.tau_pow), + ); } } /// The product of two elements of $A_{**}$. pub fn dual_mul(a: &DualElement, b: &DualElement) -> DualElement { let mut out = DualElement::new(); - for (m1, &c1) in a { - for (m2, &c2) in b { + for &(m1, c1) in a { + for &(m2, c2) in b { mul_monomials(m1, m2, c1 * c2, &mut out); } } @@ -285,19 +428,17 @@ pub fn dual_mul(a: &DualElement, b: &DualElement) -> DualElement { /// The generator $\tau_i \in A_{**}$. pub fn tau_gen(i: usize) -> DualElement { - DualElement::from([((1u32 << i, vec![]), Tau::one())]) + DualElement::from([(Monomial::new(1 << i, PPart::zero()), Tau::one())]) } /// The generator $\xi_i \in A_{**}$ (for $i \ge 1$). pub fn xi_gen(i: usize) -> DualElement { - let mut r = vec![0u32; i + 1]; - r[i] = 1; - DualElement::from([((0, r), Tau::one())]) + DualElement::from([(xi_pow_mon(i, 1), Tau::one())]) } /// The unit $1 \in A_{**}$. pub fn dual_one() -> DualElement { - DualElement::from([((0, vec![]), Tau::one())]) + DualElement::from([(Monomial::one(), Tau::one())]) } // --------------------------------------------------------------------------- @@ -306,26 +447,26 @@ pub fn dual_one() -> DualElement { /// An element of $A_{**} \otimes A_{**}$: pairs of [`Monomial`]s → $\mathbb{F}_2[\tau]$ /// coefficients. Zero coefficients are never stored. -pub type TensorElement = BTreeMap<(Monomial, Monomial), Tau>; +pub type TensorElement = SparseSum<(Monomial, Monomial)>; /// The unit $1 \otimes 1$. fn tensor_one() -> TensorElement { - TensorElement::from([(((0, vec![]), (0, vec![])), Tau::one())]) + TensorElement::from([((Monomial::one(), Monomial::one()), Tau::one())]) } /// Multiply in $A_{**} \otimes A_{**}$: $(a_L \otimes a_R)(b_L \otimes b_R) = (a_L b_L) \otimes (a_R b_R)$. fn tensor_mul(t1: &TensorElement, t2: &TensorElement) -> TensorElement { let mut out = TensorElement::new(); - for ((al, ar), &c1) in t1 { - for ((bl, br), &c2) in t2 { + for &((al, ar), c1) in t1 { + for &((bl, br), c2) in t2 { let mut left = DualElement::new(); mul_monomials(al, bl, Tau::one(), &mut left); let mut right = DualElement::new(); mul_monomials(ar, br, Tau::one(), &mut right); let c = c1 * c2; - for (ml, &cl) in &left { - for (mr, &cr) in &right { - add_term(&mut out, (ml.clone(), mr.clone()), c * (cl * cr)); + for &(ml, cl) in &left { + for &(mr, cr) in &right { + out.add_term((ml, mr), c * (cl * cr)); } } } @@ -335,23 +476,26 @@ fn tensor_mul(t1: &TensorElement, t2: &TensorElement) -> TensorElement { /// The monomial $\xi_j^p$ (with $\xi_0 = 1$). fn xi_pow_mon(j: usize, p: u32) -> Monomial { - if j == 0 || p == 0 { - (0, vec![]) - } else { - let mut r = vec![0u32; j + 1]; - r[j] = p; - (0, r) + let mut r = PPart::zero(); + if j > 0 && p > 0 { + r.set(j - 1, p); } + Monomial::new(0, r) } /// $\psi(\tau_k) = 1 \otimes \tau_k + \sum_{i=0}^{k} \tau_i \otimes \xi_{k-i}^{2^i}$. fn coprod_tau(k: usize) -> TensorElement { let mut out = TensorElement::new(); - add_term(&mut out, ((0, vec![]), (1 << k, vec![])), Tau::one()); + out.add_term( + (Monomial::one(), Monomial::new(1 << k, PPart::zero())), + Tau::one(), + ); for i in 0..=k { - add_term( - &mut out, - ((1 << i, vec![]), xi_pow_mon(k - i, 1 << i)), + out.add_term( + ( + Monomial::new(1 << i, PPart::zero()), + xi_pow_mon(k - i, 1 << i), + ), Tau::one(), ); } @@ -362,27 +506,21 @@ fn coprod_tau(k: usize) -> TensorElement { fn coprod_xi(k: usize) -> TensorElement { let mut out = TensorElement::new(); for i in 0..=k { - add_term( - &mut out, - (xi_pow_mon(i, 1), xi_pow_mon(k - i, 1 << i)), - Tau::one(), - ); + out.add_term((xi_pow_mon(i, 1), xi_pow_mon(k - i, 1 << i)), Tau::one()); } out } /// The coproduct of a basis monomial $\tau(E)\xi(R)$, computed as an algebra map: the product /// of the coproducts of its generators. -fn coproduct_monomial(e_mask: u32, r: &[u32]) -> TensorElement { +fn coproduct_monomial(m: &Monomial) -> TensorElement { let mut acc = tensor_one(); - for i in 0..u32::BITS { - if (e_mask >> i) & 1 != 0 { - acc = tensor_mul(&acc, &coprod_tau(i as usize)); - } + for i in BitflagIterator::set_bit_iterator(m.q_part as u64) { + acc = tensor_mul(&acc, &coprod_tau(i)); } - for (j, &rj) in r.iter().enumerate() { - for _ in 0..rj { - acc = tensor_mul(&acc, &coprod_xi(j)); + for (i, r) in m.p_part.iter().enumerate() { + for _ in 0..r { + acc = tensor_mul(&acc, &coprod_xi(i + 1)); } } acc @@ -391,9 +529,9 @@ fn coproduct_monomial(e_mask: u32, r: &[u32]) -> TensorElement { /// The coproduct of an arbitrary element of $A_{**}$ (extended $\mathbb{F}_2[\tau]$-linearly). pub fn coproduct(elt: &DualElement) -> TensorElement { let mut out = TensorElement::new(); - for ((e, r), &c) in elt { - for (key, cc) in coproduct_monomial(*e, r) { - add_term(&mut out, key, c * cc); + for &(m, c) in elt { + for (key, cc) in coproduct_monomial(&m) { + out.add_term(key, c * cc); } } out @@ -418,7 +556,7 @@ fn chi_xi(k: usize) -> DualElement { for i in 0..k { let term = dual_mul(&chi_xi(i), &xi_pow_elt(k - i, 1 << i)); for (mon, c) in term { - add_term(&mut acc, mon, c); + acc.add_term(mon, c); } } acc @@ -431,7 +569,7 @@ fn chi_tau(k: usize) -> DualElement { for i in 0..k { let term = dual_mul(&chi_tau(i), &xi_pow_elt(k - i, 1 << i)); for (mon, c) in term { - add_term(&mut acc, mon, c); + acc.add_term(mon, c); } } acc @@ -443,20 +581,20 @@ fn chi_tau(k: usize) -> DualElement { /// involution. pub fn antipode(elt: &DualElement) -> DualElement { let mut out = DualElement::new(); - for ((e, r), &c) in elt { - let mut m = dual_one(); + for &(mon, c) in elt { + let mut acc = dual_one(); for i in 0..u32::BITS { - if (e >> i) & 1 != 0 { - m = dual_mul(&m, &chi_tau(i as usize)); + if (mon.q_part >> i) & 1 != 0 { + acc = dual_mul(&acc, &chi_tau(i as usize)); } } - for (j, &rj) in r.iter().enumerate() { - for _ in 0..rj { - m = dual_mul(&m, &chi_xi(j)); + for (i, r) in mon.p_part.iter().enumerate() { + for _ in 0..r { + acc = dual_mul(&acc, &chi_xi(i + 1)); } } - for (mon, cc) in m { - add_term(&mut out, mon, c * cc); + for (out_mon, cc) in acc { + out.add_term(out_mon, c * cc); } } out @@ -466,24 +604,6 @@ pub fn antipode(elt: &DualElement) -> DualElement { // The product in the Steenrod algebra A_C, computed by dualizing ψ. // --------------------------------------------------------------------------- -/// The C-motivic (prime 2) bidegree `(t, w)` of the basis monomial `τ(E)ξ(R)`, in the -/// *paper's* R-convention (`r[j]` is the exponent of `ξ_j`, with `ξ_0 = 1` so `r[0]` is -/// ignored). `|τ_i| = (2^{i+1}-1, 2^i-1)`, `|ξ_j| = (2^{j+1}-2, 2^j-1)`. -fn paper_bidegree(e_mask: u32, r: &[u32]) -> (i32, i32) { - let (mut t, mut w) = (0i32, 0i32); - for i in 0..u32::BITS { - if (e_mask >> i) & 1 != 0 { - t += (1 << (i + 1)) - 1; - w += (1 << i) - 1; - } - } - for (j, &rj) in r.iter().enumerate().skip(1) { - t += rj as i32 * ((1 << (j + 1)) - 2); - w += rj as i32 * ((1 << j) - 1); - } - (t, w) -} - /// All basis monomials `τ(E)ξ(R)` of a given topological degree `target`. pub fn enum_basis(target: i32) -> Vec { // Generators available up to `target`: τ_i of degree 2^{i+1}-1, ξ_j (j≥1) of degree @@ -523,7 +643,8 @@ fn enum_basis_dfs( ) { if idx == gens.len() { if rem == 0 { - out.push((*e, trim(r.to_vec()))); + // `r` is paper-indexed: `r[0]` is the (always zero) xi_0 slot. + out.push(Monomial::from_paper(*e, r).expect("basis exponent out of packing range")); } return; } @@ -554,15 +675,15 @@ fn enum_basis_dfs( /// Computed by duality: the coefficient of `z` in `a · b` is the coefficient of /// `mon(a) ⊗ mon(b)` in `ψ(mon(z))`, summed over the basis `z` of the appropriate /// topological degree. -pub fn multiply(a: &Monomial, b: &Monomial) -> DualElement { - let key = (a.clone(), b.clone()); - let t = paper_bidegree(a.0, &a.1).0 + paper_bidegree(b.0, &b.1).0; +pub fn multiply(a: Monomial, b: Monomial) -> DualElement { + let key = (a, b); + let t = a.bidegree().0 + b.bidegree().0; let mut out = DualElement::new(); for z in enum_basis(t) { - if let Some(&c) = coproduct_monomial(z.0, &z.1).get(&key) + if let Some(&c) = coproduct_monomial(&z).get(&key) && !c.is_zero() { - out.insert(z, c); + out.add_term(z, c); } } out @@ -588,6 +709,44 @@ pub fn multiply(a: &Monomial, b: &Monomial) -> DualElement { // tightens from S(X) ≤ R₁ to S(X) = R₁ — the classical admissible matrix. What remains // over the classical product is the second matrix Y, carrying the τᵢ/Q-part interaction. +/// A list of candidate columns, flattened into one allocation. +/// +/// The recursion walks these lists constantly, so the columns are stored end to end with an +/// offset table rather than as a `Vec>`: iterating that chased a pointer per column and +/// showed up as a fifth of the product's time. +#[derive(Debug, Default)] +struct Columns { + entries: Vec, + /// `offsets[k]..offsets[k + 1]` is column `k`. Always starts with a 0. + offsets: Vec, +} + +impl Columns { + fn push(&mut self, col: &[u32]) { + if self.offsets.is_empty() { + self.offsets.push(0); + } + self.entries.extend_from_slice(col); + self.offsets.push(self.entries.len() as u32); + } + + fn iter(&self) -> impl Iterator { + self.offsets + .windows(2) + .map(|w| &self.entries[w[0] as usize..w[1] as usize]) + } +} + +impl FromIterator> for Columns { + fn from_iter>>(iter: I) -> Self { + let mut out = Self::default(); + for col in iter { + out.push(&col); + } + out + } +} + /// All columns `(v₀, v₁, …)` of non-negative integers with `Σᵢ 2ⁱ vᵢ == target`. fn columns_eq(target: u32) -> Vec> { if target == 0 { @@ -606,43 +765,39 @@ fn columns_eq(target: u32) -> Vec> { out } -/// All columns with `Σᵢ 2ⁱ vᵢ ≤ bound` (distinct weighted sums, so no duplicates). -fn columns_le(bound: u32) -> Vec> { - (0..=bound).flat_map(columns_eq).collect() -} - thread_local! { /// Per-thread memo of [`columns_eq`] / [`columns_le`], keyed by target/bound. /// The candidate lists are pure functions of one integer and are requested /// over and over across products; caching them avoids re-enumerating and /// re-allocating. (`Rc` is confined to a single product's call stack.) - static COLS_EQ: RefCell>>>> = RefCell::new(FxHashMap::default()); - static COLS_LE: RefCell>>>> = RefCell::new(FxHashMap::default()); + static COLS_EQ: RefCell>> = RefCell::new(FxHashMap::default()); + static COLS_LE: RefCell>> = RefCell::new(FxHashMap::default()); } -fn columns_eq_cached(target: u32) -> Rc>> { +fn columns_eq_cached(target: u32) -> Rc { COLS_EQ.with(|c| { Rc::clone( c.borrow_mut() .entry(target) - .or_insert_with(|| Rc::new(columns_eq(target))), + .or_insert_with(|| Rc::new(columns_eq(target).into_iter().collect())), ) }) } -fn columns_le_cached(bound: u32) -> Rc>> { +/// All columns with `Σᵢ 2ⁱ vᵢ ≤ bound` (distinct weighted sums, so no duplicates). +fn columns_le_cached(bound: u32) -> Rc { COLS_LE.with(|c| { Rc::clone( c.borrow_mut() .entry(bound) - .or_insert_with(|| Rc::new(columns_le(bound))), + .or_insert_with(|| Rc::new((0..=bound).flat_map(columns_eq).collect())), ) }) } /// The column-0 options for the `X` matrix: `x_{0,0} = 0`, and `x_{i,0} ≤ R₁[i]` /// for `i ≥ 1` (bounded because it feeds the row sum `S(X)_i ≤ R₁[i]`). -fn col0_x_options(r1: &[u32]) -> Vec> { +fn col0_x_options(r1: &[u32]) -> Columns { let mut result = vec![vec![0u32]]; for &bound in r1.iter().skip(1) { let mut next = Vec::new(); @@ -658,10 +813,13 @@ fn col0_x_options(r1: &[u32]) -> Vec> { result.into_iter().map(trim).collect() } -/// Buffer size for row / column / anti-diagonal indices. All are bounded by -/// ~2·log₂(degree) + (number of ξ generators), far below this for any feasible -/// stem (a ξ or τ at index `i` already has degree ≥ 2^i). -const NB: usize = 64; +/// Buffer size for row / column / anti-diagonal indices. +/// +/// Anti-diagonals are indexed by `i + j`, with the row index `i` bounded by a column's length +/// (at most `log2(MAX_DEGREE) + 1`) and the column index `j` by [`PPart::MAX_LEN`] + 1. That +/// caps the sum near 21, and this is the next power of two; [`Acc`] holds three of these, so the +/// size is paid on every matrix. +const NB: usize = 32; /// Immutable context for the closed-form product recursion (Theorem 5.1, ρ = 0). struct Closed<'a> { @@ -672,7 +830,7 @@ struct Closed<'a> { e2_mask: u32, sigma2_e1: i64, /// Per-column candidate columns for the `X` matrix. - x_cands: &'a [Rc>>], + x_cands: &'a [Rc], } /// Incrementally maintained matrix state: row sums `S`, a running XOR per @@ -714,6 +872,7 @@ impl<'a> Closed<'a> { if !fits { continue; } + debug_assert!(cand.len() + j <= NB, "anti-diagonal index out of range"); for (i, &v) in cand.iter().enumerate() { acc.rows[i] += v; acc.or[i + j] ^= v; @@ -733,43 +892,44 @@ impl<'a> Closed<'a> { /// A complete `X` (already `b = 1`): form `S′`, the `Y` column targets, and the /// output P-part `T(X) = acc.sum`, then enumerate `Y`. fn on_x(&self, xcols: &[&'a [u32]], acc: &Acc, out: &mut DualElement) { - // S′ = R₁ − S(X) on indices ≥ 1 (S′[0] = 0); τ-power = Σ(S′). - let sprime: Vec = (0..self.l) - .map(|j| if j == 0 { 0 } else { self.r1v[j] - acc.rows[j] }) - .collect(); - let sigma2_sprime: i64 = sprime - .iter() - .enumerate() - .map(|(j, &v)| (v as i64) << j) - .sum(); - let sigma_sprime: u32 = sprime.iter().sum(); + // S′ = R₁ − S(X) on indices ≥ 1 (S′[0] = 0); τ-power = Σ(S′). This runs once per `X` + // matrix, so the working vectors are stack arrays rather than allocations. + let mut sprime = [0u32; NB]; + let (mut sigma2_sprime, mut sigma_sprime) = (0i64, 0u32); + for (j, sp) in sprime.iter_mut().enumerate().take(self.l).skip(1) { + let v = self.r1v[j] - acc.rows[j]; + *sp = v; + sigma2_sprime += (v as i64) << j; + sigma_sprime += v; + } + let sprime = &sprime[..self.l]; // RY targets: column j ≥ 1 of Y has weighted sum R₂[j] − R(X)[j]. - let mut ry: Vec = vec![0; self.l]; - for (j, ry_j) in ry.iter_mut().enumerate().skip(1) { + let mut ry = [0i64; NB]; + for j in 1..self.l { let rxj: u32 = xcols[j].iter().enumerate().map(|(i, &v)| v << i).sum(); - *ry_j = self.r2v[j] as i64 - rxj as i64; - if *ry_j < 0 { + ry[j] = self.r2v[j] as i64 - rxj as i64; + if ry[j] < 0 { return; } } // Σ₂(S(Y)) = Σ₂(E₁) + Σ₂(S′) is forced and equals Σ_{j≥0} R(Y)[j], so column // 0 of Y has a determined weighted sum. - let ry0 = (self.sigma2_e1 + sigma2_sprime) - ry[1..].iter().sum::(); + let ry0 = (self.sigma2_e1 + sigma2_sprime) - ry[1..self.l].iter().sum::(); if ry0 < 0 { return; } - // The output P-part T(X) (index 0 dropped: ξ₀ = 1). - let out_r = trim( - std::iter::once(0) - .chain((1..NB).map(|d| acc.sum[d])) - .collect(), - ); + // The output P-part T(X). The anti-diagonal sums are paper-indexed, so entry `i` of the + // packed exponent sequence is anti-diagonal `i + 1` (index 0 is the xi_0 slot). + let mut out_r = PPart::zero(); + for i in 0..PPart::MAX_LEN { + out_r.set(i, acc.sum[i + 1]); + } - let mut y_cands: Vec>>> = Vec::with_capacity(self.l); + let mut y_cands: Vec> = Vec::with_capacity(self.l); y_cands.push(columns_eq_cached(ry0 as u32)); - for &t in ry.iter().skip(1) { + for &t in &ry[1..self.l] { y_cands.push(columns_eq_cached(t as u32)); } @@ -777,9 +937,9 @@ impl<'a> Closed<'a> { y_cands: &y_cands, e1_mask: self.e1_mask, e2_mask: self.e2_mask, - sprime: &sprime, + sprime, sigma_sprime, - out_r: &out_r, + out_r, } .enum_y(0, &mut Acc::zero(), out); } @@ -788,12 +948,12 @@ impl<'a> Closed<'a> { /// Immutable context for the `Y` half of the recursion, the counterpart of [`Closed`]. `X` is /// already fixed by the time this is built, so `sprime` and `out_r` are constants here. struct ClosedY<'a> { - y_cands: &'a [Rc>>], + y_cands: &'a [Rc], e1_mask: u32, e2_mask: u32, sprime: &'a [u32], sigma_sprime: u32, - out_r: &'a [u32], + out_r: PPart, } impl ClosedY<'_> { @@ -855,9 +1015,8 @@ impl ClosedY<'_> { out_e_mask |= 1 << i; } } - add_term( - out, - (out_e_mask, self.out_r.to_vec()), + out.add_term( + Monomial::new(out_e_mask, self.out_r), Tau::power(self.sigma_sprime), ); } @@ -865,16 +1024,21 @@ impl ClosedY<'_> { /// The product `a · b` in `A_C` via Kong–Lin Theorem 5.1 (ρ = 0). Same contract /// as [`multiply`] (the duality oracle it is validated against). -pub fn multiply_closed(a: &Monomial, b: &Monomial) -> DualElement { - let (e1_mask, r1) = (a.0, &a.1); - let (e2_mask, r2) = (b.0, &b.1); - let l = r1.len().max(r2.len()).max(1); - let r1v: Vec = (0..l).map(|j| r1.get(j).copied().unwrap_or(0)).collect(); - let r2v: Vec = (0..l).map(|j| r2.get(j).copied().unwrap_or(0)).collect(); +pub fn multiply_closed(a: Monomial, b: Monomial) -> DualElement { + let (e1_mask, e2_mask) = (a.q_part, b.q_part); + // The recursion below follows Kong–Lin's indexing, in which `R[0]` is the xi_0 slot, so the + // exponents are unpacked into paper-indexed working vectors here. + let l = a.p_part.len().max(b.p_part.len()) + 1; + let paper = |m: Monomial| -> Vec { + std::iter::once(0) + .chain((0..l - 1).map(|i| m.p_part.get(i))) + .collect() + }; + let (r1v, r2v) = (paper(a), paper(b)); // Per-column candidate columns for X: column 0 bounded by R₁ rows, columns // j ≥ 1 by R₂[j] weighted. - let mut x_cands: Vec>>> = Vec::with_capacity(l); + let mut x_cands: Vec> = Vec::with_capacity(l); x_cands.push(Rc::new(col0_x_options(&r1v))); for &bound in r2v.iter().skip(1) { x_cands.push(columns_le_cached(bound)); @@ -988,8 +1152,7 @@ impl MotivicMilnorAlgebra { /// The `(topological degree, motivic weight)` of a basis element. pub fn bidegree(&self, degree: i32, idx: usize) -> (i32, i32) { - let (e, r) = &self.basis[degree as usize][idx]; - let (t, w) = paper_bidegree(*e, r); + let (t, w) = self.basis[degree as usize][idx].bidegree(); (t, -w) } @@ -1023,25 +1186,25 @@ impl MotivicMilnorAlgebra { .ok()?; } - let degree = paper_bidegree(e_mask, &r).0; + let m = Monomial::from_paper(e_mask, &r)?; + let degree = m.bidegree().0; if degree < 0 { return None; } self.compute_basis(degree); // The basis is stored trimmed, so trimming the parsed vector makes this an exact match. - self.index_of(degree, &(e_mask, trim(r))) - .map(|idx| (degree, idx)) + self.index_of(degree, &m).map(|idx| (degree, idx)) } /// A display string for a basis element (`Q_i … P(R)`). pub fn basis_element_to_string(&self, degree: i32, idx: usize) -> String { - let (e_mask, r) = self.basis_element(degree, idx); + let m = self.basis_element(degree, idx); let mut parts = Vec::new(); - for i in fp::prime::iter::BitflagIterator::set_bit_iterator(*e_mask as u64) { + for i in fp::prime::iter::BitflagIterator::set_bit_iterator(m.q_part as u64) { parts.push(format!("Q_{i}")); } - if !r.is_empty() { - parts.push(format!("P({})", r.iter().format(", "))); + if !m.p_part.is_empty() { + parts.push(format!("P({})", m.paper_p_part().iter().format(", "))); } if parts.is_empty() { "1".to_string() @@ -1081,7 +1244,7 @@ impl MotivicMilnorAlgebra { let a = &self.basis[t1 as usize][idx1]; let b = &self.basis[t2 as usize][idx2]; let t = t1 + t2; - multiply_closed(a, b) + multiply_closed(*a, *b) .into_iter() .map(|(z, c)| { ( @@ -1146,6 +1309,8 @@ impl std::fmt::Display for MotivicMilnorAlgebra { #[cfg(test)] mod tests { + use std::collections::BTreeMap; + use super::*; /// The classical mod-2 Milnor product of `P(a) * P(b)`, as a list of exponent sequences. @@ -1200,9 +1365,14 @@ mod tests { trim(r) } - /// The inverse of [`pp_to_paper`]: drop the ξ_0 slot. - fn paper_to_pp(r: &[u32]) -> Vec { - trim(r.get(1..).unwrap_or_default().to_vec()) + /// A monomial from paper-indexed exponents, panicking on anything out of packing range. + fn mon(q_part: u32, r: &[u32]) -> Monomial { + Monomial::from_paper(q_part, r).unwrap() + } + + /// The pure-ξ monomial with the given classical exponent sequence. + fn paper(pp: &[u32]) -> Monomial { + mon(0, &pp_to_paper(pp)) } #[test] @@ -1212,11 +1382,11 @@ mod tests { // Degree 0 is the unit; degrees 1 and 2 are 1-dimensional (Q_0 and P(ξ_1)). assert_eq!(alg.dimension(0), 1); - assert_eq!(alg.basis_element(0, 0), &(0u32, vec![])); + assert_eq!(alg.basis_element(0, 0), &mon(0u32, &[])); assert_eq!(alg.dimension(1), 1); - assert_eq!(alg.basis_element(1, 0), &(0b1u32, vec![])); // Q_0 + assert_eq!(alg.basis_element(1, 0), &mon(0b1u32, &[])); // Q_0 assert_eq!(alg.dimension(2), 1); - assert_eq!(alg.basis_element(2, 0), &(0u32, vec![0, 1])); // P(ξ_1) + assert_eq!(alg.basis_element(2, 0), &mon(0u32, &[0, 1])); // P(ξ_1) // bidegree: Q_0 is (1, 0), P(ξ_1) is (2, -1) in this presentation. assert_eq!(alg.bidegree(1, 0), (1, 0)); @@ -1226,13 +1396,13 @@ mod tests { let terms: DualElement = alg .product_indexed(1, 0, 2, 0) .into_iter() - .map(|(c, idx)| (alg.basis_element(3, idx).clone(), c)) + .map(|(c, idx)| (*alg.basis_element(3, idx), c)) .collect(); assert_eq!( terms, DualElement::from([ - ((0b10, vec![]), Tau::one()), - ((0b1, vec![0, 1]), Tau::one()) + (mon(0b10, &[]), Tau::one()), + (mon(0b1, &[0, 1]), Tau::one()) ]) ); } @@ -1245,15 +1415,15 @@ mod tests { for idx1 in 0..alg.dimension(t1) { for t2 in 0..=6 { for idx2 in 0..alg.dimension(t2) { - let a = alg.basis_element(t1, idx1).clone(); - let b = alg.basis_element(t2, idx2).clone(); + let a = *alg.basis_element(t1, idx1); + let b = *alg.basis_element(t2, idx2); // Indexed product agrees with the raw monomial product. let indexed: DualElement = alg .product_indexed(t1, idx1, t2, idx2) .into_iter() - .map(|(c, idx)| (alg.basis_element(t1 + t2, idx).clone(), c)) + .map(|(c, idx)| (*alg.basis_element(t1 + t2, idx), c)) .collect(); - assert_eq!(indexed, multiply(&a, &b)); + assert_eq!(indexed, multiply(a, b)); // Weight-homogeneous: w_out - (τ-power) = w_a + w_b. let (wa, wb) = (alg.bidegree(t1, idx1).1, alg.bidegree(t2, idx2).1); @@ -1272,22 +1442,22 @@ mod tests { // τ_0^2 = τ ξ_1 (the defining relation at ρ = 0). assert_eq!( dual_mul(&tau_gen(0), &tau_gen(0)), - DualElement::from([((0, vec![0, 1]), Tau::power(1))]) // τ^1 · ξ_1 + DualElement::from([(mon(0, &[0, 1]), Tau::power(1))]) // τ^1 · ξ_1 ); // ξ_1^2 is just the monomial ξ_1^2. assert_eq!( dual_mul(&xi_gen(1), &xi_gen(1)), - DualElement::from([((0, vec![0, 2]), Tau::one())]) + DualElement::from([(mon(0, &[0, 2]), Tau::one())]) ); // Distinct τ's commute and stay square-free: τ_0 τ_1. assert_eq!( dual_mul(&tau_gen(0), &tau_gen(1)), - DualElement::from([((0b11, vec![]), Tau::one())]) + DualElement::from([(mon(0b11, &[]), Tau::one())]) ); // τ_1^2 = τ ξ_2. assert_eq!( dual_mul(&tau_gen(1), &tau_gen(1)), - DualElement::from([((0, vec![0, 0, 1]), Tau::power(1))]) + DualElement::from([(mon(0, &[0, 0, 1]), Tau::power(1))]) ); // Multiplication by the unit is the identity. assert_eq!(dual_mul(&dual_one(), &tau_gen(2)), tau_gen(2)); @@ -1306,12 +1476,12 @@ mod tests { }; let motivic = |p1: &[u32], p2: &[u32]| -> Vec> { - let mut out: Vec> = multiply(&(0, pp_to_paper(p1)), &(0, pp_to_paper(p2))) + let mut out: Vec> = multiply(paper(p1), paper(p2)) .into_iter() - .filter(|((e, _), _)| *e == 0) // keep the pure-ξ outputs - .map(|((_, r), c)| { + .filter(|(m, _)| m.q_part == 0) // keep the pure-ξ outputs + .map(|(m, c)| { assert_eq!(c, Tau::one(), "pure-ξ output carried a τ power"); - paper_to_pp(&r) + m.p_part.iter().collect() }) .collect(); out.sort(); @@ -1351,22 +1521,19 @@ mod tests { // ones. Identity checked (pure-ξ, τ-free, so over 𝔽₂): expressing both inputs and // outputs of the motivic conjugate product in the standard basis via χ reproduces the // codebase product. `A[μ][W] = coeff of conj-mon μ in std ξ(W) = [χ(mon_W)]_μ`. - let mdeg = |pp: &[u32]| paper_bidegree(0, &pp_to_paper(pp)).0; + let mdeg = |pp: &[u32]| paper(pp).bidegree().0; let xi_pps = |t: i32| -> Vec> { enum_basis(t) .into_iter() - .filter(|(e, _)| *e == 0) - .map(|(_, r)| paper_to_pp(&r)) + .filter(|m| m.q_part == 0) + .map(|m| m.p_part.iter().collect()) .collect() }; let a_coeff = |mu_pp: &[u32], w_pp: &[u32]| -> u32 { - let ap = antipode(&DualElement::from([( - (0u32, pp_to_paper(w_pp)), - Tau::one(), - )])); + let ap = antipode(&DualElement::from([(paper(w_pp), Tau::one())])); // The antipode of a pure-ξ element is pure-ξ and τ-free, so a present monomial has // coefficient τ^0. - ap.get(&(0u32, pp_to_paper(mu_pp))) + ap.get(&paper(mu_pp)) .map_or(0, |&c| u32::from(c == Tau::one())) }; let cb_mul = |a_pp: &[u32], b_pp: &[u32], acc: &mut BTreeMap, u32>| { @@ -1386,10 +1553,10 @@ mod tests { (vec![2, 1], vec![2]), ] { // LHS: motivic conjugate product (E=0), converted to std via χ. - let conj: Vec> = multiply(&(0, pp_to_paper(&r1)), &(0, pp_to_paper(&r2))) + let conj: Vec> = multiply(paper(&r1), paper(&r2)) .into_iter() - .filter(|((e, _), _)| *e == 0) - .map(|((_, r), _)| paper_to_pp(&r)) + .filter(|(m, _)| m.q_part == 0) + .map(|(m, _)| m.p_part.iter().collect()) .collect(); let mut lhs: BTreeMap, u32> = BTreeMap::new(); for w in xi_pps(mdeg(&r1) + mdeg(&r2)) { @@ -1427,9 +1594,9 @@ mod tests { // coefficient sum + zero cancellation via Tau arithmetic). fn mul_elt_basis(x: &DualElement, c: &Monomial) -> DualElement { let mut out = DualElement::new(); - for (z, &cz) in x { - for (w, cw) in multiply(z, c) { - add_term(&mut out, w, cz * cw); + for &(z, cz) in x { + for (w, cw) in multiply(z, *c) { + out.add_term(w, cz * cw); } } out @@ -1438,15 +1605,15 @@ mod tests { let basis: Vec = (0..=5).flat_map(enum_basis).collect(); for a in &basis { for b in &basis { - let ab = multiply(a, b); + let ab = multiply(*a, *b); for c in &basis { let lhs = mul_elt_basis(&ab, c); - let bc = multiply(b, c); + let bc = multiply(*b, *c); // a · (b·c): multiply basis `a` on the left of each term of bc. let mut rhs = DualElement::new(); - for (z, &cz) in &bc { - for (w, cw) in multiply(a, z) { - add_term(&mut rhs, w, cz * cw); + for &(z, cz) in &bc { + for (w, cw) in multiply(*a, z) { + rhs.add_term(w, cz * cw); } } assert_eq!(lhs, rhs, "associativity failed at {a:?},{b:?},{c:?}"); @@ -1459,13 +1626,13 @@ mod tests { fn test_product_q0_p_xi1() { // The case my reading of Kong–Lin Theorem 5.1 got wrong. By duality: // Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1) (both with coefficient τ^0). - let q0 = (0b1, vec![]); - let p_xi1 = (0, vec![0, 1]); + let q0 = mon(0b1, &[]); + let p_xi1 = mon(0, &[0, 1]); assert_eq!( - multiply(&q0, &p_xi1), + multiply(q0, p_xi1), DualElement::from([ - ((0b10, vec![]), Tau::one()), // Q_1 - ((0b1, vec![0, 1]), Tau::one()), // Q_0 P(ξ_1) + (mon(0b10, &[]), Tau::one()), // Q_1 + (mon(0b1, &[0, 1]), Tau::one()), // Q_0 P(ξ_1) ]) ); } @@ -1473,16 +1640,20 @@ mod tests { #[test] fn test_product_unit_and_squares() { // 1 · x = x. - let x = (0b101, vec![0, 2]); + let x = mon(0b101, &[0, 2]); assert_eq!( - multiply(&(0, vec![]), &x), - DualElement::from([(x.clone(), Tau::one())]) + multiply(mon(0, &[]), x), + DualElement::from([(x, Tau::one())]) ); // Q_i^2 = 0. (Q_0 = Sq^1 is the motivic Bockstein; P(ξ_1) = Sq^2 does NOT square to // zero — its square is τ Q_0 Q_1 + …, a genuine motivic feature.) for i in 0..3 { assert!( - multiply(&(1 << i, vec![]), &(1 << i, vec![])).is_empty(), + multiply( + Monomial::new(1 << i, PPart::zero()), + Monomial::new(1 << i, PPart::zero()) + ) + .is_empty(), "Q_{i}^2 ≠ 0" ); } @@ -1499,11 +1670,11 @@ mod tests { let mut saw_tau = false; for a in &basis { for b in &basis { - let (wa, wb) = (paper_bidegree(a.0, &a.1).1, paper_bidegree(b.0, &b.1).1); - for (z, c) in multiply(a, b) { + let (wa, wb) = (a.bidegree().1, b.bidegree().1); + for (z, c) in multiply(*a, *b) { // Each coefficient is a single τ power by construction (a homogeneous Tau). let k = c.valuation().expect("stored coefficients are never zero") as i32; - let wz = paper_bidegree(z.0, &z.1).1; + let wz = z.bidegree().1; assert_eq!(wa + wb - wz, k, "weight mismatch: {a:?}*{b:?} → z={z:?}"); assert!(k >= 0); saw_tau |= k > 0; @@ -1523,8 +1694,8 @@ mod tests { assert_eq!( antipode(&xi_gen(2)), DualElement::from([ - ((0, vec![0, 0, 1]), Tau::one()), - ((0, vec![0, 3]), Tau::one()) + (mon(0, &[0, 0, 1]), Tau::one()), + (mon(0, &[0, 3]), Tau::one()) ]) ); assert_eq!(antipode(&tau_gen(0)), tau_gen(0)); @@ -1551,7 +1722,7 @@ mod tests { &DualElement::from([(r, Tau::one())]), ); for (mon, cc) in prod { - add_term(&mut out, mon, c * cc); + out.add_term(mon, c * cc); } } out @@ -1576,18 +1747,18 @@ mod tests { assert_eq!( coproduct(&tau_gen(1)), TensorElement::from([ - (((0, vec![]), (0b10, vec![])), Tau::one()), // 1 ⊗ τ_1 - (((0b1, vec![]), (0, vec![0, 1])), Tau::one()), // τ_0 ⊗ ξ_1 - (((0b10, vec![]), (0, vec![])), Tau::one()), // τ_1 ⊗ 1 + ((mon(0, &[]), mon(0b10, &[])), Tau::one()), // 1 ⊗ τ_1 + ((mon(0b1, &[]), mon(0, &[0, 1])), Tau::one()), // τ_0 ⊗ ξ_1 + ((mon(0b10, &[]), mon(0, &[])), Tau::one()), // τ_1 ⊗ 1 ]) ); // ψ(ξ_2) = 1⊗ξ_2 + ξ_1⊗ξ_1^2 + ξ_2⊗1 (matches Milnor). assert_eq!( coproduct(&xi_gen(2)), TensorElement::from([ - (((0, vec![]), (0, vec![0, 0, 1])), Tau::one()), // 1 ⊗ ξ_2 - (((0, vec![0, 1]), (0, vec![0, 2])), Tau::one()), // ξ_1 ⊗ ξ_1^2 - (((0, vec![0, 0, 1]), (0, vec![])), Tau::one()), // ξ_2 ⊗ 1 + ((mon(0, &[]), mon(0, &[0, 0, 1])), Tau::one()), // 1 ⊗ ξ_2 + ((mon(0, &[0, 1]), mon(0, &[0, 2])), Tau::one()), // ξ_1 ⊗ ξ_1^2 + ((mon(0, &[0, 0, 1]), mon(0, &[])), Tau::one()), // ξ_2 ⊗ 1 ]) ); } @@ -1658,18 +1829,24 @@ mod tests { fn test_closed_form_small_cases() { // The case a naive reading of Theorem 5.1 gets wrong (the ξ₀ = 1 index // absorption): Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1). - let q0 = (0b1, vec![]); - let p_xi1 = (0, vec![0, 1]); - assert_eq!(multiply_closed(&q0, &p_xi1), multiply(&q0, &p_xi1)); + let q0 = mon(0b1, &[]); + let p_xi1 = mon(0, &[0, 1]); + assert_eq!(multiply_closed(q0, p_xi1), multiply(q0, p_xi1)); // The τ-generating case: P(ξ_1)² has a τ Q_0 Q_1 term. - assert_eq!(multiply_closed(&p_xi1, &p_xi1), multiply(&p_xi1, &p_xi1)); + assert_eq!(multiply_closed(p_xi1, p_xi1), multiply(p_xi1, p_xi1)); // Q_i² = 0. for i in 0..3 { assert_eq!( - multiply_closed(&(1 << i, vec![]), &(1 << i, vec![])), - multiply(&(1 << i, vec![]), &(1 << i, vec![])) + multiply_closed( + Monomial::new(1 << i, PPart::zero()), + Monomial::new(1 << i, PPart::zero()) + ), + multiply( + Monomial::new(1 << i, PPart::zero()), + Monomial::new(1 << i, PPart::zero()) + ) ); } } @@ -1689,8 +1866,8 @@ mod tests { continue; } assert_eq!( - multiply_closed(a, b), - multiply(a, b), + multiply_closed(*a, *b), + multiply(*a, *b), "closed form ≠ oracle for {a:?} · {b:?}" ); } From 5322999028ea7acc0ff6cf43547138a3fdd78a46 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Thu, 27 Aug 2026 22:35:22 -0400 Subject: [PATCH 04/15] algebra: extract next_disjoint from the Milnor multiplier `PPartMultiplier::next_val` steps to the next entry value that keeps the Milnor coefficient non-zero, rather than testing candidates and discarding them. At p = 2 (and not mod 4) that is one branch-free expression, and it is reusable: the motivic closed-form product needs the same predicate on its anti-diagonals, where it is currently a filter over precomputed candidates. Testing it against brute force turned up a precondition the original code satisfied implicitly and never stated: `k` must already be disjoint from `sum`. The increment may carry through `k`'s own bits but not through `sum`'s, so an overlapping `k` can come back *smaller* -- `next_disjoint(2, 2)` is 1, not 4. Every caller walks a matrix whose anti-diagonal entries are pairwise disjoint, so it holds, but it is now documented and `debug_assert`ed. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01F25ZVbsP7ULg41iY3MP6FX --- .../algebra/src/algebra/milnor_algebra.rs | 46 ++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index d62ec0793b..38d54960df 100644 --- a/ext/crates/algebra/src/algebra/milnor_algebra.rs +++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs @@ -1431,6 +1431,27 @@ impl PPartAllocation { } } +/// The least `l > k` whose binary digits are disjoint from those of `sum`, **given that `k` is +/// itself disjoint from `sum`**. +/// +/// Equivalently, the least such `l` with $\binom{\mathrm{sum} + l}{l}$ odd: adding `l` to `sum` +/// carries exactly where they share a set bit, and the 2-adic valuation of that binomial is the +/// number of carries (Kummer). So this steps straight to the next value that keeps a Milnor +/// coefficient non-zero mod 2, instead of testing candidates and discarding them. +/// +/// # Panics +/// +/// In debug builds, if `k & sum != 0`. The identity genuinely needs it: the increment is allowed +/// to carry through `k`'s own bits but not through `sum`'s, so a `k` that overlaps `sum` can come +/// back *smaller* than `k` (`next_disjoint(2, 2) == 1`). Every caller is walking a matrix whose +/// entries are already pairwise disjoint along each anti-diagonal, so the precondition holds. +/// +/// The result can exceed any bound the caller has in mind; compare it against that separately. +pub const fn next_disjoint(sum: PPartEntry, k: PPartEntry) -> PPartEntry { + debug_assert!(k & sum == 0, "next_disjoint needs k disjoint from sum"); + ((k | sum) + 1) & !sum +} + #[allow(non_snake_case)] pub struct PPartMultiplier { p: ValidPrime, @@ -1531,7 +1552,7 @@ impl PPartMultiplier { }) .unwrap_or(max + 1) } else { - ((k | sum) + 1) & !sum + next_disjoint(sum, k) } } _ => (k + 1..max + 1) @@ -2206,6 +2227,29 @@ mod tests { } } + /// `next_disjoint` is the jump-to-valid form of the "is this binomial odd" test that the + /// Milnor coefficient needs; check it against the brute-force search it replaces, over every + /// input satisfying its precondition. + #[test] + fn next_disjoint_matches_brute_force() { + for sum in 0..64u32 { + for k in (0..64u32).filter(|k| k & sum == 0) { + let expected = (k + 1..) + .find(|l| l & sum == 0) + .expect("a disjoint value always exists"); + assert_eq!(next_disjoint(sum, k), expected, "sum = {sum}, k = {k}"); + // The characterisation it is actually used for. + assert_ne!(u32::binomial2(sum + expected, expected), 0); + } + } + } + + #[test] + #[should_panic(expected = "next_disjoint needs k disjoint from sum")] + fn next_disjoint_rejects_overlapping_k() { + next_disjoint(2, 2); + } + /// The packing is only sound because each field is wide enough for every entry that can occur /// at degree at most `MAX_DEGREE`. Check that against the $\xi$-degrees directly, so that /// changing `MAX_DEGREE` or `WIDTHS` without the other fails loudly. From 1ba525f6dbec3ddd6919edb58ac2d3602915c327 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Fri, 28 Aug 2026 23:50:42 -0400 Subject: [PATCH 05/15] motivic: address the review comments on the tau and coefficient code `Tau::pow` and `Tau::shift` had no callers outside their own test, so the question of whether `0^0` should be `1` goes away with them. The `get` closure in `c_coeff` captured nothing and becomes a free function, and its Euclidean division by a power of two becomes the arithmetic shift it compiles to. Co-Authored-By: Claude Opus 5 --- .../algebra/src/algebra/motivic/milnor.rs | 14 +++++++---- ext/crates/algebra/src/algebra/motivic/mod.rs | 3 +-- ext/crates/algebra/src/algebra/motivic/tau.rs | 24 ------------------- 3 files changed, 11 insertions(+), 30 deletions(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index ce271df9cf..d92ed91c37 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -147,6 +147,11 @@ fn sigma2(r: &[u32]) -> u64 { r.iter().enumerate().map(|(i, &v)| (v as u64) << i).sum() } +/// The `i`-th entry of an exponent sequence, widened, reading past the end as $0$. +fn entry(seq: &[u32], i: usize) -> i32 { + seq.get(i).copied().unwrap_or(0) as i32 +} + /// The coefficient $c(S, R)$ of Notation 3.3 / Theorem 3.4, reduced mod 2. /// /// $c(S, R) = \prod_{n \ge 1} \binom{\lfloor \sum_{i=0}^{n-1} 2^{i-n}(s_i - r_i)\rfloor}{r_n}$ @@ -158,16 +163,17 @@ fn c_coeff(s: &[u32], r: &[u32]) -> u32 { if r.first().copied().unwrap_or(0) > 0 { return 0; } - let get = |seq: &[u32], i: usize| seq.get(i).copied().unwrap_or(0) as i32; let len = s.len().max(r.len()); let mut prod = 1u32; for n in 1..len { let mut num: i32 = 0; for i in 0..n { - num += (get(s, i) - get(r, i)) << i; + num += (entry(s, i) - entry(r, i)) << i; } - let floor = num.div_euclid(1i32 << n); - prod *= i32::binomial2(floor, get(r, n)) as u32; + // Arithmetic shift right is exactly the Euclidean quotient by a power of two, and + // stays so for the negative `num` that a large $R$ produces. + let floor = num >> n; + prod *= i32::binomial2(floor, entry(r, n)) as u32; if prod == 0 { return 0; } diff --git a/ext/crates/algebra/src/algebra/motivic/mod.rs b/ext/crates/algebra/src/algebra/motivic/mod.rs index 94f1fa1bd3..941a228b0e 100644 --- a/ext/crates/algebra/src/algebra/motivic/mod.rs +++ b/ext/crates/algebra/src/algebra/motivic/mod.rs @@ -1,5 +1,4 @@ -//! The C-motivic (over $\mathbb{C}$, prime 2) Steenrod algebra and its mod-$\tau$ -//! reduction. +//! The C-motivic prime 2 Steenrod algebra and its mod-$\tau$ reduction. //! //! This layer implements the *deformation* view of the C-motivic Adams $E_2$: the //! C-motivic dual Steenrod algebra $A_C$ over $\mathbb{F}_2[\tau]$, its mod-$\tau$ diff --git a/ext/crates/algebra/src/algebra/motivic/tau.rs b/ext/crates/algebra/src/algebra/motivic/tau.rs index 1a6e53c526..54a5b12aa1 100644 --- a/ext/crates/algebra/src/algebra/motivic/tau.rs +++ b/ext/crates/algebra/src/algebra/motivic/tau.rs @@ -53,20 +53,6 @@ impl Tau { pub fn valuation(self) -> Option { self.0 } - - /// The $n$-th power $(\tau^k)^n = \tau^{kn}$, with the convention $x^0 = 1$. - pub fn pow(self, n: u32) -> Self { - if n == 0 { - Self::ONE - } else { - Self(self.0.map(|k| k * n)) - } - } - - /// Multiply by $\tau^k$, i.e. raise the valuation by `k` (zero stays zero). - pub fn shift(self, k: u32) -> Self { - Self(self.0.map(|v| v + k)) - } } impl std::ops::Mul for Tau { @@ -129,14 +115,4 @@ mod tests { assert_eq!(Tau::zero() + Tau::power(5), Tau::power(5)); assert_eq!(Tau::zero() + Tau::zero(), Tau::zero()); } - - #[test] - fn test_pow_and_shift() { - assert_eq!(Tau::power(3).pow(2), Tau::power(6)); - assert_eq!(Tau::power(3).pow(0), Tau::one()); - assert_eq!(Tau::zero().pow(0), Tau::one()); - assert_eq!(Tau::zero().pow(3), Tau::zero()); - assert_eq!(Tau::power(2).shift(3), Tau::power(5)); - assert_eq!(Tau::zero().shift(3), Tau::zero()); - } } From 6c11c7322372349c69bea96f68209942ea12b620 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sat, 29 Aug 2026 17:30:42 -0400 Subject: [PATCH 06/15] motivic: drop the tau coefficient, which the weight already determines MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Everything the engine computes is bidegree-homogeneous and tau has bidegree (0, -1), so the topological degree is additive and the whole tau power lands in the weight: the coefficient of a term z in a product a*b can only be tau^(w_z - w_a - w_b). Nothing was free to store. So a `SparseSum` becomes a mod-2 set of keys and `product_indexed` returns bare indices. A caller that wants a coefficient asks `Grading::tau_exponent` for it. `Grading` names which of the two dual weight conventions is in play — A_C weights a basis element by the negative of the monomial it pairs with — because the exponent formula is the same on both sides once each is asked for its own weight, and a bare sign is not. `rewrite_tau` still counts the exponent it always did; `mul_monomials` now debug-asserts that it agrees with what the weights dictate, which is the invariant that licenses not keeping it. This is a wash for speed: `motivic_basis`, which the change cannot touch, moved 4% on its own, and every other group moved less than that. The point is the smaller representation — a cached structure constant halves from 16 bytes to 8 — and one less thing to thread through the resolution. Co-Authored-By: Claude Opus 5 --- ext/crates/algebra/src/algebra/mod.rs | 2 +- .../algebra/src/algebra/motivic/milnor.rs | 411 +++++++++--------- ext/crates/algebra/src/algebra/motivic/mod.rs | 24 +- ext/crates/algebra/src/algebra/motivic/tau.rs | 118 ----- 4 files changed, 228 insertions(+), 327 deletions(-) delete mode 100644 ext/crates/algebra/src/algebra/motivic/tau.rs diff --git a/ext/crates/algebra/src/algebra/mod.rs b/ext/crates/algebra/src/algebra/mod.rs index 44c0f4454a..0bb7e9fa1f 100644 --- a/ext/crates/algebra/src/algebra/mod.rs +++ b/ext/crates/algebra/src/algebra/mod.rs @@ -19,7 +19,7 @@ pub mod milnor_algebra; pub use milnor_algebra::MilnorAlgebra; pub mod motivic; -pub use motivic::{MotivicMilnorAlgebra, Tau}; +pub use motivic::MotivicMilnorAlgebra; mod steenrod_algebra; pub use steenrod_algebra::{AlgebraType, SteenrodAlgebra}; diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index d92ed91c37..8a66a3efa5 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -29,15 +29,16 @@ //! implementation; the closed-form product (Kong–Lin Theorem 5.1) can replace it //! later for speed, validated against this one. //! -//! Coefficients are homogeneous, so each one is a single power of $\tau$: a -//! [`Tau`] scalar (see [`super::tau`]). This is the $A_C$ product **engine**, -//! over $\mathbb{F}_2[\tau]$; the mod-$\tau$ reduction $A_C/\tau$ is presented to -//! the resolution engine as an ordinary $\mathbb{F}_2$-algebra in a follow-up -//! layer built on top of this engine. +//! Coefficients are never stored. Everything here is homogeneous, so each one is a single power +//! of $\tau$, and $\tau$ has weight $-1$ — which pins the exponent to the weights of the terms +//! it sits between ([`Grading::tau_exponent`]). This is the $A_C$ product **engine**, over +//! $\mathbb{F}_2[\tau]$; the mod-$\tau$ reduction $A_C/\tau$ is presented to the resolution +//! engine as an ordinary $\mathbb{F}_2$-algebra in a follow-up layer built on top of this engine. //! //! Weight convention: the motivic weight of an algebra basis element is the //! *negative* of the weight of the dual monomial it pairs with, so that products -//! are weight-homogeneous with $\tau$ (weight $-1$) absorbing the difference. +//! are weight-homogeneous with $\tau$ (weight $-1$) absorbing the difference. [`Grading`] names +//! the two sides. use std::{ cell::RefCell, @@ -51,7 +52,6 @@ use maybe_rayon::prelude::*; use once::OnceVec; use rustc_hash::FxHashMap; -use super::Tau; use crate::algebra::milnor_algebra::PPart; /// Drop the trailing zeros from an exponent vector. @@ -284,18 +284,17 @@ fn rewrite_tau_dfs( /// A sparse $\mathbb{F}_2[\tau]$-linear combination of `K`s, as a `Vec` sorted by key. /// -/// Everything we compute is homogeneous — $\psi$ and the products are graded — so every -/// coefficient is a single power of $\tau$: a [`Tau`] scalar. Zero coefficients are never stored -/// and the terms are kept sorted, so equality is canonical, and coefficient arithmetic is exactly -/// [`Tau`]'s arithmetic — `mul` adds valuations, `add` cancels equal powers mod 2 (unequal powers -/// would be inhomogeneous and cannot arise). +/// Coefficients are not stored, because the bigrading already determines them: everything we +/// compute is homogeneous, $\tau$ has weight $-1$, so the coefficient of a term is forced to be +/// $\tau^{\text{(input weight)} - \text{(term weight)}}$ — see [`Grading::tau_exponent`]. What is left is +/// a set of keys taken mod 2, and adding a term toggles it. /// /// A sorted `Vec` rather than a `BTreeMap`: these combinations are small — a product of two basis -/// elements has a handful of terms — and the keys are now `Copy` words, so one contiguous -/// allocation searched by [`slice::binary_search_by`] beats a tree that allocates per node. +/// elements has a handful of terms — and the keys are `Copy` words, so one contiguous allocation +/// searched by [`slice::binary_search`] beats a tree that allocates per node. #[derive(Debug, Clone, Default, PartialEq, Eq)] pub struct SparseSum { - terms: Vec<(K, Tau)>, + terms: Vec, } impl SparseSum { @@ -314,59 +313,46 @@ impl SparseSum { self.terms.is_empty() } - /// The coefficient of `key`, if non-zero. - pub fn get(&self, key: &K) -> Option<&Tau> { - self.position(key).ok().map(|i| &self.terms[i].1) + /// Whether `key` occurs. + pub fn contains(&self, key: &K) -> bool { + self.terms.binary_search(key).is_ok() } - /// Add `coeff * key`, dropping the term if it cancels to zero. - pub fn add_term(&mut self, key: K, coeff: Tau) { - if coeff.is_zero() { - return; - } - match self.position(&key) { + /// Add `key`, cancelling mod 2 against an occurrence already present. + pub fn add_term(&mut self, key: K) { + match self.terms.binary_search(&key) { Ok(i) => { - let sum = self.terms[i].1 + coeff; - if sum.is_zero() { - self.terms.remove(i); - } else { - self.terms[i].1 = sum; - } + self.terms.remove(i); } - Err(i) => self.terms.insert(i, (key, coeff)), + Err(i) => self.terms.insert(i, key), } } - /// The index of `key`, or where it would be inserted. - fn position(&self, key: &K) -> Result { - self.terms.binary_search_by(|(k, _)| k.cmp(key)) - } - /// The terms, in key order. - pub fn iter(&self) -> std::slice::Iter<'_, (K, Tau)> { + pub fn iter(&self) -> std::slice::Iter<'_, K> { self.terms.iter() } } -impl FromIterator<(K, Tau)> for SparseSum { - fn from_iter>(iter: I) -> Self { +impl FromIterator for SparseSum { + fn from_iter>(iter: I) -> Self { let mut out = Self::new(); - for (key, coeff) in iter { - out.add_term(key, coeff); + for key in iter { + out.add_term(key); } out } } -impl From<[(K, Tau); N]> for SparseSum { - fn from(terms: [(K, Tau); N]) -> Self { +impl From<[K; N]> for SparseSum { + fn from(terms: [K; N]) -> Self { terms.into_iter().collect() } } impl IntoIterator for SparseSum { - type IntoIter = std::vec::IntoIter<(K, Tau)>; - type Item = (K, Tau); + type IntoIter = std::vec::IntoIter; + type Item = K; fn into_iter(self) -> Self::IntoIter { self.terms.into_iter() @@ -374,16 +360,60 @@ impl IntoIterator for SparseSum { } impl<'a, K> IntoIterator for &'a SparseSum { - type IntoIter = std::slice::Iter<'a, (K, Tau)>; - type Item = &'a (K, Tau); + type IntoIter = std::slice::Iter<'a, K>; + type Item = &'a K; fn into_iter(self) -> Self::IntoIter { self.terms.iter() } } -/// An element of the dual algebra $A_{**}$: [`Monomial`]s with $\mathbb{F}_2[\tau]$ -/// coefficients. +/// Which of the two dual weight conventions a computation is expressed in. +/// +/// $A_C$ and its dual $A_{**}$ weight the same [`Monomial`] by negatives of each other: the +/// monomial's own weight grades $A_{**}$, and the basis element of $A_C$ it pairs with takes the +/// negative (see the module docs). Naming the convention keeps that sign in one place — the +/// formula for a $\tau$ exponent is the same on both sides once each is asked for its own weight. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Grading { + /// $A_{**}$, weighted by the monomial itself — the convention of [`dual_mul`]. + Dual, + /// $A_C$, weighted by the basis element the monomial is dual to — the convention of + /// [`multiply`] and [`product_indexed`](MotivicMilnorAlgebra::product_indexed). + Steenrod, +} + +impl Grading { + /// The motivic weight of `m` in this convention. + pub fn weight(self, m: Monomial) -> i32 { + let w = m.bidegree().1; + match self { + Self::Dual => w, + Self::Steenrod => -w, + } + } + + /// The exponent of $\tau$ on the term `z` of a product `a * b`. + /// + /// The bigrading forces it, which is why products store no coefficients. $\tau$ has weight + /// $-1$ and topological degree $0$, so a weight-homogeneous product can only carry + /// $\tau^{w_z - w_a - w_b}$ on `z`. + /// + /// # Panics + /// + /// In debug builds, if the weights imply a negative power — i.e. `z` cannot occur in `a * b`. + pub fn tau_exponent(self, a: Monomial, b: Monomial, z: Monomial) -> u32 { + let k = self.weight(z) - self.weight(a) - self.weight(b); + debug_assert!( + k >= 0, + "negative tau exponent: {a:?} * {b:?} → {z:?} ({self:?})" + ); + k as u32 + } +} + +/// An element of the dual algebra $A_{**}$: a mod-2 set of [`Monomial`]s, whose +/// $\mathbb{F}_2[\tau]$ coefficients are implicit in the weight ([`Grading::tau_exponent`]). pub type DualElement = SparseSum; /// Elementwise sum of two exponent sequences, i.e. the product of the two $\xi$ monomials. @@ -401,7 +431,7 @@ fn ppart_add(a: PPart, b: PPart) -> PPart { /// exponents: the exterior parts form $S = E_1 + E_2$ (entries in $\{0,1,2\}$), which is /// rewritten into the square-free basis via [`rewrite_tau`] ($\tau_i^2 = \tau\xi_{i+1}$), and /// the resulting $\xi$ exponents are added to $R_1 + R_2$. -fn mul_monomials(m1: Monomial, m2: Monomial, coeff: Tau, acc: &mut DualElement) { +fn mul_monomials(m1: Monomial, m2: Monomial, acc: &mut DualElement) { let (e1, e2) = (m1.q_part, m2.q_part); let bits = u32::BITS - (e1 | e2).leading_zeros(); let s: Vec = (0..bits) @@ -414,19 +444,20 @@ fn mul_monomials(m1: Monomial, m2: Monomial, coeff: Tau, acc: &mut DualElement) let term_r = PPart::try_from_slice(term.r.get(1..).unwrap_or_default()) .expect("rewrite_tau exponent out of packing range"); let r = ppart_add(term_r, r12); - acc.add_term( - Monomial::new(term.e_mask, r), - coeff * Tau::power(term.tau_pow), - ); + let z = Monomial::new(term.e_mask, r); + // The τ power `rewrite_tau` computed and the one the weight dictates are the same + // number; this is the invariant that lets the coefficient go unstored. + debug_assert_eq!(term.tau_pow, Grading::Dual.tau_exponent(m1, m2, z)); + acc.add_term(z); } } /// The product of two elements of $A_{**}$. pub fn dual_mul(a: &DualElement, b: &DualElement) -> DualElement { let mut out = DualElement::new(); - for &(m1, c1) in a { - for &(m2, c2) in b { - mul_monomials(m1, m2, c1 * c2, &mut out); + for &m1 in a { + for &m2 in b { + mul_monomials(m1, m2, &mut out); } } out @@ -434,17 +465,17 @@ pub fn dual_mul(a: &DualElement, b: &DualElement) -> DualElement { /// The generator $\tau_i \in A_{**}$. pub fn tau_gen(i: usize) -> DualElement { - DualElement::from([(Monomial::new(1 << i, PPart::zero()), Tau::one())]) + DualElement::from([Monomial::new(1 << i, PPart::zero())]) } /// The generator $\xi_i \in A_{**}$ (for $i \ge 1$). pub fn xi_gen(i: usize) -> DualElement { - DualElement::from([(xi_pow_mon(i, 1), Tau::one())]) + DualElement::from([xi_pow_mon(i, 1)]) } /// The unit $1 \in A_{**}$. pub fn dual_one() -> DualElement { - DualElement::from([(Monomial::one(), Tau::one())]) + DualElement::from([Monomial::one()]) } // --------------------------------------------------------------------------- @@ -457,22 +488,21 @@ pub type TensorElement = SparseSum<(Monomial, Monomial)>; /// The unit $1 \otimes 1$. fn tensor_one() -> TensorElement { - TensorElement::from([((Monomial::one(), Monomial::one()), Tau::one())]) + TensorElement::from([(Monomial::one(), Monomial::one())]) } /// Multiply in $A_{**} \otimes A_{**}$: $(a_L \otimes a_R)(b_L \otimes b_R) = (a_L b_L) \otimes (a_R b_R)$. fn tensor_mul(t1: &TensorElement, t2: &TensorElement) -> TensorElement { let mut out = TensorElement::new(); - for &((al, ar), c1) in t1 { - for &((bl, br), c2) in t2 { + for &(al, ar) in t1 { + for &(bl, br) in t2 { let mut left = DualElement::new(); - mul_monomials(al, bl, Tau::one(), &mut left); + mul_monomials(al, bl, &mut left); let mut right = DualElement::new(); - mul_monomials(ar, br, Tau::one(), &mut right); - let c = c1 * c2; - for &(ml, cl) in &left { - for &(mr, cr) in &right { - out.add_term((ml, mr), c * (cl * cr)); + mul_monomials(ar, br, &mut right); + for &ml in &left { + for &mr in &right { + out.add_term((ml, mr)); } } } @@ -492,18 +522,12 @@ fn xi_pow_mon(j: usize, p: u32) -> Monomial { /// $\psi(\tau_k) = 1 \otimes \tau_k + \sum_{i=0}^{k} \tau_i \otimes \xi_{k-i}^{2^i}$. fn coprod_tau(k: usize) -> TensorElement { let mut out = TensorElement::new(); - out.add_term( - (Monomial::one(), Monomial::new(1 << k, PPart::zero())), - Tau::one(), - ); + out.add_term((Monomial::one(), Monomial::new(1 << k, PPart::zero()))); for i in 0..=k { - out.add_term( - ( - Monomial::new(1 << i, PPart::zero()), - xi_pow_mon(k - i, 1 << i), - ), - Tau::one(), - ); + out.add_term(( + Monomial::new(1 << i, PPart::zero()), + xi_pow_mon(k - i, 1 << i), + )); } out } @@ -512,7 +536,7 @@ fn coprod_tau(k: usize) -> TensorElement { fn coprod_xi(k: usize) -> TensorElement { let mut out = TensorElement::new(); for i in 0..=k { - out.add_term((xi_pow_mon(i, 1), xi_pow_mon(k - i, 1 << i)), Tau::one()); + out.add_term((xi_pow_mon(i, 1), xi_pow_mon(k - i, 1 << i))); } out } @@ -535,9 +559,9 @@ fn coproduct_monomial(m: &Monomial) -> TensorElement { /// The coproduct of an arbitrary element of $A_{**}$ (extended $\mathbb{F}_2[\tau]$-linearly). pub fn coproduct(elt: &DualElement) -> TensorElement { let mut out = TensorElement::new(); - for &(m, c) in elt { - for (key, cc) in coproduct_monomial(&m) { - out.add_term(key, c * cc); + for &m in elt { + for key in coproduct_monomial(&m) { + out.add_term(key); } } out @@ -549,7 +573,7 @@ pub fn coproduct(elt: &DualElement) -> TensorElement { /// $\xi_j^p \in A_{**}$ (with $\xi_0 = 1$). fn xi_pow_elt(j: usize, p: u32) -> DualElement { - DualElement::from([(xi_pow_mon(j, p), Tau::one())]) + DualElement::from([xi_pow_mon(j, p)]) } /// $\chi(\xi_k)$, from the antipode axiom with this module's coproduct: @@ -560,9 +584,8 @@ fn chi_xi(k: usize) -> DualElement { } let mut acc = DualElement::new(); for i in 0..k { - let term = dual_mul(&chi_xi(i), &xi_pow_elt(k - i, 1 << i)); - for (mon, c) in term { - acc.add_term(mon, c); + for mon in dual_mul(&chi_xi(i), &xi_pow_elt(k - i, 1 << i)) { + acc.add_term(mon); } } acc @@ -573,9 +596,8 @@ fn chi_xi(k: usize) -> DualElement { fn chi_tau(k: usize) -> DualElement { let mut acc = tau_gen(k); for i in 0..k { - let term = dual_mul(&chi_tau(i), &xi_pow_elt(k - i, 1 << i)); - for (mon, c) in term { - acc.add_term(mon, c); + for mon in dual_mul(&chi_tau(i), &xi_pow_elt(k - i, 1 << i)) { + acc.add_term(mon); } } acc @@ -587,7 +609,7 @@ fn chi_tau(k: usize) -> DualElement { /// involution. pub fn antipode(elt: &DualElement) -> DualElement { let mut out = DualElement::new(); - for &(mon, c) in elt { + for &mon in elt { let mut acc = dual_one(); for i in 0..u32::BITS { if (mon.q_part >> i) & 1 != 0 { @@ -599,8 +621,8 @@ pub fn antipode(elt: &DualElement) -> DualElement { acc = dual_mul(&acc, &chi_xi(i + 1)); } } - for (out_mon, cc) in acc { - out.add_term(out_mon, c * cc); + for out_mon in acc { + out.add_term(out_mon); } } out @@ -675,8 +697,8 @@ fn enum_basis_dfs( } /// The product `a · b` in the C-motivic Steenrod algebra `A_C`, where `a`, `b` are Milnor basis -/// elements `Q(E)P(R)` given as the monomials `(E, R)` they are dual to. The result is a map -/// from basis monomials to their `𝔽₂[τ]` coefficients. +/// elements `Q(E)P(R)` given as the monomials `(E, R)` they are dual to. The result is the set of +/// basis monomials occurring, their `𝔽₂[τ]` coefficients being implicit ([`Grading::tau_exponent`]). /// /// Computed by duality: the coefficient of `z` in `a · b` is the coefficient of /// `mon(a) ⊗ mon(b)` in `ψ(mon(z))`, summed over the basis `z` of the appropriate @@ -686,10 +708,8 @@ pub fn multiply(a: Monomial, b: Monomial) -> DualElement { let t = a.bidegree().0 + b.bidegree().0; let mut out = DualElement::new(); for z in enum_basis(t) { - if let Some(&c) = coproduct_monomial(&z).get(&key) - && !c.is_zero() - { - out.add_term(z, c); + if coproduct_monomial(&z).contains(&key) { + out.add_term(z); } } out @@ -1021,10 +1041,7 @@ impl ClosedY<'_> { out_e_mask |= 1 << i; } } - out.add_term( - Monomial::new(out_e_mask, self.out_r), - Tau::power(self.sigma_sprime), - ); + out.add_term(Monomial::new(out_e_mask, self.out_r)); } } @@ -1075,8 +1092,9 @@ pub fn multiply_closed(a: Monomial, b: Monomial) -> DualElement { /// per-(topological-)degree basis indexing. /// /// This is the $A_C$ product **engine** over $\mathbb{F}_2[\tau]$: the dual-based product is -/// exposed through [`product_indexed`](MotivicMilnorAlgebra::product_indexed), whose coefficients -/// are [`Tau`] scalars. It is deliberately not an [`Algebra`](crate::algebra::Algebra) +/// exposed through [`product_indexed`](MotivicMilnorAlgebra::product_indexed), whose +/// $\tau$ coefficients are implicit in the weights ([`Grading::tau_exponent`]). It is +/// deliberately not an [`Algebra`](crate::algebra::Algebra) /// implementation, because that trait is over $\mathbb{F}_p$; the mod-$\tau$ reduction, which *is* /// such an algebra, is a follow-up layer on top of this engine, and the honest /// $\mathbb{F}_2[\tau]$ resolution is built by lifting against this engine (Phase 2). @@ -1117,7 +1135,7 @@ pub struct MotivicMilnorAlgebra { pub struct ProductBlock { /// The number of basis elements in degree `t2` (the block's row stride). dim2: usize, - entries: Vec>>, + entries: Vec>>, } impl MotivicMilnorAlgebra { @@ -1158,8 +1176,19 @@ impl MotivicMilnorAlgebra { /// The `(topological degree, motivic weight)` of a basis element. pub fn bidegree(&self, degree: i32, idx: usize) -> (i32, i32) { - let (t, w) = self.basis[degree as usize][idx].bidegree(); - (t, -w) + let m = self.basis[degree as usize][idx]; + (m.bidegree().0, Grading::Steenrod.weight(m)) + } + + /// The exponent of $\tau$ on basis element `idx` of degree `t1 + t2` in the product of the + /// elements at `(t1, idx1)` and `(t2, idx2)` — the coefficient [`Self::product_indexed`] + /// leaves implicit. See [`Grading::tau_exponent`]. + pub fn tau_exponent(&self, t1: i32, idx1: usize, t2: i32, idx2: usize, idx: usize) -> u32 { + Grading::Steenrod.tau_exponent( + self.basis[t1 as usize][idx1], + self.basis[t2 as usize][idx2], + self.basis[(t1 + t2) as usize][idx], + ) } /// Parse a basis element written as `basis_element_to_string` prints it — @@ -1246,27 +1275,25 @@ impl MotivicMilnorAlgebra { /// `(t1, idx1)` and `(t2, idx2)`, mapped into the degree-`(t1+t2)` basis. The /// bases must already be computed (the caller through [`Self::block`] ensures /// this). - fn compute_product(&self, t1: i32, idx1: usize, t2: i32, idx2: usize) -> Vec<(Tau, usize)> { + fn compute_product(&self, t1: i32, idx1: usize, t2: i32, idx2: usize) -> Vec { let a = &self.basis[t1 as usize][idx1]; let b = &self.basis[t2 as usize][idx2]; let t = t1 + t2; multiply_closed(*a, *b) .into_iter() - .map(|(z, c)| { - ( - c, - self.index_of(t, &z) - .expect("product landed outside the basis"), - ) + .map(|z| { + self.index_of(t, &z) + .expect("product landed outside the basis") }) .collect() } - /// The product of two basis elements, as an $\mathbb{F}_2[\tau]$-linear combination of basis - /// elements in degree `t1 + t2`: a list of `(coefficient, index)` pairs. The product of two - /// homogeneous basis elements is weight-homogeneous, so each coefficient is a single power of - /// $\tau$ — a [`Tau`] scalar. - pub fn product_indexed(&self, t1: i32, idx1: usize, t2: i32, idx2: usize) -> Vec<(Tau, usize)> { + /// The product of two basis elements as a list of basis indices in degree `t1 + t2`. + /// + /// The $\mathbb{F}_2[\tau]$ coefficients are not returned because they are not free: the + /// product of two homogeneous basis elements is weight-homogeneous, so index `j` carries + /// exactly $\tau^{k}$ for the `k` that [`Self::tau_exponent`] computes from the weights. + pub fn product_indexed(&self, t1: i32, idx1: usize, t2: i32, idx2: usize) -> Vec { self.cached_product(&self.block(t1, t2), t1, idx1, t2, idx2) .to_vec() } @@ -1281,7 +1308,7 @@ impl MotivicMilnorAlgebra { idx1: usize, t2: i32, idx2: usize, - ) -> &'a [(Tau, usize)] { + ) -> &'a [usize] { block.entries[idx1 * block.dim2 + idx2] .get_or_init(|| self.compute_product(t1, idx1, t2, idx2)) } @@ -1402,14 +1429,11 @@ mod tests { let terms: DualElement = alg .product_indexed(1, 0, 2, 0) .into_iter() - .map(|(c, idx)| (*alg.basis_element(3, idx), c)) + .map(|idx| *alg.basis_element(3, idx)) .collect(); assert_eq!( terms, - DualElement::from([ - (mon(0b10, &[]), Tau::one()), - (mon(0b1, &[0, 1]), Tau::one()) - ]) + DualElement::from([mon(0b10, &[]), mon(0b1, &[0, 1])]) ); } @@ -1427,15 +1451,19 @@ mod tests { let indexed: DualElement = alg .product_indexed(t1, idx1, t2, idx2) .into_iter() - .map(|(c, idx)| (*alg.basis_element(t1 + t2, idx), c)) + .map(|idx| *alg.basis_element(t1 + t2, idx)) .collect(); assert_eq!(indexed, multiply(a, b)); - // Weight-homogeneous: w_out - (τ-power) = w_a + w_b. - let (wa, wb) = (alg.bidegree(t1, idx1).1, alg.bidegree(t2, idx2).1); - for (c, idx) in alg.product_indexed(t1, idx1, t2, idx2) { - let w = alg.bidegree(t1 + t2, idx).1; - assert_eq!(w - c.valuation().unwrap() as i32, wa + wb); + // The indexed τ exponent resolves its indices to the same monomials the + // free function is given. (That the exponent is the *right* one is + // checked inside `mul_monomials`, against `rewrite_tau`'s own count.) + for idx in alg.product_indexed(t1, idx1, t2, idx2) { + let z = *alg.basis_element(t1 + t2, idx); + assert_eq!( + alg.tau_exponent(t1, idx1, t2, idx2, idx), + Grading::Steenrod.tau_exponent(a, b, z) + ); } } } @@ -1445,26 +1473,32 @@ mod tests { #[test] fn test_dual_mul_relations() { - // τ_0^2 = τ ξ_1 (the defining relation at ρ = 0). + // τ_0^2 = τ ξ_1 (the defining relation at ρ = 0). The τ^1 is read back off the weights. + let (t0, t1) = (mon(0b1, &[]), mon(0b10, &[])); assert_eq!( dual_mul(&tau_gen(0), &tau_gen(0)), - DualElement::from([(mon(0, &[0, 1]), Tau::power(1))]) // τ^1 · ξ_1 + DualElement::from([mon(0, &[0, 1])]) ); - // ξ_1^2 is just the monomial ξ_1^2. + assert_eq!(Grading::Dual.tau_exponent(t0, t0, mon(0, &[0, 1])), 1); + // ξ_1^2 is just the monomial ξ_1^2, with no τ. assert_eq!( dual_mul(&xi_gen(1), &xi_gen(1)), - DualElement::from([(mon(0, &[0, 2]), Tau::one())]) + DualElement::from([mon(0, &[0, 2])]) ); - // Distinct τ's commute and stay square-free: τ_0 τ_1. + let xi1 = mon(0, &[0, 1]); + assert_eq!(Grading::Dual.tau_exponent(xi1, xi1, mon(0, &[0, 2])), 0); + // Distinct τ's commute and stay square-free: τ_0 τ_1, again τ-free. assert_eq!( dual_mul(&tau_gen(0), &tau_gen(1)), - DualElement::from([(mon(0b11, &[]), Tau::one())]) + DualElement::from([mon(0b11, &[])]) ); + assert_eq!(Grading::Dual.tau_exponent(t0, t1, mon(0b11, &[])), 0); // τ_1^2 = τ ξ_2. assert_eq!( dual_mul(&tau_gen(1), &tau_gen(1)), - DualElement::from([(mon(0, &[0, 0, 1]), Tau::power(1))]) + DualElement::from([mon(0, &[0, 0, 1])]) ); + assert_eq!(Grading::Dual.tau_exponent(t1, t1, mon(0, &[0, 0, 1])), 1); // Multiplication by the unit is the identity. assert_eq!(dual_mul(&dual_one(), &tau_gen(2)), tau_gen(2)); } @@ -1484,9 +1518,13 @@ mod tests { let motivic = |p1: &[u32], p2: &[u32]| -> Vec> { let mut out: Vec> = multiply(paper(p1), paper(p2)) .into_iter() - .filter(|(m, _)| m.q_part == 0) // keep the pure-ξ outputs - .map(|(m, c)| { - assert_eq!(c, Tau::one(), "pure-ξ output carried a τ power"); + .filter(|m| m.q_part == 0) // keep the pure-ξ outputs + .map(|m| { + assert_eq!( + Grading::Steenrod.tau_exponent(paper(p1), paper(p2), m), + 0, + "pure-ξ output carried a τ power" + ); m.p_part.iter().collect() }) .collect(); @@ -1536,11 +1574,8 @@ mod tests { .collect() }; let a_coeff = |mu_pp: &[u32], w_pp: &[u32]| -> u32 { - let ap = antipode(&DualElement::from([(paper(w_pp), Tau::one())])); - // The antipode of a pure-ξ element is pure-ξ and τ-free, so a present monomial has - // coefficient τ^0. - ap.get(&paper(mu_pp)) - .map_or(0, |&c| u32::from(c == Tau::one())) + let ap = antipode(&DualElement::from([paper(w_pp)])); + u32::from(ap.contains(&paper(mu_pp))) }; let cb_mul = |a_pp: &[u32], b_pp: &[u32], acc: &mut BTreeMap, u32>| { for pp in classical_mul(a_pp, b_pp) { @@ -1561,8 +1596,8 @@ mod tests { // LHS: motivic conjugate product (E=0), converted to std via χ. let conj: Vec> = multiply(paper(&r1), paper(&r2)) .into_iter() - .filter(|(m, _)| m.q_part == 0) - .map(|(m, _)| m.p_part.iter().collect()) + .filter(|m| m.q_part == 0) + .map(|m| m.p_part.iter().collect()) .collect(); let mut lhs: BTreeMap, u32> = BTreeMap::new(); for w in xi_pps(mdeg(&r1) + mdeg(&r2)) { @@ -1596,13 +1631,13 @@ mod tests { #[test] fn test_product_associative() { // Associativity (a·b)·c = a·(b·c) is convention-independent and a strong global check. - // Extend the basis×basis product to element×basis, F_2[τ]-linearly (add_term handles the - // coefficient sum + zero cancellation via Tau arithmetic). + // Extend the basis×basis product to element×basis, F_2[τ]-linearly (`add_term` cancels + // mod 2, and the τ coefficients agree termwise because both sides are homogeneous). fn mul_elt_basis(x: &DualElement, c: &Monomial) -> DualElement { let mut out = DualElement::new(); - for &(z, cz) in x { - for (w, cw) in multiply(z, *c) { - out.add_term(w, cz * cw); + for &z in x { + for w in multiply(z, *c) { + out.add_term(w); } } out @@ -1617,9 +1652,9 @@ mod tests { let bc = multiply(*b, *c); // a · (b·c): multiply basis `a` on the left of each term of bc. let mut rhs = DualElement::new(); - for &(z, cz) in &bc { - for (w, cw) in multiply(*a, z) { - rhs.add_term(w, cz * cw); + for &z in &bc { + for w in multiply(*a, z) { + rhs.add_term(w); } } assert_eq!(lhs, rhs, "associativity failed at {a:?},{b:?},{c:?}"); @@ -1637,8 +1672,8 @@ mod tests { assert_eq!( multiply(q0, p_xi1), DualElement::from([ - (mon(0b10, &[]), Tau::one()), // Q_1 - (mon(0b1, &[0, 1]), Tau::one()), // Q_0 P(ξ_1) + mon(0b10, &[]), // Q_1 + mon(0b1, &[0, 1]), // Q_0 P(ξ_1) ]) ); } @@ -1647,10 +1682,7 @@ mod tests { fn test_product_unit_and_squares() { // 1 · x = x. let x = mon(0b101, &[0, 2]); - assert_eq!( - multiply(mon(0, &[]), x), - DualElement::from([(x, Tau::one())]) - ); + assert_eq!(multiply(mon(0, &[]), x), DualElement::from([x])); // Q_i^2 = 0. (Q_0 = Sq^1 is the motivic Bockstein; P(ξ_1) = Sq^2 does NOT square to // zero — its square is τ Q_0 Q_1 + …, a genuine motivic feature.) for i in 0..3 { @@ -1667,23 +1699,20 @@ mod tests { #[test] fn test_product_weight_homogeneous_and_tau_appears() { - // Every product is weight-homogeneous: term z occurs with a single power τ^k, and by - // weight-preservation of ψ, k = (w_a + w_b) - w_z ≥ 0 in *dual*-monomial weights (the - // algebra weight of a basis element is the negative of its dual monomial's weight, so - // τ, of weight -1, absorbs the difference). Also confirm τ genuinely enters at least - // one product (i.e. the τ_i^2 = τξ_{i+1} relation fires). + // The τ exponent implied by the weights must agree with the one `rewrite_tau` counts — + // the invariant that licenses not storing coefficients at all. That comparison is the + // `debug_assert` in `mul_monomials`, which this sweep exercises; here we check the + // consequences it must have: the exponent is a genuine (non-negative) power, the + // topological degree is additive, and τ really does enter — i.e. the τ_i^2 = τξ_{i+1} + // relation fires rather than the whole thing being vacuously τ-free. let basis: Vec = (0..=6).flat_map(enum_basis).collect(); let mut saw_tau = false; for a in &basis { for b in &basis { - let (wa, wb) = (a.bidegree().1, b.bidegree().1); - for (z, c) in multiply(*a, *b) { - // Each coefficient is a single τ power by construction (a homogeneous Tau). - let k = c.valuation().expect("stored coefficients are never zero") as i32; - let wz = z.bidegree().1; - assert_eq!(wa + wb - wz, k, "weight mismatch: {a:?}*{b:?} → z={z:?}"); - assert!(k >= 0); - saw_tau |= k > 0; + for z in multiply(*a, *b) { + assert_eq!(z.bidegree().0, a.bidegree().0 + b.bidegree().0); + // Panics in debug if the weights imply a negative power. + saw_tau |= Grading::Steenrod.tau_exponent(*a, *b, z) > 0; } } } @@ -1699,10 +1728,7 @@ mod tests { assert_eq!(antipode(&xi_gen(1)), xi_gen(1)); assert_eq!( antipode(&xi_gen(2)), - DualElement::from([ - (mon(0, &[0, 0, 1]), Tau::one()), - (mon(0, &[0, 3]), Tau::one()) - ]) + DualElement::from([mon(0, &[0, 0, 1]), mon(0, &[0, 3])]) ); assert_eq!(antipode(&tau_gen(0)), tau_gen(0)); @@ -1722,13 +1748,10 @@ mod tests { // Antipode axiom: m(χ ⊗ id)ψ(x) = η ε(x); for positive-degree x this is 0. let counit_via_antipode = |x: &DualElement| -> DualElement { let mut out = DualElement::new(); - for ((l, r), c) in coproduct(x) { - let prod = dual_mul( - &antipode(&DualElement::from([(l, Tau::one())])), - &DualElement::from([(r, Tau::one())]), - ); - for (mon, cc) in prod { - out.add_term(mon, c * cc); + for (l, r) in coproduct(x) { + let prod = dual_mul(&antipode(&DualElement::from([l])), &DualElement::from([r])); + for mon in prod { + out.add_term(mon); } } out @@ -1753,18 +1776,18 @@ mod tests { assert_eq!( coproduct(&tau_gen(1)), TensorElement::from([ - ((mon(0, &[]), mon(0b10, &[])), Tau::one()), // 1 ⊗ τ_1 - ((mon(0b1, &[]), mon(0, &[0, 1])), Tau::one()), // τ_0 ⊗ ξ_1 - ((mon(0b10, &[]), mon(0, &[])), Tau::one()), // τ_1 ⊗ 1 + (mon(0, &[]), mon(0b10, &[])), // 1 ⊗ τ_1 + (mon(0b1, &[]), mon(0, &[0, 1])), // τ_0 ⊗ ξ_1 + (mon(0b10, &[]), mon(0, &[])), // τ_1 ⊗ 1 ]) ); // ψ(ξ_2) = 1⊗ξ_2 + ξ_1⊗ξ_1^2 + ξ_2⊗1 (matches Milnor). assert_eq!( coproduct(&xi_gen(2)), TensorElement::from([ - ((mon(0, &[]), mon(0, &[0, 0, 1])), Tau::one()), // 1 ⊗ ξ_2 - ((mon(0, &[0, 1]), mon(0, &[0, 2])), Tau::one()), // ξ_1 ⊗ ξ_1^2 - ((mon(0, &[0, 0, 1]), mon(0, &[])), Tau::one()), // ξ_2 ⊗ 1 + (mon(0, &[]), mon(0, &[0, 0, 1])), // 1 ⊗ ξ_2 + (mon(0, &[0, 1]), mon(0, &[0, 2])), // ξ_1 ⊗ ξ_1^2 + (mon(0, &[0, 0, 1]), mon(0, &[])), // ξ_2 ⊗ 1 ]) ); } diff --git a/ext/crates/algebra/src/algebra/motivic/mod.rs b/ext/crates/algebra/src/algebra/motivic/mod.rs index 941a228b0e..17579f36c1 100644 --- a/ext/crates/algebra/src/algebra/motivic/mod.rs +++ b/ext/crates/algebra/src/algebra/motivic/mod.rs @@ -1,27 +1,23 @@ //! The C-motivic prime 2 Steenrod algebra and its mod-$\tau$ reduction. //! //! This layer implements the *deformation* view of the C-motivic Adams $E_2$: the -//! C-motivic dual Steenrod algebra $A_C$ over $\mathbb{F}_2[\tau]$, its mod-$\tau$ -//! reduction $A_C/\tau$ (a connected finite-type $\mathbb{F}_2$-algebra), and the -//! coefficient ring $\mathbb{F}_2[\tau]$ itself. +//! C-motivic dual Steenrod algebra $A_C$ over $\mathbb{F}_2[\tau]$, and its mod-$\tau$ +//! reduction $A_C/\tau$ (a connected finite-type $\mathbb{F}_2$-algebra). //! -//! This module provides the foundation layer for that computation: +//! The foundation layer is [`MotivicMilnorAlgebra`] — $A_C$, a free +//! $\mathbb{F}_2[\tau]$-module on the Milnor basis, with the Kong–Lin product +//! ([`milnor`]). This is the product **engine**; it is deliberately not an +//! [`Algebra`](crate::algebra::Algebra) (that trait is over $\mathbb{F}_p$). //! -//! - [`Tau`] — the coefficient ring $\mathbb{F}_2[\tau]$ as a small homogeneous -//! scalar ([`tau`]); the whole $\tau$-tower is carried here rather than threaded -//! through the resolution engine. -//! - [`MotivicMilnorAlgebra`] — $A_C$, a free $\mathbb{F}_2[\tau]$-module on the -//! Milnor basis, with the Kong–Lin product ([`milnor`]). This is the product -//! **engine**; it is deliberately not an [`Algebra`](crate::algebra::Algebra) -//! (that trait is over $\mathbb{F}_p$). +//! There is no representation of the coefficient ring $\mathbb{F}_2[\tau]$ here, because +//! nothing needs one: every element in sight is bidegree-homogeneous and $\tau$ has weight +//! $-1$, so a coefficient is pinned by the weights of the terms it sits between. See +//! [`milnor::Grading::tau_exponent`], which is how one is recovered when a caller wants it. //! //! The mod-$\tau$ reduction $A_C/\tau$ — the connected finite-type //! $\mathbb{F}_2$-algebra that the existing resolution engine resolves to yield the //! algebraic Novikov $E_2$ — is presented as an [`Algebra`](crate::algebra::Algebra) //! in a follow-up on top of this engine. -pub mod tau; -pub use tau::Tau; - pub mod milnor; pub use milnor::MotivicMilnorAlgebra; diff --git a/ext/crates/algebra/src/algebra/motivic/tau.rs b/ext/crates/algebra/src/algebra/motivic/tau.rs deleted file mode 100644 index 54a5b12aa1..0000000000 --- a/ext/crates/algebra/src/algebra/motivic/tau.rs +++ /dev/null @@ -1,118 +0,0 @@ -//! The coefficient ring $\mathbb{F}_2[\tau]$, represented by $\tau$-valuations. -//! -//! Every computation over the C-motivic Steenrod algebra that we perform is -//! **homogeneous** — differentials, products, and coproducts all preserve the -//! motivic weight — and a homogeneous element of $\mathbb{F}_2[\tau]$ in a fixed -//! bidegree is exactly $\tau^k$ for a single $k$ (or $0$). So we never need a -//! polynomial: a coefficient is one integer, its $\tau$-valuation, with a -//! sentinel for zero. That is [`Tau`]. -//! -//! This is a small, self-contained scalar type for the motivic layer only. It is -//! deliberately **not** an [`fp`] replacement threaded through the resolution -//! engine: the engine stays $\mathbb{F}_p$-only (see the `motivic` module docs), -//! and `Tau` is used solely by the motivic algebra and the Phase 2 lift. - -/// A homogeneous $\mathbb{F}_2[\tau]$ coefficient: either $0$ or $\tau^k$ for a -/// unique $k \ge 0$, stored as the $\tau$-valuation with a zero sentinel. -/// -/// Over $\mathbb{F}_2$ a nonzero homogeneous coefficient has unit $\mathbb{F}_2$ -/// part, so the valuation is the whole story: `None` is $0$, `Some(k)` is -/// $\tau^k$. Multiplication adds valuations; addition of two coefficients in the -/// *same* bidegree can only ever combine equal powers (homogeneity), so it -/// cancels mod 2. -#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] -pub struct Tau(Option); - -impl Tau { - /// The unit $\tau^0 = 1$. - pub const ONE: Self = Self(Some(0)); - /// The zero coefficient. - pub const ZERO: Self = Self(None); - - /// The zero coefficient (method form, for symmetry with [`Tau::one`]). - pub fn zero() -> Self { - Self::ZERO - } - - /// The unit $1 = \tau^0$. - pub fn one() -> Self { - Self::ONE - } - - /// The coefficient $\tau^k$. - pub fn power(k: u32) -> Self { - Self(Some(k)) - } - - /// Whether this coefficient is zero. - pub fn is_zero(self) -> bool { - self.0.is_none() - } - - /// The $\tau$-valuation: `Some(k)` if this is $\tau^k$, `None` if zero. - pub fn valuation(self) -> Option { - self.0 - } -} - -impl std::ops::Mul for Tau { - type Output = Self; - - /// Product: $\tau^i \cdot \tau^j = \tau^{i+j}$; zero absorbs. - // Multiplying τ-powers genuinely *adds* their valuations, so the `+` here is correct. - #[allow(clippy::suspicious_arithmetic_impl)] - fn mul(self, other: Self) -> Self { - match (self.0, other.0) { - (Some(i), Some(j)) => Self(Some(i + j)), - _ => Self::ZERO, - } - } -} - -impl std::ops::Add for Tau { - type Output = Self; - - /// Sum of two homogeneous coefficients in one bidegree. By homogeneity the - /// two nonzero powers must be equal, so $\tau^k + \tau^k = 0$ over - /// $\mathbb{F}_2$; a zero summand is absorbed. Adding *distinct* nonzero - /// powers is inhomogeneous and cannot arise (checked in debug builds). - fn add(self, other: Self) -> Self { - match (self.0, other.0) { - (None, x) | (x, None) => Self(x), - (Some(i), Some(j)) => { - debug_assert_eq!(i, j, "adding inhomogeneous tau powers τ^{i} + τ^{j}"); - Self::ZERO - } - } - } -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn test_zero_one() { - assert!(Tau::zero().is_zero()); - assert!(!Tau::one().is_zero()); - assert_eq!(Tau::one().valuation(), Some(0)); - assert_eq!(Tau::power(3).valuation(), Some(3)); - assert_eq!(Tau::zero().valuation(), None); - } - - #[test] - fn test_mul_adds_valuations() { - assert_eq!(Tau::power(2) * Tau::power(3), Tau::power(5)); - assert_eq!(Tau::one() * Tau::power(4), Tau::power(4)); - assert_eq!(Tau::zero() * Tau::power(4), Tau::zero()); - assert_eq!(Tau::power(4) * Tau::zero(), Tau::zero()); - } - - #[test] - fn test_add_cancels_equal_powers() { - assert_eq!(Tau::power(2) + Tau::power(2), Tau::zero()); - assert_eq!(Tau::power(2) + Tau::zero(), Tau::power(2)); - assert_eq!(Tau::zero() + Tau::power(5), Tau::power(5)); - assert_eq!(Tau::zero() + Tau::zero(), Tau::zero()); - } -} From 2324a95c4b164da00deeb79d34ba63351f220de1 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sat, 29 Aug 2026 18:31:34 -0400 Subject: [PATCH 07/15] motivic: give the A_C reading of a monomial its own type MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A_C and A_** are indexed by the same (E, R) data, so `Monomial` could not say which of the two it meant — and `DualElement` was the return type of both `dual_mul`, where it is an element of the dual algebra, and `multiply`, where it is not. The `Grading` enum existed only to supply, by hand and at every call, the fact that the type had lost. `Dual` carries it instead. `Bigraded::bidegree` reports each type's own weight, and a single `tau_exponent` reads whichever the values it is given belong to, so the convention is no longer something a caller chooses — it follows from what they are holding. Mixing the two in one call is a type error; reinterpreting a value from one side as the other now requires writing `Dual(..)` or `.0`, which is where such a conversion should be visible. The `debug_assert` on the exponent stays as the backstop for that deliberate case. Co-Authored-By: Claude Opus 5 --- ext/crates/algebra/benches/motivic.rs | 13 +- .../algebra/src/algebra/motivic/milnor.rs | 259 ++++++++++-------- ext/crates/algebra/src/algebra/motivic/mod.rs | 4 +- 3 files changed, 150 insertions(+), 126 deletions(-) diff --git a/ext/crates/algebra/benches/motivic.rs b/ext/crates/algebra/benches/motivic.rs index 50a39a8430..2ca72f6f26 100644 --- a/ext/crates/algebra/benches/motivic.rs +++ b/ext/crates/algebra/benches/motivic.rs @@ -15,7 +15,7 @@ use algebra::{ MotivicMilnorAlgebra, - motivic::milnor::{Monomial, enum_basis, multiply_closed}, + motivic::milnor::{Dual, Monomial, enum_basis, multiply_closed}, }; use criterion::{ BenchmarkGroup, BenchmarkId, Criterion, Throughput, criterion_group, criterion_main, @@ -25,13 +25,18 @@ use pprof::criterion::{Output, PProfProfiler}; /// `Q(E)P(R)` from the `Q` indices and the ξ exponents, in the paper's indexing where `R[0]` /// belongs to ξ_0 = 1 and is skipped. -fn elt(q: &[u32], xi: &[u32]) -> Monomial { +fn elt(q: &[u32], xi: &[u32]) -> Dual { let mut r = vec![0]; r.extend_from_slice(xi); - Monomial::from_paper(q.iter().map(|i| 1 << i).sum(), &r).unwrap() + Dual(Monomial::from_paper(q.iter().map(|i| 1 << i).sum(), &r).unwrap()) } -fn bench_product(g: &mut BenchmarkGroup, name: &str, a: Monomial, b: Monomial) { +fn bench_product( + g: &mut BenchmarkGroup, + name: &str, + a: Dual, + b: Dual, +) { g.bench_function(name, |bench| { bench.iter(|| std::hint::black_box(multiply_closed(a, b))); }); diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index 8a66a3efa5..ecb5edb55e 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -31,14 +31,14 @@ //! //! Coefficients are never stored. Everything here is homogeneous, so each one is a single power //! of $\tau$, and $\tau$ has weight $-1$ — which pins the exponent to the weights of the terms -//! it sits between ([`Grading::tau_exponent`]). This is the $A_C$ product **engine**, over +//! it sits between ([`tau_exponent`]). This is the $A_C$ product **engine**, over //! $\mathbb{F}_2[\tau]$; the mod-$\tau$ reduction $A_C/\tau$ is presented to the resolution //! engine as an ordinary $\mathbb{F}_2$-algebra in a follow-up layer built on top of this engine. //! //! Weight convention: the motivic weight of an algebra basis element is the //! *negative* of the weight of the dual monomial it pairs with, so that products -//! are weight-homogeneous with $\tau$ (weight $-1$) absorbing the difference. [`Grading`] names -//! the two sides. +//! are weight-homogeneous with $\tau$ (weight $-1$) absorbing the difference. [`Dual`] marks the +//! $A_C$ side, so the two cannot be confused. use std::{ cell::RefCell, @@ -118,23 +118,6 @@ impl Monomial { pub fn one() -> Self { Self::default() } - - /// The C-motivic (prime 2) bidegree `(t, w)`: `|τ_i| = (2^{i+1}-1, 2^i-1)` and - /// `|ξ_{i+1}| = (2^{i+2}-2, 2^{i+1}-1)`. - /// - /// [`MotivicMilnorAlgebra::bidegree`] negates `w` for the algebra's own weight convention. - pub fn bidegree(&self) -> (i32, i32) { - let (mut t, mut w) = (0i32, 0i32); - for i in BitflagIterator::set_bit_iterator(self.q_part as u64) { - t += (1 << (i + 1)) - 1; - w += (1 << i) - 1; - } - for (i, r) in self.p_part.iter().enumerate() { - t += r as i32 * ((1 << (i + 2)) - 2); - w += r as i32 * ((1 << (i + 1)) - 1); - } - (t, w) - } } /// $\Sigma(R) = \sum_i r_i$. @@ -286,7 +269,7 @@ fn rewrite_tau_dfs( /// /// Coefficients are not stored, because the bigrading already determines them: everything we /// compute is homogeneous, $\tau$ has weight $-1$, so the coefficient of a term is forced to be -/// $\tau^{\text{(input weight)} - \text{(term weight)}}$ — see [`Grading::tau_exponent`]. What is left is +/// $\tau^{\text{(input weight)} - \text{(term weight)}}$ — see [`tau_exponent`]. What is left is /// a set of keys taken mod 2, and adding a term toggles it. /// /// A sorted `Vec` rather than a `BTreeMap`: these combinations are small — a product of two basis @@ -368,54 +351,73 @@ impl<'a, K> IntoIterator for &'a SparseSum { } } -/// Which of the two dual weight conventions a computation is expressed in. +/// The basis element of $A_C$ dual to a monomial of $A_{**}$. /// -/// $A_C$ and its dual $A_{**}$ weight the same [`Monomial`] by negatives of each other: the -/// monomial's own weight grades $A_{**}$, and the basis element of $A_C$ it pairs with takes the -/// negative (see the module docs). Naming the convention keeps that sign in one place — the -/// formula for a $\tau$ exponent is the same on both sides once each is asked for its own weight. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum Grading { - /// $A_{**}$, weighted by the monomial itself — the convention of [`dual_mul`]. - Dual, - /// $A_C$, weighted by the basis element the monomial is dual to — the convention of - /// [`multiply`] and [`product_indexed`](MotivicMilnorAlgebra::product_indexed). - Steenrod, +/// $A_C$ and $A_{**}$ are indexed by the same $(E, R)$ data, so a bare [`Monomial`] cannot say +/// which of the two it means — and the distinction is not cosmetic: the two weight the same data +/// by negatives of each other, so a $\tau$ exponent read in the wrong one comes out with the +/// wrong sign. Wrapping the $A_C$ reading makes that a type error instead. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub struct Dual(pub T); + +/// A monomial carrying a motivic bidegree, in whichever of the two dual conventions its type +/// stands for. +pub trait Bigraded: Copy + std::fmt::Debug { + /// The `(topological degree, motivic weight)`. + fn bidegree(self) -> (i32, i32); } -impl Grading { - /// The motivic weight of `m` in this convention. - pub fn weight(self, m: Monomial) -> i32 { - let w = m.bidegree().1; - match self { - Self::Dual => w, - Self::Steenrod => -w, +/// The C-motivic (prime 2) bidegree of a monomial of $A_{**}$: `|τ_i| = (2^{i+1}-1, 2^i-1)` and +/// `|ξ_{i+1}| = (2^{i+2}-2, 2^{i+1}-1)`. +impl Bigraded for Monomial { + fn bidegree(self) -> (i32, i32) { + let (mut t, mut w) = (0i32, 0i32); + for i in BitflagIterator::set_bit_iterator(self.q_part as u64) { + t += (1 << (i + 1)) - 1; + w += (1 << i) - 1; + } + for (i, r) in self.p_part.iter().enumerate() { + t += r as i32 * ((1 << (i + 2)) - 2); + w += r as i32 * ((1 << (i + 1)) - 1); } + (t, w) } +} - /// The exponent of $\tau$ on the term `z` of a product `a * b`. - /// - /// The bigrading forces it, which is why products store no coefficients. $\tau$ has weight - /// $-1$ and topological degree $0$, so a weight-homogeneous product can only carry - /// $\tau^{w_z - w_a - w_b}$ on `z`. - /// - /// # Panics - /// - /// In debug builds, if the weights imply a negative power — i.e. `z` cannot occur in `a * b`. - pub fn tau_exponent(self, a: Monomial, b: Monomial, z: Monomial) -> u32 { - let k = self.weight(z) - self.weight(a) - self.weight(b); - debug_assert!( - k >= 0, - "negative tau exponent: {a:?} * {b:?} → {z:?} ({self:?})" - ); - k as u32 +/// The same data read in $A_C$, whose weight is the negative of the monomial's. +impl Bigraded for Dual { + fn bidegree(self) -> (i32, i32) { + let (t, w) = self.0.bidegree(); + (t, -w) } } +/// The exponent of $\tau$ on the term `z` of a product `a * b` — the coefficient that products +/// here leave unstored. +/// +/// The bigrading forces it, which is why nothing stores one. $\tau$ has bidegree $(0, -1)$, so +/// topological degree is additive and the whole $\tau$ power lands in the weight: a +/// weight-homogeneous product can only carry $\tau^{w_z - w_a - w_b}$ on `z`. The formula is the +/// same on both sides of the duality; which weights it reads is settled by whether `W` is +/// [`Monomial`] ($A_{**}$) or [`Dual`] ($A_C$). +/// +/// # Panics +/// +/// In debug builds, if the weights imply a negative power — i.e. `z` cannot occur in `a * b`. +pub fn tau_exponent(a: W, b: W, z: W) -> u32 { + let k = z.bidegree().1 - a.bidegree().1 - b.bidegree().1; + debug_assert!(k >= 0, "negative tau exponent: {a:?} * {b:?} → {z:?}"); + k as u32 +} + /// An element of the dual algebra $A_{**}$: a mod-2 set of [`Monomial`]s, whose -/// $\mathbb{F}_2[\tau]$ coefficients are implicit in the weight ([`Grading::tau_exponent`]). +/// $\mathbb{F}_2[\tau]$ coefficients are implicit in the weight ([`tau_exponent`]). pub type DualElement = SparseSum; +/// An element of $A_C$: a mod-2 set of basis elements, whose $\mathbb{F}_2[\tau]$ coefficients +/// are likewise implicit ([`tau_exponent`]). +pub type SteenrodElement = SparseSum>; + /// Elementwise sum of two exponent sequences, i.e. the product of the two $\xi$ monomials. fn ppart_add(a: PPart, b: PPart) -> PPart { let mut out = PPart::zero(); @@ -447,7 +449,7 @@ fn mul_monomials(m1: Monomial, m2: Monomial, acc: &mut DualElement) { let z = Monomial::new(term.e_mask, r); // The τ power `rewrite_tau` computed and the one the weight dictates are the same // number; this is the invariant that lets the coefficient go unstored. - debug_assert_eq!(term.tau_pow, Grading::Dual.tau_exponent(m1, m2, z)); + debug_assert_eq!(term.tau_pow, tau_exponent(m1, m2, z)); acc.add_term(z); } } @@ -633,7 +635,7 @@ pub fn antipode(elt: &DualElement) -> DualElement { // --------------------------------------------------------------------------- /// All basis monomials `τ(E)ξ(R)` of a given topological degree `target`. -pub fn enum_basis(target: i32) -> Vec { +pub fn enum_basis(target: i32) -> Vec> { // Generators available up to `target`: τ_i of degree 2^{i+1}-1, ξ_j (j≥1) of degree // 2^{j+1}-2. Represented as (is_tau, index, degree). let mut gens: Vec<(bool, usize, i32)> = Vec::new(); @@ -667,12 +669,14 @@ fn enum_basis_dfs( rem: i32, e: &mut u32, r: &mut [u32], - out: &mut Vec, + out: &mut Vec>, ) { if idx == gens.len() { if rem == 0 { // `r` is paper-indexed: `r[0]` is the (always zero) xi_0 slot. - out.push(Monomial::from_paper(*e, r).expect("basis exponent out of packing range")); + out.push(Dual( + Monomial::from_paper(*e, r).expect("basis exponent out of packing range"), + )); } return; } @@ -698,17 +702,18 @@ fn enum_basis_dfs( /// The product `a · b` in the C-motivic Steenrod algebra `A_C`, where `a`, `b` are Milnor basis /// elements `Q(E)P(R)` given as the monomials `(E, R)` they are dual to. The result is the set of -/// basis monomials occurring, their `𝔽₂[τ]` coefficients being implicit ([`Grading::tau_exponent`]). +/// basis monomials occurring, their `𝔽₂[τ]` coefficients being implicit ([`tau_exponent`]). /// /// Computed by duality: the coefficient of `z` in `a · b` is the coefficient of /// `mon(a) ⊗ mon(b)` in `ψ(mon(z))`, summed over the basis `z` of the appropriate /// topological degree. -pub fn multiply(a: Monomial, b: Monomial) -> DualElement { - let key = (a, b); +pub fn multiply(a: Dual, b: Dual) -> SteenrodElement { + // ψ is a statement about A_**, so the pairing is against the monomials themselves. + let key = (a.0, b.0); let t = a.bidegree().0 + b.bidegree().0; - let mut out = DualElement::new(); + let mut out = SteenrodElement::new(); for z in enum_basis(t) { - if coproduct_monomial(&z).contains(&key) { + if coproduct_monomial(&z.0).contains(&key) { out.add_term(z); } } @@ -1047,7 +1052,10 @@ impl ClosedY<'_> { /// The product `a · b` in `A_C` via Kong–Lin Theorem 5.1 (ρ = 0). Same contract /// as [`multiply`] (the duality oracle it is validated against). -pub fn multiply_closed(a: Monomial, b: Monomial) -> DualElement { +pub fn multiply_closed(a: Dual, b: Dual) -> SteenrodElement { + // Theorem 5.1 is combinatorics on the (E, R) exponents; the A_C reading is re-applied to the + // output monomials at the end. + let (a, b) = (a.0, b.0); let (e1_mask, e2_mask) = (a.q_part, b.q_part); // The recursion below follows Kong–Lin's indexing, in which `R[0]` is the xi_0 slot, so the // exponents are unpacked into paper-indexed working vectors here. @@ -1080,7 +1088,7 @@ pub fn multiply_closed(a: Monomial, b: Monomial) -> DualElement { let mut xcols: Vec<&[u32]> = Vec::with_capacity(l); let mut acc = Acc::zero(); ctx.enum_x(0, &mut xcols, &mut acc, &mut out); - out + out.into_iter().map(Dual).collect() } // --------------------------------------------------------------------------- @@ -1093,7 +1101,7 @@ pub fn multiply_closed(a: Monomial, b: Monomial) -> DualElement { /// /// This is the $A_C$ product **engine** over $\mathbb{F}_2[\tau]$: the dual-based product is /// exposed through [`product_indexed`](MotivicMilnorAlgebra::product_indexed), whose -/// $\tau$ coefficients are implicit in the weights ([`Grading::tau_exponent`]). It is +/// $\tau$ coefficients are implicit in the weights ([`tau_exponent`]). It is /// deliberately not an [`Algebra`](crate::algebra::Algebra) /// implementation, because that trait is over $\mathbb{F}_p$; the mod-$\tau$ reduction, which *is* /// such an algebra, is a follow-up layer on top of this engine, and the honest @@ -1107,7 +1115,7 @@ pub fn multiply_closed(a: Monomial, b: Monomial) -> DualElement { pub struct MotivicMilnorAlgebra { /// `basis[t]` is the $\mathbb{F}_2[\tau]$-basis in topological degree `t`, sorted for stable /// indexing. - basis: OnceVec>, + basis: OnceVec>>, /// Memoized basis-element products, one dense [`ProductBlock`] per pair of /// topological degrees `(t1, t2)`. The duality product is expensive and a /// resolution asks for the same structure constants repeatedly, so we cache @@ -1165,26 +1173,25 @@ impl MotivicMilnorAlgebra { } /// The `idx`-th basis monomial `(E, R)` in degree `degree`. - pub fn basis_element(&self, degree: i32, idx: usize) -> &Monomial { + pub fn basis_element(&self, degree: i32, idx: usize) -> &Dual { &self.basis[degree as usize][idx] } /// The index of a basis monomial in its degree, if present. - pub fn index_of(&self, degree: i32, elt: &Monomial) -> Option { + pub fn index_of(&self, degree: i32, elt: &Dual) -> Option { self.basis[degree as usize].binary_search(elt).ok() } /// The `(topological degree, motivic weight)` of a basis element. pub fn bidegree(&self, degree: i32, idx: usize) -> (i32, i32) { - let m = self.basis[degree as usize][idx]; - (m.bidegree().0, Grading::Steenrod.weight(m)) + self.basis[degree as usize][idx].bidegree() } /// The exponent of $\tau$ on basis element `idx` of degree `t1 + t2` in the product of the /// elements at `(t1, idx1)` and `(t2, idx2)` — the coefficient [`Self::product_indexed`] - /// leaves implicit. See [`Grading::tau_exponent`]. + /// leaves implicit. See [`tau_exponent`]. pub fn tau_exponent(&self, t1: i32, idx1: usize, t2: i32, idx2: usize, idx: usize) -> u32 { - Grading::Steenrod.tau_exponent( + tau_exponent( self.basis[t1 as usize][idx1], self.basis[t2 as usize][idx2], self.basis[(t1 + t2) as usize][idx], @@ -1228,12 +1235,12 @@ impl MotivicMilnorAlgebra { } self.compute_basis(degree); // The basis is stored trimmed, so trimming the parsed vector makes this an exact match. - self.index_of(degree, &m).map(|idx| (degree, idx)) + self.index_of(degree, &Dual(m)).map(|idx| (degree, idx)) } /// A display string for a basis element (`Q_i … P(R)`). pub fn basis_element_to_string(&self, degree: i32, idx: usize) -> String { - let m = self.basis_element(degree, idx); + let Dual(m) = self.basis_element(degree, idx); let mut parts = Vec::new(); for i in fp::prime::iter::BitflagIterator::set_bit_iterator(m.q_part as u64) { parts.push(format!("Q_{i}")); @@ -1408,6 +1415,16 @@ mod tests { mon(0, &pp_to_paper(pp)) } + /// The $A_C$ basis element dual to [`mon`]'s monomial. + fn dmon(q_part: u32, r: &[u32]) -> Dual { + Dual(mon(q_part, r)) + } + + /// The $A_C$ basis element dual to [`paper`]'s monomial. + fn dpaper(pp: &[u32]) -> Dual { + Dual(paper(pp)) + } + #[test] fn test_algebra_basis_and_multiply() { let alg = MotivicMilnorAlgebra::new(); @@ -1415,25 +1432,25 @@ mod tests { // Degree 0 is the unit; degrees 1 and 2 are 1-dimensional (Q_0 and P(ξ_1)). assert_eq!(alg.dimension(0), 1); - assert_eq!(alg.basis_element(0, 0), &mon(0u32, &[])); + assert_eq!(alg.basis_element(0, 0), &dmon(0u32, &[])); assert_eq!(alg.dimension(1), 1); - assert_eq!(alg.basis_element(1, 0), &mon(0b1u32, &[])); // Q_0 + assert_eq!(alg.basis_element(1, 0), &dmon(0b1u32, &[])); // Q_0 assert_eq!(alg.dimension(2), 1); - assert_eq!(alg.basis_element(2, 0), &mon(0u32, &[0, 1])); // P(ξ_1) + assert_eq!(alg.basis_element(2, 0), &dmon(0u32, &[0, 1])); // P(ξ_1) // bidegree: Q_0 is (1, 0), P(ξ_1) is (2, -1) in this presentation. assert_eq!(alg.bidegree(1, 0), (1, 0)); assert_eq!(alg.bidegree(2, 0), (2, -1)); // Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1), reconstructed from indices. - let terms: DualElement = alg + let terms: SteenrodElement = alg .product_indexed(1, 0, 2, 0) .into_iter() .map(|idx| *alg.basis_element(3, idx)) .collect(); assert_eq!( terms, - DualElement::from([mon(0b10, &[]), mon(0b1, &[0, 1])]) + SteenrodElement::from([dmon(0b10, &[]), dmon(0b1, &[0, 1])]) ); } @@ -1448,7 +1465,7 @@ mod tests { let a = *alg.basis_element(t1, idx1); let b = *alg.basis_element(t2, idx2); // Indexed product agrees with the raw monomial product. - let indexed: DualElement = alg + let indexed: SteenrodElement = alg .product_indexed(t1, idx1, t2, idx2) .into_iter() .map(|idx| *alg.basis_element(t1 + t2, idx)) @@ -1462,7 +1479,7 @@ mod tests { let z = *alg.basis_element(t1 + t2, idx); assert_eq!( alg.tau_exponent(t1, idx1, t2, idx2, idx), - Grading::Steenrod.tau_exponent(a, b, z) + tau_exponent(a, b, z) ); } } @@ -1479,26 +1496,26 @@ mod tests { dual_mul(&tau_gen(0), &tau_gen(0)), DualElement::from([mon(0, &[0, 1])]) ); - assert_eq!(Grading::Dual.tau_exponent(t0, t0, mon(0, &[0, 1])), 1); + assert_eq!(tau_exponent(t0, t0, mon(0, &[0, 1])), 1); // ξ_1^2 is just the monomial ξ_1^2, with no τ. assert_eq!( dual_mul(&xi_gen(1), &xi_gen(1)), DualElement::from([mon(0, &[0, 2])]) ); let xi1 = mon(0, &[0, 1]); - assert_eq!(Grading::Dual.tau_exponent(xi1, xi1, mon(0, &[0, 2])), 0); + assert_eq!(tau_exponent(xi1, xi1, mon(0, &[0, 2])), 0); // Distinct τ's commute and stay square-free: τ_0 τ_1, again τ-free. assert_eq!( dual_mul(&tau_gen(0), &tau_gen(1)), DualElement::from([mon(0b11, &[])]) ); - assert_eq!(Grading::Dual.tau_exponent(t0, t1, mon(0b11, &[])), 0); + assert_eq!(tau_exponent(t0, t1, mon(0b11, &[])), 0); // τ_1^2 = τ ξ_2. assert_eq!( dual_mul(&tau_gen(1), &tau_gen(1)), DualElement::from([mon(0, &[0, 0, 1])]) ); - assert_eq!(Grading::Dual.tau_exponent(t1, t1, mon(0, &[0, 0, 1])), 1); + assert_eq!(tau_exponent(t1, t1, mon(0, &[0, 0, 1])), 1); // Multiplication by the unit is the identity. assert_eq!(dual_mul(&dual_one(), &tau_gen(2)), tau_gen(2)); } @@ -1516,16 +1533,16 @@ mod tests { }; let motivic = |p1: &[u32], p2: &[u32]| -> Vec> { - let mut out: Vec> = multiply(paper(p1), paper(p2)) + let mut out: Vec> = multiply(dpaper(p1), dpaper(p2)) .into_iter() - .filter(|m| m.q_part == 0) // keep the pure-ξ outputs - .map(|m| { + .filter(|d| d.0.q_part == 0) // keep the pure-ξ outputs + .map(|d| { assert_eq!( - Grading::Steenrod.tau_exponent(paper(p1), paper(p2), m), + tau_exponent(dpaper(p1), dpaper(p2), d), 0, "pure-ξ output carried a τ power" ); - m.p_part.iter().collect() + d.0.p_part.iter().collect() }) .collect(); out.sort(); @@ -1569,8 +1586,8 @@ mod tests { let xi_pps = |t: i32| -> Vec> { enum_basis(t) .into_iter() - .filter(|m| m.q_part == 0) - .map(|m| m.p_part.iter().collect()) + .filter(|d| d.0.q_part == 0) + .map(|d| d.0.p_part.iter().collect()) .collect() }; let a_coeff = |mu_pp: &[u32], w_pp: &[u32]| -> u32 { @@ -1594,10 +1611,10 @@ mod tests { (vec![2, 1], vec![2]), ] { // LHS: motivic conjugate product (E=0), converted to std via χ. - let conj: Vec> = multiply(paper(&r1), paper(&r2)) + let conj: Vec> = multiply(dpaper(&r1), dpaper(&r2)) .into_iter() - .filter(|m| m.q_part == 0) - .map(|m| m.p_part.iter().collect()) + .filter(|d| d.0.q_part == 0) + .map(|d| d.0.p_part.iter().collect()) .collect(); let mut lhs: BTreeMap, u32> = BTreeMap::new(); for w in xi_pps(mdeg(&r1) + mdeg(&r2)) { @@ -1633,8 +1650,8 @@ mod tests { // Associativity (a·b)·c = a·(b·c) is convention-independent and a strong global check. // Extend the basis×basis product to element×basis, F_2[τ]-linearly (`add_term` cancels // mod 2, and the τ coefficients agree termwise because both sides are homogeneous). - fn mul_elt_basis(x: &DualElement, c: &Monomial) -> DualElement { - let mut out = DualElement::new(); + fn mul_elt_basis(x: &SteenrodElement, c: &Dual) -> SteenrodElement { + let mut out = SteenrodElement::new(); for &z in x { for w in multiply(z, *c) { out.add_term(w); @@ -1643,7 +1660,7 @@ mod tests { out } - let basis: Vec = (0..=5).flat_map(enum_basis).collect(); + let basis: Vec> = (0..=5).flat_map(enum_basis).collect(); for a in &basis { for b in &basis { let ab = multiply(*a, *b); @@ -1651,7 +1668,7 @@ mod tests { let lhs = mul_elt_basis(&ab, c); let bc = multiply(*b, *c); // a · (b·c): multiply basis `a` on the left of each term of bc. - let mut rhs = DualElement::new(); + let mut rhs = SteenrodElement::new(); for &z in &bc { for w in multiply(*a, z) { rhs.add_term(w); @@ -1667,13 +1684,13 @@ mod tests { fn test_product_q0_p_xi1() { // The case my reading of Kong–Lin Theorem 5.1 got wrong. By duality: // Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1) (both with coefficient τ^0). - let q0 = mon(0b1, &[]); - let p_xi1 = mon(0, &[0, 1]); + let q0 = dmon(0b1, &[]); + let p_xi1 = dmon(0, &[0, 1]); assert_eq!( multiply(q0, p_xi1), - DualElement::from([ - mon(0b10, &[]), // Q_1 - mon(0b1, &[0, 1]), // Q_0 P(ξ_1) + SteenrodElement::from([ + dmon(0b10, &[]), // Q_1 + dmon(0b1, &[0, 1]), // Q_0 P(ξ_1) ]) ); } @@ -1681,15 +1698,15 @@ mod tests { #[test] fn test_product_unit_and_squares() { // 1 · x = x. - let x = mon(0b101, &[0, 2]); - assert_eq!(multiply(mon(0, &[]), x), DualElement::from([x])); + let x = dmon(0b101, &[0, 2]); + assert_eq!(multiply(dmon(0, &[]), x), SteenrodElement::from([x])); // Q_i^2 = 0. (Q_0 = Sq^1 is the motivic Bockstein; P(ξ_1) = Sq^2 does NOT square to // zero — its square is τ Q_0 Q_1 + …, a genuine motivic feature.) for i in 0..3 { assert!( multiply( - Monomial::new(1 << i, PPart::zero()), - Monomial::new(1 << i, PPart::zero()) + Dual(Monomial::new(1 << i, PPart::zero())), + Dual(Monomial::new(1 << i, PPart::zero())) ) .is_empty(), "Q_{i}^2 ≠ 0" @@ -1705,14 +1722,14 @@ mod tests { // consequences it must have: the exponent is a genuine (non-negative) power, the // topological degree is additive, and τ really does enter — i.e. the τ_i^2 = τξ_{i+1} // relation fires rather than the whole thing being vacuously τ-free. - let basis: Vec = (0..=6).flat_map(enum_basis).collect(); + let basis: Vec> = (0..=6).flat_map(enum_basis).collect(); let mut saw_tau = false; for a in &basis { for b in &basis { for z in multiply(*a, *b) { assert_eq!(z.bidegree().0, a.bidegree().0 + b.bidegree().0); // Panics in debug if the weights imply a negative power. - saw_tau |= Grading::Steenrod.tau_exponent(*a, *b, z) > 0; + saw_tau |= tau_exponent(*a, *b, z) > 0; } } } @@ -1858,8 +1875,8 @@ mod tests { fn test_closed_form_small_cases() { // The case a naive reading of Theorem 5.1 gets wrong (the ξ₀ = 1 index // absorption): Q_0 · P(ξ_1) = Q_1 + Q_0 P(ξ_1). - let q0 = mon(0b1, &[]); - let p_xi1 = mon(0, &[0, 1]); + let q0 = dmon(0b1, &[]); + let p_xi1 = dmon(0, &[0, 1]); assert_eq!(multiply_closed(q0, p_xi1), multiply(q0, p_xi1)); // The τ-generating case: P(ξ_1)² has a τ Q_0 Q_1 term. @@ -1869,12 +1886,12 @@ mod tests { for i in 0..3 { assert_eq!( multiply_closed( - Monomial::new(1 << i, PPart::zero()), - Monomial::new(1 << i, PPart::zero()) + Dual(Monomial::new(1 << i, PPart::zero())), + Dual(Monomial::new(1 << i, PPart::zero())) ), multiply( - Monomial::new(1 << i, PPart::zero()), - Monomial::new(1 << i, PPart::zero()) + Dual(Monomial::new(1 << i, PPart::zero())), + Dual(Monomial::new(1 << i, PPart::zero())) ) ); } @@ -1886,7 +1903,7 @@ mod tests { // oracle on every ordered pair of basis elements with deg(a)+deg(b) ≤ 18, // including the τ-carrying and multi-Q cases. (The bound is kept modest // because the *oracle* is slow; the closed form is not.) - let basis: Vec<(i32, Monomial)> = (0..=14) + let basis: Vec<(i32, Dual)> = (0..=14) .flat_map(|t| enum_basis(t).into_iter().map(move |m| (t, m))) .collect(); for (ta, a) in &basis { diff --git a/ext/crates/algebra/src/algebra/motivic/mod.rs b/ext/crates/algebra/src/algebra/motivic/mod.rs index 17579f36c1..f5891b46c0 100644 --- a/ext/crates/algebra/src/algebra/motivic/mod.rs +++ b/ext/crates/algebra/src/algebra/motivic/mod.rs @@ -12,7 +12,9 @@ //! There is no representation of the coefficient ring $\mathbb{F}_2[\tau]$ here, because //! nothing needs one: every element in sight is bidegree-homogeneous and $\tau$ has weight //! $-1$, so a coefficient is pinned by the weights of the terms it sits between. See -//! [`milnor::Grading::tau_exponent`], which is how one is recovered when a caller wants it. +//! [`milnor::tau_exponent`], which is how one is recovered when a caller wants it. Which of the +//! two dual weight conventions applies is carried by the type — [`milnor::Dual`] marks the $A_C$ +//! reading of a monomial — so the two cannot be mixed up. //! //! The mod-$\tau$ reduction $A_C/\tau$ — the connected finite-type //! $\mathbb{F}_2$-algebra that the existing resolution engine resolves to yield the From 4bbc3f13c641b00975b368666fe90a3190358056 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sat, 29 Aug 2026 22:42:01 -0400 Subject: [PATCH 08/15] motivic: trim the engine after the tau removal MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Docs first: three of them still described the F_2[tau] coefficient that no longer exists, and both module headers had grown into essays that restated each other and the items below them. CLAUDE.md wants one line and the explanation on the item, so the algebra presentation moved to `Monomial`, the duality argument to `multiply` (which also stops calling the closed form future work, since it has been implemented for a while), and the weight convention to `Dual`. `SparseSum::len` and `iter` had no callers — `iter` was a second spelling of the `IntoIterator` impl next to it. `xi_gen(i)` was `xi_pow_elt(i, 1)` written out, `new()` was the derived `Default` written out, and `antipode` iterated set bits by hand where the rest of the file uses `BitflagIterator`. `on_y` bounded a loop over an `[u32; NB]` by `u32::BITS`, which only works because the two constants happen to be equal; `nb_covers_antidiagonals` now pins the bound `NB` actually needs, and its doc names the constants rather than quoting a number that would go stale. `basis_element_from_string` had no test — its round trip lives in the follow-up that consumes it, so it would have shipped unexercised. Brought down. `product_indexed_with` likewise comes from the follow-up, where copying the index list per cache hit showed up as real cost. Co-Authored-By: Claude Opus 5 --- .../algebra/src/algebra/motivic/milnor.rs | 232 +++++++++--------- ext/crates/algebra/src/algebra/motivic/mod.rs | 21 -- 2 files changed, 119 insertions(+), 134 deletions(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index ecb5edb55e..ddd7e287f4 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -1,44 +1,9 @@ -//! The C-motivic (over $\mathbb{C}$, prime 2) dual Steenrod algebra $A_{**}$ and -//! the Steenrod algebra $A_C$ in the Milnor basis, over $\mathbb{F}_2[\tau]$. +//! The C-motivic dual Steenrod algebra $A_{**}$ and Steenrod algebra $A_C$ in the Milnor basis. //! -//! This follows Kong–Lin, *Product formulas for motivic Milnor basis* -//! (arXiv:2411.12890), specialized to the C-motivic point where $\rho = 0$ (the -//! motivic cohomology of a point over $\mathbb{C}$ is $\mathbb{F}_2[\tau]$). We -//! use the *conjugate* generators of that paper, so the coproduct and -//! Milnor-matrix product match Milnor's classical formulas on the polynomial -//! (`ξ`) part. -//! -//! The dual algebra is $A_{**} = \mathbb{F}_2[\tau][\xi_1, \xi_2, \dots] \otimes -//! E(\tau_0, \tau_1, \dots)$ with $\tau_i^2 = \tau \xi_{i+1}$ (the $\rho = 0$ -//! reduction of $\tau_i^2 = \tau\xi_{i+1} + \rho\tau_{i+1}$). Bidegrees: -//! $|\xi_i| = (2(2^i-1), 2^i-1)$, $|\tau_i| = (2^{i+1}-1, 2^i-1)$, -//! $|\tau| = (0, -1)$ in (stem, weight). -//! -//! Basis elements are $Q(E)P(R)$ with $E \in \mathrm{Seq}_1$ (entries in -//! $\{0,1\}$, encoded as a bitmask: bit $i$ is $\tau_i$) and $R \in \mathrm{Seq}$ -//! (the exponent vector of the $\xi_i$). -//! -//! The product in the Steenrod algebra $A_C$ is computed by **duality**: the -//! dual algebra $A_{**}$ is a commutative $\mathbb{F}_2[\tau]$-algebra whose -//! multiplication ([`dual_mul`]) reduces $\tau_i^2 = \tau\xi_{i+1}$ via -//! [`rewrite_tau`] (Kong–Lin Theorem 3.4), and its coproduct $\psi$ -//! ([`coproduct`], Kong–Lin §2.2) is a $\tau$-free algebra map. The product -//! $a \cdot b$ in $A_C$ is then read off from $\psi$: the coefficient of $z$ in -//! $a \cdot b$ is the coefficient of $\mathrm{mon}(a) \otimes \mathrm{mon}(b)$ in -//! $\psi(\mathrm{mon}(z))$ ([`multiply`]). This is a first, correctness-oriented -//! implementation; the closed-form product (Kong–Lin Theorem 5.1) can replace it -//! later for speed, validated against this one. -//! -//! Coefficients are never stored. Everything here is homogeneous, so each one is a single power -//! of $\tau$, and $\tau$ has weight $-1$ — which pins the exponent to the weights of the terms -//! it sits between ([`tau_exponent`]). This is the $A_C$ product **engine**, over -//! $\mathbb{F}_2[\tau]$; the mod-$\tau$ reduction $A_C/\tau$ is presented to the resolution -//! engine as an ordinary $\mathbb{F}_2$-algebra in a follow-up layer built on top of this engine. -//! -//! Weight convention: the motivic weight of an algebra basis element is the -//! *negative* of the weight of the dual monomial it pairs with, so that products -//! are weight-homogeneous with $\tau$ (weight $-1$) absorbing the difference. [`Dual`] marks the -//! $A_C$ side, so the two cannot be confused. +//! Follows Kong–Lin, *Product formulas for motivic Milnor basis* (arXiv:2411.12890), at the +//! C-motivic point where $\rho = 0$ and the motivic cohomology of a point is +//! $\mathbb{F}_2[\tau]$. It uses that paper's *conjugate* generators, so the coproduct and +//! Milnor-matrix product match Milnor's classical formulas on the $\xi$ part. use std::{ cell::RefCell, @@ -65,8 +30,30 @@ fn trim(mut v: Vec) -> Vec { v } +/// The packed $\xi$ exponents of a paper-indexed sequence: drop the $\xi_0 = 1$ slot. +/// +/// Returns `None` if an exponent is too large to pack, which for the degrees this engine reaches +/// means the caller has gone out of range. The sequence may be empty (it arrives trimmed). +fn p_part_from_paper(r: &[u32]) -> Option { + debug_assert_eq!( + r.first().copied().unwrap_or(0), + 0, + "xi_0 = 1 slot must be 0" + ); + PPart::try_from_slice(r.get(1..).unwrap_or_default()) +} + /// A basis monomial $\tau(E)\xi(R)$ of the dual algebra. /// +/// $A_{**} = \mathbb{F}_2[\tau][\xi_1, \xi_2, \dots] \otimes E(\tau_0, \tau_1, \dots)$ with +/// $\tau_i^2 = \tau \xi_{i+1}$ (the $\rho = 0$ reduction of +/// $\tau_i^2 = \tau\xi_{i+1} + \rho\tau_{i+1}$), so $E \in \mathrm{Seq}_1$ has entries in +/// $\{0, 1\}$ and $R \in \mathrm{Seq}$ is unrestricted. Bidegrees, in (stem, weight): +/// $|\xi_i| = (2(2^i-1), 2^i-1)$, $|\tau_i| = (2^{i+1}-1, 2^i-1)$, $|\tau| = (0, -1)$. +/// +/// The same data indexes the $A_C$ basis element $Q(E)P(R)$ that this monomial is dual to; see +/// [`Dual`], which is that reading of it. +/// /// The derived ordering is lexicographic in `(q_part, p_part)`, which is what /// [`MotivicMilnorAlgebra`] sorts and binary-searches its per-degree bases by. #[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash)] @@ -95,15 +82,7 @@ impl Monomial { /// Returns `None` if an exponent is too large to pack, which for the degrees this engine /// reaches means the caller has gone out of range. pub fn from_paper(q_part: u32, r: &[u32]) -> Option { - debug_assert_eq!( - r.first().copied().unwrap_or(0), - 0, - "xi_0 = 1 slot must be 0" - ); - Some(Self::new( - q_part, - PPart::try_from_slice(r.get(1..).unwrap_or_default())?, - )) + Some(Self::new(q_part, p_part_from_paper(r)?)) } /// The $\xi$ exponents in the paper's indexing, with the $\xi_0$ slot restored. @@ -183,12 +162,8 @@ pub fn rewrite_tau(s: &[u32]) -> Vec { let target2 = sigma2(s); let target_sum = sigma(s); // Highest index we must consider: any nonzero e_i or r_i contributes at least 2^i, so - // 2^i <= target2. - let max_idx = if target2 == 0 { - 0 - } else { - (63 - target2.leading_zeros()) as usize - }; + // 2^i <= target2. An empty `s` gives `target2 == 0`, and index 0 still has to be walked. + let max_idx = target2.highest_one().unwrap_or(0) as usize; let mut out = Vec::new(); let mut e = vec![0u32; max_idx + 1]; @@ -265,12 +240,9 @@ fn rewrite_tau_dfs( // The dual algebra A_** (used to compute the product in A_C by duality). // --------------------------------------------------------------------------- -/// A sparse $\mathbb{F}_2[\tau]$-linear combination of `K`s, as a `Vec` sorted by key. +/// A set of `K`s taken mod 2, as a `Vec` sorted by key: adding a term toggles it. /// -/// Coefficients are not stored, because the bigrading already determines them: everything we -/// compute is homogeneous, $\tau$ has weight $-1$, so the coefficient of a term is forced to be -/// $\tau^{\text{(input weight)} - \text{(term weight)}}$ — see [`tau_exponent`]. What is left is -/// a set of keys taken mod 2, and adding a term toggles it. +/// There are no coefficients to store — see [`tau_exponent`] for why. /// /// A sorted `Vec` rather than a `BTreeMap`: these combinations are small — a product of two basis /// elements has a handful of terms — and the keys are `Copy` words, so one contiguous allocation @@ -286,11 +258,6 @@ impl SparseSum { Self { terms: Vec::new() } } - /// The number of non-zero terms. - pub fn len(&self) -> usize { - self.terms.len() - } - /// Whether this is zero. pub fn is_empty(&self) -> bool { self.terms.is_empty() @@ -310,11 +277,6 @@ impl SparseSum { Err(i) => self.terms.insert(i, key), } } - - /// The terms, in key order. - pub fn iter(&self) -> std::slice::Iter<'_, K> { - self.terms.iter() - } } impl FromIterator for SparseSum { @@ -441,10 +403,7 @@ fn mul_monomials(m1: Monomial, m2: Monomial, acc: &mut DualElement) { .collect(); let r12 = ppart_add(m1.p_part, m2.p_part); for term in rewrite_tau(&s) { - // `term.r` is paper-indexed and `term.r[0]` is zero (`c_coeff` vanishes otherwise); it - // is also trimmed, so it can be empty. - let term_r = PPart::try_from_slice(term.r.get(1..).unwrap_or_default()) - .expect("rewrite_tau exponent out of packing range"); + let term_r = p_part_from_paper(&term.r).expect("rewrite_tau exponent out of packing range"); let r = ppart_add(term_r, r12); let z = Monomial::new(term.e_mask, r); // The τ power `rewrite_tau` computed and the one the weight dictates are the same @@ -472,7 +431,7 @@ pub fn tau_gen(i: usize) -> DualElement { /// The generator $\xi_i \in A_{**}$ (for $i \ge 1$). pub fn xi_gen(i: usize) -> DualElement { - DualElement::from([xi_pow_mon(i, 1)]) + xi_pow_elt(i, 1) } /// The unit $1 \in A_{**}$. @@ -484,8 +443,8 @@ pub fn dual_one() -> DualElement { // The coproduct ψ: A_** → A_** ⊗ A_**, an algebra map (Kong–Lin §2.2). // --------------------------------------------------------------------------- -/// An element of $A_{**} \otimes A_{**}$: pairs of [`Monomial`]s → $\mathbb{F}_2[\tau]$ -/// coefficients. Zero coefficients are never stored. +/// An element of $A_{**} \otimes A_{**}$: a mod-2 set of pairs of [`Monomial`]s, coefficients +/// implicit as in [`DualElement`]. pub type TensorElement = SparseSum<(Monomial, Monomial)>; /// The unit $1 \otimes 1$. @@ -613,10 +572,8 @@ pub fn antipode(elt: &DualElement) -> DualElement { let mut out = DualElement::new(); for &mon in elt { let mut acc = dual_one(); - for i in 0..u32::BITS { - if (mon.q_part >> i) & 1 != 0 { - acc = dual_mul(&acc, &chi_tau(i as usize)); - } + for i in BitflagIterator::set_bit_iterator(mon.q_part as u64) { + acc = dual_mul(&acc, &chi_tau(i)); } for (i, r) in mon.p_part.iter().enumerate() { for _ in 0..r { @@ -700,13 +657,18 @@ fn enum_basis_dfs( } } -/// The product `a · b` in the C-motivic Steenrod algebra `A_C`, where `a`, `b` are Milnor basis -/// elements `Q(E)P(R)` given as the monomials `(E, R)` they are dual to. The result is the set of -/// basis monomials occurring, their `𝔽₂[τ]` coefficients being implicit ([`tau_exponent`]). +/// The product `a · b` in $A_C$, as the set of basis elements occurring — their +/// $\mathbb{F}_2[\tau]$ coefficients being implicit ([`tau_exponent`]). /// -/// Computed by duality: the coefficient of `z` in `a · b` is the coefficient of -/// `mon(a) ⊗ mon(b)` in `ψ(mon(z))`, summed over the basis `z` of the appropriate +/// Computed by duality against the coproduct: $A_{**}$ multiplies by reducing +/// $\tau_i^2 = \tau\xi_{i+1}$ ([`rewrite_tau`], Kong–Lin Theorem 3.4), and its coproduct $\psi$ +/// ([`coproduct`], §2.2) is a $\tau$-free algebra map, so the coefficient of `z` in `a · b` is +/// the coefficient of `mon(a) ⊗ mon(b)` in `ψ(mon(z))`, over the basis `z` of the appropriate /// topological degree. +/// +/// This is the correctness oracle, not the production path: it enumerates a whole degree and +/// coproducts each element. [`multiply_closed`] computes the same thing from Kong–Lin +/// Theorem 5.1 and is validated against this. pub fn multiply(a: Dual, b: Dual) -> SteenrodElement { // ψ is a statement about A_**, so the pairing is against the monomials themselves. let key = (a.0, b.0); @@ -753,6 +715,7 @@ struct Columns { } impl Columns { + /// Append one column. fn push(&mut self, col: &[u32]) { if self.offsets.is_empty() { self.offsets.push(0); @@ -761,6 +724,7 @@ impl Columns { self.offsets.push(self.entries.len() as u32); } + /// The columns, in insertion order. fn iter(&self) -> impl Iterator { self.offsets .windows(2) @@ -805,6 +769,7 @@ thread_local! { static COLS_LE: RefCell>> = RefCell::new(FxHashMap::default()); } +/// All columns with `Σᵢ 2ⁱ vᵢ = target`. fn columns_eq_cached(target: u32) -> Rc { COLS_EQ.with(|c| { Rc::clone( @@ -847,9 +812,9 @@ fn col0_x_options(r1: &[u32]) -> Columns { /// Buffer size for row / column / anti-diagonal indices. /// /// Anti-diagonals are indexed by `i + j`, with the row index `i` bounded by a column's length -/// (at most `log2(MAX_DEGREE) + 1`) and the column index `j` by [`PPart::MAX_LEN`] + 1. That -/// caps the sum near 21, and this is the next power of two; [`Acc`] holds three of these, so the -/// size is paid on every matrix. +/// (at most `PPart::MAX_DEGREE.ilog2() + 1`) and the column index `j` by [`PPart::MAX_LEN`] + 1. +/// This is the next power of two above that sum, checked in `nb_covers_antidiagonals`; [`Acc`] +/// holds three arrays of this size, so it is paid on every matrix. const NB: usize = 32; /// Immutable context for the closed-form product recursion (Theorem 5.1, ρ = 0). @@ -874,6 +839,7 @@ struct Acc { } impl Acc { + /// Empty accumulator, before any column has been placed. fn zero() -> Self { Self { rows: [0; NB], @@ -1037,8 +1003,8 @@ impl ClosedY<'_> { } // E_out = E₂ + T(Y) must be square-free (Seq₁). Q indices fit in a u32. let mut out_e_mask = 0u32; - for i in 0..u32::BITS as usize { - let val = ((self.e2_mask >> i) & 1) + acc.sum[i]; + for (i, &t) in acc.sum.iter().enumerate() { + let val = ((self.e2_mask >> i) & 1) + t; if val > 1 { return; } @@ -1099,18 +1065,14 @@ pub fn multiply_closed(a: Dual, b: Dual) -> SteenrodElement /// $\mathbb{F}_2[\tau]$-module on the Milnor basis $\{Q(E)P(R)\}$ with lazy /// per-(topological-)degree basis indexing. /// -/// This is the $A_C$ product **engine** over $\mathbb{F}_2[\tau]$: the dual-based product is -/// exposed through [`product_indexed`](MotivicMilnorAlgebra::product_indexed), whose -/// $\tau$ coefficients are implicit in the weights ([`tau_exponent`]). It is -/// deliberately not an [`Algebra`](crate::algebra::Algebra) -/// implementation, because that trait is over $\mathbb{F}_p$; the mod-$\tau$ reduction, which *is* -/// such an algebra, is a follow-up layer on top of this engine, and the honest -/// $\mathbb{F}_2[\tau]$ resolution is built by lifting against this engine (Phase 2). +/// The product is exposed through [`product_indexed`](MotivicMilnorAlgebra::product_indexed), +/// whose $\tau$ coefficients are implicit in the weights ([`tau_exponent`]), and +/// [`bidegree`](MotivicMilnorAlgebra::bidegree) reports them in the $A_C$ convention +/// ([`Dual`]). /// -/// Weight convention: [`bidegree`](MotivicMilnorAlgebra::bidegree) returns `(t, w)` where `t` is -/// the topological degree and `w` is the motivic weight in the presentation where $\tau$ has -/// weight $-1$ and products are weight-homogeneous (i.e. `w = -(dual monomial weight)`; see the -/// module-level note). +/// Deliberately not an [`Algebra`](crate::algebra::Algebra): that trait is over +/// $\mathbb{F}_p$. The mod-$\tau$ reduction, which *is* such an algebra, is a follow-up layer +/// on top of this engine. #[derive(Default)] pub struct MotivicMilnorAlgebra { /// `basis[t]` is the $\mathbb{F}_2[\tau]$-basis in topological degree `t`, sorted for stable @@ -1131,8 +1093,7 @@ pub struct MotivicMilnorAlgebra { /// A dense block of basis-element products for one pair of topological degrees /// `(t1, t2)`. Entry `(idx1, idx2)`, at flat position `idx1 * dim2 + idx2`, is the /// product of the `idx1`-th basis element in degree `t1` with the `idx2`-th in -/// degree `t2` — the `(Tau, index-in-degree-(t1+t2))` list [`multiply_closed`] -/// returns. +/// degree `t2`, as indices into the degree-`(t1 + t2)` basis. /// /// This is the natural **batch unit** for the product: every structure constant /// the resolution needs for a given degree pair lives in one block. Entries fill @@ -1147,11 +1108,9 @@ pub struct ProductBlock { } impl MotivicMilnorAlgebra { + /// An engine with no basis computed yet. pub fn new() -> Self { - Self { - basis: OnceVec::new(), - blocks: RwLock::new(FxHashMap::default()), - } + Self::default() } /// Compute and cache the basis in every topological degree up to and including `degree`. @@ -1297,12 +1256,25 @@ impl MotivicMilnorAlgebra { /// The product of two basis elements as a list of basis indices in degree `t1 + t2`. /// - /// The $\mathbb{F}_2[\tau]$ coefficients are not returned because they are not free: the - /// product of two homogeneous basis elements is weight-homogeneous, so index `j` carries - /// exactly $\tau^{k}$ for the `k` that [`Self::tau_exponent`] computes from the weights. + /// No coefficients: [`Self::tau_exponent`] recovers the one on any index. pub fn product_indexed(&self, t1: i32, idx1: usize, t2: i32, idx2: usize) -> Vec { - self.cached_product(&self.block(t1, t2), t1, idx1, t2, idx2) - .to_vec() + self.product_indexed_with(t1, idx1, t2, idx2, <[_]>::to_vec) + } + + /// Run `f` on the cached product **by reference**, without copying the index list. + /// + /// A cache hit is otherwise free, so a consumer that walks many structure constants — the + /// mod-$\tau$ layer's `multiply_basis_elements`, say — should prefer this to + /// [`Self::product_indexed`], which allocates a `Vec` per hit. + pub fn product_indexed_with( + &self, + t1: i32, + idx1: usize, + t2: i32, + idx2: usize, + f: impl FnOnce(&[usize]) -> R, + ) -> R { + f(self.cached_product(&self.block(t1, t2), t1, idx1, t2, idx2)) } /// One entry of `block`, computed on first request. `block` must be the block for @@ -1488,6 +1460,40 @@ mod tests { } } + #[test] + fn basis_element_from_string_round_trips() { + // The parser is the inverse of the printer, and is what lets `.json` module descriptors + // be written over this algebra. Nothing in this crate calls it yet, so without this it + // would ship untested. + let alg = MotivicMilnorAlgebra::new(); + alg.compute_basis(12); + for degree in 0..=12 { + for idx in 0..alg.dimension(degree) { + let s = alg.basis_element_to_string(degree, idx); + assert_eq!( + alg.basis_element_from_string(&s), + Some((degree, idx)), + "round trip failed on {s:?}" + ); + } + } + assert_eq!(alg.basis_element_from_string("1"), Some((0, 0))); + assert_eq!(alg.basis_element_from_string("Q_0"), Some((1, 0))); + assert_eq!(alg.basis_element_from_string("not an element"), None); + } + + #[test] + fn nb_covers_antidiagonals() { + // `NB` sizes the row / anti-diagonal arrays. An anti-diagonal index is `i + j`, so the + // bound is the widest column plus the number of columns, and nothing may exceed it. + let max_row = PPart::MAX_DEGREE.ilog2() as usize + 1; + let max_col = PPart::MAX_LEN + 1; + assert!( + max_row + max_col <= NB, + "NB = {NB} too small for {max_row} + {max_col}" + ); + } + #[test] fn test_dual_mul_relations() { // τ_0^2 = τ ξ_1 (the defining relation at ρ = 0). The τ^1 is read back off the weights. diff --git a/ext/crates/algebra/src/algebra/motivic/mod.rs b/ext/crates/algebra/src/algebra/motivic/mod.rs index f5891b46c0..7f076982f9 100644 --- a/ext/crates/algebra/src/algebra/motivic/mod.rs +++ b/ext/crates/algebra/src/algebra/motivic/mod.rs @@ -1,25 +1,4 @@ //! The C-motivic prime 2 Steenrod algebra and its mod-$\tau$ reduction. -//! -//! This layer implements the *deformation* view of the C-motivic Adams $E_2$: the -//! C-motivic dual Steenrod algebra $A_C$ over $\mathbb{F}_2[\tau]$, and its mod-$\tau$ -//! reduction $A_C/\tau$ (a connected finite-type $\mathbb{F}_2$-algebra). -//! -//! The foundation layer is [`MotivicMilnorAlgebra`] — $A_C$, a free -//! $\mathbb{F}_2[\tau]$-module on the Milnor basis, with the Kong–Lin product -//! ([`milnor`]). This is the product **engine**; it is deliberately not an -//! [`Algebra`](crate::algebra::Algebra) (that trait is over $\mathbb{F}_p$). -//! -//! There is no representation of the coefficient ring $\mathbb{F}_2[\tau]$ here, because -//! nothing needs one: every element in sight is bidegree-homogeneous and $\tau$ has weight -//! $-1$, so a coefficient is pinned by the weights of the terms it sits between. See -//! [`milnor::tau_exponent`], which is how one is recovered when a caller wants it. Which of the -//! two dual weight conventions applies is carried by the type — [`milnor::Dual`] marks the $A_C$ -//! reading of a monomial — so the two cannot be mixed up. -//! -//! The mod-$\tau$ reduction $A_C/\tau$ — the connected finite-type -//! $\mathbb{F}_2$-algebra that the existing resolution engine resolves to yield the -//! algebraic Novikov $E_2$ — is presented as an [`Algebra`](crate::algebra::Algebra) -//! in a follow-up on top of this engine. pub mod milnor; pub use milnor::MotivicMilnorAlgebra; From b1b0a2da9d106a67904a1d0c57f82e728d7c3f5d Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sat, 29 Aug 2026 19:19:01 -0400 Subject: [PATCH 09/15] motivic: speed up the closed-form product MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two changes to the same hot path, each measured on its own. `enum_y` is the innermost loop — a flamegraph puts essentially all leaf time under it — and its three loops indexed `acc.or[i + j]` and `acc.sum[i + j]` by computed index. The anti-diagonal offset defeats bounds-check elision, so every iteration paid for checks and would not vectorise. Iterating the slices instead, through a paired `Acc::place`/`unplace` that also absorbs the apply/undo bodies written out four times. The block registry then took a write lock on `RwLock>` for every new degree pair, which is the one thing its own doc claimed the design avoids. `once` has this container already: `MultiIndexed<2, V>` is a wait-free sparse map from integer coordinates, and `try_insert` gives the racing-insert loser its value back to drop. Since `get` borrows from `&self`, the `Arc` around each block goes too — `block` returns `&ProductBlock`. motivic_block/12 285 µs -> 247 -> 230 µs -18.7% motivic_block/16 1.79 ms -> 1.57 -> 1.52 ms -15.2% motivic_block/20 8.98 ms -> 7.64 -> 7.16 ms -20.3% motivic_product/xi/large 399 µs -> 308 µs -22.6% motivic_product/q/large 14.2 µs -> 11.4 µs -20.4% `motivic_basis`, which neither change can touch, moved -1% to +0.5% and is the control. The same zip rewrite in `enum_x` measures neutral — X enumeration is not the bottleneck — so it is left indexed rather than changed for symmetry. Co-Authored-By: Claude Opus 5 --- .../algebra/src/algebra/motivic/milnor.rs | 90 +++++++++++-------- 1 file changed, 54 insertions(+), 36 deletions(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index ddd7e287f4..ce862eb0d3 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -5,16 +5,12 @@ //! $\mathbb{F}_2[\tau]$. It uses that paper's *conjugate* generators, so the coproduct and //! Milnor-matrix product match Milnor's classical formulas on the $\xi$ part. -use std::{ - cell::RefCell, - rc::Rc, - sync::{Arc, OnceLock, RwLock}, -}; +use std::{cell::RefCell, rc::Rc, sync::OnceLock}; use fp::prime::{Binomial, iter::BitflagIterator}; use itertools::Itertools; use maybe_rayon::prelude::*; -use once::OnceVec; +use once::{MultiIndexed, OnceVec}; use rustc_hash::FxHashMap; use crate::algebra::milnor_algebra::PPart; @@ -839,6 +835,32 @@ struct Acc { } impl Acc { + /// Add column `cand` at column index `j`, updating rows, anti-diagonal ORs and sums. + fn place(&mut self, cand: &[u32], j: usize) { + for ((&v, r), (o, sm)) in cand + .iter() + .zip(&mut self.rows) + .zip(self.or[j..].iter_mut().zip(&mut self.sum[j..])) + { + *r += v; + *o ^= v; + *sm += v; + } + } + + /// Undo a [`Self::place`] of the same column at the same index. + fn unplace(&mut self, cand: &[u32], j: usize) { + for ((&v, r), (o, sm)) in cand + .iter() + .zip(&mut self.rows) + .zip(self.or[j..].iter_mut().zip(&mut self.sum[j..])) + { + *r -= v; + *o ^= v; + *sm -= v; + } + } + /// Empty accumulator, before any column has been placed. fn zero() -> Self { Self { @@ -962,25 +984,15 @@ impl ClosedY<'_> { self.on_y(acc, out); return; } + // Zipped rather than indexed: the anti-diagonal offset makes `acc.or[i + j]` opaque to + // the bounds-check elision, and this is the innermost loop of the whole product. for cand in self.y_cands[j].iter() { - if !cand - .iter() - .enumerate() - .all(|(i, &v)| acc.or[i + j] & v == 0) - { + if cand.iter().zip(&acc.or[j..]).any(|(&v, &o)| o & v != 0) { continue; } - for (i, &v) in cand.iter().enumerate() { - acc.rows[i] += v; - acc.or[i + j] ^= v; - acc.sum[i + j] += v; - } + acc.place(cand, j); self.enum_y(j + 1, acc, out); - for (i, &v) in cand.iter().enumerate() { - acc.rows[i] -= v; - acc.or[i + j] ^= v; - acc.sum[i + j] -= v; - } + acc.unplace(cand, j); } } @@ -1083,11 +1095,12 @@ pub struct MotivicMilnorAlgebra { /// resolution asks for the same structure constants repeatedly, so we cache /// them (the role the classical Milnor algebra's `cache-multiplication` table /// plays). Blocking by degree pair — rather than a flat `(t1, idx1, t2, idx2)` - /// map — makes the cache the natural *batch unit* (see [`ProductBlock`]): the - /// per-degree-pair registry is small and read-mostly, and each block's entries - /// fill through independent [`OnceLock`]s, so a concurrent resolution reads - /// hits lock-free instead of serializing on one global product lock. - blocks: RwLock>>, + /// map — makes the cache the natural *batch unit* (see [`ProductBlock`]). + /// + /// Nothing here takes a lock: [`MultiIndexed`] is wait-free, and each block's entries fill + /// through independent [`OnceLock`]s, so a concurrent resolution neither serializes on a + /// global product lock nor on registering a new degree pair. + blocks: MultiIndexed<2, ProductBlock>, } /// A dense block of basis-element products for one pair of topological degrees @@ -1218,9 +1231,9 @@ impl MotivicMilnorAlgebra { /// entries empty) on first request and shared thereafter. Creation is the only /// step that touches the block registry's write lock; entry computation happens /// lock-free through the block's [`OnceLock`]s. - fn block(&self, t1: i32, t2: i32) -> Arc { - if let Some(block) = self.blocks.read().unwrap().get(&(t1, t2)) { - return Arc::clone(block); + fn block(&self, t1: i32, t2: i32) -> &ProductBlock { + if let Some(block) = self.blocks.get([t1, t2]) { + return block; } // Ensure the operand and output bases exist before sizing/indexing. self.compute_basis(t1); @@ -1228,13 +1241,18 @@ impl MotivicMilnorAlgebra { self.compute_basis(t1 + t2); let dim1 = self.dimension(t1); let dim2 = self.dimension(t2); - let mut w = self.blocks.write().unwrap(); - Arc::clone(w.entry((t1, t2)).or_insert_with(|| { - Arc::new(ProductBlock { + // A racing thread may have inserted since the `get`; either block is equally good, so + // the loser drops its own and reads the winner's. + let _ = self.blocks.try_insert( + [t1, t2], + ProductBlock { dim2, entries: (0..dim1 * dim2).map(|_| OnceLock::new()).collect(), - }) - })) + }, + ); + self.blocks + .get([t1, t2]) + .expect("just inserted, and entries are never removed") } /// The closed-form product (Kong–Lin Theorem 5.1) of the two basis elements at @@ -1274,7 +1292,7 @@ impl MotivicMilnorAlgebra { idx2: usize, f: impl FnOnce(&[usize]) -> R, ) -> R { - f(self.cached_product(&self.block(t1, t2), t1, idx1, t2, idx2)) + f(self.cached_product(self.block(t1, t2), t1, idx1, t2, idx2)) } /// One entry of `block`, computed on first request. `block` must be the block for @@ -1307,7 +1325,7 @@ impl MotivicMilnorAlgebra { let dim1 = block.entries.len() / dim2; (0..dim1).into_maybe_par_iter().for_each(|idx1| { for idx2 in 0..dim2 { - self.cached_product(&block, t1, idx1, t2, idx2); + self.cached_product(block, t1, idx1, t2, idx2); } }); } From c9e1982acd5611ec270867bfe7faccca476acfca Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sun, 30 Aug 2026 01:16:58 -0400 Subject: [PATCH 10/15] motivic: bound the Y walk by the degree equation instead of testing it at the leaf MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `on_y` rejected a complete Y matrix when Σ(E₁) + 2Σ(S′) ≠ Σ(S(Y)). That target depends only on X, so it is known before the Y walk starts — and measuring the rejection showed the walk was building 1,361,550 complete Y matrices to accept 6, with 1,304,502 of them (96%) dying on that one scalar equation. Each candidate column now carries its unweighted sum, so `ClosedY` can hold the suffix bounds on what the remaining columns can still contribute and stop a branch as soon as the target is out of reach. On the a=[8,4,2,1] b=[4,2,1,1] product: Y matrices reaching on_y 1,361,550 -> 57,048 (24x fewer) rejected by the equation 1,304,502 -> 0 (subsumed by the bound) which is a 3.2x wall-clock win on the degree-138 product (10.96s -> 3.46s), and on the bench: motivic_product/xi/large -40% motivic_product/q/large -28% motivic_block/20 7.38 ms -> 5.60 ms The small cases regress 12-19%: the bounds cost a pass over the candidate lists per X matrix, which a product with a handful of columns cannot amortise. Left alone, since the shapes that regress are microseconds and the ones that gain are the ones that make large computations infeasible. Note the `motivic_basis` control also drifted +12% across this measurement, so treat anything under that as noise; the large-case wins are well clear of it. The leaf test stays as the statement of the condition, now unreachable. Co-Authored-By: Claude Opus 5 --- .../algebra/src/algebra/motivic/milnor.rs | 48 +++++++++++++++++-- 1 file changed, 44 insertions(+), 4 deletions(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index ce862eb0d3..bd3e0585cc 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -708,6 +708,11 @@ struct Columns { entries: Vec, /// `offsets[k]..offsets[k + 1]` is column `k`. Always starts with a 0. offsets: Vec, + /// Unweighted sum Σᵢ vᵢ of each column, parallel to `offsets`. + sums: Vec, + /// The smallest and largest entry of `sums`, for bounding a partial `Y`. + min_sum: u32, + max_sum: u32, } impl Columns { @@ -715,9 +720,19 @@ impl Columns { fn push(&mut self, col: &[u32]) { if self.offsets.is_empty() { self.offsets.push(0); + self.min_sum = u32::MAX; } self.entries.extend_from_slice(col); self.offsets.push(self.entries.len() as u32); + let sum: u32 = col.iter().sum(); + self.sums.push(sum); + self.min_sum = self.min_sum.min(sum); + self.max_sum = self.max_sum.max(sum); + } + + /// The columns paired with their unweighted sums. + fn iter_with_sums(&self) -> impl Iterator { + self.iter().zip(self.sums.iter().copied()) } /// The columns, in insertion order. @@ -952,15 +967,24 @@ impl<'a> Closed<'a> { y_cands.push(columns_eq_cached(t as u32)); } + // The degree equation's target and the suffix bounds, both fixed once `X` is chosen. + let (mut suffix_min, mut suffix_max) = ([0u32; NB], [0u32; NB]); + for j in (0..y_cands.len()).rev() { + suffix_min[j] = suffix_min[j + 1] + y_cands[j].min_sum; + suffix_max[j] = suffix_max[j + 1] + y_cands[j].max_sum; + } ClosedY { y_cands: &y_cands, + target_sigma_sy: self.e1_mask.count_ones() + 2 * sigma_sprime, + suffix_min, + suffix_max, e1_mask: self.e1_mask, e2_mask: self.e2_mask, sprime, sigma_sprime, out_r, } - .enum_y(0, &mut Acc::zero(), out); + .enum_y(0, 0, &mut Acc::zero(), out); } } @@ -968,6 +992,13 @@ impl<'a> Closed<'a> { /// already fixed by the time this is built, so `sprime` and `out_r` are constants here. struct ClosedY<'a> { y_cands: &'a [Rc], + /// `Σ(S(Y))` that [`Self::on_y`]'s degree equation demands. It depends only on `X`, so the + /// `Y` walk can be bounded against it instead of testing it at the leaf. + target_sigma_sy: u32, + /// `suffix_min[j]` / `suffix_max[j]`: the least and greatest `Σ(S(Y))` the columns from `j` + /// on can still contribute. + suffix_min: [u32; NB], + suffix_max: [u32; NB], e1_mask: u32, e2_mask: u32, sprime: &'a [u32], @@ -979,19 +1010,28 @@ impl ClosedY<'_> { /// Enumerate the `Y` matrix column by column (each column an exact weighted sum), /// pruning on anti-diagonal collisions (`b(Y) = 0`) and accumulating each matching /// contribution. - fn enum_y(&self, j: usize, acc: &mut Acc, out: &mut DualElement) { + fn enum_y(&self, j: usize, running: u32, acc: &mut Acc, out: &mut DualElement) { if j == self.y_cands.len() { self.on_y(acc, out); return; } + // `Σ(S(Y))` must land exactly on `target_sigma_sy`, and the columns from `j` on can only + // move it within `[suffix_min, suffix_max]`. Checking that here rather than at the leaf + // is what keeps the walk from building matrices that cannot possibly qualify. + let Some(need) = self.target_sigma_sy.checked_sub(running) else { + return; + }; + if need < self.suffix_min[j] || need > self.suffix_max[j] { + return; + } // Zipped rather than indexed: the anti-diagonal offset makes `acc.or[i + j]` opaque to // the bounds-check elision, and this is the innermost loop of the whole product. - for cand in self.y_cands[j].iter() { + for (cand, sum) in self.y_cands[j].iter_with_sums() { if cand.iter().zip(&acc.or[j..]).any(|(&v, &o)| o & v != 0) { continue; } acc.place(cand, j); - self.enum_y(j + 1, acc, out); + self.enum_y(j + 1, running + sum, acc, out); acc.unplace(cand, j); } } From 0c2e11bf9cca61c16bdb151ae4ea582b15c9d0ee Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sun, 30 Aug 2026 01:41:20 -0400 Subject: [PATCH 11/15] motivic: cut the Y walk on the square-free condition, and skip infeasible X earlier MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two more prunes of the same kind as the degree-equation bound, found by counting what the walk rejects rather than by guessing. `on_y` rejected a complete `Y` when `E₂ + T(Y)` was not square-free. Anti-diagonal sums only grow, so a diagonal that has already overflowed its cap never recovers — testing it as each column is placed cuts the subtree instead of rediscovering the failure at every leaf below it. On a=[8,4,2,1] b=[4,2,1,1] the leaves reaching `on_y` fall from 57,048 to 12, of which 6 are kept; `c_coeff`, which had been rejecting 88% of survivors, is left rejecting 6. `on_x` then built the whole `Y` candidate list before discovering the degree equation was out of reach. The window is decidable from the column targets alone — a column of weighted sum `w` has plain sum between `w.count_ones()` and `w` — so that test moves ahead of the candidate lists. degree 98 product 41.3 ms -> 13.0 ms degree 138 product 3.46 s -> 646 ms motivic_block/20 5.60 ms -> 3.14 ms Two things I tried that measured worse and are not here: carrying the same feasibility bound incrementally through the `X` walk (the popcount floor is too loose to fire — it cut 1.4% of X for a 11% slowdown), and reusing the `Y` candidate vector and accumulator across `X` matrices (threading the scratch costs more than the small `Vec` it saves). Co-Authored-By: Claude Opus 5 --- .../algebra/src/algebra/motivic/milnor.rs | 33 ++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index bd3e0585cc..c533f738df 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -954,6 +954,20 @@ impl<'a> Closed<'a> { return; } + // Whether any `Y` can hit the degree equation is decidable from the column targets alone: + // a column of weighted sum `w` has plain sum between `w.count_ones()` (all weight as high + // as it goes) and `w` (all of it in row 0). Testing that here costs a popcount per column + // and skips building the candidate lists for an `X` that cannot contribute. + let target_sigma_sy = self.e1_mask.count_ones() as i64 + 2 * sigma_sprime as i64; + let (mut lo, mut hi) = (ry0.count_ones() as i64, ry0); + for &t in &ry[1..self.l] { + lo += t.count_ones() as i64; + hi += t; + } + if target_sigma_sy < lo || target_sigma_sy > hi { + return; + } + // The output P-part T(X). The anti-diagonal sums are paper-indexed, so entry `i` of the // packed exponent sequence is anti-diagonal `i + 1` (index 0 is the xi_0 slot). let mut out_r = PPart::zero(); @@ -973,9 +987,14 @@ impl<'a> Closed<'a> { suffix_min[j] = suffix_min[j + 1] + y_cands[j].min_sum; suffix_max[j] = suffix_max[j + 1] + y_cands[j].max_sum; } + let mut e_cap = [1u32; NB]; + for (k, cap) in e_cap.iter_mut().enumerate() { + *cap = 1 - ((self.e2_mask >> k) & 1); + } ClosedY { y_cands: &y_cands, - target_sigma_sy: self.e1_mask.count_ones() + 2 * sigma_sprime, + target_sigma_sy: target_sigma_sy as u32, + e_cap, suffix_min, suffix_max, e1_mask: self.e1_mask, @@ -999,6 +1018,9 @@ struct ClosedY<'a> { /// on can still contribute. suffix_min: [u32; NB], suffix_max: [u32; NB], + /// How much anti-diagonal `k` may still accumulate: `E₂ + T(Y)` has to stay square-free, so + /// a diagonal `E₂` already occupies is capped at 0 and every other at 1. + e_cap: [u32; NB], e1_mask: u32, e2_mask: u32, sprime: &'a [u32], @@ -1030,6 +1052,15 @@ impl ClosedY<'_> { if cand.iter().zip(&acc.or[j..]).any(|(&v, &o)| o & v != 0) { continue; } + // Anti-diagonal sums only grow, so a diagonal that already overflows its cap can + // never come back: cut here rather than rediscovering it at every leaf below. + if cand + .iter() + .zip(acc.sum[j..].iter().zip(&self.e_cap[j..])) + .any(|(&v, (&t, &cap))| t + v > cap) + { + continue; + } acc.place(cand, j); self.enum_y(j + 1, running + sum, acc, out); acc.unplace(cand, j); From 0b7bb1aa00599215c6a561f6aa4ef95e3e0014f5 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sun, 30 Aug 2026 01:41:42 -0400 Subject: [PATCH 12/15] motivic: compute the mod-tau product directly instead of filtering the full one MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The deformation pipeline resolves over `A_C/τ`, not `A_C` — the follow-up's `CTauAlgebra` asks for the full `F_2[τ]` product and keeps the terms of τ-valuation 0. That pays for the whole enumeration to discard nearly all of it. Mod τ the constraint is much tighter, and it is a constraint on the *walk*, not the output: a term carries `τ^{Σ(S′)}`, so keeping only `τ⁰` forces `S′ = 0`, i.e. `S(X) = R₁` exactly rather than `≤`. That is the classical admissible-matrix condition. `multiply_closed_mod_tau` enforces it during the `X` walk, cutting a branch as soon as the remaining columns cannot fill a row. a=[8,4,2,1]·[4,2,1,1] 12.6 ms -> 56.7 µs 222x a=[16,8,2,1]·[8,4,2,1] 676 ms -> 172 µs 3922x which puts the mod-τ product about 115x the classical Milnor product on the same inputs, rather than the ~500,000x the full `A_C` product costs. `test_mod_tau_matches_the_filtered_product` pins the equivalence to filtering `multiply_closed` by `tau_exponent == 0`, which is the whole licence for constraining the walk up front. Co-Authored-By: Claude Opus 5 --- .../algebra/src/algebra/motivic/milnor.rs | 68 +++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index c533f738df..4a7f80b391 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -836,6 +836,16 @@ struct Closed<'a> { e1_mask: u32, e2_mask: u32, sigma2_e1: i64, + /// Keep only the `τ⁰` part of the product. + /// + /// A term carries `τ^{Σ(S′)}`, so dropping everything divisible by `τ` means `S′ = 0`, i.e. + /// `S(X) = R₁` exactly rather than `≤`. That is the classical admissible-matrix condition, + /// and enforcing it during the walk is far cheaper than computing the whole `A_C` product + /// and discarding the τ-divisible terms afterwards. + mod_tau: bool, + /// `row_slack[j][i]`: the most columns `j..l` can still add to row `i`, for the `mod_tau` + /// prune. + row_slack: Vec<[u32; NB]>, /// Per-column candidate columns for the `X` matrix. x_cands: &'a [Rc], } @@ -906,6 +916,17 @@ impl<'a> Closed<'a> { if !fits { continue; } + // With `S(X) = R₁` forced, a row that the remaining columns can no longer fill is + // dead: cut before recursing. + if self.mod_tau { + let slack = &self.row_slack[j + 1]; + if (1..self.l).any(|i| { + let placed = acc.rows[i] + cand.get(i).copied().unwrap_or(0); + placed + slack[i] < self.r1v[i] || placed > self.r1v[i] + }) { + continue; + } + } debug_assert!(cand.len() + j <= NB, "anti-diagonal index out of range"); for (i, &v) in cand.iter().enumerate() { acc.rows[i] += v; @@ -1102,6 +1123,19 @@ impl ClosedY<'_> { /// The product `a · b` in `A_C` via Kong–Lin Theorem 5.1 (ρ = 0). Same contract /// as [`multiply`] (the duality oracle it is validated against). pub fn multiply_closed(a: Dual, b: Dual) -> SteenrodElement { + multiply_closed_inner(a, b, false) +} + +/// The `τ⁰` part of `a · b`, i.e. their product in the mod-`τ` reduction `A_C/τ`. +/// +/// Equal to filtering [`multiply_closed`] down to the terms [`tau_exponent`] puts at `0`, but it +/// constrains the enumeration instead of the output: mod `τ` the `X` matrix must satisfy +/// `S(X) = R₁` exactly, which is the classical admissible-matrix condition. +pub fn multiply_closed_mod_tau(a: Dual, b: Dual) -> SteenrodElement { + multiply_closed_inner(a, b, true) +} + +fn multiply_closed_inner(a: Dual, b: Dual, mod_tau: bool) -> SteenrodElement { // Theorem 5.1 is combinatorics on the (E, R) exponents; the A_C reading is re-applied to the // output monomials at the end. let (a, b) = (a.0, b.0); @@ -1124,6 +1158,21 @@ pub fn multiply_closed(a: Dual, b: Dual) -> SteenrodElement x_cands.push(columns_le_cached(bound)); } + // How much each suffix of columns can still contribute to each row. + let mut row_slack = vec![[0u32; NB]; l + 1]; + for j in (0..l).rev() { + let bound = if j == 0 { None } else { Some(r2v[j]) }; + // `here[j]` is column `j`'s row of the table, `later[0]` is column `j + 1`'s. + let (here, later) = row_slack.split_at_mut(j + 1); + for (i, (slot, &prev)) in here[j].iter_mut().zip(later[0].iter()).enumerate() { + let add = match bound { + Some(b) => b >> i, + None => r1v.get(i).copied().unwrap_or(0), + }; + *slot = prev + add; + } + } + let ctx = Closed { r1v: &r1v, r2v: &r2v, @@ -1131,6 +1180,8 @@ pub fn multiply_closed(a: Dual, b: Dual) -> SteenrodElement e1_mask, e2_mask, sigma2_e1: e1_mask as i64, // Σ (bit i)·2ⁱ = e1_mask + mod_tau, + row_slack, x_cands: &x_cands, }; let mut out = DualElement::new(); @@ -1992,6 +2043,23 @@ mod tests { } } + #[test] + fn test_mod_tau_matches_the_filtered_product() { + // The specialised walk must agree with filtering the general product to tau^0 — that + // equivalence is the whole licence for constraining `S(X) = R_1` up front. + let basis: Vec> = (0..=9).flat_map(enum_basis).collect(); + for &a in &basis { + for &b in &basis { + let filtered: Vec<_> = multiply_closed(a, b) + .into_iter() + .filter(|&z| tau_exponent(a, b, z) == 0) + .collect(); + let direct: Vec<_> = multiply_closed_mod_tau(a, b).into_iter().collect(); + assert_eq!(filtered, direct, "mod-tau mismatch at {a:?} * {b:?}"); + } + } + } + #[test] fn test_closed_form_matches_duality_oracle() { // Exhaustive fuzz: the closed form (Theorem 5.1) agrees with the duality From 967b3570f23b4f12b3168cc7b235309140b01640 Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sun, 30 Aug 2026 16:39:15 -0400 Subject: [PATCH 13/15] motivic: trim the comments before merge MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The module doc was six lines restating what `multiply_closed` and `Monomial` already say; per CLAUDE.md it is one line and the facts move to the items. The conjugate generators are a fact about `Monomial`, and the Kong–Lin citation belongs with the theorem it implements. Also drop a test comment that restated its own degree bounds, so the two cannot drift. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01F25ZVbsP7ULg41iY3MP6FX --- .../algebra/src/algebra/motivic/milnor.rs | 21 ++++++++----------- ext/crates/algebra/src/algebra/motivic/mod.rs | 2 +- 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index 4a7f80b391..e32f6a6cb6 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -1,9 +1,4 @@ //! The C-motivic dual Steenrod algebra $A_{**}$ and Steenrod algebra $A_C$ in the Milnor basis. -//! -//! Follows Kong–Lin, *Product formulas for motivic Milnor basis* (arXiv:2411.12890), at the -//! C-motivic point where $\rho = 0$ and the motivic cohomology of a point is -//! $\mathbb{F}_2[\tau]$. It uses that paper's *conjugate* generators, so the coproduct and -//! Milnor-matrix product match Milnor's classical formulas on the $\xi$ part. use std::{cell::RefCell, rc::Rc, sync::OnceLock}; @@ -50,6 +45,9 @@ fn p_part_from_paper(r: &[u32]) -> Option { /// The same data indexes the $A_C$ basis element $Q(E)P(R)$ that this monomial is dual to; see /// [`Dual`], which is that reading of it. /// +/// These are Kong–Lin's *conjugate* generators, so the coproduct and the Milnor-matrix product +/// agree with Milnor's classical formulas on the $\xi$ part. +/// /// The derived ordering is lexicographic in `(q_part, p_part)`, which is what /// [`MotivicMilnorAlgebra`] sorts and binary-searches its per-degree bases by. #[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash)] @@ -1033,7 +1031,7 @@ impl<'a> Closed<'a> { struct ClosedY<'a> { y_cands: &'a [Rc], /// `Σ(S(Y))` that [`Self::on_y`]'s degree equation demands. It depends only on `X`, so the - /// `Y` walk can be bounded against it instead of testing it at the leaf. + /// `Y` walk can be bounded against it rather than reaching a leaf to test it. target_sigma_sy: u32, /// `suffix_min[j]` / `suffix_max[j]`: the least and greatest `Σ(S(Y))` the columns from `j` /// on can still contribute. @@ -1120,8 +1118,9 @@ impl ClosedY<'_> { } } -/// The product `a · b` in `A_C` via Kong–Lin Theorem 5.1 (ρ = 0). Same contract -/// as [`multiply`] (the duality oracle it is validated against). +/// The product `a · b` in `A_C` via Theorem 5.1 of Kong–Lin, *Product formulas for motivic +/// Milnor basis* (arXiv:2411.12890), at the C-motivic point $\rho = 0$. Same contract as +/// [`multiply`] (the duality oracle it is validated against). pub fn multiply_closed(a: Dual, b: Dual) -> SteenrodElement { multiply_closed_inner(a, b, false) } @@ -2062,10 +2061,8 @@ mod tests { #[test] fn test_closed_form_matches_duality_oracle() { - // Exhaustive fuzz: the closed form (Theorem 5.1) agrees with the duality - // oracle on every ordered pair of basis elements with deg(a)+deg(b) ≤ 18, - // including the τ-carrying and multi-Q cases. (The bound is kept modest - // because the *oracle* is slow; the closed form is not.) + // Exhaustive over ordered pairs of basis elements, including the τ-carrying and + // multi-Q cases. The degree bound is modest because the *oracle* is slow. let basis: Vec<(i32, Dual)> = (0..=14) .flat_map(|t| enum_basis(t).into_iter().map(move |m| (t, m))) .collect(); diff --git a/ext/crates/algebra/src/algebra/motivic/mod.rs b/ext/crates/algebra/src/algebra/motivic/mod.rs index 7f076982f9..c511ef0aaf 100644 --- a/ext/crates/algebra/src/algebra/motivic/mod.rs +++ b/ext/crates/algebra/src/algebra/motivic/mod.rs @@ -1,4 +1,4 @@ -//! The C-motivic prime 2 Steenrod algebra and its mod-$\tau$ reduction. +//! The C-motivic prime 2 Steenrod algebra. pub mod milnor; pub use milnor::MotivicMilnorAlgebra; From 1a5b4881a59bf06d5e7a898f5c632879a22d118b Mon Sep 17 00:00:00 2001 From: Joey Beauvais-Feisthauer Date: Sun, 30 Aug 2026 17:20:58 -0400 Subject: [PATCH 14/15] motivic: drop the basis-element string API `basis_element_from_string` existed so `.json` module descriptors could be written over the algebra, which is a concern of the layer above; nothing here called it, as its own test admitted. `MilnorAlgebra` already parses the same shape, so the mod-tau layer can take it from there rather than from a second parser kept alive by a round-trip test. `basis_element_to_string` took `(degree, idx)` to match a trait this type deliberately does not implement. It becomes `Display` on `Dual`, which is the type that actually has something to print. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01F25ZVbsP7ULg41iY3MP6FX --- .../algebra/src/algebra/motivic/milnor.rs | 107 +++++------------- 1 file changed, 31 insertions(+), 76 deletions(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index e32f6a6cb6..c0d62e1283 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -366,6 +366,26 @@ pub fn tau_exponent(a: W, b: W, z: W) -> u32 { k as u32 } +impl std::fmt::Display for Dual { + /// Prints the $A_C$ basis element as `1`, `Q_i`, `P(R)`, or a product `Q_i … P(R)`, with `R` + /// in the paper's indexing. + fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { + let Self(m) = self; + let mut parts = Vec::new(); + for i in BitflagIterator::set_bit_iterator(m.q_part as u64) { + parts.push(format!("Q_{i}")); + } + if !m.p_part.is_empty() { + parts.push(format!("P({})", m.paper_p_part().iter().format(", "))); + } + if parts.is_empty() { + write!(f, "1") + } else { + write!(f, "{}", parts.join(" ")) + } + } +} + /// An element of the dual algebra $A_{**}$: a mod-2 set of [`Monomial`]s, whose /// $\mathbb{F}_2[\tau]$ coefficients are implicit in the weight ([`tau_exponent`]). pub type DualElement = SparseSum; @@ -1291,63 +1311,6 @@ impl MotivicMilnorAlgebra { ) } - /// Parse a basis element written as `basis_element_to_string` prints it — - /// `1`, `Q_i`, `P(r_0, r_1, …)`, or a space-separated product `Q_i … P(R)` — - /// into its `(topological degree, index)`. The inverse of - /// [`Self::basis_element_to_string`]; returns `None` on a malformed string or a - /// monomial absent from the basis. This is what lets `.json` module descriptors - /// be written over the motivic Steenrod algebra. - pub fn basis_element_from_string(&self, elt: &str) -> Option<(i32, usize)> { - let mut e_mask: u32 = 0; - let mut r: Vec = Vec::new(); - // The `P(…)` block (if any) is last and may contain spaces after commas, so - // split it off before whitespace-tokenizing the `Q_i` factors. - let (q_part, p_part) = match elt.split_once("P(") { - Some((q, rest)) => (q, Some(rest.strip_suffix(')')?)), - None => (elt, None), - }; - for token in q_part.split_whitespace() { - if token == "1" { - continue; // the unit factor - } - let i: u32 = token.strip_prefix("Q_")?.parse().ok()?; - e_mask |= 1 << i; - } - if let Some(p) = p_part { - r = p - .split(',') - .map(|x| x.trim().parse::()) - .collect::>() - .ok()?; - } - - let m = Monomial::from_paper(e_mask, &r)?; - let degree = m.bidegree().0; - if degree < 0 { - return None; - } - self.compute_basis(degree); - // The basis is stored trimmed, so trimming the parsed vector makes this an exact match. - self.index_of(degree, &Dual(m)).map(|idx| (degree, idx)) - } - - /// A display string for a basis element (`Q_i … P(R)`). - pub fn basis_element_to_string(&self, degree: i32, idx: usize) -> String { - let Dual(m) = self.basis_element(degree, idx); - let mut parts = Vec::new(); - for i in fp::prime::iter::BitflagIterator::set_bit_iterator(m.q_part as u64) { - parts.push(format!("Q_{i}")); - } - if !m.p_part.is_empty() { - parts.push(format!("P({})", m.paper_p_part().iter().format(", "))); - } - if parts.is_empty() { - "1".to_string() - } else { - parts.join(" ") - } - } - /// The dense [`ProductBlock`] for the degree pair `(t1, t2)`, created (with all /// entries empty) on first request and shared thereafter. Creation is the only /// step that touches the block registry's write lock; entry computation happens @@ -1600,25 +1563,17 @@ mod tests { } #[test] - fn basis_element_from_string_round_trips() { - // The parser is the inverse of the printer, and is what lets `.json` module descriptors - // be written over this algebra. Nothing in this crate calls it yet, so without this it - // would ship untested. - let alg = MotivicMilnorAlgebra::new(); - alg.compute_basis(12); - for degree in 0..=12 { - for idx in 0..alg.dimension(degree) { - let s = alg.basis_element_to_string(degree, idx); - assert_eq!( - alg.basis_element_from_string(&s), - Some((degree, idx)), - "round trip failed on {s:?}" - ); - } - } - assert_eq!(alg.basis_element_from_string("1"), Some((0, 0))); - assert_eq!(alg.basis_element_from_string("Q_0"), Some((1, 0))); - assert_eq!(alg.basis_element_from_string("not an element"), None); + fn basis_elements_display() { + let one = Dual(Monomial::one()); + assert_eq!(one.to_string(), "1"); + assert_eq!( + Dual(Monomial::from_paper(1, &[]).unwrap()).to_string(), + "Q_0" + ); + assert_eq!( + Dual(Monomial::from_paper(0b101, &[0, 2, 1]).unwrap()).to_string(), + "Q_0 Q_2 P(0, 2, 1)" + ); } #[test] From 25855a4be58628f127f19535d35cc3963d8321fd Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 30 Aug 2026 23:44:55 +0000 Subject: [PATCH 15/15] motivic: fix a compute_basis race, and build the slack table only for mod-tau compute_basis read the basis length outside the write lock and then pushed, so two threads reaching an uncached degree together (fill_block runs the block walk in parallel) could both fill the same degrees. The vec ends up longer than it should be and basis[t] stops holding degree t, which then feeds wrong structure constants into every cached product above it. OnceVec::extend re-reads the length under the lock, which is what the classical algebra already uses. The added test fails against the previous code ("degree 11 misaligned after concurrent first use") and passes now. The negative-degree guard keeps the old no-op behaviour, since the loop bound tolerated a negative degree but `degree as usize` would not. Separately, row_slack is read only by the mod-tau prune, so the general product was allocating and filling (l + 1) * NB entries on every call and never reading them. It moves into row_slack_table, built only for that walk. Co-Authored-By: Claude Claude-Session: https://claude.ai/code/session_013ePtYD7Bt4iPeCtmqtqvZE --- .../algebra/src/algebra/motivic/milnor.rs | 81 ++++++++++++++----- 1 file changed, 63 insertions(+), 18 deletions(-) diff --git a/ext/crates/algebra/src/algebra/motivic/milnor.rs b/ext/crates/algebra/src/algebra/motivic/milnor.rs index c0d62e1283..b8be3e4102 100644 --- a/ext/crates/algebra/src/algebra/motivic/milnor.rs +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -1154,6 +1154,25 @@ pub fn multiply_closed_mod_tau(a: Dual, b: Dual) -> Steenrod multiply_closed_inner(a, b, true) } +/// `row_slack[j][i]`: the most columns `j..l` can still add to row `i`. Read only by the +/// mod-τ prune in [`Closed::enum_x`], so it is built only for that walk. +fn row_slack_table(l: usize, r1v: &[u32], r2v: &[u32]) -> Vec<[u32; NB]> { + let mut row_slack = vec![[0u32; NB]; l + 1]; + for j in (0..l).rev() { + let bound = if j == 0 { None } else { Some(r2v[j]) }; + // `here[j]` is column `j`'s row of the table, `later[0]` is column `j + 1`'s. + let (here, later) = row_slack.split_at_mut(j + 1); + for (i, (slot, &prev)) in here[j].iter_mut().zip(later[0].iter()).enumerate() { + let add = match bound { + Some(b) => b >> i, + None => r1v.get(i).copied().unwrap_or(0), + }; + *slot = prev + add; + } + } + row_slack +} + fn multiply_closed_inner(a: Dual, b: Dual, mod_tau: bool) -> SteenrodElement { // Theorem 5.1 is combinatorics on the (E, R) exponents; the A_C reading is re-applied to the // output monomials at the end. @@ -1177,20 +1196,12 @@ fn multiply_closed_inner(a: Dual, b: Dual, mod_tau: bool) -> x_cands.push(columns_le_cached(bound)); } - // How much each suffix of columns can still contribute to each row. - let mut row_slack = vec![[0u32; NB]; l + 1]; - for j in (0..l).rev() { - let bound = if j == 0 { None } else { Some(r2v[j]) }; - // `here[j]` is column `j`'s row of the table, `later[0]` is column `j + 1`'s. - let (here, later) = row_slack.split_at_mut(j + 1); - for (i, (slot, &prev)) in here[j].iter_mut().zip(later[0].iter()).enumerate() { - let add = match bound { - Some(b) => b >> i, - None => r1v.get(i).copied().unwrap_or(0), - }; - *slot = prev + add; - } - } + // Only the mod-τ prune reads the slack table, so the general product does not build it. + let row_slack = if mod_tau { + row_slack_table(l, &r1v, &r2v) + } else { + Vec::new() + }; let ctx = Closed { r1v: &r1v, @@ -1270,11 +1281,16 @@ impl MotivicMilnorAlgebra { /// Compute and cache the basis in every topological degree up to and including `degree`. /// Idempotent and cheap to re-call. pub fn compute_basis(&self, degree: i32) { - for t in self.basis.len() as i32..=degree { - let mut b = enum_basis(t); - b.sort(); - self.basis.push(b); + if degree < 0 { + return; } + // `extend` re-reads the length under the write lock, so two threads racing to first-use + // cannot both fill the same degree and leave `basis[t]` holding some other degree. + self.basis.extend(degree as usize, |t| { + let mut b = enum_basis(t as i32); + b.sort(); + b + }); } /// The $\mathbb{F}_2[\tau]$-rank of `A_C` in topological degree `degree`. @@ -1997,6 +2013,35 @@ mod tests { } } + #[test] + fn test_concurrent_compute_basis_keeps_degrees_aligned() { + const MAX: i32 = 12; + let reference: Vec = { + let alg = MotivicMilnorAlgebra::new(); + alg.compute_basis(MAX); + (0..=MAX).map(|t| alg.dimension(t)).collect() + }; + for _ in 0..40 { + let alg = std::sync::Arc::new(MotivicMilnorAlgebra::new()); + let hs: Vec<_> = (0..8) + .map(|_| { + let a = alg.clone(); + std::thread::spawn(move || a.compute_basis(MAX)) + }) + .collect(); + for h in hs { + h.join().unwrap(); + } + for t in 0..=MAX { + assert_eq!( + alg.dimension(t), + reference[t as usize], + "degree {t} misaligned after concurrent first use" + ); + } + } + } + #[test] fn test_mod_tau_matches_the_filtered_product() { // The specialised walk must agree with filtering the general product to tau^0 — that