From 330c51a59ccb4eb35e216faf7c8b9ced2c5357a7 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 22 Oct 2026 17:36:36 +0000 Subject: [PATCH] Round a computed value onto zero only where its terms cancelled A number an evaluation produced was rounded onto zero below half the working digits outright, 1e-50 at a hundred, whatever produced it. So e^(-155) was 0 at a hundred digits, the modulus of a complex number near 1e-136 was 0 and a quotient by it NaN, and the derivative of a right antiderivative of 1/sqrt(a + i a sinh(c + d x)) at a large argument was NaN at up to two hundred digits and 1.5329 at three hundred, where it is 2.0185e-68 (1 - i). A residual of cancellation is a difference below half the working digits of its terms, and that is now the only value rounded onto zero by size: a sum of opposite signs is measured against one of its terms, and a product, a quotient, a power, a root, an exponential or a modulus, which cancel nothing, is never rounded onto zero; a cosine's or a sine's residual is measured against the modulus it multiplies. sin(pi), e^(i pi) + 1 and sqrt(2) sqrt(3) - sqrt(6) are still 0, and #561's matrix times its inverse is still the identity. A quotient of decimals is divided to the working precision's places, as it was, and divided again to its digits where that lost them: 3.5e-68 kept thirty-three digits, and a quotient below 1e-100 was 0. The derivative is right from two hundred digits. Its terms cancel in their sixty-eighth digit, so below about twice that the difference is under half the digits and reads as a residual: at a hundred it is the real part without the imaginary one, where it was NaN. The scale is either term, which costs nothing, where comparing the two rescaled a mantissa: allocation per evaluation is master's, measured on the kernel gate's polynomial at a hundred and at fifteen digits. Part of #1944. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 28 +++++ .../Entity.Continuous.Complex.Definition.cs | 35 ++++-- .../Entity.Continuous.Real.Definition.cs | 26 +++- .../Entity/Continuous/Number/Operators.cs | 118 ++++++++++++------ ...AResidualIsZeroOnlyWhereItCancelledTest.cs | 99 +++++++++++++++ 5 files changed, 252 insertions(+), 54 deletions(-) create mode 100644 Sources/Tests/UnitTests/Core/AResidualIsZeroOnlyWhereItCancelledTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 47f838101..9eabfbfad 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -110,6 +110,34 @@ is rounded as before: `(Entity)(1m / 3)` is still `1/3` ([#1944](https://github. | The derivative of the antiderivative of `1/sqrt(a + i*a*sinh(c + d*x))` at `a, c, d, x = 21/10, 19/10, 17/10, 227283/1250`, at 2000 digits | `2.0185464…E-68` | `2.0185464…E-68 - 2.0185464…E-68i` | | `"1/3 + 1e-60".EvalNumerical()`, `"sin(pi)".EvalNumerical()`, `(Entity)(1.0 / 3)` | `1/3`, `0`, `1/3` | `1/3`, `0`, `1/3` | +### A computed value is zero only where its terms cancelled + +**Silent: wrong answers, and values where there was NaN.** A number an evaluation produced was rounded onto zero below +half the working digits outright, `1e-50` at the default hundred. That held for any number, whatever produced it. So +`e^(-155)` was `0` at a hundred digits, the modulus of a complex number near `1e-136` was `0`, and a quotient by that +modulus was `NaN`. The derivative of a right antiderivative of `1/sqrt(a + i a sinh(c + d x))` at a large argument was +therefore `NaN` at up to two hundred digits and `1.5329…` at three hundred, where it is `2.0185…E-68 - 2.0185…E-68i`. + +A number is now rounded onto zero only as a residual of the cancellation that produced it: a difference below half the +working digits of its larger term. A product, a quotient, a power, a root, an exponential or a modulus cancels nothing +and is not rounded onto zero, and a cosine's or a sine's residual is measured against the modulus it multiplies. +`sin(pi)`, `e^(i pi) + 1` and `sqrt(2) sqrt(3) - sqrt(6)` are still `0`. A quotient of decimals is rounded to the +working precision's number of decimal places, as it was, and one that has lost digits that way is divided again to the +precision's digits: a quotient of `3.5e-68` kept thirty-three digits, and one below `1e-100` was `0`. The derivative above is right from two hundred digits; its terms cancel in their sixty-eighth digit, +so at fewer than about twice that many digits the difference is below half of them and reads as a residual +([#1944](https://github.com/asc-community/AngouriMath/issues/1944)). + +| Input, at the default hundred digits unless stated | Was (2.6.0) | Now | +|---|---|---| +| `"e^(-155)".EvalNumerical()` | `0` | `4.8345416380533358…E-68` | +| `"e^(-60) * e^(-60)".EvalNumerical()` | `0` | `7.6676480737219996…E-53` | +| `"1/(10^67*sqrt(2) + 10^67*sqrt(2)*i)".EvalNumerical()` | `0` | `3.5355339059327376…E-68 - 3.5355339059327376…E-68i` | +| `"sqrt(-e^(-300))".EvalNumerical()` | `0` | `7.1750959731644104…E-66i` | +| The derivative above, at 300 digits | `1.5329402…` | `2.0185464…E-68 - 2.0185464…E-68i` | +| The derivative above, at 200 digits | `NaN` | `2.0185464…E-68 - 2.0185464…E-68i` | +| The derivative above, at 100 digits | `NaN` | `2.0185464…E-68`, without its imaginary part: a residual at that precision | +| `"sin(pi)".EvalNumerical()`, `"(1 + sqrt(2) i) (1 - sqrt(2) i)".EvalNumerical()` | `0`, `3` | `0`, `3` | + --- ## 2.6.0 — since 2.5.0 diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Complex.Definition.cs b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Complex.Definition.cs index 1413254ec..766d4027a 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Complex.Definition.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Complex.Definition.cs @@ -80,27 +80,36 @@ public static Complex Create(Real real, Real imaginary) => /// /// Creates an instance of Complex /// - public static Complex Create(EDecimal real, EDecimal imaginary) => Create(real, imaginary, fromMachineNumbers: false); + public static Complex Create(EDecimal real, EDecimal imaginary) => Create(real, imaginary, EDecimal.One, EDecimal.One, fromMachineNumbers: false); + + /// + /// A number an arithmetic operation produced, each part zero only as a residual of the cancellation that + /// produced it, as reads one. + /// https://github.com/asc-community/AngouriMath/issues/1944 + /// + internal static Complex CreateAtScales(EDecimal real, EDecimal imaginary, EDecimal realScale, EDecimal imaginaryScale) + => Create(real, imaginary, realScale, imaginaryScale, fromMachineNumbers: false); /// /// A number from a pair of s, each part rounded as /// rounds it. /// - internal static Complex FromMachineNumbers(EDecimal real, EDecimal imaginary) => Create(real, imaginary, fromMachineNumbers: true); + internal static Complex FromMachineNumbers(EDecimal real, EDecimal imaginary) => Create(real, imaginary, EDecimal.One, EDecimal.One, fromMachineNumbers: true); - private static Complex Create(EDecimal real, EDecimal imaginary, bool fromMachineNumbers) + private static Complex Create(EDecimal real, EDecimal imaginary, EDecimal realScale, EDecimal imaginaryScale, bool fromMachineNumbers) { if (!MathS.Settings.DowncastingEnabled) return new Complex(real, imaginary); - Real Part(EDecimal part) => fromMachineNumbers ? Real.FromMachineNumber(part) : Real.Create(part); + Real Part(EDecimal part, EDecimal scale) => fromMachineNumbers ? Real.FromMachineNumber(part) : Real.CreateAtScale(part, scale); if (imaginary.IsZero) - return Part(real); + return Part(real, realScale); if (real.IsNaN() || imaginary.IsNaN()) return Real.NaN; - if (imaginary.IsFinite && imaginary.Abs().LessThan(MathS.Settings.DowncastingTolerance)) - return Part(real); + if (imaginary.IsFinite && !imaginaryScale.IsZero && imaginary.Abs().LessThan( + ReferenceEquals(imaginaryScale, EDecimal.One) ? MathS.Settings.DowncastingTolerance : MathS.Settings.DowncastingTolerance.Multiply(imaginaryScale.Abs()))) + return Part(real, realScale); else - return new Complex(Part(real), Part(imaginary)); + return new Complex(Part(real, realScale), Part(imaginary, imaginaryScale)); } /// @@ -122,8 +131,10 @@ public void Deconstruct(out Real realPart, out Real imaginaryPart) => /// https://github.com/asc-community/AngouriMath/issues/1615 /// public new virtual Real Abs() - => Real.Create((RealPart.EDecimal * RealPart.EDecimal + ImaginaryPart.EDecimal * ImaginaryPart.EDecimal) - .SqrtByIntegerRoot(MathS.Settings.DecimalPrecisionContext)); + // A modulus is never a residual of cancellation: the sum of two squares cancels nothing. + // https://github.com/asc-community/AngouriMath/issues/1944 + => Real.CreateAtScale((RealPart.EDecimal * RealPart.EDecimal + ImaginaryPart.EDecimal * ImaginaryPart.EDecimal) + .SqrtByIntegerRoot(MathS.Settings.DecimalPrecisionContext), EDecimal.Zero); /// /// The phase of a complex number (aka angle) @@ -136,7 +147,9 @@ public void Deconstruct(out Real realPart, out Real imaginaryPart) => public static Complex CreatePolar(EDecimal magnitude, EDecimal phase) { var context = MathS.Settings.DecimalPrecisionContext; - return Create(magnitude.Multiply(phase.Cos(context), context), magnitude.Multiply(phase.Sin(context), context)); + // The residual of a cosine or a sine is the magnitude's: cos(pi/2) leaves 1e-100 of it. + var scale = magnitude.Abs(); + return CreateAtScales(magnitude.Multiply(phase.Cos(context), context), magnitude.Multiply(phase.Sin(context), context), scale, scale); } /// -oo + -ooi diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Real.Definition.cs b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Real.Definition.cs index ba6d381cb..c723d7a23 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Real.Definition.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Real.Definition.cs @@ -65,7 +65,18 @@ private protected Real() : base(null, null) { } /// such logarithms lost its imaginary part. /// https://github.com/asc-community/AngouriMath/issues/1944 /// - public static Real Create(EDecimal value) => Create(value, fromMachineNumber: false); + public static Real Create(EDecimal value) => Create(value, fromMachineNumber: false, EDecimal.One); + + /// + /// A value an arithmetic operation produced, which is a residual of cancellation, and zero, only below + /// half the working digits of , the size of the terms that cancelled; + /// zero for a value no cancellation produced -- a product, a quotient or a power of nonzero values -- + /// which is never rounded onto zero. With the floor at half the working digits outright, every value + /// below 1e-50 was zero at a hundred digits, e^(-155) among them, and a modulus near 1e-136 was zero, so + /// that a quotient by it was NaN. https://github.com/asc-community/AngouriMath/issues/1944 + /// + internal static Real CreateAtScale(EDecimal value, EDecimal cancelledScale) + => Create(value, fromMachineNumber: false, cancelledScale); /// /// A number from a , a or a , which carries @@ -73,11 +84,11 @@ private protected Real() : base(null, null) { } /// , ten to the minus sixteen, so that 1m/3 /// is 1/3. https://github.com/asc-community/AngouriMath/issues/1944 /// - internal static Real FromMachineNumber(EDecimal value) => Create(value, fromMachineNumber: true); + internal static Real FromMachineNumber(EDecimal value) => Create(value, fromMachineNumber: true, EDecimal.One); // The tolerance is made only past the cheap check below, which almost every value stops at: made at // the call, it was 580 more bytes on each evaluation of a cubic at a hundred digits. - private static Real Create(EDecimal value, bool fromMachineNumber) + private static Real Create(EDecimal value, bool fromMachineNumber, EDecimal cancelledScale) { if (!MathS.Settings.DowncastingEnabled) return new Real(value); @@ -91,11 +102,14 @@ private static Real Create(EDecimal value, bool fromMachineNumber) if (!Rational.MayBeASmallRational(value, MathS.Settings.FloatToRationalIterCount, MathS.Settings.DowncastingToleranceAsDouble)) return new Real(value); var (intPart, intRest) = value.SplitDecimal(); - // If the difference between value & round(value) is zero (see Number.IsZero), we consider value as an integer + // If the difference between value & round(value) is zero (see Number.IsZero), we consider value as an integer. + // Onto zero only as a residual of the cancellation that produced it, below the tolerance times the + // scale of what cancelled; a negative one reaches zero by rounding up. var tolerance = MathS.Settings.DowncastingTolerance; - if (intRest.LessThan(tolerance)) + bool IsAResidual() => value.IsZero || !cancelledScale.IsZero && value.Abs().LessThan(tolerance.Multiply(cancelledScale.Abs())); + if (intRest.LessThan(tolerance) && (!intPart.IsZero || IsAResidual())) return Integer.Create(intPart); - if (intRest.GreaterThan(1 - tolerance)) + if (intRest.GreaterThan(1 - tolerance) && (!intPart.Increment().IsZero || IsAResidual())) return Integer.Create(intPart.Increment()); var attempt = Rational.FindRationalWithin(value, fromMachineNumber ? MathS.Settings.PrecisionErrorZeroRange : tolerance); diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs index da379c455..08e992fee 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs @@ -34,13 +34,41 @@ internal static EDecimal CtxSubtract(EDecimal a, EDecimal b) => a.Subtract(b, MathS.Settings.DecimalPrecisionContext); internal static EDecimal CtxMultiply(EDecimal a, EDecimal b) => a.Multiply(b, MathS.Settings.DecimalPrecisionContext); + // To the working precision's number of decimal places, which is the precision's digits for a quotient + // near one and cheaper than dividing to digits; a quotient that has visibly lost digits that way is + // divided again to them: to places, 3.5e-68 kept thirty-three digits at a hundred, and a quotient below + // 1e-100 was zero. https://github.com/asc-community/AngouriMath/issues/1944 internal static EDecimal CtxDivide(EDecimal a, EDecimal b) - => b.IsZero ? EDecimal.NaN : a.DivideToExponent(b, -MathS.Settings.DecimalPrecisionContext.Value.Precision); + { + if (b.IsZero) + return EDecimal.NaN; + var context = MathS.Settings.DecimalPrecisionContext.Value; + var quotient = a.DivideToExponent(b, -context.Precision); + // 3.3 bits a digit, and five digits to spare; an unlimited context keeps the division to places. + if (!context.Precision.IsZero && quotient.IsFinite && !a.IsZero + && quotient.UnsignedMantissa.GetUnsignedBitLengthAsInt64() < (context.Precision.ToInt64Checked() - 5) * 3) + return a.Divide(b, context); + return quotient; + } internal static EDecimal CtxMod(EDecimal a, EDecimal b) => a.Remainder(b, MathS.Settings.DecimalPrecisionContext); internal static EDecimal CtxPow(EDecimal a, EDecimal b) => a.Pow(b, MathS.Settings.DecimalPrecisionContext); internal static T Min(T a, T b) where T : Real => a < b ? a : b; + + /// + /// The scale of what cancels in a + b: a term where their signs differ, and zero where they cannot + /// cancel, so that a sum of like signs is never a residual. + /// https://github.com/asc-community/AngouriMath/issues/1944 + /// + internal static EDecimal CancelledScale(EDecimal a, EDecimal b) + // Either term: where they nearly cancel they are alike, and where they do not the sum is near the + // larger and no residual of either. Comparing them rescaled a mantissa, 500 bytes at a hundred digits. + => a.Sign * b.Sign < 0 ? a : EDecimal.Zero; + + /// The scale of what cancels in a - b, as . + internal static EDecimal CancelledScaleOfADifference(EDecimal a, EDecimal b) + => a.Sign * b.Sign > 0 ? a : EDecimal.Zero; internal static T Max(T a, T b) where T : Real => a > b ? a : b; /// @@ -86,25 +114,29 @@ internal static T OpSum(T a, T b) where T : Number => SuperSwitch(a, b, (a, b) => Integer.Create(CtxAdd(a.EInteger, b.EInteger)), (a, b) => Rational.Create(CtxAdd(a.ERational, b.ERational)), - (a, b) => Real.Create(CtxAdd(a.EDecimal, b.EDecimal)), - (a, b) => Complex.Create(CtxAdd(a.RealPart.EDecimal, b.RealPart.EDecimal), - CtxAdd(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal)) + (a, b) => Real.CreateAtScale(CtxAdd(a.EDecimal, b.EDecimal), CancelledScale(a.EDecimal, b.EDecimal)), + (a, b) => Complex.CreateAtScales(CtxAdd(a.RealPart.EDecimal, b.RealPart.EDecimal), + CtxAdd(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal), + CancelledScale(a.RealPart.EDecimal, b.RealPart.EDecimal), + CancelledScale(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal)) ); internal static T OpSub(T a, T b) where T : Number => SuperSwitch(a, b, (a, b) => Integer.Create(CtxSubtract(a.EInteger, b.EInteger)), (a, b) => Rational.Create(CtxSubtract(a.ERational, b.ERational)), - (a, b) => Real.Create(CtxSubtract(a.EDecimal, b.EDecimal)), - (a, b) => Complex.Create(CtxSubtract(a.RealPart.EDecimal, b.RealPart.EDecimal), - CtxSubtract(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal)) + (a, b) => Real.CreateAtScale(CtxSubtract(a.EDecimal, b.EDecimal), CancelledScaleOfADifference(a.EDecimal, b.EDecimal)), + (a, b) => Complex.CreateAtScales(CtxSubtract(a.RealPart.EDecimal, b.RealPart.EDecimal), + CtxSubtract(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal), + CancelledScaleOfADifference(a.RealPart.EDecimal, b.RealPart.EDecimal), + CancelledScaleOfADifference(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal)) ); internal static T OpMul(T a, T b) where T : Number => SuperSwitch(a, b, (a, b) => Integer.Create(CtxMultiply(a.EInteger, b.EInteger)), (a, b) => Rational.Create(CtxMultiply(a.ERational, b.ERational)), - (a, b) => Real.Create(CtxMultiply(a.EDecimal, b.EDecimal)), + (a, b) => Real.CreateAtScale(CtxMultiply(a.EDecimal, b.EDecimal), EDecimal.Zero), (a, b) => { // A whole zero against a whole infinity is 0 * oo, which is NaN for a real @@ -117,9 +149,9 @@ internal static T OpMul(T a, T b) where T : Number => // Define both (oo * i) and (i * oo) to be (oo i) which is (0 + oo i) instead of (NaN + oo i) static EDecimal ModifiedMultiply(EDecimal a, EDecimal b) => a.IsInfinity() && b.IsZero || b.IsInfinity() && a.IsZero ? EDecimal.Zero : CtxMultiply(a, b); - return Complex.Create( - CtxSubtract(ModifiedMultiply(a.RealPart.EDecimal, b.RealPart.EDecimal), ModifiedMultiply(a.ImaginaryPart.EDecimal, b.ImaginaryPart.EDecimal)), - CtxAdd(ModifiedMultiply(a.RealPart.EDecimal, b.ImaginaryPart.EDecimal), ModifiedMultiply(a.ImaginaryPart.EDecimal, b.RealPart.EDecimal))); + 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)); } ); @@ -127,7 +159,7 @@ 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, (a, b) => CtxDivide(a.ERational, b.ERational) is { } n ? Rational.Create(n) : Real.NaN, - (a, b) => Real.Create(CtxDivide(a.EDecimal, b.EDecimal)), + (a, b) => Real.CreateAtScale(CtxDivide(a.EDecimal, b.EDecimal), EDecimal.Zero), (a, b) => { /* @@ -158,7 +190,7 @@ internal static Complex OpDiv(T a, T b) where T : Number => var abs2 = CtxMultiply(bAbs, bAbs); var Re = CtxDivide(conj.RealPart.EDecimal, abs2); var Im = CtxDivide(conj.ImaginaryPart.EDecimal, abs2); - var c = Complex.Create(Re, Im); + var c = Complex.CreateAtScales(Re, Im, EDecimal.Zero, EDecimal.Zero); return a * c; } ); @@ -172,28 +204,35 @@ private static Complex DivideBySmith(Complex a, Complex b) { var (ar, ai) = (a.RealPart.EDecimal, a.ImaginaryPart.EDecimal); var (c, d) = (b.RealPart.EDecimal, b.ImaginaryPart.EDecimal); - // Rounded to the context's digits at the end: CtxDivide divides to a number of - // decimal places. var context = MathS.Settings.DecimalPrecisionContext.Value; + // A part is a residual only of the sum in its numerator, at that sum's scale over the + // denominator's; the denominator is a sum of like signs, and a quotient cancels nothing. + EDecimal Over(EDecimal scale, EDecimal denominator) => scale.IsZero ? scale : scale.Abs().Divide(denominator.Abs(), context); if (d.IsZero) - return Complex.Create(CtxDivide(ar, c).RoundToPrecision(context), CtxDivide(ai, c).RoundToPrecision(context)); + return Complex.CreateAtScales(CtxDivide(ar, c), CtxDivide(ai, c), EDecimal.Zero, EDecimal.Zero); if (c.IsZero) - return Complex.Create(CtxDivide(ai, d).RoundToPrecision(context), CtxDivide(ar.Negate(), d).RoundToPrecision(context)); + return Complex.CreateAtScales(CtxDivide(ai, d), CtxDivide(ar.Negate(), d), EDecimal.Zero, EDecimal.Zero); if (c.Abs().CompareTo(d.Abs()) >= 0) { var ratio = CtxDivide(d, c); var denominator = CtxAdd(c, CtxMultiply(d, ratio)); - return Complex.Create( - CtxDivide(CtxAdd(ar, CtxMultiply(ai, ratio)), denominator).RoundToPrecision(context), - CtxDivide(CtxSubtract(ai, CtxMultiply(ar, ratio)), denominator).RoundToPrecision(context)); + var (aiRatio, arRatio) = (CtxMultiply(ai, ratio), CtxMultiply(ar, ratio)); + return Complex.CreateAtScales( + CtxDivide(CtxAdd(ar, aiRatio), denominator), + CtxDivide(CtxSubtract(ai, arRatio), denominator), + Over(CancelledScale(ar, aiRatio), denominator), + Over(CancelledScaleOfADifference(ai, arRatio), denominator)); } else { var ratio = CtxDivide(c, d); var denominator = CtxAdd(CtxMultiply(c, ratio), d); - return Complex.Create( - CtxDivide(CtxAdd(CtxMultiply(ar, ratio), ai), denominator).RoundToPrecision(context), - CtxDivide(CtxSubtract(CtxMultiply(ai, ratio), ar), denominator).RoundToPrecision(context)); + var (arRatio, aiRatio) = (CtxMultiply(ar, ratio), CtxMultiply(ai, ratio)); + return Complex.CreateAtScales( + CtxDivide(CtxAdd(arRatio, ai), denominator), + CtxDivide(CtxSubtract(aiRatio, ar), denominator), + Over(CancelledScale(arRatio, ai), denominator), + Over(CancelledScaleOfADifference(aiRatio, ar), denominator)); } } @@ -384,8 +423,8 @@ public static Complex Sqrt(Complex num) // From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,7dc9c2ee4f99814a if (num is Real { EDecimal: var real }) if (real.IsNegative) - return Complex.Create(0, real.Negate().SqrtByIntegerRoot(context)); - else return Real.Create(real.SqrtByIntegerRoot(context)); + return Complex.CreateAtScales(EDecimal.Zero, real.Negate().SqrtByIntegerRoot(context), EDecimal.Zero, EDecimal.Zero); + else return Real.CreateAtScale(real.SqrtByIntegerRoot(context), EDecimal.Zero); else { @@ -432,7 +471,8 @@ public static Complex Sqrt(Complex num) if (im.Sign < 0) y = -y; x = im.Divide(y.Multiply(2, context), context); } - return Complex.Create(x, y); + // Each part is a root of a sum of like signs or a quotient by one, and cancels nothing. + return Complex.CreateAtScales(x, y, EDecimal.Zero, EDecimal.Zero); } } @@ -443,11 +483,12 @@ public static Complex Exp(Complex x) var expReal = x.RealPart.EDecimal.Exponential(context); var imaginary = x.ImaginaryPart.EDecimal; if (imaginary.IsZero) - return expReal; + return Real.CreateAtScale(expReal, EDecimal.Zero); // From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,7dc9c2ee4f99814a var cosImaginary = expReal.Multiply(imaginary.Cos(context), context); var sinImaginary = expReal.Multiply(imaginary.Sin(context), context); - return Complex.Create(cosImaginary, sinImaginary); + // The residual of a cosine or a sine is the modulus's: e^(i pi) leaves 1e-100 i of it. + return Complex.CreateAtScales(cosImaginary, sinImaginary, expReal, expReal); } /// e.g. Pow(2, 5) = 32 @@ -490,7 +531,7 @@ static Complex BinaryIntPow(Complex num, EInteger val) // where the power was a logarithm and a series. if (@base is Real { EDecimal: var maybeE } && power is Real { EDecimal: var realExponent } && maybeE.Equals(InternalAMExtensions.ConstantCache.Lookup(MathS.Settings.DecimalPrecisionContext).E)) - return Real.Create(realExponent.Exponential(MathS.Settings.DecimalPrecisionContext)); + return Real.CreateAtScale(realExponent.Exponential(MathS.Settings.DecimalPrecisionContext), EDecimal.Zero); // TODO: make it more detailed (e. g. +oo ^ +oo = +oo) // A whole power, whether written as an integer or, for a base off the // positive axis, as a decimal with nothing after the point -- which is how `2` @@ -530,7 +571,7 @@ static Complex BinaryIntPow(Complex num, EInteger val) { var wholeContext = MathS.Settings.DecimalPrecisionContext; var raised = decimalBase.Pow(pow.Abs().ToInt32Checked(), wholeContext); - return Real.Create(pow.Sign < 0 ? EDecimal.One.Divide(raised, wholeContext) : raised); + return Real.CreateAtScale(pow.Sign < 0 ? EDecimal.One.Divide(raised, wholeContext) : raised, EDecimal.Zero); } return BinaryIntPow(@base, pow); } @@ -546,7 +587,7 @@ static Complex BinaryIntPow(Complex num, EInteger val) var root = rootBase.SqrtByIntegerRoot(halfContext); var n = halfPower.Numerator.Abs().ToInt32Checked(); var raised = n == 1 ? root : root.Pow(n, halfContext); - return Real.Create(halfPower.Numerator.Sign < 0 ? EDecimal.One.Divide(raised, halfContext) : raised); + return Real.CreateAtScale(halfPower.Numerator.Sign < 0 ? EDecimal.One.Divide(raised, halfContext) : raised, EDecimal.Zero); } // A negative real to a real power that is a whole number of halves -- a half // power above all -- is the modulus to that power times a power of i, exactly: @@ -574,20 +615,21 @@ static Complex BinaryIntPow(Complex num, EInteger val) && rootAbove.Multiply(rootAbove).Equals(exactBase.Numerator.Abs()) && rootBelow.Multiply(rootBelow).Equals(exactBase.Denominator)) modulus = Pow(Rational.Create(rootAbove, rootBelow), Integer.Create(numeratorOfPower)); else if (numeratorOfPower.IsEven) - modulus = Pow(Real.Create(negativeBase.Negate()), Integer.Create(numeratorOfPower.Divide(2))); + modulus = Pow(Real.CreateAtScale(negativeBase.Negate(), EDecimal.Zero), Integer.Create(numeratorOfPower.Divide(2))); else { var root = negativeBase.Negate().SqrtByIntegerRoot(halfContext); var n = numeratorOfPower.Abs().ToInt32Checked(); var raised = n == 1 ? root : root.Pow(n, halfContext); - modulus = Real.Create(numeratorOfPower.Sign < 0 ? EDecimal.One.Divide(raised, halfContext) : raised); + modulus = Real.CreateAtScale(numeratorOfPower.Sign < 0 ? EDecimal.One.Divide(raised, halfContext) : raised, EDecimal.Zero); } return quarterTurns.ToInt32Checked() switch { 0 => modulus, - 1 => Complex.Create(Integer.Zero, modulus.RealPart), - 2 => Real.Create(modulus.RealPart.EDecimal.Negate()), - _ => Complex.Create(Integer.Zero, Real.Create(modulus.RealPart.EDecimal.Negate())), + // A turn of the modulus, which is not a residual of anything. + 1 => Complex.CreateAtScales(EDecimal.Zero, modulus.RealPart.EDecimal, EDecimal.Zero, EDecimal.Zero), + 2 => Real.CreateAtScale(modulus.RealPart.EDecimal.Negate(), EDecimal.Zero), + _ => Complex.CreateAtScales(EDecimal.Zero, modulus.RealPart.EDecimal.Negate(), EDecimal.Zero, EDecimal.Zero), }; } if (@base.IsFinite && power is Rational r && r.ERational.Denominator.Abs() < 10 // there should be a minimal threshold to avoid long searches @@ -616,7 +658,9 @@ static Complex BinaryIntPow(Complex num, EInteger val) var newRho = powerReal.MultiplyAndAdd(theta, powerImaginary.Multiply(rho.NaturalLogarithm(context), context), context); var t = PowerOfReals(rho, powerReal, context).Multiply(powerImaginary.Multiply(-theta, context).Exponential(context), context); - return Complex.Create(t.Multiply(newRho.Cos(context), context), t.Multiply(newRho.Sin(context), context)); + // The residual of a cosine or a sine is the modulus's. + var modulusScale = t.Abs(); + return Complex.CreateAtScales(t.Multiply(newRho.Cos(context), context), t.Multiply(newRho.Sin(context), context), modulusScale, modulusScale); } /// e.g. Log(2, 32) = 5 diff --git a/Sources/Tests/UnitTests/Core/AResidualIsZeroOnlyWhereItCancelledTest.cs b/Sources/Tests/UnitTests/Core/AResidualIsZeroOnlyWhereItCancelledTest.cs new file mode 100644 index 000000000..a585ea048 --- /dev/null +++ b/Sources/Tests/UnitTests/Core/AResidualIsZeroOnlyWhereItCancelledTest.cs @@ -0,0 +1,99 @@ +// +// 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 AngouriMath.Extensions; +using PeterO.Numbers; +using Xunit; +using static AngouriMath.Entity.Number; + +namespace AngouriMath.Tests.Core +{ + /// + /// A number an evaluation produces is rounded onto zero only as a residual of the cancellation that produced it: + /// a difference below half the working digits of its larger term. A product, a quotient, a power, a root, an + /// exponential or a modulus cancels nothing and is never zero by size. Below half the working digits outright, + /// e^(-155) was 0 at a hundred digits, a modulus near 1e-136 was 0 and a quotient by it NaN, and the + /// derivative of a right antiderivative at a large hyperbolic argument was NaN at a hundred digits and wrong at + /// three hundred. References from mpmath at 120 digits. + /// #1944 + /// + [Trait("Area", "Core")] + public sealed class AResidualIsZeroOnlyWhereItCancelledTest + { + private static void AssertClose(string reference, EDecimal actual) + { + var expected = EDecimal.FromString(reference); + var error = actual.Subtract(expected, EContext.ForPrecision(200)).Abs().Divide(expected.Abs(), EContext.ForPrecision(10)); + Assert.True(error.LessThan(EDecimal.FromString("1E-40")), $"{actual} where it is {expected}"); + } + + [Theory] + [InlineData("e^(-155)", "4.83454163805333583428273878250224602960536796E-68")] + [InlineData("e^(-120)", "7.66764807372199963243420837527582768233983151E-53")] + [InlineData("e^(-60) * e^(-60)", "7.66764807372199963243420837527582768233983151E-53")] + // e^(-200) e^(i pi): the imaginary part is the residual of a sine, at the modulus's scale. + [InlineData("e^(-200 + i*pi)", "-1.38389652673673753064868145697908468540304758233947720939393E-87")] + public void AProductOrAPowerIsNotZeroBySize(string expression, string reference) + => AssertClose(reference, Assert.IsType(expression.EvalNumerical()).EDecimal); + + [Fact] + public void AQuotientKeepsItsDigitsAndItsImaginaryPart() + { + var value = Assert.IsType("1/(10^67*sqrt(2) + 10^67*sqrt(2)*i)".EvalNumerical()); + AssertClose("3.53553390593273762200422181052424519642417969E-68", value.RealPart.EDecimal); + AssertClose("-3.53553390593273762200422181052424519642417969E-68", value.ImaginaryPart.EDecimal); + } + + [Fact] + public void TheRootOfASmallNegativeIsImaginary() + { + var value = Assert.IsType("sqrt(-e^(-300))".EvalNumerical()); + Assert.True(value.RealPart.EDecimal.IsZero, $"{value}"); + AssertClose("7.17509597316441041983269290720898818840868827424412900583245E-66", value.ImaginaryPart.EDecimal); + } + + [Theory] + [InlineData("sin(pi)", "0")] + [InlineData("cos(pi/2)", "0")] + [InlineData("e^(i*pi) + 1", "0")] + [InlineData("sin(100*pi)", "0")] + [InlineData("sqrt(2)*sqrt(3) - sqrt(6)", "0")] + [InlineData("(1 + sqrt(2)*i) * (1 - sqrt(2)*i)", "3")] + // A difference of e^(-155) and a value 1e-90 of it larger: below half the digits of the terms. + [InlineData("e^(-155) - e^(-155) * (1 + 10^(-90))", "0")] + public void AResidualOfCancellationIsStillZero(string expression, string exact) + => Assert.Equal(exact.ToEntity(), expression.EvalNumerical()); + + [Theory] + [InlineData(200)] + [InlineData(300)] + public void TheDerivativeOfAnAntiderivativeAtALargeHyperbolicArgumentIsTheIntegrand(int digits) + { + // The derivative's terms cancel in their sixty-eighth digit, so below twice that it reads as a residual. + var integrand = "1/sqrt(a + i*a*sinh(c + d*x))".ToEntity(); + var derivative = integrand.Integrate("x").Differentiate("x"); + Entity Pin(Entity e) => e.Substitute("a", Rational.Create(21, 10)).Substitute("c", Rational.Create(19, 10)) + .Substitute("d", Rational.Create(17, 10)).Substitute("x", Rational.Create(227283, 1250)); + using var _ = MathS.Settings.DecimalPrecisionContext.Set(new EContext(digits, ERounding.HalfEven, -100000, 100000, false)); + var value = Assert.IsType(Pin(derivative).EvalNumerical()); + AssertClose("2.01854640568606010538100888350636282143172892E-68", value.RealPart.EDecimal); + AssertClose("-2.01854640568606010538100888350636282143172892E-68", value.ImaginaryPart.EDecimal); + } + + [Theory] + [InlineData(50)] + [InlineData(100)] + public void TheIntegrandAtALargeHyperbolicArgumentHasItsValue(int digits) + { + var integrand = "1/sqrt(21/10 + i*21/10*sinh(19/10 + 17/10*227283/1250))".ToEntity(); + using var _ = MathS.Settings.DecimalPrecisionContext.Set(new EContext(digits, ERounding.HalfEven, -100000, 100000, false)); + var value = Assert.IsType(integrand.EvalNumerical()); + AssertClose("2.01854640568606010538100888350636282143172892E-68", value.RealPart.EDecimal); + AssertClose("-2.01854640568606010538100888350636282143172892E-68", value.ImaginaryPart.EDecimal); + } + } +}