alkahest-cas 3.12.0

//! Phase 15 — Symbolic matrices and vectors.
//!
//! Provides a dense `Matrix` of `ExprId` values together with:
//! - arithmetic (`+`, `-`, `*`)
//! - `transpose()`
//! - `det()` (Bareiss fraction-free elimination when every entry is numeric,
//!   cofactor expansion when any entry is symbolic)
//! - `jacobian(f_vec, x_vec, pool)` — the `m×n` matrix `∂f_i/∂x_j`

use crate::diff::diff;
use crate::kernel::{ExprId, ExprPool};
use crate::simplify::engine::simplify;
use std::fmt;

pub mod eigen;
pub mod exp_gate;
pub mod linear_algebra;
pub mod normal_form;
pub(crate) mod poly_det;
pub(crate) mod putzer;
mod smith;
mod smith_poly;
pub mod spectrum;
pub(crate) mod zero_test;

pub use eigen::{
    characteristic_polynomial_lambda_minus_m, diagonalize, eigenvalues, eigenvectors, EigenError,
};
pub use exp_gate::{
    take_matrix_exp_refusal, take_matrix_exp_side_conditions, ExpRefusalReason, MatrixExpRefusal,
};
pub use linear_algebra::{
    cholesky, column_space_basis, jordan_form, lu_decomposition, matrix_exponential,
    matrix_inverse, minimal_polynomial, nullspace_basis, qr_decomposition, rank,
    rational_canonical_form, row_space_basis, rref, take_matrix_inverse_side_conditions,
    LinearAlgebraError, LuDecomposition, QrDecomposition,
};
pub use normal_form::{
    hermite_basis, hermite_form, hermite_form_poly, smith_form, smith_form_poly, smith_invariants,
    IntegerMatrix, NormalFormError, PolyMatrixQ, RatUniPoly,
};
pub use spectrum::{take_spectrum_refusal, SpectrumRefusal};
pub use zero_test::{take_zero_test_refusal, ZeroTestRefusal};

// ---------------------------------------------------------------------------
// Matrix type
// ---------------------------------------------------------------------------

/// A dense symbolic matrix stored in row-major order.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Matrix {
    /// Row-major flat storage of `ExprId` entries.
    data: Vec<ExprId>,
    pub rows: usize,
    pub cols: usize,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub enum MatrixError {
    DimensionMismatch {
        msg: String,
    },
    NotSquare,
    /// The matrix cannot be inverted: its determinant is zero, or could not be
    /// shown to be non-zero.
    ///
    /// Inverting via the adjugate divides by `det`, so "probably invertible" is
    /// not good enough — if the determinant is in fact zero the returned
    /// inverse is meaningless and nothing downstream can tell. Over a
    /// transcendental extension the vanishing question is not always decidable
    /// (see the `matrix::zero_test` module), so the two cases are both refusals
    /// and this variant states their disjunction rather than picking one.
    ///
    /// Which of the two it was is available from
    /// [`take_zero_test_refusal`]: `Some(..)` means the determinant was
    /// undecided (code `E-MAT-004`), `None` means it was proven zero (code
    /// `E-MAT-003`).
    SingularMatrix,
}

impl fmt::Display for MatrixError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            MatrixError::DimensionMismatch { msg } => write!(f, "dimension mismatch: {msg}"),
            MatrixError::NotSquare => write!(f, "matrix is not square"),
            MatrixError::SingularMatrix => write!(
                f,
                "matrix is singular, or its determinant could not be proven non-zero"
            ),
        }
    }
}

impl std::error::Error for MatrixError {}

impl crate::errors::AlkahestError for MatrixError {
    fn code(&self) -> &'static str {
        match self {
            MatrixError::DimensionMismatch { .. } => "E-MAT-001",
            MatrixError::NotSquare => "E-MAT-002",
            // `E-MAT-004` — the undecided-determinant refusal — is carried by
            // this same variant and read back through
            // [`take_zero_test_refusal`]; see [`MatrixError::SingularMatrix`].
            MatrixError::SingularMatrix => "E-MAT-003",
        }
    }

    fn remediation(&self) -> Option<&'static str> {
        match self {
            MatrixError::DimensionMismatch { .. } => Some(
                "ensure all rows have the same column count and operand dimensions match",
            ),
            MatrixError::NotSquare => Some(
                "determinant and inverse require a square matrix; use the pseudo-inverse for rectangular matrices",
            ),
            MatrixError::SingularMatrix => Some(
                "check your system of equations for linear dependence; if the entries are symbolic, substitute concrete values for the parameters",
            ),
        }
    }
}

impl Matrix {
    /// Create a matrix from row-major nested vectors.
    pub fn new(rows: Vec<Vec<ExprId>>) -> Result<Self, MatrixError> {
        if rows.is_empty() {
            return Ok(Matrix {
                data: vec![],
                rows: 0,
                cols: 0,
            });
        }
        let cols = rows[0].len();
        for r in &rows {
            if r.len() != cols {
                return Err(MatrixError::DimensionMismatch {
                    msg: format!("expected {cols} columns, got {}", r.len()),
                });
            }
        }
        let nrows = rows.len();
        let data: Vec<ExprId> = rows.into_iter().flatten().collect();
        Ok(Matrix {
            data,
            rows: nrows,
            cols,
        })
    }

    /// Create a zero matrix (all entries are `pool.integer(0)`).
    pub fn zeros(rows: usize, cols: usize, pool: &ExprPool) -> Self {
        let zero = pool.integer(0_i32);
        Matrix {
            data: vec![zero; rows * cols],
            rows,
            cols,
        }
    }

    /// Create an identity matrix.
    pub fn identity(n: usize, pool: &ExprPool) -> Self {
        let zero = pool.integer(0_i32);
        let one = pool.integer(1_i32);
        let mut data = vec![zero; n * n];
        for i in 0..n {
            data[i * n + i] = one;
        }
        Matrix {
            data,
            rows: n,
            cols: n,
        }
    }

    /// Get entry at row `r`, column `c` (0-indexed).
    pub fn get(&self, r: usize, c: usize) -> ExprId {
        self.data[r * self.cols + c]
    }

    /// Set entry at row `r`, column `c`.
    pub fn set(&mut self, r: usize, c: usize, val: ExprId) {
        self.data[r * self.cols + c] = val;
    }

    /// Get a row as a vector.
    pub fn row(&self, r: usize) -> Vec<ExprId> {
        self.data[r * self.cols..(r + 1) * self.cols].to_vec()
    }

    /// Get a column as a vector.
    pub fn col(&self, c: usize) -> Vec<ExprId> {
        (0..self.rows).map(|r| self.get(r, c)).collect()
    }

    /// Transpose.
    pub fn transpose(&self) -> Self {
        let mut data = Vec::with_capacity(self.rows * self.cols);
        for c in 0..self.cols {
            for r in 0..self.rows {
                data.push(self.get(r, c));
            }
        }
        Matrix {
            data,
            rows: self.cols,
            cols: self.rows,
        }
    }

    /// Element-wise addition.
    pub fn add(&self, other: &Matrix, pool: &ExprPool) -> Result<Matrix, MatrixError> {
        self.check_same_shape(other)?;
        let data = self
            .data
            .iter()
            .zip(other.data.iter())
            .map(|(&a, &b)| pool.add(vec![a, b]))
            .collect();
        Ok(Matrix {
            data,
            rows: self.rows,
            cols: self.cols,
        })
    }

    /// Element-wise subtraction.
    pub fn sub(&self, other: &Matrix, pool: &ExprPool) -> Result<Matrix, MatrixError> {
        self.check_same_shape(other)?;
        let neg_one = pool.integer(-1_i32);
        let data = self
            .data
            .iter()
            .zip(other.data.iter())
            .map(|(&a, &b)| {
                let neg_b = pool.mul(vec![neg_one, b]);
                pool.add(vec![a, neg_b])
            })
            .collect();
        Ok(Matrix {
            data,
            rows: self.rows,
            cols: self.cols,
        })
    }

    /// Matrix multiplication (`self` is m×k, `other` is k×n → result is m×n).
    pub fn mul(&self, other: &Matrix, pool: &ExprPool) -> Result<Matrix, MatrixError> {
        if self.cols != other.rows {
            return Err(MatrixError::DimensionMismatch {
                msg: format!(
                    "cannot multiply {}×{} by {}×{}",
                    self.rows, self.cols, other.rows, other.cols
                ),
            });
        }
        let m = self.rows;
        let n = other.cols;
        let k = self.cols;
        let mut data = Vec::with_capacity(m * n);
        for r in 0..m {
            for c in 0..n {
                let terms: Vec<ExprId> = (0..k)
                    .map(|i| pool.mul(vec![self.get(r, i), other.get(i, c)]))
                    .collect();
                let entry = if terms.is_empty() {
                    pool.integer(0_i32)
                } else if terms.len() == 1 {
                    terms[0]
                } else {
                    pool.add(terms)
                };
                data.push(entry);
            }
        }
        Ok(Matrix {
            data,
            rows: m,
            cols: n,
        })
    }

    /// Scalar multiplication.
    pub fn scale(&self, scalar: ExprId, pool: &ExprPool) -> Matrix {
        let data = self
            .data
            .iter()
            .map(|&e| pool.mul(vec![scalar, e]))
            .collect();
        Matrix {
            data,
            rows: self.rows,
            cols: self.cols,
        }
    }

    /// Simplify all entries.
    pub fn simplify_entries(&self, pool: &ExprPool) -> Matrix {
        let data = self.data.iter().map(|&e| simplify(e, pool).value).collect();
        Matrix {
            data,
            rows: self.rows,
            cols: self.cols,
        }
    }

    /// Every entry as an exact rational, or `None` if any entry is not a
    /// numeric literal.
    fn numeric_entries(&self, pool: &ExprPool) -> Option<Vec<rug::Rational>> {
        self.data
            .iter()
            .map(|&e| {
                pool.with(e, |d| match d {
                    crate::kernel::ExprData::Integer(i) => Some(rug::Rational::from(i.0.clone())),
                    crate::kernel::ExprData::Rational(r) => Some(r.0.clone()),
                    _ => None,
                })
            })
            .collect()
    }

    /// Determinant of a matrix whose entries are all numeric literals.
    ///
    /// `O(n³)` against the `O(n!)` of the cofactor expansion in
    /// [`det`](Matrix::det). Each row is scaled by the lcm of its
    /// denominators, the integer determinant is taken by FLINT's
    /// `fmpz_mat_det` (multimodular / fraction-free, chosen by FLINT), and the
    /// product of the row scales is divided back out. The value is exact and
    /// unique, so this is identical to the Bareiss elimination over
    /// `rug::Rational` it replaced (kept as `det_bareiss_reference` for the
    /// differential tests) — measured 0.21 ms → 0.006 ms at `n = 10` and
    /// 45 ms → 1 ms at `n = 60`, conversion included.
    fn det_numeric(&self, pool: &ExprPool) -> Option<ExprId> {
        let m = self.numeric_entries(pool)?;
        let d = det_rational_flint(&m, self.rows);
        Some(if *d.denom() == 1 {
            pool.integer(d.numer().clone())
        } else {
            pool.rational(d.numer().clone(), d.denom().clone())
        })
    }

    /// Determinant.
    ///
    /// All-numeric matrices take the `O(n³)` Bareiss route in
    /// `det_numeric` (private); symbolic entries fall back to
    /// cofactor expansion along the first row.  The expansion is memoised on
    /// the `(rows, columns)` of each sub-minor (see `DetMemo`), so it
    /// evaluates `O(n·2ⁿ)` sub-determinants rather than the `O(n!)` of the
    /// plain recursion, and returns exactly the same expression.
    pub fn det(&self, pool: &ExprPool) -> Result<ExprId, MatrixError> {
        if self.rows != self.cols {
            return Err(MatrixError::NotSquare);
        }
        let n = self.rows;
        if n == 0 {
            return Ok(pool.integer(1_i32));
        }
        if n <= DetMemo::MAX_DIM {
            let mut memo = DetMemo::new(self, pool);
            let all = DetMemo::full_mask(n);
            return Ok(memo.det(self, all, all, pool));
        }
        self.det_unmemoized(pool)
    }

    /// Determinant of the minor that removes row `skip_row` and column
    /// `skip_col`, sharing `memo` with other minors of the same matrix.
    ///
    /// Identical to `self.minor(skip_row, skip_col).det(pool)`; the adjugate
    /// asks for all `n²` minors, and with a shared memo their common
    /// sub-minors are expanded once rather than once per cofactor.
    pub(crate) fn minor_det_memo(
        &self,
        skip_row: usize,
        skip_col: usize,
        memo: &mut DetMemo,
        pool: &ExprPool,
    ) -> ExprId {
        let n = self.rows;
        debug_assert_eq!(n, self.cols);
        debug_assert!(n >= 1);
        if n - 1 == 0 {
            return pool.integer(1_i32);
        }
        if n > DetMemo::MAX_DIM {
            return self
                .minor(skip_row, skip_col)
                .det_unmemoized(pool)
                .expect("a minor of a square matrix is square");
        }
        let all = DetMemo::full_mask(n);
        memo.det(self, all & !(1 << skip_row), all & !(1 << skip_col), pool)
    }

    /// The original recursive cofactor expansion, for matrices too large for
    /// [`DetMemo`]'s bit masks (where a symbolic expansion cannot finish
    /// anyway).  [`DetMemo::det`] computes exactly the same expression.
    fn det_unmemoized(&self, pool: &ExprPool) -> Result<ExprId, MatrixError> {
        let n = self.rows;
        if n == 0 {
            return Ok(pool.integer(1_i32));
        }
        if n == 1 {
            return Ok(self.get(0, 0));
        }
        if n >= 3 {
            if let Some(d) = self.det_numeric(pool) {
                return Ok(d);
            }
        }
        if n == 2 {
            // ad - bc
            let ad = pool.mul(vec![self.get(0, 0), self.get(1, 1)]);
            let bc = pool.mul(vec![self.get(0, 1), self.get(1, 0)]);
            let neg_bc = pool.mul(vec![pool.integer(-1_i32), bc]);
            return Ok(simplify(pool.add(vec![ad, neg_bc]), pool).value);
        }
        // Cofactor expansion along first row for n >= 3
        let mut terms: Vec<ExprId> = Vec::new();
        for j in 0..n {
            let minor = self.minor(0, j);
            let minor_det = minor.det_unmemoized(pool)?;
            let sign = if j % 2 == 0 {
                pool.integer(1_i32)
            } else {
                pool.integer(-1_i32)
            };
            terms.push(pool.mul(vec![sign, self.get(0, j), minor_det]));
        }
        Ok(simplify(pool.add(terms), pool).value)
    }

    /// Submatrix obtained by removing row `r` and column `c`.
    pub(crate) fn minor(&self, skip_row: usize, skip_col: usize) -> Matrix {
        let n = self.rows;
        let mut data = Vec::with_capacity((n - 1) * (n - 1));
        for r in 0..n {
            if r == skip_row {
                continue;
            }
            for c in 0..n {
                if c == skip_col {
                    continue;
                }
                data.push(self.get(r, c));
            }
        }
        Matrix {
            data,
            rows: n - 1,
            cols: n - 1,
        }
    }

    fn check_same_shape(&self, other: &Matrix) -> Result<(), MatrixError> {
        if self.rows != other.rows || self.cols != other.cols {
            Err(MatrixError::DimensionMismatch {
                msg: format!(
                    "{}×{} vs {}×{}",
                    self.rows, self.cols, other.rows, other.cols
                ),
            })
        } else {
            Ok(())
        }
    }

    /// Return a flat reference to all entries.
    pub fn entries(&self) -> &[ExprId] {
        &self.data
    }

    /// Return entries as a nested `Vec<Vec<ExprId>>`.
    pub fn to_nested(&self) -> Vec<Vec<ExprId>> {
        (0..self.rows).map(|r| self.row(r)).collect()
    }

    /// V2-17 — `det(λI − M)` as a pooled expression plus the fresh λ symbol used.
    pub fn characteristic_polynomial_lambda_minus_m(
        &self,
        pool: &ExprPool,
    ) -> Result<(ExprId, ExprId), EigenError> {
        eigen::characteristic_polynomial_lambda_minus_m(self, pool)
    }

    /// V2-17 — Algebraic eigenvalues `(value, multiplicity)` for matrices whose characteristic
    /// polynomial factors over ℚ into linear and quadratic terms.
    pub fn eigenvalues(&self, pool: &ExprPool) -> Result<Vec<(ExprId, usize)>, EigenError> {
        eigen::eigenvalues(self, pool)
    }

    /// V2-17 — Eigenvalue tuples `(λ, multiplicity, column eigenvectors)`.
    pub fn eigenvectors(
        &self,
        pool: &ExprPool,
    ) -> Result<Vec<(ExprId, usize, Vec<Matrix>)>, EigenError> {
        eigen::eigenvectors(self, pool)
    }

    /// V2-17 — `(P, D)` with `M·P == P·D` when diagonalizable in the ℚ-splitting-field sense.
    pub fn diagonalize(&self, pool: &ExprPool) -> Result<(Matrix, Matrix), EigenError> {
        eigen::diagonalize(self, pool)
    }

    pub fn nullspace(&self, pool: &ExprPool) -> Result<Vec<Matrix>, LinearAlgebraError> {
        linear_algebra::nullspace_basis(self, pool)
    }

    pub fn rank(&self, pool: &ExprPool) -> Result<usize, LinearAlgebraError> {
        linear_algebra::rank(self, pool)
    }

    pub fn rref(&self, pool: &ExprPool) -> Result<Matrix, LinearAlgebraError> {
        linear_algebra::rref(self, pool)
    }

    pub fn column_space(&self, pool: &ExprPool) -> Result<Vec<Matrix>, LinearAlgebraError> {
        linear_algebra::column_space_basis(self, pool)
    }

    pub fn row_space(&self, pool: &ExprPool) -> Result<Vec<Matrix>, LinearAlgebraError> {
        linear_algebra::row_space_basis(self, pool)
    }

    pub fn lu(&self, pool: &ExprPool) -> Result<LuDecomposition, LinearAlgebraError> {
        linear_algebra::lu_decomposition(self, pool)
    }

    pub fn qr(&self, pool: &ExprPool) -> Result<QrDecomposition, LinearAlgebraError> {
        linear_algebra::qr_decomposition(self, pool)
    }

    pub fn cholesky(&self, pool: &ExprPool) -> Result<Matrix, LinearAlgebraError> {
        linear_algebra::cholesky(self, pool)
    }

    pub fn jordan_form(&self, pool: &ExprPool) -> Result<(Matrix, Matrix), LinearAlgebraError> {
        linear_algebra::jordan_form(self, pool)
    }

    pub fn rational_canonical_form(
        &self,
        pool: &ExprPool,
    ) -> Result<(Matrix, Matrix), LinearAlgebraError> {
        linear_algebra::rational_canonical_form(self, pool)
    }

    pub fn minimal_polynomial(
        &self,
        pool: &ExprPool,
    ) -> Result<(ExprId, ExprId), LinearAlgebraError> {
        linear_algebra::minimal_polynomial(self, pool)
    }

    pub fn matrix_exp(&self, pool: &ExprPool) -> Result<Matrix, LinearAlgebraError> {
        linear_algebra::matrix_exponential(self, pool)
    }

    pub fn inverse(&self, pool: &ExprPool) -> Result<Matrix, MatrixError> {
        linear_algebra::matrix_inverse(self, pool)
    }
}

// ---------------------------------------------------------------------------
// Jacobian
// ---------------------------------------------------------------------------

/// Compute the Jacobian matrix `J[i][j] = ∂f_i/∂x_j`.
///
/// `f_vec` is a slice of m scalar expressions; `x_vec` is a slice of n
/// variable expressions.  The result is an m×n `Matrix`.
pub fn jacobian(
    f_vec: &[ExprId],
    x_vec: &[ExprId],
    pool: &ExprPool,
) -> Result<Matrix, crate::diff::diff_impl::DiffError> {
    let m = f_vec.len();
    let n = x_vec.len();
    let mut data = Vec::with_capacity(m * n);
    for &f in f_vec {
        for &x in x_vec {
            let df = diff(f, x, pool)?.value;
            data.push(df);
        }
    }
    Ok(Matrix {
        data,
        rows: m,
        cols: n,
    })
}

// ---------------------------------------------------------------------------
// Display
// ---------------------------------------------------------------------------

impl Matrix {
    pub fn display(&self, pool: &ExprPool) -> String {
        let rows: Vec<String> = (0..self.rows)
            .map(|r| {
                let entries: Vec<String> = self
                    .row(r)
                    .into_iter()
                    .map(|e| pool.display(e).to_string())
                    .collect();
                format!("[{}]", entries.join(", "))
            })
            .collect();
        format!("[{}]", rows.join(", "))
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

/// Memo for the symbolic cofactor expansion in [`Matrix::det`].
///
/// Expanding along the first remaining row, every sub-determinant that the
/// recursion reaches is the determinant of the submatrix picked out by a set
/// of rows and a set of columns of the original matrix (both kept in their
/// original order).  The plain recursion re-derives each of them once per
/// path that reaches it — `n!` leaf calls for an `n×n` matrix — while there
/// are only `2ⁿ` distinct column sets at each depth.  Keying on the
/// `(rows, cols)` bit masks computes each once.
///
/// Every step is exactly the step the plain recursion takes on the same
/// submatrix (numeric Bareiss for an all-numeric block of size ≥ 3, `ad − bc`
/// for size 2, first-row expansion otherwise), and the pool is hash-consed,
/// so the result is the same `ExprId` the unmemoised expansion returns.
pub(crate) struct DetMemo {
    /// `(row mask, column mask) → det` of that submatrix.
    table: std::collections::HashMap<(u64, u64), ExprId>,
    /// Per row, the columns whose entry is *not* a numeric literal.
    symbolic_cols: Vec<u64>,
}

impl DetMemo {
    /// Largest dimension the `u64` masks can describe.
    const MAX_DIM: usize = 64;

    fn full_mask(n: usize) -> u64 {
        if n >= 64 {
            u64::MAX
        } else {
            (1u64 << n) - 1
        }
    }

    pub(crate) fn new(m: &Matrix, pool: &ExprPool) -> Self {
        let symbolic_cols = (0..m.rows)
            .map(|r| {
                let mut bits = 0u64;
                for c in 0..m.cols.min(Self::MAX_DIM) {
                    let numeric = pool.with(m.get(r, c), |d| {
                        matches!(
                            d,
                            crate::kernel::ExprData::Integer(_)
                                | crate::kernel::ExprData::Rational(_)
                        )
                    });
                    if !numeric {
                        bits |= 1 << c;
                    }
                }
                bits
            })
            .collect();
        DetMemo {
            table: std::collections::HashMap::new(),
            symbolic_cols,
        }
    }

    /// Ascending bit positions of `mask`.
    fn bits(mut mask: u64) -> impl Iterator<Item = usize> {
        std::iter::from_fn(move || {
            if mask == 0 {
                None
            } else {
                let i = mask.trailing_zeros() as usize;
                mask &= mask - 1;
                Some(i)
            }
        })
    }

    /// The submatrix of `m` on `rows × cols`, as a `Matrix`.
    fn submatrix(m: &Matrix, rows: u64, cols: u64) -> Matrix {
        let k = rows.count_ones() as usize;
        let mut data = Vec::with_capacity(k * k);
        for r in Self::bits(rows) {
            for c in Self::bits(cols) {
                data.push(m.get(r, c));
            }
        }
        Matrix {
            data,
            rows: k,
            cols: k,
        }
    }

    /// `det` of `m` restricted to `rows × cols` (equal popcounts, ≥ 1).
    fn det(&mut self, m: &Matrix, rows: u64, cols: u64, pool: &ExprPool) -> ExprId {
        let k = rows.count_ones() as usize;
        debug_assert_eq!(k, cols.count_ones() as usize);
        debug_assert!(k >= 1);
        if k == 1 {
            return m.get(
                rows.trailing_zeros() as usize,
                cols.trailing_zeros() as usize,
            );
        }
        if let Some(&d) = self.table.get(&(rows, cols)) {
            return d;
        }
        let d = self.det_uncached(m, rows, cols, pool);
        self.table.insert((rows, cols), d);
        d
    }

    fn det_uncached(&mut self, m: &Matrix, rows: u64, cols: u64, pool: &ExprPool) -> ExprId {
        let k = rows.count_ones() as usize;
        if k >= 3 && Self::bits(rows).all(|r| self.symbolic_cols[r] & cols == 0) {
            if let Some(d) = Self::submatrix(m, rows, cols).det_numeric(pool) {
                return d;
            }
        }
        let r0 = rows.trailing_zeros() as usize;
        let rest = rows & (rows - 1);
        if k == 2 {
            // ad - bc, exactly as `Matrix::det` builds it for a 2×2.
            let mut cs = Self::bits(cols);
            let (c0, c1) = (cs.next().unwrap(), cs.next().unwrap());
            let r1 = rest.trailing_zeros() as usize;
            let ad = pool.mul(vec![m.get(r0, c0), m.get(r1, c1)]);
            let bc = pool.mul(vec![m.get(r0, c1), m.get(r1, c0)]);
            let neg_bc = pool.mul(vec![pool.integer(-1_i32), bc]);
            return simplify(pool.add(vec![ad, neg_bc]), pool).value;
        }
        // Cofactor expansion along the first remaining row.
        let mut terms: Vec<ExprId> = Vec::with_capacity(k);
        for (j, c) in Self::bits(cols).enumerate() {
            let minor_det = self.det(m, rest, cols & !(1 << c), pool);
            let sign = if j % 2 == 0 {
                pool.integer(1_i32)
            } else {
                pool.integer(-1_i32)
            };
            terms.push(pool.mul(vec![sign, m.get(r0, c), minor_det]));
        }
        simplify(pool.add(terms), pool).value
    }
}

/// Exact determinant of the `n × n` row-major rational matrix `m`, through
/// `fmpz_mat_det` after clearing each row's denominators.
pub(crate) fn det_rational_flint(m: &[rug::Rational], n: usize) -> rug::Rational {
    use crate::flint::{mat::FlintMat, FlintInteger};
    debug_assert_eq!(m.len(), n * n);
    if n == 0 {
        return rug::Rational::from(1);
    }
    let mut fm = FlintMat::new(n, n);
    let mut scale = rug::Integer::from(1);
    for i in 0..n {
        let row = &m[i * n..(i + 1) * n];
        let mut l = rug::Integer::from(1);
        for c in row {
            if *c.denom() != 1 {
                l.lcm_mut(c.denom());
            }
        }
        for (j, c) in row.iter().enumerate() {
            let v = if l == 1 {
                c.numer().clone()
            } else {
                c.numer() * rug::Integer::from(&l / c.denom())
            };
            fm.set_entry(i, j, &FlintInteger::from_rug(&v));
        }
        scale *= l;
    }
    rug::Rational::from((fm.det().to_rug(), scale))
}

/// The Bareiss fraction-free elimination over `rug::Rational` that
/// [`det_rational_flint`] replaced — the reference its differential tests
/// compare against.
#[cfg(test)]
pub(crate) fn det_bareiss_reference(m: &[rug::Rational], n: usize) -> rug::Rational {
    if n == 0 {
        return rug::Rational::from(1);
    }
    let mut m = m.to_vec();
    let at = |i: usize, j: usize| i * n + j;
    let mut prev = rug::Rational::from(1);
    let mut sign = 1i32;
    for k in 0..n.saturating_sub(1) {
        if m[at(k, k)] == 0 {
            let Some(r) = (k + 1..n).find(|&r| m[at(r, k)] != 0) else {
                return rug::Rational::new();
            };
            for j in 0..n {
                m.swap(at(k, j), at(r, j));
            }
            sign = -sign;
        }
        for i in k + 1..n {
            for j in k + 1..n {
                let v = (m[at(i, j)].clone() * m[at(k, k)].clone()
                    - m[at(i, k)].clone() * m[at(k, j)].clone())
                    / prev.clone();
                m[at(i, j)] = v;
            }
        }
        prev = m[at(k, k)].clone();
    }
    let d = m[at(n - 1, n - 1)].clone();
    if sign < 0 {
        -d
    } else {
        d
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::kernel::{Domain, ExprPool};

    fn p() -> ExprPool {
        ExprPool::new()
    }

    /// The `O(n!)` cofactor expansion the numeric route replaced, kept here as
    /// the reference the fast path is checked against.
    fn det_cofactor(m: &Matrix, pool: &ExprPool) -> ExprId {
        let n = m.rows;
        if n == 1 {
            return m.get(0, 0);
        }
        let mut terms: Vec<ExprId> = Vec::new();
        for j in 0..n {
            let minor = m.minor(0, j);
            let minor_det = det_cofactor(&minor, pool);
            let sign = if j % 2 == 0 {
                pool.integer(1_i32)
            } else {
                pool.integer(-1_i32)
            };
            terms.push(pool.mul(vec![sign, m.get(0, j), minor_det]));
        }
        simplify(pool.add(terms), pool).value
    }

    /// A mixed symbolic/numeric `n×n` matrix: symbols, sums, powers, integer
    /// and rational literals, and zero rows/blocks so that some sub-minors
    /// take the all-numeric route and some are singular.
    fn mixed_matrix(n: usize, seed: u64, pool: &ExprPool) -> Matrix {
        let x = pool.symbol("x", Domain::Real);
        let y = pool.symbol("y", Domain::Real);
        let mut state = seed.wrapping_mul(0x9E37_79B9_7F4A_7C15) | 1;
        let mut rnd = |m: u64| {
            state ^= state << 13;
            state ^= state >> 7;
            state ^= state << 17;
            state % m
        };
        let rows: Vec<Vec<ExprId>> = (0..n)
            .map(|_| {
                (0..n)
                    .map(|_| match rnd(8) {
                        0 => x,
                        1 => pool.add(vec![x, pool.integer(rnd(5) as i64 - 2)]),
                        2 => pool.mul(vec![pool.integer(rnd(7) as i64 - 3), y]),
                        3 => pool.pow(x, pool.integer(rnd(3) as i64 + 2)),
                        4 => pool
                            .rational(rug::Integer::from(rnd(9) as i64 - 4), rug::Integer::from(3)),
                        5 => pool.integer(0_i32),
                        _ => pool.integer(rnd(9) as i64 - 4),
                    })
                    .collect()
            })
            .collect();
        Matrix::new(rows).expect("square")
    }

    #[test]
    fn memoised_det_is_the_same_expr_as_the_plain_expansion() {
        let pool = p();
        for n in 1..=7usize {
            for seed in 0..10u64 {
                let m = mixed_matrix(n, seed * 31 + n as u64, &pool);
                assert_eq!(
                    m.det(&pool).unwrap(),
                    m.det_unmemoized(&pool).unwrap(),
                    "n={n} seed={seed}"
                );
            }
        }
        // All-numeric blocks inside a symbolic matrix: a symbolic first row
        // over an integer block, so every sub-minor below row 0 is numeric.
        let x = pool.symbol("x", Domain::Real);
        let mut rows: Vec<Vec<ExprId>> = vec![(0..5).map(|_| x).collect()];
        for r in 1..5 {
            rows.push(
                (0..5)
                    .map(|c| pool.integer((r * 7 + c * 3) as i64 % 5))
                    .collect(),
            );
        }
        let m = Matrix::new(rows).unwrap();
        assert_eq!(m.det(&pool).unwrap(), m.det_unmemoized(&pool).unwrap());
    }

    #[test]
    fn shared_memo_minors_match_each_minor_det() {
        let pool = p();
        for n in 2..=6usize {
            for seed in 0..4u64 {
                let m = mixed_matrix(n, seed + 100 * n as u64, &pool);
                let mut memo = DetMemo::new(&m, &pool);
                for i in 0..n {
                    for j in 0..n {
                        assert_eq!(
                            m.minor_det_memo(j, i, &mut memo, &pool),
                            m.minor(j, i).det_unmemoized(&pool).unwrap(),
                            "n={n} seed={seed} minor ({j},{i})"
                        );
                    }
                }
            }
        }
    }

    /// The cofactor expansion is `O(n!)` unmemoised (a 10×10 never finished
    /// in practice: 6.8 s at `n = 9`, ×10 per step).  With the memo every
    /// sub-determinant is keyed by its `(rows, cols)` masks, and the rows
    /// are always a suffix of the matrix's rows, so at most `2ⁿ` entries are
    /// ever computed.  Asserted on the memo itself rather than on wall-clock
    /// time, which sanitizer builds stretch arbitrarily.
    #[test]
    fn symbolic_det_memo_is_at_most_two_to_the_n() {
        let pool = p();
        let x = pool.symbol("x", Domain::Real);
        let y = pool.symbol("y", Domain::Real);
        let n = 8;
        let rows: Vec<Vec<ExprId>> = (0..n)
            .map(|i| {
                (0..n)
                    .map(|j: usize| {
                        if i == j {
                            pool.add(vec![x, pool.integer(i as i64)])
                        } else {
                            pool.add(vec![
                                pool.pow(x, pool.integer(i.abs_diff(j) as i64)),
                                pool.mul(vec![pool.integer(i as i64 + 1), y]),
                            ])
                        }
                    })
                    .collect()
            })
            .collect();
        let m = Matrix::new(rows).unwrap();
        let mut memo = DetMemo::new(&m, &pool);
        let all = DetMemo::full_mask(n);
        let d = memo.det(&m, all, all, &pool);
        assert_eq!(d, m.det(&pool).unwrap());
        assert!(
            memo.table.len() < 1 << n,
            "{} memo entries for n = {n}",
            memo.table.len()
        );
    }

    #[test]
    #[ignore = "timing report, run with --ignored --nocapture"]
    fn symbolic_det_timing() {
        // Separate pools, so neither run profits from the other's interning
        // or simplifier memo.
        for n in 4..=9usize {
            let pa = p();
            let ma = mixed_matrix(n, 7, &pa);
            let t = std::time::Instant::now();
            let _ = ma.det(&pa).unwrap();
            let memo = t.elapsed();
            let pb = p();
            let mb = mixed_matrix(n, 7, &pb);
            let t = std::time::Instant::now();
            let _ = mb.det_unmemoized(&pb).unwrap();
            let plain = t.elapsed();
            println!("n={n}: memo {memo:?}  plain {plain:?}");
        }
    }

    #[test]
    fn numeric_det_agrees_with_cofactor_expansion() {
        let pool = p();
        let mut state = 0x2545_F491_4F6C_DD1D_u64;
        let mut rnd = |m: i64| {
            state = state
                .wrapping_mul(6_364_136_223_846_793_005)
                .wrapping_add(1);
            ((state >> 33) % (2 * m as u64 + 1)) as i64 - m
        };
        for n in 3..=6usize {
            for trial in 0..12 {
                // Integer entries, plus a rational-entry and a singular case.
                let data: Vec<ExprId> = (0..n * n)
                    .map(|k| match trial % 3 {
                        0 => pool.integer(rnd(9)),
                        1 => pool.rational(rug::Integer::from(rnd(9)), rug::Integer::from(7)),
                        _ => pool.integer(if k < n { 0_i32 } else { rnd(9) as i32 }),
                    })
                    .collect();
                let rows: Vec<Vec<ExprId>> = data.chunks(n).map(|c| c.to_vec()).collect();
                let m = Matrix::new(rows).expect("square");
                assert_eq!(
                    m.det(&pool).expect("square"),
                    det_cofactor(&m, &pool),
                    "n={n} trial={trial}: numeric Bareiss disagrees with cofactor expansion"
                );
            }
        }
    }

    #[test]
    fn numeric_det_handles_a_zero_pivot_and_a_singular_matrix() {
        let pool = p();
        // Leading zero pivot but nonsingular: det = -1·(0·1 − 1·1)·… → -1.
        let m = Matrix::new(vec![
            vec![
                pool.integer(0_i32),
                pool.integer(1_i32),
                pool.integer(0_i32),
            ],
            vec![
                pool.integer(1_i32),
                pool.integer(0_i32),
                pool.integer(0_i32),
            ],
            vec![
                pool.integer(0_i32),
                pool.integer(0_i32),
                pool.integer(1_i32),
            ],
        ])
        .expect("square");
        assert_eq!(m.det(&pool).unwrap(), pool.integer(-1_i32));

        // A zero row is singular.
        let z = Matrix::new(vec![
            vec![
                pool.integer(0_i32),
                pool.integer(0_i32),
                pool.integer(0_i32),
            ],
            vec![
                pool.integer(1_i32),
                pool.integer(2_i32),
                pool.integer(3_i32),
            ],
            vec![
                pool.integer(4_i32),
                pool.integer(5_i32),
                pool.integer(7_i32),
            ],
        ])
        .expect("square");
        assert_eq!(z.det(&pool).unwrap(), pool.integer(0_i32));
    }

    #[test]
    fn identity_2x2() {
        let pool = p();
        let id = Matrix::identity(2, &pool);
        assert_eq!(id.rows, 2);
        assert_eq!(id.cols, 2);
        assert_eq!(id.get(0, 0), pool.integer(1_i32));
        assert_eq!(id.get(0, 1), pool.integer(0_i32));
        assert_eq!(id.get(1, 0), pool.integer(0_i32));
        assert_eq!(id.get(1, 1), pool.integer(1_i32));
    }

    #[test]
    fn transpose_2x3() {
        let pool = p();
        let x = pool.symbol("x", Domain::Real);
        let y = pool.symbol("y", Domain::Real);
        let z = pool.symbol("z", Domain::Real);
        let a = pool.integer(1_i32);
        let b = pool.integer(2_i32);
        let c = pool.integer(3_i32);
        // [[x, y, z], [a, b, c]]  →  [[x,a],[y,b],[z,c]]
        let m = Matrix::new(vec![vec![x, y, z], vec![a, b, c]]).unwrap();
        let t = m.transpose();
        assert_eq!(t.rows, 3);
        assert_eq!(t.cols, 2);
        assert_eq!(t.get(0, 0), x);
        assert_eq!(t.get(1, 1), b);
    }

    #[test]
    fn add_matrices() {
        let pool = p();
        let x = pool.symbol("x", Domain::Real);
        let one = pool.integer(1_i32);
        let m1 = Matrix::new(vec![vec![x, one]]).unwrap();
        let m2 = Matrix::new(vec![vec![one, x]]).unwrap();
        let result = m1.add(&m2, &pool).unwrap();
        // result[0][0] = x + 1
        let r00_str = pool.display(result.get(0, 0)).to_string();
        assert!(
            r00_str.contains("x") && r00_str.contains("1"),
            "got: {r00_str}"
        );
    }

    #[test]
    fn mul_2x2() {
        let pool = p();
        // [[1,0],[0,1]] * [[a,b],[c,d]] = [[a,b],[c,d]]
        let id = Matrix::identity(2, &pool);
        let x = pool.symbol("x", Domain::Real);
        let y = pool.symbol("y", Domain::Real);
        let m = Matrix::new(vec![vec![x, y], vec![y, x]]).unwrap();
        let result = id.mul(&m, &pool).unwrap().simplify_entries(&pool);
        assert_eq!(result.get(0, 0), x);
        assert_eq!(result.get(0, 1), y);
    }

    #[test]
    fn det_2x2() {
        let pool = p();
        // det([[a,b],[c,d]]) = ad - bc
        let a = pool.symbol("a", Domain::Real);
        let b = pool.symbol("b", Domain::Real);
        let c = pool.symbol("c", Domain::Real);
        let d = pool.symbol("d", Domain::Real);
        let m = Matrix::new(vec![vec![a, b], vec![c, d]]).unwrap();
        let det = m.det(&pool).unwrap();
        let s = pool.display(det).to_string();
        assert!(s.contains("a") && s.contains("d"), "got: {s}");
    }

    #[test]
    fn det_3x3_identity_is_one() {
        let pool = p();
        let id = Matrix::identity(3, &pool);
        let det = id.det(&pool).unwrap();
        assert_eq!(det, pool.integer(1_i32));
    }

    #[test]
    fn jacobian_linear() {
        // f = [x + y, x - y], vars = [x, y]
        // J = [[1, 1], [1, -1]]
        let pool = p();
        let x = pool.symbol("x", Domain::Real);
        let y = pool.symbol("y", Domain::Real);
        let neg_y = pool.mul(vec![pool.integer(-1_i32), y]);
        let f1 = pool.add(vec![x, y]);
        let f2 = pool.add(vec![x, neg_y]);
        let j = jacobian(&[f1, f2], &[x, y], &pool).unwrap();
        assert_eq!(j.rows, 2);
        assert_eq!(j.cols, 2);
        assert_eq!(j.get(0, 0), pool.integer(1_i32)); // ∂f1/∂x
        assert_eq!(j.get(0, 1), pool.integer(1_i32)); // ∂f1/∂y
        assert_eq!(j.get(1, 0), pool.integer(1_i32)); // ∂f2/∂x
        assert_eq!(j.get(1, 1), pool.integer(-1_i32)); // ∂f2/∂y
    }

    #[test]
    fn jacobian_quadratic() {
        // f = [x², y²], vars = [x, y]
        // J = [[2x, 0], [0, 2y]]
        let pool = p();
        let x = pool.symbol("x", Domain::Real);
        let y = pool.symbol("y", Domain::Real);
        let f1 = pool.pow(x, pool.integer(2_i32));
        let f2 = pool.pow(y, pool.integer(2_i32));
        let j = jacobian(&[f1, f2], &[x, y], &pool).unwrap();
        // ∂f1/∂y = 0, ∂f2/∂x = 0
        assert_eq!(j.get(0, 1), pool.integer(0_i32));
        assert_eq!(j.get(1, 0), pool.integer(0_i32));
    }

    // ---- det_rational_flint vs the Bareiss reference ----------------------

    fn arb_entry() -> impl proptest::strategy::Strategy<Value = rug::Rational> {
        use proptest::prelude::*;
        prop_oneof![
            4 => (-9i64..10).prop_map(rug::Rational::from),
            2 => (-40i64..40, 1i64..9).prop_map(|(a, b)| rug::Rational::from((a, b))),
            1 => (any::<i64>(), any::<i64>(), 1u32..4).prop_map(|(a, b, e)| {
                // Huge entries: well past one limb.
                let v = rug::Integer::from(a) * rug::ops::Pow::pow(rug::Integer::from(b), e) + 1;
                rug::Rational::from(v)
            }),
            1 => Just(rug::Rational::new()),
        ]
    }

    fn arb_square() -> impl proptest::strategy::Strategy<Value = (usize, Vec<rug::Rational>)> {
        use proptest::prelude::*;
        (0usize..8).prop_flat_map(|n| (Just(n), prop::collection::vec(arb_entry(), n * n)))
    }

    proptest::proptest! {
        #[test]
        fn det_flint_matches_bareiss((n, m) in arb_square()) {
            proptest::prop_assert_eq!(det_rational_flint(&m, n), det_bareiss_reference(&m, n));
        }

        #[test]
        fn det_flint_matches_bareiss_singular((n, m) in arb_square(), k in 0usize..8, c in -3i64..4) {
            // Force a dependent row: row n-1 := c · row k.
            let mut m = m;
            if n >= 2 {
                let k = k % (n - 1);
                for j in 0..n {
                    m[(n - 1) * n + j] = rug::Rational::from(c) * m[k * n + j].clone();
                }
                proptest::prop_assert_eq!(det_rational_flint(&m, n), rug::Rational::new());
            }
            proptest::prop_assert_eq!(det_rational_flint(&m, n), det_bareiss_reference(&m, n));
        }
    }

    #[test]
    fn det_flint_small_cases() {
        let r = |a: i64, b: i64| rug::Rational::from((a, b));
        assert_eq!(det_rational_flint(&[], 0), 1);
        assert_eq!(det_rational_flint(&[r(-7, 3)], 1), r(-7, 3));
        let m = [r(1, 2), r(1, 3), r(1, 4), r(1, 5)];
        assert_eq!(det_rational_flint(&m, 2), r(1, 10) - r(1, 12));
        // Through the public entry point: 3×3 with rational entries.
        let pool = p();
        let q = |a: i64, b: i64| pool.rational(a, b);
        let mat = Matrix::new(vec![
            vec![q(1, 2), pool.integer(3), q(-5, 7)],
            vec![pool.integer(0), q(2, 9), pool.integer(4)],
            vec![q(11, 3), pool.integer(-1), q(1, 1)],
        ])
        .unwrap();
        let entries = mat.numeric_entries(&pool).unwrap();
        let expect = det_bareiss_reference(&entries, 3);
        let got = mat.det(&pool).unwrap();
        assert_eq!(
            got,
            pool.rational(expect.numer().clone(), expect.denom().clone())
        );
    }

    /// `cargo test --release -p alkahest-cas det_timing -- --ignored --nocapture`
    #[test]
    #[ignore]
    fn det_timing() {
        // Warm FLINT up so the first row does not carry its one-time setup.
        let _ = det_rational_flint(&[rug::Rational::from(1)], 1);
        for &n in &[3usize, 5, 10, 20, 40, 60] {
            let m: Vec<rug::Rational> = (0..n * n)
                .map(|k| {
                    let v = ((k as i64 * 7919 + 13) % 201) - 100;
                    if k % 5 == 0 {
                        rug::Rational::from((v, 7))
                    } else {
                        rug::Rational::from(v)
                    }
                })
                .collect();
            let reps = if n <= 5 {
                2000
            } else if n <= 20 {
                20
            } else {
                3
            };
            let t = std::time::Instant::now();
            let mut a = rug::Rational::new();
            for _ in 0..reps {
                a = det_bareiss_reference(&m, n);
            }
            let old = t.elapsed().as_secs_f64() * 1e3 / reps as f64;
            let t = std::time::Instant::now();
            let mut b = rug::Rational::new();
            for _ in 0..reps {
                b = det_rational_flint(&m, n);
            }
            let new = t.elapsed().as_secs_f64() * 1e3 / reps as f64;
            assert_eq!(a, b);
            println!("det n={n}: bareiss {old:.3} ms, flint {new:.3} ms");
        }
    }
}