symplex 0.30.0

//! Taylor / Laurent / asymptotic series expansion.
//!
//! [`series`] computes the truncated expansion of an expression around a
//! point (Maclaurin when the point is `0`, asymptotic when the point is
//! `±∞` via [`series_at_infinity`]); [`series_dir`] from one side.
//!
//! # Contract
//!
//! `series(f, x, x0, n)` is the **asymptotic expansion of `f` as `x → x0`
//! through real values** — two-sided unless a direction is given
//! ([`series_dir`]; at `±∞` the side is implied) — in the same sense as
//! `limit` and the Gruntz algorithm: the result `S` contains every term of
//! exponent `< n` in `x − x0` (in `1/x` at `±∞`), each exact, and
//! `f − S = O((x − x0)ⁿ)` as `x` approaches `x0` along the reals from the
//! side(s) in question.  It is **not** a Laurent series over the complex
//! plane:
//!
//! * `exp(−1/x²)` at `0` is `0` to every order (all real Taylor
//!   coefficients vanish; SymPy: `series(exp(-1/x**2), x, 0, 4) = O(x**4)`),
//!   although `0` is an essential singularity;
//! * `exp(−1/x)` at `0` has no two-sided expansion (it is `O(xⁿ)` from the
//!   right and unbounded from the left) and is refused; from the right it
//!   is `0`, from the left it is refused (SymPy's default `dir='+'` gives
//!   `O(x**4)`, `dir='-'` returns `exp(-1/x)` unexpanded);
//! * a two-sided kink (`|x|`) is refused unless both sides agree
//!   (`cos|x| = cos x`);
//! * expansions from the right and at `±∞` are *log-extended*: a
//!   coefficient may contain `ln(x − x0)` (`ln x` at `±∞`), as in
//!   Stirling's series `ln Γ(x) = x ln x − x − ½ ln x + ½ ln 2π +
//!   1/(12x) + …` or `Ei(x) = γ + ln x + x + …` for `x → 0⁺`; the
//!   remainder is then `O((x − x0)ⁿ·|ln(x − x0)|ᵏ)`.  Two-sided
//!   expansions refuse logarithmic singularities.
//!
//! When no such expansion exists, or it cannot be established (a
//! fractional power two-sided, an oscillation, an unknown function whose
//! derivatives are singular at the point, the Stieltjes constants that
//! `ζ(x)` needs beyond `1/(x − 1) + γ`), the result is `Err` and the public
//! API keeps a formal `Series` node.
//!
//! Consumers that need the complex-analytic (Laurent) expansion —
//! residues, Laurent coefficients, inverse transforms — must not take it
//! from [`series`] for a function that may have an essential singularity at
//! the point: [`laurent_series`] runs the engine without limits and refuses
//! an essential singularity, and residues have their own structural route
//! ([`residue`](crate::calculus::residue)).
//!
//! # Algorithm
//!
//! Expansion is performed by a **truncated Laurent-series arithmetic engine**
//! that walks the expression DAG bottom-up (explicit post-order, no
//! recursion) and combines the series of the children:
//!
//! * `Add` / `Mul` / integer `Pow` — exact truncated arithmetic (poles are
//!   represented by a negative leading exponent, so `sin x / x` needs no
//!   special handling).
//! * `exp`, `sin`, `cos`, `sinh`, `cosh`, `ln`, `atan`, `atanh`, `asin`,
//!   `asinh`, `tan`, `tanh`, `erf`, `LambertW`, `(1+u)^α` — composed from
//!   closed-form Maclaurin coefficients (Bernoulli numbers for `tan`/`tanh`,
//!   central binomials for `asin`, `(−n)^{n−1}/n!` for `W`) rather than
//!   repeated differentiation, so high orders stay fast.
//! * `Γ`, `ψ`, `ψ⁽ᵐ⁾` at their poles, `ζ` at `1`, `ln Γ` and `ψ` at `+∞`
//!   (Stirling), `Ei`, `Ci`, `Chi`, `li`, `Si`, `Shi`, `K₀`, `Y₀` at the
//!   zeros of their arguments — closed forms at the points where
//!   differentiation has nothing to evaluate (see *Special functions at
//!   their singular points* below).
//! * Any other `var`-dependent sub-expression falls back to Taylor
//!   coefficients by differentiation, evaluated at the expansion point.
//!
//! Fractional powers of a series with a zero constant term (Puiseux
//! expansions such as `√x·sin x`) and logarithmic singularities are
//! rejected — the caller then keeps an unevaluated `Series` node instead of
//! producing a wrong polynomial.  Such rejections are *definite*: the
//! differentiation fallback is never tried for them, because the low-order
//! derivatives of `x^(5/2)` or `|x²|` all vanish at `0` and would silently
//! yield the wrong polynomial `0`.
//!
//! `|g|` is expanded as `±g` when the sign of `g` near the point is known:
//! always when the leading exponent of `g` is even (`|x²| = x²`,
//! `|x − 1| = 1 − x`), and one-sidedly otherwise.  For an odd leading
//! exponent the two one-sided expansions of the *whole* expression are
//! compared and accepted only when they agree (`cos|x| = cos x`); `e^|x|`
//! and `|sin x|` have no two-sided expansion and are rejected.
//!
//! # Design
//!
//! The result is a plain expression (polynomial in `(x − a)`, possibly with
//! negative powers).  No `O(·)` term is appended — the truncation order is
//! implicit in the `order` parameter: all terms with exponent `< order` are
//! present and exact.
//!
//! A coefficient that is a sum of constants is tested for a hidden zero
//! (`asinh 2 − ln(2 + √5)`) numerically as it forms
//! ([`settle_constant`]): the valuation of a series, and so every division
//! by it, depends on knowing which leading coefficients vanish.

use num_bigint::BigInt;
use num_rational::Ratio;
use num_traits::{One, Signed, ToPrimitive, Zero};
use rustc_hash::FxHashMap;

use crate::base::arena::Arena;
use crate::base::assumptions::Props;
use crate::base::bernoulli::bernoulli;
use crate::base::combinatorics::{binomial, factorial};
use crate::base::errors::SymplexError;
use crate::base::node::{ExprId, ExprNode};
use crate::base::numeric::Q;
use crate::base::walk;
use crate::calculus::calculus_util::settle_constant;
use crate::transforms::{eval, subs};

/// Maximum number of sub-expressions expanded by differentiation before the
/// engine gives up on per-node fallbacks (the root is always tried).
const MAX_FALLBACK_NODES: usize = 4;

/// Largest derivative (tree size) the differentiation fallback evaluates.
const MAX_FALLBACK_SIZE: usize = 1_500;

/// Maximum integer exponent expanded by repeated multiplication; larger
/// exponents use the binomial series.
const MAX_INT_POWER: i64 = 64;

/// Minimum internal working precision of the engine.
const MIN_WORKING_ORDER: i64 = 4;

/// Number of precision-escalation rounds before giving up on a deep pole.
const MAX_PRECISION_ATTEMPTS: usize = 3;

// ═══════════════════════════════════════════════════════════════════════════
// Public entry points
// ═══════════════════════════════════════════════════════════════════════════

/// Compute the series of `expr` in `var` around `point` with all terms of
/// exponent `< order` (in `var − point`): the two-sided real asymptotic
/// expansion of the [contract](self#contract).
///
/// For `point = 0` this is the Maclaurin series; for `point = ±∞` an
/// asymptotic expansion in `1/var` (see [`series_at_infinity`]).  Poles at
/// the expansion point produce negative powers.
///
/// Returns `Err` when no such expansion exists or cannot be established
/// (fractional-power or logarithmic singularity, the two sides differing,
/// an unsupported sub-expression whose derivatives are singular at the
/// point).
pub(crate) fn series(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
    order: u32,
) -> Result<ExprId, SymplexError> {
    let mode = Mode {
        side: Side::Both,
        log_var: None,
        limits: true,
        partial: false,
    };
    series_with_pole_retry(arena, expr, var, point, order, mode)
}

/// [`series`] as `var → point` from one side: [`Direction::Right`]
/// (`var > point`, log-extended: the result may contain `ln(var − point)`,
/// `Ei(x) = γ + ln x + x + x²/4 + …`) or [`Direction::Left`];
/// [`Direction::Both`] is [`series`].  At `±∞` the direction is implied.
///
/// `exp(−1/x)` at `0` is `0` to every order from the right and has no
/// expansion from the left (nor a two-sided one).
///
/// [`Direction::Right`]: crate::calculus::limit::Direction::Right
/// [`Direction::Left`]: crate::calculus::limit::Direction::Left
/// [`Direction::Both`]: crate::calculus::limit::Direction::Both
pub(crate) fn series_dir(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
    order: u32,
    dir: crate::calculus::limit::Direction,
) -> Result<ExprId, SymplexError> {
    use crate::calculus::limit::Direction;
    let mode = match dir {
        Direction::Both => return series(arena, expr, var, point, order),
        Direction::Right => Mode {
            side: Side::Above,
            log_var: Some(arena.symbol(LOG_PLACEHOLDER)),
            limits: true,
            partial: false,
        },
        Direction::Left => Mode {
            side: Side::Below,
            log_var: None,
            limits: true,
            partial: false,
        },
    };
    series_with_pole_retry(arena, expr, var, point, order, mode)
}

/// Stand-in for `ln t` during a log-extended expansion (replaced before the
/// result is returned).
const LOG_PLACEHOLDER: &str = "__series_ln_t";

/// [`series_in_mode`], retried as `(x − a)⁻ᵏ·series((x − a)ᵏ·f)` for
/// `k = 1..=5` when the direct expansion fails: a pole of a function the
/// engine only knows through the differentiation fallback (`sinh(ln x) =
/// (x − 1/x)/2`) is invisible to Taylor coefficients.  For the real
/// asymptotic expansion this is exact: `xᵏ·f = S + O(xⁿ⁺ᵏ)` gives
/// `f = S/xᵏ + O(xⁿ)` for `x ≠ 0`.  (Before, this retry lived only in
/// [`laurent_series`], reached as `series`'s fallback, which is now
/// restricted to meromorphic functions.)
fn series_with_pole_retry(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
    order: u32,
    mode: Mode,
) -> Result<ExprId, SymplexError> {
    let err = match series_in_mode(arena, expr, var, point, order, mode) {
        Ok(s) => return Ok(s),
        Err(e) => e,
    };
    if matches!(
        arena.node(point),
        ExprNode::Infinity | ExprNode::NegInfinity
    ) || !matches!(arena.node(var), ExprNode::Symbol(_))
        || order == 0
    {
        return Err(err);
    }
    let x_minus_a = if arena.is_zero_structural(point) {
        var
    } else {
        arena.sub(var, point)
    };
    for k in 1u32..=5 {
        let k_id = arena.int(k as i64);
        let multiplier = arena.pow(x_minus_a, k_id);
        let modified = arena.mul(&[expr, multiplier]);
        if let Ok(ts) = series_in_mode(arena, modified, var, point, order + k, mode) {
            let neg_k = arena.int(-(k as i64));
            let divisor = arena.pow(x_minus_a, neg_k);
            let result = arena.mul(&[ts, divisor]);
            let result = crate::transforms::expand::expand(arena, result);
            return Ok(eval::eval(arena, result));
        }
    }
    Err(err)
}

/// The expansion at a finite `point` (or `±∞`) in the given [`Mode`].
fn series_in_mode(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
    order: u32,
    mode: Mode,
) -> Result<ExprId, SymplexError> {
    if order == 0 {
        return Ok(arena.zero);
    }
    match arena.node(point) {
        ExprNode::Infinity => return series_at_infinity(arena, expr, var, order, false),
        ExprNode::NegInfinity => return series_at_infinity(arena, expr, var, order, true),
        _ => {}
    }
    if !matches!(arena.node(var), ExprNode::Symbol(_)) {
        return Err(SymplexError::InvalidArgument {
            operation: "series",
            reason: "expansion variable must be a symbol".into(),
        });
    }
    let at_zero = arena.is_zero_structural(point);
    // Shift so that the expansion point becomes 0:  f(x) = f(a + t).
    let shifted = if at_zero {
        expr
    } else {
        let t_plus_a = arena.add(&[var, point]);
        subs::subs(arena, expr, var, t_plus_a)
    };
    let ts = expand_with_mode(arena, shifted, var, order as i64, mode)?;
    let mut poly = ts.to_expr(arena, var, order as i64);
    if let Some(l) = mode.log_var {
        let ln_t = arena.ln(var);
        poly = subs::subs(arena, poly, l, ln_t);
        poly = eval::eval(arena, poly);
    }
    let result = if at_zero {
        poly
    } else {
        let x_minus_a = arena.sub(var, point);
        let back = subs::subs(arena, poly, var, x_minus_a);
        eval::eval(arena, back)
    };
    if mode.log_var.is_some()
        && (contains_singular_atom(arena, result) || walk::has_unevaluated(arena, result))
    {
        return Err(SymplexError::ComputationFailed {
            operation: "series",
            reason: "a coefficient of the one-sided expansion is singular".into(),
        });
    }
    Ok(result)
}

/// Asymptotic expansion of `expr` as `var → +∞` (or `−∞` when `negative`),
/// with all terms of exponent `< order` in `1/var`.
///
/// Substitutes `var = ±1/t`, expands at `t = 0`, and substitutes back.
pub(crate) fn series_at_infinity(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    order: u32,
    negative: bool,
) -> Result<ExprId, SymplexError> {
    if order == 0 {
        return Ok(arena.zero);
    }
    let t = arena.symbol("_t");
    let one = arena.one;
    let inv_t = arena.div(one, t);
    let inv_t = if negative { arena.neg(inv_t) } else { inv_t };
    let in_t = subs::subs(arena, expr, var, inv_t);
    // t = ±1/x → 0⁺ only, so fractional powers of t^(even) are
    // single-valued, and `ln t` is real: the expansion is log-extended, with
    // `ln t = −ln(±x)` (Stirling's `x ln x − x − ½ ln x + …`).
    let log_t = arena.symbol(LOG_PLACEHOLDER);
    let ts = expand_log_extended(arena, in_t, t, order as i64, log_t, true)?;
    let poly = ts.to_expr(arena, t, order as i64);
    let ln_arg = if negative { arena.neg(var) } else { var };
    let ln_x = arena.ln(ln_arg);
    let minus_ln_x = arena.neg(ln_x);
    let poly = subs::subs(arena, poly, log_t, minus_ln_x);
    let inv_x = arena.div(one, var);
    let inv_x = if negative { arena.neg(inv_x) } else { inv_x };
    let back = subs::subs(arena, poly, t, inv_x);
    let result = eval::eval(arena, back);
    // Never return a bogus expansion: a coefficient that is infinite,
    // undefined or an unevaluated limit means the expansion failed.
    if contains_singular_atom(arena, result) || walk::has_unevaluated(arena, result) {
        return Err(SymplexError::ComputationFailed {
            operation: "series_at_infinity",
            reason: format!(
                "no asymptotic expansion in powers of 1/{}: a coefficient is singular",
                arena.display(var)
            ),
        });
    }
    Ok(result)
}

/// `true` if the expression contains `∞`, `−∞`, `zoo`, or `NaN` anywhere.
fn contains_singular_atom(arena: &Arena, id: ExprId) -> bool {
    walk::contains(arena, id, arena.infinity)
        || walk::contains(arena, id, arena.neg_infinity)
        || walk::contains(arena, id, arena.complex_infinity)
        || walk::contains(arena, id, arena.nan)
}

/// Compute a Laurent series expansion of `expr` in `var` around `point`.
///
/// Unlike [`series`] (a real asymptotic expansion) this is the
/// complex-analytic Laurent series, so it exists only where `expr` is
/// meromorphic at the point: the engine runs without taking limits in its
/// differentiation fallback, and an essential singularity is refused —
/// before, `exp(−1/x²)` had the "Laurent series" `0` (its real Taylor
/// series; the Laurent series has infinitely many negative powers).  As a
/// last resort this multiplies by `(x − a)^k`, `k = 1..=5`, and divides
/// back.
pub(crate) fn laurent_series(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    point: ExprId,
    order: u32,
) -> Result<ExprId, SymplexError> {
    let mode = Mode {
        side: Side::Both,
        log_var: None,
        limits: false,
        partial: false,
    };
    let series = |arena: &mut Arena, e: ExprId, order: u32| {
        if matches!(
            arena.node(point),
            ExprNode::Infinity | ExprNode::NegInfinity
        ) {
            return Err(SymplexError::ComputationFailed {
                operation: "laurent_series",
                reason: "a Laurent series is taken at a finite point".into(),
            });
        }
        series_in_mode(arena, e, var, point, order, mode)
    };
    if let Ok(ts) = series(arena, expr, order) {
        return Ok(ts);
    }
    let x_minus_a = if arena.is_zero_structural(point) {
        var
    } else {
        arena.sub(var, point)
    };
    for k in 1u32..=5 {
        let k_id = arena.int(k as i64);
        let multiplier = arena.pow(x_minus_a, k_id);
        let modified = arena.mul(&[expr, multiplier]);
        if let Ok(ts) = series(arena, modified, order + k) {
            let neg_k = arena.int(-(k as i64));
            let divisor = arena.pow(x_minus_a, neg_k);
            let result = arena.mul(&[ts, divisor]);
            let result = crate::transforms::expand::expand(arena, result);
            return Ok(eval::eval(arena, result));
        }
    }
    Err(SymplexError::ComputationFailed {
        operation: "laurent_series",
        reason: "could not determine pole order (tried up to order 5)".into(),
    })
}

// ═══════════════════════════════════════════════════════════════════════════
// Truncated Laurent series
// ═══════════════════════════════════════════════════════════════════════════

/// `Σ_{i} coeffs[i] · x^(shift + i)`, exact for all exponents `< known`.
///
/// Invariant: `coeffs.len() == (known − shift) as usize`.
#[derive(Clone, Debug)]
pub(crate) struct TSeries {
    shift: i64,
    known: i64,
    coeffs: Vec<ExprId>,
}

impl TSeries {
    /// Exponent bound: coefficients are exact for all exponents `< known`.
    pub(crate) fn known(&self) -> i64 {
        self.known
    }

    /// Lowest stored exponent (negative for Laurent series).
    pub(crate) fn shift(&self) -> i64 {
        self.shift
    }

    /// Coefficient of `x^e`, zero outside the stored range (public alias of
    /// [`coeff_at`](Self::coeff_at) for other modules).
    pub(crate) fn coefficient(&self, arena: &Arena, e: i64) -> ExprId {
        self.coeff_at(arena, e)
    }

    fn zero(arena: &Arena, known: i64) -> Self {
        TSeries {
            shift: 0,
            known,
            coeffs: vec![arena.zero; known.max(0) as usize],
        }
    }

    fn constant(arena: &Arena, c: ExprId, known: i64) -> Self {
        let mut s = Self::zero(arena, known);
        if known > 0 {
            s.coeffs[0] = c;
        }
        s
    }

    fn var(arena: &Arena, known: i64) -> Self {
        let mut s = Self::zero(arena, known);
        if known > 1 {
            s.coeffs[1] = arena.one;
        }
        s
    }

    /// Coefficient of `x^e` (zero outside the stored range).
    fn coeff_at(&self, arena: &Arena, e: i64) -> ExprId {
        if e < self.shift || e >= self.known {
            arena.zero
        } else {
            self.coeffs[(e - self.shift) as usize]
        }
    }

    /// Exponent of the first structurally non-zero coefficient.
    fn leading_exponent(&self, arena: &Arena) -> Option<i64> {
        self.coeffs
            .iter()
            .position(|&c| !arena.is_zero_structural(c))
            .map(|i| self.shift + i as i64)
    }

    /// Drop leading structural zeros so that `shift` is the true valuation.
    fn normalized(mut self, arena: &Arena) -> Self {
        let lead = self
            .coeffs
            .iter()
            .position(|&c| !arena.is_zero_structural(c))
            .unwrap_or(self.coeffs.len());
        if lead > 0 {
            self.coeffs.drain(0..lead);
            self.shift += lead as i64;
        }
        self
    }

    /// Do `a` and `b` agree on every coefficient of exponent `< order`?
    /// (Structural comparison of the evaluated coefficients.)
    fn same(arena: &Arena, a: &TSeries, b: &TSeries, order: i64) -> bool {
        let lo = a.shift.min(b.shift);
        let hi = a.known.min(b.known).min(order);
        (lo..hi).all(|e| a.coeff_at(arena, e) == b.coeff_at(arena, e))
    }

    /// Restrict to exponents `< n`.
    fn truncate_known(mut self, n: i64) -> Self {
        if n < self.known {
            let keep = (n - self.shift).max(0) as usize;
            self.coeffs.truncate(keep);
            self.known = n;
            if self.coeffs.is_empty() {
                self.shift = n;
            }
        }
        self
    }

    fn add(arena: &mut Arena, a: &TSeries, b: &TSeries) -> TSeries {
        let shift = a.shift.min(b.shift);
        let known = a.known.min(b.known);
        let mut coeffs = Vec::with_capacity((known - shift).max(0) as usize);
        for e in shift..known {
            let ca = a.coeff_at(arena, e);
            let cb = b.coeff_at(arena, e);
            let s = arena.add(&[ca, cb]);
            coeffs.push(
                if arena.is_zero_structural(ca) || arena.is_zero_structural(cb) {
                    eval::eval(arena, s)
                } else {
                    settle_constant(arena, s)
                },
            );
        }
        TSeries {
            shift,
            known,
            coeffs,
        }
    }

    fn scale(arena: &mut Arena, a: &TSeries, c: ExprId) -> TSeries {
        let coeffs = a
            .coeffs
            .iter()
            .map(|&x| {
                let p = arena.mul(&[c, x]);
                eval::eval(arena, p)
            })
            .collect();
        TSeries {
            shift: a.shift,
            known: a.known,
            coeffs,
        }
    }

    fn neg(arena: &mut Arena, a: &TSeries) -> TSeries {
        let m1 = arena.neg_one;
        Self::scale(arena, a, m1)
    }

    fn mul(arena: &mut Arena, a: &TSeries, b: &TSeries) -> TSeries {
        let a = a.clone().normalized(arena);
        let b = b.clone().normalized(arena);
        let shift = a.shift + b.shift;
        let known = (a.known + b.shift).min(b.known + a.shift);
        let len = (known - shift).max(0) as usize;
        let mut coeffs = Vec::with_capacity(len);
        for idx in 0..len {
            let mut terms = Vec::new();
            for (i, &ca) in a.coeffs.iter().enumerate() {
                if i > idx {
                    break;
                }
                let j = idx - i;
                if j >= b.coeffs.len() {
                    continue;
                }
                let cb = b.coeffs[j];
                if arena.is_zero_structural(ca) || arena.is_zero_structural(cb) {
                    continue;
                }
                terms.push(arena.mul(&[ca, cb]));
            }
            coeffs.push(match terms.len() {
                0 => arena.zero,
                1 => eval::eval(arena, terms[0]),
                _ => {
                    let s = arena.add(&terms);
                    settle_constant(arena, s)
                }
            });
        }
        TSeries {
            shift,
            known,
            coeffs,
        }
    }

    /// `1/a`.  Returns `None` if `a` is zero to the known precision.
    fn inverse(arena: &mut Arena, a: &TSeries) -> Option<TSeries> {
        let a = a.clone().normalized(arena);
        let v = a.leading_exponent(arena)?;
        let c0 = a.coeffs[0];
        let rel_known = a.known - v; // relative precision of 1 + w
        // w_i = a_{v+i}/c0 for i ≥ 1
        let inv_c0 = {
            let m1 = arena.neg_one;
            let p = arena.pow(c0, m1);
            eval::eval(arena, p)
        };
        let mut w: Vec<ExprId> = Vec::with_capacity(rel_known.max(0) as usize);
        w.push(arena.one);
        for i in 1..rel_known {
            let ai = a.coeff_at(arena, v + i);
            let p = arena.mul(&[ai, inv_c0]);
            w.push(eval::eval(arena, p));
        }
        // b_0 = 1, b_n = −Σ_{i=1}^{n} w_i b_{n−i}
        let mut b: Vec<ExprId> = Vec::with_capacity(w.len());
        b.push(arena.one);
        for n in 1..w.len() {
            let mut terms = Vec::new();
            for i in 1..=n {
                if arena.is_zero_structural(w[i]) || arena.is_zero_structural(b[n - i]) {
                    continue;
                }
                terms.push(arena.mul(&[w[i], b[n - i]]));
            }
            let s = match terms.len() {
                0 => arena.zero,
                1 => terms[0],
                _ => arena.add(&terms),
            };
            let ns = arena.neg(s);
            b.push(settle_constant(arena, ns));
        }
        let coeffs = b
            .iter()
            .map(|&x| {
                let p = arena.mul(&[inv_c0, x]);
                eval::eval(arena, p)
            })
            .collect();
        Some(TSeries {
            shift: -v,
            known: -v + rel_known,
            coeffs,
        })
    }

    /// `a^n` for integer `n`.
    fn pow_int(arena: &mut Arena, a: &TSeries, n: i64) -> Option<TSeries> {
        if n == 0 {
            return Some(Self::constant(arena, arena.one, a.known.max(1)));
        }
        if n.abs() > MAX_INT_POWER {
            return None;
        }
        let base = if n < 0 {
            Self::inverse(arena, a)?
        } else {
            a.clone()
        };
        let mut acc = base.clone();
        for _ in 1..n.abs() {
            acc = Self::mul(arena, &acc, &base);
        }
        Some(acc)
    }

    /// Remove the constant term, returning `(u0, w)` with `w = a − u0` (valuation ≥ 1).
    fn split_constant(&self, arena: &Arena) -> (ExprId, TSeries) {
        let u0 = self.coeff_at(arena, 0);
        let mut w = self.clone();
        if 0 >= w.shift && 0 < w.known {
            w.coeffs[(-w.shift) as usize] = arena.zero;
        }
        (u0, w.normalized(arena))
    }

    /// `Σ_n f_n · w^n` for a series `w` with valuation ≥ 1.
    fn compose(arena: &mut Arena, f: &dyn Fn(&mut Arena, usize) -> ExprId, w: &TSeries) -> TSeries {
        let known = w.known;
        if known <= 0 {
            return TSeries {
                shift: known,
                known,
                coeffs: Vec::new(),
            };
        }
        let mut acc = Self::constant(arena, arena.zero, known);
        let f0 = f(arena, 0);
        acc.coeffs[0] = f0;
        if w.coeffs.iter().all(|&c| arena.is_zero_structural(c)) {
            return acc;
        }
        let mut p = w.clone().normalized(arena);
        let mut n = 1usize;
        loop {
            if p.shift >= known || p.coeffs.is_empty() {
                break;
            }
            let fn_ = f(arena, n);
            if !arena.is_zero_structural(fn_) {
                let term = Self::scale(arena, &p, fn_);
                acc = Self::add(arena, &acc, &term);
            }
            n += 1;
            if n > known as usize + 1 {
                break;
            }
            p = Self::mul(arena, &p, w).normalized(arena);
        }
        acc
    }

    /// Build the expression `Σ coeffs[i] x^(shift+i)` for exponents `< order`.
    fn to_expr(&self, arena: &mut Arena, var: ExprId, order: i64) -> ExprId {
        let mut terms = Vec::new();
        for (i, &c) in self.coeffs.iter().enumerate() {
            let e = self.shift + i as i64;
            if e >= order {
                break;
            }
            if arena.is_zero_structural(c) {
                continue;
            }
            let term = if e == 0 {
                c
            } else if e == 1 {
                arena.mul(&[c, var])
            } else {
                let ee = arena.int(e);
                let xp = arena.pow(var, ee);
                arena.mul(&[c, xp])
            };
            terms.push(term);
        }
        let s = match terms.len() {
            0 => arena.zero,
            1 => terms[0],
            _ => arena.add(&terms),
        };
        eval::eval(arena, s)
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Closed-form Maclaurin coefficients
// ═══════════════════════════════════════════════════════════════════════════

fn rat_expr(arena: &mut Arena, r: Q) -> ExprId {
    let nid = arena.intern_num(r);
    arena.intern(ExprNode::Num(nid))
}

fn rat_i(n: i64) -> Q {
    Ratio::from_integer(BigInt::from(n))
}

/// The central binomial coefficient `C(2n, n)`.
fn central_binomial(n: u64) -> BigInt {
    binomial(2 * n, n)
}

fn pow_rat_i(r: &Q, n: i64) -> Q {
    let mut acc = Q::one();
    let base = if n < 0 { Q::one() / r } else { r.clone() };
    for _ in 0..n.unsigned_abs() {
        acc *= &base;
    }
    acc
}

/// Elementary functions with closed-form Maclaurin coefficients.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub(crate) enum FnKind {
    Exp,
    Sin,
    Cos,
    Sinh,
    Cosh,
    /// `ln(1 + u)`
    Ln1p,
    Atan,
    Atanh,
    Asin,
    Asinh,
    Tan,
    Tanh,
    Erf,
    LambertW,
}

impl FnKind {
    /// The `n`-th Maclaurin coefficient of the function, as an exact rational
    /// where possible (`erf` carries a `2/√π` factor and is built as an
    /// expression).
    pub(crate) fn coefficient(self, arena: &mut Arena, n: usize) -> ExprId {
        let r = self.rational_coefficient(n);
        match self {
            FnKind::Erf => {
                if r.is_zero() {
                    return arena.zero;
                }
                // (2/√π) · (−1)^m / (m! (2m+1))
                let re = rat_expr(arena, r);
                let two = arena.int(2);
                let pi = arena.pi;
                let sp = arena.sqrt(pi);
                let f = arena.div(two, sp);
                let v = arena.mul(&[re, f]);
                eval::eval(arena, v)
            }
            _ => rat_expr(arena, r),
        }
    }

    /// Rational part of the `n`-th coefficient.
    pub(crate) fn rational_coefficient(self, n: usize) -> Q {
        let odd = n % 2 == 1;
        let m = n / 2;
        let sign_m = if m.is_multiple_of(2) {
            Q::one()
        } else {
            -Q::one()
        };
        match self {
            FnKind::Exp => Q::new(BigInt::one(), factorial(n as u64)),
            FnKind::Sin => {
                if odd {
                    sign_m / Q::from_integer(factorial(n as u64))
                } else {
                    Q::zero()
                }
            }
            FnKind::Cos => {
                if odd {
                    Q::zero()
                } else {
                    sign_m / Q::from_integer(factorial(n as u64))
                }
            }
            FnKind::Sinh => {
                if odd {
                    Q::new(BigInt::one(), factorial(n as u64))
                } else {
                    Q::zero()
                }
            }
            FnKind::Cosh => {
                if odd {
                    Q::zero()
                } else {
                    Q::new(BigInt::one(), factorial(n as u64))
                }
            }
            FnKind::Ln1p => {
                if n == 0 {
                    Q::zero()
                } else {
                    let s = if n % 2 == 1 { Q::one() } else { -Q::one() };
                    s / rat_i(n as i64)
                }
            }
            FnKind::Atan => {
                if odd {
                    sign_m / rat_i(n as i64)
                } else {
                    Q::zero()
                }
            }
            FnKind::Atanh => {
                if odd {
                    Q::one() / rat_i(n as i64)
                } else {
                    Q::zero()
                }
            }
            FnKind::Asin | FnKind::Asinh => {
                if !odd {
                    return Q::zero();
                }
                // C(2m,m) / (4^m (2m+1))
                let c = Q::from_integer(central_binomial(m as u64));
                let d = pow_rat_i(&rat_i(4), m as i64) * rat_i(n as i64);
                let v = c / d;
                if self == FnKind::Asinh { sign_m * v } else { v }
            }
            FnKind::Tan | FnKind::Tanh => {
                // tan x = Σ_{m≥1} (−1)^{m−1} 2^{2m}(2^{2m}−1) B_{2m} x^{2m−1}/(2m)!
                // tanh x = Σ_{m≥1} 2^{2m}(2^{2m}−1) B_{2m} x^{2m−1}/(2m)!
                if !odd {
                    return Q::zero();
                }
                let mm = m + 1; // n = 2mm − 1
                let two_pow = pow_rat_i(&rat_i(2), 2 * mm as i64);
                let b = bernoulli(2 * mm);
                let v = &two_pow * (&two_pow - Q::one()) * b
                    / Q::from_integer(factorial(2 * mm as u64));
                if self == FnKind::Tan {
                    if (mm - 1).is_multiple_of(2) { v } else { -v }
                } else {
                    v
                }
            }
            FnKind::Erf => {
                if !odd {
                    return Q::zero();
                }
                sign_m / (Q::from_integer(factorial(m as u64)) * rat_i(n as i64))
            }
            FnKind::LambertW => {
                if n == 0 {
                    return Q::zero();
                }
                // (−n)^{n−1} / n!
                let base = rat_i(-(n as i64));
                pow_rat_i(&base, n as i64 - 1) / Q::from_integer(factorial(n as u64))
            }
        }
    }
}

/// Generalised binomial coefficient `C(α, n)` for a symbolic or rational `α`.
fn gen_binomial_expr(arena: &mut Arena, alpha: ExprId, n: usize) -> ExprId {
    if n == 0 {
        return arena.one;
    }
    if let Some(a) = arena.as_num(alpha).cloned() {
        let mut acc = Q::one();
        for i in 0..n {
            acc = acc * (&a - rat_i(i as i64)) / rat_i(i as i64 + 1);
        }
        return rat_expr(arena, acc);
    }
    let mut factors = Vec::with_capacity(n + 1);
    for i in 0..n {
        let ie = arena.int(-(i as i64));
        factors.push(arena.add(&[alpha, ie]));
    }
    let inv_fact = rat_expr(arena, Q::new(BigInt::one(), factorial(n as u64)));
    factors.push(inv_fact);
    let p = arena.mul(&factors);
    eval::eval(arena, p)
}

// ═══════════════════════════════════════════════════════════════════════════
// The engine
// ═══════════════════════════════════════════════════════════════════════════

/// How the variable approaches the expansion point.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Side {
    /// Two-sided (ordinary Maclaurin / Laurent expansion).
    Both,
    /// From above only (`x → 0⁺`; used for `x → ±∞` via `t = 1/x`).
    Above,
    /// From below only (`x → 0⁻`; used to cross-check `|g|` with odd
    /// valuation, see [`expand_with_mode`]).
    Below,
}

/// How one expansion is carried out.
#[derive(Clone, Copy, Debug)]
struct Mode {
    side: Side,
    /// With `Some(l)` (only for [`Side::Above`]) the expansion is
    /// *log-extended*: `l` is a `var`-free stand-in for `ln(var)`, and `ln`
    /// of a series with a non-zero valuation `v` is `v·l + ln c + ln(1 + …)`
    /// instead of a refusal.  Coefficients may then be polynomials in `l`;
    /// a term `c(l)·varᵏ` is still "of order `k`" (it is `o(var^(k−ε))`), as
    /// in Gruntz's algorithm, where `l = ln ω` belongs to a lower
    /// comparability class.  `exp`, `sinh`, `cosh` of an argument whose
    /// constant term contains `l` (`exp(2 ln t) = t²`) are refused: that
    /// coefficient would not be of lower order.
    log_var: Option<ExprId>,
    /// May the differentiation fallback take limits of derivatives that
    /// cannot be substituted?  Off inside Gruntz (the limit engine calling
    /// itself through the series fallback has no overall budget).
    limits: bool,
    /// Accept a result exact to fewer terms than requested (its `known`
    /// says how many): Gruntz only needs the leading term, and `ζ(1 + w)`
    /// is known only to `O(w)`.
    partial: bool,
}

impl Mode {
    fn with_side(self, side: Side) -> Self {
        Mode { side, ..self }
    }
}

/// Why a structural expansion step produced no series.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Obstruction {
    /// No structural rule for this node; the caller may fall back to Taylor
    /// coefficients by differentiation.
    Unknown,
    /// No Laurent expansion exists here (Puiseux exponent, `|x|`-type
    /// singularity, non-real `|·|` argument).  The differentiation fallback
    /// must **not** be tried: the low-order derivatives of `x^(5/2)` or
    /// `|x²|` (`2x·sign(x²)`) all vanish at `0`, so it would silently return
    /// the wrong polynomial `0`.
    NoExpansion,
    /// `|g|` with `g` of odd valuation in a two-sided expansion: the two
    /// one-sided expansions of `|g|` differ, but the enclosing expression may
    /// still have a two-sided expansion (`cos|x|`).
    Kink,
}

/// [`expand_with_mode`] for a function that must be *meromorphic* at `0`
/// (a Laurent series in the complex sense, for [`laurent_series`] and
/// formal power series): two-sided, and the differentiation fallback takes
/// no limits, so an essential singularity (`exp(−1/x²)`, whose real
/// Taylor coefficients all vanish) is refused.
pub(crate) fn expand_meromorphic(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    order: i64,
) -> Result<TSeries, SymplexError> {
    expand_with_mode(
        arena,
        expr,
        var,
        order,
        Mode {
            side: Side::Both,
            log_var: None,
            limits: false,
            partial: false,
        },
    )
}

/// Log-extended expansion as `var → 0⁺` (see [`Mode::log_var`]): `log_var`
/// stands for `ln(var)` and may appear (polynomially) in the coefficients.
/// `limits` allows the differentiation fallback to take limits.
///
/// Used by [`series_at_infinity`] (`var = 1/x`, `log_var` a placeholder).
pub(crate) fn expand_log_extended(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    order: i64,
    log_var: ExprId,
    limits: bool,
) -> Result<TSeries, SymplexError> {
    expand_with_mode(
        arena,
        expr,
        var,
        order,
        Mode {
            side: Side::Above,
            log_var: Some(log_var),
            limits,
            partial: false,
        },
    )
}

/// [`expand_log_extended`] for a leading term: no limits in the
/// differentiation fallback, and a result exact to fewer than `order` terms
/// is returned rather than refused (its [`TSeries::known`] tells how far
/// it is exact).  Used by Gruntz's leading-term extraction (`var = ω`,
/// `log_var = ln ω`).
pub(crate) fn expand_leading(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    order: i64,
    log_var: ExprId,
) -> Result<TSeries, SymplexError> {
    expand_with_mode(
        arena,
        expr,
        var,
        order,
        Mode {
            side: Side::Above,
            log_var: Some(log_var),
            limits: false,
            partial: true,
        },
    )
}

/// Expand `expr` around `var = 0` with all exponents `< order` exact, in the
/// given [`Mode`].
///
/// From above ([`Side::Above`], used for `x → ±∞` via `t = 1/x`)
/// `(t^v)^α = t^{vα}` holds for every integer `vα`, which is not valid
/// two-sided (`√(x²) = |x|`).
///
/// A two-sided expansion that fails only because of `|g|` with odd
/// valuation (`|x|`, `|sin x|`) is retried from each side separately and
/// accepted when both sides agree (`cos|x| = cos x`); otherwise there is
/// no two-sided expansion and the error is returned.
fn expand_with_mode(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    order: i64,
    mode: Mode,
) -> Result<TSeries, SymplexError> {
    let mut kink = false;
    match expand_from_side(arena, expr, var, order, mode, &mut kink) {
        Err(_) if kink => {
            let above = expand_from_side(
                arena,
                expr,
                var,
                order,
                mode.with_side(Side::Above),
                &mut false,
            )?;
            let below = expand_from_side(
                arena,
                expr,
                var,
                order,
                mode.with_side(Side::Below),
                &mut false,
            )?;
            if TSeries::same(arena, &above, &below, order) {
                Ok(above)
            } else {
                Err(SymplexError::ComputationFailed {
                    operation: "series",
                    reason: format!(
                        "no two-sided expansion: the expansions of {} from the left and from the right differ (|·| with odd valuation)",
                        arena.display(expr)
                    ),
                })
            }
        }
        r => r,
    }
}

/// [`expand_with_mode`] for one [`Side`], with the precision-escalation loop.
/// Sets `kink` when the failure was an odd-valuation `|g|` in a two-sided
/// expansion.
fn expand_from_side(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    order: i64,
    mode: Mode,
    kink: &mut bool,
) -> Result<TSeries, SymplexError> {
    // Work with at least a few terms so that the valuation of every
    // sub-expression is visible (at precision 1 the variable itself would
    // truncate to nothing and `1/x` could not be expanded).
    let mut working = order.max(MIN_WORKING_ORDER);
    let mut last_err = None;
    let mut best: Option<TSeries> = None;
    for _attempt in 0..MAX_PRECISION_ATTEMPTS {
        let mut hidden_valuation = false;
        match expand_with_precision(arena, expr, var, working, mode, &mut hidden_valuation, kink) {
            Ok(ts) if ts.known >= order => return Ok(ts.truncate_known(order)),
            Ok(ts) => {
                // Precision was lost through poles; increase and retry.
                working += order - ts.known + 1;
                if mode.partial && best.as_ref().is_none_or(|b| b.known < ts.known) {
                    best = Some(ts);
                }
            }
            Err(e) if hidden_valuation && !*kink => {
                // A sub-expression's leading term lay beyond the working
                // window (e.g. `1/(x⁵ + x⁶)` at low order): widen and retry.
                last_err = Some(e);
                working = working * 2 + 4;
            }
            Err(e) => return Err(e),
        }
    }
    if let Some(ts) = best {
        return Ok(ts);
    }
    Err(last_err.unwrap_or(SymplexError::ComputationFailed {
        operation: "series",
        reason: "could not reach the requested order (deep pole)".into(),
    }))
}

/// One pass of the engine at working precision `n`.  Sets `hidden_valuation`
/// when some `var`-dependent sub-expression had *no* visible term at this
/// precision, so a failure may be curable by widening the window; sets
/// `kink` when the failure was [`Obstruction::Kink`].
fn expand_with_precision(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    n: i64,
    mode: Mode,
    hidden_valuation: &mut bool,
    kink: &mut bool,
) -> Result<TSeries, SymplexError> {
    let order_ids = walk::post_order_ids(arena, expr);
    let mut cache: FxHashMap<ExprId, Option<TSeries>> = FxHashMap::default();
    let mut fallbacks_used = 0usize;
    for id in order_ids {
        if cache.contains_key(&id) {
            continue;
        }
        let ts = if !walk::contains(arena, id, var) {
            Some(TSeries::constant(arena, id, n))
        } else if id == var {
            Some(TSeries::var(arena, n))
        } else {
            match structural_series(arena, id, var, mode, &cache) {
                Ok(s) => Some(s),
                Err(Obstruction::Unknown) => {
                    if id == expr || fallbacks_used < MAX_FALLBACK_NODES {
                        fallbacks_used += 1;
                        taylor_by_differentiation(arena, id, var, n, mode)
                    } else {
                        None
                    }
                }
                Err(Obstruction::Kink) => {
                    *kink = true;
                    return Err(no_expansion_error(arena, id));
                }
                Err(Obstruction::NoExpansion) => return Err(no_expansion_error(arena, id)),
            }
        };
        if let Some(s) = &ts
            && id != expr
            && s.clone().normalized(arena).coeffs.is_empty()
        {
            *hidden_valuation = true;
        }
        cache.insert(id, ts);
    }
    match cache.remove(&expr).flatten() {
        Some(ts) => Ok(ts),
        None => Err(SymplexError::ComputationFailed {
            operation: "series",
            reason: "no Laurent expansion at this point (singularity or unsupported function)"
                .into(),
        }),
    }
}

fn no_expansion_error(arena: &Arena, id: ExprId) -> SymplexError {
    SymplexError::ComputationFailed {
        operation: "series",
        reason: format!(
            "no Laurent expansion of {} at this point (fractional power, |·| or non-real argument)",
            arena.display(id)
        ),
    }
}

/// The cached series of a child node, or [`Obstruction::Unknown`] when the
/// child itself could not be expanded.
fn child(cache: &FxHashMap<ExprId, Option<TSeries>>, id: ExprId) -> Result<TSeries, Obstruction> {
    cache
        .get(&id)
        .cloned()
        .flatten()
        .ok_or(Obstruction::Unknown)
}

/// Combine children's series according to the node type.
fn structural_series(
    arena: &mut Arena,
    id: ExprId,
    var: ExprId,
    mode: Mode,
    cache: &FxHashMap<ExprId, Option<TSeries>>,
) -> Result<TSeries, Obstruction> {
    let side = mode.side;
    let node = arena.node(id).clone();
    let analytic = match node {
        ExprNode::Add(ref ch) => {
            let mut acc: Option<TSeries> = None;
            for &c in ch.iter() {
                let s = child(cache, c)?;
                acc = Some(match acc {
                    None => s,
                    Some(a) => TSeries::add(arena, &a, &s),
                });
            }
            acc
        }
        ExprNode::Mul(ref ch) => {
            let mut acc: Option<TSeries> = None;
            for &c in ch.iter() {
                let s = child(cache, c)?;
                acc = Some(match acc {
                    None => s,
                    Some(a) => TSeries::mul(arena, &a, &s),
                });
            }
            acc
        }
        ExprNode::Neg(inner) => {
            let s = child(cache, inner)?;
            Some(TSeries::neg(arena, &s))
        }
        ExprNode::Pow(base, exp) => {
            let b = child(cache, base)?;
            if !walk::contains(arena, exp, var) {
                let e = arena.as_num(exp).cloned().ok_or(Obstruction::Unknown)?;
                if e.is_integer() {
                    let ei = e.to_integer().to_i64().ok_or(Obstruction::Unknown)?;
                    if ei.abs() <= MAX_INT_POWER {
                        return TSeries::pow_int(arena, &b, ei).ok_or(Obstruction::Unknown);
                    }
                    // Large integer exponents: binomial series with
                    // closed-form coefficients instead of repeated products.
                }
                // Rational exponent: Puiseux expansions are refused.
                return pow_rational(arena, &b, exp, side);
            }
            // b^e with var-dependent exponent: exp(e · ln b).
            let e = child(cache, exp)?;
            let lnb = apply_ln(arena, &b, mode)?;
            let prod = TSeries::mul(arena, &e, &lnb);
            apply_exp_like(arena, FnKind::Exp, &prod, mode)
        }
        ExprNode::Abs(a) => return abs_series(arena, &child(cache, a)?, side),
        ExprNode::Sign(a) => return sign_series(arena, &child(cache, a)?, side, false),
        ExprNode::Heaviside(a) => return sign_series(arena, &child(cache, a)?, side, true),
        ExprNode::Apply(f, ref args) if args.len() == 2 && is_bessel(arena, f) => {
            let u = child(cache, args[1])?;
            if let Some(s) = log_bessel_series(arena, f, args[0], &u, var, mode)? {
                return Ok(s);
            }
            return bessel_series(arena, f, args[0], &u, var);
        }
        ExprNode::Gamma(_)
        | ExprNode::LogGamma(_)
        | ExprNode::Digamma(_)
        | ExprNode::Polygamma(..)
        | ExprNode::Zeta(_)
        | ExprNode::Ei(_)
        | ExprNode::Ci(_)
        | ExprNode::Si(_)
        | ExprNode::Li(_)
        | ExprNode::Erfc(_) => return special_series(arena, &node, mode, cache),
        ExprNode::Apply(f, ref args)
            if args.len() == 1
                && (matches!(
                    arena.lib_fn(f),
                    Some(crate::base::libfn::LibFn::Shi | crate::base::libfn::LibFn::Chi)
                ) || matches!(arena.symbol_name(f), ERFCX_ASYMPTOTIC | EI_ASYMPTOTIC)) =>
        {
            return special_series(arena, &node, mode, cache);
        }
        ExprNode::Exp(a) => apply_exp_like(arena, FnKind::Exp, &child(cache, a)?, mode),
        ExprNode::Sin(a) => apply_fn(arena, FnKind::Sin, &child(cache, a)?),
        ExprNode::Cos(a) => apply_fn(arena, FnKind::Cos, &child(cache, a)?),
        ExprNode::Sinh(a) => apply_exp_like(arena, FnKind::Sinh, &child(cache, a)?, mode),
        ExprNode::Cosh(a) => apply_exp_like(arena, FnKind::Cosh, &child(cache, a)?, mode),
        ExprNode::Tan(a) => apply_fn(arena, FnKind::Tan, &child(cache, a)?),
        ExprNode::Tanh(a) => apply_fn(arena, FnKind::Tanh, &child(cache, a)?),
        ExprNode::Atan(a) => apply_fn(arena, FnKind::Atan, &child(cache, a)?),
        ExprNode::Atanh(a) => apply_fn(arena, FnKind::Atanh, &child(cache, a)?),
        ExprNode::Asin(a) => apply_fn(arena, FnKind::Asin, &child(cache, a)?),
        ExprNode::Asinh(a) => apply_fn(arena, FnKind::Asinh, &child(cache, a)?),
        ExprNode::Erf(a) => apply_fn(arena, FnKind::Erf, &child(cache, a)?),
        ExprNode::LambertW(a) => apply_fn(arena, FnKind::LambertW, &child(cache, a)?),
        ExprNode::Ln(a) => return apply_ln(arena, &child(cache, a)?, mode),
        _ => None,
    };
    analytic.ok_or(Obstruction::Unknown)
}

/// `exp`, `sinh`, `cosh` of `a`, refused in a log-extended expansion when
/// the constant term of `a` contains the stand-in for `ln(var)` (see
/// [`Mode::log_var`]).
fn apply_exp_like(arena: &mut Arena, kind: FnKind, a: &TSeries, mode: Mode) -> Option<TSeries> {
    if let Some(l) = mode.log_var {
        let u0 = a.coeff_at(arena, 0);
        if walk::contains(arena, u0, l) {
            return None;
        }
    }
    apply_fn(arena, kind, a)
}

fn is_bessel(arena: &Arena, f: crate::base::node::SymbolId) -> bool {
    use crate::base::libfn::LibFn;
    matches!(
        arena.lib_fn(f),
        Some(LibFn::BesselJ | LibFn::BesselY | LibFn::BesselI | LibFn::BesselK)
    )
}

/// A Bessel function of the series `u` at a zero of `u`: `J_n(u)` and
/// `I_n(u)` of integer order from `J_n(z) = Σ_k (−1)^k (z/2)^{2k+n}/(k!(k+n)!)`
/// (`I_n` without the signs; `J_{−n} = (−1)ⁿ J_n`, `I_{−n} = I_n`); `Y_ν`,
/// `K_ν` (logarithmic singularity) and `J_ν`, `I_ν` of any other order (a
/// fractional power, or unknown) have no Laurent expansion there.  Before,
/// the differentiation fallback returned the coefficients of `Y₀` in terms
/// of `Y₀(0) = −∞` as if they were finite.  An argument with a pole is an
/// essential singularity; elsewhere the fallback (finite values) applies.
fn bessel_series(
    arena: &mut Arena,
    f: crate::base::node::SymbolId,
    order: ExprId,
    u: &TSeries,
    var: ExprId,
) -> Result<TSeries, Obstruction> {
    use crate::base::libfn::LibFn;
    let u = u.clone().normalized(arena);
    let Some(k) = u.leading_exponent(arena) else {
        return Err(Obstruction::NoExpansion);
    };
    if k < 0 {
        return Err(Obstruction::NoExpansion);
    }
    if k == 0 {
        return Err(Obstruction::Unknown);
    }
    let lib = arena.lib_fn(f).ok_or(Obstruction::Unknown)?;
    if walk::contains(arena, order, var) || !matches!(lib, LibFn::BesselJ | LibFn::BesselI) {
        return Err(Obstruction::NoExpansion);
    }
    let n = arena
        .as_num(order)
        .filter(|r| r.is_integer())
        .and_then(|r| r.to_integer().to_i64())
        .ok_or(Obstruction::NoExpansion)?;
    let alternate = lib == LibFn::BesselJ;
    let reflect_sign = alternate && n < 0 && n % 2 != 0;
    let n = n.unsigned_abs();
    let coefficient = move |ar: &mut Arena, m: usize| -> ExprId {
        let m = m as u64;
        if m < n || (m - n) % 2 == 1 {
            return ar.zero;
        }
        let kk = (m - n) / 2;
        let denom = factorial(kk) * factorial(kk + n) * (BigInt::one() << (2 * kk + n));
        let mut c = Q::new(BigInt::one(), denom);
        if alternate && kk % 2 == 1 {
            c = -c;
        }
        if reflect_sign {
            c = -c;
        }
        rat_expr(ar, c)
    };
    Ok(TSeries::compose(arena, &coefficient, &u))
}

/// `|g|` for a series `g` with leading term `c·x^k`, `c` a real constant of
/// known sign and all visible coefficients real.
///
/// Near `0` the sign of `g` is that of `c·x^k`, so `|g| = ±g`:
///
/// * `k` even (including `k = 0`): `sign(c)·g` on both sides;
/// * `k` odd, one-sided: `sign(c)·g` from above, `−sign(c)·g` from below;
/// * `k` odd, two-sided: [`Obstruction::Kink`] — no expansion of `|g|`
///   itself, but the caller re-expands the whole expression from each side.
///
/// A constant term of unknown sign (`|a + x|`) is left to the
/// differentiation fallback (`|a| + sign(a)·x`, correct for real `a ≠ 0`);
/// an unknown sign with `k ≠ 0`, or a non-real coefficient, is a definite
/// [`Obstruction::NoExpansion`] — differentiating `|g|` at a zero of `g`
/// gives `sign(0) = 0` and hence a wrong all-zero polynomial.
fn abs_series(arena: &mut Arena, g: &TSeries, side: Side) -> Result<TSeries, Obstruction> {
    let g = g.clone().normalized(arena);
    // No visible term: the valuation is hidden at this precision.  Fail
    // definitively so that the caller widens the window instead of
    // differentiating.
    let k = g.leading_exponent(arena).ok_or(Obstruction::NoExpansion)?;
    let c = g.coeff_at(arena, k);
    let Some(c_positive) = constant_sign(arena, c) else {
        return Err(if k == 0 {
            Obstruction::Unknown
        } else {
            Obstruction::NoExpansion
        });
    };
    if !g.coeffs.iter().all(|&co| is_known_real(arena, co)) {
        return Err(Obstruction::NoExpansion);
    }
    let flip_below = k.rem_euclid(2) == 1;
    let positive = match side {
        Side::Both if flip_below => return Err(Obstruction::Kink),
        Side::Below if flip_below => !c_positive,
        Side::Both | Side::Above | Side::Below => c_positive,
    };
    Ok(if positive { g } else { TSeries::neg(arena, &g) })
}

/// `sign(g)` (or `H(g)` when `heaviside`) for a series `g` whose leading
/// term `c·x^k` has a real constant `c` of known sign and whose visible
/// coefficients are real: near `0` the sign of `g` is that of `c·x^k` (see
/// [`abs_series`] for the sides), a constant `±1` (`1` or `0`) exact to any
/// order.
///
/// Anything else is a definite [`Obstruction::NoExpansion`]: the
/// differentiation fallback evaluates `sign(g(0)) = sign(0) = 0` and
/// `d sign(g) = 0`, the wrong polynomial `0` — `atanh(x²)·sign(x²)` expanded
/// to `0` instead of `x²`.
fn sign_series(
    arena: &mut Arena,
    g: &TSeries,
    side: Side,
    heaviside: bool,
) -> Result<TSeries, Obstruction> {
    let g = g.clone().normalized(arena);
    let k = g.leading_exponent(arena).ok_or(Obstruction::NoExpansion)?;
    let c = g.coeff_at(arena, k);
    // A non-zero constant term of unknown sign: `sign(a + x)` is
    // `sign(a)` near 0, which the differentiation fallback finds.
    let Some(c_positive) = constant_sign(arena, c) else {
        return Err(if k == 0 {
            Obstruction::Unknown
        } else {
            Obstruction::NoExpansion
        });
    };
    if !g.coeffs.iter().all(|&co| is_known_real(arena, co)) {
        return Err(Obstruction::NoExpansion);
    }
    let flip_below = k.rem_euclid(2) == 1;
    let positive = match side {
        Side::Both if flip_below => return Err(Obstruction::Kink),
        Side::Below if flip_below => !c_positive,
        Side::Both | Side::Above | Side::Below => c_positive,
    };
    let value = match (positive, heaviside) {
        (true, _) => arena.one,
        (false, false) => arena.neg_one,
        (false, true) => arena.zero,
    };
    Ok(TSeries::constant(arena, value, g.known.max(1)))
}

/// `Some(true)` if the var-free constant `c` is known positive, `Some(false)`
/// if known negative, `None` if zero or of unknown sign.
fn constant_sign(arena: &mut Arena, c: ExprId) -> Option<bool> {
    if let Some(r) = arena.as_num(c) {
        return if r.is_zero() {
            None
        } else {
            Some(r.is_positive())
        };
    }
    let mut cache = crate::base::assumptions::AssumptionCache::new();
    if cache.query(arena, c, Props::POSITIVE) == Some(true) {
        return Some(true);
    }
    if cache.query(arena, c, Props::NEGATIVE) == Some(true) {
        return Some(false);
    }
    // A real constant the assumption system leaves open (`ln 3`,
    // `cos(3/2)`): its certified value decides.
    match crate::calculus::limit::const_sign(arena, c) {
        Some(1) => Some(true),
        Some(-1) => Some(false),
        _ => None,
    }
}

/// Is the var-free constant `c` known to be real?
fn is_known_real(arena: &Arena, c: ExprId) -> bool {
    arena.as_num(c).is_some() || {
        let mut cache = crate::base::assumptions::AssumptionCache::new();
        cache.query(arena, c, Props::REAL) == Some(true)
    }
}

fn is_zero_const(arena: &mut Arena, c: ExprId) -> bool {
    let v = settle_constant(arena, c);
    arena.is_zero_structural(v) || arena.as_num(v).is_some_and(|r| r.is_zero())
}

/// `f(a)` for an elementary `f` with known Maclaurin series.
fn apply_fn(arena: &mut Arena, kind: FnKind, a: &TSeries) -> Option<TSeries> {
    let a = a.clone().normalized(arena);
    if a.shift < 0 && a.leading_exponent(arena).is_some_and(|v| v < 0) {
        return None; // essential singularity
    }
    let (u0, w) = a.split_constant(arena);
    let u0_zero = is_zero_const(arena, u0);
    let compose = |arena: &mut Arena, k: FnKind, w: &TSeries| -> TSeries {
        TSeries::compose(
            arena,
            &move |ar: &mut Arena, i: usize| k.coefficient(ar, i),
            w,
        )
    };
    match kind {
        FnKind::Exp => {
            let s = compose(arena, FnKind::Exp, &w);
            if u0_zero {
                Some(s)
            } else {
                let e = arena.exp(u0);
                let e = eval::eval(arena, e);
                Some(TSeries::scale(arena, &s, e))
            }
        }
        FnKind::Sin | FnKind::Cos => {
            let sw = compose(arena, FnKind::Sin, &w);
            let cw = compose(arena, FnKind::Cos, &w);
            if u0_zero {
                return Some(if kind == FnKind::Sin { sw } else { cw });
            }
            let su = arena.sin(u0);
            let su = eval::eval(arena, su);
            let cu = arena.cos(u0);
            let cu = eval::eval(arena, cu);
            let (t1, t2) = if kind == FnKind::Sin {
                // sin(u0 + w) = sin u0 cos w + cos u0 sin w
                (
                    TSeries::scale(arena, &cw, su),
                    TSeries::scale(arena, &sw, cu),
                )
            } else {
                // cos(u0 + w) = cos u0 cos w − sin u0 sin w
                let nsu = arena.neg(su);
                (
                    TSeries::scale(arena, &cw, cu),
                    TSeries::scale(arena, &sw, nsu),
                )
            };
            Some(TSeries::add(arena, &t1, &t2))
        }
        FnKind::Sinh | FnKind::Cosh => {
            let sw = compose(arena, FnKind::Sinh, &w);
            let cw = compose(arena, FnKind::Cosh, &w);
            if u0_zero {
                return Some(if kind == FnKind::Sinh { sw } else { cw });
            }
            let su = arena.sinh(u0);
            let su = eval::eval(arena, su);
            let cu = arena.cosh(u0);
            let cu = eval::eval(arena, cu);
            let (t1, t2) = if kind == FnKind::Sinh {
                (
                    TSeries::scale(arena, &cw, su),
                    TSeries::scale(arena, &sw, cu),
                )
            } else {
                (
                    TSeries::scale(arena, &cw, cu),
                    TSeries::scale(arena, &sw, su),
                )
            };
            Some(TSeries::add(arena, &t1, &t2))
        }
        // `ln` is dispatched to `apply_ln` (it needs the constant term split
        // differently); a direct `Ln1p` request has no meaning here.
        FnKind::Ln1p => None,
        _ => {
            if !u0_zero {
                return None; // fallback: differentiate
            }
            Some(compose(arena, kind, &w))
        }
    }
}

/// `ln(a)` — requires a non-zero constant term, except in a log-extended
/// expansion ([`Mode::log_var`]), where `a = c·var^v·(1 + …)` with `c` real
/// of known sign and real visible coefficients gives
/// `ln c + v·ln(var) + ln(1 + …)` (for `var > 0` the argument of `var^v`
/// is 0, so no multiple of `2πi` is lost).
///
/// A constant term on the cut (negative real) with a non-real correction is
/// a definite [`Obstruction::NoExpansion`]: `ln(−1 + i·x)` tends to `iπ`
/// from `x > 0` and to `−iπ` from `x < 0` (before, the series was
/// `iπ − i·x + …` on both sides).
fn apply_ln(arena: &mut Arena, a: &TSeries, mode: Mode) -> Result<TSeries, Obstruction> {
    let mut a = a.clone().normalized(arena);
    let v = a.leading_exponent(arena).ok_or(Obstruction::Unknown)?;
    let mut log_shift = None;
    if v != 0 {
        // logarithmic singularity
        let l = mode.log_var.ok_or(Obstruction::Unknown)?;
        if mode.side != Side::Above {
            return Err(Obstruction::Unknown);
        }
        let c = a.coeff_at(arena, v);
        match constant_sign(arena, c) {
            Some(true) => {}
            Some(false) if a.coeffs.iter().all(|&co| is_known_real(arena, co)) => {}
            _ => return Err(Obstruction::Unknown),
        }
        a.shift -= v;
        a.known -= v;
        let vl = {
            let ve = arena.int(v);
            let p = arena.mul(&[ve, l]);
            eval::eval(arena, p)
        };
        log_shift = Some(vl);
    }
    let (u0, w) = a.split_constant(arena);
    if approaches_cut(arena, u0, &w) {
        return Err(Obstruction::NoExpansion);
    }
    // ln(u0 + w) = ln u0 + ln(1 + w/u0)
    let inv_u0 = {
        let m1 = arena.neg_one;
        let p = arena.pow(u0, m1);
        eval::eval(arena, p)
    };
    let w_over = TSeries::scale(arena, &w, inv_u0);
    let mut s = TSeries::compose(
        arena,
        &|ar: &mut Arena, i: usize| FnKind::Ln1p.coefficient(ar, i),
        &w_over,
    );
    let mut ln_u0 = ln_of_constant(arena, u0);
    if let Some(vl) = log_shift {
        let sum = arena.add(&[ln_u0, vl]);
        ln_u0 = eval::eval(arena, sum);
    }
    if !arena.is_zero_structural(ln_u0) && !s.coeffs.is_empty() && s.shift <= 0 {
        let idx = (-s.shift) as usize;
        let c = arena.add(&[s.coeffs[idx], ln_u0]);
        s.coeffs[idx] = eval::eval(arena, c);
    }
    Ok(s)
}

/// `ln c` of a constant, with `ln(p/q) = ln p − ln q` for a positive
/// rational, so that the constant cancels structurally against one spelled
/// with integers (the `ln 2` of K₀'s expansion against
/// `ln(x/2) = ln x + ln(1/2)`).
pub(crate) fn ln_of_constant(arena: &mut Arena, c: ExprId) -> ExprId {
    let l = match arena.as_num(c).cloned() {
        Some(r) if r.is_positive() && !r.is_integer() => {
            let p = rat_expr(arena, Q::from_integer(r.numer().clone()));
            let q = rat_expr(arena, Q::from_integer(r.denom().clone()));
            let lp = arena.ln(p);
            let lq = arena.ln(q);
            arena.sub(lp, lq)
        }
        _ => arena.ln(c),
    };
    eval::eval(arena, l)
}

/// Does `u0 + w` approach the negative real axis (the cut of `ln` and of
/// fractional powers) with a correction `w` that is not known to be real?
fn approaches_cut(arena: &mut Arena, u0: ExprId, w: &TSeries) -> bool {
    is_known_real(arena, u0)
        && constant_sign(arena, u0) == Some(false)
        && !w.coeffs.iter().all(|&c| is_known_real(arena, c))
}

/// `a^α` for a var-free non-integer exponent: `u0^α · Σ C(α,n) (w/u0)^n`.
///
/// A non-zero valuation `v` is accepted only when `vα` is an integer and the
/// result is single-valued: always for expansions from above, otherwise only
/// when `v / denom(α)` is even (so that no `|x|` appears).  A genuine
/// Puiseux series (`x^(5/2)`) or an `|x|` is a definite
/// [`Obstruction::NoExpansion`]: the differentiation fallback would see only
/// vanishing derivatives at low orders and return `0`.
fn pow_rational(
    arena: &mut Arena,
    a: &TSeries,
    alpha: ExprId,
    side: Side,
) -> Result<TSeries, Obstruction> {
    let mut a = a.clone().normalized(arena);
    // No visible term: the valuation is hidden at this precision; fail
    // definitively so that the caller widens the window.
    let v = a.leading_exponent(arena).ok_or(Obstruction::NoExpansion)?;
    let mut outer_shift = 0i64;
    if v != 0 {
        let ar = arena.as_num(alpha).cloned().ok_or(Obstruction::Unknown)?;
        let va = &ar * rat_i(v);
        if !va.is_integer() {
            return Err(Obstruction::NoExpansion); // genuine Puiseux series
        }
        if side != Side::Above && !ar.is_integer() {
            let q = ar.denom().to_i64().ok_or(Obstruction::Unknown)?;
            if (v / q) % 2 != 0 {
                return Err(Obstruction::NoExpansion); // would introduce |x|
            }
        }
        outer_shift = va.to_integer().to_i64().ok_or(Obstruction::Unknown)?;
        a.shift -= v;
        a.known -= v;
    }
    let (u0, w) = a.split_constant(arena);
    if arena.as_num(alpha).is_none_or(|r| !r.is_integer()) && approaches_cut(arena, u0, &w) {
        return Err(Obstruction::NoExpansion);
    }
    let inv_u0 = {
        let m1 = arena.neg_one;
        let p = arena.pow(u0, m1);
        eval::eval(arena, p)
    };
    let w_over = TSeries::scale(arena, &w, inv_u0);
    let s = TSeries::compose(
        arena,
        &move |ar: &mut Arena, i: usize| gen_binomial_expr(ar, alpha, i),
        &w_over,
    );
    let u0a = arena.pow(u0, alpha);
    let u0a = eval::eval(arena, u0a);
    let mut r = TSeries::scale(arena, &s, u0a);
    r.shift += outer_shift;
    r.known += outer_shift;
    Ok(r)
}

// ═══════════════════════════════════════════════════════════════════════════
// Special functions at their singular points
// ═══════════════════════════════════════════════════════════════════════════
//
// Where the argument of `Γ`, `ψ`, `ψ⁽ᵐ⁾`, `ζ`, `ln Γ`, `Ei`, `Ci`, `Chi`,
// `li`, `K₀`, `Y₀` tends to a pole, a logarithmic singularity or `+∞`, the
// differentiation fallback has nothing to evaluate; these rules give the
// expansion there in closed form (DLMF 5.7.1, 5.7.4, 5.11.1, 5.11.2, 5.15.1,
// 6.6.2, 6.6.5, 6.6.6, 6.2.8, 10.31.2, 10.8.2, 25.2.4; the recurrences
// Γ(z+1) = zΓ(z), ψ(z+1) = ψ(z) + 1/z).  Elsewhere they answer
// [`Obstruction::Unknown`] and the fallback applies as before.  Logarithmic
// singularities need a log-extended expansion ([`Mode::log_var`]); without
// one they are a definite [`Obstruction::NoExpansion`].

/// Most poles `Γ(−n)`, `ψ(−n)` unwound by the recurrence.
const MAX_POLE_SHIFT: i64 = 64;

/// Head of Gruntz's internal `erfcx(z) = e^{z²} erfc(z)` for `z → +∞`
/// (SymPy's `_erfs`), with the asymptotic series
/// `(1/√π) Σ (−1)ᵏ (2k−1)!!/(2ᵏ z^{2k+1})` (DLMF 7.12.1).  It exists only
/// inside a limit computation.
pub(crate) const ERFCX_ASYMPTOTIC: &str = "__erfcx_asym";

/// Head of Gruntz's internal `e^{−z} Ei(z)` for `z → ±∞` (SymPy's
/// `_eis`), with the asymptotic series `Σ k!/z^{k+1}` (DLMF 6.12.2).
pub(crate) const EI_ASYMPTOTIC: &str = "__ei_asym";

/// Most terms of an asymptotic (Stirling-type) series.
const MAX_ASYMPTOTIC_TERMS: i64 = 40;

/// The argument series `a` as `u0 + w` with `w` of valuation `≥ 1`, if `a`
/// has no pole.  `Err(NoExpansion)` when `a` has no visible term (a hidden
/// valuation: the caller widens the window).
fn split_regular(arena: &mut Arena, a: &TSeries) -> Result<Option<(ExprId, TSeries)>, Obstruction> {
    let a = a.clone().normalized(arena);
    if a.coeffs.is_empty() {
        return Err(Obstruction::NoExpansion);
    }
    if a.shift < 0 {
        return Ok(None);
    }
    let (u0, w) = a.split_constant(arena);
    Ok(Some((settle_constant(arena, u0), w)))
}

/// `n` if the constant `u0` is the integer `−n ≤ 0`.
fn nonpositive_integer(arena: &Arena, u0: ExprId) -> Option<i64> {
    if arena.is_zero_structural(u0) {
        return Some(0);
    }
    let r = arena.as_num(u0)?;
    if r.is_integer() && !r.is_positive() {
        r.to_integer().to_i64().map(|n| -n)
    } else {
        None
    }
}

fn is_one_const(arena: &Arena, u0: ExprId) -> bool {
    arena.as_num(u0).is_some_and(|r| r.is_one())
}

/// `ζ(k)` for an integer `k ≥ 2` (`π²/6`, `ζ(3)`, …).
fn zeta_int(arena: &mut Arena, k: usize) -> ExprId {
    let ke = arena.int(k as i64);
    let z = arena.zeta(ke);
    eval::eval(arena, z)
}

/// `c·s` for a rational `c`.
fn scale_rat(arena: &mut Arena, s: &TSeries, c: Q) -> TSeries {
    let ce = rat_expr(arena, c);
    TSeries::scale(arena, s, ce)
}

/// The constant series `c` with the precision of `like`.
fn constant_like(arena: &Arena, c: ExprId, like: &TSeries) -> TSeries {
    TSeries::constant(arena, c, like.known)
}

/// `(w + s)^(−power)` for an integer shift `s`.
fn recip_shifted(arena: &mut Arena, w: &TSeries, s: i64, power: i64) -> Option<TSeries> {
    let se = arena.int(s);
    let c = constant_like(arena, se, w);
    let base = TSeries::add(arena, w, &c);
    TSeries::pow_int(arena, &base, -power)
}

/// `ln Γ(1 + w) = −γw + Σ_{k≥2} (−1)ᵏ ζ(k) wᵏ/k` (DLMF 5.7.3).
fn ln_gamma_one_plus(arena: &mut Arena, w: &TSeries) -> TSeries {
    TSeries::compose(
        arena,
        &|ar: &mut Arena, k: usize| match k {
            0 => ar.zero,
            1 => {
                let g = ar.euler_gamma;
                ar.neg(g)
            }
            _ => {
                let z = zeta_int(ar, k);
                let sign = if k.is_multiple_of(2) { 1 } else { -1 };
                let c = rat_expr(ar, Q::new(BigInt::from(sign), BigInt::from(k)));
                let p = ar.mul(&[c, z]);
                eval::eval(ar, p)
            }
        },
        w,
    )
}

/// `Γ(1 + w)` from the Taylor coefficients `gₙ` of `exp(ln Γ(1 + w))`,
/// `g₀ = 1`, `n·gₙ = Σ_{k=1}^{n} k·c_k·g_{n−k}` with `c_k` those of
/// [`ln_gamma_one_plus`].  The coefficients are built once as closed forms
/// (polynomials in `γ`, `ζ(k)`; none vanishes, `gₙ ≈ (−1)ⁿ`), so the
/// series arithmetic of `exp` — with a numerical zero test of every sum of
/// `ζ` values — is avoided.
fn gamma_one_plus(arena: &mut Arena, w: &TSeries) -> TSeries {
    // `TSeries::compose` asks for coefficients up to `w.known + 1`.
    let n_max = w.known.max(1) as usize + 2;
    let mut c: Vec<ExprId> = Vec::with_capacity(n_max + 1);
    c.push(arena.zero);
    for k in 1..=n_max {
        c.push(if k == 1 {
            let g = arena.euler_gamma;
            arena.neg(g)
        } else {
            let z = zeta_int(arena, k);
            let sign = if k.is_multiple_of(2) { 1 } else { -1 };
            let q = rat_expr(arena, Q::new(BigInt::from(sign), BigInt::from(k)));
            let p = arena.mul(&[q, z]);
            eval::eval(arena, p)
        });
    }
    let mut g: Vec<ExprId> = Vec::with_capacity(n_max + 1);
    g.push(arena.one);
    for n in 1..=n_max {
        let mut terms = Vec::with_capacity(n);
        for k in 1..=n {
            let kq = rat_expr(arena, Q::new(BigInt::from(k), BigInt::from(n)));
            terms.push(arena.mul(&[kq, c[k], g[n - k]]));
        }
        let s = arena.add(&terms);
        let s = crate::transforms::expand::expand(arena, s);
        g.push(eval::eval(arena, s));
    }
    TSeries::compose(
        arena,
        &move |ar: &mut Arena, n: usize| g.get(n).copied().unwrap_or(ar.zero),
        w,
    )
}

/// `ψ⁽ᵐ⁾(1 + w) = Σ_j (−1)^{m+j+1} (m+j)!/j! ζ(m+j+1) wʲ` (`ψ(1) = −γ`),
/// the derivatives of [`ln_gamma_one_plus`].
fn polygamma_one_plus(arena: &mut Arena, m: usize, w: &TSeries) -> TSeries {
    TSeries::compose(
        arena,
        &move |ar: &mut Arena, j: usize| {
            if m == 0 && j == 0 {
                let g = ar.euler_gamma;
                return ar.neg(g);
            }
            let z = zeta_int(ar, m + j + 1);
            let mut c =
                Q::from_integer(factorial((m + j) as u64)) / Q::from_integer(factorial(j as u64));
            if (m + j + 1) % 2 == 1 {
                c = -c;
            }
            let ce = rat_expr(ar, c);
            let p = ar.mul(&[ce, z]);
            eval::eval(ar, p)
        },
        w,
    )
}

/// For an argument `a → +∞` in a log-extended expansion from above (a pole
/// with a positive leading coefficient): `(ln a, 1/a, p)` with `p` the
/// valuation of `1/a`.
fn at_plus_infinity(arena: &mut Arena, a: &TSeries, mode: Mode) -> Option<(TSeries, TSeries, i64)> {
    mode.log_var?;
    if mode.side != Side::Above {
        return None;
    }
    let a = a.clone().normalized(arena);
    let v = a.leading_exponent(arena)?;
    if v >= 0 {
        return None;
    }
    let c = a.coeff_at(arena, v);
    if constant_sign(arena, c) != Some(true) {
        return None;
    }
    let ln_a = apply_ln(arena, &a, mode).ok()?;
    let inv = TSeries::inverse(arena, &a)?;
    Some((ln_a, inv, -v))
}

/// For an argument `a → +∞` in an expansion from above (a pole with a
/// positive leading coefficient): `(1/a, p)` with `p` the valuation of
/// `1/a` (no logarithm needed).
fn pole_to_plus_infinity(arena: &mut Arena, a: &TSeries, mode: Mode) -> Option<(TSeries, i64)> {
    if mode.side != Side::Above {
        return None;
    }
    let a = a.clone().normalized(arena);
    let v = a.leading_exponent(arena)?;
    if v >= 0 || constant_sign(arena, a.coeff_at(arena, v)) != Some(true) {
        return None;
    }
    let inv = TSeries::inverse(arena, &a)?;
    Some((inv, -v))
}

/// Number of terms `K` of an asymptotic series in `1/a` (valuation `p`) so
/// that its remainder lies beyond the working precision `target`.
fn asymptotic_terms(target: i64, p: i64) -> i64 {
    (target.max(1) / (2 * p) + 2).min(MAX_ASYMPTOTIC_TERMS)
}

/// `Σ_{k=1}^{K} c_k·inv^(2k+offset)` with `c_k = coeff(k)`.
fn odd_even_sum(
    arena: &mut Arena,
    inv: &TSeries,
    k_max: i64,
    offset: i64,
    coeff: &dyn Fn(usize) -> Q,
) -> Option<TSeries> {
    let inv2 = TSeries::mul(arena, inv, inv);
    let mut power = TSeries::pow_int(arena, inv, 2 + offset)?;
    let mut acc: Option<TSeries> = None;
    for k in 1..=k_max {
        let c = coeff(k as usize);
        if !c.is_zero() {
            let term = scale_rat(arena, &power, c);
            acc = Some(match acc {
                None => term,
                Some(s) => TSeries::add(arena, &s, &term),
            });
        }
        power = TSeries::mul(arena, &power, &inv2);
    }
    acc
}

/// Expansion of a special function whose argument tends to one of its
/// singular points (see the section comment).
fn special_series(
    arena: &mut Arena,
    node: &ExprNode,
    mode: Mode,
    cache: &FxHashMap<ExprId, Option<TSeries>>,
) -> Result<TSeries, Obstruction> {
    use crate::base::libfn::LibFn;
    let unknown = Err(Obstruction::Unknown);
    match *node {
        ExprNode::Erfc(a) => {
            let e = apply_fn(arena, FnKind::Erf, &child(cache, a)?).ok_or(Obstruction::Unknown)?;
            let one = constant_like(arena, arena.one, &e);
            let ne = TSeries::neg(arena, &e);
            Ok(TSeries::add(arena, &one, &ne))
        }
        ExprNode::Gamma(a) | ExprNode::Digamma(a) => {
            let a = child(cache, a)?;
            if matches!(node, ExprNode::Digamma(_))
                && let Some((ln_a, inv, p)) = at_plus_infinity(arena, &a, mode)
            {
                // ψ(a) = ln a − 1/(2a) − Σ B₂ₖ/(2k a^{2k}) (DLMF 5.11.2)
                let k_max = asymptotic_terms(a.known - a.shift, p);
                let half = scale_rat(arena, &inv, Q::new(BigInt::from(-1), BigInt::from(2)));
                let mut acc = TSeries::add(arena, &ln_a, &half);
                let tail = odd_even_sum(arena, &inv, k_max, 0, &|k| {
                    -bernoulli(2 * k) / rat_i(2 * k as i64)
                })
                .ok_or(Obstruction::Unknown)?;
                acc = TSeries::add(arena, &acc, &tail);
                return Ok(acc.truncate_known(p * (2 * k_max + 2)));
            }
            let Some((u0, w)) = split_regular(arena, &a)? else {
                return unknown;
            };
            let Some(n) = nonpositive_integer(arena, u0) else {
                return unknown;
            };
            if n > MAX_POLE_SHIFT || w.coeffs.is_empty() {
                return Err(Obstruction::NoExpansion);
            }
            if matches!(node, ExprNode::Gamma(_)) {
                // Γ(−n + w) = Γ(1 + w) / Π_{k=0}^{n} (w − n + k)
                let mut acc = gamma_one_plus(arena, &w);
                for k in 0..=n {
                    let r = recip_shifted(arena, &w, k - n, 1).ok_or(Obstruction::NoExpansion)?;
                    acc = TSeries::mul(arena, &acc, &r);
                }
                Ok(acc)
            } else {
                // ψ(−n + w) = ψ(1 + w) − Σ_{k=0}^{n} 1/(w − n + k)
                let mut acc = polygamma_one_plus(arena, 0, &w);
                for k in 0..=n {
                    let r = recip_shifted(arena, &w, k - n, 1).ok_or(Obstruction::NoExpansion)?;
                    let nr = TSeries::neg(arena, &r);
                    acc = TSeries::add(arena, &acc, &nr);
                }
                Ok(acc)
            }
        }
        ExprNode::Polygamma(m, a) => {
            let Some(m) = arena
                .as_num(m)
                .filter(|r| r.is_integer() && r.is_positive())
                .and_then(|r| r.to_integer().to_usize())
                .filter(|&m| m <= MAX_INT_POWER as usize)
            else {
                return unknown;
            };
            let a = child(cache, a)?;
            if let Some((inv, p)) = pole_to_plus_infinity(arena, &a, mode) {
                // ψ⁽ᵐ⁾(a) = (−1)^{m+1} [(m−1)!/a^m + m!/(2a^{m+1})
                //   + Σ B₂ₖ (2k+m−1)!/((2k)! a^{2k+m})]  (DLMF 5.15.8)
                let k_max = asymptotic_terms(a.known - a.shift, p);
                let mi = m as i64;
                let sign = if m % 2 == 1 { Q::one() } else { -Q::one() };
                let lead = TSeries::pow_int(arena, &inv, mi).ok_or(Obstruction::Unknown)?;
                let lead = scale_rat(
                    arena,
                    &lead,
                    &sign * Q::from_integer(factorial(m as u64 - 1)),
                );
                let second = TSeries::pow_int(arena, &inv, mi + 1).ok_or(Obstruction::Unknown)?;
                let second = scale_rat(
                    arena,
                    &second,
                    &sign * Q::from_integer(factorial(m as u64)) / rat_i(2),
                );
                let mut acc = TSeries::add(arena, &lead, &second);
                let tail = odd_even_sum(arena, &inv, k_max, mi, &|k| {
                    &sign * bernoulli(2 * k) * Q::from_integer(factorial((2 * k + m - 1) as u64))
                        / Q::from_integer(factorial(2 * k as u64))
                })
                .ok_or(Obstruction::Unknown)?;
                acc = TSeries::add(arena, &acc, &tail);
                return Ok(acc.truncate_known(p * (2 * k_max + 2 + mi)));
            }
            let Some((u0, w)) = split_regular(arena, &a)? else {
                return unknown;
            };
            let Some(n) = nonpositive_integer(arena, u0) else {
                return unknown;
            };
            if n > MAX_POLE_SHIFT || w.coeffs.is_empty() {
                return Err(Obstruction::NoExpansion);
            }
            // ψ⁽ᵐ⁾(−n + w) = ψ⁽ᵐ⁾(1 + w) − (−1)ᵐ m! Σ_{k=0}^{n} (w − n + k)^{−m−1}
            let mut acc = polygamma_one_plus(arena, m, &w);
            let mut c = Q::from_integer(factorial(m as u64));
            if m.is_multiple_of(2) {
                c = -c;
            }
            for k in 0..=n {
                let r = recip_shifted(arena, &w, k - n, m as i64 + 1)
                    .ok_or(Obstruction::NoExpansion)?;
                let t = scale_rat(arena, &r, c.clone());
                acc = TSeries::add(arena, &acc, &t);
            }
            Ok(acc)
        }
        ExprNode::Zeta(a) => {
            let a = child(cache, a)?;
            let Some((u0, w)) = split_regular(arena, &a)? else {
                return unknown;
            };
            if !is_one_const(arena, u0) {
                return unknown;
            }
            // ζ(1 + w) = 1/w + γ + O(w): the higher coefficients are
            // Stieltjes constants, which have no representation here, so
            // the expansion is exact only below the valuation of w.
            let v = w.leading_exponent(arena).ok_or(Obstruction::NoExpansion)?;
            let inv = TSeries::inverse(arena, &w).ok_or(Obstruction::NoExpansion)?;
            let g = TSeries::constant(arena, arena.euler_gamma, v);
            Ok(TSeries::add(arena, &inv, &g).truncate_known(v))
        }
        ExprNode::LogGamma(a) => {
            let a = child(cache, a)?;
            if let Some((ln_a, inv, p)) = at_plus_infinity(arena, &a, mode) {
                // Stirling: (a − ½) ln a − a + ½ ln 2π + Σ B₂ₖ/(2k(2k−1) a^{2k−1})
                let k_max = asymptotic_terms(a.known - a.shift, p);
                let half = rat_expr(arena, Q::new(BigInt::from(-1), BigInt::from(2)));
                let hc = constant_like(arena, half, &a);
                let a_half = TSeries::add(arena, &a, &hc);
                let main = TSeries::mul(arena, &a_half, &ln_a);
                let na = TSeries::neg(arena, &a);
                // ½ ln 2π, spelled ½ ln 2 + ½ ln π (SymPy's form in limits),
                // so that its exponential simplifies to √2·√π.
                let two_pi = {
                    let two = arena.int(2);
                    let pi = arena.pi;
                    let l2 = arena.ln(two);
                    let lp = arena.ln(pi);
                    let s = arena.add(&[l2, lp]);
                    let h = rat_expr(arena, Q::new(BigInt::one(), BigInt::from(2)));
                    let p = arena.mul(&[h, s]);
                    let p = crate::transforms::expand::expand(arena, p);
                    eval::eval(arena, p)
                };
                let c = constant_like(arena, two_pi, &main);
                let mut acc = TSeries::add(arena, &main, &na);
                acc = TSeries::add(arena, &acc, &c);
                let tail = odd_even_sum(arena, &inv, k_max, -1, &|k| {
                    bernoulli(2 * k) / rat_i((2 * k * (2 * k - 1)) as i64)
                })
                .ok_or(Obstruction::Unknown)?;
                acc = TSeries::add(arena, &acc, &tail);
                return Ok(acc.truncate_known(p * (2 * k_max + 1)));
            }
            let Some((u0, _)) = split_regular(arena, &a)? else {
                return unknown;
            };
            if !arena.is_zero_structural(u0) {
                return unknown;
            }
            // ln Γ(a) = ln Γ(1 + a) − ln a for a → 0⁺.
            let a_n = a.clone().normalized(arena);
            let lead = a_n
                .leading_exponent(arena)
                .ok_or(Obstruction::NoExpansion)?;
            let c = a_n.coeff_at(arena, lead);
            if mode.log_var.is_none() || constant_sign(arena, c) != Some(true) {
                return Err(Obstruction::NoExpansion);
            }
            let ln_a = apply_ln(arena, &a_n, mode).map_err(|_| Obstruction::NoExpansion)?;
            let lg = ln_gamma_one_plus(arena, &a_n);
            let nl = TSeries::neg(arena, &ln_a);
            Ok(TSeries::add(arena, &lg, &nl))
        }
        ExprNode::Si(a) | ExprNode::Ei(a) | ExprNode::Ci(a) => {
            let a = child(cache, a)?;
            let Some((u0, _)) = split_regular(arena, &a)? else {
                return unknown;
            };
            if !arena.is_zero_structural(u0) {
                return unknown;
            }
            let kind = match node {
                ExprNode::Si(_) => IntegralKind::Si,
                ExprNode::Ei(_) => IntegralKind::Ei,
                _ => IntegralKind::Ci,
            };
            integral_at_zero(arena, kind, &a, mode)
        }
        ExprNode::Apply(f, ref args)
            if args.len() == 1
                && matches!(arena.symbol_name(f), ERFCX_ASYMPTOTIC | EI_ASYMPTOTIC) =>
        {
            let erfcx = arena.symbol_name(f) == ERFCX_ASYMPTOTIC;
            let a = child(cache, args[0])?.normalized(arena);
            let v = a.leading_exponent(arena).ok_or(Obstruction::NoExpansion)?;
            if v == 0
                && a.coeffs
                    .iter()
                    .skip(1)
                    .all(|&c| arena.is_zero_structural(c))
            {
                // A constant argument `c + O(ωᴷ)` (in a log-extended
                // expansion `ln(1/ω)` is the constant `−ln ω`): the function
                // is `f(c) + O(ωᴷ)`, its derivative being bounded away from
                // the logarithmic singularity of `e^{−z} Ei z` at `0`.
                let c0 = a.coeffs[0];
                if erfcx || !is_zero_const(arena, c0) {
                    let c = arena.intern(ExprNode::Apply(f, smallvec::smallvec![c0]));
                    return Ok(TSeries::constant(arena, c, a.known));
                }
                return unknown;
            }
            if v >= 0 || mode.side != Side::Above {
                return unknown;
            }
            let positive = constant_sign(arena, a.coeff_at(arena, v));
            if positive != Some(true) && (erfcx || positive.is_none()) {
                return unknown;
            }
            let inv = TSeries::inverse(arena, &a).ok_or(Obstruction::NoExpansion)?;
            let p = -v;
            let k_max = if erfcx {
                asymptotic_terms(a.known - a.shift, p)
            } else {
                ((a.known - a.shift).max(1) / p + 2).min(MAX_ASYMPTOTIC_TERMS)
            };
            let mut acc: Option<TSeries> = None;
            let mut power = inv.clone();
            let step = if erfcx {
                TSeries::mul(arena, &inv, &inv)
            } else {
                inv.clone()
            };
            let mut c = Q::one();
            for k in 0..=k_max {
                if k > 0 {
                    c = if erfcx {
                        -c * rat_i(2 * k - 1) / rat_i(2)
                    } else {
                        c * rat_i(k)
                    };
                }
                let term = scale_rat(arena, &power, c.clone());
                acc = Some(match acc {
                    None => term,
                    Some(s) => TSeries::add(arena, &s, &term),
                });
                power = TSeries::mul(arena, &power, &step);
            }
            let mut s = acc.ok_or(Obstruction::Unknown)?;
            // Remainder: the first omitted term.
            let omitted = if erfcx { 2 * k_max + 3 } else { k_max + 2 };
            s = s.truncate_known(p * omitted);
            if erfcx {
                let rsp = {
                    let pi = arena.pi;
                    let sp = arena.sqrt(pi);
                    let one = arena.one;
                    let q = arena.div(one, sp);
                    eval::eval(arena, q)
                };
                s = TSeries::scale(arena, &s, rsp);
            }
            Ok(s)
        }
        ExprNode::Apply(f, ref args) => {
            let kind = match arena.lib_fn(f) {
                Some(LibFn::Shi) => IntegralKind::Shi,
                Some(LibFn::Chi) => IntegralKind::Chi,
                _ => return unknown,
            };
            let a = child(cache, args[0])?;
            let Some((u0, _)) = split_regular(arena, &a)? else {
                return unknown;
            };
            if !arena.is_zero_structural(u0) {
                return unknown;
            }
            integral_at_zero(arena, kind, &a, mode)
        }
        ExprNode::Li(a) => {
            let a = child(cache, a)?;
            let Some((u0, _)) = split_regular(arena, &a)? else {
                return unknown;
            };
            if !is_one_const(arena, u0) {
                return unknown;
            }
            // li(a) = Ei(ln a), ln a → 0.
            let l = apply_ln(arena, &a, mode)?;
            integral_at_zero(arena, IntegralKind::Ei, &l, mode)
        }
        _ => unknown,
    }
}

/// The integral functions expanded at a zero of their argument.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum IntegralKind {
    Si,
    Shi,
    Ei,
    Ci,
    Chi,
}

/// `Si`, `Shi` (entire), `Ei(a) = γ + ln|a| + Σ aᵏ/(k·k!)`,
/// `Ci(a) = γ + ln a + Σ (−1)ᵐ a²ᵐ/(2m·(2m)!)`, `Chi` likewise without the
/// signs, for `a → 0` (DLMF 6.6.2, 6.6.5, 6.6.6; SymPy's `Ci(−x) = Ci(x) + iπ`
/// is the principal `ln` of a negative `a`, and `Ei` of a real argument is
/// real).  The logarithmic ones need a log-extended expansion.
fn integral_at_zero(
    arena: &mut Arena,
    kind: IntegralKind,
    a: &TSeries,
    mode: Mode,
) -> Result<TSeries, Obstruction> {
    let a = a.clone().normalized(arena);
    if a.coeffs.is_empty() {
        return Err(Obstruction::NoExpansion);
    }
    let coefficient = move |ar: &mut Arena, k: usize| -> ExprId {
        let kq = rat_i(k as i64);
        let fact = Q::from_integer(factorial(k as u64));
        let odd = k % 2 == 1;
        let m = k / 2;
        let alt = if m.is_multiple_of(2) {
            Q::one()
        } else {
            -Q::one()
        };
        let c = match kind {
            IntegralKind::Si if odd => alt / (kq * fact),
            IntegralKind::Shi if odd => Q::one() / (kq * fact),
            IntegralKind::Ei if k > 0 => Q::one() / (kq * fact),
            IntegralKind::Ci if k > 0 && !odd => alt / (kq * fact),
            IntegralKind::Chi if k > 0 && !odd => Q::one() / (kq * fact),
            IntegralKind::Ei | IntegralKind::Ci | IntegralKind::Chi if k == 0 => {
                return ar.euler_gamma;
            }
            _ => Q::zero(),
        };
        rat_expr(ar, c)
    };
    let series = TSeries::compose(arena, &coefficient, &a);
    if matches!(kind, IntegralKind::Si | IntegralKind::Shi) {
        return Ok(series);
    }
    if mode.log_var.is_none() {
        return Err(Obstruction::NoExpansion);
    }
    let log_arg = if kind == IntegralKind::Ei {
        abs_series(arena, &a, mode.side)?
    } else {
        a
    };
    let ln = apply_ln(arena, &log_arg, mode).map_err(|_| Obstruction::NoExpansion)?;
    Ok(TSeries::add(arena, &series, &ln))
}

/// `K₀(u)` and `Y₀(u)` at a zero of `u` in a log-extended expansion (DLMF
/// 10.31.2, 10.8.2):
/// `K₀(z) = −(ln(z/2) + γ) I₀(z) + Σ_{k≥1} H_k (z/2)^{2k}/(k!)²`,
/// `Y₀(z) = (2/π)[(ln(z/2) + γ) J₀(z) + Σ_{k≥1} (−1)^{k+1} H_k (z/2)^{2k}/(k!)²]`,
/// for `u → 0⁺` (positive leading coefficient).  `Ok(None)` when the rule
/// does not apply.
fn log_bessel_series(
    arena: &mut Arena,
    f: crate::base::node::SymbolId,
    order: ExprId,
    u: &TSeries,
    var: ExprId,
    mode: Mode,
) -> Result<Option<TSeries>, Obstruction> {
    use crate::base::libfn::LibFn;
    let lib = arena.lib_fn(f);
    let is_k = match lib {
        Some(LibFn::BesselK) => true,
        Some(LibFn::BesselY) => false,
        _ => return Ok(None),
    };
    if mode.log_var.is_none()
        || walk::contains(arena, order, var)
        || !arena.is_zero_structural(order)
    {
        return Ok(None);
    }
    let u = u.clone().normalized(arena);
    let Some(k) = u.leading_exponent(arena) else {
        return Err(Obstruction::NoExpansion);
    };
    if k <= 0 {
        return Ok(None);
    }
    if constant_sign(arena, u.coeff_at(arena, k)) != Some(true) {
        return Err(Obstruction::NoExpansion);
    }
    let alternate = !is_k;
    let bessel0 = move |ar: &mut Arena, m: usize| -> ExprId {
        if m % 2 == 1 {
            return ar.zero;
        }
        let kk = (m / 2) as u64;
        let f = factorial(kk);
        let mut c = Q::new(BigInt::one(), &f * &f * (BigInt::one() << (2 * kk)));
        if alternate && kk % 2 == 1 {
            c = -c;
        }
        rat_expr(ar, c)
    };
    let tail = move |ar: &mut Arena, m: usize| -> ExprId {
        if m % 2 == 1 || m == 0 {
            return ar.zero;
        }
        let kk = (m / 2) as u64;
        let h: Q = (1..=kk)
            .map(|j| Q::new(BigInt::one(), BigInt::from(j)))
            .sum();
        let f = factorial(kk);
        let mut c = h / Q::from_integer(&f * &f * (BigInt::one() << (2 * kk)));
        if alternate && kk.is_multiple_of(2) {
            c = -c;
        }
        rat_expr(ar, c)
    };
    let i0 = TSeries::compose(arena, &bessel0, &u);
    let s = TSeries::compose(arena, &tail, &u);
    let ln_u = apply_ln(arena, &u, mode).map_err(|_| Obstruction::NoExpansion)?;
    // ln(u/2) + γ
    let shift = {
        let two = arena.int(2);
        let l2 = arena.ln(two);
        let g = arena.euler_gamma;
        let nl2 = arena.neg(l2);
        let s = arena.add(&[g, nl2]);
        eval::eval(arena, s)
    };
    let sc = constant_like(arena, shift, &ln_u);
    let lg = TSeries::add(arena, &ln_u, &sc);
    let prod = TSeries::mul(arena, &lg, &i0);
    if is_k {
        let np = TSeries::neg(arena, &prod);
        Ok(Some(TSeries::add(arena, &np, &s)))
    } else {
        let sum = TSeries::add(arena, &prod, &s);
        let two_over_pi = {
            let two = arena.int(2);
            let pi = arena.pi;
            let q = arena.div(two, pi);
            eval::eval(arena, q)
        };
        Ok(Some(TSeries::scale(arena, &sum, two_over_pi)))
    }
}

/// Taylor coefficients by repeated differentiation at `0`.
///
/// When direct substitution of `0` into a derivative is singular or
/// indeterminate (`atan(1/t)` at `t = 0` gives `atan(zoo)`), the
/// coefficient is the *limit* of that derivative at `0` — one-sided
/// (`0⁺`) for expansions at `±∞`, two-sided otherwise.  A limit that does
/// not exist or is not finite means there is no Taylor expansion.
fn taylor_by_differentiation(
    arena: &mut Arena,
    expr: ExprId,
    var: ExprId,
    n: i64,
    mode: Mode,
) -> Option<TSeries> {
    let side = mode.side;
    let zero = arena.zero;
    let mut coeffs = Vec::with_capacity(n.max(0) as usize);
    let mut current = expr;
    let mut factorial = Q::one();
    // At a discontinuity of `sign`, `H`, `⌊·⌋`, … the value at the point is
    // not the one-sided limit: `tan(sign x)` substituted `sign 0 = 0` and
    // expanded to `0` from both sides.  With such nodes every coefficient
    // is the limit from the side of the expansion.
    let directional = crate::calculus::limit::has_directional_nodes(arena, expr, var);
    // Coefficients `0..n` are claimed exact up to `O(xⁿ)`, which needs the
    // `n`-th derivative bounded near the point (Taylor's remainder): its
    // limit is taken too and must be finite.  Without it `1/ln x` "expanded"
    // to `0` (`f(0⁺) = 0`, but `1/ln x` is not `O(x)`), and
    // `−2 − (x − 3/4)/ln x` to `−2`.
    for k in 0..=n {
        // A derivative that stays formal (`d re(u)/dx` for a variable not
        // declared real) has no value at the point; asking `limit` for one
        // sent Gruntz after `Subs(Derivative(…))` for seconds.
        if walk::has_unevaluated(arena, current) {
            return None;
        }
        // Derivatives of nested functions grow exponentially; the limit
        // taken of each (below) expands them, and `atanh(acosh x)³` at
        // `x = −1` never returned.
        if crate::transforms::pattern::tree_size_capped(arena, current, MAX_FALLBACK_SIZE + 1)
            > MAX_FALLBACK_SIZE
        {
            return None;
        }
        // `safe_substitute` refuses a singular sub-expression instead of
        // letting canonicalisation fold `0·sinh(ln 0)` to `0` (the series
        // of `x·sinh(ln x)` lost its constant term −1/2 that way).
        let value =
            crate::calculus::limit::safe_substitute(arena, current, var, zero).unwrap_or(arena.nan);
        if k == n {
            // Only a definitely unbounded `n`-th derivative refuses; one whose
            // limit cannot be determined (`d³(x·Γ(x))`) is given the benefit
            // of the doubt, as before.
            if !directional && is_finite_constant(arena, value, var) {
                break;
            }
            if !mode.limits {
                return None;
            }
            use crate::calculus::limit::Direction;
            let dirs: &[Direction] = match side {
                Side::Above => &[Direction::Right],
                Side::Below => &[Direction::Left],
                Side::Both => &[Direction::Right, Direction::Left],
            };
            for &dir in dirs {
                if let Ok(lim) =
                    crate::calculus::limit::limit_dir_generic(arena, current, var, zero, dir)
                    && !is_finite_constant(arena, lim, var)
                {
                    return None;
                }
            }
            break;
        }
        let value = if !directional && is_finite_constant(arena, value, var) {
            value
        } else {
            if !mode.limits {
                return None;
            }
            let dir = match side {
                Side::Above => crate::calculus::limit::Direction::Right,
                Side::Below => crate::calculus::limit::Direction::Left,
                Side::Both => crate::calculus::limit::Direction::Both,
            };
            let lim =
                crate::calculus::limit::limit_dir_generic(arena, current, var, zero, dir).ok()?;
            if !is_finite_constant(arena, lim, var) {
                return None;
            }
            lim
        };
        let coeff = if k == 0 {
            value
        } else {
            let inv = rat_expr(arena, Q::one() / &factorial);
            let p = arena.mul(&[value, inv]);
            eval::eval(arena, p)
        };
        coeffs.push(coeff);
        current = crate::transforms::diff::diff(arena, current, var);
        factorial *= rat_i(k + 1);
    }
    Some(TSeries {
        shift: 0,
        known: n,
        coeffs,
    })
}

/// A finite, fully evaluated constant (no `±∞`, `zoo`, `NaN`, hidden
/// singularity such as `ln 0`, unevaluated node, or occurrence of `var`).
fn is_finite_constant(arena: &Arena, value: ExprId, var: ExprId) -> bool {
    value != arena.infinity
        && value != arena.neg_infinity
        && crate::calculus::limit::is_valid_limit_value(arena, value, var)
}

// ═══════════════════════════════════════════════════════════════════════════
// Tests
// ═══════════════════════════════════════════════════════════════════════════

#[cfg(test)]
mod tests {
    use super::*;
    use crate::base::arena::Arena;

    fn sym(a: &mut Arena, name: &str) -> ExprId {
        a.symbol(name)
    }

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

    /// Compare a series numerically against the target function at a small point.
    fn check_close(a: &mut Arena, series: ExprId, target: ExprId, x: ExprId, at: f64, tol: f64) {
        let nid = a.intern_num(Ratio::from_float(at).unwrap());
        let pt = a.intern(ExprNode::Num(nid));
        let s = subs::subs(a, series, x, pt);
        let t = subs::subs(a, target, x, pt);
        let sv = crate::transforms::evalf::eval_const_f64(a, s).unwrap();
        let tv = crate::transforms::evalf::eval_const_f64(a, t).unwrap();
        assert!(
            (sv - tv).abs() < tol,
            "series {} vs target {} at {at}: {sv} vs {tv}",
            display(a, series),
            display(a, target)
        );
    }

    #[test]
    fn series_constant() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let five = a.int(5);
        let zero = a.zero;
        let result = series(&mut a, five, x, zero, 3).unwrap();
        assert_eq!(display(&a, result), "5");
    }

    #[test]
    fn series_x_around_zero() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let result = series(&mut a, x, x, zero, 3).unwrap();
        assert_eq!(display(&a, result), "x");
    }

    #[test]
    fn series_polynomial_is_exact() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let three = a.int(3);
        let x3 = a.pow(x, three);
        let zero = a.zero;
        let result = series(&mut a, x3, x, zero, 5).unwrap();
        assert_eq!(display(&a, result), "x^3");
    }

    #[test]
    fn series_order_zero_is_zero() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let result = series(&mut a, x, x, zero, 0).unwrap();
        assert_eq!(result, a.zero);
    }

    #[test]
    fn series_exp_sin_cos_fast_paths() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let e = a.exp(x);
        let s = series(&mut a, e, x, zero, 5).unwrap();
        assert_eq!(display(&a, s), "1/24*x^4 + 1/6*x^3 + 1/2*x^2 + x + 1");
        let sn = a.sin(x);
        let s = series(&mut a, sn, x, zero, 6).unwrap();
        assert_eq!(display(&a, s), "1/120*x^5 - 1/6*x^3 + x");
        let c = a.cos(x);
        let s = series(&mut a, c, x, zero, 5).unwrap();
        assert_eq!(display(&a, s), "1/24*x^4 - 1/2*x^2 + 1");
    }

    #[test]
    fn series_composition_sin_x_squared() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let two = a.int(2);
        let x2 = a.pow(x, two);
        let f = a.sin(x2);
        let s = series(&mut a, f, x, zero, 8).unwrap();
        assert_eq!(display(&a, s), "-1/6*x^6 + x^2");
    }

    #[test]
    fn series_sin_over_x_is_laurent_free() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let sn = a.sin(x);
        let f = a.div(sn, x);
        let s = series(&mut a, f, x, zero, 5).unwrap();
        assert_eq!(display(&a, s), "1/120*x^4 - 1/6*x^2 + 1");
    }

    #[test]
    fn series_pole_gives_laurent() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let one = a.one;
        // 1/(x(1−x)) = 1/x + 1 + x + x² + …
        let omx = a.sub(one, x);
        let den = a.mul(&[x, omx]);
        let f = a.div(one, den);
        let s = series(&mut a, f, x, zero, 3).unwrap();
        assert_eq!(display(&a, s), "x^2 + x + 1/x + 1");
    }

    #[test]
    fn series_tan_asin_erf_lambertw() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let t = a.tan(x);
        let s = series(&mut a, t, x, zero, 8).unwrap();
        assert_eq!(display(&a, s), "17/315*x^7 + 2/15*x^5 + 1/3*x^3 + x");
        let th = a.tanh(x);
        let s = series(&mut a, th, x, zero, 6).unwrap();
        assert_eq!(display(&a, s), "2/15*x^5 - 1/3*x^3 + x");
        let asn = a.asin(x);
        let s = series(&mut a, asn, x, zero, 6).unwrap();
        assert_eq!(display(&a, s), "3/40*x^5 + 1/6*x^3 + x");
        let w = a.lambertw(x);
        let s = series(&mut a, w, x, zero, 5).unwrap();
        assert_eq!(display(&a, s), "-8/3*x^4 + 3/2*x^3 - x^2 + x");
        let e = a.erf(x);
        let s = series(&mut a, e, x, zero, 4).unwrap();
        check_close(&mut a, s, e, x, 0.1, 1e-5);
    }

    #[test]
    fn series_ln_and_binomial() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let one = a.one;
        let opx = a.add(&[one, x]);
        let l = a.ln(opx);
        let s = series(&mut a, l, x, zero, 4).unwrap();
        assert_eq!(display(&a, s), "1/3*x^3 - 1/2*x^2 + x");
        let half = a.rational(1, 2);
        let sq = a.pow(opx, half);
        let s = series(&mut a, sq, x, zero, 4).unwrap();
        assert_eq!(display(&a, s), "1/16*x^3 - 1/8*x^2 + 1/2*x + 1");
        // ln(2 + x) = ln 2 + x/2 − x²/8 + …
        let two = a.int(2);
        let tpx = a.add(&[two, x]);
        let l2 = a.ln(tpx);
        let s = series(&mut a, l2, x, zero, 3).unwrap();
        check_close(&mut a, s, l2, x, 0.01, 1e-6);
    }

    #[test]
    fn series_around_nonzero_point() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let one = a.one;
        let e = a.exp(x);
        let s = series(&mut a, e, x, one, 4).unwrap();
        check_close(&mut a, s, e, x, 1.01, 1e-8);
    }

    #[test]
    fn puiseux_is_rejected() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let sq = a.sqrt(x);
        let sn = a.sin(x);
        let f = a.mul(&[sq, sn]);
        assert!(series(&mut a, f, x, zero, 5).is_err());
        let l = a.ln(x);
        assert!(series(&mut a, l, x, zero, 5).is_err());
    }

    #[test]
    fn puiseux_is_rejected_below_the_singular_derivative() {
        // x^(5/2): the first three derivatives vanish at 0, so a
        // differentiation fallback at order 3 would return the wrong
        // polynomial 0.  The refusal must be definite at every order.
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let e = a.rational(5, 2);
        let f = a.pow(x, e);
        assert!(series(&mut a, f, x, zero, 3).is_err());
        assert!(series(&mut a, f, x, zero, 6).is_err());
    }

    #[test]
    fn abs_of_even_valuation_is_signed_argument() {
        // |x²| = x²; |x − 1| = 1 − x near 0; |−x² + x³| = x² − x³.
        // sympy: series(Abs(x**2), x, 0, 6) == x**2;
        //        series(Abs(x - 1), x, 0, 6) == 1 - x
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let two = a.int(2);
        let x2 = a.pow(x, two);
        let f = a.abs(x2);
        let s = series(&mut a, f, x, zero, 6).unwrap();
        assert_eq!(s, x2);
        let one = a.one;
        let xm1 = a.sub(x, one);
        let g = a.abs(xm1);
        let s = series(&mut a, g, x, zero, 6).unwrap();
        let expected = a.sub(one, x);
        assert_eq!(s, expected);
        let three = a.int(3);
        let x3 = a.pow(x, three);
        let h_arg = a.sub(x3, x2);
        let h = a.abs(h_arg);
        let s = series(&mut a, h, x, zero, 6).unwrap();
        let expected = a.sub(x2, x3);
        assert_eq!(s, expected);
    }

    #[test]
    fn abs_of_odd_valuation_needs_both_sides_to_agree() {
        // |x| and |sin x| have a kink; cos|x| = cos x and |x|² = x² do not.
        // sympy: series(cos(Abs(x)), x, 0, 6) == x**4/24 - x**2/2 + 1
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let ax = a.abs(x);
        assert!(series(&mut a, ax, x, zero, 6).is_err());
        let sn = a.sin(x);
        let asn = a.abs(sn);
        assert!(series(&mut a, asn, x, zero, 6).is_err());
        let ex = a.exp(ax);
        assert!(series(&mut a, ex, x, zero, 6).is_err());
        let cs = a.cos(ax);
        let s = series(&mut a, cs, x, zero, 6).unwrap();
        assert_eq!(display(&a, s), "1/24*x^4 - 1/2*x^2 + 1");
        let two = a.int(2);
        let ax2 = a.pow(ax, two);
        let s = series(&mut a, ax2, x, zero, 6).unwrap();
        let x2 = a.pow(x, two);
        assert_eq!(s, x2);
    }

    #[test]
    fn abs_at_infinity_is_one_sided() {
        // sympy: series(Abs(x), x, oo, 3) == x; series(Abs(x), x, -oo, 3) == -x
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let ax = a.abs(x);
        let s = series_at_infinity(&mut a, ax, x, 3, false).unwrap();
        assert_eq!(s, x);
        let s = series_at_infinity(&mut a, ax, x, 3, true).unwrap();
        let neg_x = a.neg(x);
        assert_eq!(s, neg_x);
    }

    #[test]
    fn series_at_infinity_rational() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let one = a.one;
        let xp1 = a.add(&[x, one]);
        let f = a.div(x, xp1);
        let s = series_at_infinity(&mut a, f, x, 3, false).unwrap();
        assert_eq!(display(&a, s), "x^(-2) - 1/x + 1");
        // sqrt(x²+1) − x  ~  1/(2x) − 1/(8x³)
        let two = a.int(2);
        let x2 = a.pow(x, two);
        let x2p1 = a.add(&[x2, one]);
        let root = a.sqrt(x2p1);
        let g = a.sub(root, x);
        let s = series_at_infinity(&mut a, g, x, 4, false).unwrap();
        assert_eq!(display(&a, s), "-1/8*x^(-3) + 1/(2*x)");
    }

    #[test]
    fn exp_pow_with_variable_exponent() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let one = a.one;
        // (1+x)^x = 1 + x² − x³/2 + …
        let opx = a.add(&[one, x]);
        let f = a.pow(opx, x);
        let s = series(&mut a, f, x, zero, 4).unwrap();
        assert_eq!(display(&a, s), "-1/2*x^3 + x^2 + 1");
    }

    #[test]
    fn large_integer_powers_use_binomial_series() {
        // Beyond MAX_INT_POWER repeated multiplication is replaced by the
        // binomial series; the result must be exact and fast.
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let one = a.one;
        let opx = a.add(&[one, x]);
        let big = a.int(MAX_INT_POWER + 36); // 100
        let f = a.pow(opx, big);
        let start = std::time::Instant::now();
        let s = series(&mut a, f, x, zero, 3).unwrap();
        assert!(start.elapsed().as_secs_f64() < 1.0);
        assert_eq!(display(&a, s), "4950*x^2 + 100*x + 1");
        // negative: (1+x)^(-70) = 1 − 70x + 2485x²
        let neg = a.int(-(MAX_INT_POWER + 6));
        let f = a.pow(opx, neg);
        let s = series(&mut a, f, x, zero, 3).unwrap();
        assert_eq!(display(&a, s), "2485*x^2 - 70*x + 1");
        // with a zero at the origin: (x + x²)^70 = x^70 + 70 x^71 + …
        let two = a.int(2);
        let x2 = a.pow(x, two);
        let base = a.add(&[x, x2]);
        let e70 = a.int(70);
        let f = a.pow(base, e70);
        let s = series(&mut a, f, x, zero, 72).unwrap();
        assert_eq!(display(&a, s), "70*x^71 + x^70");
        // odd valuation with a negative integer exponent is fine two-sided:
        // (x + x²)^(-65) = x^(-65) (1 + x)^(-65) = x^(-65) − 65 x^(-64) + …
        // (the pole of order 65 forces the precision-escalation retry)
        let em65 = a.int(-65);
        let f = a.pow(base, em65);
        let ts = expand_meromorphic(&mut a, f, x, 1).unwrap();
        assert_eq!(ts.shift(), -65);
        assert!(ts.known() >= 1);
        let c = ts.coefficient(&a, -65);
        assert_eq!(display(&a, c), "1");
        let c = ts.coefficient(&a, -64);
        assert_eq!(display(&a, c), "-65");
        let c = ts.coefficient(&a, -63);
        assert_eq!(display(&a, c), "2145");
    }

    #[test]
    fn low_order_requests_still_see_the_pole() {
        // At order 1 the variable itself would be truncated away at the
        // requested precision; the engine must widen its working window.
        let mut a = Arena::new();
        let x = sym(&mut a, "x");
        let zero = a.zero;
        let one = a.one;
        let inv = a.div(one, x);
        let s = series(&mut a, inv, x, zero, 1).unwrap();
        assert_eq!(display(&a, s), "1/x");
        // 1/(x⁵ + x⁶) = x⁻⁵ − x⁻⁴ + x⁻³ − …  (valuation 5 hidden at precision 4)
        let five = a.int(5);
        let six = a.int(6);
        let x5 = a.pow(x, five);
        let x6 = a.pow(x, six);
        let d = a.add(&[x5, x6]);
        let f = a.div(one, d);
        let s = series(&mut a, f, x, zero, 1).unwrap();
        assert_eq!(
            display(&a, s),
            "1/x + x^(-3) + x^(-5) - x^(-2) - x^(-4) - 1"
        );
        // A genuine failure is still reported after the bounded retries.
        let l = a.ln(x);
        let start = std::time::Instant::now();
        assert!(series(&mut a, l, x, zero, 1).is_err());
        assert!(start.elapsed().as_secs_f64() < 1.0);
    }

    #[test]
    fn tan_coefficients_match_bernoulli_formula() {
        // 1, 1/3, 2/15, 17/315, 62/2835
        let expected = [(1, 1), (3, 1), (5, 2), (7, 17), (9, 62)];
        let denoms = [1i64, 3, 15, 315, 2835];
        for (i, (n, num)) in expected.iter().enumerate() {
            let c = FnKind::Tan.rational_coefficient(*n as usize);
            assert_eq!(
                c,
                Q::new(BigInt::from(*num), BigInt::from(denoms[i])),
                "n={n}"
            );
        }
    }
}