ocas-calc 0.27.2

//! Half-power front-end for trig radical integrands (0.27.2 Wave C6).
//!
//! Detects integrands of the shape
//!
//! ```text
//!     C · S(cos u)^{p/2} du,     p odd,  deg_`cos` S ≤ 2,
//! ```
//!
//! with `C` free of `u`, and rewrites them into the canonical algebraic form
//! `∫ W(z)/√(Q(z)) dz` with `deg Q = 4` and `Q` already in Legendre normal
//! form, then hands the result to [`super::elliptic::integrate_elliptic`] and
//! substitutes back. When the rewrite does not produce a 3rd/4th-degree
//! radicand the module declines, so elementary engines keep first claim.
//!
//! # Substitutions and their derivation
//!
//! **Weierstrass (`t = tan(u/2)`).** With `sin u = 2t/(1+t²)`,
//! `cos u = (1−t²)/(1+t²)` and `du = 2 dt/(1+t²)`, a base that is linear in
//! `cos u` becomes `S = [(c₀+c₁) + (c₀−c₁)t²]/(1+t²)`, so
//! `S^{p/2} du = 2 P(t)^{p/2} dt/(1+t²)^{(p+2)/2}`; clearing denominators
//! gives the quartic-radicand integral `2C·√(P(t)(1+t²))-shaped` forms. This
//! is the classical route and the one the wave brief specifies.
//!
//! **Half-angle refinement.** The Weierstrass quartic `(1+t²)P(t)` is *not*
//! in Legendre normal form, and its reduction needs the complementary-modulus
//! integrals for the `t²` remainder — the `u = t/√(1+t²)` composition that the
//! elliptic engine knows how to perform anyway. The engine therefore applies
//! the same substitution `z = t/√(1+t²) = sin(u/2)` first, which turns
//!
//! - `S = c₀ + c₁cos u = (c₀+c₁) − 2c₁sin²(u/2) = A(1 − M z²)` with
//!   `A = c₀+c₁`, `M = 2c₁/A`, and `du = 2 dz/√(1−z²)`;
//! - `S = c₀ + c₂cos²u = (c₀+c₂) − c₂sin²u = A(1 − M z²)` with
//!   `A = c₀+c₂`, `M = c₂/A`, and `du = dz/√(1−z²)` (here `z = sin u`),
//!
//! into the *Legendre* quartic `Q = A(1−z²)(1−M z²)`. Since
//! `sin²(u/2) = t²/(1+t²)` and `sin²u = t²/(1+t²)` hold identically, the
//! emitted amplitude `asin(z)` has derivative exactly `1/2` (resp. `1`) times
//! the required chain factor on every interval where the substitution is
//! invertible, including past `u = π/2`. This module emits that refined form;
//! it is the same substitution, just composed in the order that keeps the
//! result in the first/second kind.
//!
//! Both branches require `S` to depend on `cos u` only (no `sin u` term and,
//! for the quadratic branch, no `cos u` term): a genuine `sin u` component
//! makes the radicand a general quartic whose reduction needs a Möbius
//! normalisation this wave does not implement, and `cos u` together with
//! `cos²u` makes the Legendre form unreachable by either half-angle.
//!
//! # Known gaps (declined, never guessed)
//!
//! - Kernels that are not a polynomial in `cos u`: `sec`, `tan`, `cot`,
//!   `tanh`, `coth`, `sinh` families, and any `S` mixing `sin u`/`cos u`.
//! - Prefactors that depend on `u` (`exp(x)·√(cos x)`, `sin(x)·√(S)`).
//! - Products of two different radicals and `|p| > 5`.
//! - Arguments other than the bare integration variable (`cos(c + d·x)`).
//! - `S` whose substituted leading coefficient `A` vanishes identically
//!   (a degenerate, non-squarefree radicand).
//!
//! # Branch / parameter caveat (important)
//!
//! The emitted closed form is real — and therefore the real antiderivative —
//! **only when the radicand's leading factor `A = c₀ + c₁` (respectively
//! `c₀ + c₂`) is positive**, because it carries `1/√A`. If a parameter
//! specialization makes `A < 0`, the same expression is still a formal
//! complex antiderivative, but the branch matching the real integrand is the
//! *complementary* one and would need the imaginary-modulus transformation
//! (`F(φ, m) ↦ F(φ', 1/m)/√m` and the matching `E`/`Π` identities), which is
//! not implemented. The sign of a symbolic `A` is not decidable here, so the
//! form is emitted regardless; the corpus oracle's deterministic dummy values
//! for `a` and `b` (see the `harness_parity_for_cos_family` test) happen to
//! give `a + b < 0`, where both this module's evaluator and the harness's
//! report the point as out of domain rather than as a disagreement.

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

use super::is_constant;

/// Integrate `C·S(cos u)^{p/2}` by rewriting to a quartic radical integrand.
///
/// Returns `None` when the rewrite does not produce a quartic radicand in
/// Legendre normal form.
pub(crate) fn integrate_half_power<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<Atom<'a>> {
    if super::node_count(expr) > 200 {
        return None;
    }
    // Canonicalise first: `normalize` folds `(S^{a})^{n}` for integer `n`, so
    // `1/(S)^{3/2}` reaches this front-end as the single half power `S^{-3/2}`
    // instead of a radical wrapped in a reciprocal.
    let expr = ocas_atom::normalize::normalize(ctx, expr);
    let (coeff, core) = split_constant(ctx, expr, var)?;
    let base_atom = ctx.var(var.as_str());
    let cos_atom = ctx.fun("cos", &[base_atom]);
    // The non-constant part must be a single half power. `normalize` folds
    // `(S^{a})^{n}` for integer `n`, so what is left is either `sqrt(S)`,
    // `S^{p/2}`, or `sqrt(S)^{m}` (the `sqrt` head is not a `Pow`, so the
    // folding rule does not reach it).
    let (base, p) = match core.node() {
        AtomNode::Fun(name, args) if name.as_str() == "sqrt" && args.len() == 1 => (args[0], 1i64),
        AtomNode::Pow(b, e) => {
            if let AtomNode::Fun(fname, fargs) = b.node()
                && fname.as_str() == "sqrt"
                && fargs.len() == 1
            {
                let (pe, qe) = exp_fraction(*e)?;
                if qe != 1 {
                    return None;
                }
                (fargs[0], pe)
            } else {
                let (pe, qe) = exp_fraction(*e)?;
                if qe != 2 {
                    return None;
                }
                (*b, pe)
            }
        }
        _ => return None,
    };
    if p % 2 == 0 || p.unsigned_abs() > 5 {
        return None;
    }
    let s = poly_in_base(ctx, base, cos_atom, var, 2)?;
    let c0 = s.first().copied().unwrap_or_else(|| ctx.num(0));
    let c1 = s.get(1).copied().unwrap_or_else(|| ctx.num(0));
    let c2 = s.get(2).copied().unwrap_or_else(|| ctx.num(0));
    let zero = |a: Atom<'a>| is_zero(ctx, a);
    // Pick the half-angle branch: linear in `cos u`, or pure `cos²u`.
    //
    // The base is written as `A − β z²` (a polynomial: no compound
    // denominator), because `S = A(1 − M z²)` with `M = β/A` is the same
    // value; the engine then sees a polynomial radicand and the Legendre
    // normal form drops out directly.
    let (a_coef, beta, scale, z_repl, sign_branch) = if zero(c2) {
        if zero(c1) {
            // Constant base: no radical, elementary engines own this.
            return None;
        }
        let a = cz(ctx, ctx.add(&[c0, c1]));
        if zero(a) {
            return None;
        }
        let beta = cz(ctx, ctx.mul(&[ctx.num(2), c1]));
        let half_var = ctx.mul(&[base_atom, ctx.pow(ctx.num(2), ctx.num(-1))]);
        (a, beta, ctx.num(2), ctx.fun("sin", &[half_var]), false)
    } else {
        if !zero(c1) {
            // A `cos u` term together with `cos²u`: neither half-angle reaches
            // Legendre normal form.
            return None;
        }
        let a = cz(ctx, ctx.add(&[c0, c2]));
        if zero(a) {
            return None;
        }
        (a, c2, ctx.num(1), ctx.fun("sin", &[base_atom]), true)
    };
    let z = fresh_symbol(expr, var)?;
    let zvar = ctx.var(z.as_str());
    let radicand = ctx.add(&[
        a_coef,
        ctx.mul(&[ctx.num(-1), beta, ctx.pow(zvar, ctx.num(2))]),
    ]);
    let one_minus_z2 = ctx.add(&[
        ctx.num(1),
        ctx.mul(&[ctx.num(-1), ctx.pow(zvar, ctx.num(2))]),
    ]);
    let rewritten = ctx.mul(&[
        coeff,
        scale,
        ctx.pow(radicand, half_exp(ctx, p)),
        ctx.pow(one_minus_z2, half_exp(ctx, -1)),
    ]);
    let g = super::elliptic::integrate_elliptic(ctx, rewritten, z)?;
    let back = super::replace_symbol(ctx, g, z, z_repl);
    // Branch correction for the `z = sin u` branch.
    //
    // `dz/du = cos u` is *signed*, so the engine's `√(1 − z²) = |cos u|`
    // describes the wrong sheet wherever `cos u < 0`; multiplying the whole
    // antiderivative by `sign(cos u) = cos u/|cos u|` — written here as
    // `cos u·(1 − sin²u)^{-1/2}` so no `abs` head is needed — restores the
    // correct sheet on every interval. The `z = sin(u/2)` branch of the
    // Weierstrass substitution is only claimed on `|u| < π`, where the
    // corresponding sign is `+1`, so no factor is emitted there.
    let back = if sign_branch {
        let one_minus_sin2 = ctx.add(&[
            ctx.num(1),
            ctx.mul(&[ctx.num(-1), ctx.pow(z_repl, ctx.num(2))]),
        ]);
        ctx.mul(&[
            ctx.fun("cos", &[base_atom]),
            ctx.pow(one_minus_sin2, half_exp(ctx, -1)),
            back,
        ])
    } else {
        back
    };
    Some(ocas_atom::normalize::normalize(ctx, back))
}

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

/// `2^{-1}`.
fn half<'a>(ctx: &'a AtomArena<'a>) -> Atom<'a> {
    ctx.pow(ctx.num(2), ctx.num(-1))
}

/// The atom `p/2`.
fn half_exp<'a>(ctx: &'a AtomArena<'a>, p: i64) -> Atom<'a> {
    ctx.mul(&[ctx.num(p), half(ctx)])
}

fn cz<'a>(ctx: &'a AtomArena<'a>, a: Atom<'a>) -> Atom<'a> {
    crate::ode::util::collect_terms(ctx, a)
}

fn is_zero<'a>(ctx: &'a AtomArena<'a>, a: Atom<'a>) -> bool {
    matches!(cz(ctx, a).node(), AtomNode::Num(0))
}

/// Split a product into `(constant part, non-constant part)`.
fn split_constant<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
) -> Option<(Atom<'a>, Atom<'a>)> {
    let factors: Vec<Atom<'a>> = match expr.node() {
        AtomNode::Mul(args) => args.to_vec(),
        _ => vec![expr],
    };
    let (c, nc): (Vec<Atom<'a>>, Vec<Atom<'a>>) =
        factors.into_iter().partition(|f| is_constant(*f, var));
    if nc.is_empty() {
        return None;
    }
    let core = if nc.len() == 1 { nc[0] } else { ctx.mul(&nc) };
    let coeff = match c.len() {
        0 => ctx.num(1),
        1 => c[0],
        _ => ctx.mul(&c),
    };
    Some((coeff, core))
}

/// Parse a rational exponent atom `p/q` (same shapes as the elliptic engine).
fn exp_fraction<'a>(exp: Atom<'a>) -> Option<(i64, i64)> {
    match exp.node() {
        AtomNode::Num(n) => Some((*n, 1)),
        AtomNode::Pow(b, e) => {
            if let (AtomNode::Num(bb), AtomNode::Num(ee)) = (b.node(), e.node())
                && *ee == -1
                && *bb > 0
            {
                return Some((1, *bb));
            }
            None
        }
        AtomNode::Mul(args) => {
            let mut num: Option<i64> = None;
            let mut den: Option<i64> = None;
            for a in args.iter() {
                match a.node() {
                    AtomNode::Num(n) => {
                        if num.is_some() {
                            return None;
                        }
                        num = Some(*n);
                    }
                    AtomNode::Pow(b, e) => {
                        if let (AtomNode::Num(bb), AtomNode::Num(ee)) = (b.node(), e.node())
                            && *ee == -1
                            && *bb > 0
                        {
                            if den.is_some() {
                                return None;
                            }
                            den = Some(*bb);
                        } else {
                            return None;
                        }
                    }
                    _ => return None,
                }
            }
            match (num, den) {
                (Some(p), Some(q)) => Some((p, q)),
                (Some(p), None) => Some((p, 1)),
                (None, Some(q)) => Some((1, q)),
                (None, None) => None,
            }
        }
        _ => None,
    }
}

/// Structural occurrence test.
fn contains_atom<'a>(expr: Atom<'a>, needle: Atom<'a>) -> bool {
    if expr == needle {
        return true;
    }
    match expr.node() {
        AtomNode::Num(_) | AtomNode::Var(_) => false,
        AtomNode::Pow(b, e) => contains_atom(*b, needle) || contains_atom(*e, needle),
        AtomNode::Add(args) | AtomNode::Mul(args) | AtomNode::Fun(_, args) => {
            args.iter().any(|a| contains_atom(*a, needle))
        }
    }
}

/// A substitution variable that collides neither with `var` nor with any
/// symbol already present in `expr`.
fn fresh_symbol<'a>(expr: Atom<'a>, var: Symbol) -> Option<Symbol> {
    for name in ["_z", "_s", "_w", "_v"] {
        let s = Symbol::new(name);
        if s != var && !super::contains_symbol(expr, s) {
            return Some(s);
        }
    }
    None
}

/// Coefficients `[c₀, …, c_n]` of `expr` as a polynomial in `base`, or `None`
/// when `expr` is not such a polynomial.
///
/// Coefficients must be free of `var`; this is what rejects `√(sin u)` (whose
/// "coefficient" would be `sin u` itself) and any other base that merely
/// contains `var` without being built from `base`.
fn poly_in_base<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    base: Atom<'a>,
    var: Symbol,
    max_deg: usize,
) -> Option<Vec<Atom<'a>>> {
    let expanded = cz(ctx, expr);
    let terms: Vec<Atom<'a>> = match expanded.node() {
        AtomNode::Add(args) => args.to_vec(),
        _ => vec![expanded],
    };
    let mut out: Vec<Atom<'a>> = Vec::new();
    for t in terms {
        let (deg, coeff) = monomial(ctx, t, base, max_deg)?;
        // Coefficients must not depend on the integration variable: this is
        // what rejects a base that merely *contains* `var` (e.g. `sin x`
        // against the base `cos x`).
        if !is_constant(coeff, var) {
            return None;
        }
        if out.len() <= deg {
            out.resize(deg + 1, ctx.num(0));
        }
        out[deg] = cz(ctx, ctx.add(&[out[deg], coeff]));
    }
    while let Some(last) = out.last().copied() {
        if is_zero(ctx, last) {
            out.pop();
        } else {
            break;
        }
    }
    Some(out)
}

/// `(degree in base, coefficient)` for one additive term.
fn monomial<'a>(
    ctx: &'a AtomArena<'a>,
    t: Atom<'a>,
    base: Atom<'a>,
    max_deg: usize,
) -> Option<(usize, Atom<'a>)> {
    if t == base {
        return Some((1, ctx.num(1)));
    }
    match t.node() {
        AtomNode::Num(_) => Some((0, t)),
        AtomNode::Var(_) => {
            if contains_atom(t, base) {
                None
            } else {
                Some((0, t))
            }
        }
        AtomNode::Fun(_, _) => {
            if contains_atom(t, base) {
                None
            } else {
                Some((0, t))
            }
        }
        AtomNode::Add(_) => None,
        AtomNode::Pow(b, e) => {
            if *b == base {
                if let AtomNode::Num(n) = e.node() {
                    let n = usize::try_from(*n).ok()?;
                    if n <= max_deg {
                        return Some((n, ctx.num(1)));
                    }
                }
                return None;
            }
            if contains_atom(t, base) {
                None
            } else {
                Some((0, t))
            }
        }
        AtomNode::Mul(args) => {
            let mut deg = 0usize;
            let mut coeffs: Vec<Atom<'a>> = Vec::with_capacity(args.len());
            for a in args.iter() {
                let (d, c) = monomial(ctx, *a, base, max_deg)?;
                deg = deg.checked_add(d)?;
                if deg > max_deg {
                    return None;
                }
                coeffs.push(c);
            }
            Some((deg, ctx.mul(&coeffs)))
        }
    }
}

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

#[cfg(test)]
mod tests {
    use super::super::elliptic::testnum::{diff_local, eval_f64};
    use super::*;
    use ocas_core::arena::Arena;

    fn parse<'a>(ctx: &'a AtomArena<'a>, s: &str) -> Atom<'a> {
        ocas_parse::parse(ctx, s).expect("parse")
    }

    /// Run the front-end, require a radical-free result, and check
    /// `d/dx result == integrand` numerically at the sample points.
    fn assert_antiderivative_num<'a>(
        ctx: &'a AtomArena<'a>,
        src: &str,
        consts: &[(Symbol, f64)],
        samples: &[f64],
    ) {
        let var = Symbol::new("x");
        let integrand = parse(ctx, src);
        let result =
            integrate_half_power(ctx, integrand, var).unwrap_or_else(|| panic!("declined: {src}"));
        assert!(
            !result.to_string().contains("Integral"),
            "residue for {src}: {result}"
        );
        let d = diff_local(ctx, result, var);
        let mut checked = 0usize;
        for &xv in samples {
            let mut env = consts.to_vec();
            env.push((var, xv));
            let lhs = match eval_f64(d, &env) {
                Some(v) => v,
                None => continue,
            };
            let rhs = eval_f64(integrand, &env).expect("eval integrand");
            let tol = 1e-9 * rhs.abs().max(1.0);
            assert!(
                (lhs - rhs).abs() < tol,
                "at x={xv}: diff={lhs} integrand={rhs} (src: {src}, result: {result})"
            );
            checked += 1;
        }
        assert!(checked >= 2, "only {checked} usable samples for {src}");
    }

    fn declines<'a>(ctx: &'a AtomArena<'a>, src: &str) {
        let e = parse(ctx, src);
        let var = Symbol::new("x");
        assert!(
            integrate_half_power(ctx, e, var).is_none(),
            "expected decline for {src}"
        );
    }

    #[test]
    fn linear_in_cos_first_kind() {
        // ∫dx/√(a+b·cos x) = (2/√(a+b))·F(x/2, 2b/(a+b))
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 1.0)];
        assert_antiderivative_num(
            &ctx,
            "1/sqrt(a+b*cos(x))",
            &env,
            &[-2.0, -1.0, 0.4, 1.2, 2.4],
        );
    }

    #[test]
    fn linear_in_cos_second_kind() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 1.0)];
        // ∫√(a+b·cos x) dx = 2√(a+b)·E(x/2, 2b/(a+b))
        assert_antiderivative_num(&ctx, "sqrt(a+b*cos(x))", &env, &[-2.0, -0.8, 0.5, 1.5, 2.5]);
    }

    #[test]
    fn linear_in_cos_hermite() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 1.0)];
        assert_antiderivative_num(
            &ctx,
            "1/(a+b*cos(x))^(3/2)",
            &env,
            &[-2.2, -0.6, 0.3, 1.4, 2.6],
        );
    }

    #[test]
    fn quadratic_in_cos_first_kind() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 1.0)];
        // ∫dx/√(a+b·cos²x) = F(x, b/(a+b))/√(a+b)
        assert_antiderivative_num(
            &ctx,
            "1/sqrt(a+b*cos(x)^2)",
            &env,
            &[-2.4, -1.2, -0.4, 0.6, 1.9, 2.6],
        );
    }

    #[test]
    fn quadratic_in_cos_second_kind() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 1.0)];
        // ∫√(a+b·cos²x) dx = √(a+b)·E(x, b/(a+b))
        assert_antiderivative_num(
            &ctx,
            "sqrt(a+b*cos(x)^2)",
            &env,
            &[-2.4, -1.1, -0.3, 0.7, 1.8, 2.7],
        );
    }

    #[test]
    fn quadratic_in_cos_hermite() {
        // The brief's headline case: ∫dx/(a+b·cos²x)^{3/2}
        //   = E(x, M)/(a√(a+b)) − (b/a)·sin x cos x / ((a+b)√(a+b·cos²x)),
        //   M = b/(a+b).
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let env = [(Symbol::new("a"), 2.0), (Symbol::new("b"), 1.0)];
        assert_antiderivative_num(
            &ctx,
            "1/(a+b*cos(x)^2)^(3/2)",
            &env,
            &[-2.4, -1.3, -0.3, 0.8, 1.7, 2.8],
        );
    }

    #[test]
    fn numeric_coefficients_and_extra_powers() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        assert_antiderivative_num(&ctx, "1/sqrt(2+cos(x))", &[], &[-2.0, -0.5, 0.6, 2.0]);
        assert_antiderivative_num(&ctx, "1/(2+cos(x)^2)^(3/2)", &[], &[-2.2, -0.7, 0.5, 2.2]);
        assert_antiderivative_num(&ctx, "3/sqrt(2+cos(x))", &[], &[-1.8, -0.4, 0.9, 2.1]);
        // Symbolic coefficient on the second-kind branch.
        let env = [(Symbol::new("a"), 3.0), (Symbol::new("b"), 1.0)];
        assert_antiderivative_num(&ctx, "sqrt(a+b*cos(x)^2)", &env, &[-2.1, -0.8, 0.6, 2.2]);
    }

    /// Harness-parity check.
    ///
    /// `ocas-tests/src/integral_eval.rs` verifies a solved case with a 5-point
    /// central difference of the antiderivative at the fixed abscissae below,
    /// assigning each free parameter a deterministic dummy value keyed by its
    /// name (hash → `TABLE`, random sign). This test replicates that oracle
    /// exactly over the *whole* pipeline (`ocas_calc::integrate`) and asserts
    /// the two properties that matter:
    ///
    /// 1. With the harness's own dummy values for `a`/`b` the antiderivative
    ///    may land outside the evaluated branch (`a + b < 0` makes the emitted
    ///    `(a+b)^{-1/2}` complex, and the harness's `m = b/(a+b)` then exceeds
    ///    its `|m| ≤ 0.9` window), but it must **never** disagree where both
    ///    sides are evaluable.
    /// 2. On the same abscissae — which include `1.91` and `2.53`, i.e. past
    ///    `π/2` — with a parameter regime where the antiderivative is real,
    ///    every usable sample must agree to `1e-5` relative. This is the
    ///    regression guard for the `z = sin u` branch correction.
    #[test]
    fn harness_parity_for_cos_family() {
        const TABLE: [f64; 12] = [
            2.0,
            3.0,
            5.0,
            7.0,
            0.5,
            1.5,
            0.25,
            11.0,
            1.0 / 3.0,
            4.0,
            6.0,
            0.75,
        ];
        const SAMPLES: [f64; 8] = [-1.7, -0.9, -0.37, 0.31, 0.77, 1.23, 1.91, 2.53];
        fn param_value(name: &str) -> f64 {
            let mut h: u64 = 0xcbf2_9ce4_8422_2325;
            for b in name.bytes() {
                h ^= u64::from(b);
                h = h.wrapping_mul(0x0000_0100_0000_01b3);
            }
            let v = TABLE[(h % 12) as usize];
            if (h >> 8) & 1 == 0 { v } else { -v }
        }
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let var = Symbol::new("x");
        let cases = [
            "1/(a + b*cos(x)^2)^(3/2)",
            "1/sqrt(a + b*cos(x))",
            "1/sqrt(a + b*cos(x)^2)",
            "sqrt(a + b*cos(x))",
            "sqrt(a + b*cos(x)^2)",
        ];
        for mode in ["harness dummy parameters", "real regime a=2 b=1"] {
            for src in cases {
                let env = if mode.starts_with("harness") {
                    [
                        (Symbol::new("a"), param_value("a")),
                        (Symbol::new("b"), param_value("b")),
                    ]
                } else {
                    [(Symbol::new("a"), 2.0), (Symbol::new("b"), 1.0)]
                };
                let integrand = parse(&ctx, src);
                let result = crate::integrate(&ctx, integrand, var);
                let text = result.to_string();
                assert!(!text.contains("Integral("), "{src} fell back: {text}");
                let (mut checked, mut worst) = (0usize, 0.0f64);
                for &x in SAMPLES.iter() {
                    let h = 1e-4 * x.abs().max(1.0);
                    let eval_at = |t: f64| {
                        let mut e = env.to_vec();
                        e.push((var, t));
                        eval_f64(result, &e)
                    };
                    let (fm2, fm1, fp1, fp2) = (
                        eval_at(x - 2.0 * h),
                        eval_at(x - h),
                        eval_at(x + h),
                        eval_at(x + 2.0 * h),
                    );
                    let mut e = env.to_vec();
                    e.push((var, x));
                    let rhs = eval_f64(integrand, &e);
                    if let (Some(a), Some(b), Some(c), Some(d), Some(f)) = (fm2, fm1, fp1, fp2, rhs)
                    {
                        let deriv = (-d + 8.0 * c - 8.0 * b + a) / (12.0 * h);
                        if deriv.is_finite() {
                            checked += 1;
                            worst = worst.max((deriv - f).abs() / f.abs().max(1.0));
                        }
                    }
                }
                if mode.starts_with("harness") {
                    // Outside the evaluated branch the harness's oracle reports
                    // Domain for every sample (Indeterminate), never a
                    // mismatch; assert only the "no disagreement" half.
                    assert!(worst <= 1e-5, "{src} [{mode}]: worst rel {worst:e}");
                } else {
                    assert!(
                        checked >= 2,
                        "{src} [{mode}]: only {checked} usable samples"
                    );
                    assert!(
                        worst <= 1e-5,
                        "{src} [{mode}]: worst rel {worst:e} ({text})"
                    );
                    println!("{src}: checked={checked} worst_rel={worst:e}");
                }
            }
        }
    }

    #[test]
    fn declines_unsupported_shapes() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // `sin` component: the radicand is a general quartic.
        declines(&ctx, "sqrt(sin(x))");
        declines(&ctx, "1/sqrt(a+b*sin(x))");
        // `cos` and `cos²` together: neither half-angle reaches Legendre form.
        declines(&ctx, "1/sqrt(a+b*cos(x)+c*cos(x)^2)");
        // Non-polynomial kernels.
        declines(&ctx, "1/(a+b*sec(x))^(5/2)");
        declines(&ctx, "sqrt(tan(x))");
        declines(&ctx, "sqrt(a+b*sinh(x))");
        // Variable-dependent prefactor.
        declines(&ctx, "exp(x)*sqrt(cos(x))");
        declines(&ctx, "sin(x)*sqrt(a+b*cos(x))");
        declines(&ctx, "x/sqrt(a+b*cos(x))");
        // Radical-free / plain algebraic inputs.
        declines(&ctx, "1/sqrt(x^5+1)");
        declines(&ctx, "1/(1+cos(x))");
        declines(&ctx, "cos(x)");
        declines(&ctx, "sqrt(1+x^2)");
        // Higher half powers.
        declines(&ctx, "1/(a+b*cos(x))^(7/2)");
        // Positive third powers leave a non-constant `S(u)√Q` remainder, which
        // the elliptic engine deliberately declines (see its module docs).
        declines(&ctx, "(a+b*cos(x))^(3/2)");
        declines(&ctx, "(a+b*cos(x)^2)^(5/2)");
        // Degenerate base (A ≡ 0): 1 − cos x = 2sin²(x/2) is not our family.
        declines(&ctx, "1/sqrt(1-cos(x))");
        // Non-linear argument.
        declines(&ctx, "1/sqrt(a+b*cos(2*x))");
    }

    #[test]
    fn stress_stable_outcomes() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let var = Symbol::new("x");
        let solid = parse(&ctx, "1/sqrt(a+b*cos(x))");
        let mut first: Option<String> = None;
        for i in 0..300 {
            let r = integrate_half_power(&ctx, solid, var);
            let s = r.map(|a| a.to_string());
            if i == 0 {
                first = s.clone();
            }
            assert_eq!(s, first, "non-deterministic outcome at iteration {i}");
        }
        assert!(first.is_some(), "stable outcome must be a solve");
        let none = parse(&ctx, "sqrt(sin(x))");
        for i in 0..300 {
            assert!(
                integrate_half_power(&ctx, none, var).is_none(),
                "decline not stable at iteration {i}"
            );
        }
    }
}