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