From 552d23ae7fa20a83cd81b2b6b95825749beb5b28 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 30 Oct 2026 02:39:19 +0000 Subject: [PATCH] Compiled combinatorial functions at whole arguments are exact The compiled binomial took a whole lower argument as a product of ratios in doubles, (n - i)/(i + 1), each rounded, so binomial(54, 45) came out 5317936259.999995 where evaluation gives 5317936260. Of the 2485 binomials with 0 <= k <= n < 70, 1173 were not the integer. The factorial powers and the double factorial, products of whole factors, rounded at each step once a partial product passed 2^53, and so were off in the last places from there. At whole arguments each is now computed exactly and rounded once: - the binomial by binomial(m, i + 1) = binomial(m, i) (m - i)/(i + 1), whose division leaves no remainder, over the shorter of k and m - k, and for a negative n by (-1)^k binomial(k - n - 1, k); - the products in doubles while each partial product is below 2^53, where every step is exact, and past it in EInteger, rounded by EFloat or, for a negative order's reciprocal, by ERational. The Catalan numbers were computed exactly, in #2216, and rounded by (double)BigInteger, which truncates: 246 of the first 1000 were a unit low in the last place, C_32 = 55534064877048200 coming out 55534064877048190. They are rounded by EFloat now. Nothing else in the library rounds a BigInteger to a double for a result. A factor of 0 is checked for before the product, so a partial product past the doubles' range is known to stay past it, and the loop stops there. Checked against evaluation exactly, as doubles, in both compilers: every two-argument function on a grid of 90 whole first arguments, up to 10^12, and 58 whole second arguments from -12 to 45; the binomial past 2^53 and past the doubles, among others at (60, 30), (1100, 550), (3000, 1500) and (-1, 4001); the double factorial from -301 to 320; and the Catalan numbers from 0 to 1000. The expected value is the exact rational rounded once by ERational, since ToNumerics goes through a decimal whose exponent has a floor, and reads 1/299!! as 0. On master 9 of these tests fail. Part of #2195. Co-Authored-By: Claude Opus 5.5 (1M context) --- BREAKING-CHANGES.md | 4 +- .../Numerics/CombinatorialFunctions.cs | 114 +++++++++++++++--- .../Common/CombinatorialCompilationTest.cs | 88 ++++++++++++++ 3 files changed, 190 insertions(+), 16 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 35f5e6f2e..9aa7e4a08 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -328,13 +328,15 @@ negative even integer or off the integers, since the extension through the gamma **Not silent: what was refused is compiled.** `binomial` was refused by both compilers with `UncompilableNodeException`, and `fallingfactorial`, `risingfactorial` and `n!!`, new since 2.6.0, compile from the start. Each compiles in double precision and on the interpreter's terms: a whole order is a product, and any other goes -through the gamma function, with `NaN` at a pole of the numerator and `0` at one of the denominator. The double +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)). | 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` | ### The number-theoretic functions compile diff --git a/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs b/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs index 849ac540a..02604bc6c 100644 --- a/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs +++ b/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs @@ -6,6 +6,7 @@ // using System; +using PeterO.Numbers; using NumericsComplex = System.Numerics.Complex; namespace AngouriMath.Numerics @@ -14,13 +15,20 @@ namespace AngouriMath.Numerics /// The binomial coefficient, the falling and the rising factorial and the double factorial in /// double precision, for the compilers, on the terms the interpreter takes them: a whole order is /// a product, and any other goes through the gamma function, with NaN at a pole of the - /// numerator and 0 at one of the denominator, where the reciprocal gamma vanishes. + /// numerator and 0 at one of the denominator, where the reciprocal gamma vanishes. At whole + /// arguments each is computed exactly and rounded once, as evaluation's exact value is. /// https://github.com/asc-community/AngouriMath/issues/2195 /// internal static class CombinatorialFunctions { [AngouriMath.Core.ConstantField] private static readonly NumericsComplex NaN = new(double.NaN, 0); + // 2^53: a double holds every whole number below it, so a product of two of them below it is exact. + private const double ExactBelow = 9007199254740992; + + // Past this many bits a whole number is past the range of a double, whose largest is under 2^1024. + private const long BeyondDoubles = 1100; + /// Whether is a whole number a holds exactly. private static bool IsWhole(NumericsComplex z, out long whole) { @@ -44,6 +52,8 @@ public static NumericsComplex Binomial(NumericsComplex n, NumericsComplex k) { if (lower < 0 || IsWhole(n, out var upper) && upper >= 0 && lower > upper) return 0; + if (IsWhole(n, out var whole)) + return WholeBinomial(whole, lower); NumericsComplex product = 1; for (long i = 0; i < lower && product != 0 && !IsInfinite(product); i++) product *= (n - i) / (i + 1); @@ -105,16 +115,10 @@ public static NumericsComplex DoubleFactorial(NumericsComplex n) { if (!IsWhole(n, out var whole) || whole < 0 && whole % 2 == 0) return NaN; - NumericsComplex product = 1; - if (whole >= 0) - { - for (var k = whole; k > 1 && !IsInfinite(product); k -= 2) - product *= k; - return product; - } - for (var j = whole + 2; j <= -1; j += 2) - product *= j; - return 1 / product; + // n (n - 2) ... down to 2 or 1, and for a negative odd n the reciprocal of (n + 2) (n + 4) ... (-1). + return whole >= 0 + ? ProductOfWholeFactors(whole, -2, (whole + 1) / 2, reciprocal: false) + : ProductOfWholeFactors(whole + 2, 2, (-whole - 1) / 2, reciprocal: true); } public static double DoubleFactorial(double n) => Real(DoubleFactorial(new NumericsComplex(n, 0))); @@ -134,13 +138,19 @@ public static NumericsComplex CatalanNumber(NumericsComplex n) return 0; // Exactly, by C_(k + 1) = C_k 2 (2k + 1)/(k + 2), whose division leaves no remainder, and rounded // once, since a product of the ratios in doubles rounds at each step and C_5 would be - // 41.99999999999999. Past the index the interpreter computes exactly, the gamma ratio. + // 41.99999999999999. By EFloat, as (double)BigInteger truncates, and C_32 would be + // 55534064877048190. Past the index the interpreter computes exactly, the gamma ratio. if (whole <= 1000) { - var catalan = System.Numerics.BigInteger.One; + var catalan = EInteger.One; for (var k = 0; k < whole; k++) - catalan = catalan * (2 * (2 * k + 1)) / (k + 2); - return (double)catalan; + { + catalan = catalan.Multiply(2 * (2 * k + 1)).Divide(k + 2); + // The Catalan numbers grow, so once one is past the doubles the rest are. + if (catalan.GetUnsignedBitLengthAsInt64() > BeyondDoubles) + return double.PositiveInfinity; + } + return EFloat.FromEInteger(catalan).ToDouble(); } } if (IsWhole(n + 0.5, out var half) && half <= 0) @@ -157,6 +167,8 @@ public static NumericsComplex CatalanNumber(NumericsComplex n) private static NumericsComplex OfAWholeOrder(NumericsComplex x, long n, bool falling) { var step = falling == n > 0 ? -1 : 1; + if (IsWhole(x, out var whole)) + return ProductOfWholeFactors(n > 0 ? whole : whole + step, step, Math.Abs(n), reciprocal: n < 0); NumericsComplex product = 1; for (long i = 0; i < Math.Abs(n) && product != 0 && !IsInfinite(product); i++) { @@ -168,6 +180,78 @@ private static NumericsComplex OfAWholeOrder(NumericsComplex x, long n, bool fal return n > 0 ? product : 1 / product; } + // binomial(n, k) for a whole n and a whole k >= 0, k <= n where n >= 0, exactly and rounded once: by + // binomial(m, i + 1) = binomial(m, i) (m - i)/(i + 1), whose division leaves no remainder, over the + // shorter of k and m - k. In doubles while the product before the division is below 2^53, where it + // and the quotient are exact, and past it in EInteger. For a negative n, binomial(n, k) is + // (-1)^k binomial(k - n - 1, k). + private static double WholeBinomial(long n, long k) + { + var sign = n < 0 && k % 2 != 0 ? -1 : 1; + var top = n < 0 ? k - n - 1 : n; + var steps = Math.Min(k, top - k); + double value = 1; + long i = 0; + for (; i < steps; i++) + { + double factor = top - i; + var product = value * factor; + if (factor >= ExactBelow || product >= ExactBelow) + break; + value = product / (i + 1); + } + if (i == steps) + return sign * value; + var exact = EInteger.FromInt64((long)value); + for (; i < steps; i++) + { + exact = exact * EInteger.FromInt64(top - i) / EInteger.FromInt64(i + 1); + // The binomials grow up to the middle, so once one is past the doubles the rest are. + if (exact.GetUnsignedBitLengthAsInt64() > BeyondDoubles) + return sign * double.PositiveInfinity; + } + return sign * EFloat.FromEInteger(exact).ToDouble(); + } + + // The product of the whole factors first, first + step, ..., count of them, or its reciprocal, + // exactly and rounded once: in doubles while each partial product is below 2^53, where it is exact, + // and past it in EInteger. 0 where a factor vanishes, and NaN for its reciprocal. The callers' factors + // are consecutive, or two apart without crossing 0, so otherwise they have one sign and are each at + // least 1 in size: a partial product past the range of the doubles stays past it. + private static double ProductOfWholeFactors(long first, long step, long count, bool reciprocal) + { + if (count == 0) + return 1; + var last = first + step * (count - 1); + if (Math.Min(first, last) <= 0 && Math.Max(first, last) >= 0 && first % step == 0) + return reciprocal ? double.NaN : 0; + double product = 1; + long i = 0; + for (; i < count; i++) + { + double factor = first + step * i; + var next = product * factor; + if (Math.Abs(factor) >= ExactBelow || Math.Abs(next) >= ExactBelow) + break; + product = next; + } + if (i == count) + return reciprocal ? 1 / product : product; + var exact = EInteger.FromInt64((long)product); + for (; i < count; i++) + { + exact *= EInteger.FromInt64(first + step * i); + if (exact.GetUnsignedBitLengthAsInt64() > BeyondDoubles) + { + var sign = first < 0 && count % 2 != 0 ? -1 : 1; + return reciprocal ? sign * 0.0 : sign * double.PositiveInfinity; + } + } + return reciprocal + ? ERational.Create(exact.Sign < 0 ? EInteger.FromInt32(-1) : EInteger.One, exact.Abs()).ToDouble() + : EFloat.FromEInteger(exact).ToDouble(); + } + private static bool IsInfinite(NumericsComplex z) => double.IsInfinity(z.Real) || double.IsInfinity(z.Imaginary); // The real value of a real argument's result, and NaN where it is not real. diff --git a/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs b/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs index 146c02ca7..38d2f4825 100644 --- a/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs +++ b/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs @@ -59,6 +59,63 @@ public void ATwoArgumentFunctionIsWhatEvaluationSays(string expression) Assert.True(compared > values.Length * values.Length * 3 / 4, $"only {compared} of {values.Length * values.Length} evaluated"); } + private static readonly string[] wholeFirsts = + Enumerable.Range(-25, 86).Select(n => n.ToString()).Concat(new[] { "100", "1000", "54321", "1000000000000" }).ToArray(); + + private static readonly string[] wholeSeconds = Enumerable.Range(-12, 58).Select(n => n.ToString()).ToArray(); + + // Equal as doubles, or both NaN. + private static bool Same(Complex expected, Complex actual) + => double.IsNaN(expected.Real) ? double.IsNaN(actual.Real) : expected == actual; + + // A rational value rounded once to a double, where ToNumerics goes through a decimal whose exponent has a floor, + // so that 1/299!! would be 0. + private static Complex Exactly(Entity.Number.Complex value) + => value is Entity.Number.Rational rational ? rational.ERational.ToDouble() : value.ToNumerics(); + + /// + /// At whole arguments each is computed exactly and rounded once, so it is evaluation's exact value as a + /// double, where a product of ratios in doubles rounds at each step: binomial(54, 45) would be + /// 5317936259.999995. + /// + [Theory] + [InlineData("binomial(x, y)")] + [InlineData("fallingfactorial(x, y)")] + [InlineData("risingfactorial(x, y)")] + public void AtWholeArgumentsATwoArgumentFunctionIsExact(string expression) + { + var fast = expression.ToEntity().Compile("x", "y"); + var linq = expression.ToEntity().Compile("x", "y"); + var real = expression.ToEntity().Compile("x", "y"); + var compared = 0; + foreach (var a in wholeFirsts) + foreach (var b in wholeSeconds) + { + if (expression.ToEntity().Substitute("x", a).Substitute("y", b).Evaled is not Entity.Number.Complex evaluated) + continue; + var expected = Exactly(evaluated); + var (x, y) = (double.Parse(a), double.Parse(b)); + Assert.True(Same(expected, fast.Call(x, y)), $"{expression} at ({a}, {b}): {expected} evaluated, {fast.Call(x, y)} by FastExpression"); + Assert.True(Same(expected, linq(x, y)), $"{expression} at ({a}, {b}): {expected} evaluated, {linq(x, y)} by LINQ"); + Assert.True(Same(expected, real(x, y)), $"{expression} at ({a}, {b}): {expected} evaluated, {real(x, y)} over double"); + compared++; + } + Assert.True(compared > wholeFirsts.Length * wholeSeconds.Length * 3 / 4, $"only {compared} of {wholeFirsts.Length * wholeSeconds.Length} evaluated"); + } + + /// Past 2^53 and past the range of a double, a whole binomial is still evaluation's, rounded once. + [Theory] + [InlineData(54, 45)] + [InlineData(60, 30)] + [InlineData(-40, 25)] + [InlineData(1000000000000, 5)] + [InlineData(1100, 550)] + [InlineData(3000, 1500)] + [InlineData(-1, 4001)] + public void AWholeBinomialIsEvaluationsRoundedOnce(long n, long k) + => Assert.Equal(Exactly((Entity.Number.Complex)$"binomial({n}, {k})".ToEntity().Evaled).Real, + "binomial(x, y)".ToEntity().Compile("x", "y")(n, k)); + /// /// The Catalan numbers: exact on the whole numbers, -1/2 at -1, 0 below it, no value at a /// negative half-integer, and the gamma ratio elsewhere. @@ -80,6 +137,23 @@ public void TheCatalanNumbersAreWhatEvaluationSays() } } + /// + /// A Catalan number at a whole index is evaluation's, rounded once, past 2^53 and past the doubles: rounded by + /// truncation, as (double)BigInteger rounds, C_32 would be 55534064877048190. + /// + [Fact] + public void AtAWholeIndexTheCatalanNumberIsEvaluationsRoundedOnce() + { + var fast = "catalannumber(x)".ToEntity().Compile("x"); + var real = "catalannumber(x)".ToEntity().Compile("x"); + for (var n = 0; n <= 1000; n++) + { + var expected = Exactly(Assert.IsAssignableFrom($"catalannumber({n})".ToEntity().Evaled)); + Assert.True(Same(expected, fast.Call(n)), $"catalannumber({n}): {expected} evaluated, {fast.Call(n)} by FastExpression"); + Assert.True(Same(expected, real(n)), $"catalannumber({n}): {expected} evaluated, {real(n)} over double"); + } + } + /// A Catalan number at a whole index is the integer, exactly, as far as a double holds it. [Theory] [InlineData(5, 42)] @@ -87,6 +161,20 @@ public void TheCatalanNumbersAreWhatEvaluationSays() public void ACatalanNumberAtAWholeIndexIsExact(int n, long expected) => Assert.Equal((double)expected, "catalannumber(x)".ToEntity().Compile("x").Call(n).Real); + /// The double factorial at a whole number is evaluation's, rounded once, past 2^53 and past the doubles. + [Fact] + public void AtAWholeNumberTheDoubleFactorialIsExact() + { + var fast = "x!!".ToEntity().Compile("x"); + var real = "x!!".ToEntity().Compile("x"); + for (var n = -301; n <= 320; n++) + { + var expected = Exactly(Assert.IsAssignableFrom($"({n})!!".ToEntity().Evaled)); + Assert.True(Same(expected, fast.Call(n)), $"({n})!!: {expected} evaluated, {fast.Call(n)} by FastExpression"); + Assert.True(Same(expected, real(n)), $"({n})!!: {expected} evaluated, {real(n)} over double"); + } + } + [Fact] public void TheDoubleFactorialIsWhatEvaluationSays() {