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 9edee52d6..d051593ba 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 @@ -672,9 +672,10 @@ private static Entity FallingFactorialOver(Number n, int k) 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)), + Real real => Real.CreateAtScale(real.EDecimal.RoundToPrecision(context), PeterO.Numbers.EDecimal.Zero), + Complex complex => Complex.CreateAtScales(complex.RealPart.EDecimal.RoundToPrecision(context), + complex.ImaginaryPart.EDecimal.RoundToPrecision(context), + PeterO.Numbers.EDecimal.Zero, PeterO.Numbers.EDecimal.Zero), var other => other, }; } @@ -774,7 +775,9 @@ private static class FactorialPowers /// 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. + /// take longer than twenty seconds. The quotient is made as one that cancels nothing, and is + /// never rounded onto zero: made as any number, a value under the downcasting tolerance was + /// taken for a residual, and fallingfactorial(0.5, -100.5), 9.4e-161, was 0. /// internal static Complex? QuotientOfFactorials(Complex[] over, Complex[] under) { @@ -790,7 +793,7 @@ private static class FactorialPowers 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)); + return Real.CreateAtScale(real.RoundToPrecision(context), PeterO.Numbers.EDecimal.Zero); } Complex quotient = Integer.One; using (MathS.Settings.DowncastingEnabled.Set(false)) @@ -801,8 +804,9 @@ private static class FactorialPowers 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)); + return Complex.CreateAtScales(quotient.RealPart.EDecimal.RoundToPrecision(context), + quotient.ImaginaryPart.EDecimal.RoundToPrecision(context), + PeterO.Numbers.EDecimal.Zero, PeterO.Numbers.EDecimal.Zero); } private static Complex? FactorialOf(Complex x, PeterO.Numbers.EContext context) @@ -831,16 +835,35 @@ protected override Entity InnerSimplify(bool isExact) => Integer { EInteger: var n } when n.Equals(PeterO.Numbers.EInteger.FromInt32(-1)) => Rational.Create(-1, 2), Integer => Integer.Zero, Number value when (value + Rational.Create(1, 2)).Evaled is Integer { EInteger.Sign: <= 0 } => MathS.NaN, + // The two gammas as one quotient, and 4^n times it with no bound on the exponent: + // each gamma left the range from n near 450, where the number is near 1e267, and + // catalannumber(500.5) was +oo/+oo, which is NaN. Complex value when !isExact - && Number.Factorial((Complex)(value - Rational.Create(1, 2)).Evaled) is var top && Number.Factorial((Complex)(value + 1).Evaled) is var bottom - && !top.IsNaN && !bottom.IsNaN - => (MathS.Pow(4, value) * top / (MathS.Sqrt(MathS.pi) * bottom)).Evaled, + && FactorialPowers.QuotientOfFactorials([(Complex)(value - Rational.Create(1, 2)).Evaled], [(Complex)(value + 1).Evaled]) is { } quotient + => FourToThePowerTimes(value, quotient), _ => null }, (@this, a) => ((CatalanNumberf)@this).New(a), isExact); /// How far a Catalan number is computed exactly before the node is left as written. private const int LargestIndex = 1000; + + /// + /// 4^n quotient / √π. For a real , formed in decimals with no bound on the exponent + /// and rounded into the context: no number is made of the power itself, which for n = 10^20 has an + /// exponent near 6e19 and cannot be read as the integer a number checks it for. A complex one is formed in + /// the context, as it was. + /// + private static Entity FourToThePowerTimes(Complex n, Complex quotient) + { + if (n is not Real { EDecimal: var power } || quotient is not Real { EDecimal: var ratio }) + return (MathS.Pow(4, n) * quotient / MathS.Sqrt(MathS.pi)).Evaled; + var context = MathS.Settings.DecimalPrecisionContext.Value; + var unbounded = InternalAMExtensions.WithUnboundedExponents(context); + var value = PeterO.Numbers.EDecimal.FromInt32(4).Pow(power, unbounded).Multiply(ratio, unbounded) + .Divide(PeterO.Numbers.EDecimal.PI(unbounded).Sqrt(unbounded), unbounded); + return Real.CreateAtScale(value.RoundToPrecision(context), PeterO.Numbers.EDecimal.Zero); + } } public partial record HarmonicNumberf diff --git a/Sources/Tests/UnitTests/Convenience/ACatalanNumberPastTheExponentRangeTest.cs b/Sources/Tests/UnitTests/Convenience/ACatalanNumberPastTheExponentRangeTest.cs new file mode 100644 index 000000000..21562e0c8 --- /dev/null +++ b/Sources/Tests/UnitTests/Convenience/ACatalanNumberPastTheExponentRangeTest.cs @@ -0,0 +1,58 @@ +// +// 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 +{ + /// + /// A Catalan number of a number that is not whole is 4^n Γ(n + 1/2) / (√π Γ(n + 2)), and each gamma + /// was computed in the evaluation's exponent range: from n = 449.5, whose number is 2.5e266, both were + /// +oo and the quotient NaN, and catalannumber(-60.3) was NaN as well. And a quotient + /// of factorials under the downcasting tolerance was rounded onto zero as though it were what a cancellation left: + /// fallingfactorial(0.5, -100.5), 9.4e-161, was 0. Found by compilesweep, on + /// catalannumber(Ei(sinh(x))) at 6; no issue was filed. The values are mpmath's at 120 digits. + /// + [Trait("Area", "Discrete")] + public sealed class ACatalanNumberPastTheExponentRangeTest + { + [Theory] + [InlineData("catalannumber(10.5)", "31375.760336503176813184117125078909721937218359811")] + [InlineData("catalannumber(449.5)", "2.4958055762434131441628313609834440544021408905592e+266")] + [InlineData("catalannumber(500.5)", "1.0773809195193076398539141613347340467488925064819e+297")] + [InlineData("catalannumber(1600.25)", "2.462290001548975528409180780699097156132733259886e+958")] + [InlineData("catalannumber(-60.3)", "-8.3875179942165737514658221112507292321206364599213e-40")] + // A quotient of factorials under the downcasting tolerance, which is not a residual. + [InlineData("fallingfactorial(0.5, -100.5)", "9.4019927548161528235486105066852002441716212361456e-161")] + [InlineData("risingfactorial(100.5, -60.5)", "2.1883878150169778143163818514630537992056282322519e-111")] + public void IsItsValue(string input, string expected) + { + var value = Assert.IsAssignableFrom(Within(input)); + var reference = EDecimal.FromString(expected); + var error = value.EDecimal.Subtract(reference).Abs().Divide(reference.Abs(), EContext.Binary64); + Assert.True(error.CompareTo(EDecimal.FromString("1e-25")) < 0, $"{input} evaluated to {value}"); + } + + [Theory] + // 1.7e1199, and about 10^(6e39): past the range, where a real is +oo. + [InlineData("catalannumber(2000.5)")] + [InlineData("catalannumber(1e40 + 0.2)")] + public void IsPositiveInfinityPastTheRange(string input) + => Assert.Equal(Entity.Number.Real.PositiveInfinity, Within(input)); + + private static Entity Within(string input) + { + var task = Task.Factory.StartNew(() => input.ToEntity().EvalNumerical(), TaskCreationOptions.LongRunning); + Assert.True(task.Wait(TimeSpan.FromSeconds(30)), $"{input} did not finish in thirty seconds"); + return task.Result; + } + } +}