diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index fef54c9ea..7b8a3415d 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -360,13 +360,16 @@ start. Each compiles in double precision and on the interpreter's terms: a whole through the gamma function, with `NaN` at a pole of the numerator and `0` at one of the denominator. Where every argument is whole, the product is computed exactly and rounded once, so it is evaluation's exact value, as a double. The double factorial has no value off the whole numbers or at a negative even one -([#2195](https://github.com/asc-community/AngouriMath/issues/2195)). +([#2195](https://github.com/asc-community/AngouriMath/issues/2195)). A gamma ratio whose factorials are past the doubles, +as `216!` is, is taken in logarithms, the Catalan number's too, so `binomial(216, pi)` is about `2954596.68` where a +ratio of infinities had no value ([#2248](https://github.com/asc-community/AngouriMath/issues/2248)). | Input | Was (2.6.0) | Now | |---|---|---| | `"binomial(x, 2)".ToEntity().Compile("x").Call(5)` | `UncompilableNodeException` | `(10, 0)` | | `"binomial(x, 45)".ToEntity().Compile("x")(54)` | `UncompilableNodeException` | `5317936260` | | `"binomial(x, 1/2)".ToEntity().Compile("x")(3)` | `UncompilableNodeException` | `2.037183271576257` | +| `"binomial(x, pi)".ToEntity().Compile("x")(216)` | `UncompilableNodeException` | `2954596.6790617765` | ### The number-theoretic functions compile diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 006adb650..7d14ee5a2 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -360,6 +360,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Common/CombinatorialCompilationTest.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/GammaRatioCompilationTest.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/CombinatorialFunctions.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/Numerics/CombinatorialFunctions.cs b/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs index 02604bc6c..b72972e3b 100644 --- a/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs +++ b/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs @@ -63,7 +63,10 @@ public static NumericsComplex Binomial(NumericsComplex n, NumericsComplex k) return NaN; if (IsWhole(n - k, out var difference) && difference < 0) return 0; - return SpecialFunctions.Factorial(n) / (SpecialFunctions.Factorial(k) * SpecialFunctions.Factorial(n - k)); + var ratio = SpecialFunctions.Factorial(n) / (SpecialFunctions.Factorial(k) * SpecialFunctions.Factorial(n - k)); + return InTheDoubles(ratio) ? ratio : FromLogarithm( + SpecialFunctions.LogFactorial(n) - SpecialFunctions.LogFactorial(k) - SpecialFunctions.LogFactorial(n - k), + n.Imaginary == 0 && k.Imaginary == 0); } public static double Binomial(double n, double k) => Real(Binomial(new NumericsComplex(n, 0), new NumericsComplex(k, 0))); @@ -82,7 +85,9 @@ public static NumericsComplex FallingFactorial(NumericsComplex x, NumericsComple return NaN; if (IsWhole(x - n, out var difference) && difference < 0) return 0; - return SpecialFunctions.Factorial(x) / SpecialFunctions.Factorial(x - n); + var ratio = SpecialFunctions.Factorial(x) / SpecialFunctions.Factorial(x - n); + return InTheDoubles(ratio) ? ratio : FromLogarithm( + SpecialFunctions.LogFactorial(x) - SpecialFunctions.LogFactorial(x - n), x.Imaginary == 0 && n.Imaginary == 0); } public static double FallingFactorial(double x, double n) => Real(FallingFactorial(new NumericsComplex(x, 0), new NumericsComplex(n, 0))); @@ -101,7 +106,9 @@ public static NumericsComplex RisingFactorial(NumericsComplex x, NumericsComplex return NaN; if (IsWhole(x, out var bottom) && bottom <= 0) return 0; - return SpecialFunctions.Factorial(x + n - 1) / SpecialFunctions.Factorial(x - 1); + var ratio = SpecialFunctions.Factorial(x + n - 1) / SpecialFunctions.Factorial(x - 1); + return InTheDoubles(ratio) ? ratio : FromLogarithm( + SpecialFunctions.LogFactorial(x + n - 1) - SpecialFunctions.LogFactorial(x - 1), x.Imaginary == 0 && n.Imaginary == 0); } public static double RisingFactorial(double x, double n) => Real(RisingFactorial(new NumericsComplex(x, 0), new NumericsComplex(n, 0))); @@ -157,7 +164,10 @@ public static NumericsComplex CatalanNumber(NumericsComplex n) return NaN; var top = SpecialFunctions.Factorial(n - 0.5); var bottom = SpecialFunctions.Factorial(n + 1); - return NumericsComplex.Pow(4, n) * top / (Math.Sqrt(Math.PI) * bottom); + var ratio = NumericsComplex.Pow(4, n) * top / (Math.Sqrt(Math.PI) * bottom); + return InTheDoubles(ratio) ? ratio : FromLogarithm( + n * Math.Log(4) + SpecialFunctions.LogFactorial(n - 0.5) - SpecialFunctions.LogFactorial(n + 1) - Math.Log(Math.PI) / 2, + n.Imaginary == 0); } public static double CatalanNumber(double n) => Real(CatalanNumber(new NumericsComplex(n, 0))); @@ -254,6 +264,23 @@ private static double ProductOfWholeFactors(long first, long step, long count, b private static bool IsInfinite(NumericsComplex z) => double.IsInfinity(z.Real) || double.IsInfinity(z.Imaginary); + /// + /// Whether a ratio of factorials, at arguments that are none of its poles, came out as a number. A factorial past + /// 170 is past the doubles where a ratio of them need not be: binomial(216, pi) is about 2954596.68 + /// and 216! about 1e408, so the ratio was ∞/∞, NaN. None of the ratios is 0 away from a + /// pole, so a 0 is one whose terms were under the doubles. https://github.com/asc-community/AngouriMath/issues/2248 + /// + private static bool InTheDoubles(NumericsComplex ratio) + => ratio != 0 && !IsInfinite(ratio) && !double.IsNaN(ratio.Real) && !double.IsNaN(ratio.Imaginary); + + /// + /// The ratio whose logarithm is , where its terms are past the doubles. A real one's + /// logarithm has a multiple of pi for an imaginary part, odd where the ratio is negative, which is its sign. + /// + private static NumericsComplex FromLogarithm(NumericsComplex logarithm, bool real) + => !real ? NumericsComplex.Exp(logarithm) + : new NumericsComplex(Math.Cos(logarithm.Imaginary) < 0 ? -Math.Exp(logarithm.Real) : Math.Exp(logarithm.Real), 0); + // The real value of a real argument's result, and NaN where it is not real. private static double Real(NumericsComplex z) => z.Imaginary == 0 ? z.Real : double.NaN; } diff --git a/Sources/AngouriMath/Numerics/SpecialFunctions.cs b/Sources/AngouriMath/Numerics/SpecialFunctions.cs index 966800902..623ceb230 100644 --- a/Sources/AngouriMath/Numerics/SpecialFunctions.cs +++ b/Sources/AngouriMath/Numerics/SpecialFunctions.cs @@ -456,6 +456,17 @@ private static NumericsComplex Gamma(NumericsComplex z) ? ExpTimes(Math.Log(Math.PI) - LogSine(Math.PI * z) - LogGamma(1 - z), 1) : ExpTimes(LogGamma(z), 1); + /// + /// A logarithm of z!, Γ(z + 1), at a number that is not one of its poles, for a ratio of factorials + /// that are past the doubles where the ratio is not. At a real number its imaginary part is a multiple of + /// pi, odd where the factorial is negative. https://github.com/asc-community/AngouriMath/issues/2248 + /// + internal static NumericsComplex LogFactorial(NumericsComplex z) + { + var w = z + 1; + return w.Real < 0.5 ? Math.Log(Math.PI) - LogSine(Math.PI * w) - LogGamma(1 - w) : LogGamma(w); + } + /// A logarithm of the gamma function for Re z >= 1/2, by Lanczos's approximation. private static NumericsComplex LogGamma(NumericsComplex z) { diff --git a/Sources/Tests/UnitTests/Common/GammaRatioCompilationTest.cs b/Sources/Tests/UnitTests/Common/GammaRatioCompilationTest.cs new file mode 100644 index 000000000..1228937f0 --- /dev/null +++ b/Sources/Tests/UnitTests/Common/GammaRatioCompilationTest.cs @@ -0,0 +1,63 @@ +// +// 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 +{ + /// + /// A compiled binomial coefficient, factorial power or Catalan number whose factorials are past the doubles is its + /// value, where it was NaN: binomial(216, pi) is about 2954596.68, and 216!, about + /// 1e408, made the gamma ratio ∞/∞. The values are mpmath 1.3.0's at 50 digits. + /// #2248 + /// + [Trait("Area", "Common")] + public sealed class GammaRatioCompilationTest + { + private static bool Agree(Complex expected, Complex actual) + => Complex.Abs(expected - actual) <= 1e-10 * Complex.Abs(expected); + + [Theory] + [InlineData("binomial(216, pi)", 2954596.6790612643608, 0)] + [InlineData("binomial(10, 601/2)", 7.7046763087758266086e-22, 0)] + [InlineData("binomial(-601/2, 7/3)", 109419.23916039716406, 0)] + [InlineData("fallingfactorial(300, 1/2)", 17.327726454436284054, 0)] + [InlineData("risingfactorial(200, 1/2)", 14.133299559727925473, 0)] + [InlineData("catalannumber(2017/5)", 5.1595508663052781318e+238, 0)] + [InlineData("catalannumber(-2017/5)", -2.8923921724035808851e-247, 0)] + [InlineData("binomial(200 + i, 1/2)", 15.967717653807835931, 0.039869146616281137615)] + public void ItIsTheValue(string expression, double real, double imaginary) + { + var expected = new Complex(real, imaginary); + var compiled = $"x + {expression}".ToEntity(); + var fast = compiled.Compile("x").Call(0); + var linq = compiled.Compile("x")(0); + Assert.True(Agree(expected, fast), $"{expression}: {fast} by FastExpression, where {expected}"); + Assert.True(Agree(expected, linq), $"{expression}: {linq} by LINQ, where {expected}"); + if (imaginary == 0) + { + var overDoubles = compiled.Compile("x")(0); + Assert.True(Agree(expected, overDoubles), $"{expression}: {overDoubles} over double, where {expected}"); + } + } + + /// Past the doubles, at about 9.8e486 and 8.2e597, the value is +oo. + [Theory] + [InlineData("binomial(2000, 1001/2)")] + [InlineData("catalannumber(1001)")] + public void PastTheDoublesItIsInfinite(string expression) + { + var compiled = $"x + {expression}".ToEntity(); + Assert.Equal(new Complex(double.PositiveInfinity, 0), compiled.Compile("x").Call(0)); + Assert.Equal(new Complex(double.PositiveInfinity, 0), compiled.Compile("x")(0)); + Assert.Equal(double.PositiveInfinity, compiled.Compile("x")(0)); + } + } +}