symplex 0.39.0

//! Lazard-Rioboo-Trager `log_to_real` conversion.
//!
//! Converts algebraic logarithmic terms from the Rothstein-Trager algorithm
//! into real-valued `ln + atan` expressions.
//!
//! Given an irreducible factor `q(t)` of the Rothstein-Trager resultant
//! (degree ≥ 2) and the degree-1 PRS member `h(t, x)`, this module computes:
//!
//! ```text
//! Σ_{α: q(α)=0} α · ln(h(α, x))
//! ```
//!
//! as a real-valued sum of `ln` and `atan` terms, using the identity:
//!
//! ```text
//! α·ln(h(α,x)) + ᾱ·ln(h(ᾱ,x)) = u·ln(A²+B²) + v·log_to_atan(A, B)
//! ```
//!
//! where `α = u + iv`, `h(α,x) = A(x) + iB(x)`, and `log_to_atan`
//! converts the imaginary part to a sum of arctangents.
//!
//! The identity needs real coefficients: `q ∈ ℚ[t]` and `h ∈ ℚ(t)[x]`, so
//! that the non-real roots of `q` come in conjugate pairs and
//! `h(ᾱ, x) = conj(h(α, x))` for real `x`.  That holds for every caller:
//! the rational integrator (`try_risch_rational`) takes `ℚ[x]` numerators
//! and denominators only, so an integrand with a parameter never reaches
//! this module, and its real forms (`ln|v(x)|` with `v ∈ ℚ[x]`, `ln(A² +
//! B²)`, `atan`) are real for every real `x` — decision D4 (0.31), which
//! allows the real forms only where the argument is real under the
//! declared assumptions, asks nothing more here (audited in 0.31).
//!
//! # Algorithm
//!
//! 1. Solve `q(t) = 0` for complex roots (exact radicals via Cardano/Ferrari
//!    for degree ≤ 4).
//! 2. Decompose each root into `(Re, Im)` via [`as_real_imag`](crate::base::complex::as_real_imag).
//! 3. For each conjugate pair `(u, v)` with `v > 0`:
//!    a. Evaluate `h(u+iv, x)` and separate `A(x) + iB(x)`.
//!    b. Emit `u · ln(A² + B²) + v · log_to_atan(A, B)`.
//!
//! When the roots of `q` are `ρ·ζ` with `ρ^L` rational and `ζ` roots of
//! unity (`t^M ≡ c mod q`: the residues of `xᵏ/(xⁿ + a)`), steps 2–3 run
//! exactly in `ℚ(ζ_N)[ρ]` instead (`cyclotomic_log_to_real`), with the
//! coefficients over the cosines of the real subfield (de Moivre).
//!
//! # References
//!
//! - Lazard & Rioboo, "Integration of rational functions: rational computation
//!   of the logarithmic part", 1990
//! - Bronstein, *Symbolic Integration I*, §2.5
//! - SymPy `integrals/rationaltools.py`, `log_to_real` and `log_to_atan`

use num_bigint::BigInt;
use num_rational::Ratio;

use crate::base::arena::Arena;
use crate::base::node::{ExprId, ExprNode};
use crate::base::numeric::Q;
use crate::poly::dense::Poly;
use crate::poly::generic::GenPoly;
use crate::poly::ratfn::RationalFn;
use crate::poly::traits::Ring;

// ═══════════════════════════════════════════════════════════════════════════
// Helper: Poly → GenPoly<RationalFn>
// ═══════════════════════════════════════════════════════════════════════════

/// Convert a `Poly` (polynomial in `x` with rational coefficients) to a
/// `GenPoly<RationalFn>` where each coefficient is a constant rational
/// function (denominator = 1, independent of `t`).
pub(crate) fn poly_to_genpoly_rf(p: &Poly) -> GenPoly<RationalFn> {
    let coeffs: Vec<RationalFn> = p
        .coeffs()
        .iter()
        .map(|c| RationalFn::from_rational(c.clone()))
        .collect();
    GenPoly::from_coeffs(coeffs)
}

/// Convert a `Poly` (polynomial in `x` with rational coefficients) to a
/// `GenPoly<RationalFn>` where each coefficient is multiplied by `t`.
///
/// If `p(x) = Σ c_k · x^k`, this returns `Σ (c_k · t) · x^k` where
/// `t` is the variable of the `RationalFn` coefficient type.  Only the tests
/// build `A − t·D'` this way after 0.28.0 (the Lazard–Rioboo–Trager sequence
/// works in `ℤ[t][x]`).
#[cfg(test)]
pub(crate) fn poly_to_genpoly_rf_times_t(p: &Poly) -> GenPoly<RationalFn> {
    let t_poly = Poly::from_coeffs(vec![
        Ratio::from_integer(BigInt::from(0)),
        Ratio::from_integer(BigInt::from(1)),
    ]); // t
    let coeffs: Vec<RationalFn> = p
        .coeffs()
        .iter()
        .map(|c| {
            if c.is_zero() {
                RationalFn::zero()
            } else {
                let c_poly = Poly::constant(c.clone());
                let c_times_t = &c_poly * &t_poly;
                RationalFn::from_poly(c_times_t)
            }
        })
        .collect();
    GenPoly::from_coeffs(coeffs)
}

// ═══════════════════════════════════════════════════════════════════════════
// log_to_atan — convert imaginary log part to arctangent sum
// ═══════════════════════════════════════════════════════════════════════════

/// Convert the imaginary part of a complex logarithm to a sum of arctangents.
///
/// Given `A(x) = a₁·x + a₀` and `B(x) = b₁·x + b₀`, computes an
/// expression `F(x)` such that:
///
/// ```text
/// dF/dx = d/dx [I · ln((A + I·B) / (A - I·B))]
/// ```
///
/// For the common case where `B` is constant (`b₁ ≈ 0`), this is
/// simply `2 · atan(A / B)`.
///
/// For the general case where both `A` and `B` are linear, the result
/// is a sum of two `atan` terms computed via the Bézout identity.
///
/// Returns `None` if the inputs are degenerate (e.g., `B ≈ 0`).
pub(crate) fn log_to_atan_deg1(
    arena: &mut Arena,
    var: ExprId,
    a1: ExprId,
    a0: ExprId,
    b1: ExprId,
    b0: ExprId,
) -> Option<ExprId> {
    let two = arena.int(2);

    // Check if b1 = 0 (B is constant — the common case).  An undecided
    // zero test (`None`: the value is 0 to the evaluator's precision cap
    // without a proof) declines: taking it for nonzero divides by `b₁`.
    let b1_is_zero = crate::poly::algebraic::is_zero_checked(arena, b1)?;

    if b1_is_zero {
        tracing::debug!("log_to_atan_deg1: b1 ≈ 0, B is constant → single atan");

        // Check b0 ≠ 0 (bail unless we can confirm nonzero, since we divide by b0)
        if crate::poly::algebraic::is_zero_checked(arena, b0) != Some(false) {
            tracing::debug!("log_to_atan_deg1: b0 ≈ 0, degenerate → None");
            return None;
        }

        // F = 2 · atan((a₁·x + a₀) / b₀)
        let a_x = arena.mul(&[a1, var]);
        let a_poly = arena.add(&[a_x, a0]);
        let ratio = arena.div(a_poly, b0);
        let atan_val = arena.atan(ratio);
        let result = arena.mul(&[two, atan_val]);
        return Some(result);
    }

    tracing::debug!("log_to_atan_deg1: both A and B are degree 1 → Bézout path");

    // Both A and B are linear: A = a₁x + a₀, B = b₁x + b₀.
    //
    // Polynomial division: A / B gives quotient q = a₁/b₁, remainder r = a₀ - q·b₀.
    let q = arena.div(a1, b1);
    let q_b0 = arena.mul(&[q, b0]);
    let r = arena.sub(a0, q_b0);
    let r = crate::transforms::eval::eval(arena, r);

    let r_is_zero = crate::poly::algebraic::is_zero_checked(arena, r)?;

    if r_is_zero {
        // Division is exact: F = 2 · atan(q)
        tracing::trace!("log_to_atan_deg1: remainder is zero → single constant atan");
        let atan_q = arena.atan(q);
        return Some(arena.mul(&[two, atan_q]));
    }

    // General case: compute via explicit Bézout identity.
    //
    // For coprime degree-1 polynomials A = a₁x+a₀ and B = b₁x+b₀,
    // the extended GCD gives s, t such that s·B + t·(-A) = 1.
    // Solving the 2×2 linear system:
    //   s·b₁ - t·a₁ = 0  (coefficient of x)
    //   s·b₀ - t·a₀ = 1  (constant term)
    // yields s = a₁/D, t = b₁/D where D = a₁·b₀ - b₁·a₀.
    //
    // Then u = (A·s + B·t) / gcd = (a₁·s + b₁·t)·x + (a₀·s + b₀·t)
    // = ((a₁² + b₁²)/D)·x + ((a₀·a₁ + b₀·b₁)/D).
    //
    // The result is: F = 2·atan(u) + log_to_atan(s, t)
    // where s = a₁/D and t = b₁/D are both constants (degree 0),
    // so log_to_atan(s, t) = 2·atan(s/t) = 2·atan(a₁/b₁).

    let a1_b0 = arena.mul(&[a1, b0]);
    let b1_a0 = arena.mul(&[b1, a0]);
    let det = arena.sub(a1_b0, b1_a0);
    let det = crate::transforms::eval::eval(arena, det);

    if crate::poly::algebraic::is_zero_checked(arena, det) != Some(false) {
        tracing::debug!("log_to_atan_deg1: determinant ≈ 0 (A and B proportional) → None");
        return None;
    }

    // u(x) = ((a₁² + b₁²) / D) · x + ((a₀·a₁ + b₀·b₁) / D)
    let a1_sq = arena.mul(&[a1, a1]);
    let b1_sq = arena.mul(&[b1, b1]);
    let u_x_numer = arena.add(&[a1_sq, b1_sq]);
    let u_x_coeff = arena.div(u_x_numer, det);

    let a0_a1 = arena.mul(&[a0, a1]);
    let b0_b1 = arena.mul(&[b0, b1]);
    let u_c_numer = arena.add(&[a0_a1, b0_b1]);
    let u_const = arena.div(u_c_numer, det);

    let u_x_term = arena.mul(&[u_x_coeff, var]);
    let u_expr = arena.add(&[u_x_term, u_const]);
    let u_expr = crate::transforms::eval::eval(arena, u_expr);
    let atan_u = arena.atan(u_expr);
    let first_atan = arena.mul(&[two, atan_u]);

    // Second term: 2·atan(a₁/b₁)
    let a1_over_b1 = arena.div(a1, b1);
    let a1_over_b1 = crate::transforms::eval::eval(arena, a1_over_b1);
    let atan_a1b1 = arena.atan(a1_over_b1);
    let second_atan = arena.mul(&[two, atan_a1b1]);

    tracing::trace!("log_to_atan_deg1: Bézout path produced 2 atan terms");
    Some(arena.add(&[first_atan, second_atan]))
}

// ═══════════════════════════════════════════════════════════════════════════
// Quadratic factors: exact conversion over ℚ[x]
// ═══════════════════════════════════════════════════════════════════════════

/// `Σ_{q(α)=0} α·ln(S(α, x))` for an irreducible quadratic `q`, as a real
/// expression, with every polynomial computed exactly over `ℚ`.
///
/// Write `q = a·t² + b·t + c`, `u = −b/2a` and `w = √|b² − 4ac| / 2a`.  The
/// coefficients of `S` lie in `ℚ[t]/(q)`, so `S(t, x) = Σ (p_k + s_k·t)·x^k`
/// and at a root `α = u ± w` (real roots) or `u ± i·w` (complex roots)
/// `S(α, x) = P(x) ± w·Q(x)` resp. `P ± i·w·Q`, with
/// `P = Σ (p_k + s_k·u)·x^k` and `Q = Σ s_k·x^k` in `ℚ[x]`.  Then
///
/// - real roots: `(u + w)·ln|P + w·Q| + (u − w)·ln|P − w·Q|`;
/// - complex roots: `u·ln(P² + w²·Q²) + w·LogToAtan(P, w·Q)`
///   ([`log_to_atan_surd`]), where `P² + w²Q² ∈ ℚ[x]` has no real zero.
///
/// `S` is the log argument of the Lazard–Rioboo–Trager algorithm, monic
/// in `x`, of any degree; `None` if `q` is not quadratic or a coefficient
/// of `S` is not a polynomial in `t`.  (Bronstein, *Symbolic Integration I*,
/// §2.8; SymPy's `log_to_real` does the same with the roots found by
/// `roots`.)
pub(crate) fn quadratic_log_to_real(
    arena: &mut Arena,
    var: ExprId,
    q: &Poly,
    s: &GenPoly<RationalFn>,
) -> Option<Vec<ExprId>> {
    if q.degree() != Some(2) {
        return None;
    }
    let (c, b, a) = (q.coeff(0), q.coeff(1), q.coeff(2));
    let two = Q::from_integer(BigInt::from(2));
    let u = -&b / (&two * &a);
    let disc = &b * &b - Q::from_integer(BigInt::from(4)) * &a * &c;
    if disc.is_zero() {
        return None;
    }
    let w2 = num_traits::Signed::abs(&disc) / (Q::from_integer(BigInt::from(4)) * &a * &a);

    // P and Q from the coefficients of S, reduced modulo q.
    let mut p_coeffs = Vec::with_capacity(s.coeffs.len());
    let mut q_coeffs = Vec::with_capacity(s.coeffs.len());
    for coeff in &s.coeffs {
        if !coeff.denom().is_constant() {
            return None;
        }
        let lin = coeff.numer().rem(q).scale(&coeff.denom().coeff(0).recip());
        let (p0, p1) = (lin.coeff(0), lin.coeff(1));
        p_coeffs.push(&p0 + &p1 * &u);
        q_coeffs.push(p1);
    }
    let big_p = Poly::from_coeffs(p_coeffs);
    let big_q = Poly::from_coeffs(q_coeffs);
    if big_q.is_zero() {
        // S(α, x) would not depend on the root: impossible for the log
        // argument of distinct residues.
        return None;
    }

    let u_id = arena.num_ratio(u.clone());
    let w = surd(arena, &w2);
    let p_expr = crate::poly::polybridge::poly_to_expr(arena, &big_p, var);
    let q_expr = crate::poly::polybridge::poly_to_expr(arena, &big_q, var);
    let mut terms = Vec::new();

    if num_traits::Signed::is_positive(&disc) {
        let w_q = arena.mul(&[w, q_expr]);
        let neg_w = arena.neg(w);
        for (sign_w, sign_wq) in [(w, w_q), (neg_w, arena.neg(w_q))] {
            let alpha = arena.add(&[u_id, sign_w]);
            let arg = arena.add(&[p_expr, sign_wq]);
            let abs_arg = arena.abs(arg);
            let ln = arena.ln(abs_arg);
            terms.push(arena.mul(&[alpha, ln]));
        }
    } else {
        if !u.is_zero() {
            let norm = &(&big_p * &big_p) + &(&big_q * &big_q).scale(&w2);
            let norm_expr = crate::poly::polybridge::poly_to_expr(arena, &norm, var);
            let ln = arena.ln(norm_expr);
            terms.push(arena.mul(&[u_id, ln]));
        }
        let atans = log_to_atan_surd(arena, var, &big_p, &big_q, &w2)?;
        terms.push(arena.mul(&[w, atans]));
    }
    Some(terms)
}

/// `√r` for a positive rational `r`, canonicalised by the arena (`√(1/64)`
/// is `1/8`, `√(3/4)` is `√3/2`).
fn surd(arena: &mut Arena, r: &Q) -> ExprId {
    let r_id = arena.num_ratio(r.clone());
    let half = arena.rational(1, 2);
    arena.pow(r_id, half)
}

/// `2·atan(u(x)/√r2)`, written `±2·atan((c/√r2)·p(x))` with `c > 0` and
/// `p` primitive over ℤ with a positive leading coefficient (atan is odd):
/// `2·atan(√3/3·(2x² − 1))` rather than `2·atan(4√3·(x²/6 − 1/12))`.
fn two_atan_poly_over_surd(arena: &mut Arena, u: &Poly, var: ExprId, r2: &Q) -> ExprId {
    let c = u.content();
    if c.is_zero() {
        return arena.zero;
    }
    let negative = u
        .leading_coeff()
        .is_some_and(num_traits::Signed::is_negative);
    let p = u.scale(&if negative { -c.recip() } else { c.recip() });
    let p_expr = crate::poly::polybridge::poly_to_expr(arena, &p, var);
    let c_id = arena.num_ratio(c);
    let inv_r = surd(arena, &r2.recip());
    let arg = arena.mul(&[c_id, inv_r, p_expr]);
    let atan = arena.atan(arg);
    let two = arena.int(if negative { -2 } else { 2 });
    arena.mul(&[two, atan])
}

/// Bronstein's `LogToAtan(A, r·B)` for `A, B ∈ ℚ[x]` and a real constant
/// `r = √r2` (`r2 ∈ ℚ`, positive): a sum `F` of arctangents of real
/// polynomials with `F' = (i·ln((A + i·rB)/(A − i·rB)))'`.
///
/// Bronstein's recursion (*Symbolic Integration I*, §2.8) stays in `ℚ[x]`
/// because only `r²` enters the polynomial arithmetic: with
/// `s·B − t·A = h` over `ℚ` (so `(s/r)·(rB) − t·A = h`),
///
/// - `rB | A`: `2·atan(A/(rB))`;
/// - `deg A < deg B`: `LogToAtan(−rB, A)`, which is `LogToAtan(−B, A/r)` up
///   to the common real factor `r` (the value is invariant under scaling
///   both arguments), i.e. `B ← A`, `A ← −B`, `r2 ← 1/r2`;
/// - otherwise `2·atan((A·s + r2·B·t)/(h·r)) + LogToAtan(s/r, t)`, the
///   last again `LogToAtan(s, r·t)` up to the factor `r`.
///
/// `None` if `B` is zero.
pub(crate) fn log_to_atan_surd(
    arena: &mut Arena,
    var: ExprId,
    a: &Poly,
    b: &Poly,
    r2: &Q,
) -> Option<ExprId> {
    if b.is_zero() {
        return None;
    }
    let (mut a, mut b, mut r2) = (a.clone(), b.clone(), r2.clone());
    let mut atans = Vec::new();
    // Each step lowers deg A + deg B, as in Euclid's algorithm.
    let max_steps = 2 * (a.degree().unwrap_or(0) + b.degree().unwrap_or(0)) + 4;
    for _ in 0..max_steps {
        if b.is_zero() {
            break;
        }
        if a.is_zero() {
            // i·ln(i·rB / (−i·rB)) is constant.
            break;
        }
        if a.degree() < b.degree() {
            (a, b) = (-&b, a);
            r2 = r2.recip();
        }
        let (quot, rem) = a.div_rem(&b);
        if rem.is_zero() {
            atans.push(two_atan_poly_over_surd(arena, &quot, var, &r2));
            break;
        }
        let eg = Poly::extended_gcd(&b, &(-&a));
        let (s, t, h) = (eg.x, eg.y, eg.gcd);
        let u_num = &(&a * &s) + &(&b * &t).scale(&r2);
        let (u, u_rem) = u_num.div_rem(&h);
        if !u_rem.is_zero() {
            return None;
        }
        atans.push(two_atan_poly_over_surd(arena, &u, var, &r2));
        (a, b) = (s, t);
    }
    Some(match atans.len() {
        0 => arena.zero,
        1 => atans[0],
        _ => arena.add(&atans),
    })
}

// ═══════════════════════════════════════════════════════════════════════════
// log_to_real — main conversion function
// ═══════════════════════════════════════════════════════════════════════════

/// Convert an algebraic logarithmic term to real `ln + atan` expressions.
///
/// Given an irreducible factor `q(t)` of the Rothstein-Trager resultant
/// and the degree-1 PRS member `h(t, x)` (as a `GenPoly<RationalFn>`),
/// computes the real-valued contribution:
///
/// ```text
/// Σ_{α: q(α)=0} α · ln(h(α, x))
/// = Σ_j [ u_j · ln(A_j² + B_j²) + v_j · log_to_atan(A_j, B_j) ]
/// ```
///
/// where `(u_j, v_j)` are the real/imaginary parts of the roots of `q`
/// with `v_j > 0` (one from each conjugate pair).
///
/// Returns `None` if the roots of `q` cannot be found (degree ≥ 5 with
/// non-solvable Galois group) or if any intermediate computation fails.
pub(crate) fn log_to_real(
    arena: &mut Arena,
    var: ExprId,
    q: &Poly,
    h_prs: &GenPoly<RationalFn>,
) -> Option<Vec<ExprId>> {
    let q_deg = q.degree().unwrap_or(0);
    tracing::debug!(q_degree = q_deg, "log_to_real: starting conversion");

    if q_deg < 2 {
        tracing::debug!("log_to_real: q has degree < 2, should be handled as rational term");
        return None;
    }

    // ── Step 1: Solve q(t) = 0 for complex roots ──────────────────
    let t_sym = arena.symbol("__lrt_t");
    let q_expr = crate::poly::polybridge::poly_to_expr(arena, q, t_sym);
    let roots = crate::transforms::solve::solve(arena, q_expr, t_sym);

    tracing::debug!(n_roots = roots.len(), "log_to_real: solved q(t) = 0");

    if roots.is_empty() {
        tracing::debug!("log_to_real: no roots found (degree ≥ 5 non-solvable?)");
        return None;
    }

    // Roots `ρ·ζ` with `ρ` rational and `ζ` roots of unity (the residues
    // of `xᵏ/(xⁿ ± 1)`): every coefficient exactly in `ℚ(ζ_N)`.  For a
    // factor of degree ≤ 4 (radical roots, as SymPy's `roots`) first; for a
    // higher degree only where the radical roots below give a real form —
    // otherwise the answer stays a `RootSum`, as in SymPy (`∫ 1/(x⁷ + 1)`).
    if q_deg <= 4
        && let Some(terms) = cyclotomic_log_to_real(arena, var, q, h_prs)
    {
        return Some(terms);
    }

    // ── Step 2: Decompose roots into (u, v) = (Re, Im) pairs ──────
    let i_unit = arena.i_unit;
    let mut pairs: Vec<(ExprId, ExprId)> = Vec::new();
    let mut reals: Vec<ExprId> = Vec::new();
    let mut used = vec![false; roots.len()];

    for (idx, root) in roots.iter().enumerate() {
        if used[idx] {
            continue;
        }

        // The decomposition must be exact (no opaque `re(…)`/`im(…)` nodes):
        // nested Cardano/Ferrari radicals of complex numbers cannot be split
        // into real algebraic parts, and using them would only produce an
        // unusable `ln`/`atan` form.  Bail out so the caller can fall back.
        let parts = crate::base::complex::decompose(arena, root.value);
        if !parts.exact {
            tracing::debug!(idx, "log_to_real: root has no exact Re/Im decomposition");
            return None;
        }
        let (re_raw, im_raw) = (parts.re, parts.im);
        let u_val = crate::transforms::eval::eval(arena, re_raw);
        let v_val = crate::transforms::eval::eval(arena, im_raw);

        // Check imaginary part: cross-checked zero test (undecided: decline).
        let v_is_zero = crate::poly::algebraic::is_zero_checked(arena, v_val)?;

        if v_is_zero {
            // Real (irrational) root α: contributes α·ln|h(α, x)|.  Before
            // 0.25 it was skipped, so `∫ 1/(x³ − 2) dx` lost its
            // `ln|x − ∛2|` term and was rejected.
            tracing::trace!(idx, "log_to_real: real root");
            reals.push(u_val);
            used[idx] = true;
            continue;
        }

        // Sign test: keep roots with positive imaginary part.
        // A root listed before its conjugate may have the negative
        // imaginary part; the pair is represented by the conjugate.
        match crate::poly::algebraic::sign_checked(arena, v_val) {
            Some(1) => {
                tracing::trace!(idx, "log_to_real: found root with positive Im");
                pairs.push((u_val, v_val));
            }
            Some(-1) => {
                let neg_v = arena.neg(v_val);
                let neg_v = crate::transforms::eval::eval(arena, neg_v);
                pairs.push((u_val, neg_v));
            }
            _ => return None,
        }

        // Mark this root as used.
        used[idx] = true;

        // Find and mark the conjugate root (same Re, opposite Im).
        // Conjugate matching uses f64 sum-to-zero heuristic.
        let v_f64 = crate::transforms::evalf::eval_const_f64(arena, v_val)?;
        for j in (idx + 1)..roots.len() {
            if used[j] {
                continue;
            }
            let (_, im_j_raw) = crate::base::complex::as_real_imag(arena, roots[j].value);
            let v_j = crate::transforms::eval::eval(arena, im_j_raw);
            let v_j_f64 = crate::transforms::evalf::eval_const_f64(arena, v_j);
            if let Some(vj) = v_j_f64
                && (vj + v_f64).abs() < 1e-10
            {
                // Conjugate found.
                tracing::trace!(j, "log_to_real: conjugate root found and marked");
                used[j] = true;
                break;
            }
        }
    }

    tracing::debug!(
        n_pairs = pairs.len(),
        "log_to_real: conjugate pairs identified"
    );

    if pairs.is_empty() && reals.is_empty() {
        tracing::debug!("log_to_real: no roots to convert");
        return None;
    }
    if pairs.len() * 2 + reals.len() != q_deg {
        tracing::debug!(
            n_pairs = pairs.len(),
            n_reals = reals.len(),
            "log_to_real: roots not all classified"
        );
        return None;
    }

    // ── Step 3: Extract h(t,x) coefficients ────────────────────────
    // h(t,x) is degree 1 in x (monic): h = 1·x + h₀(t)
    // where h₀(t) is a RationalFn.
    if h_prs.degree() != Some(1) {
        tracing::debug!(
            h_degree = ?h_prs.degree(),
            "log_to_real: h(t,x) is not degree 1 in x, cannot proceed"
        );
        return None;
    }
    let h_coeff_1: RationalFn = h_prs.coeff(1);
    let h_coeff_0: RationalFn = h_prs.coeff(0);

    // Convert to arena expressions in t_sym.
    let h1_expr = crate::poly::polybridge::ratfn_to_expr(arena, &h_coeff_1, t_sym);
    let h0_expr = crate::poly::polybridge::ratfn_to_expr(arena, &h_coeff_0, t_sym);

    tracing::trace!("log_to_real: h(t,x) coefficients converted to arena");

    // ── Step 4: For each (u_j, v_j), build ln + atan terms ────────
    let mut result_terms: Vec<ExprId> = Vec::new();

    for &alpha in &reals {
        let h1_at = crate::transforms::subs::subs(arena, h1_expr, t_sym, alpha);
        let h1_at = crate::transforms::eval::eval(arena, h1_at);
        let h0_at = crate::transforms::subs::subs(arena, h0_expr, t_sym, alpha);
        let h0_at = crate::transforms::eval::eval(arena, h0_at);
        let h1_x = arena.mul(&[h1_at, var]);
        let h_at = arena.add(&[h1_x, h0_at]);
        let abs_h = arena.abs(h_at);
        let ln_h = arena.ln(abs_h);
        result_terms.push(arena.mul(&[alpha, ln_h]));
    }

    for (pair_idx, (u_j, v_j)) in pairs.iter().enumerate() {
        tracing::debug!(pair_idx, "log_to_real: processing conjugate pair");

        // Construct t_val = u_j + I·v_j
        let i_v = arena.mul(&[i_unit, *v_j]);
        let t_val = arena.add(&[*u_j, i_v]);
        let t_val = crate::transforms::eval::eval(arena, t_val);

        // Evaluate h coefficients at t = t_val
        let h1_at = crate::transforms::subs::subs(arena, h1_expr, t_sym, t_val);
        let h1_at = crate::transforms::eval::eval(arena, h1_at);
        let h0_at = crate::transforms::subs::subs(arena, h0_expr, t_sym, t_val);
        let h0_at = crate::transforms::eval::eval(arena, h0_at);

        // h(α, x) = h1_at · x + h0_at
        // Separate Re/Im for each coefficient.
        let (re_h1, im_h1) = crate::base::complex::as_real_imag(arena, h1_at);
        let re_h1 = crate::transforms::eval::eval(arena, re_h1);
        let im_h1 = crate::transforms::eval::eval(arena, im_h1);

        let (re_h0, im_h0) = crate::base::complex::as_real_imag(arena, h0_at);
        let re_h0 = crate::transforms::eval::eval(arena, re_h0);
        let im_h0 = crate::transforms::eval::eval(arena, im_h0);

        // A(x) = Re(h1)·x + Re(h0)
        // B(x) = Im(h1)·x + Im(h0)
        tracing::trace!(
            pair_idx,
            "log_to_real: A(x) = Re(h1)·x + Re(h0), B(x) = Im(h1)·x + Im(h0)"
        );

        // ── Build ln(A² + B²) ──────────────────────────────────────
        //
        // The norm |h(α,x)|² = h(α,x)·h(ᾱ,x) is always a polynomial in x
        // with rational coefficients (the norm over the extension field).
        //
        // For h(t,x) = x + c(t):  |h|² = (x + c(α))(x + c(ᾱ))
        //   = x² + (c(α)+c(ᾱ))·x + c(α)·c(ᾱ)
        //   = x² + 2·Re(c(α))·x + |c(α)|²
        //
        // We compute this symbolically: A² + B² where A and B are arena
        // expressions, then expand + eval to simplify.
        let re_h1_x = arena.mul(&[re_h1, var]);
        let a_expr = arena.add(&[re_h1_x, re_h0]);
        let im_h1_x = arena.mul(&[im_h1, var]);
        let b_expr = arena.add(&[im_h1_x, im_h0]);

        let a_sq = arena.mul(&[a_expr, a_expr]);
        let b_sq = arena.mul(&[b_expr, b_expr]);
        let norm_sq = arena.add(&[a_sq, b_sq]);

        // Expand and eval to simplify: (x+1/2)² + (√3/2)² → x²+x+1
        let norm_sq = crate::transforms::expand::expand(arena, norm_sq);
        let norm_sq = crate::transforms::eval::eval(arena, norm_sq);

        // Apply radical simplification (powdenest + powsimp_base) to
        // clean up terms like (√3)² → 3.
        let norm_sq = crate::simplify::powsimp::powdenest(arena, norm_sq);
        let norm_sq = crate::simplify::powsimp::powsimp_base(arena, norm_sq);
        let norm_sq = crate::transforms::eval::eval(arena, norm_sq);

        tracing::trace!(pair_idx, "log_to_real: norm |h|² computed and simplified");

        // ln term: u_j · ln(|h|²)
        // ── Check if u_j = 0 (pure imaginary root) ────────────────
        // (An undecided zero test declines: dropping or keeping a term on a
        // guess would give a wrong antiderivative.)
        let u_is_zero = crate::poly::algebraic::is_zero_checked(arena, *u_j)?;

        // ── Check if B(x) = 0 (imaginary part vanishes) ───────────
        let im_h1_is_zero = crate::poly::algebraic::is_zero_checked(arena, im_h1)?;
        let im_h0_is_zero = crate::poly::algebraic::is_zero_checked(arena, im_h0)?;
        let b_is_zero = im_h1_is_zero && im_h0_is_zero;

        if u_is_zero && b_is_zero {
            // Both u ≈ 0 and B ≈ 0: this pair contributes nothing.
            // This happens for integrands like x/(x⁴+x²+1) where the
            // Rothstein-Trager roots are pure imaginary and h evaluates
            // to a real expression.  Return None to let the algebraic
            // remainder fallback handle this case (via apart).
            tracing::debug!(
                pair_idx,
                "log_to_real: u ≈ 0 and B ≈ 0 — pair contributes nothing, bailing"
            );
            return None;
        }

        // Build ln term: u_j · ln(|h|²)  (skip if u_j ≈ 0)
        if !u_is_zero {
            let ln_norm = arena.ln(norm_sq);
            let ln_term = arena.mul(&[*u_j, ln_norm]);
            result_terms.push(ln_term);
        }

        if b_is_zero {
            // B ≈ 0 but u ≠ 0: only the ln term contributes (already pushed above).
            tracing::debug!(pair_idx, "log_to_real: B ≈ 0, ln-only term (no atan)");
        } else {
            // ── Build atan via log_to_atan_deg1 ────────────────────────
            let atan_result = log_to_atan_deg1(arena, var, re_h1, re_h0, im_h1, im_h0);

            match atan_result {
                Some(atan_expr) => {
                    // atan term: v_j · log_to_atan(A, B)
                    // (log_to_atan already includes the factor of 2)
                    let atan_term = arena.mul(&[*v_j, atan_expr]);
                    result_terms.push(atan_term);
                    tracing::debug!(pair_idx, "log_to_real: emitted atan term");
                }
                None => {
                    tracing::debug!(
                        pair_idx,
                        "log_to_real: log_to_atan_deg1 failed, skipping pair"
                    );
                    return None;
                }
            }
        }
    }

    if result_terms.is_empty() {
        None
    } else {
        tracing::debug!(
            n_terms = result_terms.len(),
            "log_to_real: conversion complete"
        );
        // The same terms exactly in `ℚ(ζ_N)` when the roots allow (above).
        if q_deg > 4
            && let Some(terms) = cyclotomic_log_to_real(arena, var, q, h_prs)
        {
            return Some(terms);
        }
        Some(result_terms)
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Roots of unity: log_to_real in ℚ(ζ_N)
// ═══════════════════════════════════════════════════════════════════════════

/// Largest `M` tried for `t^M ≡ c (mod q)`.
const MAX_CYCLOTOMIC_ORDER: usize = 96;

/// The cyclotomic field `ℚ(ζ) = ℚ[z]/(Φ_N)`, `ζ = e^{2πi/N}`, `4 | N` (so
/// that `i = ζ^{N/4}` is in it).  Elements are polynomials of degree
/// `< φ(N)`, which is canonical (`Φ_N` is irreducible).
struct Cyclotomic {
    n: usize,
    phi: Poly,
}

impl Cyclotomic {
    /// `Φ_N = (x^N − 1)/∏_{d | N, d < N} Φ_d`, the divisors in increasing
    /// order (each `Φ_d` from the smaller ones).
    fn new(n: usize) -> Option<Cyclotomic> {
        let divisors: Vec<usize> = (1..=n).filter(|d| n.is_multiple_of(*d)).collect();
        let mut phis: Vec<(usize, Poly)> = Vec::with_capacity(divisors.len());
        for &d in &divisors {
            let mut p = monomial(d, Q::from_integer(BigInt::from(1)));
            p = &p - &Poly::constant(Q::from_integer(BigInt::from(1)));
            for (e, pe) in &phis {
                if d % e == 0 {
                    let (quot, rem) = p.div_rem(pe);
                    if !rem.is_zero() {
                        return None;
                    }
                    p = quot;
                }
            }
            phis.push((d, p));
        }
        let (_, phi) = phis.pop()?;
        Some(Cyclotomic { n, phi })
    }

    fn reduce(&self, p: &Poly) -> Poly {
        p.rem(&self.phi)
    }

    fn mul(&self, a: &Poly, b: &Poly) -> Poly {
        self.reduce(&(a * b))
    }

    /// Complex conjugation `ζ^k ↦ ζ^{N−k}`.
    fn conj(&self, a: &Poly) -> Poly {
        let mut c = vec![Q::from_integer(BigInt::from(0)); self.n];
        for (k, v) in a.coeffs().iter().enumerate() {
            c[(self.n - k) % self.n] += v;
        }
        self.reduce(&Poly::from_coeffs(c))
    }

    /// Real and imaginary parts, `(a + ā)/2` and `−i·(a − ā)/2`.
    fn re_im(&self, a: &Poly) -> [Poly; 2] {
        let half = Q::new(BigInt::from(1), BigInt::from(2));
        let ab = self.conj(a);
        let re = (a + &ab).scale(&half);
        let neg_i = monomial(3 * self.n / 4, -half);
        let im = self.mul(&(a - &ab), &neg_i);
        [re, im]
    }

    /// The `f64` value of the real part of `a` (for a sign that only
    /// chooses between two equal forms).
    fn approx(&self, a: &Poly) -> f64 {
        let mut s = 0.0;
        for (k, v) in a.coeffs().iter().enumerate() {
            let angle = 2.0 * std::f64::consts::PI * k as f64 / self.n as f64;
            s += num_traits::ToPrimitive::to_f64(v).unwrap_or(0.0) * angle.cos();
        }
        s
    }

    /// A real element `r = Σ aₖ·cos(2πk/N)` over the basis
    /// `k = 0 … φ(N)/2 − 1` of the real subfield (`2cos(2πk/N) = ζ^k + ζ^{−k}`
    /// is a monic polynomial of degree `k` in `2cos(2π/N)`, of degree
    /// `φ(N)/2`): the coefficients `aₖ`, by elimination over ℚ.
    fn cos_coordinates(&self, r: &Poly) -> Option<Vec<Q>> {
        let zero = Q::from_integer(BigInt::from(0));
        let dim = self.phi.degree()?;
        let d = dim / 2;
        let half = Q::new(BigInt::from(1), BigInt::from(2));
        // Columns: the basis vectors, then r.
        let mut cols: Vec<Vec<Q>> = Vec::with_capacity(d + 1);
        for k in 0..d {
            let b = &monomial(k, half.clone()) + &monomial((self.n - k) % self.n, half.clone());
            let b = self.reduce(&b);
            cols.push((0..dim).map(|i| b.coeff(i)).collect());
        }
        cols.push((0..dim).map(|i| r.coeff(i)).collect());
        // Row-reduce the dim × (d+1) system.
        let mut rows: Vec<Vec<Q>> = (0..dim)
            .map(|i| cols.iter().map(|c| c[i].clone()).collect())
            .collect();
        let mut piv_row = 0usize;
        for col in 0..d {
            let p = (piv_row..dim).find(|&i| rows[i][col] != zero)?;
            rows.swap(piv_row, p);
            let inv = rows[piv_row][col].recip();
            for v in rows[piv_row].iter_mut() {
                *v *= &inv;
            }
            let pivot = rows[piv_row].clone();
            for (i, row) in rows.iter_mut().enumerate() {
                if i != piv_row && row[col] != zero {
                    let f = row[col].clone();
                    for (x, pv) in row.iter_mut().zip(&pivot) {
                        *x -= &f * pv;
                    }
                }
            }
            piv_row += 1;
        }
        // Consistent (r real): the rows below the rank are zero.
        if rows[piv_row..].iter().any(|row| row[d] != zero) {
            return None;
        }
        Some((0..d).map(|k| rows[k][d].clone()).collect())
    }

    /// `r` real as an expression: `Σ aₖ·cos(2πk/N)`, each `cos` of an
    /// angle above `π/4` written as the `sin` of its complement, or the same
    /// with the radical forms of `cos(π/5)`, `cos(π/8)`, `cos(π/10)`,
    /// `cos(π/12)`, … — whichever is smaller after `eval` (which folds the
    /// multiples of `π/4` and `π/6`).
    fn real_expr(&self, arena: &mut Arena, r: &Poly) -> Option<ExprId> {
        let coords = self.cos_coordinates(r)?;
        let mut trig: Vec<ExprId> = Vec::new();
        let mut radical: Vec<ExprId> = Vec::new();
        let mut any_radical = false;
        for (k, a) in coords.iter().enumerate() {
            if a == &Q::from_integer(BigInt::from(0)) {
                continue;
            }
            let a_id = arena.num_ratio(a.clone());
            if k == 0 {
                trig.push(a_id);
                radical.push(a_id);
                continue;
            }
            // cos(π·angle), angle = 2k/N ∈ (0, 1/2)
            let angle = Q::new(BigInt::from(2 * k), BigInt::from(self.n));
            let quarter = Q::new(BigInt::from(1), BigInt::from(4));
            let atom = if angle > quarter {
                let comp = Q::new(BigInt::from(1), BigInt::from(2)) - &angle;
                let c = arena.num_ratio(comp);
                let arg = arena.mul(&[c, arena.pi()]);
                arena.sin(arg)
            } else {
                let c = arena.num_ratio(angle.clone());
                let arg = arena.mul(&[c, arena.pi()]);
                arena.cos(arg)
            };
            trig.push(arena.mul(&[a_id, atom]));
            let rad = match cos_radical(arena, &angle) {
                Some(v) => {
                    any_radical = true;
                    v
                }
                None => atom,
            };
            radical.push(arena.mul(&[a_id, rad]));
        }
        let build = |arena: &mut Arena, parts: &[ExprId]| -> ExprId {
            let s = match parts.len() {
                0 => arena.zero,
                1 => parts[0],
                _ => arena.add(parts),
            };
            crate::transforms::eval::eval(arena, s)
        };
        let t = build(arena, &trig);
        if !any_radical {
            return Some(t);
        }
        let rd = build(arena, &radical);
        let cap = 10_000;
        let size = |arena: &Arena, e| crate::transforms::pattern::tree_size_capped(arena, e, cap);
        Some(if size(arena, rd) < size(arena, t) {
            rd
        } else {
            t
        })
    }
}

/// `ℚ(ζ)[ρ]/(ρ^L − κ)` for the real `ρ = κ^{1/L} > 0` (`κ` rational, `L`
/// the least power of `ρ` that is rational): an element is `Σ ρ^r·aᵣ`
/// (`r < L`, `aᵣ ∈ ℚ(ζ)`).  A zero representation is the number 0; for
/// `L > 1` the converse can fail (`ρ⁴ = √2 = ζ₈ + ζ₈⁻¹`), so
/// [`is_zero`](Self::is_zero) confirms a nonzero one numerically.
struct RadCyc {
    cyc: Cyclotomic,
    l: usize,
    kappa: Q,
}

impl RadCyc {
    fn zero_elem(&self) -> Vec<Poly> {
        vec![Poly::from_coeffs(Vec::new()); self.l]
    }

    /// `p(ρ·ζ^e)`.
    fn eval_at(&self, p: &Poly, e: usize) -> Vec<Poly> {
        let n = self.cyc.n;
        let zero = Q::from_integer(BigInt::from(0));
        let mut comps = vec![vec![zero; n]; self.l];
        let mut kpow = Q::from_integer(BigInt::from(1));
        for (m, v) in p.coeffs().iter().enumerate() {
            if m > 0 && m % self.l == 0 {
                kpow *= &self.kappa;
            }
            comps[m % self.l][(e * m) % n] += v * &kpow;
        }
        comps
            .into_iter()
            .map(|c| self.cyc.reduce(&Poly::from_coeffs(c)))
            .collect()
    }

    /// `ρ·ζ^e`
    fn root(&self, e: usize) -> Vec<Poly> {
        let mut out = self.zero_elem();
        if self.l == 1 {
            out[0] = self.cyc.reduce(&monomial(e, self.kappa.clone()));
        } else {
            out[1] = self
                .cyc
                .reduce(&monomial(e, Q::from_integer(BigInt::from(1))));
        }
        out
    }

    fn mul(&self, a: &[Poly], b: &[Poly]) -> Vec<Poly> {
        let mut out = self.zero_elem();
        for (i, x) in a.iter().enumerate() {
            for (j, y) in b.iter().enumerate() {
                if x.is_zero() || y.is_zero() {
                    continue;
                }
                let mut t = self.cyc.mul(x, y);
                if i + j >= self.l {
                    t = t.scale(&self.kappa);
                }
                let k = (i + j) % self.l;
                out[k] = &out[k] + &t;
            }
        }
        out
    }

    fn conj(&self, a: &[Poly]) -> Vec<Poly> {
        a.iter().map(|x| self.cyc.conj(x)).collect()
    }

    fn scale(a: &[Poly], c: &Q) -> Vec<Poly> {
        a.iter().map(|x| x.scale(c)).collect()
    }

    fn re_im(&self, a: &[Poly]) -> [Vec<Poly>; 2] {
        let parts: Vec<[Poly; 2]> = a.iter().map(|x| self.cyc.re_im(x)).collect();
        [
            parts.iter().map(|[re, _]| re.clone()).collect(),
            parts.iter().map(|[_, im]| im.clone()).collect(),
        ]
    }

    fn rho(&self) -> f64 {
        num_traits::ToPrimitive::to_f64(&self.kappa)
            .unwrap_or(f64::NAN)
            .powf(1.0 / self.l as f64)
    }

    /// The value of a real element (numerically).
    fn approx(&self, a: &[Poly]) -> f64 {
        let rho = self.rho();
        a.iter()
            .enumerate()
            .map(|(r, x)| rho.powi(r as i32) * self.cyc.approx(x))
            .sum()
    }

    /// `Some(true)` for a zero representation; `Some(false)` when the value
    /// is clearly nonzero; `None` when a nonzero representation (`L > 1`)
    /// has a value too small to tell.
    fn is_zero(&self, a: &[Poly]) -> Option<bool> {
        if a.iter().all(Poly::is_zero) {
            return Some(true);
        }
        if self.l == 1 {
            return Some(false);
        }
        // The size of the terms ρ^r·c·ζ^k (a cancellation below 10⁻⁹ of it
        // is left undecided).
        let rho = self.rho();
        let scale: f64 = a
            .iter()
            .enumerate()
            .flat_map(|(r, x)| {
                x.coeffs().iter().map(move |c| {
                    rho.powi(r as i32)
                        * num_traits::ToPrimitive::to_f64(c).unwrap_or(f64::NAN).abs()
                })
            })
            .fold(0.0, f64::max);
        (self.approx(a).abs() > 1e-9 * scale).then_some(false)
    }

    /// `Σ ρ^r·aᵣ` as an expression (`ρ^r = κ^{r/L}`).
    fn real_expr(&self, arena: &mut Arena, a: &[Poly]) -> Option<ExprId> {
        let mut parts = Vec::new();
        for (r, x) in a.iter().enumerate() {
            if x.is_zero() {
                continue;
            }
            let v = self.cyc.real_expr(arena, x)?;
            if r == 0 {
                parts.push(v);
            } else {
                let k = arena.num_ratio(self.kappa.clone());
                let ex = arena.rational(r as i64, self.l as i64);
                let rho_r = arena.pow(k, ex);
                parts.push(arena.mul(&[rho_r, v]));
            }
        }
        let s = match parts.len() {
            0 => arena.zero,
            1 => parts[0],
            _ => arena.add(&parts),
        };
        Some(crate::transforms::eval::eval(arena, s))
    }
}

/// `h₀/h₁ mod q` over ℚ (`q` irreducible, so `ℚ[t]/(q)` is a field).
fn ratfn_quotient_mod(h0: &RationalFn, h1: &RationalFn, q: &Poly) -> Option<Poly> {
    let num = (h0.numer() * h1.denom()).rem(q);
    let den = (h0.denom() * h1.numer()).rem(q);
    if den.is_zero() {
        return None;
    }
    let eg = Poly::extended_gcd(&den, q);
    if eg.gcd.degree() != Some(0) {
        return None;
    }
    let inv = eg.x.scale(&eg.gcd.coeff(0).recip());
    Some((&num * &inv).rem(q))
}

/// `c·x^e` as a [`Poly`].
fn monomial(e: usize, c: Q) -> Poly {
    let mut v = vec![Q::from_integer(BigInt::from(0)); e + 1];
    v[e] = c;
    Poly::from_coeffs(v)
}

/// `cos(π·angle)` in radicals for the angles of denominator 5, 8, 10 and
/// 12 in `(0, 1/2)`; `None` otherwise.
fn cos_radical(arena: &mut Arena, angle: &Q) -> Option<ExprId> {
    let (p, d) = (
        num_traits::ToPrimitive::to_i64(angle.numer())?,
        num_traits::ToPrimitive::to_i64(angle.denom())?,
    );
    let half = arena.rational(1, 2);
    let sqrt = |arena: &mut Arena, e: ExprId| arena.pow(e, half);
    // (a + b·√s)/c or √(a + b·√s)/c
    let (a, b, s, c, nested) = match (p, d) {
        (1, 5) => (1, 1, 5, 4, false),
        (2, 5) => (-1, 1, 5, 4, false),
        (1, 8) => (2, 1, 2, 2, true),
        (3, 8) => (2, -1, 2, 2, true),
        (1, 10) => (10, 2, 5, 4, true),
        (3, 10) => (10, -2, 5, 4, true),
        (1, 12) | (5, 12) => {
            // (√6 ± √2)/4
            let six = arena.int(6);
            let two = arena.int(2);
            let (r6, r2) = (sqrt(arena, six), sqrt(arena, two));
            let r2 = if p == 1 { r2 } else { arena.neg(r2) };
            let sum = arena.add(&[r6, r2]);
            let quarter = arena.rational(1, 4);
            return Some(arena.mul(&[quarter, sum]));
        }
        _ => return None,
    };
    let s_id = arena.int(s);
    let rs = sqrt(arena, s_id);
    let b_id = arena.int(b);
    let a_id = arena.int(a);
    let brs = arena.mul(&[b_id, rs]);
    let mut inner = arena.add(&[a_id, brs]);
    if nested {
        inner = sqrt(arena, inner);
    }
    let inv_c = arena.rational(1, c);
    Some(arena.mul(&[inv_c, inner]))
}

/// [`log_to_real`] exactly in a cyclotomic field when the roots of `q` are
/// `ρ·ζ` with `ρ > 0` a radical and `ζ` roots of unity — `t^M ≡ c (mod q)`
/// for a rational `c`, `ρ = |c|^{1/M}` — as for the Rothstein–Trager
/// residues of `xᵏ/(xⁿ + a)`.  The roots are the `ρ·ζ^e`, `ζ = e^{2πi/N}`,
/// at which `q` vanishes in `ℚ(ζ)[ρ]` ([`RadCyc`]); the log argument
/// `x + g(α)`, `g = h₀/h₁ mod q` (a constant multiple of `h(α, x)`), its
/// real and imaginary parts, `|g|²` and `α = u + iv` are all elements of
/// `ℚ(ζ)[ρ]`, written as `Σ ρ^r·(…)` over the basis `cos(2πk/N)` of the
/// real subfield (or in radicals, whichever is smaller).  Before 0.32 they were the
/// expansions of powers of `cos(π/n) + i·sin(π/n)` from the solver's
/// roots, de Moivre not applied (`∫ x⁶/(x⁸ + 1)` printed 18 KB).  The
/// terms are those of [`log_to_real`]: `α·ln|x + g|` for a real root and
/// `u·ln(x² + 2·Re(g)·x + |g|²) + 2v·atan((x + Re g)/Im g)` for each pair
/// `α, ᾱ` with `v > 0`.  `None` when `q` is not of this form.
fn cyclotomic_log_to_real(
    arena: &mut Arena,
    var: ExprId,
    q: &Poly,
    h_prs: &GenPoly<RationalFn>,
) -> Option<Vec<ExprId>> {
    use num_traits::Signed;
    let deg = q.degree()?;
    if deg < 3 || h_prs.degree() != Some(1) {
        return None;
    }
    // t^M mod q for M = 1, 2, …: the first constant one.
    let t = Poly::x();
    let mut p = t.rem(q);
    let mut found: Option<(usize, Q)> = None;
    for m in 2..=MAX_CYCLOTOMIC_ORDER.min(8 * deg) {
        p = (&p * &t).rem(q);
        if p.degree() == Some(0) && m >= deg {
            found = Some((m, p.coeff(0)));
            break;
        }
    }
    let (m, c) = found?;
    // ρ = |c|^{1/M}; L = the least divisor of M with ρ^L rational.
    let exact_root = |z: &BigInt, k: u32| -> Option<BigInt> {
        let r = z.nth_root(k);
        (num_traits::Pow::pow(&r, k) == *z).then_some(r)
    };
    let abs_c = c.abs();
    let (l, kappa) = (1..=m).filter(|l| m.is_multiple_of(*l)).find_map(|l| {
        let k = u32::try_from(m / l).ok()?;
        Some((
            l,
            Q::new(exact_root(abs_c.numer(), k)?, exact_root(abs_c.denom(), k)?),
        ))
    })?;
    let m_full = if c.is_positive() { m } else { 2 * m };
    let n = num_integer::lcm(m_full, 4);
    let rc = RadCyc {
        cyc: Cyclotomic::new(n)?,
        l,
        kappa,
    };
    let step = n / m_full;
    let candidates = (0..m_full)
        .filter(|j| c.is_positive() || j % 2 == 1)
        .map(|j| j * step);
    // All the roots of t^M = c when q is that binomial; otherwise those at
    // which q vanishes (exactly: a zero representation is a zero).
    let exps: Vec<usize> = if deg == m {
        candidates.collect()
    } else {
        candidates
            .filter(|&e| rc.eval_at(q, e).iter().all(Poly::is_zero))
            .collect()
    };
    if exps.len() != deg {
        return None;
    }
    let g_t = ratfn_quotient_mod(&h_prs.coeff(0), &h_prs.coeff(1), q)?;
    let mut terms: Vec<ExprId> = Vec::new();
    let two = arena.int(2);
    for &e in &exps {
        // One root of each conjugate pair (v > 0 ⇔ 0 < e < N/2) and the real ones.
        if e > n / 2 {
            continue;
        }
        let g = rc.eval_at(&g_t, e);
        let [re_g, im_g] = rc.re_im(&g);
        let [u, v] = rc.re_im(&rc.root(e));
        let re_g_id = rc.real_expr(arena, &re_g)?;
        let x_plus = arena.add(&[var, re_g_id]);
        let x_plus = crate::transforms::eval::eval(arena, x_plus);
        let im_zero = rc.is_zero(&im_g)?;
        if e == 0 || 2 * e == n {
            // Real root α = ±ρ: α·ln|x + g| (g real).
            if !im_zero {
                return None;
            }
            let a_id = rc.real_expr(arena, &u)?;
            let abs = arena.abs(x_plus);
            let ln = arena.ln(abs);
            terms.push(arena.mul(&[a_id, ln]));
            continue;
        }
        let u_zero = rc.is_zero(&u)?;
        if u_zero && im_zero {
            return None;
        }
        if !u_zero {
            // x² + 2·Re(g)·x + |g|²
            let norm_g = rc.mul(&g, &rc.conj(&g));
            let [norm_g, _] = rc.re_im(&norm_g);
            let n_id = rc.real_expr(arena, &norm_g)?;
            let lin_c = rc.real_expr(
                arena,
                &RadCyc::scale(&re_g, &Q::from_integer(BigInt::from(2))),
            )?;
            let x2 = arena.pow(var, two);
            let lin = arena.mul(&[lin_c, var]);
            let norm = arena.add(&[x2, lin, n_id]);
            let norm = crate::transforms::eval::eval(arena, norm);
            let u_id = rc.real_expr(arena, &u)?;
            let ln = arena.ln(norm);
            terms.push(arena.mul(&[u_id, ln]));
        }
        if !im_zero {
            // atan is odd: a negative Im g goes to the coefficient.
            let flip = rc.approx(&im_g) < 0.0;
            let sgn = Q::from_integer(BigInt::from(if flip { -1 } else { 1 }));
            let b_id = rc.real_expr(arena, &RadCyc::scale(&im_g, &sgn))?;
            let ratio = arena.div(x_plus, b_id);
            let ratio = crate::transforms::eval::eval(arena, ratio);
            let atan = arena.atan(ratio);
            let two_v = RadCyc::scale(&v, &(Q::from_integer(BigInt::from(2)) * &sgn));
            let two_v = rc.real_expr(arena, &two_v)?;
            terms.push(arena.mul(&[two_v, atan]));
        }
    }
    tracing::debug!(
        n,
        deg,
        n_terms = terms.len(),
        "log_to_real: roots of unity, exact in Q(zeta_N)"
    );
    (!terms.is_empty()).then_some(terms)
}

// ═══════════════════════════════════════════════════════════════════════════
// Vieta's formulas — symmetric functions of roots without root-finding
// ═══════════════════════════════════════════════════════════════════════════

/// Extract the elementary symmetric polynomials `e_k` from a polynomial's
/// coefficients using Vieta's formulas.
///
/// For a monic polynomial `t^n + a_{n-1}·t^{n-1} + ... + a_1·t + a_0`,
/// the elementary symmetric polynomials of its roots are:
///
/// ```text
/// e_1 = -a_{n-1}          (sum of roots)
/// e_2 =  a_{n-2}          (sum of products of pairs)
/// e_k = (-1)^k · a_{n-k}  (k-th elementary symmetric polynomial)
/// ```
///
/// Returns `e_1, e_2, ..., e_n` as `Ratio<BigInt>` values.
/// The input polynomial must be monic (leading coefficient = 1).
///
/// Returns `None` if the polynomial is zero, constant, or not monic.
pub(crate) fn vieta_elementary_symmetric(poly: &Poly) -> Option<Vec<Q>> {
    let n = poly.degree()?;
    if n == 0 {
        return None;
    }

    // Check monic (leading coefficient = 1).
    let lc = poly.leading_coeff()?;
    if !num_traits::One::is_one(lc) {
        tracing::trace!("vieta_elementary_symmetric: polynomial is not monic");
        return None;
    }

    let one: Q = Ratio::from_integer(BigInt::from(1));
    let neg_one: Q = Ratio::from_integer(BigInt::from(-1));

    let mut e = Vec::with_capacity(n);
    for k in 1..=n {
        // e_k = (-1)^k · a_{n-k}
        let coeff = poly.coeff(n - k);
        let sign = if k % 2 == 0 { &one } else { &neg_one };
        e.push(sign * &coeff);
    }

    Some(e)
}

/// Compute the power sum `p_k = Σ_{i=1}^{n} α_i^k` using Newton's identities,
/// given the elementary symmetric polynomials `e_1, ..., e_n`.
///
/// Newton's identities:
/// ```text
/// p_1 = e_1
/// p_2 = e_1·p_1 - 2·e_2
/// p_k = Σ_{i=1}^{k-1} (-1)^{i-1}·e_i·p_{k-i} + (-1)^{k-1}·k·e_k   (k ≤ n)
/// p_k = Σ_{i=1}^{n}   (-1)^{i-1}·e_i·p_{k-i}                        (k > n)
/// ```
///
/// Returns `p_1, p_2, ..., p_max_k`.
pub(crate) fn vieta_power_sums(elementary: &[Q], max_k: usize) -> Vec<Q> {
    let n = elementary.len(); // degree of the polynomial
    let one: Q = Ratio::from_integer(BigInt::from(1));
    let neg_one: Q = Ratio::from_integer(BigInt::from(-1));
    let mut p: Vec<Q> = Vec::with_capacity(max_k);

    for k in 1..=max_k {
        let mut pk = Ratio::from_integer(BigInt::from(0));

        let upper = if k <= n { k - 1 } else { n };
        for i in 1..=upper {
            // (-1)^{i-1} · e_i · p_{k-i}
            let sign = if (i - 1) % 2 == 0 { &one } else { &neg_one };
            let e_i = &elementary[i - 1];
            let p_km = if k - i >= 1 {
                &p[k - i - 1] // p_{k-i} (0-indexed)
            } else {
                // k - i == 0 → this shouldn't happen since i ≤ k-1
                continue;
            };
            pk += sign * e_i * p_km;
        }

        // For k ≤ n: add the (-1)^{k-1} · k · e_k term
        if k <= n {
            let sign = if (k - 1) % 2 == 0 { &one } else { &neg_one };
            let k_rat = Ratio::from_integer(BigInt::from(k));
            pk += sign * &k_rat * &elementary[k - 1];
        }

        p.push(pk);
    }

    p
}

/// Try to evaluate `RootSum(poly, body, sumvar)` when the body is a
/// polynomial in `sumvar` (no other variables), using Vieta's formulas
/// and Newton's identities.
///
/// For `body = c_m·t^m + ... + c_1·t + c_0`, the sum over all roots is:
/// ```text
/// Σ body(α_i) = c_m·p_m + ... + c_1·p_1 + n·c_0
/// ```
/// where `p_k = Σ α_i^k` (power sums) and `n` is the polynomial degree.
///
/// Returns `Some(rational_value)` if the body is a polynomial in sumvar
/// with rational coefficients and no other free variables.
/// Returns `None` otherwise.
pub(crate) fn vieta_rootsum_poly_body(
    arena: &Arena,
    poly_id: ExprId,
    body_id: ExprId,
    sumvar_id: ExprId,
) -> Option<Q> {
    // Extract the polynomial as Poly.
    let poly = crate::poly::polybridge::expr_to_poly(arena, poly_id, sumvar_id)?;
    let n = poly.degree()?;

    // Make monic for Vieta's formulas.
    let monic = poly.make_monic();

    // Extract body as Poly in sumvar.
    let body_poly = crate::poly::polybridge::expr_to_poly(arena, body_id, sumvar_id)?;
    let body_deg = body_poly.degree().unwrap_or(0);

    tracing::debug!(
        poly_degree = n,
        body_degree = body_deg,
        "vieta_rootsum_poly_body: attempting Vieta evaluation"
    );

    // Get elementary symmetric polynomials.
    let elementary = vieta_elementary_symmetric(&monic)?;

    // Compute power sums up to the body degree.
    let power_sums = vieta_power_sums(&elementary, body_deg);

    // Evaluate: Σ body(α_i) = Σ_k c_k · p_k + n · c_0
    let n_rat = Ratio::from_integer(BigInt::from(n));
    let mut result = &n_rat * body_poly.coeff(0); // n · c_0

    for k in 1..=body_deg {
        let c_k = body_poly.coeff(k);
        if !c_k.is_zero() && k <= power_sums.len() {
            result += &c_k * &power_sums[k - 1];
        }
    }

    tracing::debug!(%result, "vieta_rootsum_poly_body: computed via Vieta");
    Some(result)
}

// ═══════════════════════════════════════════════════════════════════════════
// Numeric RootSum with a rational body: exact trace
// ═══════════════════════════════════════════════════════════════════════════

/// A Gaussian rational `re + im·i`: the coefficients of a numeric
/// `RootSum` once the parameters of an antiderivative are bound to the
/// complex values of the integrator's self-check.
#[derive(Clone, Debug, PartialEq)]
struct GaussQ {
    re: Q,
    im: Q,
}

impl GaussQ {
    fn real(re: Q) -> Self {
        GaussQ {
            re,
            im: Q::from_integer(BigInt::from(0)),
        }
    }

    fn from_usize(n: usize) -> Self {
        GaussQ::real(Q::from_integer(BigInt::from(n)))
    }
}

impl Ring for GaussQ {
    fn zero() -> Self {
        GaussQ::real(Q::from_integer(BigInt::from(0)))
    }
    fn one() -> Self {
        GaussQ::real(Q::from_integer(BigInt::from(1)))
    }
    fn is_zero(&self) -> bool {
        num_traits::Zero::is_zero(&self.re) && num_traits::Zero::is_zero(&self.im)
    }
    fn add(&self, rhs: &Self) -> Self {
        GaussQ {
            re: &self.re + &rhs.re,
            im: &self.im + &rhs.im,
        }
    }
    fn sub(&self, rhs: &Self) -> Self {
        GaussQ {
            re: &self.re - &rhs.re,
            im: &self.im - &rhs.im,
        }
    }
    fn mul(&self, rhs: &Self) -> Self {
        GaussQ {
            re: &self.re * &rhs.re - &self.im * &rhs.im,
            im: &self.re * &rhs.im + &self.im * &rhs.re,
        }
    }
    fn neg(&self) -> Self {
        GaussQ {
            re: -&self.re,
            im: -&self.im,
        }
    }
}

impl crate::poly::traits::EuclideanDomain for GaussQ {
    fn div_rem(&self, other: &Self) -> (Self, Self) {
        (crate::poly::traits::Field::div(self, other), GaussQ::zero())
    }
}

impl crate::poly::traits::Field for GaussQ {
    fn div(&self, other: &Self) -> Self {
        self.mul(&other.inv())
    }
    /// `(a + b·i)⁻¹ = (a − b·i)/(a² + b²)`; callers only invert non-zero
    /// elements (leading coefficients, a gcd's constant), and a zero one
    /// gives zero rather than a panic.
    fn inv(&self) -> Self {
        let norm = &self.re * &self.re + &self.im * &self.im;
        if num_traits::Zero::is_zero(&norm) {
            return GaussQ::zero();
        }
        GaussQ {
            re: &self.re / &norm,
            im: -&self.im / &norm,
        }
    }
}

type GaussPoly = GenPoly<GaussQ>;

/// Largest degree in the root variable [`gauss_fraction`] builds.
const MAX_GAUSS_DEGREE: usize = 64;

/// `e` as a quotient `(P, Q)` of polynomials in `t` over `ℚ(i)`, when it
/// is built from rational numbers, `i` and `t` by sums, products and
/// integer powers (`None` for anything else: another symbol, `√2`, `π`).
/// Nothing is cancelled.
fn gauss_fraction(arena: &Arena, e: ExprId, t: ExprId) -> Option<(GaussPoly, GaussPoly)> {
    let mut cache: rustc_hash::FxHashMap<ExprId, (GaussPoly, GaussPoly)> =
        rustc_hash::FxHashMap::default();
    let one = GaussPoly::one();
    for id in crate::base::walk::post_order_ids(arena, e) {
        let value = match arena.node(id) {
            ExprNode::Num(_) => {
                let q = arena.as_num(id)?.clone();
                (GaussPoly::constant(GaussQ::real(q)), one.clone())
            }
            ExprNode::ImaginaryUnit => {
                let i = GaussQ {
                    re: Q::from_integer(BigInt::from(0)),
                    im: Q::from_integer(BigInt::from(1)),
                };
                (GaussPoly::constant(i), one.clone())
            }
            ExprNode::Symbol(_) if id == t => (GaussPoly::x(), one.clone()),
            ExprNode::Add(children) => {
                let mut acc = (GaussPoly::zero(), one.clone());
                for c in children.iter() {
                    let (p, q) = cache.get(c)?;
                    acc = if *q == acc.1 {
                        (acc.0.add(p), acc.1)
                    } else {
                        (acc.0.mul(q).add(&p.mul(&acc.1)), acc.1.mul(q))
                    };
                }
                acc
            }
            ExprNode::Mul(children) => {
                let mut acc = (one.clone(), one.clone());
                for c in children.iter() {
                    let (p, q) = cache.get(c)?;
                    acc = (acc.0.mul(p), acc.1.mul(q));
                }
                acc
            }
            ExprNode::Neg(inner) => {
                let (p, q) = cache.get(inner)?;
                (p.neg(), q.clone())
            }
            ExprNode::Pow(base, exp) => {
                let n = arena.as_num(*exp)?;
                if !n.is_integer() {
                    return None;
                }
                let n = i64::try_from(n.to_integer()).ok()?;
                let k = usize::try_from(n.unsigned_abs()).ok()?;
                let (p, q) = cache.get(base)?;
                if k.saturating_mul(p.degree().unwrap_or(0).max(q.degree().unwrap_or(0)))
                    > MAX_GAUSS_DEGREE
                {
                    return None;
                }
                if n >= 0 {
                    (p.pow(k), q.pow(k))
                } else {
                    if p.is_zero() {
                        return None;
                    }
                    (q.pow(k), p.pow(k))
                }
            }
            _ => return None,
        };
        if value.0.degree().unwrap_or(0) > MAX_GAUSS_DEGREE
            || value.1.degree().unwrap_or(0) > MAX_GAUSS_DEGREE
            || value.1.is_zero()
        {
            return None;
        }
        cache.insert(id, value);
    }
    cache.remove(&e)
}

/// The exact value of `RootSum(poly, body, sumvar)` = `Σ_{poly(r) = 0}
/// body(r)` when `poly` and `body` have no free symbol but `sumvar`, their
/// coefficients are Gaussian rationals, `poly` is square-free and `body` is
/// a rational function of the root without a pole at one: `body ≡ S mod
/// poly` with `deg S < deg poly` (the inverse of the body's denominator
/// modulo `poly` by the extended Euclidean algorithm), and `Σ S(r) =
/// Σₖ sₖ·pₖ` with the power sums `pₖ = Σ rᵏ` of the roots from
/// Newton's identities.  `None` otherwise.
///
/// The integrator's self-check meets such sums at every sample point: the
/// derivative of `RootSum(e, c(t)·ln(x − t))` is `RootSum(e, c(t)/(x − t))`,
/// numeric once `x` and the parameters are bound.  Solving `e` in radicals
/// (as `rootsum_doit` does for `eval`) and evaluating the radical
/// expression to 30 digits took most of the check's time (0.31: 8 ms per
/// check of a `RootSum` answer, most of the cost of the Rubi entries
/// `∫ P(x)/(a + b·x² + c·x⁴)² dx`); the trace is exact rational arithmetic.
pub(crate) fn numeric_rootsum_rational_body(
    arena: &mut Arena,
    poly_id: ExprId,
    body_id: ExprId,
    sumvar_id: ExprId,
) -> Option<ExprId> {
    use crate::poly::traits::Field;
    let (e_num, e_den) = gauss_fraction(arena, poly_id, sumvar_id)?;
    if e_den.degree()? != 0 {
        return None;
    }
    let n = e_num.degree()?;
    if n == 0 {
        return None;
    }
    let e = e_num.make_monic();
    if GaussPoly::gcd(&e, &e.derivative()).degree()? != 0 {
        return None;
    }
    let (p, q) = gauss_fraction(arena, body_id, sumvar_id)?;
    let q_mod = q.rem(&e);
    if q_mod.is_zero() {
        return None;
    }
    let eg = GaussPoly::extended_gcd(&q_mod, &e);
    if eg.gcd.degree()? != 0 {
        return None;
    }
    let inv = eg.gcd.leading_coeff()?.inv();
    let s = p.mul(&eg.x.scale(&inv)).rem(&e);
    // Newton: pₖ = −k·a_{n−k} − Σ_{i<k} a_{n−i}·p_{k−i} for the monic
    // e = tⁿ + a_{n−1}·tⁿ⁻¹ + … + a₀.
    let a = e.coeffs();
    let mut power_sums: Vec<GaussQ> = Vec::with_capacity(n);
    power_sums.push(GaussQ::from_usize(n));
    for k in 1..n {
        let mut pk = a[n - k].mul(&GaussQ::from_usize(k)).neg();
        for i in 1..k {
            pk = pk.sub(&a[n - i].mul(&power_sums[k - i]));
        }
        power_sums.push(pk);
    }
    let mut value = GaussQ::zero();
    for (k, sk) in s.coeffs().iter().enumerate() {
        value = value.add(&sk.mul(power_sums.get(k)?));
    }
    let re = arena.num_ratio(value.re);
    if num_traits::Zero::is_zero(&value.im) {
        return Some(re);
    }
    let im = arena.num_ratio(value.im);
    let i = arena.i_unit();
    let im_i = arena.mul(&[im, i]);
    Some(arena.add(&[re, im_i]))
}

// ═══════════════════════════════════════════════════════════════════════════
// rootsum_doit — expand RootSum when the polynomial is solvable
// ═══════════════════════════════════════════════════════════════════════════

/// Try to expand a `RootSum(poly, body, sumvar)` into an explicit sum
/// by finding the roots of the polynomial via [`solve`](crate::transforms::solve::solve).
///
/// For degree ≤ 4 polynomials, `solve` produces exact radical roots
/// (quadratic formula, Cardano, Ferrari).  Each root is substituted
/// into the body, and the results are summed.
///
/// Returns `Some(expanded_sum)` if all roots are found, `None` if
/// `solve` can't find them (degree ≥ 5 non-solvable).
///
/// # Example
///
/// ```text
/// RootSum(9t²+3t+1, t -> t·ln(x-3t))
///   → (-1/6+i√3/6)·ln(x-3(-1/6+i√3/6)) + (-1/6-i√3/6)·ln(x-3(-1/6-i√3/6))
/// ```
///
/// The result can then be simplified via `eval` / `as_real_imag` / etc.
///
/// A polynomial with parameters (the parametric rational integrator's
/// `RootSum(t³ + a, …)`) is not expanded: its radical roots (`∛(−a)`,
/// Cardano's formula in the parameters) multiply the size by the degree,
/// and by parts then integrated such a sum term by term — a 300-term,
/// 6 s antiderivative of `x¹⁹/(a + b·x⁵)` while the `RootSum` route was
/// being added (0.31).  SymPy keeps it as well; once the parameters have
/// values the sum is expanded like any other.
pub(crate) fn rootsum_doit(
    arena: &mut Arena,
    poly_id: ExprId,
    body_id: ExprId,
    sumvar_id: ExprId,
) -> Option<ExprId> {
    let ExprNode::Symbol(sumvar_sym) = *arena.node(sumvar_id) else {
        return None;
    };
    let parametric = crate::base::walk::free_symbols(arena, poly_id)
        .iter()
        .any(|&s| !matches!(arena.node(s), ExprNode::Symbol(t) if *t == sumvar_sym));
    if parametric {
        tracing::debug!("rootsum_doit: the polynomial has parameters, keeping RootSum");
        return None;
    }
    // Solve the polynomial for roots.
    let roots = crate::transforms::solve::solve(arena, poly_id, sumvar_id);

    if roots.is_empty() {
        tracing::debug!("rootsum_doit: solve returned no roots, cannot expand");
        return None;
    }
    // Expanding over `RootOf` placeholders gains nothing and multiplies the
    // size by the degree: since 0.22 `solve` returns one `RootOf` per root of
    // an irreducible factor of degree ≥ 5, and `∫ atan(√x − x⁹) dx` became a
    // 486 KB sum of twelve copies of the body.  Keep the compact `RootSum`.
    let has_rootof = roots.iter().any(|r| {
        crate::base::walk::post_order_ids(arena, r.value)
            .into_iter()
            .any(|id| matches!(arena.node(id), ExprNode::RootOf(..)))
    });
    if has_rootof {
        tracing::debug!("rootsum_doit: roots are RootOf placeholders, keeping RootSum");
        return None;
    }

    // Check that we got the expected number of roots (= degree of poly).
    let poly_obj = crate::poly::polybridge::expr_to_poly(arena, poly_id, sumvar_id);
    if let Some(ref p) = poly_obj
        && let Some(deg) = p.degree()
        && roots.len() != deg
    {
        tracing::debug!(
            expected = deg,
            found = roots.len(),
            "rootsum_doit: solve found fewer roots than polynomial degree, cannot expand fully"
        );
        return None;
    }

    tracing::debug!(
        n_roots = roots.len(),
        "rootsum_doit: expanding RootSum by substituting each root"
    );

    // For each root α_k: substitute sumvar = α_k into body, eval.
    let mut terms: Vec<ExprId> = Vec::with_capacity(roots.len());
    for root in &roots {
        let substituted = crate::transforms::subs::subs(arena, body_id, sumvar_id, root.value);
        let evaluated = crate::transforms::eval::eval(arena, substituted);
        terms.push(evaluated);
    }

    // Sum all terms.
    let result = if terms.len() == 1 {
        terms[0]
    } else {
        arena.add(&terms)
    };

    tracing::debug!("rootsum_doit: expansion complete");
    Some(result)
}

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

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

    // ── Vieta tests ────────────────────────────────────────────────

    #[test]
    fn vieta_elementary_symmetric_quadratic() {
        // q(t) = t² + 3t + 2 = (t+1)(t+2)
        // Roots: -1, -2
        // e_1 = -(-1-2) = 3  (but Vieta: e_1 = -a_{n-1} = -3)
        // Wait: for t² + 3t + 2, a_1=3, a_0=2
        // e_1 = -a_1 = -3  (sum of roots = -1 + -2 = -3) ✓
        // e_2 = a_0 = 2   (product of roots = (-1)(-2) = 2) ✓
        let q = Poly::from_coeffs(vec![r(2, 1), r(3, 1), r(1, 1)]);
        let e = vieta_elementary_symmetric(&q).unwrap();
        assert_eq!(e.len(), 2);
        assert_eq!(e[0], r(-3, 1), "e_1 = sum of roots = -3");
        assert_eq!(e[1], r(2, 1), "e_2 = product of roots = 2");
    }

    #[test]
    fn vieta_power_sums_quadratic() {
        // Roots: -1, -2
        // p_1 = (-1) + (-2) = -3
        // p_2 = (-1)² + (-2)² = 1 + 4 = 5
        // p_3 = (-1)³ + (-2)³ = -1 + (-8) = -9
        let e = vec![r(-3, 1), r(2, 1)];
        let p = vieta_power_sums(&e, 3);
        assert_eq!(p.len(), 3);
        assert_eq!(p[0], r(-3, 1), "p_1 = -3");
        assert_eq!(p[1], r(5, 1), "p_2 = 5");
        assert_eq!(p[2], r(-9, 1), "p_3 = -9");
    }

    #[test]
    fn vieta_power_sums_cubic() {
        // q(t) = t³ - 6t² + 11t - 6 = (t-1)(t-2)(t-3)
        // Roots: 1, 2, 3
        // e_1 = 6, e_2 = 11, e_3 = 6
        let q = Poly::from_coeffs(vec![r(-6, 1), r(11, 1), r(-6, 1), r(1, 1)]);
        let e = vieta_elementary_symmetric(&q).unwrap();
        assert_eq!(e[0], r(6, 1), "e_1 = 1+2+3 = 6");
        assert_eq!(e[1], r(11, 1), "e_2 = 1·2+1·3+2·3 = 11");
        assert_eq!(e[2], r(6, 1), "e_3 = 1·2·3 = 6");

        // p_1 = 6, p_2 = 1+4+9 = 14, p_3 = 1+8+27 = 36
        let p = vieta_power_sums(&e, 3);
        assert_eq!(p[0], r(6, 1), "p_1 = 6");
        assert_eq!(p[1], r(14, 1), "p_2 = 14");
        assert_eq!(p[2], r(36, 1), "p_3 = 36");
    }

    #[test]
    fn vieta_rootsum_poly_body_sum_of_roots() {
        // RootSum(t²+3t+2, t -> t) = sum of roots = -3
        let mut arena = Arena::new();
        let t = arena.symbol("t");
        let two = arena.int(2);
        let three = arena.int(3);
        let t_sq = arena.pow(t, two);
        let poly_expr = {
            let three_t = arena.mul(&[three, t]);
            arena.add(&[t_sq, three_t, two])
        };

        let result = vieta_rootsum_poly_body(&arena, poly_expr, t, t);
        assert_eq!(
            result,
            Some(r(-3, 1)),
            "sum of roots of t²+3t+2 should be -3"
        );
    }

    #[test]
    fn vieta_rootsum_poly_body_sum_of_squares() {
        // RootSum(t²+3t+2, t -> t²) = sum of squares of roots = 5
        let mut arena = Arena::new();
        let t = arena.symbol("t");
        let two = arena.int(2);
        let three = arena.int(3);
        let t_sq = arena.pow(t, two);
        let poly_expr = {
            let three_t = arena.mul(&[three, t]);
            arena.add(&[t_sq, three_t, two])
        };
        let body = arena.pow(t, two); // t²

        let result = vieta_rootsum_poly_body(&arena, poly_expr, body, t);
        assert_eq!(
            result,
            Some(r(5, 1)),
            "sum of squares of roots of t²+3t+2 should be 5"
        );
    }

    /// The exact trace of a numeric `RootSum` with a rational body, over
    /// `ℚ` and over `ℚ(i)`; `None` for a body with a pole at a root.
    #[test]
    fn numeric_rootsum_with_a_rational_body() {
        let mut arena = Arena::new();
        let t = arena.symbol("t");
        let two = arena.int(2);
        let three = arena.int(3);
        let t2 = arena.pow(t, two);
        let three_t = arena.mul(&[three, t]);
        // t² + 3t + 2 = (t + 1)(t + 2): Σ 1/(r + 3) = 1/2 + 1 = 3/2.
        let q = arena.add(&[t2, three_t, two]);
        let t_plus_3 = arena.add(&[t, three]);
        let neg_one = arena.int(-1);
        let body = arena.pow(t_plus_3, neg_one);
        let v = numeric_rootsum_rational_body(&mut arena, q, body, t).unwrap();
        assert_eq!(display(&arena, v), "3/2");
        // A pole at a root (1/(t + 1)): no value.
        let one = arena.one;
        let t_plus_1 = arena.add(&[t, one]);
        let pole = arena.pow(t_plus_1, neg_one);
        assert!(numeric_rootsum_rational_body(&mut arena, q, pole, t).is_none());
        // t² + 1 (roots ±i): Σ 1/(r − i·2) = 1/(−i) + 1/(−3i) = 4i/3.
        let t2_plus_1 = arena.add(&[t2, one]);
        let i = arena.i_unit();
        let two_i = arena.mul(&[two, i]);
        let neg_two_i = arena.neg(two_i);
        let shifted = arena.add(&[t, neg_two_i]);
        let body = arena.pow(shifted, neg_one);
        let v = numeric_rootsum_rational_body(&mut arena, t2_plus_1, body, t).unwrap();
        assert_eq!(display(&arena, v), "4/3*I");
        // A parameter: no value.
        let a = arena.symbol("a");
        let t_plus_a = arena.add(&[t, a]);
        let body = arena.pow(t_plus_a, neg_one);
        assert!(numeric_rootsum_rational_body(&mut arena, q, body, t).is_none());
    }

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

    #[test]
    fn vieta_rootsum_constant_body() {
        // RootSum(t²+3t+2, t -> 7) = 2 * 7 = 14 (n roots, each contributing 7)
        let mut arena = Arena::new();
        let t = arena.symbol("t");
        let two = arena.int(2);
        let three = arena.int(3);
        let seven = arena.int(7);
        let t_sq = arena.pow(t, two);
        let poly_expr = {
            let three_t = arena.mul(&[three, t]);
            arena.add(&[t_sq, three_t, two])
        };

        let result = vieta_rootsum_poly_body(&arena, poly_expr, seven, t);
        assert_eq!(
            result,
            Some(r(14, 1)),
            "RootSum with constant body 7 over degree-2 poly = 14"
        );
    }

    fn r(n: i64, d: i64) -> Q {
        Ratio::new(BigInt::from(n), BigInt::from(d))
    }

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

    // ── poly_to_genpoly_rf tests ───────────────────────────────────

    #[test]
    fn poly_to_genpoly_rf_constant() {
        let p = Poly::from_int(5);
        let gp = poly_to_genpoly_rf(&p);
        assert_eq!(gp.degree(), Some(0));
        let c = gp.coeff(0);
        assert_eq!(c.to_rational(), Some(r(5, 1)));
    }

    #[test]
    fn poly_to_genpoly_rf_linear() {
        // p(x) = 3x + 2
        let p = Poly::from_coeffs(vec![r(2, 1), r(3, 1)]);
        let gp = poly_to_genpoly_rf(&p);
        assert_eq!(gp.degree(), Some(1));
        assert_eq!(gp.coeff(0).to_rational(), Some(r(2, 1)));
        assert_eq!(gp.coeff(1).to_rational(), Some(r(3, 1)));
    }

    #[test]
    fn poly_to_genpoly_rf_times_t_linear() {
        // p(x) = 3x + 2 → 3t·x + 2t
        let p = Poly::from_coeffs(vec![r(2, 1), r(3, 1)]);
        let gp = poly_to_genpoly_rf_times_t(&p);
        assert_eq!(gp.degree(), Some(1));
        // coeff(0) = 2t → numer should be 2t polynomial
        let c0 = gp.coeff(0);
        assert_eq!(c0.numer().degree(), Some(1));
        assert_eq!(c0.numer().coeff(0), r(0, 1));
        assert_eq!(c0.numer().coeff(1), r(2, 1));
    }

    // ── log_to_atan_deg1 tests ─────────────────────────────────────

    #[test]
    fn log_to_atan_constant_b() {
        // A(x) = x + 1/2, B(x) = -√3/2
        // → 2·atan((x + 1/2) / (-√3/2)) = 2·atan(-(2x+1)/√3)
        let mut arena = Arena::new();
        let x = sym(&mut arena, "x");
        let a1 = arena.one;
        let a0 = arena.rational(1, 2);
        let b1 = arena.zero;
        let three = arena.int(3);
        let half = arena.rational(1, 2);
        let sqrt3 = arena.pow(three, half);
        let half2 = arena.rational(1, 2);
        let half_sqrt3 = arena.mul(&[half2, sqrt3]);
        let b0 = arena.neg(half_sqrt3);

        let result = log_to_atan_deg1(&mut arena, x, a1, a0, b1, b0);
        assert!(result.is_some(), "should produce atan for constant B");

        // Verify numerically at x = 1
        let one = arena.int(1);
        let at_1 = crate::transforms::subs::subs(&mut arena, result.unwrap(), x, one);
        let at_1 = crate::transforms::eval::eval(&mut arena, at_1);
        let val = crate::transforms::evalf::eval_const_f64(&mut arena, at_1);
        assert!(val.is_some(), "should evaluate to f64");

        // Expected: 2·atan((1 + 0.5) / (-√3/2)) = 2·atan(1.5 / (-0.866)) = 2·atan(-1.7321)
        let expected = 2.0 * (1.5_f64 / (-(3.0_f64.sqrt()) / 2.0)).atan();
        let v = val.unwrap();
        assert!(
            (v - expected).abs() < 1e-8,
            "log_to_atan at x=1: got {v}, expected {expected}"
        );
    }

    #[test]
    fn log_to_atan_degenerate_b_zero() {
        // B(x) = 0 → should return None
        let mut arena = Arena::new();
        let x = sym(&mut arena, "x");
        let one = arena.one;
        let zero = arena.zero;
        let result = log_to_atan_deg1(&mut arena, x, one, zero, zero, zero);
        assert!(result.is_none(), "B = 0 should be degenerate");
    }

    // ── log_to_real tests ──────────────────────────────────────────

    #[test]
    fn log_to_real_x3_minus_1_quadratic_factor() {
        // For 1/(x³-1), the quadratic factor of R(t) is q(t) = 9t²+3t+1.
        // The PRS degree-1 member is h(t,x) = x - 3t.
        //
        // log_to_real should produce:
        //   u·ln(x²+x+1) + v·(2·atan(...))
        // where u = -1/6 and v = √3/6.
        let mut arena = Arena::new();
        let x = sym(&mut arena, "x");

        // q(t) = 9t² + 3t + 1
        let q = Poly::from_coeffs(vec![r(1, 1), r(3, 1), r(9, 1)]);

        // h(t,x) = x - 3t as GenPoly<RationalFn>
        let neg_3t = RationalFn::from_poly(Poly::from_coeffs(vec![r(0, 1), r(-3, 1)]));
        let h_prs: GenPoly<RationalFn> =
            GenPoly::from_coeffs(vec![neg_3t, RationalFn::from_rational(r(1, 1))]);

        let terms = log_to_real(&mut arena, x, &q, &h_prs);
        assert!(
            terms.is_some(),
            "log_to_real should succeed for x³-1 quadratic factor"
        );

        let terms = terms.unwrap();
        assert!(
            terms.len() >= 2,
            "should produce at least 2 terms (ln + atan)"
        );

        // Sum all terms and evaluate the definite integral [2, 3].
        let sum = arena.add(&terms);
        let val_3 = arena.int(3);
        let f3 = crate::transforms::subs::subs(&mut arena, sum, x, val_3);
        let f3 = crate::transforms::eval::eval(&mut arena, f3);
        let val_2 = arena.int(2);
        let f2 = crate::transforms::subs::subs(&mut arena, sum, x, val_2);
        let f2 = crate::transforms::eval::eval(&mut arena, f2);

        let f3_f64 = crate::transforms::evalf::eval_const_f64(&mut arena, f3);
        let f2_f64 = crate::transforms::evalf::eval_const_f64(&mut arena, f2);

        if let (Some(f3v), Some(f2v)) = (f3_f64, f2_f64) {
            let integral = f3v - f2v;
            // This is ONLY the quadratic-factor contribution.
            // Full answer = 1/3·ln|x-1| + (this).
            // Ground truth for the quadratic part on [2,3]: ≈ -0.1557
            let expected = {
                let gt_full = 0.07539_f64;
                let gt_linear = (1.0 / 3.0) * (2.0_f64.ln()); // 1/3 · ln(2)
                gt_full - gt_linear
            };
            assert!(
                (integral - expected).abs() < 0.01,
                "quadratic factor ∫₂³: got {integral}, expected ≈ {expected}"
            );
        } else {
            panic!("could not evaluate log_to_real result to f64");
        }
    }

    /// `f(x0)` as an `f64`.
    fn value_at(arena: &mut Arena, f: ExprId, x: ExprId, x0: Q) -> f64 {
        let x0 = arena.num_ratio(x0);
        let v = crate::transforms::subs::subs(arena, f, x, x0);
        let v = crate::transforms::eval::eval(arena, v);
        crate::transforms::evalf::eval_const_f64(arena, v).expect("numeric value")
    }

    fn poly_q(coeffs: &[(i64, i64)]) -> Poly {
        Poly::from_coeffs(coeffs.iter().map(|&(n, d)| r(n, d)).collect())
    }

    fn ratfn_poly(coeffs: &[(i64, i64)]) -> RationalFn {
        RationalFn::from_poly(poly_q(coeffs))
    }

    #[test]
    fn quadratic_factor_with_a_degree_4_log_argument() {
        // Σ_{64α²+1=0} α·ln(x⁴ + 8α) = atan(x⁴)/4 (the log part of
        // ∫ x³/(x⁸ + 1) dx; SymPy 1.14 `integrate(x**3/(x**8+1), x)`).
        let mut arena = Arena::new();
        let x = sym(&mut arena, "x");
        let q = poly_q(&[(1, 1), (0, 1), (64, 1)]);
        let zero = ratfn_poly(&[]);
        let s: GenPoly<RationalFn> = GenPoly::from_coeffs(vec![
            ratfn_poly(&[(0, 1), (8, 1)]),
            zero.clone(),
            zero.clone(),
            zero,
            ratfn_poly(&[(1, 1)]),
        ]);
        let terms = quadratic_log_to_real(&mut arena, x, &q, &s).expect("real form");
        let sum = arena.add(&terms);
        assert_eq!(arena.display(sum).to_string(), "1/4*atan(x^4)");
    }

    #[test]
    fn quadratic_factor_with_real_roots() {
        // Σ_{8α²−1=0} α·ln(x − 4α) = (√2/4)·(ln|x − √2| − ln|x + √2|), the
        // log part of ∫ 1/(x² − 2) dx (SymPy 1.14 `integrate(1/(x**2-2), x)`).
        let mut arena = Arena::new();
        let x = sym(&mut arena, "x");
        let q = poly_q(&[(-1, 1), (0, 1), (8, 1)]);
        let s: GenPoly<RationalFn> =
            GenPoly::from_coeffs(vec![ratfn_poly(&[(0, 1), (-4, 1)]), ratfn_poly(&[(1, 1)])]);
        let terms = quadratic_log_to_real(&mut arena, x, &q, &s).expect("real form");
        let sum = arena.add(&terms);
        for x0 in [r(3, 1), r(1, 3), r(-7, 2)] {
            let got = value_at(&mut arena, sum, x, x0.clone());
            let xf = num_traits::ToPrimitive::to_f64(&x0).expect("f64");
            let s2 = 2f64.sqrt();
            let want = s2 / 4.0 * ((xf - s2).abs().ln() - (xf + s2).abs().ln());
            assert!((got - want).abs() < 1e-12, "x = {xf}: {got} vs {want}");
        }
    }

    #[test]
    fn log_to_atan_surd_satisfies_its_defining_derivative() {
        // F' = 2(A'·rB − A·rB')/(A² + r²B²) for F = LogToAtan(A, rB), here with
        // A = x³ + 2x, B = x² − 3, r = √3: B ∤ A, so the Bézout step and the
        // recursion both run.
        let mut arena = Arena::new();
        let x = sym(&mut arena, "x");
        let a = poly_q(&[(0, 1), (2, 1), (0, 1), (1, 1)]);
        let b = poly_q(&[(-3, 1), (0, 1), (1, 1)]);
        let r2 = r(3, 1);
        let f = log_to_atan_surd(&mut arena, x, &a, &b, &r2).expect("atan sum");
        let df = crate::transforms::diff::diff(&mut arena, f, x);
        for x0 in [r(1, 2), r(5, 3), r(-9, 4)] {
            let got = value_at(&mut arena, df, x, x0.clone());
            let xf = num_traits::ToPrimitive::to_f64(&x0).expect("f64");
            let rt = 3f64.sqrt();
            let (av, dav) = (xf.powi(3) + 2.0 * xf, 3.0 * xf * xf + 2.0);
            let (bv, dbv) = (rt * (xf * xf - 3.0), rt * 2.0 * xf);
            let want = 2.0 * (dav * bv - av * dbv) / (av * av + bv * bv);
            assert!((got - want).abs() < 1e-10, "x = {xf}: {got} vs {want}");
        }
    }

    #[test]
    fn log_to_real_keeps_the_real_roots_of_a_cubic() {
        // ∫ 1/(x³ − 2) dx: q = 108t³ − 1 has the real root ∛2/6, whose term
        // (∛2/6)·ln|x − ∛2| was dropped before 0.25.  Checked through the
        // whole integrator: F' = f at three points.
        let mut arena = Arena::new();
        let x = sym(&mut arena, "x");
        let three = arena.int(3);
        let x3 = arena.pow(x, three);
        let m2 = arena.int(-2);
        let den = arena.add(&[x3, m2]);
        let m1 = arena.int(-1);
        let f = arena.pow(den, m1);
        let big_f =
            crate::calculus::risch::try_risch_rational(&mut arena, f, x).expect("closed form");
        assert!(!crate::base::walk::has_unevaluated(&arena, big_f));
        let df = crate::transforms::diff::diff(&mut arena, big_f, x);
        for x0 in [r(3, 1), r(1, 2), r(-5, 2)] {
            let got = value_at(&mut arena, df, x, x0.clone());
            let xf = num_traits::ToPrimitive::to_f64(&x0).expect("f64");
            let want = 1.0 / (xf.powi(3) - 2.0);
            assert!((got - want).abs() < 1e-10, "x = {xf}: {got} vs {want}");
        }
    }

    #[test]
    fn log_to_real_degree_less_than_2_returns_none() {
        let mut a = Arena::new();
        let x = sym(&mut a, "x");

        // Linear q(t) = 3t - 1 → should return None
        let q = Poly::from_coeffs(vec![r(-1, 1), r(3, 1)]);
        let h_prs: GenPoly<RationalFn> = GenPoly::from_coeffs(vec![
            RationalFn::from_rational(r(0, 1)),
            RationalFn::from_rational(r(1, 1)),
        ]);

        let result = log_to_real(&mut a, x, &q, &h_prs);
        assert!(result.is_none(), "degree-1 q should return None");
    }
}