Skip to main content

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