Divide a quotient of factorials before rounding it into the exponent range

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

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;
             }
+
+            /// <summary>
+            /// The product of the factorials of <paramref name="over"/> over that of <paramref name="under"/>,
+            /// 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 <c>binomial(10000, 1.43)</c> was <c>+oo/+oo</c>, which is <c>NaN</c>, 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 <c>binomial(10000000, 1.43)</c>
+            /// take longer than twenty seconds.
+            /// </summary>
+            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

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<EContext, EContext> fiveMoreDigits = new();
         [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable<EContext, EContext> twelveMoreDigits = new();
 
+        /// <summary>
+        /// <paramref name="context"/> with no bound on the exponent, one per context instance for the
+        /// reason <see cref="WithGuardDigits"/> gives: for a quotient whose terms leave the caller's
+        /// range when the quotient does not.
+        /// </summary>
+        internal static EContext WithUnboundedExponents(EContext context)
+            => unboundedExponents.GetValue(context, static c => c.WithUnlimitedExponents());
+        [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable<EContext, EContext> unboundedExponents = new();
+
         /// <summary>
         /// The sine and the cosine of <paramref name="x"/> together, each to the precision of
         /// <paramref name="context"/>: the argument reduced to <c>[-pi, pi]</c> 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

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
+{
+    /// <summary>
+    /// 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, <c>binomial(450.5, 3.25)</c> was <c>+oo</c> and <c>binomial(10000, 1.43)</c> was <c>+oo/+oo</c>,
+    /// which is <c>NaN</c>, where it is about 414190; and a whole <c>n</c> had its factorial computed exactly, so
+    /// <c>binomial(479001600, 1.43)</c> built a number of four billion digits and did not stop when cancelled. With
+    /// a whole <c>k</c> and a decimal <c>n</c>, the falling factorial over <c>k!</c> overflowed the same way. Found
+    /// by compilesweep; no issue was filed. The values are mpmath's at 50 digits.
+    /// </summary>
+    [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<Entity.Number.Real>(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}");
+        }
+    }
+}