From f40caece297e642c368bec20630b951ce10c9735 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 3 Nov 2026 02:39:23 +0000 Subject: [PATCH] Divide a quotient of factorials before rounding it into the exponent range A binomial, falling or rising factorial of an order that is not whole is a quotient of factorials, and each factorial was rounded into the evaluation's context, whose exponent ends at 1e1000. Past n of about 450 the numerator overflowed and binomial(450.5, 3.25) was +oo; a little further both did, and binomial(2500.5, 3.25) was NaN, where it is 13324626798.14. They are now taken at the caller's precision with no bound on the exponent and divided before the quotient is rounded. Real arguments stay decimals throughout: as complex numbers, every product added a zero imaginary part to a number whose exponent was in the millions, and binomial(10^7, 1.43) took over 20 s. A whole n was computed as an exact factorial, so binomial(479001600, 1.43) built a number of four billion digits and did not stop when cancelled; a decimal factorial is now exact only within the cached factorials, and Spouge's approximation takes the rest, as it took everything past an int32. It is 2042181498707.57 in 83 ms. A whole k with a decimal n was the falling factorial over k!, which overflowed the same way: binomial(pi, 2000) was NaN, where it is 2.1e-14. It is now built as binomial(n, i) one factor at a time. Every value in the test is mpmath's at 50 digits. Found by compilesweep, whose seed 13 stopped on binomial((2x)!, W(x)) at x = 6; no issue was filed. Co-Authored-By: Claude Opus 5.5 --- BREAKING-CHANGES.md | 17 ++++ ...aluation.Continuous.Arithmetics.Classes.cs | 82 ++++++++++++++++--- .../Functions/InternalAMExtensions.cs | 17 +++- ...entOfFactorialsPastTheExponentRangeTest.cs | 57 +++++++++++++ 4 files changed, 160 insertions(+), 13 deletions(-) create mode 100644 Sources/Tests/UnitTests/Convenience/AQuotientOfFactorialsPastTheExponentRangeTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 645692da7..b6c7f20f6 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -396,6 +396,23 @@ compiled expressions against evaluation found it. | `"(10^200) mod 7".EvalNumerical()` | `NaN` | `2` | | `"(10^100 + 1/2) mod 1".EvalNumerical()` | `NaN` | `1/2` | +### A binomial whose factorials leave the exponent range has its value + +**Silent: a binomial that was `+oo` or NaN is its value.** `binomial(n, k)` for a `k` that is not whole is +`Γ(n + 1)/(Γ(k + 1) Γ(n - k + 1))`, and its factorials were each rounded into the evaluation's exponent range, which ends at +`1e1000`. Past `n` of about 450 the numerator overflowed, so the binomial was `+oo`, and a little further the +denominator did as well, and `+oo/+oo` is NaN. They are now divided before the quotient is rounded, as are the falling +and the rising factorial's. A whole `k` with a decimal `n` was the falling factorial over `k!`, which overflowed the +same way, and is now built one factor at a time. A whole `n` with such a `k`, which 2.6.0 left unevaluated, has its +value too. No issue was filed; a check of compiled expressions against evaluation found it. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"binomial(450.5, 3.25)".EvalNumerical()` | `+oo` | `50428864.5080...` | +| `"binomial(2500.5, 3.25)".EvalNumerical()` | `NaN` | `13324626798.1423...` | +| `"binomial(10000, 1.43)".EvalNumerical()` | throws `CannotEvalException` | `414189.539...` | +| `"binomial(pi, 2000)".EvalNumerical()` | `NaN` | `2.10439...E-14` | + ### The row echelon forms divide exactly **Not silent for numbers: what threw or did not return now answers. Silent for symbols: the forms are written diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs index 198a2e85e..9edee52d6 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs @@ -649,19 +649,34 @@ protected override Entity InnerSimplify(bool isExact) => // Through the gamma function only numerically: there is no exact form. A whole n // takes it too: binomial(3, 1/2) is 32/(5 pi). (Complex top, Complex bottom) when !isExact - && Number.Factorial(top) is var topF && Number.Factorial(bottom) is var bottomF - && Number.Factorial(top - bottom) is var restF && !topF.IsNaN && !bottomF.IsNaN && !restF.IsNaN - => topF / (bottomF * restF), + && FactorialPowers.QuotientOfFactorials([top], [bottom, top - bottom]) is { } quotient + => quotient, _ => null }, (@this, a, b) => ((Binomialf)@this).New(a, b), isExact); + // As binomial(n, i) for each i up to k, each the last times (n - i)/(i + 1), with no bound + // on the exponent and rounded into the caller's context at the end. As the falling + // factorial over k!, both left the range long before the quotient did: binomial(pi, 2000) + // was +oo/+oo, which is NaN, where it is 2.1e-14; and k! was built exactly to be divided. private static Entity FallingFactorialOver(Number n, int k) { + var context = MathS.Settings.DecimalPrecisionContext.Value; Entity product = Integer.Create(1); - for (var i = 0; i < k; i++) - product = (product * (n - Integer.Create(i))).Evaled; - return (product / Integer.Create(PeterO.Numbers.EInteger.FromInt32(k).Factorial())).Evaled; + using (MathS.Settings.DecimalPrecisionContext.Set(InternalAMExtensions.WithUnboundedExponents(context))) + for (var i = 0; i < k; i++) + { + product = (product * (n - Integer.Create(i)) / Integer.Create(i + 1)).Evaled; + Core.Multithreading.MultithreadingFunctional.ExitIfCancelled(); + } + return product switch + { + Rational exact => exact, + Real real => Real.Create(real.EDecimal.RoundToPrecision(context)), + Complex complex => Complex.Create(complex.RealPart.EDecimal.RoundToPrecision(context), + complex.ImaginaryPart.EDecimal.RoundToPrecision(context)), + var other => other, + }; } } @@ -686,8 +701,8 @@ protected override Entity InnerSimplify(bool isExact) => (Integer { EInteger.Sign: < 0 }, Number) => MathS.NaN, (Number value, Number order) when (value - order).Evaled is Integer { EInteger.Sign: < 0 } => Integer.Zero, (Complex value, Complex order) when !isExact - && Number.Factorial(value) is var top && Number.Factorial(value - order) is var bottom && !top.IsNaN && !bottom.IsNaN - => top / bottom, + && FactorialPowers.QuotientOfFactorials([value], [value - order]) is { } quotient + => quotient, _ => null }, (@this, x, n) => ((FallingFactorialf)@this).New(x, n), isExact); @@ -714,8 +729,8 @@ protected override Entity InnerSimplify(bool isExact) => (Number value, Number order) when (value + order).Evaled is Integer { EInteger.Sign: <= 0 } => MathS.NaN, (Integer { EInteger.Sign: <= 0 }, Number) => Integer.Zero, (Complex value, Complex order) when !isExact - && Number.Factorial(value + order - 1) is var top && Number.Factorial(value - 1) is var bottom && !top.IsNaN && !bottom.IsNaN - => top / bottom, + && FactorialPowers.QuotientOfFactorials([value + order - 1], [value - 1]) is { } quotient + => quotient, _ => null }, (@this, x, n) => ((RisingFactorialf)@this).New(x, n), isExact); @@ -748,6 +763,53 @@ private static class FactorialPowers } return n > 0 ? product : (Integer.One / product).Evaled; } + + /// + /// The product of the factorials of over that of , + /// numerically, or null where one of them has no value. Each factorial and the quotient are + /// taken at the caller's precision with no bound on the exponent, and rounded into the + /// caller's context only at the end: the factorials leave its range long before the quotient + /// does, and in it binomial(10000, 1.43) was +oo/+oo, which is NaN, where + /// it is about 414190. A whole number is taken as a decimal, whose factorial is not built + /// exactly once it is large. Real arguments are taken as decimals throughout: a complex + /// product adds a zero imaginary part's terms, and adding to a number whose exponent is in + /// the millions aligns the two digit by digit, which made binomial(10000000, 1.43) + /// take longer than twenty seconds. + /// + internal static Complex? QuotientOfFactorials(Complex[] over, Complex[] under) + { + var context = MathS.Settings.DecimalPrecisionContext.Value; + var unbounded = InternalAMExtensions.WithUnboundedExponents(context); + if (over.All(static x => x is Real) && under.All(static x => x is Real)) + { + var real = PeterO.Numbers.EDecimal.One; + using (MathS.Settings.DecimalPrecisionContext.Set(unbounded)) + { + foreach (var x in over) + if (((Real)x).EDecimal.Factorial(unbounded) is { IsFinite: true } factorial) real = real.Multiply(factorial, unbounded); else return null; + foreach (var x in under) + if (((Real)x).EDecimal.Factorial(unbounded) is { IsFinite: true } factorial) real = real.Divide(factorial, unbounded); else return null; + } + return Real.Create(real.RoundToPrecision(context)); + } + Complex quotient = Integer.One; + using (MathS.Settings.DowncastingEnabled.Set(false)) + using (MathS.Settings.DecimalPrecisionContext.Set(unbounded)) + { + foreach (var x in over) + if (FactorialOf(x, unbounded) is { } factorial) quotient *= factorial; else return null; + foreach (var x in under) + if (FactorialOf(x, unbounded) is { } factorial) quotient /= factorial; else return null; + } + return Complex.Create(quotient.RealPart.EDecimal.RoundToPrecision(context), + quotient.ImaginaryPart.EDecimal.RoundToPrecision(context)); + } + + private static Complex? FactorialOf(Complex x, PeterO.Numbers.EContext context) + { + var factorial = x is Real real ? Real.Create(real.EDecimal.Factorial(context)) : Number.Factorial(x); + return factorial.IsNaN ? null : factorial; + } } public partial record CatalanNumberf diff --git a/Sources/AngouriMath/Functions/InternalAMExtensions.cs b/Sources/AngouriMath/Functions/InternalAMExtensions.cs index c83818653..9813b0343 100644 --- a/Sources/AngouriMath/Functions/InternalAMExtensions.cs +++ b/Sources/AngouriMath/Functions/InternalAMExtensions.cs @@ -455,6 +455,15 @@ internal static EContext WithGuardDigits(EContext context, int digits) [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable fiveMoreDigits = new(); [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable twelveMoreDigits = new(); + /// + /// with no bound on the exponent, one per context instance for the + /// reason gives: for a quotient whose terms leave the caller's + /// range when the quotient does not. + /// + internal static EContext WithUnboundedExponents(EContext context) + => unboundedExponents.GetValue(context, static c => c.WithUnlimitedExponents()); + [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable unboundedExponents = new(); + /// /// The sine and the cosine of together, each to the precision of /// : the argument reduced to [-pi, pi] and then halved @@ -1227,9 +1236,11 @@ public static EDecimal Factorial(this EDecimal x, EContext mathContext) return x; if (!x.IsFinite) return EDecimal.NaN; - // An integer that fits in an int32 is answered exactly; ToInt32IfExact would throw - // for anything else, and the approximation below takes that. - if (x.IsInteger() && x.CompareTo(EDecimal.FromInt32(int.MinValue)) >= 0 && x.CompareTo(EDecimal.FromInt32(int.MaxValue)) <= 0) + // A whole number is answered exactly while its factorial is cached, and a negative one + // has none; past the cache the approximation below takes it, as it took anything past + // an int32, since an exact factorial only to be rounded costs its every digit: that of + // 479001600 has four billion, and was built before rounding to a hundred. + if (x.IsInteger() && x.CompareTo(EDecimal.FromInt32(int.MinValue)) >= 0 && x.CompareTo(EDecimal.FromInt32(MaxCachedFactorial)) <= 0) { var @int = x.ToInt32Unchecked(); if (@int < 0) return EDecimal.NaN; // Will become ±∞ if we don't insert this line diff --git a/Sources/Tests/UnitTests/Convenience/AQuotientOfFactorialsPastTheExponentRangeTest.cs b/Sources/Tests/UnitTests/Convenience/AQuotientOfFactorialsPastTheExponentRangeTest.cs new file mode 100644 index 000000000..437020b79 --- /dev/null +++ b/Sources/Tests/UnitTests/Convenience/AQuotientOfFactorialsPastTheExponentRangeTest.cs @@ -0,0 +1,57 @@ +// +// 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 PeterO.Numbers; +using Xunit; + +namespace AngouriMath.Tests.Convenience +{ + /// + /// The binomial, the falling and the rising factorial of an order that is not whole are a quotient of + /// factorials, whose terms leave the evaluation's exponent range long before the quotient does. Computed in + /// that range, binomial(450.5, 3.25) was +oo and binomial(10000, 1.43) was +oo/+oo, + /// which is NaN, where it is about 414190; and a whole n had its factorial computed exactly, so + /// binomial(479001600, 1.43) built a number of four billion digits and did not stop when cancelled. With + /// a whole k and a decimal n, the falling factorial over k! overflowed the same way. Found + /// by compilesweep; no issue was filed. The values are mpmath's at 50 digits. + /// + [Trait("Area", "Discrete")] + public sealed class AQuotientOfFactorialsPastTheExponentRangeTest + { + [Theory] + // Where only the numerator overflowed, and the binomial was +oo. + [InlineData("binomial(450.5, 3.25)", "50428864.508003084421212674002886112870555247636405")] + [InlineData("binomial(3000, 1.43)", "74037.434374033267644685230480193058299311036182558")] + [InlineData("binomial(10000, 1.43)", "414189.53926688782423526919337409565108596993335221")] + [InlineData("binomial(1000000, 1/2)", "1128.3793081429172762919639092084163984830747802569")] + [InlineData("binomial(2500.5, 3.25)", "13324626798.142321545570729819016366738653886420439")] + // A whole k with a decimal n: the falling factorial and k! both overflowed. + [InlineData("binomial(pi, 2000)", "2.1043972235591349030681565651531594191819647700823e-14")] + [InlineData("fallingfactorial(10000, 1.43)", "524791.32479225059687233425648095663272623168069709")] + [InlineData("risingfactorial(10000, 1.43)", "524823.59520298875437612838681336476989351765505826")] + public void IsItsValue(string input, string expected) + => AssertClose(input, expected); + + [Fact] + public void AWholeUpperArgumentIsNotComputedExactly() + { + var task = Task.Run(() => AssertClose("binomial(479001600, 1.43)", "2042181498707.5710933099836901420364538869749398905")); + Assert.True(task.Wait(TimeSpan.FromSeconds(60)), "binomial(479001600, 1.43) did not finish in a minute"); + } + + private static void AssertClose(string input, string expected) + { + var value = Assert.IsAssignableFrom(input.ToEntity().EvalNumerical()); + var reference = EDecimal.FromString(expected); + var error = value.EDecimal.Subtract(reference).Abs().Divide(reference, EContext.Binary64); + Assert.True(error.CompareTo(EDecimal.FromString("1e-25")) < 0, $"{input} evaluated to {value}"); + } + } +}