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.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/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_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs index 07b3bd35f8..cd72a42a56 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,230 @@ 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 `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 +/// 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 { + /// `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; + /// `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]; + let mut i = 0; + while i < Self::MAX_LEN { + shifts[i + 1] = shifts[i] + Self::WIDTHS[i]; + i += 1; + } + shifts + }; + /// 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 { + ((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::packed_signature` in `ext`). + pub const fn bits(self) -> u64 { + self.0 + } + + /// 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.entry(i) + } + + /// 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 +350,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 +359,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 +382,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 +412,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, @@ -257,7 +431,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>, @@ -323,8 +504,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 { @@ -371,7 +565,7 @@ impl Algebra for MilnorAlgebra { MilnorBasisElement { degree: 1, q_part: 1, - p_part: vec![], + p_part: PPart::zero(), }, )); } @@ -383,7 +577,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 +596,7 @@ impl Algebra for MilnorAlgebra { MilnorBasisElement { degree, q_part: 0, - p_part: vec![1 << i], + p_part: PPart::from_iter([1 << i]), }, )); } @@ -417,6 +611,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() { @@ -428,11 +629,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() { - map.insert(b.clone(), i); + 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 }); @@ -453,8 +655,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 }) @@ -474,7 +676,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"))] @@ -585,19 +791,34 @@ 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 >= combinatorics::xi_degrees(p).len() { + return None; + } + 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(); + 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 +829,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) @@ -637,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) } } @@ -715,7 +940,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 +959,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 +978,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 +1058,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 +1088,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 +1142,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, }), @@ -933,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.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)); - } + table.sort_by_cached_key(|e| e.excess(fp::prime::TWO)); table }); } @@ -971,13 +1197,21 @@ 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 { + 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)) } @@ -987,8 +1221,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.clone())]; + 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) { @@ -1011,23 +1245,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, @@ -1051,8 +1283,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) @@ -1063,8 +1295,8 @@ impl MilnorAlgebra { &self, mut res: FpSliceMut, coef: u32, - m1: &MilnorBasisElement, - m2: &MilnorBasisElement, + m1: MilnorBasisElement, + m2: MilnorBasisElement, excess: i32, mut allocation: PPartAllocation, ) -> PPartAllocation { @@ -1074,8 +1306,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 +1324,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, @@ -1114,7 +1346,7 @@ impl MilnorAlgebra { &self, res: FpSliceMut, coef: u32, - m1: &MilnorBasisElement, + m1: MilnorBasisElement, s_deg: i32, s: FpSlice, ) { @@ -1127,7 +1359,7 @@ impl MilnorAlgebra { &self, mut res: FpSliceMut, coef: u32, - m1: &MilnorBasisElement, + m1: MilnorBasisElement, s_deg: i32, s: FpSlice, mut allocation: PPartAllocation, @@ -1149,7 +1381,7 @@ impl MilnorAlgebra { #[derive(Debug, Default)] struct Matrix2D { cols: usize, - inner: PPart, + inner: Vec, } impl std::fmt::Display for Matrix2D { @@ -1202,8 +1434,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 +1449,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 +1460,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 +1482,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 +1504,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.entry(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.entry(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 +1534,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 +1617,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.entry(row - 1); for col in 1..self.cols { self.M[0][col] += self.M[row][col]; self.M[row][col] = 0; @@ -1416,13 +1641,12 @@ 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(); let mut coef = 1; if self.init { @@ -1443,27 +1667,22 @@ 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. 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.entry(i); + if i + 1 < self.cols { + entry += self.M[0][i + 1]; + } + 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() { - self.ans.p_part.reserve(self.diag_num); + 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); @@ -1510,12 +1729,17 @@ 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 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 { @@ -1543,7 +1767,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 +1775,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 +1791,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 +1831,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 +1849,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 +1878,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 +1895,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 +1928,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 +1957,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 +1965,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, @@ -1906,6 +2124,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); @@ -1920,13 +2167,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 +2252,216 @@ 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); + } + } + } + + /// 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. + #[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(); @@ -2012,34 +2473,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 +2536,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/milnor_rank.rs b/ext/crates/algebra/src/algebra/milnor_rank.rs new file mode 100644 index 0000000000..65c6f0a7e4 --- /dev/null +++ b/ext/crates/algebra/src/algebra/milnor_rank.rs @@ -0,0 +1,267 @@ +//! 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. +//! +//! 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; + +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..6627432dcd 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}; diff --git a/ext/crates/algebra/src/algebra/pair_algebra.rs b/ext/crates/algebra/src/algebra/pair_algebra.rs index e080d34116..34815baaa0 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() { @@ -204,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, ); @@ -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,16 +263,17 @@ 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(), coeff * c, - &r, + r, 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 @@ -364,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), @@ -377,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 { @@ -385,19 +387,19 @@ 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); } } }) } /// 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(); +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, - p_part: vec![], + p_part: PPart::zero(), degree: 0, }; @@ -410,17 +412,15 @@ 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; // 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); } @@ -445,7 +445,7 @@ mod tests { MilnorElt { q_part: 0, - p_part: p_part.into(), + p_part: PPart::from_slice(p_part), degree, } } @@ -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 6d42a1f9f8..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 { @@ -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..9c847781ea 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,40 @@ 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. + /// 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); } - true + Some((mask, value)) } fn zero_signature(&self) -> Vec { @@ -125,25 +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: _, - }| { - algebra - .ppart_table(degree - gen_deg) - .iter() - .enumerate() - .filter_map(move |(n, op)| { - if self.has_signature(op, signature) { - 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 @@ -1325,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:?}" + ); + } + } + } + } + } } 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 }