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..2ca72f6f26 --- /dev/null +++ b/ext/crates/algebra/benches/motivic.rs @@ -0,0 +1,106 @@ +//! 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::{Dual, 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]) -> Dual { + let mut r = vec![0]; + r.extend_from_slice(xi); + Dual(Monomial::from_paper(q.iter().map(|i| 1 << i).sum(), &r).unwrap()) +} + +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))); + }); +} + +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); diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index afbe59bf8c..38d54960df 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 { @@ -1425,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, @@ -1525,7 +1552,7 @@ impl PPartMultiplier { }) .unwrap_or(max + 1) } else { - ((k | sum) + 1) & !sum + next_disjoint(sum, k) } } _ => (k + 1..max + 1) @@ -2200,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. diff --git a/ext/crates/algebra/src/algebra/mod.rs b/ext/crates/algebra/src/algebra/mod.rs index 67baed2b04..0bb7e9fa1f 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; + 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..b8be3e4102 --- /dev/null +++ b/ext/crates/algebra/src/algebra/motivic/milnor.rs @@ -0,0 +1,2120 @@ +//! The C-motivic dual Steenrod algebra $A_{**}$ and Steenrod algebra $A_C$ in the Milnor basis. + +use std::{cell::RefCell, rc::Rc, sync::OnceLock}; + +use fp::prime::{Binomial, iter::BitflagIterator}; +use itertools::Itertools; +use maybe_rayon::prelude::*; +use once::{MultiIndexed, OnceVec}; +use rustc_hash::FxHashMap; + +use crate::algebra::milnor_algebra::PPart; + +/// 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 +} + +/// 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. +/// +/// 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)] +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 { + Some(Self::new(q_part, p_part_from_paper(r)?)) + } + + /// 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() + } +} + +/// $\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 `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}$ +/// 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 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 += (entry(s, i) - entry(r, i)) << i; + } + // 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; + } + } + 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. 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]; + 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). +// --------------------------------------------------------------------------- + +/// A set of `K`s taken mod 2, as a `Vec` sorted by key: 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 +/// searched by [`slice::binary_search`] beats a tree that allocates per node. +#[derive(Debug, Clone, Default, PartialEq, Eq)] +pub struct SparseSum { + terms: Vec, +} + +impl SparseSum { + /// The zero combination. + pub fn new() -> Self { + Self { terms: Vec::new() } + } + + /// Whether this is zero. + pub fn is_empty(&self) -> bool { + self.terms.is_empty() + } + + /// Whether `key` occurs. + pub fn contains(&self, key: &K) -> bool { + self.terms.binary_search(key).is_ok() + } + + /// 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) => { + self.terms.remove(i); + } + Err(i) => self.terms.insert(i, key), + } + } +} + +impl FromIterator for SparseSum { + fn from_iter>(iter: I) -> Self { + let mut out = Self::new(); + for key in iter { + out.add_term(key); + } + out + } +} + +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; + type Item = K; + + fn into_iter(self) -> Self::IntoIter { + self.terms.into_iter() + } +} + +impl<'a, K> IntoIterator for &'a SparseSum { + type IntoIter = std::slice::Iter<'a, K>; + type Item = &'a K; + + fn into_iter(self) -> Self::IntoIter { + self.terms.iter() + } +} + +/// The basis element of $A_C$ dual to a monomial of $A_{**}$. +/// +/// $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); +} + +/// 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 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 +} + +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; + +/// 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(); + 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, 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 = ppart_add(m1.p_part, m2.p_part); + for term in rewrite_tau(&s) { + 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 + // number; this is the invariant that lets the coefficient go unstored. + debug_assert_eq!(term.tau_pow, 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 in a { + for &m2 in b { + mul_monomials(m1, m2, &mut out); + } + } + out +} + +/// The generator $\tau_i \in A_{**}$. +pub fn tau_gen(i: usize) -> DualElement { + 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 { + xi_pow_elt(i, 1) +} + +/// The unit $1 \in A_{**}$. +pub fn dual_one() -> DualElement { + DualElement::from([Monomial::one()]) +} + +// --------------------------------------------------------------------------- +// The coproduct ψ: A_** → A_** ⊗ A_**, an algebra map (Kong–Lin §2.2). +// --------------------------------------------------------------------------- + +/// 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$. +fn tensor_one() -> TensorElement { + 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) in t1 { + for &(bl, br) in t2 { + let mut left = DualElement::new(); + mul_monomials(al, bl, &mut left); + let mut right = DualElement::new(); + mul_monomials(ar, br, &mut right); + for &ml in &left { + for &mr in &right { + out.add_term((ml, mr)); + } + } + } + } + out +} + +/// The monomial $\xi_j^p$ (with $\xi_0 = 1$). +fn xi_pow_mon(j: usize, p: u32) -> Monomial { + 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(); + 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), + )); + } + 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 { + out.add_term((xi_pow_mon(i, 1), xi_pow_mon(k - i, 1 << i))); + } + 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(m: &Monomial) -> TensorElement { + let mut acc = tensor_one(); + for i in BitflagIterator::set_bit_iterator(m.q_part as u64) { + acc = tensor_mul(&acc, &coprod_tau(i)); + } + for (i, r) in m.p_part.iter().enumerate() { + for _ in 0..r { + acc = tensor_mul(&acc, &coprod_xi(i + 1)); + } + } + 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 &m in elt { + for key in coproduct_monomial(&m) { + out.add_term(key); + } + } + 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)]) +} + +/// $\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 { + for mon in dual_mul(&chi_xi(i), &xi_pow_elt(k - i, 1 << i)) { + acc.add_term(mon); + } + } + 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 { + for mon in dual_mul(&chi_tau(i), &xi_pow_elt(k - i, 1 << i)) { + acc.add_term(mon); + } + } + 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 &mon in elt { + let mut acc = dual_one(); + 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 { + acc = dual_mul(&acc, &chi_xi(i + 1)); + } + } + for out_mon in acc { + out.add_term(out_mon); + } + } + out +} + +// --------------------------------------------------------------------------- +// The product in the Steenrod algebra A_C, computed by dualizing ψ. +// --------------------------------------------------------------------------- + +/// 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 { + // `r` is paper-indexed: `r[0]` is the (always zero) xi_0 slot. + out.push(Dual( + Monomial::from_paper(*e, r).expect("basis exponent out of packing range"), + )); + } + 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 $A_C$, as the set of basis elements occurring — their +/// $\mathbb{F}_2[\tau]$ coefficients being implicit ([`tau_exponent`]). +/// +/// 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); + let t = a.bidegree().0 + b.bidegree().0; + let mut out = SteenrodElement::new(); + for z in enum_basis(t) { + if coproduct_monomial(&z.0).contains(&key) { + out.add_term(z); + } + } + 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. + +/// 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, + /// 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 { + /// Append one column. + 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. + 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 { + 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 +} + +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()); +} + +/// All columns with `Σᵢ 2ⁱ vᵢ = target`. +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).into_iter().collect())), + ) + }) +} + +/// 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((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]) -> Columns { + 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. +/// +/// Anti-diagonals are indexed by `i + j`, with the row index `i` bounded by a column's length +/// (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). +struct Closed<'a> { + r1v: &'a [u32], + r2v: &'a [u32], + l: usize, + 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], +} + +/// 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 { + /// 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 { + 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; + } + // 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; + 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′). 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 = [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 { + 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..self.l].iter().sum::(); + if ry0 < 0 { + 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(); + 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); + y_cands.push(columns_eq_cached(ry0 as u32)); + for &t in &ry[1..self.l] { + 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; + } + 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: target_sigma_sy as u32, + e_cap, + suffix_min, + suffix_max, + e1_mask: self.e1_mask, + e2_mask: self.e2_mask, + sprime, + sigma_sprime, + out_r, + } + .enum_y(0, 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], + /// `Σ(S(Y))` that [`Self::on_y`]'s degree equation demands. It depends only on `X`, so the + /// `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. + 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], + sigma_sprime: u32, + out_r: PPart, +} + +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, 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, sum) in self.y_cands[j].iter_with_sums() { + 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); + } + } + + /// 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, &t) in acc.sum.iter().enumerate() { + let val = ((self.e2_mask >> i) & 1) + t; + if val > 1 { + return; + } + if val == 1 { + out_e_mask |= 1 << i; + } + } + out.add_term(Monomial::new(out_e_mask, self.out_r)); + } +} + +/// 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) +} + +/// 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) +} + +/// `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. + 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. + 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); + x_cands.push(Rc::new(col0_x_options(&r1v))); + for &bound in r2v.iter().skip(1) { + x_cands.push(columns_le_cached(bound)); + } + + // 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, + r2v: &r2v, + l, + 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(); + 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.into_iter().map(Dual).collect() +} + +// --------------------------------------------------------------------------- +// 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. +/// +/// 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`]). +/// +/// 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 + /// 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`]). + /// + /// 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 +/// `(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`, 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 +/// 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 { + /// An engine with no basis computed yet. + pub fn new() -> Self { + Self::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) { + 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`. + 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) -> &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: &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) { + 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 [`tau_exponent`]. + pub fn tau_exponent(&self, t1: i32, idx1: usize, t2: i32, idx2: usize, idx: usize) -> u32 { + tau_exponent( + self.basis[t1 as usize][idx1], + self.basis[t2 as usize][idx2], + self.basis[(t1 + t2) as usize][idx], + ) + } + + /// 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) -> &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); + self.compute_basis(t2); + self.compute_basis(t1 + t2); + let dim1 = self.dimension(t1); + let dim2 = self.dimension(t2); + // 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 + /// `(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 { + 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| { + self.index_of(t, &z) + .expect("product landed outside the basis") + }) + .collect() + } + + /// The product of two basis elements as a list of basis indices in degree `t1 + t2`. + /// + /// 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.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 + /// `(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 [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 std::collections::BTreeMap; + + 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) + } + + /// 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)) + } + + /// 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(); + 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), &dmon(0u32, &[])); + assert_eq!(alg.dimension(1), 1); + assert_eq!(alg.basis_element(1, 0), &dmon(0b1u32, &[])); // Q_0 + assert_eq!(alg.dimension(2), 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: SteenrodElement = alg + .product_indexed(1, 0, 2, 0) + .into_iter() + .map(|idx| *alg.basis_element(3, idx)) + .collect(); + assert_eq!( + terms, + SteenrodElement::from([dmon(0b10, &[]), dmon(0b1, &[0, 1])]) + ); + } + + #[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); + let b = *alg.basis_element(t2, idx2); + // Indexed product agrees with the raw monomial product. + let indexed: SteenrodElement = alg + .product_indexed(t1, idx1, t2, idx2) + .into_iter() + .map(|idx| *alg.basis_element(t1 + t2, idx)) + .collect(); + assert_eq!(indexed, multiply(a, b)); + + // 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), + tau_exponent(a, b, z) + ); + } + } + } + } + } + } + + #[test] + 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] + 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. + let (t0, t1) = (mon(0b1, &[]), mon(0b10, &[])); + assert_eq!( + dual_mul(&tau_gen(0), &tau_gen(0)), + DualElement::from([mon(0, &[0, 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!(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!(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!(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)); + } + + #[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(dpaper(p1), dpaper(p2)) + .into_iter() + .filter(|d| d.0.q_part == 0) // keep the pure-ξ outputs + .map(|d| { + assert_eq!( + tau_exponent(dpaper(p1), dpaper(p2), d), + 0, + "pure-ξ output carried a τ power" + ); + d.0.p_part.iter().collect() + }) + .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(pp).bidegree().0; + let xi_pps = |t: i32| -> Vec> { + enum_basis(t) + .into_iter() + .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 { + 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) { + *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(dpaper(&r1), dpaper(&r2)) + .into_iter() + .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)) { + 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` cancels + // mod 2, and the τ coefficients agree termwise because both sides are homogeneous). + 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); + } + } + 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 = SteenrodElement::new(); + for &z in &bc { + for w in multiply(*a, z) { + rhs.add_term(w); + } + } + 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 = dmon(0b1, &[]); + let p_xi1 = dmon(0, &[0, 1]); + assert_eq!( + multiply(q0, p_xi1), + SteenrodElement::from([ + dmon(0b10, &[]), // Q_1 + dmon(0b1, &[0, 1]), // Q_0 P(ξ_1) + ]) + ); + } + + #[test] + fn test_product_unit_and_squares() { + // 1 · x = 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( + Dual(Monomial::new(1 << i, PPart::zero())), + Dual(Monomial::new(1 << i, PPart::zero())) + ) + .is_empty(), + "Q_{i}^2 ≠ 0" + ); + } + } + + #[test] + fn test_product_weight_homogeneous_and_tau_appears() { + // 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 { + 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 |= tau_exponent(*a, *b, z) > 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([mon(0, &[0, 0, 1]), mon(0, &[0, 3])]) + ); + 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) in coproduct(x) { + let prod = dual_mul(&antipode(&DualElement::from([l])), &DualElement::from([r])); + for mon in prod { + out.add_term(mon); + } + } + 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([ + (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])), // 1 ⊗ ξ_2 + (mon(0, &[0, 1]), mon(0, &[0, 2])), // ξ_1 ⊗ ξ_1^2 + (mon(0, &[0, 0, 1]), mon(0, &[])), // ξ_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 = 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. + 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( + Dual(Monomial::new(1 << i, PPart::zero())), + Dual(Monomial::new(1 << i, PPart::zero())) + ), + multiply( + Dual(Monomial::new(1 << i, PPart::zero())), + Dual(Monomial::new(1 << i, PPart::zero())) + ) + ); + } + } + + #[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 + // 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 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(); + 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..c511ef0aaf --- /dev/null +++ b/ext/crates/algebra/src/algebra/motivic/mod.rs @@ -0,0 +1,4 @@ +//! The C-motivic prime 2 Steenrod algebra. + +pub mod milnor; +pub use milnor::MotivicMilnorAlgebra;