feat: const-ify Lu/Ldlt det + solve_vec and Matrix inf_norm + det_errbound

Make six public methods `const fn`, completing const-evaluation parity
with `Matrix::det_direct` and `Vector::dot` now that MSRV 1.95 exposes
`f64::mul_add` as const (1.94) and `core::hint::cold_path` as const
(1.95):

- `Lu::det`, `Lu::solve_vec`
- `Ldlt::det`, `Ldlt::solve_vec`
- `Matrix::inf_norm`, `Matrix::det_errbound`

Iterator chains (`.iter().map().sum()`, `.enumerate().take(i)`,
`.enumerate().skip(i + 1)`, `.iter_mut().enumerate().take(D)`) were
rewritten as `while` loops since they are not const-stable.

Fix error-variant correctness in both solve_vec paths:

- A corrupt stored `U` / `D` diagonal at `(i, i)` now surfaces as
  `LaError::NonFinite { row: Some(i), col: i }`, matching the
  convention used by `Matrix::det`, `Lu::factor`, and `Ldlt::factor`.
- A computed-intermediate overflow keeps `row: None, col: i`.
- Previously both were conflated into `row: None, col: i`, defeating
  debuggability for callers who construct factorizations directly.

Sharpen `LaError::NonFinite` variant docs and the `# Errors` sections of
`Lu::solve_vec` / `Ldlt::solve_vec` to spell out the `(row, col)`
contract for each failure mode.

Cleanup:

- Import `ERR_COEFF_{2,3,4}` at `src/matrix.rs` top; drop `crate::`
  prefixes inside `Matrix::det_errbound`.
- Rename intermediate bindings `q` / `v` → `quotient` to stay under
  `clippy::many_single_char_names`.
- Rename const-evaluability tests from `*_is_const_evaluable_*` to
  `*_const_eval_*` for a concise, consistent style.

Tests:

- Added `{lu,ldlt,inf_norm,det_errbound}_const_eval_*` tests that force
  compile-time evaluation inside `const { … }` initializers.
- Added `solve_vec_defensive_non_finite_diagonal_{2,3,4,5}d` in
  `src/lu.rs` covering the new split error path (previously only the
  `Singular` defensive path was exercised).
- Updated matching `Ldlt` defensive tests to assert the new
  `row: Some(D - 1), col: D - 1` payload.

Validation: `just ci` passes — 164 unit tests, 324 integration/feature
tests, 28 doctests; clippy pedantic + nursery as errors; all examples.

Co-Authored-By: Oz <oz-agent@warp.dev>

Author: acgetchell

src/ldlt.rs

@@ -124,10 +124,12 @@ impl<const D: usize> Ldlt<D> {
     /// ```
     #[inline]
     #[must_use]
-    pub fn det(&self) -> f64 {
+    pub const fn det(&self) -> f64 {
         let mut det = 1.0;
-        for i in 0..D {
+        let mut i = 0;
+        while i < D {
             det *= self.factors.rows[i][i];
+            i += 1;
         }
         det
     }
@@ -154,57 +156,82 @@ impl<const D: usize> Ldlt<D> {
     /// # Errors
     /// Returns [`LaError::Singular`] if a diagonal entry `d = D[i,i]` satisfies `d <= tol`
     /// (non-positive or too small), where `tol` is the tolerance that was used during factorization.
-    /// Returns [`LaError::NonFinite`] if NaN/∞ is detected.
+    ///
+    /// Returns [`LaError::NonFinite`] if NaN/∞ is detected. The `row`/`col` coordinates
+    /// follow the convention documented on [`LaError::NonFinite`]:
+    ///
+    /// - `row: Some(i), col: i` — the stored `D` diagonal at `(i, i)` is non-finite
+    ///   (only reachable via direct `Ldlt` construction; [`Matrix::ldlt`](crate::Matrix::ldlt)
+    ///   rejects such factorizations).
+    /// - `row: None, col: i` — a computed intermediate (forward/back-substitution
+    ///   accumulator or the quotient `x[i] / diag`) overflowed to NaN/∞ at step `i`.
     #[inline]
-    pub fn solve_vec(&self, b: Vector<D>) -> Result<Vector<D>, LaError> {
+    pub const fn solve_vec(&self, b: Vector<D>) -> Result<Vector<D>, LaError> {
         let mut x = b.data;
 
         // Forward substitution: L y = b (L has unit diagonal).
-        for i in 0..D {
+        let mut i = 0;
+        while i < D {
             let mut sum = x[i];
             let row = self.factors.rows[i];
-            for (j, x_j) in x.iter().enumerate().take(i) {
-                sum = (-row[j]).mul_add(*x_j, sum);
+            let mut j = 0;
+            while j < i {
+                sum = (-row[j]).mul_add(x[j], sum);
+                j += 1;
             }
             if !sum.is_finite() {
                 cold_path();
                 return Err(LaError::NonFinite { row: None, col: i });
             }
             x[i] = sum;
+            i += 1;
         }
 
         // Diagonal solve: D z = y.
-        for (i, x_i) in x.iter_mut().enumerate().take(D) {
+        let mut i = 0;
+        while i < D {
             let diag = self.factors.rows[i][i];
+            // A corrupt stored diagonal is a specific matrix cell (i, i),
+            // distinct from a computed overflow — report it with
+            // `row: Some(i)` per the `LaError::NonFinite` convention used by
+            // `Matrix::det`, `Lu::factor`, and `Ldlt::factor`.
             if !diag.is_finite() {
                 cold_path();
-                return Err(LaError::NonFinite { row: None, col: i });
+                return Err(LaError::NonFinite {
+                    row: Some(i),
+                    col: i,
+                });
             }
             if diag <= self.tol {
                 cold_path();
                 return Err(LaError::Singular { pivot_col: i });
             }
 
-            let v = *x_i / diag;
-            if !v.is_finite() {
+            let quotient = x[i] / diag;
+            if !quotient.is_finite() {
                 cold_path();
                 return Err(LaError::NonFinite { row: None, col: i });
             }
-            *x_i = v;
+            x[i] = quotient;
+            i += 1;
         }
 
         // Back substitution: Lᵀ x = z.
-        for ii in 0..D {
+        let mut ii = 0;
+        while ii < D {
             let i = D - 1 - ii;
             let mut sum = x[i];
-            for (j, x_j) in x.iter().enumerate().skip(i + 1) {
-                sum = (-self.factors.rows[j][i]).mul_add(*x_j, sum);
+            let mut j = i + 1;
+            while j < D {
+                sum = (-self.factors.rows[j][i]).mul_add(x[j], sum);
+                j += 1;
             }
             if !sum.is_finite() {
                 cold_path();
                 return Err(LaError::NonFinite { row: None, col: i });
             }
             x[i] = sum;
+            ii += 1;
         }
 
         Ok(Vector::new(x))
@@ -487,7 +514,9 @@ mod tests {
             paste! {
                 /// `solve_vec` must surface `NonFinite` when a stored
                 /// diagonal is NaN, even though `factor` cannot produce
-                /// such a factorization.
+                /// such a factorization. The error must pinpoint the
+                /// corrupt cell at `(D-1, D-1)` per the
+                /// [`LaError::NonFinite`] convention.
                 #[test]
                 fn [<solve_vec_defensive_non_finite_diagonal_ $d d>]() {
                     let mut factors = Matrix::<$d>::identity();
@@ -498,7 +527,13 @@ mod tests {
                     };
                     let b = Vector::<$d>::new([1.0; $d]);
                     let err = ldlt.solve_vec(b).unwrap_err();
-                    assert_eq!(err, LaError::NonFinite { row: None, col: $d - 1 });
+                    assert_eq!(
+                        err,
+                        LaError::NonFinite {
+                            row: Some($d - 1),
+                            col: $d - 1,
+                        }
+                    );
                 }
 
                 /// `solve_vec` must surface `Singular` when a stored
@@ -524,4 +559,65 @@ mod tests {
     gen_solve_vec_defensive_tests!(3);
     gen_solve_vec_defensive_tests!(4);
     gen_solve_vec_defensive_tests!(5);
+
+    // -----------------------------------------------------------------------
+    // Const-evaluability tests.
+    //
+    // These prove that `Ldlt::det` and `Ldlt::solve_vec` are truly `const fn`
+    // by forcing the compiler to evaluate them inside a `const` initializer.
+    // `Ldlt::factor` is not (yet) `const fn` because the rank-1 update loop
+    // uses array indexing patterns that still require non-const helpers on
+    // some toolchains; we therefore construct `Ldlt<D>` directly.
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn ldlt_det_const_eval_d2() {
+        const DET: f64 = {
+            // Diagonal D = [4.0, 0.25] ⇒ det = 1.0.
+            let mut factors = Matrix::<2>::identity();
+            factors.rows[0][0] = 4.0;
+            factors.rows[1][1] = 0.25;
+            let ldlt = Ldlt::<2> {
+                factors,
+                tol: DEFAULT_SINGULAR_TOL,
+            };
+            ldlt.det()
+        };
+        assert!((DET - 1.0).abs() <= 1e-12);
+    }
+
+    #[test]
+    fn ldlt_det_const_eval_d3() {
+        const DET: f64 = {
+            // Diagonal D = [2.0, 3.0, 5.0] ⇒ det = 30.0.
+            let mut factors = Matrix::<3>::identity();
+            factors.rows[0][0] = 2.0;
+            factors.rows[1][1] = 3.0;
+            factors.rows[2][2] = 5.0;
+            let ldlt = Ldlt::<3> {
+                factors,
+                tol: DEFAULT_SINGULAR_TOL,
+            };
+            ldlt.det()
+        };
+        assert!((DET - 30.0).abs() <= 1e-12);
+    }
+
+    #[test]
+    fn ldlt_solve_vec_const_eval_d2() {
+        // Identity factors ⇒ solve_vec returns the RHS untouched.
+        const X: [f64; 2] = {
+            let ldlt = Ldlt::<2> {
+                factors: Matrix::<2>::identity(),
+                tol: DEFAULT_SINGULAR_TOL,
+            };
+            let b = Vector::<2>::new([1.0, 2.0]);
+            match ldlt.solve_vec(b) {
+                Ok(v) => v.into_array(),
+                Err(_) => [0.0, 0.0],
+            }
+        };
+        assert!((X[0] - 1.0).abs() <= 1e-12);
+        assert!((X[1] - 2.0).abs() <= 1e-12);
+    }
 }

src/lib.rs

@@ -119,11 +119,22 @@ pub enum LaError {
         pivot_col: usize,
     },
     /// A non-finite value (NaN/∞) was encountered.
+    ///
+    /// The `(row, col)` coordinate follows a consistent convention across the crate:
+    ///
+    /// - `row: Some(r), col: c` — a *stored* matrix cell at `(r, c)` is non-finite.
+    ///   Used by `Matrix::det`, `Lu::factor`, `Ldlt::factor`, and the `solve_vec`
+    ///   paths when they detect a corrupt stored factor (only reachable via
+    ///   direct struct construction; `factor` itself rejects such inputs).
+    /// - `row: None, col: c` — the non-finite value is either a *vector input*
+    ///   entry at index `c`, or a *computed intermediate* at step `c`
+    ///   (e.g. an accumulator that overflowed during forward/back substitution).
     NonFinite {
-        /// Row of the non-finite entry (for matrix inputs), or `None` when
-        /// the error originates from a vector input or a computed intermediate.
+        /// Row of the non-finite entry for a stored matrix cell, or `None` for
+        /// a vector-input entry or a computed intermediate. See the variant
+        /// docs for the full convention.
         row: Option<usize>,
-        /// Column index (for matrix inputs), vector index, or factorization
+        /// Column index (stored cell), vector index, or factorization/solve
         /// step where the non-finite value was detected.
         col: usize,
     },