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}"); + } + } +}