From bd6baab14d5a9549daadb4b9fa706947f8460815 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 3 Nov 2026 10:05:51 +0000 Subject: [PATCH 1/2] Form a complex product past the exponent range before rounding it A complex number's parts are decimals in the precision context, whose exponents end at 1e1000. 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: (3+4i)^2000, (1+i)^(10^6) and (1e600+1e600i)^2, and Simplify returned the NaN as well. 2.6.0 does the same. A complex product whose parts leave the range is formed with unbounded exponents and each part rounded into the context, and a complex power that does is raised the same way. Each part past the range is the infinity of its sign, as a real past it is +oo; within the range nothing changes, since the bounded result is used whenever it is finite. (1+i)^(10^6) is +oo and (1e600+1e600i)^2 is +oo i. The guard against an infinity only rounding made now covers any sum, product or power of exact finite operands (a power's base not zero), which is exact and finite, so an exact simplification leaves (3+4i)^2000 as written. An operand that is already a decimal is not covered: the integrator reads a decimal sample past the range as +oo, and taking that away made one Rubi antiderivative 23,030 characters long. Co-Authored-By: Claude Opus 5.5 --- BREAKING-CHANGES.md | 16 ++++ .../Entity/Continuous/Number/Operators.cs | 35 ++++++- .../Evaluation/Evaluation.Classes.cs | 23 +++-- .../Core/AComplexOverflowIsNotNaNTest.cs | 91 +++++++++++++++++++ 4 files changed, 154 insertions(+), 11 deletions(-) create mode 100644 Sources/Tests/UnitTests/Core/AComplexOverflowIsNotNaNTest.cs 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..6519ee03b 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs @@ -151,10 +151,30 @@ 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)); } ); + /// + /// 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 +635,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 (built.IsFinite || @base is Real) + 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 diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs index 487682612..3e80684d4 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 { IsFinite: false } && !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; + } + } +} From dc07fd1abdd825e6cab44170c4938f89da44f24e Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 3 Nov 2026 10:47:37 +0000 Subject: [PATCH 2/2] Read finiteness from a number's parts where every evaluation passes IsFinite is a cache made on an entity's first read of it, and each number an evaluation or a power produces is a new entity. Asked of every exact value in the exact-evaluation shortcut (the guard #2347 added) and of every power built (the check this branch added), it cost about 250 KB and 237 KB on each solve of the kernel benchmark's SolveHard, +5.1% against the gate's baseline. The parts answer the same question without the cache, a rational part being finite as it is. The size bound #2347 added is computed in longs as well; it measured free as EIntegers, but it is asked of every exact power. SolveHard allocates 9,549,209 B, against 9,546,434 B before #2347 and 10,027,870 B before this commit. Co-Authored-By: Claude Opus 5.5 --- .../Entity/Continuous/Number/Operators.cs | 39 ++++++++++++------- .../Evaluation/Evaluation.Classes.cs | 2 +- 2 files changed, 26 insertions(+), 15 deletions(-) diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs index 6519ee03b..1b9651440 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs @@ -158,6 +158,16 @@ static EDecimal ModifiedMultiply(EDecimal a, EDecimal b) => } ); + /// + /// 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, @@ -636,7 +646,7 @@ static Complex BinaryIntPow(Complex num, EInteger val) return Real.CreateAtScale(pow.Sign < 0 ? EDecimal.One.Divide(raised, wholeContext) : raised, EDecimal.Zero); } var built = BinaryIntPow(@base, pow); - if (built.IsFinite || @base is Real) + 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 @@ -746,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 @@ -755,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 3e80684d4..4f8643b92 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Classes.cs @@ -65,7 +65,7 @@ private static bool IsAZeroOnlyRoundingMade(Entity node, Entity value) /// 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(IsFiniteRoundedOntoAnInfinity); + => 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