From 9fc90829960ce639a1c9a3646bba2d047f924ec0 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 22 Oct 2026 18:37:24 +0000 Subject: [PATCH] Measure a function's residual against its argument A function's value was rounded onto zero below half the working digits outright. That floor is for the residual a function leaves at its own zeros: sin(pi) evaluates to about 1e-100. That residual is relative to the argument, though, and a function of an exact small argument cancels nothing, so sin(10^(-60)), arctan(10^(-60)), erf(10^(-60)) and Si(10^(-60)) were all 0 at a hundred digits. The sine, cosine, tangent, logarithm, arcsine, arccosine and arctangent of a real argument, the error functions, and Si and Shi now create their value at the argument's scale, as #2072 made arithmetic do with what cancelled. sin(pi), sin(100 pi), cos(pi/2), tan(pi) and arccos(1) are still 0. The argument is passed as it is and its modulus taken only on the path that tests for a residual, so the kernel workloads allocate as before. sinh(10^(-60)) is still 0: it is evaluated from e^(10^(-60)), which is 1 within half the digits, a cancellation the working precision cannot resolve. Fixes #2074. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 17 +++++++ .../Continuous/Number/ErrorFunctions.cs | 7 ++- .../Entity/Continuous/Number/Operators.cs | 16 +++--- .../Number/TrigonometricIntegrals.cs | 7 +++ .../Core/AFunctionOfASmallArgumentTest.cs | 51 +++++++++++++++++++ 5 files changed, 89 insertions(+), 9 deletions(-) create mode 100644 Sources/Tests/UnitTests/Core/AFunctionOfASmallArgumentTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index fc715e3b5..ad8f9a9fd 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -172,6 +172,23 @@ so at fewer than about twice that many digits the difference is below half of th | 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` | +### A function of a small argument is its value, not zero + +**Silent: wrong answers.** A function's value was rounded onto zero below half the working digits outright, for the +residual a function leaves at its own zeros, such as `sin(pi)`'s. That residual is relative to the argument, and a +function of an exact small argument cancels nothing, so `sin(10^(-60))` was `0` at the default hundred digits. The +sine, cosine, tangent, logarithm, arcsine, arccosine and arctangent of a real argument, and the error functions and +`Si` and `Shi`, now measure their residual against their argument. `sin(pi)`, `cos(pi/2)` and `tan(pi)` are still `0` +([#2074](https://github.com/asc-community/AngouriMath/issues/2074)). + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"sin(10^(-60))".EvalNumerical()` | `0` | `1E-60`, to the hundredth digit | +| `"arctan(10^(-60))".EvalNumerical()` | `0` | `1E-60`, to the hundredth digit | +| `"erf(10^(-60))".EvalNumerical()` | `0` | `1.1283791670955125…E-60` | +| `"Si(10^(-60))".EvalNumerical()` | `0` | `1E-60`, to the hundredth digit | +| `"sin(pi)".EvalNumerical()` | `0` | `0` | + --- ## 2.6.0 — since 2.5.0 diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/ErrorFunctions.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/ErrorFunctions.cs index 013a24b83..d2483619f 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/ErrorFunctions.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/ErrorFunctions.cs @@ -53,7 +53,10 @@ public static Complex Erf(Complex z) if (re.IsZero && im.IsZero) return Integer.Zero; var (erfRe, erfIm) = Erf(re, im, context, extraDigits: 0); - return Complex.Create(erfRe.RoundToPrecision(context), erfIm.RoundToPrecision(context)); + // Its residual is its argument's: erf(1e-60) is 1.13e-60 and no residual of anything. + // https://github.com/asc-community/AngouriMath/issues/2074 + var scale = im.IsZero ? re : Magnitude((re, im)); + return Complex.CreateAtScales(erfRe.RoundToPrecision(context), erfIm.RoundToPrecision(context), scale, scale); } /// @@ -115,7 +118,7 @@ public static Complex Erfi(Complex z) if (x.IsInfinity()) return x.IsNegative ? Real.NegativeInfinity : Real.PositiveInfinity; var context = MathS.Settings.DecimalPrecisionContext.Value; - return ErfiOfAReal(x, context).RoundToPrecision(context); + return Real.CreateAtScale(ErfiOfAReal(x, context).RoundToPrecision(context), x); } var erf = Erf(Complex.Create(z.ImaginaryPart.EDecimal.Negate(), z.RealPart.EDecimal)); // -i (a + i b) = b - i a diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs index f7ebfd984..30c5620be 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs @@ -782,8 +782,10 @@ public static Complex Ln(Complex x) var context = MathS.Settings.DecimalPrecisionContext; if (LostToExponentRange(x) && x is Rational { ERational: { Sign: > 0 } ratio }) return LnOfRatio(ratio, context); + // A function's residual is its argument's: sin(pi) leaves about 1e-100 pi, where sin(1e-60) + // cancels nothing and is 1e-60. https://github.com/asc-community/AngouriMath/issues/2074 if (x is Real { EDecimal: { IsNegative: false } real }) - return real.NaturalLogarithm(context); + return Real.CreateAtScale(real.NaturalLogarithm(context), real); // From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1 return Complex.Create(x.Abs().EDecimal.NaturalLogarithm(context), x.ImaginaryPart.EDecimal.Arctan2(x.RealPart.EDecimal, context)); } @@ -793,7 +795,7 @@ public static Complex Sin(Complex num) { var context = MathS.Settings.DecimalPrecisionContext; if (num is Real { EDecimal: var real }) - return real.Sin(context); + return Real.CreateAtScale(real.Sin(context), real); // From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1 var (re, im) = (num.RealPart.EDecimal, num.ImaginaryPart.EDecimal); // We need both sinh and cosh of imaginary part. @@ -857,7 +859,7 @@ public static Complex Cos(Complex num) { var context = MathS.Settings.DecimalPrecisionContext; if (num is Real { EDecimal: var real }) - return real.Cos(context); + return Real.CreateAtScale(real.Cos(context), real); // From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1 var (re, im) = (num.RealPart.EDecimal, num.ImaginaryPart.EDecimal); var p = im.Exponential(context); @@ -874,7 +876,7 @@ public static Complex Tan(Complex num) { var context = MathS.Settings.DecimalPrecisionContext; if (num is Real { EDecimal: var real }) - return real.Tan(context); + return Real.CreateAtScale(real.Tan(context), real); // From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1 // tan z = sin z / cos z, but to avoid unnecessary repeated trig computations, use @@ -954,7 +956,7 @@ public static Complex Arcsin(Complex num) { var context = MathS.Settings.DecimalPrecisionContext; if (num is Real { EDecimal: var real } && !(real.GreaterThan(EDecimal.One) || real.LessThan(-EDecimal.One))) - return real.Arcsin(context); + return Real.CreateAtScale(real.Arcsin(context), real); var (beta, v) = ArcSinCosInner(num, context); return Complex.Create(beta.Arcsin(context), v); } @@ -964,7 +966,7 @@ public static Complex Arccos(Complex num) { var context = MathS.Settings.DecimalPrecisionContext; if (num is Real { EDecimal: var real } && !(real.GreaterThan(EDecimal.One) || real.LessThan(-EDecimal.One))) - return real.Acos(context); + return Real.CreateAtScale(real.Acos(context), real); var (beta, v) = ArcSinCosInner(num, context); return Complex.Create(beta.Acos(context), -v); } @@ -974,7 +976,7 @@ public static Complex Arctan(Complex num) { var context = MathS.Settings.DecimalPrecisionContext; if (num is Real { EDecimal: var real }) - return real.Arctan(context); + return Real.CreateAtScale(real.Arctan(context), real); // From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1 var one = Integer.One; var two = Integer.Create(2); diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/TrigonometricIntegrals.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/TrigonometricIntegrals.cs index 24c91d697..0174ac64c 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/TrigonometricIntegrals.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/TrigonometricIntegrals.cs @@ -93,6 +93,13 @@ private static Complex TrigonometricIntegral(Complex z, bool hyperbolic, bool od else valueIm = im.IsNegative ? pi.Divide(-2, context) : pi.Divide(2, context); } + // Si and Shi vanish at 0 like their argument, and their residual is its: Si(1e-60) is 1e-60. + // https://github.com/asc-community/AngouriMath/issues/2074 + if (odd) + { + var scale = im.IsZero ? re : re.IsZero ? im : Magnitude((re, im)); + return Complex.CreateAtScales(valueRe.RoundToPrecision(context), valueIm.RoundToPrecision(context), scale, scale); + } return Complex.Create(valueRe.RoundToPrecision(context), valueIm.RoundToPrecision(context)); } diff --git a/Sources/Tests/UnitTests/Core/AFunctionOfASmallArgumentTest.cs b/Sources/Tests/UnitTests/Core/AFunctionOfASmallArgumentTest.cs new file mode 100644 index 000000000..b8262a0c3 --- /dev/null +++ b/Sources/Tests/UnitTests/Core/AFunctionOfASmallArgumentTest.cs @@ -0,0 +1,51 @@ +// +// 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 function's value is rounded onto zero only as the residual it leaves at its own zeros, which is relative + /// to its argument: sin(pi) leaves about 1e-100 pi, where sin(1e-60) cancels nothing. + /// Every one of these was 0 at the default hundred digits. References from mpmath at 80 digits. + /// #2074 + /// + [Trait("Area", "Core")] + public sealed class AFunctionOfASmallArgumentTest + { + [Theory] + [InlineData("sin(10^(-60))", "1E-60")] + [InlineData("sin(-10^(-60))", "-1E-60")] + [InlineData("tan(10^(-60))", "1E-60")] + [InlineData("arcsin(10^(-60))", "1E-60")] + [InlineData("arctan(10^(-60))", "1E-60")] + [InlineData("erf(10^(-60))", "1.128379167095512573896158903121545171688101258658E-60")] + [InlineData("erfi(10^(-60))", "1.128379167095512573896158903121545171688101258658E-60")] + [InlineData("Si(10^(-60))", "1E-60")] + [InlineData("Shi(10^(-60))", "1E-60")] + public void IsItsValue(string expression, string reference) + { + var expected = EDecimal.FromString(reference); + var actual = Assert.IsType(expression.EvalNumerical()).EDecimal; + var error = actual.Subtract(expected, EContext.ForPrecision(200)).Abs().Divide(expected.Abs(), EContext.ForPrecision(10)); + Assert.True(error.LessThan(EDecimal.FromString("1E-40")), $"{expression} is {actual}, where it is {expected}"); + } + + [Theory] + [InlineData("sin(pi)")] + [InlineData("sin(100*pi)")] + [InlineData("cos(pi/2)")] + [InlineData("tan(pi)")] + [InlineData("arccos(1)")] + public void TheResidualAtAZeroIsStillZero(string expression) + => Assert.Equal(Integer.Zero, expression.EvalNumerical()); + } +}