From 90768f52ec86b1a49196a1287cb85995f08f3119 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 3 Nov 2026 04:05:45 +0000 Subject: [PATCH] Take a compiled binomial of a large whole order through a ratio of gammas MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A compiled binomial(n, k) with a whole k is a product of k factors while k fits a long, and past that it went through the gamma function. Past the longs the logarithms of Γ(k + 1) and Γ(n - k + 1), both near k ln k, cancelled to their rounding error, and all three compilers answered 1: binomial(-5/2, 10^20) is 7.52e29, and binomial(277.7, 36!) is 2e-11027, which is 0 in a double. Within the longs, an n whose factors never take the product to 0 ran k steps: binomial(-1/2, 10^12) did not return. From an order of a million, with |n| at most half of it and n not a whole number at least 0, the binomial is now (-1)^k Γ(k - n)/(Γ(k + 1) Γ(-n)), with ln Γ(k - n) - ln Γ(k + 1) from Stirling's series written as differences, so that nothing near k ln k is subtracted, and ln(1 + z) by its series where z is small. Every double past 2^53 is even. A whole n at least 0 below such a k, and a negative k, are 0, as they are within the longs. 2.6.0 did not compile the binomial. The values in the test are mpmath's. Found by compilesweep at seed 13; no issue was filed. Co-Authored-By: Claude Opus 5.5 --- .../Numerics/CombinatorialFunctions.cs | 59 +++++++++++++++++- ...ACompiledBinomialOfALargeWholeOrderTest.cs | 61 +++++++++++++++++++ 2 files changed, 119 insertions(+), 1 deletion(-) create mode 100644 Sources/Tests/UnitTests/Common/ACompiledBinomialOfALargeWholeOrderTest.cs 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); + } +}