A compiled gamma ratio past the doubles is taken in logarithms

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<double, double>("x")(54)` | `UncompilableNodeException` | `5317936260` |
 | `"binomial(x, 1/2)".ToEntity().Compile<double, double>("x")(3)` | `UncompilableNodeException` | `2.037183271576257` |
+| `"binomial(x, pi)".ToEntity().Compile<double, double>("x")(216)` | `UncompilableNodeException` | `2954596.6790617765` |
 
 ### The number-theoretic functions compile
 

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
 

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);
 
+        /// <summary>
+        /// 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: <c>binomial(216, pi)</c> is about <c>2954596.68</c>
+        /// and <c>216!</c> about <c>1e408</c>, so the ratio was <c>∞/∞</c>, <c>NaN</c>. 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
+        /// </summary>
+        private static bool InTheDoubles(NumericsComplex ratio)
+            => ratio != 0 && !IsInfinite(ratio) && !double.IsNaN(ratio.Real) && !double.IsNaN(ratio.Imaginary);
+
+        /// <summary>
+        /// The ratio whose logarithm is <paramref name="logarithm"/>, where its terms are past the doubles. A real one's
+        /// logarithm has a multiple of <c>pi</c> for an imaginary part, odd where the ratio is negative, which is its sign.
+        /// </summary>
+        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;
     }

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);
 
+        /// <summary>
+        /// A logarithm of <c>z!</c>, <c>Γ(z + 1)</c>, 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
+        /// <c>pi</c>, odd where the factorial is negative. https://github.com/asc-community/AngouriMath/issues/2248
+        /// </summary>
+        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);
+        }
+
         /// <summary>A logarithm of the gamma function for <c>Re z &gt;= 1/2</c>, by Lanczos's approximation.</summary>
         private static NumericsComplex LogGamma(NumericsComplex z)
         {

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
+{
+    /// <summary>
+    /// A compiled binomial coefficient, factorial power or Catalan number whose factorials are past the doubles is its
+    /// value, where it was <c>NaN</c>: <c>binomial(216, pi)</c> is about <c>2954596.68</c>, and <c>216!</c>, about
+    /// <c>1e408</c>, made the gamma ratio <c>∞/∞</c>. The values are mpmath 1.3.0's at 50 digits.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/2248">#2248</a>
+    /// </summary>
+    [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<Complex, Complex>("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<double, double>("x")(0);
+                Assert.True(Agree(expected, overDoubles), $"{expression}: {overDoubles} over double, where {expected}");
+            }
+        }
+
+        /// <summary>Past the doubles, at about <c>9.8e486</c> and <c>8.2e597</c>, the value is <c>+oo</c>.</summary>
+        [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<Complex, Complex>("x")(0));
+            Assert.Equal(double.PositiveInfinity, compiled.Compile<double, double>("x")(0));
+        }
+    }
+}