From 2e088bd95aa8d4e1d605f5b948e762b657422aa8 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 22 Oct 2026 17:48:57 +0000 Subject: [PATCH] Integrate the polygamma functions, and give an odd order its exact values Each order integrates to the one before, int psi^(n)(u) = psi^(n - 1)(u), and the digamma function to ln Gamma(u) = ln((u - 1)!), for an argument linear in the variable; the integrator reads them from the table it reads the other special functions from. An odd order at a whole or half-whole argument is exact: psi^(n)(1) = n! zeta(n + 1) and psi^(n)(1/2) = (2^(n + 1) - 1) n! zeta(n + 1) (DLMF 5.15.2, 5.15.3), with zeta(2m) a rational times pi^(2m) from the exact Bernoulli numbers, and the recurrence moves the argument to either: polygamma(1, 1) is pi^2/6 and polygamma(1, 3) is pi^2/6 - 5/4. An even order wants zeta at an odd number or Euler's constant, and stays as written. Each exact value is checked against the numerical one. Part of #171. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 4 ++ .../Entity/Continuous/Number/Polygamma.cs | 50 +++++++++++++++++++ Sources/AngouriMath/Docs/Usage/Syntax.md | 7 ++- .../Integration/IntegralPatterns.cs | 13 +++++ ...aluation.Continuous.Arithmetics.Classes.cs | 1 + .../UnitTests/Convenience/PolygammaTest.cs | 41 +++++++++++++++ 6 files changed, 114 insertions(+), 2 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 1b777ab77..e2290d4e3 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -103,6 +103,8 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa - **New functions.** `polygamma(n, z)` and `digamma(z)` are the polygamma functions `psi^(n)(z)`, the `(n + 1)`-th derivatives of `ln Gamma(z)`, and the digamma function, which is `polygamma(0, z)`. Neither name was recognised, so each was read as a product with an undeclared variable. - **Evaluation.** Each is evaluated to the working precision anywhere in the complex plane but at `0, -1, -2, ...`, its poles, where it is `NaN`. - **Derivatives.** Each order differentiates to the next. The derivative of `x!`, which was left unevaluated, is `x! * digamma(x + 1)` ([#171](https://github.com/asc-community/AngouriMath/issues/171)). +- **Integrals.** Each order integrates to the one before, and `digamma(x)` to `ln((x - 1)!)`. +- **Exact values.** An odd order at a whole or half-whole argument is a rational times a power of `pi` plus a rational: `polygamma(1, 1)` is `1/6 * pi ^ 2`. | Input | Was (2.6.0) | Now | |---|---|---| @@ -111,6 +113,8 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa | `"polygamma(1, 1)".ToEntity()` | an `UnhandledParseException` | `polygamma(1, 1)`, which evaluates to `1.6449340668482264…` | | `"x!".ToEntity().Differentiate("x")` | `derivative(x!, x)`, unevaluated | `x! * digamma(x + 1)` | | `"gamma(x)".ToEntity().Differentiate("x")` | `derivative((x - 1)!, x)` | `(x - 1)! * digamma(x - 1 + 1)` | +| `"polygamma(1, 3)".ToEntity().InnerSimplified` | an `UnhandledParseException` | `1/6 * pi ^ 2 - 5/4` | +| `"digamma(x)".ToEntity().Integrate("x")` | `digamma * x ^ 2 / 2 + C`, of the product | `ln((x - 1)!) + C` | ### Numerical equality is to half the evaluation's digits diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/Polygamma.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/Polygamma.cs index c505ed5c3..e0421ef1d 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/Polygamma.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/Polygamma.cs @@ -59,6 +59,56 @@ public static Complex Polygamma(EInteger order, Complex z) return Complex.Create(valueRe.RoundToPrecision(context), valueIm.RoundToPrecision(context)); } + /// + /// psi^(n)(q) exactly, for an odd order n and a whole or half-whole q other than a pole, + /// or : a rational times pi^(n + 1), plus a rational. + /// + /// + /// psi^(n)(1) = n! zeta(n + 1) and psi^(n)(1/2) = (2^(n + 1) - 1) n! zeta(n + 1) for an odd + /// n (DLMF 5.15.2, 5.15.3), where zeta(2m) = (-1)^(m + 1) B_(2m) (2 pi)^(2m)/(2 (2m)!) is a + /// rational times pi^(2m); the recurrence psi^(n)(z + 1) = psi^(n)(z) - n!/z^(n + 1) moves + /// the argument to either. An even order wants zeta at an odd number, or Euler's constant, which + /// have no closed form here. https://github.com/asc-community/AngouriMath/issues/171 + /// + internal static Entity? ExactPolygamma(EInteger order, ERational argument) + { + if (order.Sign <= 0 || order.IsEven || order.CompareTo(100) > 0) + return null; + var half = ERational.Create(EInteger.One, EInteger.FromInt32(2)); + var start = argument.Denominator.Equals(EInteger.One) ? ERational.One + : argument.Denominator.Equals(EInteger.FromInt32(2)) ? half : null; + if (start is null || argument.Denominator.Equals(EInteger.One) && argument.Sign <= 0) + return null; + var steps = argument.Subtract(start).ToEInteger(); + if (steps.Abs().CompareTo(1000) > 0) + return null; + var n = order.ToInt32Checked(); + var m = (n + 1) / 2; + var factorial = Factorial(n); + // zeta(2m) / pi^(2m) = (-1)^(m + 1) B_(2m) 2^(2m) / (2 (2m)!) + var zetaOverPiPower = Bernoulli(m).Multiply(ERational.FromEInteger(EInteger.One.ShiftLeft(2 * m))) + .Divide(ERational.FromEInteger(Factorial(2 * m).Multiply(2))); + if (m % 2 == 0) + zetaOverPiPower = zetaOverPiPower.Negate(); + var piCoefficient = zetaOverPiPower.Multiply(ERational.FromEInteger(factorial)); + if (start.Equals(half)) + piCoefficient = piCoefficient.Multiply(ERational.FromEInteger(EInteger.One.ShiftLeft(n + 1).Subtract(EInteger.One))); + // The recurrence, up from the start or down to the argument. + var rest = ERational.Zero; + var power = n + 1; + if (steps.Sign > 0) + for (var j = EInteger.Zero; j.CompareTo(steps) < 0; j = j.Add(EInteger.One)) + rest = rest.Subtract(ERational.FromEInteger(factorial).Divide(Pow(start.Add(ERational.FromEInteger(j)), power))); + else + for (var j = EInteger.One; j.CompareTo(steps.Negate()) <= 0; j = j.Add(EInteger.One)) + rest = rest.Add(ERational.FromEInteger(factorial).Divide(Pow(start.Subtract(ERational.FromEInteger(j)), power))); + var piTerm = Rational.Create(piCoefficient) * MathS.Pow(MathS.pi, Integer.Create(power)); + return (rest.Sign < 0 ? piTerm - Rational.Create(rest.Negate()) : piTerm + Rational.Create(rest)).InnerSimplified; + + static ERational Pow(ERational x, int exponent) + => ERational.Create(x.Numerator.Pow(exponent), x.Denominator.Pow(exponent)); + } + private static (EDecimal Re, EDecimal Im) Polygamma(int n, EDecimal re, EDecimal im, EContext work) { if (re.CompareTo(EDecimal.FromString("0.5")) < 0) diff --git a/Sources/AngouriMath/Docs/Usage/Syntax.md b/Sources/AngouriMath/Docs/Usage/Syntax.md index 7ae45cc54..85aff5bd8 100644 --- a/Sources/AngouriMath/Docs/Usage/Syntax.md +++ b/Sources/AngouriMath/Docs/Usage/Syntax.md @@ -332,8 +332,11 @@ undefined at `0`; `Si(+oo)` is `pi / 2` and `Ci(+oo)` is `0`. to the working precision anywhere in the complex plane but at `0, -1, -2, ...`, its poles, and each order is the derivative of the one before: the derivative of `polygamma(n, x)` is `polygamma(n + 1, x)`. The derivative of `x!` is `x! * digamma(x + 1)`, and of `gamma(x)` it is -`(x - 1)! * digamma(x)`. Special values, such as `digamma(1) = -gamma`, stay as written, since the -library has no symbol for Euler's constant. +`(x - 1)! * digamma(x)`. Integrating goes the other way: `polygamma(n, x)` integrates to +`polygamma(n - 1, x)`, and `digamma(x)` to `ln((x - 1)!)`. An odd order at a whole or half-whole +argument is exact, a rational times a power of `pi` plus a rational: `polygamma(1, 1)` is +`1/6 * pi ^ 2` and `polygamma(1, 3)` is `1/6 * pi ^ 2 - 5/4`. Other special values, such as +`digamma(1) = -gamma`, stay as written, since the library has no symbol for Euler's constant. **Calculus** — `derivative(expr, var, order)`, `integral(expr, var)`, `integral(expr, var, from, to)`, `limit(expr, var, dest)`, `limitleft(...)`, `limitright(...)`; diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index d3ff34dfb..4d845c94b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -67,6 +67,7 @@ internal static Entity AntiderivativeLog(Entity arg) Entity.Cif(var arg) => arg, Entity.Shif(var arg) => arg, Entity.Chif(var arg) => arg, + Entity.Polygammaf(var order, var arg) when !order.ContainsNode(x) => arg, Entity.Arcsinf(var arg) => arg, Entity.Arccosf(var arg) => arg, Entity.Arctanf(var arg) => arg, @@ -384,6 +385,18 @@ _ when GaussianMoment(expr, x) is { } moment => moment, TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => (arg * MathS.Chi(arg) - MathS.Hyperbolic.Sinh(arg)) / a, + // The polygamma functions, each the derivative of the one before: int psi^(n)(u) = + // psi^(n - 1)(u) for n >= 1, and int psi(u) = ln Gamma(u) = ln((u - 1)!), whose value is + // i pi more where Gamma is negative, a constant on each interval between the poles. + // https://github.com/asc-community/AngouriMath/issues/171 + Entity.Polygammaf(Entity.Number.Integer { EInteger: var order }, var arg) when order.Sign > 0 && + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + MathS.Polygamma(Entity.Number.Integer.Create(order - 1), arg) / a, + + Entity.Polygammaf(Entity.Number.Integer { IsZero: true }, var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + MathS.Ln(MathS.Factorial(arg - 1)) / a, + Entity.Absf(var arg) when TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => // ∫ |ax + b| dx = sgn(ax + b) * (ax + b)^2 / (2a) MathS.Signum(arg) * MathS.Pow(arg, 2) / (2 * a), diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs index 797005a64..e0245a716 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs @@ -602,6 +602,7 @@ protected override Entity InnerSimplify(bool isExact) => ExpandOnTwoArguments(Order, Argument, (order, argument) => (order, argument) switch { + (Integer { EInteger: var n }, Rational q) when isExact && Number.ExactPolygamma(n, q.ERational) is { } exact => exact, (Integer { EInteger: var n }, Complex z) when !isExact && n.Sign >= 0 => Number.Polygamma(n, z), (Integer { EInteger: var n }, Real r) when n.Sign >= 0 && r.EDecimal.IsPositiveInfinity() => n.IsZero ? Real.PositiveInfinity : Integer.Zero, _ => null diff --git a/Sources/Tests/UnitTests/Convenience/PolygammaTest.cs b/Sources/Tests/UnitTests/Convenience/PolygammaTest.cs index 85f5938bf..9e3322c86 100644 --- a/Sources/Tests/UnitTests/Convenience/PolygammaTest.cs +++ b/Sources/Tests/UnitTests/Convenience/PolygammaTest.cs @@ -78,6 +78,47 @@ public void EachOrderIsTheDerivativeOfTheOneBefore() Assert.Equal((MathS.Polygamma(3, "x^2") * (2 * MathS.Var("x"))).InnerSimplified, "polygamma(2, x^2)".ToEntity().Differentiate("x").InnerSimplified); } + /// + /// An odd order at a whole or half-whole argument is a rational times pi^(n + 1) plus a rational + /// (DLMF 5.15.2, 5.15.3); each is checked against its numerical value. + /// + [Theory] + [InlineData("polygamma(1, 1)", "1/6 * pi ^ 2")] + [InlineData("polygamma(1, 1/2)", "1/2 * pi ^ 2")] + [InlineData("polygamma(3, 1)", "1/15 * pi ^ 4")] + [InlineData("polygamma(1, 3)", "1/6 * pi ^ 2 - 5/4")] + [InlineData("polygamma(1, -1/2)", "1/2 * pi ^ 2 + 4")] + [InlineData("polygamma(3, 5/2)", "pi ^ 4 - 2624/27")] + public void AnOddOrderAtAWholeOrHalfWholeArgumentIsExact(string expression, string exact) + { + var simplified = expression.ToEntity().InnerSimplified; + Assert.Equal(exact.ToEntity().InnerSimplified, simplified); + AssertClose(expression.EvalNumerical().Stringize(), Assert.IsType(simplified.EvalNumerical()).EDecimal); + } + + /// An even order wants zeta at an odd number, or Euler's constant: it stays as written. + [Theory] + [InlineData("digamma(1)")] + [InlineData("polygamma(2, 1)")] + [InlineData("polygamma(1, 1/3)")] + public void OtherwiseItStaysAsWritten(string expression) + => Assert.IsType(expression.ToEntity().InnerSimplified); + + /// Each order integrates to the one before, and the digamma function to ln Gamma. + [Theory] + [InlineData("polygamma(2, x)", "polygamma(1, x)")] + [InlineData("polygamma(1, 3 * x + 2)", "digamma(3 * x + 2) / 3")] + [InlineData("digamma(x)", "ln((x - 1)!)")] + [InlineData("digamma(2 * x)", "ln((2 * x - 1)!) / 2")] + public void IntegratesToTheOrderBefore(string integrand, string antiderivative) + { + var integral = integrand.ToEntity().Integrate("x").Substitute("C", 0).Simplify(); + Assert.Equal(antiderivative.ToEntity().Simplify(), integral); + var back = integral.Differentiate("x").Substitute("x", Rational.Create(7, 3)).EvalNumerical(); + var expected = integrand.ToEntity().Substitute("x", Rational.Create(7, 3)).EvalNumerical(); + AssertClose(expected.Stringize(), Assert.IsType(back).EDecimal); + } + [Theory] [InlineData("digamma(x)", "digamma(x)", @"\psi\left(x\right)", "sympy.digamma(x)")] [InlineData("polygamma(2, x)", "polygamma(2, x)", @"\psi^{(2)}\left(x\right)", "sympy.polygamma(2, x)")]