symplex 0.43.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::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]` with adaptive
    /// Gauss–Kronrod (G7/K15) quadrature using the default [`QuadOpts`].
    ///
    /// 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.
    ///
    /// # 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`] — the quadrature did not reach
    ///   the requested accuracy (e.g. a divergent integral).
    ///
    /// # 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 > 1e3 * 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 error estimate is returned even if the tolerance was not met
    /// within `max_subdivisions`; check it before trusting the value.
    ///
    /// # 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.  Before 0.43 `∫₋₁¹ δ(x) dx` came out `0`.
        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))
            })
        {
            return Err(SymplexError::Unevaluable {
                reason: "a definite integral of DiracDelta of the integration variable has no \
                         numerical value (integrate it symbolically)"
                    .into(),
            });
        }

        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)
            ))
        })?;
        drop(inner);

        let f = |t: f64| func(&[t]);
        definite::quadrature(&f, a, b, opts)
    }

    /// 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)
            }
        }
    }
}