From 5bb9b71417addce8366a399fdbd031c1c40396a9 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 30 Oct 2026 14:39:41 +0000 Subject: [PATCH] A compiled gamma ratio past the doubles is taken in logarithms A compiled binomial coefficient, falling or rising factorial or Catalan number at an argument that is not whole is a ratio of gamma functions. Past 170! a factorial is past the doubles where the ratio need not be, so the ratio was oo/oo, NaN: binomial(216, pi) is about 2954596.68, fallingfactorial(300, 1/2) about 17.33, catalannumber(2017/5) about 5.16e238, and each was NaN. A ratio whose terms are under the doubles was 0/0, as binomial(-601/2, 7/3) was, which is about 109419.24. Where the ratio does not come out as a number, it is now the exponential of the difference of the factorials' logarithms. SpecialFunctions.LogFactorial is the logarithm of Gamma(z + 1), by the Lanczos approximation the compiled gamma function already uses, and by reflection below 1/2. At a real argument its imaginary part is a multiple of pi, odd where the factorial is negative, which gives the ratio's sign. None of the ratios is 0 away from a pole, which the callers answer first, so a ratio that came out as a number is kept as it was, to the last place. GammaRatioCompilationTest compares ten values with mpmath 1.3.0's at 50 digits, under both compilers and over double, and all ten fail on master. Two of them are past the doubles, binomial(2000, 1001/2) and catalannumber(1001), and are +oo. Evaluation is NaN for three of the ten, binomial(10, 601/2), binomial(2000, 1001/2) and catalannumber(-2017/5), and EvalNumerical throws CannotEvalException for catalannumber(1001). The last two of the three are past 1e200 or under 1e-200, the decimal range #1498 is about. Those are evaluation's, and not changed here. Found by work/compilesweep. Fixes #2248. Co-Authored-By: Claude Opus 5.5 (1M context) --- BREAKING-CHANGES.md | 5 +- Sources/.editorconfig | 3 + .../Numerics/CombinatorialFunctions.cs | 35 +++++++++-- .../AngouriMath/Numerics/SpecialFunctions.cs | 11 ++++ .../Common/GammaRatioCompilationTest.cs | 63 +++++++++++++++++++ 5 files changed, 112 insertions(+), 5 deletions(-) create mode 100644 Sources/Tests/UnitTests/Common/GammaRatioCompilationTest.cs 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)); + } + } +}