symplex 0.39.0

//! Compute residues of functions at poles of any order.
//!
//! The residue of `f(x)` at `x = a` is the coefficient of `1/(x−a)` in the
//! Laurent expansion of `f` around `a`.
//!
//! # Algorithm
//!
//! 1. Shift: `g(t) = f(a + t)` and split `g = N(t)/D(t)`.
//! 2. If `D` is a polynomial in `t`, its multiplicity `m` at `t = 0` is the
//!    index of the first non-zero coefficient.  With `D̃(t) = D(t)/tᵐ`
//!    (exact, non-vanishing at 0), the residue is the Taylor coefficient
//!    `Res = 1/(m−1)! · dᵐ⁻¹/dtᵐ⁻¹ [N(t)/D̃(t)]` at `t = 0`, which is a plain
//!    substitution — no limits needed.  This handles algebraic points such
//!    as `a = i` where polynomial cancellation over ℚ is unavailable.
//! 3. Otherwise (transcendental denominators such as `1/sin t`), the order
//!    is probed with `lim tᵏ g(t)` for `k = 1, …, 6` and the same derivative
//!    formula is evaluated as a limit.
//! 4. Both formulas presuppose a pole of finite order: `N` (and, in step 3,
//!    `D`) must be *analytic* at `t = 0`.  That is certified structurally
//!    ([`analytic_at_zero`]); a finite value `N(0)` is not enough —
//!    `t·e^{1/t}` is `0` at `t = 0` after `0·x → 0`, and `t²·sin(1/t)` has
//!    real limits although `t = 0` is an essential singularity.
//! 5. When the certificate fails, `t = 0` may be an essential singularity.
//!    If `g` is single-valued and analytic on the Riemann sphere except at
//!    `0` and finitely many poles — rational functions of `t` times entire
//!    functions (`exp`, `sin`, `cos`, `sinh`, `cosh`, `erf`, `erfc`, `c^w`)
//!    of Laurent polynomials in `1/t` ([`singular_polys_off_origin`]) —
//!    the residue theorem on the sphere gives
//!    `Res₀ g = −Res_∞ g − Σ_{p ≠ 0} Res_p g`, where the other poles `p` are
//!    the (explicitly solved) roots of the polynomial denominators and
//!    `Res_∞`, `Res_p` are residues at poles of finite order (steps 2–3).
//!    This is the coefficient of `t⁻¹` of the product of the Laurent
//!    series: for `g = P(t)·h(1/t)` with `P` a polynomial it is the finite
//!    sum `Σₙ pₙ₋₁ hₙ` (`Res(z²e^{1/z}, 0) = 1/6`, `Res(z²sin(1/z), 0) =
//!    −1/6`); for a rational `P` with other poles the infinite series sums
//!    to the closed form (`Res(e^{1/z}/(z − 1), 0) = 1 − e`).  When the
//!    series does not close this way (`e^z·e^{1/z}`: `∞` is essential too,
//!    the value is `I₁(2)`), the residue is refused.
//!
//! If nothing above succeeds, an error is returned and the caller keeps a
//! formal `Residue` node.

use num_bigint::BigInt;
use num_traits::One;

use crate::base::arena::Arena;
use crate::base::errors::SymplexError;
use crate::base::node::{ExprId, ExprNode};
use crate::base::walk;
use crate::calculus::calculus_util;
use crate::transforms::{diff, eval, subs};

/// Maximum pole order probed by the limit-based fallback.
const MAX_PROBE_ORDER: usize = 6;

fn failed(reason: impl Into<String>) -> SymplexError {
    SymplexError::computation_failed("residue", reason)
}

/// Compute the residue of `expr` at `var = point`.
///
/// See the module documentation for the algorithm.  Returns `Err` when the
/// residue cannot be determined (essential singularity, unknown pole
/// order, or a non-symbol `var`).
pub(crate) fn residue(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
) -> Result<ExprId, SymplexError> {
    residue_with(arena, expr, var, point, true)
}

/// [`residue`], with the essential-singularity route (module step 5)
/// enabled or not.  The route itself computes residues at poles only
/// (`essential = false`), so the mutual calls stop after one level.
fn residue_with(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
    essential: bool,
) -> Result<ExprId, SymplexError> {
    if !matches!(arena.node(var), ExprNode::Symbol(_)) {
        return Err(failed("residue variable must be a symbol"));
    }
    if point == arena.infinity() {
        return residue_at_infinity_with(arena, expr, var, essential);
    }
    let result = residue_inner(arena, expr, var, point, essential)?;
    // The limit engine can leak internal dummy symbols into a half-finished
    // result; a residue may only mention symbols of the input.
    let mut allowed: Vec<ExprId> = walk::free_symbols(arena, expr);
    allowed.extend(walk::free_symbols(arena, point));
    if walk::free_symbols(arena, result)
        .iter()
        .any(|s| !allowed.contains(s))
    {
        return Err(failed("limit engine returned an incomplete result"));
    }
    Ok(result)
}

/// [`residue`] without the final symbol-leak check.
fn residue_inner(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
    essential: bool,
) -> Result<ExprId, SymplexError> {
    // Shift the pole to the origin: g(t) = f(point + t).
    let t = arena.symbol("__res_t");
    let shifted_var = arena.add(&[point, t]);
    let g = subs::subs(arena, expr, var, shifted_var);
    let g = eval::eval(arena, g);

    if !walk::contains(arena, g, t) {
        // Constant in t: analytic, residue 0 (unless it is itself non-finite).
        return Ok(arena.zero());
    }

    match residue_at_pole(arena, g, t) {
        Err(pole_err) if essential => essential_residue(arena, g, t).map_err(|e| {
            failed(format!(
                "{}; as an essential singularity: {}",
                pole_reason(&pole_err),
                pole_reason(&e)
            ))
        }),
        other => other,
    }
}

/// The reason text of a [`failed`] error (the whole message otherwise).
fn pole_reason(e: &SymplexError) -> String {
    match e {
        SymplexError::ComputationFailed { reason, .. } => reason.clone(),
        other => other.to_string(),
    }
}

/// Residue of `g` at `t = 0` when it is a pole of finite order (module
/// steps 2–4); an error for anything else.
fn residue_at_pole(arena: &mut Arena, g: ExprId, t: ExprId) -> Result<ExprId, SymplexError> {
    // Clear nested fractions (e.g. from f(1/t)) so that the denominator is a
    // genuine polynomial whenever f is rational.
    let g = crate::transforms::integrate::clear_nested_fractions(arena, g, t);
    let (numer, denom) = crate::poly::polybridge::as_numer_denom(arena, g);

    if let Some(coeffs) = calculus_util::poly_coeffs_symbolic(arena, denom, t) {
        return residue_polynomial_denominator(arena, numer, &coeffs, t);
    }

    // The real-axis limits below determine the order of a *pole*; they say
    // nothing about an essential singularity (`sin(1/t)` is bounded on the
    // real axis).  Require a quotient of functions analytic at 0.
    if !analytic_at_zero(arena, numer, t) || !analytic_at_zero(arena, denom, t) {
        return Err(failed(
            "not a quotient of functions analytic at the point (essential singularity or branch point?)",
        ));
    }
    residue_by_limits(arena, g, t)
}

/// Residue when the shifted denominator is a polynomial in `t` with the
/// given coefficients.
fn residue_polynomial_denominator(
    arena: &mut Arena,
    numer: ExprId,
    coeffs: &[ExprId],
    t: ExprId,
) -> Result<ExprId, SymplexError> {
    // Multiplicity of t = 0: index of the first coefficient that is not
    // provably zero.
    let mut m = coeffs.len();
    for (k, &c) in coeffs.iter().enumerate() {
        if !is_provably_zero(arena, c) {
            m = k;
            break;
        }
    }
    if m == coeffs.len() {
        return Err(failed("denominator is identically zero"));
    }

    // Numerator must be analytic at t = 0 (a finite value N(0) is not
    // enough: t·e^{1/t} evaluates to 0 there).
    let n0 = subs::subs(arena, numer, t, arena.zero());
    let n0 = eval::eval(arena, n0);
    if !is_finite_constant(arena, n0, t) || !analytic_at_zero(arena, numer, t) {
        return Err(failed(
            "numerator is not analytic at the point (essential singularity or branch point?)",
        ));
    }

    if m == 0 {
        // No pole: the function is analytic here.
        return Ok(arena.zero());
    }

    // D̃(t) = Σ_{k ≥ m} c_k t^{k−m}
    let mut terms: Vec<ExprId> = Vec::with_capacity(coeffs.len() - m);
    for (k, &c) in coeffs.iter().enumerate().skip(m) {
        let e = arena.int((k - m) as i64);
        let tp = arena.pow(t, e);
        terms.push(arena.mul(&[c, tp]));
    }
    let d_tilde = arena.add(&terms);

    // h(t) = N(t)/D̃(t); Res = h^{(m−1)}(0) / (m−1)!
    let mut h = arena.div(numer, d_tilde);
    for _ in 1..m {
        h = diff::diff(arena, h, t);
    }
    let h0 = subs::subs(arena, h, t, arena.zero());
    let h0 = eval::eval(arena, h0);
    if !is_finite_constant(arena, h0, t) {
        return Err(failed("could not evaluate the Laurent coefficient"));
    }
    let res = divide_by_factorial(arena, h0, m - 1);
    Ok(eval::eval(arena, res))
}

/// Limit-based fallback for non-polynomial denominators.
fn residue_by_limits(arena: &mut Arena, g: ExprId, t: ExprId) -> Result<ExprId, SymplexError> {
    let zero = arena.zero();
    for m in 1..=MAX_PROBE_ORDER {
        let e = arena.int(m as i64);
        let tm = arena.pow(t, e);
        let scaled = arena.mul(&[tm, g]);
        let scaled = eval::eval(arena, scaled);
        let Ok(l) = crate::calculus::limit::limit(arena, scaled, t, zero) else {
            continue;
        };
        if !is_finite_constant(arena, l, t) {
            continue;
        }
        if m == 1 {
            return Ok(eval::eval(arena, l));
        }
        // Order ≤ m established: Res = lim d^{m−1}/dt^{m−1}[tᵐ g] / (m−1)!
        let mut h = scaled;
        for _ in 1..m {
            h = diff::diff(arena, h, t);
        }
        let h = eval::eval(arena, h);
        let l2 = crate::calculus::limit::limit(arena, h, t, zero)?;
        if !is_finite_constant(arena, l2, t) {
            return Err(failed("Laurent coefficient limit is not finite"));
        }
        let res = divide_by_factorial(arena, l2, m - 1);
        return Ok(eval::eval(arena, res));
    }
    Err(failed(
        "could not determine the pole order (essential singularity or order > 6)",
    ))
}

/// Residue at an essential singularity at `t = 0` by the residue theorem
/// on the Riemann sphere (module step 5):
/// `Res₀ g = −Res_∞ g − Σ_{p ≠ 0} Res_p g`.
///
/// Valid because [`singular_polys_off_origin`] certifies that `g` is
/// analytic on `ℂ ∖ {0}` except at the roots of the returned polynomials,
/// all of which are found explicitly ([`nonzero_roots`]); the residues at
/// those roots and at `∞` are computed with the pole formulas only, which
/// refuse anything but a pole of finite order (so a second essential
/// singularity, at `∞` or elsewhere, makes the whole computation refuse).
fn essential_residue(arena: &mut Arena, g: ExprId, t: ExprId) -> Result<ExprId, SymplexError> {
    let polys = singular_polys_off_origin(arena, g, t).ok_or_else(|| {
        failed("not a rational function times entire functions of polynomials in 1/(z − z0)")
    })?;
    let roots = nonzero_roots(arena, &polys, t)
        .ok_or_else(|| failed("the other poles could not all be found exactly"))?;
    let mut parts: Vec<ExprId> = Vec::with_capacity(roots.len() + 1);
    parts.push(residue_at_infinity_with(arena, g, t, false)?);
    for p in roots {
        parts.push(residue_inner(arena, g, t, p, false)?);
    }
    let sum = arena.add(&parts);
    let neg = arena.neg(sum);
    let res = eval::eval(arena, neg);
    if !is_finite_constant(arena, res, t) {
        return Err(failed(
            "the residues of the other singularities are not finite",
        ));
    }
    Ok(res)
}

/// Certify that `g` is single-valued and analytic on `ℂ ∖ {0}` except at
/// the roots of finitely many polynomials in `t`, and return those
/// polynomials (the bases of negative integer powers).
///
/// Every subexpression depending on `t` is checked with an explicit stack:
/// sums, products and non-negative integer powers; the entire functions
/// `exp`, `sin`, `cos`, `sinh`, `cosh`, `erf`, `erfc` and `c^w` (`c ≠ 0`
/// free of `t`); negative integer powers of products of polynomials in `t`
/// (poles at their roots) and of exponentials (which never vanish).
/// Anything else — fractional powers of `t`, `ln`, `tan`, `1/(e^{1/t} − 1)`
/// (poles accumulating at 0) — gives `None`.
fn singular_polys_off_origin(arena: &mut Arena, g: ExprId, t: ExprId) -> Option<Vec<ExprId>> {
    let mut polys: Vec<ExprId> = Vec::new();
    let mut stack = vec![g];
    let mut seen = rustc_hash::FxHashSet::default();
    while let Some(id) = stack.pop() {
        if !seen.insert(id) || !walk::contains(arena, id, t) {
            continue;
        }
        let node = arena.node(id).clone();
        let ok = match node {
            ExprNode::Symbol(_) | ExprNode::Add(_) | ExprNode::Mul(_) | ExprNode::Neg(_) => true,
            ExprNode::Exp(_)
            | ExprNode::Sin(_)
            | ExprNode::Cos(_)
            | ExprNode::Sinh(_)
            | ExprNode::Cosh(_)
            | ExprNode::Erf(_)
            | ExprNode::Erfc(_) => true,
            ExprNode::Pow(b, x) => {
                let int_exp = arena.as_num(x).filter(|q| q.is_integer()).cloned();
                match int_exp {
                    Some(q) if !num_traits::Signed::is_negative(&q) => true,
                    Some(_) => {
                        let factors: Vec<ExprId> = match arena.node(b) {
                            ExprNode::Mul(fs) => fs.to_vec(),
                            _ => vec![b],
                        };
                        let mut all = true;
                        for f in factors {
                            if !walk::contains(arena, f, t)
                                || matches!(arena.node(f), ExprNode::Exp(_))
                            {
                                continue;
                            }
                            if calculus_util::poly_coeffs_symbolic(arena, f, t).is_some() {
                                polys.push(f);
                            } else {
                                all = false;
                                break;
                            }
                        }
                        all
                    }
                    // c^w = exp(w·ln c) is entire in w for a constant c ≠ 0.
                    None => {
                        !walk::contains(arena, b, t) && {
                            let pv = point_value(arena, b);
                            pv.nonzero
                        }
                    }
                }
            }
            _ => false,
        };
        if !ok {
            return None;
        }
        node.for_each_child(|c| stack.push(c));
    }
    Some(polys)
}

/// All distinct nonzero roots of the polynomials `polys` in `t`, or `None`
/// unless they are certainly all found.
///
/// Polynomials with rational coefficients are multiplied together and made
/// square-free, so a root shared by two factors or repeated is counted
/// once; the solver's distinct roots are complete exactly when their
/// number is the degree.  A polynomial with symbolic coefficients (at most
/// one) must have as many distinct roots found as its degree.
fn nonzero_roots(arena: &mut Arena, polys: &[ExprId], t: ExprId) -> Option<Vec<ExprId>> {
    use crate::poly::Poly;
    let mut rational = Poly::from_coeffs(vec![num_rational::Ratio::one()]);
    let mut symbolic: Vec<ExprId> = Vec::new();
    for &p in polys {
        let coeffs = calculus_util::poly_coeffs_symbolic(arena, p, t)?;
        // Divide out the root at 0: p̃(t) = p(t)/t^m.
        let m = coeffs.iter().position(|&c| !arena.is_zero_structural(c))?;
        let tail = &coeffs[m..];
        if tail.len() <= 1 {
            continue;
        }
        let nums: Option<Vec<num_rational::Ratio<BigInt>>> =
            tail.iter().map(|&c| arena.as_num(c).cloned()).collect();
        if let Some(nums) = nums {
            rational = rational.mul(&Poly::from_coeffs(nums));
        } else {
            let mut terms: Vec<ExprId> = Vec::with_capacity(tail.len());
            for (k, &c) in tail.iter().enumerate() {
                let e = arena.int(k as i64);
                let tk = arena.pow(t, e);
                terms.push(arena.mul(&[c, tk]));
            }
            let q = arena.add(&terms);
            if !symbolic.contains(&q) {
                symbolic.push(q);
            }
        }
    }
    if symbolic.len() > 1 {
        return None;
    }
    let mut roots: Vec<ExprId> = Vec::new();
    let mut solve_all = |arena: &mut Arena, q: ExprId, degree: usize| -> Option<()> {
        let found = crate::transforms::solve::solve(arena, q, t);
        if found.len() != degree {
            return None;
        }
        for s in found {
            let v = eval::eval(arena, s.value);
            if arena.is_zero_structural(v) || walk::contains(arena, v, t) {
                return None;
            }
            if !roots.contains(&v) {
                roots.push(v);
            }
        }
        Some(())
    };
    let sqf = rational.square_free_part();
    let degree = sqf.degree().unwrap_or(0);
    if degree > 0 {
        let q = crate::poly::polybridge::poly_to_expr(arena, &sqf, t);
        solve_all(arena, q, degree)?;
    }
    if let Some(&q) = symbolic.first() {
        let degree = calculus_util::poly_coeffs_symbolic(arena, q, t)?.len() - 1;
        solve_all(arena, q, degree)?;
    }
    Some(roots)
}

/// Residue at infinity: `Res_{z=∞} f = −Res_{t=0} f(1/t)/t²`.
pub(crate) fn residue_at_infinity(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
) -> Result<ExprId, SymplexError> {
    residue_at_infinity_with(arena, expr, var, true)
}

/// [`residue_at_infinity`] with the essential-singularity route enabled or
/// not (see [`residue_with`]).
fn residue_at_infinity_with(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    essential: bool,
) -> Result<ExprId, SymplexError> {
    if !matches!(arena.node(var), ExprNode::Symbol(_)) {
        return Err(failed("residue variable must be a symbol"));
    }
    let t = arena.symbol("__res_inf_t");
    let one = arena.one();
    let inv_t = arena.div(one, t);
    let f_inv = subs::subs(arena, expr, var, inv_t);
    let neg_two = arena.int(-2);
    let t_m2 = arena.pow(t, neg_two);
    let h = arena.mul(&[f_inv, t_m2]);
    let h = eval::eval(arena, h);
    let zero = arena.zero();
    let r = residue_with(arena, h, t, zero, essential)?;
    let neg = arena.neg(r);
    Ok(eval::eval(arena, neg))
}

/// `v / n!`
fn divide_by_factorial(arena: &mut Arena, v: ExprId, n: usize) -> ExprId {
    let mut fact = BigInt::one();
    for k in 2..=n {
        fact *= BigInt::from(k);
    }
    let f = arena.big_int(fact);
    arena.div(v, f)
}

/// Is `c` provably zero (structurally after `eval`, after full
/// simplification, or numerically when constant)?
fn is_provably_zero(arena: &mut Arena, c: ExprId) -> bool {
    let c = eval::eval(arena, c);
    if arena.is_zero_structural(c) {
        return true;
    }
    let expanded = crate::transforms::expand::expand(arena, c);
    let expanded = eval::eval(arena, expanded);
    if arena.is_zero_structural(expanded) {
        return true;
    }
    if walk::free_symbols(arena, expanded).is_empty() {
        // Constant: decide numerically (real or complex).
        if let Ok(s) = crate::transforms::evalf::evalf(arena, expanded, 20) {
            let s = s.trim();
            if s == "0" {
                return true;
            }
            if let Ok(v) = s.parse::<f64>() {
                return v.abs() < 1e-18;
            }
            return false;
        }
    }
    let opts = crate::simplify::simplify_engine::SimplifyOpts::default();
    let r = crate::simplify::simplify_engine::unified_simplify(arena, expanded, &opts);
    arena.is_zero_structural(r.expr)
}

/// What is known about a constant `v` (the value of a subexpression at the
/// expansion point).
struct PointValue {
    /// Provably nonzero (or, with free parameters, not structurally zero —
    /// the generic case).
    nonzero: bool,
    /// The complex value, when `v` is a finite number.
    value: Option<num_complex::Complex64>,
}

fn point_value(arena: &mut Arena, v: ExprId) -> PointValue {
    let v = eval::eval(arena, v);
    if arena.is_zero_structural(v) {
        return PointValue {
            nonzero: false,
            value: Some(num_complex::Complex64::new(0.0, 0.0)),
        };
    }
    if !walk::free_symbols(arena, v).is_empty() {
        // Symbolic parameter: generic position, as for every other
        // parameter-dependent residue.
        return PointValue {
            nonzero: true,
            value: None,
        };
    }
    match crate::transforms::evalf::evalf_complex64(arena, v) {
        Ok(z) if z.re.is_finite() && z.im.is_finite() => PointValue {
            nonzero: !is_provably_zero(arena, v),
            value: Some(z),
        },
        _ => PointValue {
            nonzero: false,
            value: None,
        },
    }
}

/// Is `e` provably analytic (holomorphic) in a neighbourhood of `t = 0`?
///
/// Every subexpression depending on `t` is checked with an explicit stack
/// (no recursion over the tree): sums, products and non-negative integer
/// powers of analytic functions are analytic, as are the entire functions
/// (`exp`, `sin`, `cos`, `sinh`, `cosh`, `erf`, `erfc`) of analytic
/// arguments.  Negative powers need a nonzero base at 0; fractional powers,
/// logarithms and the inverse functions additionally need the argument at 0
/// off their branch cuts; `tan`/`tanh`/`Γ`/`ψ` need to avoid their poles.
/// Anything else (`|·|`, `re`, `conj`, `sign`, `H`, `floor`, …) is not
/// certified.  With symbolic parameters the value at 0 is taken in generic
/// position (nonzero unless structurally zero).
pub(crate) fn analytic_at_zero(arena: &mut Arena, e: ExprId, t: ExprId) -> bool {
    use num_complex::Complex64;
    let zero = arena.zero();
    let mut stack = vec![e];
    let mut seen = rustc_hash::FxHashSet::default();
    // Value of a child at t = 0.
    let at0 = |arena: &mut Arena, c: ExprId| -> PointValue {
        let v = subs::subs(arena, c, t, zero);
        point_value(arena, v)
    };
    // Off the negative real axis (principal branch cut of ln / powers).
    let off_negative_axis = |z: Option<Complex64>| match z {
        Some(z) => z.re > 0.0 || z.im.abs() > 1e-12 * z.norm(),
        None => true,
    };
    while let Some(id) = stack.pop() {
        if !seen.insert(id) || !walk::contains(arena, id, t) {
            continue;
        }
        let node = arena.node(id).clone();
        let ok = match node {
            ExprNode::Symbol(_) => true,
            ExprNode::Add(_) | ExprNode::Mul(_) | ExprNode::Neg(_) => true,
            ExprNode::Exp(_)
            | ExprNode::Sin(_)
            | ExprNode::Cos(_)
            | ExprNode::Sinh(_)
            | ExprNode::Cosh(_)
            | ExprNode::Erf(_)
            | ExprNode::Erfc(_) => true,
            ExprNode::Pow(b, x) => {
                let nonneg_int = arena
                    .as_num(x)
                    .is_some_and(|q| q.is_integer() && !num_traits::Signed::is_negative(q));
                let neg_int = arena
                    .as_num(x)
                    .is_some_and(|q| q.is_integer() && num_traits::Signed::is_negative(q));
                if nonneg_int {
                    true
                } else if !walk::contains(arena, b, t) {
                    // c^x(t) = exp(x(t)·ln c): entire in x for c ≠ 0.
                    at0(arena, b).nonzero
                } else {
                    let pv = at0(arena, b);
                    pv.nonzero && (neg_int || off_negative_axis(pv.value))
                }
            }
            ExprNode::Ln(a) => {
                let pv = at0(arena, a);
                pv.nonzero && off_negative_axis(pv.value)
            }
            ExprNode::Tan(a) => {
                let c = arena.cos(a);
                at0(arena, c).nonzero
            }
            ExprNode::Tanh(a) => {
                let c = arena.cosh(a);
                at0(arena, c).nonzero
            }
            ExprNode::Asin(a) | ExprNode::Acos(a) | ExprNode::Atanh(a) => {
                // cuts: real |x| ≥ 1
                match at0(arena, a).value {
                    Some(z) => z.im.abs() > 1e-12 * z.norm() || z.re.abs() < 1.0,
                    None => true,
                }
            }
            ExprNode::Atan(a) | ExprNode::Asinh(a) => {
                // cuts: imaginary |y| ≥ 1
                match at0(arena, a).value {
                    Some(z) => z.re.abs() > 1e-12 * z.norm() || z.im.abs() < 1.0,
                    None => true,
                }
            }
            ExprNode::Acosh(a) => {
                // cut: real x ≤ 1
                match at0(arena, a).value {
                    Some(z) => z.im.abs() > 1e-12 * z.norm() || z.re > 1.0,
                    None => true,
                }
            }
            ExprNode::Gamma(a) | ExprNode::Digamma(a) => match at0(arena, a).value {
                Some(z) => {
                    !(z.im.abs() <= 1e-12 * z.norm().max(1.0)
                        && z.re <= 0.5
                        && (z.re - z.re.round()).abs() < 1e-12)
                }
                None => true,
            },
            _ => false,
        };
        if !ok {
            return false;
        }
        arena.node(id).clone().for_each_child(|c| stack.push(c));
    }
    true
}

/// Is `v` a finite value free of `t`, infinities and formal nodes?
fn is_finite_constant(arena: &mut Arena, v: ExprId, t: ExprId) -> bool {
    if walk::contains(arena, v, t) || walk::has_unevaluated(arena, v) {
        return false;
    }
    let mut stack = vec![v];
    let mut seen = rustc_hash::FxHashSet::default();
    while let Some(id) = stack.pop() {
        if !seen.insert(id) {
            continue;
        }
        match arena.node(id) {
            ExprNode::Infinity
            | ExprNode::NegInfinity
            | ExprNode::ComplexInfinity
            | ExprNode::NaN => {
                return false;
            }
            ExprNode::Ln(inner) if arena.is_zero_structural(*inner) => return false,
            ExprNode::Pow(base, exp)
                if arena.is_zero_structural(*base)
                    && arena
                        .as_num(*exp)
                        .is_some_and(num_traits::Signed::is_negative) =>
            {
                return false;
            }
            node => node.for_each_child(|c| stack.push(c)),
        }
    }
    if walk::free_symbols(arena, v).is_empty() && arena.as_num(v).is_none() {
        // Constant: make sure it evaluates to a finite (possibly complex) number.
        return match crate::transforms::evalf::evalf(arena, v, 16) {
            Ok(s) => !s.contains("inf") && !s.contains("nan") && !s.contains("NaN"),
            Err(_) => false,
        };
    }
    true
}

#[cfg(test)]
mod tests {
    use super::*;

    fn show(a: &Arena, id: ExprId) -> String {
        a.display(id).to_string()
    }

    #[test]
    fn simple_pole() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let f = a.pow(z, a.neg_one());
        let zero = a.zero();
        let r = residue(&mut a, f, z, zero).unwrap();
        assert_eq!(show(&a, r), "1");
    }

    #[test]
    fn double_pole_exp() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let m2 = a.int(-2);
        let zm2 = a.pow(z, m2);
        let ez = a.exp(z);
        let f = a.mul(&[zm2, ez]);
        let zero = a.zero();
        let r = residue(&mut a, f, z, zero).unwrap();
        assert_eq!(show(&a, r), "1");
    }

    #[test]
    fn triple_pole_zero_residue() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let one = a.one();
        let zm1 = a.sub(z, one);
        let m3 = a.int(-3);
        let f = a.pow(zm1, m3);
        let r = residue(&mut a, f, z, one).unwrap();
        assert_eq!(show(&a, r), "0");
    }

    #[test]
    fn cos_over_z_cubed() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let c = a.cos(z);
        let m3 = a.int(-3);
        let zm3 = a.pow(z, m3);
        let f = a.mul(&[c, zm3]);
        let zero = a.zero();
        let r = residue(&mut a, f, z, zero).unwrap();
        assert_eq!(show(&a, r), "-1/2");
    }

    #[test]
    fn double_pole_at_i() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let two = a.int(2);
        let z2 = a.pow(z, two);
        let one = a.one();
        let d = a.add(&[z2, one]);
        let m2 = a.int(-2);
        let f = a.pow(d, m2);
        let i = a.i_unit();
        let r = residue(&mut a, f, z, i).unwrap();
        // −i/4
        let s = show(&a, r);
        assert!(s == "-1/4*I" || s == "-1/4*i", "{s}");
    }

    #[test]
    fn analytic_point_is_zero() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let two = a.int(2);
        let f = a.pow(z, two);
        let zero = a.zero();
        let r = residue(&mut a, f, z, zero).unwrap();
        assert_eq!(show(&a, r), "0");
    }

    #[test]
    fn essential_singularity_by_the_residue_theorem() {
        // Res(e^{1/z}, 0) = 1, the coefficient of 1/z in Σ z^{-n}/n! (this
        // was refused; before that it was 0).
        let mut a = Arena::new();
        let z = a.symbol("z");
        let inv = a.pow(z, a.neg_one());
        let f = a.exp(inv);
        let zero = a.zero();
        let r = residue(&mut a, f, z, zero).unwrap();
        assert_eq!(show(&a, r), "1");
        // e^z·e^{1/z}: ∞ is essential too (the residue is I₁(2)): refused.
        let ez = a.exp(z);
        let g = a.mul(&[ez, f]);
        assert!(residue(&mut a, g, z, zero).is_err());
    }

    #[test]
    fn transcendental_denominator() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let s = a.sin(z);
        let f = a.pow(s, a.neg_one());
        let zero = a.zero();
        let r = residue(&mut a, f, z, zero).unwrap();
        assert_eq!(show(&a, r), "1");
    }

    #[test]
    fn at_infinity() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let f = a.pow(z, a.neg_one());
        let r = residue_at_infinity(&mut a, f, z).unwrap();
        assert_eq!(show(&a, r), "-1");
    }
}