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)")]