diff --git a/CLAUDE.md b/CLAUDE.md index f83c164..5d7df3d 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -41,6 +41,7 @@ dotnet run -c Release --project PreciseNumber.Benchmarks -- --filter '*' --job s - Exponentials, logarithms and powers (`PreciseNumber/PreciseNumber.Exponentials.cs`, satisfying `IExponentialFunctions`, `ILogarithmicFunctions` and `IPowerFunctions`) follow the same precision rule and do not route through `double` either. `ln(m · 10^k)` is `ln m + k · ln 10` against the stored `Ln10`, with the mantissa centred on `[1/√10, √10)` and fed to the atanh series; `exp(v)` factors out `10^round(v / ln 10)` as an exponent shift and halves what is left before a Taylor sum. Nothing here is a free-standing decision: `Exp10`/`Log10` must not route through the natural log, because the exponent is the whole answer for a power of ten, and the `…M1`/`…P1` variants must not be computed as `Exp(x) - 1`/`Log(1 + x)`, because that cancels away the precision near zero they exist to keep - A fractional `Pow` is `exp(y · ln x)`, carried wider by the integer digits of `y · ln x` because `Exp`'s range reduction consumes them. The integer path stays exponentiation by squaring and is exact; tests pin that exactness rather than a tolerance - Hyperbolics (`PreciseNumber/PreciseNumber.Hyperbolics.cs`, satisfying `IHyperbolicFunctions`) are the exponentials and logarithms under other names and add no transcendental machinery of their own. Which of them needs its textbook form rearranged is narrower than floating-point habit suggests, and the reason is worth keeping straight: addition, subtraction and multiplication here are *exact*, so cancelling two nearly equal values costs nothing by itself. Digits are lost only by cancelling against something `Exp`, `Log`, `Sqrt` or `Divide` has **already** rounded to a working width. That is why `sinh` sums its own series below one half instead of taking `(e^x - e^-x)/2`, and why `asinh` and `atanh` subtract their one analytically and go through `LogP1` — each of those three otherwise returns about thirty correct digits from a fifty-digit type at `1e-30`, and the tests fail outright on them rather than drifting in the last place. `cosh`, `tanh` and `acosh` do **not** need it: `cosh` sums two positive terms, and the other two cancel only against operands nothing has rounded yet. `tanh` is still written as `-t/(2 + t)` with `t = expm1(-2x)`, but for range rather than precision — the negative exponent decays instead of growing, so it saturates to `±1` where `e^2x` would overflow — and `acosh` still factors the difference of squares, to keep a `2n`-digit intermediate out of the root. Don't simplify the first three back, and don't defend the last two on precision grounds +- `Add` and `Subtract` are exact, which is a trap when one operand sits millions of digits below the other: aligning them builds a significand that wide only for the caller to round it away. Where a sum is only ever going to be rounded (the `±1` in `ExpM1` and `LogP1`, the reciprocal in `Sinh` and `Cosh`, the ones in `Asinh` and `Acosh`), use `AddToPrecision`/`SubtractToPrecision`, which collapse an operand wholly below the kept digits into one sticky unit that rounds the same way. `ExpM1` also returns `-1` outright once `e^x` is below the working width, since `Exp` of a very negative argument overflows its exponent - The `sanitize` constructor parameter controls whether trailing zeros are removed (default: true) - Constants (`Zero`, `One`, `Pi`, `E`, `Tau`) are pre-computed static instances - `ConstantPrecision` is a ceiling, not just a width. `PiTo`/`ETo`/`Ln2To`/`Ln10To` serve at most that many digits and return the capped constant rather than failing, which is deliberate and pinned by `TestConstantAccessorsReturnTheWholeConstantWhenAskedForMore` — so a caller that *needs* what it asked for has to check, because the accessor will not. The circular functions are where this bites: reducing modulo `π/2` spends one digit of π per integer digit of the argument, so only about `ConstantPrecision - integerDigits` are left for the answer. `RequireReducibleArgument` enforces exactly that sum and nothing wider. Do not tighten it to `reductionDigits`: that pads the requirement with `TrigonometricGuardDigits` twice over as margin, margin is allowed to be unavailable, and bounding it would refuse `Sin(1000000, 130)` — a call whose every reported digit is correct. The half-turn family (`SinPi` and the rest) has no such ceiling, because it reduces before it multiplies, and `Log` has none either, because it carries the magnitude in the decimal exponent rather than cancelling it against a constant @@ -49,7 +50,7 @@ dotnet run -c Release --project PreciseNumber.Benchmarks -- --filter '*' --job s ### Test Structure -Tests use MSTest. `PreciseNumber.Test/PreciseNumberTests.cs` covers arithmetic, parsing, and formatting, `PreciseNumberConversionTests.cs` covers generic math conversion in every mode, `PreciseNumberRootTests.cs` pins the roots against published digits and against squaring back, `PreciseNumberExponentialTests.cs` does the same for the exponentials and logarithms and additionally pins the cases a `double` fallback cannot reach — fifty published digits of a fractional power, and `ExpM1`/`LogP1` of `1e-30` not collapsing to zero — `PreciseNumberHyperbolicTests.cs` pins the hyperbolics against published digits and against `cosh²x - sinh²x = 1`, and separates the small-argument assertions that actually discriminate (`sinh`, `asinh`, `atanh`) from the two that read like they do and don't (`tanh`, `acosh`) — the class remark records which is which, so the distinction survives the next person to read it — and `PreciseNumberValueTypeTests.cs` pins `default` as zero and asserts that small-value addition, subtraction, multiplication, and comparison allocate nothing. The test project targets only .NET 10.0 while the main library multi-targets net7.0, net8.0, net9.0, and net10.0. +Tests use MSTest. `PreciseNumber.Test/PreciseNumberTests.cs` covers arithmetic, parsing, and formatting, `PreciseNumberConversionTests.cs` covers generic math conversion in every mode, `PreciseNumberRootTests.cs` pins the roots against published digits and against squaring back, `PreciseNumberExponentialTests.cs` does the same for the exponentials and logarithms and additionally pins the cases a `double` fallback cannot reach — fifty published digits of a fractional power, and `ExpM1`/`LogP1` of `1e-30` not collapsing to zero — `PreciseNumberHyperbolicTests.cs` pins the hyperbolics against published digits and against `cosh²x - sinh²x = 1`, and separates the small-argument assertions that actually discriminate (`sinh`, `asinh`, `atanh`) from the two that read like they do and don't (`tanh`, `acosh`) — the class remark records which is which, so the distinction survives the next person to read it — `PreciseNumberLargeArgumentTests.cs` bounds the time of `ExpM1`/`LogP1`/hyperbolics on huge arguments and pins digits the exact sums produced, and `PreciseNumberValueTypeTests.cs` pins `default` as zero and asserts that small-value addition, subtraction, multiplication, and comparison allocate nothing. The test project targets only .NET 10.0 while the main library multi-targets net7.0, net8.0, net9.0, and net10.0. ### Benchmarks diff --git a/PreciseNumber.Test/PreciseNumberLargeArgumentTests.cs b/PreciseNumber.Test/PreciseNumberLargeArgumentTests.cs new file mode 100644 index 0000000..12f221a --- /dev/null +++ b/PreciseNumber.Test/PreciseNumberLargeArgumentTests.cs @@ -0,0 +1,105 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Test; + +using System; +using System.Globalization; +using System.Threading.Tasks; + +/// +/// Covers the functions that add or subtract one, or a reciprocal, against an intermediate whose +/// exponent lies far from zero: ExpM1, LogP1, Sinh, Cosh, Asinh +/// and Acosh. +/// +/// +/// Exact addition aligns both operands to the smaller exponent, so e^-1e7 - 1 used to build a +/// four-million-digit significand for an answer that rounds to -1, and ExpM1(-1e10) +/// overflowed the exponent of e^-1e10 before it got that far. The pinned digits below were +/// produced by the exact path before it was bounded, so they show the bounded path rounds the same. +/// +[TestClass] +public class PreciseNumberLargeArgumentTests +{ + /// + /// Longer than any of these calls takes once bounded, and far shorter than any took before. + /// + private static readonly TimeSpan Budget = TimeSpan.FromSeconds(5); + + private static PreciseNumber Parse(string text) => + PreciseNumber.Parse(text, CultureInfo.InvariantCulture); + + /// + /// Runs a call on the thread pool and fails as soon as it overruns , rather + /// than holding the test run until a call that takes minutes finishes. + /// + private static PreciseNumber WithinBudget(Func call, string description) + { + Task task = Task.Run(call); + Assert.IsTrue(task.Wait(Budget), $"{description} took longer than {Budget.TotalSeconds} seconds"); + return task.Result; + } + + [TestMethod] + public void TestExpM1OfAVeryNegativeArgumentIsMinusOne() + { + Assert.AreEqual(PreciseNumber.NegativeOne, PreciseNumber.ExpM1(Parse("-1E10")), "ExpM1(-1e10) is not -1"); + Assert.AreEqual(PreciseNumber.NegativeOne, PreciseNumber.Exp2M1(Parse("-1E10")), "Exp2M1(-1e10) is not -1"); + Assert.AreEqual(PreciseNumber.NegativeOne, PreciseNumber.Exp10M1(Parse("-1E10")), "Exp10M1(-1e10) is not -1"); + Assert.AreEqual(PreciseNumber.NegativeOne, PreciseNumber.ExpM1(Parse("-1E10"), 200), "ExpM1(-1e10) at 200 digits is not -1"); + } + + [TestMethod] + public void TestExpM1OfALargeArgumentFinishesQuickly() + { + Assert.AreEqual(PreciseNumber.NegativeOne, WithinBudget(() => PreciseNumber.ExpM1(Parse("-1E7")), "ExpM1(-1e7)"), "ExpM1(-1e7) is not -1"); + Assert.AreEqual(PreciseNumber.NegativeOne, WithinBudget(() => PreciseNumber.Exp10M1(Parse("-1E7")), "Exp10M1(-1e7)"), "Exp10M1(-1e7) is not -1"); + + PreciseNumber expected = PreciseNumber.Exp(Parse("1E6")); + PreciseNumber actual = WithinBudget(() => PreciseNumber.ExpM1(Parse("1E6")), "ExpM1(1e6)"); + Assert.AreEqual(expected, actual, "ExpM1(1e6) differs from Exp(1e6), though the one is far below its last digit"); + } + + [TestMethod] + public void TestLogP1OfAHugeArgumentFinishesQuickly() + { + PreciseNumber x = Parse("1E1000000"); + Assert.AreEqual(PreciseNumber.Log(x), WithinBudget(() => PreciseNumber.LogP1(x), "LogP1(1e1000000)"), "LogP1(1e1000000) differs from Log(1e1000000)"); + } + + [TestMethod] + public void TestSinhAndCoshOfALargeArgumentFinishQuickly() + { + PreciseNumber sinh = WithinBudget(() => PreciseNumber.Sinh(Parse("1E7")), "Sinh(1e7)"); + PreciseNumber cosh = WithinBudget(() => PreciseNumber.Cosh(Parse("-1E7")), "Cosh(-1e7)"); + + // At this size e^-x is millions of digits below e^x, so both are e^x / 2 to every digit kept. + Assert.AreEqual(sinh, cosh, "Sinh(1e7) and Cosh(-1e7) disagree"); + Assert.AreEqual(-sinh, WithinBudget(() => PreciseNumber.Sinh(Parse("-1E7")), "Sinh(-1e7)"), "Sinh is not odd at 1e7"); + } + + [TestMethod] + public void TestAsinhAndAcoshOfAHugeArgumentFinishQuickly() + { + PreciseNumber x = Parse("1E1000000"); + + // Both are ln(2x) once 1/x is far below the last digit kept. + PreciseNumber expected = PreciseNumber.Log(Parse("2E1000000")); + Assert.AreEqual(expected, WithinBudget(() => PreciseNumber.Asinh(x), "Asinh(1e1000000)"), "Asinh(1e1000000) is not ln(2e1000000)"); + Assert.AreEqual(-expected, WithinBudget(() => PreciseNumber.Asinh(-x), "Asinh(-1e1000000)"), "Asinh(-1e1000000) is not -ln(2e1000000)"); + Assert.AreEqual(expected, WithinBudget(() => PreciseNumber.Acosh(x), "Acosh(1e1000000)"), "Acosh(1e1000000) is not ln(2e1000000)"); + } + + [TestMethod] + public void TestBoundedSumsKeepTheDigitsTheExactSumsGave() + { + Assert.AreEqual(Parse("-0.99999999999999999999999999999999999999999996279924"), PreciseNumber.ExpM1(Parse("-100"), 50), "ExpM1(-100) changed"); + Assert.AreEqual(Parse("1.9424263952412559365842088360176992193662086219516E130"), PreciseNumber.ExpM1(Parse("300"), 50), "ExpM1(300) changed"); + Assert.AreEqual(Parse("3.1622776601683793319988935444327185337195551393252E100"), PreciseNumber.Exp10M1(Parse("100.5"), 50), "Exp10M1(100.5) changed"); + Assert.AreEqual(Parse("692.08386071786388396574754062165634523662942268717"), PreciseNumber.LogP1(Parse("3.7E300"), 50), "LogP1(3.7e300) changed"); + Assert.AreEqual(Parse("9.712131976206279682921044180088496096831043109758E129"), PreciseNumber.Sinh(Parse("300"), 50), "Sinh(300) changed"); + Assert.AreEqual(Parse("1.60125985730181641832516084312451811087207753341E130"), PreciseNumber.Cosh(Parse("-300.5"), 50), "Cosh(-300.5) changed"); + Assert.AreEqual(Parse("-692.3849658106478055799981957385354499198560479429"), PreciseNumber.Asinh(Parse("-2.5E300"), 50), "Asinh(-2.5e300) changed"); + Assert.AreEqual(Parse("693.44967654764023392316347631621299731774945936819"), PreciseNumber.Acosh(Parse("7.25E300"), 50), "Acosh(7.25e300) changed"); + Assert.AreEqual(Parse("92.796550900321772670136890308832744872119559679511"), PreciseNumber.Acosh(Parse("1E40"), 50), "Acosh(1e40) changed"); + } +} diff --git a/PreciseNumber/PreciseNumber.Exponentials.cs b/PreciseNumber/PreciseNumber.Exponentials.cs index b4a3d5f..ece0703 100644 --- a/PreciseNumber/PreciseNumber.Exponentials.cs +++ b/PreciseNumber/PreciseNumber.Exponentials.cs @@ -317,7 +317,7 @@ public static PreciseNumber LogP1(PreciseNumber x, int significantDigits) // better conditioned than an atanh argument approaching one. if (Abs(x) > DirectSeriesLimit) { - return Log(Add(One, x), significantDigits); + return Log(AddToPrecision(One, x, significantDigits + ExponentialGuardDigits), significantDigits); } int working = significantDigits + ExponentialGuardDigits; @@ -486,7 +486,16 @@ public static PreciseNumber ExpM1(PreciseNumber x, int significantDigits) if (Abs(x) > DirectSeriesLimit) { int wide = significantDigits + ExponentialGuardDigits; - return Subtract(Exp(x, wide), One).ReduceSignificance(significantDigits); + + // Below -3 · (wide + 1), e^x is under 10^-(wide + 1) because ln 10 is under three, so it + // cannot reach the digits kept and the answer is -1. Exp itself would build that tiny + // value first, or overflow its exponent trying to. + if (x < new PreciseNumber(0, -3 * ((BigInteger)wide + 1))) + { + return NegativeOne; + } + + return SubtractToPrecision(Exp(x, wide), One, wide).ReduceSignificance(significantDigits); } return ExpSeriesWithoutLeadingOne(x, significantDigits + ExponentialGuardDigits) diff --git a/PreciseNumber/PreciseNumber.Hyperbolics.cs b/PreciseNumber/PreciseNumber.Hyperbolics.cs index b544cce..c60c5b5 100644 --- a/PreciseNumber/PreciseNumber.Hyperbolics.cs +++ b/PreciseNumber/PreciseNumber.Hyperbolics.cs @@ -107,7 +107,7 @@ public static PreciseNumber Sinh(PreciseNumber x, int significantDigits) PreciseNumber raised = Exp(x, working); PreciseNumber lowered = Divide(One, raised, working); - return Divide(Subtract(raised, lowered), Two, working).ReduceSignificance(significantDigits); + return Divide(SubtractToPrecision(raised, lowered, working), Two, working).ReduceSignificance(significantDigits); } /// @@ -150,7 +150,7 @@ public static PreciseNumber Cosh(PreciseNumber x, int significantDigits) int working = significantDigits + ExponentialGuardDigits; PreciseNumber raised = Exp(x, working); PreciseNumber lowered = Divide(One, raised, working); - return Divide(Add(raised, lowered), Two, working).ReduceSignificance(significantDigits); + return Divide(AddToPrecision(raised, lowered, working), Two, working).ReduceSignificance(significantDigits); } /// @@ -256,8 +256,8 @@ public static PreciseNumber Asinh(PreciseNumber x, int significantDigits) int working = significantDigits + ExponentialGuardDigits; PreciseNumber magnitude = Abs(x); PreciseNumber square = Multiply(magnitude, magnitude); - PreciseNumber root = Sqrt(Add(One, square), working); - PreciseNumber excess = Divide(square, Add(One, root), working); + PreciseNumber root = Sqrt(AddToPrecision(One, square, working), working); + PreciseNumber excess = Divide(square, AddToPrecision(One, root, working), working); PreciseNumber result = LogP1(Add(magnitude, excess), significantDigits); return x.Significand.Sign > 0 ? result : -result; @@ -320,8 +320,8 @@ public static PreciseNumber Acosh(PreciseNumber x, int significantDigits) } int working = significantDigits + ExponentialGuardDigits; - PreciseNumber below = Subtract(x, One); - PreciseNumber root = Sqrt(Multiply(below, Add(x, One)), working); + PreciseNumber below = SubtractToPrecision(x, One, working); + PreciseNumber root = Sqrt(Multiply(below, AddToPrecision(x, One, working)), working); return LogP1(Add(below, root), significantDigits); } diff --git a/PreciseNumber/PreciseNumber.cs b/PreciseNumber/PreciseNumber.cs index 621bd81..066ea33 100644 --- a/PreciseNumber/PreciseNumber.cs +++ b/PreciseNumber/PreciseNumber.cs @@ -1554,6 +1554,69 @@ public static PreciseNumber Add(PreciseNumber left, PreciseNumber right) return new PreciseNumber(commonExponent, commonLeft + commonRight); } + /// + /// Digits kept below the last one a precision-bounded sum is wanted to, before a far smaller + /// operand is collapsed into a single sticky digit. + /// + private const int StickyGuardDigits = 3; + + /// + /// Adds two numbers whose sum is only ever going to be rounded to a given number of significant + /// digits. + /// + /// The first number to add. + /// The second number to add. + /// The precision the caller will round the sum to. + /// + /// The exact sum when the operands are close in magnitude. Otherwise the larger operand plus a + /// single unit of the smaller one's sign, placed below every digit of the larger operand and + /// below the last digit the caller keeps. + /// + /// + /// is exact, so 1 + 1e-4000000 scales the + /// one by 10^4000000 and builds a four-million-digit significand only for the caller to + /// round almost all of it away. When the smaller operand lies wholly below the larger one's last + /// digit and below the caller's rounding position, the exact sum and the sticky sum lie strictly + /// inside the same gap between neighbouring multiples of that digit. Every rounding boundary at + /// the caller's precision, or one digit coarser after a halving, is such a multiple, so both + /// round the same way. Only the sign of the small operand can still matter, and the sticky unit + /// keeps it. + /// + internal static PreciseNumber AddToPrecision(PreciseNumber left, PreciseNumber right, int significantDigits) + { + if (left.Significand.IsZero || right.Significand.IsZero) + { + return Add(left, right); + } + + long leftMagnitude = (long)left.Exponent + left.SignificantDigits; + long rightMagnitude = (long)right.Exponent + right.SignificantDigits; + (PreciseNumber larger, PreciseNumber smaller, long largerMagnitude, long smallerMagnitude) = leftMagnitude >= rightMagnitude + ? (left, right, leftMagnitude, rightMagnitude) + : (right, left, rightMagnitude, leftMagnitude); + + // The smaller operand is below 10^smallerMagnitude, so it is wholly beneath 10^floor once + // smallerMagnitude is at most floor. + long floor = Math.Min(larger.Exponent, largerMagnitude - significantDigits - StickyGuardDigits); + if (smallerMagnitude > floor || floor - 1 < int.MinValue) + { + return Add(left, right); + } + + return Add(larger, new PreciseNumber((int)(floor - 1), smaller.Significand.Sign)); + } + + /// + /// Subtracts one number from another where the difference is only ever going to be rounded to a + /// given number of significant digits. + /// + /// The number to subtract from. + /// The number to subtract. + /// The precision the caller will round the difference to. + /// The difference, as forms it. + internal static PreciseNumber SubtractToPrecision(PreciseNumber left, PreciseNumber right, int significantDigits) => + AddToPrecision(left, -right, significantDigits); + /// /// Multiplies two numbers. ///