symplex 0.34.0

//! Summation, products, and series extension methods on [`Ex`].
//!
//! * [`Ex::summation`] / [`Ex::try_summation`] — symbolic `Σ`
//! * [`Ex::product_over`] / [`Ex::try_product_over`] — symbolic `Π`
//! * [`Ex::hypergeometric_ratio`] — the term ratio `t(k+1)/t(k)`
//! * [`Ex::is_absolutely_convergent`] — absolute convergence of `Σ t(k)`
//! * [`Ex::series_at_infinity`] / [`Ex::series_at_neg_infinity`] — asymptotic expansions
//! * [`Ex::series_dir`] / [`Ex::try_series_dir`] — one-sided expansions
//!
//! # What a series is
//!
//! [`Ex::series`] and the methods here compute the **asymptotic expansion
//! as `var → point` through real values** (two-sided for `series`, one-sided
//! for [`Ex::series_dir`] and at `±∞`), in the same sense as
//! [`Ex::limit`]: every term of exponent `< order` is present and exact, and
//! the remainder is `O((var − point)^order)` along the reals.  It is not a
//! Laurent series over the complex plane: `exp(−1/x²)` expands to `0` at
//! `0` (SymPy: `O(x**4)`), `exp(−1/x)` has no two-sided expansion there
//! (`0` from the right, none from the left).  One-sided expansions from the
//! right and those at `±∞` may contain logarithms (`ln Γ(x)` at `∞` is
//! Stirling's series, `Ei(x)` at `0⁺` is `γ + ln x + x + …`).  Where no
//! such expansion exists or it cannot be established, a formal `Series`
//! node is returned (`try_*` variants: `Err`).  Use
//! [`Ex::residue`] for residues: a series is not a Laurent expansion at an
//! essential singularity.

use tracing::debug_span;

use crate::api::expr::{Ex, Expr, Numeric};
use crate::base::errors::SymplexError;
use crate::base::node::ExprNode;
use crate::calculus::limit::Direction;
use crate::calculus::summation::{self, SumOutcome};

impl Expr<Numeric> {
    // ── Summation ──────────────────────────────────────────────────

    /// Evaluate `Σ_{var=lower}^{upper} self` symbolically.
    ///
    /// Bounds may be concrete integers, symbolic expressions, or
    /// `ctx.infinity()` / `ctx.neg_infinity()`.  Strategies include exact
    /// enumeration for small concrete ranges, Faulhaber's formula (any
    /// degree, Bernoulli numbers), partial-fraction telescoping, harmonic
    /// numbers, geometric and arithmetico-geometric series, binomial
    /// identities (`Σ P(k)·C(n,k)·xᵏ` for any polynomial `P`), Gosper's
    /// algorithm, and — for infinite sums — p-series (`ζ(2m)` in closed
    /// form, `ζ(2m+1)` as a `zeta` node, `Σ (−1)ᵏ/(2k+1)² = G`), alternating
    /// series, and a table of classical power series (`Σ xᵏ/k! = eˣ`, `sin`,
    /// `cos`, `atan`, …).
    ///
    /// Convergence of a geometric series with a *symbolic* ratio cannot be
    /// decided (there is no `|r| < 1` assumption), so `Σ_{k≥0} rᵏ` with
    /// symbolic `r` stays unevaluated; numeric ratios are decided exactly.
    ///
    /// When no closed form is known the formal `Sum` node is returned.
    /// Sums that provably diverge to `±∞` evaluate to `oo` / `-oo`; for
    /// oscillating divergence the formal `Sum` is kept.  Use
    /// [`try_summation`](Self::try_summation) to get an error instead.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let k = ctx.symbol("k");
    /// let n = ctx.symbol("n");
    ///
    /// // Σ_{k=1}^{n} k² = n³/3 + n²/2 + n/6
    /// let s = k.powi(2).summation(&k, &ctx.int(1), &n);
    /// assert_eq!(s.subs_i64(&n, 10).eval().to_string(), "385");
    ///
    /// // Σ_{k=1}^{∞} 1/k² = π²/6
    /// let basel = k.powi(-2).summation(&k, &ctx.int(1), &ctx.infinity());
    /// assert_eq!(basel.to_string(), "1/6*pi^2");
    ///
    /// // Σ_{k=0}^{∞} x^k/k! = e^x
    /// let x = ctx.symbol("x");
    /// let e = (x.pow(&k) / k.factorial()).summation(&k, &ctx.int(0), &ctx.infinity());
    /// assert_eq!(e.to_string(), "exp(x)");
    ///
    /// // Σ_{k=1}^{∞} 1/k³ = ζ(3)
    /// let apery = k.powi(-3).summation(&k, &ctx.int(1), &ctx.infinity());
    /// assert_eq!(apery.to_string(), "zeta(3)");
    /// ```
    #[must_use = "returns the evaluated sum; does not modify in place"]
    pub fn summation(&self, var: &Ex, lower: &Ex, upper: &Ex) -> Ex {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lower);
        let hi_id = self.checked_id(upper);
        let _span = debug_span!("summation", body = ?self.raw_id()).entered();
        let mut inner = self.inner.write();
        let outcome = summation::summation(&mut inner.arena, self.raw_id(), var_id, lo_id, hi_id);
        let id = match outcome {
            SumOutcome::Closed(id) => id,
            SumOutcome::Divergent(Some(inf)) => inf,
            SumOutcome::Divergent(None) | SumOutcome::Unevaluated => inner
                .arena
                .intern(ExprNode::Sum(self.raw_id(), var_id, lo_id, hi_id)),
        };
        drop(inner);
        self.wrap(id)
    }

    /// Like [`summation`](Self::summation), but returns `Err` when the sum
    /// diverges ([`SymplexError::Divergent`]) or no closed form was found
    /// ([`SymplexError::ComputationFailed`]).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let k = ctx.symbol("k");
    /// let harmonic = k.powi(-1).try_summation(&k, &ctx.int(1), &ctx.infinity());
    /// assert!(matches!(harmonic, Err(SymplexError::Divergent { .. })));
    ///
    /// let geometric = ctx.rational(1, 2).pow(&k).try_summation(&k, &ctx.int(0), &ctx.infinity());
    /// assert_eq!(geometric.unwrap().to_string(), "2");
    /// ```
    pub fn try_summation(&self, var: &Ex, lower: &Ex, upper: &Ex) -> Result<Ex, SymplexError> {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lower);
        let hi_id = self.checked_id(upper);
        let _span = debug_span!("try_summation", body = ?self.raw_id()).entered();
        let mut inner = self.inner.write();
        let outcome = summation::summation(&mut inner.arena, self.raw_id(), var_id, lo_id, hi_id);
        drop(inner);
        match outcome {
            SumOutcome::Closed(id) => {
                let result = self.wrap(id);
                if result.has_unevaluated() {
                    Err(SymplexError::ComputationFailed {
                        operation: "summation",
                        reason: "no closed form for part of the sum".into(),
                    })
                } else {
                    Ok(result)
                }
            }
            SumOutcome::Divergent(Some(inf)) => Err(SymplexError::Divergent {
                operation: "summation",
                reason: format!("the series diverges to {}", self.wrap(inf)),
            }),
            SumOutcome::Divergent(None) => Err(SymplexError::Divergent {
                operation: "summation",
                reason: "the series diverges (terms do not tend to zero)".into(),
            }),
            SumOutcome::Unevaluated => Err(SymplexError::ComputationFailed {
                operation: "summation",
                reason: "no closed form found".into(),
            }),
        }
    }

    // ── Products ───────────────────────────────────────────────────

    /// Evaluate `Π_{var=lower}^{upper} self` symbolically.
    ///
    /// Handles constant factors (`cⁿ`), `Π k = n!`, `Π (k+a)` as Gamma /
    /// factorial ratios, `Π aᶠ⁽ᵏ⁾ = a^{Σ f(k)}`, products of factors, and
    /// rational functions whose numerator and denominator factor into
    /// linear factors over ℚ (`Π (1 − 1/k²) = (n+1)/(2n)`).  Infinite
    /// products are evaluated through the limit of the finite closed form
    /// when that limit is exactly computable.
    ///
    /// Returns the formal `Product` node when no closed form is known.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let k = ctx.symbol("k");
    /// let n = ctx.symbol("n");
    ///
    /// assert_eq!(k.product_over(&k, &ctx.int(1), &n).to_string(), "n!");
    ///
    /// // Π_{k=1}^{n} (1 + 1/k) = n + 1
    /// let p = (&ctx.int(1) + &k.powi(-1)).product_over(&k, &ctx.int(1), &n);
    /// assert_eq!(p.to_string(), "n + 1");
    ///
    /// // Π_{k=2}^{∞} (1 − 1/k²) = 1/2
    /// let p = (&ctx.int(1) - &k.powi(-2)).product_over(&k, &ctx.int(2), &ctx.infinity());
    /// assert_eq!(p.to_string(), "1/2");
    /// ```
    #[must_use = "returns the evaluated product; does not modify in place"]
    pub fn product_over(&self, var: &Ex, lower: &Ex, upper: &Ex) -> Ex {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lower);
        let hi_id = self.checked_id(upper);
        let _span = debug_span!("product_over", body = ?self.raw_id()).entered();
        let mut inner = self.inner.write();
        let outcome = summation::product(&mut inner.arena, self.raw_id(), var_id, lo_id, hi_id);
        let id = match outcome {
            SumOutcome::Closed(id) => id,
            SumOutcome::Divergent(Some(inf)) => inf,
            SumOutcome::Divergent(None) | SumOutcome::Unevaluated => inner
                .arena
                .intern(ExprNode::Product_(self.raw_id(), var_id, lo_id, hi_id)),
        };
        drop(inner);
        self.wrap(id)
    }

    /// Like [`product_over`](Self::product_over), but returns `Err` when the
    /// product diverges or no closed form was found.
    pub fn try_product_over(&self, var: &Ex, lower: &Ex, upper: &Ex) -> Result<Ex, SymplexError> {
        let var_id = self.checked_id(var);
        let lo_id = self.checked_id(lower);
        let hi_id = self.checked_id(upper);
        let mut inner = self.inner.write();
        let outcome = summation::product(&mut inner.arena, self.raw_id(), var_id, lo_id, hi_id);
        drop(inner);
        match outcome {
            SumOutcome::Closed(id) => {
                let result = self.wrap(id);
                if result.has_unevaluated() {
                    Err(SymplexError::ComputationFailed {
                        operation: "product",
                        reason: "no closed form for part of the product".into(),
                    })
                } else {
                    Ok(result)
                }
            }
            SumOutcome::Divergent(_) => Err(SymplexError::Divergent {
                operation: "product",
                reason: "the infinite product diverges".into(),
            }),
            SumOutcome::Unevaluated => Err(SymplexError::ComputationFailed {
                operation: "product",
                reason: "no closed form found".into(),
            }),
        }
    }

    // ── Hypergeometric terms ───────────────────────────────────────

    /// If `self` is a hypergeometric term in `var`, return the ratio
    /// `self(var+1) / self(var)` as a rational function of `var`.
    ///
    /// Returns `None` when the ratio is not a rational function of `var`
    /// (e.g. for `sin(k)` or `k^k`).  Symbolic parameters other than `var`
    /// are allowed.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let k = ctx.symbol("k");
    /// let n = ctx.symbol("n");
    ///
    /// let r = k.factorial().hypergeometric_ratio(&k).unwrap();
    /// assert_eq!(r.to_string(), "k + 1");
    ///
    /// let r = n.binomial(&k).hypergeometric_ratio(&k).unwrap();
    /// assert_eq!(r.subs_i64(&n, 5).subs_i64(&k, 2).eval().to_string(), "1");
    ///
    /// assert!(k.sin().hypergeometric_ratio(&k).is_none());
    /// ```
    #[must_use]
    pub fn hypergeometric_ratio(&self, var: &Ex) -> Option<Ex> {
        let var_id = self.checked_id(var);
        let mut inner = self.inner.write();
        let r = crate::calculus::gosper::is_hypergeometric(&mut inner.arena, self.raw_id(), var_id);
        drop(inner);
        r.map(|id| self.wrap(id))
    }

    // ── Convergence ────────────────────────────────────────────────

    /// Test whether `Σ |self(var)|` converges (absolute convergence).
    ///
    /// Sign-alternating factors such as `(−1)^k` are stripped before the
    /// convergence tests run, so conditionally convergent series like
    /// `Σ (−1)^k/k` return `Some(false)` here but `Some(true)` from
    /// [`is_convergent`](Self::is_convergent).
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let k = ctx.symbol("k");
    /// let alt = ctx.int(-1).pow(&k) / &k;
    /// assert_eq!(alt.is_convergent(&k), Some(true));
    /// assert_eq!(alt.is_absolutely_convergent(&k), Some(false));
    /// ```
    #[must_use]
    pub fn is_absolutely_convergent(&self, var: &Ex) -> Option<bool> {
        let var_id = self.checked_id(var);
        let _span = debug_span!("is_absolutely_convergent").entered();
        let mut inner = self.inner.write();
        crate::calculus::convergence::is_absolutely_convergent(
            &mut inner.arena,
            self.raw_id(),
            var_id,
        )
    }

    // ── Asymptotic expansions ──────────────────────────────────────

    /// Asymptotic expansion of `self` as `var → +∞`, with every term of
    /// exponent `< n_terms` in `1/var` (what a series is: see
    /// [`series_dir`](Self::series_dir)).
    ///
    /// Internally substitutes `var = 1/t`, expands at `t → 0⁺`, and
    /// substitutes back; coefficients may contain `ln(var)` (Stirling's
    /// series `ln Γ(x) = x ln x − x − ½ ln x + ½ ln 2π + 1/(12x) + …`).
    /// Returns a formal `Series` node when the expansion cannot be
    /// computed.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    /// // x/(x+1) = 1 − 1/x + 1/x² − …
    /// let s = (&x / &(&x + 1)).series_at_infinity(&x, 3);
    /// let at_10 = s.subs_i64(&x, 10).eval_f64().unwrap();
    /// assert!((at_10 - 0.91).abs() < 1e-12);
    /// ```
    #[must_use = "returns the expansion; does not modify in place"]
    pub fn series_at_infinity(&self, var: &Ex, n_terms: u32) -> Ex {
        self.series_at_infinity_impl(var, n_terms, false)
    }

    /// Asymptotic expansion of `self` as `var → −∞` (see
    /// [`series_at_infinity`](Self::series_at_infinity)).
    #[must_use = "returns the expansion; does not modify in place"]
    pub fn series_at_neg_infinity(&self, var: &Ex, n_terms: u32) -> Ex {
        self.series_at_infinity_impl(var, n_terms, true)
    }

    /// Like [`series_at_infinity`](Self::series_at_infinity), but returns
    /// `Err` if the expansion could not be computed.
    pub fn try_series_at_infinity(&self, var: &Ex, n_terms: u32) -> Result<Ex, SymplexError> {
        let r = self.series_at_infinity(var, n_terms);
        if r.has_unevaluated() {
            Err(SymplexError::ComputationFailed {
                operation: "series_at_infinity",
                reason: "could not compute asymptotic expansion".into(),
            })
        } else {
            Ok(r)
        }
    }

    /// The expansion of `self` as `var → point` from the side `dir`, with
    /// every term of exponent `< order` in `var − point`.
    ///
    /// Like [`series`](Self::series) (which is `series_dir` with
    /// [`Direction::Both`]) this is the **asymptotic expansion as
    /// `var → point` through real values**, in the same sense as
    /// [`limit_dir`](Self::limit_dir): every term of exponent `< order` is
    /// present and exact, and the remainder is `O((var − point)^order)` as
    /// `var` approaches `point` from the side(s) in question.  It is not a
    /// Laurent series over the complex plane: `exp(−1/x²)` expands to `0` at
    /// `0` (SymPy: `O(x**4)`), `exp(−1/x)` has no two-sided expansion there
    /// (`0` from the right, none from the left).  From either side, and at
    /// `±∞` (where the side is implied and `dir` is ignored), the expansion
    /// may contain logarithms — `ln(var − point)` from the right,
    /// `ln(point − var)` from the left: `Ei(x) = γ + ln x + x + …` for
    /// `x → 0⁺`, `x^x = 1 + x·(ln(−x) + iπ) + …` for `x → 0⁻`, Stirling's
    /// series for `ln Γ(x)` at `∞`.  At a branch point a one-sided
    /// expansion may be a Puiseux series, in powers of `√(var − point)` or
    /// `∛(var − point)` (of `√(point − var)` from the left), as SymPy gives:
    /// `asin(x²)` at `x → 1⁺` is `π/2 − 2i·√(x − 1) − …`.  Where no such
    /// expansion exists or it cannot be established, a formal `Series` node
    /// is returned.  For residues use [`residue`](Self::residue): at an
    /// essential singularity a series is not the Laurent expansion.
    ///
    /// # Examples
    ///
    /// ```
    /// use symplex::prelude::*;
    ///
    /// let ctx = Context::new();
    /// let x = ctx.symbol("x");
    /// let f = (-(1 / &x)).exp(); // exp(−1/x)
    /// let zero = ctx.int(0);
    /// // O(x⁴) from the right, no expansion from the left or two-sided
    /// assert_eq!(f.series_dir(&x, &zero, 4, Direction::Right).to_string(), "0");
    /// assert!(f.series_dir(&x, &zero, 4, Direction::Left).has_unevaluated());
    /// assert!(f.series(&x, &zero, 4).has_unevaluated());
    /// // Ei(x) = γ + ln x + x + x²/4 + … for x → 0⁺
    /// let ei = x.ei().series_dir(&x, &zero, 3, Direction::Right);
    /// assert!(ei.to_string().contains("ln(x)"), "{ei}");
    /// // A branch point: asin(x²) for x → 1⁻ is π/2 − 2√(1 − x) + …
    /// let one = ctx.int(1);
    /// let a = x.powi(2).asin();
    /// assert!(a.series(&x, &one, 2).has_unevaluated()); // no two-sided one
    /// let left = a.series_dir(&x, &one, 1, Direction::Left);
    /// assert_eq!(left.to_string(), "-2*sqrt(-x + 1) + 1/2*pi");
    /// ```
    #[must_use = "returns the expansion; does not modify in place"]
    pub fn series_dir(&self, var: &Ex, point: &Ex, order: u32, dir: Direction) -> Ex {
        let var_id = self.checked_id(var);
        let point_id = self.checked_id(point);
        let _span = debug_span!("series_dir", expr = ?self.raw_id(), order = order).entered();
        let mut inner = self.inner.write();
        let id = match crate::calculus::series::series_dir(
            &mut inner.arena,
            self.raw_id(),
            var_id,
            point_id,
            order,
            dir,
        ) {
            Ok(id) => id,
            Err(_) => {
                let order_id = inner.arena.int(order as i64);
                inner
                    .arena
                    .intern(ExprNode::Series(self.raw_id(), var_id, point_id, order_id))
            }
        };
        drop(inner);
        self.wrap(id)
    }

    /// Like [`series_dir`](Self::series_dir), but returns `Err` if there
    /// is no expansion.
    pub fn try_series_dir(
        &self,
        var: &Ex,
        point: &Ex,
        order: u32,
        dir: Direction,
    ) -> Result<Ex, SymplexError> {
        let r = self.series_dir(var, point, order, dir);
        if r.has_unevaluated() {
            Err(SymplexError::ComputationFailed {
                operation: "series_dir",
                reason: "no expansion from this side".into(),
            })
        } else {
            Ok(r)
        }
    }

    fn series_at_infinity_impl(&self, var: &Ex, n_terms: u32, negative: bool) -> Ex {
        let var_id = self.checked_id(var);
        let _span = debug_span!("series_at_infinity", expr = ?self.raw_id()).entered();
        let mut inner = self.inner.write();
        let id = crate::calculus::series::series_at_infinity(
            &mut inner.arena,
            self.raw_id(),
            var_id,
            n_terms,
            negative,
        );
        let id = match id {
            Ok(id) => id,
            Err(_) => {
                let point = if negative {
                    inner.arena.neg_infinity
                } else {
                    inner.arena.infinity
                };
                let order = inner.arena.int(n_terms as i64);
                inner
                    .arena
                    .intern(ExprNode::Series(self.raw_id(), var_id, point, order))
            }
        };
        drop(inner);
        self.wrap(id)
    }
}