la_stack/polynomial.rs
1#![forbid(unsafe_code)]
2
3//! Bounded exact determinants over sparse univariate rational polynomials.
4//!
5//! The subset recurrence reorganizes the Leibniz identity (REFERENCES.md \[12\]).
6//! See the [invariant and resource contract](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#exact-sparse-polynomial-determinants).
7
8use std::collections::{BTreeMap, btree_map::Entry};
9
10use num_rational::BigRational;
11use num_traits::{One, Zero};
12
13use crate::rational::canonicalize_rational;
14use crate::{LaError, NonFiniteLocation, NonFiniteOrigin, PolynomialWork};
15
16/// Largest supported order of [`polynomial_determinant`].
17pub const MAX_POLYNOMIAL_MATRIX_DIM: usize = 8;
18
19/// An exact polynomial `Σ coefficient[e] * t^e` with nonnegative `u128` exponents.
20///
21/// Terms are stored in increasing exponent order. Coefficients are reduced
22/// rational numbers with positive denominators; zero coefficients are absent.
23/// Duplicate input exponents are summed exactly, including cancellation.
24/// [`Self::leading_term`] means the *least* exponent with a nonzero coefficient,
25/// appropriate to an infinitesimal `t → 0+`, rather than the highest degree.
26///
27/// Storage and arbitrary-precision coefficients allocate on the heap. Construction
28/// consumes owned terms; queries borrow them; cloning makes a deep copy. The
29/// determinant budget does not cover construction. Multivariate encodings,
30/// perturbation-variable assignment and interpretation of coefficient signs
31/// belong to callers, who must prove their encoding preserves the desired order
32/// and exponent addition. No multivariate ordering is inferred here.
33#[derive(Clone, Debug, Default, Eq, PartialEq)]
34#[must_use]
35pub struct SparsePolynomial {
36 terms: BTreeMap<u128, BigRational>,
37}
38
39impl SparsePolynomial {
40 /// Construct the zero polynomial without allocating.
41 pub const fn zero() -> Self {
42 Self {
43 terms: BTreeMap::new(),
44 }
45 }
46
47 /// Consume terms, normalize rational coefficients and combine equal exponents.
48 ///
49 /// Input order is immaterial to a successful result. Raw unreduced or
50 /// negative-denominator rationals are accepted. Validation precedes zero
51 /// removal, so even `0/0` is rejected. The iterator is consumed only up to
52 /// the first invalid coefficient; this constructor has no work budget.
53 ///
54 /// # Errors
55 /// Returns [`LaError::NonFinite`] with input origin and a
56 /// [`NonFiniteLocation::PolynomialTerm`] for a zero denominator.
57 pub fn try_from_terms(
58 terms: impl IntoIterator<Item = (u128, BigRational)>,
59 ) -> Result<Self, LaError> {
60 let mut polynomial = Self::zero();
61 for (index, (exponent, coefficient)) in terms.into_iter().enumerate() {
62 if coefficient.denom().is_zero() {
63 return Err(LaError::NonFinite {
64 location: NonFiniteLocation::PolynomialTerm { index, exponent },
65 origin: NonFiniteOrigin::Input,
66 });
67 }
68 polynomial.add_term(exponent, canonicalize_rational(coefficient));
69 }
70 Ok(polynomial)
71 }
72
73 /// Borrow all nonzero terms in increasing exponent order.
74 #[must_use]
75 pub const fn terms(&self) -> &BTreeMap<u128, BigRational> {
76 &self.terms
77 }
78
79 /// Whether every coefficient is exactly zero.
80 #[must_use]
81 pub fn is_zero(&self) -> bool {
82 self.terms.is_empty()
83 }
84
85 /// Borrow the first nonzero coefficient in increasing exponent order.
86 ///
87 /// `None` proves this is the zero polynomial. A failed determinant evaluation
88 /// never returns a polynomial, so it cannot be confused with this result.
89 #[must_use]
90 pub fn leading_term(&self) -> Option<(u128, &BigRational)> {
91 self.terms
92 .first_key_value()
93 .map(|(&exponent, coefficient)| (exponent, coefficient))
94 }
95
96 fn add_term(&mut self, exponent: u128, coefficient: BigRational) {
97 if coefficient.is_zero() {
98 return;
99 }
100 match self.terms.entry(exponent) {
101 Entry::Vacant(entry) => {
102 entry.insert(coefficient);
103 }
104 Entry::Occupied(mut entry) => {
105 *entry.get_mut() += coefficient;
106 if entry.get().is_zero() {
107 entry.remove();
108 }
109 }
110 }
111 }
112}
113
114/// Caller-selected work limits for one [`polynomial_determinant`] invocation.
115///
116/// No default policy is supplied. All limits, including zero, are valid.
117/// These count combinatorial work, not elapsed time, bytes or integer bit
118/// operations: a single rational product can be expensive for wide inputs.
119#[derive(Clone, Copy, Debug, Eq, PartialEq)]
120#[must_use]
121pub struct PolynomialDeterminantBudget {
122 /// Maximum structural subset transitions. The complete `D * 2^(D-1)`
123 /// envelope is checked before allocation, even for an all-zero matrix.
124 /// Order zero requires zero transitions; order eight requires 1,024.
125 pub max_subset_transitions: usize,
126 /// Maximum individual nonzero coefficient-pair multiplications. Each is
127 /// charged before exponent addition or rational arithmetic; zero partials
128 /// and zero entries cost no products. Cancellation does not refund work.
129 pub max_term_products: usize,
130}
131
132/// Compute the complete exact determinant polynomial of a borrowed square array.
133///
134/// Supports orders `0..=8`; the empty determinant is the constant one. Singular
135/// polynomial matrices return the zero polynomial successfully. All cancellation
136/// terms are retained, so [`SparsePolynomial::leading_term`] is valid only after
137/// the complete calculation succeeds. There is no early sign fallback, truncation,
138/// division, rounding or conditioning-dependent numerical error.
139///
140/// The kernel allocates `2^D` sparse states (at most 256) and consumes each partial
141/// once all its predecessors have contributed. Inputs are neither cloned nor
142/// changed. At most one new stored term is created per charged product, plus the
143/// initial constant one. Rational sizes, allocator overhead, and caller-owned
144/// inputs are outside these combinatorial bounds; allocation failure follows
145/// Rust's allocator behavior. Limits do not promise a byte or wall-clock bound.
146///
147/// # Errors
148/// In order: [`LaError::UnsupportedDimension`] for `D > 8`,
149/// [`LaError::PolynomialBudgetExceeded`] for insufficient structural transitions,
150/// then product-budget or [`LaError::PolynomialExponentOverflow`] failures during
151/// evaluation. An overflowing intermediate exponent is an error even if later
152/// contributions could cancel it. Errors discard partial results and leave the
153/// input reusable. Traversal is increasing subset mask, column, left exponent,
154/// then right exponent, making budget accounting deterministic.
155///
156/// # Examples
157/// ```
158/// use la_stack::prelude::*;
159/// # fn main() -> Result<(), LaError> {
160/// let one = SparsePolynomial::try_from_terms([(0, BigRational::from_integer(1.into()))])?;
161/// let perturbed = SparsePolynomial::try_from_terms([
162/// (0, BigRational::from_integer(1.into())),
163/// (3, BigRational::new(1.into(), 2.into())),
164/// ])?;
165/// let matrix = [[one.clone(), one.clone()], [one, perturbed]];
166/// let determinant = polynomial_determinant(&matrix, PolynomialDeterminantBudget {
167/// max_subset_transitions: 4,
168/// max_term_products: 5,
169/// })?;
170/// assert_eq!(determinant.leading_term(), Some((3, &BigRational::new(1.into(), 2.into()))));
171/// # Ok(())
172/// # }
173/// ```
174pub fn polynomial_determinant<const D: usize>(
175 matrix: &[[SparsePolynomial; D]; D],
176 budget: PolynomialDeterminantBudget,
177) -> Result<SparsePolynomial, LaError> {
178 if D > MAX_POLYNOMIAL_MATRIX_DIM {
179 return Err(LaError::unsupported_dimension(D, MAX_POLYNOMIAL_MATRIX_DIM));
180 }
181 let state_count = 1usize << D;
182 let transitions = D * (state_count / 2);
183 if transitions > budget.max_subset_transitions {
184 return Err(LaError::PolynomialBudgetExceeded {
185 work: PolynomialWork::SubsetTransitions,
186 limit: budget.max_subset_transitions,
187 completed: 0,
188 requested: transitions,
189 });
190 }
191 let mut partials = vec![SparsePolynomial::zero(); state_count];
192 partials[0].add_term(0, BigRational::one());
193 let mut products = 0;
194 for mask in 0..state_count - 1 {
195 let row = mask.count_ones() as usize;
196 let partial = std::mem::take(&mut partials[mask]);
197 if partial.is_zero() {
198 continue;
199 }
200 for (col, entry) in matrix[row].iter().enumerate() {
201 if mask & (1 << col) != 0 {
202 continue;
203 }
204 let negative = (mask >> (col + 1)).count_ones() % 2 != 0;
205 let next = &mut partials[mask | (1 << col)];
206 for (&left, left_coefficient) in &partial.terms {
207 for (&right, right_coefficient) in &entry.terms {
208 if products == budget.max_term_products {
209 return Err(LaError::PolynomialBudgetExceeded {
210 work: PolynomialWork::TermProducts,
211 limit: budget.max_term_products,
212 completed: products,
213 requested: 1,
214 });
215 }
216 products += 1;
217 let exponent =
218 left.checked_add(right)
219 .ok_or(LaError::PolynomialExponentOverflow {
220 row,
221 col,
222 left,
223 right,
224 })?;
225 let coefficient = left_coefficient * right_coefficient;
226 next.add_term(exponent, if negative { -coefficient } else { coefficient });
227 }
228 }
229 }
230 }
231 Ok(std::mem::take(&mut partials[state_count - 1]))
232}