alkahest-cas 3.12.0

//! Symbolic Laplace transform `L{f(t)}(s)` and inverse `L⁻¹{F(s)}(t)`.
//!
//! # Forward transform
//!
//! [`laplace_transform`] is a rule/table-based structural recursion over `f(t)`:
//!
//! | `f(t)`                 | `L{f}(s)`                       | rule              |
//! |------------------------|---------------------------------|-------------------|
//! | `c` (const)            | `c/s`                           | constant          |
//! | `t^n` (`n ∈ ℤ₊`)       | `n! / s^{n+1}`                  | power             |
//! | `e^{a t}`              | `1/(s−a)`                       | exponential       |
//! | `sin(b t)`             | `b/(s²+b²)`                     | sine              |
//! | `cos(b t)`             | `s/(s²+b²)`                     | cosine            |
//! | `sinh(b t)`            | `b/(s²−b²)`                     | hyperbolic sine   |
//! | `cosh(b t)`            | `s/(s²−b²)`                     | hyperbolic cosine |
//! | `α·f + β·g`            | `α·L{f} + β·L{g}`               | linearity         |
//! | `e^{a t}·f(t)`         | `F(s−a)`                        | s-shift theorem   |
//! | `t^n·f(t)`             | `(−1)^n F^{(n)}(s)`             | frequency-diff    |
//! | `θ(t−a)·g(t−a)`        | `e^{−a s} G(s)`                | t-shift (`a ≥ 0`, reported) |
//! | `θ(t−a)`               | `e^{−a s}/s`                    | shifted step      |
//! | `δ(t−a)`               | `e^{−a s}`  (`δ(t) ↦ 1`)        | impulse           |
//!
//! The derivative rule `L{f^{(n)}} = sⁿF − sⁿ⁻¹f(0) − … − f^{(n−1)}(0)` is
//! exposed separately as [`laplace_derivative_rule`] for the ODE workflow (it
//! operates on the *symbol* `F = L{f}` plus initial values, since `f` itself is
//! an unknown function).
//!
//! # Inverse transform
//!
//! [`inverse_laplace_transform`] inverts a **proper rational** `F(s) = p(s)/q(s)`
//! by partial fractions (via [`crate::poly::apart`]) and a per-term table:
//!
//! | term in `F(s)`              | `L⁻¹` term                       |
//! |-----------------------------|----------------------------------|
//! | `A/(s−a)^n`                 | `A·t^{n−1} e^{a t}/(n−1)!`        |
//! | `(B s + C)/((s−p)²+ω²)`     | damped sin/cos (`n = 1`, `ω² > 0`) |
//! | `(B s + C)/((s−p)²−κ²)`     | damped sinh/cosh (`n = 1`, `κ² > 0`) |
//! | `(B s + C)/((s−p)²+ω²)²`    | `t`-weighted damped sin/cos      |
//! | `(B s + C)/((s−p)²−κ²)²`    | `t`-weighted damped sinh/cosh    |
//! | `e^{−a s} F(s)`             | `θ(t−a)·(L⁻¹F)(t−a)`             |
//!
//! A leading polynomial part of `F` (improper rational) maps back to derivatives
//! of `δ(t)`, which we decline rather than fabricate.
//!
//! ## Symbolic parameters
//!
//! `F`'s coefficients need not be rational numbers.  `apart` decomposes over
//! `ℚ(params)` when other symbols are present, so the two workhorse forms of
//! applied mathematics invert:
//!
//! | `F(s)` | `L⁻¹` | hypotheses |
//! |---|---|---|
//! | `K/(s² + 2ζωs + ω²)` | `K·e^{−ζωt}·sin(ω√(1−ζ²)·t)/(ω√(1−ζ²))` | `ω ≠ 0`, `ζ ≠ ±1`, `ω²(1−ζ²) > 0` |
//! | `D·ka/((s+ka)(s+ke))` | `D·ka·(e^{−ka t} − e^{−ke t})/(ke − ka)` | `ka ≠ ke` |
//!
//! The hypotheses are not decoration.  At `ka = ke` the Bateman answer is
//! `0/0` and the true inverse is `D·ka·t·e^{−ka t}`; at `ζ > 1` the
//! second-order answer is a *hyperbolic* function that the printed form
//! expresses only through an imaginary `√(1−ζ²)`.  They are returned by
//! [`inverse_laplace_transform_with_conditions`], discharged where the caller's
//! [`crate::simplify::assumptions::AssumptionContext`] proves them
//! ([`inverse_laplace_transform_with_assumptions`]), and — for the historical
//! [`inverse_laplace_transform`] signature — carried out of band on
//! [`super::take_transform_side_conditions`].
//!
//! ## The unilateral shift
//!
//! The second-shift rule is the transform of `g` only for `a ≥ 0`: at `a < 0`
//! the step edge lies before the origin, outside the range `∫₀^∞` sees, and the
//! answer is a different function (`L{θ(t+1)} = 1/s`, not `e^{s}/s`).  A
//! *literal* negative shift refutes the rule and is refused; a *symbolic* one
//! is a hypothesis, recorded as
//! [`crate::deriv::SideCondition::InDomain`]`(a, NonNegative)` and returned by
//! [`laplace_transform_with_conditions`] — or, for the historical signature, on
//! [`super::take_transform_side_conditions`].  A symbol whose static
//! [`crate::kernel::Domain`] is `Positive`/`NonNegative` discharges it.
//!
//! The same fact governs the inverse's delay factor `e^{−a s}`, with one
//! difference: an *advance* `e^{+a s}` with a literal `a > 0` is not the
//! transform of any causal function, so it is declined rather than answered
//! under a hypothesis that cannot hold.
//!
//! # Caveats
//!
//! Both directions are otherwise **formal** — no convergence region is computed
//! and no existence side-conditions are attached.  Unrecognised forms return
//! [`LaplaceError::NoRule`] (forward) / [`LaplaceError::NotInvertible`] (inverse)
//! rather than guessing.

use rug::Integer;

use crate::deriv::SideCondition;
use crate::kernel::{ExprData, ExprId, ExprPool};
use crate::simplify::assumptions::AssumptionContext;

use super::{stash_transform_side_conditions, Genericity};

/// Errors from the Laplace transform routines.
///
/// The two string-carrying variants each cover *two* codes, split by a message
/// tag rather than by a new variant — the enum is public and exhaustive, so a
/// fourth variant would be a major semver break. The split is not cosmetic:
/// [`Self::NoRule`] / [`Self::NotInvertible`] normally mean "this table does
/// not reach that form yet", a fact about the implementation that a caller
/// answers by rewriting the input or waiting for a wider table, whereas
/// `E-TRANSFORM-004` means the *unilateral transform's own hypothesis is
/// refuted* — there is no rule to add, because no causal function has that
/// transform. Build the tagged form only through
/// [`Self::forward_not_causal`] / [`Self::inverse_not_causal`] so the tag and
/// the code cannot drift apart.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum LaplaceError {
    /// No forward rule matched `f(t)` (`E-TRANSFORM-001`), or — when tagged by
    /// [`Self::forward_not_causal`] — the unilateral hypothesis is refuted
    /// (`E-TRANSFORM-004`).
    NoRule(String),
    /// The inverse-transform input is not a form the table can invert
    /// (`E-TRANSFORM-002`), or — when tagged by [`Self::inverse_not_causal`] —
    /// it is not the transform of any causal function (`E-TRANSFORM-004`).
    NotInvertible(String),
    /// The frequency variable `s` and time variable `t` must be distinct
    /// symbols (`E-TRANSFORM-003`).
    SameVariable,
}

impl LaplaceError {
    /// Message tag for the `E-TRANSFORM-004` refusal.
    const CAUSALITY_TAG: &'static str = "unilateral hypothesis refuted: ";

    /// The forward rule's `a ≥ 0` hypothesis is *refuted*, not merely unproven:
    /// the step/impulse edge lies strictly before the origin, so `∫₀^∞` never
    /// sees it and `L{θ(t−a)g(t−a)} = e^{−as}G(s)` is false here.
    pub fn forward_not_causal(detail: impl std::fmt::Display) -> Self {
        LaplaceError::NoRule(format!("{}{detail}", Self::CAUSALITY_TAG))
    }

    /// The input is not the transform of any causal function — an *advance*
    /// factor `e^{+as}` with literal `a > 0`. Nothing to invert to, rather than
    /// something this table cannot reach.
    pub fn inverse_not_causal(detail: impl std::fmt::Display) -> Self {
        LaplaceError::NotInvertible(format!("{}{detail}", Self::CAUSALITY_TAG))
    }

    /// Was this decline a refuted causality hypothesis (`E-TRANSFORM-004`)?
    pub fn is_causality_refutation(&self) -> bool {
        match self {
            LaplaceError::NoRule(m) | LaplaceError::NotInvertible(m) => {
                m.starts_with(Self::CAUSALITY_TAG)
            }
            LaplaceError::SameVariable => false,
        }
    }
}

impl std::fmt::Display for LaplaceError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            // A tagged message already says what went wrong; "no rule for"
            // would misdescribe it as a gap in the table.
            LaplaceError::NoRule(m) if self.is_causality_refutation() => {
                write!(f, "laplace_transform: {m}")
            }
            LaplaceError::NotInvertible(m) if self.is_causality_refutation() => {
                write!(f, "inverse_laplace_transform: {m}")
            }
            LaplaceError::NoRule(m) => {
                write!(f, "laplace_transform: no rule for {m}")
            }
            LaplaceError::NotInvertible(m) => {
                write!(f, "inverse_laplace_transform: cannot invert {m}")
            }
            LaplaceError::SameVariable => write!(
                f,
                "laplace_transform: time and frequency variables must differ"
            ),
        }
    }
}

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

impl crate::errors::AlkahestError for LaplaceError {
    fn code(&self) -> &'static str {
        match self {
            _ if self.is_causality_refutation() => "E-TRANSFORM-004",
            LaplaceError::NoRule(_) => "E-TRANSFORM-001",
            LaplaceError::NotInvertible(_) => "E-TRANSFORM-002",
            LaplaceError::SameVariable => "E-TRANSFORM-003",
        }
    }

    fn remediation(&self) -> Option<&'static str> {
        match self {
            _ if self.is_causality_refutation() => Some(
                "the unilateral transform integrates over t ≥ 0 only, so a step or impulse \
                 edge at a < 0 — or an advance factor e^{+a·s} on the inverse — is outside \
                 what it can see; no wider table fixes this. Shift the edge to a ≥ 0, or \
                 use a bilateral transform",
            ),
            LaplaceError::NoRule(_) => Some(
                "laplace_transform is table-based: reduce f(t) to a sum of constants, \
                 t^n, e^{a·t}, sin/cos/sinh/cosh of b·t, Heaviside/Dirac shifts, and \
                 products of those with e^{a·t} or t^n",
            ),
            LaplaceError::NotInvertible(_) => Some(
                "inverse_laplace_transform inverts a proper rational F(s) whose \
                 denominator factors into poles of degree ≤ 2 (times an optional delay \
                 e^{−a·s}); an improper part would invert to derivatives of δ, which are \
                 declined rather than fabricated",
            ),
            LaplaceError::SameVariable => {
                Some("pass distinct symbols for the time and frequency variables")
            }
        }
    }
}

// ===========================================================================
// Small helpers
// ===========================================================================

fn is_free_of(expr: ExprId, var: ExprId, pool: &ExprPool) -> bool {
    crate::integrate::risch::poly_rde::is_free_of_var(expr, var, pool)
}

/// `n!` as an interned integer.
fn factorial(n: u64, pool: &ExprPool) -> ExprId {
    let mut acc = Integer::from(1);
    for k in 2..=n {
        acc *= Integer::from(k);
    }
    pool.integer(acc)
}

/// Extract `a` from an expression that is `a·var + b` with `a, b` free of `var`,
/// returning `(a, b)`.  Returns `None` when the expression is not affine in
/// `var`.
fn as_affine(expr: ExprId, var: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId)> {
    if expr == var {
        return Some((pool.integer(1_i32), pool.integer(0_i32)));
    }
    if is_free_of(expr, var, pool) {
        return Some((pool.integer(0_i32), expr));
    }
    match pool.get(expr) {
        ExprData::Mul(_) => {
            let (a, b) = as_affine_term(expr, var, pool)?;
            // a single Mul term has no constant part
            if b == pool.integer(0_i32) {
                Some((a, pool.integer(0_i32)))
            } else {
                None
            }
        }
        ExprData::Add(args) => {
            let mut a_acc: Vec<ExprId> = Vec::new();
            let mut b_acc: Vec<ExprId> = Vec::new();
            for arg in args {
                if is_free_of(arg, var, pool) {
                    b_acc.push(arg);
                } else {
                    let (a, b) = as_affine_term(arg, var, pool)?;
                    if b != pool.integer(0_i32) {
                        return None;
                    }
                    a_acc.push(a);
                }
            }
            let a = match a_acc.len() {
                0 => pool.integer(0_i32),
                1 => a_acc[0],
                _ => pool.add(a_acc),
            };
            let b = match b_acc.len() {
                0 => pool.integer(0_i32),
                1 => b_acc[0],
                _ => pool.add(b_acc),
            };
            Some((a, b))
        }
        _ => None,
    }
}

/// Affine analysis of a single (non-Add) term: returns `(coeff, 0)` for a term
/// `coeff·var`, else `None` unless it is free of `var` (`(0, term)`).
fn as_affine_term(expr: ExprId, var: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId)> {
    if expr == var {
        return Some((pool.integer(1_i32), pool.integer(0_i32)));
    }
    if is_free_of(expr, var, pool) {
        return Some((pool.integer(0_i32), expr));
    }
    if let ExprData::Mul(args) = pool.get(expr) {
        let pos = args.iter().position(|&a| a == var)?;
        let others: Vec<ExprId> = args
            .iter()
            .enumerate()
            .filter(|&(i, _)| i != pos)
            .map(|(_, &a)| a)
            .collect();
        if others.iter().all(|&o| is_free_of(o, var, pool)) {
            let coeff = match others.len() {
                0 => pool.integer(1_i32),
                1 => others[0],
                _ => pool.mul(others),
            };
            return Some((coeff, pool.integer(0_i32)));
        }
    }
    None
}

/// Simplify and return the bare `ExprId`.
fn simp(expr: ExprId, pool: &ExprPool) -> ExprId {
    crate::simplify::simplify(expr, pool).value
}

fn neg(expr: ExprId, pool: &ExprPool) -> ExprId {
    pool.mul(vec![pool.integer(-1_i32), expr])
}

fn recip(expr: ExprId, pool: &ExprPool) -> ExprId {
    pool.pow(expr, pool.integer(-1_i32))
}

// ===========================================================================
// Forward transform
// ===========================================================================

/// Compute the Laplace transform `L{f(t)}(s) = ∫₀^∞ f(t) e^{−s t} dt`.
///
/// `t` is the time variable, `s` the frequency variable; both must be distinct
/// symbols.  This is a *formal* transform — see the [module docs](self) for the
/// rule table, caveats, and declines.
///
/// # Errors
///
/// - [`LaplaceError::SameVariable`] if `t == s`.
/// - [`LaplaceError::NoRule`] if no table rule matches `f(t)`.
///
/// # Examples
///
/// ```
/// use alkahest_cas::kernel::{Domain, ExprPool};
/// use alkahest_cas::transform::laplace_transform;
///
/// let pool = ExprPool::new();
/// let t = pool.symbol("t", Domain::Real);
/// let s = pool.symbol("s", Domain::Real);
/// // L{1}(s) = 1/s
/// let one = pool.integer(1_i32);
/// let f = laplace_transform(one, t, s, &pool).unwrap();
/// assert_eq!(pool.display(f).to_string(), pool.display(pool.pow(s, pool.integer(-1_i32))).to_string());
/// ```
pub fn laplace_transform(
    f: ExprId,
    t: ExprId,
    s: ExprId,
    pool: &ExprPool,
) -> Result<ExprId, LaplaceError> {
    if t == s {
        return Err(LaplaceError::SameVariable);
    }
    stash_transform_side_conditions(Vec::new());
    let (out, conds) = laplace_transform_with_conditions(f, t, s, pool)?;
    stash_transform_side_conditions(conds);
    Ok(out)
}

/// [`laplace_transform`] with its genericity hypotheses returned in band rather
/// than through [`super::take_transform_side_conditions`].
///
/// The list is empty for every input whose shifts are literal numbers.  It is
/// non-empty exactly when a table entry was selected on a fact about a symbolic
/// parameter that the input does not settle — the unilateral second-shift rule
/// `L{θ(t−a)·g(t−a)} = e^{−a s}·G(s)`, which is that transform only for
/// `a ≥ 0`.  At `a < 0` the step edge (or the impulse) sits *before* `t = 0`,
/// outside the range the unilateral integral sees, and the true transform is a
/// different function: `L{θ(t+1)} = 1/s`, not `e^{s}/s`.
///
/// A hypothesis a symbol's static [`crate::kernel::Domain`] already proves
/// (`Domain::Positive` / `Domain::NonNegative`) is discharged rather than
/// reported.  A shift that is a *literal* negative number is not a hypothesis
/// at all and is refused by [`LaplaceError::NoRule`].
pub fn laplace_transform_with_conditions(
    f: ExprId,
    t: ExprId,
    s: ExprId,
    pool: &ExprPool,
) -> Result<(ExprId, Vec<SideCondition>), LaplaceError> {
    if t == s {
        return Err(LaplaceError::SameVariable);
    }
    let mut gen = Genericity::new(None);
    let out = laplace_inner(f, t, s, pool, 0, &mut gen)?;
    Ok((simp(out, pool), gen.into_conditions()))
}

const MAX_DEPTH: usize = 32;

fn laplace_inner(
    f: ExprId,
    t: ExprId,
    s: ExprId,
    pool: &ExprPool,
    depth: usize,
    gen: &mut Genericity<'_>,
) -> Result<ExprId, LaplaceError> {
    if depth > MAX_DEPTH {
        return Err(LaplaceError::NoRule("recursion depth exceeded".into()));
    }

    // Constant (free of t): L{c} = c/s.
    if is_free_of(f, t, pool) {
        return Ok(pool.mul(vec![f, recip(s, pool)]));
    }

    // Bare t: L{t} = 1/s².
    if f == t {
        return Ok(recip(pool.pow(s, pool.integer(2_i32)), pool));
    }

    match pool.get(f) {
        // Linearity over sums.
        ExprData::Add(args) => {
            let mut terms = Vec::with_capacity(args.len());
            for a in args {
                terms.push(laplace_inner(a, t, s, pool, depth + 1, gen)?);
            }
            Ok(pool.add(terms))
        }

        // Products: split off the t-free scalar (linearity), then dispatch the
        // remaining t-dependent factor through the structural product rules.
        ExprData::Mul(args) => laplace_mul(&args, t, s, pool, depth, gen),

        // Pure power t^n.
        ExprData::Pow { base, exp } if base == t => {
            if let Some(n) = nonneg_int_exp(exp, pool) {
                // L{t^n} = n! / s^{n+1}
                let fact = factorial(n, pool);
                let denom = pool.pow(s, pool.integer(Integer::from(n + 1)));
                Ok(pool.mul(vec![fact, recip(denom, pool)]))
            } else {
                Err(LaplaceError::NoRule(format!(
                    "t^e with non-integer exponent: {}",
                    pool.display(f)
                )))
            }
        }

        ExprData::Func { name, args } if args.len() == 1 => {
            laplace_func(&name, args[0], t, s, pool, depth, gen)
        }

        _ => Err(LaplaceError::NoRule(pool.display(f).to_string())),
    }
}

/// Non-negative integer value of an exponent `ExprId`, if it is one.
fn nonneg_int_exp(exp: ExprId, pool: &ExprPool) -> Option<u64> {
    if let ExprData::Integer(n) = pool.get(exp) {
        let n = n.0;
        if n >= 0 {
            return n.to_u64();
        }
    }
    None
}

/// Laplace transform of a product `∏ args`.
fn laplace_mul(
    args: &[ExprId],
    t: ExprId,
    s: ExprId,
    pool: &ExprPool,
    depth: usize,
    gen: &mut Genericity<'_>,
) -> Result<ExprId, LaplaceError> {
    // Pull out the constant (t-free) scalar prefactor.
    let (consts, rest): (Vec<ExprId>, Vec<ExprId>) =
        args.iter().partition(|&&a| is_free_of(a, t, pool));
    let scalar = match consts.len() {
        0 => None,
        1 => Some(consts[0]),
        _ => Some(pool.mul(consts.clone())),
    };

    let inner = match rest.len() {
        0 => {
            // Wholly constant — handled by caller, but be safe.
            let c = scalar.unwrap_or_else(|| pool.integer(1_i32));
            return Ok(pool.mul(vec![c, recip(s, pool)]));
        }
        1 => rest[0],
        _ => pool.mul(rest.clone()),
    };

    let transformed = laplace_product_body(inner, t, s, pool, depth, gen)?;
    Ok(match scalar {
        Some(c) => pool.mul(vec![c, transformed]),
        None => transformed,
    })
}

/// Transform a t-dependent product with no constant scalar factor, applying the
/// structural product theorems (frequency-diff `t^n·f`, s-shift `e^{at}·f`,
/// t-shift `θ(t−a)·g(t−a)`).
fn laplace_product_body(
    body: ExprId,
    t: ExprId,
    s: ExprId,
    pool: &ExprPool,
    depth: usize,
    gen: &mut Genericity<'_>,
) -> Result<ExprId, LaplaceError> {
    let factors: Vec<ExprId> = match pool.get(body) {
        ExprData::Mul(a) => a,
        _ => vec![body],
    };

    // (1) e^{a t} · g(t)  →  G(s − a)   [s-shift theorem]
    for (i, &fac) in factors.iter().enumerate() {
        if let Some(a) = match_exp_linear(fac, t, pool) {
            let rest = remove_index(&factors, i, pool);
            let g_transform = laplace_inner(rest, t, s, pool, depth + 1, gen)?;
            let s_minus_a = simp(pool.add(vec![s, neg(a, pool)]), pool);
            return Ok(subs_one(g_transform, s, s_minus_a, pool));
        }
    }

    // (2) θ(t − a) · g(t − a)  →  e^{−a s} G(s)   [t-shift / second shift]
    if let Some(res) = try_time_shift(&factors, t, s, pool, depth, gen)? {
        return Ok(res);
    }

    // (3) t^n · g(t)  →  (−1)^n F^{(n)}(s)   [frequency differentiation]
    for (i, &fac) in factors.iter().enumerate() {
        if let Some(n) = match_t_power(fac, t, pool) {
            let rest = remove_index(&factors, i, pool);
            let mut g_transform = laplace_inner(rest, t, s, pool, depth + 1, gen)?;
            for _ in 0..n {
                g_transform = crate::diff::diff(g_transform, s, pool)
                    .map_err(|_| LaplaceError::NoRule("frequency-diff failed".into()))?
                    .value;
            }
            let sign = if n % 2 == 0 {
                pool.integer(1_i32)
            } else {
                pool.integer(-1_i32)
            };
            return Ok(pool.mul(vec![sign, g_transform]));
        }
    }

    // No product theorem applied.  If `body` was not actually a product (a lone
    // factor, e.g. a bare `cos(b t)` whose t-free scalar was already peeled off
    // by `laplace_mul`), fall back to the structural table.  This cannot recurse
    // forever: `laplace_inner` only re-enters `laplace_product_body` for a `Mul`,
    // and here `body` is not a `Mul`.
    if !matches!(pool.get(body), ExprData::Mul(_)) {
        return laplace_inner(body, t, s, pool, depth + 1, gen);
    }

    Err(LaplaceError::NoRule(pool.display(body).to_string()))
}

/// If `fac` is `e^{a·t}` (or `e^{a·t+b}` collapsed), return `a` (the linear
/// coefficient), with the additive constant folded into the coefficient via
/// `e^{b}` only when `b` is free of `t` — here we require pure `a·t`.
fn match_exp_linear(fac: ExprId, t: ExprId, pool: &ExprPool) -> Option<ExprId> {
    if let ExprData::Func { name, args } = pool.get(fac) {
        if name == "exp" && args.len() == 1 {
            let (a, b) = as_affine(args[0], t, pool)?;
            if b == pool.integer(0_i32) && a != pool.integer(0_i32) {
                return Some(a);
            }
        }
    }
    // exp written as a power e^(…) is normalised to Func("exp", …) in this CAS.
    None
}

/// If `fac` is `t^n` with `n ∈ ℤ₊`, or bare `t`, return `n`.
fn match_t_power(fac: ExprId, t: ExprId, pool: &ExprPool) -> Option<u64> {
    if fac == t {
        return Some(1);
    }
    if let ExprData::Pow { base, exp } = pool.get(fac) {
        if base == t {
            return nonneg_int_exp(exp, pool).filter(|&n| n >= 1);
        }
    }
    None
}

/// Recognise `θ(t − a) · g(t − a)` (any number of g-factors), returning
/// `e^{−a s} · L{g(t)}` when every t-dependent non-Heaviside factor is a
/// function of `(t − a)` with the *same* shift `a ≥ 0` as the Heaviside.
fn try_time_shift(
    factors: &[ExprId],
    t: ExprId,
    s: ExprId,
    pool: &ExprPool,
    depth: usize,
    gen: &mut Genericity<'_>,
) -> Result<Option<ExprId>, LaplaceError> {
    // Find a Heaviside factor and its shift a (from θ(t − a)).
    let mut heaviside_idx = None;
    let mut shift = None;
    for (i, &fac) in factors.iter().enumerate() {
        if let ExprData::Func { name, args } = pool.get(fac) {
            if name == "heaviside" && args.len() == 1 {
                // arg = t − a  ⇒  (coeff 1)·t + (−a)
                if let Some((coeff, b)) = as_affine(args[0], t, pool) {
                    if coeff == pool.integer(1_i32) {
                        heaviside_idx = Some(i);
                        shift = Some(neg(b, pool)); // a = −b
                        break;
                    }
                }
            }
        }
    }
    let (hi, a) = match (heaviside_idx, shift) {
        (Some(hi), Some(a)) => (hi, simp(a, pool)),
        _ => return Ok(None),
    };
    require_nonneg_shift(a, pool, gen)?;

    // The remaining factors form g(t − a).  Substitute u = t + a (i.e. shift the
    // argument back) and transform g(u), then attach e^{−a s}.
    let rest = remove_index(factors, hi, pool);
    let exp_neg_as = pool.func("exp", vec![simp(neg(pool.mul(vec![a, s]), pool), pool)]);

    // g(t − a): replace t by (t + a) to recover g(t), require the result be a
    // genuine function of t (the standard second-shift form).  If `rest` is just
    // the constant 1 (bare shifted step), G(s) = 1/s.
    if rest == pool.integer(1_i32) {
        return Ok(Some(pool.mul(vec![exp_neg_as, recip(s, pool)])));
    }
    let t_plus_a = simp(pool.add(vec![t, a]), pool);
    let g_of_t = subs_one(rest, t, t_plus_a, pool);
    let g_transform = laplace_inner(simp(g_of_t, pool), t, s, pool, depth + 1, gen)?;
    Ok(Some(pool.mul(vec![exp_neg_as, g_transform])))
}

/// Single-argument primitive functions: sin/cos/sinh/cosh/exp/heaviside/dirac.
#[allow(clippy::too_many_arguments)]
fn laplace_func(
    name: &str,
    arg: ExprId,
    t: ExprId,
    s: ExprId,
    pool: &ExprPool,
    _depth: usize,
    gen: &mut Genericity<'_>,
) -> Result<ExprId, LaplaceError> {
    // exp(a t) → 1/(s − a)
    if name == "exp" {
        let (a, b) = as_affine(arg, t, pool).ok_or_else(|| {
            LaplaceError::NoRule(format!("exp of non-affine argument: {}", pool.display(arg)))
        })?;
        if b != pool.integer(0_i32) {
            return Err(LaplaceError::NoRule(
                "exp(a t + b): nonzero constant offset".into(),
            ));
        }
        let denom = pool.add(vec![s, neg(a, pool)]);
        return Ok(recip(denom, pool));
    }

    // For sin/cos/sinh/cosh the argument must be a pure linear b·t.
    let trig = matches!(name, "sin" | "cos" | "sinh" | "cosh");
    if trig {
        let (b, off) = as_affine(arg, t, pool).ok_or_else(|| {
            LaplaceError::NoRule(format!(
                "{name} of non-affine argument: {}",
                pool.display(arg)
            ))
        })?;
        if off != pool.integer(0_i32) || b == pool.integer(0_i32) {
            return Err(LaplaceError::NoRule(format!(
                "{name}(b t): argument must be a nonzero multiple of t"
            )));
        }
        let b2 = square_of_freq(b, pool);
        let s2 = pool.pow(s, pool.integer(2_i32));
        return Ok(match name {
            // sin(bt) = b/(s²+b²)
            "sin" => {
                let denom = pool.add(vec![s2, b2]);
                pool.mul(vec![b, recip(denom, pool)])
            }
            // cos(bt) = s/(s²+b²)
            "cos" => {
                let denom = pool.add(vec![s2, b2]);
                pool.mul(vec![s, recip(denom, pool)])
            }
            // sinh(bt) = b/(s²−b²)
            "sinh" => {
                let denom = pool.add(vec![s2, neg(b2, pool)]);
                pool.mul(vec![b, recip(denom, pool)])
            }
            // cosh(bt) = s/(s²−b²)
            "cosh" => {
                let denom = pool.add(vec![s2, neg(b2, pool)]);
                pool.mul(vec![s, recip(denom, pool)])
            }
            _ => unreachable!(),
        });
    }

    // θ(t − a) → e^{−a s}/s   (requires a ≥ 0 when `a` is a literal).
    if name == "heaviside" {
        let (coeff, b) = as_affine(arg, t, pool).ok_or_else(|| {
            LaplaceError::NoRule(format!(
                "heaviside of non-affine argument: {}",
                pool.display(arg)
            ))
        })?;
        if coeff != pool.integer(1_i32) {
            return Err(LaplaceError::NoRule(
                "heaviside(c·t − a): coefficient of t must be 1".into(),
            ));
        }
        let a = simp(neg(b, pool), pool); // a = −b
        require_nonneg_shift(a, pool, gen)?;
        let exp_neg_as = pool.func("exp", vec![simp(neg(pool.mul(vec![a, s]), pool), pool)]);
        return Ok(pool.mul(vec![exp_neg_as, recip(s, pool)]));
    }

    // δ(t − a) → e^{−a s}   (δ(t) ↦ 1; requires a ≥ 0 when `a` is a literal).
    if name == "diracdelta" {
        let (coeff, b) = as_affine(arg, t, pool).ok_or_else(|| {
            LaplaceError::NoRule(format!(
                "diracdelta of non-affine argument: {}",
                pool.display(arg)
            ))
        })?;
        if coeff != pool.integer(1_i32) {
            return Err(LaplaceError::NoRule(
                "diracdelta(c·t − a): coefficient of t must be 1".into(),
            ));
        }
        let a = simp(neg(b, pool), pool);
        require_nonneg_shift(a, pool, gen)?;
        return Ok(pool.func("exp", vec![simp(neg(pool.mul(vec![a, s]), pool), pool)]));
    }

    Err(LaplaceError::NoRule(format!("{name}(...)")))
}

/// `b²` for a frequency coefficient, folding `(√c)² → c` so the forward
/// sinh/cosh table emits a ℚ-rational denominator (needed for inverse `apart`).
fn square_of_freq(b: ExprId, pool: &ExprPool) -> ExprId {
    let half = pool.rational(1_i32, 2_i32);
    if let ExprData::Pow { base, exp } = pool.get(b) {
        if exp == half {
            return base;
        }
    }
    simp(pool.pow(b, pool.integer(2_i32)), pool)
}

/// Substitute every occurrence of `from` with `to` in `expr`.
fn subs_one(expr: ExprId, from: ExprId, to: ExprId, pool: &ExprPool) -> ExprId {
    let mut map = std::collections::HashMap::new();
    map.insert(from, to);
    crate::kernel::subs(expr, &map, pool)
}

/// Remove the factor at `idx` from `factors`, returning the product of the rest
/// (or `1` if none remain).
fn remove_index(factors: &[ExprId], idx: usize, pool: &ExprPool) -> ExprId {
    let rest: Vec<ExprId> = factors
        .iter()
        .enumerate()
        .filter(|&(i, _)| i != idx)
        .map(|(_, &f)| f)
        .collect();
    match rest.len() {
        0 => pool.integer(1_i32),
        1 => rest[0],
        _ => pool.mul(rest),
    }
}

// ===========================================================================
// Derivative rule (for the ODE workflow)
// ===========================================================================

/// The Laplace transform of the `order`-th derivative `f^{(order)}(t)` in terms
/// of `F = L{f}(s)` and the initial values `f(0), f'(0), …, f^{(order−1)}(0)`:
///
/// ```text
///   L{f^{(n)}} = sⁿ F − sⁿ⁻¹ f(0) − sⁿ⁻² f'(0) − … − f^{(n−1)}(0).
/// ```
///
/// `initial_values[k]` must be `f^{(k)}(0)`.  Because `f` is an *unknown*
/// function, this operates on the placeholder symbol `f_transform = F(s)` rather
/// than a concrete `f`; it is the building block consumed by an
/// (algebraic-solve) ODE-via-Laplace workflow.  Missing trailing initial values
/// default to `0`.
///
/// Returns the simplified expression for `L{f^{(order)}}(s)`.
pub fn laplace_derivative_rule(
    f_transform: ExprId,
    s: ExprId,
    order: u32,
    initial_values: &[ExprId],
    pool: &ExprPool,
) -> ExprId {
    // sⁿ F
    let s_n = pool.pow(s, pool.integer(order as i32));
    let mut terms = vec![pool.mul(vec![s_n, f_transform])];
    // − Σ_{k=0}^{n−1} s^{n−1−k} f^{(k)}(0)
    for k in 0..order {
        let f_k0 = initial_values
            .get(k as usize)
            .copied()
            .unwrap_or_else(|| pool.integer(0_i32));
        if f_k0 == pool.integer(0_i32) {
            continue;
        }
        let power = (order - 1 - k) as i32;
        let s_pow = pool.pow(s, pool.integer(power));
        terms.push(pool.mul(vec![pool.integer(-1_i32), s_pow, f_k0]));
    }
    simp(pool.add(terms), pool)
}

// ===========================================================================
// Inverse transform
// ===========================================================================

/// Compute the inverse Laplace transform `L⁻¹{F(s)}(t)` for a **proper rational**
/// `F(s)` (optionally times a delay factor `e^{−a s}`).
///
/// Strategy: factor out any `e^{−a s}` delay (→ Heaviside shift), partial-fraction
/// the rational part via [`crate::poly::apart`], then map each term through the
/// inverse table.  See the [module docs](self) for the table and caveats.
///
/// # Errors
///
/// - [`LaplaceError::SameVariable`] if `s == t`.
/// - [`LaplaceError::NotInvertible`] for non-rational `F`, an improper rational
///   (polynomial part ≠ 0 ⇒ derivatives of `δ`), or an irreducible denominator
///   factor of degree > 2.
pub fn inverse_laplace_transform(
    big_f: ExprId,
    s: ExprId,
    t: ExprId,
    pool: &ExprPool,
) -> Result<ExprId, LaplaceError> {
    // Clear first: on the `?` below the call returns without reaching the
    // stash, and a caller that catches the error and then reads the channel
    // must not be handed an *earlier* call's hypotheses.
    stash_transform_side_conditions(Vec::new());
    let (out, conds) = inverse_laplace_transform_with_conditions(big_f, s, t, pool)?;
    stash_transform_side_conditions(conds);
    Ok(out)
}

/// [`inverse_laplace_transform`] with its genericity hypotheses returned in
/// band rather than through
/// [`super::take_transform_side_conditions`].
///
/// The list is empty for every input whose coefficients are rational numbers.
/// It is non-empty exactly when the answer depends on a fact about a symbolic
/// parameter that the input does not settle — `ka ≠ ke` for the two-compartment
/// denominator `(s+ka)(s+ke)`, `ζ < 1` for the under-damped second-order form.
pub fn inverse_laplace_transform_with_conditions(
    big_f: ExprId,
    s: ExprId,
    t: ExprId,
    pool: &ExprPool,
) -> Result<(ExprId, Vec<SideCondition>), LaplaceError> {
    let mut g = Genericity::new(None);
    let out = inverse_laplace_inner(big_f, s, t, pool, &mut g)?;
    Ok((out, g.into_conditions()))
}

/// [`inverse_laplace_transform_with_conditions`] under an explicit
/// [`AssumptionContext`].
///
/// A hypothesis the context **proves** is discharged rather than reported;
/// anything else is still reported. An assumption context makes conditions
/// disappear only by settling them, never by being present.
///
/// # How much it can settle
///
/// Discharge goes through, in order: a literal value; a matching
/// [`crate::deriv::SideCondition`] fact from the context or from a symbol's
/// static [`crate::kernel::Domain`]; the structural rules "a product is
/// non-zero when every factor is" and "a positive quantity is non-zero"; and
/// finally [`crate::logic::satisfiable`] on `assumptions ∧ ¬goal`.
///
/// That last step is an **interval** decision procedure, not real quantifier
/// elimination, so it answers `Unknown` for anything nonlinear — and `Unknown`
/// is not a proof. Concretely, declaring `0 < ζ < 1 ∧ ω > 0` discharges
/// `ω ≠ 0` but leaves `ζ ≠ ±1` and `ω²(1 − ζ²) > 0` on the wire, even though
/// both follow. Over-reporting is the safe direction and it is the one taken;
/// closing the gap is a job for [`crate::real::cad`], whose cost is not
/// something a transform call should pay unasked.
pub fn inverse_laplace_transform_with_assumptions(
    big_f: ExprId,
    s: ExprId,
    t: ExprId,
    pool: &ExprPool,
    assumptions: &AssumptionContext,
) -> Result<(ExprId, Vec<SideCondition>), LaplaceError> {
    let mut g = Genericity::new(Some(assumptions));
    let out = inverse_laplace_inner(big_f, s, t, pool, &mut g)?;
    Ok((out, g.into_conditions()))
}

fn inverse_laplace_inner(
    big_f: ExprId,
    s: ExprId,
    t: ExprId,
    pool: &ExprPool,
    gen: &mut Genericity<'_>,
) -> Result<ExprId, LaplaceError> {
    if s == t {
        return Err(LaplaceError::SameVariable);
    }

    // Peel a delay factor e^{−a s}: L⁻¹{e^{−a s} G(s)} = θ(t−a)·g(t−a).  The
    // unilateral rule needs `a ≥ 0`; a literal negative `a` is refused by the
    // forward table, and a *symbolic* `a` is a hypothesis, not a fact.
    if let Some((a, g)) = split_delay(big_f, s, pool) {
        require_nonneg_shift(a, pool, gen).map_err(|_| {
            LaplaceError::inverse_not_causal(format!(
                "advance factor e^{{{}·s}} — the unilateral inverse needs a delay, not an \
                 advance, and no causal f has this transform",
                pool.display(simp(neg(a, pool), pool))
            ))
        })?;
        let g_inv = inverse_laplace_inner(g, s, t, pool, gen)?;
        let t_minus_a = simp(pool.add(vec![t, neg(a, pool)]), pool);
        let shifted = subs_one(g_inv, t, t_minus_a, pool);
        let heaviside = pool.func("heaviside", vec![t_minus_a]);
        return Ok(simp(pool.mul(vec![heaviside, shifted]), pool));
    }

    // Peel s-free scalar factors (e.g. √2 in √2/(s²−2) from L{sinh(√2 t)})
    // so `apart` sees a rational function of `s`.  Over ℚ(params) most such
    // factors are ordinary coefficients, but a genuine algebraic number (√2)
    // or a transcendental (π) still has to come out first.
    if let Some((scalar, rest)) = split_s_free_scalar(big_f, s, pool) {
        let rest_inv = inverse_laplace_inner(rest, s, t, pool, gen)?;
        return Ok(simp(pool.mul(vec![scalar, rest_inv]), pool));
    }

    // Rational route: partial fractions over ℚ — or, when `F` mentions symbolic
    // parameters, over ℚ(params) — then the table per term.  `apart` reports
    // any genericity its decomposition rests on; those are the caller's
    // hypotheses too.
    let (pf, apart_conds) = crate::poly::apart_with_conditions(big_f, s, pool)
        .map_err(|e| LaplaceError::NotInvertible(format!("apart failed: {e}")))?;
    for cond in apart_conds {
        gen.adopt(cond, pool);
    }

    let terms: Vec<ExprId> = match pool.get(pf) {
        ExprData::Add(args) => args,
        _ => vec![pf],
    };

    let mut out = Vec::with_capacity(terms.len());
    for term in terms {
        out.push(invert_term(term, s, t, pool, gen)?);
    }
    Ok(simp(pool.add(out), pool))
}

/// If `F` is a product with at least one factor free of `s`, return
/// `(product of s-free factors, product of the rest)`.  Used so algebraic
/// amplitudes from the forward sinh/cosh table do not block `apart`.
fn split_s_free_scalar(big_f: ExprId, s: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId)> {
    let factors: Vec<ExprId> = match pool.get(big_f) {
        ExprData::Mul(a) => a,
        _ => return None,
    };
    let mut free = Vec::new();
    let mut rest = Vec::new();
    for &fac in &factors {
        if is_free_of(fac, s, pool) {
            free.push(fac);
        } else {
            rest.push(fac);
        }
    }
    if free.is_empty() || rest.is_empty() {
        return None;
    }
    let scalar = match free.len() {
        1 => free[0],
        _ => pool.mul(free),
    };
    let body = match rest.len() {
        1 => rest[0],
        _ => pool.mul(rest),
    };
    Some((scalar, body))
}

/// Split `F = e^{−a s} · G(s)` returning `(a, G)`, or `None` if there is no
/// exponential delay factor.
fn split_delay(big_f: ExprId, s: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId)> {
    let factors: Vec<ExprId> = match pool.get(big_f) {
        ExprData::Mul(a) => a,
        _ => vec![big_f],
    };
    for (i, &fac) in factors.iter().enumerate() {
        if let ExprData::Func { name, args } = pool.get(fac) {
            if name == "exp" && args.len() == 1 {
                // arg should be −a·s (linear in s, no constant)
                if let Some((coeff, b)) = as_affine(args[0], s, pool) {
                    if b == pool.integer(0_i32) && coeff != pool.integer(0_i32) {
                        let a = simp(neg(coeff, pool), pool); // arg = −a s ⇒ a = −coeff
                        let g = remove_index(&factors, i, pool);
                        return Some((a, g));
                    }
                }
            }
        }
    }
    None
}

/// Invert a single partial-fraction term `A·(s−a)^{−n}` or a quadratic-denominator
/// term `(B s + C)/((s−p)² + ω²)`.
fn invert_term(
    term: ExprId,
    s: ExprId,
    t: ExprId,
    pool: &ExprPool,
    gen: &mut Genericity<'_>,
) -> Result<ExprId, LaplaceError> {
    // Split term = numer · denom_pow, where denom_pow = D(s)^{−n}.
    let (numer, base, n) = split_rational_term(term, s, pool)
        .ok_or_else(|| LaplaceError::NotInvertible(pool.display(term).to_string()))?;

    // A polynomial part (n == 0) is the transform of `δ(t)` (constant term) or
    // its derivatives (higher terms) — outside the rational table.  We decline
    // it rather than emit distributional `δ` output, keeping the inverse on
    // ordinary functions (matching the "proper rational" scope).
    if n == 0 {
        return Err(LaplaceError::NotInvertible(format!(
            "polynomial part {} (δ / derivatives of δ — improper rational)",
            pool.display(term)
        )));
    }

    match poly_degree(base, s, pool) {
        Some(1) => invert_linear_pole(numer, base, n, s, t, pool),
        Some(2) => invert_quadratic(numer, base, n, s, t, pool, gen),
        _ => Err(LaplaceError::NotInvertible(format!(
            "denominator factor of degree > 2: {}",
            pool.display(base)
        ))),
    }
}

/// Decompose a term into `(numerator, denom_base, n)` with `term = numerator ·
/// denom_base^{−n}` and `n ≥ 0`, `numerator` free of any negative power of `s`.
///
/// A negative power whose base is free of `s` — `(ka − ke)^{−1}` in the
/// ℚ(params) decomposition of the two-compartment model — is a *coefficient*,
/// not a pole, and belongs in the numerator. Reading it as a second pole is
/// what used to make every parametric term look like "two distinct denominator
/// factors" and fail.
fn split_rational_term(term: ExprId, s: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId, u64)> {
    let factors: Vec<ExprId> = match pool.get(term) {
        ExprData::Mul(a) => a,
        _ => vec![term],
    };
    let mut numer_parts: Vec<ExprId> = Vec::new();
    let mut base: Option<ExprId> = None;
    let mut n: u64 = 0;

    for &fac in &factors {
        if let ExprData::Pow { base: b, exp } = pool.get(fac) {
            if let ExprData::Integer(e) = pool.get(exp) {
                let ev = e.0;
                if ev < 0 && !is_free_of(b, s, pool) {
                    // negative power in `s` → part of denominator
                    if base.is_some() && base != Some(b) {
                        // Two distinct denominator factors — not a single PF term.
                        return None;
                    }
                    base = Some(b);
                    n = (-ev).to_u64()?;
                    continue;
                }
            }
        }
        numer_parts.push(fac);
    }

    let numer = match numer_parts.len() {
        0 => pool.integer(1_i32),
        1 => numer_parts[0],
        _ => pool.mul(numer_parts),
    };
    match base {
        Some(b) => Some((numer, b, n)),
        None => Some((numer, pool.integer(1_i32), 0)),
    }
}

/// Degree of `base` as a polynomial in `s` (only handles `s`, `s ± c`, and
/// `s² + …` forms via structural inspection).  Returns `None` if not obviously
/// polynomial of degree 1 or 2.
fn poly_degree(base: ExprId, s: ExprId, pool: &ExprPool) -> Option<u64> {
    if base == s {
        return Some(1);
    }
    match pool.get(base) {
        ExprData::Add(args) => {
            let mut deg = 0u64;
            for a in args {
                deg = deg.max(monomial_degree(a, s, pool)?);
            }
            Some(deg)
        }
        ExprData::Pow { .. } | ExprData::Mul(_) => monomial_degree(base, s, pool),
        _ if is_free_of(base, s, pool) => Some(0),
        _ => None,
    }
}

fn monomial_degree(term: ExprId, s: ExprId, pool: &ExprPool) -> Option<u64> {
    if term == s {
        return Some(1);
    }
    if is_free_of(term, s, pool) {
        return Some(0);
    }
    match pool.get(term) {
        ExprData::Pow { base, exp } if base == s => nonneg_int_exp(exp, pool),
        ExprData::Mul(args) => {
            let mut deg = 0u64;
            for a in args {
                deg += monomial_degree(a, s, pool)?;
            }
            Some(deg)
        }
        _ => None,
    }
}

/// `L⁻¹{A/(s−a)^n} = A·t^{n−1} e^{a t}/(n−1)!`.
fn invert_linear_pole(
    numer: ExprId,
    base: ExprId,
    n: u64,
    s: ExprId,
    t: ExprId,
    pool: &ExprPool,
) -> Result<ExprId, LaplaceError> {
    // numer must be free of s (deg < deg(base) = 1).
    if !is_free_of(numer, s, pool) {
        return Err(LaplaceError::NotInvertible(format!(
            "linear-pole numerator depends on s: {}",
            pool.display(numer)
        )));
    }
    // base = s − a (monic).  Extract a from (coeff·s + b): a = −b/coeff, require coeff = 1.
    let (coeff, b) = as_affine(base, s, pool)
        .ok_or_else(|| LaplaceError::NotInvertible(pool.display(base).to_string()))?;
    if coeff != pool.integer(1_i32) {
        return Err(LaplaceError::NotInvertible(
            "non-monic linear denominator".into(),
        ));
    }
    let a = simp(neg(b, pool), pool); // a = −b

    let exp_at = pool.func("exp", vec![pool.mul(vec![a, t])]);
    let mut parts = vec![numer, exp_at];
    if n >= 2 {
        let t_pow = pool.pow(t, pool.integer(Integer::from(n - 1)));
        parts.push(t_pow);
        let fact = factorial(n - 1, pool);
        parts.push(recip(fact, pool));
    }
    Ok(pool.mul(parts))
}

/// Invert `(B s + C)/((s−p)² ± λ²)^n` for `n ∈ {1, 2}`.
///
/// Completing the square yields `ω² = γ − β²/4`.  The sign selects the table:
///
/// ```text
///   ω² > 0:  oscillatory — sin/cos with ω = √(ω²)
///   ω² < 0:  hyperbolic  — sinh/cosh with κ = √(−ω²)
///   ω² = 0:  declined (degenerate; should have been a linear factor)
///
///   n = 1, oscillatory:
///     e^{p t} ( B cos(ω t) + ((C + B p)/ω) sin(ω t) )
///   n = 1, hyperbolic:
///     e^{p t} ( B cosh(κ t) + ((C + B p)/κ) sinh(κ t) )
///
///   n = 2: analogous t-weighted forms (sin−ωt cos ↔ κt cosh−sinh).
/// ```
///
/// Higher powers (`n ≥ 3`) are declined.  The `n = 2` oscillatory case is
/// required for round-trips of `t·sin` / `t·cos`.
fn invert_quadratic(
    numer: ExprId,
    base: ExprId,
    n: u64,
    s: ExprId,
    t: ExprId,
    pool: &ExprPool,
    gen: &mut Genericity<'_>,
) -> Result<ExprId, LaplaceError> {
    if n != 1 && n != 2 {
        return Err(LaplaceError::NotInvertible(
            "repeated irreducible quadratic pole (n ≥ 3) not in table".into(),
        ));
    }
    // Write base = s² + β s + γ. Complete the square: (s + β/2)² + (γ − β²/4).
    let (alpha, beta, gamma) = quadratic_coeffs(base, s, pool)
        .ok_or_else(|| LaplaceError::NotInvertible(pool.display(base).to_string()))?;
    if alpha != pool.integer(1_i32) {
        return Err(LaplaceError::NotInvertible(
            "non-monic quadratic denominator".into(),
        ));
    }
    // p = −β/2 ; ω² = γ − β²/4.
    let half = pool.rational(1_i32, 2_i32);
    let p = simp(pool.mul(vec![neg(beta, pool), half]), pool);
    let beta2 = pool.pow(beta, pool.integer(2_i32));
    let quarter = pool.rational(1_i32, 4_i32);
    let omega_sq = simp(
        pool.add(vec![gamma, neg(pool.mul(vec![beta2, quarter]), pool)]),
        pool,
    );

    // The sign of ω² selects sin/cos vs sinh/cosh, and ω² = 0 is a third,
    // genuinely different closed form (`t·e^{pt}`, a double real pole).  A
    // literal decides itself.  A *parametric* ω² decides nothing on its own:
    //
    //   • if the caller's assumptions settle the sign, use it and say nothing;
    //   • otherwise take the oscillatory branch — which is the analytic
    //     continuation of both, since `cos ωt` and `sin(ωt)/ω` are entire in
    //     ω² — and record what makes it the *real-valued* answer, namely
    //     ω² > 0, plus ω² ≠ 0 which rules out the double-pole form.
    //
    // Silently taking the branch, as this did before, returns an expression
    // that is imaginary at ζ > 1 with nothing on the wire to say so.
    let hyperbolic = match literal_rational(omega_sq, pool) {
        Some(r) if r < 0 => true,
        Some(r) if r == 0 => {
            return Err(LaplaceError::NotInvertible(
                "degenerate quadratic pole (ω² = 0)".into(),
            ));
        }
        Some(_) => false,
        None => {
            if gen.known_negative(omega_sq, pool) {
                true
            } else {
                if !gen.known_positive(omega_sq, pool) {
                    gen.need_nonzero(omega_sq, pool);
                    gen.need_positive(omega_sq, pool);
                }
                false
            }
        }
    };

    let freq_sq = if hyperbolic {
        simp(neg(omega_sq, pool), pool) // κ² = −ω²
    } else {
        omega_sq
    };
    // ω or κ.  When ω² is a perfect square — `w²` from `L{sin(w t)}` — take the
    // root symbolically rather than leaving `(w²)^{1/2}`.  This is sound for
    // *either* branch of the root because every form emitted below is an even
    // function of the frequency: `cos ωt`, `sin(ωt)/ω`, `t·sin(ωt)/(2ω)` and
    // `(sin ωt − ωt cos ωt)/(2ω³)` are all unchanged by `ω ↦ −ω`.  So
    // `sin(√(w²)·t)/√(w²)` and `sin(w t)/w` are the same function of `w`.
    let freq = match perfect_square_root(freq_sq, pool) {
        Some(r) => r,
        None => simp(pool.pow(freq_sq, half), pool),
    };

    // numerator B s + C.
    let (bb, cc) = as_affine(numer, s, pool)
        .ok_or_else(|| LaplaceError::NotInvertible(pool.display(numer).to_string()))?;

    let exp_pt = pool.func("exp", vec![pool.mul(vec![p, t])]);
    let freq_t = pool.mul(vec![freq, t]);
    let (odd_fn, even_fn) = if hyperbolic {
        (
            pool.func("sinh", vec![freq_t]),
            pool.func("cosh", vec![freq_t]),
        )
    } else {
        (
            pool.func("sin", vec![freq_t]),
            pool.func("cos", vec![freq_t]),
        )
    };

    if n == 1 {
        // e^{p t} [ B · even(freq·t) + ((C + B p)/freq) · odd(freq·t) ]
        let even_term = pool.mul(vec![bb, even_fn]);
        let bp = pool.mul(vec![bb, p]);
        let odd_coeff = pool.mul(vec![pool.add(vec![cc, bp]), recip(freq, pool)]);
        let odd_term = pool.mul(vec![odd_coeff, odd_fn]);
        return Ok(pool.mul(vec![exp_pt, pool.add(vec![even_term, odd_term])]));
    }

    // n = 2.
    // Oscillatory: B·(t/(2ω)) sin + (Bp+C)/(2ω³)·(sin − ωt cos)
    // Hyperbolic:  B·(t/(2κ)) sinh + (Bp+C)/(2κ³)·(κt cosh − sinh)
    let two = pool.integer(2_i32);
    let two_freq = pool.mul(vec![two, freq]);
    let bp_plus_c = pool.add(vec![pool.mul(vec![bb, p]), cc]);

    let t_odd = pool.mul(vec![bb, t, recip(two_freq, pool), odd_fn]);
    let combo = if hyperbolic {
        pool.add(vec![pool.mul(vec![freq, t, even_fn]), neg(odd_fn, pool)])
    } else {
        pool.add(vec![odd_fn, neg(pool.mul(vec![freq, t, even_fn]), pool)])
    };
    let freq3 = pool.mul(vec![freq, freq, freq]);
    let two_freq3 = pool.mul(vec![two, freq3]);
    let second = pool.mul(vec![bp_plus_c, recip(two_freq3, pool), combo]);
    Ok(pool.mul(vec![exp_pt, pool.add(vec![t_odd, second])]))
}

/// A symbolic square root of `e` when one is evident: `w²` → `w`, `9·a⁴` →
/// `3·a²`, a perfect-square literal → its root.  `None` when no exact root is
/// apparent, in which case the caller keeps the formal `e^{1/2}`.
///
/// The sign of the root is not determined and does not need to be — see the
/// evenness argument at the call site.
fn perfect_square_root(e: ExprId, pool: &ExprPool) -> Option<ExprId> {
    if let Some(r) = literal_rational(e, pool) {
        let (n, d) = r.into_numer_denom();
        if n < 0 {
            return None;
        }
        let (ns, nrem) = n.clone().sqrt_rem(Integer::new());
        let (ds, drem) = d.clone().sqrt_rem(Integer::new());
        if nrem != 0 || drem != 0 {
            return None;
        }
        return Some(if ds == 1 {
            pool.integer(ns)
        } else {
            pool.rational(ns, ds)
        });
    }
    match pool.get(e) {
        ExprData::Pow { base, exp } => {
            let ExprData::Integer(k) = pool.get(exp) else {
                return None;
            };
            let k = k.0.to_i64()?;
            if k <= 0 || k % 2 != 0 {
                return None;
            }
            let half = k / 2;
            Some(if half == 1 {
                base
            } else {
                pool.pow(base, pool.integer(half as i32))
            })
        }
        ExprData::Mul(args) => {
            let mut parts = Vec::with_capacity(args.len());
            for a in args {
                parts.push(perfect_square_root(a, pool)?);
            }
            Some(pool.mul(parts))
        }
        _ => None,
    }
}

/// The unilateral second-shift rule `θ(t−a)·g(t−a) ↦ e^{−a s}·G(s)` needs
/// `a ≥ 0`, and the two ways that can fail are not the same failure.
///
/// A *literal* negative `a` refutes the hypothesis outright: `θ(t+1) ≡ 1` on
/// `t ≥ 0`, so its transform is `1/s`, and `e^{s}/s` is not a conditional
/// answer but a wrong one.  That is refused.
///
/// A *symbolic* `a` settles nothing either way, so it is recorded as a
/// hypothesis — unless the symbol's static [`crate::kernel::Domain`] or the
/// ambient facts already prove it.  Before this was recorded,
/// `L{θ(t+a)} = e^{a s}/s` came back with nothing on the wire to say it is the
/// transform only for `a ≤ 0`: at `a = 1.5, s = 2.5` it reads `17.008` against
/// a true `1/s = 0.4`.
fn require_nonneg_shift(
    a: ExprId,
    pool: &ExprPool,
    gen: &mut Genericity<'_>,
) -> Result<(), LaplaceError> {
    if let Some(r) = literal_rational(a, pool) {
        if r < 0 {
            return Err(LaplaceError::forward_not_causal(format!(
                "shift a = {} must be ≥ 0 for unilateral Heaviside/Dirac",
                pool.display(a)
            )));
        }
        return Ok(());
    }
    let cond = SideCondition::InDomain(a, crate::kernel::Domain::NonNegative);
    if gen.refuted(&cond, pool) {
        return Err(LaplaceError::forward_not_causal(format!(
            "shift a = {} is known to be negative; the unilateral Heaviside/Dirac \
             rule needs a ≥ 0",
            pool.display(a)
        )));
    }
    gen.adopt(cond, pool);
    Ok(())
}

/// If `expr` is a literal rational (integer or ratio), return it.
fn literal_rational(expr: ExprId, pool: &ExprPool) -> Option<rug::Rational> {
    match pool.get(expr) {
        ExprData::Integer(n) => Some(rug::Rational::from(n.0.clone())),
        ExprData::Rational(r) => Some(r.0.clone()),
        _ => None,
    }
}

/// Coefficients `(α, β, γ)` of a monic-or-scaled quadratic `α s² + β s + γ`.
fn quadratic_coeffs(base: ExprId, s: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId, ExprId)> {
    let args: Vec<ExprId> = match pool.get(base) {
        ExprData::Add(a) => a,
        _ => vec![base],
    };
    let mut alpha = pool.integer(0_i32);
    let mut beta = pool.integer(0_i32);
    let mut gamma_parts: Vec<ExprId> = Vec::new();
    for term in args {
        match monomial_degree(term, s, pool)? {
            2 => alpha = monomial_coeff(term, s, 2, pool)?,
            1 => beta = monomial_coeff(term, s, 1, pool)?,
            0 => gamma_parts.push(term),
            _ => return None,
        }
    }
    let gamma = match gamma_parts.len() {
        0 => pool.integer(0_i32),
        1 => gamma_parts[0],
        _ => pool.add(gamma_parts),
    };
    Some((alpha, beta, gamma))
}

/// Coefficient of `s^deg` in a single monomial `c·s^deg`.
fn monomial_coeff(term: ExprId, s: ExprId, deg: u64, pool: &ExprPool) -> Option<ExprId> {
    if deg == 0 {
        return Some(term);
    }
    if deg == 1 && term == s {
        return Some(pool.integer(1_i32));
    }
    if let ExprData::Pow { base, exp } = pool.get(term) {
        if base == s && nonneg_int_exp(exp, pool) == Some(deg) {
            return Some(pool.integer(1_i32));
        }
    }
    if let ExprData::Mul(args) = pool.get(term) {
        let mut coeff_parts: Vec<ExprId> = Vec::new();
        let mut found = false;
        for a in args {
            if a == s && deg == 1 {
                found = true;
                continue;
            }
            if let ExprData::Pow { base, exp } = pool.get(a) {
                if base == s && nonneg_int_exp(exp, pool) == Some(deg) {
                    found = true;
                    continue;
                }
            }
            coeff_parts.push(a);
        }
        if found {
            return Some(match coeff_parts.len() {
                0 => pool.integer(1_i32),
                1 => coeff_parts[0],
                _ => pool.mul(coeff_parts),
            });
        }
    }
    None
}

#[cfg(test)]
mod tests;