diff --git a/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs b/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs
index b72972e3b..0c3129450 100644
--- a/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs
+++ b/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 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
///
+ // 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;
+
+ ///
+ /// binomial(n, k) for a whole of at least a million and an
+ /// at most half of it that is not a whole number at least 0: (-1)^k Γ(k - n)/(Γ(k + 1) Γ(-n)), with
+ /// the logarithm of Γ(k - n)/Γ(k + 1) by ; null for any other. As a product
+ /// it took k steps, which for binomial(-1/2, 10^12), whose factors never take it to 0, is 10^12; and past
+ /// the longs it went through the gamma function, whose two logarithms near k ln k cancelled to their
+ /// rounding error: binomial(-5/2, 10^20) was 1, where it is 7.52e29.
+ ///
+ 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;
+ }
+
+ ///
+ /// ln Γ(k + a) - ln Γ(k + b) for a of at least a million and an
+ /// and at most half of it, by Stirling's series written as differences, so that nothing
+ /// near k ln k is subtracted: (a - b) ln k + (k + a - 1/2) ln(1 + a/k) - (k + b - 1/2) ln(1 + b/k)
+ /// - (a - b) and the series' next two terms, whose remainder at k/2 is below 10^-30.
+ ///
+ 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;
+ }
+
+ /// ln(1 + z), by its series where z is small, which ln of 1 + z would round away.
+ 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);
diff --git a/Sources/Tests/UnitTests/Common/ACompiledBinomialOfALargeWholeOrderTest.cs b/Sources/Tests/UnitTests/Common/ACompiledBinomialOfALargeWholeOrderTest.cs
new file mode 100644
index 000000000..91edb72e9
--- /dev/null
+++ b/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
+{
+ ///
+ /// A compiled binomial(n, k) for a whole k far past n. Past the longs it went through
+ /// the gamma function, whose two logarithms near k ln k cancelled to their rounding error, and was 1:
+ /// binomial(-5/2, 10^20) is 7.52e29 and binomial(277.7, 36!) is 2e-11027. Within
+ /// them it was a product of k factors, which for binomial(-1/2, 10^12) never reach 0. Found by
+ /// compilesweep; no issue was filed. The values are mpmath's at 30 digits, to 20.
+ ///
+ [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("x", "y");
+ var fast = "binomial(x, y)".ToEntity().Compile("x", "y");
+ var linq = "binomial(x, y)".ToEntity().Compile("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);
+ }
+}