symplex 0.42.0

//! Additional matrix decompositions and structure tests.
//!
//! Everything here is a method on [`Matrix`] (or a free function taking
//! matrices); the module split is purely organisational.
//!
//! * **Decompositions:** [`Matrix::qr`] (Gram–Schmidt, exact with
//!   radicals), [`Matrix::cholesky`], [`Matrix::ldl`], the free function
//!   [`gram_schmidt`].
//! * **Structure tests** (`Option<bool>`, `None` = undecidable):
//!   [`Matrix::is_symmetric`], [`Matrix::is_hermitian`],
//!   [`Matrix::is_orthogonal`], [`Matrix::is_unitary`],
//!   [`Matrix::is_positive_definite`], [`Matrix::is_positive_semidefinite`],
//!   [`Matrix::is_upper_triangular`], [`Matrix::is_lower_triangular`],
//!   [`Matrix::is_diagonal`], [`Matrix::is_identity`], [`Matrix::is_zero`],
//!   [`Matrix::is_skew_symmetric`], [`Matrix::is_nilpotent`].
//! * **Norms:** [`Matrix::norm_1`], [`Matrix::norm_inf`],
//!   [`Matrix::norm_p`] (vectors), plus
//!   [`norm_frobenius`](Matrix::norm_frobenius) in the core module.
//! * **Functions of diagonalizable matrices:**
//!   [`Matrix::matrix_pow_symbolic`], [`Matrix::matrix_sqrt`],
//!   [`Matrix::matrix_log`] (any Jordan form with non-zero eigenvalues).
//! * **Calculus:** [`hessian`], [`wronskian`], [`Matrix::casoratian`].

use crate::api::expr::Ex;
use crate::base::combinatorics::factorial;
use crate::base::errors::SymplexError;
use crate::base::numeric::Q;
use crate::domains::decompositions::{Diagonalization, JordanForm, Ldl, Qr};
use crate::domains::matrix::{
    Matrix, all3, budget_check, ex_is_nonnegative, ex_is_positive, ex_is_zero, is_swell, reop,
    simplify_function_entry,
};

/// Orthogonal basis vectors and the upper-triangular coefficient matrix
/// produced by [`gram_schmidt_cols`], both as raw row-major `Vec`s.
type GramSchmidtParts = (Vec<Vec<Ex>>, Vec<Vec<Ex>>);

/// The product a Cholesky-type factorisation uses ([`Matrix::cholesky`],
/// [`Matrix::ldl`]): `A = L·Lᵀ` or `A = L·Lᴴ`.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Symmetry {
    Symmetric,
    Hermitian,
}

fn invalid(operation: &'static str, reason: impl Into<String>) -> SymplexError {
    SymplexError::invalid_argument(operation, reason)
}

fn failed(operation: &'static str, reason: impl Into<String>) -> SymplexError {
    SymplexError::computation_failed(operation, reason)
}

/// Dot product of two equal-length column vectors given as `Vec<Ex>`:
/// `Σ āₖ·bₖ` (Hermitian) or `Σ aₖ·bₖ`.
fn dot_vec(a: &[Ex], b: &[Ex], hermitian: bool) -> Ex {
    let left = |x: &Ex| if hermitian { x.conjugate() } else { x.clone() };
    let mut acc = &left(&a[0]) * &b[0];
    for (x, y) in a.iter().zip(b).skip(1) {
        acc += &left(x) * y;
    }
    acc
}

/// `√λ` on the principal branch, as `i·√(−λ)` for a negative eigenvalue
/// of a Hermitian matrix (see [`negative_real_eigenvalue`]).
fn principal_sqrt(lam: &Ex, real_spectrum: bool) -> Ex {
    match negative_real_eigenvalue(lam, real_spectrum) {
        Some(minus) => &lam.context().i_unit() * &minus.sqrt(),
        None => lam.sqrt(),
    }
}

/// `ln λ` on the principal branch, as `ln(−λ) + iπ` for a negative
/// eigenvalue of a Hermitian matrix (see [`negative_real_eigenvalue`]).
fn principal_ln(lam: &Ex, real_spectrum: bool) -> Ex {
    match negative_real_eigenvalue(lam, real_spectrum) {
        Some(minus) => {
            let ctx = lam.context();
            minus.ln() + &ctx.i_unit() * &ctx.pi()
        }
        None => lam.ln(),
    }
}

/// Did `diagonalize` fail because the matrix is defective (as opposed to
/// eigenvalues not found or expression swell)?
fn not_diagonalizable(e: &SymplexError) -> bool {
    matches!(e, SymplexError::ComputationFailed { reason, .. }
        if reason.starts_with("matrix is not diagonalizable"))
}

impl Matrix {
    /// `f(A)` through the Jordan form `A = P·J·P⁻¹`: a block `Jₖ(λ)` maps
    /// to the upper triangular Toeplitz matrix with `f⁽ᵈ⁾(λ)/d!`
    /// (`taylor(λ, d)`) on its `d`-th superdiagonal — SymPy's
    /// `_matrix_pow_by_jordan_blocks`.  Used by `matrix_sqrt` and
    /// `matrix_pow_symbolic` for defective matrices, which before 0.38 were
    /// refused ("requires a diagonalizable matrix").
    fn jordan_function(
        &self,
        op: &'static str,
        taylor: &dyn Fn(&Ex, usize) -> Result<Ex, SymplexError>,
    ) -> Result<Matrix, SymplexError> {
        let n = self.nrows();
        let JordanForm { p, j } = self
            .jordan_form()
            .map_err(|e| failed(op, format!("Jordan form unavailable ({e})")))?;
        let zero = self.context().zero();
        let mut fj: Vec<Vec<Ex>> = vec![vec![zero; n]; n];
        let mut col = 0;
        while col < n {
            let lambda = j[(col, col)].clone();
            let mut block = 1;
            while col + block < n
                && j[(col + block - 1, col + block)].is_one_structural()
                && j[(col + block, col + block)] == lambda
            {
                block += 1;
            }
            for d in 0..block {
                let v = taylor(&lambda, d)?;
                for i in 0..(block - d) {
                    fj[col + i][col + i + d] = v.clone();
                }
            }
            col += block;
        }
        let p_inv = p.inv()?;
        let result = p.matmul(&Matrix::new(fj)?)?.matmul(&p_inv)?;
        budget_check(result.iter(), op)?;
        Ok(result.map(simplify_function_entry))
    }
}

/// `−λ` for an eigenvalue `λ` known to be real (`real_spectrum`: the
/// matrix is Hermitian) that is negative and in Cardano form, so that
/// `matrix_sqrt` / `matrix_log` can take the principal values
/// `√λ = i·√(−λ)` and `ln λ = ln(−λ) + iπ`.  Written `√λ`, such a
/// radicand is exactly real but its imaginary part cancels only
/// numerically, so it sits on the branch cut with the sign of rounding
/// noise and does not evaluate: before 0.38 every entry of `√A` of such a
/// 3×3 Hermitian matrix over `ℚ(√2, √3, i)` was unevaluable.
fn negative_real_eigenvalue(lam: &Ex, real_spectrum: bool) -> Option<Ex> {
    if !real_spectrum || !crate::domains::matrix::has_nested_higher_root(lam) {
        return None;
    }
    let z = lam.eval_complex64().ok()?;
    (z.re < 0.0 && z.im.abs() <= 1e-9 * z.re.abs()).then(|| (-lam).eval())
}

fn dependent_error(op: &'static str, j: usize) -> SymplexError {
    failed(
        op,
        format!(
            "vectors are linearly dependent (vector {j} lies in the span of the previous ones)"
        ),
    )
}

// ═══════════════════════════════════════════════════════════════════════════
// Gram–Schmidt / QR
// ═══════════════════════════════════════════════════════════════════════════

/// Gram–Schmidt orthogonalisation of a list of column vectors.
///
/// Returns vectors spanning the same space, pairwise orthogonal (and of
/// unit length if `normalize` is set) for the Hermitian inner product
/// `⟨x, y⟩ = Σ x̄ₖ·yₖ`, as SymPy's `GramSchmidt` (before 0.32 the bilinear
/// `Σ xₖ·yₖ`; the two agree for real vectors).  Radicals are kept exact: a
/// normalised entry is `uᵢ · ‖u‖⁻¹` with the norm left as `√(‖u‖²)`, so
/// that dot products of the results cancel structurally (`qᵢ·qⱼ` is
/// `0`, `qᵢ·qᵢ` is `1` without further simplification).  Symbolic entries
/// are simplified; constant entries are only constant-folded.
///
/// # Errors
///
/// - [`SymplexError::InvalidArgument`] if `vectors` is empty, any vector
///   is not a column vector, or lengths differ.
/// - [`SymplexError::ComputationFailed`] if the vectors are linearly
///   dependent (a residual vector is provably zero), or if symbolic
///   entries swell beyond
///   [`EXPRESSION_BUDGET`](crate::domains::matrix::EXPRESSION_BUDGET).
///   Symbolic residuals whose zero-ness cannot be decided are assumed
///   non-zero.
///
/// # Examples
///
/// ```
/// use symplex::prelude::*;
/// use symplex::matrix_decomp::gram_schmidt;
///
/// let ctx = Context::new();
/// let v1 = matrix![ctx, [1], [1]];
/// let v2 = matrix![ctx, [1], [0]];
/// let q = gram_schmidt(&[v1, v2], true).unwrap();
/// // q0 = (1/√2, 1/√2), q1 = (1/√2, −1/√2)
/// assert_eq!(q[0][(0, 0)].powi(2), ctx.rational(1, 2));
/// assert_eq!(symplex::matrix::dot(&q[0], &q[1]).unwrap(), ctx.int(0));
/// assert_eq!(symplex::matrix::dot(&q[1], &q[1]).unwrap(), ctx.int(1));
/// ```
pub fn gram_schmidt(vectors: &[Matrix], normalize: bool) -> Result<Vec<Matrix>, SymplexError> {
    if vectors.is_empty() {
        return Err(invalid("gram_schmidt", "need at least one vector"));
    }
    let n = vectors[0].nrows();
    for (idx, v) in vectors.iter().enumerate() {
        if v.ncols() != 1 {
            return Err(invalid(
                "gram_schmidt",
                format!(
                    "vector {idx} is {}×{}, expected a column vector",
                    v.nrows(),
                    v.ncols()
                ),
            ));
        }
        if v.nrows() != n {
            return Err(invalid(
                "gram_schmidt",
                format!("vector {idx} has length {}, expected {n}", v.nrows()),
            ));
        }
    }
    let cols: Vec<Vec<Ex>> = vectors.iter().map(|v| v.col(0)).collect();
    let (basis, _) = gram_schmidt_cols(&cols, normalize, "gram_schmidt")?;
    basis.into_iter().map(Matrix::col_vector).collect()
}

/// Core Gram–Schmidt on raw columns.
///
/// Returns `(orthogonal vectors, R)` where `R` is the `k×k` upper
/// triangular coefficient matrix (`R[i][j] = ⟨q_i, a_j⟩` for `i < j`,
/// `R[j][j] = ‖u_j‖`) when `normalize` is true.
///
/// The inner product is the Hermitian `⟨x, y⟩ = Σ x̄ₖ·yₖ`, as SymPy's
/// `QRdecomposition`/`GramSchmidt` (`u.dot(v, hermitian=True)`).  Before
/// 0.32 it was the bilinear `Σ xₖ·yₖ`, so `Q` satisfied `QᵀQ = I` rather
/// than `QᴴQ = I` for complex entries, and a column such as `[1, i]` (with
/// `1 + i² = 0`) was "linearly dependent".  When every entry is provably
/// real the two agree and the real path ([`gram_schmidt_expr`]) runs
/// unchanged, falling back to the exact path only where it fails (a 4×4
/// matrix over `ℚ(√2, √3)` exceeded the expression budget); otherwise the
/// columns are orthogonalised exactly over `ℚ(√p…, i)(params,
/// conj(params))` ([`hermitian_gram_schmidt`](crate::domains::linalg::hermitian_gram_schmidt)),
/// or, for entries outside that field, on expressions with conjugates.
fn gram_schmidt_cols(
    cols: &[Vec<Ex>],
    normalize: bool,
    op: &'static str,
) -> Result<GramSchmidtParts, SymplexError> {
    use crate::domains::linalg::{HermitianGs, hermitian_gram_schmidt};
    let exact = || match hermitian_gram_schmidt(cols, normalize) {
        Some(HermitianGs::Factors(f)) => Some(Ok((f.q, f.r))),
        Some(HermitianGs::Dependent(j)) => Some(Err(dependent_error(op, j))),
        None => None,
    };
    let hermitian = !cols.iter().flatten().all(|e| e.is_real() == Some(true));
    if hermitian {
        if let Some(res) = exact() {
            return res;
        }
        return gram_schmidt_expr(cols, normalize, op, true);
    }
    match gram_schmidt_expr(cols, normalize, op, false) {
        Err(e) => exact().unwrap_or(Err(e)),
        ok => ok,
    }
}

/// Gram–Schmidt on expressions (see [`gram_schmidt_cols`]), with the
/// Hermitian or the bilinear inner product.
fn gram_schmidt_expr(
    cols: &[Vec<Ex>],
    normalize: bool,
    op: &'static str,
    hermitian: bool,
) -> Result<GramSchmidtParts, SymplexError> {
    let k = cols.len();
    let zero = cols[0][0].context().zero();
    let mut q: Vec<Vec<Ex>> = Vec::with_capacity(k);
    let mut r: Vec<Vec<Ex>> = vec![vec![zero.clone(); k]; k];
    // Constant entries only need folding; `simplify` would rewrite the
    // radicals inconsistently (`√(2/3)` vs `2^(-1/2)·√3`) and break the
    // structural cancellation of dot products.
    let tidy = |e: Ex| {
        if e.is_constant() {
            e.eval()
        } else {
            e.simplify()
        }
    };

    for j in 0..k {
        let mut u = cols[j].clone();
        for i in 0..j {
            // r_ij = q_i · a_j  (q_i normalised) or (q_i · a_j)/(q_i · q_i) otherwise
            let proj = if normalize {
                tidy(dot_vec(&q[i], &cols[j], hermitian))
            } else {
                let qq = dot_vec(&q[i], &q[i], hermitian);
                tidy(&dot_vec(&q[i], &cols[j], hermitian) / &qq)
            };
            for (u_e, q_e) in u.iter_mut().zip(q[i].iter()) {
                *u_e = &*u_e - &(&proj * q_e);
            }
            r[i][j] = proj;
        }
        budget_check(u.iter(), op)?;
        let u: Vec<Ex> = u.into_iter().map(tidy).collect();
        let norm_sq = tidy(dot_vec(&u, &u, hermitian));
        // The exact test of a rational function over ℚ(radicals, i) first:
        // `ex_is_zero` would `simplify` a norm that is not identically zero,
        // which swelled for 3×3 symbolic matrices with radicals (> 10 s).
        let dependent = match crate::domains::matrix::rational_function_is_zero(&norm_sq) {
            Some(b) => b,
            None => ex_is_zero(&norm_sq) == Some(true),
        };
        if dependent {
            return Err(dependent_error(op, j));
        }
        if normalize {
            let norm = norm_sq.sqrt();
            let one = cols[0][0].context().one();
            // `√(1/n)` displays nicely for rational `n > 0` and cancels
            // against `√n` structurally (canonical numeric radicals); for
            // symbolic `n` only `n^(-1/2)` is guaranteed to cancel.  For
            // `n < 0` the principal `√(1/n)` is `−1/√n` (before 0.31 the
            // column of `Q·R` was then `−aⱼ`: `[0, 1/2, −i]` has `n = −3/4`).
            let positive_rational = norm_sq
                .as_rational()
                .is_some_and(|q| num_traits::Signed::is_positive(&q));
            let inv_norm = if positive_rational {
                (&one / &norm_sq).sqrt()
            } else {
                &one / &norm
            };
            r[j][j] = norm;
            q.push(u.iter().map(|e| tidy(e * &inv_norm)).collect());
        } else {
            r[j][j] = cols[0][0].context().one();
            q.push(u);
        }
    }
    Ok((q, r))
}

impl Matrix {
    /// QR decomposition via Gram–Schmidt: [`Qr`]`{ q, r }` with `A = Q·R`,
    /// `Q` (m×n) having orthonormal columns (`QᴴQ = I`) and `R` (n×n) upper
    /// triangular with a real positive diagonal.
    ///
    /// As SymPy's `QRdecomposition`, the inner product is the Hermitian
    /// `⟨x, y⟩ = Σ x̄ₖ·yₖ`, so `Q` is unitary for complex entries (a symbol
    /// without a real assumption is complex: its entries carry
    /// `conjugate(x)`).  Before 0.32 the bilinear `Σ xₖ·yₖ` gave `QᵀQ = I`
    /// instead, and `[1, i]ᵀ` was refused as linearly dependent.  Real
    /// matrices are unaffected.  Complex and symbolic entries that are
    /// rational functions over square roots of rationals and `i` are
    /// orthogonalised in exact field arithmetic, simplified once.
    ///
    /// Works exactly with radicals (entries like `1/√2`).
    ///
    /// A rank-deficient `A` (rank `r < n`, decided by [`rref`](Matrix::rref))
    /// gets SymPy's reduced form: `Q` (m×r) is orthonormalised from the
    /// pivot columns, the linearly independent columns from the left, and
    /// `R = QᴴA` (r×n) is upper trapezoidal in echelon form (row `k` is
    /// zero left of the `k`-th pivot column).  Before 0.38 it was refused.
    ///
    /// # Errors
    ///
    /// Returns [`SymplexError::ComputationFailed`] for the zero matrix (`Q`
    /// and `R` would be empty) or if symbolic entries swell beyond
    /// [`EXPRESSION_BUDGET`](crate::domains::matrix::EXPRESSION_BUDGET).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let a = matrix![ctx, [1, 1], [0, 1]];
    /// let Qr { q, r } = a.qr().unwrap();
    /// assert_eq!((&q * &r).simplify(), a);
    /// assert_eq!((&q.transpose() * &q).simplify(), Matrix::identity(&ctx, 2).unwrap());
    /// assert!(r[(1, 0)].is_zero_structural());
    ///
    /// // Complex: SymPy's Matrix([1, I]).QRdecomposition() is
    /// // Q = [sqrt(2)/2, sqrt(2)*I/2], R = [sqrt(2)].
    /// let v = Matrix::col_vector(vec![ctx.int(1), ctx.i_unit()]).unwrap();
    /// let Qr { q, r } = v.qr().unwrap();
    /// assert_eq!(q[(1, 0)], (ctx.rational(1, 2).sqrt() * ctx.i_unit()).eval());
    /// assert_eq!(r[(0, 0)], ctx.int(2).sqrt());
    /// assert_eq!(q.adjoint().matmul(&q).unwrap().eval(), Matrix::identity(&ctx, 1).unwrap());
    /// ```
    pub fn qr(&self) -> Result<Qr<Matrix>, SymplexError> {
        budget_check(self.iter(), "qr")?;
        let (_, pivots) = self.rref();
        if pivots.is_empty() {
            return Err(failed(
                "qr",
                "matrix is zero (rank 0): Q (m×0) and R (0×n) would be empty",
            ));
        }
        let n = self.ncols();
        let cols: Vec<Vec<Ex>> = pivots.iter().map(|&j| self.col(j)).collect();
        let (q_cols, mut r_rows) = gram_schmidt_cols(&cols, true, "qr")?;
        if pivots.len() < n {
            r_rows = self.qr_spread_r(&pivots, &q_cols, &r_rows);
        }
        let q_mats: Vec<Matrix> = q_cols
            .into_iter()
            .map(Matrix::col_vector)
            .collect::<Result<_, _>>()?;
        let q_refs: Vec<&Matrix> = q_mats.iter().collect();
        let q = Matrix::hstack(&q_refs)?;
        let r = Matrix::new(r_rows)?;
        Ok(Qr { q, r })
    }

    /// The `r×n` factor `R = QᴴA` of a rank-deficient [`qr`](Self::qr) from
    /// the `r×r` factor `r_piv` of the pivot columns: column `pivots[c]` is
    /// column `c` of `r_piv`, and a non-pivot column `j` has `⟨q_k, a_j⟩`
    /// in the rows `k` with `pivots[k] < j` — the others are zero, since
    /// `a_j` lies in the span of the pivot columns left of it, which is
    /// that of `q_0, …, q_{k−1}`.
    fn qr_spread_r(&self, pivots: &[usize], q_cols: &[Vec<Ex>], r_piv: &[Vec<Ex>]) -> Vec<Vec<Ex>> {
        let n = self.ncols();
        let zero = self.context().zero();
        let hermitian = !q_cols.iter().flatten().all(|e| e.is_real() == Some(true));
        let tidy = |e: Ex| {
            if e.is_constant() {
                crate::domains::linalg::algebraic_constant_normal_form(&e)
                    .unwrap_or_else(|| e.eval())
            } else {
                e.simplify()
            }
        };
        let mut full = vec![vec![zero; n]; pivots.len()];
        for (k, row) in full.iter_mut().enumerate() {
            for (c, &pj) in pivots.iter().enumerate() {
                row[pj] = r_piv[k][c].clone();
            }
            for j in (pivots[k] + 1..n).filter(|j| !pivots.contains(j)) {
                row[j] = tidy(dot_vec(&q_cols[k], &self.col(j), hermitian));
            }
        }
        full
    }

    // ── Cholesky / LDLᵀ ────────────────────────────────────────────────

    /// Cholesky decomposition `A = L·Lᵀ` for a symmetric positive-definite
    /// matrix, and `A = L·Lᴴ` (conjugate transpose) for a Hermitian one
    /// that is not symmetric (`L` lower triangular with positive diagonal).
    ///
    /// As SymPy's `cholesky` (`hermitian=True`, its default), a Hermitian
    /// matrix takes `lᵢⱼ = (aᵢⱼ − Σₖ lᵢₖ·l̄ⱼₖ)/lⱼⱼ` and the pivots `aⱼⱼ − Σₖ |lⱼₖ|²`:
    /// `[[2, i], [−i, 2]]` gives `L = [[√2, 0], [−i·√2/2, √6/2]]` (refused as
    /// "not symmetric" before 0.36).  A matrix provably symmetric — every
    /// real one, and a complex symmetric one such as `[[2, i], [i, 2]]`,
    /// which SymPy takes only with `hermitian=False` — keeps `A = L·Lᵀ`.
    ///
    /// Positive-definiteness is decided pivot by pivot: each diagonal
    /// pivot must be provably positive (assumption system, or numeric
    /// evaluation for constant entries).
    ///
    /// # Errors
    ///
    /// - [`SymplexError::InvalidArgument`] if the matrix is not square, or
    ///   provably not symmetric and not provably Hermitian.
    /// - [`SymplexError::ComputationFailed`] if a pivot is provably
    ///   non-positive (not positive definite) **or** its sign cannot be
    ///   decided symbolically.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let a = matrix![ctx, [4, 2], [2, 3]];
    /// let l = a.cholesky().unwrap();
    /// assert_eq!((&l * &l.transpose()).simplify(), a);
    /// assert!(matrix![ctx, [1, 2], [2, 1]].cholesky().is_err()); // indefinite
    ///
    /// // Hermitian: SymPy's Matrix([[9, 3*I], [-3*I, 5]]).cholesky() is
    /// // [[3, 0], [-I, 2]].
    /// let i = ctx.i_unit();
    /// let h = Matrix::new(vec![vec![ctx.int(9), &i * 3], vec![&i * -3, ctx.int(5)]]).unwrap();
    /// let l = h.cholesky().unwrap();
    /// assert_eq!(l, Matrix::new(vec![vec![ctx.int(3), ctx.int(0)], vec![-&i, ctx.int(2)]]).unwrap());
    /// assert_eq!(l.matmul(&l.adjoint()).unwrap().eval(), h);
    /// ```
    pub fn cholesky(&self) -> Result<Matrix, SymplexError> {
        if !self.is_square() {
            return Err(invalid(
                "cholesky",
                format!(
                    "requires a square matrix, got {}×{}",
                    self.nrows(),
                    self.ncols()
                ),
            ));
        }
        let hermitian = self.symmetry_kind("cholesky")? == Symmetry::Hermitian;
        let n = self.nrows();
        let zero = self.context().zero();
        let mut l: Vec<Vec<Ex>> = vec![vec![zero.clone(); n]; n];
        // `x̄` in the Hermitian mode, `x` otherwise.
        let conj = |x: &Ex| if hermitian { x.conjugate() } else { x.clone() };

        for j in 0..n {
            let mut sum_sq = zero.clone();
            for item in l[j].iter().take(j) {
                if hermitian {
                    sum_sq += item * &item.conjugate();
                } else {
                    sum_sq += item.powi(2);
                }
            }
            let diag = (self.get(j, j) - &sum_sq).simplify();
            // Over one denominator a pivot such as `2 − x²/(x² + 1)` (real `x`)
            // is `(x² + 2)/(x² + 1)`, whose sign the assumptions decide.
            match ex_is_positive(&diag).or_else(|| ex_is_positive(&diag.together())) {
                Some(true) => {}
                Some(false) => {
                    return Err(failed(
                        "cholesky",
                        format!("matrix is not positive definite (pivot {j} is {diag})"),
                    ));
                }
                None => {
                    return Err(failed(
                        "cholesky",
                        format!(
                            "cannot decide the sign of pivot {j} = {diag}; add assumptions or use ldl()"
                        ),
                    ));
                }
            }
            l[j][j] = diag.sqrt().simplify();
            for i in (j + 1)..n {
                let mut sum_prod = zero.clone();
                for (l_ik, l_jk) in l[i].iter().zip(l[j].iter()).take(j) {
                    sum_prod += &(l_ik * &conj(l_jk));
                }
                let num = self.get(i, j) - &sum_prod;
                l[i][j] = (&num / &l[j][j]).simplify();
            }
        }
        Matrix::new(l)
    }

    /// Which product [`cholesky`](Self::cholesky) and [`ldl`](Self::ldl)
    /// factor `A` with: `Lᵀ` unless `A` is provably Hermitian and not
    /// provably symmetric (a matrix that is both is real), `Lᴴ` then.
    /// `InvalidArgument` when `A` is provably not symmetric and not
    /// provably Hermitian; `ComputationFailed` when neither can be decided.
    ///
    /// Before 0.38 an undecided symmetry was taken as symmetric: the
    /// factors were computed from the lower triangle alone, so
    /// `[[x², x + y], [2x, 0]]` came back with `L·D·Lᵀ = [[x², 2x], [2x, 0]]`
    /// (SymPy: "Matrix must be Hermitian.").
    fn symmetry_kind(&self, op: &'static str) -> Result<Symmetry, SymplexError> {
        let symmetric = self.is_symmetric();
        if symmetric == Some(true) {
            return Ok(Symmetry::Symmetric);
        }
        if self.is_hermitian() == Some(true) {
            return Ok(Symmetry::Hermitian);
        }
        if symmetric == Some(false) {
            return Err(invalid(op, "matrix is neither symmetric nor Hermitian"));
        }
        Err(failed(
            op,
            "cannot decide whether the matrix is symmetric or Hermitian \
             (an entry differs from its transposed entry by an expression not provably 0)",
        ))
    }

    /// LDLᵀ decomposition [`Ldl`]`{ l, d }` with `A = L·D·Lᵀ` for a
    /// symmetric matrix, and `A = L·D·Lᴴ` for a Hermitian one that is not
    /// symmetric (`L` unit lower triangular, `D` diagonal) — SymPy's
    /// `LDLdecomposition` with `hermitian=True`, its default:
    /// `[[9, 3i], [−3i, 5]]` gives `L = [[1, 0], [−i/3, 1]]`, `D = diag(9, 4)`
    /// (refused as "not symmetric" before 0.36).  See
    /// [`cholesky`](Self::cholesky) for which matrices take which form.
    ///
    /// Unlike [`cholesky`](Self::cholesky) this needs no square roots and
    /// no sign information, so it works for symbolic symmetric matrices
    /// and for indefinite ones (SymPy refuses an indefinite Hermitian one).
    ///
    /// # Errors
    ///
    /// - [`SymplexError::InvalidArgument`] if the matrix is not square, or
    ///   provably not symmetric and not provably Hermitian.
    /// - [`SymplexError::ComputationFailed`] if a zero pivot is
    ///   encountered with a nonzero entry below it (no pivoting is
    ///   performed).  A zero pivot whose column is zero below it needs no
    ///   division and is kept (`d_j = 0`, `l_ij = 0`): `[[0]]` and
    ///   `[[1, 1], [1, 1]]` decompose (before 0.30 they were errors).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let a = matrix![ctx, [4, 2], [2, 3]];
    /// let Ldl { l, d } = a.ldl().unwrap();
    /// assert_eq!((&(&l * &d) * &l.transpose()).eval(), a);
    /// assert_eq!(d, matrix![ctx, [4, 0], [0, 2]]);
    /// ```
    pub fn ldl(&self) -> Result<Ldl<Matrix>, SymplexError> {
        if !self.is_square() {
            return Err(invalid(
                "ldl",
                format!(
                    "requires a square matrix, got {}×{}",
                    self.nrows(),
                    self.ncols()
                ),
            ));
        }
        let hermitian = self.symmetry_kind("ldl")? == Symmetry::Hermitian;
        let n = self.nrows();
        let ctx = self.context();
        let zero = ctx.zero();
        let one = ctx.one();
        let mut l: Vec<Vec<Ex>> = vec![vec![zero.clone(); n]; n];
        let mut d: Vec<Ex> = vec![zero.clone(); n];
        // `x̄` in the Hermitian mode, `x` otherwise.
        let conj = |x: &Ex| if hermitian { x.conjugate() } else { x.clone() };

        for j in 0..n {
            l[j][j] = one.clone();
            // d_j = a_jj − Σ_{k<j} l_jk² d_k  (l_jk·l̄_jk Hermitian)
            let mut acc = zero.clone();
            for k in 0..j {
                if hermitian {
                    acc += &(&l[j][k] * &l[j][k].conjugate()) * &d[k];
                } else {
                    acc += &l[j][k].powi(2) * &d[k];
                }
            }
            let dj = (self.get(j, j) - &acc).simplify();
            // l_ij = (a_ij − Σ_{k<j} l_ik l_jk d_k) / d_j  (l̄_jk Hermitian)
            let numerators: Vec<Ex> = ((j + 1)..n)
                .map(|i| {
                    let mut acc = zero.clone();
                    for k in 0..j {
                        acc += &(&l[i][k] * &conj(&l[j][k])) * &d[k];
                    }
                    self.get(i, j) - &acc
                })
                .collect();
            if ex_is_zero(&dj) == Some(true) {
                if numerators.iter().all(|u| ex_is_zero(u) == Some(true)) {
                    // Nothing to divide: the column below is zero too.
                    d[j] = zero.clone();
                    continue;
                }
                return Err(failed(
                    "ldl",
                    format!("zero pivot at position {j}; matrix needs pivoting or is singular"),
                ));
            }
            d[j] = dj;
            for (i, u) in ((j + 1)..n).zip(numerators) {
                l[i][j] = (&u / &d[j]).simplify();
            }
        }
        Ok(Ldl {
            l: Matrix::new(l)?,
            d: Matrix::diag(&d)?,
        })
    }

    // ── Structure tests ────────────────────────────────────────────────

    /// Is `A = Aᵀ`?  Three-valued; `Some(false)` for non-square.
    pub fn is_symmetric(&self) -> Option<bool> {
        if !self.is_square() {
            return Some(false);
        }
        let n = self.nrows();
        all3((0..n).flat_map(|i| {
            ((i + 1)..n).map(move |j| ex_is_zero(&(self.get(i, j) - self.get(j, i))))
        }))
    }

    /// Is `A = −Aᵀ`?  Three-valued; `Some(false)` for non-square.
    pub fn is_skew_symmetric(&self) -> Option<bool> {
        if !self.is_square() {
            return Some(false);
        }
        let n = self.nrows();
        all3(
            (0..n)
                .flat_map(|i| (i..n).map(move |j| ex_is_zero(&(self.get(i, j) + self.get(j, i))))),
        )
    }

    /// Is `A = Aᴴ` (conjugate transpose)?  Three-valued.
    pub fn is_hermitian(&self) -> Option<bool> {
        if !self.is_square() {
            return Some(false);
        }
        let adj = self.adjoint();
        self.equals(&adj)
    }

    /// Is `AᵀA = I`?  Three-valued; `Some(false)` for non-square.
    pub fn is_orthogonal(&self) -> Option<bool> {
        if !self.is_square() {
            return Some(false);
        }
        let prod = self.transpose().matmul(self).ok()?;
        prod.is_identity()
    }

    /// Is `AᴴA = I`?  Three-valued; `Some(false)` for non-square.
    pub fn is_unitary(&self) -> Option<bool> {
        if !self.is_square() {
            return Some(false);
        }
        let prod = self.adjoint().matmul(self).ok()?;
        prod.is_identity()
    }

    /// Are all entries below the main diagonal zero?  Three-valued.
    pub fn is_upper_triangular(&self) -> Option<bool> {
        all3(
            (0..self.nrows())
                .flat_map(|i| (0..i.min(self.ncols())).map(move |j| ex_is_zero(self.get(i, j)))),
        )
    }

    /// Are all entries above the main diagonal zero?  Three-valued.
    pub fn is_lower_triangular(&self) -> Option<bool> {
        all3(
            (0..self.nrows())
                .flat_map(|i| ((i + 1)..self.ncols()).map(move |j| ex_is_zero(self.get(i, j)))),
        )
    }

    /// Are all off-diagonal entries zero?  Three-valued (works for
    /// rectangular matrices too).
    pub fn is_diagonal(&self) -> Option<bool> {
        all3((0..self.nrows()).flat_map(|i| {
            (0..self.ncols())
                .filter(move |&j| j != i)
                .map(move |j| ex_is_zero(self.get(i, j)))
        }))
    }

    /// Is this the identity matrix?  Three-valued; `Some(false)` for
    /// non-square.
    pub fn is_identity(&self) -> Option<bool> {
        if !self.is_square() {
            return Some(false);
        }
        let one = self.context().one();
        let one = &one;
        all3((0..self.nrows()).flat_map(|i| {
            (0..self.ncols()).map(move |j| {
                if i == j {
                    ex_is_zero(&(self.get(i, j) - one))
                } else {
                    ex_is_zero(self.get(i, j))
                }
            })
        }))
    }

    /// Are all entries zero?  Three-valued.
    pub fn is_zero(&self) -> Option<bool> {
        all3(self.iter().map(ex_is_zero))
    }

    /// Is `Aⁿ = 0` (n = dimension)?  Three-valued; `Some(false)` for
    /// non-square.
    ///
    /// A square matrix is nilpotent iff `Aⁿ = 0`, so exactly one matrix
    /// power is examined.
    pub fn is_nilpotent(&self) -> Option<bool> {
        if !self.is_square() {
            return Some(false);
        }
        let p = self.powi(self.nrows() as u32).ok()?;
        p.expand().is_zero()
    }

    /// Positive-definiteness via Sylvester's criterion: symmetric and all
    /// leading principal minors strictly positive.  Three-valued.
    ///
    /// Returns `Some(false)` for non-square or provably non-symmetric
    /// matrices and `None` when a minor's sign cannot be decided
    /// (symbolic entries without assumptions).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// assert_eq!(matrix![ctx, [2, -1], [-1, 2]].is_positive_definite(), Some(true));
    /// assert_eq!(matrix![ctx, [1, 2], [2, 1]].is_positive_definite(), Some(false));
    /// let x = ctx.symbol("x");
    /// let m = Matrix::new(vec![vec![x.clone(), ctx.int(0)], vec![ctx.int(0), ctx.int(1)]]).unwrap();
    /// assert_eq!(m.is_positive_definite(), None);
    /// ```
    pub fn is_positive_definite(&self) -> Option<bool> {
        match self.is_symmetric() {
            Some(true) => {}
            Some(false) => return Some(false),
            None => return None,
        }
        let n = self.nrows();
        all3((1..=n).map(|k| {
            let minor = self.submatrix(0..k, 0..k).ok()?.det().ok()?;
            ex_is_positive(&minor)
        }))
    }

    /// Positive-semidefiniteness: symmetric and **all** principal minors
    /// non-negative (leading minors alone are not sufficient).
    /// Three-valued.
    ///
    /// Examines `2ⁿ − 1` minors, so this is intended for small matrices.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// assert_eq!(matrix![ctx, [1, 1], [1, 1]].is_positive_semidefinite(), Some(true));
    /// assert_eq!(matrix![ctx, [0, 0], [0, -1]].is_positive_semidefinite(), Some(false));
    /// ```
    pub fn is_positive_semidefinite(&self) -> Option<bool> {
        match self.is_symmetric() {
            Some(true) => {}
            Some(false) => return Some(false),
            None => return None,
        }
        let n = self.nrows();
        all3((1u32..(1u32 << n)).map(|mask| {
            let idx: Vec<usize> = (0..n).filter(|&i| mask & (1 << i) != 0).collect();
            let rows: Vec<Vec<Ex>> = idx
                .iter()
                .map(|&i| idx.iter().map(|&j| self.get(i, j).clone()).collect())
                .collect();
            let minor = Matrix::new(rows).ok()?.det().ok()?;
            ex_is_nonnegative(&minor)
        }))
    }

    // ── Norms ──────────────────────────────────────────────────────────

    /// Induced 1-norm: maximum absolute column sum.
    ///
    /// For symbolic entries the result is a `max(…)` of `abs(…)` sums;
    /// call `.eval()` to fold constants.
    pub fn norm_1(&self) -> Ex {
        let ctx = self.context();
        let sums = (0..self.ncols()).map(|j| {
            let mut acc = self.get(0, j).abs();
            for i in 1..self.nrows() {
                acc += self.get(i, j).abs();
            }
            acc
        });
        Ex::max_of(&ctx, sums).eval()
    }

    /// Induced ∞-norm: maximum absolute row sum.
    pub fn norm_inf(&self) -> Ex {
        let ctx = self.context();
        let sums = (0..self.nrows()).map(|i| {
            let mut acc = self.get(i, 0).abs();
            for j in 1..self.ncols() {
                acc += self.get(i, j).abs();
            }
            acc
        });
        Ex::max_of(&ctx, sums).eval()
    }

    /// Vector p-norm `(Σ |xᵢ|ᵖ)^(1/p)` for a row or column vector.
    ///
    /// # Errors
    ///
    /// Returns [`SymplexError::InvalidArgument`] if the matrix is not a
    /// row or column vector.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let v = matrix![ctx, [3], [-4]];
    /// assert_eq!(v.norm_p(&ctx.int(2)).unwrap().eval(), ctx.int(5));
    /// assert_eq!(v.norm_p(&ctx.int(1)).unwrap().eval(), ctx.int(7));
    /// ```
    pub fn norm_p(&self, p: &Ex) -> Result<Ex, SymplexError> {
        if self.nrows() != 1 && self.ncols() != 1 {
            return Err(invalid(
                "norm_p",
                format!(
                    "requires a row or column vector, got {}×{}",
                    self.nrows(),
                    self.ncols()
                ),
            ));
        }
        let ctx = self.context();
        let mut acc = ctx.zero();
        for e in self.iter() {
            acc += e.abs().pow(p);
        }
        let inv_p = &ctx.one() / p;
        Ok(acc.pow(&inv_p))
    }

    // ── Functions of diagonalizable matrices ───────────────────────────

    /// Symbolic power `Aⁿ` for a diagonalizable matrix via `P·Dⁿ·P⁻¹`.
    ///
    /// `n` may be any expression (a symbol, a rational, …); each diagonal
    /// entry becomes `λᵢⁿ`.  Entries are simplified.  A defective matrix
    /// goes through its Jordan form, a block `Jₖ(λ)` taking
    /// `binomial(n, d)·λⁿ⁻ᵈ` on its `d`-th superdiagonal (SymPy's
    /// `_matrix_pow_by_jordan_blocks`; refused before 0.38).  When the
    /// eigenvector matrix is too large to invert and the eigenvalues are
    /// distinct, Sylvester's formula is used instead.
    ///
    /// # Errors
    ///
    /// Same conditions as [`jordan_form`](Matrix::jordan_form); also a
    /// nilpotent Jordan block of size > 1 (SymPy leaves `A**n` unevaluated).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let n = ctx.symbol("n");
    /// let a = matrix![ctx, [2, 0], [0, 3]];
    /// let an = a.matrix_pow_symbolic(&n).unwrap();
    /// assert_eq!(an[(0, 0)], ctx.int(2).pow(&n));
    /// assert_eq!(an[(1, 1)], ctx.int(3).pow(&n));
    /// ```
    pub fn matrix_pow_symbolic(&self, n: &Ex) -> Result<Matrix, SymplexError> {
        let op = "matrix_pow_symbolic";
        let Diagonalization { p, d } = match self.diagonalize() {
            Ok(dg) => dg,
            Err(e) if not_diagonalizable(&e) => {
                let ctx = self.context();
                return self.jordan_function(op, &|lam, k| {
                    if k == 0 {
                        return Ok(lam.pow(n));
                    }
                    if ex_is_zero(lam) == Some(true) {
                        return Err(failed(
                            op,
                            "a nilpotent Jordan block (eigenvalue 0, size > 1) has no closed-form power for a symbolic exponent",
                        ));
                    }
                    // binomial(n, k)·λ^(n−k), as SymPy's `jordan_cell_power`
                    let mut binom = ctx.one();
                    for i in 0..k {
                        binom = &binom * &(n - &ctx.int(i as i64));
                    }
                    let fact = ctx.from_bigint(factorial(k as u64));
                    Ok(&(&binom / &fact) * &lam.pow(&(n - &ctx.int(k as i64))))
                });
            }
            Err(e) if is_swell(&e) => return self.sylvester_or(op, e, &|l| l.pow(n)),
            Err(e) => return Err(failed(op, format!("requires a diagonalizable matrix: {e}"))),
        };
        let dn = d.map_indexed(|i, j, e| if i == j { e.pow(n) } else { e.clone() });
        let p_inv = match p.inv() {
            Ok(pi) => pi,
            Err(e) if is_swell(&e) => return self.sylvester_or(op, e, &|l| l.pow(n)),
            Err(e) => return Err(e),
        };
        Ok(p.matmul(&dn)?.matmul(&p_inv)?.map(simplify_function_entry))
    }

    /// `f(A)` by [`Matrix::sylvester_function`] after the eigenvector route
    /// failed with `e` on the expression budget; `e` again when the
    /// eigenvalues are not distinct.
    fn sylvester_or(
        &self,
        op: &'static str,
        e: SymplexError,
        f: &dyn Fn(&Ex) -> Ex,
    ) -> Result<Matrix, SymplexError> {
        self.sylvester_function(op, f)?.ok_or(e)
    }

    /// Principal square root `√A` of a diagonalizable matrix via
    /// `P·√D·P⁻¹` (principal branch on each eigenvalue).  A defective
    /// matrix goes through its Jordan form, a block `Jₖ(λ)` taking
    /// `binomial(1/2, d)·√λ/λᵈ` on its `d`-th superdiagonal (refused
    /// before 0.38).  When the eigenvector matrix is too large to invert and
    /// the eigenvalues are distinct, Sylvester's formula is used instead.
    ///
    /// # Errors
    ///
    /// Same conditions as [`jordan_form`](Matrix::jordan_form); also a
    /// nilpotent Jordan block of size > 1, which has no square root of this
    /// form.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let a = matrix![ctx, [4, 0], [0, 9]];
    /// let s = a.matrix_sqrt().unwrap();
    /// assert_eq!((&s * &s).simplify(), a);
    /// ```
    pub fn matrix_sqrt(&self) -> Result<Matrix, SymplexError> {
        let op = "matrix_sqrt";
        let real_spectrum = self.is_hermitian() == Some(true);
        let Diagonalization { p, d } = match self.diagonalize() {
            Ok(dg) => dg,
            Err(e) if not_diagonalizable(&e) => {
                let ctx = self.context();
                return self.jordan_function(op, &|lam, k| {
                    if k == 0 {
                        return Ok(lam.sqrt());
                    }
                    if ex_is_zero(lam) == Some(true) {
                        return Err(failed(
                            op,
                            "a nilpotent Jordan block (eigenvalue 0, size > 1) has no square root of this form",
                        ));
                    }
                    // binomial(1/2, k)·√λ/λᵏ
                    let mut binom = Q::from_integer(1.into());
                    for i in 0..k {
                        binom = binom * (Q::new(1.into(), 2.into()) - Q::from_integer(i.into()))
                            / Q::from_integer((i + 1).into());
                    }
                    Ok(&(&ctx.from_ratio(binom) * &lam.sqrt()) / &lam.powi(k as i64))
                });
            }
            Err(e) if is_swell(&e) => {
                return self.sylvester_or(op, e, &|l| principal_sqrt(l, real_spectrum));
            }
            Err(e) => return Err(failed(op, format!("requires a diagonalizable matrix: {e}"))),
        };
        let sd = d.map_indexed(|i, j, e| {
            if i == j {
                principal_sqrt(e, real_spectrum)
            } else {
                e.clone()
            }
        });
        let p_inv = match p.inv() {
            Ok(pi) => pi,
            Err(e) if is_swell(&e) => {
                return self.sylvester_or(op, e, &|l| principal_sqrt(l, real_spectrum));
            }
            Err(e) => return Err(e),
        };
        Ok(p.matmul(&sd)?.matmul(&p_inv)?.map(simplify_function_entry))
    }

    /// Principal matrix logarithm `log A` via the Jordan decomposition
    /// `A = P·J·P⁻¹`: `log A = P·log(J)·P⁻¹`.  SymPy: `Matrix.log()`.
    ///
    /// For a Jordan block `J_k(λ) = λI + N`,
    /// `log J = ln(λ)·I + Σ_{d=1}^{k−1} (−1)^{d+1} Nᵈ / (d·λᵈ)`, i.e. the
    /// entry `d` places above the diagonal is `(−1)^{d+1} / (d·λᵈ)`.  For
    /// diagonalizable matrices this is `P·diag(ln λᵢ)·P⁻¹`.  Each
    /// eigenvalue uses the principal branch of `ln`, so `matrix_exp` of
    /// the result reproduces `A` but `matrix_log(matrix_exp(A))` need not
    /// reproduce `A` when eigenvalues of `A` have imaginary part outside
    /// `(−π, π]`.  Entries are simplified.
    ///
    /// # Errors
    ///
    /// - [`SymplexError::InvalidArgument`] if the matrix is not square.
    /// - [`SymplexError::ComputationFailed`] if the matrix is singular
    ///   (an eigenvalue is provably zero: the logarithm does not exist),
    ///   or the Jordan form cannot be computed (see
    ///   [`jordan_form`](Matrix::jordan_form)).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// // SymPy: Matrix([[2, 0], [0, 3]]).log() == [[log(2), 0], [0, log(3)]]
    /// let l = matrix![ctx, [2, 0], [0, 3]].matrix_log().unwrap();
    /// assert_eq!(l, Matrix::diag(&[ctx.int(2).ln(), ctx.int(3).ln()]).unwrap());
    /// // Defective: log [[1, 1], [0, 1]] = [[0, 1], [0, 0]]
    /// let n = matrix![ctx, [1, 1], [0, 1]];
    /// assert_eq!(n.matrix_log().unwrap(), matrix![ctx, [0, 1], [0, 0]]);
    /// assert!(matrix![ctx, [1, 2], [2, 4]].matrix_log().is_err());   // singular
    /// ```
    pub fn matrix_log(&self) -> Result<Matrix, SymplexError> {
        if !self.is_square() {
            return Err(invalid(
                "matrix_log",
                format!(
                    "requires a square matrix, got {}×{}",
                    self.nrows(),
                    self.ncols()
                ),
            ));
        }
        let real_spectrum = self.is_hermitian() == Some(true);
        match self.matrix_log_jordan(real_spectrum) {
            Err(e) if is_swell(&e) => {
                if ex_is_zero(&self.det()?) == Some(true) {
                    return Err(failed(
                        "matrix_log",
                        "matrix is singular (eigenvalue 0); the logarithm does not exist",
                    ));
                }
                self.sylvester_or("matrix_log", e, &|l| principal_ln(l, real_spectrum))
            }
            r => r,
        }
    }

    /// [`matrix_log`](Self::matrix_log) through the Jordan form.
    fn matrix_log_jordan(&self, real_spectrum: bool) -> Result<Matrix, SymplexError> {
        let n = self.nrows();
        let ctx = self.context();
        let JordanForm { p, j } = self.jordan_form().map_err(|e| match e {
            SymplexError::ComputationFailed { reason, .. }
                if reason.starts_with("expression swell") =>
            {
                failed("matrix_log", reason)
            }
            e => failed("matrix_log", format!("Jordan form unavailable ({e})")),
        })?;

        let zero = ctx.zero();
        let mut log_j: Vec<Vec<Ex>> = vec![vec![zero; n]; n];
        let mut col = 0;
        while col < n {
            let lambda = j[(col, col)].clone();
            let mut block_size = 1;
            while col + block_size < n
                && j[(col + block_size - 1, col + block_size)].is_one_structural()
                && j[(col + block_size, col + block_size)] == lambda
            {
                block_size += 1;
            }
            if ex_is_zero(&lambda) == Some(true) {
                return Err(failed(
                    "matrix_log",
                    "matrix is singular (eigenvalue 0); the logarithm does not exist",
                ));
            }
            let ln_lambda = principal_ln(&lambda, real_spectrum);
            for i in 0..block_size {
                log_j[col + i][col + i] = ln_lambda.clone();
                for d in 1..(block_size - i) {
                    // (−1)^{d+1} / (d · λ^d)
                    let sign = if d % 2 == 1 { 1 } else { -1 };
                    let coef = ctx.rational(sign, d as i64);
                    log_j[col + i][col + i + d] = &coef / &lambda.powi(d as i64);
                }
            }
            col += block_size;
        }
        let log_j = Matrix::new(log_j)?;
        let p_inv = p.inv().map_err(|e| match e {
            SymplexError::ComputationFailed { reason, .. }
                if reason.starts_with("expression swell") =>
            {
                failed("matrix_log", reason)
            }
            _ => failed(
                "matrix_log",
                "eigenvector matrix is singular (internal inconsistency)",
            ),
        })?;
        let result = p.matmul(&log_j)?.matmul(&p_inv)?;
        budget_check(result.iter(), "matrix_log")?;
        Ok(result.map(simplify_function_entry))
    }

    /// Casoratian (discrete Wronskian) of the sequences `f₁(n), …, fₖ(n)`:
    /// `det[ fⱼ(n + i) ]` for `i, j = 0..k`.  A non-zero Casoratian
    /// proves the sequences linearly independent.  SymPy:
    /// `casoratian(seqs, n, zero=False)`; SymPy's default `zero=True`
    /// evaluates at `n = 0`, i.e. `Matrix::casoratian(seqs, n)?.subs(n, 0)`.
    ///
    /// # Errors
    ///
    /// - [`SymplexError::InvalidArgument`] if `seqs` is empty.
    /// - [`SymplexError::ComputationFailed`] if the determinant exceeds
    ///   [`EXPRESSION_BUDGET`](crate::matrix::EXPRESSION_BUDGET).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let n = ctx.symbol("n");
    /// // SymPy: casoratian([2**n, 3**n], n, zero=False) == 6**n, and == 1 at n = 0
    /// let w = Matrix::casoratian(&[ctx.int(2).pow(&n), ctx.int(3).pow(&n)], &n).unwrap();
    /// assert_eq!(w.simplify(), ctx.int(6).pow(&n));
    /// assert_eq!(w.subs(&n, &ctx.int(0)).eval(), ctx.int(1));
    /// // SymPy: casoratian([1, n, n**2], n) == 2
    /// let w2 = Matrix::casoratian(&[ctx.one(), n.clone(), n.powi(2)], &n).unwrap();
    /// assert_eq!(w2.expand(), ctx.int(2));
    /// ```
    pub fn casoratian(seqs: &[Ex], n: &Ex) -> Result<Ex, SymplexError> {
        if seqs.is_empty() {
            return Err(invalid("casoratian", "seqs must be non-empty"));
        }
        let k = seqs.len();
        let ctx = n.context();
        let rows: Vec<Vec<Ex>> = (0..k)
            .map(|i| {
                let shifted = n + &ctx.int(i as i64);
                seqs.iter().map(|f| f.subs(n, &shifted)).collect()
            })
            .collect();
        // k rows of k entries, k ≥ 1.
        Matrix::from_rows_unchecked(rows)
            .det()
            .map_err(|e| reop(e, "casoratian"))
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Calculus: Hessian & Wronskian
// ═══════════════════════════════════════════════════════════════════════════

/// Hessian matrix `H[i][j] = ∂²f / ∂vars[i] ∂vars[j]`.
///
/// # Errors
///
/// Returns [`SymplexError::InvalidArgument`] if `vars` is empty.
///
/// # Examples
///
/// ```
/// use symplex::prelude::*;
/// use symplex::matrix_decomp::hessian;
///
/// let ctx = Context::new();
/// let (x, y) = (ctx.symbol("x"), ctx.symbol("y"));
/// let f = &x.powi(2) * &y + &y.powi(3);
/// let h = hessian(&f, &[&x, &y]).unwrap();
/// assert_eq!(h[(0, 0)], &y * 2);
/// assert_eq!(h[(0, 1)], &x * 2);
/// assert_eq!(h[(1, 1)], &y * 6);
/// assert_eq!(h.is_symmetric(), Some(true));
/// assert!(hessian(&f, &[]).is_err());
/// ```
pub fn hessian(f: &Ex, vars: &[&Ex]) -> Result<Matrix, SymplexError> {
    if vars.is_empty() {
        return Err(invalid("hessian", "vars must be non-empty"));
    }
    let firsts: Vec<Ex> = vars.iter().map(|v| f.diff(v)).collect();
    Matrix::from_fn(vars.len(), vars.len(), |i, j| firsts[i].diff(vars[j]))
}

/// Wronskian `W(f₁, …, fₙ)(x) = det[ fⱼ^(i) ]` of a list of functions.
///
/// Row `i` holds the `i`-th derivatives.  A non-zero Wronskian proves
/// linear independence of the functions.
///
/// Returns NaN when the determinant cannot be formed: `funcs` is empty, or
/// the derivatives exceed [`EXPRESSION_BUDGET`](crate::matrix::EXPRESSION_BUDGET).
/// Use [`try_wronskian`] to get the error instead.
///
/// # Examples
///
/// ```
/// use symplex::prelude::*;
/// use symplex::matrix_decomp::wronskian;
///
/// let ctx = Context::new();
/// let x = ctx.symbol("x");
/// // W(sin x, cos x) = −sin² − cos² = −1
/// let w = wronskian(&[&x.sin(), &x.cos()], &x);
/// assert_eq!(w.simplify(), ctx.int(-1));
/// // W(x, x²) = x²
/// let w2 = wronskian(&[&x, &x.powi(2)], &x);
/// assert_eq!(w2.simplify(), x.powi(2));
/// ```
pub fn wronskian(funcs: &[&Ex], var: &Ex) -> Ex {
    try_wronskian(funcs, var).unwrap_or_else(|_| var.context().nan())
}

/// Wronskian `W(f₁, …, fₙ)(x)` of a list of functions; see [`wronskian`].
///
/// # Errors
///
/// - [`SymplexError::InvalidArgument`] if `funcs` is empty.
/// - [`SymplexError::ComputationFailed`] if the derivatives exceed
///   [`EXPRESSION_BUDGET`](crate::matrix::EXPRESSION_BUDGET).
///
/// # Examples
///
/// ```
/// use symplex::prelude::*;
/// use symplex::matrix_decomp::try_wronskian;
///
/// let ctx = Context::new();
/// let x = ctx.symbol("x");
/// let w = try_wronskian(&[&x.exp(), &(&x * 2).exp()], &x).unwrap();
/// assert_eq!(w.simplify(), (&x * 3).exp());
/// assert!(try_wronskian(&[], &x).is_err());
/// ```
pub fn try_wronskian(funcs: &[&Ex], var: &Ex) -> Result<Ex, SymplexError> {
    if funcs.is_empty() {
        return Err(invalid("wronskian", "funcs must be non-empty"));
    }
    let n = funcs.len();
    let mut rows: Vec<Vec<Ex>> = Vec::with_capacity(n);
    let mut current: Vec<Ex> = funcs.iter().map(|f| (*f).clone()).collect();
    for i in 0..n {
        if i > 0 {
            current = current.iter().map(|f| f.diff(var)).collect();
        }
        rows.push(current.clone());
    }
    // n rows of n entries, n ≥ 1.
    Matrix::from_rows_unchecked(rows).det()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::api::context::Context;
    use crate::base::assumptions::Assumption;

    fn ctxi(ctx: &Context, rows: &[&[i64]]) -> Matrix {
        Matrix::from_i64(ctx, rows).unwrap()
    }

    #[test]
    fn qr_reconstructs_and_is_orthonormal() {
        let ctx = Context::new();
        let a = ctxi(&ctx, &[&[1, 2], &[3, 4], &[5, 6]]);
        let Qr { q, r } = a.qr().unwrap();
        assert_eq!(q.shape(), (3, 2));
        assert_eq!(r.shape(), (2, 2));
        assert_eq!((&q * &r).simplify(), a);
        assert_eq!(
            (&q.transpose() * &q).simplify(),
            Matrix::identity(&ctx, 2).unwrap()
        );
        assert_eq!(r.is_upper_triangular(), Some(true));
        assert_eq!(ex_is_positive(r.get(0, 0)), Some(true));
    }

    #[test]
    fn qr_of_dependent_columns_is_reduced() {
        // SymPy: Matrix([[1, 2], [2, 4]]).QRdecomposition() has a 2×1 Q and
        // R = [[sqrt(5), 2*sqrt(5)]]; only the zero matrix is refused.
        let ctx = Context::new();
        let a = ctxi(&ctx, &[&[1, 2], &[2, 4]]);
        let Qr { q, r } = a.qr().unwrap();
        assert_eq!((q.shape(), r.shape()), ((2, 1), (1, 2)));
        assert_eq!((&q * &r).simplify(), a);
        assert!(ctxi(&ctx, &[&[0, 0], &[0, 0]]).qr().is_err());
    }

    #[test]
    fn gram_schmidt_unnormalized() {
        let ctx = Context::new();
        let v1 = ctxi(&ctx, &[&[1], &[1], &[0]]);
        let v2 = ctxi(&ctx, &[&[1], &[0], &[1]]);
        let g = gram_schmidt(&[v1.clone(), v2], false).unwrap();
        assert_eq!(g[0], v1);
        assert_eq!(
            crate::domains::matrix::dot(&g[0], &g[1]).unwrap().eval(),
            ctx.int(0)
        );
        assert!(gram_schmidt(&[v1.clone(), &v1 * 2], false).is_err());
        assert!(gram_schmidt(&[], false).is_err());
    }

    #[test]
    fn cholesky_and_ldl() {
        let ctx = Context::new();
        let a = ctxi(&ctx, &[&[4, 12, -16], &[12, 37, -43], &[-16, -43, 98]]);
        let l = a.cholesky().unwrap();
        assert_eq!(l, ctxi(&ctx, &[&[2, 0, 0], &[6, 1, 0], &[-8, 5, 3]]));
        let Ldl { l: l2, d } = a.ldl().unwrap();
        assert_eq!((&(&l2 * &d) * &l2.transpose()).eval(), a);
        assert_eq!(d.diagonal(), vec![ctx.int(4), ctx.int(1), ctx.int(9)]);
    }

    #[test]
    fn cholesky_errors() {
        let ctx = Context::new();
        assert!(ctxi(&ctx, &[&[-1, 0], &[0, 1]]).cholesky().is_err());
        assert!(ctxi(&ctx, &[&[1, 2], &[3, 4]]).cholesky().is_err()); // not symmetric
        assert!(ctxi(&ctx, &[&[1, 2, 3]]).cholesky().is_err());
        let x = ctx.symbol("x");
        let sym = Matrix::new(vec![
            vec![x.clone(), ctx.int(0)],
            vec![ctx.int(0), ctx.int(1)],
        ])
        .unwrap();
        assert!(sym.cholesky().is_err(), "undecidable pivot must be Err");
        // With a positivity assumption the symbolic factorization works.
        let p = ctx.symbol_with("p", &[Assumption::Positive]).unwrap();
        let sym_p = Matrix::new(vec![
            vec![p.clone(), ctx.int(0)],
            vec![ctx.int(0), ctx.int(1)],
        ])
        .unwrap();
        let l = sym_p.cholesky().unwrap();
        assert_eq!(l.get(0, 0), &p.sqrt());
    }

    #[test]
    fn ldl_symbolic_indefinite() {
        let ctx = Context::new();
        let (a, b, c) = (ctx.symbol("a"), ctx.symbol("b"), ctx.symbol("c"));
        let m = Matrix::new(vec![vec![a.clone(), b.clone()], vec![b.clone(), c.clone()]]).unwrap();
        let Ldl { l, d } = m.ldl().unwrap();
        let back = (&(&l * &d) * &l.transpose()).simplify();
        assert_eq!(back.equals(&m), Some(true));
        assert!(ctxi(&ctx, &[&[0, 1], &[1, 0]]).ldl().is_err());
    }

    #[test]
    fn structure_predicates() {
        let ctx = Context::new();
        let x = ctx.symbol("x");
        let sym = ctxi(&ctx, &[&[1, 2], &[2, 1]]);
        assert_eq!(sym.is_symmetric(), Some(true));
        assert_eq!(ctxi(&ctx, &[&[1, 2], &[3, 4]]).is_symmetric(), Some(false));
        assert_eq!(ctxi(&ctx, &[&[1, 2, 3]]).is_symmetric(), Some(false));
        let trig = Matrix::new(vec![
            vec![ctx.int(1), &x.sin().powi(2) + &x.cos().powi(2)],
            vec![ctx.int(1), ctx.int(0)],
        ])
        .unwrap();
        assert_eq!(trig.is_symmetric(), Some(true));
        let unknown = Matrix::new(vec![
            vec![ctx.int(1), x.clone()],
            vec![ctx.int(1), ctx.int(0)],
        ])
        .unwrap();
        assert_eq!(unknown.is_symmetric(), None);

        assert_eq!(
            ctxi(&ctx, &[&[0, 1], &[-1, 0]]).is_skew_symmetric(),
            Some(true)
        );
        assert_eq!(sym.is_skew_symmetric(), Some(false));
        assert_eq!(sym.is_hermitian(), Some(true));
        let i = ctx.i_unit();
        let herm = Matrix::new(vec![vec![ctx.int(1), i.clone()], vec![-&i, ctx.int(2)]]).unwrap();
        assert_eq!(herm.is_hermitian(), Some(true));
        assert_eq!(herm.is_symmetric(), Some(false));

        let rot = Matrix::new(vec![vec![x.cos(), -&x.sin()], vec![x.sin(), x.cos()]]).unwrap();
        assert_eq!(rot.is_orthogonal(), Some(true));
        assert_eq!(ctxi(&ctx, &[&[1, 1], &[0, 1]]).is_orthogonal(), Some(false));
        let u = Matrix::new(vec![
            vec![ctx.int(0), i.clone()],
            vec![i.clone(), ctx.int(0)],
        ])
        .unwrap();
        assert_eq!(u.is_unitary(), Some(true));

        let up = ctxi(&ctx, &[&[1, 2], &[0, 3]]);
        assert_eq!(up.is_upper_triangular(), Some(true));
        assert_eq!(up.is_lower_triangular(), Some(false));
        assert_eq!(up.transpose().is_lower_triangular(), Some(true));
        assert_eq!(Matrix::identity(&ctx, 3).unwrap().is_diagonal(), Some(true));
        assert_eq!(Matrix::identity(&ctx, 3).unwrap().is_identity(), Some(true));
        assert_eq!(up.is_identity(), Some(false));
        assert_eq!(Matrix::zeros(&ctx, 2, 3).unwrap().is_zero(), Some(true));
        assert_eq!(up.is_zero(), Some(false));
        assert_eq!(ctxi(&ctx, &[&[0, 1], &[0, 0]]).is_nilpotent(), Some(true));
        assert_eq!(up.is_nilpotent(), Some(false));
        assert_eq!(ctxi(&ctx, &[&[1, 2, 3]]).is_nilpotent(), Some(false));
    }

    #[test]
    fn definiteness() {
        let ctx = Context::new();
        assert_eq!(
            ctxi(&ctx, &[&[2, -1, 0], &[-1, 2, -1], &[0, -1, 2]]).is_positive_definite(),
            Some(true)
        );
        assert_eq!(
            ctxi(&ctx, &[&[1, 2], &[2, 1]]).is_positive_definite(),
            Some(false)
        );
        assert_eq!(
            ctxi(&ctx, &[&[1, 2], &[3, 4]]).is_positive_definite(),
            Some(false)
        );
        assert_eq!(
            ctxi(&ctx, &[&[1, 1], &[1, 1]]).is_positive_definite(),
            Some(false)
        );
        assert_eq!(
            ctxi(&ctx, &[&[1, 1], &[1, 1]]).is_positive_semidefinite(),
            Some(true)
        );
        // Leading minors are all zero but the matrix is not PSD:
        assert_eq!(
            ctxi(&ctx, &[&[0, 0], &[0, -1]]).is_positive_semidefinite(),
            Some(false)
        );
        let p = ctx.symbol_with("p", &[Assumption::Positive]).unwrap();
        let m = Matrix::new(vec![
            vec![p.clone(), ctx.int(0)],
            vec![ctx.int(0), p.clone()],
        ])
        .unwrap();
        assert_eq!(m.is_positive_definite(), Some(true));
    }

    #[test]
    fn norms() {
        let ctx = Context::new();
        let a = ctxi(&ctx, &[&[1, -2], &[3, 4]]);
        assert_eq!(a.norm_1(), ctx.int(6));
        assert_eq!(a.norm_inf(), ctx.int(7));
        assert_eq!(a.norm_frobenius().eval(), ctx.int(30).sqrt().eval());
        assert_eq!(a.norm(), a.norm_frobenius());
        let v = ctxi(&ctx, &[&[3], &[-4]]);
        assert_eq!(v.norm_p(&ctx.int(2)).unwrap().eval(), ctx.int(5));
        assert!(a.norm_p(&ctx.int(2)).is_err());
    }

    #[test]
    fn symbolic_power_and_sqrt() {
        let ctx = Context::new();
        let n = ctx.symbol("n");
        let a = ctxi(&ctx, &[&[2, 1], &[0, 3]]);
        let an = a.matrix_pow_symbolic(&n).unwrap();
        for k in 0..4u32 {
            let direct = a.powi(k).unwrap();
            let via = an.subs(&n, &ctx.int(k as i64)).eval().simplify();
            assert_eq!(via, direct, "A^{k}");
        }
        let s = ctxi(&ctx, &[&[4, 0], &[0, 9]]).matrix_sqrt().unwrap();
        assert_eq!(s, ctxi(&ctx, &[&[2, 0], &[0, 3]]));
        let spd = ctxi(&ctx, &[&[2, 1], &[1, 2]]);
        let r = spd.matrix_sqrt().unwrap();
        assert_eq!((&r * &r).simplify(), spd);
        // Defective: through the Jordan blocks (SymPy:
        // Matrix([[1, 1], [0, 1]])**(1/2) == [[1, 1/2], [0, 1]]); a nilpotent
        // block has no square root.
        let half = ctx.rational(1, 2);
        let s = ctxi(&ctx, &[&[1, 1], &[0, 1]]).matrix_sqrt().unwrap();
        assert_eq!(
            s,
            Matrix::new(vec![vec![ctx.one(), half], vec![ctx.zero(), ctx.one()]]).unwrap()
        );
        assert!(ctxi(&ctx, &[&[0, 1], &[0, 0]]).matrix_sqrt().is_err());
    }

    #[test]
    fn hessian_and_wronskian() {
        let ctx = Context::new();
        let (x, y) = (ctx.symbol("x"), ctx.symbol("y"));
        let f = &x.powi(3) + &(&x * &y.powi(2));
        let h = hessian(&f, &[&x, &y]).unwrap();
        assert_eq!(h.get(0, 0), &(&x * 6));
        assert_eq!(h.get(0, 1), &(&y * 2));
        assert_eq!(h.get(1, 0), &(&y * 2));
        assert_eq!(h.get(1, 1), &(&x * 2));
        let w = wronskian(&[&x.exp(), &(&x * 2).exp()], &x).simplify();
        // W(e^x, e^{2x}) = 2e^{3x} − e^{3x} = e^{3x}
        assert_eq!(w, (&x * 3).exp());
        // Dependent functions → zero Wronskian
        let w0 = wronskian(&[&x, &(&x * 2)], &x).simplify();
        assert!(w0.is_zero_structural());
    }
}