Integrate the polygamma functions, and give an odd order its exact values

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
 

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));
             }
 
+            /// <summary>
+            /// <c>psi^(n)(q)</c> exactly, for an odd order <c>n</c> and a whole or half-whole <c>q</c> other than a pole,
+            /// or <see langword="null"/>: a rational times <c>pi^(n + 1)</c>, plus a rational.
+            /// </summary>
+            /// <remarks>
+            /// <c>psi^(n)(1) = n! zeta(n + 1)</c> and <c>psi^(n)(1/2) = (2^(n + 1) - 1) n! zeta(n + 1)</c> for an odd
+            /// <c>n</c> (DLMF 5.15.2, 5.15.3), where <c>zeta(2m) = (-1)^(m + 1) B_(2m) (2 pi)^(2m)/(2 (2m)!)</c> is a
+            /// rational times <c>pi^(2m)</c>; the recurrence <c>psi^(n)(z + 1) = psi^(n)(z) - n!/z^(n + 1)</c> moves
+            /// the argument to either. An even order wants <c>zeta</c> at an odd number, or Euler's constant, which
+            /// have no closed form here. https://github.com/asc-community/AngouriMath/issues/171
+            /// </remarks>
+            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)

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(...)`;

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),

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

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);
         }
 
+        /// <summary>
+        /// An odd order at a whole or half-whole argument is a rational times <c>pi^(n + 1)</c> plus a rational
+        /// (DLMF 5.15.2, 5.15.3); each is checked against its numerical value.
+        /// </summary>
+        [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<Real>(simplified.EvalNumerical()).EDecimal);
+        }
+
+        /// <summary>An even order wants zeta at an odd number, or Euler's constant: it stays as written.</summary>
+        [Theory]
+        [InlineData("digamma(1)")]
+        [InlineData("polygamma(2, 1)")]
+        [InlineData("polygamma(1, 1/3)")]
+        public void OtherwiseItStaysAsWritten(string expression)
+            => Assert.IsType<Entity.Polygammaf>(expression.ToEntity().InnerSimplified);
+
+        /// <summary>Each order integrates to the one before, and the digamma function to <c>ln Gamma</c>.</summary>
+        [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<Real>(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)")]