symplex 0.53.0

//! Definite / improper / numeric integration methods on [`Ex`].

use tracing::debug_span;

use crate::api::expr::{Ex, Expr, Numeric};
use crate::base::errors::SymplexError;
use crate::base::interval::Interval;
use crate::base::node::ExprNode;
use crate::calculus::definite::{self, QuadOpts, QuadResult, QuadRule};

impl Expr<Numeric> {
    /// Compute the definite integral `∫_lo^hi self dvar`.
    ///
    /// Unlike a naive `F(hi) − F(lo)`, this locates singularities of the
    /// integrand inside the interval, treats infinite bounds and endpoint
    /// singularities as improper integrals via one-sided limits, applies
    /// symmetry shortcuts, resolves `Abs`/`Sign`/`Heaviside`/`DiracDelta`/
    /// `Piecewise` integrands, and consults a table of classical improper
    /// integrals when no elementary antiderivative exists.
    ///
    /// If the integral is divergent or cannot be evaluated the result is an
    /// **unevaluated** definite-integral node `Integral(f, x, lo, hi)`
    /// (see [`definite_integral_node`](Self::definite_integral_node)) —
    /// never a wrong finite number.  Use
    /// [`try_integrate_definite`](Self::try_integrate_definite) to
    /// distinguish "diverges" from "could not compute".
    ///
    /// With parameters, a closed form is returned only where it holds for
    /// every allowed value of them: no principal branch in the
    /// antiderivative (`ln`, `atan`, powers) can jump along `[lo, hi]`, and
    /// no pole of the integrand can lie on it.  A parameter is what its
    /// assumptions say, so an unassumed one may be complex (decision D4):
    /// `∫₀¹ dx/(x + a − i)` stays unevaluated, since the pole `x = −a + i`
    /// is on `[0, 1]` for `a = −1/2 + i`, while for a declared-real `a` it
    /// is `ln(a − i + 1) − ln(a − i)`.  Where that cannot be decided the
    /// integral is left unevaluated rather than split into cases.
    ///
    /// If `self` is, or contains, such a node, the inner integrals are
    /// evaluated first (innermost out), so nested integrals can be built up
    /// with repeated calls.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    ///
    /// // ∫₀¹ x² dx = 1/3
    /// let v = x.powi(2).integrate_definite(&x, &ctx.int(0), &ctx.int(1));
    /// assert_eq!(format!("{v}"), "1/3");
    ///
    /// // ∫₀^∞ e^{−x} dx = 1
    /// let v = (-&x).exp().integrate_definite(&x, &ctx.int(0), &ctx.infinity());
    /// assert_eq!(format!("{v}"), "1");
    ///
    /// // ∫₋∞^∞ e^{−x²} dx = √π
    /// let v = (-x.powi(2)).exp().integrate_definite(&x, &ctx.neg_infinity(), &ctx.infinity());
    /// assert_eq!(format!("{v}"), "sqrt(pi)");
    ///
    /// // ∫₋₁¹ x⁻² dx diverges: the result stays unevaluated, bounds intact.
    /// let v = x.powi(-2).integrate_definite(&x, &ctx.int(-1), &ctx.int(1));
    /// assert!(v.has_unevaluated());
    /// assert!(v.is_definite_integral());
    /// assert_eq!(format!("{v}"), "Integral(x^(-2), x, -1, 1)");
    ///
    /// // ∫₀¹ xˣ dx has no closed form; the node still has a numeric value.
    /// let v = x.pow(&x).integrate_definite(&x, &ctx.int(0), &ctx.int(1));
    /// assert!(v.is_definite_integral());
    /// assert!((v.eval_f64().unwrap() - 0.7834305107).abs() < 1e-8);
    /// ```
    #[must_use = "returns the definite integral value"]
    pub fn integrate_definite(&self, var: &Ex, lo: &Ex, hi: &Ex) -> Ex {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lo);
        let hi_id = self.checked_id(hi);
        let _span =
            debug_span!("integrate_definite", expr = ?self.raw_id(), var = ?var_id).entered();
        let mut inner = self.inner.write();
        let id = match definite::integrate_definite(
            &mut inner.arena,
            self.raw_id(),
            var_id,
            lo_id,
            hi_id,
        ) {
            Ok(id) => id,
            Err(_) => definite::unevaluated_definite(
                &mut inner.arena,
                self.raw_id(),
                var_id,
                lo_id,
                hi_id,
            ),
        };
        drop(inner);
        self.wrap(id)
    }

    /// Like [`integrate_definite`](Self::integrate_definite), but returns
    /// `Err` instead of an unevaluated form.
    ///
    /// * [`SymplexError::Divergent`] — the integral was proven to diverge,
    /// * [`SymplexError::ComputationFailed`] — no closed form could be
    ///   established, or the result still contains an unevaluated form
    ///   (a `DefiniteIntegral`, `Integral`, `Limit`, … node),
    /// * [`SymplexError::InvalidArgument`] — `var` is not a symbol.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    ///
    /// let r = x.powi(-2).try_integrate_definite(&x, &ctx.int(-1), &ctx.int(1));
    /// assert!(matches!(r, Err(SymplexError::Divergent { .. })));
    ///
    /// let v = x.ln().try_integrate_definite(&x, &ctx.int(0), &ctx.int(1)).unwrap();
    /// assert_eq!(format!("{v}"), "-1");
    /// ```
    pub fn try_integrate_definite(&self, var: &Ex, lo: &Ex, hi: &Ex) -> Result<Ex, SymplexError> {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lo);
        let hi_id = self.checked_id(hi);
        let _span =
            debug_span!("try_integrate_definite", expr = ?self.raw_id(), var = ?var_id).entered();
        let mut inner = self.inner.write();
        let id =
            definite::integrate_definite(&mut inner.arena, self.raw_id(), var_id, lo_id, hi_id)?;
        let unevaluated = crate::base::walk::has_unevaluated(&inner.arena, id);
        drop(inner);
        if unevaluated {
            return Err(SymplexError::ComputationFailed {
                operation: "integrate_definite",
                reason: "result contains unevaluated forms".into(),
            });
        }
        Ok(self.wrap(id))
    }

    /// Build the formal, unevaluated definite integral `∫_lo^hi self dvar`
    /// **without** attempting to evaluate it.
    ///
    /// Useful for display, LaTeX, and formal manipulation (differentiation
    /// by the Leibniz rule, substitution into the bounds, numeric
    /// evaluation by quadrature).  Only the cheap structural folds are
    /// applied: `lo == hi` gives `0`, an integrand free of `var` over a
    /// finite interval gives `self · (hi − lo)`, and numeric bounds with
    /// `lo > hi` are reordered with a sign flip.  To evaluate the node
    /// later, call [`eval_integrals`](Self::eval_integrals) on it (or on any
    /// expression containing it).
    ///
    /// The integration variable is bound inside the integrand: it is not
    /// reported by [`free_symbols`](Self::free_symbols) and is not touched
    /// by [`subs`](Self::subs); the bounds are in the enclosing scope.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    /// let t = ctx.symbol("t");
    ///
    /// let node = x.sin().definite_integral_node(&x, &ctx.int(0), &t);
    /// assert!(node.is_definite_integral());
    /// assert!(node.has_unevaluated());
    /// assert_eq!(format!("{node}"), "Integral(sin(x), x, 0, t)");
    /// assert_eq!(node.to_latex(), r"\int_{0}^{t} \sin\left(x\right)\, dx");
    ///
    /// // Only `t` is free; `x` is bound.
    /// assert_eq!(node.free_symbols().len(), 1);
    ///
    /// // Leibniz rule: d/dt ∫₀ᵗ sin(x) dx = sin(t).
    /// assert_eq!(format!("{}", node.diff(&t)), "sin(t)");
    ///
    /// // Evaluating recovers the closed form 1 − cos(t).
    /// let v = node.eval_integrals();
    /// assert_eq!(format!("{v}"), "-cos(t) + 1");
    /// ```
    #[must_use = "returns the formal definite integral node"]
    pub fn definite_integral_node(&self, var: &Ex, lo: &Ex, hi: &Ex) -> Ex {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lo);
        let hi_id = self.checked_id(hi);
        let mut inner = self.inner.write();
        let id = inner
            .arena
            .definite_integral(self.raw_id(), var_id, lo_id, hi_id);
        drop(inner);
        self.wrap(id)
    }

    /// Is this expression an unevaluated definite integral node
    /// (`Integral(f, x, lo, hi)`)?
    ///
    /// Only the root node is inspected; use
    /// [`has_unevaluated`](Self::has_unevaluated) to search the whole tree.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    /// // No closed form → the node comes back.
    /// let v = x.pow(&x).integrate_definite(&x, &ctx.int(0), &ctx.int(1));
    /// assert!(v.is_definite_integral());
    /// // Closed form → a number.
    /// let v = x.integrate_definite(&x, &ctx.int(0), &ctx.int(1));
    /// assert!(!v.is_definite_integral());
    /// ```
    #[must_use]
    pub fn is_definite_integral(&self) -> bool {
        let inner = self.inner.read();
        matches!(
            inner.arena.node(self.raw_id()),
            ExprNode::DefiniteIntegral(..)
        )
    }

    /// Evaluate every formal definite-integral node in this expression,
    /// innermost first.
    ///
    /// This is the "doit" operation for `Integral(f, x, lo, hi)` nodes (the
    /// analogue of [`eval_derivatives`](Self::eval_derivatives) for
    /// `Derivative`): each node is run through the definite integrator and
    /// replaced by its closed form.  Nodes that still cannot be evaluated —
    /// no closed form, or a proven divergence — are left in place, so the
    /// result is never a wrong finite number; check with
    /// [`has_unevaluated`](Self::has_unevaluated).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    /// let t = ctx.symbol("t");
    ///
    /// // d/dt ∫₀¹ sin(t·x) dx = ∫₀¹ x·cos(t·x) dx, then evaluate it.
    /// let node = (&t * &x).sin().definite_integral_node(&x, &ctx.int(0), &ctx.int(1));
    /// let d = node.diff(&t);
    /// assert!(d.is_definite_integral());
    /// let v = d.eval_integrals();
    /// assert!(!v.has_unevaluated(), "{v}");
    ///
    /// // No closed form: the node is returned unchanged.
    /// let n = x.pow(&x).definite_integral_node(&x, &ctx.int(0), &ctx.int(1));
    /// assert!(n.eval_integrals().is_definite_integral());
    /// ```
    #[must_use = "returns the expression with definite integrals evaluated"]
    pub fn eval_integrals(&self) -> Ex {
        let _span = debug_span!("eval_integrals", expr = ?self.raw_id()).entered();
        let id = {
            let mut inner = self.inner.write();
            definite::evaluate_inner_definite(&mut inner.arena, self.raw_id())
        };
        self.wrap(id)
    }

    /// Numerically integrate `self` over `[lo, hi]` to the default
    /// [`QuadOpts`] tolerance, `max(10⁻¹², 10⁻¹⁰·|value|)`: adaptive
    /// Gauss–Kronrod (G7/K15) quadrature in `f64`, and where its error
    /// estimate exceeds the tolerance the arbitrary-precision quadrature of
    /// [`integrate_numeric_with`](Self::integrate_numeric_with).
    ///
    /// Infinite bounds are supported.  The integrand is compiled with
    /// [`compile`](Self::compile), so it must contain no free symbols other
    /// than `var` and only nodes the compiler supports.
    ///
    /// A value is returned only when its error estimate is within the
    /// tolerance.  Before 0.44 an estimate up to 1000 times the tolerance
    /// was accepted, and next to a singularity at a nonzero end the estimate
    /// itself was too small: `∫₀¹ dx/√(1 − x)` was `1.9999999898808887`,
    /// `∫₀¹ log²(1 − x)/√(1 − x) dx` `15.99998` (both now exact to the last
    /// digit or two).
    ///
    /// # Errors
    ///
    /// * [`SymplexError::FreeSymbol`] — another free symbol is present,
    /// * [`SymplexError::NotImplemented`] — the integrand contains a node
    ///   that cannot be compiled to `f64` code,
    /// * [`SymplexError::Unevaluable`] — a bound is not a real number, or the
    ///   integrand contains `DiracDelta` of `var` (a distribution, 0 at every
    ///   point a quadrature samples; before 0.43 `∫₋₁¹ δ(x) dx` was `0`),
    /// * [`SymplexError::ComputationFailed`] — neither quadrature reached
    ///   the tolerance (e.g. a divergent integral), the integrand has more
    ///   than 10,000 poles, kinks or steps in the range (`Γ(9·10¹⁸x)` on
    ///   `[−1, 1]`: up to 0.49 a search of minutes), or a kink or step the
    ///   breakpoint scan cannot place could not be located numerically,
    /// * [`SymplexError::Divergent`] — the arbitrary-precision quadrature
    ///   met terms that do not decay and a comparison test proves the
    ///   divergence (`(x − c)·f(x)`, `x·f(x)` at `±∞`, or with a factor
    ///   `ln`, tends to a nonzero limit).
    ///
    /// Bounds that round to the same `f64`, or lie within `2⁻²⁰` of each
    /// other relatively, are integrated by the arbitrary-precision quadrature
    /// with the bounds as they are, and the rounding of the bounds to `f64`
    /// counts in the error estimate otherwise: up to 0.49 `∫|sin x| dx` over
    /// `[10²⁰, 10²⁰ + 4]` was `Ok(0.0)` (it is `2.2483…`).  The integrand is
    /// split at its breakpoints (poles, kinks, steps), the periodic ones
    /// enumerated between the exact bounds; before 0.50 only an integrand
    /// with a kink was split, and `∫₋₂² tan x dx`, which diverges, was
    /// `Ok(0.0)`.
    ///
    /// On an infinite range the `f64` value is returned only when the
    /// integrand's tail beyond `|x| = 2¹⁶`, evaluated with certified bounds
    /// out to `2²⁰⁴⁸`, carries less than the tolerance or was seen by the
    /// compiled integrand.  Before 0.47 an overflow in the compiled integrand
    /// cut a divergent integral off: `∫_{−∞}^{∞} |x|/((x² + 1)π) dx` was
    /// `225.93` (`|x|/(x² + 1)` is 0 in `f64` once `x²` overflows).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    ///
    /// let v = x.sin().integrate_numeric(&x, &ctx.int(0), &ctx.pi()).unwrap();
    /// assert!((v - 2.0).abs() < 1e-10);
    ///
    /// // ∫₀^∞ e^{−x²} dx = √π/2
    /// let v = (-x.powi(2)).exp().integrate_numeric(&x, &ctx.int(0), &ctx.infinity()).unwrap();
    /// assert!((v - std::f64::consts::PI.sqrt() / 2.0).abs() < 1e-10);
    ///
    /// // A divergent integral is reported as an error, not a number.
    /// assert!(x.powi(-2).integrate_numeric(&x, &ctx.int(-1), &ctx.int(1)).is_err());
    /// ```
    ///
    /// Conditionally convergent or slowly decaying oscillatory tails (e.g.
    /// `sin(x)/x` on `[0, ∞)`) are beyond plain adaptive quadrature and
    /// also return an error; use [`integrate_definite`](Self::integrate_definite)
    /// for those.
    pub fn integrate_numeric(&self, var: &Ex, lo: &Ex, hi: &Ex) -> Result<f64, SymplexError> {
        let opts = QuadOpts::default();
        let QuadResult { value, error: err } = self.integrate_numeric_with(var, lo, hi, &opts)?;
        let tol = opts.abs_tol.max(opts.rel_tol * value.abs());
        if err > tol {
            return Err(SymplexError::ComputationFailed {
                operation: "integrate_numeric",
                reason: format!(
                    "quadrature did not converge: estimate {value} with error {err:e} (integral may diverge)"
                ),
            });
        }
        Ok(value)
    }

    /// Numerically integrate with explicit [`QuadOpts`], returning the
    /// estimate and its error estimate as a [`QuadResult`].
    ///
    /// The `f64` adaptive Gauss–Kronrod quadrature runs first (see
    /// [`QuadResult::error`] for what its estimate covers), split at the
    /// centres of narrow peaks — complex poles next to the range, narrow
    /// Gaussians — and, for an integrand with a kink or a step (`|·|`,
    /// `sign`, `Heaviside`, `floor`, `Piecewise`, …), at its breakpoints.
    /// When its estimate exceeds the tolerance — an endpoint or interior singularity
    /// the `f64` bisection cannot resolve, a tolerance beyond `f64` — the
    /// integral is evaluated as [`integrate_numeric_decimal`](Self::integrate_numeric_decimal)
    /// evaluates it (tanh-sinh quadrature over arbitrary-precision floats,
    /// breakpoints split off, the bounds exact) to 17 certified digits, and
    /// that value is returned with an error of two units in its last place
    /// (`2⁻⁶⁰` for 0).  Otherwise the `f64` estimate is returned with its
    /// error estimate even if the tolerance was not met; check it before
    /// trusting the value.
    ///
    /// Before 0.44 the error estimate next to a singularity at a nonzero
    /// end was 30 to 60 times smaller than the error: `∫₀¹ dx/√(1 − x)` was
    /// `1.9999999898808887` with an estimate of `3.2·10⁻¹⁰`, the Beta moment
    /// `∫₀¹ x⁴e^{x/10}/(√(1 − x)·B(5, ½)) dx` `1.0952385480452325` with
    /// `4.4·10⁻¹⁰` (truly `1.0952385618059876`, mpmath); and a peak no node
    /// saw was missed: `exp(−10⁶(x − 2/3)²)·(x + 1)²cos(4x/3)` on `[0, 2]`
    /// integrated to `0` with an estimate of 0 (truly `0.0031`).
    ///
    /// # Errors
    ///
    /// As [`integrate_numeric`](Self::integrate_numeric), except that a
    /// result short of the tolerance is returned with its estimate.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    /// use symplex::definite::QuadOpts;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    /// let opts = QuadOpts { rel_tol: 1e-6, ..QuadOpts::default() };
    /// let r = (-x.powi(2)).exp()
    ///     .integrate_numeric_with(&x, &ctx.neg_infinity(), &ctx.infinity(), &opts)
    ///     .unwrap();
    /// assert!((r.value - std::f64::consts::PI.sqrt()).abs() < 1e-6);
    /// assert!(r.error < 1e-4);
    /// ```
    pub fn integrate_numeric_with(
        &self,
        var: &Ex,
        lo: &Ex,
        hi: &Ex,
        opts: &QuadOpts,
    ) -> Result<QuadResult, SymplexError> {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lo);
        let hi_id = self.checked_id(hi);

        let mut inner = self.inner.write();
        let arena = &mut inner.arena;

        let var_name = match arena.node(var_id) {
            ExprNode::Symbol(sid) => arena.symbol_name(*sid).to_string(),
            _ => {
                return Err(SymplexError::InvalidArgument {
                    operation: "integrate_numeric",
                    reason: "integration variable must be a symbol".into(),
                });
            }
        };

        // Bounds → f64 (infinities allowed).
        let bound_f64 = |arena: &mut crate::base::arena::Arena, id| -> Result<f64, SymplexError> {
            if id == arena.infinity() {
                return Ok(f64::INFINITY);
            }
            if id == arena.neg_infinity() {
                return Ok(f64::NEG_INFINITY);
            }
            crate::transforms::evalf::eval_const_f64(arena, id).ok_or_else(|| {
                SymplexError::Unevaluable {
                    reason: format!(
                        "integration bound {} is not a real number",
                        arena.display(id)
                    ),
                }
            })
        };
        let a = bound_f64(arena, lo_id)?;
        let b = bound_f64(arena, hi_id)?;

        // Free symbols other than the variable.
        let evaled = crate::transforms::eval::eval(arena, self.raw_id());
        for s in crate::base::walk::free_symbols(arena, evaled) {
            if s != var_id
                && let ExprNode::Symbol(sid) = arena.node(s)
            {
                return Err(SymplexError::FreeSymbol {
                    name: arena.symbol_name(*sid).to_string(),
                });
            }
        }

        // `DiracDelta` compiles to 0 at every point: its integral is not a
        // sum of values, and only its closed form has one (`evalf`'s exact
        // route; refused without).  Before 0.43 `∫₋₁¹ δ(x) dx` came out `0`,
        // and up to 0.51 it was refused.
        if crate::base::walk::post_order_ids(arena, evaled)
            .iter()
            .any(|&id| {
                matches!(arena.node(id), ExprNode::DiracDelta(inner)
                if crate::base::walk::contains(arena, *inner, var_id))
            })
        {
            drop(inner);
            return self.certified_quadrature(var, lo, hi);
        }

        let func = crate::output::lambdify::compile(arena, evaled, &[&var_name]).map_err(|e| {
            SymplexError::NotImplemented(format!(
                "integrand {} cannot be compiled for numeric quadrature: {e}",
                arena.display(evaled)
            ))
        })?;
        let (lo_f, hi_f) = if a <= b { (a, b) } else { (b, a) };
        let lo_x = crate::transforms::eval::eval(arena, lo_id);
        let hi_x = crate::transforms::eval::eval(arena, hi_id);
        // How far the `f64` bounds are from the bounds themselves.
        let moved =
            definite::bound_rounding(arena, lo_x, a) + definite::bound_rounding(arena, hi_x, b);
        let exact = (a.is_finite() && b.is_finite()).then(|| {
            if a <= b {
                Interval::closed(lo_x, hi_x)
            } else {
                Interval::closed(hi_x, lo_x)
            }
        });
        let points = definite::quadrature_split_points(arena, evaled, var_id, lo_f, hi_f, exact)?;
        drop(inner);

        // Distinct bounds with one `f64` value (`[10²⁰, 10²⁰ + 4]`), or a
        // few thousand ulps apart: the `f64` quadrature has no interval, or
        // too few `f64` nodes, to integrate over.  Up to 0.49 it returned
        // `Ok(0.0)` for the first; the arbitrary-precision quadrature takes
        // the bounds as they are.
        if definite::f64_narrow(lo_f, hi_f) && lo_x != hi_x {
            return self.certified_quadrature(var, lo, hi);
        }

        let peak = std::cell::Cell::new(0.0f64);
        let f = |t: f64| {
            let v = func(&[t]);
            if v.is_finite() {
                peak.set(peak.get().max(v.abs()));
            }
            v
        };
        // The bounds' rounding to `f64` moves the range: its effect is in the
        // error estimate (at most the largest value seen times the shift).
        let with_bounds = |r: QuadResult| QuadResult {
            error: r.error + 2.0 * moved * peak.get(),
            ..r
        };
        let tolerance = |r: &QuadResult| opts.abs_tol.max(opts.rel_tol * r.value.abs());
        // On an infinite range the compiled body is cut off where an
        // intermediate value overflows (`|x|/(x² + 1)` is 0 beyond `x ≈
        // 1.3·10¹⁵⁴`) and no node lies beyond `f64::MAX`: the integrand's
        // far tail is checked by certified evaluations.
        let tail_seen = |r: &QuadResult| {
            let guard = self.inner.read();
            crate::transforms::evalf::f64_quadrature_tail_trusted(
                &guard.arena,
                evaled,
                var_id,
                a,
                b,
                &f,
                tolerance(r),
            )
        };
        let report = if a <= b {
            definite::quadrature_split(&f, a, b, &points, opts)
        } else {
            definite::quadrature_split(&f, b, a, &points, opts).map(|(r, ok)| {
                (
                    QuadResult {
                        value: -r.value,
                        ..r
                    },
                    ok,
                )
            })
        }
        .map(|(r, ok)| (with_bounds(r), ok));
        if let Ok((r, _)) = &report
            && r.error <= tolerance(r)
            && tail_seen(r)
        {
            return Ok(*r);
        }
        match self.certified_quadrature(var, lo, hi) {
            Ok(r) => return Ok(r),
            Err(e @ SymplexError::Divergent { .. }) => return Err(e),
            Err(_) => {}
        }
        match report {
            Ok((r, _)) if !tail_seen(&r) => Err(SymplexError::ComputationFailed {
                operation: "integrate_numeric",
                reason: format!(
                    "quadrature did not converge: the estimate {} misses the integrand's tail \
                     beyond |x| = 2^16, which does not decay fast enough (the integral may diverge)",
                    r.value
                ),
            }),
            report => report.map(|(r, _)| r),
        }
    }

    /// `∫_lo^hi self dvar` by the arbitrary-precision quadrature of
    /// [`integrate_numeric_decimal`](Self::integrate_numeric_decimal) to 17
    /// digits, read as an `f64` with an error of two units in its last place
    /// (`2⁻⁶⁰`, the radius `evalf` returns 0 within, for 0); its error when
    /// it fails, [`SymplexError::Unevaluable`] when the value is not real.
    fn certified_quadrature(&self, var: &Ex, lo: &Ex, hi: &Ex) -> Result<QuadResult, SymplexError> {
        let node = self.definite_integral_node(var, lo, hi).eval();
        let printed = {
            let guard = node.inner.read();
            crate::transforms::evalf::evalf_quadrature(&guard.arena, node.raw_id(), 17, false)
        }?;
        let not_real = || SymplexError::Unevaluable {
            reason: format!("the integral is not a real number: {printed}"),
        };
        let value: f64 = printed.parse().map_err(|_| not_real())?;
        if !value.is_finite() {
            return Err(not_real());
        }
        Ok(QuadResult {
            value,
            error: (2.0 * f64::EPSILON * value.abs()).max((-60f64).exp2()),
        })
    }

    /// Numerically integrate `self` in `var` from the first of `points` to
    /// the last, split at the points in between, to `digits` correct
    /// significant digits — mpmath's `quad(f, [a, p₁, …, b])` with
    /// symplex's certification.  The bounds may be expressions or `±∞`.
    ///
    /// Each piece is an unevaluated definite integral evaluated as
    /// [`eval_decimal`](Self::eval_decimal) evaluates one: by
    /// arbitrary-precision quadrature with `rule` (tanh-sinh by default; see
    /// [`QuadRule`]), its own breakpoints split off, two successive levels
    /// required to agree to the working precision, and only the digits the
    /// resulting bound covers printed.  Supply the points where the
    /// integrand is not smooth when the breakpoint scanner cannot find them.
    ///
    /// # Errors
    ///
    /// * [`SymplexError::InvalidArgument`] — fewer than two points, or `var`
    ///   is not a symbol,
    /// * as [`eval_decimal`](Self::eval_decimal): a free symbol other than
    ///   `var`, [`SymplexError::PrecisionExhausted`] when the digits cannot
    ///   be certified, [`SymplexError::ComputationFailed`] when the
    ///   quadrature does not converge (a divergent integral, an oscillating
    ///   integrand on an infinite range).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    /// use symplex::definite::QuadRule;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    /// // ∫₀¹ x^(−9/10) dx = 10: an endpoint singularity at full precision.
    /// let f = x.pow(&ctx.rational(-9, 10));
    /// let v = f.integrate_numeric_decimal(&x, &[&ctx.int(0), &ctx.int(1)], 30, QuadRule::TanhSinh);
    /// assert_eq!(v.unwrap(), "10");
    /// // ∫₀² |x − 1| dx = 1, the kink supplied.
    /// let g = (&x - 1).abs();
    /// let pts = [&ctx.int(0), &ctx.int(1), &ctx.int(2)];
    /// assert_eq!(g.integrate_numeric_decimal(&x, &pts, 20, QuadRule::default()).unwrap(), "1");
    /// ```
    pub fn integrate_numeric_decimal(
        &self,
        var: &Ex,
        points: &[&Ex],
        digits: u32,
        rule: QuadRule,
    ) -> Result<String, SymplexError> {
        let _span =
            debug_span!("integrate_numeric_decimal", expr = ?self.raw_id(), digits).entered();
        if points.len() < 2 {
            return Err(SymplexError::InvalidArgument {
                operation: "integrate_numeric_decimal",
                reason: "at least two points (the bounds) are required".into(),
            });
        }
        {
            let var_id = self.checked_id(var);
            let inner = self.inner.read();
            if !matches!(inner.arena.node(var_id), ExprNode::Symbol(_)) {
                return Err(SymplexError::InvalidArgument {
                    operation: "integrate_numeric_decimal",
                    reason: "integration variable must be a symbol".into(),
                });
            }
        }
        let mut total: Option<Ex> = None;
        for w in points.windows(2) {
            let piece = self.definite_integral_node(var, w[0], w[1]);
            total = Some(match total {
                None => piece,
                Some(t) => &t + &piece,
            });
        }
        let Some(total) = total else {
            return Err(SymplexError::InvalidArgument {
                operation: "integrate_numeric_decimal",
                reason: "at least two points (the bounds) are required".into(),
            });
        };
        let evaled = total.eval();
        let guard = evaled.inner.read();
        crate::transforms::evalf::evalf_quadrature(
            &guard.arena,
            evaled.raw_id(),
            digits,
            rule == QuadRule::GaussLegendre,
        )
    }

    /// Residue of `self` at `var = ∞`, defined as
    /// `Res_{z=∞} f(z) = −Res_{t=0} f(1/t)/t²`.
    ///
    /// The sum of all finite residues plus the residue at infinity is
    /// zero for a function meromorphic on the extended plane.  If the
    /// residue cannot be computed, a formal `Residue` node is returned.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let z = ctx.symbol("z");
    /// // f = 1/z has residue 1 at 0, hence −1 at ∞.
    /// let f = &ctx.int(1) / &z;
    /// assert_eq!(format!("{}", f.residue_at_infinity(&z)), "-1");
    /// // Polynomials: Res_{∞} z = −Res_{t=0} 1/t³ = 0.
    /// assert_eq!(format!("{}", z.residue_at_infinity(&z)), "0");
    /// ```
    #[must_use = "returns the residue value; does not modify in place"]
    pub fn residue_at_infinity(&self, var: &Ex) -> Ex {
        let var_id = self.checked_id(var);
        let _span =
            debug_span!("residue_at_infinity", expr = ?self.raw_id(), var = ?var_id).entered();
        let mut inner = self.inner.write();
        match crate::calculus::residue::residue_at_infinity(&mut inner.arena, self.raw_id(), var_id)
        {
            Ok(id) => {
                drop(inner);
                self.wrap(id)
            }
            Err(_) => {
                let inf = inner.arena.infinity();
                let id = inner
                    .arena
                    .intern(ExprNode::Residue(self.raw_id(), var_id, inf));
                drop(inner);
                self.wrap(id)
            }
        }
    }
}