From d1186d2b8892124c91749e80408897780d6c4bec Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 30 Oct 2026 00:40:51 +0000 Subject: [PATCH] The polygamma functions and the harmonic numbers compile polygamma, digamma and harmonicnumber compile in both compilers, through PolygammaFunction, a double-precision kernel by the interpreter's algorithm: reflected left of Re z = 1/2 near the axis, moved right by the recurrence to Re w >= 16 + n, and summed by the asymptotic series with the Bernoulli numbers from the tangent numbers. Off the axis the recurrence is taken without the reflection, whose cotangent polynomial cancels there to all but a few digits, and cot is taken with its argument reduced, since at a half-whole z its rounding is magnified by pi^(n + 1) and the polynomial's coefficients. The harmonic numbers are the sum for a whole n and the polygamma form elsewhere, as the interpreter has them. Both compilers agree with evaluation on a grid of orders and arguments, and both, with evaluation, agree with mpmath at 7,170 polygamma points and 2,997 harmonic ones. Part of #2195 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 4 +- Sources/.editorconfig | 6 + .../Functions/Compilation/IntoFE/Compiler.cs | 25 ++ .../Compilation/IntoFE/FastExpression.cs | 8 + .../IntoLinq/CompilationProtocol.cs | 4 + .../Compilation/IntoLinq/MathAllMethods.cs | 11 + .../Compilation/IntoLinq/MathAllMethods.tt | 9 + .../AngouriMath/Numerics/PolygammaFunction.cs | 229 ++++++++++++++++++ Sources/Tests/AotSmokeTest/Program.cs | 2 + .../Common/PolygammaCompilationTest.cs | 69 ++++++ .../UnitTests/Core/UserInvalidExceptions.cs | 5 +- 11 files changed, 369 insertions(+), 3 deletions(-) create mode 100644 Sources/AngouriMath/Numerics/PolygammaFunction.cs create mode 100644 Sources/Tests/UnitTests/Common/PolygammaCompilationTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 023f26be6..423f1e962 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -654,7 +654,8 @@ other a parse error.** They are `1 + 1/2^m + ... + 1/n^m` for a whole `n >= 0`, `ζ(m) - ζ(m, n + 1)` elsewhere, as Mathematica's `HarmonicNumber` is: `ψ(n + 1) - ψ(1)` for `m = 1`, evaluated numerically through `digamma` and `polygamma` for a whole `m >= 1`, and exactly at a half-whole `n` for `m = 1` or an even `m`. A negative integer `n` has no value but for a whole `m <= 0` -([#2182](https://github.com/asc-community/AngouriMath/issues/2182)). +([#2182](https://github.com/asc-community/AngouriMath/issues/2182)). They compile, on the same terms, in both compilers +([#2195](https://github.com/asc-community/AngouriMath/issues/2195)). | Input | Was (2.6.0) | Now | |---|---|---| @@ -866,6 +867,7 @@ no others. Where that proof does not go through, or a bound is symbolic, the int - **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`. The digamma function there is a rational less Euler's constant, and at a half-whole argument less `2 ln 2` as well: `digamma(1)` is `-eulergamma`. +- **Compiled forms.** Each compiles, in double precision and by the same reflection, recurrence and asymptotic series, in both compilers ([#2195](https://github.com/asc-community/AngouriMath/issues/2195)). | Input | Was (2.6.0) | Now | |---|---|---| diff --git a/Sources/.editorconfig b/Sources/.editorconfig index a24681bd6..14650bf62 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -369,6 +369,12 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Common/LambertWCompilationTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[AngouriMath/Numerics/PolygammaFunction.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[Tests/UnitTests/Common/PolygammaCompilationTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [Tests/UnitTests/Common/CompileMatrixTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs index 1f036fcc2..2b56e8204 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs @@ -617,6 +617,31 @@ private protected override void CompileNode(Compiler compiler) } } + // The polygamma functions and the harmonic numbers. https://github.com/asc-community/AngouriMath/issues/2195 + public partial record Polygammaf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + Order.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_POLYGAMMA)); + } + } + + public partial record HarmonicNumberf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Order.InnerCompile(compiler); + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_HARMONIC_NUMBER)); + } + } + public partial record Maxf { private protected override bool CompilesItself => true; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs index a2cf38559..ba6963ad6 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs @@ -82,6 +82,8 @@ internal enum InstructionType CALL_FALLING_FACTORIAL, CALL_RISING_FACTORIAL, CALL_LAMBERTW, + CALL_POLYGAMMA, + CALL_HARMONIC_NUMBER, // 3-arg: a value where its predicate holds, and the one below it elsewhere. CALL_SELECT = 150, @@ -359,6 +361,12 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value case InstructionType.CALL_LAMBERTW: stack.Push(AngouriMath.Numerics.LambertWFunction.LambertW(stack.Pop(), stack.Pop())); break; + case InstructionType.CALL_POLYGAMMA: + stack.Push(AngouriMath.Numerics.PolygammaFunction.Polygamma(stack.Pop(), stack.Pop())); + break; + case InstructionType.CALL_HARMONIC_NUMBER: + stack.Push(AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(stack.Pop(), stack.Pop())); + break; case InstructionType.CALL_DOUBLE_FACTORIAL: stack.Push(AngouriMath.Numerics.CombinatorialFunctions.DoubleFactorial(stack.Pop())); break; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs index 61cea63ea..93bfbeac6 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs @@ -244,6 +244,10 @@ RisingFactorialf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLe => Expression.Call(GetDef("RisingFactorial", 2, newLeft.Type), newLeft, newRight), LambertWf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(right) is var newRight => Expression.Call(GetDef("LambertW", 2, newLeft.Type), newLeft, newRight), + Polygammaf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(right) is var newRight + => Expression.Call(GetDef("Polygamma", 2, newLeft.Type), newLeft, newRight), + HarmonicNumberf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(right) is var newRight + => Expression.Call(GetDef("HarmonicNumber", 2, newLeft.Type), newLeft, newRight), Maxf => Expression.Call(GetDef("Max", 2, left.Type), left, right), Minf => Expression.Call(GetDef("Min", 2, left.Type), left, right), diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs index ef23af0b7..bfa9ad7eb 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs @@ -347,6 +347,13 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a) public static System.Numerics.Complex LambertW(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.LambertWFunction.LambertW(a, b); public static double LambertW(double a, double b) => AngouriMath.Numerics.LambertWFunction.LambertW(a, b); + // The polygamma functions of a whole order, and the harmonic numbers. + // https://github.com/asc-community/AngouriMath/issues/2195 + public static System.Numerics.Complex Polygamma(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.PolygammaFunction.Polygamma(a, b); + public static double Polygamma(double a, double b) => AngouriMath.Numerics.PolygammaFunction.Polygamma(a, b); + public static System.Numerics.Complex HarmonicNumber(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); + public static double HarmonicNumber(double a, double b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); + /// /// True where System.Numerics.Complex.IsNaN is, spelled out because that overload is /// .NET 7 and later and this assembly also targets netstandard2.0. A complex number is NaN @@ -540,7 +547,11 @@ public static bool IsNaN(System.Numerics.Complex a) { ("RisingFactorial", 2, typeof(System.Numerics.Complex)), Def((a, b) => RisingFactorial(a, b)) }, { ("RisingFactorial", 2, typeof(double)), Def((a, b) => RisingFactorial(a, b)) }, { ("LambertW", 2, typeof(System.Numerics.Complex)), Def((a, b) => LambertW(a, b)) }, + { ("Polygamma", 2, typeof(System.Numerics.Complex)), Def((a, b) => Polygamma(a, b)) }, + { ("HarmonicNumber", 2, typeof(System.Numerics.Complex)), Def((a, b) => HarmonicNumber(a, b)) }, { ("LambertW", 2, typeof(double)), Def((a, b) => LambertW(a, b)) }, + { ("Polygamma", 2, typeof(double)), Def((a, b) => Polygamma(a, b)) }, + { ("HarmonicNumber", 2, typeof(double)), Def((a, b) => HarmonicNumber(a, b)) }, }; // An expression tree is how C# spells "the MethodInfo of this method" without a string: diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt index 4443f45db..d8f28e101 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt @@ -237,6 +237,13 @@ namespace AngouriMath.Core.Compilation.IntoLinq public static System.Numerics.Complex LambertW(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.LambertWFunction.LambertW(a, b); public static double LambertW(double a, double b) => AngouriMath.Numerics.LambertWFunction.LambertW(a, b); + // The polygamma functions of a whole order, and the harmonic numbers. + // https://github.com/asc-community/AngouriMath/issues/2195 + public static System.Numerics.Complex Polygamma(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.PolygammaFunction.Polygamma(a, b); + public static double Polygamma(double a, double b) => AngouriMath.Numerics.PolygammaFunction.Polygamma(a, b); + public static System.Numerics.Complex HarmonicNumber(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); + public static double HarmonicNumber(double a, double b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); + /// /// True where System.Numerics.Complex.IsNaN is, spelled out because that overload is /// .NET 7 and later and this assembly also targets netstandard2.0. A complex number is NaN @@ -293,6 +300,8 @@ namespace AngouriMath.Core.Compilation.IntoLinq <# } #> <# foreach (var type in new[]{ "System.Numerics.Complex", "double" }) { #> { ("LambertW", 2, typeof(<#= type #>)), Def<<#= type #>>((a, b) => LambertW(a, b)) }, + { ("Polygamma", 2, typeof(<#= type #>)), Def<<#= type #>>((a, b) => Polygamma(a, b)) }, + { ("HarmonicNumber", 2, typeof(<#= type #>)), Def<<#= type #>>((a, b) => HarmonicNumber(a, b)) }, <# } #> }; diff --git a/Sources/AngouriMath/Numerics/PolygammaFunction.cs b/Sources/AngouriMath/Numerics/PolygammaFunction.cs new file mode 100644 index 000000000..557171135 --- /dev/null +++ b/Sources/AngouriMath/Numerics/PolygammaFunction.cs @@ -0,0 +1,229 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using PeterO.Numbers; +using NumericsComplex = System.Numerics.Complex; + +namespace AngouriMath.Numerics +{ + /// + /// The polygamma functions and the harmonic numbers in double precision, for the compilers, on the + /// interpreter's terms and by its algorithm: psi^(n)(z) for a whole n >= 0, reflected + /// left of Re z = 1/2 (DLMF 5.15.6), moved right by the recurrence (DLMF 5.15.5) and summed by the + /// asymptotic series (DLMF 5.11.2, 5.15.8); NaN at its poles 0, -1, -2, ... and at an order + /// that is not whole. + /// https://github.com/asc-community/AngouriMath/issues/2195 + /// + internal static class PolygammaFunction + { + [AngouriMath.Core.ConstantField] private static readonly NumericsComplex NaN = new(double.NaN, 0); + + /// The most terms of the asymptotic series, and so the Bernoulli numbers kept. + private const int Terms = 40; + + /// B_2, B_4, ..., B_(2 Terms), from the tangent numbers exactly and then rounded. + [AngouriMath.Core.ConstantField] private static readonly double[] bernoulli = Bernoulli(); + + public static NumericsComplex Polygamma(NumericsComplex order, NumericsComplex z) + { + if (order.Imaginary != 0 || double.IsNaN(order.Real) || double.IsInfinity(order.Real) + || Math.Floor(order.Real) != order.Real || order.Real < 0 || order.Real > 170) + return NaN; + var n = (int)order.Real; + if (double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) + return NaN; + if (double.IsInfinity(z.Real) || double.IsInfinity(z.Imaginary)) + // Toward +oo on the real line the digamma function grows as ln x and the others vanish. + return z.Imaginary == 0 && double.IsPositiveInfinity(z.Real) ? (n == 0 ? double.PositiveInfinity : 0) : NaN; + if (z.Imaginary == 0 && z.Real <= 0 && Math.Floor(z.Real) == z.Real) + return NaN; + return Of(n, z); + } + + public static double Polygamma(double order, double z) + { + var value = Polygamma(new NumericsComplex(order, 0), new NumericsComplex(z, 0)); + return value.Imaginary == 0 ? value.Real : double.NaN; + } + + private static NumericsComplex Of(int n, NumericsComplex z) + { + // Reflected only near the axis, where moving right passes its poles, and far left of it, where it + // would take thousands of steps: off the axis cot(pi z) is near i or -i, where 1 + cot^2 vanishes, + // and Q_n's terms cancel to all but a few digits. + if (z.Real < 0.5 && (Math.Abs(z.Imaginary) < 1 || z.Real < -1000)) + { + // psi^(n)(z) = (-1)^n psi^(n)(1 - z) - pi^(n + 1) Q_n(cot(pi z)) + var reflected = Of(n, 1 - z); + if (n % 2 == 1) + reflected = -reflected; + var cot = Cot(z); + return reflected - Math.Pow(Math.PI, n + 1) * Evaluate(CotangentDerivative(n), cot); + } + var radius = 16 + n; + var shift = Math.Max(0, (int)Math.Ceiling(radius - z.Real)); + // sum_(k < shift) 1/(z + k)^(n + 1) + NumericsComplex sum = 0; + for (var k = 0; k < shift; k++) + sum += 1 / NumericsComplex.Pow(z + k, n + 1); + var tail = Asymptotic(n, z + shift); + // psi^(n)(z) = psi^(n)(z + shift) - (-1)^n n! sum + var scaled = SpecialFunctions.Factorial((double)n) * sum; + return n % 2 == 0 ? tail - scaled : tail + scaled; + } + + // ln w - 1/(2 w) - sum B_(2k)/(2k w^(2k)) for the digamma function, and + // (-1)^(n + 1) ((n - 1)!/w^n + n!/(2 w^(n + 1)) + sum B_(2k) (2k + n - 1)!/((2k)! w^(2k + n))) otherwise. + private static NumericsComplex Asymptotic(int n, NumericsComplex w) + { + var inverse = 1 / w; + var inverseSquared = inverse * inverse; + NumericsComplex sum, power; + if (n == 0) + { + sum = NumericsComplex.Log(w) - inverse / 2; + power = inverseSquared; + } + else + { + var toTheN = NumericsComplex.Pow(inverse, n); + sum = SpecialFunctions.Factorial((double)(n - 1)) * toTheN + SpecialFunctions.Factorial((double)n) / 2 * toTheN * inverse; + power = toTheN * inverseSquared; + } + for (var k = 1; k <= Terms; k++) + { + // B_(2k)/(2k) for the digamma function, B_(2k) (2k + n - 1)!/(2k)! otherwise. + var coefficient = bernoulli[k - 1]; + if (n == 0) + coefficient /= 2 * k; + else + for (var i = 0; i < n - 1; i++) + coefficient *= 2 * k + 1 + i; + var term = coefficient * power; + sum = n == 0 ? sum - term : sum + term; + if (NumericsComplex.Abs(term) <= 1e-17 * NumericsComplex.Abs(sum)) + break; + power *= inverseSquared; + } + return n == 0 || n % 2 == 1 ? sum : -sum; + } + + /// + /// The harmonic number H_n^(m): the sum 1 + 1/2^m + ... + 1/n^m for a whole + /// n >= 0 up to 10000, no value at a negative whole n but for a whole m <= 0, + /// and elsewhere, for a whole m from 1 to 64, + /// (-1)^m/(m - 1)! (psi^(m - 1)(1) - psi^(m - 1)(n + 1)). + /// + public static NumericsComplex HarmonicNumber(NumericsComplex n, NumericsComplex m) + { + var wholeN = IsWhole(n, out var number); + var wholeM = IsWhole(m, out var order); + if (wholeN && number >= 0 && number <= 10000) + { + NumericsComplex sum = 0; + for (var k = 1; k <= number; k++) + sum += NumericsComplex.Pow(k, -m); + return sum; + } + if (wholeN && number < 0 && !(wholeM && order <= 0)) + return NaN; + if (wholeM && order >= 1 && order <= 64) + { + var sign = order % 2 == 0 ? 1 : -1; + var o = (int)order; + return sign / SpecialFunctions.Factorial((double)(o - 1)) + * (Polygamma(o - 1, 1) - Polygamma(o - 1, n + 1)); + } + return NaN; + } + + public static double HarmonicNumber(double n, double m) + { + var value = HarmonicNumber(new NumericsComplex(n, 0), new NumericsComplex(m, 0)); + return value.Imaginary == 0 ? value.Real : double.NaN; + } + + private static bool IsWhole(NumericsComplex z, out long whole) + { + whole = 0; + if (z.Imaginary != 0 || double.IsNaN(z.Real) || double.IsInfinity(z.Real) + || Math.Floor(z.Real) != z.Real || Math.Abs(z.Real) > 9e15) + return false; + whole = (long)z.Real; + return true; + } + + /// + /// cot(pi z), accurate where it vanishes: Q_n's coefficients and pi^(n + 1) magnify its + /// error, and at a half-whole z cos(pi z) is a rounding away from 0. The real part is reduced + /// by its nearest integer, the period, and near a half cot(pi r) is tan(pi (1/2 - r)), whose + /// argument is then exact. Complex.Tan is taken only of a small argument: of one near pi/2 it + /// takes the imaginary part as sinh(2y)/(cos(2x) + cosh(2y)), 0/0 on the axis. + /// + private static NumericsComplex Cot(NumericsComplex z) + { + var x = z.Real - Math.Round(z.Real); + if (Math.Abs(x) <= 0.25) + { + var r = Math.PI * new NumericsComplex(x, z.Imaginary); + return NumericsComplex.Cos(r) / NumericsComplex.Sin(r); + } + var half = x > 0 ? 0.5 : -0.5; + return NumericsComplex.Tan(Math.PI * new NumericsComplex(half - x, -z.Imaginary)); + } + + private static NumericsComplex Evaluate(double[] coefficients, NumericsComplex x) + { + NumericsComplex value = 0; + for (var i = coefficients.Length - 1; i >= 0; i--) + value = value * x + coefficients[i]; + return value; + } + + /// The coefficients, lowest first, of Q_n, with d^n/dz^n cot(pi z) = pi^n Q_n(cot(pi z)). + private static double[] CotangentDerivative(int n) + { + var q = new double[] { 0, 1 }; + for (var step = 0; step < n; step++) + { + // -(1 + c^2) q'(c) + var derivative = new double[q.Length - 1]; + for (var i = 1; i < q.Length; i++) + derivative[i - 1] = q[i] * i; + var next = new double[derivative.Length + 2]; + for (var i = 0; i < derivative.Length; i++) + { + next[i] -= derivative[i]; + next[i + 2] -= derivative[i]; + } + q = next; + } + return q; + } + + // B_(2k) = (-1)^(k - 1) 2k T_k/(4^k (4^k - 1)), the tangent numbers by Brent and Harvey's recurrence. + private static double[] Bernoulli() + { + var t = new EInteger[Terms + 1]; + t[1] = EInteger.One; + for (var i = 2; i <= Terms; i++) + t[i] = t[i - 1].Multiply(i - 1); + for (var i = 2; i <= Terms; i++) + for (var j = i; j <= Terms; j++) + t[j] = t[j - 1].Multiply(j - i).Add(t[j].Multiply(j - i + 2)); + var values = new double[Terms]; + for (var k = 1; k <= Terms; k++) + { + var fourToTheK = EInteger.One.ShiftLeft(2 * k); + var value = ERational.Create(t[k].Multiply(2 * k), fourToTheK.Multiply(fourToTheK.Subtract(EInteger.One))).ToDouble(); + values[k - 1] = k % 2 == 1 ? value : -value; + } + return values; + } + } +} diff --git a/Sources/Tests/AotSmokeTest/Program.cs b/Sources/Tests/AotSmokeTest/Program.cs index 23563ef3f..f1dc291f8 100644 --- a/Sources/Tests/AotSmokeTest/Program.cs +++ b/Sources/Tests/AotSmokeTest/Program.cs @@ -311,6 +311,8 @@ private static void Compilation() $"{fastBinomial.Call(5)}, {linqBinomial(5)}"); var compiledW = "lambertw(-1, x)".Compile(x); CheckClose("compiled lambertw(-1, x) at x = -0.3", compiledW(-0.3), -1.7813370234216276); + var compiledPolygamma = "polygamma(1, x) + harmonicnumber(x, 2)".Compile(x); + CheckClose("compiled polygamma(1, x) + harmonicnumber(x, 2) at x = 1/2", compiledPolygamma(0.5), 4.934802200544679 + 0.7101318663035472); // FastExpression: a statement is its truth value, and a piecewise the case that holds. var fastPiecewise = "piecewise(x provided x > 0, -x)".Compile(x); var fastStatement = "x > 1 or x < 0".Compile(x); diff --git a/Sources/Tests/UnitTests/Common/PolygammaCompilationTest.cs b/Sources/Tests/UnitTests/Common/PolygammaCompilationTest.cs new file mode 100644 index 000000000..6a2d3cfcc --- /dev/null +++ b/Sources/Tests/UnitTests/Common/PolygammaCompilationTest.cs @@ -0,0 +1,69 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System.Numerics; +using AngouriMath; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Common +{ + /// + /// The polygamma functions and the harmonic numbers compile, in FastExpression and in the LINQ + /// compiler, to what evaluation gives: at its poles, either side of them, left of 1/2 where + /// the reflection is taken, far out, and off the axis. Each was refused. + /// #2195 + /// + [Trait("Area", "Common")] + public sealed class PolygammaCompilationTest + { + private static readonly string[] orders = { "0", "1", "2", "3", "7", "1/2" }; + private static readonly string[] arguments = + { "1", "2", "1/2", "-1/2", "-5/2", "0", "-3", "3/10", "-7/10", "25", "1000", "i", "1 + i", "-2 - 3i", "1/2 + 10i" }; + + private static Complex Numeric(string value) => ((Entity.Number.Complex)value.ToEntity().Evaled).ToNumerics(); + + // Both NaN, or within a part in a billion. + private static bool Agree(Complex expected, Complex actual) + => double.IsNaN(expected.Real) || double.IsNaN(expected.Imaginary) + ? double.IsNaN(actual.Real) || double.IsNaN(actual.Imaginary) + : Complex.Abs(expected - actual) <= 1e-9 * System.Math.Max(1, Complex.Abs(expected)); + + [Theory] + [InlineData("polygamma(n, z)")] + [InlineData("harmonicnumber(z, n)")] + public void ItIsWhatEvaluationSays(string expression) + { + var fast = expression.ToEntity().Compile("n", "z"); + var linq = expression.ToEntity().Compile("n", "z"); + var real = expression.ToEntity().Compile("n", "z"); + var compared = 0; + foreach (var order in orders) + foreach (var argument in arguments) + { + if (expression.ToEntity().Substitute("n", order).Substitute("z", argument).Evaled is not Entity.Number.Complex evaluated) + continue; + var expected = evaluated.ToNumerics(); + var (n, z) = (Numeric(order), Numeric(argument)); + Assert.True(Agree(expected, fast.Call(n, z)), $"{expression} at ({order}, {argument}): {expected} evaluated, {fast.Call(n, z)} by FastExpression"); + Assert.True(Agree(expected, linq(n, z)), $"{expression} at ({order}, {argument}): {expected} evaluated, {linq(n, z)} by LINQ"); + if (z.Imaginary == 0) + { + var overTheReals = System.Math.Abs(expected.Imaginary) <= 1e-12 * System.Math.Max(1, Complex.Abs(expected)) + ? new Complex(expected.Real, 0) : new Complex(double.NaN, 0); + Assert.True(Agree(overTheReals, real(n.Real, z.Real)), $"{expression} at ({order}, {argument}): {overTheReals} over the reals, {real(n.Real, z.Real)} over double"); + } + compared++; + } + Assert.True(compared > orders.Length * arguments.Length * 3 / 4, $"only {compared} of {orders.Length * arguments.Length} evaluated"); + } + + [Fact] + public void TheDigammaFunctionCompiles() + => Assert.Equal(-0.5772156649015329, "digamma(x)".ToEntity().Compile("x")(1), 14); + } +} diff --git a/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs b/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs index 068f32a7e..687c1f163 100644 --- a/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs +++ b/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs @@ -57,9 +57,10 @@ [Fact] public void UncompilableNode4() => // of these were nodes 2.0 added and never taught the compiler; floor, ceil and round have // been taught since (https://github.com/asc-community/AngouriMath/issues/1603), phi, // gamma and the factorial (https://github.com/asc-community/AngouriMath/issues/1607), and - // the binomial coefficient (https://github.com/asc-community/AngouriMath/issues/2195). + // the binomial coefficient and the polygamma functions + // (https://github.com/asc-community/AngouriMath/issues/2195). [Theory] - [InlineData("polygamma(1, x)")] + [InlineData("gamma(2, x)")] [InlineData("gcd(x, 6)")] [InlineData("lcm(x, 6)")] [InlineData("prime(x)")]