Take a compiled binomial of a large whole order through a ratio of gammas

Sources/AngouriMath/Numerics/CombinatorialFunctions.cs

@@ -34,7 +34,7 @@ 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)
+                || Math.Floor(z.Real) != z.Real || Math.Abs(z.Real) > LongestWhole)
                 return false;
             whole = (long)z.Real;
             return true;
@@ -54,11 +54,21 @@ public static NumericsComplex Binomial(NumericsComplex n, NumericsComplex k)
                     return 0;
                 if (IsWhole(n, out var whole))
                     return WholeBinomial(whole, lower);
+                if (OfALargeWholeOrder(n, k.Real) is { } large)
+                    return large;
                 NumericsComplex product = 1;
                 for (long i = 0; i < lower && product != 0 && !IsInfinite(product); i++)
                     product *= (n - i) / (i + 1);
                 return product;
             }
+            // Every double past the longs is whole: the cases above, for an order there.
+            if (k.Imaginary == 0 && Math.Abs(k.Real) > LongestWhole && !double.IsInfinity(k.Real))
+            {
+                if (k.Real < 0 || n.Imaginary == 0 && Math.Floor(n.Real) == n.Real && n.Real >= 0 && k.Real > n.Real)
+                    return 0;
+                if (OfALargeWholeOrder(n, k.Real) is { } large)
+                    return large;
+            }
             if (IsWhole(n, out var negative) && negative < 0)
                 return NaN;
             if (IsWhole(n - k, out var difference) && difference < 0)
@@ -270,6 +280,53 @@ private static double ProductOfWholeFactors(long first, long step, long count, b
         /// 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>
+        // The largest whole number IsWhole takes.
+        private const double LongestWhole = 9e15;
+
+        // From this order on, a binomial whose n is far below it is taken through the logarithm of a ratio of gammas.
+        private const double LargeOrder = 1e6;
+
+        /// <summary>
+        /// <c>binomial(n, k)</c> for a whole <paramref name="k"/> of at least a million and an <paramref name="n"/>
+        /// at most half of it that is not a whole number at least 0: <c>(-1)^k Γ(k - n)/(Γ(k + 1) Γ(-n))</c>, with
+        /// the logarithm of <c>Γ(k - n)/Γ(k + 1)</c> by <see cref="LogGammaRatio"/>; null for any other. As a product
+        /// it took k steps, which for <c>binomial(-1/2, 10^12)</c>, whose factors never take it to 0, is 10^12; and past
+        /// the longs it went through the gamma function, whose two logarithms near <c>k ln k</c> cancelled to their
+        /// rounding error: <c>binomial(-5/2, 10^20)</c> was 1, where it is <c>7.52e29</c>.
+        /// </summary>
+        private static NumericsComplex? OfALargeWholeOrder(NumericsComplex n, double k)
+        {
+            if (k < LargeOrder || n.Magnitude > k / 2 || n.Imaginary == 0 && Math.Floor(n.Real) == n.Real && n.Real >= 0)
+                return null;
+            var value = FromLogarithm(LogGammaRatio(-n, 1, k) - SpecialFunctions.LogFactorial(-n - 1), n.Imaginary == 0);
+            // Past 2^53 every double is even.
+            return k < ExactBelow && ((long)k & 1) == 1 ? -value : value;
+        }
+
+        /// <summary>
+        /// <c>ln Γ(k + a) - ln Γ(k + b)</c> for a <paramref name="k"/> of at least a million and an <paramref name="a"/>
+        /// and <paramref name="b"/> at most half of it, by Stirling's series written as differences, so that nothing
+        /// near <c>k ln k</c> is subtracted: <c>(a - b) ln k + (k + a - 1/2) ln(1 + a/k) - (k + b - 1/2) ln(1 + b/k)
+        /// - (a - b)</c> and the series' next two terms, whose remainder at <c>k/2</c> is below <c>10^-30</c>.
+        /// </summary>
+        private static NumericsComplex LogGammaRatio(NumericsComplex a, NumericsComplex b, double k)
+        {
+            NumericsComplex za = k + a, zb = k + b;
+            return (a - b) * Math.Log(k) + (za - 0.5) * LogOnePlus(a / k) - (zb - 0.5) * LogOnePlus(b / k) - (a - b)
+                + (1 / za - 1 / zb) / 12 - (1 / (za * za * za) - 1 / (zb * zb * zb)) / 360;
+        }
+
+        /// <summary><c>ln(1 + z)</c>, by its series where <c>z</c> is small, which <c>ln</c> of <c>1 + z</c> would round away.</summary>
+        private static NumericsComplex LogOnePlus(NumericsComplex z)
+        {
+            if (z.Magnitude > 1e-3)
+                return NumericsComplex.Log(1 + z);
+            NumericsComplex sum = 0, power = z;
+            for (var i = 1; i <= 7; i++, power *= -z)
+                sum += power / i;
+            return sum;
+        }
+
         private static bool InTheDoubles(NumericsComplex ratio)
             => ratio != 0 && !IsInfinite(ratio) && !double.IsNaN(ratio.Real) && !double.IsNaN(ratio.Imaginary);
 

Sources/Tests/UnitTests/Common/ACompiledBinomialOfALargeWholeOrderTest.cs

@@ -0,0 +1,61 @@
+//
+// 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 System.Threading.Tasks;
+using AngouriMath.Extensions;
+using Xunit;
+using Complex = System.Numerics.Complex;
+
+namespace AngouriMath.Tests.Common
+{
+    /// <summary>
+    /// A compiled <c>binomial(n, k)</c> for a whole <c>k</c> far past <c>n</c>. Past the longs it went through
+    /// the gamma function, whose two logarithms near <c>k ln k</c> cancelled to their rounding error, and was 1:
+    /// <c>binomial(-5/2, 10^20)</c> is <c>7.52e29</c> and <c>binomial(277.7, 36!)</c> is <c>2e-11027</c>. Within
+    /// them it was a product of k factors, which for <c>binomial(-1/2, 10^12)</c> never reach 0. Found by
+    /// compilesweep; no issue was filed. The values are mpmath's at 30 digits, to 20.
+    /// </summary>
+    [Trait("Area", "Compilation")]
+    public sealed class ACompiledBinomialOfALargeWholeOrderTest
+    {
+        [Theory]
+        [InlineData(-2.5, 1e20, 7.5225277806367504928e+29)]
+        [InlineData(277.7, 1e20, 0)]
+        [InlineData(277.7, 3.7199332678990125e41, 0)]
+        [InlineData(-0.5, 1e12, 5.6418958354768576325e-7)]
+        [InlineData(-0.5, 1e12 + 1, -5.6418958354740366846e-7)]
+        [InlineData(3.25, 2e6, 3.0992465282944267377e-27)]
+        [InlineData(-7.5, 1e7 + 1, -1.6899292703910971418e+42)]
+        // A whole n at least 0 below k has nothing to choose, and nor has a negative k.
+        [InlineData(5, 1e20, 0)]
+        [InlineData(-2.5, -1e20, 0)]
+        public void IsItsValue(double n, double k, double expected)
+        {
+            var real = "binomial(x, y)".ToEntity().Compile<double, double, double>("x", "y");
+            var fast = "binomial(x, y)".ToEntity().Compile("x", "y");
+            var linq = "binomial(x, y)".ToEntity().Compile<Complex, Complex, Complex>("x", "y");
+            var task = Task.Run(() => (real(n, k), fast.Call(n, k), linq(n, k)));
+            Assert.True(task.Wait(TimeSpan.FromSeconds(30)), $"binomial({n}, {k}) did not finish");
+            var (overDouble, byFast, byLinq) = task.Result;
+            Assert.True(Close(expected, overDouble), $"binomial({n}, {k}) over double is {overDouble}");
+            Assert.True(Close(expected, byFast.Real) && byFast.Imaginary == 0, $"binomial({n}, {k}) by FastExpression is {byFast}");
+            Assert.True(Close(expected, byLinq.Real) && byLinq.Imaginary == 0, $"binomial({n}, {k}) by LINQ is {byLinq}");
+        }
+
+        [Fact]
+        public void AComplexUpperArgumentIsItsValue()
+        {
+            var value = "binomial(x, y)".ToEntity().Compile("x", "y").Call(new Complex(0.5, 1), 1e7);
+            var expected = new Complex(6.5666206043323881952e-11, -1.7327461680588210534e-11);
+            Assert.True((value - expected).Magnitude < 1e-12 * expected.Magnitude, $"binomial(1/2 + i, 10^7) is {value}");
+        }
+
+        private static bool Close(double expected, double actual)
+            => expected == 0 ? actual == 0 : Math.Abs(actual - expected) < 1e-12 * Math.Abs(expected);
+    }
+}