docs: add det_direct example, hoist 4×4 minors, update docs

- Add examples/const_det_4x4.rs showing compile-time 4×4 determinant
- Hoist 6 unique 2×2 minors in D=4 path (was computing 12, 6 duplicated)
- Fix det_direct doc to include D=0 (empty product → 1.0) with doctest
- Add const fn design goal and compile-time determinants section to README
- Update AGENTS.md with det_direct and new example
- Add const_det_4x4 to justfile examples recipe

Author: acgetchell

AGENTS.md

@@ -28,7 +28,7 @@ When making changes in this repo, prioritize (in order):
 - Pre-commit validation: `just ci`
 - Python tests: `just test-python`
 - Run a single test (by name filter): `cargo test solve_2x2_basic` (or the full path: `cargo test lu::tests::solve_2x2_basic`)
-- Run examples: `just examples` (or `cargo run --example det_5x5` / `cargo run --example solve_5x5`)
+- Run examples: `just examples` (or `cargo run --example det_5x5` / `cargo run --example solve_5x5` / `cargo run --example const_det_4x4`)
 - Spell check: `just spell-check` (uses `typos.toml` at repo root; add false positives to `[default.extend-words]`)
 
 ## Code structure (big picture)
@@ -37,7 +37,7 @@ When making changes in this repo, prioritize (in order):
 - The linear algebra implementation is split across:
   - `src/lib.rs`: crate root + shared items (`LaError`, `DEFAULT_SINGULAR_TOL`, `DEFAULT_PIVOT_TOL`) + re-exports
   - `src/vector.rs`: `Vector<const D: usize>` (`[f64; D]`)
-  - `src/matrix.rs`: `Matrix<const D: usize>` (`[[f64; D]; D]`) + helpers (`get`, `set`, `inf_norm`, `det`)
+  - `src/matrix.rs`: `Matrix<const D: usize>` (`[[f64; D]; D]`) + helpers (`get`, `set`, `inf_norm`, `det`, `det_direct`)
   - `src/lu.rs`: `Lu<const D: usize>` factorization with partial pivoting (`solve_vec`, `det`)
   - `src/ldlt.rs`: `Ldlt<const D: usize>` factorization without pivoting for symmetric SPD/PSD matrices (`solve_vec`, `det`)
 - A minimal `justfile` exists for common workflows (see `just --list`).

README.md

@@ -29,6 +29,7 @@ while keeping the API intentionally small and explicit.
 
 - ✅ `Copy` types where possible
 - ✅ Const-generic dimensions (no dynamic sizes)
+- ✅ `const fn` where possible (compile-time evaluation of determinants, dot products, etc.)
 - ✅ Explicit algorithms (LU, solve, determinant)
 - ✅ No runtime dependencies (dev-dependencies are for contributors only)
 - ✅ Stack storage only (no heap allocation in core types)
@@ -103,12 +104,36 @@ let det = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap().det();
 assert!((det - 1.0).abs() <= 1e-12);
 ```
 
+## ⚡ Compile-time determinants (D ≤ 4)
+
+`det_direct()` is a `const fn` that computes closed-form determinants for D=1–4
+using fused multiply-add, bypassing LU factorization entirely. This enables
+compile-time evaluation when inputs are known at compile time:
+
+```rust
+use la_stack::prelude::*;
+
+// Evaluated entirely at compile time — no runtime cost.
+const DET: Option<f64> = {
+    let m = Matrix::<3>::from_rows([
+        [2.0, 0.0, 0.0],
+        [0.0, 3.0, 0.0],
+        [0.0, 0.0, 5.0],
+    ]);
+    m.det_direct()
+};
+assert_eq!(DET, Some(30.0));
+```
+
+The public `det()` method automatically dispatches through the closed-form path
+for D ≤ 4 and falls back to LU for D ≥ 5 — no API change needed.
+
 ## 🧩 API at a glance
 
 | Type | Storage | Purpose | Key methods |
 |---|---|---|---|
 | `Vector<D>` | `[f64; D]` | Fixed-length vector | `new`, `zero`, `dot`, `norm2_sq` |
-| `Matrix<D>` | `[[f64; D]; D]` | Fixed-size square matrix | `from_rows`, `zero`, `identity`, `lu`, `ldlt`, `det` |
+| `Matrix<D>` | `[[f64; D]; D]` | Fixed-size square matrix | `from_rows`, `zero`, `identity`, `lu`, `ldlt`, `det`, `det_direct` |
 | `Lu<D>` | `Matrix<D>` + pivot array | Factorization for solves/det | `solve_vec`, `det` |
 | `Ldlt<D>` | `Matrix<D>` | Factorization for symmetric SPD/PSD solves/det | `solve_vec`, `det` |
 
@@ -123,6 +148,7 @@ just examples
 # or:
 cargo run --example solve_5x5
 cargo run --example det_5x5
+cargo run --example const_det_4x4
 ```
 
 ## 🤝 Contributing

examples/const_det_4x4.rs

@@ -0,0 +1,36 @@
+//! Compile-time 4×4 determinant via `det_direct()`.
+//!
+//! Because `det_direct` is a `const fn` (Rust 1.94+), the determinant is
+//! evaluated entirely at compile time — zero runtime cost.
+
+use la_stack::prelude::*;
+
+/// A 4×4 Hilbert-like matrix with exact rational entries scaled to integers.
+const MAT: Matrix<4> = Matrix::<4>::from_rows([
+    [1.0, 2.0, 3.0, 4.0],
+    [5.0, 6.0, 7.0, 8.0],
+    [2.0, 6.0, 1.0, 5.0],
+    [3.0, 8.0, 2.0, 9.0],
+]);
+
+/// Determinant computed at compile time.
+const DET: f64 = match MAT.det_direct() {
+    Some(d) => d,
+    None => panic!("det_direct only supports D <= 4"),
+};
+
+fn main() {
+    println!("4×4 matrix:");
+    for r in 0..4 {
+        print!("  [");
+        for c in 0..4 {
+            if c > 0 {
+                print!(", ");
+            }
+            print!("{:5.1}", MAT.get(r, c).unwrap());
+        }
+        println!("]");
+    }
+    println!();
+    println!("det (computed at compile time) = {DET}");
+}

justfile

@@ -187,6 +187,7 @@ doc-check:
 examples:
     cargo run --quiet --example det_5x5
     cargo run --quiet --example solve_5x5
+    cargo run --quiet --example const_det_4x4
 
 # Fix (mutating): apply formatters/auto-fixes
 fix: toml-fmt fmt python-fix shell-fmt markdown-fix yaml-fix

src/matrix.rs

@@ -196,9 +196,10 @@ impl<const D: usize> Matrix<D> {
         Ldlt::factor(self, tol)
     }
 
-    /// Closed-form determinant for dimensions 1–4, bypassing LU factorization.
+    /// Closed-form determinant for dimensions 0–4, bypassing LU factorization.
     ///
-    /// Returns `Some(det)` for `D` ∈ {1, 2, 3, 4}, `None` for larger matrices.
+    /// Returns `Some(det)` for `D` ∈ {0, 1, 2, 3, 4}, `None` for D ≥ 5.
+    /// `D = 0` returns `Some(1.0)` (empty product).
     /// This is a `const fn` (Rust 1.94+) and uses fused multiply-add (`mul_add`)
     /// for improved accuracy and performance.
     ///
@@ -212,6 +213,9 @@ impl<const D: usize> Matrix<D> {
     /// let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
     /// assert!((m.det_direct().unwrap() - (-2.0)).abs() <= 1e-12);
     ///
+    /// // D = 0 is the empty product.
+    /// assert_eq!(Matrix::<0>::zero().det_direct(), Some(1.0));
+    ///
     /// // D ≥ 5 returns None.
     /// assert!(Matrix::<5>::identity().det_direct().is_none());
     /// ```
@@ -239,33 +243,23 @@ impl<const D: usize> Matrix<D> {
                 )
             }
             4 => {
-                // Cofactor expansion on first row → four 3×3 sub-determinants,
-                // each computed inline (closures are not const-compatible).
+                // Cofactor expansion on first row → four 3×3 sub-determinants.
+                // Hoist the 6 unique 2×2 minors from rows 2–3 (each used twice).
                 let r = &self.rows;
 
-                // Minor M00: rows 1-3, cols 1-3
-                let m00_0 = r[2][2].mul_add(r[3][3], -(r[2][3] * r[3][2]));
-                let m00_1 = r[2][1].mul_add(r[3][3], -(r[2][3] * r[3][1]));
-                let m00_2 = r[2][1].mul_add(r[3][2], -(r[2][2] * r[3][1]));
-                let c00 = r[1][1].mul_add(m00_0, (-r[1][2]).mul_add(m00_1, r[1][3] * m00_2));
-
-                // Minor M01: rows 1-3, cols 0,2,3
-                let m01_0 = r[2][2].mul_add(r[3][3], -(r[2][3] * r[3][2]));
-                let m01_1 = r[2][0].mul_add(r[3][3], -(r[2][3] * r[3][0]));
-                let m01_2 = r[2][0].mul_add(r[3][2], -(r[2][2] * r[3][0]));
-                let c01 = r[1][0].mul_add(m01_0, (-r[1][2]).mul_add(m01_1, r[1][3] * m01_2));
-
-                // Minor M02: rows 1-3, cols 0,1,3
-                let m02_0 = r[2][1].mul_add(r[3][3], -(r[2][3] * r[3][1]));
-                let m02_1 = r[2][0].mul_add(r[3][3], -(r[2][3] * r[3][0]));
-                let m02_2 = r[2][0].mul_add(r[3][1], -(r[2][1] * r[3][0]));
-                let c02 = r[1][0].mul_add(m02_0, (-r[1][1]).mul_add(m02_1, r[1][3] * m02_2));
-
-                // Minor M03: rows 1-3, cols 0,1,2
-                let m03_0 = r[2][1].mul_add(r[3][2], -(r[2][2] * r[3][1]));
-                let m03_1 = r[2][0].mul_add(r[3][2], -(r[2][2] * r[3][0]));
-                let m03_2 = r[2][0].mul_add(r[3][1], -(r[2][1] * r[3][0]));
-                let c03 = r[1][0].mul_add(m03_0, (-r[1][1]).mul_add(m03_1, r[1][2] * m03_2));
+                // 2×2 minors: s_ij = r[2][i]*r[3][j] - r[2][j]*r[3][i]
+                let s23 = r[2][2].mul_add(r[3][3], -(r[2][3] * r[3][2])); // cols 2,3
+                let s13 = r[2][1].mul_add(r[3][3], -(r[2][3] * r[3][1])); // cols 1,3
+                let s12 = r[2][1].mul_add(r[3][2], -(r[2][2] * r[3][1])); // cols 1,2
+                let s03 = r[2][0].mul_add(r[3][3], -(r[2][3] * r[3][0])); // cols 0,3
+                let s02 = r[2][0].mul_add(r[3][2], -(r[2][2] * r[3][0])); // cols 0,2
+                let s01 = r[2][0].mul_add(r[3][1], -(r[2][1] * r[3][0])); // cols 0,1
+
+                // 3×3 cofactors via row 1 expansion using hoisted minors.
+                let c00 = r[1][1].mul_add(s23, (-r[1][2]).mul_add(s13, r[1][3] * s12));
+                let c01 = r[1][0].mul_add(s23, (-r[1][2]).mul_add(s03, r[1][3] * s02));
+                let c02 = r[1][0].mul_add(s13, (-r[1][1]).mul_add(s03, r[1][3] * s01));
+                let c03 = r[1][0].mul_add(s12, (-r[1][1]).mul_add(s02, r[1][2] * s01));
 
                 Some(r[0][0].mul_add(
                     c00,