alkahest_cas/ode/dsolve/

system.rs

//! Linear constant-coefficient **systems** `y' = A·y + f(t)`.
//!
//! [`super::dsolve`] is scalar by construction — its input is one equation in
//! one unknown — so a two-compartment pharmacokinetic model
//! (`x' = −ka·x`, `y' = ka·x − ke·y`) or any state-space plant had no route
//! through it at all.  This module takes the [`ODE`] type the rest of the crate
//! already uses to model `d(state)/dt = rhs` and solves it in closed form when
//! the right-hand side is linear in the state with coefficients free of `t`.
//!
//! # Why a separate entry point rather than an overload of `dsolve`
//!
//! [`super::dsolve`]'s [`OdeInput`](super::OdeInput) names *one* unknown `y`
//! and a list of derivative symbols `y', y'', …` of that one unknown.  A system
//! has `n` unknowns and one derivative each; there is no reading of `OdeInput`
//! under which the two coincide, and widening it would mean a second
//! interpretation of `derivs` selected by a flag.  [`ODE`] already models
//! exactly the shape needed (`state_vars`, `rhs`, `time_var`), is already
//! public, and is already what the numeric integrators consume — so a system
//! solved here and a system integrated by `ode_integrate_rk45` are the same
//! object, which they would not be under an overload.
//!
//! # Method: Putzer, not Jordan
//!
//! The fundamental matrix is `Φ(t) = e^{At}`, computed by **Putzer's
//! algorithm**:
//!
//! ```text
//! e^{At} = Σ_{k=1}^{n} r_k(t)·P_{k−1},
//!     P_0 = I,  P_k = (A − λ_k I)·P_{k−1},
//!     r_1' = λ_1 r_1,  r_1(0) = 1,
//!     r_{k+1}' = λ_{k+1} r_{k+1} + r_k,  r_{k+1}(0) = 0.
//! ```
//!
//! It needs the eigenvalues **and nothing else** — no eigenvectors, no
//! nullspaces, no similarity transform.  That matters twice over:
//!
//! * A **defective** `A` costs nothing extra.  For `[[2,1],[0,2]]` the
//!   eigenvalue list is `2, 2`, the recurrence gives `r_1 = e^{2t}`,
//!   `r_2 = t·e^{2t}`, and `e^{At} = e^{2t}I + t·e^{2t}(A − 2I)` is
//!   `[[e^{2t}, t·e^{2t}], [0, e^{2t}]]` — the Jordan answer, with no Jordan
//!   machinery and no `jordan_form` refusal to route around.
//! * Every eigenvector routine in the crate has to decide whether a candidate
//!   pivot vanishes, and over symbolic entries that question is undecidable
//!   (see `crate::matrix::zero_test`).  Putzer moves the only such question
//!   to one place — `λ_i − λ_j`, in the recurrence for `r_{k+1}` — where it can
//!   be reported instead of guessed.
//!
//! `r_k` is kept as an explicit `Σ c·t^j·e^{λt}` list rather than being handed
//! to the integration engine, so the recurrence is exact term by term and the
//! `λ_i = λ_j` confluence is a visible branch rather than a failed integral.
//!
//! # The branch problem
//!
//! When `λ_i − λ_j` is neither provably zero nor a non-zero literal — `ka − ke`
//! in the compartment model — the generic formula divides by it.  The result is
//! then **undefined**, not merely non-general, at the confluence, so the
//! condition travels with the answer in
//! [`SystemSolution::side_conditions`] and the confluent limit is named in
//! [`SystemSolution::notes`].  This is a stronger caveat than the scalar
//! module's: there the two-exponential form still solves the equation at the
//! repeated root and only stops spanning, here the expression has a removable
//! singularity that this solver does not remove.
//!
//! # Verification
//!
//! Every component of the returned solution is substituted back into *its own*
//! equation and required to reduce to zero — symbolically, or numerically over
//! ℂ at sampled `t`, integration constants and free parameters, with the
//! finite / non-finite / unevaluable classification the scalar gate uses.  A
//! candidate that does not verify is withheld.

// Matrix and per-equation loops index several parallel structures at once
// (`A`, `P`, the candidate, its derivative and the system's own right-hand
// side), so an index is the clearer spelling.  `matrix::eigen` carries the same
// allow for the same reason.
#![allow(clippy::needless_range_loop)]

use super::verify::{eval_complex, C64, CONST_SETS, PARAM_SETS, X_SAMPLES};
use super::{contains, ddx, integrate_or_decline, simp, simp_plain, sub, AssumedSign, SolveCtx};
use crate::deriv::SideCondition;
use crate::eval::symbols::collect_free_symbols;
use crate::kernel::eval_const::try_expr_f64;
use crate::kernel::{ExprData, ExprId, ExprPool};
use crate::matrix::putzer::{self, Gap, GapPolicy};
use crate::matrix::zero_test::{zero_status, ZeroStatus};
use crate::matrix::Matrix;
use crate::ode::ODE;
use crate::simplify::assumptions::AssumptionContext;
use std::collections::{HashMap, HashSet};
use std::fmt;

// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------

/// The general solution of a linear constant-coefficient system.
#[derive(Clone, Debug)]
pub struct SystemSolution {
    /// The state variables, in the order the [`ODE`] listed them.
    pub state_vars: Vec<ExprId>,
    /// `y_of_t[i]` is the closed form of `state_vars[i]` as a function of `t`.
    pub y_of_t: Vec<ExprId>,
    /// The fresh constants `C1 … Cn` appearing in [`Self::y_of_t`].
    pub constants: Vec<ExprId>,
    /// The fundamental matrix `e^{At}`, row-major, `n × n`.
    pub fundamental_matrix: Vec<Vec<ExprId>>,
    /// Short label of the route taken.
    pub method: &'static str,
    /// Conditions under which [`Self::y_of_t`] is defined and general.
    pub side_conditions: Vec<SideCondition>,
    /// Prose for the cases [`Self::side_conditions`] excludes.
    pub notes: Vec<String>,
}

/// Errors / declines from [`dsolve_system`].
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum DsolveSystemError {
    /// A right-hand side is not affine in the state variables.
    NotLinear(String),
    /// A coefficient of the state depends on the independent variable, so the
    /// system is `y' = A(t)·y` and `e^{At}` is not its fundamental matrix.
    NonConstantCoefficient(String),
    /// The eigenvalues of `A` are not available in closed form.
    UnsupportedSpectrum(String),
    /// The forcing term could not be closed (a required integral is not
    /// elementary), or some other step declined.
    Unsupported(String),
    /// A candidate was produced and failed the substitution gate, so it is
    /// withheld rather than returned.
    VerificationFailed(String),
}

impl fmt::Display for DsolveSystemError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            DsolveSystemError::NotLinear(m) => {
                write!(f, "dsolve_system: system is not linear in the state: {m}")
            }
            DsolveSystemError::NonConstantCoefficient(m) => write!(
                f,
                "dsolve_system: coefficient depends on the independent variable: {m}"
            ),
            DsolveSystemError::UnsupportedSpectrum(m) => {
                write!(f, "dsolve_system: no closed-form spectrum: {m}")
            }
            DsolveSystemError::Unsupported(m) => write!(f, "dsolve_system: unsupported: {m}"),
            DsolveSystemError::VerificationFailed(m) => {
                write!(f, "dsolve_system: candidate failed verification: {m}")
            }
        }
    }
}

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

impl crate::errors::AlkahestError for DsolveSystemError {
    fn code(&self) -> &'static str {
        match self {
            DsolveSystemError::NotLinear(_) => "E-ODE-030",
            DsolveSystemError::NonConstantCoefficient(_) => "E-ODE-031",
            DsolveSystemError::UnsupportedSpectrum(_) => "E-ODE-032",
            DsolveSystemError::Unsupported(_) => "E-ODE-033",
            DsolveSystemError::VerificationFailed(_) => "E-ODE-034",
        }
    }

    fn remediation(&self) -> Option<&'static str> {
        match self {
            DsolveSystemError::NotLinear(_) => Some(
                "each right-hand side must be affine in the state variables; \
                 nonlinear systems have no closed-form matrix exponential",
            ),
            DsolveSystemError::NonConstantCoefficient(_) => Some(
                "y' = A(t)·y needs a Magnus/Peano series, not e^{At}; freeze the \
                 coefficients or use a numeric integrator",
            ),
            DsolveSystemError::UnsupportedSpectrum(_) => Some(
                "the characteristic polynomial of A has no closed-form roots here; \
                 substitute concrete values for the symbolic entries, or reduce the \
                 system to 2×2/3×3 or a triangular shape",
            ),
            DsolveSystemError::Unsupported(_) => Some(
                "a required integral of the forcing term did not close in elementary \
                 form; no unevaluated-integral output is returned",
            ),
            DsolveSystemError::VerificationFailed(_) => Some(
                "the solver found a candidate that did not verify by substitution; \
                 this is reported rather than returned as a (possibly wrong) answer",
            ),
        }
    }
}

/// Solve `y' = A·y + f(t)` in closed form.
///
/// See the [module docs](self) for the method and for how a parameter-dependent
/// eigenvalue confluence is reported.
///
/// # Errors
///
/// See [`DsolveSystemError`].  Every returned solution has been verified by
/// substitution into every equation of the system.
pub fn dsolve_system(ode: &ODE, pool: &ExprPool) -> Result<SystemSolution, DsolveSystemError> {
    dsolve_system_with(ode, &AssumptionContext::new(), pool)
}

/// [`dsolve_system`], consulting caller-stated facts before assuming a branch.
///
/// A stated `λ_i − λ_j ≠ 0` (or any fact pinning it up to a non-zero constant
/// multiple) removes the corresponding side condition.
///
/// # Errors
///
/// See [`DsolveSystemError`].
pub fn dsolve_system_with(
    ode: &ODE,
    assumptions: &AssumptionContext,
    pool: &ExprPool,
) -> Result<SystemSolution, DsolveSystemError> {
    let mut ctx = SolveCtx {
        assumptions,
        conds: super::Conditions::default(),
    };
    solve_inner(ode, &mut ctx, pool)
}

// ---------------------------------------------------------------------------
// Extraction: y' = A·y + f(t)
// ---------------------------------------------------------------------------

/// Split each `rhs_i` into `Σ_j A[i][j]·y_j + f_i(t)`.
fn extract_linear_system(
    ode: &ODE,
    pool: &ExprPool,
) -> Result<(Matrix, Vec<ExprId>), DsolveSystemError> {
    let n = ode.state_vars.len();
    let t = ode.time_var;
    let mut rows: Vec<Vec<ExprId>> = Vec::with_capacity(n);
    let mut forcing: Vec<ExprId> = Vec::with_capacity(n);

    for &r in &ode.rhs {
        let mut row = Vec::with_capacity(n);
        for &v in &ode.state_vars {
            let c = ddx(r, v, pool).map_err(|e| DsolveSystemError::NotLinear(e.to_string()))?;
            // A coefficient still mentioning a state variable means the row is
            // not affine in the state.
            for &w in &ode.state_vars {
                if contains(c, w, pool) {
                    return Err(DsolveSystemError::NotLinear(format!(
                        "∂({})/∂({}) still depends on {}",
                        pool.display(r),
                        pool.display(v),
                        pool.display(w)
                    )));
                }
            }
            if contains(c, t, pool) {
                return Err(DsolveSystemError::NonConstantCoefficient(format!(
                    "∂({})/∂({}) = {} depends on {}",
                    pool.display(r),
                    pool.display(v),
                    pool.display(c),
                    pool.display(t)
                )));
            }
            row.push(simp(c, pool));
        }
        // f = rhs − Σ A[i][j]·y_j
        let mut acc = r;
        for (c, &v) in row.iter().zip(ode.state_vars.iter()) {
            acc = sub(acc, simp(pool.mul(vec![*c, v]), pool), pool);
        }
        for &w in &ode.state_vars {
            if contains(acc, w, pool) {
                return Err(DsolveSystemError::NotLinear(format!(
                    "the state-independent remainder {} still depends on {}",
                    pool.display(acc),
                    pool.display(w)
                )));
            }
        }
        rows.push(row);
        forcing.push(simp(acc, pool));
    }
    let a = Matrix::new(rows).map_err(|e| DsolveSystemError::NotLinear(e.to_string()))?;
    Ok((a, forcing))
}

// ---------------------------------------------------------------------------
// Spectrum
// ---------------------------------------------------------------------------

/// The `n` eigenvalues of `A`, listed with multiplicity.
///
/// A triangular `A` is answered from its diagonal directly.  That is not only a
/// shortcut: for a compartment model — lower triangular by construction, and
/// the case this module exists for — the general routine would build
/// `det(λI − A)`, fail to clear it to ℤ\[λ\] because the entries are symbolic,
/// and fall through to the Cardano formula, returning three nested-radical
/// expressions for what the diagonal states in three symbols.
fn spectrum(a: &Matrix, pool: &ExprPool) -> Result<Vec<ExprId>, DsolveSystemError> {
    let n = a.rows;
    if a.cols != n {
        return Err(DsolveSystemError::NotLinear(
            "coefficient matrix is not square".to_string(),
        ));
    }
    if is_triangular(a, pool) {
        return Ok((0..n).map(|i| simp(a.get(i, i), pool)).collect());
    }
    // A symbolic, non-triangular 3×3 has a cubic characteristic polynomial with
    // free-symbol coefficients, and the only closed form for its roots is
    // Cardano's.  Those radicals are correct only on a coordinated choice of
    // three cube-root branches, they are not what a caller wants to read, and
    // `simplify` grows rather than shrinks them: the resulting `e^{At}` is large
    // enough that building it, let alone verifying it, does not terminate in any
    // useful time.  Refusing here is cheaper and more honest than timing out and
    // then declining anyway.
    if n >= 3 && a.entries().iter().any(|&e| try_expr_f64(e, pool).is_none()) {
        return Err(DsolveSystemError::UnsupportedSpectrum(format!(
            "a {n}×{n} coefficient matrix with symbolic entries that is not triangular \
             needs Cardano (or higher) radicals for its spectrum; those do not simplify, \
             do not verify, and are not returned.  Substitute concrete values for the \
             symbolic entries, or write the system in triangular (compartment) form"
        )));
    }
    let eigs = crate::matrix::eigenvalues(a, pool)
        .map_err(|e| DsolveSystemError::UnsupportedSpectrum(e.to_string()))?;
    let mut out = Vec::with_capacity(n);
    for (lam, mult) in eigs {
        for _ in 0..mult {
            out.push(simp(lam, pool));
        }
    }
    if out.len() != n {
        return Err(DsolveSystemError::UnsupportedSpectrum(format!(
            "eigenvalue multiplicities sum to {} for an {n}×{n} matrix; the \
             characteristic polynomial did not split",
            out.len()
        )));
    }
    Ok(out)
}

/// Cross-check the eigenvalue list before anything is built on it.
///
/// Delegates to `matrix::spectrum::confirm_spectrum`, which is the
/// canonical form of this check and lives beside the routine that produces the
/// spectrum — one copy, checked by one set of tests, rather than a second that
/// drifts.  What stays here is the *policy*.
///
/// Putzer's algorithm turns a wrong eigenvalue list into a page-long `e^{At}`
/// candidate that the substitution gate then refuses — correctly, but slowly,
/// and under an error saying "verification failed" when the truth is "these are
/// not the eigenvalues".  So a system is solved only on a spectrum that was
/// positively **confirmed**, and an unevaluable one is refused here.
/// `matrix::eigenvalues` is deliberately laxer: it has no expansion to protect,
/// and refusing an unevaluable spectrum there would discard every correct
/// symbolic answer whose radicals the numeric evaluator happens not to model.
fn confirm_spectrum(
    a: &Matrix,
    lambdas: &[ExprId],
    pool: &ExprPool,
) -> Result<(), DsolveSystemError> {
    match crate::matrix::spectrum::confirm_spectrum(a, lambdas, pool) {
        Ok(check) if check.is_confirmed() => Ok(()),
        Ok(_) => Err(DsolveSystemError::UnsupportedSpectrum(
            "the eigenvalues of the coefficient matrix could not be evaluated at any \
             sample, so nothing confirms they are its spectrum"
                .to_string(),
        )),
        Err(refusal) => Err(DsolveSystemError::UnsupportedSpectrum(refusal.to_string())),
    }
}

fn is_triangular(a: &Matrix, pool: &ExprPool) -> bool {
    let n = a.rows;
    let upper =
        (0..n).all(|i| (0..i).all(|j| matches!(zero_status(pool, a.get(i, j)), ZeroStatus::Zero)));
    let lower = (0..n)
        .all(|i| ((i + 1)..n).all(|j| matches!(zero_status(pool, a.get(i, j)), ZeroStatus::Zero)));
    upper || lower
}

// ---------------------------------------------------------------------------
// Putzer
// ---------------------------------------------------------------------------

/// This module's policy for the one undecidable question Putzer asks.
///
/// A [`SystemSolution`] has a `side_conditions` field, so an undecided
/// `λ_i − λ_j` is a *hypothesis carried with the answer*: the two-compartment
/// model's `ka ≠ ke` is exactly this, together with a note naming the confluent
/// limit.  [`crate::matrix::matrix_exponential`] divides by the same gaps but
/// returns a bare `Matrix`, so it reports them out of band through
/// [`putzer::RecordedGaps`] instead — same algorithm, different place to put
/// the hypothesis, which is why the policy is a parameter.
struct DsolveGaps<'a, 'b> {
    ctx: &'a mut SolveCtx<'b>,
}

impl GapPolicy for DsolveGaps<'_, '_> {
    fn classify(&mut self, lambda: ExprId, mu: ExprId, pool: &ExprPool) -> Gap {
        let d = putzer::normalised_gap(lambda, mu, pool);
        if lambda == mu || matches!(zero_status(pool, d), ZeroStatus::Zero) {
            return Gap::Confluent;
        }
        require_nonzero(d, self.ctx, pool);
        Gap::Distinct(d)
    }
}

/// `e^{At}` by Putzer's algorithm, row-major.
///
/// The expansion itself lives in [`crate::matrix::putzer`], shared with
/// [`crate::matrix::matrix_exponential`]: one `e^{At}` in the crate rather than
/// two, after the other one was found to be wrong for every defective matrix.
fn matrix_exponential(
    a: &Matrix,
    lambdas: &[ExprId],
    t: ExprId,
    ctx: &mut SolveCtx<'_>,
    pool: &ExprPool,
) -> Vec<Vec<ExprId>> {
    let n = a.rows;
    let phi = putzer::matrix_exponential(a, lambdas, t, &mut DsolveGaps { ctx }, pool)
        .expect("a system with no state variables is rejected before this point");
    (0..n)
        .map(|i| (0..n).map(|j| phi.get(i, j)).collect())
        .collect()
}

// ---------------------------------------------------------------------------
// Assembly
// ---------------------------------------------------------------------------

fn solve_inner(
    ode: &ODE,
    ctx: &mut SolveCtx<'_>,
    pool: &ExprPool,
) -> Result<SystemSolution, DsolveSystemError> {
    let n = ode.state_vars.len();
    if n == 0 {
        return Err(DsolveSystemError::NotLinear(
            "the system has no state variables".to_string(),
        ));
    }
    if ode.rhs.len() != n {
        return Err(DsolveSystemError::NotLinear(
            "state variable and right-hand-side counts differ".to_string(),
        ));
    }
    let t = ode.time_var;
    let (a, forcing) = extract_linear_system(ode, pool)?;
    let lambdas = spectrum(&a, pool)?;
    confirm_spectrum(&a, &lambdas, pool)?;
    let phi = matrix_exponential(&a, &lambdas, t, ctx, pool);

    // Homogeneous part: y = Φ(t)·C.
    let mut used: HashSet<String> = HashSet::new();
    for &e in ode.rhs.iter().chain(ode.state_vars.iter()) {
        collect_symbol_names(e, pool, &mut used);
    }
    let mut gen = super::ConstGen::with_used(used);
    let constants: Vec<ExprId> = (0..n).map(|_| gen.fresh(pool)).collect();
    let mut y: Vec<ExprId> = (0..n)
        .map(|i| {
            let terms: Vec<ExprId> = (0..n)
                .map(|j| pool.mul(vec![phi[i][j], constants[j]]))
                .collect();
            simp(pool.add(terms), pool)
        })
        .collect();

    // Forcing: y_p = Φ(t)·∫ Φ(t)⁻¹ f(t) dt, using Φ(t)⁻¹ = Φ(−t).
    let forced = forcing.iter().any(|&f| !is_zero(f, pool));
    let method = if forced {
        let phi_inv = negate_time(&phi, t, pool);
        let mut integrated = Vec::with_capacity(n);
        for row in &phi_inv {
            let integrand = simp(
                pool.add(
                    row.iter()
                        .zip(forcing.iter())
                        .map(|(&p, &f)| pool.mul(vec![p, f]))
                        .collect(),
                ),
                pool,
            );
            let g = integrate_or_decline(integrand, t, pool)
                .map_err(|e| DsolveSystemError::Unsupported(e.to_string()))?;
            integrated.push(g);
        }
        for (i, yi) in y.iter_mut().enumerate() {
            let terms: Vec<ExprId> = (0..n)
                .map(|j| pool.mul(vec![phi[i][j], integrated[j]]))
                .collect();
            *yi = simp(pool.add(vec![*yi, pool.add(terms)]), pool);
        }
        "linear_system_putzer_variation_of_parameters"
    } else {
        "linear_system_putzer"
    };

    verify_system(ode, &y, &constants, pool)?;

    Ok(SystemSolution {
        state_vars: ode.state_vars.clone(),
        y_of_t: y,
        constants,
        fundamental_matrix: phi,
        method,
        side_conditions: ctx.conds.side.clone(),
        notes: ctx.conds.notes.clone(),
    })
}

/// `Φ(−t)`: substitute `t → −t` in every entry.
fn negate_time(phi: &[Vec<ExprId>], t: ExprId, pool: &ExprPool) -> Vec<Vec<ExprId>> {
    let neg_t = pool.mul(vec![pool.integer(-1_i32), t]);
    let mut m = HashMap::new();
    m.insert(t, neg_t);
    phi.iter()
        .map(|row| {
            row.iter()
                .map(|&e| simp(crate::kernel::subs::subs(e, &m, pool), pool))
                .collect()
        })
        .collect()
}

fn is_zero(e: ExprId, pool: &ExprPool) -> bool {
    matches!(zero_status(pool, e), ZeroStatus::Zero)
}

fn collect_symbol_names(expr: ExprId, pool: &ExprPool, out: &mut HashSet<String>) {
    pool.with(expr, |d| match d {
        ExprData::Symbol { name, .. } => {
            out.insert(name.clone());
        }
        ExprData::Add(args) | ExprData::Mul(args) | ExprData::Func { args, .. } => {
            for &a in args {
                collect_symbol_names(a, pool, out);
            }
        }
        ExprData::Pow { base, exp } => {
            collect_symbol_names(*base, pool, out);
            collect_symbol_names(*exp, pool, out);
        }
        _ => {}
    });
}

/// Record `e ≠ 0` and the confluence it excludes, unless already settled.
fn require_nonzero(e: ExprId, ctx: &mut SolveCtx<'_>, pool: &ExprPool) {
    if super::is_settled_nonzero(e, pool) || ctx.assumed_sign(e, pool) != AssumedSign::Unknown {
        return;
    }
    ctx.conds.require_nonzero(e);
    ctx.conds.note(format!(
        "two eigenvalues of the coefficient matrix differ by {}; where that vanishes \
         the returned expression divides by zero.  Its limit there exists and is the \
         confluent (t·e^{{λt}}) form, but this solver does not take it — pass concrete \
         values, or assert the difference is non-zero, to get an answer valid at that \
         point",
        pool.display(e)
    ));
}

// ---------------------------------------------------------------------------
// Verification
// ---------------------------------------------------------------------------

/// Absolute floor of the "this residual is zero" band.
const ZERO_TOL: f64 = 1e-9;
/// Relative part of the band; see [`super::verify`] on why it is relative.
const REL_TOL: f64 = 1e-7;
/// Resolved, agreeing `(sample, equation)` pairs required for a numeric
/// certificate.
const MIN_AGREEING: usize = 12;

/// What one `(parameters, constants, t)` sample of the candidate produced.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum SampleOutcome {
    /// Every component and every derivative is a finite complex number.
    Resolved,
    /// Every component evaluated, and at least one is not finite.
    CandidateNonFinite,
    /// Some component contains a construct the evaluator does not know.
    Unevaluable,
}

/// Is the system's own right-hand side finite at this `t` and these parameters,
/// candidate aside?  Probes several finite states, as the scalar gate's
/// `ode_is_regular_at` does; one finite probe is enough, since the question is
/// whether the *equations* have a singularity here.
fn system_is_regular_at(ode: &ODE, env: &HashMap<ExprId, C64>, pool: &ExprPool) -> bool {
    const PROBES: [f64; 3] = [1.0, 2.5, -1.5];
    PROBES.iter().any(|&p| {
        let mut e = env.clone();
        for (k, &v) in ode.state_vars.iter().enumerate() {
            e.insert(v, C64::real(p + 0.25 * k as f64));
        }
        ode.rhs
            .iter()
            .all(|&r| matches!(eval_complex(r, &e, pool), Some(v) if v.is_finite()))
    })
}

/// Substitute the candidate into **every** equation and require zero.
fn verify_system(
    ode: &ODE,
    y: &[ExprId],
    constants: &[ExprId],
    pool: &ExprPool,
) -> Result<(), DsolveSystemError> {
    let t = ode.time_var;
    let mut subst: HashMap<ExprId, ExprId> = HashMap::new();
    for (&v, &sol) in ode.state_vars.iter().zip(y.iter()) {
        subst.insert(v, sol);
    }
    let mut residuals = Vec::with_capacity(y.len());
    let mut derivs = Vec::with_capacity(y.len());
    for (i, &sol) in y.iter().enumerate() {
        let dy = ddx(sol, t, pool).map_err(|e| DsolveSystemError::Unsupported(e.to_string()))?;
        derivs.push(dy);
        let rhs = simp(crate::kernel::subs::subs(ode.rhs[i], &subst, pool), pool);
        residuals.push(simp(sub(dy, rhs, pool), pool));
    }
    if residuals
        .iter()
        .all(|&r| is_literal_zero(r, pool) || is_literal_zero(simp_plain(r, pool), pool))
    {
        return Ok(());
    }

    // Numeric fallback over ℂ: `t`, the constants, and any free parameters.
    let mut bound: Vec<ExprId> = vec![t];
    bound.extend_from_slice(&ode.state_vars);
    bound.extend_from_slice(constants);
    // `pi` and the imaginary unit are symbols in this crate but they are not
    // parameters: `collect_free_symbols` leaves them out and `eval_complex`
    // gives them their own values.  Sampling `pi` at 1.7 is what refused the
    // casus-irreducibilis companion matrix, whose eigenvalues are the correct
    // `2√(−p/3)·cos((acos c + 2πk)/3)`.
    let mut params: Vec<ExprId> = Vec::new();
    for &r in &residuals {
        collect_free_symbols(r, pool, &mut params);
    }
    params.retain(|s| !bound.contains(s));
    params.sort_by_key(|&s| pool.display(s).to_string());

    // The scalar gate's tables, not a second copy of them: a parameter must be
    // sampled the same way whichever gate is asking, and the sign coverage of
    // `PARAM_SETS` is a soundness property that has to hold in both.
    let param_sets = PARAM_SETS;
    let const_sets = CONST_SETS;
    let t_samples = X_SAMPLES;

    let (mut agree, mut disagree, mut skipped) = (0usize, 0usize, 0usize);
    let mut sets_agreeing = 0usize;
    for ps in param_sets {
        let mut env: HashMap<ExprId, C64> = HashMap::new();
        for (i, &p) in params.iter().enumerate() {
            env.insert(p, C64::real(ps[i % ps.len()]));
        }
        let before = agree;
        for cs in const_sets {
            for (i, &c) in constants.iter().enumerate() {
                env.insert(c, C64::real(cs[i % cs.len()]));
            }
            for &tv in &t_samples {
                env.insert(t, C64::real(tv));
                // Split evaluation: the candidate on one side, the system's own
                // right-hand side on the other, so a non-finite sample can be
                // attributed rather than guessed at.
                //
                // The two failures are *not* symmetric, and treating them alike
                // is the mistake `integrate::gate` documents.  A construct the
                // evaluator does not know is a property of the expression, not
                // of the point, and carries no information.  A value it does
                // know and that comes back non-finite where the system itself is
                // perfectly regular is evidence the candidate is wrong — a
                // blanket skip there is how a domain hole clears a grid built to
                // catch one.
                let mut yv = Vec::with_capacity(y.len());
                let mut dv = Vec::with_capacity(y.len());
                let mut outcome = SampleOutcome::Resolved;
                for i in 0..y.len() {
                    match (
                        eval_complex(y[i], &env, pool),
                        eval_complex(derivs[i], &env, pool),
                    ) {
                        (Some(a), Some(b)) if a.is_finite() && b.is_finite() => {
                            yv.push(a);
                            dv.push(b);
                        }
                        (Some(_), Some(_)) => {
                            outcome = SampleOutcome::CandidateNonFinite;
                            break;
                        }
                        _ => {
                            outcome = SampleOutcome::Unevaluable;
                            break;
                        }
                    }
                }
                match outcome {
                    SampleOutcome::Resolved => {}
                    SampleOutcome::Unevaluable => {
                        skipped += 1;
                        continue;
                    }
                    SampleOutcome::CandidateNonFinite => {
                        // Is the *system* regular here, candidate aside?  A
                        // constant-coefficient right-hand side is finite at any
                        // finite state unless the forcing term has a pole at
                        // this `t`, which is the one case that carries no
                        // information.
                        if system_is_regular_at(ode, &env, pool) {
                            disagree += 1;
                        } else {
                            skipped += 1;
                        }
                        continue;
                    }
                }
                let mut rhs_env = env.clone();
                for (k, &v) in ode.state_vars.iter().enumerate() {
                    rhs_env.insert(v, yv[k]);
                }
                let scale = yv
                    .iter()
                    .chain(dv.iter())
                    .fold(1.0_f64, |m, v| m.max(v.abs()));
                for i in 0..y.len() {
                    match eval_complex(ode.rhs[i], &rhs_env, pool) {
                        Some(v) if v.is_finite() => {
                            let resid = dv[i].sub(v).abs();
                            if resid <= ZERO_TOL.max(REL_TOL * scale) {
                                agree += 1;
                            } else {
                                disagree += 1;
                            }
                        }
                        // The system itself is not finite at this state — no
                        // information about the candidate.
                        _ => skipped += 1,
                    }
                }
            }
        }
        if agree > before {
            sets_agreeing += 1;
        }
    }
    if disagree == 0 && agree >= MIN_AGREEING && sets_agreeing >= 2 {
        return Ok(());
    }
    let mut shown = pool.display(residuals[0]).to_string();
    // A refused candidate can be pages long; the message is a diagnostic, not
    // the expression itself.
    if shown.len() > 400 {
        shown.truncate(400);
        shown.push_str(" …");
    }
    Err(DsolveSystemError::VerificationFailed(format!(
        "samples: {agree} agreeing, {disagree} disagreeing, {skipped} unresolved \
         ({sets_agreeing} parameter sets agreeing); first residual: {shown}"
    )))
}

fn is_literal_zero(e: ExprId, pool: &ExprPool) -> bool {
    matches!(pool.get(e), ExprData::Integer(n) if n.0 == 0)
}

#[cfg(test)]
mod tests;