la_stack/
error.rs

1#![forbid(unsafe_code)]
2
3//! Typed error categories and helpers for linear algebra operations.
4
5use core::fmt;
6
7/// Arithmetic operation associated with a computation failure.
8///
9/// This identifies both non-finite intermediates/results through
10/// [`NonFiniteOrigin::Computation`] and finite-endpoint range exhaustion through
11/// [`LaError::IntervalRangeExhausted`].
12///
13/// # Examples
14/// ```
15/// use la_stack::ArithmeticOperation;
16///
17/// assert_eq!(ArithmeticOperation::LuSolve.to_string(), "LU solve");
18/// ```
19#[derive(Clone, Copy, Debug, PartialEq, Eq)]
20#[non_exhaustive]
21pub enum ArithmeticOperation {
22    /// Matrix infinity-norm calculation.
23    MatrixInfinityNorm,
24    /// Matrix symmetry validation.
25    SymmetryCheck,
26    /// LU factorization.
27    LuFactorization,
28    /// LDLT factorization.
29    LdltFactorization,
30    /// Forward or backward substitution with an LU factorization.
31    LuSolve,
32    /// Forward, diagonal, or backward substitution with an LDLT factorization.
33    LdltSolve,
34    /// Determinant calculation.
35    Determinant,
36    /// Determinant error-bound calculation.
37    DeterminantErrorBound,
38    /// Outward-rounded interval addition.
39    IntervalAddition,
40    /// Exact-real subtraction enclosed by an outward-rounded interval.
41    IntervalSubtraction,
42    /// Outward-rounded interval multiplication.
43    IntervalMultiplication,
44    /// Outward-rounded interval square.
45    IntervalSquare,
46    /// Division-free interval determinant calculation.
47    IntervalDeterminant,
48    /// Vector dot-product calculation.
49    VectorDotProduct,
50    /// `axis · (left - right)` without first rounding the vector difference.
51    VectorDotDifference,
52    /// Vector squared-norm calculation.
53    VectorSquaredNorm,
54    /// Overflow- and underflow-safe vector Euclidean-norm calculation.
55    VectorNorm,
56    /// Unsigned angle between two vectors.
57    VectorAngle,
58}
59
60impl fmt::Display for ArithmeticOperation {
61    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
62        f.write_str(match self {
63            Self::MatrixInfinityNorm => "matrix infinity norm",
64            Self::SymmetryCheck => "symmetry check",
65            Self::LuFactorization => "LU factorization",
66            Self::LdltFactorization => "LDLT factorization",
67            Self::LuSolve => "LU solve",
68            Self::LdltSolve => "LDLT solve",
69            Self::Determinant => "determinant",
70            Self::DeterminantErrorBound => "determinant error bound",
71            Self::IntervalAddition => "interval addition",
72            Self::IntervalSubtraction => "interval subtraction",
73            Self::IntervalMultiplication => "interval multiplication",
74            Self::IntervalSquare => "interval square",
75            Self::IntervalDeterminant => "interval determinant",
76            Self::VectorDotProduct => "vector dot product",
77            Self::VectorDotDifference => "vector dot difference",
78            Self::VectorSquaredNorm => "vector squared norm",
79            Self::VectorNorm => "vector Euclidean norm",
80            Self::VectorAngle => "vector angle",
81        })
82    }
83}
84
85/// Factorization whose pivot policy rejected a matrix as numerically singular.
86///
87/// # Examples
88/// ```
89/// use la_stack::FactorizationKind;
90///
91/// assert_eq!(FactorizationKind::Ldlt.to_string(), "LDLT");
92/// ```
93#[derive(Clone, Copy, Debug, PartialEq, Eq)]
94#[non_exhaustive]
95pub enum FactorizationKind {
96    /// LU factorization with partial pivoting.
97    Lu,
98    /// LDLT factorization without pivoting.
99    Ldlt,
100}
101
102impl fmt::Display for FactorizationKind {
103    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
104        f.write_str(match self {
105            Self::Lu => "LU",
106            Self::Ldlt => "LDLT",
107        })
108    }
109}
110
111/// Reason a raw tolerance cannot become a [`crate::Tolerance`].
112///
113/// # Examples
114/// ```
115/// use la_stack::prelude::*;
116///
117/// match LaError::invalid_tolerance(-1.0) {
118///     LaError::InvalidTolerance {
119///         reason: InvalidToleranceReason::Negative,
120///         ..
121///     } => {}
122///     _ => unreachable!("a finite negative tolerance has a negative reason"),
123/// }
124/// ```
125#[derive(Clone, Copy, Debug, PartialEq, Eq)]
126#[non_exhaustive]
127pub enum InvalidToleranceReason {
128    /// The tolerance is finite but negative.
129    Negative,
130    /// The tolerance is NaN or positive/negative infinity.
131    NotFinite,
132}
133
134/// Endpoint of an interval constructor input.
135#[derive(Clone, Copy, Debug, PartialEq, Eq)]
136#[non_exhaustive]
137pub enum IntervalBound {
138    /// Lower endpoint.
139    Lower,
140    /// Upper endpoint.
141    Upper,
142}
143
144/// Operand of a binary scalar operation used to construct an interval.
145#[derive(Clone, Copy, Debug, PartialEq, Eq)]
146#[non_exhaustive]
147pub enum IntervalOperand {
148    /// Left operand.
149    Left,
150    /// Right operand.
151    Right,
152}
153
154/// Operand of a binary vector operation.
155///
156/// Identifies which input needs correction when [`LaError::ZeroVector`] or
157/// [`NonFiniteLocation::VectorOperandEntry`] reports an invalid direction or
158/// coordinate, including errors from [`crate::VectorAngle::angle`].
159///
160/// # Examples
161/// ```
162/// use la_stack::prelude::*;
163///
164/// match [1.0, 0.0].as_slice().angle(&[0.0, -0.0]) {
165///     Err(LaError::ZeroVector { operand, .. }) => {
166///         assert_eq!(operand, VectorOperand::Right);
167///     }
168///     _ => unreachable!("the right input is an all-zero vector"),
169/// }
170/// ```
171#[derive(Clone, Copy, Debug, PartialEq, Eq)]
172#[non_exhaustive]
173pub enum VectorOperand {
174    /// Left operand (`self` for a method).
175    Left,
176    /// Right operand (the method argument).
177    Right,
178}
179
180impl fmt::Display for VectorOperand {
181    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
182        f.write_str(match self {
183            Self::Left => "left",
184            Self::Right => "right",
185        })
186    }
187}
188
189/// Location at which a non-finite value was observed.
190///
191/// # Examples
192/// ```
193/// use la_stack::prelude::*;
194///
195/// match LaError::non_finite_input_matrix(1, 2) {
196///     LaError::NonFinite {
197///         location: NonFiniteLocation::MatrixCell { row, col, .. },
198///         ..
199///     } => assert_eq!((row, col), (1, 2)),
200///     _ => unreachable!("constructor returns a matrix-cell location"),
201/// }
202/// ```
203#[derive(Clone, Copy, Debug, PartialEq, Eq)]
204#[non_exhaustive]
205pub enum NonFiniteLocation {
206    /// Cell `(row, col)` in matrix-shaped storage or computation.
207    #[non_exhaustive]
208    MatrixCell {
209        /// Matrix row.
210        row: usize,
211        /// Matrix column.
212        col: usize,
213    },
214    /// Entry in a vector input.
215    #[non_exhaustive]
216    VectorEntry {
217        /// Vector index.
218        index: usize,
219    },
220    /// Entry in one of two borrowed vector operands.
221    #[non_exhaustive]
222    VectorOperandEntry {
223        /// Operand containing the non-finite coordinate.
224        operand: VectorOperand,
225        /// Coordinate index within that operand.
226        index: usize,
227    },
228    /// Indexed step in a factorization, solve, or reduction.
229    #[non_exhaustive]
230    Step {
231        /// Step index.
232        index: usize,
233    },
234    /// Endpoint supplied to an interval constructor.
235    #[non_exhaustive]
236    IntervalBound {
237        /// Bound at which the non-finite value was observed.
238        bound: IntervalBound,
239    },
240    /// Operand supplied to a binary scalar interval constructor.
241    #[non_exhaustive]
242    IntervalOperand {
243        /// Operand at which the non-finite value was observed.
244        operand: IntervalOperand,
245    },
246    /// Scalar value without a meaningful matrix or vector coordinate.
247    Scalar,
248}
249
250/// Provenance of a non-finite value.
251///
252/// # Examples
253/// ```
254/// use la_stack::prelude::*;
255///
256/// let err = LaError::non_finite_computation_scalar(ArithmeticOperation::Determinant);
257/// match err {
258///     LaError::NonFinite {
259///         origin: NonFiniteOrigin::Computation { operation, .. },
260///         ..
261///     } => assert_eq!(operation, ArithmeticOperation::Determinant),
262///     _ => unreachable!("constructor preserves computation provenance"),
263/// }
264/// ```
265#[derive(Clone, Copy, Debug, PartialEq, Eq)]
266#[non_exhaustive]
267pub enum NonFiniteOrigin {
268    /// The caller supplied a non-finite input.
269    Input,
270    /// Finite inputs produced a non-finite arithmetic result.
271    #[non_exhaustive]
272    Computation {
273        /// Operation that produced the value.
274        operation: ArithmeticOperation,
275    },
276}
277
278/// Computed LDLT condition that violates the no-pivot positive-semidefinite
279/// factorization requirements.
280///
281/// These values are computed in binary64. They explain why LDLT rejected the
282/// matrix, but do not prove that the stored matrix is exactly indefinite.
283///
284/// # Examples
285/// ```
286/// use la_stack::prelude::*;
287///
288/// match LaError::not_positive_semidefinite_negative(1, -3.0) {
289///     LaError::NotPositiveSemidefinite {
290///         violation: PositiveSemidefiniteViolation::NegativePivot { value, .. },
291///         ..
292///     } => assert_eq!(value, -3.0),
293///     _ => unreachable!("constructor preserves the PSD violation"),
294/// }
295/// ```
296#[derive(Clone, Copy, Debug, PartialEq)]
297#[non_exhaustive]
298pub enum PositiveSemidefiniteViolation {
299    /// LDLT produced a strictly negative diagonal pivot.
300    #[non_exhaustive]
301    NegativePivot {
302        /// Observed negative pivot value.
303        value: f64,
304    },
305    /// A zero diagonal pivot still has a non-zero coupling below it.
306    #[non_exhaustive]
307    ZeroPivotCoupling {
308        /// Row containing the non-zero coupling.
309        row: usize,
310        /// Observed coupling value.
311        value: f64,
312    },
313}
314
315/// Mathematical or numerical reason a matrix was classified as singular.
316///
317/// # Examples
318/// ```
319/// use la_stack::prelude::*;
320///
321/// match LaError::singular_numerical(1, FactorizationKind::Lu, 0.0, 1e-12) {
322///     LaError::Singular {
323///         reason: SingularityReason::Numerical { factorization, .. }, ..
324///     } => assert_eq!(factorization, FactorizationKind::Lu),
325///     _ => unreachable!("constructor preserves the singularity reason"),
326/// }
327/// ```
328#[derive(Clone, Copy, Debug, PartialEq)]
329#[non_exhaustive]
330pub enum SingularityReason {
331    /// The algorithm proved that the pivot is exactly zero.
332    Exact,
333    /// A finite pivot was rejected by the factorization tolerance.
334    #[non_exhaustive]
335    Numerical {
336        /// Factorization applying the tolerance policy.
337        factorization: FactorizationKind,
338        /// Absolute magnitude of the rejected pivot.
339        pivot_magnitude: f64,
340        /// Tolerance against which the magnitude was compared.
341        tolerance: f64,
342    },
343}
344
345/// Reason an exact result cannot satisfy an exact-to-`f64` conversion contract.
346///
347/// `RequiresRounding` is recoverable when the caller is willing to opt into a
348/// rounded exact-to-`f64` API. `NotFinite` means no finite `f64` can represent
349/// the result even after rounding.
350///
351/// # Examples
352/// ```
353/// use la_stack::prelude::*;
354///
355/// let err = LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding);
356/// assert!(err.requires_rounding());
357///
358/// let err = LaError::unrepresentable(None, UnrepresentableReason::NotFinite);
359/// assert_eq!(
360///     err.unrepresentable_reason(),
361///     Some(UnrepresentableReason::NotFinite)
362/// );
363/// ```
364#[derive(Clone, Copy, Debug, PartialEq, Eq)]
365#[non_exhaustive]
366pub enum UnrepresentableReason {
367    /// A finite `f64` exists only after rounding, but the requested conversion
368    /// requires an exact binary64 representation.
369    RequiresRounding,
370    /// No finite `f64` can represent the exact value after rounding.
371    NotFinite,
372}
373
374/// Linear algebra errors.
375///
376/// This enum and each struct-style variant are `#[non_exhaustive]` so downstream
377/// matches must retain a wildcard for future error categories and fields.
378///
379/// # Examples
380/// ```
381/// use la_stack::prelude::*;
382///
383/// match LaError::singular_exact(2) {
384///     LaError::Singular {
385///         pivot_col,
386///         reason: SingularityReason::Exact,
387///         ..
388///     } => assert_eq!(pivot_col, 2),
389///     _ => unreachable!("constructor returns an exact singularity"),
390/// }
391/// ```
392#[derive(Clone, Copy, Debug, PartialEq)]
393#[non_exhaustive]
394pub enum LaError {
395    /// Two vector operands have unequal ambient lengths.
396    #[non_exhaustive]
397    DimensionMismatch {
398        /// Left operand length.
399        left: usize,
400        /// Right operand length.
401        right: usize,
402    },
403    /// An operation requiring a direction received empty vector storage.
404    EmptyVector,
405    /// An operation requiring a direction received an all-zero vector.
406    #[non_exhaustive]
407    ZeroVector {
408        /// Operand whose coordinates are all zero (of either sign).
409        operand: VectorOperand,
410    },
411    /// A matrix is exactly or numerically singular.
412    #[non_exhaustive]
413    Singular {
414        /// Factorization column or step where a usable pivot was unavailable.
415        pivot_col: usize,
416        /// Typed reason for the singularity classification.
417        reason: SingularityReason,
418    },
419    /// A caller input or arithmetic intermediate/result is NaN or infinite.
420    #[non_exhaustive]
421    NonFinite {
422        /// Typed location of the value.
423        location: NonFiniteLocation,
424        /// Whether the value came from input or a particular computation.
425        origin: NonFiniteOrigin,
426    },
427    /// An exact-real interval intermediate or result cannot be enclosed by
428    /// finite binary64 endpoints.
429    #[non_exhaustive]
430    IntervalRangeExhausted {
431        /// Operation whose mathematical intermediate or result exceeded the
432        /// finite interval endpoint domain.
433        operation: ArithmeticOperation,
434    },
435    /// An exact result cannot satisfy the requested finite-`f64` conversion.
436    #[non_exhaustive]
437    Unrepresentable {
438        /// Failed vector component, or `None` for a scalar result.
439        index: Option<usize>,
440        /// Reason the conversion contract cannot be satisfied.
441        reason: UnrepresentableReason,
442    },
443    /// A finite interval's lower bound is greater than its upper bound.
444    #[non_exhaustive]
445    InvertedInterval {
446        /// Rejected lower bound.
447        lower: f64,
448        /// Rejected upper bound.
449        upper: f64,
450    },
451    /// Exact determinant scaling overflowed the internal exponent representation.
452    #[non_exhaustive]
453    DeterminantScaleOverflow {
454        /// Matrix dimension `D`.
455        dim: usize,
456        /// Minimum decomposed binary64 exponent among non-zero entries.
457        min_exponent: i32,
458    },
459    /// A matrix algorithm or runtime dispatch helper does not support a dimension.
460    #[non_exhaustive]
461    UnsupportedDimension {
462        /// Matrix dimension requested by the caller.
463        requested: usize,
464        /// Largest matrix dimension supported by the operation.
465        max: usize,
466    },
467    /// A matrix index is outside the `D×D` storage domain.
468    #[non_exhaustive]
469    IndexOutOfBounds {
470        /// Requested row.
471        row: usize,
472        /// Requested column.
473        col: usize,
474        /// Matrix dimension `D`; valid indices are less than this value.
475        dim: usize,
476    },
477    /// A raw tolerance is negative or non-finite.
478    #[non_exhaustive]
479    InvalidTolerance {
480        /// Raw value supplied by the caller.
481        value: f64,
482        /// Typed reason the value violates the tolerance invariant.
483        reason: InvalidToleranceReason,
484    },
485    /// A matrix required to be symmetric has an asymmetric off-diagonal pair.
486    #[non_exhaustive]
487    Asymmetric {
488        /// Row of the upper-triangular entry.
489        row: usize,
490        /// Column of the upper-triangular entry.
491        col: usize,
492        /// Matrix dimension `D`.
493        dim: usize,
494        /// Observed entry at `(row, col)`.
495        upper: f64,
496        /// Observed entry at `(col, row)`.
497        lower: f64,
498        /// Maximum absolute difference allowed by the symmetry check.
499        allowed_abs_diff: f64,
500    },
501    /// A computed LDLT pivot violates the no-pivot positive-semidefinite
502    /// factorization requirements.
503    ///
504    /// This diagnoses a binary64 factorization rejection; it is not an exact
505    /// certificate that the stored matrix is indefinite.
506    #[non_exhaustive]
507    NotPositiveSemidefinite {
508        /// LDLT pivot column or step where the violation was detected.
509        pivot_col: usize,
510        /// Typed PSD-domain violation.
511        violation: PositiveSemidefiniteViolation,
512    },
513}
514
515impl LaError {
516    /// Construct a [`LaError::Singular`] error proving that the pivot at
517    /// `pivot_col` is exactly zero.
518    #[inline]
519    #[must_use]
520    pub const fn singular_exact(pivot_col: usize) -> Self {
521        Self::Singular {
522            pivot_col,
523            reason: SingularityReason::Exact,
524        }
525    }
526
527    /// Construct a [`LaError::Singular`] error for a pivot rejected by a
528    /// factorization tolerance, preserving the factorization, pivot magnitude,
529    /// and tolerance in [`SingularityReason::Numerical`].
530    #[inline]
531    #[must_use]
532    pub const fn singular_numerical(
533        pivot_col: usize,
534        factorization: FactorizationKind,
535        pivot_magnitude: f64,
536        tolerance: f64,
537    ) -> Self {
538        Self::Singular {
539            pivot_col,
540            reason: SingularityReason::Numerical {
541                factorization,
542                pivot_magnitude,
543                tolerance,
544            },
545        }
546    }
547
548    /// Construct a [`LaError::NonFinite`] input error located at matrix cell
549    /// `(row, col)`.
550    #[inline]
551    #[must_use]
552    pub const fn non_finite_input_matrix(row: usize, col: usize) -> Self {
553        Self::NonFinite {
554            location: NonFiniteLocation::MatrixCell { row, col },
555            origin: NonFiniteOrigin::Input,
556        }
557    }
558
559    /// Construct a [`LaError::NonFinite`] input error located at vector entry
560    /// `index`.
561    #[inline]
562    #[must_use]
563    pub const fn non_finite_input_vector(index: usize) -> Self {
564        Self::NonFinite {
565            location: NonFiniteLocation::VectorEntry { index },
566            origin: NonFiniteOrigin::Input,
567        }
568    }
569
570    /// Construct a [`LaError::NonFinite`] input error for a scalar without an
571    /// index or matrix coordinate.
572    #[inline]
573    #[must_use]
574    pub const fn non_finite_input_scalar() -> Self {
575        Self::NonFinite {
576            location: NonFiniteLocation::Scalar,
577            origin: NonFiniteOrigin::Input,
578        }
579    }
580
581    /// Construct a [`LaError::NonFinite`] input error for a specific interval
582    /// endpoint.
583    #[inline]
584    #[must_use]
585    pub const fn non_finite_input_interval_bound(bound: IntervalBound) -> Self {
586        Self::NonFinite {
587            location: NonFiniteLocation::IntervalBound { bound },
588            origin: NonFiniteOrigin::Input,
589        }
590    }
591
592    /// Construct a [`LaError::NonFinite`] input error for a specific binary
593    /// scalar operand.
594    #[inline]
595    #[must_use]
596    pub const fn non_finite_input_interval_operand(operand: IntervalOperand) -> Self {
597        Self::NonFinite {
598            location: NonFiniteLocation::IntervalOperand { operand },
599            origin: NonFiniteOrigin::Input,
600        }
601    }
602
603    /// Construct a [`LaError::NonFinite`] computation error at matrix cell
604    /// `(row, col)`, retaining the originating `operation`.
605    #[inline]
606    #[must_use]
607    pub const fn non_finite_computation_matrix(
608        operation: ArithmeticOperation,
609        row: usize,
610        col: usize,
611    ) -> Self {
612        Self::NonFinite {
613            location: NonFiniteLocation::MatrixCell { row, col },
614            origin: NonFiniteOrigin::Computation { operation },
615        }
616    }
617
618    /// Construct a [`LaError::NonFinite`] computation error at `index`,
619    /// retaining the originating `operation`.
620    #[inline]
621    #[must_use]
622    pub const fn non_finite_computation_step(operation: ArithmeticOperation, index: usize) -> Self {
623        Self::NonFinite {
624            location: NonFiniteLocation::Step { index },
625            origin: NonFiniteOrigin::Computation { operation },
626        }
627    }
628
629    /// Construct a scalar [`LaError::NonFinite`] computation error retaining
630    /// the originating `operation`.
631    #[inline]
632    #[must_use]
633    pub const fn non_finite_computation_scalar(operation: ArithmeticOperation) -> Self {
634        Self::NonFinite {
635            location: NonFiniteLocation::Scalar,
636            origin: NonFiniteOrigin::Computation { operation },
637        }
638    }
639
640    /// Construct an [`LaError::IntervalRangeExhausted`] failure retaining the
641    /// responsible interval operation.
642    #[inline]
643    #[must_use]
644    pub const fn interval_range_exhausted(operation: ArithmeticOperation) -> Self {
645        Self::IntervalRangeExhausted { operation }
646    }
647
648    /// Construct a [`LaError::Unrepresentable`] conversion failure for a scalar
649    /// (`index = None`) or vector component (`index = Some(_)`).
650    #[inline]
651    #[must_use]
652    pub const fn unrepresentable(index: Option<usize>, reason: UnrepresentableReason) -> Self {
653        Self::Unrepresentable { index, reason }
654    }
655
656    /// Construct a [`LaError::InvertedInterval`] error preserving both rejected
657    /// finite bounds.
658    #[inline]
659    #[must_use]
660    pub const fn inverted_interval(lower: f64, upper: f64) -> Self {
661        Self::InvertedInterval { lower, upper }
662    }
663
664    /// Return the typed exact-to-`f64` conversion reason, or `None` for every
665    /// other error variant.
666    #[inline]
667    #[must_use]
668    pub const fn unrepresentable_reason(&self) -> Option<UnrepresentableReason> {
669        match self {
670            Self::Unrepresentable { reason, .. } => Some(*reason),
671            _ => None,
672        }
673    }
674
675    /// Return whether this is a `RequiresRounding` conversion failure for which
676    /// retrying with an explicit rounded API may succeed.
677    #[inline]
678    #[must_use]
679    pub const fn requires_rounding(&self) -> bool {
680        matches!(
681            self,
682            Self::Unrepresentable {
683                reason: UnrepresentableReason::RequiresRounding,
684                ..
685            }
686        )
687    }
688
689    /// Construct a [`LaError::DeterminantScaleOverflow`] retaining the matrix
690    /// dimension and minimum decomposed entry exponent.
691    #[inline]
692    #[must_use]
693    pub const fn determinant_scale_overflow(dim: usize, min_exponent: i32) -> Self {
694        Self::DeterminantScaleOverflow { dim, min_exponent }
695    }
696
697    /// Construct a [`LaError::UnsupportedDimension`] retaining the requested
698    /// and maximum supported dimensions.
699    #[inline]
700    #[must_use]
701    pub const fn unsupported_dimension(requested: usize, max: usize) -> Self {
702        Self::UnsupportedDimension { requested, max }
703    }
704
705    /// Construct a [`LaError::IndexOutOfBounds`] retaining the requested matrix
706    /// coordinates and dimension.
707    #[inline]
708    #[must_use]
709    pub const fn index_out_of_bounds(row: usize, col: usize, dim: usize) -> Self {
710        Self::IndexOutOfBounds { row, col, dim }
711    }
712
713    /// Construct an invalid-tolerance error and classify its typed reason.
714    ///
715    /// This low-level constructor assumes `value` has already failed the
716    /// tolerance invariant; raw caller input should normally be parsed through
717    /// [`crate::Tolerance::try_new`].
718    /// Non-finiteness takes precedence over negativity, so negative infinity is
719    /// classified as [`InvalidToleranceReason::NotFinite`].
720    #[inline]
721    #[must_use]
722    pub const fn invalid_tolerance(value: f64) -> Self {
723        let reason = if value.is_finite() {
724            InvalidToleranceReason::Negative
725        } else {
726            InvalidToleranceReason::NotFinite
727        };
728        Self::InvalidTolerance { value, reason }
729    }
730
731    /// Construct a [`LaError::Asymmetric`] error for the pair `(row, col)` and
732    /// `(col, row)`, retaining both observed values and the effective absolute
733    /// difference bound.
734    #[inline]
735    #[must_use]
736    pub const fn asymmetric(
737        row: usize,
738        col: usize,
739        dim: usize,
740        upper: f64,
741        lower: f64,
742        allowed_abs_diff: f64,
743    ) -> Self {
744        Self::Asymmetric {
745            row,
746            col,
747            dim,
748            upper,
749            lower,
750            allowed_abs_diff,
751        }
752    }
753
754    /// Construct a [`LaError::NotPositiveSemidefinite`] error for a computed
755    /// negative LDLT diagonal pivot.
756    #[inline]
757    #[must_use]
758    pub const fn not_positive_semidefinite_negative(pivot_col: usize, value: f64) -> Self {
759        Self::NotPositiveSemidefinite {
760            pivot_col,
761            violation: PositiveSemidefiniteViolation::NegativePivot { value },
762        }
763    }
764
765    /// Construct a [`LaError::NotPositiveSemidefinite`] error for a computed
766    /// zero LDLT diagonal with a non-zero coupling at `row`, distinguishing the
767    /// observed factorization violation from an uncoupled singular pivot.
768    #[inline]
769    #[must_use]
770    pub const fn not_positive_semidefinite_zero_coupling(
771        pivot_col: usize,
772        row: usize,
773        value: f64,
774    ) -> Self {
775        Self::NotPositiveSemidefinite {
776            pivot_col,
777            violation: PositiveSemidefiniteViolation::ZeroPivotCoupling { row, value },
778        }
779    }
780}
781
782/// Write the structured location portion of [`LaError::NonFinite`]'s public
783/// display contract without allocating an intermediate string.
784///
785fn write_non_finite_location(
786    f: &mut fmt::Formatter<'_>,
787    location: NonFiniteLocation,
788) -> fmt::Result {
789    match location {
790        NonFiniteLocation::MatrixCell { row, col } => {
791            write!(f, "matrix cell ({row}, {col})")
792        }
793        NonFiniteLocation::VectorEntry { index } => write!(f, "vector entry {index}"),
794        NonFiniteLocation::VectorOperandEntry { operand, index } => {
795            write!(f, "{operand} vector entry {index}")
796        }
797        NonFiniteLocation::Step { index } => write!(f, "step {index}"),
798        NonFiniteLocation::IntervalBound {
799            bound: IntervalBound::Lower,
800        } => f.write_str("interval lower bound"),
801        NonFiniteLocation::IntervalBound {
802            bound: IntervalBound::Upper,
803        } => f.write_str("interval upper bound"),
804        NonFiniteLocation::IntervalOperand {
805            operand: IntervalOperand::Left,
806        } => f.write_str("left interval operand"),
807        NonFiniteLocation::IntervalOperand {
808            operand: IntervalOperand::Right,
809        } => f.write_str("right interval operand"),
810        NonFiniteLocation::Scalar => f.write_str("scalar value"),
811    }
812}
813
814/// Write a [`LaError::NonFinite`] message from its structured location and
815/// origin while distinguishing scalar input from a computed scalar result.
816fn write_non_finite(
817    f: &mut fmt::Formatter<'_>,
818    location: NonFiniteLocation,
819    origin: NonFiniteOrigin,
820) -> fmt::Result {
821    match (location, origin) {
822        (NonFiniteLocation::Scalar, NonFiniteOrigin::Input) => {
823            f.write_str("non-finite scalar input")
824        }
825        (NonFiniteLocation::Scalar, NonFiniteOrigin::Computation { operation }) => {
826            write!(f, "non-finite scalar result computed during {operation}")
827        }
828        (location, NonFiniteOrigin::Input) => {
829            f.write_str("non-finite input value at ")?;
830            write_non_finite_location(f, location)
831        }
832        (location, NonFiniteOrigin::Computation { operation }) => {
833            write!(f, "non-finite value computed during {operation} at ")?;
834            write_non_finite_location(f, location)
835        }
836    }
837}
838
839/// Write exact-conversion failures with their optional component context.
840fn write_unrepresentable(
841    f: &mut fmt::Formatter<'_>,
842    index: Option<usize>,
843    reason: UnrepresentableReason,
844) -> fmt::Result {
845    f.write_str(match reason {
846        UnrepresentableReason::RequiresRounding => {
847            "exact result requires rounding to fit finite f64"
848        }
849        UnrepresentableReason::NotFinite => {
850            "exact result has no finite f64 representation after rounding"
851        }
852    })?;
853    if let Some(index) = index {
854        write!(f, " at index {index}")?;
855    }
856    Ok(())
857}
858
859impl fmt::Display for LaError {
860    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
861        match *self {
862            Self::DimensionMismatch { left, right } => {
863                write!(
864                    f,
865                    "vector dimension mismatch: left length {left}, right length {right}"
866                )
867            }
868            Self::EmptyVector => f.write_str("an empty vector has no direction"),
869            Self::ZeroVector { operand } => write!(f, "the {operand} zero vector has no direction"),
870            Self::Singular {
871                pivot_col,
872                reason: SingularityReason::Exact,
873            } => write!(f, "matrix is exactly singular at pivot column {pivot_col}"),
874            Self::Singular {
875                pivot_col,
876                reason:
877                    SingularityReason::Numerical {
878                        factorization,
879                        pivot_magnitude,
880                        tolerance,
881                    },
882            } => write!(
883                f,
884                "matrix is numerically singular during {factorization} factorization at pivot column {pivot_col}: pivot magnitude {pivot_magnitude} <= tolerance {tolerance}"
885            ),
886            Self::NonFinite { location, origin } => write_non_finite(f, location, origin),
887            Self::IntervalRangeExhausted { operation } => write!(
888                f,
889                "exact-real intermediate or result of {operation} has no enclosure with finite binary64 endpoints"
890            ),
891            Self::Unrepresentable { index, reason } => write_unrepresentable(f, index, reason),
892            Self::InvertedInterval { lower, upper } => write!(
893                f,
894                "invalid interval bounds [{lower}, {upper}]; expected lower <= upper"
895            ),
896            Self::DeterminantScaleOverflow { dim, min_exponent } => write!(
897                f,
898                "exact determinant scale exponent overflows for dimension {dim} with minimum entry exponent {min_exponent}"
899            ),
900            Self::UnsupportedDimension { requested, max } => write!(
901                f,
902                "unsupported matrix dimension {requested}; maximum supported dimension is {max}"
903            ),
904            Self::IndexOutOfBounds { row, col, dim } => write!(
905                f,
906                "matrix index ({row}, {col}) is out of bounds for dimension {dim}"
907            ),
908            Self::InvalidTolerance {
909                value,
910                reason: InvalidToleranceReason::Negative,
911            } => write!(f, "invalid tolerance {value}; expected value >= 0"),
912            Self::InvalidTolerance {
913                value,
914                reason: InvalidToleranceReason::NotFinite,
915            } => write!(f, "invalid tolerance {value}; expected a finite value"),
916            Self::Asymmetric {
917                row,
918                col,
919                dim,
920                upper,
921                lower,
922                allowed_abs_diff,
923            } => write!(
924                f,
925                "matrix is not symmetric for dimension {dim}: entry ({row}, {col}) = {upper} and entry ({col}, {row}) = {lower} differ by more than allowed absolute difference {allowed_abs_diff}"
926            ),
927            Self::NotPositiveSemidefinite {
928                pivot_col,
929                violation: PositiveSemidefiniteViolation::NegativePivot { value },
930            } => write!(
931                f,
932                "LDLT rejected the matrix at pivot column {pivot_col}: computed diagonal value {value} < 0"
933            ),
934            Self::NotPositiveSemidefinite {
935                pivot_col,
936                violation: PositiveSemidefiniteViolation::ZeroPivotCoupling { row, value },
937            } => write!(
938                f,
939                "LDLT rejected the matrix at pivot column {pivot_col}: computed zero diagonal has non-zero coupling at row {row} with value {value}"
940            ),
941        }
942    }
943}
944
945impl std::error::Error for LaError {}
946
947#[cfg(test)]
948mod tests {
949    use std::error::Error;
950
951    use super::*;
952    use crate::MAX_STACK_MATRIX_DISPATCH_DIM;
953
954    #[test]
955    fn category_displays_are_concise() {
956        assert_eq!(FactorizationKind::Lu.to_string(), "LU");
957        assert_eq!(FactorizationKind::Ldlt.to_string(), "LDLT");
958        assert_eq!(
959            ArithmeticOperation::MatrixInfinityNorm.to_string(),
960            "matrix infinity norm"
961        );
962        assert_eq!(
963            ArithmeticOperation::SymmetryCheck.to_string(),
964            "symmetry check"
965        );
966        assert_eq!(
967            ArithmeticOperation::LuFactorization.to_string(),
968            "LU factorization"
969        );
970        assert_eq!(
971            ArithmeticOperation::LdltFactorization.to_string(),
972            "LDLT factorization"
973        );
974        assert_eq!(ArithmeticOperation::LuSolve.to_string(), "LU solve");
975        assert_eq!(ArithmeticOperation::LdltSolve.to_string(), "LDLT solve");
976        assert_eq!(ArithmeticOperation::Determinant.to_string(), "determinant");
977        assert_eq!(
978            ArithmeticOperation::DeterminantErrorBound.to_string(),
979            "determinant error bound"
980        );
981        assert_eq!(
982            ArithmeticOperation::VectorDotProduct.to_string(),
983            "vector dot product"
984        );
985        assert_eq!(
986            ArithmeticOperation::VectorDotDifference.to_string(),
987            "vector dot difference"
988        );
989        assert_eq!(
990            ArithmeticOperation::VectorSquaredNorm.to_string(),
991            "vector squared norm"
992        );
993        assert_eq!(
994            ArithmeticOperation::VectorNorm.to_string(),
995            "vector Euclidean norm"
996        );
997        assert_eq!(
998            ArithmeticOperation::IntervalAddition.to_string(),
999            "interval addition"
1000        );
1001        assert_eq!(
1002            ArithmeticOperation::IntervalSubtraction.to_string(),
1003            "interval subtraction"
1004        );
1005        assert_eq!(
1006            ArithmeticOperation::IntervalMultiplication.to_string(),
1007            "interval multiplication"
1008        );
1009        assert_eq!(
1010            ArithmeticOperation::IntervalSquare.to_string(),
1011            "interval square"
1012        );
1013        assert_eq!(
1014            ArithmeticOperation::IntervalDeterminant.to_string(),
1015            "interval determinant"
1016        );
1017    }
1018
1019    #[test]
1020    fn singular_constructors_and_displays_preserve_reason() {
1021        let exact = LaError::singular_exact(3);
1022        assert_eq!(
1023            exact,
1024            LaError::Singular {
1025                pivot_col: 3,
1026                reason: SingularityReason::Exact,
1027            }
1028        );
1029        assert_eq!(
1030            exact.to_string(),
1031            "matrix is exactly singular at pivot column 3"
1032        );
1033
1034        let numerical = LaError::singular_numerical(2, FactorizationKind::Lu, 1e-14, 1e-12);
1035        assert_eq!(
1036            numerical,
1037            LaError::Singular {
1038                pivot_col: 2,
1039                reason: SingularityReason::Numerical {
1040                    factorization: FactorizationKind::Lu,
1041                    pivot_magnitude: 1e-14,
1042                    tolerance: 1e-12,
1043                },
1044            }
1045        );
1046        assert_eq!(
1047            numerical.to_string(),
1048            "matrix is numerically singular during LU factorization at pivot column 2: pivot magnitude 0.00000000000001 <= tolerance 0.000000000001"
1049        );
1050    }
1051
1052    #[test]
1053    fn non_finite_constructors_preserve_location_and_origin() {
1054        assert_eq!(
1055            LaError::non_finite_input_matrix(1, 2),
1056            LaError::NonFinite {
1057                location: NonFiniteLocation::MatrixCell { row: 1, col: 2 },
1058                origin: NonFiniteOrigin::Input,
1059            }
1060        );
1061        assert_eq!(
1062            LaError::non_finite_input_vector(3).to_string(),
1063            "non-finite input value at vector entry 3"
1064        );
1065        assert_eq!(
1066            LaError::non_finite_input_scalar().to_string(),
1067            "non-finite scalar input"
1068        );
1069        assert_eq!(
1070            LaError::non_finite_input_interval_bound(IntervalBound::Upper).to_string(),
1071            "non-finite input value at interval upper bound"
1072        );
1073        assert_eq!(
1074            LaError::non_finite_input_interval_bound(IntervalBound::Lower).to_string(),
1075            "non-finite input value at interval lower bound"
1076        );
1077        assert_eq!(
1078            LaError::non_finite_input_interval_operand(IntervalOperand::Left).to_string(),
1079            "non-finite input value at left interval operand"
1080        );
1081        assert_eq!(
1082            LaError::non_finite_input_interval_operand(IntervalOperand::Right).to_string(),
1083            "non-finite input value at right interval operand"
1084        );
1085        assert_eq!(
1086            LaError::non_finite_computation_matrix(ArithmeticOperation::LuFactorization, 2, 1)
1087                .to_string(),
1088            "non-finite value computed during LU factorization at matrix cell (2, 1)"
1089        );
1090        assert_eq!(
1091            LaError::non_finite_computation_step(ArithmeticOperation::LuSolve, 1).to_string(),
1092            "non-finite value computed during LU solve at step 1"
1093        );
1094        assert_eq!(
1095            LaError::non_finite_computation_scalar(ArithmeticOperation::Determinant).to_string(),
1096            "non-finite scalar result computed during determinant"
1097        );
1098    }
1099
1100    #[test]
1101    fn unrepresentable_helpers_preserve_recovery_reason() {
1102        let rounding = LaError::unrepresentable(Some(2), UnrepresentableReason::RequiresRounding);
1103        let scalar_rounding =
1104            LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding);
1105        let indexed_not_finite =
1106            LaError::unrepresentable(Some(2), UnrepresentableReason::NotFinite);
1107        let not_finite = LaError::unrepresentable(None, UnrepresentableReason::NotFinite);
1108        assert_eq!(
1109            rounding.unrepresentable_reason(),
1110            Some(UnrepresentableReason::RequiresRounding)
1111        );
1112        assert!(rounding.requires_rounding());
1113        assert_eq!(
1114            rounding.to_string(),
1115            "exact result requires rounding to fit finite f64 at index 2"
1116        );
1117        assert_eq!(
1118            scalar_rounding.to_string(),
1119            "exact result requires rounding to fit finite f64"
1120        );
1121        assert_eq!(
1122            indexed_not_finite.to_string(),
1123            "exact result has no finite f64 representation after rounding at index 2"
1124        );
1125        assert_eq!(
1126            not_finite.to_string(),
1127            "exact result has no finite f64 representation after rounding"
1128        );
1129        assert!(!not_finite.requires_rounding());
1130        assert_eq!(LaError::singular_exact(0).unrepresentable_reason(), None);
1131    }
1132
1133    #[test]
1134    fn unrepresentable_display_propagates_writer_failure() {
1135        struct LimitedWrites {
1136            remaining: usize,
1137        }
1138
1139        impl fmt::Write for LimitedWrites {
1140            fn write_str(&mut self, _: &str) -> fmt::Result {
1141                if self.remaining == 0 {
1142                    return Err(fmt::Error);
1143                }
1144                self.remaining -= 1;
1145                Ok(())
1146            }
1147        }
1148
1149        for reason in [
1150            UnrepresentableReason::RequiresRounding,
1151            UnrepresentableReason::NotFinite,
1152        ] {
1153            for index in [None, Some(2)] {
1154                let error = LaError::unrepresentable(index, reason);
1155                for remaining in [0, 1] {
1156                    // One write accepts the base message, then rejects any
1157                    // optional index suffix. Zero rejects the base message.
1158                    let expected = if remaining == 1 && index.is_none() {
1159                        Ok(())
1160                    } else {
1161                        Err(fmt::Error)
1162                    };
1163                    assert_eq!(
1164                        fmt::write(&mut LimitedWrites { remaining }, format_args!("{error}")),
1165                        expected,
1166                        "reason={reason:?}, index={index:?}, remaining={remaining}"
1167                    );
1168                }
1169            }
1170        }
1171    }
1172
1173    #[test]
1174    fn invalid_tolerance_classifies_non_finite_before_negative() {
1175        assert_eq!(
1176            LaError::invalid_tolerance(-1.0),
1177            LaError::InvalidTolerance {
1178                value: -1.0,
1179                reason: InvalidToleranceReason::Negative,
1180            }
1181        );
1182        assert_eq!(
1183            LaError::invalid_tolerance(f64::NEG_INFINITY),
1184            LaError::InvalidTolerance {
1185                value: f64::NEG_INFINITY,
1186                reason: InvalidToleranceReason::NotFinite,
1187            }
1188        );
1189        assert_eq!(
1190            LaError::invalid_tolerance(-1.0).to_string(),
1191            "invalid tolerance -1; expected value >= 0"
1192        );
1193        assert_eq!(
1194            LaError::invalid_tolerance(f64::NEG_INFINITY).to_string(),
1195            "invalid tolerance -inf; expected a finite value"
1196        );
1197    }
1198
1199    #[test]
1200    fn inverted_interval_error_preserves_both_bounds() {
1201        let error = LaError::inverted_interval(2.0, 1.0);
1202        assert_eq!(
1203            error,
1204            LaError::InvertedInterval {
1205                lower: 2.0,
1206                upper: 1.0,
1207            }
1208        );
1209        assert_eq!(
1210            error.to_string(),
1211            "invalid interval bounds [2, 1]; expected lower <= upper"
1212        );
1213    }
1214
1215    #[test]
1216    fn interval_range_error_preserves_operation() {
1217        let error = LaError::interval_range_exhausted(ArithmeticOperation::IntervalSquare);
1218        assert_eq!(
1219            error,
1220            LaError::IntervalRangeExhausted {
1221                operation: ArithmeticOperation::IntervalSquare,
1222            }
1223        );
1224        assert_eq!(
1225            error.to_string(),
1226            "exact-real intermediate or result of interval square has no enclosure with finite binary64 endpoints"
1227        );
1228    }
1229
1230    #[test]
1231    fn asymmetric_error_retains_observed_values_and_bound() {
1232        let err = LaError::asymmetric(0, 2, 3, 1.0, 1.5, 1e-12);
1233        assert_eq!(
1234            err,
1235            LaError::Asymmetric {
1236                row: 0,
1237                col: 2,
1238                dim: 3,
1239                upper: 1.0,
1240                lower: 1.5,
1241                allowed_abs_diff: 1e-12,
1242            }
1243        );
1244        assert_eq!(
1245            err.to_string(),
1246            "matrix is not symmetric for dimension 3: entry (0, 2) = 1 and entry (2, 0) = 1.5 differ by more than allowed absolute difference 0.000000000001"
1247        );
1248    }
1249
1250    #[test]
1251    fn positive_semidefinite_errors_preserve_distinct_violations() {
1252        assert_eq!(
1253            LaError::not_positive_semidefinite_negative(1, -3.0).to_string(),
1254            "LDLT rejected the matrix at pivot column 1: computed diagonal value -3 < 0"
1255        );
1256        assert_eq!(
1257            LaError::not_positive_semidefinite_zero_coupling(0, 1, 2.0).to_string(),
1258            "LDLT rejected the matrix at pivot column 0: computed zero diagonal has non-zero coupling at row 1 with value 2"
1259        );
1260    }
1261
1262    #[test]
1263    fn remaining_helpers_and_displays_preserve_fields() {
1264        assert_eq!(
1265            LaError::determinant_scale_overflow(3, -1074).to_string(),
1266            "exact determinant scale exponent overflows for dimension 3 with minimum entry exponent -1074"
1267        );
1268        assert_eq!(
1269            LaError::unsupported_dimension(8, MAX_STACK_MATRIX_DISPATCH_DIM).to_string(),
1270            "unsupported matrix dimension 8; maximum supported dimension is 7"
1271        );
1272        assert_eq!(
1273            LaError::index_out_of_bounds(3, 0, 3).to_string(),
1274            "matrix index (3, 0) is out of bounds for dimension 3"
1275        );
1276    }
1277
1278    #[test]
1279    fn is_std_error_with_no_source() {
1280        let err = LaError::singular_exact(0);
1281        let error: &dyn Error = &err;
1282        assert!(error.source().is_none());
1283    }
1284}
la_stack/error.rs

la_stack/
error.rs
