Skip to content

Compiled combinatorial functions at whole arguments are exact #2218

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Oct 30, 2026
Merged

Compiled combinatorial functions at whole arguments are exact #2218

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion BREAKING-CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double, double>("x")(54)` | `UncompilableNodeException` | `5317936260` |
| `"binomial(x, 1/2)".ToEntity().Compile<double, double>("x")(3)` | `UncompilableNodeException` | `2.037183271576257` |

### The number-theoretic functions compile
Expand Down
114 changes: 99 additions & 15 deletions Sources/AngouriMath/Numerics/CombinatorialFunctions.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
//

using System;
using PeterO.Numbers;
using NumericsComplex = System.Numerics.Complex;

namespace AngouriMath.Numerics
Expand All @@ -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 <c>NaN</c> at a pole of the
/// numerator and <c>0</c> at one of the denominator, where the reciprocal gamma vanishes.
/// numerator and <c>0</c> 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
/// </summary>
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;

/// <summary>Whether <paramref name="z"/> is a whole number a <see cref="long"/> holds exactly.</summary>
private static bool IsWhole(NumericsComplex z, out long whole)
{
Expand All @@ -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);
Expand Down Expand Up @@ -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)));
Expand All @@ -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)
Expand All @@ -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++)
{
Expand All @@ -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.
Expand Down
88 changes: 88 additions & 0 deletions Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();

/// <summary>
/// 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: <c>binomial(54, 45)</c> would be
/// <c>5317936259.999995</c>.
/// </summary>
[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<Complex, Complex, Complex>("x", "y");
var real = expression.ToEntity().Compile<double, double, double>("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");
}

/// <summary>Past 2^53 and past the range of a double, a whole binomial is still evaluation's, rounded once.</summary>
[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<double, double, double>("x", "y")(n, k));

/// <summary>
/// The Catalan numbers: exact on the whole numbers, <c>-1/2</c> at <c>-1</c>, <c>0</c> below it, no value at a
/// negative half-integer, and the gamma ratio elsewhere.
Expand All @@ -80,13 +137,44 @@ public void TheCatalanNumbersAreWhatEvaluationSays()
}
}

/// <summary>
/// A Catalan number at a whole index is evaluation's, rounded once, past 2^53 and past the doubles: rounded by
/// truncation, as <c>(double)BigInteger</c> rounds, <c>C_32</c> would be <c>55534064877048190</c>.
/// </summary>
[Fact]
public void AtAWholeIndexTheCatalanNumberIsEvaluationsRoundedOnce()
{
var fast = "catalannumber(x)".ToEntity().Compile("x");
var real = "catalannumber(x)".ToEntity().Compile<double, double>("x");
for (var n = 0; n <= 1000; n++)
{
var expected = Exactly(Assert.IsAssignableFrom<Entity.Number.Complex>($"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");
}
}

/// <summary>A Catalan number at a whole index is the integer, exactly, as far as a double holds it.</summary>
[Theory]
[InlineData(5, 42)]
[InlineData(30, 3814986502092304)]
public void ACatalanNumberAtAWholeIndexIsExact(int n, long expected)
=> Assert.Equal((double)expected, "catalannumber(x)".ToEntity().Compile("x").Call(n).Real);

/// <summary>The double factorial at a whole number is evaluation's, rounded once, past 2^53 and past the doubles.</summary>
[Fact]
public void AtAWholeNumberTheDoubleFactorialIsExact()
{
var fast = "x!!".ToEntity().Compile("x");
var real = "x!!".ToEntity().Compile<double, double>("x");
for (var n = -301; n <= 320; n++)
{
var expected = Exactly(Assert.IsAssignableFrom<Entity.Number.Complex>($"({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()
{
Expand Down
Loading