la_stack/matrix.rs
1//! Fixed-size, stack-allocated square matrices.
2
3use crate::LaError;
4use crate::ldlt::Ldlt;
5use crate::lu::Lu;
6
7/// Fixed-size square matrix `D×D`, stored inline.
8#[must_use]
9#[derive(Clone, Copy, Debug, PartialEq)]
10pub struct Matrix<const D: usize> {
11 pub(crate) rows: [[f64; D]; D],
12}
13
14impl<const D: usize> Matrix<D> {
15 /// Construct from row-major storage.
16 ///
17 /// # Examples
18 /// ```
19 /// use la_stack::prelude::*;
20 ///
21 /// let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
22 /// assert_eq!(m.get(0, 1), Some(2.0));
23 /// ```
24 #[inline]
25 pub const fn from_rows(rows: [[f64; D]; D]) -> Self {
26 Self { rows }
27 }
28
29 /// All-zeros matrix.
30 ///
31 /// # Examples
32 /// ```
33 /// use la_stack::prelude::*;
34 ///
35 /// let z = Matrix::<2>::zero();
36 /// assert_eq!(z.get(1, 1), Some(0.0));
37 /// ```
38 #[inline]
39 pub const fn zero() -> Self {
40 Self {
41 rows: [[0.0; D]; D],
42 }
43 }
44
45 /// Identity matrix.
46 ///
47 /// # Examples
48 /// ```
49 /// use la_stack::prelude::*;
50 ///
51 /// let i = Matrix::<3>::identity();
52 /// assert_eq!(i.get(0, 0), Some(1.0));
53 /// assert_eq!(i.get(0, 1), Some(0.0));
54 /// assert_eq!(i.get(2, 2), Some(1.0));
55 /// ```
56 #[inline]
57 pub const fn identity() -> Self {
58 let mut m = Self::zero();
59
60 let mut i = 0;
61 while i < D {
62 m.rows[i][i] = 1.0;
63 i += 1;
64 }
65
66 m
67 }
68
69 /// Get an element with bounds checking.
70 ///
71 /// # Examples
72 /// ```
73 /// use la_stack::prelude::*;
74 ///
75 /// let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
76 /// assert_eq!(m.get(1, 0), Some(3.0));
77 /// assert_eq!(m.get(2, 0), None);
78 /// ```
79 #[inline]
80 #[must_use]
81 pub const fn get(&self, r: usize, c: usize) -> Option<f64> {
82 if r < D && c < D {
83 Some(self.rows[r][c])
84 } else {
85 None
86 }
87 }
88
89 /// Set an element with bounds checking.
90 ///
91 /// Returns `true` if the index was in-bounds.
92 ///
93 /// # Examples
94 /// ```
95 /// use la_stack::prelude::*;
96 ///
97 /// let mut m = Matrix::<2>::zero();
98 /// assert!(m.set(0, 1, 2.5));
99 /// assert_eq!(m.get(0, 1), Some(2.5));
100 /// assert!(!m.set(10, 0, 1.0));
101 /// ```
102 #[inline]
103 pub const fn set(&mut self, r: usize, c: usize, value: f64) -> bool {
104 if r < D && c < D {
105 self.rows[r][c] = value;
106 true
107 } else {
108 false
109 }
110 }
111
112 /// Infinity norm (maximum absolute row sum).
113 ///
114 /// # Examples
115 /// ```
116 /// use la_stack::prelude::*;
117 ///
118 /// let m = Matrix::<2>::from_rows([[1.0, -2.0], [3.0, 4.0]]);
119 /// assert!((m.inf_norm() - 7.0).abs() <= 1e-12);
120 /// ```
121 #[inline]
122 #[must_use]
123 pub fn inf_norm(&self) -> f64 {
124 let mut max_row_sum: f64 = 0.0;
125
126 for row in &self.rows {
127 let row_sum: f64 = row.iter().map(|&x| x.abs()).sum();
128 if row_sum > max_row_sum {
129 max_row_sum = row_sum;
130 }
131 }
132
133 max_row_sum
134 }
135
136 /// Compute an LU decomposition with partial pivoting.
137 ///
138 /// # Examples
139 /// ```
140 /// use la_stack::prelude::*;
141 ///
142 /// # fn main() -> Result<(), LaError> {
143 /// let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
144 /// let lu = a.lu(DEFAULT_PIVOT_TOL)?;
145 ///
146 /// let b = Vector::<2>::new([5.0, 11.0]);
147 /// let x = lu.solve_vec(b)?.into_array();
148 ///
149 /// assert!((x[0] - 1.0).abs() <= 1e-12);
150 /// assert!((x[1] - 2.0).abs() <= 1e-12);
151 /// # Ok(())
152 /// # }
153 /// ```
154 ///
155 /// # Errors
156 /// Returns [`LaError::Singular`] if, for some column `k`, the largest-magnitude candidate pivot
157 /// in that column satisfies `|pivot| <= tol` (so no numerically usable pivot exists).
158 /// Returns [`LaError::NonFinite`] if NaN/∞ is detected during factorization.
159 #[inline]
160 pub fn lu(self, tol: f64) -> Result<Lu<D>, LaError> {
161 Lu::factor(self, tol)
162 }
163
164 /// Compute an LDLT factorization (`A = L D Lᵀ`) without pivoting.
165 ///
166 /// This is intended for symmetric positive definite (SPD) and positive semi-definite (PSD)
167 /// matrices such as Gram matrices.
168 ///
169 /// # Examples
170 /// ```
171 /// use la_stack::prelude::*;
172 ///
173 /// # fn main() -> Result<(), LaError> {
174 /// let a = Matrix::<2>::from_rows([[4.0, 2.0], [2.0, 3.0]]);
175 /// let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?;
176 ///
177 /// // det(A) = 8
178 /// assert!((ldlt.det() - 8.0).abs() <= 1e-12);
179 ///
180 /// // Solve A x = b
181 /// let b = Vector::<2>::new([1.0, 2.0]);
182 /// let x = ldlt.solve_vec(b)?.into_array();
183 /// assert!((x[0] - (-0.125)).abs() <= 1e-12);
184 /// assert!((x[1] - 0.75).abs() <= 1e-12);
185 /// # Ok(())
186 /// # }
187 /// ```
188 ///
189 /// # Errors
190 /// Returns [`LaError::Singular`] if, for some step `k`, the required diagonal entry `d = D[k,k]`
191 /// is `<= tol` (non-positive or too small). This treats PSD degeneracy (and indefinite inputs)
192 /// as singular/degenerate.
193 /// Returns [`LaError::NonFinite`] if NaN/∞ is detected during factorization.
194 #[inline]
195 pub fn ldlt(self, tol: f64) -> Result<Ldlt<D>, LaError> {
196 Ldlt::factor(self, tol)
197 }
198
199 /// Closed-form determinant for dimensions 0–4, bypassing LU factorization.
200 ///
201 /// Returns `Some(det)` for `D` ∈ {0, 1, 2, 3, 4}, `None` for D ≥ 5.
202 /// `D = 0` returns `Some(1.0)` (empty product).
203 /// This is a `const fn` (Rust 1.94+) and uses fused multiply-add (`mul_add`)
204 /// for improved accuracy and performance.
205 ///
206 /// For a determinant that works for any dimension (falling back to LU for D ≥ 5),
207 /// use [`det`](Self::det).
208 ///
209 /// # Examples
210 /// ```
211 /// use la_stack::prelude::*;
212 ///
213 /// let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
214 /// assert!((m.det_direct().unwrap() - (-2.0)).abs() <= 1e-12);
215 ///
216 /// // D = 0 is the empty product.
217 /// assert_eq!(Matrix::<0>::zero().det_direct(), Some(1.0));
218 ///
219 /// // D ≥ 5 returns None.
220 /// assert!(Matrix::<5>::identity().det_direct().is_none());
221 /// ```
222 #[inline]
223 #[must_use]
224 pub const fn det_direct(&self) -> Option<f64> {
225 match D {
226 0 => Some(1.0),
227 1 => Some(self.rows[0][0]),
228 2 => {
229 // ad - bc
230 Some(self.rows[0][0].mul_add(self.rows[1][1], -(self.rows[0][1] * self.rows[1][0])))
231 }
232 3 => {
233 // Cofactor expansion on first row.
234 let m00 =
235 self.rows[1][1].mul_add(self.rows[2][2], -(self.rows[1][2] * self.rows[2][1]));
236 let m01 =
237 self.rows[1][0].mul_add(self.rows[2][2], -(self.rows[1][2] * self.rows[2][0]));
238 let m02 =
239 self.rows[1][0].mul_add(self.rows[2][1], -(self.rows[1][1] * self.rows[2][0]));
240 Some(
241 self.rows[0][0]
242 .mul_add(m00, (-self.rows[0][1]).mul_add(m01, self.rows[0][2] * m02)),
243 )
244 }
245 4 => {
246 // Cofactor expansion on first row → four 3×3 sub-determinants.
247 // Hoist the 6 unique 2×2 minors from rows 2–3 (each used twice).
248 let r = &self.rows;
249
250 // 2×2 minors: s_ij = r[2][i]*r[3][j] - r[2][j]*r[3][i]
251 let s23 = r[2][2].mul_add(r[3][3], -(r[2][3] * r[3][2])); // cols 2,3
252 let s13 = r[2][1].mul_add(r[3][3], -(r[2][3] * r[3][1])); // cols 1,3
253 let s12 = r[2][1].mul_add(r[3][2], -(r[2][2] * r[3][1])); // cols 1,2
254 let s03 = r[2][0].mul_add(r[3][3], -(r[2][3] * r[3][0])); // cols 0,3
255 let s02 = r[2][0].mul_add(r[3][2], -(r[2][2] * r[3][0])); // cols 0,2
256 let s01 = r[2][0].mul_add(r[3][1], -(r[2][1] * r[3][0])); // cols 0,1
257
258 // 3×3 cofactors via row 1 expansion using hoisted minors.
259 let c00 = r[1][1].mul_add(s23, (-r[1][2]).mul_add(s13, r[1][3] * s12));
260 let c01 = r[1][0].mul_add(s23, (-r[1][2]).mul_add(s03, r[1][3] * s02));
261 let c02 = r[1][0].mul_add(s13, (-r[1][1]).mul_add(s03, r[1][3] * s01));
262 let c03 = r[1][0].mul_add(s12, (-r[1][1]).mul_add(s02, r[1][2] * s01));
263
264 Some(r[0][0].mul_add(
265 c00,
266 (-r[0][1]).mul_add(c01, r[0][2].mul_add(c02, -(r[0][3] * c03))),
267 ))
268 }
269 _ => None,
270 }
271 }
272
273 /// Determinant, using closed-form formulas for D ≤ 4 and LU decomposition for D ≥ 5.
274 ///
275 /// For D ∈ {1, 2, 3, 4}, this bypasses LU factorization entirely for a significant
276 /// speedup (see [`det_direct`](Self::det_direct)). The `tol` parameter is only used
277 /// by the LU fallback path for D ≥ 5.
278 ///
279 /// # Examples
280 /// ```
281 /// use la_stack::prelude::*;
282 ///
283 /// # fn main() -> Result<(), LaError> {
284 /// let det = Matrix::<3>::identity().det(DEFAULT_PIVOT_TOL)?;
285 /// assert!((det - 1.0).abs() <= 1e-12);
286 /// # Ok(())
287 /// # }
288 /// ```
289 ///
290 /// # Errors
291 /// Returns [`LaError::NonFinite`] if the result contains NaN or infinity.
292 /// For D ≥ 5, propagates LU factorization errors (e.g. [`LaError::Singular`]).
293 #[inline]
294 pub fn det(self, tol: f64) -> Result<f64, LaError> {
295 if let Some(d) = self.det_direct() {
296 return if d.is_finite() {
297 Ok(d)
298 } else {
299 // Scan for the first non-finite entry to preserve coordinates.
300 for r in 0..D {
301 for c in 0..D {
302 if !self.rows[r][c].is_finite() {
303 return Err(LaError::NonFinite {
304 row: Some(r),
305 col: c,
306 });
307 }
308 }
309 }
310 // All entries are finite but the determinant overflowed.
311 Err(LaError::NonFinite { row: None, col: 0 })
312 };
313 }
314 self.lu(tol).map(|lu| lu.det())
315 }
316
317 /// Conservative absolute error bound for `det_direct()`.
318 ///
319 /// Returns `Some(bound)` such that `|det_direct() - det_exact| ≤ bound`,
320 /// or `None` for D ≥ 5 where no fast bound is available.
321 ///
322 /// For D ≤ 4, the bound is derived from the matrix permanent using
323 /// Shewchuk-style error analysis. For D = 0 or 1, returns `Some(0.0)`
324 /// since the determinant computation is exact (no arithmetic).
325 ///
326 /// This method does NOT require the `exact` feature — the bounds use
327 /// pure f64 arithmetic and are useful for custom adaptive-precision logic.
328 ///
329 /// # When to use
330 ///
331 /// Use this to build adaptive-precision logic: if `|det_direct()| > bound`,
332 /// the f64 sign is provably correct. Otherwise fall back to exact arithmetic.
333 ///
334 /// # Examples
335 /// ```
336 /// use la_stack::prelude::*;
337 ///
338 /// let m = Matrix::<3>::from_rows([
339 /// [1.0, 2.0, 3.0],
340 /// [4.0, 5.0, 6.0],
341 /// [7.0, 8.0, 9.0],
342 /// ]);
343 /// let bound = m.det_errbound().unwrap();
344 /// let det_approx = m.det_direct().unwrap();
345 /// // If |det_approx| > bound, the sign is guaranteed correct.
346 /// ```
347 ///
348 /// # Adaptive precision pattern (requires `exact` feature)
349 /// ```ignore
350 /// use la_stack::prelude::*;
351 ///
352 /// let m = Matrix::<3>::identity();
353 /// if let Some(bound) = m.det_errbound() {
354 /// let det = m.det_direct().unwrap();
355 /// if det.abs() > bound {
356 /// // f64 sign is guaranteed correct
357 /// let sign = det.signum() as i8;
358 /// } else {
359 /// // Fall back to exact arithmetic (requires `exact` feature)
360 /// let sign = m.det_sign_exact().unwrap();
361 /// }
362 /// } else {
363 /// // D ≥ 5: no fast filter, use exact directly
364 /// let sign = m.det_sign_exact().unwrap();
365 /// }
366 /// ```
367 #[must_use]
368 #[inline]
369 pub fn det_errbound(&self) -> Option<f64> {
370 match D {
371 0 | 1 => Some(0.0), // No arithmetic — result is exact.
372 2 => {
373 let r = &self.rows;
374 let permanent = (r[0][0] * r[1][1]).abs() + (r[0][1] * r[1][0]).abs();
375 Some(crate::ERR_COEFF_2 * permanent)
376 }
377 3 => {
378 let r = &self.rows;
379 let pm00 = (r[1][1] * r[2][2]).abs() + (r[1][2] * r[2][1]).abs();
380 let pm01 = (r[1][0] * r[2][2]).abs() + (r[1][2] * r[2][0]).abs();
381 let pm02 = (r[1][0] * r[2][1]).abs() + (r[1][1] * r[2][0]).abs();
382 let permanent = r[0][2]
383 .abs()
384 .mul_add(pm02, r[0][1].abs().mul_add(pm01, r[0][0].abs() * pm00));
385 Some(crate::ERR_COEFF_3 * permanent)
386 }
387 4 => {
388 let r = &self.rows;
389 // 2×2 minor permanents from rows 2–3.
390 let sp23 = (r[2][2] * r[3][3]).abs() + (r[2][3] * r[3][2]).abs();
391 let sp13 = (r[2][1] * r[3][3]).abs() + (r[2][3] * r[3][1]).abs();
392 let sp12 = (r[2][1] * r[3][2]).abs() + (r[2][2] * r[3][1]).abs();
393 let sp03 = (r[2][0] * r[3][3]).abs() + (r[2][3] * r[3][0]).abs();
394 let sp02 = (r[2][0] * r[3][2]).abs() + (r[2][2] * r[3][0]).abs();
395 let sp01 = (r[2][0] * r[3][1]).abs() + (r[2][1] * r[3][0]).abs();
396 // 3×3 cofactor permanents from row 1.
397 let pc0 = r[1][3]
398 .abs()
399 .mul_add(sp12, r[1][2].abs().mul_add(sp13, r[1][1].abs() * sp23));
400 let pc1 = r[1][3]
401 .abs()
402 .mul_add(sp02, r[1][2].abs().mul_add(sp03, r[1][0].abs() * sp23));
403 let pc2 = r[1][3]
404 .abs()
405 .mul_add(sp01, r[1][1].abs().mul_add(sp03, r[1][0].abs() * sp13));
406 let pc3 = r[1][2]
407 .abs()
408 .mul_add(sp01, r[1][1].abs().mul_add(sp02, r[1][0].abs() * sp12));
409 let permanent = r[0][3].abs().mul_add(
410 pc3,
411 r[0][2]
412 .abs()
413 .mul_add(pc2, r[0][1].abs().mul_add(pc1, r[0][0].abs() * pc0)),
414 );
415 Some(crate::ERR_COEFF_4 * permanent)
416 }
417 _ => None,
418 }
419 }
420}
421
422impl<const D: usize> Default for Matrix<D> {
423 #[inline]
424 fn default() -> Self {
425 Self::zero()
426 }
427}
428
429#[cfg(test)]
430mod tests {
431 use super::*;
432 use crate::DEFAULT_PIVOT_TOL;
433
434 use approx::assert_abs_diff_eq;
435 use pastey::paste;
436 use std::hint::black_box;
437
438 macro_rules! gen_public_api_matrix_tests {
439 ($d:literal) => {
440 paste! {
441 #[test]
442 fn [<public_api_matrix_from_rows_get_set_bounds_checked_ $d d>]() {
443 let mut rows = [[0.0f64; $d]; $d];
444 rows[0][0] = 1.0;
445 rows[$d - 1][$d - 1] = -2.0;
446
447 let mut m = Matrix::<$d>::from_rows(rows);
448
449 assert_eq!(m.get(0, 0), Some(1.0));
450 assert_eq!(m.get($d - 1, $d - 1), Some(-2.0));
451
452 // Out-of-bounds is None.
453 assert_eq!(m.get($d, 0), None);
454
455 // Out-of-bounds set fails.
456 assert!(!m.set($d, 0, 3.0));
457
458 // In-bounds set works.
459 assert!(m.set(0, $d - 1, 3.0));
460 assert_eq!(m.get(0, $d - 1), Some(3.0));
461 }
462
463 #[test]
464 fn [<public_api_matrix_zero_and_default_are_zero_ $d d>]() {
465 let z = Matrix::<$d>::zero();
466 assert_abs_diff_eq!(z.inf_norm(), 0.0, epsilon = 0.0);
467
468 let d = Matrix::<$d>::default();
469 assert_abs_diff_eq!(d.inf_norm(), 0.0, epsilon = 0.0);
470 }
471
472 #[test]
473 fn [<public_api_matrix_inf_norm_max_row_sum_ $d d>]() {
474 let mut rows = [[0.0f64; $d]; $d];
475
476 // Row 0 has absolute row sum = D.
477 for c in 0..$d {
478 rows[0][c] = -1.0;
479 }
480
481 // Row 1 has smaller absolute row sum.
482 for c in 0..$d {
483 rows[1][c] = 0.5;
484 }
485
486 let m = Matrix::<$d>::from_rows(rows);
487 assert_abs_diff_eq!(m.inf_norm(), f64::from($d), epsilon = 0.0);
488 }
489
490 #[test]
491 fn [<public_api_matrix_identity_lu_det_solve_vec_ $d d>]() {
492 let m = Matrix::<$d>::identity();
493
494 // Identity has ones on diag and zeros off diag.
495 for r in 0..$d {
496 for c in 0..$d {
497 let expected = if r == c { 1.0 } else { 0.0 };
498 assert_abs_diff_eq!(m.get(r, c).unwrap(), expected, epsilon = 0.0);
499 }
500 }
501
502 // Determinant is 1.
503 let det = m.det(DEFAULT_PIVOT_TOL).unwrap();
504 assert_abs_diff_eq!(det, 1.0, epsilon = 1e-12);
505
506 // LU solve on identity returns the RHS.
507 let lu = m.lu(DEFAULT_PIVOT_TOL).unwrap();
508
509 let b_arr = {
510 let mut arr = [0.0f64; $d];
511 let values = [1.0f64, 2.0, 3.0, 4.0, 5.0];
512 for (dst, src) in arr.iter_mut().zip(values.iter()) {
513 *dst = *src;
514 }
515 arr
516 };
517
518 let b = crate::Vector::<$d>::new(b_arr);
519 let x = lu.solve_vec(b).unwrap().into_array();
520
521 for (x_i, b_i) in x.iter().zip(b_arr.iter()) {
522 assert_abs_diff_eq!(*x_i, *b_i, epsilon = 1e-12);
523 }
524 }
525 }
526 };
527 }
528
529 // Mirror delaunay-style multi-dimension tests.
530 gen_public_api_matrix_tests!(2);
531 gen_public_api_matrix_tests!(3);
532 gen_public_api_matrix_tests!(4);
533 gen_public_api_matrix_tests!(5);
534
535 // === det_direct tests ===
536
537 #[test]
538 fn det_direct_d0_is_one() {
539 assert_eq!(Matrix::<0>::zero().det_direct(), Some(1.0));
540 }
541
542 #[test]
543 fn det_direct_d1_returns_element() {
544 let m = Matrix::<1>::from_rows([[42.0]]);
545 assert_eq!(m.det_direct(), Some(42.0));
546 }
547
548 #[test]
549 fn det_direct_d2_known_value() {
550 // [[1,2],[3,4]] → det = 1*4 - 2*3 = -2
551 // black_box prevents compile-time constant folding of the const fn.
552 let m = black_box(Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]));
553 assert_abs_diff_eq!(m.det_direct().unwrap(), -2.0, epsilon = 1e-15);
554 }
555
556 #[test]
557 fn det_direct_d3_known_value() {
558 // Classic 3×3: det = 0
559 let m = black_box(Matrix::<3>::from_rows([
560 [1.0, 2.0, 3.0],
561 [4.0, 5.0, 6.0],
562 [7.0, 8.0, 9.0],
563 ]));
564 assert_abs_diff_eq!(m.det_direct().unwrap(), 0.0, epsilon = 1e-12);
565 }
566
567 #[test]
568 fn det_direct_d3_nonsingular() {
569 // [[2,1,0],[0,3,1],[1,0,2]] → det = 2*(6-0) - 1*(0-1) + 0 = 13
570 let m = black_box(Matrix::<3>::from_rows([
571 [2.0, 1.0, 0.0],
572 [0.0, 3.0, 1.0],
573 [1.0, 0.0, 2.0],
574 ]));
575 assert_abs_diff_eq!(m.det_direct().unwrap(), 13.0, epsilon = 1e-12);
576 }
577
578 #[test]
579 fn det_direct_d4_identity() {
580 let m = black_box(Matrix::<4>::identity());
581 assert_abs_diff_eq!(m.det_direct().unwrap(), 1.0, epsilon = 1e-15);
582 }
583
584 #[test]
585 fn det_direct_d4_known_value() {
586 // Diagonal matrix: det = product of diagonal entries.
587 let mut rows = [[0.0f64; 4]; 4];
588 rows[0][0] = 2.0;
589 rows[1][1] = 3.0;
590 rows[2][2] = 5.0;
591 rows[3][3] = 7.0;
592 let m = black_box(Matrix::<4>::from_rows(rows));
593 assert_abs_diff_eq!(m.det_direct().unwrap(), 210.0, epsilon = 1e-12);
594 }
595
596 #[test]
597 fn det_direct_d5_returns_none() {
598 assert_eq!(Matrix::<5>::identity().det_direct(), None);
599 }
600
601 #[test]
602 fn det_direct_d8_returns_none() {
603 assert_eq!(Matrix::<8>::zero().det_direct(), None);
604 }
605
606 macro_rules! gen_det_direct_agrees_with_lu {
607 ($d:literal) => {
608 paste! {
609 #[test]
610 #[allow(clippy::cast_precision_loss)] // r, c, D are tiny integers
611 fn [<det_direct_agrees_with_lu_ $d d>]() {
612 // Well-conditioned matrix: diagonally dominant.
613 let mut rows = [[0.0f64; $d]; $d];
614 for r in 0..$d {
615 for c in 0..$d {
616 rows[r][c] = if r == c {
617 (r as f64) + f64::from($d) + 1.0
618 } else {
619 0.1 / ((r + c + 1) as f64)
620 };
621 }
622 }
623 let m = Matrix::<$d>::from_rows(rows);
624 let direct = m.det_direct().unwrap();
625 let lu_det = m.lu(DEFAULT_PIVOT_TOL).unwrap().det();
626 let eps = lu_det.abs().mul_add(1e-12, 1e-12);
627 assert_abs_diff_eq!(direct, lu_det, epsilon = eps);
628 }
629 }
630 };
631 }
632
633 gen_det_direct_agrees_with_lu!(1);
634 gen_det_direct_agrees_with_lu!(2);
635 gen_det_direct_agrees_with_lu!(3);
636 gen_det_direct_agrees_with_lu!(4);
637
638 #[test]
639 fn det_direct_identity_all_dims() {
640 assert_abs_diff_eq!(
641 Matrix::<1>::identity().det_direct().unwrap(),
642 1.0,
643 epsilon = 0.0
644 );
645 assert_abs_diff_eq!(
646 Matrix::<2>::identity().det_direct().unwrap(),
647 1.0,
648 epsilon = 0.0
649 );
650 assert_abs_diff_eq!(
651 Matrix::<3>::identity().det_direct().unwrap(),
652 1.0,
653 epsilon = 0.0
654 );
655 assert_abs_diff_eq!(
656 Matrix::<4>::identity().det_direct().unwrap(),
657 1.0,
658 epsilon = 0.0
659 );
660 }
661
662 #[test]
663 fn det_direct_zero_matrix() {
664 assert_abs_diff_eq!(
665 Matrix::<2>::zero().det_direct().unwrap(),
666 0.0,
667 epsilon = 0.0
668 );
669 assert_abs_diff_eq!(
670 Matrix::<3>::zero().det_direct().unwrap(),
671 0.0,
672 epsilon = 0.0
673 );
674 assert_abs_diff_eq!(
675 Matrix::<4>::zero().det_direct().unwrap(),
676 0.0,
677 epsilon = 0.0
678 );
679 }
680
681 #[test]
682 fn det_returns_nonfinite_error_for_nan_d2() {
683 let m = Matrix::<2>::from_rows([[f64::NAN, 1.0], [1.0, 1.0]]);
684 assert_eq!(
685 m.det(DEFAULT_PIVOT_TOL),
686 Err(LaError::NonFinite {
687 row: Some(0),
688 col: 0
689 })
690 );
691 }
692
693 #[test]
694 fn det_returns_nonfinite_error_for_inf_d3() {
695 let m =
696 Matrix::<3>::from_rows([[f64::INFINITY, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]);
697 assert_eq!(
698 m.det(DEFAULT_PIVOT_TOL),
699 Err(LaError::NonFinite {
700 row: Some(0),
701 col: 0
702 })
703 );
704 }
705
706 #[test]
707 fn det_direct_is_const_evaluable_d2() {
708 // Const evaluation proves the function is truly const fn.
709 const DET: Option<f64> = {
710 let m = Matrix::<2>::from_rows([[1.0, 0.0], [0.0, 1.0]]);
711 m.det_direct()
712 };
713 assert_eq!(DET, Some(1.0));
714 }
715
716 #[test]
717 fn det_direct_is_const_evaluable_d3() {
718 const DET: Option<f64> = {
719 let m = Matrix::<3>::from_rows([[2.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 5.0]]);
720 m.det_direct()
721 };
722 assert_eq!(DET, Some(30.0));
723 }
724
725 // === det_errbound tests (no `exact` feature required) ===
726
727 #[test]
728 fn det_errbound_available_without_exact_feature() {
729 // Verify det_errbound is accessible without exact feature
730 let m = Matrix::<3>::identity();
731 let bound = m.det_errbound();
732 assert!(bound.is_some());
733 assert!(bound.unwrap() > 0.0);
734 }
735
736 #[test]
737 fn det_errbound_d5_returns_none() {
738 // D=5 has no fast filter
739 assert_eq!(Matrix::<5>::identity().det_errbound(), None);
740 }
741}