la_stack/lib.rs
1#![forbid(unsafe_code)]
2#![warn(missing_docs)]
3#![doc = include_str!("../README.md")]
4
5#[cfg(doc)]
6mod readme_doctests {
7 //! Executable versions of README examples.
8 /// ```rust
9 /// use la_stack::prelude::*;
10 ///
11 /// # fn main() -> Result<(), LaError> {
12 /// // This system requires pivoting (a[0][0] = 0), so it's a good LU demo.
13 /// let a = Matrix::<5>::try_from_rows([
14 /// [0.0, 1.0, 1.0, 1.0, 1.0],
15 /// [1.0, 0.0, 1.0, 1.0, 1.0],
16 /// [1.0, 1.0, 0.0, 1.0, 1.0],
17 /// [1.0, 1.0, 1.0, 0.0, 1.0],
18 /// [1.0, 1.0, 1.0, 1.0, 0.0],
19 /// ])?;
20 ///
21 /// let b = Vector::<5>::try_new([14.0, 13.0, 12.0, 11.0, 10.0])?;
22 ///
23 /// let lu = a.lu(DEFAULT_PIVOT_TOL)?;
24 /// let x = lu.solve_vec(b)?.into_array();
25 ///
26 /// // Floating-point rounding is expected; compare with a tolerance.
27 /// let expected = [1.0, 2.0, 3.0, 4.0, 5.0];
28 /// for (x_i, e_i) in x.iter().zip(expected.iter()) {
29 /// assert!((*x_i - *e_i).abs() <= 1e-12);
30 /// }
31 /// # Ok(())
32 /// # }
33 /// ```
34 fn solve_5x5_example() {}
35
36 /// ```rust
37 /// use la_stack::prelude::*;
38 ///
39 /// # fn main() -> Result<(), LaError> {
40 /// // This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting.
41 /// let a = Matrix::<5>::try_from_rows([
42 /// [1.0, 1.0, 0.0, 0.0, 0.0],
43 /// [1.0, 2.0, 1.0, 0.0, 0.0],
44 /// [0.0, 1.0, 2.0, 1.0, 0.0],
45 /// [0.0, 0.0, 1.0, 2.0, 1.0],
46 /// [0.0, 0.0, 0.0, 1.0, 2.0],
47 /// ])?;
48 ///
49 /// let det = a.ldlt(DEFAULT_SINGULAR_TOL)?.det()?;
50 /// assert!((det - 1.0).abs() <= 1e-12);
51 /// # Ok(())
52 /// # }
53 /// ```
54 fn det_5x5_ldlt_example() {}
55}
56
57#[cfg(feature = "exact")]
58mod exact;
59#[cfg(feature = "exact")]
60pub use num_bigint::BigInt;
61#[cfg(feature = "exact")]
62pub use num_rational::BigRational;
63#[cfg(feature = "exact")]
64pub use num_traits::{FromPrimitive, Signed, ToPrimitive};
65
66mod ldlt;
67mod lu;
68mod matrix;
69mod vector;
70
71use core::fmt;
72
73// ---------------------------------------------------------------------------
74// Error-bound constants for `Matrix::det_errbound()`.
75//
76// For `D ∈ {2, 3, 4}`, `Matrix::det_direct()` evaluates the Leibniz expansion
77// of the determinant as a tree of f64 multiplies and fused multiply-adds
78// (FMAs). Following Shewchuk's error-analysis methodology (REFERENCES.md
79// [8]), the absolute error of that computation is bounded by
80//
81// |det_direct(A) - det_exact(A)| ≤ ERR_COEFF_D · p(|A|)
82//
83// where `p(|A|)` is the **absolute Leibniz sum**
84//
85// p(|A|) = Σ_σ ∏ᵢ |A[i, σ(i)]|,
86//
87// i.e. the same cofactor-expansion tree as `det_direct` but with each
88// entry replaced by its magnitude. Note that `p(|A|)` is *not* the
89// combinatorial matrix permanent — the name "permanent" appears in the
90// source for brevity and to match the cited literature.
91//
92// Each constant has the shape `a · EPS + b · EPS²`: the linear term bounds
93// the first-order rounding and the quadratic term absorbs the interaction
94// of errors in nested FMAs. The coefficients `a` and `b` are conservative
95// over-estimates derived from the longest dependency chain of `det_direct`
96// at that dimension.
97//
98// These constants are NOT feature-gated — they rely only on f64 arithmetic
99// and are useful for adaptive-precision logic even without the `exact`
100// feature. Most callers should prefer `Matrix::det_errbound()`, which
101// applies these constants to the actual matrix; the raw constants are
102// exposed for advanced use cases (composing the bound with a pre-reduced
103// permanent, rolling a custom adaptive filter, etc.). See
104// `Matrix::det_sign_exact()` (behind the `exact` feature) for the
105// reference adaptive-filter that consumes these internally.
106// ---------------------------------------------------------------------------
107
108const EPS: f64 = f64::EPSILON; // 2^-52
109
110/// Absolute error coefficient for [`Matrix::<2>::det_direct`](crate::Matrix::det_direct).
111///
112/// This constant is not a caller-tuned tolerance. It is the dimension-specific
113/// multiplier that turns the matrix's absolute Leibniz sum into a conservative
114/// bound on floating-point roundoff in the closed-form 2×2 determinant formula.
115///
116/// For any 2×2 matrix `A = [[a, b], [c, d]]` with finite f64 entries,
117///
118/// ```text
119/// |A.det_direct() - det_exact(A)| ≤ ERR_COEFF_2 · (|a·d| + |b·c|)
120/// ```
121///
122/// `det_direct` evaluates `a·d - b·c` as one multiply followed by one FMA
123/// (2 rounding events); the linear `3·EPS` term bounds those roundings
124/// and the quadratic `16·EPS²` term is a conservative cushion for their
125/// interaction. Derivation follows Shewchuk's framework; see
126/// `REFERENCES.md` \[8\].
127///
128/// Prefer [`Matrix::det_errbound`](crate::Matrix::det_errbound) unless
129/// you already have the absolute-Leibniz sum available; see
130/// `Matrix::det_sign_exact` (under the `exact` feature) for the reference
131/// adaptive-precision filter.
132///
133/// # Example
134/// ```
135/// use la_stack::prelude::*;
136///
137/// # fn main() -> Result<(), LaError> {
138/// let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
139/// let Some(det) = m.det_direct()? else {
140/// return Ok(());
141/// };
142/// assert_eq!(det, -2.0);
143/// // Compute the bound from the raw constant for illustration; most
144/// // callers would match on `m.det_errbound()?` instead.
145/// let p = (1.0_f64 * 4.0).abs() + (2.0_f64 * 3.0).abs();
146/// let bound = ERR_COEFF_2 * p;
147/// if det.abs() > bound {
148/// // The f64 sign is provably correct without exact arithmetic.
149/// }
150/// # Ok(())
151/// # }
152/// ```
153pub const ERR_COEFF_2: f64 = 3.0 * EPS + 16.0 * EPS * EPS;
154
155/// Absolute error coefficient for [`Matrix::<3>::det_direct`](crate::Matrix::det_direct).
156///
157/// This constant is not a caller-tuned tolerance. It is the dimension-specific
158/// multiplier that turns the matrix's absolute Leibniz sum into a conservative
159/// bound on floating-point roundoff in the closed-form 3×3 determinant formula.
160///
161/// For any 3×3 matrix `A` with finite f64 entries,
162///
163/// ```text
164/// |A.det_direct() - det_exact(A)| ≤ ERR_COEFF_3 · p(|A|)
165/// ```
166///
167/// where `p(|A|)` is the absolute Leibniz sum (the same cofactor
168/// expansion as `det_direct` but with `|·|` at every leaf).
169/// `det_direct` for D=3 uses three 2×2 FMA minors combined by a nested
170/// FMA, yielding the `8·EPS + 64·EPS²` bound. See `REFERENCES.md`
171/// \[8\] for the Shewchuk framework these bounds follow.
172///
173/// Prefer [`Matrix::det_errbound`](crate::Matrix::det_errbound) over this
174/// constant for typical use; see [`ERR_COEFF_2`] for a worked example.
175pub const ERR_COEFF_3: f64 = 8.0 * EPS + 64.0 * EPS * EPS;
176
177/// Absolute error coefficient for [`Matrix::<4>::det_direct`](crate::Matrix::det_direct).
178///
179/// This constant is not a caller-tuned tolerance. It is the dimension-specific
180/// multiplier that turns the matrix's absolute Leibniz sum into a conservative
181/// bound on floating-point roundoff in the closed-form 4×4 determinant formula.
182///
183/// For any 4×4 matrix `A` with finite f64 entries,
184///
185/// ```text
186/// |A.det_direct() - det_exact(A)| ≤ ERR_COEFF_4 · p(|A|)
187/// ```
188///
189/// where `p(|A|)` is the absolute Leibniz sum. `det_direct` for D=4
190/// hoists six 2×2 minors, combines them into four 3×3 cofactors, then
191/// reduces those with an FMA row combination, yielding the
192/// `12·EPS + 128·EPS²` bound. See `REFERENCES.md` \[8\] for the
193/// Shewchuk framework these bounds follow.
194///
195/// Prefer [`Matrix::det_errbound`](crate::Matrix::det_errbound) over this
196/// constant for typical use; see [`ERR_COEFF_2`] for a worked example.
197pub const ERR_COEFF_4: f64 = 12.0 * EPS + 128.0 * EPS * EPS;
198
199/// Largest dimension supported by [`try_with_stack_matrix!`].
200///
201/// The crate can represent `Matrix<D>` for any compile-time `D`, but runtime
202/// dispatch must enumerate a finite set of concrete stack types. Dimensions
203/// `0..=7` cover downstream geometric predicate matrices while keeping the
204/// dispatch surface explicit.
205pub const MAX_STACK_MATRIX_DISPATCH_DIM: usize = 7;
206
207/// Finite, non-negative tolerance used by numerical predicates and factorizations.
208///
209/// Construct with [`Tolerance::new`] when accepting raw caller input. Once
210/// constructed, the stored value is guaranteed to be finite and `>= 0`, so
211/// downstream algorithms do not need to revalidate the tolerance.
212///
213/// This is the crate-wide tolerance contract: raw negative, NaN, and infinite
214/// values are rejected with [`LaError::InvalidTolerance`] at construction time.
215#[must_use]
216#[derive(Clone, Copy, Debug, PartialEq)]
217pub struct Tolerance {
218 value: f64,
219}
220
221impl Tolerance {
222 /// Construct a tolerance without checking the raw value.
223 ///
224 /// This crate-internal escape hatch is only for constants whose finite,
225 /// non-negative value is visible at the call site. Public callers should
226 /// use [`Tolerance::new`] so the returned value carries the validation
227 /// proof.
228 pub(crate) const fn new_unchecked(value: f64) -> Self {
229 Self { value }
230 }
231
232 /// Construct a finite, non-negative tolerance.
233 ///
234 /// # Examples
235 /// ```
236 /// use la_stack::prelude::*;
237 ///
238 /// # fn main() -> Result<(), LaError> {
239 /// let tol = Tolerance::new(1e-12)?;
240 /// assert_eq!(tol.get(), 1e-12);
241 /// # Ok(())
242 /// # }
243 /// ```
244 ///
245 /// # Errors
246 /// Returns [`LaError::InvalidTolerance`] when `value` is NaN, infinite, or
247 /// negative.
248 #[inline]
249 pub const fn new(value: f64) -> Result<Self, LaError> {
250 if value >= 0.0 && value.is_finite() {
251 Ok(Self::new_unchecked(value))
252 } else {
253 Err(LaError::invalid_tolerance(value))
254 }
255 }
256
257 /// Return the raw finite, non-negative tolerance value.
258 ///
259 /// # Examples
260 /// ```
261 /// use la_stack::prelude::*;
262 ///
263 /// # fn main() -> Result<(), LaError> {
264 /// let tol = Tolerance::new(0.0)?;
265 /// assert_eq!(tol.get(), 0.0);
266 /// # Ok(())
267 /// # }
268 /// ```
269 #[inline]
270 #[must_use]
271 pub const fn get(self) -> f64 {
272 self.value
273 }
274}
275
276/// Default absolute threshold used for singularity/degeneracy detection.
277///
278/// This is intentionally conservative for geometric predicates and small systems.
279///
280/// Conceptually, this is an absolute bound for deciding when a scalar should be treated
281/// as "numerically zero" (e.g. LU pivots, LDLT diagonal entries).
282pub const DEFAULT_SINGULAR_TOL: Tolerance = Tolerance::new_unchecked(1e-12);
283
284/// Default absolute pivot magnitude threshold used for LU pivot selection / singularity detection.
285///
286/// This name is kept for backwards compatibility; prefer [`DEFAULT_SINGULAR_TOL`] when the
287/// tolerance is not specifically about pivot selection.
288pub const DEFAULT_PIVOT_TOL: Tolerance = DEFAULT_SINGULAR_TOL;
289
290/// Relative tolerance used to validate matrices for LDLT factorization.
291pub(crate) const LDLT_SYMMETRY_REL_TOL: Tolerance = Tolerance::new_unchecked(1e-12);
292
293/// Linear algebra errors.
294///
295/// This enum is `#[non_exhaustive]` — downstream `match` arms must include a
296/// wildcard (`_`) pattern to compile, allowing new variants to be added in
297/// future minor releases without breaking existing code.
298#[derive(Clone, Copy, Debug, PartialEq)]
299#[non_exhaustive]
300pub enum LaError {
301 /// The matrix is (numerically) singular.
302 Singular {
303 /// The factorization column/step where a suitable pivot/diagonal could not be found.
304 pivot_col: usize,
305 },
306 /// A non-finite value (NaN/∞) was encountered.
307 ///
308 /// The `(row, col)` coordinate follows a consistent convention across the crate:
309 ///
310 /// - `row: Some(r), col: c` — the non-finite value is tied to a matrix/factor
311 /// cell at `(r, c)`, either because a stored input/factor cell is already
312 /// non-finite or because factorization computed a non-finite value for
313 /// that cell before storing it.
314 /// - `row: None, col: c` — the non-finite value is tied to a vector entry,
315 /// determinant product, tolerance-scale accumulator, solve accumulator, or
316 /// other scalar/intermediate that has no matrix row coordinate.
317 NonFinite {
318 /// Row of the non-finite entry for a stored matrix cell, or `None` for
319 /// a vector-input entry or a computed intermediate. See the variant
320 /// docs for the full convention.
321 row: Option<usize>,
322 /// Column index (stored cell), vector index, or factorization/solve
323 /// step where the non-finite value was detected.
324 col: usize,
325 },
326 /// The exact result overflows the target representation (e.g. `f64`).
327 ///
328 /// Returned by `Matrix::det_exact_f64` and `Matrix::solve_exact_f64`
329 /// (requires `exact` feature) when an exact value is too large to
330 /// represent as a finite `f64`.
331 Overflow {
332 /// For vector results (e.g. `solve_exact_f64`), the index of the
333 /// component that overflowed. `None` for scalar results.
334 index: Option<usize>,
335 },
336 /// A requested runtime matrix dimension has no stack-dispatch arm.
337 UnsupportedDimension {
338 /// Runtime dimension requested by the caller.
339 requested: usize,
340 /// Largest runtime dimension supported by the dispatch helper.
341 max: usize,
342 },
343 /// A matrix index is outside the `D×D` storage domain.
344 IndexOutOfBounds {
345 /// Requested row index.
346 row: usize,
347 /// Requested column index.
348 col: usize,
349 /// Matrix dimension `D`; valid row and column indices are `< dim`.
350 dim: usize,
351 },
352 /// A tolerance value is not finite and non-negative.
353 InvalidTolerance {
354 /// Raw tolerance supplied by the caller.
355 value: f64,
356 },
357 /// A matrix required to be symmetric has an asymmetric off-diagonal pair.
358 Asymmetric {
359 /// Row index of the first asymmetric pair.
360 row: usize,
361 /// Column index of the first asymmetric pair.
362 col: usize,
363 /// Matrix dimension `D`.
364 dim: usize,
365 },
366 /// A symmetric matrix failed the positive-semidefinite LDLT domain check.
367 NotPositiveSemidefinite {
368 /// Factorization column/step where a negative LDLT diagonal was found.
369 pivot_col: usize,
370 /// Negative diagonal value observed at that step.
371 value: f64,
372 },
373}
374
375impl LaError {
376 /// Construct a [`LaError::NonFinite`] pinpointing a stored matrix cell at `(row, col)`.
377 ///
378 /// Use this for non-finite values read from a stored `Matrix` entry or
379 /// factorization cell, and for non-finite factorization updates that would
380 /// be stored at `(row, col)` if accepted. The resulting error has
381 /// `row: Some(row), col`, matching the matrix/factor-cell convention
382 /// documented on [`NonFinite`](Self::NonFinite). For vector-input entries
383 /// or scalar intermediates without a matrix row coordinate, use
384 /// [`non_finite_at`](Self::non_finite_at).
385 ///
386 /// # Examples
387 /// ```
388 /// use la_stack::prelude::*;
389 ///
390 /// assert_eq!(
391 /// LaError::non_finite_cell(1, 2),
392 /// LaError::NonFinite {
393 /// row: Some(1),
394 /// col: 2,
395 /// }
396 /// );
397 /// ```
398 #[inline]
399 #[must_use]
400 pub const fn non_finite_cell(row: usize, col: usize) -> Self {
401 Self::NonFinite {
402 row: Some(row),
403 col,
404 }
405 }
406
407 /// Construct a [`LaError::NonFinite`] pinpointing a vector-input entry or
408 /// computed scalar/intermediate at index `col`.
409 ///
410 /// Use this for non-finite values in a `Vector` input, determinant scalar,
411 /// tolerance-scale accumulator, or solve accumulator that overflowed during
412 /// forward/back substitution. The resulting error has `row: None, col`,
413 /// matching the vector/scalar-intermediate convention documented on
414 /// [`NonFinite`](Self::NonFinite). For stored matrix cells or computed
415 /// factorization updates tied to a matrix cell, use
416 /// [`non_finite_cell`](Self::non_finite_cell).
417 ///
418 /// # Examples
419 /// ```
420 /// use la_stack::prelude::*;
421 ///
422 /// assert_eq!(
423 /// LaError::non_finite_at(2),
424 /// LaError::NonFinite { row: None, col: 2 }
425 /// );
426 /// ```
427 #[inline]
428 #[must_use]
429 pub const fn non_finite_at(col: usize) -> Self {
430 Self::NonFinite { row: None, col }
431 }
432
433 /// Construct a [`LaError::UnsupportedDimension`] for runtime stack dispatch.
434 ///
435 /// # Examples
436 /// ```
437 /// use la_stack::prelude::*;
438 ///
439 /// assert_eq!(
440 /// LaError::unsupported_dimension(8, MAX_STACK_MATRIX_DISPATCH_DIM),
441 /// LaError::UnsupportedDimension {
442 /// requested: 8,
443 /// max: MAX_STACK_MATRIX_DISPATCH_DIM,
444 /// }
445 /// );
446 /// ```
447 #[inline]
448 #[must_use]
449 pub const fn unsupported_dimension(requested: usize, max: usize) -> Self {
450 Self::UnsupportedDimension { requested, max }
451 }
452
453 /// Construct a [`LaError::IndexOutOfBounds`] for a `D×D` matrix index.
454 ///
455 /// # Examples
456 /// ```
457 /// use la_stack::prelude::*;
458 ///
459 /// assert_eq!(
460 /// LaError::index_out_of_bounds(2, 0, 2),
461 /// LaError::IndexOutOfBounds {
462 /// row: 2,
463 /// col: 0,
464 /// dim: 2,
465 /// }
466 /// );
467 /// ```
468 #[inline]
469 #[must_use]
470 pub const fn index_out_of_bounds(row: usize, col: usize, dim: usize) -> Self {
471 Self::IndexOutOfBounds { row, col, dim }
472 }
473
474 /// Construct a [`LaError::InvalidTolerance`] for a raw tolerance value.
475 ///
476 /// # Examples
477 /// ```
478 /// use la_stack::prelude::*;
479 ///
480 /// assert_eq!(
481 /// LaError::invalid_tolerance(-1.0),
482 /// LaError::InvalidTolerance { value: -1.0 }
483 /// );
484 /// ```
485 #[inline]
486 #[must_use]
487 pub const fn invalid_tolerance(value: f64) -> Self {
488 Self::InvalidTolerance { value }
489 }
490
491 /// Construct a [`LaError::Asymmetric`] for a `D×D` matrix.
492 ///
493 /// # Examples
494 /// ```
495 /// use la_stack::prelude::*;
496 ///
497 /// assert_eq!(
498 /// LaError::asymmetric(0, 1, 3),
499 /// LaError::Asymmetric {
500 /// row: 0,
501 /// col: 1,
502 /// dim: 3,
503 /// }
504 /// );
505 /// ```
506 #[inline]
507 #[must_use]
508 pub const fn asymmetric(row: usize, col: usize, dim: usize) -> Self {
509 Self::Asymmetric { row, col, dim }
510 }
511
512 /// Construct a [`LaError::NotPositiveSemidefinite`] for LDLT factorization.
513 ///
514 /// # Examples
515 /// ```
516 /// use la_stack::prelude::*;
517 ///
518 /// assert_eq!(
519 /// LaError::not_positive_semidefinite(1, -3.0),
520 /// LaError::NotPositiveSemidefinite {
521 /// pivot_col: 1,
522 /// value: -3.0,
523 /// }
524 /// );
525 /// ```
526 #[inline]
527 #[must_use]
528 pub const fn not_positive_semidefinite(pivot_col: usize, value: f64) -> Self {
529 Self::NotPositiveSemidefinite { pivot_col, value }
530 }
531
532 /// Parse a raw tolerance into a finite, non-negative [`Tolerance`].
533 ///
534 /// # Examples
535 /// ```
536 /// use la_stack::prelude::*;
537 ///
538 /// assert_eq!(LaError::validate_tolerance(1e-12)?.get(), 1e-12);
539 ///
540 /// let raw = 0.0;
541 /// let tol = LaError::validate_tolerance(raw)?;
542 /// let _lu = Matrix::<2>::identity().lu(tol)?;
543 ///
544 /// assert_eq!(
545 /// LaError::validate_tolerance(-1.0),
546 /// Err(LaError::InvalidTolerance { value: -1.0 })
547 /// );
548 /// # Ok::<(), LaError>(())
549 /// ```
550 ///
551 /// # Errors
552 /// Returns [`LaError::InvalidTolerance`] when `value` is NaN, infinite, or
553 /// negative.
554 #[inline]
555 pub const fn validate_tolerance(value: f64) -> Result<Tolerance, Self> {
556 Tolerance::new(value)
557 }
558}
559
560impl fmt::Display for LaError {
561 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
562 match *self {
563 Self::Singular { pivot_col } => {
564 write!(f, "singular matrix at pivot column {pivot_col}")
565 }
566 Self::NonFinite { row: Some(r), col } => {
567 write!(f, "non-finite value at ({r}, {col})")
568 }
569 Self::NonFinite { row: None, col } => {
570 write!(f, "non-finite value at index {col}")
571 }
572 Self::Overflow { index: Some(i) } => {
573 write!(
574 f,
575 "exact result overflows the target representation at index {i}"
576 )
577 }
578 Self::Overflow { index: None } => {
579 write!(f, "exact result overflows the target representation")
580 }
581 Self::UnsupportedDimension { requested, max } => {
582 write!(
583 f,
584 "unsupported matrix dimension {requested}; maximum stack-dispatch dimension is {max}"
585 )
586 }
587 Self::IndexOutOfBounds { row, col, dim } => {
588 write!(
589 f,
590 "matrix index ({row}, {col}) is out of bounds for dimension {dim}"
591 )
592 }
593 Self::InvalidTolerance { value } => {
594 write!(f, "invalid tolerance {value}; expected finite value >= 0")
595 }
596 Self::Asymmetric { row, col, dim } => {
597 write!(
598 f,
599 "matrix is not symmetric for dimension {dim}: asymmetric pair ({row}, {col})"
600 )
601 }
602 Self::NotPositiveSemidefinite { pivot_col, value } => {
603 write!(
604 f,
605 "matrix is not positive semidefinite at LDLT pivot column {pivot_col}: diagonal value {value} < 0"
606 )
607 }
608 }
609 }
610}
611
612impl std::error::Error for LaError {}
613
614pub use ldlt::Ldlt;
615pub use lu::Lu;
616pub use matrix::Matrix;
617pub use vector::Vector;
618
619/// Fallibly dispatch a runtime dimension to a concrete stack-allocated matrix.
620///
621/// The macro creates a zero matrix with type `Matrix<N>` for the selected
622/// runtime dimension `N`, then evaluates the supplied closure body. Supported
623/// runtime dimensions run from `0` through [`MAX_STACK_MATRIX_DISPATCH_DIM`].
624/// Unsupported dimensions return
625/// `Err(LaError::UnsupportedDimension { requested, max })` converted with
626/// `From<LaError>`, so downstream crates can use their own public error type.
627///
628/// # Errors
629/// Returns [`LaError::UnsupportedDimension`] (converted through `From<LaError>`)
630/// when the requested runtime dimension is greater than
631/// [`MAX_STACK_MATRIX_DISPATCH_DIM`]. The closure body may return any other
632/// error representable by its declared `Result` type.
633///
634/// # Examples
635/// ```
636/// use la_stack::prelude::*;
637///
638/// # fn main() -> Result<(), LaError> {
639/// let requested = 2usize;
640/// let det = try_with_stack_matrix!(requested, |mut m| -> Result<f64, LaError> {
641/// m.set_checked(0, 0, 1.0)?;
642/// m.set_checked(1, 1, 1.0)?;
643/// m.det()
644/// })?;
645///
646/// assert_eq!(det, 1.0);
647/// # Ok(())
648/// # }
649/// ```
650#[macro_export]
651macro_rules! try_with_stack_matrix {
652 ($dim:expr, |$matrix:ident| -> $ret:ty $body:block $(,)?) => {{
653 let __la_stack_requested_dim: usize = $dim;
654 match __la_stack_requested_dim {
655 0 => $crate::try_with_stack_matrix!(@arm 0, $matrix, $ret, $body),
656 1 => $crate::try_with_stack_matrix!(@arm 1, $matrix, $ret, $body),
657 2 => $crate::try_with_stack_matrix!(@arm 2, $matrix, $ret, $body),
658 3 => $crate::try_with_stack_matrix!(@arm 3, $matrix, $ret, $body),
659 4 => $crate::try_with_stack_matrix!(@arm 4, $matrix, $ret, $body),
660 5 => $crate::try_with_stack_matrix!(@arm 5, $matrix, $ret, $body),
661 6 => $crate::try_with_stack_matrix!(@arm 6, $matrix, $ret, $body),
662 7 => $crate::try_with_stack_matrix!(@arm 7, $matrix, $ret, $body),
663 requested => Err(::core::convert::From::from(
664 $crate::LaError::unsupported_dimension(
665 requested,
666 $crate::MAX_STACK_MATRIX_DISPATCH_DIM,
667 ),
668 )),
669 }
670 }};
671 ($dim:expr, |mut $matrix:ident| -> $ret:ty $body:block $(,)?) => {{
672 let __la_stack_requested_dim: usize = $dim;
673 match __la_stack_requested_dim {
674 0 => $crate::try_with_stack_matrix!(@arm_mut 0, $matrix, $ret, $body),
675 1 => $crate::try_with_stack_matrix!(@arm_mut 1, $matrix, $ret, $body),
676 2 => $crate::try_with_stack_matrix!(@arm_mut 2, $matrix, $ret, $body),
677 3 => $crate::try_with_stack_matrix!(@arm_mut 3, $matrix, $ret, $body),
678 4 => $crate::try_with_stack_matrix!(@arm_mut 4, $matrix, $ret, $body),
679 5 => $crate::try_with_stack_matrix!(@arm_mut 5, $matrix, $ret, $body),
680 6 => $crate::try_with_stack_matrix!(@arm_mut 6, $matrix, $ret, $body),
681 7 => $crate::try_with_stack_matrix!(@arm_mut 7, $matrix, $ret, $body),
682 requested => Err(::core::convert::From::from(
683 $crate::LaError::unsupported_dimension(
684 requested,
685 $crate::MAX_STACK_MATRIX_DISPATCH_DIM,
686 ),
687 )),
688 }
689 }};
690 (@arm $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
691 let __la_stack_body = |$matrix: $crate::Matrix<$d>| -> $ret { $body };
692 __la_stack_body($crate::Matrix::<$d>::zero())
693 }};
694 (@arm_mut $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
695 let __la_stack_body = |mut $matrix: $crate::Matrix<$d>| -> $ret { $body };
696 __la_stack_body($crate::Matrix::<$d>::zero())
697 }};
698}
699
700/// Common imports for ergonomic usage.
701///
702/// This prelude re-exports the primary types and constants: [`Matrix`],
703/// [`Vector`], [`Lu`], [`Ldlt`], [`Tolerance`], [`LaError`],
704/// [`DEFAULT_PIVOT_TOL`], [`DEFAULT_SINGULAR_TOL`], and the determinant error
705/// bound coefficients [`ERR_COEFF_2`], [`ERR_COEFF_3`], and [`ERR_COEFF_4`].
706/// It also re-exports [`MAX_STACK_MATRIX_DISPATCH_DIM`] and
707/// [`try_with_stack_matrix!`] for runtime-to-const matrix dispatch.
708///
709/// When the `exact` feature is enabled, `BigInt` and `BigRational` are also
710/// re-exported so callers can construct exact values (e.g. as the expected
711/// result of `Matrix::det_exact`) without adding `num-bigint` / `num-rational`
712/// to their own dependencies. The most commonly needed `num-traits` items are
713/// re-exported alongside them: `FromPrimitive` for `BigRational::from_f64` /
714/// `from_i64`, `ToPrimitive` for `BigRational::to_f64` / `to_i64`, and `Signed`
715/// for `.is_positive()` / `.is_negative()` / `.abs()`.
716pub mod prelude {
717 pub use crate::{
718 DEFAULT_PIVOT_TOL, DEFAULT_SINGULAR_TOL, ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4, LaError,
719 Ldlt, Lu, MAX_STACK_MATRIX_DISPATCH_DIM, Matrix, Tolerance, Vector, try_with_stack_matrix,
720 };
721
722 #[cfg(feature = "exact")]
723 pub use crate::{BigInt, BigRational, FromPrimitive, Signed, ToPrimitive};
724}
725
726#[cfg(test)]
727mod tests {
728 use core::assert_matches;
729
730 use super::*;
731
732 use approx::assert_abs_diff_eq;
733
734 #[test]
735 fn default_singular_tol_is_expected() {
736 assert_abs_diff_eq!(DEFAULT_SINGULAR_TOL.get(), 1e-12, epsilon = 0.0);
737 assert_abs_diff_eq!(
738 DEFAULT_PIVOT_TOL.get(),
739 DEFAULT_SINGULAR_TOL.get(),
740 epsilon = 0.0
741 );
742 }
743
744 #[test]
745 fn tolerance_new_accepts_finite_non_negative_values() {
746 assert_eq!(
747 Tolerance::new(0.0).unwrap().get().to_bits(),
748 0.0f64.to_bits()
749 );
750 assert_eq!(
751 Tolerance::new(1e-12).unwrap().get().to_bits(),
752 1e-12f64.to_bits()
753 );
754 assert_eq!(
755 Tolerance::new(f64::MAX).unwrap().get().to_bits(),
756 f64::MAX.to_bits()
757 );
758 }
759
760 #[test]
761 fn tolerance_new_rejects_negative_nan_and_infinity() {
762 assert_eq!(
763 Tolerance::new(-1.0),
764 Err(LaError::InvalidTolerance { value: -1.0 })
765 );
766 assert_matches!(
767 Tolerance::new(f64::NAN),
768 Err(LaError::InvalidTolerance { value }) if value.is_nan()
769 );
770 assert_eq!(
771 Tolerance::new(f64::INFINITY),
772 Err(LaError::InvalidTolerance {
773 value: f64::INFINITY,
774 })
775 );
776 assert_eq!(
777 Tolerance::new(f64::NEG_INFINITY),
778 Err(LaError::InvalidTolerance {
779 value: f64::NEG_INFINITY,
780 })
781 );
782 }
783
784 #[test]
785 fn validate_tolerance_matches_tolerance_new() {
786 for value in [0.0, 1e-12, f64::MAX] {
787 assert_eq!(LaError::validate_tolerance(value), Tolerance::new(value));
788 }
789
790 assert_eq!(
791 LaError::validate_tolerance(-1.0),
792 Err(LaError::InvalidTolerance { value: -1.0 })
793 );
794 assert_matches!(
795 LaError::validate_tolerance(f64::NAN),
796 Err(LaError::InvalidTolerance { value }) if value.is_nan()
797 );
798 assert_eq!(
799 LaError::validate_tolerance(f64::INFINITY),
800 Err(LaError::InvalidTolerance {
801 value: f64::INFINITY,
802 })
803 );
804 assert_eq!(
805 LaError::validate_tolerance(f64::NEG_INFINITY),
806 Err(LaError::InvalidTolerance {
807 value: f64::NEG_INFINITY,
808 })
809 );
810 }
811
812 #[test]
813 fn laerror_display_formats_singular() {
814 let err = LaError::Singular { pivot_col: 3 };
815 assert_eq!(err.to_string(), "singular matrix at pivot column 3");
816 }
817
818 #[test]
819 fn laerror_display_formats_nonfinite_with_row() {
820 let err = LaError::NonFinite {
821 row: Some(1),
822 col: 2,
823 };
824 assert_eq!(err.to_string(), "non-finite value at (1, 2)");
825 }
826
827 #[test]
828 fn laerror_display_formats_nonfinite_without_row() {
829 let err = LaError::NonFinite { row: None, col: 3 };
830 assert_eq!(err.to_string(), "non-finite value at index 3");
831 }
832
833 #[test]
834 fn laerror_display_formats_overflow() {
835 let err = LaError::Overflow { index: None };
836 assert_eq!(
837 err.to_string(),
838 "exact result overflows the target representation"
839 );
840 }
841
842 #[test]
843 fn laerror_display_formats_overflow_with_index() {
844 let err = LaError::Overflow { index: Some(2) };
845 assert_eq!(
846 err.to_string(),
847 "exact result overflows the target representation at index 2"
848 );
849 }
850
851 #[test]
852 fn laerror_display_formats_unsupported_dimension() {
853 let err = LaError::UnsupportedDimension {
854 requested: 8,
855 max: MAX_STACK_MATRIX_DISPATCH_DIM,
856 };
857 assert_eq!(
858 err.to_string(),
859 "unsupported matrix dimension 8; maximum stack-dispatch dimension is 7"
860 );
861 }
862
863 #[test]
864 fn laerror_display_formats_index_out_of_bounds() {
865 let err = LaError::IndexOutOfBounds {
866 row: 3,
867 col: 0,
868 dim: 3,
869 };
870 assert_eq!(
871 err.to_string(),
872 "matrix index (3, 0) is out of bounds for dimension 3"
873 );
874 }
875
876 #[test]
877 fn laerror_display_formats_invalid_tolerance() {
878 let err = LaError::InvalidTolerance { value: -1.0 };
879 assert_eq!(
880 err.to_string(),
881 "invalid tolerance -1; expected finite value >= 0"
882 );
883 }
884
885 #[test]
886 fn laerror_display_formats_asymmetric() {
887 let err = LaError::Asymmetric {
888 row: 0,
889 col: 2,
890 dim: 3,
891 };
892 assert_eq!(
893 err.to_string(),
894 "matrix is not symmetric for dimension 3: asymmetric pair (0, 2)"
895 );
896 }
897
898 #[test]
899 fn laerror_display_formats_not_positive_semidefinite() {
900 let err = LaError::NotPositiveSemidefinite {
901 pivot_col: 1,
902 value: -3.0,
903 };
904 assert_eq!(
905 err.to_string(),
906 "matrix is not positive semidefinite at LDLT pivot column 1: diagonal value -3 < 0"
907 );
908 }
909
910 #[test]
911 fn laerror_is_std_error_with_no_source() {
912 let err = LaError::Singular { pivot_col: 0 };
913 let e: &dyn std::error::Error = &err;
914 assert!(e.source().is_none());
915 }
916
917 #[test]
918 fn prelude_reexports_compile_and_work() {
919 use crate::prelude::*;
920
921 // Use the items so we know they are in scope and usable.
922 let m = Matrix::<2>::identity();
923 let v = Vector::<2>::new([1.0, 2.0]);
924 let _ = m.lu(DEFAULT_PIVOT_TOL).unwrap().solve_vec(v).unwrap();
925 let _ = m.ldlt(DEFAULT_SINGULAR_TOL).unwrap().solve_vec(v).unwrap();
926 assert_eq!(MAX_STACK_MATRIX_DISPATCH_DIM, 7);
927 }
928
929 macro_rules! gen_stack_matrix_dispatch_tests {
930 ($d:literal) => {
931 pastey::paste! {
932 #[test]
933 fn [<try_with_stack_matrix_dispatches_ $d d>]() {
934 let requested = $d;
935 let got = try_with_stack_matrix!(requested, |mut m| -> Result<usize, LaError> {
936 if $d > 0 {
937 m.set_checked($d - 1, $d - 1, f64::from($d))?;
938 assert_abs_diff_eq!(
939 m.get_checked($d - 1, $d - 1)?,
940 f64::from($d),
941 epsilon = 0.0
942 );
943 }
944 Ok($d)
945 });
946
947 assert_eq!(got, Ok($d));
948 }
949 }
950 };
951 }
952
953 gen_stack_matrix_dispatch_tests!(2);
954 gen_stack_matrix_dispatch_tests!(3);
955 gen_stack_matrix_dispatch_tests!(4);
956 gen_stack_matrix_dispatch_tests!(5);
957 gen_stack_matrix_dispatch_tests!(6);
958 gen_stack_matrix_dispatch_tests!(7);
959
960 #[test]
961 fn try_with_stack_matrix_supports_zero_dimension() {
962 let got = try_with_stack_matrix!(0usize, |m| -> Result<Option<f64>, LaError> {
963 m.det_direct()
964 });
965
966 assert_eq!(got, Ok(Some(1.0)));
967 }
968
969 #[test]
970 fn try_with_stack_matrix_reports_unsupported_dimension() {
971 let got = try_with_stack_matrix!(8usize, |m| -> Result<f64, LaError> { m.det() });
972
973 assert_eq!(
974 got,
975 Err(LaError::UnsupportedDimension {
976 requested: 8,
977 max: MAX_STACK_MATRIX_DISPATCH_DIM,
978 })
979 );
980 }
981
982 #[derive(Debug, PartialEq)]
983 struct DownstreamError(LaError);
984
985 impl From<LaError> for DownstreamError {
986 fn from(err: LaError) -> Self {
987 Self(err)
988 }
989 }
990
991 #[test]
992 fn try_with_stack_matrix_converts_unsupported_dimension_error() {
993 let got = try_with_stack_matrix!(9usize, |m| -> Result<usize, DownstreamError> {
994 assert_abs_diff_eq!(m.inf_norm().unwrap(), 0.0, epsilon = 0.0);
995 Ok(0)
996 });
997
998 assert_eq!(
999 got,
1000 Err(DownstreamError(LaError::UnsupportedDimension {
1001 requested: 9,
1002 max: MAX_STACK_MATRIX_DISPATCH_DIM,
1003 }))
1004 );
1005 }
1006
1007 /// Exercise every exact-feature re-export via the prelude so a future
1008 /// refactor that drops one (e.g. removing `Signed` from the prelude
1009 /// list) fails to compile rather than silently breaking downstream.
1010 #[cfg(feature = "exact")]
1011 #[test]
1012 fn prelude_exact_reexports_compile_and_work() {
1013 use crate::prelude::*;
1014
1015 // `BigInt` and `BigRational` constructors.
1016 let n = BigInt::from(7);
1017 let r = BigRational::from_integer(n.clone());
1018 assert_eq!(*r.numer(), n);
1019
1020 // `FromPrimitive::from_f64` / `from_i64` on `BigRational`.
1021 let half = BigRational::from_f64(0.5).unwrap();
1022 let two = BigRational::from_i64(2).unwrap();
1023 assert_eq!(
1024 half.clone() + half.clone(),
1025 BigRational::from_integer(BigInt::from(1))
1026 );
1027
1028 // `Signed::is_positive` / `is_negative` / `abs`.
1029 assert!(half.is_positive());
1030 assert!(!half.is_negative());
1031 let neg = -half.clone();
1032 assert!(neg.is_negative());
1033 assert_eq!(neg.abs(), half);
1034
1035 // `ToPrimitive::to_f64` / `to_i64`.
1036 assert!((half.to_f64().unwrap() - 0.5).abs() <= f64::EPSILON);
1037 assert_eq!(two.to_i64().unwrap(), 2);
1038 }
1039}