Stop ExpM1, LogP1 and the hyperbolics building million-digit sums for large arguments [patch]

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
 

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;
+
+/// <summary>
+/// Covers the functions that add or subtract one, or a reciprocal, against an intermediate whose
+/// exponent lies far from zero: <c>ExpM1</c>, <c>LogP1</c>, <c>Sinh</c>, <c>Cosh</c>, <c>Asinh</c>
+/// and <c>Acosh</c>.
+/// </summary>
+/// <remarks>
+/// Exact addition aligns both operands to the smaller exponent, so <c>e^-1e7 - 1</c> used to build a
+/// four-million-digit significand for an answer that rounds to <c>-1</c>, and <c>ExpM1(-1e10)</c>
+/// overflowed the exponent of <c>e^-1e10</c> 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.
+/// </remarks>
+[TestClass]
+public class PreciseNumberLargeArgumentTests
+{
+	/// <summary>
+	/// Longer than any of these calls takes once bounded, and far shorter than any took before.
+	/// </summary>
+	private static readonly TimeSpan Budget = TimeSpan.FromSeconds(5);
+
+	private static PreciseNumber Parse(string text) =>
+		PreciseNumber.Parse(text, CultureInfo.InvariantCulture);
+
+	/// <summary>
+	/// Runs a call on the thread pool and fails as soon as it overruns <see cref="Budget"/>, rather
+	/// than holding the test run until a call that takes minutes finishes.
+	/// </summary>
+	private static PreciseNumber WithinBudget(Func<PreciseNumber> call, string description)
+	{
+		Task<PreciseNumber> 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");
+	}
+}

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)

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);
 	}
 
 	/// <summary>
@@ -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);
 	}
 
 	/// <summary>
@@ -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);
 	}
 

PreciseNumber/PreciseNumber.cs

@@ -731,7 +731,7 @@
 	/// A fixed precision such as <c>E15</c> rounds values that need 17 digits, which turns
 	/// <see cref="double.MaxValue"/> into a number that converts back to infinity.
 	/// </remarks>
 	internal static string GetStringFormatForFloatType<TFloat>()
 		where TFloat : INumber<TFloat>
 		=> "R";
 
@@ -1189,7 +1189,7 @@
 	public static PreciseNumber MinMagnitudeNumber(PreciseNumber x, PreciseNumber y) => MinMagnitude(x, y);
 
 	/// <inheritdoc/>
 	public static PreciseNumber Parse(ReadOnlySpan<char> s, NumberStyles style, IFormatProvider? provider)
 	{
 		// Read the separator and signs from the same culture TryFormat writes them with, so that
 		// Parse(x.ToString(p), p) round-trips for every provider.
@@ -1239,14 +1239,14 @@
 					}
 
 					hasDecimal = true;
 					i += decimalSeparator.Length - 1;
 					continue;
 				}
 
 				// Group separators may only sit between integral digits.
 				if (allowThousands && !hasDecimal && digitCount > 0 && s[i..].StartsWith(groupSeparator, StringComparison.Ordinal))
 				{
 					i += groupSeparator.Length - 1;
 					continue;
 				}
 
@@ -1554,6 +1554,69 @@
 		return new PreciseNumber(commonExponent, commonLeft + commonRight);
 	}
 
+	/// <summary>
+	/// 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.
+	/// </summary>
+	private const int StickyGuardDigits = 3;
+
+	/// <summary>
+	/// Adds two numbers whose sum is only ever going to be rounded to a given number of significant
+	/// digits.
+	/// </summary>
+	/// <param name="left">The first number to add.</param>
+	/// <param name="right">The second number to add.</param>
+	/// <param name="significantDigits">The precision the caller will round the sum to.</param>
+	/// <returns>
+	/// 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
+	/// <see cref="StickyGuardDigits"/> below the last digit the caller keeps.
+	/// </returns>
+	/// <remarks>
+	/// <see cref="Add(PreciseNumber, PreciseNumber)"/> is exact, so <c>1 + 1e-4000000</c> scales the
+	/// one by <c>10^4000000</c> 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.
+	/// </remarks>
+	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));
+	}
+
+	/// <summary>
+	/// Subtracts one number from another where the difference is only ever going to be rounded to a
+	/// given number of significant digits.
+	/// </summary>
+	/// <param name="left">The number to subtract from.</param>
+	/// <param name="right">The number to subtract.</param>
+	/// <param name="significantDigits">The precision the caller will round the difference to.</param>
+	/// <returns>The difference, as <see cref="AddToPrecision(PreciseNumber, PreciseNumber, int)"/> forms it.</returns>
+	internal static PreciseNumber SubtractToPrecision(PreciseNumber left, PreciseNumber right, int significantDigits) =>
+		AddToPrecision(left, -right, significantDigits);
+
 	/// <summary>
 	/// Multiplies two numbers.
 	/// </summary>