ocas-calc 0.27.3

//! Trig-denominator power reduction, linear-numerator decomposition,
//! polynomial×trig closed forms, phase-shift normalization, and
//! polynomial×kernel ratio reduction.
//!
//! Five mechanisms, each declining honestly (`None`) when inapplicable:
//!
//! - **T1 [`integrate_trig_denom`]** — `∫ dx/(a + b·T(u))^n` for
//!   `T ∈ {sin, cos}`, `u = c + d·x` linear (symbolic coefficients allowed),
//!   `a, b` constant w.r.t. `x`, integer `2 ≤ n ≤ 6`. The classical
//!   recurrence, derived by differentiating `S(u)/D^(n−1)` with
//!   `D = a + b·T(u)` and `S` the complementary trig function:
//!
//!   ```text
//!   J_n = [ ε·b·S(u)/(d·D^(n−1)) + a(2n−3)·J_{n−1} − (n−2)·J_{n−2} ] / [ (n−1)(a²−b²) ]
//!   with  ε = −1, S = sin  when T = cos;   ε = +1, S = cos  when T = sin.
//!   ```
//!
//!   The `n = 1` base re-enters the chain via [`integrate_raw`] (Weierstrass
//!   and the symbolic-rational backend solve it); any `Integral` residue
//!   there declines the whole reduction. The singular case `a² = b²`
//!   declines — exactly when `a, b` are numeric, syntactically via
//!   `collect_terms` when symbolic.
//!
//! - **T2 [`integrate_trig_num_linear`]** — `(A + B·T(u))/(a + b·T(u))^n`
//!   splits into `A·J_n + B·K_n` with `K_n = ∫ T/D^n = (J_{n−1} − a·J_n)/b`,
//!   reusing T1. The numerator's trig function and argument must match the
//!   denominator's exactly; a product of two denominator powers declines.
//!
//! - **T3 [`integrate_poly_trig`]** — `∫ x^m·T(a·x + b) dx` (integer
//!   `1 ≤ m ≤ 6`) and `∫ (c + d·x)^k·T(a·x + b) dx` (`k ≤ 4`, binomially
//!   expanded into monomial terms) as direct closed forms by tabular
//!   reduction:
//!
//!   ```text
//!   ∫ x^m sin(ax+b) dx = Σ_{k=0..m} fall(m,k)·x^(m−k)·a^(−k−1)·s_k,
//!       s_k = −cos, +sin, +cos, −sin, …  (k mod 4 = 0,1,2,3)
//!   ∫ x^m cos(ax+b) dx = Σ_{k=0..m} fall(m,k)·x^(m−k)·a^(−k−1)·c_k,
//!       c_k = +sin, +cos, −sin, −cos, …
//!   ```
//!
//!   The closed form never re-enters the chain and contains no `Integral`,
//!   so it cannot ping-pong with the parts heuristic (the 0.27 lesson).
//!
//! - **P1 [`integrate_phase_shift`]** — *phase-shift normalization*. A
//!   denominator that mixes two kernels at one argument,
//!   `a + b·cos(u) + c·sin(u)` (`b·sin(u) + c·cos(u)` in either order and
//!   either order of the two summands), is the polar form of
//!   `a + R·cos(u − φ)` with
//!
//!   ```text
//!   R = √(b² + c²),   φ = atan(c/b).
//!   ```
//!
//!   The whole integrand is rewritten by replacing the denominator base with
//!   the shifted one and delegating to T1/T2/P2 with the argument
//!   `v = u − φ`; `v` is linear in `x` again because `φ` is a constant
//!   w.r.t. `x`. Branch conditions relied on: the identity
//!   `R·cos(u − atan(c/b)) = b·cos(u) + c·sin(u)` holds for **all** real
//!   `b, c` (a negative `b` is absorbed into `cos(φ) = b/R < 0`), so no
//!   quadrant split is needed and `atan`'s principal branch is sound; the
//!   branch `φ + π` (equivalently `R → −R`, same discriminant `a² − R²`) is
//!   only a fallback when the principal branch leaves an `Integral` residue
//!   behind. The `b < 0` case needs no separate treatment for the same
//!   reason. `a² = R²` (a squared singularity of the shifted denominator)
//!   still declines through the usual discriminant test, and the shifted
//!   denominator's discriminant is developed symbolically
//!   (`a² − (b² + c²)`), so no branch of the original integrand can slip
//!   through a numeric special case.
//!
//! - **P2 [`integrate_poly_kernel_ratio`]** —//!   `∫ C·P(x)·T(u)/(a + b·T(u))^n dx` for a polynomial `P` with
//!   `0 < deg P ≤ n − 1` and `deg P ≤ 4`, `T ∈ {sin, cos}` at the *same*
//!   linear argument as the denominator. With `D = a + b·T(u)` and `ε = +1`
//!   for `T = sin`, `ε = −1` for `T = cos` (so `S' = ε·T·du` for the
//!   complementary kernel `S`), one integration by parts
//!
//!   ```text
//!   ∫ P·T·D^(−n) dx
//!       = C·P·D^(1−n) + ∫ [ (C·a/b)·P − C·P'/(b·d·(n−1)) ]·T·D^(1−n) dx,
//!       C = −ε/b,
//!   ```
//!
//!   lowers the denominator power by one and `deg P` by one, so `deg P` steps
//!   leave a constant residual `∫ T·D^(−n+deg P)`. That base case has the
//!   closed form `C·S(u)/(d·(n−1))·D^(1−n)`; when the terminal power is
//!   exactly `1` the residual is instead T1's own `n = 1` base, which
//!   re-enters the chain and declines on a residue. The emitted result is a
//!   plain sum of products, so it cannot ping-pong with the parts heuristic.

use ocas_atom::normalize::normalize;
use ocas_atom::{Atom, AtomArena, AtomNode, Symbol};

use super::{
    contains_integral, integrate_raw, is_constant, is_fallback, linear_form, node_count, rat_atom,
};

/// Cap on the denominator power `n` in T1/T2.
const MAX_DENOM_POW: i64 = 6;
/// Cap on the monomial degree `m` in T3.
const MAX_POLY_DEG: i64 = 6;
/// Cap on `k` for `(c + d·x)^k` expansion in T3.
const MAX_LIN_POW: i64 = 4;
/// Cap on `deg P` in P2.
const MAX_RATIO_DEG: i64 = 4;
/// Node budget for the input integrand.
const MAX_NODES: usize = 200;
/// Master switch for the two 0.27.2 additions (P1 phase shift and P2
/// polynomial×kernel ratio).
///
/// Both are implemented and reachable, but their *emitted atoms* are
/// still being debugged: the closed-form base [`weierstrass_base`] and
/// the P2 recurrence have each been re-derived and checked against
/// numerical quadrature, yet the composed atom that reaches
/// `crate::diff` does not reproduce the integrand (the measured errors
/// are large and structural, not tolerance-level). Rather than let a
/// wrong antiderivative leave this module, the phase is disabled:
/// T1/T2/T3 keep exactly their 0.27.1 behaviour and P1/P2 decline.
/// Flip to `true` to re-enable.
const PHASE_RATIO_ENABLED: bool = false;
/// Node budget for a P2 result: the recurrence nests, so the emitted sum
/// grows like `2^deg P` and must be bounded independently of the input.
const MAX_RESULT_NODES: usize = 4000;
/// Cap on the T1 base re-entries performed by one P1/P2 reduction.
const MAX_SUB_CALLS: u32 = 2;

/// Combined entry: T1 (denominator power) → P1 (phase shift) → T2 (linear
/// numerator) → P2 (polynomial×kernel ratio) → T3 (polynomial×trig closed
/// form).
///
/// P1 sits after T1 so a plain `1/(a + b·T(u))^n` never pays for the
/// mixed-kernel scan, and before T2/P2 because the shifted denominator is a
/// *single* kernel, which is exactly what those phases consume.
pub(crate) fn integrate_trig_reduction<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    if !PHASE_RATIO_ENABLED {
        return integrate_trig_denom(ctx, expr, var)
            .or_else(|| integrate_trig_num_linear(ctx, expr, var))
            .or_else(|| integrate_poly_trig(ctx, expr, var));
    }
    integrate_trig_denom(ctx, expr, var)
        .or_else(|| integrate_phase_shift(ctx, expr, var))
        .or_else(|| integrate_trig_num_linear(ctx, expr, var))
        .or_else(|| integrate_poly_kernel_ratio(ctx, expr, var))
        .or_else(|| integrate_poly_trig(ctx, expr, var))
}

// =========================================================================
// Shared matcher: one `(a + b·T(u))^(−n)` factor in a product
// =========================================================================

/// A product split around a single trig-denominator power.
struct DenomMatch<'a> {
    /// Constant (w.r.t. `var`) leftover factors.
    rest: Vec<Atom<'a>>,
    /// Non-constant factors that are not the denominator power.
    other: Vec<Atom<'a>>,
    /// `a + b·T(u)` as found.
    base: Atom<'a>,
    a: Atom<'a>,
    b: Atom<'a>,
    sin: bool,
    u: Atom<'a>,
    /// `du/dx` (nonzero).
    du: Atom<'a>,
    /// The denominator power `n ≥ 1`.
    n: i64,
    /// `a² − b²`, verified nonzero.
    disc: Atom<'a>,
}

fn match_trig_denom<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<DenomMatch<'a>> {
    let factors: Vec<Atom<'a>> = match expr.node() {
        AtomNode::Mul(args) => args.to_vec(),
        _ => vec![expr],
    };
    let mut rest = Vec::new();
    let mut other = Vec::new();
    let mut found: Option<(Atom<'a>, Atom<'a>, Atom<'a>, bool, Atom<'a>, Atom<'a>, i64)> = None;
    for f in factors {
        if is_constant(f, var) {
            rest.push(f);
            continue;
        }
        if let Some((base, n)) = as_recip_pow(f)
            && let Some((a, b, sin, u, du)) = split_trig_base(ctx, base, var)
        {
            // A second trig-denominator power is out of scope (decline).
            if found.is_some() {
                return None;
            }
            found = Some((base, a, b, sin, u, du, n));
            continue;
        }
        other.push(f);
    }
    let (base, a, b, sin, u, du, n) = found?;
    // Fold the leftover factors into one atom: consumers (T2/P2) want the
    // whole non-denominator part as a single product, and a constant-only
    // leftover must read as empty.
    let other: Vec<Atom<'a>> = if other.is_empty() {
        Vec::new()
    } else {
        let folded = normalize(ctx, ctx.mul(&other));
        if is_constant(folded, var) {
            Vec::new()
        } else {
            vec![folded]
        }
    };
    let disc = discriminant(ctx, a, b)?;
    Some(DenomMatch {
        rest,
        other,
        base,
        a,
        b,
        sin,
        u,
        du,
        n,
        disc,
    })
}

/// `(base, n)` for `(base^n)^(−1)` or `base^(−n)` with integer `n ≥ 1` —
/// the two shapes the parser/normalizer produce for `1/base^n`. The nested
/// form is tried first: `(base^2)^(−1)` must read as `n = 2`, not as
/// `Pow(base^2, −1)` with `n = 1`.
fn as_recip_pow(f: Atom<'_>) -> Option<(Atom<'_>, i64)> {
    let AtomNode::Pow(b, e) = f.node() else {
        return None;
    };
    if let (AtomNode::Pow(bb, ee), AtomNode::Num(-1)) = (b.node(), e.node())
        && let AtomNode::Num(n) = ee.node()
        && *n >= 1
    {
        return Some((*bb, *n));
    }
    if let AtomNode::Num(n) = e.node()
        && *n <= -1
    {
        return Some((*b, n.checked_neg()?));
    }
    None
}

/// Decomposition of one additive term.
enum TermSplit<'a> {
    Constant,
    /// `coeff · T(u)` with `coeff` constant and `T ∈ {sin, cos}`.
    Trig {
        coeff: Atom<'a>,
        sin: bool,
        u: Atom<'a>,
    },
    Other,
}

/// Decompose one additive term: constant, or constants × a single bare
/// `sin(u)`/`cos(u)` factor.
fn split_trig_term<'a>(ctx: &'a AtomArena<'a>, term: Atom<'a>, var: Symbol) -> TermSplit<'a> {
    if is_constant(term, var) {
        return TermSplit::Constant;
    }
    let factors: Vec<Atom<'a>> = match term.node() {
        AtomNode::Mul(args) => args.to_vec(),
        // A bare `sin(u)`/`cos(u)`: the coefficient is the implicit unit.
        AtomNode::Fun(name, args)
            if args.len() == 1 && (name.as_str() == "sin" || name.as_str() == "cos") =>
        {
            return TermSplit::Trig {
                coeff: ctx.num(1),
                sin: name.as_str() == "sin",
                u: args[0],
            };
        }
        _ => vec![term],
    };
    let mut coeff: Vec<Atom<'a>> = Vec::new();
    let mut trig: Option<(bool, Atom<'a>)> = None;
    for f in factors {
        if is_constant(f, var) {
            coeff.push(f);
            continue;
        }
        match f.node() {
            AtomNode::Fun(name, args)
                if args.len() == 1 && (name.as_str() == "sin" || name.as_str() == "cos") =>
            {
                if trig.is_some() {
                    return TermSplit::Other;
                }
                trig = Some((name.as_str() == "sin", args[0]));
            }
            _ => return TermSplit::Other,
        }
    }
    let Some((sin, u)) = trig else {
        return TermSplit::Other;
    };
    let coeff = if coeff.is_empty() {
        ctx.num(1)
    } else {
        normalize(ctx, ctx.mul(&coeff))
    };
    TermSplit::Trig { coeff, sin, u }
}

/// One bare trig term of a base: `(coefficient, is_sin, argument)`.
///
type TrigTerm<'a> = (Atom<'a>, bool, Atom<'a>);

/// The base as a constant plus its bare trig terms: `(a, [(coeff, sin, u)])`.
/// The caller decides how many trig terms its class accepts.
fn split_base_trigs<'a>(
    ctx: &'a AtomArena<'a>,
    base: Atom<'a>,
    var: Symbol,
) -> Option<(Atom<'a>, Vec<TrigTerm<'a>>)> {
    let AtomNode::Add(args) = base.node() else {
        return None;
    };
    let mut consts: Vec<Atom<'a>> = Vec::new();
    let mut trigs: Vec<TrigTerm<'a>> = Vec::new();
    for t in args.iter() {
        match split_trig_term(ctx, *t, var) {
            TermSplit::Constant => consts.push(*t),
            TermSplit::Trig { coeff, sin, u } => trigs.push((coeff, sin, u)),
            TermSplit::Other => return None,
        }
    }
    let a = if consts.is_empty() {
        ctx.num(0)
    } else {
        normalize(ctx, ctx.add(&consts))
    };
    Some((a, trigs))
}

/// Match `base` as `a + b·T(u)`: constant terms plus exactly one trig term
/// whose argument is linear in `var` with nonzero slope. Returns
/// `(a, b, sin, u, du/dx)`. A two-kernel base is outside T1/T2's class.
fn split_trig_base<'a>(
    ctx: &'a AtomArena<'a>,
    base: Atom<'a>,
    var: Symbol,
) -> Option<(Atom<'a>, Atom<'a>, bool, Atom<'a>, Atom<'a>)> {
    let (a, trigs) = split_base_trigs(ctx, base, var)?;
    if trigs.len() != 1 {
        return None;
    }
    let (b, sin, u) = trigs[0];
    let (du, _phase) = linear_form(ctx, u, var)?;
    if matches!(du.node(), AtomNode::Num(0)) {
        return None;
    }
    Some((a, b, sin, u, du))
}

/// Match `base` as the two-kernel sum `a + b·cos(u) + c·sin(u)`: exactly two
/// trig terms, one of each kernel, sharing one linear argument. Returns
/// `(a, b, c, u, du/dx)` with `b` the `cos` coefficient and `c` the `sin`
/// one.
fn split_mixed_base<'a>(
    ctx: &'a AtomArena<'a>,
    base: Atom<'a>,
    var: Symbol,
) -> Option<(Atom<'a>, Atom<'a>, Atom<'a>, Atom<'a>, Atom<'a>)> {
    let (a, trigs) = split_base_trigs(ctx, base, var)?;
    if trigs.len() != 2 {
        return None;
    }
    let (mut b, mut c) = (None, None);
    let mut u_norm = None;
    for (coeff, sin, u) in trigs {
        if !is_constant(coeff, var) {
            return None;
        }
        let slot = if sin { &mut c } else { &mut b };
        if slot.is_some() {
            return None;
        }
        *slot = Some(coeff);
        let un = normalize(ctx, u);
        match u_norm {
            None => u_norm = Some(un),
            Some(prev) if prev == un => {}
            Some(_) => return None,
        }
    }
    let (b, c) = (b?, c?);
    let u = u_norm?;
    // Both coefficients vanishing leaves nothing to normalize
    // (`φ = atan(0/0)` is not a number).
    if matches!(b.node(), AtomNode::Num(0)) && matches!(c.node(), AtomNode::Num(0)) {
        return None;
    }
    let (du, _phase) = linear_form(ctx, u, var)?;
    if matches!(du.node(), AtomNode::Num(0)) {
        return None;
    }
    Some((a, b, c, u, du))
}

/// `a² − b²`, folded for numeric `a`/`b`; `None` when it is zero (exactly
/// for numeric values, syntactically via `collect_terms` for symbolic
/// ones) — the singular case the recurrence cannot divide by.
fn discriminant<'a>(ctx: &'a AtomArena<'a>, a: Atom<'a>, b: Atom<'a>) -> Option<Atom<'a>> {
    if let (Some(av), Some(bv)) = (numeric(a), numeric(b)) {
        let d = (av as i128)
            .checked_mul(av as i128)?
            .checked_sub((bv as i128).checked_mul(bv as i128)?)?;
        if d == 0 {
            return None;
        }
        return i64::try_from(d).ok().map(|v| ctx.num(v));
    }
    let sq_a = square(ctx, a)?;
    let sq_b = square(ctx, b)?;
    let sum = normalize(
        ctx,
        fold_sqrt_squares(ctx, ctx.add(&[sq_a, ctx.mul(&[ctx.num(-1), sq_b])])),
    );
    let d = crate::ode::util::collect_terms(ctx, sum);
    if matches!(d.node(), AtomNode::Num(0)) {
        return None;
    }
    Some(d)
}

/// Fold `sqrt(e)² → e` (and the equivalent `e^(1/2)` power) bottom-up, so a
/// discriminant built from a radical coefficient depends on the radicand
/// rather than on a nested radical. `(√5)² → 5` lets `9 − (√5)²` collapse to
/// `4` and a numeric half-power fold to an exact rational.
fn fold_sqrt_squares<'a>(ctx: &'a AtomArena<'a>, expr: Atom<'a>) -> Atom<'a> {
    match expr.node() {
        AtomNode::Num(_) | AtomNode::Var(_) => expr,
        AtomNode::Add(args) => {
            let kids: Vec<Atom<'a>> = args.iter().map(|a| fold_sqrt_squares(ctx, *a)).collect();
            normalize(ctx, ctx.add(&kids))
        }
        AtomNode::Mul(args) => {
            let kids: Vec<Atom<'a>> = args.iter().map(|a| fold_sqrt_squares(ctx, *a)).collect();
            normalize(ctx, ctx.mul(&kids))
        }
        AtomNode::Fun(_, _) => expr,
        AtomNode::Pow(b, e) => {
            let base = fold_sqrt_squares(ctx, *b);
            let exp = fold_sqrt_squares(ctx, *e);
            if let AtomNode::Num(k) = exp.node()
                && *k % 2 == 0
                && let AtomNode::Fun(name, fargs) = base.node()
                && name.as_str() == "sqrt"
                && fargs.len() == 1
                && let Some(folded) = int_pow(ctx, fargs[0], *k)
            {
                return folded;
            }
            normalize(ctx, ctx.pow(base, exp))
        }
    }
}

fn square<'a>(ctx: &'a AtomArena<'a>, e: Atom<'a>) -> Option<Atom<'a>> {
    match e.node() {
        AtomNode::Num(v) => Some(ctx.num(v.checked_mul(*v)?)),
        _ => Some(ctx.pow(e, ctx.num(2))),
    }
}

// =========================================================================
// T1: ∫ dx/(a + b·T(u))^n
// =========================================================================

/// T1 entry: `∫ C·dx/(a + b·T(u))^n` for `2 ≤ n ≤ 6` (see module docs).
pub(crate) fn integrate_trig_denom<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    if is_constant(expr, var) || node_count(expr) > MAX_NODES {
        return None;
    }
    let m = match_trig_denom(ctx, expr, var)?;
    if !(2..=MAX_DENOM_POW).contains(&m.n) || !m.other.is_empty() {
        return None;
    }
    let core = denom_integral(ctx, &m, m.n, var)?;
    let mut factors = m.rest.clone();
    factors.push(core);
    Some(normalize(ctx, ctx.mul(&factors)))
}

/// `J_n = ∫ dx/(a + b·T(u))^n` by the recurrence in the module docs.
///
/// The `n = 1` base first re-enters the chain (Weierstrass / symbolic
/// rational), which is what 0.27.1 shipped and what the T1/T2 families are
/// verified against. The chain declines whenever the kernel coefficient is a
/// radical — precisely the phase-shifted P1 base — so the module's own
/// Weierstrass closed form [`weierstrass_base`] takes over as the fallback.
fn denom_integral<'a>(
    ctx: &'a AtomArena<'a>,
    m: &DenomMatch<'a>,
    n: i64,
    var: Symbol,
) -> Option<Atom<'a>> {
    debug_assert!((1..=MAX_DENOM_POW).contains(&n));
    if n == 1 {
        // The chain first (it owns the `a² < b²` and negative-`a` shapes and
        // its answer is what 0.27.1 shipped), then the module's closed form.
        // The closed form is the *only* route when the kernel coefficient is
        // a radical, which is what the phase shift P1 feeds it.
        let g = ctx.pow(m.base, ctx.num(-1));
        let r = integrate_raw(ctx, g, var, 0, true, 0, 0);
        if contains_integral(r) || is_fallback(&r) {
            return weierstrass_base(ctx, m, var);
        }
        return Some(r);
    }
    let j1 = denom_integral(ctx, m, n - 1, var)?;
    let j2 = if n == 2 {
        None
    } else {
        Some(denom_integral(ctx, m, n - 2, var)?)
    };
    // Boundary term ε·b·S(u)/(du·D^(n−1)): ε = −1, S = sin for T = cos;
    // ε = +1, S = cos for T = sin.
    let (sign, other_sin) = if m.sin { (1, false) } else { (-1, true) };
    let boundary = ctx.mul(&[
        ctx.num(sign),
        m.b,
        ctx.fun(if other_sin { "sin" } else { "cos" }, &[m.u]),
        ctx.pow(m.base, ctx.num(-(n - 1))),
        ctx.pow(m.du, ctx.num(-1)),
    ]);
    let mut terms = vec![boundary, ctx.mul(&[m.a, ctx.num(2 * n - 3), j1])];
    if let Some(j2) = j2 {
        terms.push(ctx.mul(&[ctx.num(-(n - 2)), j2]));
    }
    let bracket = normalize(ctx, ctx.add(&terms));
    let denom = normalize(ctx, ctx.mul(&[ctx.num(n - 1), m.disc]));
    Some(normalize(
        ctx,
        ctx.mul(&[bracket, ctx.pow(denom, ctx.num(-1))]),
    ))
}

// =========================================================================
// Base case: the Weierstrass closed form of `∫ dx/(a + b·T(u))`
// =========================================================================

/// Closed form of `∫ dx/(a + b·T(u))`, the `n = 1` base of the T1
/// recurrence, by the Weierstrass half-angle substitution `t = tan(v/2)`
/// with `v = d·x + e` (`d` may be symbolic, and for P1 it is `u − φ`):
///
/// ```text
/// a² > b², T = cos:  2·atan( t·√((a−b)/(a+b)) ) / √(a²−b²)
/// a² > b², T = sin:  2·atan( (b + a·t)/√(a²−b²) ) / √(a²−b²)
/// ```
///
/// Both forms were checked by numerically differentiating the closed form
/// against the integrand (the `sin` argument `(a+b)·t` is **wrong** — it was
/// the first form tried and it fails that check; `(b + a·t)` passes). The
/// `atan` branch constant is absorbed into the integration constant, which
/// is what makes the antiderivative continuous on `|v| < π`. The caller has
/// already established `a² ≠ b²`.
///
/// The `a² < b²` family **declines**: its real form is a logarithm whose
/// bookkeeping did not survive the same numerical check, and a wrong closed
/// form is worse than a decline (the T1 `n = 1` base then re-enters the
/// chain, which owns those shapes).
fn weierstrass_base<'a>(
    ctx: &'a AtomArena<'a>,
    m: &DenomMatch<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    let (d, e) = linear_form(ctx, m.u, var)?;
    if matches!(d.node(), AtomNode::Num(0)) {
        return None;
    }
    // A negative discriminant means `a² < b²`: decline (see the docs).
    if matches!(m.disc.node(), AtomNode::Num(v) if *v < 0) {
        return None;
    }
    // The verified `cos` argument is the antiderivative only for `a > 0`:
    // numerically, `a = -5, b = 3` makes it the negative of the integrand,
    // and `sign(a)` repairs it exactly. A symbolic `a` has no known sign, so
    // that family declines rather than risk a wrong answer; `a = 0` is already
    // excluded by the discriminant.
    // The verified `sin` argument `(b + a·t)/√(a²−b²)` is sign-robust (checked
    // for both signs of `a` and `b`), so it serves a symbolic `a` too. The
    // `cos` argument `t·√((a−b)/(a+b))` is the antiderivative only for `a > 0`:
    // numerically `a = -5, b = 3` makes it the negative of the integrand, and
    // `sign(a)` repairs it exactly. A symbolic `a` therefore declines for `cos`
    // only, rather than risk a wrong sign.
    let sign_a = if m.sin {
        ctx.num(1)
    } else {
        match numeric(m.a) {
            Some(v) if v > 0 => ctx.num(1),
            Some(v) if v < 0 => ctx.num(-1),
            Some(_) => return None,
            None => return None,
        }
    };

    // `-1 * 2^(-1)` silently yields `disc^(-1)` instead of `disc^(-1/2)`.
    let half = rat_atom(ctx, 1, 2);
    let inv_half = rat_atom(ctx, -1, 2);
    let v = ctx.add(&[ctx.mul(&[d, ctx.var(var.as_str())]), e]);
    let t = ctx.fun("tan", &[normalize(ctx, ctx.mul(&[half, v]))]);
    let inv_root = ctx.pow(m.disc, inv_half);
    // Verified closed forms (see the doc comment): the `cos` argument is
    // `t·√((a−b)/(a+b))` with **no** `1/√(a²−b²)` factor, while the `sin`
    // argument is `(b + a·t)/√(a²−b²)`. Multiplying the `cos` argument by the
    // reciprocal root as well produced a wrong antiderivative.
    let atan_arg = if m.sin {
        ctx.mul(&[ctx.add(&[m.b, ctx.mul(&[m.a, t])]), inv_root])
    } else {
        let ratio = normalize(
            ctx,
            ctx.mul(&[
                ctx.add(&[m.a, ctx.mul(&[ctx.num(-1), m.b])]),
                inv(ctx, ctx.add(&[m.a, m.b])),
            ]),
        );
        ctx.mul(&[t, ctx.fun("sqrt", &[ratio])])
    };
    let atan = ctx.fun("atan", &[normalize(ctx, atan_arg)]);
    Some(normalize(
        ctx,
        ctx.mul(&[ctx.num(2), sign_a, inv_root, atan]),
    ))
}

// =========================================================================
// P1: phase-shift normalization of a two-kernel denominator
// =========================================================================
/// P1 entry: `a + b·cos(u) + c·sin(u)` → `a + R·cos(u − φ)` with
/// `R = √(b² + c²)` and `φ = atan(c/b)`, then delegate to T1/T2/P2 on the
/// shifted (still linear) argument. See the module docs for the branch
/// discussion.
pub(crate) fn integrate_phase_shift<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    if is_constant(expr, var) || node_count(expr) > MAX_NODES {
        return None;
    }
    let factors: Vec<Atom<'a>> = match expr.node() {
        AtomNode::Mul(args) => args.to_vec(),
        _ => vec![expr],
    };
    let mut rest: Vec<Atom<'a>> = Vec::new();
    let mut found: Option<(Atom<'a>, Atom<'a>, Atom<'a>, Atom<'a>)> = None;
    for f in &factors {
        if is_constant(*f, var) {
            rest.push(*f);
            continue;
        }
        let Some((base, _n)) = as_recip_pow(*f) else {
            continue;
        };
        let Some((a, b, c, u, _du)) = split_mixed_base(ctx, base, var) else {
            continue;
        };
        // Two mixed-kernel denominator powers are out of scope.
        if found.is_some() {
            return None;
        }
        found = Some((a, b, c, u));
    }
    let (a, b, c, u) = found?;
    let phi = ctx.fun("atan", &[normalize(ctx, ctx.mul(&[c, inv(ctx, b)]))]);
    // `R²` is developed as the algebraic `b² + c²` rather than as the square
    // of the radical `R`, so the shifted discriminant stays free of nested
    // radicals and a numeric `a² − R²` folds exactly.
    let r_sq = normalize(ctx, ctx.add(&[square(ctx, b)?, square(ctx, c)?]));
    let r = ctx.fun("sqrt", &[r_sq]);
    // The shifted denominator `a + R·cos(u − φ)` is a *single*-kernel base,
    // so it must survive the same discriminant test T1 applies.
    discriminant(ctx, a, r)?;
    for flip in [false, true] {
        let base = shifted_base(ctx, a, r, u, phi, flip);
        let shifted = rebuild(ctx, &factors, &rest, base, var);
        let Some(out) = integrate_trig_reduction(ctx, shifted, var) else {
            continue;
        };
        if contains_integral(out) {
            continue;
        }
        return Some(out);
    }
    None
}

/// `a + R·cos(u − φ)`; `flip` selects the branch `φ + π`, i.e. `R → −R`,
/// which has the same discriminant `a² − R²`.
fn shifted_base<'a>(
    ctx: &'a AtomArena<'a>,
    a: Atom<'a>,
    r: Atom<'a>,
    u: Atom<'a>,
    phi: Atom<'a>,
    flip: bool,
) -> Atom<'a> {
    let r_signed = if flip { ctx.mul(&[ctx.num(-1), r]) } else { r };
    let v = normalize(ctx, ctx.add(&[u, ctx.mul(&[ctx.num(-1), phi])]));
    let kernel = ctx.mul(&[r_signed, ctx.fun("cos", &[v])]);
    if matches!(a.node(), AtomNode::Num(0)) {
        normalize(ctx, kernel)
    } else {
        normalize(ctx, ctx.add(&[a, kernel]))
    }
}

/// Rebuild the integrand with the mixed-kernel denominator base replaced by
/// `base`, leaving every other factor (including the denominator's power)
/// untouched.
fn rebuild<'a>(
    ctx: &'a AtomArena<'a>,
    factors: &[Atom<'a>],
    rest: &[Atom<'a>],
    base: Atom<'a>,
    var: Symbol,
) -> Atom<'a> {
    let mut out: Vec<Atom<'a>> = rest.to_vec();
    for f in factors {
        if is_constant(*f, var) {
            continue;
        }
        if let Some((old_base, n)) = as_recip_pow(*f)
            && split_mixed_base(ctx, old_base, var).is_some()
        {
            out.push(ctx.pow(base, ctx.num(-n)));
            continue;
        }
        out.push(*f);
    }
    if out.is_empty() {
        return ctx.num(1);
    }
    normalize(ctx, ctx.mul(&out))
}

// =========================================================================
// P2: ∫ P(x)·T(u)/(a + b·T(u))^n dx
// =========================================================================

/// P2 entry: `∫ C·P(x)·T(u)/(a + b·T(u))^n dx` for
/// `deg P ≤ min(n − 1, 4)` (see the module docs for the reduction). The
/// numerator kernel and the denominator kernel must share one argument and
/// one function.
pub(crate) fn integrate_poly_kernel_ratio<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    if is_constant(expr, var) || node_count(expr) > MAX_NODES {
        return None;
    }
    let m = match_trig_denom(ctx, expr, var)?;
    if !(1..=MAX_DENOM_POW).contains(&m.n) {
        return None;
    }
    let (coef, deg) = split_numerator_monomial(ctx, *m.other.first()?, var, m.sin, m.u)?;
    if deg == 0 || deg > MAX_RATIO_DEG || deg >= m.n {
        return None;
    }
    let mut calls = 0u32;
    let core = reduce_ratio(ctx, &m, coef, deg, var, &mut calls)?;
    let result = multiply_all(ctx, &m.rest, core);
    if contains_integral(result) || node_count(result) > MAX_RESULT_NODES {
        return None;
    }
    Some(result)
}

/// `factors[0]·…·factors[n]·extra`, tolerating an empty factor list (the
/// `AtomArena::mul` builder rejects an empty slice).
fn multiply_all<'a>(ctx: &'a AtomArena<'a>, factors: &[Atom<'a>], extra: Atom<'a>) -> Atom<'a> {
    match factors {
        [] => extra,
        [f] => normalize(ctx, ctx.mul(&[*f, extra])),
        _ => {
            let mut all: Vec<Atom<'a>> = factors.to_vec();
            all.push(extra);
            normalize(ctx, ctx.mul(&all))
        }
    }
}

/// `∫ P(x)·T(u)·D^(−n) dx` by the recurrence in the module docs; `P` is the
/// monomial `coef·x^deg` and `deg ≤ n − 1`, so the loop terminates at a
/// constant residual.
fn reduce_ratio<'a>(
    ctx: &'a AtomArena<'a>,
    m: &DenomMatch<'a>,
    coef: Atom<'a>,
    deg: i64,
    var: Symbol,
    calls: &mut u32,
) -> Option<Atom<'a>> {
    debug_assert!(deg >= 0 && deg < m.n);
    let mut p = coef;
    let mut power = m.n;
    let mut out: Vec<Atom<'a>> = Vec::new();
    // `deg ≤ n − 1`, so at most `deg` steps of the recurrence are taken and
    // the walk stops at the constant residual `∫ T·D^(−power)`: a closed
    // form for `power ≥ 2`, T1's `n = 1` base for `power = 1`, and the
    // `deg = n` case (rejected upstream) for `power = 0`.
    loop {
        debug_assert!(power >= 1);
        out.push(ratio_boundary(ctx, m, p, power)?);
        if power == 1 {
            *calls += 1;
            if *calls > MAX_SUB_CALLS {
                return None;
            }
            let j1 = denom_integral(ctx, m, 1, var)?;
            if contains_integral(j1) || is_fallback(&j1) {
                return None;
            }
            out.push(ctx.mul(&[ratio_const(ctx, m)?, p, j1]));
            return Some(normalize(ctx, ctx.add(&out)));
        }
        p = ratio_next(ctx, m, &p, power, var)?;
        if !is_polynomial_in_var(p, var) {
            return None;
        }
        power -= 1;
    }
}

/// `C = −ε/b` with `ε = +1` for `T = sin` and `−1` for `T = cos`.
fn ratio_const<'a>(ctx: &'a AtomArena<'a>, m: &DenomMatch<'a>) -> Option<Atom<'a>> {
    let eps = if m.sin { 1 } else { -1 };
    let signed = ctx.mul(&[ctx.num(eps), inv(ctx, m.b)]);
    Some(normalize(ctx, ctx.mul(&[ctx.num(-1), signed])))
}

/// The boundary term `C·P·D^(1−n)` of one integration-by-parts step.
fn ratio_boundary<'a>(
    ctx: &'a AtomArena<'a>,
    m: &DenomMatch<'a>,
    p: Atom<'a>,
    n: i64,
) -> Option<Atom<'a>> {
    let c = ratio_const(ctx, m)?;
    Some(normalize(
        ctx,
        ctx.mul(&[c, p, ctx.pow(m.base, ctx.num(1 - n))]),
    ))
}

/// The polynomial of the next residual:
/// `R = (a/b)·P + P'/b`, i.e.
/// `K(P, n) = C·P·D^(1−n) + K(R, n−1)` with `C = 1/(b·d)`. The boundary
/// term does not depend on `n`, so the `P'` coefficient carries no `1/(n−1)`
/// factor.
fn ratio_next<'a>(
    ctx: &'a AtomArena<'a>,
    m: &DenomMatch<'a>,
    p: &Atom<'a>,
    _n: i64,
    var: Symbol,
) -> Option<Atom<'a>> {
    let shift = normalize(ctx, ctx.mul(&[m.a, inv(ctx, m.b), *p]));
    let derivative = crate::diff(ctx, *p, var);
    let slope = normalize(ctx, ctx.mul(&[inv(ctx, m.b), inv(ctx, m.du), derivative]));
    Some(normalize(ctx, ctx.add(&[shift, slope])))
}

/// Split a P2 numerator as `P(x)·T(u)`: returns the polynomial's
/// `(coeff, degree)` with `T` the same kernel and argument as the
/// denominator's. A bare kernel (constant `P`) and any non-polynomial
/// leftover (a second kernel, a half power, …) decline.
fn split_numerator_monomial<'a>(
    ctx: &'a AtomArena<'a>,
    num: Atom<'a>,
    var: Symbol,
    sin: bool,
    u_den: Atom<'a>,
) -> Option<(Atom<'a>, i64)> {
    let factors: Vec<Atom<'a>> = match num.node() {
        AtomNode::Mul(args) => args.to_vec(),
        _ => vec![num],
    };
    let mut coef: Vec<Atom<'a>> = Vec::new();
    let mut deg: Option<i64> = None;
    let mut kernel_seen = false;
    for f in factors {
        if let AtomNode::Fun(name, args) = f.node()
            && args.len() == 1
            && (name.as_str() == "sin" || name.as_str() == "cos")
        {
            let is_sin = name.as_str() == "sin";
            let arg_eq = normalize(ctx, args[0]) == normalize(ctx, u_den);
            eprintln!();
            if kernel_seen || is_sin != sin || !arg_eq {
                return None;
            }
            kernel_seen = true;
            continue;
        }
        if deg.is_none() {
            match f.node() {
                AtomNode::Var(v) if *v == var => {
                    deg = Some(1);
                    continue;
                }
                AtomNode::Pow(b, e)
                    if matches!(b.node(), AtomNode::Var(v) if *v == var)
                        && matches!(e.node(), AtomNode::Num(k) if *k >= 1) =>
                {
                    let AtomNode::Num(k) = e.node() else {
                        return None;
                    };
                    deg = Some(*k);
                    continue;
                }
                _ => {}
            }
        }
        if !is_constant(f, var) {
            return None;
        }
        coef.push(f);
    }
    if !kernel_seen {
        return None;
    }
    let deg = deg?;
    let c = if coef.is_empty() {
        ctx.num(1)
    } else {
        normalize(ctx, ctx.mul(&coef))
    };
    Some((c, deg))
}

/// True when `expr` is a polynomial in `var` of degree `≤ MAX_RATIO_DEG`
/// (a constant, or a sum of monomials with constant coefficients). Every
/// intermediate residual of the P2 recurrence has this shape by
/// construction; the check is a cheap guard on the arithmetic.
fn is_polynomial_in_var(expr: Atom<'_>, var: Symbol) -> bool {
    let terms: Vec<Atom<'_>> = match expr.node() {
        AtomNode::Add(args) => args.to_vec(),
        _ => vec![expr],
    };
    terms.iter().all(|t| {
        let factors: Vec<Atom<'_>> = match t.node() {
            AtomNode::Mul(args) => args.to_vec(),
            _ => vec![*t],
        };
        let mut deg = 0i64;
        for f in factors {
            match f.node() {
                AtomNode::Var(v) if *v == var => deg += 1,
                AtomNode::Pow(b, e) => {
                    if matches!(b.node(), AtomNode::Var(v) if *v == var)
                        && let AtomNode::Num(k) = e.node()
                    {
                        deg += *k;
                    } else if !is_constant(f, var) {
                        return false;
                    }
                }
                _ => {
                    if !is_constant(f, var) {
                        return false;
                    }
                }
            }
        }
        deg <= MAX_RATIO_DEG
    })
}

// =========================================================================
// T2: (A + B·T(u))/(a + b·T(u))^n
// =========================================================================

/// T2 entry: `(A + B·T(u))/(a + b·T(u))^n` for `2 ≤ n ≤ 6`, with the
/// numerator's trig function and argument matching the denominator's.
pub(crate) fn integrate_trig_num_linear<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    if is_constant(expr, var) || node_count(expr) > MAX_NODES {
        return None;
    }
    let m = match_trig_denom(ctx, expr, var)?;
    if !(2..=MAX_DENOM_POW).contains(&m.n) || m.other.len() != 1 {
        return None;
    }
    let u_norm = normalize(ctx, m.u);
    let (num_a, num_b) = split_numerator(ctx, m.other[0], var, m.sin, u_norm)?;
    if matches!(num_b.node(), AtomNode::Num(0)) {
        return None; // constant numerator: T1's shape
    }
    let jn = denom_integral(ctx, &m, m.n, var)?;
    let jn1 = denom_integral(ctx, &m, m.n - 1, var)?;
    // (A + B·T)/D^n = A·J_n + (B/b)·(J_{n−1} − a·J_n), from T = (D − a)/b.
    let term_a = ctx.mul(&[num_a, jn]);
    let inner = normalize(ctx, ctx.add(&[jn1, ctx.mul(&[ctx.num(-1), m.a, jn])]));
    let term_b = ctx.mul(&[num_b, inv(ctx, m.b), inner]);
    let core = normalize(ctx, ctx.add(&[term_a, term_b]));
    let mut factors = m.rest.clone();
    factors.push(core);
    Some(normalize(ctx, ctx.mul(&factors)))
}

/// Split the numerator as `A + B·T(u)`: constant terms vs. coefficients of
/// trig factors matching (`sin`, `u_norm`). Every term must classify.
fn split_numerator<'a>(
    ctx: &'a AtomArena<'a>,
    num: Atom<'a>,
    var: Symbol,
    sin: bool,
    u_norm: Atom<'a>,
) -> Option<(Atom<'a>, Atom<'a>)> {
    let terms: Vec<Atom<'a>> = match num.node() {
        AtomNode::Add(args) => args.to_vec(),
        _ => vec![num],
    };
    let mut a_terms = Vec::new();
    let mut b_terms = Vec::new();
    for t in terms {
        match split_trig_term(ctx, t, var) {
            TermSplit::Constant => a_terms.push(t),
            TermSplit::Trig { coeff, sin: s2, u } => {
                if s2 != sin || normalize(ctx, u) != u_norm {
                    return None;
                }
                b_terms.push(coeff);
            }
            TermSplit::Other => return None,
        }
    }
    let a = if a_terms.is_empty() {
        ctx.num(0)
    } else {
        normalize(ctx, ctx.add(&a_terms))
    };
    let b = if b_terms.is_empty() {
        ctx.num(0)
    } else {
        normalize(ctx, ctx.add(&b_terms))
    };
    Some((a, b))
}

// =========================================================================
// T3: polynomial × trig closed form
// =========================================================================

/// T3 entry: `∫ C·x^m·T(a·x + b) dx` (`1 ≤ m ≤ 6`) and
/// `∫ C·(c + d·x)^k·T(a·x + b) dx` (`1 ≤ k ≤ 4`) as direct closed forms.
pub(crate) fn integrate_poly_trig<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    if is_constant(expr, var) || node_count(expr) > MAX_NODES {
        return None;
    }
    let factors: Vec<Atom<'a>> = match expr.node() {
        AtomNode::Mul(args) => args.to_vec(),
        _ => vec![expr],
    };
    let mut rest = Vec::new();
    let mut poly: Option<Vec<(Atom<'a>, i64)>> = None;
    let mut trig: Option<(bool, Atom<'a>, Atom<'a>)> = None;
    for f in factors {
        if is_constant(f, var) {
            rest.push(f);
            continue;
        }
        match f.node() {
            AtomNode::Fun(name, args)
                if args.len() == 1 && (name.as_str() == "sin" || name.as_str() == "cos") =>
            {
                if trig.is_some() {
                    return None;
                }
                let (slope, _phase) = linear_form(ctx, args[0], var)?;
                if matches!(slope.node(), AtomNode::Num(0)) {
                    return None;
                }
                trig = Some((name.as_str() == "sin", args[0], slope));
            }
            AtomNode::Var(v) if *v == var => {
                if poly.is_some() {
                    return None;
                }
                poly = Some(vec![(ctx.num(1), 1)]);
            }
            AtomNode::Pow(b, e) => {
                if poly.is_some() {
                    return None;
                }
                let AtomNode::Num(k) = e.node() else {
                    return None;
                };
                if matches!(b.node(), AtomNode::Var(v) if *v == var) {
                    if !(1..=MAX_POLY_DEG).contains(k) {
                        return None;
                    }
                    poly = Some(vec![(ctx.num(1), *k)]);
                } else {
                    if !(1..=MAX_LIN_POW).contains(k) {
                        return None;
                    }
                    poly = Some(expand_linear_pow(ctx, *b, *k, var)?);
                }
            }
            _ => return None,
        }
    }
    let poly = poly?;
    let (sin, u, slope) = trig?;
    let mut terms = Vec::new();
    for (coef, deg) in poly {
        let cf = closed_form(ctx, deg, sin, u, slope, var)?;
        match cf.node() {
            AtomNode::Add(args) => {
                for t in args.iter() {
                    terms.push(ctx.mul(&[coef, *t]));
                }
            }
            _ => terms.push(ctx.mul(&[coef, cf])),
        }
    }
    let core = normalize(ctx, ctx.add(&terms));
    let mut factors = rest;
    factors.push(core);
    Some(normalize(ctx, ctx.mul(&factors)))
}

/// `(c + d·x)^k` as `(coeff, degree)` monomial pairs via the binomial
/// theorem (zero coefficients dropped).
fn expand_linear_pow<'a>(
    ctx: &'a AtomArena<'a>,
    base: Atom<'a>,
    k: i64,
    var: Symbol,
) -> Option<Vec<(Atom<'a>, i64)>> {
    let (d, c) = linear_form(ctx, base, var)?;
    if matches!(d.node(), AtomNode::Num(0)) {
        return None;
    }
    let mut out = Vec::new();
    for j in 0..=k {
        let coef = normalize(
            ctx,
            ctx.mul(&[
                ctx.num(binom(k, j)),
                int_pow(ctx, c, k - j)?,
                int_pow(ctx, d, j)?,
            ]),
        );
        if matches!(coef.node(), AtomNode::Num(0)) {
            continue;
        }
        out.push((coef, j));
    }
    if out.is_empty() { None } else { Some(out) }
}

/// Closed form of `∫ x^m·T(u) dx` (`T` = sin when `sin`, else cos) with
/// `u = slope·x + phase`, by m-fold tabular reduction (see module docs).
/// Emits the full sum directly — no chain re-entry.
fn closed_form<'a>(
    ctx: &'a AtomArena<'a>,
    m: i64,
    sin: bool,
    u: Atom<'a>,
    slope: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    debug_assert!((0..=MAX_POLY_DEG).contains(&m));
    let x = ctx.var(var.as_str());
    let mut terms = Vec::new();
    let mut fall: i64 = 1; // falling factorial m!/(m−k)!
    for k in 0..=m {
        if k > 0 {
            fall = fall.checked_mul(m - k + 1)?;
        }
        let out_sin = if sin { k % 2 == 1 } else { k % 2 == 0 };
        let neg = if sin {
            matches!(k % 4, 0 | 3)
        } else {
            k % 4 >= 2
        };
        let acoef = match slope.node() {
            AtomNode::Num(av) => {
                let p = av.checked_pow((k + 1) as u32)?;
                ctx.pow(ctx.num(p), ctx.num(-1))
            }
            _ => ctx.pow(slope, ctx.num(-(k + 1))),
        };
        let mut fs: Vec<Atom<'a>> = Vec::new();
        if neg {
            fs.push(ctx.num(-1));
        }
        if fall != 1 {
            fs.push(ctx.num(fall));
        }
        let deg = m - k;
        if deg > 0 {
            fs.push(int_pow(ctx, x, deg)?);
        }
        fs.push(acoef);
        fs.push(ctx.fun(if out_sin { "sin" } else { "cos" }, &[u]));
        terms.push(normalize(ctx, ctx.mul(&fs)));
    }
    Some(normalize(ctx, ctx.add(&terms)))
}

/// `C(n, k)` for small non-negative integers (0 outside the range).
fn binom(n: i64, mut k: i64) -> i64 {
    if k < 0 || k > n {
        return 0;
    }
    if k > n - k {
        k = n - k;
    }
    let mut r: i64 = 1;
    for i in 0..k {
        r = r.saturating_mul(n - i) / (i + 1);
    }
    r
}

// =========================================================================
// Helpers
// =========================================================================

/// Integer value of a `Num` atom.
fn numeric(e: Atom<'_>) -> Option<i64> {
    if let AtomNode::Num(n) = e.node() {
        Some(*n)
    } else {
        None
    }
}

/// `a^(−1)`, folding the trivial units.
fn inv<'a>(ctx: &'a AtomArena<'a>, a: Atom<'a>) -> Atom<'a> {
    match a.node() {
        AtomNode::Num(1) => ctx.num(1),
        AtomNode::Num(-1) => ctx.num(-1),
        _ => ctx.pow(a, ctx.num(-1)),
    }
}

/// `base^e` for integer `e ≥ 0`, folding numeric bases and `e ∈ {0, 1}`.
/// Negative exponents (never produced here) fall back to a `Pow` node.
fn int_pow<'a>(ctx: &'a AtomArena<'a>, base: Atom<'a>, e: i64) -> Option<Atom<'a>> {
    if e < 0 {
        return Some(ctx.pow(base, ctx.num(e)));
    }
    match e {
        0 => Some(ctx.num(1)),
        1 => Some(base),
        _ => {
            if let AtomNode::Num(v) = base.node() {
                Some(ctx.num(v.checked_pow(e as u32)?))
            } else {
                Some(ctx.pow(base, ctx.num(e)))
            }
        }
    }
}

// =========================================================================
// Tests
// =========================================================================

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

    /// Parse and normalize (mirroring what the chain feeds the mechanism).
    fn parse_norm<'a>(ctx: &'a AtomArena<'a>, s: &str) -> Atom<'a> {
        normalize(ctx, ocas_parse::parse(ctx, s).unwrap())
    }

    /// Numeric f64 evaluator for test verification (handles the operators
    /// the reduction results produce, incl. the Weierstrass atan/log base).
    fn eval_f64(expr: Atom<'_>, env: &[(Symbol, f64)]) -> Option<f64> {
        match expr.node() {
            AtomNode::Num(n) => Some(*n as f64),
            AtomNode::Var(v) => env.iter().find(|(s, _)| s == v).map(|(_, val)| *val),
            AtomNode::Add(args) => args
                .iter()
                .try_fold(0.0, |acc, a| Some(acc + eval_f64(*a, env)?)),
            AtomNode::Mul(args) => args
                .iter()
                .try_fold(1.0, |acc, a| Some(acc * eval_f64(*a, env)?)),
            AtomNode::Pow(b, e) => Some(eval_f64(*b, env)?.powf(eval_f64(*e, env)?)),
            AtomNode::Fun(name, args) => {
                let v = eval_f64(*args.first()?, env)?;
                Some(match name.as_str() {
                    "sin" => v.sin(),
                    "cos" => v.cos(),
                    "tan" => v.tan(),
                    "sec" => v.cos().recip(),
                    "csc" => v.sin().recip(),
                    "cot" => v.tan().recip(),
                    "exp" => v.exp(),
                    "log" => v.ln(),
                    "sqrt" => v.sqrt(),
                    "atan" => v.atan(),
                    _ => return None,
                })
            }
        }
    }

    /// Run the mechanism, require `Some` without residue, and check
    /// `diff(result) == integrand` numerically at the given sample points.
    /// Sample points must keep `a + b·T(u)` away from zero and `u/2` away
    /// from the `tan` poles (the Weierstrass base answer is a sawtooth).
    fn assert_antiderivative_num<'a>(
        ctx: &'a AtomArena<'a>,
        integrand: Atom<'a>,
        var: Symbol,
        consts: &[(Symbol, f64)],
        samples: &[f64],
    ) {
        let result = integrate_trig_reduction(ctx, integrand, var).expect("mechanism declined");
        assert!(
            !result.to_string().contains("Integral"),
            "residue: {result}"
        );
        let d = crate::diff(ctx, result, var);
        for &xv in samples {
            let mut env = consts.to_vec();
            env.push((var, xv));
            let lhs = eval_f64(d, &env).expect("eval diff");
            let rhs = eval_f64(integrand, &env).expect("eval integrand");
            let tol = 1e-6 * rhs.abs().max(1.0);
            assert!(
                (lhs - rhs).abs() < tol,
                "at x={xv}: diff={lhs} integrand={rhs} (result: {result})"
            );
        }
    }

    fn assert_declined<'a>(ctx: &'a AtomArena<'a>, input: &str) {
        let expr = parse_norm(ctx, input);
        assert!(
            integrate_trig_reduction(ctx, expr, Symbol::new("x")).is_none(),
            "expected None for {input}"
        );
    }

    // ------------------------- T1 -------------------------

    #[test]
    fn t1_cos_n2_numeric() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(2 + cos(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t1_cos_n3_corpus() {
        // Corpus sample: symbolic linear argument, numeric a/b.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(-5 + 3*cos(c + d*x))^3");
        let env = [(Symbol::new("c"), 0.4), (Symbol::new("d"), 1.3)];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.2, 0.5, 0.9]);
    }

    #[test]
    fn t1_sin_n2_numeric() {
        // a = 2, b = 1: the Weierstrass base 1/(2+sin x) is chain-correct.
        // (a ∤ 2b bases like 1/(3+2 sin x) hit the rational.rs
        // sqrt_positive_rational missing-/q bug — main-line fix.)
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(2 + sin(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t1_cos_n4_symbolic() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(a + b*cos(c + d*x))^4");
        let env = [
            (Symbol::new("a"), 2.0),
            (Symbol::new("b"), 0.5),
            (Symbol::new("c"), 0.3),
            (Symbol::new("d"), 0.8),
        ];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.2, 0.6, 1.0]);
    }

    #[test]
    fn t1_singular_declines() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // a² = b² exactly (numeric) and syntactically (symbolic a = b).
        assert_declined(&ctx, "1/(2 + 2*cos(x))^2");
        assert_declined(&ctx, "1/(a + a*sin(c + d*x))^2");
    }

    #[test]
    fn t1_out_of_scope_declines() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // n = 1 is the chain's (Weierstrass); n = 7 exceeds the budget;
        // a Mul base (no additive constant term) is not the T1 shape.
        assert_declined(&ctx, "1/(2 + cos(x))");
        assert_declined(&ctx, "1/(2 + cos(x))^7");
        assert_declined(&ctx, "1/(2*cos(x))^2");
    }

    // ------------------------- T2 -------------------------

    #[test]
    fn t2_cos_linear_num() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "(1 + 3*cos(x))/(2 + cos(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t2_sin_symbolic() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "(A + B*cos(c + d*x))/(a + b*cos(c + d*x))^3");
        let env = [
            (Symbol::new("A"), 1.5),
            (Symbol::new("B"), -0.5),
            (Symbol::new("a"), 2.0),
            (Symbol::new("b"), 0.5),
            (Symbol::new("c"), 0.3),
            (Symbol::new("d"), 1.1),
        ];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.2, 0.5, 0.8]);
    }

    #[test]
    fn t2_bare_trig_numerator() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "sin(x)/(2 + sin(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t2_two_denominator_powers_decline() {
        // Corpus sample: a product of two denominator powers is out of
        // scope (and a = b in the first factor is singular anyway).
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        assert_declined(
            &ctx,
            "(A + B*sin(e + f*x))/((a + a*sin(e + f*x))^3*(c - c*sin(e + f*x))^2)",
        );
    }

    #[test]
    fn t2_mismatched_trig_declines() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // sin numerator over a cos denominator is a u-substitution shape,
        // not this decomposition.
        assert_declined(&ctx, "sin(x)/(2 + cos(x))^2");
    }

    // ------------------------- P1: phase shift -------------------------

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p1_numeric_two_kernels_n2() {
        // 1/(3 + 2·cos x + sin x)² → 1/(3 + √5·cos(x − atan(1/2)))².
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(3 + 2*cos(x) + sin(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p1_symbolic_coeff_n2() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(a + 2*cos(x) + 3*sin(x))^2");
        let env = [(Symbol::new("a"), 4.0)];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.3, 0.7, 1.1]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p1_fully_symbolic_n1() {
        // The n = 1 base: R = √(b²+c²), φ = atan(c/b), a symbolic.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(a + b*cos(x) + c*sin(x))");
        let env = [
            (Symbol::new("a"), 3.0),
            (Symbol::new("b"), 1.5),
            (Symbol::new("c"), 0.7),
        ];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.3, 0.7, 1.1]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p1_two_kernels_n3() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(5 + 2*cos(x) + 2*sin(x))^3");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p1_negative_and_swapped_kernels() {
        // b < 0 (absorbed into cos φ < 0), kernels in the other order and a
        // linear argument.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(4 + 3*sin(1 + 2*x) - 4*cos(1 + 2*x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.1, 0.4, 0.8]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p1_bare_kernels_no_constant() {
        // a = 0: the denominator is a pure phase-shifted cosine.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(cos(x) + sin(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn p1_singular_declines() {
        // a² = b² + c²: the shifted denominator vanishes somewhere.
        assert_declined_ctx("1/(sqrt(2)*cos(x) + sqrt(2)*sin(x))^2");
        assert_declined_ctx("1/(2*cos(x))^2");
    }

    // ------------------------- P2: P(x)·T(u)/D^n -------------------------

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p2_x2_sin_over_sin_cubed() {
        // The numerator kernel must equal the denominator kernel: with
        // `D = a + b·sin u`, the reduction's complement is `cos u`.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "x^2*sin(x)/(a + b*sin(x))^3");
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 0.7)];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.3, 0.7, 1.1]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p2_x2_cos_over_cos_cubed() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "x^2*cos(x)/(a + b*cos(x))^3");
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 0.7)];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.3, 0.7, 1.1]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p2_linear_poly_corpus_shape() {
        // The corpus shape, symbolic throughout.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "(e + f*x)^2*cos(c + d*x)/(a + b*cos(c + d*x))^3");
        let env = [
            (Symbol::new("e"), 1.2),
            (Symbol::new("f"), 0.5),
            (Symbol::new("c"), 0.4),
            (Symbol::new("d"), 1.3),
            (Symbol::new("a"), 2.0),
            (Symbol::new("b"), 0.7),
        ];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.2, 0.5, 0.9]);
    }

    #[test]
    // P1/P2 are gated off (see PHASE_RATIO_ENABLED): the reduction is
    // derived and numerically checked, but the emitted atom is still wrong.
    #[ignore = "P1/P2 atom construction still being debugged"]
    fn p2_linear_over_quadratic() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "x*cos(x)/(3 + 2*cos(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn p2_out_of_scope_declines() {
        // deg P = n: the terminal residual would still carry an `x`.
        assert_declined_ctx("x^2*sin(x)/(2 + sin(x))^2");
        // Degree beyond the P2 budget.
        assert_declined_ctx("x^5*sin(x)/(2 + sin(x))^6");
        // Kernel mismatch: the numerator kernel differs from the
        // denominator's.
        assert_declined_ctx("x*sin(x)/(2 + cos(x))^3");
        // Half-integer denominator power.
        assert_declined_ctx("x*cos(x)/(2 + cos(x))^(3/2)");
    }

    // ------------------------- declines -------------------------

    #[test]
    fn declines_owned_by_sibling_agents() {
        // Hyperbolic twin (sibling agent).
        assert_declined_ctx("1/(a + b*tanh(x))");
        // Half-power / elliptic routing (sibling agent).
        assert_declined_ctx("cos(x)^(7/2)/(a + b*cos(x)^2)^(3/2)");
        // Exponential factor: not a trig-rational shape at all.
        assert_declined_ctx("exp(x)*cos(x)");
        // Hyperbolic two-kernel denominator.
        assert_declined_ctx("1/(a + b*cosh(x) + c*sinh(x))^2");
    }

    fn assert_declined_ctx(input: &str) {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        assert_declined(&ctx, input);
    }

    // ------------------------- stress -------------------------

    #[test]
    fn stress_repeat_is_stable() {
        // The `n = 1` base prefers the chain, whose availability depends on
        // the *global* chain budget, so a deep reduction's printed form can
        // legitimately differ between rounds. What must hold is that every
        // round either declines or emits a *correct* antiderivative, that the
        // outcome kind never flips, and that no round grows the expression.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let var = Symbol::new("x");
        let env = [
            (Symbol::new("a"), 2.0),
            (Symbol::new("b"), 0.7),
            (Symbol::new("A"), 1.5),
            (Symbol::new("B"), -0.5),
        ];
        let cases = [
            "1/(3 + 2*cos(x) + sin(x))^2",
            "1/(5 + 2*cos(x) + 2*sin(x))^3",
            "1/(a + b*cos(x) + c*sin(x))",
            "x^2*cos(x)/(a + b*sin(x))^3",
            "(1 + 2*x)^2*sin(x)/(3 + cos(x))^3",
            "1/(2 + cos(x))^2",
            "1/(3 + 2*sin(x))^2",
            "(A + B*cos(x))/(2 + cos(x))^3",
            "x^2*cos(x)",
            "sin(x)/(2 + cos(x))^2",
            "1/(a + b*tanh(x))",
            "cos(x)^(7/2)/(a + b*cos(x)^2)^(3/2)",
        ];
        let exprs: Vec<Atom<'_>> = cases.iter().map(|s| parse_norm(&ctx, s)).collect();
        let mut kinds: Vec<bool> = Vec::with_capacity(exprs.len());
        for (i, e) in exprs.iter().enumerate() {
            let r = integrate_trig_reduction(&ctx, *e, var);
            let nodes = r.map(node_count).unwrap_or(0);
            assert!(
                nodes <= MAX_RESULT_NODES,
                "case {i} produced {nodes} nodes (> {MAX_RESULT_NODES})"
            );
            kinds.push(r.is_some());
        }
        for round in 0..300 {
            for (i, e) in exprs.iter().enumerate() {
                let r = integrate_trig_reduction(&ctx, *e, var);
                assert_eq!(
                    r.is_some(),
                    kinds[i],
                    "case {i} flipped solved/declined at round {round}"
                );
                if let Some(a) = r {
                    let nodes = node_count(a);
                    assert!(nodes <= MAX_RESULT_NODES, "case {i} grew to {nodes} nodes");
                    // Every accepted answer must differentiate back.
                    let d = crate::diff(&ctx, a, var);
                    for &xv in &[0.37, 0.91] {
                        let mut sample = env.to_vec();
                        sample.push((var, xv));
                        if let (Some(lhs), Some(rhs)) =
                            (eval_f64(d, &sample), eval_f64(*e, &sample))
                        {
                            assert!(
                                (lhs - rhs).abs() < 1e-5 * rhs.abs().max(1.0),
                                "case {i} round {round} wrong at x={xv}: {lhs} vs {rhs}"
                            );
                        }
                    }
                }
            }
        }
    }

    // ------------------------- T3 -------------------------

    #[test]
    fn t3_x2_cos() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "x^2*cos(x)");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t3_x4_sin_symbolic_slope() {
        // Beyond the parts budget (PARTS_MAX_DEPTH = 2): the 0.27.1 gap.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "x^4*sin(a + b*x)");
        let env = [(Symbol::new("a"), 0.4), (Symbol::new("b"), 1.3)];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t3_x6_cos_numeric() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "x^6*cos(2*x - 1)");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t3_linpow4_cos() {
        // The post-trig_reduce corpus term shape: `(c+d*x)^k·T(a+b*x)`.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "(c + d*x)^4*cos(a + b*x)");
        let env = [
            (Symbol::new("a"), 1.1),
            (Symbol::new("b"), 0.6),
            (Symbol::new("c"), 0.7),
            (Symbol::new("d"), -0.4),
        ];
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &env, &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t3_linpow3_sin_numeric() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "(1 + 2*x)^3*sin(3*x)");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t3_out_of_scope_declines() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // Two trig factors: trig_reduce's product-to-sum owns this shape.
        assert_declined(&ctx, "x*sin(x)*cos(x)");
        // A trig power: rules C12/C14 and trig_reduce own this shape.
        assert_declined(&ctx, "x*sin(x)^2");
        // Degree beyond the budget.
        assert_declined(&ctx, "x^7*sin(x)");
    }
    #[test]
    fn t1_closed_form_sin_arm() {
        // `a² − b² = 5` is not a perfect square, so the chain's Weierstrass
        // route declines and the module's `sin` arm answers. This is the arm
        // that was wrong (`(a+b)·t` instead of `b + a·t`) and produced corpus
        // wrong answers.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(3 + 2*sin(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t1_closed_form_cos_arm() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(3 + 2*cos(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t1_closed_form_declines_where_unverified() {
        // `a² < b²`: the real branch of the closed form did not survive the
        // numeric check, so the base declines rather than emit a wrong
        // logarithm.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let var = Symbol::new("x");
        let e = parse_norm(&ctx, "1/(2 + 3*sin(x))^2");
        let m = match_trig_denom(&ctx, e, var).unwrap();
        assert!(weierstrass_base(&ctx, &m, var).is_none());
    }

    #[test]
    fn t1_closed_form_negative_a_uses_the_sign_rule() {
        // `sign(a)` repairs the `cos` argument for `a < 0`; the numeric check
        // is the guard.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let expr = parse_norm(&ctx, "1/(-5 + 3*cos(x))^2");
        assert_antiderivative_num(&ctx, expr, Symbol::new("x"), &[], &[0.3, 0.7, 1.1]);
    }

    #[test]
    fn t1_closed_form_symbolic_a() {
        // The `sin` arm is sign-robust and serves a symbolic `a`; the `cos`
        // arm declines for a symbolic `a` because its correctness depends on
        // `sign(a)`, which a symbol does not carry.
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let var = Symbol::new("x");
        let env = [(Symbol::new("a"), 3.0), (Symbol::new("b"), 2.0)];
        let es = parse_norm(&ctx, "1/(a + b*sin(x))^2");
        assert_antiderivative_num(&ctx, es, var, &env, &[0.3, 0.7, 1.1]);
        let ec = parse_norm(&ctx, "1/(a + b*cos(x))^2");
        let m = match_trig_denom(&ctx, ec, var).unwrap();
        assert!(weierstrass_base(&ctx, &m, var).is_none());
    }
}