diff --git a/CLAUDE.md b/CLAUDE.md index 6ab5460..58d57f2 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -40,6 +40,7 @@ dotnet run -c Release --project PreciseNumber.Benchmarks -- --filter '*' --job s - Roots (`PreciseNumber/PreciseNumber.Roots.cs`, satisfying `IRootFunctions`) follow `Divide`'s precision rule and do not route through `double`. Each scales the significand by a power of ten until the degree divides the exponent, then takes an integer Newton root of the significand, so an exact root stops on the exact answer rather than on a tolerance and no seed has to survive a value outside `double`'s range - 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 - The `sanitize` constructor parameter controls whether trailing zeros are removed (default: true) - Constants (`Zero`, `One`, `Pi`, `E`, `Tau`) are pre-computed static instances - As a value type it can't be null or inherited. Don't add null checks for `PreciseNumber` parameters, and don't reintroduce `protected` members @@ -47,7 +48,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 — 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 — 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.Benchmarks/HyperbolicBenchmarks.cs b/PreciseNumber.Benchmarks/HyperbolicBenchmarks.cs new file mode 100644 index 0000000..6740bb5 --- /dev/null +++ b/PreciseNumber.Benchmarks/HyperbolicBenchmarks.cs @@ -0,0 +1,105 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Benchmarks; + +using BenchmarkDotNet.Attributes; + +/// +/// Measures the hyperbolic functions and their inverses. +/// +/// +/// None of these carries a series of its own except Sinh below one half, so read them against +/// rather than against each other in isolation: each should cost +/// about what the exponential or logarithm underneath it costs, plus the arithmetic named below. A +/// row that is several times its underlying transcendental is doing work it should not be. +/// +/// against is the pair worth watching. +/// The small case sums its own series and the larger one takes an exponential and a reciprocal, so +/// the two are different algorithms rather than the same one at different arguments. The series +/// should win at every point on the Digits axis — it exists for precision rather than speed, +/// but it would be worth knowing if it cost more than the path it replaces. +/// +/// +/// is one exponential, one reciprocal and one halving, which makes it the +/// floor for this file. is the opposite end: it returns +/// ±1 without taking an exponential at all, so it should be near free and should not move +/// with the Digits axis. Anything else there means the saturation test is not firing. +/// +/// +[MemoryDiagnoser] +public class HyperbolicBenchmarks +{ + private PreciseNumber value = PreciseNumber.Zero; + private PreciseNumber smallValue = PreciseNumber.Zero; + private PreciseNumber unitInterval = PreciseNumber.Zero; + private PreciseNumber aboveOne = PreciseNumber.Zero; + private PreciseNumber saturated = PreciseNumber.Zero; + + /// + /// Gets or sets the number of significant digits in the operand, and so in the answer. + /// + [Params(8, 30, 200)] + public int Digits { get; set; } + + /// + /// Prepares the operands. + /// + [GlobalSetup] + public void Setup() + { + // A little over one, so Sinh takes the exponential path rather than the series. + value = Operands.Number(Digits, -(Digits - 1)); + + // Well under one half, so Sinh sums its series instead. + smallValue = Operands.Number(Digits, -(Digits + 1), offset: 13); + + // In (0, 1), the domain of Atanh. + unitInterval = Operands.Number(Digits, -Digits, offset: 29); + + // Above one, the domain of Acosh. + aboveOne = PreciseNumber.One + Operands.Number(Digits, -(Digits - 1), offset: 5); + + // Far enough out that Tanh returns without taking an exponential. + saturated = Operands.Number(Digits, -(Digits - 4), offset: 17); + } + + /// Takes the hyperbolic sine, by way of an exponential and a reciprocal. + /// The result. + [Benchmark(Baseline = true)] + public PreciseNumber SinhOfAValue() => PreciseNumber.Sinh(value); + + /// Takes the hyperbolic sine of a small value, which sums its own series. + /// The result. + [Benchmark] + public PreciseNumber SinhOfASmallValue() => PreciseNumber.Sinh(smallValue); + + /// Takes the hyperbolic cosine, which is the cheapest thing here. + /// The result. + [Benchmark] + public PreciseNumber CoshOfAValue() => PreciseNumber.Cosh(value); + + /// Takes the hyperbolic tangent, which is one exponential and a division. + /// The result. + [Benchmark] + public PreciseNumber TanhOfAValue() => PreciseNumber.Tanh(value); + + /// Takes the hyperbolic tangent of a value past saturation, which takes no exponential. + /// The result. + [Benchmark] + public PreciseNumber TanhOfASaturatedValue() => PreciseNumber.Tanh(saturated); + + /// Takes the inverse hyperbolic sine, which is a root and a logarithm. + /// The result. + [Benchmark] + public PreciseNumber AsinhOfAValue() => PreciseNumber.Asinh(value); + + /// Takes the inverse hyperbolic cosine, which is a root and a logarithm. + /// The result. + [Benchmark] + public PreciseNumber AcoshOfAValue() => PreciseNumber.Acosh(aboveOne); + + /// Takes the inverse hyperbolic tangent, which is a division and a logarithm. + /// The result. + [Benchmark] + public PreciseNumber AtanhOfAValue() => PreciseNumber.Atanh(unitInterval); +} diff --git a/PreciseNumber.Test/PreciseNumberHyperbolicTests.cs b/PreciseNumber.Test/PreciseNumberHyperbolicTests.cs new file mode 100644 index 0000000..87421aa --- /dev/null +++ b/PreciseNumber.Test/PreciseNumberHyperbolicTests.cs @@ -0,0 +1,431 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Test; + +using System.Globalization; +using System.Linq; +using System.Numerics; + +/// +/// Covers on . +/// +/// +/// The digit-for-digit assertions carry values computed independently of this library, at eighty +/// digits, from the constants each function is defined against — so a change that makes the +/// implementation agree with itself but not with mathematics still fails. +/// +/// Three of the small-argument assertions are load-bearing, and it is worth recording which, since +/// the set is narrower than the floating-point intuition suggests. Substituting the textbook form +/// back in was tried, one function at a time, and only sinh's (e^x - e^-x)/2, +/// asinh's ln(x + √(x² + 1)) and atanh's ½ ln((1 + x)/(1 - x)) actually +/// lose digits: each cancels against a value that Exp, Sqrt or Divide had +/// already rounded to the working width, and a 1e-30 argument comes back with about thirty +/// of the fifty digits asked for. , +/// and +/// fail outright on those forms rather than drifting in the last place, which is what stops a later +/// simplification quietly undoing the rearrangements. +/// +/// +/// and +/// read like members of that set and are not. +/// Their textbook forms pass, because this type's addition, subtraction and multiplication are +/// exact and nothing has been rounded before the cancellation. They are kept as ordinary regression +/// assertions — the values they pin are still the right ones — but they should not be relied on to +/// catch a rewrite of those two functions. +/// +/// +[TestClass] +public class PreciseNumberHyperbolicTests +{ + /// + /// The first fifty significant digits of the hyperbolic functions of one. + /// + private const string Sinh1Digits = "11752011936438014568823818505956008151557179813341"; + private const string Cosh1Digits = "15430806348152437784779056207570616826015291123659"; + private const string Tanh1Digits = "76159415595576488811945828260479359041276859725794"; + + /// + /// The first fifty significant digits of the inverse hyperbolic functions, each of which is a + /// logarithm in closed form: asinh 1 = ln(1 + √2), acosh 2 = ln(2 + √3) and + /// atanh ½ = ½ ln 3. + /// + private const string Asinh1Digits = "88137358701954302523260932497979230902816032826164"; + private const string Acosh2Digits = "13169578969248167086250463473079684440269819714675"; + private const string AtanhHalfDigits = "54930614433405484569762261846126285232374527891137"; + + /// + /// The first fifty significant digits of acosh(1 + 1e-30), which is 1.414…e-15. + /// + /// + /// Close to √(2 · 1e-30) but not equal to it: acosh(1 + d) = √(2d) · (1 - d/12 + …), + /// so this and √2 share their first twenty-nine digits and then diverge. Comparing against + /// the true value rather than against √2 is what makes the assertion mean something. + /// + private const string AcoshJustAboveOneDigits = "14142135623730950488016887242095802274394741174562"; + + private static PreciseNumber Parse(string text) => + PreciseNumber.Parse(text, CultureInfo.InvariantCulture); + + private static string Digits(PreciseNumber value) => + value.Significand.ToString(CultureInfo.InvariantCulture); + + /// + /// Asserts that two values agree to a number of significant digits, comparing relative to the + /// expected magnitude so the assertion means the same thing at every exponent. + /// + private static void AssertAgreesTo(PreciseNumber expected, PreciseNumber actual, int digits, string message) + { + PreciseNumber difference = PreciseNumber.Abs(actual - expected); + PreciseNumber tolerance = PreciseNumber.Abs(expected) * Parse($"1E-{digits.ToString(CultureInfo.InvariantCulture)}"); + + Assert.IsLessThanOrEqualTo( + tolerance, + difference, + $"{message}: expected {expected}, got {actual}, which differs by {difference}"); + } + + /// + /// A spread that straddles the magnitude at which Sinh stops summing its own series and + /// defers to Exp, and includes both signs of each. + /// + private static PreciseNumber[] Sweep() => + [ + Parse("1E-30"), + Parse("-1E-30"), + Parse("0.25"), + Parse("-0.25"), + Parse("0.5"), + Parse("1"), + Parse("-1"), + Parse("2.75"), + Parse("-2.75"), + ]; + + [TestMethod] + public void TestSinhMatchesPublishedDigits() => + Assert.AreEqual(Sinh1Digits, Digits(PreciseNumber.Sinh(PreciseNumber.One, 50)), "Sinh(1) is wrong"); + + [TestMethod] + public void TestCoshMatchesPublishedDigits() => + Assert.AreEqual(Cosh1Digits, Digits(PreciseNumber.Cosh(PreciseNumber.One, 50)), "Cosh(1) is wrong"); + + [TestMethod] + public void TestTanhMatchesPublishedDigits() => + Assert.AreEqual(Tanh1Digits, Digits(PreciseNumber.Tanh(PreciseNumber.One, 50)), "Tanh(1) is wrong"); + + [TestMethod] + public void TestAsinhMatchesPublishedDigits() => + Assert.AreEqual(Asinh1Digits, Digits(PreciseNumber.Asinh(PreciseNumber.One, 50)), "Asinh(1) is wrong"); + + [TestMethod] + public void TestAcoshMatchesPublishedDigits() => + Assert.AreEqual(Acosh2Digits, Digits(PreciseNumber.Acosh(2.ToPreciseNumber(), 50)), "Acosh(2) is wrong"); + + [TestMethod] + public void TestAtanhMatchesPublishedDigits() => + Assert.AreEqual(AtanhHalfDigits, Digits(PreciseNumber.Atanh(Parse("0.5"), 50)), "Atanh(0.5) is wrong"); + + /// + /// Pins the series path against the identity that defines it. + /// + /// + /// sinh below one half is summed rather than taken from two exponentials, so this is the + /// assertion that the two paths agree on the answer and not merely on the neighbourhood. + /// + [TestMethod] + public void TestCoshSquaredLessSinhSquaredIsOneAcrossASweep() + { + foreach (PreciseNumber x in Sweep()) + { + PreciseNumber sinh = PreciseNumber.Sinh(x, 50); + PreciseNumber cosh = PreciseNumber.Cosh(x, 50); + AssertAgreesTo( + PreciseNumber.One, + (cosh * cosh) - (sinh * sinh), + 45, + $"cosh²x - sinh²x is not one at x = {x}"); + } + } + + [TestMethod] + public void TestTanhIsSinhOverCoshAcrossASweep() + { + foreach (PreciseNumber x in Sweep()) + { + PreciseNumber expected = PreciseNumber.Divide( + PreciseNumber.Sinh(x, 60), + PreciseNumber.Cosh(x, 60), + 50); + + AssertAgreesTo(expected, PreciseNumber.Tanh(x, 50), 45, $"Tanh disagrees with Sinh/Cosh at x = {x}"); + } + } + + [TestMethod] + public void TestSinhAndAsinhRoundTripAcrossASweep() + { + foreach (PreciseNumber x in Sweep()) + { + AssertAgreesTo(x, PreciseNumber.Asinh(PreciseNumber.Sinh(x, 60), 50), 45, $"Asinh(Sinh(x)) lost x = {x}"); + } + } + + [TestMethod] + public void TestTanhAndAtanhRoundTripAcrossASweep() + { + foreach (PreciseNumber x in Sweep()) + { + AssertAgreesTo(x, PreciseNumber.Atanh(PreciseNumber.Tanh(x, 60), 50), 45, $"Atanh(Tanh(x)) lost x = {x}"); + } + } + + /// + /// acosh is the inverse of cosh only on the non-negative half, since cosh is + /// even, so the round trip returns the magnitude rather than the value. + /// + /// + /// The sweep's smallest values are left out, and the reason is the function rather than the + /// implementation: cosh is flat at zero, so cosh(1e-30) - 1 is 5e-61 and + /// falls below any working precision short of sixty-one digits. Once cosh has rounded to + /// exactly one there is no acosh that can recover the argument, and asserting otherwise + /// would be asserting against arithmetic rather than against this code. + /// + [TestMethod] + public void TestCoshAndAcoshRoundTripToTheMagnitudeAcrossASweep() + { + foreach (PreciseNumber x in Sweep().Where(x => PreciseNumber.Abs(x) >= Parse("0.25"))) + { + PreciseNumber recovered = PreciseNumber.Acosh(PreciseNumber.Cosh(x, 60), 50); + AssertAgreesTo(PreciseNumber.Abs(x), recovered, 40, $"Acosh(Cosh(x)) lost |x| at x = {x}"); + } + } + + /// + /// The assertion that fails on (e^x - e^-x) / 2. + /// + /// + /// Both exponentials are one to thirty digits at this argument, so the textbook difference keeps + /// only the digits that survive it — about thirty of the fifty asked for. The series keeps all + /// fifty, and forty-five is comfortably on the far side of the gap between them. + /// + [TestMethod] + public void TestSinhKeepsItsDigitsNearZero() + { + PreciseNumber tiny = Parse("1E-30"); + AssertAgreesTo(tiny, PreciseNumber.Sinh(tiny, 50), 45, "Sinh(1e-30) has lost its precision"); + Assert.AreNotEqual(PreciseNumber.Zero, PreciseNumber.Sinh(tiny, 50), "Sinh(1e-30) collapsed to zero"); + } + + /// + /// The assertion that fails on ln(x + √(x² + 1)). + /// + [TestMethod] + public void TestAsinhKeepsItsDigitsNearZero() + { + PreciseNumber tiny = Parse("1E-30"); + AssertAgreesTo(tiny, PreciseNumber.Asinh(tiny, 50), 45, "Asinh(1e-30) has lost its precision"); + } + + /// + /// The assertion that fails on ½ ln((1 + x)/(1 - x)). + /// + [TestMethod] + public void TestAtanhKeepsItsDigitsNearZero() + { + PreciseNumber tiny = Parse("1E-30"); + AssertAgreesTo(tiny, PreciseNumber.Atanh(tiny, 50), 45, "Atanh(1e-30) has lost its precision"); + } + + /// + /// Pins tanh near zero. Not a discriminating assertion — see the note on this class. + /// + /// + /// A tanh built as (1 - e^-2x) / (1 + e^-2x) against a literal one passes this too, + /// because the subtraction is exact and the exponential has not been rounded against anything + /// first. The value is still worth pinning; it just does not police the implementation the way + /// does. + /// + [TestMethod] + public void TestTanhKeepsItsDigitsNearZero() + { + PreciseNumber tiny = Parse("1E-30"); + AssertAgreesTo(tiny, PreciseNumber.Tanh(tiny, 50), 45, "Tanh(1e-30) has lost its precision"); + } + + /// + /// Pins acosh just above one, where the answer is most easily got wrong. + /// + /// + /// The unfactored ln(x + √(x² - 1)) passes this as well — squaring and subtracting are + /// both exact here, so the cancellation costs nothing but width. What the assertion does catch is + /// an answer that is merely plausible: the expected value is not √(2e-30), since + /// acosh(1 + d) = √(2d) · (1 - d/12 + …), so the two share their first twenty-nine digits + /// and then part. A rearrangement that kept the precision but dropped the correction term would + /// look right to any shorter comparison and fails this one. + /// + [TestMethod] + public void TestAcoshKeepsItsDigitsJustAboveOne() + { + PreciseNumber justAbove = Parse("1.000000000000000000000000000001"); + Assert.AreEqual( + AcoshJustAboveOneDigits, + Digits(PreciseNumber.Acosh(justAbove, 50)), + "Acosh just above one is wrong"); + } + + [TestMethod] + public void TestHyperbolicFunctionsOfZeroAreExact() + { + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Sinh(PreciseNumber.Zero)); + Assert.AreEqual(PreciseNumber.One, PreciseNumber.Cosh(PreciseNumber.Zero)); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Tanh(PreciseNumber.Zero)); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Asinh(PreciseNumber.Zero)); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Acosh(PreciseNumber.One)); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Atanh(PreciseNumber.Zero)); + } + + [TestMethod] + public void TestSinhTanhAsinhAndAtanhAreOddAndCoshIsEven() + { + foreach (PreciseNumber x in Sweep()) + { + Assert.AreEqual(-PreciseNumber.Sinh(x, 50), PreciseNumber.Sinh(-x, 50), $"Sinh is not odd at x = {x}"); + Assert.AreEqual(-PreciseNumber.Tanh(x, 50), PreciseNumber.Tanh(-x, 50), $"Tanh is not odd at x = {x}"); + Assert.AreEqual(-PreciseNumber.Asinh(x, 50), PreciseNumber.Asinh(-x, 50), $"Asinh is not odd at x = {x}"); + Assert.AreEqual(PreciseNumber.Cosh(x, 50), PreciseNumber.Cosh(-x, 50), $"Cosh is not even at x = {x}"); + } + } + + /// + /// Pins the saturation branch, and with it the claim that Tanh cannot overflow. + /// + /// + /// Tanh takes the exponential of minus twice its argument, so a large argument sends that + /// exponential towards zero rather than towards an unrepresentable magnitude. Past the point + /// where it falls below the last digit being carried the answer is ±1 exactly, and is + /// returned without the exponential being taken at all — which is what keeps 1e15 below + /// from throwing, since e^-2e15 needs an exponent far outside an . + /// + [TestMethod] + public void TestTanhSaturatesRatherThanOverflowing() + { + Assert.AreEqual(PreciseNumber.One, PreciseNumber.Tanh(1000.ToPreciseNumber(), 50)); + Assert.AreEqual(PreciseNumber.NegativeOne, PreciseNumber.Tanh((-1000).ToPreciseNumber(), 50)); + Assert.AreEqual(PreciseNumber.One, PreciseNumber.Tanh(Parse("1E15"), 50)); + Assert.AreEqual(PreciseNumber.NegativeOne, PreciseNumber.Tanh(Parse("-1E15"), 50)); + } + + /// + /// Just below saturation the answer is still strictly inside (-1, 1), so the threshold is + /// pinned from both sides rather than only from the far one. + /// + [TestMethod] + public void TestTanhBelowSaturationIsStrictlyInsideItsRange() + { + PreciseNumber tanh = PreciseNumber.Tanh(20.ToPreciseNumber(), 50); + + Assert.IsLessThan(PreciseNumber.One, tanh, "Tanh(20) reached one"); + AssertAgreesTo(PreciseNumber.One, tanh, 16, "Tanh(20) is not close to one"); + } + + [TestMethod] + public void TestHyperbolicFunctionsRejectValuesOutsideTheirDomain() + { + // There is no NaN and no infinity to return, so the domain is enforced rather than encoded. + Assert.ThrowsExactly(() => PreciseNumber.Acosh(Parse("0.5"))); + Assert.ThrowsExactly(() => PreciseNumber.Acosh(PreciseNumber.Zero)); + Assert.ThrowsExactly(() => PreciseNumber.Acosh(PreciseNumber.NegativeOne)); + Assert.ThrowsExactly(() => PreciseNumber.Atanh(PreciseNumber.One)); + Assert.ThrowsExactly(() => PreciseNumber.Atanh(PreciseNumber.NegativeOne)); + Assert.ThrowsExactly(() => PreciseNumber.Atanh(2.ToPreciseNumber())); + } + + [TestMethod] + public void TestHyperbolicFunctionsRejectAPrecisionBelowOne() + { + Assert.ThrowsExactly(() => PreciseNumber.Sinh(PreciseNumber.One, 0)); + Assert.ThrowsExactly(() => PreciseNumber.Cosh(PreciseNumber.One, 0)); + Assert.ThrowsExactly(() => PreciseNumber.Tanh(PreciseNumber.One, 0)); + Assert.ThrowsExactly(() => PreciseNumber.Asinh(PreciseNumber.One, 0)); + Assert.ThrowsExactly(() => PreciseNumber.Acosh(2.ToPreciseNumber(), 0)); + Assert.ThrowsExactly(() => PreciseNumber.Atanh(PreciseNumber.Zero, 0)); + } + + /// + /// The cheap regression net: agreement with the runtime's own implementations to fifteen digits, + /// which is about all a carries. + /// + [TestMethod] + public void TestHyperbolicFunctionsAgreeWithTheRuntimeToFifteenDigits() + { + foreach (PreciseNumber x in Sweep()) + { + double value = x.To(); + + AssertAgreesTo(Math.Sinh(value).ToPreciseNumber(), PreciseNumber.Sinh(x, 50), 14, $"Sinh disagrees at x = {x}"); + AssertAgreesTo(Math.Cosh(value).ToPreciseNumber(), PreciseNumber.Cosh(x, 50), 14, $"Cosh disagrees at x = {x}"); + AssertAgreesTo(Math.Tanh(value).ToPreciseNumber(), PreciseNumber.Tanh(x, 50), 14, $"Tanh disagrees at x = {x}"); + AssertAgreesTo(Math.Asinh(value).ToPreciseNumber(), PreciseNumber.Asinh(x, 50), 14, $"Asinh disagrees at x = {x}"); + } + } + + /// + /// The inverses have narrower domains than the sweep, so they get their own comparison against + /// the runtime over values each of them accepts. + /// + [TestMethod] + public void TestInverseHyperbolicFunctionsAgreeWithTheRuntimeToFifteenDigits() + { + foreach (PreciseNumber x in new[] { Parse("0.125"), Parse("-0.125"), Parse("0.75"), Parse("-0.75"), Parse("0.9") }) + { + double value = x.To(); + AssertAgreesTo(Math.Atanh(value).ToPreciseNumber(), PreciseNumber.Atanh(x, 50), 14, $"Atanh disagrees at x = {x}"); + } + + foreach (PreciseNumber x in new[] { Parse("1.5"), Parse("2"), Parse("10"), Parse("1000") }) + { + double value = x.To(); + AssertAgreesTo(Math.Acosh(value).ToPreciseNumber(), PreciseNumber.Acosh(x, 50), 14, $"Acosh disagrees at x = {x}"); + } + } + + /// + /// A large argument still produces a result, rather than reaching the exponent limit on the way. + /// + /// + /// The identity is checked at a hundred and forty digits rather than at fifty, and the reason is + /// worth stating because the obvious version of this test is wrong. sinh and cosh + /// of a hundred are both about 1.34e43 and differ by e^-100, so their squares are + /// about 1.8e86 and a difference of one first becomes visible at the eighty-seventh + /// significant digit. At fifty digits the two are the same number and cosh² - sinh² is + /// exactly zero — correctly, not defectively. A hundred and forty leaves the identity fifty-odd + /// digits of room to hold in. + /// + [TestMethod] + public void TestSinhAndCoshOfALargeArgument() + { + PreciseNumber hundred = 100.ToPreciseNumber(); + PreciseNumber sinh = PreciseNumber.Sinh(hundred, 140); + PreciseNumber cosh = PreciseNumber.Cosh(hundred, 140); + + AssertAgreesTo(PreciseNumber.One, (cosh * cosh) - (sinh * sinh), 40, "cosh²x - sinh²x is not one at x = 100"); + AssertAgreesTo(Math.Sinh(100.0).ToPreciseNumber(), sinh, 14, "Sinh(100) disagrees with the runtime"); + } + + /// + /// The default overload produces at least + /// digits, and a wider argument carries its own width through. + /// + [TestMethod] + public void TestDefaultPrecisionFollowsTheArgument() + { + Assert.AreEqual(50, PreciseNumber.Sinh(PreciseNumber.One).SignificantDigits); + Assert.AreEqual(50, PreciseNumber.Cosh(PreciseNumber.One).SignificantDigits); + + PreciseNumber wide = Parse("1.00000000000000000000000000000000000000000000000000000000001"); + Assert.IsGreaterThanOrEqualTo( + wide.SignificantDigits, + PreciseNumber.Sinh(wide).SignificantDigits, + "Sinh narrowed a wider argument"); + } +} diff --git a/PreciseNumber/PreciseNumber.Hyperbolics.cs b/PreciseNumber/PreciseNumber.Hyperbolics.cs new file mode 100644 index 0000000..b544cce --- /dev/null +++ b/PreciseNumber/PreciseNumber.Hyperbolics.cs @@ -0,0 +1,445 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber; + +using System; +using System.Numerics; + +/// +/// Hyperbolic functions and their inverses, none of which route through . +/// +/// +/// Every one of these is the exponential or the logarithm wearing a different name, so none of them +/// carries a series of its own except where the identity it is named for would cancel away the +/// answer. The definitions are reached through , +/// , and +/// rather than restated. +/// +/// Precision follows the exponentials: a result carries the significant digits of its argument, and +/// never fewer than . Every function has an overload taking +/// that count, and the work is carried digits beyond it. +/// +/// +/// Three of these need their textbook form rearranged, and which three is worth being exact about, +/// because the usual floating-point reasoning does not transfer. Addition, subtraction and +/// multiplication are exact here, so a difference of two nearly equal values loses nothing +/// by itself. What loses digits is cancelling against a value that has already been rounded +/// to a working width — which is what , +/// , and +/// all return. +/// +/// +/// That is the test each form has to pass. sinh x = (e^x - e^-x) / 2 fails it: both +/// exponentials come back rounded to the working width and then cancel down to something of the +/// order of x, so a 1e-30 argument keeps about thirty of the fifty digits asked for. +/// So does asinh x = ln(x + √(x² + 1)), where the root arrives rounded, and +/// atanh x = ½ ln((1 + x)/(1 - x)), where the quotient does. Near zero +/// therefore sums its own series, and both inverses subtract +/// their one analytically and hand the remainder to LogP1. Sinh(1e-30) is +/// 1e-30, not zero, and so are Asinh and Atanh of the same. +/// +/// +/// , and +/// pass it as written, and none of them is rearranged for +/// precision. Cosh sums two positive terms and has nothing to cancel at all; the other two +/// are written the way they are for reasons that are not about lost digits, and each gives its own +/// under its own remarks. +/// +/// +/// has no NaN and no infinity, so Acosh below one and +/// Atanh at or beyond one throw where a would quietly return one of +/// those and carry on, the same way does. +/// +/// +public readonly partial record struct PreciseNumber + : IHyperbolicFunctions +{ + /// The message carried by the exception thrown by Acosh below one. + private const string AcoshDomainMessage = "Acosh is only defined for a value of at least one."; + + /// The message carried by the exception thrown by Atanh at or beyond one. + private const string AtanhDomainMessage = "Atanh is only defined for a value strictly between negative one and one."; + + /// + /// Returns the hyperbolic sine of a value. + /// + /// The value. + /// The hyperbolic sine of . + /// Thrown when the result needs an exponent outside the range of an . + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. sinh 0 is exactly zero. + /// + public static PreciseNumber Sinh(PreciseNumber x) => + Sinh(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the hyperbolic sine of a value, to a chosen number of significant digits. + /// + /// The value. + /// The number of significant digits to produce. + /// The hyperbolic sine of . + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// Near zero this is Σ x^(2n+1)/(2n+1)! rather than (e^x - e^-x) / 2, because that + /// difference cancels away exactly the digits the function is being asked for: both terms + /// approach one while their difference approaches 2x, so a 1e-30 argument would + /// lose thirty digits before the halving. Away from zero the two exponentials differ by enough + /// that the identity costs nothing, and one reciprocal is cheaper than a second series. + /// + public static PreciseNumber Sinh(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + + if (Abs(x) <= DirectSeriesLimit) + { + return SinhSeries(x, working).ReduceSignificance(significantDigits); + } + + PreciseNumber raised = Exp(x, working); + PreciseNumber lowered = Divide(One, raised, working); + return Divide(Subtract(raised, lowered), Two, working).ReduceSignificance(significantDigits); + } + + /// + /// Returns the hyperbolic cosine of a value. + /// + /// The value. + /// The hyperbolic cosine of . + /// Thrown when the result needs an exponent outside the range of an . + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. cosh 0 is exactly one. + /// + public static PreciseNumber Cosh(PreciseNumber x) => + Cosh(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the hyperbolic cosine of a value, to a chosen number of significant digits. + /// + /// The value. + /// The number of significant digits to produce. + /// The hyperbolic cosine of . + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// (e^x + e^-x) / 2 everywhere, with no series of its own and no small-argument case. This + /// is a sum of two positive terms, so unlike there is + /// nothing for it to cancel against: the answer approaches one as the argument approaches zero, + /// which is where the precision is wanted relative to, and the identity delivers it there. + /// + public static PreciseNumber Cosh(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return One; + } + + int working = significantDigits + ExponentialGuardDigits; + PreciseNumber raised = Exp(x, working); + PreciseNumber lowered = Divide(One, raised, working); + return Divide(Add(raised, lowered), Two, working).ReduceSignificance(significantDigits); + } + + /// + /// Returns the hyperbolic tangent of a value. + /// + /// The value. + /// The hyperbolic tangent of , in (-1, 1). + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. tanh 0 is exactly zero. + /// + public static PreciseNumber Tanh(PreciseNumber x) => + Tanh(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the hyperbolic tangent of a value, to a chosen number of significant digits. + /// + /// The value. + /// The number of significant digits to produce. + /// The hyperbolic tangent of , in (-1, 1). + /// Thrown when is less than one. + /// + /// tanh x = -t / (2 + t) with t = expm1(-2x), which is the quotient + /// (1 - e^-2x) / (1 + e^-2x) with the numerator's subtraction done by ExpM1 rather + /// than against a literal one. + /// + /// The part that earns its keep is the sign of the exponent, not the ExpM1. Taking the + /// exponential of minus twice the magnitude means it decays towards zero for a large + /// argument instead of growing, so this form saturates where a tanh written the obvious + /// way round would need e^2x and overflow. Past the point where e^-2|x| falls below + /// the requested precision the answer is ±1 to every digit asked for, and is returned + /// without taking an exponential that would only confirm it. + /// + /// + /// The ExpM1 is the smaller point: 1 - e^-2x would in fact hold its digits here, + /// because the subtraction is exact and nothing has been rounded before it. Reaching for + /// ExpM1 costs nothing and keeps the numerator from depending on that, but unlike + /// it is not repairing a real loss. + /// + /// + public static PreciseNumber Tanh(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + + if (IsBeyondTanhSaturation(x, working)) + { + return x.Significand.Sign > 0 ? One : -One; + } + + PreciseNumber shifted = ExpM1(Multiply(new(0, -2), x), working); + return Divide(-shifted, Add(Two, shifted), working).ReduceSignificance(significantDigits); + } + + /// + /// Returns the inverse hyperbolic sine of a value. + /// + /// The value. + /// The value whose hyperbolic sine is . + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. asinh 0 is exactly zero. + /// + public static PreciseNumber Asinh(PreciseNumber x) => + Asinh(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the inverse hyperbolic sine of a value, to a chosen number of significant digits. + /// + /// The value. + /// The number of significant digits to produce. + /// The value whose hyperbolic sine is . + /// Thrown when is less than one. + /// + /// asinh x = ln(x + √(x² + 1)), rearranged so the one is subtracted analytically: + /// √(1 + x²) - 1 is x² / (1 + √(1 + x²)), so the whole logarithm is + /// LogP1( x + x² / (1 + √(1 + x²)) ). Written the first way, a small argument asks for the + /// logarithm of a number indistinguishable from one at the precision it was given; written this + /// way, the argument handed to LogP1 is of the order of x itself. + /// + /// The function is odd, and is evaluated on the magnitude with the sign restored afterwards. + /// That is not only economy: x + √(x² + 1) for a large negative x is a difference + /// of two nearly equal numbers, and taking the magnitude first is what avoids it. + /// + /// + public static PreciseNumber Asinh(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + 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 result = LogP1(Add(magnitude, excess), significantDigits); + + return x.Significand.Sign > 0 ? result : -result; + } + + /// + /// Returns the inverse hyperbolic cosine of a value. + /// + /// The value, which must be at least one. + /// The non-negative value whose hyperbolic cosine is . + /// Thrown when is less than one. + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. acosh 1 is exactly zero. + /// + public static PreciseNumber Acosh(PreciseNumber x) => + Acosh(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the inverse hyperbolic cosine of a value, to a chosen number of significant digits. + /// + /// The value, which must be at least one. + /// The number of significant digits to produce. + /// The non-negative value whose hyperbolic cosine is . + /// + /// Thrown when is less than one, or when + /// is less than one. + /// + /// + /// acosh x = ln(x + √(x² - 1)), rearranged as + /// LogP1( (x - 1) + √((x - 1)(x + 1)) ). + /// + /// In fixed-precision arithmetic that rearrangement is a precision fix, because x² - 1 + /// just above one cancels the digits it is about to take the root of. Here it is not: squaring is + /// exact and so is the subtraction, so the unfactored form keeps its digits too. What the + /// factored form saves is width. x² has twice the significand of x, and an argument + /// a hundred digits wide would carry a two-hundred-digit intermediate into the root for a result + /// wanted at fifty — which on a type whose digits live in a is paid for + /// in allocation. Factoring the difference of squares also keeps the result independent of + /// exactness rather than resting on it. + /// + /// + /// acosh of a value below one is a value this type cannot represent, and throws rather + /// than returning something wrong. + /// + /// + public static PreciseNumber Acosh(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x < One) + { + throw new ArgumentOutOfRangeException(nameof(x), x, AcoshDomainMessage); + } + + if (x.IsUnit) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + PreciseNumber below = Subtract(x, One); + PreciseNumber root = Sqrt(Multiply(below, Add(x, One)), working); + return LogP1(Add(below, root), significantDigits); + } + + /// + /// Returns the inverse hyperbolic tangent of a value. + /// + /// The value, which must lie strictly between negative one and one. + /// The value whose hyperbolic tangent is . + /// Thrown when the magnitude of is at least one. + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. atanh 0 is exactly zero. + /// + public static PreciseNumber Atanh(PreciseNumber x) => + Atanh(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the inverse hyperbolic tangent of a value, to a chosen number of significant digits. + /// + /// The value, which must lie strictly between negative one and one. + /// The number of significant digits to produce. + /// The value whose hyperbolic tangent is . + /// + /// Thrown when the magnitude of is at least one, or when + /// is less than one. + /// + /// + /// atanh x = ½ ln((1 + x)/(1 - x)), written as ½ LogP1( 2x / (1 - x) ) because + /// (1 + x)/(1 - x) is 1 + 2x/(1 - x) and the quotient approaches one near zero. + /// The rearranged argument approaches 2x instead, which is what LogP1 is for. + /// + /// At ±1 the value is unbounded, and has no infinity, so the + /// endpoints throw along with everything past them. + /// + /// + public static PreciseNumber Atanh(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (Abs(x) >= One) + { + throw new ArgumentOutOfRangeException(nameof(x), x, AtanhDomainMessage); + } + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + PreciseNumber ratio = Divide(Multiply(Two, x), Subtract(One, x), working); + return Divide(LogP1(ratio, working), Two, working).ReduceSignificance(significantDigits); + } + + /// + /// Sums the hyperbolic sine series. + /// + /// The argument. + /// The significant digits to carry through the sum. + /// sinh x. + /// Thrown when the series does not converge. + /// + /// Σ x^(2n+1)/(2n+1)! from zero, each term built from its predecessor by multiplying in + /// x² and dividing by 2n(2n + 1). The leading term is x itself, so the sum + /// is of the order of x and its digits are kept relative to x rather than to one — + /// the same reason omits its leading one. + /// + /// Only called with an argument no larger than , where the + /// factorial outruns the power immediately and the sum is a handful of terms. + /// + /// + private static PreciseNumber SinhSeries(PreciseNumber x, int workingDigits) + { + if (x.Significand.IsZero) + { + return Zero; + } + + PreciseNumber square = Multiply(x, x).ReduceSignificance(workingDigits); + PreciseNumber term = x; + PreciseNumber sum = x; + + for (int n = 1; n <= SeriesIterationAllowance(workingDigits); n++) + { + // The divisor outgrows an int well before the allowance does, so it is formed as a long. + long even = 2L * n; + term = Divide(Multiply(term, square), new(0, even * (even + 1)), workingDigits); + if (term.Significand.IsZero) + { + return sum; + } + + PreciseNumber next = Add(sum, term).ReduceSignificance(workingDigits); + if (next == sum) + { + return sum; + } + + sum = next; + } + + throw new ArithmeticException( + $"The hyperbolic sine series did not converge to {workingDigits.ToString(InvariantCulture)} significant digits."); + } + + /// + /// Reports whether a hyperbolic tangent has saturated at the requested precision. + /// + /// The argument. + /// The significant digits being carried. + /// when tanh of is ±1 to every digit asked for. + /// + /// tanh differs from ±1 by about 2·e^-2|x|, so once 2|x| passes + /// (workingDigits + 1) · ln 10 the difference is below the last digit being carried. The + /// threshold only has to be recognised, not resolved, so the constant behind it is read at the + /// stored precision rather than at the working one. + /// + private static bool IsBeyondTanhSaturation(PreciseNumber x, int workingDigits) => + Multiply(Two, Abs(x)) > Multiply(new(0, workingDigits + 1), Ln10To(MinimumDivisionPrecision)); +}