symplex 0.9.0

//! Arbitrary-precision numerical evaluation with complex number support.
//!
//! This module implements [`evalf`], which evaluates a symbolic expression
//! to a decimal string with a specified number of significant digits,
//! using the [`astro_float`] crate for arbitrary-precision arithmetic.
//!
//! # Design
//!
//! Evaluation is performed bottom-up using an explicit post-order
//! traversal (no recursion).  Each sub-expression is evaluated to a
//! [`Complex`] (a pair of [`BigFloat`] values for the real and imaginary
//! parts) and cached, so shared sub-expressions (common in a hash-consed
//! DAG) are only evaluated once.
//!
//! The working precision is set higher than the requested output
//! precision to absorb rounding errors from intermediate computations.
//! If a sub-expression cannot be evaluated (e.g., it contains free
//! symbols), the function returns an error.
//!
//! # Definite integrals
//!
//! An unevaluated `DefiniteIntegral(body, var, lo, hi)` node is evaluated
//! by compiling `body` to an `f64` function and running the adaptive
//! Gauss–Kronrod quadrature from [`crate::calculus::definite`].  That
//! delivers double precision at best, so requests for more than
//! [`QUADRATURE_MAX_DIGITS`] digits are refused with
//! [`SymplexError::NotImplemented`] rather than padded with digits that
//! carry no information.  `Ex::eval_f64` (16 digits) is served; a 30-digit
//! `eval_decimal` is not.

use num_bigint::BigInt;
use num_rational::Ratio;
use num_traits::Zero;
use rustc_hash::FxHashMap;

use astro_float::{BigFloat, Consts, Radix, RoundingMode, Sign};

use crate::base::arena::Arena;
use crate::base::errors::SymplexError;
use crate::base::node::{ExprId, ExprNode};
use crate::base::walk;
use tracing::debug;

// ═══════════════════════════════════════════════════════════════════════════
// Complex type
// ═══════════════════════════════════════════════════════════════════════════

/// A complex number represented as (real_part, imaginary_part).
type Complex = (BigFloat, BigFloat);

/// Largest number of decimal digits for which an expression containing a
/// `DefiniteIntegral` is evaluated (by `f64` quadrature) instead of
/// refused.  `Ex::eval_f64` requests exactly this many.
pub(crate) const QUADRATURE_MAX_DIGITS: u32 = 16;

// ═══════════════════════════════════════════════════════════════════════════
// Public entry point
// ═══════════════════════════════════════════════════════════════════════════

/// Evaluate `expr` numerically to `digits` decimal digits of precision.
///
/// Returns the decimal string representation of the result, or an error
/// if the expression contains free symbols, infinities, or NaN.
pub(crate) fn evalf(arena: &Arena, expr: ExprId, digits: u32) -> Result<String, SymplexError> {
    // Convert decimal digits to binary precision with guard bits.
    // log2(10) ≈ 3.3219, so we use digits * 3.4 + 64 extra bits.
    let binary_prec = (digits as usize) * 34 / 10 + 64;
    // Ensure a reasonable minimum.
    let prec = binary_prec.max(128);

    debug!(
        digits = digits,
        binary_precision = prec,
        "evalf: starting numerical evaluation"
    );

    // Enforce the configured maximum precision.
    let max_prec = arena.config.max_evalf_precision as usize;
    if prec > max_prec {
        return Err(SymplexError::PrecisionExhausted {
            requested: digits,
            achieved: (max_prec * 10 / 34).saturating_sub(6) as u32,
        });
    }

    // Free symbols make the whole expression unevaluable; report the
    // offending symbol by name up front instead of surfacing an internal
    // cache miss from some parent node later.  Bound variables (a `Sum`
    // index, a `RootOf` polynomial variable, …) are not free and are
    // handled by the sub-tree evaluators below.
    let free = walk::free_symbols(arena, expr);
    if !free.is_empty() {
        let mut names: Vec<&str> = free
            .iter()
            .filter_map(|&id| match arena.node(id) {
                ExprNode::Symbol(sid) => Some(arena.symbol_name(*sid)),
                _ => None,
            })
            .collect();
        names.sort_unstable();
        if let Some(name) = names.first() {
            debug!(symbol = name, "evalf: expression has free symbols");
            return Err(SymplexError::FreeSymbol {
                name: (*name).to_owned(),
            });
        }
    }

    // Definite integrals are evaluated by `f64` quadrature; refuse rather
    // than print digits that are not there.
    if digits > QUADRATURE_MAX_DIGITS && contains_definite_integral(arena, expr) {
        return Err(SymplexError::NotImplemented(format!(
            "arbitrary-precision evaluation of a definite integral is not implemented \
             ({digits} digits requested, at most {QUADRATURE_MAX_DIGITS} available from \
             f64 Gauss–Kronrod quadrature); use eval_f64 or integrate_definite"
        )));
    }

    let rm = RoundingMode::ToEven;
    let mut cc = Consts::new().map_err(|e| {
        SymplexError::NotImplemented(format!("astro-float constants init failed: {e:?}"))
    })?;

    // Post-order traversal — evaluate children before parents.
    let post_order = walk::post_order_ids(arena, expr);
    let mut cache: FxHashMap<ExprId, Complex> = FxHashMap::default();

    for &id in &post_order {
        match eval_node(arena, id, &cache, prec, rm, &mut cc) {
            Ok(value) => {
                cache.insert(id, value);
            }
            Err(e) => {
                if id == expr {
                    return Err(e);
                }
                // Non-root node failed — skip it.  Parent nodes
                // (Sum, Product, Piecewise) will handle their own
                // sub-tree evaluation with substitution.
            }
        }
    }

    let result = cache.get(&expr).ok_or_else(|| {
        SymplexError::NotImplemented("evalf: expression not found in cache".into())
    })?;

    if result.0.is_nan() || result.1.is_nan() {
        return Err(SymplexError::PrecisionExhausted {
            requested: digits,
            achieved: 0,
        });
    }

    format_complex(result, digits, prec, rm, &mut cc)
}

/// Evaluate a constant (variable-free) arena expression to `f64`.
///
/// Returns `None` if the expression contains free symbols, produces
/// a complex result, or evaluation fails for any reason.
///
/// The expression is first simplified via [`eval`](crate::transforms::eval::eval)
/// to reduce exact values (e.g., `sin(π) → 0`, `Γ(5) → 24`) before
/// numerical evaluation.
pub(crate) fn eval_const_f64(arena: &mut Arena, expr: ExprId) -> Option<f64> {
    // Simplify exact values first.
    let evaled = crate::transforms::eval::eval(arena, expr);

    // Fast path: if the result is already a rational number, convert directly
    // without invoking the expensive arbitrary-precision machinery.
    if let Some(r) = arena.as_num(evaled) {
        let n: f64 = r.numer().to_string().parse().ok()?;
        let d: f64 = r.denom().to_string().parse().ok()?;
        if d == 0.0 {
            tracing::debug!("eval_const_f64: rational with zero denominator");
            return None;
        }
        let result = n / d;
        tracing::trace!(result, "eval_const_f64: rational fast path");
        return Some(result);
    }

    // Fall back to arbitrary-precision evaluation to 16 decimal digits,
    // then parse the resulting string.  The `evalf` call takes `&Arena`
    // (immutable), which Rust allows via automatic reborrowing of our
    // `&mut Arena`.
    match evalf(arena, evaled, 16) {
        Ok(s) => {
            let result = s.parse::<f64>().ok();
            tracing::trace!(?result, decimal_str = %s, "eval_const_f64: evalf path");
            result
        }
        Err(e) => {
            tracing::debug!(?e, "eval_const_f64: evalf failed");
            None
        }
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Node evaluation
// ═══════════════════════════════════════════════════════════════════════════

/// Evaluate a single node given its children's values in the cache.
fn eval_node(
    arena: &Arena,
    id: ExprId,
    cache: &FxHashMap<ExprId, Complex>,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<Complex, SymplexError> {
    let node = arena.node(id);

    match node {
        // ── Atoms ──────────────────────────────────────────────────
        ExprNode::Num(nid) => {
            let r = arena.num(*nid);
            Ok((ratio_to_bigfloat(r, prec, rm), BigFloat::new(prec)))
        }

        ExprNode::Pi => Ok((cc.pi(prec, rm).clone(), BigFloat::new(prec))),

        ExprNode::E => Ok((cc.e(prec, rm).clone(), BigFloat::new(prec))),

        ExprNode::ImaginaryUnit => Ok((BigFloat::new(prec), BigFloat::from_i32(1, prec))),

        // ── Named constants (arbitrary precision) ───────────────────────────
        ExprNode::EulerGamma => {
            debug!(prec, "evalf: EulerGamma via Brent–McMillan");
            let g = arb_euler_gamma(prec, rm, cc)?;
            Ok((g, BigFloat::new(prec)))
        }
        ExprNode::Catalan => {
            debug!(prec, "evalf: Catalan via Ramanujan series");
            let g = arb_catalan(prec, rm, cc)?;
            Ok((g, BigFloat::new(prec)))
        }
        ExprNode::GoldenRatio => {
            // φ = (1 + √5)/2, exactly via a single square root.
            let wp = prec + 16;
            let five = BigFloat::from_i32(5, wp);
            let sqrt5 = five.sqrt(wp, rm);
            let one = BigFloat::from_i32(1, wp);
            let two = BigFloat::from_i32(2, wp);
            let phi = one.add(&sqrt5, wp, rm).div(&two, prec, rm);
            Ok((phi, BigFloat::new(prec)))
        }

        // ── Complex analysis ────────────────────────────────────────────
        ExprNode::Re(inner) => {
            let val = get_cached(cache, *inner)?;
            Ok((val.0.clone(), BigFloat::new(prec)))
        }
        ExprNode::Im(inner) => {
            let val = get_cached(cache, *inner)?;
            Ok((val.1.clone(), BigFloat::new(prec)))
        }
        ExprNode::Conjugate(inner) => {
            let val = get_cached(cache, *inner)?;
            Ok((val.0.clone(), val.1.neg()))
        }
        ExprNode::Arg(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.0.is_zero() && val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "arg(0) is undefined".into(),
                });
            }
            Ok((atan2_bf(&val.1, &val.0, prec, rm, cc), BigFloat::new(prec)))
        }

        // ── Special functions (0.2) ──────────────────────────────────────
        ExprNode::Si(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Si of complex argument not yet supported in evalf".into(),
                });
            }
            debug!(prec, "evalf: Si via series/asymptotic");
            let (si, _ci) = arb_si_ci(&val.0, true, prec, rm, cc)?;
            Ok((si, BigFloat::new(prec)))
        }
        ExprNode::Ci(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Ci of complex argument not yet supported in evalf".into(),
                });
            }
            if val.0.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Ci(0) is -∞".into(),
                });
            }
            debug!(prec, "evalf: Ci via series/asymptotic");
            let x_abs = val.0.abs();
            let (_si, ci) = arb_si_ci(&x_abs, false, prec, rm, cc)?;
            if val.0.is_negative() {
                // Ci(-x) = Ci(x) + iπ (principal branch of the logarithm).
                Ok((ci, cc.pi(prec, rm).clone()))
            } else {
                Ok((ci, BigFloat::new(prec)))
            }
        }
        ExprNode::Ei(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Ei of complex argument not yet supported in evalf".into(),
                });
            }
            if val.0.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Ei(0) is -∞".into(),
                });
            }
            debug!(prec, "evalf: Ei via series/asymptotic");
            let r = arb_ei(&val.0, prec, rm, cc)?;
            Ok((r, BigFloat::new(prec)))
        }
        ExprNode::Li(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() || val.0.is_negative() {
                return Err(SymplexError::Unevaluable {
                    reason: "li of complex or negative argument not yet supported in evalf".into(),
                });
            }
            if val.0.is_zero() {
                return Ok(c_zero(prec));
            }
            let wp = prec + 32;
            let ln_x = val.0.ln(wp, rm, cc);
            if ln_x.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "li(1) is -∞".into(),
                });
            }
            debug!(prec, "evalf: li via Ei(ln x)");
            let r = arb_ei(&ln_x, prec, rm, cc)?;
            Ok((r, BigFloat::new(prec)))
        }
        ExprNode::Zeta(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "zeta of complex argument not yet supported in evalf".into(),
                });
            }
            debug!(prec, "evalf: zeta via Borwein / functional equation");
            let r = arb_zeta(&val.0, prec, rm, cc)?;
            Ok((r, BigFloat::new(prec)))
        }
        ExprNode::Polygamma(n_id, x_id) => {
            let n_val = get_cached(cache, *n_id)?;
            let x_val = get_cached(cache, *x_id)?;
            if !n_val.1.is_zero() || !x_val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "polygamma of complex arguments not yet supported in evalf".into(),
                });
            }
            let n_f = bigfloat_to_f64(&n_val.0, rm, cc)?;
            let n_round = n_f.round();
            if (n_f - n_round).abs() > 1e-12 || !(0.0..=10_000.0).contains(&n_round) {
                return Err(SymplexError::Unevaluable {
                    reason: "polygamma order must be a non-negative integer".into(),
                });
            }
            let n = n_round as u32;
            if n == 0 {
                let r = arb_digamma(&x_val.0, prec, rm, cc)?;
                return Ok((r, BigFloat::new(prec)));
            }
            debug!(
                prec,
                n, "evalf: polygamma via recurrence + asymptotic series"
            );
            let r = arb_polygamma(n, &x_val.0, prec, rm, cc)?;
            Ok((r, BigFloat::new(prec)))
        }
        ExprNode::KroneckerDelta(i_id, j_id) => {
            let iv = get_cached(cache, *i_id)?;
            let jv = get_cached(cache, *j_id)?;
            let d_re = iv.0.sub(&jv.0, prec, rm);
            let d_im = iv.1.sub(&jv.1, prec, rm);
            if d_re.is_zero() && d_im.is_zero() {
                Ok(c_one(prec))
            } else {
                Ok(c_zero(prec))
            }
        }

        ExprNode::PhysicalConstant(_, value_id) => {
            // Recursively evaluate the stored exact value to a float.
            eval_node_or_subtree(arena, *value_id, cache, prec, rm, cc)
        }

        ExprNode::Infinity | ExprNode::NegInfinity | ExprNode::ComplexInfinity => {
            Err(SymplexError::Unevaluable {
                reason: "cannot evaluate infinity to finite precision".into(),
            })
        }

        ExprNode::Factorial(_) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate symbolic factorial; call eval() first to reduce"
                .into(),
        }),

        ExprNode::Binomial(n_id, k_id) => {
            debug!(prec, "evalf: Binomial via arbitrary-precision Gamma");
            let n_val = get_cached(cache, *n_id)?;
            let k_val = get_cached(cache, *k_id)?;
            if !n_val.1.is_zero() || !k_val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Binomial of complex arguments not yet supported in evalf".into(),
                });
            }
            // C(n, k) = Gamma(n+1) / (Gamma(k+1) * Gamma(n-k+1))
            let one_bf = BigFloat::from_i32(1, prec);
            let n_plus_1 = n_val.0.add(&one_bf, prec, rm);
            let k_plus_1 = k_val.0.add(&one_bf, prec, rm);
            let n_minus_k = n_val.0.sub(&k_val.0, prec, rm);
            let n_minus_k_plus_1 = n_minus_k.add(&one_bf, prec, rm);

            let g_numer = arb_gamma_real(&n_plus_1, prec, rm, cc)?;
            let g_k = arb_gamma_real(&k_plus_1, prec, rm, cc)?;
            let g_nmk = arb_gamma_real(&n_minus_k_plus_1, prec, rm, cc)?;
            let g_denom = g_k.mul(&g_nmk, prec, rm);
            if g_denom.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Binomial denominator is zero (pole in Gamma)".into(),
                });
            }
            let result = g_numer.div(&g_denom, prec, rm);
            Ok((result, BigFloat::new(prec)))
        }

        ExprNode::Gamma(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Gamma of complex argument not yet supported in evalf".into(),
                });
            }
            debug!(prec, "evalf: Gamma via Stirling series");
            let result = arb_gamma_real(&val.0, prec, rm, cc)?;
            Ok((result, BigFloat::new(prec)))
        }

        ExprNode::LogGamma(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "LogGamma of complex argument not yet supported in evalf".into(),
                });
            }
            // Use arbitrary-precision Stirling series for log Γ(x).
            // For positive x, compute directly; for negative x, use
            // reflection: log Γ(x) = log(π) − log(sin(πx)) − log Γ(1−x).
            let x_approx = bigfloat_to_f64(&val.0, rm, cc)?;
            if x_approx <= 0.0 {
                let rounded = x_approx.round();
                if (rounded - x_approx).abs() < 1e-12 && rounded <= 0.0 {
                    return Err(SymplexError::Unevaluable {
                        reason: "LogGamma at non-positive integer pole".into(),
                    });
                }
                // Reflection: ln Γ(x) = ln π − ln|sin(πx)| − ln Γ(1−x)
                let one = BigFloat::from_i32(1, prec);
                let one_minus_x = one.sub(&val.0, prec, rm);
                let log_gamma_1mx = stirling_log_gamma(&one_minus_x, prec, rm, cc)?;
                let pi_val = cc.pi(prec, rm).clone();
                let pi_x = pi_val.mul(&val.0, prec, rm);
                let sin_pi_x = pi_x.sin(prec, rm, cc);
                let abs_sin = sin_pi_x.abs();
                if abs_sin.is_zero() {
                    return Err(SymplexError::Unevaluable {
                        reason: "LogGamma at non-positive integer pole".into(),
                    });
                }
                let ln_pi = cc.pi(prec, rm).clone().ln(prec, rm, cc);
                let ln_abs_sin = abs_sin.ln(prec, rm, cc);
                let result = ln_pi
                    .sub(&ln_abs_sin, prec, rm)
                    .sub(&log_gamma_1mx, prec, rm);
                Ok((result, BigFloat::new(prec)))
            } else {
                debug!(prec, "evalf: LogGamma via Stirling series");
                let result = stirling_log_gamma(&val.0, prec, rm, cc)?;
                Ok((result, BigFloat::new(prec)))
            }
        }

        ExprNode::Digamma(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Digamma of complex argument not yet supported in evalf".into(),
                });
            }
            let result = arb_digamma(&val.0, prec, rm, cc)?;
            Ok((result, BigFloat::new(prec)))
        }

        ExprNode::Erf(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "erf of complex argument not yet supported in evalf".into(),
                });
            }
            debug!(prec, "evalf: erf via arbitrary-precision series");
            let result = arb_erf(&val.0, prec, rm, cc)?;
            Ok((result, BigFloat::new(prec)))
        }

        ExprNode::Erfc(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "erfc of complex argument not yet supported in evalf".into(),
                });
            }
            debug!(prec, "evalf: erfc via arbitrary-precision series");
            let erf_val = arb_erf(&val.0, prec, rm, cc)?;
            let one = BigFloat::from_i32(1, prec);
            let result = one.sub(&erf_val, prec, rm);
            Ok((result, BigFloat::new(prec)))
        }

        ExprNode::LambertW(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "LambertW of complex argument not yet supported in evalf".into(),
                });
            }
            debug!(prec, "evalf: LambertW via Halley iteration");
            let result = arb_lambert_w(&val.0, prec, rm, cc)?;
            Ok((result, BigFloat::new(prec)))
        }

        ExprNode::Beta(a_id, b_id) => {
            let a_val = get_cached(cache, *a_id)?;
            let b_val = get_cached(cache, *b_id)?;
            if !a_val.1.is_zero() || !b_val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Beta of complex arguments not yet supported in evalf".into(),
                });
            }
            // Beta(a,b) = Gamma(a)*Gamma(b)/Gamma(a+b) at arbitrary precision.
            debug!(prec, "evalf: Beta via arbitrary-precision Gamma");
            let ga = arb_gamma_real(&a_val.0, prec, rm, cc)?;
            let gb = arb_gamma_real(&b_val.0, prec, rm, cc)?;
            let a_plus_b = a_val.0.add(&b_val.0, prec, rm);
            let gab = arb_gamma_real(&a_plus_b, prec, rm, cc)?;
            if gab.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Beta function: Gamma(a+b) is zero".into(),
                });
            }
            let numer = ga.mul(&gb, prec, rm);
            let result = numer.div(&gab, prec, rm);
            Ok((result, BigFloat::new(prec)))
        }

        ExprNode::NaN => Err(SymplexError::Unevaluable {
            reason: "cannot evaluate NaN".into(),
        }),

        ExprNode::Symbol(sid) => {
            let name = arena.symbol_name(*sid);
            Err(SymplexError::FreeSymbol {
                name: name.to_owned(),
            })
        }

        // ── Add ────────────────────────────────────────────────────
        ExprNode::Add(children) => {
            let mut sum = c_zero(prec);
            for &child in children.iter() {
                let val = get_cached(cache, child)?;
                sum = c_add(&sum, val, prec, rm);
            }
            Ok(sum)
        }

        // ── Mul ────────────────────────────────────────────────────
        ExprNode::Mul(children) => {
            let mut product = c_one(prec);
            for &child in children.iter() {
                let val = get_cached(cache, child)?;
                product = c_mul(&product, val, prec, rm);
            }
            Ok(product)
        }

        // ── Pow ────────────────────────────────────────────────────
        ExprNode::Pow(base, exp) => {
            let b = get_cached(cache, *base)?;
            let e = get_cached(cache, *exp)?;

            let b_is_real = b.1.is_zero();
            let e_is_real = e.1.is_zero();

            // Optimization: use sqrt for x^(1/2) with non-negative real base.
            if b_is_real
                && !b.0.is_negative()
                && let ExprNode::Num(nid) = arena.node(*exp)
            {
                let r = arena.num(*nid);
                if *r == Ratio::new(BigInt::from(1), BigInt::from(2)) {
                    return Ok((b.0.sqrt(prec, rm), BigFloat::new(prec)));
                }
            }

            // Optimization: small integer exponents.
            if let Some(n) = try_as_small_int(arena, *exp) {
                if b_is_real {
                    // Real base with integer exponent: use real powi.
                    if n >= 0 {
                        return Ok((b.0.powi(n as usize, prec, rm), BigFloat::new(prec)));
                    } else {
                        let pow_pos = b.0.powi((-n) as usize, prec, rm);
                        let one_bf = BigFloat::from_i32(1, prec);
                        return Ok((one_bf.div(&pow_pos, prec, rm), BigFloat::new(prec)));
                    }
                } else {
                    // Complex base with integer exponent: use c_powi.
                    if n >= 0 {
                        return Ok(c_powi(b, n as usize, prec, rm));
                    } else {
                        let pow_pos = c_powi(b, (-n) as usize, prec, rm);
                        let one_c = c_one(prec);
                        return Ok(c_div(&one_c, &pow_pos, prec, rm));
                    }
                }
            }

            // Real positive base with real exponent: exp(e·ln b).
            // (`BigFloat::pow` is avoided — it hangs on exactly
            // representable results such as 4^(1/2); see `bf_pow`.)
            if b_is_real && e_is_real && b.0.is_positive() {
                return Ok((bf_pow(&b.0, &e.0, prec, rm, cc), BigFloat::new(prec)));
            }

            // General complex power: b^e = exp(e * ln(b)).
            Ok(c_pow(b, e, prec, rm, cc))
        }

        // ── Neg ────────────────────────────────────────────────────
        ExprNode::Neg(inner) => {
            let val = get_cached(cache, *inner)?;
            Ok(c_neg(val, prec, rm))
        }

        // ── Trig ───────────────────────────────────────────────────
        ExprNode::Sin(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.sin(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_sin(val, prec, rm, cc))
            }
        }

        ExprNode::Cos(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.cos(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_cos(val, prec, rm, cc))
            }
        }

        ExprNode::Tan(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.tan(prec, rm, cc), BigFloat::new(prec)))
            } else {
                let s = c_sin(val, prec, rm, cc);
                let c = c_cos(val, prec, rm, cc);
                Ok(c_div(&s, &c, prec, rm))
            }
        }

        // ── Exp / Ln ───────────────────────────────────────────────
        ExprNode::Exp(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.exp(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_exp(val, prec, rm, cc))
            }
        }

        ExprNode::Ln(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() && val.0.is_positive() {
                Ok((val.0.ln(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_ln(val, prec, rm, cc))
            }
        }

        // ── Abs ────────────────────────────────────────────────────
        ExprNode::Abs(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.abs(), BigFloat::new(prec)))
            } else {
                Ok((c_abs(val, prec, rm), BigFloat::new(prec)))
            }
        }

        // ── Inverse trig ──────────────────────────────────────────
        ExprNode::Asin(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.asin(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_asin(val, prec, rm, cc))
            }
        }

        ExprNode::Acos(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.acos(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_acos(val, prec, rm, cc))
            }
        }

        ExprNode::Atan(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.atan(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_atan(val, prec, rm, cc))
            }
        }

        // ── Atan2 ─────────────────────────────────────────────────
        ExprNode::Atan2(y_id, x_id) => {
            let yv = get_cached(cache, *y_id)?;
            let xv = get_cached(cache, *x_id)?;
            // Both args must be real for a real atan2 result.
            if yv.1.is_zero() && xv.1.is_zero() {
                Ok((atan2_bf(&yv.0, &xv.0, prec, rm, cc), BigFloat::new(prec)))
            } else {
                // Complex fallback: atan2 is not standard for complex args;
                // compute as -i * ln((x + iy) / sqrt(x² + y²))
                // but for now just use the real parts as a best-effort.
                Err(SymplexError::Unevaluable {
                    reason: "atan2 with complex arguments is not supported".into(),
                })
            }
        }

        // ── Hyperbolic ────────────────────────────────────────────
        ExprNode::Sinh(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.sinh(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_sinh(val, prec, rm, cc))
            }
        }

        ExprNode::Cosh(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.cosh(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_cosh(val, prec, rm, cc))
            }
        }

        ExprNode::Tanh(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.tanh(prec, rm, cc), BigFloat::new(prec)))
            } else {
                let s = c_sinh(val, prec, rm, cc);
                let c = c_cosh(val, prec, rm, cc);
                Ok(c_div(&s, &c, prec, rm))
            }
        }

        // ── Inverse hyperbolic ────────────────────────────────────
        ExprNode::Asinh(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.asinh(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_asinh(val, prec, rm, cc))
            }
        }

        ExprNode::Acosh(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.acosh(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_acosh(val, prec, rm, cc))
            }
        }

        ExprNode::Atanh(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                Ok((val.0.atanh(prec, rm, cc), BigFloat::new(prec)))
            } else {
                Ok(c_atanh(val, prec, rm, cc))
            }
        }

        // ── Sign ───────────────────────────────────────────────────
        ExprNode::Sign(inner) => {
            let val = get_cached(cache, *inner)?;
            if val.1.is_zero() {
                // Real case: sign returns -1, 0, or 1
                if val.0.is_zero() {
                    Ok(c_zero(prec))
                } else if val.0.is_negative() {
                    Ok((BigFloat::from_i32(-1, prec), BigFloat::new(prec)))
                } else {
                    Ok((BigFloat::from_i32(1, prec), BigFloat::new(prec)))
                }
            } else {
                // Complex case: sign(z) = z / |z|
                let modulus = c_abs(val, prec, rm);
                if modulus.is_zero() {
                    Ok(c_zero(prec))
                } else {
                    let denom = (modulus, BigFloat::new(prec));
                    Ok(c_div(val, &denom, prec, rm))
                }
            }
        }

        // ── Floor ──────────────────────────────────────────────────
        ExprNode::Floor(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "floor of complex number is not defined".into(),
                });
            }
            // BigFloat::int() returns the integer part (truncation toward zero).
            // floor(x) = x if x is integer, else truncate toward -infinity.
            let truncated = val.0.int();
            let result = if val.0.is_negative() && truncated != val.0 {
                // For negative non-integers, floor = trunc - 1
                truncated.sub(&BigFloat::from_i32(1, prec), prec, rm)
            } else {
                truncated
            };
            Ok((result, BigFloat::new(prec)))
        }

        // ── Ceiling ────────────────────────────────────────────────
        ExprNode::Ceiling(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "ceiling of complex number is not defined".into(),
                });
            }
            let truncated = val.0.int();
            let result = if val.0.is_positive() && truncated != val.0 {
                // For positive non-integers, ceil = trunc + 1
                truncated.add(&BigFloat::from_i32(1, prec), prec, rm)
            } else {
                truncated
            };
            Ok((result, BigFloat::new(prec)))
        }

        // ── Min ────────────────────────────────────────────────────
        ExprNode::Min(children) => {
            if children.is_empty() {
                return Err(SymplexError::Unevaluable {
                    reason: "min of empty set".into(),
                });
            }
            let mut best = get_cached(cache, children[0])?.clone();
            if !best.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "min requires real arguments".into(),
                });
            }
            for &child in &children[1..] {
                let val = get_cached(cache, child)?;
                if !val.1.is_zero() {
                    return Err(SymplexError::Unevaluable {
                        reason: "min requires real arguments".into(),
                    });
                }
                if val.0.sub(&best.0, prec, rm).is_negative() {
                    best = val.clone();
                }
            }
            Ok(best)
        }

        // ── Max ────────────────────────────────────────────────────
        ExprNode::Max(children) => {
            if children.is_empty() {
                return Err(SymplexError::Unevaluable {
                    reason: "max of empty set".into(),
                });
            }
            let mut best = get_cached(cache, children[0])?.clone();
            if !best.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "max requires real arguments".into(),
                });
            }
            for &child in &children[1..] {
                let val = get_cached(cache, child)?;
                if !val.1.is_zero() {
                    return Err(SymplexError::Unevaluable {
                        reason: "max requires real arguments".into(),
                    });
                }
                if val.0.sub(&best.0, prec, rm).is_positive() {
                    best = val.clone();
                }
            }
            Ok(best)
        }

        // ── Heaviside ──────────────────────────────────────────────
        ExprNode::Heaviside(inner) => {
            let val = get_cached(cache, *inner)?;
            if !val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Heaviside of complex argument not supported in evalf".into(),
                });
            }
            if val.0.is_zero() {
                // H(0) = 0.5
                Ok((BigFloat::from_f64(0.5, prec), BigFloat::new(prec)))
            } else if val.0.is_negative() {
                Ok(c_zero(prec))
            } else {
                Ok((BigFloat::from_i32(1, prec), BigFloat::new(prec)))
            }
        }

        // ── DiracDelta ─────────────────────────────────────────────
        // Distributional: zero everywhere except at a single point of measure zero.
        ExprNode::DiracDelta(inner) => {
            let _val = get_cached(cache, *inner)?;
            // Numerically, δ(x) = 0 for all representable floats
            Ok(c_zero(prec))
        }

        // ── Apply-based special functions ───────────────────────────
        ExprNode::Apply(sid, args) => {
            let name = arena.symbol_name(*sid);
            match name {
                "besselj" if args.len() == 2 => {
                    let order = get_cached(cache, args[0])?;
                    let arg = get_cached(cache, args[1])?;
                    if !order.1.is_zero() || !arg.1.is_zero() {
                        return Err(SymplexError::Unevaluable {
                            reason: "Bessel of complex argument not yet supported in evalf".into(),
                        });
                    }
                    tracing::debug!(prec, "evalf: BesselJ via series/asymptotic");
                    let result = arb_bessel_j(&order.0, &arg.0, prec, rm, cc)?;
                    Ok((result, BigFloat::new(prec)))
                }
                "bessely" if args.len() == 2 => {
                    let order = get_cached(cache, args[0])?;
                    let arg = get_cached(cache, args[1])?;
                    if !order.1.is_zero() || !arg.1.is_zero() {
                        return Err(SymplexError::Unevaluable {
                            reason: "Bessel of complex argument not yet supported in evalf".into(),
                        });
                    }
                    tracing::debug!(prec, "evalf: BesselY via series/asymptotic");
                    let result = arb_bessel_y(&order.0, &arg.0, prec, rm, cc)?;
                    Ok((result, BigFloat::new(prec)))
                }
                "besseli" if args.len() == 2 => {
                    let order = get_cached(cache, args[0])?;
                    let arg = get_cached(cache, args[1])?;
                    if !order.1.is_zero() || !arg.1.is_zero() {
                        return Err(SymplexError::Unevaluable {
                            reason: "Bessel of complex argument not yet supported in evalf".into(),
                        });
                    }
                    tracing::debug!(prec, "evalf: BesselI via ascending series");
                    let result = arb_bessel_i(&order.0, &arg.0, prec, rm, cc)?;
                    Ok((result, BigFloat::new(prec)))
                }
                "besselk" if args.len() == 2 => {
                    let order = get_cached(cache, args[0])?;
                    let arg = get_cached(cache, args[1])?;
                    if !order.1.is_zero() || !arg.1.is_zero() {
                        return Err(SymplexError::Unevaluable {
                            reason: "Bessel of complex argument not yet supported in evalf".into(),
                        });
                    }
                    tracing::debug!(prec, "evalf: BesselK via series/asymptotic");
                    let result = arb_bessel_k(&order.0, &arg.0, prec, rm, cc)?;
                    Ok((result, BigFloat::new(prec)))
                }
                "legendre" | "chebyshev_t" | "chebyshev_u" | "hermite" | "laguerre"
                    if args.len() == 2 =>
                {
                    let n_val = get_cached(cache, args[0])?;
                    let x_val = get_cached(cache, args[1])?;
                    if !n_val.1.is_zero() || !x_val.1.is_zero() {
                        return Err(SymplexError::Unevaluable {
                            reason: "orthogonal polynomial of complex arguments not supported"
                                .into(),
                        });
                    }
                    let n_f = bigfloat_to_f64(&n_val.0, rm, cc)?;
                    let n_round = n_f.round();
                    if (n_f - n_round).abs() > 1e-12 || !(0.0..=1.0e7).contains(&n_round) {
                        return Err(SymplexError::Unevaluable {
                            reason: format!("{name}: degree must be a non-negative integer"),
                        });
                    }
                    tracing::debug!(
                        prec,
                        n = n_round,
                        "evalf: orthogonal polynomial via recurrence"
                    );
                    let kind = match name {
                        "legendre" => OrthoPoly::Legendre,
                        "chebyshev_t" => OrthoPoly::ChebyshevT,
                        "chebyshev_u" => OrthoPoly::ChebyshevU,
                        "hermite" => OrthoPoly::Hermite,
                        _ => OrthoPoly::Laguerre,
                    };
                    let result = arb_orthopoly(kind, n_round as u64, &x_val.0, prec, rm);
                    Ok((result, BigFloat::new(prec)))
                }
                // ── More special functions (0.9) ──
                n if crate::transforms::eval::is_special_09(n) => {
                    eval_special_09(n, args, cache, prec, rm, cc)
                }
                _ => Err(SymplexError::Unevaluable {
                    reason: format!("cannot evaluate function '{name}'"),
                }),
            }
        }

        ExprNode::Derivative(_, _) => Err(SymplexError::Unevaluable {
            reason: "cannot evaluate unevaluated derivative".into(),
        }),

        ExprNode::Integral(_, _) => Err(SymplexError::Unevaluable {
            reason: "cannot evaluate unevaluated integral".into(),
        }),

        // ── DefiniteIntegral: f64 Gauss–Kronrod quadrature ───────────────
        // The body is compiled to a stack-VM function of the integration
        // variable (no arena mutation), the bounds come from the cache
        // (±∞ allowed), and the adaptive G7/K15 rule integrates it.  The
        // result is an `f64` widened to the working precision — see the
        // module docs for the precision contract.
        ExprNode::DefiniteIntegral(body_id, var_id, lo_id, hi_id) => {
            debug!("evalf: DefiniteIntegral — f64 Gauss–Kronrod quadrature");
            let var_name = match arena.node(*var_id) {
                ExprNode::Symbol(sid) => arena.symbol_name(*sid).to_string(),
                _ => {
                    return Err(SymplexError::Unevaluable {
                        reason: "integration variable of a definite integral must be a symbol"
                            .into(),
                    });
                }
            };
            let a = definite_bound_f64(arena, cache, *lo_id, rm, cc)?;
            let b = definite_bound_f64(arena, cache, *hi_id, rm, cc)?;
            let func = crate::output::lambdify::compile_raw(arena, *body_id, &[&var_name])
                .map_err(|e| SymplexError::Unevaluable {
                    reason: format!(
                        "definite integral body cannot be compiled for quadrature: {e}"
                    ),
                })?;
            let f = |t: f64| func(&[t]);
            let opts = crate::calculus::definite::QuadOpts::default();
            let (value, err) = crate::calculus::definite::quadrature(&f, a, b, &opts)?;
            let tol = opts.abs_tol.max(opts.rel_tol * value.abs());
            if err > 1e3 * tol {
                return Err(SymplexError::ComputationFailed {
                    operation: "evalf",
                    reason: format!(
                        "quadrature of the definite integral did not converge: estimate {value} \
                         with error {err:e} (integral may diverge)"
                    ),
                });
            }
            debug!(value, err, "evalf: DefiniteIntegral evaluated");
            Ok((BigFloat::from_f64(value, prec), BigFloat::new(prec)))
        }

        ExprNode::Sum(body_id, var_id, lo_id, hi_id) => {
            debug!("evalf: Sum — attempting finite evaluation");
            let lo_val = get_cached(cache, *lo_id)?;
            let hi_val = get_cached(cache, *hi_id)?;
            if !lo_val.1.is_zero() || !hi_val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Sum bounds must be real".into(),
                });
            }
            let lo_f = bigfloat_to_f64(&lo_val.0, rm, cc)?;
            let hi_f = bigfloat_to_f64(&hi_val.0, rm, cc)?;
            let lo_i = lo_f.round() as i64;
            let hi_i = hi_f.round() as i64;
            if (lo_f - lo_i as f64).abs() > 1e-9 || (hi_f - hi_i as f64).abs() > 1e-9 {
                return Err(SymplexError::Unevaluable {
                    reason: "Sum bounds are not integers".into(),
                });
            }
            if hi_i - lo_i > 10_000 {
                return Err(SymplexError::Unevaluable {
                    reason: "Sum range too large for numerical evaluation (> 10000 terms)".into(),
                });
            }
            let mut acc = c_zero(prec);
            for i in lo_i..=hi_i {
                let sub_val = (BigFloat::from_f64(i as f64, prec), BigFloat::new(prec));
                let term =
                    evalf_subtree_with_sub(arena, *body_id, *var_id, &sub_val, prec, rm, cc)?;
                acc = c_add(&acc, &term, prec, rm);
            }
            debug!(lo = lo_i, hi = hi_i, "evalf: Sum evaluated");
            Ok(acc)
        }

        ExprNode::Product_(body_id, var_id, lo_id, hi_id) => {
            debug!("evalf: Product — attempting finite evaluation");
            let lo_val = get_cached(cache, *lo_id)?;
            let hi_val = get_cached(cache, *hi_id)?;
            if !lo_val.1.is_zero() || !hi_val.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "Product bounds must be real".into(),
                });
            }
            let lo_f = bigfloat_to_f64(&lo_val.0, rm, cc)?;
            let hi_f = bigfloat_to_f64(&hi_val.0, rm, cc)?;
            let lo_i = lo_f.round() as i64;
            let hi_i = hi_f.round() as i64;
            if (lo_f - lo_i as f64).abs() > 1e-9 || (hi_f - hi_i as f64).abs() > 1e-9 {
                return Err(SymplexError::Unevaluable {
                    reason: "Product bounds are not integers".into(),
                });
            }
            if hi_i - lo_i > 10_000 {
                return Err(SymplexError::Unevaluable {
                    reason: "Product range too large for numerical evaluation (> 10000 terms)"
                        .into(),
                });
            }
            let mut acc = c_one(prec);
            for i in lo_i..=hi_i {
                let sub_val = (BigFloat::from_f64(i as f64, prec), BigFloat::new(prec));
                let term =
                    evalf_subtree_with_sub(arena, *body_id, *var_id, &sub_val, prec, rm, cc)?;
                acc = c_mul(&acc, &term, prec, rm);
            }
            debug!(lo = lo_i, hi = hi_i, "evalf: Product evaluated");
            Ok(acc)
        }

        ExprNode::Piecewise(pairs) => {
            debug!(
                branches = pairs.len(),
                "evalf: Piecewise — evaluating conditions"
            );
            // Branches are examined in order.  The first condition that is
            // decidedly true selects its value; a decidedly false condition
            // is skipped.  An *undecided* condition is an error — falling
            // through to a later `True` branch would be silently wrong.
            for &(value_id, cond_id) in pairs.iter() {
                match decide_condition(arena, cond_id, cache, prec, rm) {
                    Some(true) => {
                        return eval_node_or_subtree(arena, value_id, cache, prec, rm, cc);
                    }
                    Some(false) => continue,
                    None => {
                        return Err(SymplexError::Unevaluable {
                            reason: format!(
                                "cannot evaluate piecewise: condition `{}` is undecided",
                                arena.display(cond_id)
                            ),
                        });
                    }
                }
            }
            Err(SymplexError::Unevaluable {
                reason: "cannot evaluate piecewise: every condition is false".into(),
            })
        }

        ExprNode::BoolTrue
        | ExprNode::BoolFalse
        | ExprNode::Gt(_, _)
        | ExprNode::Ge(_, _)
        | ExprNode::Eq_(_, _)
        | ExprNode::Ne(_, _)
        | ExprNode::And(_)
        | ExprNode::Or(_)
        | ExprNode::Not(_) => Err(SymplexError::Unevaluable {
            reason: "boolean expression".into(),
        }),

        ExprNode::EmptySet
        | ExprNode::UniversalSet
        | ExprNode::Interval(_, _, _)
        | ExprNode::FiniteSet(_)
        | ExprNode::SetUnion(_)
        | ExprNode::SetIntersection(_)
        | ExprNode::SetComplement(_, _) => Err(SymplexError::Unevaluable {
            reason: "set-valued expressions cannot be numerically evaluated".into(),
        }),

        ExprNode::Limit(_, _, _) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate unevaluated Limit".into(),
        }),

        ExprNode::Series(_, _, _, _) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate unevaluated Series".into(),
        }),

        ExprNode::LaplaceTransform(_, _, _) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate unevaluated LaplaceTransform".into(),
        }),

        ExprNode::InverseLaplaceTransform(_, _, _) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate unevaluated InverseLaplaceTransform".into(),
        }),

        ExprNode::Residue(_, _, _) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate unevaluated Residue".into(),
        }),

        ExprNode::RootOf(poly_id, idx_id) => {
            // ── Extract the index as a non-negative integer ────────
            let idx: usize = match arena.node(*idx_id) {
                ExprNode::Num(nid) => {
                    let r = arena.num(*nid);
                    if r.is_integer() {
                        let n = r.to_integer();
                        // Convert BigInt → i64 → usize safely.
                        let n_i64: i64 = n.try_into().map_err(|_| SymplexError::Unevaluable {
                            reason: "RootOf index out of range".into(),
                        })?;
                        if n_i64 < 0 {
                            return Err(SymplexError::Unevaluable {
                                reason: "RootOf index must be non-negative".into(),
                            });
                        }
                        n_i64 as usize
                    } else {
                        return Err(SymplexError::Unevaluable {
                            reason: "RootOf index must be an integer".into(),
                        });
                    }
                }
                _ => {
                    return Err(SymplexError::Unevaluable {
                        reason: "RootOf index must be numeric".into(),
                    });
                }
            };

            // ── Identify the variable in the polynomial ───────────
            let syms = walk::free_symbols(arena, *poly_id);
            if syms.is_empty() {
                return Err(SymplexError::Unevaluable {
                    reason: "RootOf polynomial has no variables".into(),
                });
            }
            let var_id = syms[0];

            // ── Convert expression → dense Poly over ℚ ───────────
            let poly = match crate::poly::polybridge::expr_to_poly(arena, *poly_id, var_id) {
                Some(p) => p,
                None => {
                    return Err(SymplexError::Unevaluable {
                        reason: "could not convert RootOf expression to polynomial".into(),
                    });
                }
            };

            // ── Find all roots via Aberth's method ────────────────
            // Aberth handles both real and complex roots simultaneously
            // with cubic convergence.
            let roots = crate::poly::roots::aberth_roots(&poly, prec + 64, 200);

            if idx >= roots.len() {
                return Err(SymplexError::Unevaluable {
                    reason: format!(
                        "RootOf index {} exceeds the {} root(s) found",
                        idx,
                        roots.len()
                    ),
                });
            }

            let (re, im) = &roots[idx];
            Ok((re.clone(), im.clone()))
        }

        // ── RootSum: numerical evaluation via root-finding + summation ──
        //
        // RootSum(poly, body, sumvar) = Σ_{α: poly(α)=0} body(α, x).
        // We find the roots of poly numerically, substitute each into body,
        // evaluate, and sum.  If the body contains free variables other than
        // sumvar, this will fail (those variables must be substituted first).
        ExprNode::RootSum(poly, body, sumvar) => {
            // RootSum(poly, body, sumvar) = Σ_{α: poly(α)=0} body(α, x).
            //
            // Strategy:
            // 1. Convert the polynomial expression to a Poly via expr_to_poly
            //    (which takes &Arena — no mutation needed).
            // 2. Find ALL roots numerically via Aberth's method.
            // 3. For each root α_k, evaluate body with sumvar = α_k using
            //    evalf_subtree_with_sub (same mechanism as Sum evaluation).
            // 4. Sum all contributions.
            //
            // The imaginary parts should cancel for real-valued integrals.

            let poly_id = *poly;
            let body_id = *body;
            let sumvar_id = *sumvar;

            tracing::debug!("evalf: RootSum — attempting numerical evaluation via Aberth roots");

            // Step 1: Convert polynomial expression to Poly.
            let poly_obj = crate::poly::polybridge::expr_to_poly(arena, poly_id, sumvar_id)
                .ok_or_else(|| SymplexError::Unevaluable {
                    reason: "RootSum: cannot convert polynomial expression to Poly \
                             (may contain free symbols)"
                        .into(),
                })?;

            let poly_deg = poly_obj.degree().unwrap_or(0);
            tracing::debug!(
                degree = poly_deg,
                "evalf: RootSum — polynomial extracted, finding roots"
            );

            // Step 2: Find all roots via Aberth's method.
            // Use the same working precision as the rest of the evalf computation.
            let roots = crate::poly::roots::aberth_roots(&poly_obj, prec, 200);

            if roots.len() != poly_deg {
                tracing::debug!(
                    expected = poly_deg,
                    found = roots.len(),
                    "evalf: RootSum — Aberth returned fewer roots than expected"
                );
            }

            // Step 3+4: Evaluate body at each root and sum.
            let mut sum = c_zero(prec);
            for (k, root) in roots.iter().enumerate() {
                let term = evalf_subtree_with_sub(arena, body_id, sumvar_id, root, prec, rm, cc)?;
                tracing::trace!(root_idx = k, "evalf: RootSum — evaluated body at root");
                sum = c_add(&sum, &term, prec, rm);
            }

            tracing::debug!(
                n_roots = roots.len(),
                "evalf: RootSum — numerical evaluation complete"
            );
            Ok(sum)
        }

        ExprNode::DSolve(_, _, _) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate unevaluated DSolve".into(),
        }),

        ExprNode::ConditionSet(_, _) => Err(SymplexError::Unevaluable {
            reason: "cannot numerically evaluate ConditionSet".into(),
        }),
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Definite-integral helpers
// ═══════════════════════════════════════════════════════════════════════════

/// Does the tree contain a `DefiniteIntegral` node?
fn contains_definite_integral(arena: &Arena, root: ExprId) -> bool {
    let mut stack = vec![root];
    let mut seen = rustc_hash::FxHashSet::default();
    while let Some(id) = stack.pop() {
        if !seen.insert(id) {
            continue;
        }
        let node = arena.node(id);
        if matches!(node, ExprNode::DefiniteIntegral(..)) {
            return true;
        }
        node.for_each_child(|c| stack.push(c));
    }
    false
}

/// An integration bound as `f64`: `±∞` nodes directly, anything else from
/// the evaluated cache (must be real).
fn definite_bound_f64(
    arena: &Arena,
    cache: &FxHashMap<ExprId, Complex>,
    id: ExprId,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<f64, SymplexError> {
    match arena.node(id) {
        ExprNode::Infinity => Ok(f64::INFINITY),
        ExprNode::NegInfinity => Ok(f64::NEG_INFINITY),
        _ => {
            let v = get_cached(cache, id)?;
            if !v.1.is_zero() {
                return Err(SymplexError::Unevaluable {
                    reason: "integration bounds of a definite integral must be real".into(),
                });
            }
            bigfloat_to_f64(&v.0, rm, cc)
        }
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Sub-tree evaluation helpers (Sum, Product, Piecewise)
// ═══════════════════════════════════════════════════════════════════════════

/// Evaluate a sub-expression with one variable replaced by a given numeric
/// value.  Used by `Sum` and `Product_` to iterate over finite ranges.
fn evalf_subtree_with_sub(
    arena: &Arena,
    expr: ExprId,
    var_id: ExprId,
    var_value: &Complex,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<Complex, SymplexError> {
    let post_order = walk::post_order_ids(arena, expr);
    let mut sub_cache: FxHashMap<ExprId, Complex> = FxHashMap::default();
    // Pre-load the substitution so `var_id` resolves to the numeric value.
    sub_cache.insert(var_id, var_value.clone());

    for &id in &post_order {
        if sub_cache.contains_key(&id) {
            continue; // already present (e.g. the substituted variable)
        }
        let value = eval_node(arena, id, &sub_cache, prec, rm, cc)?;
        sub_cache.insert(id, value);
    }

    sub_cache
        .get(&expr)
        .cloned()
        .ok_or_else(|| SymplexError::Unevaluable {
            reason: "subtree evaluation with substitution failed".into(),
        })
}

/// Decide a boolean condition numerically using already-evaluated operands.
///
/// Relational nodes compare the cached real values of their operands;
/// `And`/`Or`/`Not` are combined with three-valued logic.  Returns `None`
/// when any needed operand is missing from the cache (free symbol, complex
/// value, unsupported node), so the caller can refuse rather than guess.
fn decide_condition(
    arena: &Arena,
    cond: ExprId,
    cache: &FxHashMap<ExprId, Complex>,
    prec: usize,
    rm: RoundingMode,
) -> Option<bool> {
    // Real-valued difference `a - b` when both operands are cached reals.
    let real_diff = |a: &ExprId, b: &ExprId| -> Option<BigFloat> {
        let (av, bv) = (cache.get(a)?, cache.get(b)?);
        if av.1.is_zero() && bv.1.is_zero() {
            Some(av.0.sub(&bv.0, prec, rm))
        } else {
            None
        }
    };

    let order = walk::post_order_ids(arena, cond);
    let mut truth: FxHashMap<ExprId, Option<bool>> = FxHashMap::default();
    for &id in &order {
        let v: Option<bool> = match arena.node(id) {
            ExprNode::BoolTrue => Some(true),
            ExprNode::BoolFalse => Some(false),
            ExprNode::Gt(a, b) => real_diff(a, b).map(|d| d.is_positive()),
            ExprNode::Ge(a, b) => real_diff(a, b).map(|d| d.is_positive() || d.is_zero()),
            ExprNode::Eq_(a, b) => real_diff(a, b).map(|d| d.is_zero()),
            ExprNode::Ne(a, b) => real_diff(a, b).map(|d| !d.is_zero()),
            ExprNode::Not(inner) => truth.get(inner).copied().flatten().map(|b| !b),
            ExprNode::And(kids) => {
                let vals: Vec<Option<bool>> = kids
                    .iter()
                    .map(|k| truth.get(k).copied().flatten())
                    .collect();
                if vals.contains(&Some(false)) {
                    Some(false)
                } else if vals.iter().all(|v| *v == Some(true)) {
                    Some(true)
                } else {
                    None
                }
            }
            ExprNode::Or(kids) => {
                let vals: Vec<Option<bool>> = kids
                    .iter()
                    .map(|k| truth.get(k).copied().flatten())
                    .collect();
                if vals.contains(&Some(true)) {
                    Some(true)
                } else if vals.iter().all(|v| *v == Some(false)) {
                    Some(false)
                } else {
                    None
                }
            }
            // Numeric operands and anything else carry no truth value.
            _ => None,
        };
        truth.insert(id, v);
    }
    truth.get(&cond).copied().flatten()
}

/// Return the cached value for `id`, or fall back to calling `eval_node`.
fn eval_node_or_subtree(
    arena: &Arena,
    id: ExprId,
    cache: &FxHashMap<ExprId, Complex>,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<Complex, SymplexError> {
    if let Some(val) = cache.get(&id) {
        return Ok(val.clone());
    }
    eval_node(arena, id, cache, prec, rm, cc)
}

// ═══════════════════════════════════════════════════════════════════════════
// Complex arithmetic helpers
// ═══════════════════════════════════════════════════════════════════════════

pub(crate) fn c_zero(prec: usize) -> Complex {
    (BigFloat::new(prec), BigFloat::new(prec))
}

pub(crate) fn c_one(prec: usize) -> Complex {
    (BigFloat::from_i32(1, prec), BigFloat::new(prec))
}

fn c_i(prec: usize) -> Complex {
    (BigFloat::new(prec), BigFloat::from_i32(1, prec))
}

pub(crate) fn c_from_real(r: BigFloat, prec: usize) -> Complex {
    (r, BigFloat::new(prec))
}

pub(crate) fn c_add(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    (a.0.add(&b.0, prec, rm), a.1.add(&b.1, prec, rm))
}

pub(crate) fn c_sub(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    (a.0.sub(&b.0, prec, rm), a.1.sub(&b.1, prec, rm))
}

pub(crate) fn c_mul(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    // (a+bi)(c+di) = (ac-bd) + (ad+bc)i
    let ac = a.0.mul(&b.0, prec, rm);
    let bd = a.1.mul(&b.1, prec, rm);
    let ad = a.0.mul(&b.1, prec, rm);
    let bc = a.1.mul(&b.0, prec, rm);
    (ac.sub(&bd, prec, rm), ad.add(&bc, prec, rm))
}

pub(crate) fn c_div(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    // (a+bi)/(c+di) = ((ac+bd) + (bc-ad)i) / (c²+d²)
    let ac = a.0.mul(&b.0, prec, rm);
    let bd = a.1.mul(&b.1, prec, rm);
    let bc = a.1.mul(&b.0, prec, rm);
    let ad = a.0.mul(&b.1, prec, rm);
    let denom =
        b.0.mul(&b.0, prec, rm)
            .add(&b.1.mul(&b.1, prec, rm), prec, rm);
    let re = ac.add(&bd, prec, rm).div(&denom, prec, rm);
    let im = bc.sub(&ad, prec, rm).div(&denom, prec, rm);
    (re, im)
}

#[allow(unused_variables)]
pub(crate) fn c_neg(a: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    (a.0.neg(), a.1.neg())
}

pub(crate) fn c_abs(a: &Complex, prec: usize, rm: RoundingMode) -> BigFloat {
    // |z| = sqrt(re² + im²)
    let re2 = a.0.mul(&a.0, prec, rm);
    let im2 = a.1.mul(&a.1, prec, rm);
    let sum = re2.add(&im2, prec, rm);
    sum.sqrt(prec, rm)
}

fn c_exp(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // exp(a+bi) = exp(a)(cos(b) + i·sin(b))
    let exp_a = z.0.exp(prec, rm, cc);
    let cos_b = z.1.cos(prec, rm, cc);
    let sin_b = z.1.sin(prec, rm, cc);
    (exp_a.mul(&cos_b, prec, rm), exp_a.mul(&sin_b, prec, rm))
}

fn c_ln(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // ln(z) = ln|z| + i·arg(z)
    let r = c_abs(z, prec, rm);
    let theta = atan2_bf(&z.1, &z.0, prec, rm, cc);
    (r.ln(prec, rm, cc), theta)
}

fn c_sin(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // sin(a+bi) = sin(a)cosh(b) + i·cos(a)sinh(b)
    let sin_a = z.0.sin(prec, rm, cc);
    let cos_a = z.0.cos(prec, rm, cc);
    let cosh_b = z.1.cosh(prec, rm, cc);
    let sinh_b = z.1.sinh(prec, rm, cc);
    (sin_a.mul(&cosh_b, prec, rm), cos_a.mul(&sinh_b, prec, rm))
}

fn c_cos(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // cos(a+bi) = cos(a)cosh(b) - i·sin(a)sinh(b)
    let cos_a = z.0.cos(prec, rm, cc);
    let sin_a = z.0.sin(prec, rm, cc);
    let cosh_b = z.1.cosh(prec, rm, cc);
    let sinh_b = z.1.sinh(prec, rm, cc);
    (
        cos_a.mul(&cosh_b, prec, rm),
        sin_a.mul(&sinh_b, prec, rm).neg(),
    )
}

fn c_sinh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // sinh(a+bi) = sinh(a)cos(b) + i·cosh(a)sin(b)
    let sinh_a = z.0.sinh(prec, rm, cc);
    let cosh_a = z.0.cosh(prec, rm, cc);
    let cos_b = z.1.cos(prec, rm, cc);
    let sin_b = z.1.sin(prec, rm, cc);
    (sinh_a.mul(&cos_b, prec, rm), cosh_a.mul(&sin_b, prec, rm))
}

fn c_cosh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // cosh(a+bi) = cosh(a)cos(b) + i·sinh(a)sin(b)
    let cosh_a = z.0.cosh(prec, rm, cc);
    let sinh_a = z.0.sinh(prec, rm, cc);
    let cos_b = z.1.cos(prec, rm, cc);
    let sin_b = z.1.sin(prec, rm, cc);
    (cosh_a.mul(&cos_b, prec, rm), sinh_a.mul(&sin_b, prec, rm))
}

fn c_pow(base: &Complex, exp: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // base^exp = exp(exp * ln(base))
    // Special case: base is zero.
    if base.0.is_zero() && base.1.is_zero() {
        // 0^0 = 1 by convention; 0^(positive) = 0.
        if exp.0.is_zero() && exp.1.is_zero() {
            return c_one(prec);
        }
        return c_zero(prec);
    }
    let ln_base = c_ln(base, prec, rm, cc);
    let product = c_mul(exp, &ln_base, prec, rm);
    c_exp(&product, prec, rm, cc)
}

fn c_powi(base: &Complex, n: usize, prec: usize, rm: RoundingMode) -> Complex {
    if n == 0 {
        return c_one(prec);
    }
    let mut result = c_one(prec);
    let mut b = base.clone();
    let mut exp = n;
    while exp > 0 {
        if exp & 1 == 1 {
            result = c_mul(&result, &b, prec, rm);
        }
        b = c_mul(&b, &b, prec, rm);
        exp >>= 1;
    }
    result
}

fn c_sqrt(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    if z.0.is_zero() && z.1.is_zero() {
        return c_zero(prec);
    }
    let half = c_from_real(
        BigFloat::from_i32(1, prec).div(&BigFloat::from_i32(2, prec), prec, rm),
        prec,
    );
    c_pow(z, &half, prec, rm, cc)
}

// ── Inverse trig (complex) ─────────────────────────────────────────

fn c_asin(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // asin(z) = -i * ln(iz + sqrt(1 - z²))
    let i_unit = c_i(prec);
    let one = c_one(prec);
    let z_sq = c_mul(z, z, prec, rm);
    let one_minus_z_sq = c_sub(&one, &z_sq, prec, rm);
    let sqrt_term = c_sqrt(&one_minus_z_sq, prec, rm, cc);
    let iz = c_mul(&i_unit, z, prec, rm);
    let sum = c_add(&iz, &sqrt_term, prec, rm);
    let ln_sum = c_ln(&sum, prec, rm, cc);
    let neg_i = c_neg(&i_unit, prec, rm);
    c_mul(&neg_i, &ln_sum, prec, rm)
}

fn c_acos(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // acos(z) = pi/2 - asin(z)
    let half_pi_val = cc.pi(prec, rm).div(&BigFloat::from_i32(2, prec), prec, rm);
    let half_pi = c_from_real(half_pi_val, prec);
    let asin_z = c_asin(z, prec, rm, cc);
    c_sub(&half_pi, &asin_z, prec, rm)
}

fn c_atan(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // atan(z) = (i/2) * ln((i+z)/(i-z))
    let i_unit = c_i(prec);
    let two = c_from_real(BigFloat::from_i32(2, prec), prec);
    let i_over_2 = c_div(&i_unit, &two, prec, rm);
    let i_plus_z = c_add(&i_unit, z, prec, rm);
    let i_minus_z = c_sub(&i_unit, z, prec, rm);
    let ratio = c_div(&i_plus_z, &i_minus_z, prec, rm);
    let ln_ratio = c_ln(&ratio, prec, rm, cc);
    c_mul(&i_over_2, &ln_ratio, prec, rm)
}

// ── Inverse hyperbolic (complex) ──────────────────────────────────

fn c_asinh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // asinh(z) = ln(z + sqrt(z² + 1))
    let one = c_one(prec);
    let z_sq = c_mul(z, z, prec, rm);
    let z_sq_plus_one = c_add(&z_sq, &one, prec, rm);
    let sqrt_term = c_sqrt(&z_sq_plus_one, prec, rm, cc);
    let sum = c_add(z, &sqrt_term, prec, rm);
    c_ln(&sum, prec, rm, cc)
}

fn c_acosh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // acosh(z) = ln(z + sqrt(z² - 1))
    let one = c_one(prec);
    let z_sq = c_mul(z, z, prec, rm);
    let z_sq_minus_one = c_sub(&z_sq, &one, prec, rm);
    let sqrt_term = c_sqrt(&z_sq_minus_one, prec, rm, cc);
    let sum = c_add(z, &sqrt_term, prec, rm);
    c_ln(&sum, prec, rm, cc)
}

fn c_atanh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
    // atanh(z) = (1/2) * ln((1+z)/(1-z))
    let one = c_one(prec);
    let two = c_from_real(BigFloat::from_i32(2, prec), prec);
    let half = c_div(&one, &two, prec, rm);
    let one_plus_z = c_add(&one, z, prec, rm);
    let one_minus_z = c_sub(&one, z, prec, rm);
    let ratio = c_div(&one_plus_z, &one_minus_z, prec, rm);
    let ln_ratio = c_ln(&ratio, prec, rm, cc);
    c_mul(&half, &ln_ratio, prec, rm)
}

// ═══════════════════════════════════════════════════════════════════════════
// Real-valued helpers
// ═══════════════════════════════════════════════════════════════════════════

/// Compute atan2(y, x) using BigFloat arithmetic.
fn atan2_bf(
    y: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> BigFloat {
    if x.is_zero() && y.is_zero() {
        return BigFloat::new(prec);
    }

    let pi_val = cc.pi(prec, rm);
    let two = BigFloat::from_i32(2, prec);
    let half_pi = pi_val.div(&two, prec, rm);

    if x.is_zero() {
        return if y.is_positive() {
            half_pi
        } else {
            half_pi.neg()
        };
    }

    let ratio = y.div(x, prec, rm);
    let atan_val = ratio.atan(prec, rm, cc);

    if x.is_positive() {
        atan_val
    } else {
        let pi_val2 = cc.pi(prec, rm);
        if y.is_negative() {
            atan_val.sub(&pi_val2, prec, rm)
        } else {
            atan_val.add(&pi_val2, prec, rm)
        }
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Special function helpers (f64-based)
// ═══════════════════════════════════════════════════════════════════════════

/// Convert a `BigFloat` to `f64` via decimal radix conversion.
fn bigfloat_to_f64(bf: &BigFloat, rm: RoundingMode, cc: &mut Consts) -> Result<f64, SymplexError> {
    if bf.is_zero() {
        return Ok(0.0);
    }
    if bf.is_inf_pos() {
        return Ok(f64::INFINITY);
    }
    if bf.is_inf_neg() {
        return Ok(f64::NEG_INFINITY);
    }
    if bf.is_nan() {
        return Err(SymplexError::Unevaluable {
            reason: "NaN in BigFloat to f64 conversion".into(),
        });
    }

    let (sign, mantissa, exponent) =
        bf.convert_to_radix(Radix::Dec, rm, cc)
            .map_err(|e| SymplexError::Unevaluable {
                reason: format!("BigFloat to f64 conversion failed: {e:?}"),
            })?;

    // mantissa is [d1, d2, ...] representing 0.d1d2d3... × 10^exponent
    let mut s = String::with_capacity(mantissa.len() + 8);
    if sign == Sign::Neg {
        s.push('-');
    }
    s.push_str("0.");
    for &d in &mantissa {
        s.push((b'0' + d) as char);
    }
    s.push('e');
    s.push_str(&exponent.to_string());

    s.parse::<f64>().map_err(|_| SymplexError::Unevaluable {
        reason: "failed to parse BigFloat decimal representation as f64".into(),
    })
}

/// Lanczos approximation for the Gamma function (g=7, 9 coefficients).
#[allow(dead_code)]
#[allow(clippy::excessive_precision)]
fn lanczos_gamma_f64(x: f64) -> Result<f64, SymplexError> {
    if x.is_nan() || x.is_infinite() {
        return Err(SymplexError::Unevaluable {
            reason: "Gamma of special float value".into(),
        });
    }

    // Reflection formula for x < 0.5:  Gamma(x) = pi / (sin(pi*x) * Gamma(1-x))
    if x < 0.5 {
        let reflected = lanczos_gamma_f64(1.0 - x)?;
        let sin_pi_x = (std::f64::consts::PI * x).sin();
        if sin_pi_x.abs() < 1e-300 {
            return Err(SymplexError::Unevaluable {
                reason: "Gamma at non-positive integer pole".into(),
            });
        }
        let result = std::f64::consts::PI / (sin_pi_x * reflected);
        if result.is_finite() {
            return Ok(result);
        }
        return Err(SymplexError::Unevaluable {
            reason: "Gamma at pole".into(),
        });
    }

    // Lanczos approximation with g=7
    const COEFFICIENTS: [f64; 9] = [
        0.99999999999980993,
        676.5203681218851,
        -1259.1392167224028,
        771.32342877765313,
        -176.61502916214059,
        12.507343278686905,
        -0.13857109526572012,
        9.9843695780195716e-6,
        1.5056327351493116e-7,
    ];

    let x = x - 1.0; // Lanczos uses Gamma(x+1) = x!
    let t = x + 7.0 + 0.5; // g = 7

    let mut sum = COEFFICIENTS[0];
    for (i, &c) in COEFFICIENTS[1..].iter().enumerate() {
        sum += c / (x + i as f64 + 1.0);
    }

    let result = (2.0 * std::f64::consts::PI).sqrt() * t.powf(x + 0.5) * (-t).exp() * sum;
    Ok(result)
}

/// Error function via Taylor series (small |x|) or asymptotic expansion (large |x|).
#[allow(dead_code)]
fn erf_f64(x: f64) -> f64 {
    if x.abs() < 4.0 {
        // Taylor series: erf(x) = 2/sqrt(pi) * sum_{n=0}^{inf} (-1)^n * x^(2n+1) / (n! * (2n+1))
        let mut sum = 0.0;
        let mut term = x; // first term: x
        sum += term;
        for n in 1..50 {
            term *= -x * x / n as f64;
            sum += term / (2 * n + 1) as f64;
        }
        sum * 2.0 / std::f64::consts::PI.sqrt()
    } else {
        // For large |x|, use complementary: erf(x) = 1 - erfc(x)
        // erfc(x) ~ exp(-x^2)/(x*sqrt(pi)) * sum_{n=0} (-1)^n * (2n-1)!! / (2x^2)^n
        let sign = x.signum();
        let ax = x.abs();
        let mut sum = 1.0;
        let mut term = 1.0;
        for n in 1..20 {
            term *= -(2 * n - 1) as f64 / (2.0 * ax * ax);
            if term.abs() < 1e-16 {
                break;
            }
            sum += term;
        }
        let erfc = (-ax * ax).exp() / (ax * std::f64::consts::PI.sqrt()) * sum;
        sign * (1.0 - erfc)
    }
}

/// Arbitrary-precision digamma (psi) function via recurrence + asymptotic series.
///
/// Algorithm:
/// 1. For x < 0, use reflection: ψ(x) = ψ(1−x) − π·cot(πx)
/// 2. Use recurrence ψ(x+1) = ψ(x) + 1/x to shift x to a large value
/// 3. Asymptotic expansion: ψ(x) ~ ln(x) − 1/(2x) − Σ B_{2k}/(2k · x^{2k})
fn arb_digamma(
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let guard = 20;
    let wp = prec + guard;

    if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() {
        return Err(SymplexError::Unevaluable {
            reason: "Digamma of special float value".into(),
        });
    }

    let one = BigFloat::from_i32(1, wp);

    // Handle negative x via reflection: ψ(x) = ψ(1−x) − π·cot(πx)
    if x.is_negative() {
        let pi_val = cc.pi(wp, rm).clone();
        let pi_x = pi_val.mul(x, wp, rm);
        let sin_val = pi_x.sin(wp, rm, cc);

        // Check for pole at non-positive integers
        if sin_val.is_zero() {
            return Err(SymplexError::Unevaluable {
                reason: "Digamma at non-positive integer pole".into(),
            });
        }
        if let Some(s_exp) = sin_val.exponent()
            && (s_exp as i64) < -(wp as i64 / 2)
        {
            return Err(SymplexError::Unevaluable {
                reason: "Digamma at non-positive integer pole".into(),
            });
        }

        let cos_val = pi_x.cos(wp, rm, cc);
        let cot_val = cos_val.div(&sin_val, wp, rm);
        let one_minus_x = one.sub(x, wp, rm);
        let psi_1mx = arb_digamma(&one_minus_x, prec, rm, cc)?;
        let pi_cot = cc.pi(wp, rm).clone().mul(&cot_val, wp, rm);
        return Ok(psi_1mx.sub(&pi_cot, wp, rm));
    }

    // Use recurrence ψ(x+1) = ψ(x) + 1/x to shift x >= threshold
    let threshold_val = (wp as i32 / 3).max(10);
    let threshold = BigFloat::from_i32(threshold_val, wp);
    let mut result = BigFloat::new(wp); // 0
    let mut x = x.clone();

    while x.sub(&threshold, wp, rm).is_negative() {
        // Check for pole at zero
        if x.is_zero() {
            return Err(SymplexError::Unevaluable {
                reason: "Digamma at non-positive integer pole".into(),
            });
        }
        if let Some(x_exp) = x.exponent()
            && (x_exp as i64) < -(wp as i64 / 2)
        {
            return Err(SymplexError::Unevaluable {
                reason: "Digamma at non-positive integer pole".into(),
            });
        }
        let inv_x = one.div(&x, wp, rm);
        result = result.sub(&inv_x, wp, rm);
        x = x.add(&one, wp, rm);
    }

    // Asymptotic expansion: ψ(x) ~ ln(x) − 1/(2x) − Σ_{k=1}^{N} B_{2k}/(2k · x^{2k})
    let ln_x = x.ln(wp, rm, cc);
    result = result.add(&ln_x, wp, rm);

    let two = BigFloat::from_i32(2, wp);
    let half_inv_x = one.div(&x.mul(&two, wp, rm), wp, rm);
    result = result.sub(&half_inv_x, wp, rm);

    let x2 = x.mul(&x, wp, rm);
    let mut x_pow = x2.clone(); // x^2

    // Bernoulli numbers B_{2k} for k=1..12 as (numerator, denominator):
    // B2=1/6, B4=−1/30, B6=1/42, B8=−1/30, B10=5/66, B12=−691/2730
    // B14=7/6, B16=−3617/510, B18=43867/798, B20=−174611/330
    // B22=854513/138, B24=−236364091/2730
    let bernoulli_nums: &[(i128, i128)] = &[
        (1, 6),
        (-1, 30),
        (1, 42),
        (-1, 30),
        (5, 66),
        (-691, 2730),
        (7, 6),
        (-3617, 510),
        (43867, 798),
        (-174611, 330),
        (854513, 138),
        (-236364091, 2730),
    ];

    let n_terms = (prec / 6 + 2).min(bernoulli_nums.len());

    for (k_idx, &(bn, bd)) in bernoulli_nums[..n_terms].iter().enumerate() {
        let k = (k_idx + 1) as i128;
        let two_k = 2 * k;
        // Term = B_{2k} / (2k · x^{2k})
        let coeff_n = BigFloat::from_i128(bn, wp);
        let coeff_d = BigFloat::from_i128(bd * two_k, wp);
        let coeff = coeff_n.div(&coeff_d, wp, rm);
        let inv_xpow = one.div(&x_pow, wp, rm);
        let term = coeff.mul(&inv_xpow, wp, rm);
        result = result.sub(&term, wp, rm);

        // Convergence check via binary exponents
        if let (Some(t_exp), Some(r_exp)) = (term.exponent(), result.exponent())
            && (r_exp as i64 - t_exp as i64) > wp as i64
        {
            break;
        }

        x_pow = x_pow.mul(&x2, wp, rm);
    }

    Ok(result)
}

// ═══════════════════════════════════════════════════════════════════════════
// Arbitrary-precision Gamma via Stirling series
// ═══════════════════════════════════════════════════════════════════════════

/// Compute log Γ(z) for real z > 0 using the Stirling asymptotic series.
///
/// Algorithm:
/// 1. Argument reduction — shift z by integer r until z + r ≥ 0.2 · p.
/// 2. Stirling series:
///    log Γ(z) = (z − ½) ln(z) − z + ½ ln(2π) + Σ B₂ₖ / [2k(2k−1) z^{2k−1}].
/// 3. Undo shift:
///    log Γ(z) = log Γ(z + r) − Σ_{i=0}^{r−1} ln(z + i).
fn stirling_log_gamma(
    z: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    tracing::debug!(prec, "stirling_log_gamma: entry");

    let guard = 20; // extra guard bits
    let wp = prec + guard;

    // Determine how far to shift z for convergence.
    let z_approx = bigfloat_to_f64(z, rm, cc)?;
    let threshold = ((wp as f64) * 0.2).ceil() as i64;
    let shift = if z_approx < threshold as f64 {
        (threshold - z_approx.floor() as i64).max(0) as usize
    } else {
        0
    };

    tracing::trace!(
        z_approx,
        threshold,
        shift,
        "stirling_log_gamma: argument reduction"
    );

    // z_shifted = z + shift
    let shift_bf = BigFloat::from_i128(shift as i128, wp);
    let z_shifted = z.add(&shift_bf, wp, rm);

    // Main Stirling formula:
    //   log Γ(z) = (z − 1/2)·ln(z) − z + (1/2)·ln(2π) + Σ ...
    let half = BigFloat::from_f64(0.5, wp);
    let z_minus_half = z_shifted.sub(&half, wp, rm);
    let log_z = z_shifted.ln(wp, rm, cc);

    let mut result = z_minus_half.mul(&log_z, wp, rm);
    result = result.sub(&z_shifted, wp, rm);

    // (1/2) · ln(2π) at working precision.
    let two = BigFloat::from_i32(2, wp);
    let pi_val = cc.pi(wp, rm).clone();
    let two_pi = two.mul(&pi_val, wp, rm);
    let half_log_2pi = two_pi.ln(wp, rm, cc).mul(&half, wp, rm);
    result = result.add(&half_log_2pi, wp, rm);

    // Series: Σ_{k=1}^{N} B_{2k} / [2k(2k−1) · z^{2k−1}]
    let z_sq = z_shifted.mul(&z_shifted, wp, rm);
    let mut z_pow = z_shifted.clone(); // z^1 → z^3 → z^5 → …

    // Upper bound on useful terms: safely below the divergence point πz.
    let max_terms = ((wp as f64) * 0.35) as usize + 10;
    let mut prev_term_exp: Option<i32> = None;

    for k in 1..=max_terms {
        let b2k = crate::base::bernoulli::bernoulli(2 * k);
        if b2k.is_zero() {
            continue;
        }

        // Advance z_pow: z^1, z^3, z^5, …
        if k > 1 {
            z_pow = z_pow.mul(&z_sq, wp, rm);
        }

        // term = B_{2k} / (2k · (2k−1) · z^{2k−1})
        let b2k_bf = ratio_to_bigfloat(&b2k, wp, rm);
        let denom_int = (2 * k * (2 * k - 1)) as i128;
        let denom = BigFloat::from_i128(denom_int, wp);
        let denom_full = denom.mul(&z_pow, wp, rm);
        let term = b2k_bf.div(&denom_full, wp, rm);

        // Divergence / convergence detection via binary exponents.
        if let Some(t_exp) = term.exponent() {
            if let Some(prev) = prev_term_exp
                && t_exp > prev + 10
            {
                tracing::trace!(k, "stirling series: diverging, stopping");
                break;
            }
            prev_term_exp = Some(t_exp);

            // Term is negligible relative to accumulated result.
            if let Some(r_exp) = result.exponent()
                && (r_exp as i64 - t_exp as i64) > wp as i64
            {
                tracing::trace!(k, "stirling series: converged");
                break;
            }
        }

        result = result.add(&term, wp, rm);
    }

    // Undo argument reduction:
    //   log Γ(z) = log Γ(z+shift) − Σ_{i=0}^{shift−1} ln(z + i)
    for i in 0..shift {
        let z_plus_i = z.add(&BigFloat::from_i128(i as i128, wp), wp, rm);
        let log_zi = z_plus_i.ln(wp, rm, cc);
        result = result.sub(&log_zi, wp, rm);
    }

    Ok(result)
}

/// Compute Γ(x) at arbitrary precision for real x.
///
/// For positive x, evaluates exp(stirling_log_gamma(x)).
/// For negative non-integer x, applies the reflection formula
/// Γ(z) = π / (sin(πz) · Γ(1 − z)).
fn arb_gamma_real(
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    // Pole at zero.
    if x.is_zero() {
        return Err(SymplexError::Unevaluable {
            reason: "Gamma at non-positive integer pole".into(),
        });
    }

    let x_approx = bigfloat_to_f64(x, rm, cc)?;

    // Detect non-positive integer poles via f64 approximation.
    if x_approx <= 0.0 {
        let rounded = x_approx.round();
        if (rounded - x_approx).abs() < 1e-12 && rounded <= 0.0 {
            return Err(SymplexError::Unevaluable {
                reason: "Gamma at non-positive integer pole".into(),
            });
        }
    }

    if x_approx > 0.0 {
        // Direct Stirling path.
        tracing::debug!(x_approx, "arb_gamma_real: positive argument");
        let log_gamma = stirling_log_gamma(x, prec, rm, cc)?;
        Ok(log_gamma.exp(prec, rm, cc))
    } else {
        // Reflection: Γ(z) = π / (sin(πz) · Γ(1 − z))
        tracing::debug!(x_approx, "arb_gamma_real: reflection formula");
        let one = BigFloat::from_i32(1, prec);
        let one_minus_x = one.sub(x, prec, rm);

        let log_gamma_1mx = stirling_log_gamma(&one_minus_x, prec, rm, cc)?;
        let gamma_1mx = log_gamma_1mx.exp(prec, rm, cc);

        let pi_val = cc.pi(prec, rm).clone();
        let pi_x = pi_val.mul(x, prec, rm);
        let sin_pi_x = pi_x.sin(prec, rm, cc);

        // Guard against the pole (sin(πz) ≈ 0 for integer z).
        if sin_pi_x.is_zero() {
            return Err(SymplexError::Unevaluable {
                reason: "Gamma at non-positive integer pole".into(),
            });
        }
        if let Some(s_exp) = sin_pi_x.exponent()
            && (s_exp as i64) < -(prec as i64 / 2)
        {
            return Err(SymplexError::Unevaluable {
                reason: "Gamma at non-positive integer pole".into(),
            });
        }

        let denom = sin_pi_x.mul(&gamma_1mx, prec, rm);
        let pi_val2 = cc.pi(prec, rm).clone();
        Ok(pi_val2.div(&denom, prec, rm))
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Arbitrary-precision erf via Taylor series / asymptotic expansion
// ═══════════════════════════════════════════════════════════════════════════

/// Compute erf(x) at arbitrary precision for real x.
///
/// For small-to-moderate |x|, uses the Taylor series:
///   erf(x) = (2/√π) Σ_{n=0}^{N} (-1)^n x^{2n+1} / (n! (2n+1))
///
/// For large |x|, uses the asymptotic expansion of erfc:
///   erfc(x) = exp(-x²) / (x√π) · Σ_{n=0}^{N} (-1)^n (2n-1)!! / (2x²)^n
///   erf(x) = sign(x) · (1 − erfc(|x|))
fn arb_erf(
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_zero() {
        return Ok(BigFloat::new(prec));
    }

    let guard = 32;
    let wp = prec + guard;

    let x_approx = bigfloat_to_f64(x, rm, cc)?;
    // Threshold: for |x| beyond this, the asymptotic expansion converges
    // faster than the Taylor series. Roughly √(wp * ln(2) / 2).
    let threshold = ((wp as f64) * 0.35).sqrt() + 2.0;

    if x_approx.abs() < threshold {
        // ── Taylor series ──────────────────────────────────────────
        // erf(x) = (2/√π) · Σ_{n=0}^{N} prod_n / (2n+1)
        // where prod_0 = x, prod_{n+1} = prod_n · (-x²) / (n+1)
        let neg_x_sq = x.mul(x, wp, rm).neg();
        let mut prod = x.clone(); // (-x²)^n · x / n!
        let mut sum = x.clone(); // accumulator (first term = x)

        let max_terms = (wp as f64 * 0.6) as usize + 60;
        for n in 1..=max_terms {
            // prod *= -x² / n
            prod = prod.mul(&neg_x_sq, wp, rm);
            prod = prod.div(&BigFloat::from_i32(n as i32, wp), wp, rm);

            // term = prod / (2n+1)
            let divisor = BigFloat::from_i32((2 * n + 1) as i32, wp);
            let term = prod.div(&divisor, wp, rm);

            // Convergence check
            if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
                && (s_exp as i64 - t_exp as i64) > wp as i64
            {
                break;
            }

            sum = sum.add(&term, wp, rm);
        }

        // Multiply by 2/√π
        let two = BigFloat::from_i32(2, wp);
        let pi_val = cc.pi(wp, rm).clone();
        let sqrt_pi = pi_val.sqrt(wp, rm);
        let two_over_sqrt_pi = two.div(&sqrt_pi, wp, rm);

        Ok(sum.mul(&two_over_sqrt_pi, wp, rm))
    } else {
        // ── Asymptotic expansion for large |x| ────────────────────
        // erfc(x) = exp(-x²)/(x√π) · Σ_{n=0}^{N} (-1)^n (2n-1)!! / (2x²)^n
        let is_neg = x.is_negative();
        let ax = x.abs();
        let x_sq = ax.mul(&ax, wp, rm);
        let two_x_sq = x_sq.mul(&BigFloat::from_i32(2, wp), wp, rm);

        let mut sum = BigFloat::from_i32(1, wp);
        let mut term = BigFloat::from_i32(1, wp);
        let max_terms = wp / 2 + 30;

        for n in 1..=max_terms {
            // term *= -(2n-1) / (2x²)
            let factor = BigFloat::from_i32(2 * n as i32 - 1, wp);
            term = term.mul(&factor, wp, rm);
            term = term.div(&two_x_sq, wp, rm);
            term = term.neg();

            // Divergence check: if |term| starts growing, stop
            if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
                && t_exp > s_exp
            {
                break;
            }

            sum = sum.add(&term, wp, rm);
        }

        let exp_neg_x_sq = x_sq.neg().exp(wp, rm, cc);
        let sqrt_pi = cc.pi(wp, rm).clone().sqrt(wp, rm);
        let x_sqrt_pi = ax.mul(&sqrt_pi, wp, rm);
        let erfc_val = exp_neg_x_sq.mul(&sum, wp, rm).div(&x_sqrt_pi, wp, rm);

        let one = BigFloat::from_i32(1, wp);
        let erf_val = one.sub(&erfc_val, wp, rm);

        if is_neg {
            Ok(erf_val.neg())
        } else {
            Ok(erf_val)
        }
    }
}

/// Arbitrary-precision Lambert W function (principal branch) via Halley iteration.
///
/// Solves w·exp(w) = x for w, using cubic-convergent Halley steps.
fn arb_lambert_w(
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    // Work with extra guard bits for intermediate rounding.
    let wp = prec + 32;

    // Handle x = 0 exactly.
    if x.is_zero() {
        return Ok(BigFloat::new(prec));
    }

    // Handle negative x near -1/e boundary.
    // The principal branch is defined for x >= -1/e.
    let neg_inv_e = {
        let e_val = cc.e(wp, rm).clone();
        let one = BigFloat::from_i32(1, wp);
        one.div(&e_val, wp, rm).neg()
    };

    let diff = x.sub(&neg_inv_e, wp, rm);
    if diff.is_negative() {
        return Err(SymplexError::Unevaluable {
            reason: "LambertW: argument < -1/e, outside principal branch domain".into(),
        });
    }

    // Initial guess.
    let mut w = if x.is_negative() {
        // Near -1/e: use series w ≈ -1 + sqrt(2(1 + ex))
        let e_val = cc.e(wp, rm).clone();
        let ex = e_val.mul(x, wp, rm);
        let one = BigFloat::from_i32(1, wp);
        let two = BigFloat::from_i32(2, wp);
        let inner = one.add(&ex, wp, rm).mul(&two, wp, rm);
        if inner.is_negative() || inner.is_zero() {
            BigFloat::from_i32(-1, wp)
        } else {
            let sq = inner.sqrt(wp, rm);
            BigFloat::from_i32(-1, wp).add(&sq, wp, rm)
        }
    } else {
        // For small positive x, w ≈ x is a good start.
        // For large x, w ≈ ln(x) - ln(ln(x)).
        let threshold = BigFloat::from_f64(2.5, wp);
        if x.sub(&threshold, wp, rm).is_negative() {
            x.clone()
        } else {
            let ln_x = x.ln(wp, rm, cc);
            let ln_ln_x = ln_x.ln(wp, rm, cc);
            ln_x.sub(&ln_ln_x, wp, rm)
        }
    };

    // Halley iteration: cubically convergent.
    //
    // Given f(w) = w·e^w − x, f'(w) = (w+1)·e^w, f''(w) = (w+2)·e^w,
    // the Halley step is:
    //   δ = f / (f' − f·f'' / (2·f'))
    //     = (w·e^w − x) / ((w+1)·e^w − (w+2)·(w·e^w − x) / (2·(w+1)))
    let max_iter = 100;
    for _ in 0..max_iter {
        let ew = w.exp(wp, rm, cc);
        let wew = w.mul(&ew, wp, rm);
        let residual = wew.sub(x, wp, rm); // w·e^w − x

        let w_plus_1 = w.add(&BigFloat::from_i32(1, wp), wp, rm);
        let denom_base = w_plus_1.mul(&ew, wp, rm); // (w+1)·e^w

        // Halley denominator: (w+1)·e^w − (w+2)·residual / (2·(w+1))
        let w_plus_2 = w.add(&BigFloat::from_i32(2, wp), wp, rm);
        let two_w_plus_1 = w_plus_1.mul(&BigFloat::from_i32(2, wp), wp, rm);

        let correction_numer = w_plus_2.mul(&residual, wp, rm);
        let correction = if two_w_plus_1.is_zero() {
            BigFloat::new(wp)
        } else {
            correction_numer.div(&two_w_plus_1, wp, rm)
        };
        let denom = denom_base.sub(&correction, wp, rm);

        if denom.is_zero() {
            break;
        }

        let delta = residual.div(&denom, wp, rm);
        w = w.sub(&delta, wp, rm);

        // Check convergence: |delta| exponent is far below working precision.
        if delta.is_zero() {
            break;
        }
        if let (Some(d_exp), Some(w_e)) = (delta.exponent(), w.exponent())
            && (d_exp as i64) < (w_e as i64) - (wp as i64)
        {
            break;
        }
    }

    Ok(w)
}

// ═══════════════════════════════════════════════════════════════════════════
// Helpers
// ═══════════════════════════════════════════════════════════════════════════

// ═══════════════════════════════════════════════════════════════════════════
// Euler-Mascheroni constant at arbitrary precision (Brent-McMillan B1)
// ═══════════════════════════════════════════════════════════════════════════

/// Compute the Euler-Mascheroni constant γ at `prec` bits of precision
/// using the Brent-McMillan B1 algorithm.
///
/// This is the fastest known algorithm for computing γ.  It uses the identity
/// involving modified Bessel functions, simplified to a fixed-point summation
/// with convergence rate O(e^{-4n}) where n = 2^p.
///
/// Reference: Brent & McMillan (1980), mpmath `euler_fixed`.
fn arb_euler_gamma(
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    // We need ln(2) at working precision.
    let extra = 30;
    let wp = prec + extra;

    // Choose p such that e^{-4·2^p} < 2^{-wp}, i.e. 4·2^p > wp·ln(2),
    // i.e. p > log2(wp·ln(2)/4).
    let p = ((wp as f64 / 4.0) * std::f64::consts::LN_2).log2().ceil() as u64 + 1;
    let n: i128 = 1i128 << p;
    let n_sq = n * n;

    // ln(2) at working precision.
    let two_bf = BigFloat::from_i32(2, wp);
    let ln2 = two_bf.ln(wp, rm, cc);

    // A = U = -p · ln(2)   (= -ln(n) since n = 2^p)
    let p_bf = BigFloat::from_i128(p as i128, wp);
    let neg_p_ln2 = p_bf.mul(&ln2, wp, rm).neg();

    let mut a = neg_p_ln2.clone();
    let mut u = neg_p_ln2;
    // B = V = 1
    let one = BigFloat::from_i32(1, wp);
    let mut b = one.clone();
    let mut v = one.clone();

    let n_sq_bf = BigFloat::from_i128(n_sq, wp);

    let mut k: i128 = 1;
    loop {
        let k_bf = BigFloat::from_i128(k, wp);
        let k_sq_bf = k_bf.mul(&k_bf, wp, rm);

        // B = B · n² / k²
        b = b.mul(&n_sq_bf, wp, rm).div(&k_sq_bf, wp, rm);

        // A = (A · n² / k + B) / k
        a = a.mul(&n_sq_bf, wp, rm).div(&k_bf, wp, rm);
        a = a.add(&b, wp, rm).div(&k_bf, wp, rm);

        u = u.add(&a, wp, rm);
        v = v.add(&b, wp, rm);

        // Convergence: both A and B are negligibly small.
        let a_small = a.exponent().is_none_or(|e| e < -(wp as i32) + 10);
        let b_small = b.exponent().is_none_or(|e| e < -(wp as i32) + 10);
        if a_small && b_small {
            tracing::trace!(k, "arb_euler_gamma: converged");
            break;
        }

        k += 1;
        if k > 10 * (wp as i128) {
            return Err(SymplexError::ComputationFailed {
                operation: "euler_gamma",
                reason: "Brent-McMillan did not converge".into(),
            });
        }
    }

    // γ = U / V
    Ok(u.div(&v, prec, rm))
}

// ═══════════════════════════════════════════════════════════════════════════
// Hankel P/Q series for Bessel asymptotic expansion
// ═══════════════════════════════════════════════════════════════════════════

/// Compute the Hankel auxiliary functions P_ν(z) and Q_ν(z) for the
/// Bessel asymptotic expansion (DLMF §10.17).
///
/// P_ν(z) = Σ_{k=0}^{N} (-1)^k · a_{2k}(ν) / z^{2k}
/// Q_ν(z) = Σ_{k=0}^{N} (-1)^k · a_{2k+1}(ν) / z^{2k+1}
///
/// where a_k(ν) = (4ν²−1²)(4ν²−3²)⋯(4ν²−(2k−1)²) / (k! · 8^k)
/// (DLMF 10.17.1), so a_k = a_{k−1} · (4ν² − (2k−1)²) / (8k).
///
/// Uses optimal truncation: stops when terms start increasing (divergent
/// series).  The error is bounded by the first omitted term for real
/// positive z (DLMF §10.17(iii)).  The returned flag is `true` when the
/// truncation happened because the terms fell below `2^{-prec}` relative
/// to the sums — i.e. the expansion actually delivered `prec` bits.  When
/// it is `false` the series started diverging before reaching the target
/// and the caller must fall back to the power series.
fn hankel_pq(
    order: &BigFloat,
    z: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    _cc: &mut Consts,
) -> (BigFloat, BigFloat, bool) {
    let one = BigFloat::from_i32(1, prec);
    let four = BigFloat::from_i32(4, prec);
    let eight = BigFloat::from_i32(8, prec);

    let four_nu_sq = four.mul(&order.mul(order, prec, rm), prec, rm);

    let z_inv = one.div(z, prec, rm);

    // P and Q accumulators; the k = 0 term (a_0 = 1) is added by the loop.
    let mut p_sum = BigFloat::new(prec);
    let mut q_sum = BigFloat::new(prec);

    let mut a_k = one.clone(); // a_0 = 1
    let mut z_power = one.clone(); // z^0 = 1

    // The terms reach their minimum near k ≈ 2|z|; allow a little slack.
    let max_terms = (3 * prec + 100).min(100_000);

    let mut prev_a_abs: Option<BigFloat> = None;
    let mut converged = false;

    for k in 0..max_terms {
        if k > 0 {
            // a_k = a_{k-1} · (4ν² - (2k-1)²) / (8k)
            let two_km1 = BigFloat::from_i32((2 * k as i32) - 1, prec);
            let two_km1_sq = two_km1.mul(&two_km1, prec, rm);
            let numer = four_nu_sq.sub(&two_km1_sq, prec, rm);
            let k_bf = BigFloat::from_i32(k as i32, prec);
            let denom = eight.mul(&k_bf, prec, rm);
            a_k = a_k.mul(&numer, prec, rm).div(&denom, prec, rm);

            // Update z_power: multiply by 1/z each step.
            z_power = z_power.mul(&z_inv, prec, rm);
        }

        // term = a_k / z^k
        let term = a_k.mul(&z_power, prec, rm);

        // Check for divergence: if |term| > |prev_term|, stop.
        let term_abs = term.abs();
        if let Some(ref prev) = prev_a_abs
            && term_abs.partial_cmp(prev) == Some(std::cmp::Ordering::Greater)
            && k > 2
        {
            tracing::trace!(k, "hankel_pq: optimal truncation (terms diverging)");
            break;
        }

        // Convergence: term negligible relative to accumulated sums.
        let p_exp = p_sum.exponent().unwrap_or(0);
        let q_exp = q_sum.exponent().unwrap_or(0);
        let ref_exp = p_exp.max(q_exp);
        let t_exp = term.exponent().unwrap_or(i32::MIN / 2);
        if k > 0 && (ref_exp as i64 - t_exp as i64) > prec as i64 {
            tracing::trace!(k, "hankel_pq: converged (term negligible)");
            converged = true;
            break;
        }

        prev_a_abs = Some(term_abs);

        // Even k → contributes to P with sign (-1)^(k/2)
        // Odd k → contributes to Q with sign (-1)^((k-1)/2)
        if k % 2 == 0 {
            // P term: (-1)^(k/2) · a_k / z^k
            if (k / 2) % 2 == 0 {
                p_sum = p_sum.add(&term, prec, rm);
            } else {
                p_sum = p_sum.sub(&term, prec, rm);
            }
        } else {
            // Q term: (-1)^((k-1)/2) · a_k / z^k
            if ((k - 1) / 2) % 2 == 0 {
                q_sum = q_sum.add(&term, prec, rm);
            } else {
                q_sum = q_sum.sub(&term, prec, rm);
            }
        }
    }

    (p_sum, q_sum, converged)
}

/// Smallest `|x|` at which the Hankel expansion is expected to reach `wp`
/// bits for order `ν`.
///
/// The smallest term of the (divergent) expansion is `~e^{−2|x|}` times a
/// factor that grows like `e^{2|ν|}` (the first `|ν|` coefficients are
/// inflated by `4ν²/(2j−1)²`), so `2|x| ≳ (wp + 3|ν|)·ln 2`.  This is
/// only a heuristic for choosing the algorithm; [`hankel_pq`] reports
/// whether it actually converged.
fn hankel_threshold(wp: usize, order_abs: f64) -> f64 {
    (wp as f64 + 3.0 * order_abs) * std::f64::consts::LN_2 / 2.0 + 4.0
}

/// Evaluate `J_ν(x)` (`want_j == true`) or `Y_ν(x)` via the Hankel
/// asymptotic expansion at `x > 0`,
///
/// ```text
/// J_ν(x) = √(2/(πx)) · [cos(ω)·P_ν(x) − sin(ω)·Q_ν(x)]
/// Y_ν(x) = √(2/(πx)) · [sin(ω)·P_ν(x) + cos(ω)·Q_ν(x)],   ω = x − νπ/2 − π/4
/// ```
///
/// Returns `None` when the expansion does not reach `wp` bits at this
/// `x` (the caller then uses the power series).
fn bessel_jy_hankel(
    order: &BigFloat,
    x: &BigFloat,
    want_j: bool,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Option<BigFloat> {
    // The phase ω = x − νπ/2 − π/4 loses the leading `log2 x` bits of x to
    // cancellation, so the trigonometric part needs that many extra bits.
    let wp = wp + x.exponent().unwrap_or(0).max(0) as usize;
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);
    let mut nu = order.clone();
    let _ = nu.set_precision(wp, rm);

    let (p_val, q_val, converged) = hankel_pq(&nu, &xw, wp, rm, cc);
    if !converged {
        tracing::trace!("bessel_jy_hankel: expansion did not reach target precision");
        return None;
    }
    let pi = cc.pi(wp, rm).clone();
    let two = BigFloat::from_i32(2, wp);
    let four = BigFloat::from_i32(4, wp);

    let two_over_pi_x = two.div(&pi.mul(&xw, wp, rm), wp, rm);
    let amplitude = two_over_pi_x.sqrt(wp, rm);

    let nu_pi_half = nu.mul(&pi, wp, rm).div(&two, wp, rm);
    let pi_quarter = pi.div(&four, wp, rm);
    let phase = xw.sub(&nu_pi_half, wp, rm).sub(&pi_quarter, wp, rm);

    let cos_phase = phase.cos(wp, rm, cc);
    let sin_phase = phase.sin(wp, rm, cc);

    let combo = if want_j {
        cos_phase
            .mul(&p_val, wp, rm)
            .sub(&sin_phase.mul(&q_val, wp, rm), wp, rm)
    } else {
        sin_phase
            .mul(&p_val, wp, rm)
            .add(&cos_phase.mul(&q_val, wp, rm), wp, rm)
    };
    Some(amplitude.mul(&combo, wp, rm))
}

// ═══════════════════════════════════════════════════════════════════════════
// Bessel function evaluation (arbitrary precision)
// ═══════════════════════════════════════════════════════════════════════════

/// Arbitrary-precision Bessel function of the first kind J_ν(x).
///
/// * `|x|` below [`hankel_threshold`] (or whenever the asymptotic expansion
///   cannot reach the target precision): ascending power series
///   `J_ν(x) = Σ_k (-1)^k (x/2)^(ν+2k) / (k! Γ(ν+k+1))`.  Its terms grow
///   to `~e^{|x|}` before the alternating sum cancels down to `O(1)`, so
///   the working precision is raised by [`cancellation_guard_bits`] — this
///   makes the series correct at *any* `x`, merely slower for huge `x`.
/// * Otherwise the Hankel expansion with the full `P/Q` series
///   ([`bessel_jy_hankel`]), whose optimal-truncation error is `~e^{−2|x|}`.
fn arb_bessel_j(
    order: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let guard = 32;
    let wp = prec + guard;

    // Special case: x = 0.
    if x.is_zero() {
        let order_f64 = bigfloat_to_f64(order, rm, cc)?;
        return if order_f64.abs() < 1e-15 {
            // J_0(0) = 1
            Ok(BigFloat::from_i32(1, prec))
        } else if order_f64 > 0.0 {
            // J_n(0) = 0 for n > 0
            Ok(BigFloat::new(prec))
        } else {
            Err(SymplexError::Unevaluable {
                reason: "BesselJ at x=0 with negative order".into(),
            })
        };
    }

    let x_f64 = bigfloat_to_f64(x, rm, cc)?;
    let order_f64 = bigfloat_to_f64(order, rm, cc)?;

    // Check for integer order (most common case).
    let order_int = order_f64.round() as i64;
    let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;

    // Negative integer order: J_{-n}(x) = (-1)^n J_n(x).
    if is_int_order && order_int < 0 {
        let pos_order = BigFloat::from_i128((-order_int) as i128, wp);
        let j = arb_bessel_j(&pos_order, x, prec, rm, cc)?;
        return Ok(if order_int % 2 == 0 { j } else { j.neg() });
    }

    // Large |x|: Hankel asymptotic expansion (J is even/odd in x for
    // integer order; non-integer order with x < 0 is complex and rejected).
    if x_f64.abs() >= hankel_threshold(wp, order_f64.abs()) {
        if x.is_negative() {
            if !is_int_order {
                return Err(SymplexError::Unevaluable {
                    reason: "BesselJ of negative argument with non-integer order is complex".into(),
                });
            }
            let j = arb_bessel_j(order, &x.abs(), prec, rm, cc)?;
            return Ok(if order_int % 2 == 0 { j } else { j.neg() });
        }
        if let Some(j) = bessel_jy_hankel(order, x, true, wp, rm, cc) {
            return Ok(round_to(j, prec, rm));
        }
    }

    // ── Ascending series with cancellation guard bits ──────────────────
    // J_ν(x) = (x/2)^ν · Σ_{k=0}^N (-1)^k · (x/2)^{2k} / (k! · Γ(ν+k+1))
    {
        let wp = wp + cancellation_guard_bits(x_f64.abs());
        let mut xw = x.clone();
        let _ = xw.set_precision(wp, rm);
        let mut order_w = order.clone();
        let _ = order_w.set_precision(wp, rm);
        let order = &order_w;
        let two = BigFloat::from_i32(2, wp);
        let x_half = xw.div(&two, wp, rm);
        let x_half_sq = x_half.mul(&x_half, wp, rm);
        let neg_x_half_sq = x_half_sq.neg();

        // Compute (x/2)^ν.  For integer order, use repeated multiplication.
        let prefix = if is_int_order && order_int >= 0 {
            let mut p = BigFloat::from_i32(1, wp);
            for _ in 0..order_int {
                p = p.mul(&x_half, wp, rm);
            }
            p
        } else {
            // General: (x/2)^ν = exp(ν · ln(x/2))
            let ln_xh = x_half.abs().ln(wp, rm, cc);
            let nu_ln = order.mul(&ln_xh, wp, rm);
            nu_ln.exp(wp, rm, cc)
        };

        // Series: sum = Σ (-1)^k · (x/2)^{2k} / (k! · Γ(ν+k+1))
        // We track term = (-x²/4)^k / (k! · Γ(ν+k+1)) incrementally.
        // term_{k+1} = term_k · (-x²/4) / ((k+1) · (ν+k+1))
        let mut sum = BigFloat::new(wp); // will add 1/Γ(ν+1) as first term

        // First term (k=0): 1 / Γ(ν+1)
        let gamma_nu1 = arb_gamma_real(&order.add(&BigFloat::from_i32(1, wp), wp, rm), wp, rm, cc)?;
        let mut term = BigFloat::from_i32(1, wp).div(&gamma_nu1, wp, rm);
        sum = sum.add(&term, wp, rm);

        // The terms peak near k ≈ |x|/2 and then decay super-exponentially;
        // `|x| + wp` comfortably covers the tail down to 2^{-wp}.
        let max_terms = (x_f64.abs() * 1.5) as usize + wp + 40;
        for k in 1..=max_terms {
            // term *= (-x²/4) / (k · (ν + k))
            let k_bf = BigFloat::from_i32(k as i32, wp);
            let nu_plus_k = order.add(&k_bf, wp, rm);
            let denom = k_bf.mul(&nu_plus_k, wp, rm);
            term = term.mul(&neg_x_half_sq, wp, rm);
            term = term.div(&denom, wp, rm);

            // Convergence check (only once past the peak of the terms).
            if (k as f64) > x_f64.abs() / 2.0
                && let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
                && (s_exp as i64 - t_exp as i64) > wp as i64
            {
                tracing::trace!(k, "arb_bessel_j: series converged");
                break;
            }

            sum = sum.add(&term, wp, rm);
        }

        Ok(round_to(prefix.mul(&sum, wp, rm), prec, rm))
    }
}

/// Arbitrary-precision Bessel function of the second kind Y_ν(x), `x > 0`.
///
/// * Large `x` (see [`hankel_threshold`]): Hankel expansion
///   ([`bessel_jy_hankel`]), falling back to the series if it cannot reach
///   the target precision.
/// * Otherwise, integer order `n` only: the Neumann series (A&S 9.1.11)
///   with guard bits for the `e^{x}` cancellation.
fn arb_bessel_y(
    order: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let guard = 32;
    let wp = prec + guard;

    if x.is_zero() || x.is_negative() {
        return Err(SymplexError::Unevaluable {
            reason: "BesselY undefined at x ≤ 0".into(),
        });
    }

    let x_f64 = bigfloat_to_f64(x, rm, cc)?;
    let order_f64 = bigfloat_to_f64(order, rm, cc)?;
    let order_int = order_f64.round() as i64;
    let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;

    // Negative integer order: Y_{-n}(x) = (-1)^n Y_n(x).
    if is_int_order && order_int < 0 {
        let pos_order = BigFloat::from_i128((-order_int) as i128, wp);
        let y = arb_bessel_y(&pos_order, x, prec, rm, cc)?;
        return Ok(if order_int % 2 == 0 { y } else { y.neg() });
    }

    if x_f64 >= hankel_threshold(wp, order_f64.abs())
        && let Some(y) = bessel_jy_hankel(order, x, false, wp, rm, cc)
    {
        tracing::trace!("arb_bessel_y: used Hankel P/Q expansion");
        return Ok(round_to(y, prec, rm));
    }

    // ── Small |x|, integer order n ≥ 0: Neumann series (A&S 9.1.11) ─────
    //
    //   Y_n(x) = −(1/π)(x/2)^{−n} Σ_{k=0}^{n−1} (n−k−1)!/k! (x²/4)^k
    //            + (2/π) ln(x/2) J_n(x)
    //            − (1/π)(x/2)^n Σ_{k≥0} [ψ(k+1) + ψ(n+k+1)] (−x²/4)^k / (k!(n+k)!)
    //
    // with ψ(m+1) = −γ + H_m.  Computed directly for every n (no forward
    // recurrence or Wronskian division, which is unstable near zeros of J).
    // The log term and the series cancel to O(e^{x}) relative size, so a
    // few extra guard bits are added.

    if !is_int_order {
        return Err(SymplexError::Unevaluable {
            reason: format!(
                "BesselY for non-integer order {order_f64} not yet implemented in small-|x| regime"
            ),
        });
    }
    let n = order_int as usize;
    let wp = wp + cancellation_guard_bits(x_f64);
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);

    let pi = cc.pi(wp, rm).clone();
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    let inv_pi = one.div(&pi, wp, rm);
    let x_half = xw.div(&two, wp, rm);
    let x_half_sq = x_half.mul(&x_half, wp, rm);
    let ln_x_half = x_half.ln(wp, rm, cc);
    let euler_gamma = arb_euler_gamma(wp, rm, cc)?;
    let n_bf = BigFloat::from_i128(n as i128, wp);
    let j_n = arb_bessel_j(&n_bf, &xw, wp, rm, cc)?;

    // Finite sum: −(1/π)(x/2)^{−n} Σ_{k<n} (n−k−1)!/k! (x²/4)^k
    let mut finite = BigFloat::new(wp);
    if n > 0 {
        let mut f = Ratio::<BigInt>::from_integer(BigInt::from(1)); // (n−1)!
        for i in 2..n {
            f *= Ratio::from_integer(BigInt::from(i as u64));
        }
        let mut pow = one.clone();
        for k in 0..n {
            if k > 0 {
                f /= Ratio::from_integer(BigInt::from(((n - k) * k) as u64));
                pow = pow.mul(&x_half_sq, wp, rm);
            }
            let coeff = ratio_to_bigfloat(&f, wp, rm);
            finite = finite.add(&coeff.mul(&pow, wp, rm), wp, rm);
        }
        let x_half_pow_neg_n = one.div(&x_half.powi(n, wp, rm), wp, rm);
        finite = inv_pi
            .mul(&x_half_pow_neg_n, wp, rm)
            .mul(&finite, wp, rm)
            .neg();
    }

    // Log term: (2/π) ln(x/2) J_n(x)
    let log_term = two
        .mul(&inv_pi, wp, rm)
        .mul(&ln_x_half, wp, rm)
        .mul(&j_n, wp, rm);

    // Series: −(1/π)(x/2)^n Σ_k [ψ(k+1) + ψ(n+k+1)] (−x²/4)^k / (k!(n+k)!)
    let mut h_k = BigFloat::new(wp);
    let mut h_nk = BigFloat::new(wp);
    for m in 1..=n {
        h_nk = h_nk.add(
            &one.div(&BigFloat::from_i128(m as i128, wp), wp, rm),
            wp,
            rm,
        );
    }
    let mut n_fact = Ratio::<BigInt>::from_integer(BigInt::from(1));
    for i in 2..=n {
        n_fact *= Ratio::from_integer(BigInt::from(i as u64));
    }
    let mut inv_fact = one.div(&ratio_to_bigfloat(&n_fact, wp, rm), wp, rm);
    let neg_x_half_sq = x_half_sq.neg();
    let mut pow = one.clone();
    let mut series = BigFloat::new(wp);
    let two_gamma = euler_gamma.mul(&two, wp, rm);
    let max_terms = (x_f64.abs() * 2.0) as usize + wp + 40;
    for k in 0..=max_terms {
        if k > 0 {
            let k_bf = BigFloat::from_i128(k as i128, wp);
            let nk_bf = BigFloat::from_i128((n + k) as i128, wp);
            h_k = h_k.add(&one.div(&k_bf, wp, rm), wp, rm);
            h_nk = h_nk.add(&one.div(&nk_bf, wp, rm), wp, rm);
            inv_fact = inv_fact.div(&k_bf.mul(&nk_bf, wp, rm), wp, rm);
            pow = pow.mul(&neg_x_half_sq, wp, rm);
        }
        let psi_sum = h_k.add(&h_nk, wp, rm).sub(&two_gamma, wp, rm);
        let term = psi_sum.mul(&pow, wp, rm).mul(&inv_fact, wp, rm);
        series = series.add(&term, wp, rm);
        if k > 2
            && let (Some(t_exp), Some(s_exp)) = (term.exponent(), series.exponent())
            && (s_exp as i64 - t_exp as i64) > wp as i64
        {
            tracing::trace!(k, "arb_bessel_y: Neumann series converged");
            break;
        }
    }
    let series_term = inv_pi
        .mul(&x_half.powi(n, wp, rm), wp, rm)
        .mul(&series, wp, rm)
        .neg();

    let y = finite.add(&log_term, wp, rm).add(&series_term, wp, rm);
    Ok(round_to(y, prec, rm))
}

// ═══════════════════════════════════════════════════════════════════════════
// Catalan's constant
// ═══════════════════════════════════════════════════════════════════════════

/// Catalan's constant `G` at `prec` bits via Ramanujan's formula
///
/// ```text
/// G = (π/8)·ln(2 + √3) + (3/8)·Σ_{n≥0} 1 / ((2n+1)² · C(2n, n))
/// ```
///
/// The series gains ~2 bits per term (`1/C(2n,n) ~ 4^{-n}`), so about
/// `prec/2` terms are needed; every term is positive (no cancellation).
fn arb_catalan(prec: usize, rm: RoundingMode, cc: &mut Consts) -> Result<BigFloat, SymplexError> {
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    let three = BigFloat::from_i32(3, wp);
    let eight = BigFloat::from_i32(8, wp);

    // (π/8)·ln(2 + √3)
    let sqrt3 = three.sqrt(wp, rm);
    let ln_term = two.add(&sqrt3, wp, rm).ln(wp, rm, cc);
    let pi = cc.pi(wp, rm).clone();
    let first = pi.div(&eight, wp, rm).mul(&ln_term, wp, rm);

    // Σ b_n / (2n+1)² with b_n = 1/C(2n,n), b_{n+1} = b_n·(n+1)/(2(2n+1)).
    let mut b = one.clone();
    let mut sum = BigFloat::new(wp);
    let max_terms = wp / 2 + 64;
    for n in 0..max_terms {
        let two_n_plus_1 = BigFloat::from_i128(2 * n as i128 + 1, wp);
        let denom = two_n_plus_1.mul(&two_n_plus_1, wp, rm);
        let term = b.div(&denom, wp, rm);
        sum = sum.add(&term, wp, rm);
        if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
            && (s_exp as i64 - t_exp as i64) > wp as i64
        {
            break;
        }
        let n_plus_1 = BigFloat::from_i128(n as i128 + 1, wp);
        let two_2n_plus_1 = two.mul(&two_n_plus_1, wp, rm);
        b = b.mul(&n_plus_1, wp, rm).div(&two_2n_plus_1, wp, rm);
    }
    let second = three.div(&eight, wp, rm).mul(&sum, wp, rm);
    let mut g = first.add(&second, wp, rm);
    g.set_precision(prec, rm)
        .map_err(|e| SymplexError::ComputationFailed {
            operation: "catalan",
            reason: format!("precision adjustment failed: {e:?}"),
        })?;
    Ok(g)
}

// ═══════════════════════════════════════════════════════════════════════════
// Trigonometric / exponential / logarithmic integrals
// ═══════════════════════════════════════════════════════════════════════════

/// Number of guard bits needed to absorb the cancellation in an
/// alternating power series whose largest term is `~e^{|x|}`.
fn cancellation_guard_bits(x_abs: f64) -> usize {
    (x_abs * std::f64::consts::LOG2_E).ceil() as usize + 32
}

/// Auxiliary functions `f(x)`, `g(x)` of the asymptotic expansions
///
/// ```text
/// Si(x) = π/2 − f(x) cos x − g(x) sin x,   Ci(x) = f(x) sin x − g(x) cos x
/// f(x) ~ (1/x) Σ (−1)^k (2k)!/x^{2k},   g(x) ~ (1/x²) Σ (−1)^k (2k+1)!/x^{2k}
/// ```
///
/// for `x > 0`.  Uses optimal truncation (stop when terms grow); the
/// caller guarantees `x` is large enough that the smallest term is below
/// the working precision.
fn si_ci_asymptotic_fg(x: &BigFloat, wp: usize, rm: RoundingMode) -> (BigFloat, BigFloat) {
    let one = BigFloat::from_i32(1, wp);
    let inv_x = one.div(x, wp, rm);
    let inv_x2 = inv_x.mul(&inv_x, wp, rm);

    let mut f_sum = BigFloat::new(wp);
    let mut g_sum = BigFloat::new(wp);
    // t_f = (2k)!/x^{2k}, t_g = (2k+1)!/x^{2k}
    let mut t_f = one.clone();
    let mut t_g = one.clone();
    let mut prev_f: Option<BigFloat> = None;
    let max_terms = wp * 2 + 100;
    for k in 0..max_terms {
        if k > 0 {
            // (2k)! / (2k-2)! = (2k-1)(2k);  (2k+1)!/(2k-1)! = (2k)(2k+1)
            let a = BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp);
            let b = BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp);
            t_f = t_f.mul(&a, wp, rm).mul(&inv_x2, wp, rm);
            t_g = t_g.mul(&b, wp, rm).mul(&inv_x2, wp, rm);
        }
        // Optimal truncation: stop once terms start growing.
        if let Some(ref p) = prev_f
            && t_f.abs().cmp(p).is_some_and(|c| c > 0)
        {
            break;
        }
        if k % 2 == 0 {
            f_sum = f_sum.add(&t_f, wp, rm);
            g_sum = g_sum.add(&t_g, wp, rm);
        } else {
            f_sum = f_sum.sub(&t_f, wp, rm);
            g_sum = g_sum.sub(&t_g, wp, rm);
        }
        if let (Some(t_exp), Some(s_exp)) = (t_f.exponent(), f_sum.exponent())
            && (s_exp as i64 - t_exp as i64) > wp as i64
        {
            break;
        }
        prev_f = Some(t_f.abs());
    }
    (f_sum.mul(&inv_x, wp, rm), g_sum.mul(&inv_x2, wp, rm))
}

/// Sine and cosine integrals for real `x`.
///
/// * `want_si == true`: returns `(Si(x), 0)` — valid for all real `x`
///   (`Si` is odd).
/// * `want_si == false`: returns `(0, Ci(x))` — requires `x > 0`.
///
/// Power series for moderate `|x|` (with guard bits for the `e^{|x|}`
/// cancellation) and the Hankel-type asymptotic expansion when
/// `|x| > wp·ln 2` (where its optimal-truncation error `~e^{−x}` is below
/// the working precision).
fn arb_si_ci(
    x: &BigFloat,
    want_si: bool,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<(BigFloat, BigFloat), SymplexError> {
    if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() {
        return Err(SymplexError::Unevaluable {
            reason: "Si/Ci of special float value".into(),
        });
    }
    let x_f64 = bigfloat_to_f64(x, rm, cc)?;
    let x_abs = x_f64.abs();
    let base_wp = prec + 32;

    if x.is_zero() {
        return if want_si {
            Ok((BigFloat::new(prec), BigFloat::new(prec)))
        } else {
            Err(SymplexError::Unevaluable {
                reason: "Ci(0) is -∞".into(),
            })
        };
    }
    if !want_si && x.is_negative() {
        return Err(SymplexError::Unevaluable {
            reason: "Ci requires a positive argument here".into(),
        });
    }

    let switch = (base_wp as f64) * std::f64::consts::LN_2;
    if x_abs > switch {
        // ── Asymptotic expansion ──
        let wp = base_wp;
        let ax = x.abs();
        let (f, g) = si_ci_asymptotic_fg(&ax, wp, rm);
        let sin_x = ax.sin(wp, rm, cc);
        let cos_x = ax.cos(wp, rm, cc);
        let pi = cc.pi(wp, rm).clone();
        let two = BigFloat::from_i32(2, wp);
        let half_pi = pi.div(&two, wp, rm);
        let mut si =
            half_pi
                .sub(&f.mul(&cos_x, wp, rm), wp, rm)
                .sub(&g.mul(&sin_x, wp, rm), wp, rm);
        if x.is_negative() {
            si = si.neg();
        }
        let ci = f.mul(&sin_x, wp, rm).sub(&g.mul(&cos_x, wp, rm), wp, rm);
        return Ok((round_to(si, prec, rm), round_to(ci, prec, rm)));
    }

    // ── Power series ──
    let wp = base_wp + cancellation_guard_bits(x_abs);
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);
    let x2 = xw.mul(&xw, wp, rm);
    let max_terms = (x_abs * 1.5) as usize + wp + 40;

    let si = if want_si {
        // Si(x) = Σ (−1)^k t_k/(2k+1),  t_k = x^{2k+1}/(2k+1)!
        let mut t = xw.clone();
        let mut sum = BigFloat::new(wp);
        for k in 0..max_terms {
            if k > 0 {
                let d = BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp);
                t = t.mul(&x2, wp, rm).div(&d, wp, rm);
            }
            let term = t.div(&BigFloat::from_i128(2 * k as i128 + 1, wp), wp, rm);
            if k % 2 == 0 {
                sum = sum.add(&term, wp, rm);
            } else {
                sum = sum.sub(&term, wp, rm);
            }
            if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
                && (s_exp as i64 - t_exp as i64) > wp as i64
            {
                break;
            }
        }
        round_to(sum, prec, rm)
    } else {
        BigFloat::new(prec)
    };

    let ci = if want_si {
        BigFloat::new(prec)
    } else {
        // Ci(x) = γ + ln x + Σ_{k≥1} (−1)^k u_k/(2k),  u_k = x^{2k}/(2k)!
        let gamma = arb_euler_gamma(wp, rm, cc)?;
        let ln_x = xw.ln(wp, rm, cc);
        let mut sum = BigFloat::new(wp);
        let mut u = BigFloat::from_i32(1, wp);
        for k in 1..max_terms {
            let d = BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp);
            u = u.mul(&x2, wp, rm).div(&d, wp, rm);
            let term = u.div(&BigFloat::from_i128(2 * k as i128, wp), wp, rm);
            if k % 2 == 0 {
                sum = sum.add(&term, wp, rm);
            } else {
                sum = sum.sub(&term, wp, rm);
            }
            if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
                && (s_exp as i64 - t_exp as i64) > wp as i64
            {
                break;
            }
        }
        let ci = gamma.add(&ln_x, wp, rm).add(&sum, wp, rm);
        round_to(ci, prec, rm)
    };

    Ok((si, ci))
}

/// Exponential integral `Ei(x)` for real `x ≠ 0`.
///
/// * `|x| ≤ wp·ln 2`: `Ei(x) = γ + ln|x| + Σ_{k≥1} xᵏ/(k·k!)`, with guard
///   bits for the alternating cancellation when `x < 0`;
/// * otherwise the asymptotic expansion `Ei(x) ~ (eˣ/x) Σ k!/xᵏ` with
///   optimal truncation (relative error `~e^{−|x|}`).
fn arb_ei(
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() || x.is_zero() {
        return Err(SymplexError::Unevaluable {
            reason: "Ei of zero or special float value".into(),
        });
    }
    let x_f64 = bigfloat_to_f64(x, rm, cc)?;
    let x_abs = x_f64.abs();
    let base_wp = prec + 32;
    let switch = (base_wp as f64) * std::f64::consts::LN_2;

    if x_abs > switch {
        // ── Asymptotic: Ei(x) ~ (e^x / x) Σ_{k≥0} k! / x^k ──
        let wp = base_wp;
        let one = BigFloat::from_i32(1, wp);
        let inv_x = one.div(x, wp, rm);
        let mut term = one.clone();
        let mut sum = BigFloat::new(wp);
        let mut prev: Option<BigFloat> = None;
        let max_terms = wp * 2 + 100;
        for k in 0..max_terms {
            if k > 0 {
                term = term
                    .mul(&BigFloat::from_i128(k as i128, wp), wp, rm)
                    .mul(&inv_x, wp, rm);
            }
            if let Some(ref p) = prev
                && term.abs().cmp(p).is_some_and(|c| c > 0)
            {
                break;
            }
            sum = sum.add(&term, wp, rm);
            if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
                && (s_exp as i64 - t_exp as i64) > wp as i64
            {
                break;
            }
            prev = Some(term.abs());
        }
        let ex = x.exp(wp, rm, cc);
        let r = ex.mul(&inv_x, wp, rm).mul(&sum, wp, rm);
        return Ok(round_to(r, prec, rm));
    }

    // ── Power series ──
    let wp = if x.is_negative() {
        base_wp + cancellation_guard_bits(x_abs)
    } else {
        base_wp
    };
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);
    let gamma = arb_euler_gamma(wp, rm, cc)?;
    let ln_abs_x = xw.abs().ln(wp, rm, cc);
    let mut v = BigFloat::from_i32(1, wp); // x^k/k!
    let mut sum = BigFloat::new(wp);
    let max_terms = (x_abs * 3.0) as usize + wp + 40;
    for k in 1..max_terms {
        let k_bf = BigFloat::from_i128(k as i128, wp);
        v = v.mul(&xw, wp, rm).div(&k_bf, wp, rm);
        let term = v.div(&k_bf, wp, rm);
        sum = sum.add(&term, wp, rm);
        if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
            && (s_exp as i64 - t_exp as i64) > wp as i64
            && (k as f64) > x_abs
        {
            break;
        }
    }
    let r = gamma.add(&ln_abs_x, wp, rm).add(&sum, wp, rm);
    Ok(round_to(r, prec, rm))
}

/// Round a `BigFloat` down to `prec` bits (no-op on failure).
fn round_to(mut v: BigFloat, prec: usize, rm: RoundingMode) -> BigFloat {
    let _ = v.set_precision(prec, rm);
    v
}

/// `b^e` for real `b > 0` (or integer `e`), computed as `powi` for integer
/// exponents and `exp(e·ln b)` otherwise.
///
/// This deliberately avoids [`BigFloat::pow`], whose correct-rounding loop
/// never terminates when the exact result is representable in binary
/// (e.g. `4^{1/2} = 2`), which is exactly what happens for the integer
/// bases used in series like Borwein's `ζ` algorithm.
fn bf_pow(b: &BigFloat, e: &BigFloat, wp: usize, rm: RoundingMode, cc: &mut Consts) -> BigFloat {
    if e.is_zero() {
        return BigFloat::from_i32(1, wp);
    }
    if e.is_int()
        && let Ok(ef) = bigfloat_to_f64(e, rm, cc)
        && ef.abs() < 1.0e9
    {
        let n = ef.abs() as usize;
        let p = b.powi(n, wp, rm);
        return if ef < 0.0 {
            BigFloat::from_i32(1, wp).div(&p, wp, rm)
        } else {
            p
        };
    }
    let ln_b = b.ln(wp + 32, rm, cc);
    ln_b.mul(e, wp + 32, rm).exp(wp, rm, cc)
}

// ═══════════════════════════════════════════════════════════════════════════
// Riemann zeta function
// ═══════════════════════════════════════════════════════════════════════════

/// `ζ(s)` for real `s ≠ 1`.
///
/// * `s > 0`: Borwein's algorithm (P. Borwein, *An efficient algorithm for
///   the Riemann zeta function*, 1991, Algorithm 2) — an accelerated
///   alternating series with error `≤ 3/(3+√8)ⁿ`, so `n ≈ 0.39·wp` terms.
/// * `s < 0`: functional equation
///   `ζ(s) = 2ˢ π^{s−1} sin(πs/2) Γ(1−s) ζ(1−s)`; trivial zeros at the
///   negative even integers are returned exactly.
/// * `s = 0`: `−1/2`.
fn arb_zeta(
    s: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if s.is_nan() || s.is_inf_neg() {
        return Err(SymplexError::Unevaluable {
            reason: "zeta of special float value".into(),
        });
    }
    if s.is_inf_pos() {
        return Ok(BigFloat::from_i32(1, prec));
    }
    if s.is_zero() {
        return Ok(BigFloat::from_f64(-0.5, prec));
    }
    let s_f64 = bigfloat_to_f64(s, rm, cc)?;
    if (s_f64 - 1.0).abs() < 1e-300 {
        return Err(SymplexError::Unevaluable {
            reason: "zeta(1) is a pole".into(),
        });
    }
    let wp = prec + 32;

    if s.is_negative() {
        // Trivial zeros.
        let s_round = s_f64.round();
        if (s_f64 - s_round).abs() < 1e-14 && (s_round as i64) % 2 == 0 {
            let mut sw = s.clone();
            let _ = sw.set_precision(wp, rm);
            let diff = sw.sub(&BigFloat::from_f64(s_round, wp), wp, rm);
            if diff.is_zero() {
                return Ok(BigFloat::new(prec));
            }
        }
        // ζ(s) = 2^s π^{s−1} sin(π s/2) Γ(1−s) ζ(1−s)
        let mut sw = s.clone();
        let _ = sw.set_precision(wp, rm);
        let one = BigFloat::from_i32(1, wp);
        let two = BigFloat::from_i32(2, wp);
        let pi = cc.pi(wp, rm).clone();
        let one_minus_s = one.sub(&sw, wp, rm);
        let two_pow_s = bf_pow(&two, &sw, wp, rm, cc);
        let s_minus_1 = sw.sub(&one, wp, rm);
        let pi_pow = bf_pow(&pi, &s_minus_1, wp, rm, cc);
        let half_pi_s = pi.mul(&sw, wp, rm).div(&two, wp, rm);
        let sin_term = half_pi_s.sin(wp, rm, cc);
        let gamma_term = arb_gamma_real(&one_minus_s, wp, rm, cc)?;
        let zeta_term = arb_zeta_borwein(&one_minus_s, wp, rm, cc)?;
        let r = two_pow_s
            .mul(&pi_pow, wp, rm)
            .mul(&sin_term, wp, rm)
            .mul(&gamma_term, wp, rm)
            .mul(&zeta_term, wp, rm);
        return Ok(round_to(r, prec, rm));
    }

    let r = arb_zeta_borwein(s, wp, rm, cc)?;
    Ok(round_to(r, prec, rm))
}

/// Borwein's Algorithm 2 for `ζ(s)`, real `s > 0`, `s ≠ 1`.
///
/// ```text
/// d_k = n Σ_{i=0}^{k} (n+i−1)! 4ⁱ / ((n−i)! (2i)!)
/// ζ(s) = − 1/(d_n (1 − 2^{1−s})) · Σ_{k=0}^{n−1} (−1)ᵏ (d_k − d_n)/(k+1)ˢ
/// ```
fn arb_zeta_borwein(
    s: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    // n such that (3+√8)^{-n} < 2^{-wp}:  n > wp·ln2/ln(3+√8) ≈ 0.393·wp
    let n = ((wp as f64) * std::f64::consts::LN_2 / (3.0 + 8f64.sqrt()).ln()).ceil() as usize + 4;

    // Exact integer coefficients d_k via rational arithmetic.
    let mut d: Vec<Ratio<BigInt>> = Vec::with_capacity(n + 1);
    let mut acc = Ratio::<BigInt>::zero();
    // term_i = n · (n+i−1)! · 4^i / ((n−i)! (2i)!)
    // term_0 = n · (n−1)!/n! = 1
    let mut term = Ratio::from_integer(BigInt::from(1));
    for i in 0..=n {
        if i > 0 {
            // term_i / term_{i−1} = (n+i−1)·4·(n−i+1) / ((2i−1)(2i))
            let numer = BigInt::from((n + i - 1) as u64)
                * BigInt::from(4u64)
                * BigInt::from((n - i + 1) as u64);
            let denom = BigInt::from((2 * i - 1) as u64) * BigInt::from((2 * i) as u64);
            term *= Ratio::new(numer, denom);
        }
        acc += &term;
        d.push(acc.clone());
    }
    let d_n = ratio_to_bigfloat(&d[n], wp, rm);

    // S = Σ_{k=0}^{n−1} (−1)^k (d_k − d_n) / (k+1)^s
    let mut sum = BigFloat::new(wp);
    let s_is_int = s.is_int();
    let s_int = if s_is_int {
        bigfloat_to_f64(s, rm, cc)?.round() as usize
    } else {
        0
    };
    for k in 0..n {
        let coeff = ratio_to_bigfloat(&(&d[k] - &d[n]), wp, rm);
        let kp1 = BigFloat::from_i128(k as i128 + 1, wp);
        let kp1_pow_s = if s_is_int {
            kp1.powi(s_int, wp, rm)
        } else {
            bf_pow(&kp1, s, wp, rm, cc)
        };
        let term = coeff.div(&kp1_pow_s, wp, rm);
        if k % 2 == 0 {
            sum = sum.add(&term, wp, rm);
        } else {
            sum = sum.sub(&term, wp, rm);
        }
    }

    // 1 − 2^{1−s}
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    let one_minus_s = one.sub(s, wp, rm);
    let two_pow = bf_pow(&two, &one_minus_s, wp, rm, cc);
    let denom_factor = one.sub(&two_pow, wp, rm);
    if denom_factor.is_zero() {
        return Err(SymplexError::Unevaluable {
            reason: "zeta(1) is a pole".into(),
        });
    }
    let denom = d_n.mul(&denom_factor, wp, rm);
    Ok(sum.div(&denom, wp, rm).neg())
}

// ═══════════════════════════════════════════════════════════════════════════
// Polygamma
// ═══════════════════════════════════════════════════════════════════════════

/// `ψ⁽ⁿ⁾(x)` for integer `n ≥ 1` and real `x` (not a non-positive integer).
///
/// 1. Recurrence `ψ⁽ⁿ⁾(x) = ψ⁽ⁿ⁾(x+N) − (−1)ⁿ n! Σ_{k<N} 1/(x+k)^{n+1}`
///    shifts the argument to `x + N ≥ max(wp/3, n) + 10`.
/// 2. Asymptotic series (DLMF 5.15.8):
///    `ψ⁽ⁿ⁾(x) ~ (−1)^{n+1} [ (n−1)!/xⁿ + n!/(2x^{n+1}) + Σ_{k≥1} B₂ₖ (2k+n−1)!/((2k)! x^{2k+n}) ]`
///    with exact Bernoulli numbers, truncated when the terms fall below
///    the working precision (or start to diverge).
fn arb_polygamma(
    n: u32,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    _cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() {
        return Err(SymplexError::Unevaluable {
            reason: "polygamma of special float value".into(),
        });
    }
    let guard = 32;
    let wp = prec + guard;
    let n_us = n as usize;
    let one = BigFloat::from_i32(1, wp);

    // n! and (n−1)! as BigFloats.
    let mut n_fact = Ratio::<BigInt>::from_integer(BigInt::from(1));
    for i in 2..=n_us {
        n_fact *= Ratio::from_integer(BigInt::from(i as u64));
    }
    let nm1_fact = &n_fact / Ratio::from_integer(BigInt::from(n_us as u64));
    let n_fact_bf = ratio_to_bigfloat(&n_fact, wp, rm);
    let nm1_fact_bf = ratio_to_bigfloat(&nm1_fact, wp, rm);
    // (−1)^{n+1}
    let outer_sign_neg = n_us.is_multiple_of(2);

    // ── Step 1: shift x upward ──
    let threshold_val = ((wp / 3).max(n_us) + 10) as i128;
    let threshold = BigFloat::from_i128(threshold_val, wp);
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);
    let mut shift_sum = BigFloat::new(wp); // Σ 1/(x+k)^{n+1}
    let mut steps: usize = 0;
    while xw.sub(&threshold, wp, rm).is_negative() {
        if xw.is_zero() {
            return Err(SymplexError::Unevaluable {
                reason: "polygamma at non-positive integer pole".into(),
            });
        }
        if let Some(x_exp) = xw.exponent()
            && (x_exp as i64) < -(wp as i64 / 2)
        {
            return Err(SymplexError::Unevaluable {
                reason: "polygamma at non-positive integer pole".into(),
            });
        }
        let p = xw.powi(n_us + 1, wp, rm);
        let inv = one.div(&p, wp, rm);
        shift_sum = shift_sum.add(&inv, wp, rm);
        xw = xw.add(&one, wp, rm);
        steps += 1;
        if steps > 10 * wp + 1_000_000 {
            return Err(SymplexError::ComputationFailed {
                operation: "polygamma",
                reason: "argument shift did not terminate".into(),
            });
        }
    }
    // ψ⁽ⁿ⁾(x) = ψ⁽ⁿ⁾(x+N) + (−1)^{n+1} n! Σ 1/(x+k)^{n+1}
    let shift_term = n_fact_bf.mul(&shift_sum, wp, rm);

    // ── Step 2: asymptotic series at xw ──
    // bracket = (n−1)!/x^n + n!/(2 x^{n+1}) + Σ B_{2k} (2k+n−1)!/((2k)! x^{2k+n})
    let x_pow_n = xw.powi(n_us, wp, rm);
    let x_pow_np1 = x_pow_n.mul(&xw, wp, rm);
    let two = BigFloat::from_i32(2, wp);
    let mut bracket = nm1_fact_bf.div(&x_pow_n, wp, rm);
    bracket = bracket.add(&n_fact_bf.div(&two.mul(&x_pow_np1, wp, rm), wp, rm), wp, rm);

    let x2 = xw.mul(&xw, wp, rm);
    let mut x_pow = x_pow_n.mul(&x2, wp, rm); // x^{n+2} for k = 1
    // (2k+n−1)!/(2k)! as a running rational: start k=1 → (n+1)!/2!
    let mut ratio_fact = &n_fact * Ratio::from_integer(BigInt::from(n_us as u64 + 1))
        / Ratio::from_integer(BigInt::from(2));
    let mut prev_term_exp: Option<i32> = None;
    let max_terms = wp / 2 + 20;
    for k in 1..=max_terms {
        if k > 1 {
            // multiply by (2k+n−2)(2k+n−1) / ((2k−1)(2k))
            let a =
                BigInt::from((2 * k + n_us - 2) as u64) * BigInt::from((2 * k + n_us - 1) as u64);
            let b = BigInt::from((2 * k - 1) as u64) * BigInt::from((2 * k) as u64);
            ratio_fact *= Ratio::new(a, b);
            x_pow = x_pow.mul(&x2, wp, rm);
        }
        let b2k = crate::base::bernoulli::bernoulli(2 * k);
        let coeff = ratio_to_bigfloat(&(&b2k * &ratio_fact), wp, rm);
        let term = coeff.div(&x_pow, wp, rm);
        if let Some(t_exp) = term.exponent() {
            if let Some(prev) = prev_term_exp
                && t_exp > prev + 10
            {
                break; // diverging
            }
            prev_term_exp = Some(t_exp);
            if let Some(b_exp) = bracket.exponent()
                && (b_exp as i64 - t_exp as i64) > wp as i64
            {
                bracket = bracket.add(&term, wp, rm);
                break;
            }
        }
        bracket = bracket.add(&term, wp, rm);
    }

    let total = bracket.add(&shift_term, wp, rm);
    let result = if outer_sign_neg { total.neg() } else { total };
    Ok(round_to(result, prec, rm))
}

// ═══════════════════════════════════════════════════════════════════════════
// Modified Bessel functions
// ═══════════════════════════════════════════════════════════════════════════

/// Ascending series for `I_ν(x)` at working precision `wp`:
/// `I_ν(x) = (x/2)^ν Σ_{k≥0} (x²/4)^k / (k! Γ(ν+k+1))`.
///
/// All terms are positive for `x > 0`, so no cancellation occurs; the
/// number of terms is `O(|x| + wp)`.
fn bessel_i_series(
    order: &BigFloat,
    x: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let two = BigFloat::from_i32(2, wp);
    let x_half = x.div(&two, wp, rm);
    let x_half_sq = x_half.mul(&x_half, wp, rm);
    let order_f64 = bigfloat_to_f64(order, rm, cc)?;
    let order_int = order_f64.round() as i64;
    let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;

    if x.is_zero() {
        return if is_int_order && order_int == 0 {
            Ok(BigFloat::from_i32(1, wp))
        } else if order_f64 > 0.0 {
            Ok(BigFloat::new(wp))
        } else {
            Err(SymplexError::Unevaluable {
                reason: "BesselI at x=0 with negative order".into(),
            })
        };
    }

    // (x/2)^ν
    let prefix = if is_int_order && order_int >= 0 {
        x_half.powi(order_int as usize, wp, rm)
    } else if is_int_order {
        // Negative integer order: I_{-n} = I_n.
        x_half.powi((-order_int) as usize, wp, rm)
    } else {
        if x.is_negative() {
            return Err(SymplexError::Unevaluable {
                reason: "BesselI of negative argument with non-integer order is complex".into(),
            });
        }
        let ln_xh = x_half.ln(wp, rm, cc);
        order.mul(&ln_xh, wp, rm).exp(wp, rm, cc)
    };
    let eff_order = if is_int_order && order_int < 0 {
        BigFloat::from_i128((-order_int) as i128, wp)
    } else {
        order.clone()
    };

    // term_0 = 1/Γ(ν+1); term_{k+1} = term_k · (x²/4) / ((k+1)(ν+k+1))
    let one = BigFloat::from_i32(1, wp);
    let gamma_nu1 = arb_gamma_real(&eff_order.add(&one, wp, rm), wp, rm, cc)?;
    let mut term = one.div(&gamma_nu1, wp, rm);
    let mut sum = term.clone();
    let x_f64 = bigfloat_to_f64(x, rm, cc)?.abs();
    let max_terms = (x_f64 * 2.0) as usize + wp + 40;
    for k in 1..=max_terms {
        let k_bf = BigFloat::from_i128(k as i128, wp);
        let nu_plus_k = eff_order.add(&k_bf, wp, rm);
        let denom = k_bf.mul(&nu_plus_k, wp, rm);
        term = term.mul(&x_half_sq, wp, rm).div(&denom, wp, rm);
        sum = sum.add(&term, wp, rm);
        if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
            && (s_exp as i64 - t_exp as i64) > wp as i64
        {
            break;
        }
    }
    Ok(prefix.mul(&sum, wp, rm))
}

/// Modified Bessel function of the first kind `I_ν(x)` (real `ν`, `x`).
fn arb_bessel_i(
    order: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let wp = prec + 32;
    let r = bessel_i_series(order, x, wp, rm, cc)?;
    Ok(round_to(r, prec, rm))
}

/// Asymptotic expansion for `K_ν(x)`, large `x > 0` (DLMF 10.40.2):
/// `K_ν(x) ~ √(π/(2x)) e^{−x} Σ_k a_k(ν)/x^k`,
/// `a_k(ν) = ∏_{j=1}^{k} (4ν² − (2j−1)²) / (k! 8^k)`.
fn bessel_k_asymptotic(
    order: &BigFloat,
    x: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> BigFloat {
    let one = BigFloat::from_i32(1, wp);
    let four = BigFloat::from_i32(4, wp);
    let eight = BigFloat::from_i32(8, wp);
    let four_nu_sq = four.mul(&order.mul(order, wp, rm), wp, rm);
    let inv_x = one.div(x, wp, rm);
    let mut a_k = one.clone();
    let mut sum = one.clone();
    let mut prev: Option<BigFloat> = None;
    let max_terms = wp * 2 + 100;
    for k in 1..max_terms {
        let two_km1 = BigFloat::from_i32((2 * k as i32) - 1, wp);
        let numer = four_nu_sq.sub(&two_km1.mul(&two_km1, wp, rm), wp, rm);
        let denom = eight.mul(&BigFloat::from_i32(k as i32, wp), wp, rm);
        a_k = a_k
            .mul(&numer, wp, rm)
            .div(&denom, wp, rm)
            .mul(&inv_x, wp, rm);
        if let Some(ref p) = prev
            && a_k.abs().cmp(p).is_some_and(|c| c > 0)
        {
            break;
        }
        sum = sum.add(&a_k, wp, rm);
        if let (Some(t_exp), Some(s_exp)) = (a_k.exponent(), sum.exponent())
            && (s_exp as i64 - t_exp as i64) > wp as i64
        {
            break;
        }
        prev = Some(a_k.abs());
    }
    let pi = cc.pi(wp, rm).clone();
    let two = BigFloat::from_i32(2, wp);
    let amp = pi.div(&two.mul(x, wp, rm), wp, rm).sqrt(wp, rm);
    let e_neg_x = x.neg().exp(wp, rm, cc);
    amp.mul(&e_neg_x, wp, rm).mul(&sum, wp, rm)
}

/// Modified Bessel function of the second kind `K_ν(x)`, `x > 0`.
///
/// * Large `x` (`2x > wp·ln 2`): asymptotic expansion.
/// * Integer order `n` (A&S 9.6.11):
///   `K_n(x) = ½(x/2)^{−n} Σ_{k<n} (n−k−1)!/k! (−x²/4)^k + (−1)^{n+1} ln(x/2) I_n(x)
///           + (−1)^n ½ (x/2)^n Σ_{k≥0} [ψ(k+1) + ψ(n+k+1)] (x²/4)^k/(k!(n+k)!)`
///   with `ψ(m+1) = −γ + H_m`.
/// * Non-integer order: `K_ν = (π/2)(I_{−ν} − I_ν)/sin(νπ)`.
///
/// Both series formulas suffer `e^{x}`-scale cancellation (the result is
/// `~e^{−x}`), which is absorbed by extra guard bits.
fn arb_bessel_k(
    order: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_zero() || x.is_negative() || x.is_nan() {
        return Err(SymplexError::Unevaluable {
            reason: "BesselK undefined at x ≤ 0".into(),
        });
    }
    let base_wp = prec + 32;
    let x_f64 = bigfloat_to_f64(x, rm, cc)?;
    let order_f64 = bigfloat_to_f64(order, rm, cc)?;
    let order_int = order_f64.round() as i64;
    let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;
    let order_abs = if is_int_order {
        BigFloat::from_i128(order_int.unsigned_abs() as i128, base_wp)
    } else {
        order.abs()
    };

    // Large x: asymptotic (relative truncation error ~ e^{-2x}).
    if 2.0 * x_f64 > (base_wp as f64) * std::f64::consts::LN_2 {
        let r = bessel_k_asymptotic(&order_abs, x, base_wp, rm, cc);
        return Ok(round_to(r, prec, rm));
    }

    // Series: guard against e^{x} cancellation (two exponentially large
    // terms of opposite sign nearly cancel).
    let wp = base_wp + 2 * cancellation_guard_bits(x_f64);
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);

    if !is_int_order {
        // K_ν = (π/2)(I_{−ν} − I_ν)/sin(νπ)
        let neg_order = order.neg();
        let i_neg = bessel_i_series(&neg_order, &xw, wp, rm, cc)?;
        let i_pos = bessel_i_series(order, &xw, wp, rm, cc)?;
        let pi = cc.pi(wp, rm).clone();
        let two = BigFloat::from_i32(2, wp);
        let sin_nu_pi = order.mul(&pi, wp, rm).sin(wp, rm, cc);
        let r = pi
            .div(&two, wp, rm)
            .mul(&i_neg.sub(&i_pos, wp, rm), wp, rm)
            .div(&sin_nu_pi, wp, rm);
        return Ok(round_to(r, prec, rm));
    }

    let n = order_int.unsigned_abs() as usize;
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    let half = one.div(&two, wp, rm);
    let x_half = xw.div(&two, wp, rm);
    let x_half_sq = x_half.mul(&x_half, wp, rm);
    let ln_x_half = x_half.ln(wp, rm, cc);
    let gamma = arb_euler_gamma(wp, rm, cc)?;
    let n_bf = BigFloat::from_i128(n as i128, wp);
    let i_n = bessel_i_series(&n_bf, &xw, wp, rm, cc)?;

    // Finite sum: ½ (x/2)^{-n} Σ_{k=0}^{n-1} (n−k−1)!/k! (−x²/4)^k
    let mut finite = BigFloat::new(wp);
    if n > 0 {
        // f_k = (n−k−1)!/k! ; f_0 = (n−1)!
        let mut f = Ratio::<BigInt>::from_integer(BigInt::from(1));
        for i in 2..n {
            f *= Ratio::from_integer(BigInt::from(i as u64));
        }
        let mut pow = one.clone(); // (−x²/4)^k
        let neg_x_half_sq = x_half_sq.neg();
        for k in 0..n {
            if k > 0 {
                // f_k = f_{k−1} / ((n−k) · k)
                f /= Ratio::from_integer(BigInt::from(((n - k) * k) as u64));
                pow = pow.mul(&neg_x_half_sq, wp, rm);
            }
            let coeff = ratio_to_bigfloat(&f, wp, rm);
            finite = finite.add(&coeff.mul(&pow, wp, rm), wp, rm);
        }
        let x_half_pow_neg_n = one.div(&x_half.powi(n, wp, rm), wp, rm);
        finite = half.mul(&x_half_pow_neg_n, wp, rm).mul(&finite, wp, rm);
    }

    // Log term: (−1)^{n+1} ln(x/2) I_n(x)
    let log_term = ln_x_half.mul(&i_n, wp, rm);
    let log_term = if n.is_multiple_of(2) {
        log_term.neg()
    } else {
        log_term
    };

    // Infinite sum: (−1)^n ½ (x/2)^n Σ_k [ψ(k+1) + ψ(n+k+1)] (x²/4)^k/(k!(n+k)!)
    // ψ(m+1) = −γ + H_m.
    let mut h_k = BigFloat::new(wp); // H_0 = 0
    let mut h_nk = BigFloat::new(wp); // H_n
    for m in 1..=n {
        h_nk = h_nk.add(
            &one.div(&BigFloat::from_i128(m as i128, wp), wp, rm),
            wp,
            rm,
        );
    }
    // 1/(k!(n+k)!) start: 1/n!
    let mut n_fact = Ratio::<BigInt>::from_integer(BigInt::from(1));
    for i in 2..=n {
        n_fact *= Ratio::from_integer(BigInt::from(i as u64));
    }
    let mut inv_fact = one.div(&ratio_to_bigfloat(&n_fact, wp, rm), wp, rm);
    let mut pow = one.clone();
    let mut series = BigFloat::new(wp);
    let two_gamma = gamma.mul(&two, wp, rm);
    let max_terms = (x_f64 * 2.0) as usize + wp + 40;
    for k in 0..=max_terms {
        if k > 0 {
            let k_bf = BigFloat::from_i128(k as i128, wp);
            let nk_bf = BigFloat::from_i128((n + k) as i128, wp);
            h_k = h_k.add(&one.div(&k_bf, wp, rm), wp, rm);
            h_nk = h_nk.add(&one.div(&nk_bf, wp, rm), wp, rm);
            inv_fact = inv_fact.div(&k_bf.mul(&nk_bf, wp, rm), wp, rm);
            pow = pow.mul(&x_half_sq, wp, rm);
        }
        let psi_sum = h_k.add(&h_nk, wp, rm).sub(&two_gamma, wp, rm);
        let term = psi_sum.mul(&pow, wp, rm).mul(&inv_fact, wp, rm);
        series = series.add(&term, wp, rm);
        if k > 2
            && let (Some(t_exp), Some(s_exp)) = (term.exponent(), series.exponent())
            && (s_exp as i64 - t_exp as i64) > wp as i64
        {
            break;
        }
    }
    let series_term = half
        .mul(&x_half.powi(n, wp, rm), wp, rm)
        .mul(&series, wp, rm);
    let series_term = if n.is_multiple_of(2) {
        series_term
    } else {
        series_term.neg()
    };

    let r = finite.add(&log_term, wp, rm).add(&series_term, wp, rm);
    Ok(round_to(r, prec, rm))
}

// ═══════════════════════════════════════════════════════════════════════════
// Orthogonal polynomials
// ═══════════════════════════════════════════════════════════════════════════

/// Families of classical orthogonal polynomials evaluated numerically.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum OrthoPoly {
    Legendre,
    ChebyshevT,
    ChebyshevU,
    Hermite,
    Laguerre,
}

/// Evaluate `P_n(x)` (or `T_n`, `U_n`, `H_n`, `L_n`) for any non-negative
/// integer `n` via the three-term recurrence:
///
/// * Legendre: `(k+1) P_{k+1} = (2k+1) x P_k − k P_{k−1}`
/// * Chebyshev T/U: `T_{k+1} = 2x T_k − T_{k−1}` (`T_1 = x`, `U_1 = 2x`)
/// * Hermite (physicists'): `H_{k+1} = 2x H_k − 2k H_{k−1}`
/// * Laguerre: `(k+1) L_{k+1} = (2k+1−x) L_k − k L_{k−1}`
fn arb_orthopoly(kind: OrthoPoly, n: u64, x: &BigFloat, prec: usize, rm: RoundingMode) -> BigFloat {
    let wp = prec + 32 + (n as f64).log2().ceil().max(0.0) as usize;
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);
    if n == 0 {
        return BigFloat::from_i32(1, prec);
    }
    let mut p_prev = one.clone();
    let mut p_curr = match kind {
        OrthoPoly::Legendre | OrthoPoly::ChebyshevT => xw.clone(),
        OrthoPoly::ChebyshevU | OrthoPoly::Hermite => two.mul(&xw, wp, rm),
        OrthoPoly::Laguerre => one.sub(&xw, wp, rm),
    };
    for k in 1..n {
        let k_bf = BigFloat::from_i128(k as i128, wp);
        let kp1 = BigFloat::from_i128(k as i128 + 1, wp);
        let next = match kind {
            OrthoPoly::Legendre => {
                let two_k_plus_1 = BigFloat::from_i128(2 * k as i128 + 1, wp);
                let a = two_k_plus_1.mul(&xw, wp, rm).mul(&p_curr, wp, rm);
                let b = k_bf.mul(&p_prev, wp, rm);
                a.sub(&b, wp, rm).div(&kp1, wp, rm)
            }
            OrthoPoly::ChebyshevT | OrthoPoly::ChebyshevU => two
                .mul(&xw, wp, rm)
                .mul(&p_curr, wp, rm)
                .sub(&p_prev, wp, rm),
            OrthoPoly::Hermite => {
                let a = two.mul(&xw, wp, rm).mul(&p_curr, wp, rm);
                let b = two.mul(&k_bf, wp, rm).mul(&p_prev, wp, rm);
                a.sub(&b, wp, rm)
            }
            OrthoPoly::Laguerre => {
                let two_k_plus_1 = BigFloat::from_i128(2 * k as i128 + 1, wp);
                let a = two_k_plus_1.sub(&xw, wp, rm).mul(&p_curr, wp, rm);
                let b = k_bf.mul(&p_prev, wp, rm);
                a.sub(&b, wp, rm).div(&kp1, wp, rm)
            }
        };
        p_prev = p_curr;
        p_curr = next;
    }
    round_to(p_curr, prec, rm)
}

// ═══════════════════════════════════════════════════════════════════════════
// More special functions (0.9)
// ═══════════════════════════════════════════════════════════════════════════
//
// All routines take real arguments and return real results at `prec` bits;
// complex arguments and points outside the real domain of a function are
// refused with `Unevaluable`.  Working precision is `prec + 32` plus
// problem-dependent guard bits wherever an alternating series or a
// difference of nearly equal quantities loses digits.

fn unevaluable(reason: impl Into<String>) -> SymplexError {
    SymplexError::Unevaluable {
        reason: reason.into(),
    }
}

/// The cached values of `args`, all required to be real.
fn real_args(
    name: &str,
    args: &[ExprId],
    cache: &FxHashMap<ExprId, Complex>,
) -> Result<Vec<BigFloat>, SymplexError> {
    let mut out = Vec::with_capacity(args.len());
    for &a in args {
        let v = get_cached(cache, a)?;
        if !v.1.is_zero() {
            return Err(unevaluable(format!(
                "{name} of complex argument not yet supported in evalf"
            )));
        }
        out.push(v.0.clone());
    }
    Ok(out)
}

/// `|a| < 2^{-bits} |b|` — has the series term `a` become negligible
/// against the accumulated sum `b`?
fn negligible(a: &BigFloat, b: &BigFloat, bits: usize) -> bool {
    if a.is_zero() {
        return true;
    }
    match (a.exponent(), b.exponent()) {
        (Some(ae), Some(be)) => (be as i64 - ae as i64) > bits as i64,
        _ => false,
    }
}

/// `a > b`?
fn bf_gt(a: &BigFloat, b: &BigFloat) -> bool {
    a.cmp(b).is_some_and(|c| c > 0)
}

/// `a < b`?
fn bf_lt(a: &BigFloat, b: &BigFloat) -> bool {
    a.cmp(b).is_some_and(|c| c < 0)
}

/// Nearest integer to a `BigFloat`, if it is one (to within `1e-12`) and
/// fits in an `i64`.
fn bf_as_int(x: &BigFloat, rm: RoundingMode, cc: &mut Consts) -> Result<Option<i64>, SymplexError> {
    let f = bigfloat_to_f64(x, rm, cc)?;
    let r = f.round();
    if (f - r).abs() < 1e-12 && r.abs() < 9.0e15 {
        Ok(Some(r as i64))
    } else {
        Ok(None)
    }
}

/// Dispatch for the 0.9 special functions inside `eval_node`.
fn eval_special_09(
    name: &str,
    args: &[ExprId],
    cache: &FxHashMap<ExprId, Complex>,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<Complex, SymplexError> {
    use crate::base::arena::{
        FN_AIRYAI, FN_AIRYAIPRIME, FN_AIRYBI, FN_AIRYBIPRIME, FN_ASSOC_LAGUERRE, FN_ASSOC_LEGENDRE,
        FN_CHI, FN_DIRICHLET_ETA, FN_ELLIPTIC_E, FN_ELLIPTIC_F, FN_ELLIPTIC_K, FN_ELLIPTIC_PI,
        FN_ERFCINV, FN_ERFI, FN_ERFINV, FN_EXPINT, FN_FRESNELC, FN_FRESNELS, FN_GEGENBAUER,
        FN_JACOBI, FN_LOWERGAMMA, FN_POLYLOG, FN_SHI, FN_UPPERGAMMA,
    };
    let v = real_args(name, args, cache)?;
    let real = |r: BigFloat| Ok((r, BigFloat::new(prec)));
    debug!(prec, name, "evalf: special function (0.9)");
    match (name, v.len()) {
        (FN_ERFI, 1) => real(arb_erfi(&v[0], prec, rm, cc)?),
        (FN_ERFINV, 1) => real(arb_erfinv(&v[0], prec, rm, cc)?),
        (FN_ERFCINV, 1) => real(arb_erfcinv(&v[0], prec, rm, cc)?),
        (FN_EXPINT, 2) => real(arb_expint(&v[0], &v[1], prec, rm, cc)?),
        (FN_SHI, 1) => real(arb_shi_chi(&v[0], true, prec, rm, cc)?),
        (FN_CHI, 1) => {
            if v[0].is_zero() {
                return Err(unevaluable("Chi(0) is -∞"));
            }
            let chi = arb_shi_chi(&v[0].abs(), false, prec, rm, cc)?;
            if v[0].is_negative() {
                // Chi(−x) = Chi(x) + iπ
                Ok((chi, cc.pi(prec, rm).clone()))
            } else {
                real(chi)
            }
        }
        (FN_FRESNELS, 1) => real(arb_fresnel(&v[0], true, prec, rm, cc)?),
        (FN_FRESNELC, 1) => real(arb_fresnel(&v[0], false, prec, rm, cc)?),
        (FN_LOWERGAMMA, 2) => real(arb_lowergamma(&v[0], &v[1], prec, rm, cc)?),
        (FN_UPPERGAMMA, 2) => real(arb_uppergamma(&v[0], &v[1], prec, rm, cc)?),
        (FN_POLYLOG, 2) => real(arb_polylog(&v[0], &v[1], prec, rm, cc)?),
        (FN_DIRICHLET_ETA, 1) => real(arb_dirichlet_eta(&v[0], prec, rm, cc)?),
        (FN_AIRYAI, 1) => real(arb_airy(&v[0], AiryKind::Ai, prec, rm, cc)?),
        (FN_AIRYBI, 1) => real(arb_airy(&v[0], AiryKind::Bi, prec, rm, cc)?),
        (FN_AIRYAIPRIME, 1) => real(arb_airy(&v[0], AiryKind::AiPrime, prec, rm, cc)?),
        (FN_AIRYBIPRIME, 1) => real(arb_airy(&v[0], AiryKind::BiPrime, prec, rm, cc)?),
        (FN_ELLIPTIC_K, 1) => real(arb_elliptic_k(&v[0], prec, rm, cc)?),
        (FN_ELLIPTIC_E, 1) => real(arb_elliptic_e(&v[0], prec, rm, cc)?),
        (FN_ELLIPTIC_F, 2) => real(arb_elliptic_f(&v[0], &v[1], prec, rm, cc)?),
        (FN_ELLIPTIC_PI, 2) => real(arb_elliptic_pi(&v[0], &v[1], prec, rm, cc)?),
        (FN_GEGENBAUER, 3) => real(arb_gegenbauer(&v[0], &v[1], &v[2], prec, rm, cc)?),
        (FN_JACOBI, 4) => real(arb_jacobi(&v[0], &v[1], &v[2], &v[3], prec, rm, cc)?),
        (FN_ASSOC_LEGENDRE, 3) => real(arb_assoc_legendre(&v[0], &v[1], &v[2], prec, rm, cc)?),
        (FN_ASSOC_LAGUERRE, 3) => real(arb_assoc_laguerre(&v[0], &v[1], &v[2], prec, rm, cc)?),
        _ => Err(unevaluable(format!(
            "{name} called with {} argument(s)",
            v.len()
        ))),
    }
}

// ── Error-function family ──────────────────────────────────────────────────────

/// `erfi(x) = (2/√π) Σ x^{2n+1}/(n!(2n+1))` (all terms of one sign, so no
/// cancellation) for `x² ≤ wp·ln 2`, otherwise the asymptotic expansion
/// `erfi(x) ~ e^{x²}/(x√π) Σ (2k−1)!!/(2x²)^k` (optimal truncation error
/// `~e^{−x²}` relative to the leading term).
fn arb_erfi(
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_zero() {
        return Ok(BigFloat::new(prec));
    }
    if x.is_nan() || x.is_inf() {
        return Err(unevaluable("erfi of special float value"));
    }
    let wp = prec + 32;
    let x_f = bigfloat_to_f64(x, rm, cc)?;
    let ax = x_f.abs();
    let two = BigFloat::from_i32(2, wp);
    let sqrt_pi = cc.pi(wp, rm).clone().sqrt(wp, rm);
    let x_sq = x.mul(x, wp, rm);

    if ax * ax > (wp as f64) * std::f64::consts::LN_2 + 4.0 {
        let two_x_sq = x_sq.mul(&two, wp, rm);
        let mut sum = BigFloat::from_i32(1, wp);
        let mut term = BigFloat::from_i32(1, wp);
        let mut prev: Option<BigFloat> = None;
        for k in 1..(wp * 2 + 100) {
            let f = BigFloat::from_i128(2 * k as i128 - 1, wp);
            term = term.mul(&f, wp, rm).div(&two_x_sq, wp, rm);
            if let Some(ref p) = prev
                && bf_gt(&term, p)
            {
                break;
            }
            sum = sum.add(&term, wp, rm);
            if negligible(&term, &sum, wp) {
                break;
            }
            prev = Some(term.clone());
        }
        let e = x_sq.exp(wp, rm, cc);
        let denom = x.abs().mul(&sqrt_pi, wp, rm);
        let r = e.mul(&sum, wp, rm).div(&denom, wp, rm);
        let r = if x.is_negative() { r.neg() } else { r };
        return Ok(round_to(r, prec, rm));
    }

    let mut prod = x.clone();
    let _ = prod.set_precision(wp, rm);
    let mut sum = prod.clone();
    let max_terms = (3.0 * ax * ax) as usize + wp + 60;
    for n in 1..=max_terms {
        let n_bf = BigFloat::from_i128(n as i128, wp);
        prod = prod.mul(&x_sq, wp, rm).div(&n_bf, wp, rm);
        let d = BigFloat::from_i128(2 * n as i128 + 1, wp);
        let term = prod.div(&d, wp, rm);
        sum = sum.add(&term, wp, rm);
        if negligible(&term, &sum, wp) {
            break;
        }
    }
    let r = sum.mul(&two, wp, rm).div(&sqrt_pi, wp, rm);
    Ok(round_to(r, prec, rm))
}

/// `erfc(x)` with *relative* precision `prec` for `x > 0` (absolute
/// precision, via `1 − erf`, for `x ≤ 1/2`): the working precision is
/// raised by `x² log₂ e` bits to survive the cancellation in `1 − erf(x)`.
fn arb_erfc(
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let x_f = bigfloat_to_f64(x, rm, cc)?;
    let extra = if x_f > 0.5 {
        (x_f * x_f * std::f64::consts::LOG2_E).ceil() as usize + 8
    } else {
        0
    };
    let wp = prec + 32 + extra;
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);
    let erf = arb_erf(&xw, wp, rm, cc)?;
    let one = BigFloat::from_i32(1, wp);
    Ok(round_to(one.sub(&erf, wp, rm), prec, rm))
}

/// Halley iteration for `erf(x) = y` (`solve_erfc == false`) or
/// `erfc(x) = y` (`solve_erfc == true`) from the starting point `x0`.
///
/// With `g = f/f'`, the Halley step for both equations is
/// `x ← x − g / (1 + x g)` (since `f'' = −2x f'`).
fn erf_inverse_halley(
    y: &BigFloat,
    x0: f64,
    solve_erfc: bool,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let wp = prec + 32;
    let mut yw = y.clone();
    let _ = yw.set_precision(wp, rm);
    let half_sqrt_pi = cc
        .pi(wp, rm)
        .clone()
        .sqrt(wp, rm)
        .div(&BigFloat::from_i32(2, wp), wp, rm);
    let one = BigFloat::from_i32(1, wp);
    let mut x = BigFloat::from_f64(x0, wp);
    for _ in 0..200 {
        let f = if solve_erfc {
            arb_erfc(&x, wp, rm, cc)?.sub(&yw, wp, rm)
        } else {
            arb_erf(&x, wp, rm, cc)?.sub(&yw, wp, rm)
        };
        // f' = ±(2/√π) e^{−x²}  ⇒  g = f/f' = ±f (√π/2) e^{x²}
        let x_sq = x.mul(&x, wp, rm);
        let g = f
            .mul(&half_sqrt_pi, wp, rm)
            .mul(&x_sq.exp(wp, rm, cc), wp, rm);
        let g = if solve_erfc { g.neg() } else { g };
        let denom = one.add(&x.mul(&g, wp, rm), wp, rm);
        let delta = if denom.is_zero() {
            g.clone()
        } else {
            g.div(&denom, wp, rm)
        };
        x = x.sub(&delta, wp, rm);
        if negligible(&delta, &x, wp) {
            break;
        }
    }
    Ok(round_to(x, prec, rm))
}

/// `erfinv(y)` for `|y| < 1`: Winitzki's closed-form approximation refined
/// by `f64` Newton steps, then Halley iteration on `erf` (`|y| ≤ 1/2`) or
/// on `erfc(x) = 1 − |y|` (`|y| > 1/2`, where `erf` alone would lose the
/// digits of `1 − y`).
fn arb_erfinv(
    y: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if y.is_zero() {
        return Ok(BigFloat::new(prec));
    }
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let ay = y.abs();
    let one_minus_ay = one.sub(&ay, wp, rm);
    if one_minus_ay.is_zero() || one_minus_ay.is_negative() {
        return Err(unevaluable("erfinv: argument must satisfy |y| < 1"));
    }
    let y_f = bigfloat_to_f64(y, rm, cc)?;
    let r = if y_f.abs() > 0.5 {
        arb_erfcinv_small(&one_minus_ay, prec, rm, cc)?
    } else {
        // Winitzki (a = 0.147), |relative error| < 2·10⁻³.
        let a = 0.147_f64;
        let l = (1.0 - y_f * y_f).ln();
        let t = 2.0 / (std::f64::consts::PI * a) + l / 2.0;
        let mut x0 = ((t * t - l / a).sqrt() - t).sqrt().copysign(y_f);
        for _ in 0..3 {
            let f = erf_f64(x0) - y_f;
            x0 -= f * (std::f64::consts::PI.sqrt() / 2.0) * (x0 * x0).exp();
        }
        erf_inverse_halley(y, x0, false, prec, rm, cc)?
    };
    Ok(if y.is_negative() && !r.is_negative() && y_f.abs() > 0.5 {
        r.neg()
    } else {
        r
    })
}

/// `erfcinv(u)` for `0 < u ≤ 1/2`, Halley iteration on `erfc`.  The initial
/// guess inverts the leading asymptotic `erfc(x) ≈ e^{−x²}/(x√π)`, which
/// also covers `u` far below the `f64` range.
fn arb_erfcinv_small(
    u: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let wp = prec + 32;
    // L = −ln u  (computed in BigFloat so tiny u is fine)
    let l = u.ln(wp, rm, cc).neg();
    let l_f = bigfloat_to_f64(&l, rm, cc)?;
    let mut x0 = l_f.max(0.5).sqrt();
    for _ in 0..4 {
        // x² = L − ln(x√π)
        let v = l_f - (x0 * std::f64::consts::PI.sqrt()).ln();
        x0 = v.max(0.25).sqrt();
    }
    let u_f = bigfloat_to_f64(u, rm, cc)?;
    if u_f > 1e-300 {
        for _ in 0..3 {
            let f = (1.0 - erf_f64(x0)) - u_f;
            x0 += f * (std::f64::consts::PI.sqrt() / 2.0) * (x0 * x0).exp();
        }
    }
    erf_inverse_halley(u, x0, true, prec, rm, cc)
}

/// `erfcinv(u)` for `0 < u < 2`: `erfcinv(u) = −erfcinv(2 − u)` maps to
/// `u ≤ 1`, then `erfinv(1 − u)` for `u ≥ 1/2` and the `erfc` Halley
/// iteration for smaller `u`.
fn arb_erfcinv(
    u: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    if u.is_zero() || u.is_negative() || !bf_lt(u, &two) {
        return Err(unevaluable("erfcinv: argument must satisfy 0 < u < 2"));
    }
    if bf_gt(u, &one) {
        let v = two.sub(u, wp, rm);
        return Ok(arb_erfcinv(&v, prec, rm, cc)?.neg());
    }
    let u_f = bigfloat_to_f64(u, rm, cc)?;
    if u_f >= 0.5 {
        let y = one.sub(u, wp, rm);
        return arb_erfinv(&y, prec, rm, cc);
    }
    arb_erfcinv_small(u, prec, rm, cc)
}

// ── Incomplete gamma / exponential integrals ───────────────────────────────────

/// Legendre's continued fraction for `Γ(s, x)`, `x > 0` (modified Lentz):
///
/// ```text
/// Γ(s, x) = e^{−x} x^s / (x + 1 − s − 1·(1−s)/(x + 3 − s − 2·(2−s)/(x + 5 − s − …)))
/// ```
///
/// Converges for every `x > 0` and every real `s`; it is the method of
/// choice for `x ≥ max(1, s + 1)`.
fn uppergamma_cf(
    s: &BigFloat,
    x: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    // "tiny" replaces exact zeros in Lentz's algorithm.
    let mut tiny = BigFloat::from_i32(1, wp);
    tiny.set_exponent(-(4 * wp as i32));
    let mut b = x.add(&one, wp, rm).sub(s, wp, rm);
    let mut c = one.div(&tiny, wp, rm);
    let mut d = if b.is_zero() {
        one.div(&tiny, wp, rm)
    } else {
        one.div(&b, wp, rm)
    };
    let mut h = d.clone();
    let max_iter = 40 * wp + 2000;
    let mut converged = false;
    for i in 1..=max_iter {
        let i_bf = BigFloat::from_i128(i as i128, wp);
        // a_i = −i (i − s)
        let an = i_bf.mul(&i_bf.sub(s, wp, rm), wp, rm).neg();
        b = b.add(&two, wp, rm);
        d = an.mul(&d, wp, rm).add(&b, wp, rm);
        if d.is_zero() {
            d = tiny.clone();
        }
        c = b.add(&an.div(&c, wp, rm), wp, rm);
        if c.is_zero() {
            c = tiny.clone();
        }
        d = one.div(&d, wp, rm);
        let del = d.mul(&c, wp, rm);
        h = h.mul(&del, wp, rm);
        let dev = del.sub(&one, wp, rm);
        if negligible(&dev, &one, wp) {
            converged = true;
            break;
        }
    }
    if !converged {
        return Err(SymplexError::ComputationFailed {
            operation: "evalf",
            reason: "continued fraction for the incomplete gamma function did not converge".into(),
        });
    }
    let e_neg_x = x.neg().exp(wp, rm, cc);
    let x_pow_s = bf_pow(x, s, wp, rm, cc);
    Ok(e_neg_x.mul(&x_pow_s, wp, rm).mul(&h, wp, rm))
}

/// Series for `γ(s, x) = x^s e^{−x} Σ_{k≥0} x^k / (s (s+1) … (s+k))`,
/// `s > 0`, `x ≥ 0` (all terms positive).
fn lowergamma_series(
    s: &BigFloat,
    x: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_zero() {
        return Ok(BigFloat::new(wp));
    }
    let one = BigFloat::from_i32(1, wp);
    let x_f = bigfloat_to_f64(x, rm, cc)?;
    let mut term = one.div(s, wp, rm);
    let mut sum = term.clone();
    let max_terms = (2.0 * x_f) as usize + wp + 50;
    for k in 1..=max_terms {
        let s_plus_k = s.add(&BigFloat::from_i128(k as i128, wp), wp, rm);
        term = term.mul(x, wp, rm).div(&s_plus_k, wp, rm);
        sum = sum.add(&term, wp, rm);
        if negligible(&term, &sum, wp) {
            break;
        }
    }
    let e_neg_x = x.neg().exp(wp, rm, cc);
    let x_pow_s = bf_pow(x, s, wp, rm, cc);
    Ok(x_pow_s.mul(&e_neg_x, wp, rm).mul(&sum, wp, rm))
}

/// `Γ(s, x)` at working precision `wp` for real `s`, `x > 0`.
///
/// * `x ≥ max(1, s+1)`: continued fraction.
/// * `s > 0`, smaller `x`: `Γ(s) − γ(s, x)` (series), with guard bits for
///   the cancellation when `Γ(s, x) ≪ Γ(s)`.
/// * `s = 0`, `x < 1`: `Γ(0, x) = E₁(x) = −Ei(−x)`.
/// * `s < 0`, `x < 1`: shift `s` into `(0, 1]` (or to `0` for integers) by
///   `Γ(s, x) = (Γ(s+1, x) − x^s e^{−x}) / s`, downward from the base value.
fn uppergamma_wp(
    s: &BigFloat,
    x: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let s_f = bigfloat_to_f64(s, rm, cc)?;
    let x_f = bigfloat_to_f64(x, rm, cc)?;
    if x_f >= (s_f + 1.0).max(1.0) {
        return uppergamma_cf(s, x, wp, rm, cc);
    }
    if s_f > 0.0 {
        // Γ(s, x)/Γ(s) is not tiny here (x < s + 1), so a fixed guard suffices.
        let wp2 = wp + 32 + cancellation_guard_bits(x_f);
        let mut sw = s.clone();
        let _ = sw.set_precision(wp2, rm);
        let mut xw = x.clone();
        let _ = xw.set_precision(wp2, rm);
        let g = arb_gamma_real(&sw, wp2, rm, cc)?;
        let lower = lowergamma_series(&sw, &xw, wp2, rm, cc)?;
        return Ok(round_to(g.sub(&lower, wp2, rm), wp, rm));
    }
    // s ≤ 0, x < 1.  Each downward step `Γ(s−1) = (Γ(s) − x^{s−1}e^{−x})/(s−1)`
    // may cancel a few bits; reserve 16 per step.
    let s_int = bf_as_int(s, rm, cc)?;
    let steps = match s_int {
        Some(k) => (-k) as usize,
        None => (-s_f).ceil() as usize,
    };
    let wp2 = wp + 16 * steps + 32;
    let mut sw = s.clone();
    let _ = sw.set_precision(wp2, rm);
    let mut xw = x.clone();
    let _ = xw.set_precision(wp2, rm);
    let e_neg_x = xw.neg().exp(wp2, rm, cc);
    // Base point: Γ(0, x) = E₁(x) = −Ei(−x) for integer s; otherwise
    // Γ(s₀, x) = Γ(s₀) − γ(s₀, x) with s₀ = s + steps ∈ (0, 1).
    let (mut cur_s, mut val) = match s_int {
        Some(_) => (BigFloat::new(wp2), arb_ei(&xw.neg(), wp2, rm, cc)?.neg()),
        None => {
            let s0 = sw.add(&BigFloat::from_i128(steps as i128, wp2), wp2, rm);
            let g = arb_gamma_real(&s0, wp2, rm, cc)?;
            let lower = lowergamma_series(&s0, &xw, wp2, rm, cc)?;
            (s0, g.sub(&lower, wp2, rm))
        }
    };
    let one = BigFloat::from_i32(1, wp2);
    for _ in 0..steps {
        let next_s = cur_s.sub(&one, wp2, rm);
        let x_pow = bf_pow(&xw, &next_s, wp2, rm, cc);
        let t = x_pow.mul(&e_neg_x, wp2, rm);
        val = val.sub(&t, wp2, rm).div(&next_s, wp2, rm);
        cur_s = next_s;
    }
    Ok(round_to(val, wp, rm))
}

/// Upper incomplete gamma `Γ(s, x)` for real `s`, `x ≥ 0`.
fn arb_uppergamma(
    s: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || s.is_nan() || x.is_inf_neg() || s.is_inf() {
        return Err(unevaluable("uppergamma of special float value"));
    }
    if x.is_inf_pos() {
        return Ok(BigFloat::new(prec));
    }
    if x.is_negative() {
        return Err(unevaluable(
            "uppergamma of negative argument not yet supported in evalf",
        ));
    }
    let wp = prec + 32;
    if x.is_zero() {
        if !s.is_positive() {
            return Err(unevaluable("uppergamma(s, 0) diverges for s ≤ 0"));
        }
        return Ok(round_to(arb_gamma_real(s, wp, rm, cc)?, prec, rm));
    }
    let r = uppergamma_wp(s, x, wp, rm, cc)?;
    Ok(round_to(r, prec, rm))
}

/// Lower incomplete gamma `γ(s, x)` for real `s > 0`, `x ≥ 0`.
fn arb_lowergamma(
    s: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || s.is_nan() || s.is_inf() || x.is_inf_neg() {
        return Err(unevaluable("lowergamma of special float value"));
    }
    if !s.is_positive() {
        return Err(unevaluable(
            "lowergamma requires s > 0 for numerical evaluation",
        ));
    }
    if x.is_negative() {
        return Err(unevaluable(
            "lowergamma of negative argument not yet supported in evalf",
        ));
    }
    let wp = prec + 32;
    if x.is_inf_pos() {
        return Ok(round_to(arb_gamma_real(s, wp, rm, cc)?, prec, rm));
    }
    let s_f = bigfloat_to_f64(s, rm, cc)?;
    let x_f = bigfloat_to_f64(x, rm, cc)?;
    if x_f < s_f + 1.0 {
        let r = lowergamma_series(s, x, wp, rm, cc)?;
        return Ok(round_to(r, prec, rm));
    }
    // γ = Γ(s) − Γ(s, x); Γ(s, x) ≤ Γ(s) here so the difference is benign.
    let g = arb_gamma_real(s, wp, rm, cc)?;
    let upper = uppergamma_cf(s, x, wp, rm, cc)?;
    Ok(round_to(g.sub(&upper, wp, rm), prec, rm))
}

/// Generalised exponential integral `E_n(x) = x^{n−1} Γ(1−n, x)` for real
/// `n` and `x > 0` (`E_n(0) = 1/(n−1)` for `n > 1`).
fn arb_expint(
    n: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || n.is_nan() || n.is_inf() {
        return Err(unevaluable("expint of special float value"));
    }
    if x.is_inf_pos() {
        return Ok(BigFloat::new(prec));
    }
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    if x.is_zero() {
        let n_minus_1 = n.sub(&one, wp, rm);
        if !n_minus_1.is_positive() {
            return Err(unevaluable("expint(n, 0) diverges for n ≤ 1"));
        }
        return Ok(round_to(one.div(&n_minus_1, wp, rm), prec, rm));
    }
    if x.is_negative() {
        return Err(unevaluable(
            "expint of negative argument is complex; not yet supported in evalf",
        ));
    }
    let s = one.sub(n, wp, rm);
    let g = uppergamma_wp(&s, x, wp, rm, cc)?;
    let x_pow = bf_pow(x, &n.sub(&one, wp, rm), wp, rm, cc);
    Ok(round_to(x_pow.mul(&g, wp, rm), prec, rm))
}

/// `Shi(x)` (`want_shi`) or `Chi(x)` for `x > 0` (`Shi` is odd and is
/// returned for any real `x`).
///
/// * `|x| ≤ 40`: `Shi(x) = Σ x^{2k+1}/((2k+1)(2k+1)!)`,
///   `Chi(x) = γ + ln x + Σ_{k≥1} x^{2k}/((2k)(2k)!)` — all terms positive.
/// * larger `|x|`: `Shi = (Ei(x) − Ei(−x))/2`, `Chi = (Ei(x) + Ei(−x))/2`.
fn arb_shi_chi(
    x: &BigFloat,
    want_shi: bool,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || x.is_inf() {
        return Err(unevaluable("Shi/Chi of special float value"));
    }
    if x.is_zero() {
        return if want_shi {
            Ok(BigFloat::new(prec))
        } else {
            Err(unevaluable("Chi(0) is -∞"))
        };
    }
    let wp = prec + 32;
    let ax = x.abs();
    let x_f = bigfloat_to_f64(&ax, rm, cc)?;
    let r = if x_f > 40.0 {
        let ei_pos = arb_ei(&ax, wp, rm, cc)?;
        let ei_neg = arb_ei(&ax.neg(), wp, rm, cc)?;
        let two = BigFloat::from_i32(2, wp);
        if want_shi {
            ei_pos.sub(&ei_neg, wp, rm).div(&two, wp, rm)
        } else {
            ei_pos.add(&ei_neg, wp, rm).div(&two, wp, rm)
        }
    } else {
        let mut xw = ax.clone();
        let _ = xw.set_precision(wp, rm);
        let x2 = xw.mul(&xw, wp, rm);
        let max_terms = (x_f * 1.5) as usize + wp + 40;
        if want_shi {
            // t_k = x^{2k+1}/(2k+1)!
            let mut t = xw.clone();
            let mut sum = BigFloat::new(wp);
            for k in 0..max_terms {
                if k > 0 {
                    let d = BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp);
                    t = t.mul(&x2, wp, rm).div(&d, wp, rm);
                }
                let term = t.div(&BigFloat::from_i128(2 * k as i128 + 1, wp), wp, rm);
                sum = sum.add(&term, wp, rm);
                if negligible(&term, &sum, wp) {
                    break;
                }
            }
            sum
        } else {
            let gamma = arb_euler_gamma(wp, rm, cc)?;
            let ln_x = xw.ln(wp, rm, cc);
            let mut u = BigFloat::from_i32(1, wp); // x^{2k}/(2k)!
            let mut sum = BigFloat::new(wp);
            for k in 1..max_terms {
                let d = BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp);
                u = u.mul(&x2, wp, rm).div(&d, wp, rm);
                let term = u.div(&BigFloat::from_i128(2 * k as i128, wp), wp, rm);
                sum = sum.add(&term, wp, rm);
                if negligible(&term, &sum, wp) {
                    break;
                }
            }
            gamma.add(&ln_x, wp, rm).add(&sum, wp, rm)
        }
    };
    let r = if want_shi && x.is_negative() {
        r.neg()
    } else {
        r
    };
    Ok(round_to(r, prec, rm))
}

/// Fresnel integrals `S(x)` (`want_s`) and `C(x)` for real `x` (both odd).
///
/// * `πx²/2 ≤ wp·ln 2`: Maclaurin series
///   `S = Σ (−1)^k (π/2)^{2k+1} x^{4k+3}/((2k+1)!(4k+3))`,
///   `C = Σ (−1)^k (π/2)^{2k} x^{4k+1}/((2k)!(4k+1))`, with guard bits for
///   the `e^{πx²/2}` cancellation.
/// * beyond: DLMF 7.12 auxiliary functions
///   `S = 1/2 − f cos(πx²/2) − g sin(πx²/2)`, `C = 1/2 + f sin(πx²/2) − g cos(πx²/2)`,
///   `f ~ (1/πx) Σ (−1)^m (4m−1)!!/(πx²)^{2m}`, `g ~ (1/πx) Σ (−1)^m (4m+1)!!/(πx²)^{2m+1}`
///   (optimal truncation error `~e^{−πx²/2}`).
fn arb_fresnel(
    x: &BigFloat,
    want_s: bool,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || x.is_inf() {
        return Err(unevaluable("Fresnel integral of special float value"));
    }
    if x.is_zero() {
        return Ok(BigFloat::new(prec));
    }
    let base_wp = prec + 32;
    let ax = x.abs();
    let x_f = bigfloat_to_f64(&ax, rm, cc)?;
    let arg_f = std::f64::consts::FRAC_PI_2 * x_f * x_f;

    let r = if arg_f > (base_wp as f64) * std::f64::consts::LN_2 {
        let wp = base_wp;
        let mut xw = ax.clone();
        let _ = xw.set_precision(wp, rm);
        let one = BigFloat::from_i32(1, wp);
        let pi = cc.pi(wp, rm).clone();
        let pi_x2 = pi.mul(&xw, wp, rm).mul(&xw, wp, rm);
        let inv = one.div(&pi_x2, wp, rm); // 1/(πx²)
        let inv2 = inv.mul(&inv, wp, rm);
        // f-series: (4m−1)!!/(πx²)^{2m}; g-series: (4m+1)!!/(πx²)^{2m+1}
        let mut tf = one.clone();
        let mut tg = inv.clone();
        let mut f_sum = BigFloat::new(wp);
        let mut g_sum = BigFloat::new(wp);
        let mut prev: Option<BigFloat> = None;
        for m in 0..(wp * 2 + 100) {
            if m > 0 {
                // (4m−1)!!/(4m−5)!! = (4m−3)(4m−1); (4m+1)!!/(4m−3)!! = (4m−1)(4m+1)
                let a = BigFloat::from_i128((4 * m as i128 - 3) * (4 * m as i128 - 1), wp);
                let b = BigFloat::from_i128((4 * m as i128 - 1) * (4 * m as i128 + 1), wp);
                tf = tf.mul(&a, wp, rm).mul(&inv2, wp, rm);
                tg = tg.mul(&b, wp, rm).mul(&inv2, wp, rm);
            }
            if let Some(ref p) = prev
                && bf_gt(&tf, p)
            {
                break;
            }
            if m % 2 == 0 {
                f_sum = f_sum.add(&tf, wp, rm);
                g_sum = g_sum.add(&tg, wp, rm);
            } else {
                f_sum = f_sum.sub(&tf, wp, rm);
                g_sum = g_sum.sub(&tg, wp, rm);
            }
            if negligible(&tf, &one, wp) {
                break;
            }
            prev = Some(tf.clone());
        }
        let pi_x = pi.mul(&xw, wp, rm);
        let f = f_sum.div(&pi_x, wp, rm);
        let g = g_sum.div(&pi_x, wp, rm);
        let two = BigFloat::from_i32(2, wp);
        let theta = pi_x2.div(&two, wp, rm);
        let s = theta.sin(wp, rm, cc);
        let c = theta.cos(wp, rm, cc);
        let half = one.div(&two, wp, rm);
        if want_s {
            half.sub(&f.mul(&c, wp, rm), wp, rm)
                .sub(&g.mul(&s, wp, rm), wp, rm)
        } else {
            half.add(&f.mul(&s, wp, rm), wp, rm)
                .sub(&g.mul(&c, wp, rm), wp, rm)
        }
    } else {
        let wp = base_wp + (arg_f * std::f64::consts::LOG2_E).ceil() as usize + 16;
        let mut xw = ax.clone();
        let _ = xw.set_precision(wp, rm);
        let pi = cc.pi(wp, rm).clone();
        let two = BigFloat::from_i32(2, wp);
        let half_pi = pi.div(&two, wp, rm);
        let x2 = xw.mul(&xw, wp, rm);
        // u = (π/2) x²;  terms: S_k = (−1)^k u^{2k+1} x /((2k+1)! (4k+3)),
        //                          C_k = (−1)^k u^{2k} x /((2k)! (4k+1))
        let u = half_pi.mul(&x2, wp, rm);
        let u2 = u.mul(&u, wp, rm);
        let mut t = if want_s {
            u.mul(&xw, wp, rm)
        } else {
            xw.clone()
        };
        let mut sum = BigFloat::new(wp);
        let max_terms = (arg_f * 2.0) as usize + wp + 40;
        for k in 0..max_terms {
            if k > 0 {
                let d = if want_s {
                    BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp)
                } else {
                    BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp)
                };
                t = t.mul(&u2, wp, rm).div(&d, wp, rm);
            }
            let div = BigFloat::from_i128(4 * k as i128 + if want_s { 3 } else { 1 }, wp);
            let term = t.div(&div, wp, rm);
            if k % 2 == 0 {
                sum = sum.add(&term, wp, rm);
            } else {
                sum = sum.sub(&term, wp, rm);
            }
            if negligible(&term, &sum, wp) {
                break;
            }
        }
        sum
    };
    let r = if x.is_negative() { r.neg() } else { r };
    Ok(round_to(r, prec, rm))
}

// ── Polylogarithm / Dirichlet eta ───────────────────────────────────────────────

/// `ζ(m)` for an integer `m ≤ 0`, exactly: `ζ(0) = −1/2`,
/// `ζ(−k) = −B_{k+1}/(k+1)`.
fn zeta_nonpositive_int(m: i64, wp: usize, rm: RoundingMode) -> BigFloat {
    if m == 0 {
        return BigFloat::from_f64(-0.5, wp);
    }
    let k = (-m) as usize;
    let b = crate::base::bernoulli::bernoulli(k + 1);
    if b.is_zero() {
        return BigFloat::new(wp);
    }
    let v = -b / Ratio::from_integer(BigInt::from(k as i64 + 1));
    ratio_to_bigfloat(&v, wp, rm)
}

/// `k^s` for integer `k ≥ 1`.
fn int_pow_s(
    k: usize,
    s: &BigFloat,
    s_int: Option<i64>,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> BigFloat {
    let kb = BigFloat::from_i128(k as i128, wp);
    match s_int {
        Some(n) if n >= 0 => kb.powi(n as usize, wp, rm),
        Some(n) => BigFloat::from_i32(1, wp).div(&kb.powi((-n) as usize, wp, rm), wp, rm),
        None => bf_pow(&kb, s, wp, rm, cc),
    }
}

/// Direct series `Li_s(z) = Σ_{k≥1} z^k/k^s` for `|z| < 1` (used up to
/// [`polylog_series_limit`]); about `wp·ln 2 / (−ln|z|)` terms.
fn polylog_series(
    s: &BigFloat,
    s_int: Option<i64>,
    z: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> BigFloat {
    let z_f = bigfloat_to_f64(z, rm, cc).unwrap_or(0.5).abs();
    let max_terms = if z_f > 0.0 && z_f < 1.0 {
        ((wp as f64) * std::f64::consts::LN_2 / (-z_f.ln())).ceil() as usize + 20
    } else {
        wp * 2 + 60
    };
    let mut zk = z.clone();
    let mut sum = BigFloat::new(wp);
    for k in 1..max_terms {
        if k > 1 {
            zk = zk.mul(z, wp, rm);
        }
        let term = zk.div(&int_pow_s(k, s, s_int, wp, rm, cc), wp, rm);
        sum = sum.add(&term, wp, rm);
        if negligible(&term, &sum, wp) {
            break;
        }
    }
    sum
}

/// `Li_{−n}(z)` for integer `n ≥ 0` and any real `z ≠ 1`:
/// `Σ_{k=0}^{n} k! S(n+1, k+1) (z/(1−z))^{k+1}`.
fn polylog_nonpositive(
    n: usize,
    z: &BigFloat,
    wp: usize,
    rm: RoundingMode,
) -> Result<BigFloat, SymplexError> {
    let one = BigFloat::from_i32(1, wp);
    let one_minus_z = one.sub(z, wp, rm);
    if one_minus_z.is_zero() {
        return Err(unevaluable("polylog(s, 1) diverges for s ≤ 1"));
    }
    let w = z.div(&one_minus_z, wp, rm);
    let mut k_fact = BigInt::from(1);
    let mut w_pow = w.clone();
    let mut sum = BigFloat::new(wp);
    for k in 0..=n {
        if k > 0 {
            k_fact *= BigInt::from(k as u64);
            w_pow = w_pow.mul(&w, wp, rm);
        }
        let s2 = crate::domains::combinatorics::stirling2(n as u64 + 1, k as u64 + 1)
            .ok_or_else(|| unevaluable("polylog: Stirling number overflow"))?;
        let c = ratio_to_bigfloat(&Ratio::from_integer(&k_fact * s2), wp, rm);
        sum = sum.add(&c.mul(&w_pow, wp, rm), wp, rm);
    }
    Ok(sum)
}

/// `Li_s(z)` for `1/2 < z < 1` via the expansion in `μ = ln z` (`|μ| < 2π`):
///
/// * integer `s = n ≥ 1`:
///   `Li_n(e^μ) = μ^{n−1}/(n−1)! (H_{n−1} − ln(−μ)) + Σ_{k≥0, k≠n−1} ζ(n−k) μ^k/k!`
/// * non-integer `s`: `Li_s(e^μ) = Γ(1−s)(−μ)^{s−1} + Σ_{k≥0} ζ(s−k) μ^k/k!`
fn polylog_near_one(
    s: &BigFloat,
    s_int: Option<i64>,
    z: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let one = BigFloat::from_i32(1, wp);
    let mu = z.ln(wp, rm, cc); // negative
    let neg_mu = mu.neg();
    let max_terms = wp * 2 + 100;
    let mut sum = BigFloat::new(wp);
    let mut mu_pow_over_fact = one.clone(); // μ^k/k!
    match s_int {
        Some(n) if n >= 1 => {
            let n = n as usize;
            for k in 0..max_terms {
                if k > 0 {
                    mu_pow_over_fact = mu_pow_over_fact.mul(&mu, wp, rm).div(
                        &BigFloat::from_i128(k as i128, wp),
                        wp,
                        rm,
                    );
                }
                if k + 1 == n {
                    // μ^{n−1}/(n−1)! (H_{n−1} − ln(−μ))
                    let mut h = BigFloat::new(wp);
                    for j in 1..n {
                        h = h.add(
                            &one.div(&BigFloat::from_i128(j as i128, wp), wp, rm),
                            wp,
                            rm,
                        );
                    }
                    let ln_neg_mu = neg_mu.ln(wp, rm, cc);
                    let t = mu_pow_over_fact.mul(&h.sub(&ln_neg_mu, wp, rm), wp, rm);
                    sum = sum.add(&t, wp, rm);
                    continue;
                }
                let arg = n as i64 - k as i64;
                let zeta = if arg <= 0 {
                    zeta_nonpositive_int(arg, wp, rm)
                } else {
                    arb_zeta(&BigFloat::from_i128(arg as i128, wp), wp, rm, cc)?
                };
                let term = zeta.mul(&mu_pow_over_fact, wp, rm);
                sum = sum.add(&term, wp, rm);
                if k > n && negligible(&mu_pow_over_fact, &one, wp + 8) {
                    break;
                }
            }
            Ok(sum)
        }
        _ => {
            let one_minus_s = one.sub(s, wp, rm);
            let g = arb_gamma_real(&one_minus_s, wp, rm, cc)?;
            let s_minus_1 = s.sub(&one, wp, rm);
            let lead = g.mul(&bf_pow(&neg_mu, &s_minus_1, wp, rm, cc), wp, rm);
            for k in 0..max_terms {
                if k > 0 {
                    mu_pow_over_fact = mu_pow_over_fact.mul(&mu, wp, rm).div(
                        &BigFloat::from_i128(k as i128, wp),
                        wp,
                        rm,
                    );
                }
                let arg = s.sub(&BigFloat::from_i128(k as i128, wp), wp, rm);
                let zeta = arb_zeta(&arg, wp, rm, cc)?;
                let term = zeta.mul(&mu_pow_over_fact, wp, rm);
                sum = sum.add(&term, wp, rm);
                if k > 2 && negligible(&mu_pow_over_fact, &one, wp + 8) {
                    break;
                }
            }
            Ok(lead.add(&sum, wp, rm))
        }
    }
}

/// Largest `|z|` for which the direct series is used.  For integer `s` the
/// `μ`-expansion is cheap (Bernoulli numbers), so it takes over at `1/2`;
/// for non-integer `s` every term of that expansion costs a `ζ` evaluation,
/// and the direct series stays cheaper until `|z|` is very close to `1`.
fn polylog_series_limit(s_int: Option<i64>, wp: usize) -> BigFloat {
    if s_int.is_some() {
        BigFloat::from_f64(0.5, wp)
    } else {
        BigFloat::from_f64(0.95, wp)
    }
}

/// `Li_s(z)` for `0 < |z| < 1`.
fn polylog_unit_interval(
    s: &BigFloat,
    s_int: Option<i64>,
    z: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let limit = polylog_series_limit(s_int, wp);
    let az = z.abs();
    if !bf_gt(&az, &limit) {
        return Ok(polylog_series(s, s_int, z, wp, rm, cc));
    }
    if z.is_negative() {
        // Li_s(−x) = 2^{1−s} Li_s(x²) − Li_s(x)
        let one = BigFloat::from_i32(1, wp);
        let x2 = az.mul(&az, wp, rm);
        let a = polylog_unit_interval(s, s_int, &x2, wp, rm, cc)?;
        let b = polylog_near_one(s, s_int, &az, wp, rm, cc)?;
        let two = BigFloat::from_i32(2, wp);
        let one_minus_s = one.sub(s, wp, rm);
        let scale = bf_pow(&two, &one_minus_s, wp, rm, cc);
        return Ok(scale.mul(&a, wp, rm).sub(&b, wp, rm));
    }
    polylog_near_one(s, s_int, z, wp, rm, cc)
}

/// Polylogarithm `Li_s(z)` for real `s` and real `−1 ≤ z ≤ 1` (any real
/// `z ≠ 1` when `s` is a non-positive integer, where `Li_s` is rational).
fn arb_polylog(
    s: &BigFloat,
    z: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if s.is_nan() || z.is_nan() || s.is_inf() || z.is_inf() {
        return Err(unevaluable("polylog of special float value"));
    }
    if z.is_zero() {
        return Ok(BigFloat::new(prec));
    }
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let s_int = bf_as_int(s, rm, cc)?;
    if let Some(n) = s_int
        && n <= 0
    {
        let r = polylog_nonpositive((-n) as usize, z, wp, rm)?;
        return Ok(round_to(r, prec, rm));
    }
    let az = z.abs();
    if bf_gt(&az, &one) {
        return Err(unevaluable(
            "polylog outside the unit disc not yet supported in evalf",
        ));
    }
    if az.sub(&one, wp, rm).is_zero() {
        // z = ±1
        if !bf_gt(s, &one) {
            return Err(unevaluable("polylog(s, ±1) diverges for s ≤ 1"));
        }
        let zeta = arb_zeta(s, wp, rm, cc)?;
        if z.is_negative() {
            // Li_s(−1) = −η(s) = −(1 − 2^{1−s}) ζ(s)
            let two = BigFloat::from_i32(2, wp);
            let one_minus_s = one.sub(s, wp, rm);
            let coeff = one.sub(&bf_pow(&two, &one_minus_s, wp, rm, cc), wp, rm);
            return Ok(round_to(coeff.mul(&zeta, wp, rm).neg(), prec, rm));
        }
        return Ok(round_to(zeta, prec, rm));
    }
    let r = polylog_unit_interval(s, s_int, z, wp, rm, cc)?;
    Ok(round_to(r, prec, rm))
}

/// Dirichlet eta `η(s) = (1 − 2^{1−s}) ζ(s)`, with `η(1) = ln 2`.
fn arb_dirichlet_eta(
    s: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if s.is_nan() || s.is_inf_neg() {
        return Err(unevaluable("dirichlet_eta of special float value"));
    }
    if s.is_inf_pos() {
        return Ok(BigFloat::from_i32(1, prec));
    }
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    let one_minus_s = one.sub(s, wp, rm);
    if one_minus_s.is_zero() {
        return Ok(round_to(two.ln(wp, rm, cc), prec, rm));
    }
    let coeff = one.sub(&bf_pow(&two, &one_minus_s, wp, rm, cc), wp, rm);
    let zeta = arb_zeta(s, wp, rm, cc)?;
    Ok(round_to(coeff.mul(&zeta, wp, rm), prec, rm))
}

// ── Airy functions ─────────────────────────────────────────────────────────────

#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum AiryKind {
    Ai,
    Bi,
    AiPrime,
    BiPrime,
}

/// Airy functions for real `x`.
///
/// * `|x|^{3/2} ≤ (3/4)·wp·ln 2`: Maclaurin series (DLMF 9.4.1–9.4.4)
///   `Ai = c₁ f − c₂ g`, `Bi = √3 (c₁ f + c₂ g)` with
///   `f = Σ x^{3k}/∏_{j≤k}(3j−1)(3j)`, `g = Σ x^{3k+1}/∏_{j≤k}(3j)(3j+1)`,
///   `c₁ = 3^{−2/3}/Γ(2/3)`, `c₂ = 3^{−1/3}/Γ(1/3)`; the working precision
///   carries `(4/3)|x|^{3/2} log₂ e` guard bits for the cancellation.
/// * beyond, with `ζ = (2/3)|x|^{3/2}`: the asymptotic expansions
///   DLMF 9.7.5–9.7.12 (`u_k`, `v_k = −(6k+1)/(6k−1) u_k`), whose optimal
///   truncation error is `~e^{−2ζ}`.
fn arb_airy(
    x: &BigFloat,
    kind: AiryKind,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if x.is_nan() || x.is_inf() {
        return Err(unevaluable("Airy function of special float value"));
    }
    let base_wp = prec + 32;
    let x_f = bigfloat_to_f64(x, rm, cc)?;
    let ax_f = x_f.abs();
    let zeta_f = 2.0 / 3.0 * ax_f.powf(1.5);

    if 2.0 * zeta_f > (base_wp as f64) * std::f64::consts::LN_2 {
        return Ok(round_to(
            airy_asymptotic(x, kind, base_wp, rm, cc),
            prec,
            rm,
        ));
    }

    let wp = base_wp + (2.0 * zeta_f * std::f64::consts::LOG2_E).ceil() as usize + 16;
    let mut xw = x.clone();
    let _ = xw.set_precision(wp, rm);
    let one = BigFloat::from_i32(1, wp);
    let three = BigFloat::from_i32(3, wp);
    let third = one.div(&three, wp, rm);
    let two_thirds = third.mul(&BigFloat::from_i32(2, wp), wp, rm);
    let g13 = arb_gamma_real(&third, wp, rm, cc)?;
    let g23 = arb_gamma_real(&two_thirds, wp, rm, cc)?;
    let c1 = one.div(
        &bf_pow(&three, &two_thirds, wp, rm, cc).mul(&g23, wp, rm),
        wp,
        rm,
    );
    let c2 = one.div(
        &bf_pow(&three, &third, wp, rm, cc).mul(&g13, wp, rm),
        wp,
        rm,
    );

    // Series with term recurrences; the derivatives use
    // f'_k = 3k f_k / x, g'_k = (3k+1) g_k / x (k = 0 handled directly).
    let x3 = xw.mul(&xw, wp, rm).mul(&xw, wp, rm);
    let derivative = matches!(kind, AiryKind::AiPrime | AiryKind::BiPrime);
    let (mut f_sum, mut g_sum) = if derivative {
        (BigFloat::new(wp), one.clone()) // f'(0) = 0, g'(0) = 1
    } else {
        (one.clone(), xw.clone()) // f(0) = 1, g_0 = x
    };
    let mut fk = one.clone(); // x^{3k}/∏(3j−1)(3j)
    let mut gk = xw.clone(); // x^{3k+1}/∏(3j)(3j+1)
    let max_terms = (ax_f * 2.0) as usize + wp + 40;
    for k in 1..=max_terms {
        let k3 = 3 * k as i128;
        fk = fk
            .mul(&x3, wp, rm)
            .div(&BigFloat::from_i128((k3 - 1) * k3, wp), wp, rm);
        gk = gk
            .mul(&x3, wp, rm)
            .div(&BigFloat::from_i128(k3 * (k3 + 1), wp), wp, rm);
        let (tf, tg) = if derivative {
            // 3k x^{3k−1}/∏ = (3k/x) f_k ;  (3k+1) x^{3k}/∏ = ((3k+1)/x) g_k
            (
                fk.mul(&BigFloat::from_i128(k3, wp), wp, rm)
                    .div(&xw, wp, rm),
                gk.mul(&BigFloat::from_i128(k3 + 1, wp), wp, rm)
                    .div(&xw, wp, rm),
            )
        } else {
            (fk.clone(), gk.clone())
        };
        f_sum = f_sum.add(&tf, wp, rm);
        g_sum = g_sum.add(&tg, wp, rm);
        if negligible(&tf, &one, wp + 4) && negligible(&tg, &one, wp + 4) {
            break;
        }
    }
    let a = c1.mul(&f_sum, wp, rm);
    let b = c2.mul(&g_sum, wp, rm);
    let r = match kind {
        AiryKind::Ai | AiryKind::AiPrime => a.sub(&b, wp, rm),
        AiryKind::Bi | AiryKind::BiPrime => three.sqrt(wp, rm).mul(&a.add(&b, wp, rm), wp, rm),
    };
    Ok(round_to(r, prec, rm))
}

/// Asymptotic expansions of the Airy functions for large `|x|`
/// (see [`arb_airy`]).  Returns the value at working precision `wp`.
fn airy_asymptotic(
    x: &BigFloat,
    kind: AiryKind,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> BigFloat {
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    let ax = x.abs();
    let sqrt_ax = ax.sqrt(wp, rm);
    let x_quarter = sqrt_ax.sqrt(wp, rm); // |x|^{1/4}
    let zeta = two
        .mul(&ax, wp, rm)
        .mul(&sqrt_ax, wp, rm)
        .div(&BigFloat::from_i32(3, wp), wp, rm);
    let inv_zeta = one.div(&zeta, wp, rm);
    let sqrt_pi = cc.pi(wp, rm).clone().sqrt(wp, rm);
    let derivative = matches!(kind, AiryKind::AiPrime | AiryKind::BiPrime);

    // Coefficients c_k = u_k (or v_k) / ζ^k, u_k = u_{k−1}(6k−5)(6k−3)(6k−1)/(216 k (2k−1)),
    // v_k = −(6k+1)/(6k−1) u_k.  Optimal truncation: stop when |c_k| grows.
    let mut coeffs: Vec<BigFloat> = Vec::new();
    let mut u = one.clone();
    let mut zeta_pow = one.clone();
    let mut prev: Option<BigFloat> = None;
    for k in 0..(wp * 2 + 100) {
        if k > 0 {
            let k_i = k as i128;
            let num = BigFloat::from_i128((6 * k_i - 5) * (6 * k_i - 3) * (6 * k_i - 1), wp);
            let den = BigFloat::from_i128(216 * k_i * (2 * k_i - 1), wp);
            u = u.mul(&num, wp, rm).div(&den, wp, rm);
            zeta_pow = zeta_pow.mul(&inv_zeta, wp, rm);
        }
        let mut c = u.mul(&zeta_pow, wp, rm);
        if derivative {
            let k_i = k as i128;
            let ratio = BigFloat::from_i128(6 * k_i + 1, wp).div(
                &BigFloat::from_i128(6 * k_i - 1, wp),
                wp,
                rm,
            );
            c = c.mul(&ratio, wp, rm).neg();
        }
        if let Some(ref p) = prev
            && bf_gt(&c.abs(), p)
        {
            break;
        }
        let done = negligible(&c, &one, wp);
        prev = Some(c.abs());
        coeffs.push(c);
        if done {
            break;
        }
    }

    if !x.is_negative() {
        // Ai ~ e^{−ζ}/(2√π x^{1/4}) Σ (−1)^k u_k/ζ^k,   Bi ~ e^{ζ}/(√π x^{1/4}) Σ u_k/ζ^k
        // Ai' ~ −x^{1/4} e^{−ζ}/(2√π) Σ (−1)^k v_k/ζ^k, Bi' ~ x^{1/4} e^{ζ}/√π Σ v_k/ζ^k
        let mut sum = BigFloat::new(wp);
        let alternate = matches!(kind, AiryKind::Ai | AiryKind::AiPrime);
        for (k, c) in coeffs.iter().enumerate() {
            if alternate && k % 2 == 1 {
                sum = sum.sub(c, wp, rm);
            } else {
                sum = sum.add(c, wp, rm);
            }
        }
        let amp = match kind {
            AiryKind::Ai => zeta.neg().exp(wp, rm, cc).div(
                &two.mul(&sqrt_pi, wp, rm).mul(&x_quarter, wp, rm),
                wp,
                rm,
            ),
            AiryKind::Bi => zeta
                .exp(wp, rm, cc)
                .div(&sqrt_pi.mul(&x_quarter, wp, rm), wp, rm),
            AiryKind::AiPrime => x_quarter
                .mul(&zeta.neg().exp(wp, rm, cc), wp, rm)
                .div(&two.mul(&sqrt_pi, wp, rm), wp, rm)
                .neg(),
            AiryKind::BiPrime => x_quarter
                .mul(&zeta.exp(wp, rm, cc), wp, rm)
                .div(&sqrt_pi, wp, rm),
        };
        return amp.mul(&sum, wp, rm);
    }

    // x < 0: with θ = ζ + π/4, P = Σ (−1)^k c_{2k}, Q = Σ (−1)^k c_{2k+1}:
    // Ai  ~ (sinθ P − cosθ Q)/(√π x^{1/4}),  Bi  ~ (cosθ P + sinθ Q)/(√π x^{1/4})
    // Ai' ~ −x^{1/4}(cosθ P + sinθ Q)/√π,  Bi' ~  x^{1/4}(sinθ P − cosθ Q)/√π
    let mut p = BigFloat::new(wp);
    let mut q = BigFloat::new(wp);
    for (k, c) in coeffs.iter().enumerate() {
        let sign_neg = (k / 2) % 2 == 1;
        let target = if k % 2 == 0 { &mut p } else { &mut q };
        *target = if sign_neg {
            target.sub(c, wp, rm)
        } else {
            target.add(c, wp, rm)
        };
    }
    let pi = cc.pi(wp, rm).clone();
    let theta = zeta.add(&pi.div(&BigFloat::from_i32(4, wp), wp, rm), wp, rm);
    let s = theta.sin(wp, rm, cc);
    let c = theta.cos(wp, rm, cc);
    let sp_minus_cq = s.mul(&p, wp, rm).sub(&c.mul(&q, wp, rm), wp, rm);
    let cp_plus_sq = c.mul(&p, wp, rm).add(&s.mul(&q, wp, rm), wp, rm);
    match kind {
        AiryKind::Ai => sp_minus_cq.div(&sqrt_pi.mul(&x_quarter, wp, rm), wp, rm),
        AiryKind::Bi => cp_plus_sq.div(&sqrt_pi.mul(&x_quarter, wp, rm), wp, rm),
        AiryKind::AiPrime => x_quarter
            .mul(&cp_plus_sq, wp, rm)
            .div(&sqrt_pi, wp, rm)
            .neg(),
        AiryKind::BiPrime => x_quarter.mul(&sp_minus_cq, wp, rm).div(&sqrt_pi, wp, rm),
    }
}

// ── Elliptic integrals (Carlson symmetric forms) ────────────────────────────────

/// Are `a`, `b`, `c` (and `d`) equal to within `2^{−bits}` relative?
fn all_close(vals: &[&BigFloat], bits: usize, wp: usize, rm: RoundingMode) -> bool {
    let first = vals[0];
    vals[1..]
        .iter()
        .all(|v| negligible(&v.sub(first, wp, rm), first, bits))
}

/// Convergence threshold for the Carlson duplication iterations: the
/// remainder of `R_F`, `R_D`, `R_J` at arguments within relative spread `δ`
/// of their mean is `O(δ²)` (the Taylor expansions have no linear term), so
/// a spread of `2^{−wp/2}` already gives `wp` correct bits — and stopping
/// there avoids stalling on rounding noise in the last bits.
fn carlson_bits(wp: usize) -> usize {
    wp / 2 + 8
}

/// Carlson `R_C(x, y) = R_F(x, y, y)`, `x ≥ 0`, `y > 0`.
///
/// With `w = (y − x)/y`, `R_C = y^{−1/2} arcsin(√w)/√w = y^{−1/2} Σ_k C(2k,k) w^k / (4^k (2k+1))`
/// (the series is analytic through `w = 0`, covering `x > y` as well).
/// The series is used for `|w| < 1/10`: the closed forms
/// `arccos(√(x/y))/√(y−x)` / `arccosh(√(x/y))/√(x−y)` lose relative
/// precision like `1/w` there, which the `R_J` duplication (where
/// `y − x` shrinks quadratically) would amplify into a first-order error.
fn carlson_rc(
    x: &BigFloat,
    y: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> BigFloat {
    let one = BigFloat::from_i32(1, wp);
    let diff = y.sub(x, wp, rm);
    let inv_sqrt_y = one.div(&y.sqrt(wp, rm), wp, rm);
    if diff.is_zero() {
        return inv_sqrt_y;
    }
    let w = diff.div(y, wp, rm);
    let w_f = bigfloat_to_f64(&w, rm, cc).unwrap_or(1.0);
    if w_f.abs() < 0.1 {
        // c_k = c_{k−1} (2k−1)² / ((2k)(2k+1)):  1, 1/6, 3/40, 5/112, …
        let mut sum = one.clone();
        let mut c = one.clone();
        let mut wk = one.clone();
        for k in 1..(wp * 2 + 100) {
            let k_i = k as i128;
            let num = BigFloat::from_i128((2 * k_i - 1) * (2 * k_i - 1), wp);
            let den = BigFloat::from_i128(2 * k_i * (2 * k_i + 1), wp);
            c = c.mul(&num, wp, rm).div(&den, wp, rm);
            wk = wk.mul(&w, wp, rm);
            let term = c.mul(&wk, wp, rm);
            sum = sum.add(&term, wp, rm);
            if negligible(&term, &sum, wp) {
                break;
            }
        }
        return inv_sqrt_y.mul(&sum, wp, rm);
    }
    let ratio = x.div(y, wp, rm).sqrt(wp, rm); // √(x/y)
    if diff.is_positive() {
        // y > x: arccos(√(x/y)) / √(y − x)
        ratio.acos(wp, rm, cc).div(&diff.sqrt(wp, rm), wp, rm)
    } else {
        // x > y: arccosh(√(x/y)) / √(x − y)
        ratio
            .acosh(wp, rm, cc)
            .div(&diff.neg().sqrt(wp, rm), wp, rm)
    }
}

/// Carlson `R_F(x, y, z)` (`x, y, z ≥ 0`, at most one zero) by the
/// duplication theorem iterated until the arguments coincide to working
/// precision, where `R_F(μ, μ, μ) = μ^{−1/2}`.
fn carlson_rf(
    x: &BigFloat,
    y: &BigFloat,
    z: &BigFloat,
    wp: usize,
    rm: RoundingMode,
) -> Result<BigFloat, SymplexError> {
    let four = BigFloat::from_i32(4, wp);
    let (mut x, mut y, mut z) = (x.clone(), y.clone(), z.clone());
    for _ in 0..(wp * 2 + 200) {
        if all_close(&[&x, &y, &z], carlson_bits(wp), wp, rm) {
            // R_F(μ, μ, μ) = μ^{−1/2}, μ = (x + y + z)/3
            let one = BigFloat::from_i32(1, wp);
            let mu = x
                .add(&y, wp, rm)
                .add(&z, wp, rm)
                .div(&BigFloat::from_i32(3, wp), wp, rm);
            return Ok(one.div(&mu.sqrt(wp, rm), wp, rm));
        }
        let (sx, sy, sz) = (x.sqrt(wp, rm), y.sqrt(wp, rm), z.sqrt(wp, rm));
        let lambda =
            sx.mul(&sy, wp, rm)
                .add(&sy.mul(&sz, wp, rm), wp, rm)
                .add(&sz.mul(&sx, wp, rm), wp, rm);
        x = x.add(&lambda, wp, rm).div(&four, wp, rm);
        y = y.add(&lambda, wp, rm).div(&four, wp, rm);
        z = z.add(&lambda, wp, rm).div(&four, wp, rm);
    }
    Err(SymplexError::ComputationFailed {
        operation: "evalf",
        reason: "Carlson R_F did not converge".into(),
    })
}

/// Carlson `R_D(x, y, z) = R_J(x, y, z, z)` (`z > 0`) by duplication:
/// `R_D(x,y,z) = R_D(x',y',z')/4 + 3/(√z (z + λ))`, `R_D(μ,μ,μ) = μ^{−3/2}`.
fn carlson_rd(
    x: &BigFloat,
    y: &BigFloat,
    z: &BigFloat,
    wp: usize,
    rm: RoundingMode,
) -> Result<BigFloat, SymplexError> {
    let one = BigFloat::from_i32(1, wp);
    let three = BigFloat::from_i32(3, wp);
    let four = BigFloat::from_i32(4, wp);
    let (mut x, mut y, mut z) = (x.clone(), y.clone(), z.clone());
    let mut acc = BigFloat::new(wp);
    let mut scale = one.clone(); // 4^{−k}
    for _ in 0..(wp * 2 + 200) {
        if all_close(&[&x, &y, &z], carlson_bits(wp), wp, rm) {
            // R_D(μ, μ, μ) = μ^{−3/2}, μ = (x + y + 3z)/5
            let mu = x.add(&y, wp, rm).add(&three.mul(&z, wp, rm), wp, rm).div(
                &BigFloat::from_i32(5, wp),
                wp,
                rm,
            );
            let tail = one.div(&mu.mul(&mu.sqrt(wp, rm), wp, rm), wp, rm);
            return Ok(acc.add(&scale.mul(&tail, wp, rm), wp, rm));
        }
        let (sx, sy, sz) = (x.sqrt(wp, rm), y.sqrt(wp, rm), z.sqrt(wp, rm));
        let lambda =
            sx.mul(&sy, wp, rm)
                .add(&sy.mul(&sz, wp, rm), wp, rm)
                .add(&sz.mul(&sx, wp, rm), wp, rm);
        let term = three.div(&sz.mul(&z.add(&lambda, wp, rm), wp, rm), wp, rm);
        acc = acc.add(&scale.mul(&term, wp, rm), wp, rm);
        scale = scale.div(&four, wp, rm);
        x = x.add(&lambda, wp, rm).div(&four, wp, rm);
        y = y.add(&lambda, wp, rm).div(&four, wp, rm);
        z = z.add(&lambda, wp, rm).div(&four, wp, rm);
    }
    Err(SymplexError::ComputationFailed {
        operation: "evalf",
        reason: "Carlson R_D did not converge".into(),
    })
}

/// Carlson `R_J(x, y, z, p)` for `p > 0` by duplication:
/// `R_J = R_J(x',y',z',p')/4 + 3 R_C(α², β²)` with
/// `α = p(√x + √y + √z) + √x√y√z`, `β = √p (p + λ)`, `R_J(μ,μ,μ,μ) = μ^{−3/2}`.
fn carlson_rj(
    x: &BigFloat,
    y: &BigFloat,
    z: &BigFloat,
    p: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if !p.is_positive() {
        return Err(unevaluable(
            "elliptic_pi: the characteristic n must satisfy n < 1 for a real value",
        ));
    }
    let one = BigFloat::from_i32(1, wp);
    let three = BigFloat::from_i32(3, wp);
    let four = BigFloat::from_i32(4, wp);
    let (mut x, mut y, mut z, mut p) = (x.clone(), y.clone(), z.clone(), p.clone());
    let mut acc = BigFloat::new(wp);
    let mut scale = one.clone();
    for _ in 0..(wp * 2 + 200) {
        if all_close(&[&x, &y, &z, &p], carlson_bits(wp), wp, rm) {
            // R_J(μ, μ, μ, μ) = μ^{−3/2}, μ = (x + y + z + 2p)/5
            let two_p = p.mul(&BigFloat::from_i32(2, wp), wp, rm);
            let mu = x.add(&y, wp, rm).add(&z, wp, rm).add(&two_p, wp, rm).div(
                &BigFloat::from_i32(5, wp),
                wp,
                rm,
            );
            let tail = one.div(&mu.mul(&mu.sqrt(wp, rm), wp, rm), wp, rm);
            return Ok(acc.add(&scale.mul(&tail, wp, rm), wp, rm));
        }
        let (sx, sy, sz, sp) = (
            x.sqrt(wp, rm),
            y.sqrt(wp, rm),
            z.sqrt(wp, rm),
            p.sqrt(wp, rm),
        );
        let lambda =
            sx.mul(&sy, wp, rm)
                .add(&sy.mul(&sz, wp, rm), wp, rm)
                .add(&sz.mul(&sx, wp, rm), wp, rm);
        let alpha = p.mul(&sx.add(&sy, wp, rm).add(&sz, wp, rm), wp, rm).add(
            &sx.mul(&sy, wp, rm).mul(&sz, wp, rm),
            wp,
            rm,
        );
        let beta = sp.mul(&p.add(&lambda, wp, rm), wp, rm);
        let rc = carlson_rc(
            &alpha.mul(&alpha, wp, rm),
            &beta.mul(&beta, wp, rm),
            wp,
            rm,
            cc,
        );
        acc = acc.add(&scale.mul(&three.mul(&rc, wp, rm), wp, rm), wp, rm);
        scale = scale.div(&four, wp, rm);
        x = x.add(&lambda, wp, rm).div(&four, wp, rm);
        y = y.add(&lambda, wp, rm).div(&four, wp, rm);
        z = z.add(&lambda, wp, rm).div(&four, wp, rm);
        p = p.add(&lambda, wp, rm).div(&four, wp, rm);
    }
    Err(SymplexError::ComputationFailed {
        operation: "evalf",
        reason: "Carlson R_J did not converge".into(),
    })
}

/// `1 − m` with a domain check (`m < 1` required for a real `K`, `E`, `Π`).
fn one_minus_m(
    m: &BigFloat,
    wp: usize,
    rm: RoundingMode,
    what: &str,
) -> Result<BigFloat, SymplexError> {
    let one = BigFloat::from_i32(1, wp);
    let v = one.sub(m, wp, rm);
    if v.is_negative() {
        return Err(unevaluable(format!(
            "{what}: parameter m must satisfy m ≤ 1 for a real value"
        )));
    }
    Ok(v)
}

/// Complete elliptic integral of the first kind `K(m) = R_F(0, 1−m, 1)`.
fn arb_elliptic_k(
    m: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if m.is_nan() || m.is_inf() {
        return Err(unevaluable("elliptic_k of special float value"));
    }
    let wp = prec + 32;
    let y = one_minus_m(m, wp, rm, "elliptic_k")?;
    if y.is_zero() {
        return Err(unevaluable("elliptic_k(1) is infinite"));
    }
    let _ = cc;
    let r = carlson_rf(&BigFloat::new(wp), &y, &BigFloat::from_i32(1, wp), wp, rm)?;
    Ok(round_to(r, prec, rm))
}

/// Complete elliptic integral of the second kind
/// `E(m) = R_F(0, 1−m, 1) − (m/3) R_D(0, 1−m, 1)`.
fn arb_elliptic_e(
    m: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if m.is_nan() || m.is_inf() {
        return Err(unevaluable("elliptic_e of special float value"));
    }
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let y = one_minus_m(m, wp, rm, "elliptic_e")?;
    if y.is_zero() {
        return Ok(BigFloat::from_i32(1, prec));
    }
    let _ = cc;
    let zero = BigFloat::new(wp);
    let rf = carlson_rf(&zero, &y, &one, wp, rm)?;
    let rd = carlson_rd(&zero, &y, &one, wp, rm)?;
    let three = BigFloat::from_i32(3, wp);
    let r = rf.sub(&m.mul(&rd, wp, rm).div(&three, wp, rm), wp, rm);
    Ok(round_to(r, prec, rm))
}

/// Incomplete elliptic integral of the first kind
/// `F(φ | m) = sinφ · R_F(cos²φ, 1 − m sin²φ, 1)` for `|φ| ≤ π/2`, extended
/// by `F(φ + kπ | m) = F(φ | m) + 2k K(m)`.
fn arb_elliptic_f(
    phi: &BigFloat,
    m: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if phi.is_nan() || phi.is_inf() || m.is_nan() || m.is_inf() {
        return Err(unevaluable("elliptic_f of special float value"));
    }
    if phi.is_zero() {
        return Ok(BigFloat::new(prec));
    }
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let pi = cc.pi(wp, rm).clone();
    // Reduce φ to [−π/2, π/2]: φ = φ₀ + kπ.
    let k_bf = phi
        .div(&pi, wp, rm)
        .add(&BigFloat::from_f64(0.5, wp), wp, rm)
        .floor();
    let phi0 = phi.sub(&k_bf.mul(&pi, wp, rm), wp, rm);
    let s = phi0.sin(wp, rm, cc);
    let c = phi0.cos(wp, rm, cc);
    let s2 = s.mul(&s, wp, rm);
    let c2 = c.mul(&c, wp, rm);
    let y = one.sub(&m.mul(&s2, wp, rm), wp, rm);
    if !y.is_positive() {
        return Err(unevaluable(
            "elliptic_f: 1 − m sin²φ must be positive for a real value",
        ));
    }
    let mut r = s.mul(&carlson_rf(&c2, &y, &one, wp, rm)?, wp, rm);
    if !k_bf.is_zero() {
        let k = arb_elliptic_k(m, wp, rm, cc)?;
        let two = BigFloat::from_i32(2, wp);
        r = r.add(&two.mul(&k_bf, wp, rm).mul(&k, wp, rm), wp, rm);
    }
    Ok(round_to(r, prec, rm))
}

/// Complete elliptic integral of the third kind
/// `Π(n | m) = R_F(0, 1−m, 1) + (n/3) R_J(0, 1−m, 1, 1−n)` for `n < 1`, `m < 1`.
fn arb_elliptic_pi(
    n: &BigFloat,
    m: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    if n.is_nan() || n.is_inf() || m.is_nan() || m.is_inf() {
        return Err(unevaluable("elliptic_pi of special float value"));
    }
    let wp = prec + 32;
    let one = BigFloat::from_i32(1, wp);
    let zero = BigFloat::new(wp);
    let y = one_minus_m(m, wp, rm, "elliptic_pi")?;
    if y.is_zero() {
        return Err(unevaluable("elliptic_pi(n, 1) is infinite"));
    }
    let p = one.sub(n, wp, rm);
    let rf = carlson_rf(&zero, &y, &one, wp, rm)?;
    if n.is_zero() {
        return Ok(round_to(rf, prec, rm));
    }
    let rj = carlson_rj(&zero, &y, &one, &p, wp, rm, cc)?;
    let three = BigFloat::from_i32(3, wp);
    let r = rf.add(&n.mul(&rj, wp, rm).div(&three, wp, rm), wp, rm);
    Ok(round_to(r, prec, rm))
}

// ── Orthogonal polynomials with parameters ───────────────────────────────────────

/// The degree of an orthogonal polynomial as a non-negative integer.
fn poly_degree(
    n: &BigFloat,
    what: &str,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<u64, SymplexError> {
    match bf_as_int(n, rm, cc)? {
        Some(d) if (0..=10_000_000).contains(&d) => Ok(d as u64),
        _ => Err(unevaluable(format!(
            "{what}: degree must be a non-negative integer"
        ))),
    }
}

/// Gegenbauer `C_n^{(a)}(x)`: `(k+1) C_{k+1} = 2(k+a) x C_k − (k+2a−1) C_{k−1}`.
fn arb_gegenbauer(
    n: &BigFloat,
    a: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let deg = poly_degree(n, "gegenbauer", rm, cc)?;
    let wp = prec + 32 + (deg as f64).log2().ceil().max(0.0) as usize;
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    if deg == 0 {
        return Ok(BigFloat::from_i32(1, prec));
    }
    let mut prev = one.clone();
    let mut curr = two.mul(a, wp, rm).mul(x, wp, rm);
    for k in 1..deg {
        let k_bf = BigFloat::from_i128(k as i128, wp);
        let t1 = two
            .mul(&k_bf.add(a, wp, rm), wp, rm)
            .mul(x, wp, rm)
            .mul(&curr, wp, rm);
        let coeff = k_bf.add(&two.mul(a, wp, rm), wp, rm).sub(&one, wp, rm);
        let t2 = coeff.mul(&prev, wp, rm);
        let next = t1
            .sub(&t2, wp, rm)
            .div(&BigFloat::from_i128(k as i128 + 1, wp), wp, rm);
        prev = curr;
        curr = next;
    }
    Ok(round_to(curr, prec, rm))
}

/// Jacobi `P_n^{(a,b)}(x)` via the explicit sum
/// `Σ_s C(n+a, n−s) C(n+b, s) ((x−1)/2)^s ((x+1)/2)^{n−s}`.
fn arb_jacobi(
    n: &BigFloat,
    a: &BigFloat,
    b: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let deg = poly_degree(n, "jacobi", rm, cc)? as usize;
    if deg > 100_000 {
        return Err(unevaluable("jacobi: degree too large for evalf"));
    }
    let wp = prec + 32 + 2 * deg;
    let one = BigFloat::from_i32(1, wp);
    let two = BigFloat::from_i32(2, wp);
    if deg == 0 {
        return Ok(BigFloat::from_i32(1, prec));
    }
    let lo = x.sub(&one, wp, rm).div(&two, wp, rm);
    let hi = x.add(&one, wp, rm).div(&two, wp, rm);
    // binomial-like coefficients ∏_{j=lo}^{hi}(p + j) / len!
    let shifted = |p: &BigFloat, lo_j: usize, hi_j: usize| -> BigFloat {
        let mut acc = one.clone();
        for j in lo_j..=hi_j {
            acc = acc.mul(&p.add(&BigFloat::from_i128(j as i128, wp), wp, rm), wp, rm);
        }
        acc
    };
    let mut sum = BigFloat::new(wp);
    for s in 0..=deg {
        let pa = shifted(a, s + 1, deg);
        let pb = shifted(b, deg - s + 1, deg);
        let denom = ratio_to_bigfloat(&(factorial_bigint(deg - s) * factorial_bigint(s)), wp, rm);
        let term = pa
            .mul(&pb, wp, rm)
            .div(&denom, wp, rm)
            .mul(&lo.powi(s, wp, rm), wp, rm)
            .mul(&hi.powi(deg - s, wp, rm), wp, rm);
        sum = sum.add(&term, wp, rm);
    }
    Ok(round_to(sum, prec, rm))
}

/// `n!` as a rational.
fn factorial_bigint(n: usize) -> Ratio<BigInt> {
    let mut f = BigInt::from(1);
    for i in 2..=n {
        f *= BigInt::from(i as u64);
    }
    Ratio::from_integer(f)
}

/// Associated Legendre `P_n^m(x)` (Condon–Shortley phase) for integer
/// `n ≥ 0`, integer `m`, real `x` (`|x| ≤ 1` when `m` is odd).
fn arb_assoc_legendre(
    n: &BigFloat,
    m: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let deg = poly_degree(n, "assoc_legendre", rm, cc)? as i64;
    let order = bf_as_int(m, rm, cc)?
        .ok_or_else(|| unevaluable("assoc_legendre: order must be an integer"))?;
    if order.abs() > deg {
        return Ok(BigFloat::new(prec));
    }
    let wp = prec + 32 + 2 * (deg as f64).log2().ceil().max(0.0) as usize;
    let one = BigFloat::from_i32(1, wp);
    let mu = order.unsigned_abs() as usize;
    let x2 = x.mul(x, wp, rm);
    let one_minus_x2 = one.sub(&x2, wp, rm);
    // (1 − x²)^{μ/2}
    let root_pow = if mu.is_multiple_of(2) {
        one_minus_x2.powi(mu / 2, wp, rm)
    } else {
        if one_minus_x2.is_negative() {
            return Err(unevaluable(
                "assoc_legendre with odd order is complex for |x| > 1",
            ));
        }
        one_minus_x2
            .sqrt(wp, rm)
            .mul(&one_minus_x2.powi(mu / 2, wp, rm), wp, rm)
    };
    // P_μ^μ = (−1)^μ (2μ−1)!! (1−x²)^{μ/2}
    let mut dfact = BigInt::from(1);
    for j in 1..=mu {
        dfact *= BigInt::from((2 * j - 1) as u64);
    }
    if mu % 2 == 1 {
        dfact = -dfact;
    }
    let mut prev = ratio_to_bigfloat(&Ratio::from_integer(dfact), wp, rm).mul(&root_pow, wp, rm);
    let deg_u = deg as usize;
    let mut curr = if deg_u > mu {
        BigFloat::from_i128(2 * mu as i128 + 1, wp)
            .mul(x, wp, rm)
            .mul(&prev, wp, rm)
    } else {
        prev.clone()
    };
    for k in (mu + 1)..deg_u {
        let t1 = BigFloat::from_i128(2 * k as i128 + 1, wp)
            .mul(x, wp, rm)
            .mul(&curr, wp, rm);
        let t2 = BigFloat::from_i128((k + mu) as i128, wp).mul(&prev, wp, rm);
        let next = t1
            .sub(&t2, wp, rm)
            .div(&BigFloat::from_i128((k - mu + 1) as i128, wp), wp, rm);
        prev = curr;
        curr = next;
    }
    if order < 0 {
        // P_n^{−μ} = (−1)^μ (n−μ)!/(n+μ)! P_n^μ
        let mut ratio = factorial_bigint(deg_u - mu) / factorial_bigint(deg_u + mu);
        if mu % 2 == 1 {
            ratio = -ratio;
        }
        curr = curr.mul(&ratio_to_bigfloat(&ratio, wp, rm), wp, rm);
    }
    Ok(round_to(curr, prec, rm))
}

/// Generalised Laguerre `L_n^{(a)}(x)`:
/// `(k+1) L_{k+1} = (2k+1+a−x) L_k − (k+a) L_{k−1}`.
fn arb_assoc_laguerre(
    n: &BigFloat,
    a: &BigFloat,
    x: &BigFloat,
    prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
    let deg = poly_degree(n, "assoc_laguerre", rm, cc)?;
    let wp = prec + 32 + (deg as f64).log2().ceil().max(0.0) as usize;
    let one = BigFloat::from_i32(1, wp);
    if deg == 0 {
        return Ok(BigFloat::from_i32(1, prec));
    }
    let mut prev = one.clone();
    let mut curr = one.add(a, wp, rm).sub(x, wp, rm);
    for k in 1..deg {
        let k_bf = BigFloat::from_i128(k as i128, wp);
        let coeff = BigFloat::from_i128(2 * k as i128 + 1, wp)
            .add(a, wp, rm)
            .sub(x, wp, rm);
        let t1 = coeff.mul(&curr, wp, rm);
        let t2 = k_bf.add(a, wp, rm).mul(&prev, wp, rm);
        let next = t1
            .sub(&t2, wp, rm)
            .div(&BigFloat::from_i128(k as i128 + 1, wp), wp, rm);
        prev = curr;
        curr = next;
    }
    Ok(round_to(curr, prec, rm))
}

fn get_cached(cache: &FxHashMap<ExprId, Complex>, id: ExprId) -> Result<&Complex, SymplexError> {
    cache.get(&id).ok_or_else(|| SymplexError::Unevaluable {
        reason: format!("sub-expression {id:?} not in cache (likely contains free symbols)"),
    })
}

/// Convert a `Ratio<BigInt>` to a `BigFloat` with the given precision in bits.
fn ratio_to_bigfloat(r: &Ratio<BigInt>, prec: usize, rm: RoundingMode) -> BigFloat {
    if r.is_zero() {
        return BigFloat::from_i32(0, prec);
    }

    let numer = r.numer();
    let denom = r.denom();

    let n_bf = bigint_to_bigfloat(numer, prec);

    if denom == &BigInt::from(1) {
        return n_bf;
    }

    let d_bf = bigint_to_bigfloat(denom, prec);
    n_bf.div(&d_bf, prec, rm)
}

/// Convert a `BigInt` to a `BigFloat`, exactly up to the final rounding
/// to `prec` bits.
///
/// Values that fit in `i128` are converted directly.  Larger integers are
/// accumulated limb-by-limb (`acc = acc·2⁶⁴ + limb`) at a working
/// precision wide enough to hold every bit, then rounded once.
fn bigint_to_bigfloat(n: &BigInt, prec: usize) -> BigFloat {
    if let Ok(v) = <BigInt as TryInto<i128>>::try_into(n.clone()) {
        return BigFloat::from_i128(v, prec);
    }
    let (sign, limbs) = n.to_u64_digits();
    let wp = (limbs.len() * 64 + 64).max(prec);
    let rm = RoundingMode::ToEven;
    let base = BigFloat::from_u64(1u64 << 32, wp).powi(2, wp, rm); // 2^64
    let mut acc = BigFloat::new(wp);
    for &limb in limbs.iter().rev() {
        acc = acc
            .mul(&base, wp, rm)
            .add(&BigFloat::from_u64(limb, wp), wp, rm);
    }
    if sign == num_bigint::Sign::Minus {
        acc = acc.neg();
    }
    let _ = acc.set_precision(prec, rm);
    acc
}

/// Try to extract a small integer exponent (fits in i32) from an ExprId.
fn try_as_small_int(arena: &Arena, id: ExprId) -> Option<i32> {
    let r = arena.as_num(id)?;
    if !r.is_integer() {
        return None;
    }
    let n = r.to_integer();
    let val: i64 = n.try_into().ok()?;
    if val >= i32::MIN as i64 && val <= i32::MAX as i64 {
        Some(val as i32)
    } else {
        None
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Decimal / complex formatting
// ═══════════════════════════════════════════════════════════════════════════

/// Determine whether `part` is negligibly small compared to `other`.
///
/// Returns `true` if `part` is zero, or if its binary exponent is more
/// than `digits * log2(10)` bits below `other`'s exponent (meaning it
/// is below the requested output precision).
fn is_negligible_part(part: &BigFloat, other: &BigFloat, digits: u32) -> bool {
    if part.is_zero() {
        return true;
    }
    if other.is_zero() {
        return false;
    }
    match (part.exponent(), other.exponent()) {
        (Some(p_exp), Some(o_exp)) => {
            let bit_threshold = (digits as i64) * 34 / 10 + 4;
            (o_exp as i64 - p_exp as i64) > bit_threshold
        }
        _ => false,
    }
}

/// Format a complex result as a string.
///
/// If the imaginary part is negligible, formats as a real number.
/// If the real part is negligible, formats as a pure imaginary number.
/// Otherwise formats as `a + b*i` or `a - b*i`.
fn format_complex(
    z: &Complex,
    digits: u32,
    _prec: usize,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<String, SymplexError> {
    let re = &z.0;
    let im = &z.1;

    let re_negligible = is_negligible_part(re, im, digits);
    let im_negligible = is_negligible_part(im, re, digits);

    if re_negligible && im_negligible {
        return Ok("0".to_string());
    }

    if im_negligible {
        return format_decimal(re, digits, rm, cc);
    }

    if re_negligible {
        // Pure imaginary.
        let im_str = format_decimal(im, digits, rm, cc)?;
        return if im_str == "1" {
            Ok("i".to_string())
        } else if im_str == "-1" {
            Ok("-i".to_string())
        } else {
            Ok(format!("{im_str}*i"))
        };
    }

    // Both parts present.
    let re_str = format_decimal(re, digits, rm, cc)?;
    if im.is_negative() {
        let im_abs_str = format_decimal(&im.abs(), digits, rm, cc)?;
        if im_abs_str == "1" {
            Ok(format!("{re_str} - i"))
        } else {
            Ok(format!("{re_str} - {im_abs_str}*i"))
        }
    } else {
        let im_str = format_decimal(im, digits, rm, cc)?;
        if im_str == "1" {
            Ok(format!("{re_str} + i"))
        } else {
            Ok(format!("{re_str} + {im_str}*i"))
        }
    }
}

/// Round a decimal digit string `0.d1d2d3…` to `digits` significant digits.
///
/// Rounds to nearest, ties to even (the digits beyond the requested
/// count are exact decimal digits of the binary value, so an exact tie is
/// a genuine `…5000…`).  Carries propagate leftwards; a carry out of the
/// leading digit (`9.9996 → 10.00`) is reported as an exponent increment
/// so the caller can re-place the decimal point.
///
/// Returns `(rounded_digits, exponent_increment)`; when `mantissa` has no
/// more than `digits` digits it is returned unchanged.
fn round_decimal_digits(mantissa: &[u8], digits: usize) -> (Vec<u8>, i64) {
    if digits == 0 || mantissa.len() <= digits {
        return (mantissa.to_vec(), 0);
    }
    let mut kept: Vec<u8> = mantissa[..digits].to_vec();
    let next = mantissa[digits];
    let rest_nonzero = mantissa[digits + 1..].iter().any(|&d| d != 0);
    let round_up = match next.cmp(&5) {
        std::cmp::Ordering::Greater => true,
        std::cmp::Ordering::Less => false,
        // Exactly half: ties to even.
        std::cmp::Ordering::Equal => rest_nonzero || kept[digits - 1] % 2 == 1,
    };
    if !round_up {
        return (kept, 0);
    }
    // Propagate the carry from the last kept digit leftwards.
    for d in kept.iter_mut().rev() {
        if *d == 9 {
            *d = 0;
        } else {
            *d += 1;
            return (kept, 0);
        }
    }
    // Carried out of the leading digit: 0.999… → 1.000… × 10.
    kept[0] = 1;
    (kept, 1)
}

/// Format a `BigFloat` as a decimal string with `digits` significant digits.
///
/// Uses astro-float's `convert_to_radix` for reliable base-10 conversion,
/// rounds to `digits` significant digits, then formats the result with
/// proper decimal point placement.
fn format_decimal(
    val: &BigFloat,
    digits: u32,
    rm: RoundingMode,
    cc: &mut Consts,
) -> Result<String, SymplexError> {
    if val.is_zero() {
        return Ok("0".to_string());
    }

    if val.is_inf_pos() {
        return Ok("oo".to_string());
    }
    if val.is_inf_neg() {
        return Ok("-oo".to_string());
    }

    // Use the library's decimal conversion.
    let (sign, mantissa, exponent) = val
        .convert_to_radix(Radix::Dec, rm, cc)
        .map_err(|e| SymplexError::NotImplemented(format!("decimal conversion failed: {e:?}")))?;

    let sign_str = if sign == Sign::Neg { "-" } else { "" };

    // mantissa is Vec<u8> of decimal digits [d1, d2, d3, ...]
    // representing the value 0.d1d2d3... × 10^exponent.
    // So the actual value is d1.d2d3... × 10^(exponent-1).

    // Round (not truncate) to the requested digit count.
    let (rounded, exp_carry) = round_decimal_digits(&mantissa, digits as usize);
    let n = rounded.len();
    if n == 0 {
        return Ok("0".to_string());
    }

    let digit_chars: Vec<char> = rounded.iter().map(|&d| (b'0' + d) as char).collect();

    // The "adjusted exponent" is exponent-1 (shifting from 0.ddd to d.ddd notation).
    let adj_exp = exponent as i64 - 1 + exp_carry;

    // Decide between plain decimal and scientific notation.
    if adj_exp >= 0 && (adj_exp as usize) < n {
        // The decimal point falls within the digit string.
        // E.g., digits=[3,1,4,1,5], adj_exp=0 → "3.1415"
        let dot_pos = adj_exp as usize + 1;
        let mut s = String::with_capacity(n + 2);
        s.push_str(sign_str);
        for (i, &ch) in digit_chars.iter().enumerate() {
            if i == dot_pos && dot_pos < n {
                s.push('.');
            }
            s.push(ch);
        }
        // Trim trailing zeros after the decimal point.
        if s.contains('.') {
            let trimmed = s.trim_end_matches('0').trim_end_matches('.');
            Ok(trimmed.to_string())
        } else {
            Ok(s)
        }
    } else if adj_exp >= 0 && (adj_exp as usize) < n + 6 {
        // Small integer — pad with trailing zeros.
        let mut s = String::with_capacity(adj_exp as usize + 2);
        s.push_str(sign_str);
        for &ch in &digit_chars {
            s.push(ch);
        }
        let needed_zeros = (adj_exp as usize + 1).saturating_sub(digit_chars.len());
        for _ in 0..needed_zeros {
            s.push('0');
        }
        Ok(s)
    } else if (-4..0).contains(&adj_exp) {
        // Small fraction like 0.00123.
        let mut s = String::with_capacity(n + (-adj_exp) as usize + 3);
        s.push_str(sign_str);
        s.push_str("0.");
        let leading_zeros = (-adj_exp - 1) as usize;
        for _ in 0..leading_zeros {
            s.push('0');
        }
        for &ch in &digit_chars {
            s.push(ch);
        }
        // Trim trailing zeros.
        let trimmed = s.trim_end_matches('0');
        Ok(trimmed.to_string())
    } else {
        // Scientific notation: d1.d2d3...dn × 10^adj_exp.
        let mut s = String::with_capacity(n + 10);
        s.push_str(sign_str);
        s.push(digit_chars[0]);
        if n > 1 {
            s.push('.');
            for &ch in &digit_chars[1..] {
                s.push(ch);
            }
            // Trim trailing zeros in the fractional part.
            while s.ends_with('0') {
                s.pop();
            }
            if s.ends_with('.') {
                s.pop();
            }
        }
        s.push_str(&format!("e{adj_exp}"));
        Ok(s)
    }
}

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

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

    /// Helper: evaluate and assert the result starts with expected prefix.
    fn assert_evalf_starts_with(arena: &Arena, expr: ExprId, digits: u32, prefix: &str) {
        let result = evalf(arena, expr, digits).unwrap();
        assert!(
            result.starts_with(prefix),
            "expected to start with '{prefix}', got: '{result}'"
        );
    }

    /// Helper: evaluate and assert exact string match.
    fn assert_evalf_eq(arena: &Arena, expr: ExprId, digits: u32, expected: &str) {
        let result = evalf(arena, expr, digits).unwrap();
        assert_eq!(result, expected, "evalf mismatch");
    }

    // ── Integers ────────────────────────────────────────────────────

    #[test]
    fn integer_zero() {
        let a = Arena::new();
        assert_evalf_eq(&a, a.zero, 10, "0");
    }

    #[test]
    fn integer_one() {
        let a = Arena::new();
        assert_evalf_eq(&a, a.one, 10, "1");
    }

    #[test]
    fn integer_negative_one() {
        let a = Arena::new();
        assert_evalf_eq(&a, a.neg_one, 10, "-1");
    }

    #[test]
    fn integer_five() {
        let mut a = Arena::new();
        let five = a.int(5);
        assert_evalf_eq(&a, five, 10, "5");
    }

    #[test]
    fn integer_large() {
        let mut a = Arena::new();
        let n = a.int(123456);
        assert_evalf_eq(&a, n, 10, "123456");
    }

    // ── Rationals ───────────────────────────────────────────────────

    #[test]
    fn rational_half() {
        let mut a = Arena::new();
        let half = a.rational(1, 2);
        assert_evalf_eq(&a, half, 10, "0.5");
    }

    #[test]
    fn rational_one_third() {
        let mut a = Arena::new();
        let third = a.rational(1, 3);
        assert_evalf_starts_with(&a, third, 10, "0.333333333");
    }

    #[test]
    fn rational_negative() {
        let mut a = Arena::new();
        let r = a.rational(-3, 4);
        assert_evalf_starts_with(&a, r, 10, "-0.75");
    }

    // ── Constants ───────────────────────────────────────────────────

    #[test]
    fn constant_pi_15_digits() {
        let a = Arena::new();
        assert_evalf_starts_with(&a, a.pi, 15, "3.14159265358979");
    }

    #[test]
    fn constant_pi_50_digits() {
        let a = Arena::new();
        let result = evalf(&a, a.pi, 50).unwrap();
        assert!(
            result.starts_with("3.1415926535897932384626433832795"),
            "pi to 50 digits: {result}"
        );
    }

    #[test]
    fn constant_e_15_digits() {
        let a = Arena::new();
        assert_evalf_starts_with(&a, a.e_const, 15, "2.71828182845905");
    }

    // ── Arithmetic ──────────────────────────────────────────────────

    #[test]
    fn add_integers() {
        let mut a = Arena::new();
        let two = a.int(2);
        let three = a.int(3);
        let expr = a.add(&[two, three]);
        assert_evalf_eq(&a, expr, 10, "5");
    }

    #[test]
    fn mul_integers() {
        let mut a = Arena::new();
        let two = a.int(2);
        let three = a.int(3);
        let expr = a.mul(&[two, three]);
        assert_evalf_eq(&a, expr, 10, "6");
    }

    #[test]
    fn pow_integer() {
        let mut a = Arena::new();
        let two = a.int(2);
        let ten = a.int(10);
        let expr = a.pow(two, ten);
        assert_evalf_eq(&a, expr, 10, "1024");
    }

    #[test]
    fn pow_negative_exponent() {
        let mut a = Arena::new();
        let two = a.int(2);
        let neg_one = a.int(-1);
        let expr = a.pow(two, neg_one);
        assert_evalf_eq(&a, expr, 10, "0.5");
    }

    #[test]
    fn pi_plus_one() {
        let mut a = Arena::new();
        let pi = a.pi;
        let one = a.one;
        let expr = a.add(&[pi, one]);
        assert_evalf_starts_with(&a, expr, 15, "4.14159265358979");
    }

    // ── Transcendental functions ─────────────────────────────────────

    #[test]
    fn sin_zero() {
        let mut a = Arena::new();
        let zero = a.zero;
        let expr = a.sin(zero);
        assert_evalf_eq(&a, expr, 10, "0");
    }

    #[test]
    fn cos_zero() {
        let mut a = Arena::new();
        let zero = a.zero;
        let expr = a.cos(zero);
        assert_evalf_eq(&a, expr, 10, "1");
    }

    #[test]
    fn exp_zero() {
        let mut a = Arena::new();
        let zero = a.zero;
        let expr = a.exp(zero);
        assert_evalf_eq(&a, expr, 10, "1");
    }

    #[test]
    fn exp_one() {
        let mut a = Arena::new();
        let one = a.one;
        let expr = a.exp(one);
        assert_evalf_starts_with(&a, expr, 15, "2.71828182845905");
    }

    #[test]
    fn ln_one() {
        let mut a = Arena::new();
        let one = a.one;
        let expr = a.ln(one);
        assert_evalf_eq(&a, expr, 10, "0");
    }

    #[test]
    fn ln_e() {
        let mut a = Arena::new();
        let e = a.e_const;
        let expr = a.ln(e);
        // ln(e) should be 1, or very close.
        let result = evalf(&a, expr, 15).unwrap();
        let val: f64 = result.parse().unwrap_or(999.0);
        assert!(
            (val - 1.0).abs() < 1e-10,
            "ln(e) should be ~1, got: {result}"
        );
    }

    #[test]
    fn sqrt_four() {
        let mut a = Arena::new();
        let four = a.int(4);
        let expr = a.sqrt(four);
        assert_evalf_eq(&a, expr, 10, "2");
    }

    #[test]
    fn sqrt_two_30_digits() {
        let mut a = Arena::new();
        let two = a.int(2);
        let expr = a.sqrt(two);
        assert_evalf_starts_with(&a, expr, 30, "1.4142135623730950488");
    }

    #[test]
    fn sin_pi_is_near_zero() {
        let mut a = Arena::new();
        let pi = a.pi;
        let expr = a.sin(pi);
        let result = evalf(&a, expr, 20).unwrap();
        // sin(pi) is mathematically 0, numerically very small.
        let val: f64 = result.parse().unwrap_or(999.0);
        assert!(val.abs() < 1e-15, "sin(pi) should be ~0, got: {result}");
    }

    // ── Neg and Abs ─────────────────────────────────────────────────

    #[test]
    fn neg_pi() {
        let mut a = Arena::new();
        let pi = a.pi;
        let expr = a.neg(pi);
        assert_evalf_starts_with(&a, expr, 10, "-3.14159265");
    }

    #[test]
    fn abs_neg_pi() {
        let mut a = Arena::new();
        let pi = a.pi;
        let neg_pi = a.neg(pi);
        let expr = a.abs(neg_pi);
        assert_evalf_starts_with(&a, expr, 10, "3.14159265");
    }

    #[test]
    fn abs_negative_integer() {
        let mut a = Arena::new();
        let neg_five = a.int(-5);
        let expr = a.abs(neg_five);
        assert_evalf_eq(&a, expr, 10, "5");
    }

    // ── Error cases ─────────────────────────────────────────────────

    #[test]
    fn free_symbol_errors() {
        let mut a = Arena::new();
        let x = a.symbol("x");
        let result = evalf(&a, x, 10);
        assert!(result.is_err(), "evalf of free symbol should error");
    }

    #[test]
    fn infinity_errors() {
        let a = Arena::new();
        let result = evalf(&a, a.infinity, 10);
        assert!(result.is_err(), "evalf of infinity should error");
    }

    #[test]
    fn nan_errors() {
        let a = Arena::new();
        let result = evalf(&a, a.nan, 10);
        assert!(result.is_err(), "evalf of NaN should error");
    }

    #[test]
    fn imaginary_unit_works() {
        let a = Arena::new();
        let result = evalf(&a, a.i_unit, 10);
        assert!(result.is_ok(), "evalf of i should work now");
        let s = result.unwrap();
        assert!(
            s.contains("i") || s.contains("I"),
            "should display as imaginary: {s}"
        );
    }

    // ── Complex number tests ────────────────────────────────────────

    #[test]
    fn evalf_one_plus_i() {
        let mut a = Arena::new();
        let one = a.one;
        let i = a.i_unit;
        let expr = a.add(&[one, i]);
        let result = evalf(&a, expr, 10);
        assert!(result.is_ok(), "evalf(1+i) should work: {:?}", result.err());
        let s = result.unwrap();
        assert!(
            s.contains("i") || s.contains("I"),
            "1+i should contain 'i': {s}"
        );
    }

    #[test]
    fn evalf_exp_i_pi() {
        // e^(i*pi) ≈ -1
        let mut a = Arena::new();
        let i = a.i_unit;
        let pi = a.pi;
        let i_times_pi = a.mul(&[i, pi]);
        let expr = a.exp(i_times_pi);
        let result = evalf(&a, expr, 15);
        assert!(
            result.is_ok(),
            "evalf(exp(i*pi)) should work: {:?}",
            result.err()
        );
        let s = result.unwrap();
        let val: f64 = s.parse().unwrap_or(999.0);
        assert!(
            (val - (-1.0)).abs() < 1e-10,
            "exp(i*pi) should be ~-1, got: {s}"
        );
    }

    #[test]
    fn evalf_i_squared() {
        // i^2 = -1 (may be canonicalized by the arena)
        let mut a = Arena::new();
        let i = a.i_unit;
        let two = a.int(2);
        let expr = a.pow(i, two);
        let result = evalf(&a, expr, 10);
        assert!(result.is_ok(), "evalf(i^2) should work: {:?}", result.err());
        let s = result.unwrap();
        let val: f64 = s.parse().unwrap_or(999.0);
        assert!((val - (-1.0)).abs() < 1e-10, "i^2 should be -1, got: {s}");
    }

    #[test]
    fn evalf_sqrt_neg_one() {
        // sqrt(-1) = (-1)^(1/2) should give i
        let mut a = Arena::new();
        let neg_one = a.neg_one;
        let expr = a.sqrt(neg_one);
        let result = evalf(&a, expr, 10);
        assert!(
            result.is_ok(),
            "evalf(sqrt(-1)) should work: {:?}",
            result.err()
        );
        let s = result.unwrap();
        assert!(
            s.contains("i") || s.contains("I"),
            "sqrt(-1) should be imaginary: {s}"
        );
    }

    #[test]
    fn evalf_abs_3_plus_4i() {
        // |3 + 4i| = 5
        let mut a = Arena::new();
        let three = a.int(3);
        let four = a.int(4);
        let i = a.i_unit;
        let four_i = a.mul(&[four, i]);
        let sum = a.add(&[three, four_i]);
        let expr = a.abs(sum);
        let result = evalf(&a, expr, 10);
        assert!(
            result.is_ok(),
            "evalf(|3+4i|) should work: {:?}",
            result.err()
        );
        let s = result.unwrap();
        let val: f64 = s.parse().unwrap_or(999.0);
        assert!((val - 5.0).abs() < 1e-10, "|3+4i| should be 5, got: {s}");
    }

    // ── Composite expressions ───────────────────────────────────────

    #[test]
    fn polynomial_at_value() {
        // Evaluate x^2 + 2*x + 1 at x=3 → 16.
        // We substitute first, then evalf.
        let mut a = Arena::new();
        let x = a.symbol("x");
        let two = a.int(2);
        let x_sq = a.pow(x, two);
        let two_x = a.mul(&[two, x]);
        let expr = a.add(&[x_sq, two_x, a.one]);

        // Substitute x → 3.
        let three = a.int(3);
        let evaluated = crate::transforms::subs::subs(&mut a, expr, x, three);
        assert_evalf_eq(&a, evaluated, 10, "16");
    }

    #[test]
    fn sin_cos_identity() {
        // sin^2(1) + cos^2(1) should be very close to 1.
        let mut a = Arena::new();
        let one = a.one;
        let two = a.int(2);
        let sin_1 = a.sin(one);
        let cos_1 = a.cos(one);
        let sin_sq = a.pow(sin_1, two);
        let cos_sq = a.pow(cos_1, two);
        let expr = a.add(&[sin_sq, cos_sq]);
        let result = evalf(&a, expr, 15).unwrap();
        let val: f64 = result.parse().unwrap_or(0.0);
        assert!(
            (val - 1.0).abs() < 1e-12,
            "sin^2(1) + cos^2(1) should be ~1, got: {result}"
        );
    }

    // ── Idempotence ─────────────────────────────────────────────────

    #[test]
    fn same_precision_same_result() {
        let a = Arena::new();
        let r1 = evalf(&a, a.pi, 20).unwrap();
        let r2 = evalf(&a, a.pi, 20).unwrap();
        assert_eq!(r1, r2, "same precision should give same result");
    }

    // ── Bessel function evaluation ─────────────────────────────────

    #[test]
    fn bessel_j0_at_zero() {
        // J_0(0) = 1
        let mut a = Arena::new();
        let zero = a.zero;
        let j0_0 = a.besselj(zero, zero);
        let result = evalf(&a, j0_0, 15).unwrap();
        let val: f64 = result.parse().unwrap_or(f64::NAN);
        assert!(
            (val - 1.0).abs() < 1e-10,
            "J_0(0) should be 1.0, got {result}"
        );
    }

    #[test]
    fn bessel_j0_at_one() {
        // J_0(1) ≈ 0.7651976865579666
        let mut a = Arena::new();
        let zero = a.zero;
        let one = a.one;
        let j0_1 = a.besselj(zero, one);
        let result = evalf(&a, j0_1, 15).unwrap();
        let val: f64 = result.parse().unwrap_or(f64::NAN);
        assert!(
            (val - 0.7651976865579666).abs() < 1e-8,
            "J_0(1) should be ~0.7652, got {result}"
        );
    }

    #[test]
    fn bessel_j1_at_zero() {
        // J_1(0) = 0
        let mut a = Arena::new();
        let one = a.one;
        let zero = a.zero;
        let j1_0 = a.besselj(one, zero);
        let result = evalf(&a, j1_0, 15).unwrap();
        let val: f64 = result.parse().unwrap_or(f64::NAN);
        assert!(val.abs() < 1e-10, "J_1(0) should be 0, got {result}");
    }

    #[test]
    fn bessel_j1_at_one() {
        // J_1(1) ≈ 0.44005058574493355
        let mut a = Arena::new();
        let one = a.one;
        let j1_1 = a.besselj(one, one);
        let result = evalf(&a, j1_1, 15).unwrap();
        let val: f64 = result.parse().unwrap_or(f64::NAN);
        assert!(
            (val - 0.44005058574493355).abs() < 1e-8,
            "J_1(1) should be ~0.4401, got {result}"
        );
    }

    #[test]
    fn bessel_y0_at_one() {
        // Y_0(1) ≈ 0.08825696421567696
        let mut a = Arena::new();
        let zero = a.zero;
        let one = a.one;
        let y0_1 = a.bessely(zero, one);
        let result = evalf(&a, y0_1, 15).unwrap();
        let val: f64 = result.parse().unwrap_or(f64::NAN);
        assert!(
            (val - 0.08825696421567696).abs() < 1e-6,
            "Y_0(1) should be ~0.0883, got {result}"
        );
    }

    // ── 0.2: named constants ───────────────────────────────────────────────

    #[test]
    fn euler_gamma_and_catalan_high_precision() {
        let a = Arena::new();
        assert_evalf_starts_with(
            &a,
            a.euler_gamma,
            60,
            "0.57721566490153286060651209008240243104215933593992359880576",
        );
        assert_evalf_starts_with(
            &a,
            a.catalan,
            60,
            "0.91596559417721901505460351493238411077414937428167213426649",
        );
        assert_evalf_starts_with(
            &a,
            a.golden_ratio,
            40,
            "1.61803398874989484820458683436563811772",
        );
    }

    #[test]
    fn bigint_to_bigfloat_is_exact_for_huge_integers() {
        // 2^200 + 1 is far beyond i128; its BigFloat must keep all bits.
        let n = (BigInt::from(1) << 200) + BigInt::from(1);
        let bf = bigint_to_bigfloat(&n, 256);
        let two200 = BigFloat::from_i32(2, 256).powi(200, 256, RoundingMode::ToEven);
        let diff = bf.sub(&two200, 256, RoundingMode::ToEven);
        assert_eq!(diff, BigFloat::from_i32(1, 256));
    }

    #[test]
    fn bf_pow_handles_exactly_representable_results() {
        // BigFloat::pow would hang on 4^(1/2) = 2; bf_pow must not.
        let mut cc = Consts::new().unwrap();
        let rm = RoundingMode::ToEven;
        let four = BigFloat::from_i32(4, 160);
        let half = BigFloat::from_f64(0.5, 160);
        let r = bf_pow(&four, &half, 160, rm, &mut cc);
        let diff = r.sub(&BigFloat::from_i32(2, 160), 160, rm);
        assert!(
            diff.is_zero() || diff.exponent().is_none_or(|e| e < -140),
            "4^0.5 = {r}"
        );
        // integer exponent path
        let three = BigFloat::from_i32(3, 160);
        let r = bf_pow(&three, &BigFloat::from_i32(5, 160), 160, rm, &mut cc);
        assert_eq!(r, BigFloat::from_i32(243, 160));
        let r = bf_pow(&three, &BigFloat::from_i32(-1, 160), 160, rm, &mut cc);
        let back = r.mul(&three, 160, rm);
        let d = back.sub(&BigFloat::from_i32(1, 160), 160, rm);
        assert!(d.is_zero() || d.exponent().is_none_or(|e| e < -140));
    }

    // ── 0.2: special functions ──────────────────────────────────────────────

    #[test]
    fn si_ci_ei_li_at_one_and_two() {
        let mut a = Arena::new();
        let one = a.one;
        let two = a.int(2);
        let si = a.si(one);
        assert_evalf_starts_with(&a, si, 20, "0.94608307036718301494");
        let ci = a.ci(one);
        assert_evalf_starts_with(&a, ci, 20, "0.33740392290096813466");
        let ei = a.ei(one);
        assert_evalf_starts_with(&a, ei, 20, "1.8951178163559367555");
        let li = a.li(two);
        assert_evalf_starts_with(&a, li, 20, "1.0451637801174927848");
    }

    #[test]
    fn si_asymptotic_and_series_regimes_agree() {
        // Evaluate Si(40) at two precisions: at low precision the asymptotic
        // expansion is used (40 > 160·ln2 ≈ 111 is false …), so force the
        // regimes by comparing against the same value at higher precision
        // where the series is used.
        let mut a = Arena::new();
        let x = a.int(200);
        let si = a.si(x);
        let lo = evalf(&a, si, 12).unwrap(); // asymptotic (200 > 128·ln2)
        let hi = evalf(&a, si, 120).unwrap(); // series (200 < 472·ln2)
        assert!(hi.starts_with(&lo[..12]), "{lo} vs {hi}");
    }

    #[test]
    fn zeta_borwein_and_functional_equation() {
        let mut a = Arena::new();
        let three = a.int(3);
        let z3 = a.zeta(three);
        assert_evalf_starts_with(&a, z3, 30, "1.20205690315959428539973816151");
        let half = a.rational(1, 2);
        let zh = a.zeta(half);
        assert_evalf_starts_with(&a, zh, 20, "-1.4603545088095868129");
        // ζ(−5/2) > 0 (ζ is positive on (−4, −2); ζ(−3) = 1/120).
        let neg = a.rational(-5, 2);
        let zn = a.zeta(neg);
        assert_evalf_starts_with(&a, zn, 15, "0.00851692877785");
        // ζ(−1/2) < 0
        let neg_half = a.rational(-1, 2);
        let znh = a.zeta(neg_half);
        assert_evalf_starts_with(&a, znh, 15, "-0.20788622497735");
    }

    #[test]
    fn polygamma_matches_trigamma_closed_forms() {
        let mut a = Arena::new();
        // Build the node structurally so the exact folding does not kick in.
        let one = a.one;
        let x = a.rational(7, 3);
        let pg = a.intern(ExprNode::Polygamma(one, x));
        assert_evalf_starts_with(&a, pg, 20, "0.53309712542709408179");
        // ψ'(1/2) = π²/2 through the numerical path
        let half = a.rational(1, 2);
        let pg_half = a.intern(ExprNode::Polygamma(one, half));
        assert_evalf_starts_with(&a, pg_half, 20, "4.9348022005446793094");
    }

    #[test]
    fn bessel_i_k_reference() {
        let mut a = Arena::new();
        let zero = a.zero;
        let one = a.one;
        let i0 = a.besseli(zero, one);
        assert_evalf_starts_with(&a, i0, 20, "1.2660658777520083356");
        let k0 = a.besselk(zero, one);
        assert_evalf_starts_with(&a, k0, 20, "0.42102443824070833334");
        let k1 = a.besselk(one, one);
        assert_evalf_starts_with(&a, k1, 20, "0.60190723019723457474");
    }

    #[test]
    fn orthogonal_polynomials_large_degree() {
        let mut a = Arena::new();
        let n = a.int(100);
        let third = a.rational(1, 3);
        let t = a.chebyshev_t(n, third);
        let s = evalf(&a, t, 15).unwrap();
        let v: f64 = s.parse().unwrap();
        let expected = (100.0 * (1.0f64 / 3.0).acos()).cos();
        assert!((v - expected).abs() < 1e-12, "{v} vs {expected}");
        let p = a.legendre(n, a.one);
        assert_evalf_starts_with(&a, p, 10, "1");
    }

    #[test]
    fn complex_nodes_evaluate_numerically() {
        let mut a = Arena::new();
        let z = a.symbol("z");
        let re_z = a.intern(ExprNode::Re(z));
        let arg_z = a.intern(ExprNode::Arg(z));
        let conj_z = a.intern(ExprNode::Conjugate(z));
        let three = a.int(3);
        let four = a.int(4);
        let four_i = a.mul(&[four, a.i_unit]);
        let w = a.add(&[three, four_i]);
        // Substitute structurally *without* refolding: intern the nodes
        // directly on w.
        let _ = re_z;
        let re_w = a.intern(ExprNode::Re(w));
        let arg_w = a.intern(ExprNode::Arg(w));
        let conj_w = a.intern(ExprNode::Conjugate(w));
        let _ = (arg_z, conj_z);
        assert_eq!(evalf(&a, re_w, 10).unwrap(), "3");
        assert_evalf_starts_with(&a, arg_w, 15, "0.92729521800161");
        let s = evalf(&a, conj_w, 10).unwrap();
        assert!(s.contains('3') && s.contains('4') && s.contains('-'), "{s}");
    }
}