diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 13ca3f066..2072155fb 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -412,6 +412,22 @@ compiled expressions against evaluation found it. | `"(10^200) mod 7".EvalNumerical()` | `NaN` | `2` | | `"(10^100 + 1/2) mod 1".EvalNumerical()` | `NaN` | `1/2` | +### A complex product past the exponent range has its parts, not NaN + +**Silent: a different answer.** A complex number's parts are decimals in the precision context, whose exponents end at +a thousand. A product or a whole power of finite numbers, a part of which passed that end on the way, had `oo - oo` in +it and was `NaN`, which says the value does not exist; `Simplify` returned that `NaN` as well. Each part past the range +is the infinity of its sign now, as a real past it is `+oo`, and a part within it keeps its value. An exact +simplification leaves the expression as written, since the infinity is the context's and not the value's. No issue was +filed. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"(3+4i)^2000".EvalNumerical()` | `NaN` | `+oo + +ooi` | +| `"(1+i)^(10^6)".EvalNumerical()` | `NaN` | `+oo` | +| `"(1e600+1e600i)*(1e600+1e600i)".EvalNumerical()` | `NaN` | `+ooi` | +| `"(3+4i)^2000".Simplify()` | `NaN` | `(3 + 4i) ^ 2000` | + ### 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 diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs index 1070164b5..1b9651440 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs @@ -151,10 +151,40 @@ static EDecimal ModifiedMultiply(EDecimal a, EDecimal b) => a.IsInfinity() && b.IsZero || b.IsInfinity() && a.IsZero ? EDecimal.Zero : CtxMultiply(a, b); var (ac, bd) = (ModifiedMultiply(a.RealPart.EDecimal, b.RealPart.EDecimal), ModifiedMultiply(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal)); var (ad, bc) = (ModifiedMultiply(a.RealPart.EDecimal, b.ImaginaryPart.EDecimal), ModifiedMultiply(a.ImaginaryPart.EDecimal, b.RealPart.EDecimal)); - return Complex.CreateAtScales(CtxSubtract(ac, bd), CtxAdd(ad, bc), CancelledScaleOfADifference(ac, bd), CancelledScale(ad, bc)); + var (real, imaginary) = (CtxSubtract(ac, bd), CtxAdd(ad, bc)); + if (!(real.IsFinite && imaginary.IsFinite) && a.IsFinite && b.IsFinite) + return ProductPastTheRange(a, b); + return Complex.CreateAtScales(real, imaginary, CancelledScaleOfADifference(ac, bd), CancelledScale(ad, bc)); } ); + /// + /// Whether both parts are finite, read from the parts: is a cache made on + /// an entity's first read of it, and asked of every power built, and of every exact value an evaluation + /// takes, it was half a megabyte on each solve of the kernel benchmark's SolveHard. A rational part is + /// finite and is not written out to a decimal to say so. + /// + internal static bool PartsAreFinite(Complex number) => PartIsFinite(number.RealPart) && PartIsFinite(number.ImaginaryPart); + + private static bool PartIsFinite(Real part) => part is Rational || part.EDecimal.IsFinite; + + /// + /// The product of two finite numbers a part of which leaves the context's exponent range. Within it, + /// (1e600 + 1e600 i)^2 had 1e1200 - 1e1200 for its real part, which is oo - oo there, + /// and the product was NaN; the parts are formed past the range and rounded into it, so that it is + /// 2e1200 i rounded, oo i, as a real product past the range is +oo. + /// + private static Complex ProductPastTheRange(Complex a, Complex b) + { + var context = MathS.Settings.DecimalPrecisionContext.Value; + var unbounded = InternalAMExtensions.WithUnboundedExponents(context); + var (ar, ai, br, bi) = (a.RealPart.EDecimal, a.ImaginaryPart.EDecimal, b.RealPart.EDecimal, b.ImaginaryPart.EDecimal); + var (ac, bd) = (ar.Multiply(br, unbounded), ai.Multiply(bi, unbounded)); + var (ad, bc) = (ar.Multiply(bi, unbounded), ai.Multiply(br, unbounded)); + return Complex.CreateAtScales(ac.Subtract(bd, unbounded).RoundToPrecision(context), ad.Add(bc, unbounded).RoundToPrecision(context), + CancelledScaleOfADifference(ac, bd), CancelledScale(ad, bc)); + } + internal static Complex OpDiv(T a, T b) where T : Number => SuperSwitch(a, b, (a, b) => CtxDivide(a.ERational, b.ERational) is { } n ? Rational.Create(n) : Real.NaN, @@ -615,7 +645,18 @@ static Complex BinaryIntPow(Complex num, EInteger val) var raised = decimalBase.Pow(pow.Abs().ToInt32Checked(), wholeContext); return Real.CreateAtScale(pow.Sign < 0 ? EDecimal.One.Divide(raised, wholeContext) : raised, EDecimal.Zero); } - return BinaryIntPow(@base, pow); + var built = BinaryIntPow(@base, pow); + if (@base is Real || PartsAreFinite(built)) + return built; + // A finite complex base whose power leaves the exponent range on the way: each squaring past + // it was oo - oo in a part, and (3 + 4i)^2000 and (1 + i)^(10^6) were NaN. The squarings are + // made past the range and the parts rounded into it at the end. + var rangeContext = MathS.Settings.DecimalPrecisionContext.Value; + Complex pastTheRange; + using (MathS.Settings.DecimalPrecisionContext.Set(InternalAMExtensions.WithUnboundedExponents(rangeContext))) + pastTheRange = BinaryIntPow(@base, pow); + return Complex.Create(pastTheRange.RealPart.EDecimal.RoundToPrecision(rangeContext), + pastTheRange.ImaginaryPart.EDecimal.RoundToPrecision(rangeContext)); } // A half power of a nonnegative decimal is its square root, not its exponential @@ -715,7 +756,7 @@ private static Complex ByThePolarForm(EDecimal baseReal, EDecimal baseImaginary, /// The most bits an exact power is built to: 3^(10^7), at 15.8 million, takes three seconds, /// and each doubling of the exponent doubles the size and more than doubles the time. /// - [ConstantField] private static readonly EInteger MaxExactPowerBits = EInteger.FromInt32(1 << 24); + private const long MaxExactPowerBits = 1 << 24; /// /// About how many bits an exact to the whole has: for @@ -724,20 +765,21 @@ private static Complex ByThePolarForm(EDecimal baseReal, EDecimal baseImaginary, /// and a power of one is as small as it is. A base that is not exact is 0, being a decimal. /// internal static bool IsTooLargeToBuild(Complex @base, EInteger power) - => ExactPowerBits(@base, power).CompareTo(MaxExactPowerBits) > 0; - - private static EInteger ExactPowerBits(Complex @base, EInteger power) { + // In longs: the bound is asked of every exact power, and as EIntegers it was a quarter of a + // megabyte on each solve of the kernel benchmark's SolveHard. if (@base.RealPart is not Rational { ERational: var re } || @base.ImaginaryPart is not Rational { ERational: var im }) - return EInteger.Zero; - static EInteger Bits(ERational part) - => part.IsZero ? EInteger.Zero - : part.Numerator.Abs().GetUnsignedBitLengthAsEInteger().Subtract(1) - .Add(part.Denominator.GetUnsignedBitLengthAsEInteger().Subtract(1)); - var bits = Bits(re).CompareTo(Bits(im)) >= 0 ? Bits(re) : Bits(im); - if (!re.IsZero && !im.IsZero) - bits = bits.Add(1); - return bits.Multiply(power.Abs()); + return false; + static long Bits(ERational part) + => part.IsZero ? 0 : part.Numerator.GetUnsignedBitLengthAsInt64() - 1 + part.Denominator.GetUnsignedBitLengthAsInt64() - 1; + var bits = Math.Max(Bits(re), Bits(im)) + (!re.IsZero && !im.IsZero ? 1 : 0); + if (bits == 0) + return false; + var powerBits = power.GetUnsignedBitLengthAsInt64(); + if (powerBits > 40) + return true; + var magnitude = power.Sign < 0 ? -power.ToInt64Unchecked() : power.ToInt64Unchecked(); + return magnitude > MaxExactPowerBits / bits; } /// e.g. Log(2, 32) = 5 diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs index 487682612..4f8643b92 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs @@ -56,17 +56,22 @@ private static bool IsAZeroOnlyRoundingMade(Entity node, Entity value) && (node is Mulf or Divf or Powf && CannotBeZero(node) || node is Sumf or Minusf && SignOfANonzero(node) is not null); /// - /// Whether is an infinity that only the decimals made: an exact power too - /// large to build is evaluated as a decimal, which is +oo past the context's range, and an - /// infinity counts as exact, so that 2^(10^10), and every expression holding it, was - /// simplified to +oo. Asked only of a value that is not finite. + /// Whether is an infinity that only the decimals made. A sum, a product or a + /// power of exact finite numbers, the power's base not zero, is exact and finite, and is an infinity only + /// where the context's exponent range rounds it onto one: an exact power too large to build is the + /// decimal, and a complex product's parts are decimals, so 2^(10^10) and (3 + 4i)^2000 are + /// infinite as evaluated. An infinity counts as exact, and every expression holding one was simplified to + /// it. An operand that is a decimal already is not asked about: its power past the range is the + /// decimal's +oo, as it was. Asked only of a value that is not finite. /// private static bool IsAnInfinityOnlyRoundingMade(Entity node, Entity value) - => value is Complex { IsFinite: false } && !value.IsNaN - && node.Nodes.Any(static part => part is Powf(var @base, var exponent) - && @base.Evaled is Complex { IsExact: true, IsFinite: true } exactBase - && exponent.Evaled is Integer { EInteger: var power } - && Number.IsTooLargeToBuild(exactBase, power)); + => value is Complex number && !Number.PartsAreFinite(number) && !value.IsNaN && node.Nodes.Any(IsFiniteRoundedOntoAnInfinity); + + private static bool IsFiniteRoundedOntoAnInfinity(Entity part) + => part is Sumf or Minusf or Mulf or Powf + && part.Evaled is Complex { IsFinite: false } evaluated && !evaluated.IsNaN + && part.DirectChildren.All(static operand => operand.Evaled is Complex { IsFinite: true, IsExact: true }) + && !(part is Powf(var @base, _) && @base.Evaled is Complex { IsZero: true }); /// /// Whether is a number or a constant that is not zero, or a product, diff --git a/Sources/Tests/UnitTests/Core/AComplexOverflowIsNotNaNTest.cs b/Sources/Tests/UnitTests/Core/AComplexOverflowIsNotNaNTest.cs new file mode 100644 index 000000000..786bb599b --- /dev/null +++ b/Sources/Tests/UnitTests/Core/AComplexOverflowIsNotNaNTest.cs @@ -0,0 +1,91 @@ +// +// 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; +using static AngouriMath.Entity; +using static AngouriMath.Entity.Number; + +namespace AngouriMath.Tests.Core +{ + /// + /// A complex number's parts are decimals in the precision context, whose exponents end at a thousand. A + /// product or a power of finite numbers a part of which passed that end had oo - oo in it and was + /// NaN, which says the value does not exist: (3 + 4 * i)^2000, (1 + i)^(10^6) and + /// (1e600 + 1e600 * i)^2 among them. Each part past the range is the infinity of its sign now, as a + /// real past it is +oo, and an exact simplification leaves the expression as written. No issue was + /// filed. + /// + [Trait("Area", "Core")] + public sealed class AComplexOverflowIsNotNaNTest + { + [Theory] + [InlineData("(1e600 + 1e600 * i)^2")] + [InlineData("(1e600 + 1e600 * i) * (1e600 + 1e600 * i)")] + public void AProductPastTheRangeHasItsParts(string input) + { + var value = Within(() => input.ToEntity().Evaled, input); + var complex = Assert.IsAssignableFrom(value); + Assert.True(complex.RealPart.EDecimal.IsZero, $"{input} has the real part {complex.RealPart}"); + Assert.True(complex.ImaginaryPart.EDecimal.IsPositiveInfinity(), $"{input} has the imaginary part {complex.ImaginaryPart}"); + } + + [Fact] + public void APowerOfOnePlusIPastTheRangeIsPositiveInfinity() + => Assert.Equal(Real.PositiveInfinity, Within(() => "(1 + i)^(10^6)".ToEntity().Evaled, "(1 + i)^(10^6)")); + + [Theory] + [InlineData(3, 4, 2000)] + [InlineData(3, 4, 2001)] + [InlineData(1, 2, 3000)] + [InlineData(-5, 7, 2500)] + public void AGaussianPowerPastTheRangeHasTheSignsOfItsParts(int real, int imaginary, int power) + { + var input = $"({real} + {imaginary} * i)^{power}"; + var value = Assert.IsAssignableFrom(Within(() => input.ToEntity().Evaled, input)); + var (exactReal, exactImaginary) = (EInteger.One, EInteger.Zero); + for (var k = 0; k < power; k++) + (exactReal, exactImaginary) = (exactReal * real - exactImaginary * imaginary, exactReal * imaginary + exactImaginary * real); + AssertRoundedInto(exactReal, value.RealPart.EDecimal, input); + AssertRoundedInto(exactImaginary, value.ImaginaryPart.EDecimal, input); + } + + [Theory] + [InlineData("(3 + 4 * i)^1400")] + [InlineData("(3 + 4 * i)^(-2000)")] + public void APowerWithinTheRangeIsFinite(string input) + => Assert.True(Within(() => input.ToEntity().Evaled, input).IsFinite, $"{input} is not finite"); + + [Theory] + [InlineData("(3 + 4 * i)^2000")] + [InlineData("(1 + i)^(10^6)")] + [InlineData("(3 + 4 * i)^2000 + x")] + public void AnExactSimplificationLeavesItAsWritten(string input) + { + Assert.True(Within(() => input.ToEntity().InnerSimplified, input) is not Number { IsFinite: false }, $"{input} was simplified to a number that is not finite"); + Assert.True(Within(() => input.ToEntity().Simplify(), input) is not Number { IsFinite: false }, $"{input} was simplified to a number that is not finite"); + } + + /// The exact part rounded into the context: the same infinity past it, and the same sign within. + private static void AssertRoundedInto(EInteger exact, EDecimal evaluated, string input) + { + var rounded = EDecimal.FromEInteger(exact).RoundToPrecision(MathS.Settings.DecimalPrecisionContext); + Assert.True(rounded.IsInfinity() == evaluated.IsInfinity() && rounded.Sign == evaluated.Sign, + $"{input} has the part {evaluated}, where the exact part is {rounded}"); + } + + private static Entity Within(Func evaluation, string input) + { + var task = Task.Factory.StartNew(evaluation, TaskCreationOptions.LongRunning); + Assert.True(task.Wait(TimeSpan.FromSeconds(30)), $"{input} did not finish in thirty seconds"); + return task.Result; + } + } +}