From e379b9cbbe27746c1685377e88197daa056abee9 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 31 Jul 2026 05:59:15 +0000 Subject: [PATCH 1/6] Bit-pack Milnor basis elements into a u64 The p-part of a Milnor basis element was a `Vec`, costing a heap allocation and a pointer chase per element. At p = 2 the internal degree of P(R) is sum_i r_i (2^i - 1) with non-negative terms, so r_i <= deg/(2^i - 1); sizing each field by that bound packs the whole exponent sequence into 64 bits for every degree up to 2045. At odd primes the same bound applies divided by q = 2(p-1), so one layout serves every prime. `MilnorBasisElement` is now 16 bytes, `Copy`, and entirely inline. Measured over degrees 0..=300 at p = 2, `basis_table` drops from 252 MiB in 5,036,688 allocations to 77 MiB in none. Three things fall out of the packing: - The packed value is a canonical key, so the hand-rolled `MilnorHashMap` specialization for `not(odd-primes)` is gone; a plain `HashMap` now hashes a single word on every path. That code also assumed a degree bound of 1536 without enforcing it. `compute_basis` now asserts the bound up front, which is what lets everything downstream skip range checks. - Trailing zeros are not represented, so the "pop trailing zeros" loops after building a product disappear. - `PPartMultiplier` no longer borrows its inputs, so its lifetime parameter is gone, and `PPartAllocation` loses the buffer it existed to recycle. In `ext`, `MilnorSubalgebra`'s signature test becomes one masked comparison on the packed word instead of a loop over entries, with the mask hoisted out of `signature_mask`'s inner loop. Two behaviour changes worth noting: - `basis_element_from_string("P0")` and `("Sq0")` now return the identity rather than `None`. P(0) is the identity, and `AdemAlgebra::try_beps_pn` already special-cases `x == 0` this way; the old `None` came from `vec![0]` and `vec![]` hashing differently, an artifact of the representation. - `increment_p_part` now carries before incrementing. The old order transiently stored `max[i] + 1`, which need not fit a field whose width is exactly saturated by `max[i]`. The enumeration is unchanged. The observation that every Milnor exponent sequence up to degree 512 fits in 64 bits is due to Lixiong Wu; this implementation works out the widths, finds that the same layout holds all the way to degree 2045, and carries it through the algebra. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV --- ext/crates/algebra/benches/milnor.rs | 32 +- .../algebra/src/algebra/milnor_algebra.rs | 733 ++++++++++++------ .../algebra/src/algebra/pair_algebra.rs | 38 +- ext/crates/algebra/src/module/rpn.rs | 2 +- ext/crates/algebra/src/steenrod_evaluator.rs | 29 +- ext/crates/algebra/src/steenrod_parser.rs | 4 +- ext/examples/bruner.rs | 7 +- ext/examples/sq0.rs | 9 +- ext/src/nassau.rs | 31 +- ext/src/yoneda.rs | 2 +- 10 files changed, 588 insertions(+), 299 deletions(-) diff --git a/ext/crates/algebra/benches/milnor.rs b/ext/crates/algebra/benches/milnor.rs index 6a38188e59..696bb7325e 100644 --- a/ext/crates/algebra/benches/milnor.rs +++ b/ext/crates/algebra/benches/milnor.rs @@ -1,6 +1,6 @@ //! Benchmarks for the low-level Milnor `PPartMultiplier` kernel. -use algebra::milnor_algebra::{PPartAllocation, PPartEntry, PPartMultiplier}; +use algebra::milnor_algebra::{PPart, PPartAllocation, PPartMultiplier}; use criterion::{ BenchmarkGroup, Criterion, criterion_group, criterion_main, measurement::WallTime, }; @@ -11,8 +11,8 @@ fn bench_ppart( g: &mut BenchmarkGroup, name: &str, p: u32, - r: Vec, - s: Vec, + r: PPart, + s: PPart, ) { let p = ValidPrime::new(p); g.bench_function(name, |bench| { @@ -21,7 +21,7 @@ fn bench_ppart( bench.iter_batched( PPartAllocation::default, |alloc| { - let m = PPartMultiplier::::new_from_allocation(p, &r, &s, alloc, 0, 0); + let m = PPartMultiplier::::new_from_allocation(p, r, s, alloc, 0, 0); for c in m { std::hint::black_box(c); } @@ -38,30 +38,30 @@ fn ppart(c: &mut Criterion) { &mut g, "ppart_2/a", 2, - vec![60, 30, 8, 2, 1], - vec![20, 30, 20, 4, 1, 2], + PPart::from_slice(&[60, 30, 8, 2, 1]), + PPart::from_slice(&[20, 30, 20, 4, 1, 2]), ); bench_ppart::( &mut g, "ppart_2/b", 2, - vec![35, 12, 20, 14, 1, 3], - vec![60, 30, 0, 2, 1], + PPart::from_slice(&[35, 12, 20, 14, 1, 3]), + PPart::from_slice(&[60, 30, 0, 2, 1]), ); bench_ppart::( &mut g, "ppart_4/a", 2, - vec![60, 30, 8, 2, 1], - vec![20, 30, 20, 4, 1, 2], + PPart::from_slice(&[60, 30, 8, 2, 1]), + PPart::from_slice(&[20, 30, 20, 4, 1, 2]), ); bench_ppart::( &mut g, "ppart_4/b", 2, - vec![35, 12, 20, 14, 1, 3], - vec![60, 30, 0, 2, 1], + PPart::from_slice(&[35, 12, 20, 14, 1, 3]), + PPart::from_slice(&[60, 30, 0, 2, 1]), ); #[cfg(feature = "odd-primes")] @@ -70,15 +70,15 @@ fn ppart(c: &mut Criterion) { &mut g, "ppart_3/a", 3, - vec![120, 70, 40, 2], - vec![60, 35, 21, 6], + PPart::from_slice(&[120, 70, 40, 2]), + PPart::from_slice(&[60, 35, 21, 6]), ); bench_ppart::( &mut g, "ppart_3/b", 3, - vec![30, 12, 35, 24], - vec![100, 80, 16, 2, 3], + PPart::from_slice(&[30, 12, 35, 24]), + PPart::from_slice(&[100, 80, 16, 2, 3]), ); } diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index 07b3bd35f8..bbf780041c 100644 --- a/ext/crates/algebra/src/algebra/milnor_algebra.rs +++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs @@ -24,8 +24,11 @@ pub struct MilnorProfile { #[serde(default = "q_part_default")] pub q_part: u32, /// The profile function for the Q part. + /// + /// Unlike the exponent sequence of a basis element (see [`PPart`]), these are *exponents* of + /// the profile function and use [`PPartEntry::MAX`] to mean infinity, so they stay unpacked. #[serde(default)] - pub p_part: PPart, + pub p_part: Vec, } impl MilnorProfile { @@ -99,9 +102,207 @@ impl Default for MilnorProfile { } pub type PPartEntry = u32; -pub type PPart = Vec; -#[derive(Debug, Clone, Default)] +/// The exponent sequence $(r_1, r_2, \ldots)$ of a Milnor basis element $P(r_1, r_2, \ldots)$, +/// bit-packed into a single `u64`. +/// +/// Entry $r_{i+1}$ occupies [`Self::WIDTHS`]`[i]` bits starting at bit [`Self::SHIFTS`]`[i]`. The +/// widths are forced by the degree bound: at $p = 2$ the internal degree of $P(R)$ is +/// $\sum_i r_i (2^i - 1)$ and every term is non-negative, so an element of degree at most +/// [`Self::MAX_DEGREE`] has $r_i \le \mathrm{MAX\\_DEGREE}/(2^i - 1)$. At an odd prime the same +/// argument bounds $r_i$ by that quantity divided by $q = 2(p-1)$, so the $p = 2$ widths are valid +/// for every prime and this type is prime-agnostic. +/// +/// Trailing zeros are not represented: $P(2, 1)$ and $P(2, 1, 0)$ have the same packed value. That +/// is what makes the packed value a canonical key, and it makes [`Self::len`] the position of the +/// highest non-zero entry rather than a stored field. +/// +/// # Invariant +/// +/// Every entry fits in its field. This holds for any element of degree at most +/// [`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)] +pub struct PPart(u64); + +impl PPart { + /// The largest internal degree whose exponent sequences are guaranteed to fit. + /// + /// This is the largest bound for which [`Self::WIDTHS`] sums to at most 64. It is far beyond + /// anything reachable — the Milnor algebra already has over 5 million basis elements below + /// degree 300 — and exceeds the degree 1536 the previous hand-rolled packing assumed. + pub const MAX_DEGREE: i32 = 2045; + + /// The number of entries that can be stored. This equals `fp`'s `MAX_MULTINOMIAL_LEN`, which + /// already bounds the length of the $\xi$-degree table, so it is not a new restriction. + pub const MAX_LEN: usize = 10; + + /// `WIDTHS[i]` is the number of bits holding $r_{i+1}$: the number of bits needed to represent + /// `MAX_DEGREE / (2^(i+1) - 1)`. + const WIDTHS: [u32; Self::MAX_LEN] = [11, 10, 9, 8, 7, 6, 5, 4, 3, 1]; + + /// `SHIFTS[i]` is the bit offset of entry `i`; `SHIFTS[MAX_LEN]` is the total width, 64. + const SHIFTS: [u32; Self::MAX_LEN + 1] = { + let mut shifts = [0; Self::MAX_LEN + 1]; + let mut i = 0; + while i < Self::MAX_LEN { + shifts[i + 1] = shifts[i] + Self::WIDTHS[i]; + i += 1; + } + shifts + }; + + /// `FIELD_OF_BIT[b]` is the index of the entry owning bit `b`, letting [`Self::len`] turn a + /// `leading_zeros` into an entry index without looping. + const FIELD_OF_BIT: [u8; 64] = { + let mut table = [0; 64]; + let mut i = 0; + while i < Self::MAX_LEN { + let mut b = Self::SHIFTS[i]; + while b < Self::SHIFTS[i + 1] { + table[b as usize] = i as u8; + b += 1; + } + i += 1; + } + table + }; + + /// The largest value entry `i` can hold. + pub const fn max_entry(i: usize) -> PPartEntry { + ((1u64 << Self::WIDTHS[i]) - 1) as PPartEntry + } + + /// The number of bits holding entry `i`. Together with [`Self::shift`] this lets callers build + /// a mask over [`Self::bits`] directly, e.g. to test many entries in one comparison. + pub const fn width(i: usize) -> u32 { + Self::WIDTHS[i] + } + + /// The bit offset of entry `i` within [`Self::bits`]. + pub const fn shift(i: usize) -> u32 { + Self::SHIFTS[i] + } + + const fn mask(i: usize) -> u64 { + ((1u64 << Self::WIDTHS[i]) - 1) << Self::SHIFTS[i] + } + + pub const fn zero() -> Self { + Self(0) + } + + /// The raw packed value. Two exponent sequences are equal exactly when their bits are, so this + /// is a complete hash key, and it can be compared against a packed mask in one operation (see + /// `MilnorSubalgebra::has_signature` in `ext`). + pub const fn bits(self) -> u64 { + self.0 + } + + /// Entry `i`, or 0 if `i` is past the end. + #[inline] + pub const fn get(self, i: usize) -> PPartEntry { + if i >= Self::MAX_LEN { + return 0; + } + ((self.0 >> Self::SHIFTS[i]) & ((1 << Self::WIDTHS[i]) - 1)) as PPartEntry + } + + /// Set entry `i` to `v`. + /// + /// # Panics + /// + /// If `i >= MAX_LEN`, or `v` does not fit in entry `i`. Both are unreachable for elements of + /// degree at most [`Self::MAX_DEGREE`]. + #[inline] + pub fn set(&mut self, i: usize, v: PPartEntry) { + assert!(i < Self::MAX_LEN, "p-part index {i} out of range"); + assert!( + v <= Self::max_entry(i), + "p-part entry {v} does not fit in the {} bits at index {i}", + Self::WIDTHS[i], + ); + self.0 = (self.0 & !Self::mask(i)) | ((v as u64) << Self::SHIFTS[i]); + } + + /// The number of entries up to and including the last non-zero one. + #[inline] + pub const fn len(self) -> usize { + if self.0 == 0 { + 0 + } else { + Self::FIELD_OF_BIT[63 - self.0.leading_zeros() as usize] as usize + 1 + } + } + + #[inline] + pub const fn is_empty(self) -> bool { + self.0 == 0 + } + + /// Zero every entry from `n` onwards, i.e. the packed form of `self[..n]`. + #[inline] + pub const fn truncate(self, n: usize) -> Self { + if n >= Self::MAX_LEN { + self + } else { + Self(self.0 & ((1 << Self::SHIFTS[n]) - 1)) + } + } + + pub fn iter(self) -> impl DoubleEndedIterator + ExactSizeIterator { + (0..self.len()).map(move |i| self.get(i)) + } + + /// Pack `entries`, returning `None` if they do not fit. Use this for anything derived from + /// user input; use [`Self::from_slice`] when the degree bound already guarantees a fit. + pub fn try_from_slice(entries: &[PPartEntry]) -> Option { + let mut result = Self::zero(); + for (i, &entry) in entries.iter().enumerate() { + // A zero past the end is just padding, which the packed form drops anyway. + if entry == 0 { + continue; + } + if i >= Self::MAX_LEN || entry > Self::max_entry(i) { + return None; + } + result.set(i, entry); + } + Some(result) + } + + /// Pack `entries`, panicking if they do not fit. + pub fn from_slice(entries: &[PPartEntry]) -> Self { + Self::try_from_slice(entries).unwrap_or_else(|| { + panic!( + "p-part {entries:?} exceeds the degree {} bound", + Self::MAX_DEGREE + ) + }) + } +} + +impl FromIterator for PPart { + fn from_iter>(iter: I) -> Self { + let mut result = Self::zero(); + for (i, entry) in iter.into_iter().enumerate() { + if entry != 0 { + result.set(i, entry); + } + } + result + } +} + +impl std::fmt::Debug for PPart { + fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { + f.debug_list().entries(self.iter()).finish() + } +} + +/// A Milnor basis element. This is `Copy` and entirely inline: 16 bytes, no heap. +#[derive(Debug, Clone, Copy, Default)] pub struct MilnorBasisElement { pub q_part: u32, pub p_part: PPart, @@ -126,10 +327,7 @@ impl MilnorBasisElement { } pub fn clone_into(&self, other: &mut Self) { - other.q_part = self.q_part; - other.degree = self.degree; - other.p_part.clear(); - other.p_part.extend_from_slice(&self.p_part); + *other = *self; } /// Update the degree component to the correct degree @@ -138,8 +336,8 @@ impl MilnorBasisElement { let xi_degrees = combinatorics::xi_degrees(p); let tau_degrees = combinatorics::tau_degrees(p); - self.degree = q * std::iter::zip(xi_degrees, &self.p_part) - .map(|(&a, &b)| a * b as i32) + self.degree = q * std::iter::zip(xi_degrees, self.p_part.iter()) + .map(|(&a, b)| a * b as i32) .sum::() + BitflagIterator::set_bit_iterator(self.q_part as u64) .map(|k| tau_degrees[k]) @@ -161,7 +359,9 @@ impl std::cmp::Eq for MilnorBasisElement {} impl std::hash::Hash for MilnorBasisElement { fn hash(&self, state: &mut H) { - self.p_part.hash(state); + // The p-part is a single `u64`, so this is one hasher round rather than a pointer chase + // plus a variable-length slice hash. + self.p_part.bits().hash(state); #[cfg(feature = "odd-primes")] self.q_part.hash(state); } @@ -189,62 +389,13 @@ impl std::fmt::Display for MilnorBasisElement { } } -/// A version of `HashMap` that is more efficient at the prime 2. -#[cfg(feature = "odd-primes")] +/// A map from the basis elements of a single degree to their indices. +/// +/// [`MilnorBasisElement`] hashes and compares on its p-part (and, at odd primes, its q-part), both +/// of which are now single machine words, so a plain `HashMap` is already the specialised form +/// this used to hand-roll for `p = 2`. type MilnorHashMap = HashMap; -#[cfg(not(feature = "odd-primes"))] -struct MilnorHashMap { - degree: i32, - inner: HashMap, -} - -#[cfg(not(feature = "odd-primes"))] -impl Default for MilnorHashMap { - fn default() -> Self { - Self { - degree: -1, - inner: HashMap::default(), - } - } -} - -#[cfg(not(feature = "odd-primes"))] -impl MilnorHashMap { - /// Encode a [`MilnorBasisElement`] of a known degree into a `u64`. This is achieved by packing - /// the PPart into a single `u64`, where we omit the first entry since it can be derived from - /// the degree. This currently supports elements up to degree 2^9 * 3 = 1536. - fn code(x: &MilnorBasisElement) -> u64 { - let mut counter = 0; - let mut shift = 0; - for (idx, &entry) in x.p_part.iter().skip(1).enumerate() { - counter += (entry as u64) << shift; - shift += 9 - idx; - } - counter - } - - fn reserve(&mut self, additional: usize) { - self.inner.reserve(additional); - } - - fn insert(&mut self, k: MilnorBasisElement, v: V) { - if self.degree == -1 { - self.degree = k.degree; - } - assert_eq!(k.degree, self.degree); - assert!( - self.inner.insert(Self::code(&k), v).is_none(), - "Duplicate entry for {k}" - ); - } - - fn get(&self, k: &MilnorBasisElement) -> Option<&V> { - assert_eq!(k.degree, self.degree); - self.inner.get(&Self::code(k)) - } -} - pub struct MilnorAlgebra { profile: MilnorProfile, p: ValidPrime, @@ -371,7 +522,7 @@ impl Algebra for MilnorAlgebra { MilnorBasisElement { degree: 1, q_part: 1, - p_part: vec![], + p_part: PPart::zero(), }, )); } @@ -383,7 +534,7 @@ impl Algebra for MilnorAlgebra { MilnorBasisElement { degree: (2 * self.prime() - 2) as i32, q_part: 0, - p_part: vec![1], + p_part: PPart::from_iter([1]), }, )); } @@ -402,7 +553,7 @@ impl Algebra for MilnorAlgebra { MilnorBasisElement { degree, q_part: 0, - p_part: vec![1 << i], + p_part: PPart::from_iter([1 << i]), }, )); } @@ -417,6 +568,13 @@ impl Algebra for MilnorAlgebra { } fn compute_basis(&self, max_degree: i32) { + // This is the single gate that makes [`PPart`]'s packing safe: past this degree an + // exponent could outgrow its field. Everything downstream may then assume entries fit. + assert!( + max_degree <= PPart::MAX_DEGREE, + "Milnor basis elements are only supported up to degree {}, got {max_degree}", + PPart::MAX_DEGREE, + ); self.compute_ppart(max_degree); if self.generic() { @@ -432,7 +590,7 @@ impl Algebra for MilnorAlgebra { let mut map = MilnorHashMap::default(); map.reserve(basis.len()); for (i, b) in basis.iter().enumerate() { - map.insert(b.clone(), i); + assert!(map.insert(*b, i).is_none(), "Duplicate entry for {b}"); } map }); @@ -585,19 +743,29 @@ impl Algebra for MilnorAlgebra { map(char('1'), |_| Some((0, 0))), map(char('b'), |_| Some((1, 0))), map(preceded(p_or_sq, digits), |i| self.try_beps_pn(0, i)), - map((tag("P^"), digits, char('_'), digits), |(_, s, _, t)| { - let entry = p.pow(s); - let degree = entry as i32 * self.q() * combinatorics::xi_degrees(p)[t]; - let mut elt = MilnorBasisElement { - degree, - q_part: 0, - p_part: vec![0; t], - }; - elt.p_part[t - 1] = entry as PPartEntry; - self.compute_basis(degree); - self.try_basis_element_to_index(&elt) - .map(|idx| (degree, idx)) - }), + map( + (tag("P^"), digits, char('_'), digits::), + |(_, s, _, t)| { + if t == 0 || t > PPart::MAX_LEN { + return None; + } + let entry = p.pow(s) as PPartEntry; + let degree = entry as i32 * self.q() * combinatorics::xi_degrees(p)[t]; + if degree > PPart::MAX_DEGREE || entry > PPart::max_entry(t - 1) { + return None; + } + let mut p_part = PPart::zero(); + p_part.set(t - 1, entry); + let elt = MilnorBasisElement { + degree, + q_part: 0, + p_part, + }; + self.compute_basis(degree); + self.try_basis_element_to_index(&elt) + .map(|idx| (degree, idx)) + }, + ), map( ( many0(preceded(tag("Q_"), digits::)), @@ -608,12 +776,16 @@ impl Algebra for MilnorAlgebra { ), |(q_list, p_list)| { let q_part = q_list.into_iter().fold(0, |acc, q| acc + (1 << q)); + let p_part = PPart::try_from_slice(&p_list.unwrap_or_default())?; let mut elt = MilnorBasisElement { degree: 0, q_part, - p_part: p_list.unwrap_or_default(), + p_part, }; elt.compute_degree(p); + if elt.degree > PPart::MAX_DEGREE { + return None; + } self.compute_basis(elt.degree); self.try_basis_element_to_index(&elt) @@ -715,7 +887,7 @@ impl GeneratedAlgebra for MilnorAlgebra { return vec![self.basis_element_to_index(&MilnorBasisElement { degree, q_part, - p_part: vec![], + p_part: PPart::zero(), })]; } } @@ -734,7 +906,7 @@ impl GeneratedAlgebra for MilnorAlgebra { return vec![self.basis_element_to_index(&MilnorBasisElement { degree, q_part: 0, - p_part: vec![(degree as u32 / q) as PPartEntry], + p_part: PPart::from_iter([(degree as u32 / q) as PPartEntry]), })]; } vec![] @@ -753,8 +925,8 @@ impl GeneratedAlgebra for MilnorAlgebra { if self.profile.get_p_part(j as usize - 1) <= k as PPartEntry { return vec![]; } - let mut p_part = vec![0; j as usize]; - p_part[j as usize - 1] = p.pow(k) as PPartEntry; + let mut p_part = PPart::zero(); + p_part.set(j as usize - 1, p.pow(k) as PPartEntry); return vec![self.basis_element_to_index(&MilnorBasisElement { degree, q_part: 0, @@ -833,7 +1005,7 @@ impl GeneratedAlgebra for MilnorAlgebra { // Compute basis functions impl MilnorAlgebra { fn compute_ppart(&self, max_degree: i32) { - self.ppart_table.extend(0, |_| vec![Vec::new()]); + self.ppart_table.extend(0, |_| vec![PPart::zero()]); let p = self.prime().as_i32(); let q = if p == 2 { 1 } else { 2 * p - 2 }; @@ -863,20 +1035,19 @@ impl MilnorAlgebra { } let rem = (d - xi_degrees[i]) as usize; - for old in &self.ppart_table[rem] { + for &old in &self.ppart_table[rem] { // ppart_table[rem] is arranged in increasing order of highest // xi_i. If we get something too large, we may abort; if old.len() > i + 1 { break; } - if old.len() == i + 1 && old[i] == profile_list[i] { + // `profile_list[i]` is non-zero here, so `old.get(i) == profile_list[i]` + // already implies `old.len() == i + 1`. + if old.get(i) == profile_list[i] { continue; } - let mut new = old.clone(); - if new.len() < i + 1 { - new.resize(i + 1, 0); - } - new[i] += 1; + let mut new = old; + new.set(i, old.get(i) + 1); new_row.push(new); } } @@ -918,8 +1089,8 @@ impl MilnorAlgebra { table.extend( self.ppart_table[(d - q_degree as usize) / q as usize] .iter() - .map(|p_part| MilnorBasisElement { - p_part: p_part.clone(), + .map(|&p_part| MilnorBasisElement { + p_part, q_part, degree: d as i32, }), @@ -936,7 +1107,7 @@ impl MilnorAlgebra { self.basis_table.extend(max_degree as usize, |d| { let mut table: Vec<_> = self.ppart_table[d] .iter() - .map(|p| MilnorBasisElement::from_p(p.clone(), d as i32)) + .map(|&p| MilnorBasisElement::from_p(p, d as i32)) .collect(); if self.unstable_enabled { table.sort_by_cached_key(|e| e.excess(fp::prime::TWO)); @@ -973,11 +1144,14 @@ impl MilnorAlgebra { pub fn try_beps_pn(&self, e: u32, x: PPartEntry) -> Option<(i32, usize)> { let q = self.q() as u32; let degree = (q * x + e) as i32; + if degree > PPart::MAX_DEGREE || x > PPart::max_entry(0) { + return None; + } self.compute_basis(degree); self.try_basis_element_to_index(&MilnorBasisElement { degree, q_part: e, - p_part: vec![x as PPartEntry], + p_part: PPart::from_iter([x]), }) .map(|index| (degree, index)) } @@ -988,7 +1162,7 @@ impl MilnorAlgebra { } fn multiply_qpart(&self, m1: &MilnorBasisElement, f: u32) -> Vec<(u32, MilnorBasisElement)> { - let mut new_result: Vec<(u32, MilnorBasisElement)> = vec![(1, m1.clone())]; + let mut new_result: Vec<(u32, MilnorBasisElement)> = vec![(1, *m1)]; let mut old_result: Vec<(u32, MilnorBasisElement)> = Vec::new(); for k in BitflagIterator::set_bit_iterator(f as u64) { @@ -1011,23 +1185,21 @@ impl MilnorAlgebra { if term.q_part & (1 << (k + i as u32)) != 0 { continue; } - // Check if R - p^k e_i < 0. Only do this from the first term onwards. - if i > 0 && term.p_part[i - 1] < pk { - continue; - } - - let mut new_p = term.p_part.clone(); + let mut new_p = term.p_part; if i > 0 { - new_p[i - 1] -= pk; + // Check if R - p^k e_i < 0. Only do this from the first term onwards. + let entry = new_p.get(i - 1); + if entry < pk { + continue; + } + new_p.set(i - 1, entry - pk); } // Now calculate the number of Q's we are moving past let larger_q = (term.q_part >> (k + i as u32 + 1)).count_ones(); - // If new_p ends with 0, drop them - while let Some(0) = new_p.last() { - new_p.pop(); - } + // Trailing zeros are not represented in a packed p-part, so there is nothing + // to trim here. // Now put everything together let m = MilnorBasisElement { p_part: new_p, @@ -1074,8 +1246,8 @@ impl MilnorAlgebra { for (cc, basis) in m1f { let mut multiplier = PPartMultiplier::::new_from_allocation( self.prime(), - &basis.p_part, - &m2.p_part, + basis.p_part, + m2.p_part, allocation, basis.q_part, target_deg, @@ -1092,8 +1264,8 @@ impl MilnorAlgebra { } else { let mut multiplier = PPartMultiplier::::new_from_allocation( self.prime(), - &m1.p_part, - &m2.p_part, + m1.p_part, + m2.p_part, allocation, 0, target_deg, @@ -1149,7 +1321,7 @@ impl MilnorAlgebra { #[derive(Debug, Default)] struct Matrix2D { cols: usize, - inner: PPart, + inner: Vec, } impl std::fmt::Display for Matrix2D { @@ -1202,8 +1374,7 @@ impl std::ops::IndexMut for Matrix2D { pub struct PPartAllocation { m: Matrix2D, #[cfg(feature = "odd-primes")] - diagonal: PPart, - p_part: PPart, + diagonal: Vec, } thread_local! { @@ -1218,9 +1389,6 @@ impl PPartAllocation { m: Matrix2D::with_capacity(n + 1, n), #[cfg(feature = "odd-primes")] diagonal: Vec::with_capacity(n), - // This size should be the number of diagonals. Even though the answer cannot be that - // long, we still insert zeros then pop them out later. - p_part: Vec::with_capacity(2 * n), } } @@ -1232,21 +1400,21 @@ impl PPartAllocation { } #[allow(non_snake_case)] -pub struct PPartMultiplier<'a, const MOD4: bool> { +pub struct PPartMultiplier { p: ValidPrime, M: Matrix2D, - r: &'a PPart, + r: PPart, rows: usize, cols: usize, diag_num: usize, init: bool, pub ans: MilnorBasisElement, #[cfg(feature = "odd-primes")] - diagonal: PPart, + diagonal: Vec, } #[allow(non_snake_case)] -impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> { +impl PPartMultiplier { fn prime(&self) -> ValidPrime { self.p } @@ -1254,8 +1422,8 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> { #[allow(unused_mut)] // Mut is only used with odd primes pub fn new_from_allocation( p: ValidPrime, - r: &'a PPart, - s: &'a PPart, + r: PPart, + s: PPart, mut allocation: PPartAllocation, q_part: u32, degree: i32, @@ -1276,20 +1444,18 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> { M.reset(rows, cols); for i in 1..rows { - M[i][0] = r[i - 1]; + M[i][0] = r.get(i - 1); } - // This is somehow quite significantly faster than copy_from_slice - #[allow(clippy::manual_memcpy)] for k in 1..cols { - M[0][k] = s[k - 1]; + M[0][k] = s.get(k - 1); } let ans = MilnorBasisElement { q_part, - p_part: allocation.p_part, + p_part: PPart::zero(), degree, }; - PPartMultiplier { + Self { #[cfg(feature = "odd-primes")] diagonal: allocation.diagonal, p, @@ -1308,7 +1474,6 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> { m: self.M, #[cfg(feature = "odd-primes")] diagonal: self.diagonal, - p_part: self.ans.p_part, } } @@ -1392,7 +1557,7 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> { if inc <= max_inc { // If so, we found our next matrix. for row in 1..i { - self.M[row][0] = self.r[row - 1]; + self.M[row][0] = self.r.get(row - 1); for col in 1..self.cols { self.M[0][col] += self.M[row][col]; self.M[row][col] = 0; @@ -1416,13 +1581,13 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> { } } -impl Iterator for PPartMultiplier<'_, MOD4> { +impl Iterator for PPartMultiplier { type Item = u32; fn next(&mut self) -> Option { let p = self.prime().as_u32() as PPartEntry; 'outer: loop { - self.ans.p_part.clear(); + self.ans.p_part = PPart::zero(); let mut coef = 1; if self.init { @@ -1443,27 +1608,19 @@ impl Iterator for PPartMultiplier<'_, MOD4> { continue 'outer; } } - self.ans - .p_part - .reserve(std::cmp::max(self.cols, self.rows) - 1); - self.ans.p_part.extend(&self.M[0][1..self.cols]); - - if self.rows > self.cols { - self.ans.p_part.resize(self.r.len(), 0); - } - self.ans - .p_part - .iter_mut() - .zip(self.r.iter()) - .for_each(|(l, r)| *l += r); - - // If new_p ends with 0, drop them - while let Some(0) = self.ans.p_part.last() { - self.ans.p_part.pop(); + // The answer is the top row of the matrix plus `r`, entrywise. Writing a zero + // is a no-op on a packed p-part, so trailing zeros need no trimming. + for i in 0..std::cmp::max(self.cols, self.rows) - 1 { + let mut entry = self.r.get(i); + if i + 1 < self.cols { + entry += self.M[0][i + 1]; + } + if entry != 0 { + self.ans.p_part.set(i, entry); + } } return Some(coef); } else if self.update() { - self.ans.p_part.reserve(self.diag_num); for diag_idx in 1..=self.diag_num { let i_min = (diag_idx + 1).saturating_sub(self.cols); let i_max = std::cmp::min(diag_idx + 1, self.rows); @@ -1510,11 +1667,12 @@ impl Iterator for PPartMultiplier<'_, MOD4> { } } } - self.ans.p_part.push(sum); - } - // If new_p ends with 0, drop them - while let Some(0) = self.ans.p_part.last() { - self.ans.p_part.pop(); + // `diag_num` counts diagonals of the working matrix, which can exceed the + // number of entries a p-part of this degree can have; those trailing + // diagonals are necessarily zero and need not be stored. + if sum != 0 { + self.ans.p_part.set(diag_idx - 1, sum); + } } return Some(coef); @@ -1543,7 +1701,7 @@ impl MilnorAlgebra { let p_idx = self .basis_element_to_index(&MilnorBasisElement::from_p( - vec![ppow as PPartEntry], + PPart::from_iter([ppow as PPartEntry]), p_degree, )) .to_owned(); @@ -1551,7 +1709,7 @@ impl MilnorAlgebra { let q_idx = self .basis_element_to_index(&MilnorBasisElement { q_part: 1 << (i - 1), - p_part: Vec::new(), + p_part: PPart::zero(), degree: q_degree, }) .to_owned(); @@ -1567,13 +1725,13 @@ impl MilnorAlgebra { let first_idx = self.basis_element_to_index(&MilnorBasisElement { q_part: 1 << i, - p_part: Vec::new(), + p_part: PPart::zero(), degree: first_degree, }); let second_idx = self.basis_element_to_index(&MilnorBasisElement { q_part: basis.q_part ^ (1 << i), - p_part: basis.p_part.clone(), + p_part: basis.p_part, degree: second_degree, }); @@ -1607,9 +1765,9 @@ impl MilnorAlgebra { let b = self.basis_element_from_index(degree, idx); let len = b.p_part.len(); - if b.p_part[0..len - 1].iter().all(|&x| x == 0) { + if b.p_part.truncate(len - 1).is_empty() { // There is only one entry - let entry = b.p_part[len - 1]; + let entry = b.p_part.get(len - 1); let (k, m) = factor_pk(p, entry); // This is a power of p @@ -1625,12 +1783,12 @@ impl MilnorAlgebra { let l_degree = l_entry as i32 * self.q(); let l_index = self.basis_element_to_index(&MilnorBasisElement { q_part: 0, - p_part: vec![l_entry], + p_part: PPart::from_iter([l_entry]), degree: l_degree, }); - let mut r_p_part = vec![0; len - 1]; - r_p_part[len - 2] = r_entry; + let mut r_p_part = PPart::zero(); + r_p_part.set(len - 2, r_entry); let r_degree = r_entry as i32 * combinatorics::xi_degrees(p)[len - 2] * self.q(); @@ -1654,15 +1812,15 @@ impl MilnorAlgebra { let mut elt = MilnorBasisElement { q_part: 0, degree: 0, - p_part: vec![0; len], + p_part: PPart::zero(), }; - elt.p_part[len - 1] = pk; - elt.degree = entry_deg * elt.p_part[len - 1] as i32; + elt.p_part.set(len - 1, pk); + elt.degree = entry_deg * pk as i32; let first = (elt.degree, self.basis_element_to_index(&elt)); - elt.p_part[len - 1] = rem_entry; - elt.degree = entry_deg * elt.p_part[len - 1] as i32; + elt.p_part.set(len - 1, rem_entry); + elt.degree = entry_deg * rem_entry as i32; let second = (elt.degree, self.basis_element_to_index(&elt)); let coef = @@ -1671,22 +1829,18 @@ impl MilnorAlgebra { } } else { // There is more than one entry. Just separate out the last entry. - let last_entry = b.p_part[len - 1]; + let last_entry = b.p_part.get(len - 1); let last_deg = combinatorics::xi_degrees(p)[len - 1] * self.q() * last_entry as i32; let mut elt = MilnorBasisElement { q_part: 0, - p_part: vec![0; len], + p_part: PPart::zero(), degree: last_deg, }; - elt.p_part[len - 1] = last_entry; + elt.p_part.set(len - 1, last_entry); let first = (elt.degree, self.basis_element_to_index(&elt)); elt.degree = degree - last_deg; - elt.p_part.clear(); - elt.p_part.extend_from_slice(&b.p_part[0..len - 1]); - while let Some(0) = elt.p_part.last() { - elt.p_part.pop(); - } + elt.p_part = b.p_part.truncate(len - 1); let second = (elt.degree, self.basis_element_to_index(&elt)); buffer.extend([(p - c, first, second)]); }; @@ -1708,16 +1862,23 @@ impl MilnorAlgebra { } impl MilnorAlgebra { - /// Returns `true` if the new element is not within the bounds - fn increment_p_part(element: &mut PPart, max: &[PPartEntry]) -> bool { - element[0] += 1; - for i in 0..element.len() - 1 { - if element[i] > max[i] { - element[i] = 0; - element[i + 1] += 1; + /// Advance `element` to the next p-part bounded entrywise by `max`, in odometer order. + /// + /// Returns `true` once the odometer wraps, i.e. when `element` was already `max`. + /// + /// This carries *before* incrementing rather than after. The two orders enumerate the same + /// sequence, but incrementing first would transiently store `max[i] + 1`, which need not fit + /// in a packed field whose width is exactly saturated by `max[i]`. + fn increment_p_part(element: &mut PPart, max: PPart) -> bool { + for i in 0..max.len() { + let entry = element.get(i); + if entry < max.get(i) { + element.set(i, entry + 1); + return false; } + element.set(i, 0); } - element.last().unwrap() > max.last().unwrap() + true } } @@ -1730,7 +1891,7 @@ impl Bialgebra for MilnorAlgebra { let xi_degrees = combinatorics::xi_degrees(self.prime()); let mut len = 1; - let p_part = &self.basis_element_from_index(op_deg, op_idx).p_part; + let p_part = self.basis_element_from_index(op_deg, op_idx).p_part; for i in p_part.iter() { len *= i + 1; @@ -1738,32 +1899,23 @@ impl Bialgebra for MilnorAlgebra { let len = len as usize; let mut result = Vec::with_capacity(len); - let mut cur_ppart: PPart = vec![0; p_part.len()]; + let n = p_part.len(); + let mut cur_ppart = PPart::zero(); loop { let mut left_degree: i32 = 0; - for i in 0..cur_ppart.len() { - left_degree += cur_ppart[i] as i32 * xi_degrees[i]; + let mut right_ppart = PPart::zero(); + for (i, &xi_degree) in xi_degrees.iter().enumerate().take(n) { + let entry = cur_ppart.get(i); + left_degree += entry as i32 * xi_degree; + // Trailing zeros are dropped by the packing, so no trimming is needed. + right_ppart.set(i, p_part.get(i) - entry); } let right_degree: i32 = op_deg - left_degree; - let mut left_ppart = cur_ppart.clone(); - while let Some(0) = left_ppart.last() { - left_ppart.pop(); - } - - let mut right_ppart = cur_ppart - .iter() - .enumerate() - .map(|(i, v)| p_part[i] - *v) - .collect::>(); - while let Some(0) = right_ppart.last() { - right_ppart.pop(); - } - let left_idx = self.basis_element_to_index(&MilnorBasisElement { degree: left_degree, q_part: 0, - p_part: left_ppart, + p_part: cur_ppart, }); let right_idx = self.basis_element_to_index(&MilnorBasisElement { degree: right_degree, @@ -1920,13 +2072,17 @@ mod tests { assert_eq!(algebra.basis_element_to_string(d, i), name); } - // Syntactically-valid names that name no basis element must return `None` + // "P0"/"Sq0" name the identity. A packed p-part does not represent trailing zeros, so + // `P(0)` and `P()` are the same value, and `try_beps_pn(0, 0)` finds the degree-0 basis + // element. This matches `AdemAlgebra::try_beps_pn`, which special-cases `x == 0` to + // `Some((0, 0))`; the previous `None` here came from `vec![0]` and `vec![]` hashing + // differently, which was an artifact of the unpacked representation. + assert_eq!(algebra.basis_element_from_string("P0"), Some((0, 0))); + assert_eq!(algebra.basis_element_from_string("Sq0"), Some((0, 0))); + + // Syntactically-valid names that name no basis element must still return `None` // (they previously panicked in `basis_element_to_index`). // - // "P0"/"Sq0" parse via `try_beps_pn(0, 0)`, building the element - // {q_part: 0, p_part: [0]} in degree 0, which is not a basis element. - assert_eq!(algebra.basis_element_from_string("P0"), None); - assert_eq!(algebra.basis_element_from_string("Sq0"), None); // "Q_5" parses via the Q/P branch into a candidate element (degree 63) // whose basis lookup finds nothing at p = 2. assert_eq!(algebra.basis_element_from_string("Q_5"), None); @@ -2001,6 +2157,133 @@ mod tests { } } + /// 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. + #[test] + fn ppart_widths_cover_max_degree() { + let xi_degrees = combinatorics::xi_degrees(fp::prime::TWO); + for (i, &xi_degree) in xi_degrees.iter().enumerate().take(PPart::MAX_LEN) { + // deg P(R) = sum_i r_i (2^i - 1) with non-negative terms, so r_i <= deg / (2^i - 1). + let bound = PPart::MAX_DEGREE / xi_degree; + assert!( + bound <= PPart::max_entry(i) as i32, + "entry {i} needs to hold {bound} but only holds up to {}", + PPart::max_entry(i), + ); + } + // There is no entry beyond `MAX_LEN` to store: the xi-degree table itself stops there, so + // `compute_ppart` cannot produce a longer p-part. If `fp` ever raises + // `MAX_MULTINOMIAL_LEN`, this fires and `WIDTHS` has to be revisited. + assert_eq!(xi_degrees.len(), PPart::MAX_LEN); + // ... and the layout uses the whole word, so `MAX_DEGREE` is as large as it can be. + assert_eq!( + PPart::shift(PPart::MAX_LEN - 1) + PPart::width(PPart::MAX_LEN - 1), + 64 + ); + } + + #[test] + fn ppart_accessors() { + let mut p = PPart::from_slice(&[3, 0, 5]); + assert_eq!(p.len(), 3); + assert_eq!(p.iter().collect::>(), vec![3, 0, 5]); + assert_eq!(p.get(1), 0); + assert_eq!(p.get(2), 5); + // Reading past the end is zero, not a panic. + assert_eq!(p.get(7), 0); + assert_eq!(p.get(PPart::MAX_LEN), 0); + + // Trailing zeros are not represented, so they do not affect equality, length or hashing. + assert_eq!(PPart::from_slice(&[3, 0, 5, 0, 0]), p); + assert_eq!(PPart::from_slice(&[]), PPart::zero()); + assert_eq!(PPart::from_slice(&[0, 0]), PPart::zero()); + assert_eq!(PPart::zero().len(), 0); + assert!(PPart::zero().is_empty()); + + assert_eq!(p.truncate(2), PPart::from_slice(&[3])); + assert_eq!(p.truncate(0), PPart::zero()); + assert_eq!(p.truncate(PPart::MAX_LEN + 3), p); + + p.set(1, 7); + assert_eq!(p, PPart::from_slice(&[3, 7, 5])); + p.set(2, 0); + assert_eq!(p, PPart::from_slice(&[3, 7])); + } + + #[test] + fn ppart_rejects_out_of_range() { + // Too many entries, and an entry too large for its field. + assert_eq!(PPart::try_from_slice(&[1; PPart::MAX_LEN + 1]), None); + assert_eq!(PPart::try_from_slice(&[0, PPart::max_entry(1) + 1]), None); + // ... but a zero past the end is only padding. + let mut padded = vec![0; PPart::MAX_LEN + 4]; + padded[0] = 2; + assert_eq!( + PPart::try_from_slice(&padded), + Some(PPart::from_slice(&[2])) + ); + } + + #[test] + #[should_panic(expected = "does not fit")] + fn ppart_set_out_of_range_panics() { + PPart::zero().set(0, PPart::max_entry(0) + 1); + } + + /// `increment_p_part` walks up to and including `max`, whose top entry may saturate its field. + /// Incrementing before carrying would overflow there. + #[test] + fn ppart_odometer_handles_saturated_field() { + let top = PPart::MAX_LEN - 1; + let mut max = PPart::from_slice(&[2]); + max.set(top, PPart::max_entry(top)); + + let mut count = 0; + let mut cur = PPart::zero(); + loop { + count += 1; + if MilnorAlgebra::increment_p_part(&mut cur, max) { + break; + } + } + assert_eq!(count, 3 * (PPart::max_entry(top) as usize + 1)); + // Wrapping leaves the odometer back at zero. + assert_eq!(cur, PPart::zero()); + } + + /// Pack every basis element the algebra actually produces and check nothing collides or is + /// lost. This is the property the whole representation rests on. + #[rstest] + #[case(2, 120)] + #[case(3, 200)] + fn ppart_packing_is_faithful(#[case] p: u32, #[case] max_degree: i32) { + let algebra = MilnorAlgebra::new(ValidPrime::new(p), false); + algebra.compute_basis(max_degree); + + for t in 0..=max_degree { + let mut seen = HashMap::default(); + for i in 0..algebra.dimension(t) { + let elt = algebra.basis_element_from_index(t, i); + // The packed value plus the q-part identifies the element within its degree. + assert!( + seen.insert((elt.p_part.bits(), elt.q_part), i).is_none(), + "collision at degree {t} for {elt}" + ); + // Round-trip through a slice, and back through the index map. + assert_eq!( + PPart::from_slice(&elt.p_part.iter().collect::>()), + elt.p_part + ); + assert_eq!(algebra.basis_element_to_index(elt), i); + // The degree really is recoverable from the entries. + let mut recomputed = *elt; + recomputed.compute_degree(ValidPrime::new(p)); + assert_eq!(recomputed.degree, t); + } + } + } + #[test] fn test_clone_into() { let mut other = MilnorBasisElement::default(); @@ -2012,34 +2295,34 @@ mod tests { check(&MilnorBasisElement { q_part: 3, - p_part: vec![3, 2], + p_part: PPart::from_slice(&[3, 2]), degree: 12, }); check(&MilnorBasisElement { q_part: 1, - p_part: vec![3], + p_part: PPart::from_slice(&[3]), degree: 11, }); check(&MilnorBasisElement { q_part: 5, - p_part: vec![1, 3, 5, 2], + p_part: PPart::from_slice(&[1, 3, 5, 2]), degree: 7, }); check(&MilnorBasisElement { q_part: 0, - p_part: vec![], + p_part: PPart::zero(), degree: 2, }); } #[test] fn test_ppart_multiplier_2() { - let r = vec![1, 4]; - let s = vec![2, 4]; + let r = PPart::from_slice(&[1, 4]); + let s = PPart::from_slice(&[2, 4]); let mut m = PPartMultiplier::::new_from_allocation( fp::prime::TWO, - &r, - &s, + r, + s, PPartAllocation::default(), 0, 0, @@ -2075,12 +2358,12 @@ mod tests { #[test] fn test_ppart_multiplier_3() { - let r = vec![3, 4]; - let s = vec![1, 4]; + let r = PPart::from_slice(&[3, 4]); + let s = PPart::from_slice(&[1, 4]); let mut m = PPartMultiplier::::new_from_allocation( ValidPrime::new(3), - &r, - &s, + r, + s, PPartAllocation::default(), 0, 0, diff --git a/ext/crates/algebra/src/algebra/pair_algebra.rs b/ext/crates/algebra/src/algebra/pair_algebra.rs index e080d34116..23ab93ae6e 100644 --- a/ext/crates/algebra/src/algebra/pair_algebra.rs +++ b/ext/crates/algebra/src/algebra/pair_algebra.rs @@ -95,16 +95,17 @@ use std::cell::RefCell; use crate::{ MilnorAlgebra, - milnor_algebra::{MilnorBasisElement as MilnorElt, PPartAllocation, PPartMultiplier}, + milnor_algebra::{MilnorBasisElement as MilnorElt, PPart, PPartAllocation, PPartMultiplier}, }; macro_rules! sub { ($elt:ident, $k:expr, $n:expr) => { if $k > 0 { - if $elt.p_part[$k - 1] < (1 << $n) { + let entry = $elt.p_part.get($k - 1); + if entry < (1 << $n) { continue; } - $elt.p_part[$k - 1] -= 1 << $n; + $elt.p_part.set($k - 1, entry - (1 << $n)); $elt.degree -= combinatorics::xi_degrees(TWO)[$k - 1] * (1 << $n); } }; @@ -112,7 +113,7 @@ macro_rules! sub { macro_rules! unsub { ($elt:ident, $k:expr, $n:expr) => { if $k > 0 { - $elt.p_part[$k - 1] += 1 << $n; + $elt.p_part.set($k - 1, $elt.p_part.get($k - 1) + (1 << $n)); $elt.degree += combinatorics::xi_degrees(TWO)[$k - 1] * (1 << $n); } }; @@ -191,8 +192,8 @@ impl PairAlgebra for MilnorAlgebra { assert_eq!(r_degree + s_degree, result.degree); // First write the Y terms - let mut r = self.basis_element_from_index(r_degree, r_idx).clone(); - let mut s = self.basis_element_from_index(s_degree, s_idx).clone(); + let mut r = *self.basis_element_from_index(r_degree, r_idx); + let mut s = *self.basis_element_from_index(s_degree, s_idx); PPartAllocation::with_local(|mut allocation| { for k in 0..s.p_part.len() { @@ -219,8 +220,8 @@ impl PairAlgebra for MilnorAlgebra { // Now the product terms let mut multiplier = PPartMultiplier::::new_from_allocation( TWO, - &r.p_part, - &s.p_part, + r.p_part, + s.p_part, allocation, 0, r.degree + s.degree, @@ -262,7 +263,7 @@ impl PairAlgebra for MilnorAlgebra { // The twos terms for (r_idx, c) in r.iter_nonzero() { - let mut r = self.basis_element_from_index(r_degree, r_idx).clone(); + let mut r = *self.basis_element_from_index(r_degree, r_idx); sub!(r, 1, 0); self.multiply_basis_by_element( result.copy(), @@ -271,7 +272,8 @@ impl PairAlgebra for MilnorAlgebra { s_degree, s.twos.as_slice(), ); - unsub!(r, 1, 0); + // No matching `unsub!`: unlike the loops above, `r` is a fresh copy of the basis + // element on each iteration, so there is nothing to restore. } // The Y terms @@ -385,7 +387,7 @@ fn a_y_cached( None => { let v = a_y_inner(algebra, a, k, l); f(&v); - cache.insert((a.clone(), (k, l)), v); + cache.insert((*a, (k, l)), v); } } }) @@ -393,11 +395,11 @@ fn a_y_cached( /// Actually computes $A(a, Y_{k, l})$ and returns the result. fn a_y_inner(algebra: &MilnorAlgebra, a: &MilnorElt, k: usize, l: usize) -> FpVector { - let mut a = a.clone(); + let mut a = *a; let mut result = FpVector::new(TWO, algebra.dimension(a.degree + (1 << k) + (1 << l) - 2)); let mut t = MilnorElt { q_part: 0, - p_part: vec![], + p_part: PPart::zero(), degree: 0, }; @@ -410,11 +412,9 @@ fn a_y_inner(algebra: &MilnorAlgebra, a: &MilnorElt, k: usize, l: usize) -> FpVe for j in 0..=std::cmp::min(i + k - l, a.p_part.len()) { sub!(a, j, l); - t.p_part.clear(); - t.p_part.resize(k + i, 0); - - t.p_part[k + i - 1] += 1; - t.p_part[l + j - 1] += 1; + t.p_part = PPart::zero(); + t.p_part.set(k + i - 1, 1); + t.p_part.set(l + j - 1, t.p_part.get(l + j - 1) + 1); t.degree = (1 << (k + i)) + (1 << (l + j)) - 2; @@ -445,7 +445,7 @@ mod tests { MilnorElt { q_part: 0, - p_part: p_part.into(), + p_part: PPart::from_slice(p_part), degree, } } diff --git a/ext/crates/algebra/src/module/rpn.rs b/ext/crates/algebra/src/module/rpn.rs index 6d42a1f9f8..5d2b3f43f0 100644 --- a/ext/crates/algebra/src/module/rpn.rs +++ b/ext/crates/algebra/src/module/rpn.rs @@ -170,7 +170,7 @@ fn coef_milnor(algebra: &MilnorAlgebra, op_deg: i32, op_idx: usize, mut mod_degr let mut list = Vec::with_capacity(elt.p_part.len() + 1); list.push(mod_degree - sum); - list.extend_from_slice(&elt.p_part); + list.extend(elt.p_part.iter()); PPartEntry::multinomial2(&list) == 1 } diff --git a/ext/crates/algebra/src/steenrod_evaluator.rs b/ext/crates/algebra/src/steenrod_evaluator.rs index c5b9c56135..27239cb626 100644 --- a/ext/crates/algebra/src/steenrod_evaluator.rs +++ b/ext/crates/algebra/src/steenrod_evaluator.rs @@ -8,7 +8,7 @@ use fp::{ use crate::{ algebra::{AdemAlgebra, Algebra, MilnorAlgebra, adem_algebra::AdemBasisElement}, - milnor_algebra::{MilnorBasisElement, PPartEntry}, + milnor_algebra::{MilnorBasisElement, PPart, PPartEntry}, steenrod_parser::*, }; @@ -157,7 +157,7 @@ impl SteenrodEvaluator { * q; let elt = MilnorBasisElement { degree, - p_part: p_list, + p_part: PPart::from_slice(&p_list), q_part: 0, }; @@ -270,9 +270,9 @@ impl SteenrodEvaluator { return; } let mut t: Vec = vec![0; elt.p_part.len()]; - t[elt.p_part.len() - 1] = elt.p_part[elt.p_part.len() - 1]; + t[elt.p_part.len() - 1] = elt.p_part.get(elt.p_part.len() - 1); for i in (0..elt.p_part.len() - 1).rev() { - t[i] = elt.p_part[i] + 2 * t[i + 1]; + t[i] = elt.p_part.get(i) + 2 * t[i + 1]; } let t_idx = self.adem.basis_element_to_index(&AdemBasisElement { degree, @@ -307,19 +307,10 @@ impl SteenrodEvaluator { (31u32.saturating_sub(elt.q_part.leading_zeros())) as usize, ); let mut t = vec![0; t_len]; - let last_p_part = if t_len <= elt.p_part.len() { - elt.p_part[t_len - 1] - } else { - 0 - }; - t[t_len - 1] = last_p_part + ((elt.q_part >> (t_len)) & 1); + // `PPart::get` already reads past the end as zero. + t[t_len - 1] = elt.p_part.get(t_len - 1) + ((elt.q_part >> (t_len)) & 1); for i in (0..t_len - 1).rev() { - let p_part = if i < elt.p_part.len() { - elt.p_part[i] - } else { - 0 - }; - t[i] = p_part + ((elt.q_part >> (i + 1)) & 1) + p * t[i + 1]; + t[i] = elt.p_part.get(i) + ((elt.q_part >> (i + 1)) & 1) + p * t[i + 1]; } let t_idx = self.adem.basis_element_to_index(&AdemBasisElement { degree, @@ -344,11 +335,11 @@ impl SteenrodEvaluator { MilnorBasisElement { degree, q_part: 1 << qi, - p_part: vec![], + p_part: PPart::zero(), } } else { - let mut p_part = vec![0; qi as usize + 1]; - p_part[qi as usize] = 1; + let mut p_part = PPart::zero(); + p_part.set(qi as usize, 1); MilnorBasisElement { degree, q_part: 0, diff --git a/ext/crates/algebra/src/steenrod_parser.rs b/ext/crates/algebra/src/steenrod_parser.rs index 2d03718653..a5f719571f 100644 --- a/ext/crates/algebra/src/steenrod_parser.rs +++ b/ext/crates/algebra/src/steenrod_parser.rs @@ -15,14 +15,14 @@ use nom::{ sequence::{delimited, pair, preceded}, }; -use crate::{adem_algebra::AdemBasisElement, algebra::milnor_algebra::PPart}; +use crate::{adem_algebra::AdemBasisElement, algebra::milnor_algebra::PPartEntry}; type IResult = IResultBase>; #[derive(Debug, Clone)] pub enum AlgebraBasisElt { AList(Vec), // Admissible list. - PList(PPart), + PList(Vec), P(u32), Q(u32), } diff --git a/ext/examples/bruner.rs b/ext/examples/bruner.rs index 71c5c0e22a..6e4411b1fd 100644 --- a/ext/examples/bruner.rs +++ b/ext/examples/bruner.rs @@ -25,7 +25,7 @@ use std::{ use algebra::{ Algebra, MilnorAlgebra, - milnor_algebra::MilnorBasisElement, + milnor_algebra::{MilnorBasisElement, PPartEntry}, module::{FreeModule as FM, Module, homomorphism::FreeModuleHomomorphism as FMH}, }; use anyhow::{Context, Error, Result}; @@ -95,7 +95,10 @@ fn get_algebra_element<'a>( let entry = &entry[1..]; let elt = MilnorBasisElement { q_part: 0, - p_part: entry.split(',').map(|x| x.parse().unwrap()).collect(), + p_part: entry + .split(',') + .map(|x| x.parse::().unwrap()) + .collect(), degree: t, }; a.basis_element_to_index(&elt) diff --git a/ext/examples/sq0.rs b/ext/examples/sq0.rs index dc671c159b..9836e18c00 100644 --- a/ext/examples/sq0.rs +++ b/ext/examples/sq0.rs @@ -77,8 +77,9 @@ mod double { mod double_algebra { use algebra::{ - AdemAlgebra, Algebra, MilnorAlgebra, SteenrodAlgebra, adem_algebra::AdemBasisElement, - milnor_algebra::MilnorBasisElement, + AdemAlgebra, Algebra, MilnorAlgebra, SteenrodAlgebra, + adem_algebra::AdemBasisElement, + milnor_algebra::{MilnorBasisElement, PPart}, }; pub trait DoubleAlgebra: Algebra { @@ -92,14 +93,14 @@ mod double { let p_part = elt .p_part .iter() - .map(|&x| { + .map(|x| { if x.is_multiple_of(2) { Some(x / 2) } else { None } }) - .collect::>>()?; + .collect::>()?; Some(self.basis_element_to_index(&MilnorBasisElement { degree: degree / 2, p_part, diff --git a/ext/src/nassau.rs b/ext/src/nassau.rs index 8d1919136a..5135aabeec 100644 --- a/ext/src/nassau.rs +++ b/ext/src/nassau.rs @@ -20,7 +20,7 @@ use std::{ use algebra::{ Algebra, combinatorics, - milnor_algebra::{MilnorAlgebra, PPartEntry}, + milnor_algebra::{MilnorAlgebra, PPart, PPartEntry}, module::{ FreeModule, GeneratorData, Module, ZeroModule, homomorphism::{FreeModuleHomomorphism, FullModuleHomomorphism, ModuleHomomorphism}, @@ -100,15 +100,23 @@ impl MilnorSubalgebra { Self { profile: vec![] } } - /// Computes the signature of an element - fn has_signature(&self, ppart: &[PPartEntry], signature: &[PPartEntry]) -> bool { - for (i, (&profile, &signature)) in self.profile.iter().zip(signature).enumerate() { - let ppart = ppart.get(i).copied().unwrap_or(0); - if ppart & ((1 << profile) - 1) != signature { - return false; - } + /// The test "does this element have this signature" compiled into a `(mask, value)` pair to + /// match against the packed p-part. + /// + /// The per-entry test is `ppart[i] & ((1 << profile[i]) - 1) == signature[i]`. Because each + /// entry occupies a fixed field of the packed word, the low `profile[i]` bits of entry `i` are + /// a fixed bit range of that word, so the whole conjunction is a single `&` and `==`. Entries + /// past the end of the p-part read as zero, which the packing already gives us for free. + fn packed_signature(&self, signature: &[PPartEntry]) -> (u64, u64) { + let mut mask = 0; + let mut value = 0; + for (i, (&profile, &entry)) in self.profile.iter().zip(signature).enumerate() { + // A profile wider than the field constrains the whole field. + let width = std::cmp::min(profile as u32, PPart::width(i)); + mask |= ((1u64 << width) - 1) << PPart::shift(i); + value |= (entry as u64) << PPart::shift(i); } - true + (mask, value) } fn zero_signature(&self) -> Vec { @@ -131,12 +139,15 @@ impl MilnorSubalgebra { start: [offset], end: _, }| { + // Hoist the mask out of the inner loop: every element in this block is tested + // against the same signature. + let (mask, value) = self.packed_signature(signature); algebra .ppart_table(degree - gen_deg) .iter() .enumerate() .filter_map(move |(n, op)| { - if self.has_signature(op, signature) { + if op.bits() & mask == value { Some(offset + n) } else { None diff --git a/ext/src/yoneda.rs b/ext/src/yoneda.rs index 2fffe713c9..c2abd2b23f 100644 --- a/ext/src/yoneda.rs +++ b/ext/src/yoneda.rs @@ -60,7 +60,7 @@ fn rate_milnor_operation(algebra: &MilnorAlgebra, deg: i32, idx: usize) -> i32 { elt.p_part .iter() .enumerate() - .map(|(i, &r)| r.count_ones() << i) + .map(|(i, r)| r.count_ones() << i) .sum::() as i32 } From 87d7bcc4a44a1fa5bd1fd105547dd3dd42703fc9 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 31 Jul 2026 05:59:15 +0000 Subject: [PATCH 2/6] Cut per-entry overhead out of the Milnor multiplier The first packing pass regressed `milnor_ppart` by up to 8% at odd primes and mod 4, because assembling the answer went from a memcpy plus a vectorized add to a per-entry read-modify-write through the checked `PPart::set`, and because `PPart::get`'s range branch landed in `update`'s inner loop. Two changes, both confined to the kernel: - Assemble the answer in a plain `u64` and store it once. Entries are written in increasing index order into a value that starts at zero, so a shift and an `or` suffice; the range checks become debug assertions backed by `compute_basis`'s degree gate. - Pad the layout tables to 16 entries so the private `PPart::entry` can mask its index rather than branch on it. Padded entries have width zero and so read as zero, which is the answer `get` would have returned anyway. The public `get` keeps its explicit check, since callers outside the multiplier index it with a q-part-derived length that is not bounded by `MAX_LEN`. This recovers the regression (`ppart_4/a` and `ppart_3/a` back to baseline, `ppart_4/b` -8%) and improves the Nassau regime further. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV --- .../algebra/src/algebra/milnor_algebra.rs | 70 ++++++++++++++----- 1 file changed, 52 insertions(+), 18 deletions(-) diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index bbf780041c..3359b1b3e3 100644 --- a/ext/crates/algebra/src/algebra/milnor_algebra.rs +++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs @@ -138,13 +138,18 @@ impl PPart { /// already bounds the length of the $\xi$-degree table, so it is not a new restriction. pub const MAX_LEN: usize = 10; + /// The length of the layout tables. This is [`Self::MAX_LEN`] rounded up to a power of two so + /// that [`Self::entry`] can mask its index instead of bounds-checking it; entries at or past + /// `MAX_LEN` are given width 0, so they read as zero. + const TABLE_LEN: usize = 16; + /// `WIDTHS[i]` is the number of bits holding $r_{i+1}$: the number of bits needed to represent /// `MAX_DEGREE / (2^(i+1) - 1)`. - const WIDTHS: [u32; Self::MAX_LEN] = [11, 10, 9, 8, 7, 6, 5, 4, 3, 1]; + const WIDTHS: [u32; Self::TABLE_LEN] = [11, 10, 9, 8, 7, 6, 5, 4, 3, 1, 0, 0, 0, 0, 0, 0]; /// `SHIFTS[i]` is the bit offset of entry `i`; `SHIFTS[MAX_LEN]` is the total width, 64. - const SHIFTS: [u32; Self::MAX_LEN + 1] = { - let mut shifts = [0; Self::MAX_LEN + 1]; + const SHIFTS: [u32; Self::TABLE_LEN] = { + let mut shifts = [0; Self::TABLE_LEN]; let mut i = 0; while i < Self::MAX_LEN { shifts[i + 1] = shifts[i] + Self::WIDTHS[i]; @@ -195,18 +200,41 @@ impl PPart { /// The raw packed value. Two exponent sequences are equal exactly when their bits are, so this /// is a complete hash key, and it can be compared against a packed mask in one operation (see - /// `MilnorSubalgebra::has_signature` in `ext`). + /// `MilnorSubalgebra::packed_signature` in `ext`). pub const fn bits(self) -> u64 { self.0 } - /// Entry `i`, or 0 if `i` is past the end. + /// Reinterpret a raw packed value. + /// + /// Callers that assemble entries by shifting must uphold the type invariant themselves: each + /// entry must lie within its field, which holds for any element of degree at most + /// [`Self::MAX_DEGREE`]. This exists so hot loops can accumulate into a plain `u64` and store + /// once, rather than read-modify-write through [`Self::set`] per entry. + pub(crate) const fn from_bits(bits: u64) -> Self { + Self(bits) + } + + /// Entry `i`, for `i < TABLE_LEN`, with no bounds check. + /// + /// Masking the index keeps the table lookups in range without a branch. Entries in + /// `MAX_LEN..TABLE_LEN` have width 0 and so read as zero, which is the right answer; an index + /// at or beyond `TABLE_LEN` would silently wrap, which is why this is private and + /// `debug_assert`ed. Callers in the multiplier are all bounded by `MAX_LEN`. + #[inline] + const fn entry(self, i: usize) -> PPartEntry { + debug_assert!(i < Self::TABLE_LEN); + let i = i & (Self::TABLE_LEN - 1); + ((self.0 >> Self::SHIFTS[i]) & ((1 << Self::WIDTHS[i]) - 1)) as PPartEntry + } + + /// Entry `i`, or 0 if `i` is past the end. Accepts any index. #[inline] pub const fn get(self, i: usize) -> PPartEntry { if i >= Self::MAX_LEN { return 0; } - ((self.0 >> Self::SHIFTS[i]) & ((1 << Self::WIDTHS[i]) - 1)) as PPartEntry + self.entry(i) } /// Set entry `i` to `v`. @@ -1444,10 +1472,10 @@ impl PPartMultiplier { M.reset(rows, cols); for i in 1..rows { - M[i][0] = r.get(i - 1); + M[i][0] = r.entry(i - 1); } for k in 1..cols { - M[0][k] = s.get(k - 1); + M[0][k] = s.entry(k - 1); } let ans = MilnorBasisElement { @@ -1557,7 +1585,7 @@ impl PPartMultiplier { if inc <= max_inc { // If so, we found our next matrix. for row in 1..i { - self.M[row][0] = self.r.get(row - 1); + self.M[row][0] = self.r.entry(row - 1); for col in 1..self.cols { self.M[0][col] += self.M[row][col]; self.M[row][col] = 0; @@ -1587,7 +1615,6 @@ impl Iterator for PPartMultiplier { fn next(&mut self) -> Option { let p = self.prime().as_u32() as PPartEntry; 'outer: loop { - self.ans.p_part = PPart::zero(); let mut coef = 1; if self.init { @@ -1608,19 +1635,22 @@ impl Iterator for PPartMultiplier { continue 'outer; } } - // The answer is the top row of the matrix plus `r`, entrywise. Writing a zero - // is a no-op on a packed p-part, so trailing zeros need no trimming. + // The answer is the top row of the matrix plus `r`, entrywise. Accumulate into + // a plain word and store once; trailing zeros contribute nothing, so there is no + // trimming to do. + let mut ans = 0; for i in 0..std::cmp::max(self.cols, self.rows) - 1 { - let mut entry = self.r.get(i); + let mut entry = self.r.entry(i); if i + 1 < self.cols { entry += self.M[0][i + 1]; } - if entry != 0 { - self.ans.p_part.set(i, entry); - } + debug_assert!(entry <= PPart::max_entry(i)); + ans |= (entry as u64) << PPart::shift(i); } + self.ans.p_part = PPart::from_bits(ans); return Some(coef); } else if self.update() { + let mut ans = 0; for diag_idx in 1..=self.diag_num { let i_min = (diag_idx + 1).saturating_sub(self.cols); let i_max = std::cmp::min(diag_idx + 1, self.rows); @@ -1670,10 +1700,14 @@ impl Iterator for PPartMultiplier { // `diag_num` counts diagonals of the working matrix, which can exceed the // number of entries a p-part of this degree can have; those trailing // diagonals are necessarily zero and need not be stored. - if sum != 0 { - self.ans.p_part.set(diag_idx - 1, sum); + if diag_idx <= PPart::MAX_LEN { + debug_assert!(sum <= PPart::max_entry(diag_idx - 1)); + ans |= (sum as u64) << PPart::shift(diag_idx - 1); + } else { + debug_assert_eq!(sum, 0); } } + self.ans.p_part = PPart::from_bits(ans); return Some(coef); } else { From accef7c06dc7cb34e30c7d1393035e0d38bfa104 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 31 Jul 2026 17:32:41 +0000 Subject: [PATCH 3/6] Derive the p=2 Milnor basis instead of storing it `basis_table` held a `MilnorBasisElement` per basis element. At p = 2 with unstable support off, that element is exactly `from_p(ppart_table[t][i], t)` -- the p-part again, with a q-part that is always zero and a degree that is the index. It was a redundant copy. Deriving it on demand costs nothing now that `MilnorBasisElement` is `Copy` and 16 bytes: `basis_element_from_index` returns by value and builds it in registers rather than handing out a reference into a table. The multiply family takes the element by value for the same reason. The table is still built at odd primes, where the q-part varies within a degree, and when unstable support is on, where the basis is re-sorted by excess. Neither is a re-wrapping of `ppart_table`. Measured over degrees 0..=250 at p = 2 (1,958,958 elements), RSS growth from `compute_basis` drops 125.0 MB -> 95.0 MB, i.e. 66.9 -> 50.8 bytes per element. Projected to degree 500 that is 5.24 GB -> 3.95 GB. Unlike the ranker, this needs no basis renumbering and costs nothing at lookup time. A test verifies the derivation matches what the table used to hold, for every element, so the redundancy is asserted rather than assumed. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV --- .../algebra/src/algebra/milnor_algebra.rs | 124 ++++++++++++++---- .../algebra/src/algebra/pair_algebra.rs | 26 ++-- ext/crates/algebra/src/module/rpn.rs | 2 +- 3 files changed, 114 insertions(+), 38 deletions(-) diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index 3359b1b3e3..ad55e988b8 100644 --- a/ext/crates/algebra/src/algebra/milnor_algebra.rs +++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs @@ -436,7 +436,14 @@ pub struct MilnorAlgebra { /// degree `q * i`. ppart_table: OnceVec>, - /// A list of all basis elements of each degree, constructed from [`Self::ppart_table`] + /// A list of all basis elements of each degree, constructed from [`Self::ppart_table`]. + /// + /// Only populated when [`Self::stores_basis_table`] holds. At `p = 2` with unstable support + /// off, the basis element at index `i` of degree `t` is exactly + /// `MilnorBasisElement::from_p(ppart_table[t][i], t)`, so storing it repeats the p-part with a + /// known q-part and degree bolted on -- 16 bytes per element, about a quarter of the algebra's + /// footprint. [`Self::basis_element_from_index`] reconstructs it instead, which is free now + /// that the type is `Copy` and fits in registers. basis_table: OnceVec>, excess_table: OnceVec>, @@ -502,8 +509,21 @@ impl MilnorAlgebra { &self.profile } - pub fn basis_element_from_index(&self, degree: i32, idx: usize) -> &MilnorBasisElement { - &self.basis_table[degree as usize][idx] + /// Whether the basis of each degree has to be stored rather than derived. + /// + /// At odd primes the q-part varies within a degree, and with unstable support enabled the + /// basis is re-sorted by excess; in both cases the basis is not a re-wrapping of + /// [`Self::ppart_table`] and must be kept. + fn stores_basis_table(&self) -> bool { + self.generic() || self.unstable_enabled + } + + pub fn basis_element_from_index(&self, degree: i32, idx: usize) -> MilnorBasisElement { + if self.stores_basis_table() { + self.basis_table[degree as usize][idx] + } else { + MilnorBasisElement::from_p(self.ppart_table[degree as usize][idx], degree) + } } pub fn try_basis_element_to_index(&self, elt: &MilnorBasisElement) -> Option { @@ -614,11 +634,12 @@ impl Algebra for MilnorAlgebra { // Populate hash map self.basis_element_to_index_map .extend(max_degree as usize, |d| { - let basis = &self.basis_table[d]; let mut map = MilnorHashMap::default(); - map.reserve(basis.len()); - for (i, b) in basis.iter().enumerate() { - assert!(map.insert(*b, i).is_none(), "Duplicate entry for {b}"); + let dim = self.dimension(d as i32); + map.reserve(dim); + for i in 0..dim { + let b = self.basis_element_from_index(d as i32, i); + assert!(map.insert(b, i).is_none(), "Duplicate entry for {b}"); } map }); @@ -639,8 +660,8 @@ impl Algebra for MilnorAlgebra { self.multiply( res.as_slice_mut(), 1, - &self.basis_table[d][i], - &self.basis_table[e][j], + &self.basis_element_from_index(d as i32, i), + &self.basis_element_from_index(e as i32, j), ); res }) @@ -660,7 +681,11 @@ impl Algebra for MilnorAlgebra { if degree < 0 { return 0; } - self.basis_table[degree as usize].len() + if self.stores_basis_table() { + self.basis_table[degree as usize].len() + } else { + self.ppart_table[degree as usize].len() + } } #[cfg(not(feature = "cache-multiplication"))] @@ -837,7 +862,7 @@ impl UnstableAlgebra for MilnorAlgebra { } else if excess < degree { self.excess_table[degree as usize][excess as usize] } else { - self.basis_table[degree as usize].len() + self.dimension(degree) } } @@ -1132,14 +1157,16 @@ impl MilnorAlgebra { } fn generate_basis_2(&self, max_degree: i32) { + if !self.stores_basis_table() { + // Derived on demand from `ppart_table`; see the field docs. + return; + } self.basis_table.extend(max_degree as usize, |d| { let mut table: Vec<_> = self.ppart_table[d] .iter() .map(|&p| MilnorBasisElement::from_p(p, d as i32)) .collect(); - if self.unstable_enabled { - table.sort_by_cached_key(|e| e.excess(fp::prime::TWO)); - } + table.sort_by_cached_key(|e| e.excess(fp::prime::TWO)); table }); } @@ -1189,8 +1216,8 @@ impl MilnorAlgebra { self.try_beps_pn(e, x).unwrap() } - fn multiply_qpart(&self, m1: &MilnorBasisElement, f: u32) -> Vec<(u32, MilnorBasisElement)> { - let mut new_result: Vec<(u32, MilnorBasisElement)> = vec![(1, *m1)]; + fn multiply_qpart(&self, m1: MilnorBasisElement, f: u32) -> Vec<(u32, MilnorBasisElement)> { + let mut new_result: Vec<(u32, MilnorBasisElement)> = vec![(1, m1)]; let mut old_result: Vec<(u32, MilnorBasisElement)> = Vec::new(); for k in BitflagIterator::set_bit_iterator(f as u64) { @@ -1251,8 +1278,8 @@ impl MilnorAlgebra { &self, res: FpSliceMut, coef: u32, - m1: &MilnorBasisElement, - m2: &MilnorBasisElement, + m1: MilnorBasisElement, + m2: MilnorBasisElement, ) { PPartAllocation::with_local(|allocation| { self.multiply_with_allocation(res, coef, m1, m2, i32::MAX, allocation) @@ -1263,8 +1290,8 @@ impl MilnorAlgebra { &self, mut res: FpSliceMut, coef: u32, - m1: &MilnorBasisElement, - m2: &MilnorBasisElement, + m1: MilnorBasisElement, + m2: MilnorBasisElement, excess: i32, mut allocation: PPartAllocation, ) -> PPartAllocation { @@ -1314,7 +1341,7 @@ impl MilnorAlgebra { &self, res: FpSliceMut, coef: u32, - m1: &MilnorBasisElement, + m1: MilnorBasisElement, s_deg: i32, s: FpSlice, ) { @@ -1327,7 +1354,7 @@ impl MilnorAlgebra { &self, mut res: FpSliceMut, coef: u32, - m1: &MilnorBasisElement, + m1: MilnorBasisElement, s_deg: i32, s: FpSlice, mut allocation: PPartAllocation, @@ -2309,15 +2336,64 @@ mod tests { PPart::from_slice(&elt.p_part.iter().collect::>()), elt.p_part ); - assert_eq!(algebra.basis_element_to_index(elt), i); + assert_eq!(algebra.basis_element_to_index(&elt), i); // The degree really is recoverable from the entries. - let mut recomputed = *elt; + let mut recomputed = elt; recomputed.compute_degree(ValidPrime::new(p)); assert_eq!(recomputed.degree, t); } } } + /// At `p = 2` with unstable support off, the basis is not stored: it is derived from + /// `ppart_table`. Check the derivation reproduces exactly what the table used to hold, so the + /// redundancy this relies on is asserted rather than assumed. + #[test] + fn basis_is_derived_at_p2() { + let p = fp::prime::TWO; + let algebra = MilnorAlgebra::new(p, false); + algebra.compute_basis(120); + assert!( + !algebra.stores_basis_table(), + "p = 2 stable should not be storing the basis" + ); + + for t in 0..=120 { + let pparts = algebra.ppart_table(t); + assert_eq!(algebra.dimension(t), pparts.len()); + for (i, &p_part) in pparts.iter().enumerate() { + // This is precisely what `generate_basis_2` used to store. + let expected = MilnorBasisElement { + q_part: 0, + p_part, + degree: t, + }; + let actual = algebra.basis_element_from_index(t, i); + assert_eq!(actual.p_part, expected.p_part, "degree {t}, index {i}"); + assert_eq!(actual.q_part, expected.q_part, "degree {t}, index {i}"); + assert_eq!(actual.degree, expected.degree, "degree {t}, index {i}"); + } + } + } + + /// The two configurations that still need the table really do differ from `ppart_table`, so + /// the exemption in `stores_basis_table` is not over-broad. + #[rstest] + #[case(3, false)] + #[case(2, true)] + fn basis_is_stored_when_it_must_be(#[case] p: u32, #[case] unstable: bool) { + let algebra = MilnorAlgebra::new(ValidPrime::new(p), unstable); + algebra.compute_basis(60); + assert!(algebra.stores_basis_table()); + // Every stored element still round-trips through the index map. + for t in 0..=60 { + for i in 0..algebra.dimension(t) { + let elt = algebra.basis_element_from_index(t, i); + assert_eq!(algebra.basis_element_to_index(&elt), i); + } + } + } + #[test] fn test_clone_into() { let mut other = MilnorBasisElement::default(); diff --git a/ext/crates/algebra/src/algebra/pair_algebra.rs b/ext/crates/algebra/src/algebra/pair_algebra.rs index 23ab93ae6e..34815baaa0 100644 --- a/ext/crates/algebra/src/algebra/pair_algebra.rs +++ b/ext/crates/algebra/src/algebra/pair_algebra.rs @@ -192,8 +192,8 @@ impl PairAlgebra for MilnorAlgebra { assert_eq!(r_degree + s_degree, result.degree); // First write the Y terms - let mut r = *self.basis_element_from_index(r_degree, r_idx); - let mut s = *self.basis_element_from_index(s_degree, s_idx); + let mut r = self.basis_element_from_index(r_degree, r_idx); + let mut s = self.basis_element_from_index(s_degree, s_idx); PPartAllocation::with_local(|mut allocation| { for k in 0..s.p_part.len() { @@ -205,8 +205,8 @@ impl PairAlgebra for MilnorAlgebra { allocation = self.multiply_with_allocation( result.ys[m + k][n + k].as_slice_mut(), coeff, - &r, - &s, + r, + s, i32::MAX, allocation, ); @@ -263,12 +263,12 @@ impl PairAlgebra for MilnorAlgebra { // The twos terms for (r_idx, c) in r.iter_nonzero() { - let mut r = *self.basis_element_from_index(r_degree, r_idx); + let mut r = self.basis_element_from_index(r_degree, r_idx); sub!(r, 1, 0); self.multiply_basis_by_element( result.copy(), coeff * c, - &r, + r, s_degree, s.twos.as_slice(), ); @@ -366,7 +366,7 @@ thread_local! { /// [`a_y_inner`] if not available. fn a_y_cached( algebra: &MilnorAlgebra, - a: &MilnorElt, + a: MilnorElt, k: usize, l: usize, f: impl FnOnce(&FpVector), @@ -379,7 +379,7 @@ fn a_y_cached( let raw_entry = cache.raw_entry(); let result = raw_entry - .from_hash(hasher.finish(), |v| &v.0 == a && v.1 == (k, l)) + .from_hash(hasher.finish(), |v| v.0 == a && v.1 == (k, l)) .map(|(_, y)| y); match result { @@ -387,15 +387,15 @@ fn a_y_cached( None => { let v = a_y_inner(algebra, a, k, l); f(&v); - cache.insert((*a, (k, l)), v); + cache.insert((a, (k, l)), v); } } }) } /// Actually computes $A(a, Y_{k, l})$ and returns the result. -fn a_y_inner(algebra: &MilnorAlgebra, a: &MilnorElt, k: usize, l: usize) -> FpVector { - let mut a = *a; +fn a_y_inner(algebra: &MilnorAlgebra, a: MilnorElt, k: usize, l: usize) -> FpVector { + let mut a = a; let mut result = FpVector::new(TWO, algebra.dimension(a.degree + (1 << k) + (1 << l) - 2)); let mut t = MilnorElt { q_part: 0, @@ -420,7 +420,7 @@ fn a_y_inner(algebra: &MilnorAlgebra, a: &MilnorElt, k: usize, l: usize) -> FpVe // We can just read off the value of the product instead of passing through the // algorithm, but this is cached so problem for another day... - algebra.multiply(result.as_slice_mut(), 1, &t, &a); + algebra.multiply(result.as_slice_mut(), 1, t, a); unsub!(a, j, l); } @@ -462,7 +462,7 @@ mod tests { let target_deg = a.degree + (1 << k) + (1 << l) - 2; algebra.compute_basis(target_deg + 1); result.set_scratch_vector_size(algebra.dimension(target_deg)); - a_y_cached(&algebra, &a, k, l, |v| result.add(v, 1)); + a_y_cached(&algebra, a, k, l, |v| result.add(v, 1)); ans.assert_eq(&algebra.element_to_string(target_deg, result.as_slice())); }; diff --git a/ext/crates/algebra/src/module/rpn.rs b/ext/crates/algebra/src/module/rpn.rs index 5d2b3f43f0..1c4aab1d0b 100644 --- a/ext/crates/algebra/src/module/rpn.rs +++ b/ext/crates/algebra/src/module/rpn.rs @@ -157,7 +157,7 @@ fn coef_milnor(algebra: &MilnorAlgebra, op_deg: i32, op_idx: usize, mut mod_degr return false; } - let elt: &MilnorBasisElement = algebra.basis_element_from_index(op_deg, op_idx); + let elt: MilnorBasisElement = algebra.basis_element_from_index(op_deg, op_idx); let sum: PPartEntry = elt.p_part.iter().sum(); if mod_degree < 0 { From 6334813614ebe09a737301c5266fe3e4f589c66d Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 31 Jul 2026 17:32:41 +0000 Subject: [PATCH 4/6] Add PPartRanker behind an opt-in feature, not wired in `basis_element_to_index` runs once per term of every product, and is a hash map storing an entry per basis element. The index it returns is a position in an enumeration, so a canonical key alone cannot replace the map -- but the position can be computed. Let counts[i][d] be the number of exponent sequences of degree d using only xi_1..xi_i. Splitting on whether r_i is zero gives the coin-change recurrence counts[i][d] = counts[i-1][d] + counts[i][d - xi_i]. Ranking needs the number of sequences with r_i > v, and substituting r_i -> r_i - (v+1) is a bijection onto all sequences of degree d - (v+1)*xi_i, so that count is a single table lookup rather than a sum. Walking the entries downward ranks a p-part in one lookup each, against a table covering every degree at once, where the map it would replace grows with the basis. Whether that is worth it depends entirely on scale, which took some measuring to see. Against the map, at p = 2: degree per-degree map hashmap ranker ratio 120 0.10 MB 11.6us 26.7us 0.43x 300 3.12 MB 792us 1384us 0.57x 400 12.50 MB 4490us 5310us 0.85x 500 37.50 MB 33118us 15865us 2.09x A lookup probes only its own degree's map. While that fits in cache the map wins easily: one hash round and one probe, against six to ten dependent table reads. Once it does not -- the map is 37 MB in degree 500 -- every probe misses to DRAM at ~33 ns, whereas the ranker's table is ~43 KB, stays in L1, and costs ~16 ns regardless of degree. Benchmarking only up to degree 120, where the map is 0.1 MB, shows a 2x loss and hides all of this; the sweep here deliberately spans the crossover. Tuning does not move the small-degree end: nested vs flat table, a zero-padded prefix to drop the branch, and one- vs two-pass to break the dependency chain were all measured, and the padded variant was worst, because doubling the table pushed it out of L1. So the two suit opposite ends of the range, and the ranker is on the right side of the end where the algebra's memory is the problem worth solving: replacing the map there is 3.3 GB smaller and 2x faster. It stays off by default and unwired even when enabled, because it numbers the basis in colex order rather than the order compute_ppart emits, which would invalidate saved resolutions. That order is rankable in principle, but its natural recursion has depth equal to the sum of the entries, which is worse than hashing. The unstable path, which re-sorts each degree by excess, is not modelled either. Tests verify the table reproduces the algebra's own p-part counts and that the rank is a bijection onto 0..dim in every degree, at p = 2 and p = 3, plus one pinning down that it really does disagree with the current basis order. Also measured and rejected: an `unrank` recovering the p-part at a given index, which would let basis_element_from_index drop ppart_table entirely. It ran ~15x slower than the array read it would replace, at every degree, with none of the crossover above -- ppart_table is 8 bytes per element against ~43 for the map, so it stays cache-resident. The bit-packing that makes rank worth having is the same thing that makes unrank not. The likelier route, if it is ever revisited, is enumerating the basis in index order, which is O(1) amortised and matches how callers actually walk it. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV --- ext/crates/algebra/Cargo.toml | 9 + ext/crates/algebra/benches/milnor_rank.rs | 100 +++++++ ext/crates/algebra/src/algebra/milnor_rank.rs | 265 ++++++++++++++++++ ext/crates/algebra/src/algebra/mod.rs | 5 + 4 files changed, 379 insertions(+) create mode 100644 ext/crates/algebra/benches/milnor_rank.rs create mode 100644 ext/crates/algebra/src/algebra/milnor_rank.rs diff --git a/ext/crates/algebra/Cargo.toml b/ext/crates/algebra/Cargo.toml index 3b98825cb5..f40a7c1cc5 100644 --- a/ext/crates/algebra/Cargo.toml +++ b/ext/crates/algebra/Cargo.toml @@ -37,6 +37,10 @@ rstest = "0.25.0" [features] default = ["odd-primes"] cache-multiplication = [] +# An arithmetic replacement for the Milnor basis index map. Off by default and not wired in: +# adopting it would renumber the basis and invalidate saved resolutions. See +# `algebra::milnor_rank`. +milnor-rank = [] concurrent = ["fp/concurrent", "maybe-rayon/concurrent"] odd-primes = ["fp/odd-primes"] @@ -55,3 +59,8 @@ harness = false [[bench]] name = "nassau_milnor" harness = false + +[[bench]] +name = "milnor_rank" +harness = false +required-features = ["milnor-rank"] diff --git a/ext/crates/algebra/benches/milnor_rank.rs b/ext/crates/algebra/benches/milnor_rank.rs new file mode 100644 index 0000000000..bd063f2017 --- /dev/null +++ b/ext/crates/algebra/benches/milnor_rank.rs @@ -0,0 +1,100 @@ +//! Compares [`PPartRanker`] against the hash map lookup it would replace. +//! +//! `MilnorAlgebra::basis_element_to_index` is called once per term of every product, so it is one +//! of the hottest operations in a resolution. It is currently a hash map from the (packed) basis +//! element to its position in `basis_table`. The ranker computes that position arithmetically +//! instead, from a table that covers every degree at once. +//! +//! The two are compared on the same workload: recover the index of every basis element of a +//! degree. Note that they do not agree on *which* index — see [`PPartRanker`] — so this measures +//! the cost of the two strategies, not a drop-in substitution. +//! +//! Each is measured under two access orders, because the choice flatters the map: +//! +//! - **sequential** — sweep the basis in order. This is the map's insertion order, so every probe +//! walks memory linearly and prefetches perfectly. Flattering, and not what callers do. +//! - **scattered** — the same elements in a fixed pseudo-random permutation. This is closer to +//! real use, where `basis_element_to_index` is called on multiplication *outputs*, which arrive +//! in no particular order. It matters because the two structures scale differently: the map +//! stores an entry per basis element and leaves cache as the basis grows (~60 KiB in degree 120 +//! alone), whereas the ranker's table is a few KiB covering every degree at once. +//! +//! [`PPartRanker`]: algebra::milnor_rank::PPartRanker + +use std::hint::black_box; + +use algebra::{Algebra, MilnorAlgebra, milnor_rank::PPartRanker}; +use criterion::{Criterion, Throughput, criterion_group, criterion_main}; +use fp::prime::TWO; +use pprof::criterion::{Output, PProfProfiler}; + +/// Degrees to sweep. +/// +/// The range matters more than it looks, because the two structures live in different parts of the +/// memory hierarchy and the crossover is inside this range. A lookup probes only its own degree's +/// map, which is ~0.1 MB in degree 120 (L2-resident) but ~3 MB in degree 300 and ~12 MB in degree +/// 400 — well past L3, so every probe is a DRAM miss. The ranker's table is ~35 KB for *all* +/// degrees and stays in L1 throughout. Measuring only the small degrees answers a question nobody +/// is asking; the large ones are where expanding the algebra actually hurts. +const DEGREES: &[i32] = &[120, 300, 400, 500]; + +/// A fixed permutation of `0..n`, from a Fisher-Yates shuffle driven by a small LCG. Deterministic +/// so the two variants see exactly the same access order. +fn scattered(n: usize) -> Vec { + let mut order: Vec = (0..n).collect(); + let mut state = 0x2545_f491_4f6c_dd1d_u64; + for i in (1..n).rev() { + state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + order.swap(i, (state >> 33) as usize % (i + 1)); + } + order +} + +fn milnor_rank(c: &mut Criterion) { + let algebra = MilnorAlgebra::new(TWO, false); + let max_degree = *DEGREES.iter().max().unwrap(); + algebra.compute_basis(max_degree); + let ranker = PPartRanker::new(TWO, max_degree); + + let mut g = c.benchmark_group("milnor_rank"); + for °ree in DEGREES { + let dim = algebra.dimension(degree); + g.throughput(Throughput::Elements(dim as u64)); + + // Collect the elements once so neither variant pays for the table walk itself. + let elements: Vec<_> = (0..dim) + .map(|i| algebra.basis_element_from_index(degree, i)) + .collect(); + let shuffled: Vec<_> = scattered(dim).into_iter().map(|i| elements[i]).collect(); + + for (order, elements) in [("seq", &elements), ("scattered", &shuffled)] { + g.bench_function(format!("hashmap_{order}/deg{degree}"), |b| { + b.iter(|| { + for elt in elements { + black_box(algebra.basis_element_to_index(elt)); + } + }); + }); + + g.bench_function(format!("ranker_{order}/deg{degree}"), |b| { + b.iter(|| { + for elt in elements { + black_box(ranker.rank(elt.p_part, degree)); + } + }); + }); + } + } + 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 = milnor_rank +} +criterion_main!(benches); diff --git a/ext/crates/algebra/src/algebra/milnor_rank.rs b/ext/crates/algebra/src/algebra/milnor_rank.rs new file mode 100644 index 0000000000..e0b592e19b --- /dev/null +++ b/ext/crates/algebra/src/algebra/milnor_rank.rs @@ -0,0 +1,265 @@ +//! An arithmetic replacement for the `MilnorBasisElement -> index` hash map. +//! +//! Behind the off-by-default `milnor-rank` feature, and **not wired into [`MilnorAlgebra`]** even +//! when enabled -- `basis_element_to_index` still goes through the hash map. It is here so the +//! design and its measurements survive, ready to switch on when the renumbering below is worth +//! taking on. +//! +//! **This is not wired into [`MilnorAlgebra`].** It computes a *different* numbering of each +//! degree's basis than [`MilnorAlgebra::compute_basis`] produces, so adopting it would renumber +//! the basis and invalidate every saved resolution. See [`PPartRanker`] for why the numbering +//! cannot simply be made to match, and the crate benchmarks (`milnor_rank`) for what it costs +//! relative to the hash map it would replace. + +use fp::prime::ValidPrime; + +use crate::algebra::{combinatorics, milnor_algebra::PPart}; + +/// Computes a p-part's index within its degree arithmetically, instead of looking it up. +/// +/// # How it works +/// +/// `counts[i][d]` is the number of exponent sequences of degree `d` (in units of `q`) that use +/// only $\xi_1, \ldots, \xi_i$. Splitting on whether $r_i$ is zero gives the coin-change +/// recurrence +/// +/// ```text +/// counts[i][d] = counts[i - 1][d] + counts[i][d - xi_i] +/// ``` +/// +/// so the table costs `O(MAX_LEN * max_degree)` to build, and `counts[MAX_LEN][d]` is the +/// dimension of the algebra in degree `d`. +/// +/// Ranking then rests on one identity. Among the sequences of degree `d` using +/// $\xi_1, \ldots, \xi_i$, those with $r_i \ge u$ are in bijection with *all* sequences of degree +/// `d - u * xi_i` using $\xi_1, \ldots, \xi_i$, via $r_i \mapsto r_i - u$. So the number of them +/// with $r_i > v$ is `counts[i][d - (v + 1) * xi_i]`: a single lookup, with no summation. Walking +/// the entries from the top down therefore ranks a p-part in [`PPart::MAX_LEN`] lookups and adds, +/// against a table of a few hundred KiB that serves *every* degree at once — where the hash map +/// it replaces stores one entry per basis element. +/// +/// # Speed: it depends entirely on scale +/// +/// Measured against the hash map it would replace (`cargo bench --bench milnor_rank`), at p = 2: +/// +/// ```text +/// degree per-degree map hashmap ranker ratio +/// 120 0.10 MB 11.6us 26.7us 0.43x +/// 300 3.12 MB 792us 1384us 0.57x +/// 400 12.50 MB 4490us 5310us 0.85x +/// 500 37.50 MB 33118us 15865us 2.09x +/// ``` +/// +/// The crossover is a cache effect, and it is not marginal. A lookup probes only its own degree's +/// map. While that fits in cache the map wins easily: one hash round and one probe, against six to +/// ten dependent table reads for a rank. Once it does not — the map is 37 MB in degree 500 — every +/// probe misses to DRAM at ~33 ns, whereas the ranker's whole table is ~43 KB, stays in L1, and +/// costs ~16 ns regardless of degree. The map's cost grows with the basis; the ranker's does not. +/// +/// So the ranker is the wrong tool for small degrees and the right one for large. That matters +/// because large degrees are exactly where the algebra's memory becomes the problem worth solving. +/// +/// The reverse direction does not share this, and is deliberately absent. An `unrank` -- recovering +/// the p-part at a given index, which would let `basis_element_from_index` drop `ppart_table` +/// altogether -- was written, tested and benchmarked, and removed again: it ran ~15x slower than +/// the array read it would replace, at every degree measured, with no sign of the crossover that +/// makes `rank` worthwhile. The reason is that `ppart_table` is only 8 bytes per element, against +/// ~43 for the hash map, so it stays cache-resident where the map does not. The bit-packing that +/// makes `rank` worth having is the same thing that makes `unrank` not. (See the history around +/// "Speed up unrank 1.5x" if it needs revisiting; the likelier route is enumerating the basis in +/// index order, which is O(1) amortised and matches how callers actually walk it.) +/// +/// Tuning does not move this. Five arrangements were measured — nested vs flat count table, with +/// and without a zero-padded prefix to drop the branch, and one- vs two-pass to break the +/// dependency chain. None changed the small-degree verdict; the padded variant was worst, because +/// doubling the table pushed it out of L1. +/// +/// # Why it is still not wired in +/// +/// Not speed, but numbering. This ranks in colex order on $(r_{10}, \ldots, r_1)$. `compute_ppart` emits a different order: +/// it groups by the highest non-zero entry, and recurses by decrementing one entry at a time. +/// That order *is* rankable in principle, but its natural recursion has depth $\sum_i r_i$ — up to +/// `MAX_DEGREE` — which is far worse than hashing. Getting the `O(MAX_LEN)` cost requires adopting +/// the colex order, i.e. renumbering the basis. +/// +/// A renumbering is not intrinsically hard — [`crate::Algebra::magic`] already exists to +/// discriminate save files — but it invalidates stored resolutions, so it is a migration rather +/// than a drop-in change. Note also that [`MilnorAlgebra`] re-sorts each degree by excess when +/// unstable support is enabled, which this does not model. +/// +/// [`MilnorAlgebra`]: crate::MilnorAlgebra +/// [`MilnorAlgebra::compute_basis`]: crate::Algebra::compute_basis +pub struct PPartRanker { + /// `counts[i][d]` flattened to `counts[i * stride + d]`, for `i` in `0..=PPart::MAX_LEN`. + /// + /// Flat rather than `Vec>`, so a lookup is not a dependent pointer chase. + counts: Vec, + stride: usize, + /// `effective_len[d]` is the number of $\xi_i$ of degree at most `d`. + /// + /// Entries beyond it cannot contribute to a rank in degree `d`: such an entry must be zero, and + /// its `cut` is then `d - xi_i < 0`. At degree 120 this is 6 rather than 10, so it removes + /// roughly a third of the work. + effective_len: Vec, + /// `xi[i]` is the degree of $\xi_{i+1}$, divided by `q`. + xi: [i32; PPart::MAX_LEN], + max_degree: i32, +} + +impl PPartRanker { + /// Build the table for degrees `0..=max_degree`, where `max_degree` is measured in units of + /// `q` (so it is the internal degree at `p = 2`, and the internal degree divided by + /// `2(p - 1)` otherwise). + pub fn new(p: ValidPrime, max_degree: i32) -> Self { + assert!(max_degree >= 0); + let mut xi = [0; PPart::MAX_LEN]; + xi.copy_from_slice(&combinatorics::xi_degrees(p)[..PPart::MAX_LEN]); + + let stride = max_degree as usize + 1; + + let mut counts = vec![0; (PPart::MAX_LEN + 1) * stride]; + // The empty sequence is the unique sequence of degree 0 using no generators. + counts[0] = 1; + for i in 1..=PPart::MAX_LEN { + for d in 0..stride { + // Either r_i is zero, or we can subtract one from it. + counts[i * stride + d] = counts[(i - 1) * stride + d]; + if d >= xi[i - 1] as usize { + counts[i * stride + d] += counts[i * stride + d - xi[i - 1] as usize]; + } + } + } + + let effective_len = (0..=max_degree) + .map(|d| xi.iter().filter(|&&x| x <= d).count() as u8) + .collect(); + + Self { + counts, + stride, + effective_len, + xi, + max_degree, + } + } + + /// The number of p-parts of degree `degree`, i.e. what `MilnorAlgebra::dimension` returns for + /// the p-part factor of the basis. + pub fn dimension(&self, degree: i32) -> u64 { + if degree < 0 || degree > self.max_degree { + 0 + } else { + self.counts[PPart::MAX_LEN * self.stride + degree as usize] + } + } + + /// The index of `p_part` among the p-parts of degree `degree`, in the colex order described on + /// [`PPartRanker`]. + /// + /// `degree` must be the degree of `p_part` (in units of `q`), and at most the `max_degree` + /// this was built with. + /// + /// # Cost + /// + /// One table read per entry, serialised through the running `remaining`. That is the reason + /// this loses to the hash map it was meant to replace, and no arrangement of the table fixes + /// it — see the module docs. + #[inline] + pub fn rank(&self, p_part: PPart, degree: i32) -> usize { + debug_assert!(degree >= 0 && degree <= self.max_degree); + let mut rank = 0; + let mut remaining = degree; + // Only the entries with `xi_i <= degree` can contribute; the rest are zero with a negative + // cut. At degree 120 that is 6 iterations rather than 10. + for i in (0..self.effective_len[degree as usize] as usize).rev() { + let entry = p_part.get(i) as i32; + // Everything with a larger entry here sorts earlier, and there are exactly + // `counts[i + 1][remaining - (entry + 1) * xi_i]` of them. + let cut = remaining - (entry + 1) * self.xi[i]; + if cut >= 0 { + rank += self.counts[(i + 1) * self.stride + cut as usize]; + } + remaining -= entry * self.xi[i]; + } + debug_assert_eq!(remaining, 0, "degree does not match the p-part"); + rank as usize + } +} + +#[cfg(test)] +mod tests { + use fp::prime::Prime; + use rstest::rstest; + + use super::*; + use crate::{Algebra, MilnorAlgebra}; + + /// `counts[MAX_LEN]` must agree with the algebra's own count of p-parts in each degree. + #[rstest] + #[case(2, 120)] + #[case(3, 40)] + fn table_matches_ppart_table(#[case] p: u32, #[case] max_degree: i32) { + let p = ValidPrime::new(p); + let algebra = MilnorAlgebra::new(p, false); + let q = if p == 2 { 1 } else { 2 * (p.as_i32() - 1) }; + algebra.compute_basis(max_degree * q); + + let ranker = PPartRanker::new(p, max_degree); + for d in 0..=max_degree { + assert_eq!( + ranker.dimension(d), + algebra.ppart_table(d).len() as u64, + "dimension mismatch in degree {d}" + ); + } + } + + /// The whole point: `rank` must be a bijection from the p-parts of each degree onto + /// `0..dimension`. If it is, it is a valid numbering and could replace the hash map. + #[rstest] + #[case(2, 120)] + #[case(3, 40)] + fn rank_is_a_bijection(#[case] p: u32, #[case] max_degree: i32) { + let p = ValidPrime::new(p); + let algebra = MilnorAlgebra::new(p, false); + let q = if p == 2 { 1 } else { 2 * (p.as_i32() - 1) }; + algebra.compute_basis(max_degree * q); + + let ranker = PPartRanker::new(p, max_degree); + for d in 0..=max_degree { + let table = algebra.ppart_table(d); + let mut seen = vec![false; table.len()]; + for &p_part in table { + let rank = ranker.rank(p_part, d); + assert!(rank < table.len(), "rank {rank} out of range in degree {d}"); + assert!(!seen[rank], "rank {rank} hit twice in degree {d}"); + seen[rank] = true; + } + } + } + + /// Ranking in colex order really is a different numbering than the one the algebra uses. This + /// is the reason the ranker is not wired in, so pin it down rather than leave it to prose. + #[test] + fn rank_disagrees_with_the_current_basis_order() { + let p = ValidPrime::new(2); + let algebra = MilnorAlgebra::new(p, false); + algebra.compute_basis(120); + let ranker = PPartRanker::new(p, 120); + + let mut agree = 0; + let mut total = 0; + for d in 0..=120 { + for (i, &p_part) in algebra.ppart_table(d).iter().enumerate() { + total += 1; + if ranker.rank(p_part, d) == i { + agree += 1; + } + } + } + assert!( + agree * 100 < total, + "expected the orders to differ on almost everything, but {agree}/{total} agreed" + ); + } +} diff --git a/ext/crates/algebra/src/algebra/mod.rs b/ext/crates/algebra/src/algebra/mod.rs index 67baed2b04..c2425320af 100644 --- a/ext/crates/algebra/src/algebra/mod.rs +++ b/ext/crates/algebra/src/algebra/mod.rs @@ -18,6 +18,11 @@ pub use field::Field; pub mod milnor_algebra; pub use milnor_algebra::MilnorAlgebra; +/// Opt-in: an arithmetic alternative to the Milnor basis index map. Not wired in; see the module +/// docs for what it costs and what it would take to adopt. +#[cfg(feature = "milnor-rank")] +pub mod milnor_rank; + mod steenrod_algebra; pub use steenrod_algebra::{AlgebraType, SteenrodAlgebra}; From d09141e0f629576e3cb6e4e200cd0fc8633523ce Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 1 Aug 2026 18:06:52 +0000 Subject: [PATCH 5/6] Reject inputs the packing cannot represent, instead of panicking Three paths computed with unvalidated input before checking it against the packing bounds, so the intermediate arithmetic went wrong first. All three are reachable from public, non-panicking entry points. - `basis_element_from_string("P^s_t")` indexed the xi-degree table with `t`, which has exactly `MAX_LEN` entries, so `t = MAX_LEN` was out of bounds. `p^s` and the degree product could also overflow. Now `t` is bounded by the table itself and both are computed with checked arithmetic. - `try_beps_pn` computed `q * x + e` before bounding `x`, which overflows for a large `x`. The bound moves above the computation. - `MilnorSubalgebra::packed_signature` assumed the profile was no longer than `PPart::MAX_LEN` and that each signature entry fit its field. Neither holds: `SubalgebraIterator` grows a profile without limit and `from_bytes` reads whatever length a file gives. Out of range, `PPart::shift` returns 64 and the shift overflowed; an oversized entry silently spilled into the neighbouring field, which could select unrelated basis elements. It now returns `None` for a signature no element can have, and `signature_mask` yields nothing. `basis_element_from_string` is documented as total and `try_beps_pn` is the non-panicking half of `beps_pn`, so these were contract violations rather than merely untidy. Tests cover each. The signature test checks the packed mask against the per-entry comparison it replaced, over every element up to degree 60, for profiles that are narrower than their fields, wider than their fields, and longer than a p-part can be. Also adds `basis_order_at_p2_is_stable`, which pins the first nine degrees to fixed element names. The basis order is a wire format -- saved resolutions store coefficients by index -- so it needs a guard that does not read from `ppart_table`, which is the thing being guarded. Verified separately that the order is unchanged from the base commit: identical for all 4156 elements in degrees 0..=60. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV --- .../algebra/src/algebra/milnor_algebra.rs | 91 ++++++++++++-- ext/crates/algebra/src/algebra/milnor_rank.rs | 12 +- ext/crates/algebra/src/algebra/mod.rs | 4 +- ext/src/nassau.rs | 119 ++++++++++++++---- 4 files changed, 186 insertions(+), 40 deletions(-) diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index ad55e988b8..956a015d2b 100644 --- a/ext/crates/algebra/src/algebra/milnor_algebra.rs +++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs @@ -106,7 +106,7 @@ pub type PPartEntry = u32; /// The exponent sequence $(r_1, r_2, \ldots)$ of a Milnor basis element $P(r_1, r_2, \ldots)$, /// bit-packed into a single `u64`. /// -/// Entry $r_{i+1}$ occupies [`Self::WIDTHS`]`[i]` bits starting at bit [`Self::SHIFTS`]`[i]`. The +/// Entry $r_{i+1}$ occupies `WIDTHS[i]` bits starting at bit `SHIFTS[i]` (both private). The /// widths are forced by the degree bound: at $p = 2$ the internal degree of $P(R)$ is /// $\sum_i r_i (2^i - 1)$ and every term is non-negative, so an element of degree at most /// [`Self::MAX_DEGREE`] has $r_i \le \mathrm{MAX\\_DEGREE}/(2^i - 1)$. At an odd prime the same @@ -129,7 +129,7 @@ pub struct PPart(u64); impl PPart { /// The largest internal degree whose exponent sequences are guaranteed to fit. /// - /// This is the largest bound for which [`Self::WIDTHS`] sums to at most 64. It is far beyond + /// This is the largest bound for which the field widths sum to at most 64. It is far beyond /// anything reachable — the Milnor algebra already has over 5 million basis elements below /// degree 300 — and exceeds the degree 1536 the previous hand-rolled packing assumed. pub const MAX_DEGREE: i32 = 2045; @@ -660,8 +660,8 @@ impl Algebra for MilnorAlgebra { self.multiply( res.as_slice_mut(), 1, - &self.basis_element_from_index(d as i32, i), - &self.basis_element_from_index(e as i32, j), + self.basis_element_from_index(d as i32, i), + self.basis_element_from_index(e as i32, j), ); res }) @@ -799,12 +799,17 @@ impl Algebra for MilnorAlgebra { map( (tag("P^"), digits, char('_'), digits::), |(_, s, _, t)| { - if t == 0 || t > PPart::MAX_LEN { + if t == 0 || t >= combinatorics::xi_degrees(p).len() { return None; } - let entry = p.pow(s) as PPartEntry; - let degree = entry as i32 * self.q() * combinatorics::xi_degrees(p)[t]; - if degree > PPart::MAX_DEGREE || entry > PPart::max_entry(t - 1) { + let entry: PPartEntry = p.as_u32().checked_pow(s)?; + if entry > PPart::max_entry(t - 1) { + return None; + } + let degree = (entry as i32) + .checked_mul(self.q())? + .checked_mul(combinatorics::xi_degrees(p)[t])?; + if degree > PPart::MAX_DEGREE { return None; } let mut p_part = PPart::zero(); @@ -1197,9 +1202,14 @@ impl MilnorAlgebra { /// Return the degree and index of $Q_1^e P(x)$, or `None` if the element is not present /// (e.g. out of range or excluded by the profile). pub fn try_beps_pn(&self, e: u32, x: PPartEntry) -> Option<(i32, usize)> { + // Bound `x` first: `q * x + e` overflows for a large `x`, so the degree cannot be + // computed before it has been rejected. + if x > PPart::max_entry(0) { + return None; + } let q = self.q() as u32; let degree = (q * x + e) as i32; - if degree > PPart::MAX_DEGREE || x > PPart::max_entry(0) { + if degree > PPart::MAX_DEGREE { return None; } self.compute_basis(degree); @@ -2119,6 +2129,35 @@ mod tests { assert_eq!(a2.try_beps_pn(0, 8), None); } + /// `basis_element_from_string` is documented to be total. Inputs whose exponents overflow + /// intermediate arithmetic must return `None`, not panic. + #[test] + fn basis_element_from_string_rejects_overflowing_exponents() { + let algebra = MilnorAlgebra::new(fp::prime::TWO, false); + algebra.compute_basis(8); + + // `t` indexes the xi-degree table, which has exactly `MAX_LEN` entries. + assert_eq!(algebra.basis_element_from_string("P^1_10"), None); + assert_eq!(algebra.basis_element_from_string("P^1_99"), None); + // `p^s` overflows for large `s`. + assert_eq!(algebra.basis_element_from_string("P^64_2"), None); + assert_eq!(algebra.basis_element_from_string("P^4294967295_2"), None); + // ... and so does `q * x` in the `Sq`/`P` path. + assert_eq!(algebra.basis_element_from_string("Sq4294967295"), None); + } + + /// `try_beps_pn` is the non-panicking half of `beps_pn`; an out-of-range `x` must not trip + /// overflow on the way to the bounds check. + #[test] + fn try_beps_pn_rejects_overflowing_x() { + for p in [2, 3] { + let algebra = MilnorAlgebra::new(ValidPrime::new(p), false); + assert_eq!(algebra.try_beps_pn(0, PPartEntry::MAX), None); + assert_eq!(algebra.try_beps_pn(0, PPartEntry::MAX / 2), None); + assert_eq!(algebra.try_beps_pn(1, PPartEntry::MAX), None); + } + } + #[test] fn basis_element_from_string_total_milnor() { let p = ValidPrime::new(2); @@ -2345,6 +2384,40 @@ mod tests { } } + /// The basis *order* at `p = 2` is a wire format: saved resolutions store coefficients by + /// index, so reordering silently invalidates them without `magic()` changing. Deriving the + /// basis from `ppart_table` preserves the order `generate_basis_2` produced, since the stable + /// path never sorted. Pin that down against fixed expected names, so the check does not depend + /// on `ppart_table` -- the very thing it is guarding. + #[test] + fn basis_order_at_p2_is_stable() { + let algebra = MilnorAlgebra::new(fp::prime::TWO, false); + algebra.compute_basis(8); + + let expected: [&[&str]; 9] = [ + &["1"], + &["P(1)"], + &["P(2)"], + &["P(3)", "P(0, 1)"], + &["P(4)", "P(1, 1)"], + &["P(5)", "P(2, 1)"], + &["P(6)", "P(3, 1)", "P(0, 2)"], + &["P(7)", "P(4, 1)", "P(1, 2)", "P(0, 0, 1)"], + &["P(8)", "P(5, 1)", "P(2, 2)", "P(1, 0, 1)"], + ]; + for (t, names) in expected.iter().enumerate() { + let t = t as i32; + assert_eq!(algebra.dimension(t), names.len(), "dimension in degree {t}"); + for (i, name) in names.iter().enumerate() { + assert_eq!( + &algebra.basis_element_to_string(t, i), + name, + "degree {t}, index {i}" + ); + } + } + } + /// At `p = 2` with unstable support off, the basis is not stored: it is derived from /// `ppart_table`. Check the derivation reproduces exactly what the table used to hold, so the /// redundancy this relies on is asserted rather than assumed. diff --git a/ext/crates/algebra/src/algebra/milnor_rank.rs b/ext/crates/algebra/src/algebra/milnor_rank.rs index e0b592e19b..65c6f0a7e4 100644 --- a/ext/crates/algebra/src/algebra/milnor_rank.rs +++ b/ext/crates/algebra/src/algebra/milnor_rank.rs @@ -5,11 +5,13 @@ //! design and its measurements survive, ready to switch on when the renumbering below is worth //! taking on. //! -//! **This is not wired into [`MilnorAlgebra`].** It computes a *different* numbering of each -//! degree's basis than [`MilnorAlgebra::compute_basis`] produces, so adopting it would renumber -//! the basis and invalidate every saved resolution. See [`PPartRanker`] for why the numbering -//! cannot simply be made to match, and the crate benchmarks (`milnor_rank`) for what it costs -//! relative to the hash map it would replace. +//! It computes a *different* numbering of each degree's basis than [`MilnorAlgebra::compute_basis`] +//! produces, so adopting it would renumber the basis and invalidate every saved resolution. See +//! [`PPartRanker`] for why the numbering cannot simply be made to match, and the crate benchmarks +//! (`milnor_rank`) for what it costs relative to the hash map it would replace. +//! +//! [`MilnorAlgebra`]: crate::MilnorAlgebra +//! [`MilnorAlgebra::compute_basis`]: crate::Algebra::compute_basis use fp::prime::ValidPrime; diff --git a/ext/crates/algebra/src/algebra/mod.rs b/ext/crates/algebra/src/algebra/mod.rs index c2425320af..6627432dcd 100644 --- a/ext/crates/algebra/src/algebra/mod.rs +++ b/ext/crates/algebra/src/algebra/mod.rs @@ -18,8 +18,8 @@ pub use field::Field; pub mod milnor_algebra; pub use milnor_algebra::MilnorAlgebra; -/// Opt-in: an arithmetic alternative to the Milnor basis index map. Not wired in; see the module -/// docs for what it costs and what it would take to adopt. +// Opt-in: an arithmetic alternative to the Milnor basis index map. Not wired in; see the module +// docs for what it costs and what it would take to adopt. #[cfg(feature = "milnor-rank")] pub mod milnor_rank; diff --git a/ext/src/nassau.rs b/ext/src/nassau.rs index 5135aabeec..9c847781ea 100644 --- a/ext/src/nassau.rs +++ b/ext/src/nassau.rs @@ -107,16 +107,33 @@ impl MilnorSubalgebra { /// entry occupies a fixed field of the packed word, the low `profile[i]` bits of entry `i` are /// a fixed bit range of that word, so the whole conjunction is a single `&` and `==`. Entries /// past the end of the p-part read as zero, which the packing already gives us for free. - fn packed_signature(&self, signature: &[PPartEntry]) -> (u64, u64) { + /// Returns `None` if no element can have this signature, which the caller turns into an empty + /// result. A profile is not bounded by [`PPart::MAX_LEN`] -- `SubalgebraIterator` grows one + /// without limit and `from_bytes` reads whatever length a file gives -- so both ways a + /// signature can fail to be representable have to be handled here rather than assumed away. + fn packed_signature(&self, signature: &[PPartEntry]) -> Option<(u64, u64)> { let mut mask = 0; let mut value = 0; for (i, (&profile, &entry)) in self.profile.iter().zip(signature).enumerate() { + if i >= PPart::MAX_LEN { + // No p-part of a representable degree has an entry this far out, so it reads as + // zero: a non-zero constraint is unsatisfiable and a zero one is vacuous. + if entry != 0 { + return None; + } + continue; + } // A profile wider than the field constrains the whole field. let width = std::cmp::min(profile as u32, PPart::width(i)); + // The masked entry has only `width` bits, so a signature wanting more matches nothing. + // Packing it anyway would spill into the neighbouring field. + if (entry as u64) >> width != 0 { + return None; + } mask |= ((1u64 << width) - 1) << PPart::shift(i); value |= (entry as u64) << PPart::shift(i); } - (mask, value) + Some((mask, value)) } fn zero_signature(&self) -> Vec { @@ -133,28 +150,31 @@ impl MilnorSubalgebra { degree: i32, signature: &'a [PPartEntry], ) -> impl Iterator + 'a { - module.iter_gen_offsets([degree]).flat_map( - move |GeneratorData { - gen_deg, - start: [offset], - end: _, - }| { - // Hoist the mask out of the inner loop: every element in this block is tested - // against the same signature. - let (mask, value) = self.packed_signature(signature); - algebra - .ppart_table(degree - gen_deg) - .iter() - .enumerate() - .filter_map(move |(n, op)| { - if op.bits() & mask == value { - Some(offset + n) - } else { - None - } - }) - }, - ) + // The mask depends only on the signature, so compute it once for the whole sweep. An + // unrepresentable signature yields no elements at all. + self.packed_signature(signature) + .into_iter() + .flat_map(move |(mask, value)| { + module.iter_gen_offsets([degree]).flat_map( + move |GeneratorData { + gen_deg, + start: [offset], + end: _, + }| { + algebra + .ppart_table(degree - gen_deg) + .iter() + .enumerate() + .filter_map(move |(n, op)| { + if op.bits() & mask == value { + Some(offset + n) + } else { + None + } + }) + }, + ) + }) } /// Get the matrix of a free module homomorphism when restricted to the subquotient given by @@ -1336,4 +1356,55 @@ mod tests { vec![vec![0, 1, 0, 0], vec![0, 2, 0, 0], vec![0, 0, 1, 0],] ); } + + /// The packed signature test must agree with the per-entry comparison it replaced, including + /// on signatures that no element can have. Packing those naively would spill bits into the + /// neighbouring field and select unrelated elements. + #[test] + fn packed_signature_matches_per_entry_test() { + // The comparison the packed mask replaced, kept here as the reference. + fn has_signature(profile: &[u8], ppart: PPart, signature: &[PPartEntry]) -> bool { + for (i, (&profile, &signature)) in profile.iter().zip(signature).enumerate() { + if ppart.get(i) & ((1u64 << profile) - 1) as PPartEntry != signature { + return false; + } + } + true + } + + let algebra = MilnorAlgebra::new(TWO, false); + algebra.compute_basis(60); + + for profile in [ + vec![1u8, 1, 1], + vec![4, 3, 2, 1], + vec![2, 0, 3], + // Wider than the fields they constrain. + vec![9, 9, 9, 9], + // Longer than a p-part can be, so the tail entries can never be non-zero. + vec![1; PPart::MAX_LEN + 3], + ] { + let subalgebra = MilnorSubalgebra::new(profile.clone()); + for signature in [ + vec![0; profile.len()], + (0..profile.len()).map(|i| (i % 3) as PPartEntry).collect(), + // An entry too wide for its field, which must match nothing. + (0..profile.len()) + .map(|i| if i == profile.len() - 1 { 255 } else { 0 }) + .collect(), + ] { + let packed = subalgebra.packed_signature(&signature); + for t in 0..=60 { + for &op in algebra.ppart_table(t) { + let expected = has_signature(&profile, op, &signature); + let actual = packed.is_some_and(|(mask, value)| op.bits() & mask == value); + assert_eq!( + actual, expected, + "profile {profile:?}, signature {signature:?}, element {op:?}" + ); + } + } + } + } + } } From b4097d7390c2673e1d646fdf9293008dd7f0c76c Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 1 Aug 2026 18:11:26 +0000 Subject: [PATCH 6/6] Apply nightly rustfmt to the packed p-part constants CI lints with the nightly toolchain, where the unstable options in `rustfmt.toml` -- `reorder_impl_items` among them -- actually take effect. Stable rustfmt skips them with a warning, so `cargo fmt --check` passed locally and failed in CI. Formatting only: the constants are sorted and the blank lines between them dropped. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV --- .../algebra/src/algebra/milnor_algebra.rs | 49 +++++++++---------- 1 file changed, 22 insertions(+), 27 deletions(-) diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index 956a015d2b..cd72a42a56 100644 --- a/ext/crates/algebra/src/algebra/milnor_algebra.rs +++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs @@ -127,26 +127,30 @@ pub type PPartEntry = u32; pub struct PPart(u64); impl PPart { + /// `FIELD_OF_BIT[b]` is the index of the entry owning bit `b`, letting [`Self::len`] turn a + /// `leading_zeros` into an entry index without looping. + const FIELD_OF_BIT: [u8; 64] = { + let mut table = [0; 64]; + let mut i = 0; + while i < Self::MAX_LEN { + let mut b = Self::SHIFTS[i]; + while b < Self::SHIFTS[i + 1] { + table[b as usize] = i as u8; + b += 1; + } + i += 1; + } + table + }; /// The largest internal degree whose exponent sequences are guaranteed to fit. /// /// This is the largest bound for which the field widths sum to at most 64. It is far beyond /// anything reachable — the Milnor algebra already has over 5 million basis elements below /// degree 300 — and exceeds the degree 1536 the previous hand-rolled packing assumed. pub const MAX_DEGREE: i32 = 2045; - /// The number of entries that can be stored. This equals `fp`'s `MAX_MULTINOMIAL_LEN`, which /// already bounds the length of the $\xi$-degree table, so it is not a new restriction. pub const MAX_LEN: usize = 10; - - /// The length of the layout tables. This is [`Self::MAX_LEN`] rounded up to a power of two so - /// that [`Self::entry`] can mask its index instead of bounds-checking it; entries at or past - /// `MAX_LEN` are given width 0, so they read as zero. - const TABLE_LEN: usize = 16; - - /// `WIDTHS[i]` is the number of bits holding $r_{i+1}$: the number of bits needed to represent - /// `MAX_DEGREE / (2^(i+1) - 1)`. - const WIDTHS: [u32; Self::TABLE_LEN] = [11, 10, 9, 8, 7, 6, 5, 4, 3, 1, 0, 0, 0, 0, 0, 0]; - /// `SHIFTS[i]` is the bit offset of entry `i`; `SHIFTS[MAX_LEN]` is the total width, 64. const SHIFTS: [u32; Self::TABLE_LEN] = { let mut shifts = [0; Self::TABLE_LEN]; @@ -157,22 +161,13 @@ impl PPart { } shifts }; - - /// `FIELD_OF_BIT[b]` is the index of the entry owning bit `b`, letting [`Self::len`] turn a - /// `leading_zeros` into an entry index without looping. - const FIELD_OF_BIT: [u8; 64] = { - let mut table = [0; 64]; - let mut i = 0; - while i < Self::MAX_LEN { - let mut b = Self::SHIFTS[i]; - while b < Self::SHIFTS[i + 1] { - table[b as usize] = i as u8; - b += 1; - } - i += 1; - } - table - }; + /// The length of the layout tables. This is [`Self::MAX_LEN`] rounded up to a power of two so + /// that [`Self::entry`] can mask its index instead of bounds-checking it; entries at or past + /// `MAX_LEN` are given width 0, so they read as zero. + const TABLE_LEN: usize = 16; + /// `WIDTHS[i]` is the number of bits holding $r_{i+1}$: the number of bits needed to represent + /// `MAX_DEGREE / (2^(i+1) - 1)`. + const WIDTHS: [u32; Self::TABLE_LEN] = [11, 10, 9, 8, 7, 6, 5, 4, 3, 1, 0, 0, 0, 0, 0, 0]; /// The largest value entry `i` can hold. pub const fn max_entry(i: usize) -> PPartEntry {