From 943f4bf555f38513af295a0dfdf7e5a6c9b7b2a6 Mon Sep 17 00:00:00 2001 From: Matthew Edmondson Date: Mon, 21 Sep 2026 16:42:56 +0000 Subject: [PATCH 1/2] Implement IExponentialFunctions, ILogarithmicFunctions and IPowerFunctions [minor] Exp and non-integer Pow computed Math.Exp(Math.Log(x.To()) * y), so a fractional power of a fifty-digit value came back with about fifteen correct digits wearing a fifty-digit type, with nothing in the signature, the return value or the documentation to say so. Both are now computed on the significand, following the precision rule Divide and the roots already follow: the significant digits of the argument, never fewer than MinimumDivisionPrecision, with an overload to choose that count. The base-10 representation does most of the work. ln(m . 10^k) splits into ln m + k . ln 10 against the stored Ln10, with the mantissa centred on [1/root10, root10) and fed to the atanh series. exp(v) factors out 10^round(v / ln 10) as an exponent shift, then halves what is left until it is under 1/64 and squares the sum back afterwards. Two constraints in the design are load-bearing rather than stylistic, and each has a test that fails on the naive implementation: - Exp10 and Log10 do not route through the natural log. Exp10 of an integer is an exponent and no series at all; Log10 of a power of ten returns that exponent exactly. - The M1 and P1 variants are computed directly, not as Exp(x) - 1 and Log(1 + x), which would cancel away the precision near zero they exist to keep. ExpM1(1e-30) is 1e-30, not 0. The integer Pow path is untouched and stays exact; the tests pin that with equality against Squared() and Cubed() rather than with a tolerance. A fractional power of a negative value now throws ArgumentOutOfRangeException instead of returning a NaN round trip. Four existing tests pinned the old double-precision answers, one of them documenting the 17th digit as coming from binary rounding. They now carry the published digits. Fixes #81 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01JVVMyN8uWMrqzKf3iCd84T --- CLAUDE.md | 6 +- .../ExponentialBenchmarks.cs | 106 ++ .../PreciseNumberExponentialTests.cs | 340 ++++++ PreciseNumber.Test/PreciseNumberTests.cs | 15 +- PreciseNumber/PreciseNumber.Exponentials.cs | 1010 +++++++++++++++++ PreciseNumber/PreciseNumber.cs | 35 +- README.md | 15 +- 7 files changed, 1498 insertions(+), 29 deletions(-) create mode 100644 PreciseNumber.Benchmarks/ExponentialBenchmarks.cs create mode 100644 PreciseNumber.Test/PreciseNumberExponentialTests.cs create mode 100644 PreciseNumber/PreciseNumber.Exponentials.cs diff --git a/CLAUDE.md b/CLAUDE.md index 03dc40e..6ab5460 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -36,8 +36,10 @@ dotnet run -c Release --project PreciseNumber.Benchmarks -- --filter '*' --job s - Factory methods `CreateFromInteger()` and `CreateFromFloatingPoint()` handle type-specific conversion logic - Addition, subtraction and modulus align exponents before calculating; multiplication and division work on the significands directly -- `Divide` is exact when the quotient terminates, and otherwise rounds to a precision that never falls below the wider operand or `MinimumDivisionPrecision`. `Exp` and non-integer `Pow` still route through `double` +- `Divide` is exact when the quotient terminates, and otherwise rounds to a precision that never falls below the wider operand or `MinimumDivisionPrecision` - 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 - 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 @@ -45,7 +47,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, 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 — 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/ExponentialBenchmarks.cs b/PreciseNumber.Benchmarks/ExponentialBenchmarks.cs new file mode 100644 index 0000000..d689623 --- /dev/null +++ b/PreciseNumber.Benchmarks/ExponentialBenchmarks.cs @@ -0,0 +1,106 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Benchmarks; + +using BenchmarkDotNet.Attributes; + +/// +/// Measures the exponentials, the logarithms, and the fractional power built on them. +/// +/// +/// These replaced a round trip, so the cost of correctness belongs on the +/// record rather than being discovered later. and +/// are that record: they are what Exp and a fractional +/// Pow used to do, and they answer in about fifteen correct digits whatever the +/// Digits axis says, so read them as a floor on the measurement rather than as an +/// alternative. +/// +/// The Digits axis drives both the operand and the digits asked of the answer, since the +/// default precision follows the operand. Both series are summed at the requested precision, so +/// the cost grows with the digit count twice over: more terms are needed, and each term is a +/// wider . Read it against +/// 's division, which is the operation inside both loops. +/// +/// +/// and are the cases the +/// representation answers for free — an exponent shift and an exponent read. They should not move +/// with the Digits axis at all, and should allocate nothing beyond the one significand. A +/// regression there means a series is being run where none is needed. +/// +/// +[MemoryDiagnoser] +public class ExponentialBenchmarks +{ + private PreciseNumber value = PreciseNumber.Zero; + private PreciseNumber exponent = PreciseNumber.Zero; + private PreciseNumber fractionalPower = PreciseNumber.Zero; + private PreciseNumber integerPower = PreciseNumber.Zero; + private PreciseNumber powerOfTen = PreciseNumber.Zero; + private PreciseNumber tenExponent = PreciseNumber.Zero; + private double valueAsDouble; + private double exponentAsDouble; + + /// + /// 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() + { + // Kept near one so that the logarithm's mantissa reduction and the exponential's power-of-ten + // reduction are both exercised without either dominating. + value = Operands.Number(Digits, -(Digits - 1)); + exponent = Operands.Number(Digits, -Digits, offset: 11); + fractionalPower = Operands.Number(Digits, -Digits, offset: 23); + integerPower = 64.ToPreciseNumber(); + powerOfTen = PreciseNumber.Parse("1E50", System.Globalization.CultureInfo.InvariantCulture); + tenExponent = 50.ToPreciseNumber(); + valueAsDouble = value.To(); + exponentAsDouble = exponent.To(); + } + + /// Takes the natural logarithm. + /// The result. + [Benchmark] + public PreciseNumber Log() => PreciseNumber.Log(value); + + /// Raises e to a power. + /// The result. + [Benchmark] + public PreciseNumber Exp() => PreciseNumber.Exp(exponent); + + /// Raises a value to a fractional power, which is an exponential of a logarithm. + /// The result. + [Benchmark] + public PreciseNumber PowFractional() => value.Pow(fractionalPower); + + /// Raises a value to an integer power, which is exponentiation by squaring and exact. + /// The result. + [Benchmark] + public PreciseNumber PowInteger() => value.Pow(integerPower); + + /// Raises ten to an integer power, which is an exponent shift and no series. + /// The result. + [Benchmark] + public PreciseNumber Exp10OfAnInteger() => PreciseNumber.Exp10(tenExponent); + + /// Takes the base-10 logarithm of a power of ten, which is an exponent read and no series. + /// The result. + [Benchmark] + public PreciseNumber Log10OfAPowerOfTen() => PreciseNumber.Log10(powerOfTen); + + /// The logarithm these replaced, as a floor on the measurement. + /// The result. + [Benchmark(Baseline = true)] + public double DoubleLogBaseline() => Math.Log(valueAsDouble); + + /// The exponential these replaced, as a floor on the measurement. + /// The result. + [Benchmark] + public double DoubleExpBaseline() => Math.Exp(exponentAsDouble); +} diff --git a/PreciseNumber.Test/PreciseNumberExponentialTests.cs b/PreciseNumber.Test/PreciseNumberExponentialTests.cs new file mode 100644 index 0000000..625c325 --- /dev/null +++ b/PreciseNumber.Test/PreciseNumberExponentialTests.cs @@ -0,0 +1,340 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Test; + +using System.Globalization; +using System.Numerics; + +/// +/// Covers , and +/// on . +/// +/// +/// The digit-for-digit assertions carry published values rather than values this library produced, +/// so a change that makes the series agree with themselves but not with mathematics still fails. +/// +/// Several of these assertions exist specifically to pin the type's precision claim. Before these +/// functions were implemented on the significand, a fractional power went out to +/// and and came back with about +/// fifteen correct digits wearing a fifty-digit type. Any assertion here that compares fifty +/// published digits against a fractional power or an exponential fails outright on that +/// implementation — it does not merely drift in the last place. +/// +/// +[TestClass] +public class PreciseNumberExponentialTests +{ + /// + /// The first fifty significant digits of the natural logarithm of two, ten and three halves. + /// + private const string Ln2Digits = "69314718055994530941723212145817656807550013436026"; + private const string Ln10Digits = "23025850929940456840179914546843642076011014886288"; + private const string Ln1Point5Digits = "40546510810816438197801311546434913657199042346249"; + + /// + /// The first fifty significant digits of e squared, and of e to the power of minus 3.7. + /// + private const string ESquaredDigits = "73890560989306502272304274605750078131803155705518"; + private const string ExpNegative3Point7Digits = "24723526470339391202757382983402629344505070337871"; + + /// + /// The first fifty significant digits of two to the half, ten to the third, and two to the 3.5. + /// + private const string Root2Digits = "14142135623730950488016887242096980785696718753769"; + private const string CubeRoot10Digits = "21544346900318837217592935665193504952593449421921"; + private const string TwoToThreeAndAHalfDigits = "11313708498984760390413509793677584628557375003016"; + + /// + /// The first fifty significant digits of the base-10 logarithm of two, and the base-2 logarithm of ten. + /// + private const string Log10Of2Digits = "30102999566398119521373889472449302676818988146211"; + private const string Log2Of10Digits = "33219280948873623478703194294893901758648313930246"; + + 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.IsTrue( + difference <= tolerance, + $"{message}: expected {expected}, got {actual}, which differs by {difference}"); + } + + /// + /// A spread of magnitudes and digit counts, including values whose significands are far wider + /// than a can hold. + /// + private static PreciseNumber[] Sweep() => + [ + Parse("1.0000000000000000000000000000000000000000000000001"), + Parse("1.5"), + Parse("2"), + Parse("3.1622776601683793319988935444327185337195551393252"), + Parse("7"), + Parse("9.9999999999999999999999999999999999999999999999999"), + Parse("123.456789012345678901234567890123456789012345"), + Parse("0.000000000000000000000000000000000000000000001234567"), + Parse("6.02214076E23"), + Parse("1E-300"), + Parse("1E300"), + ]; + + [TestMethod] + public void TestLogMatchesPublishedDigits() + { + Assert.AreEqual(Ln2Digits, Digits(PreciseNumber.Log(2.ToPreciseNumber(), 50)), "Log(2) is wrong"); + Assert.AreEqual(Ln10Digits, Digits(PreciseNumber.Log(10.ToPreciseNumber(), 50)), "Log(10) is wrong"); + Assert.AreEqual(Ln1Point5Digits, Digits(PreciseNumber.Log(Parse("1.5"), 50)), "Log(1.5) is wrong"); + } + + [TestMethod] + public void TestLogAgreesWithTheStoredConstants() + { + // The constants were sourced independently of the series, so agreeing with them at fifty + // digits is a check on the series rather than on itself. + Assert.AreEqual(PreciseNumber.Ln2To(50), PreciseNumber.Log(2.ToPreciseNumber(), 50), "Log(2) disagrees with Ln2"); + Assert.AreEqual(PreciseNumber.Ln10To(50), PreciseNumber.Log(10.ToPreciseNumber(), 50), "Log(10) disagrees with Ln10"); + } + + [TestMethod] + public void TestLogPlacesTheDecimalPoint() + { + Assert.AreEqual("0.69314718055994530941723212145817656807550013436026", PreciseNumber.Log(2.ToPreciseNumber(), 50).ToString()); + Assert.AreEqual("0.6931471806", PreciseNumber.Log(2.ToPreciseNumber(), 10).ToString()); + } + + [TestMethod] + public void TestLogOfOneIsExactlyZero() => + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Log(PreciseNumber.One)); + + [TestMethod] + public void TestLogRejectsValuesOutsideItsDomain() + { + // There is no NaN and no infinity to return, so the domain is enforced rather than encoded. + Assert.ThrowsExactly(() => PreciseNumber.Log(PreciseNumber.Zero)); + Assert.ThrowsExactly(() => PreciseNumber.Log(PreciseNumber.NegativeOne)); + Assert.ThrowsExactly(() => PreciseNumber.Log(2.ToPreciseNumber(), 0)); + } + + [TestMethod] + public void TestExpMatchesPublishedDigits() + { + Assert.AreEqual(ESquaredDigits, Digits(PreciseNumber.Exp(2.ToPreciseNumber(), 50)), "Exp(2) is wrong"); + Assert.AreEqual(ExpNegative3Point7Digits, Digits(PreciseNumber.Exp(Parse("-3.7"), 50)), "Exp(-3.7) is wrong"); + } + + [TestMethod] + public void TestExpOfZeroAndOne() + { + Assert.AreEqual(PreciseNumber.One, PreciseNumber.Exp(PreciseNumber.Zero)); + Assert.AreEqual(PreciseNumber.ETo(50), PreciseNumber.Exp(PreciseNumber.One, 50)); + } + + [TestMethod] + public void TestExpAndLogRoundTripAcrossASweep() + { + // The logarithm is carried wider than the answer is wanted, because Exp inverts an absolute + // error into a relative one: ln(1E-300) has three integer digits, so a fifty-digit logarithm + // only pins forty-seven digits of the value it came from. + foreach (PreciseNumber value in Sweep()) + { + PreciseNumber roundTripped = PreciseNumber.Exp(PreciseNumber.Log(value, 60), 60); + AssertAgreesTo(value, roundTripped, 49, $"Exp(Log({value})) did not return its argument"); + } + } + + [TestMethod] + public void TestLogAndExpRoundTripAcrossASweep() + { + foreach (PreciseNumber value in Sweep()) + { + // Log of an exponential, rather than the other way round, so the range reduction inside + // Exp is the thing being inverted. + PreciseNumber exponent = PreciseNumber.Log(value, 50); + AssertAgreesTo(exponent, PreciseNumber.Log(PreciseNumber.Exp(exponent, 50), 50), 48, $"Log(Exp({exponent})) did not return its argument"); + } + } + + [TestMethod] + public void TestPowWithAnIntegerExponentStaysExact() + { + // The integer path is exact today and must not be traded away for the fractional one. An + // exact answer has no tolerance to compare against, so this is equality, not agreement. + foreach (PreciseNumber value in Sweep()) + { + Assert.AreEqual(value.Squared(), value.Pow(2.ToPreciseNumber()), $"Pow({value}, 2) is not exactly its square"); + Assert.AreEqual(value.Cubed(), value.Pow(3.ToPreciseNumber()), $"Pow({value}, 3) is not exactly its cube"); + } + + Assert.AreEqual(Parse("1024"), 2.ToPreciseNumber().Pow(10.ToPreciseNumber()), "Pow(2, 10) is not exactly 1024"); + } + + [TestMethod] + public void TestPowWithAFractionalExponentMatchesPublishedDigits() + { + Assert.AreEqual(Root2Digits, Digits(2.ToPreciseNumber().Pow(Parse("0.5"))), "Pow(2, 0.5) is wrong"); + Assert.AreEqual(CubeRoot10Digits, Digits(PreciseNumber.Pow(10.ToPreciseNumber(), PreciseNumber.Divide(PreciseNumber.One, 3.ToPreciseNumber(), 55)).ReduceSignificance(50)), "Pow(10, 1/3) is wrong"); + Assert.AreEqual(TwoToThreeAndAHalfDigits, Digits(2.ToPreciseNumber().Pow(Parse("3.5"))), "Pow(2, 3.5) is wrong"); + } + + [TestMethod] + public void TestPowWithAFractionalExponentAgreesWithTheRoots() + { + // Two independent routes to the same answer: the integer Newton root on the significand, and + // exp(y · ln x). Fifteen of these digits would agree under a double fallback; forty-eight + // only agree if neither route went near one. + foreach (PreciseNumber value in Sweep()) + { + AssertAgreesTo(PreciseNumber.Sqrt(value, 50), value.Pow(Parse("0.5")), 48, $"Pow({value}, 0.5) disagrees with Sqrt"); + } + } + + [TestMethod] + public void TestPowRejectsAFractionalPowerOfANegativeValue() + { + // An odd integer power of a negative value is real and still allowed. + Assert.AreEqual(Parse("-8"), Parse("-2").Pow(3.ToPreciseNumber())); + Assert.ThrowsExactly(() => Parse("-2").Pow(Parse("0.5"))); + } + + [TestMethod] + public void TestExpM1KeepsTheDigitsOfASmallArgument() + { + // Computed as Exp(x) - 1 this is exactly zero: every digit of the answer lies below the last + // digit the subtraction kept. + PreciseNumber result = PreciseNumber.ExpM1(Parse("1E-30"), 50); + + Assert.AreNotEqual(PreciseNumber.Zero, result, "ExpM1 of a small argument collapsed to zero"); + AssertAgreesTo(Parse("1.0000000000000000000000000000005E-30"), result, 40, "ExpM1(1E-30) is wrong"); + } + + [TestMethod] + public void TestLogP1KeepsTheDigitsOfASmallArgument() + { + PreciseNumber result = PreciseNumber.LogP1(Parse("1E-30"), 50); + + Assert.AreNotEqual(PreciseNumber.Zero, result, "LogP1 of a small argument collapsed to zero"); + AssertAgreesTo(Parse("9.999999999999999999999999999995E-31"), result, 40, "LogP1(1E-30) is wrong"); + } + + [TestMethod] + public void TestExpM1AndLogP1InvertEachOther() + { + foreach (string text in new[] { "1E-30", "-1E-30", "0.25", "-0.25", "3", "-3" }) + { + PreciseNumber value = Parse(text); + AssertAgreesTo(value, PreciseNumber.LogP1(PreciseNumber.ExpM1(value, 55), 55), 45, $"LogP1(ExpM1({value})) did not return its argument"); + } + } + + [TestMethod] + public void TestExp10OfAnIntegerIsAnExponentAndNothingElse() + { + PreciseNumber result = PreciseNumber.Exp10(50.ToPreciseNumber()); + + Assert.AreEqual(BigInteger.One, result.Significand, "Exp10 of an integer ran a series instead of shifting the exponent"); + Assert.AreEqual(50, result.Exponent); + Assert.AreEqual(1, result.SignificantDigits); + Assert.AreEqual(Parse("1E-30"), PreciseNumber.Exp10(Parse("-30"))); + } + + [TestMethod] + public void TestExp10AndLog10InvertEachOther() + { + AssertAgreesTo(Parse("2"), PreciseNumber.Exp10(PreciseNumber.Log10(2.ToPreciseNumber(), 50), 50), 48, "Exp10(Log10(2)) did not return two"); + Assert.AreEqual(Log10Of2Digits, Digits(PreciseNumber.Log10(2.ToPreciseNumber(), 50)), "Log10(2) is wrong"); + } + + [TestMethod] + public void TestLog10OfAPowerOfTenIsExact() + { + // The exponent is the whole answer, so no series runs and nothing is rounded. + Assert.AreEqual(50.ToPreciseNumber(), PreciseNumber.Log10(Parse("1E50"))); + Assert.AreEqual(Parse("-30"), PreciseNumber.Log10(Parse("1E-30"))); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Log10(PreciseNumber.One)); + Assert.AreEqual(3.ToPreciseNumber(), PreciseNumber.Log10(Parse("1000"))); + } + + [TestMethod] + public void TestLog2MatchesPublishedDigits() + { + Assert.AreEqual(Log2Of10Digits, Digits(PreciseNumber.Log2(10.ToPreciseNumber(), 50)), "Log2(10) is wrong"); + AssertAgreesTo(10.ToPreciseNumber(), PreciseNumber.Log2(Parse("1024"), 50), 48, "Log2(1024) is wrong"); + } + + [TestMethod] + public void TestLogInAChosenBase() + { + AssertAgreesTo(3.ToPreciseNumber(), PreciseNumber.Log(Parse("1000"), 10.ToPreciseNumber(), 50), 48, "Log(1000, 10) is wrong"); + AssertAgreesTo(Parse("0.5"), PreciseNumber.Log(3.ToPreciseNumber(), 9.ToPreciseNumber(), 50), 48, "Log(3, 9) is wrong"); + } + + [TestMethod] + public void TestExp2MatchesPublishedDigits() + { + Assert.AreEqual(TwoToThreeAndAHalfDigits, Digits(PreciseNumber.Exp2(Parse("3.5"), 50)), "Exp2(3.5) is wrong"); + + // An integer exponent goes through repeated squaring, which is exact. + Assert.AreEqual(Parse("1024"), PreciseNumber.Exp2(10.ToPreciseNumber())); + Assert.AreEqual(PreciseNumber.One, PreciseNumber.Exp2(PreciseNumber.Zero)); + } + + [TestMethod] + public void TestExp2M1AndExp10M1KeepTheDigitsOfASmallArgument() + { + // Published values, not the first-order approximations x·ln2 and x·ln10 — those only agree + // to thirty digits, which would let a second-order error through unnoticed. + AssertAgreesTo( + Parse("6.931471805599453094172321214584167945824592350726E-31"), + PreciseNumber.Exp2M1(Parse("1E-30"), 50), + 48, + "Exp2M1(1E-30) is wrong"); + + AssertAgreesTo( + Parse("2.3025850929940456840179914546870151566563406876341E-30"), + PreciseNumber.Exp10M1(Parse("1E-30"), 50), + 48, + "Exp10M1(1E-30) is wrong"); + } + + [TestMethod] + public void TestLog2P1AndLog10P1KeepTheDigitsOfASmallArgument() + { + AssertAgreesTo( + Parse("1.4426950408889634073599246810011707899062014724493E-30"), + PreciseNumber.Log2P1(Parse("1E-30"), 50), + 48, + "Log2P1(1E-30) is wrong"); + + AssertAgreesTo( + Parse("4.3429448190325182765112891891638793505344537988984E-31"), + PreciseNumber.Log10P1(Parse("1E-30"), 50), + 48, + "Log10P1(1E-30) is wrong"); + } + + [TestMethod] + public void TestPrecisionFollowsTheWiderOperandRatherThanDouble() + { + // A fifty-digit answer is the point. A double round trip caps at about seventeen, so this + // counts digits rather than comparing them. + PreciseNumber result = PreciseNumber.Log(Parse("1.234567890123456789012345678901234567890123456789"), 50); + + Assert.AreEqual(50, result.SignificantDigits, "Log did not produce the digits it was asked for"); + Assert.AreEqual(50, PreciseNumber.Exp(Parse("1.5"), 50).SignificantDigits, "Exp did not produce the digits it was asked for"); + } + + [TestMethod] + public void TestExponentialRejectsAnArgumentItCannotRepresent() => + Assert.ThrowsExactly(() => PreciseNumber.Exp(Parse("1E30"))); +} diff --git a/PreciseNumber.Test/PreciseNumberTests.cs b/PreciseNumber.Test/PreciseNumberTests.cs index e53650f..c200c19 100644 --- a/PreciseNumber.Test/PreciseNumberTests.cs +++ b/PreciseNumber.Test/PreciseNumberTests.cs @@ -1617,8 +1617,10 @@ public void PowShouldReturnCorrectValue() Assert.AreEqual(PreciseNumber.One, PreciseNumber.One.Pow(10.ToPreciseNumber())); Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Zero.Pow(10.ToPreciseNumber())); + // A fractional power is exp(y · ln x) on the significand, so it carries the type's own + // precision rather than the seventeen digits the double fallback this replaced returned. result = number.Pow(2.5.ToPreciseNumber()); - expected = 5.656854249492381.ToPreciseNumber(); + expected = PreciseNumber.Parse("5.6568542494923801952067548968387923142786875015078", CultureInfo.InvariantCulture); Assert.AreEqual(expected, result); } @@ -1671,9 +1673,10 @@ public void TestExpWithNegativePower() { PreciseNumber result = PreciseNumber.Exp(-1.ToPreciseNumber()); - // Exp routes through a double, and the result keeps every digit that double needs to round-trip. - // Its 17th digit comes from binary rounding, so it's 3 where 1/e continues 0.36787944117144232159. - PreciseNumber expected = PreciseNumber.Parse("0.36787944117144233", CultureInfo.InvariantCulture); + // Exp is computed on the significand, so these are 1/e's own digits rather than the ones a + // double needs to round-trip. The 17th is 2, where the double fallback this replaced rounded + // it to 3. + PreciseNumber expected = PreciseNumber.Parse("0.36787944117144232159552377016146086744581113103177", CultureInfo.InvariantCulture); Assert.AreEqual(expected, result); Assert.AreEqual(Math.Exp(-1), result.To()); @@ -1683,7 +1686,7 @@ public void TestExpWithNegativePower() public void TestExpWithLargePositivePower() { PreciseNumber result = PreciseNumber.Exp(5.ToPreciseNumber()); - PreciseNumber expected = 148.4131591025766m.ToPreciseNumber(); + PreciseNumber expected = PreciseNumber.Parse("148.41315910257660342111558004055227962348766759388", CultureInfo.InvariantCulture); Assert.AreEqual(expected, result); } @@ -1691,7 +1694,7 @@ public void TestExpWithLargePositivePower() public void TestExpWithLargeNegativePower() { PreciseNumber result = PreciseNumber.Exp(-5.ToPreciseNumber()); - PreciseNumber expected = 0.006737946999085467m.ToPreciseNumber(); + PreciseNumber expected = PreciseNumber.Parse("0.0067379469990854670966360484231484242488495850273551", CultureInfo.InvariantCulture); Assert.AreEqual(expected, result); } diff --git a/PreciseNumber/PreciseNumber.Exponentials.cs b/PreciseNumber/PreciseNumber.Exponentials.cs new file mode 100644 index 0000000..0a3cdf6 --- /dev/null +++ b/PreciseNumber/PreciseNumber.Exponentials.cs @@ -0,0 +1,1010 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber; + +using System; +using System.Numerics; + +/// +/// Exponentials, logarithms and powers, none of which route through . +/// +/// +/// The base-10 representation does half the work. A value is significand × 10^exponent, so +/// splitting it into a mantissa in [1, 10) and a decimal exponent costs nothing: the exponent +/// contributes e · ln 10 to a logarithm by one multiplication against a stored constant, and +/// consumes a whole factor of 10^k in an exponential by shifting the exponent field. Only the +/// mantissa needs a series, and it is always centred on one before it gets there. +/// +/// Precision follows and the roots: a result +/// carries the significant digits of its argument, and never fewer than +/// . Every function has an overload taking that count, and +/// every series runs at digits beyond it so the digit the final +/// rounding decision is made on is itself correct. +/// +/// +/// has no NaN and no infinity, so the logarithm of zero or of a negative +/// value throws where a would quietly return one of those and carry on, the +/// same way does. +/// +/// +public readonly partial record struct PreciseNumber + : IExponentialFunctions, + ILogarithmicFunctions, + IPowerFunctions +{ + /// + /// Digits computed past the ones the caller asked for, so that the digit the final rounding + /// decision is made on is itself correct. + /// + /// + /// Wider than the two a root needs. A root makes one rounding decision at the end; a series makes + /// one per term, and the exponential squares its result back up to eight times, each of which + /// doubles whatever relative error it was handed. + /// + private const int ExponentialGuardDigits = 10; + + /// + /// Times the exponential series may halve its argument before summing. + /// + /// + /// The argument reaching the series is at most ln(10) / 2, so seven halvings always bring + /// it under ; the eighth is slack. + /// + private const int MaximumHalvings = 8; + + /// + /// Integer digits allowed in the argument of an exponential before it is rejected as overflowing. + /// + /// + /// The power of ten an exponential factors out lives in the field, so it + /// has to fit an . This bound only keeps the intermediate arithmetic sane; the + /// range itself is checked exactly once the power of ten is known. + /// + private const int ExponentialArgumentDigitLimit = 20; + + /// The message carried by the exception thrown for the logarithm of a non-positive value. + private const string NonPositiveLogarithmMessage = "A logarithm is only defined for a positive value."; + + /// The message carried by the exception thrown by LogP1 below negative one. + private const string LogP1DomainMessage = "LogP1 is only defined for a value greater than negative one."; + + /// The message carried by the exception thrown for a fractional power of a negative value. + private const string NegativeBaseMessage = "A negative value has no real power with a fractional exponent."; + + /// Two, as the divisor of the halving step and the multiplier of the atanh series. + private static PreciseNumber Two { get; } = new(0, 2); + + /// Ten, as the bound the mantissa of a logarithm is centred against. + private static PreciseNumber Ten { get; } = new(1, 1); + + /// One half, added before flooring to round to the nearest integer. + private static PreciseNumber Half { get; } = new(-1, 5); + + /// + /// The largest power of ten an exponential may factor into the field. + /// + /// + /// Bounded symmetrically, so the one extra value an holds below zero is given + /// up. An exponent at the very edge of the range is unusable for anything that follows anyway. + /// + private static BigInteger MaximumExponentShift { get; } = int.MaxValue; + + /// + /// The magnitude at or below which a series is summed directly rather than reduced first. + /// + /// + /// One half. Below it the M1 and P1 variants have no cancellation to avoid and the + /// series is short; above it the range reduction is worth more than the series it replaces. + /// + private static PreciseNumber DirectSeriesLimit { get; } = Half; + + /// + /// The magnitude the exponential series halves its argument down to before summing. + /// + /// + /// One sixty-fourth. Each term then shrinks by at least a further factor of 64n, so fifty + /// digits take about twenty terms instead of about fifty. + /// + private static PreciseNumber SeriesArgumentLimit { get; } = new(-6, 15625); + + /// + /// Returns the natural logarithm of a value. + /// + /// The value to take the logarithm of, which must be positive. + /// The natural logarithm of . + /// Thrown when is zero or negative. + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. ln 1 is exactly zero. + /// + public static PreciseNumber Log(PreciseNumber x) => + Log(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the natural logarithm of a value, to a chosen number of significant digits. + /// + /// The value to take the logarithm of, which must be positive. + /// The number of significant digits to produce. + /// The natural logarithm of . + /// + /// Thrown when is zero or negative, or when + /// is less than one. + /// + /// + /// ln(m · 10^k) = ln m + k · ln 10, with ln 10 read from + /// rather than computed. The mantissa is centred on + /// [1/√10, √10) and fed to the atanh series, whose argument is then never larger than + /// about 0.52. + /// + public static PreciseNumber Log(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.Sign <= 0) + { + throw new ArgumentOutOfRangeException(nameof(x), x, NonPositiveLogarithmMessage); + } + + if (x.IsUnit) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + (PreciseNumber mantissa, long decimalExponent) = SplitAroundRootTen(x); + PreciseNumber result = LogMantissa(mantissa, working); + + if (decimalExponent != 0) + { + // The exponent's own digits are consumed by the product, so the constant has to be read + // wider than the answer is wanted. + int exponentDigits = CountDigits(new BigInteger(decimalExponent)); + result = Add(result, Multiply(new(0, decimalExponent), Ln10To(working + exponentDigits))); + } + + return result.ReduceSignificance(significantDigits); + } + + /// + /// Returns the logarithm of a value in a chosen base. + /// + /// The value to take the logarithm of, which must be positive. + /// The base of the logarithm, which must be positive and not one. + /// The logarithm of in base . + /// Thrown when either argument is zero or negative. + /// Thrown when is one, whose logarithm is zero. + public static PreciseNumber Log(PreciseNumber x, PreciseNumber newBase) => + Log(x, newBase, Math.Max(DefaultExponentialPrecision(x), DefaultExponentialPrecision(newBase))); + + /// + /// Returns the logarithm of a value in a chosen base, to a chosen number of significant digits. + /// + /// The value to take the logarithm of, which must be positive. + /// The base of the logarithm, which must be positive and not one. + /// The number of significant digits to produce. + /// The logarithm of in base . + /// + /// Thrown when either value is zero or negative, or when is + /// less than one. + /// + /// Thrown when is one, whose logarithm is zero. + public static PreciseNumber Log(PreciseNumber x, PreciseNumber newBase, int significantDigits) + { + RequireSignificantDigits(significantDigits); + int working = significantDigits + ExponentialGuardDigits; + return Divide(Log(x, working), Log(newBase, working), significantDigits); + } + + /// + /// Returns the base-2 logarithm of a value. + /// + /// The value to take the logarithm of, which must be positive. + /// The base-2 logarithm of . + /// Thrown when is zero or negative. + /// + /// Unlike , no part of this is free in a base-10 + /// representation, so an exact power of two still arrives through a series and is correct to the + /// digits asked for rather than exact. + /// + public static PreciseNumber Log2(PreciseNumber x) => + Log2(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the base-2 logarithm of a value, to a chosen number of significant digits. + /// + /// The value to take the logarithm of, which must be positive. + /// The number of significant digits to produce. + /// The base-2 logarithm of . + /// + /// Thrown when is zero or negative, or when + /// is less than one. + /// + public static PreciseNumber Log2(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + int working = significantDigits + ExponentialGuardDigits; + return Divide(Log(x, working), Ln2To(working), significantDigits); + } + + /// + /// Returns the base-10 logarithm of a value. + /// + /// The value to take the logarithm of, which must be positive. + /// The base-10 logarithm of . + /// Thrown when is zero or negative. + /// + /// Does not route through the natural logarithm for the part of the answer the representation + /// already holds. log10(m · 10^k) = k + log10 m, so an exact power of ten returns its own + /// exponent exactly, with no series run at all, and everything else pays for one mantissa. + /// + public static PreciseNumber Log10(PreciseNumber x) => + Log10(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the base-10 logarithm of a value, to a chosen number of significant digits. + /// + /// The value to take the logarithm of, which must be positive. + /// The number of significant digits to produce. + /// The base-10 logarithm of . + /// + /// Thrown when is zero or negative, or when + /// is less than one. + /// + public static PreciseNumber Log10(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.Sign <= 0) + { + throw new ArgumentOutOfRangeException(nameof(x), x, NonPositiveLogarithmMessage); + } + + (PreciseNumber mantissa, long decimalExponent) = SplitAroundRootTen(x); + PreciseNumber wholePart = new(0, decimalExponent); + + // A power of ten is entirely exponent, so there is nothing left for a series to do. + if (mantissa.IsUnit) + { + return wholePart; + } + + int working = significantDigits + ExponentialGuardDigits; + PreciseNumber fraction = Divide(LogMantissa(mantissa, working), Ln10To(working), working); + return Add(wholePart, fraction).ReduceSignificance(significantDigits); + } + + /// + /// Returns the natural logarithm of one plus a value. + /// + /// The value to add to one, which must be greater than negative one. + /// ln(1 + x). + /// Thrown when is negative one or less. + /// + /// Near zero this is computed as 2 · atanh(x / (x + 2)), which never forms 1 + x + /// and so keeps every digit of a small argument. Computing it as Log(One + x) would throw + /// away precisely the precision this function exists to preserve: + /// LogP1(1e-30) is 1e-30, not zero. + /// + public static PreciseNumber LogP1(PreciseNumber x) => + LogP1(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the natural logarithm of one plus a value, to a chosen number of significant digits. + /// + /// The value to add to one, which must be greater than negative one. + /// The number of significant digits to produce. + /// ln(1 + x). + /// + /// Thrown when is negative one or less, or when + /// is less than one. + /// + public static PreciseNumber LogP1(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + if (x <= NegativeOne) + { + throw new ArgumentOutOfRangeException(nameof(x), x, LogP1DomainMessage); + } + + // Away from zero there is nothing to cancel, so the general logarithm is both simpler and + // better conditioned than an atanh argument approaching one. + if (Abs(x) > DirectSeriesLimit) + { + return Log(Add(One, x), significantDigits); + } + + int working = significantDigits + ExponentialGuardDigits; + PreciseNumber z = Divide(x, Add(x, Two), working); + return Multiply(Two, AtanhSeries(z, working)).ReduceSignificance(significantDigits); + } + + /// + /// Returns the base-2 logarithm of one plus a value. + /// + /// The value to add to one, which must be greater than negative one. + /// log2(1 + x). + /// Thrown when is negative one or less. + public static PreciseNumber Log2P1(PreciseNumber x) => + Log2P1(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the base-2 logarithm of one plus a value, to a chosen number of significant digits. + /// + /// The value to add to one, which must be greater than negative one. + /// The number of significant digits to produce. + /// log2(1 + x). + /// + /// Thrown when is negative one or less, or when + /// is less than one. + /// + public static PreciseNumber Log2P1(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + return Divide(LogP1(x, working), Ln2To(working), significantDigits); + } + + /// + /// Returns the base-10 logarithm of one plus a value. + /// + /// The value to add to one, which must be greater than negative one. + /// log10(1 + x). + /// Thrown when is negative one or less. + public static PreciseNumber Log10P1(PreciseNumber x) => + Log10P1(x, DefaultExponentialPrecision(x)); + + /// + /// Returns the base-10 logarithm of one plus a value, to a chosen number of significant digits. + /// + /// The value to add to one, which must be greater than negative one. + /// The number of significant digits to produce. + /// log10(1 + x). + /// + /// Thrown when is negative one or less, or when + /// is less than one. + /// + public static PreciseNumber Log10P1(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + return Divide(LogP1(x, working), Ln10To(working), significantDigits); + } + + /// + /// Returns e raised to a power. + /// + /// The power to raise e to. + /// e^x. + /// 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. + /// + /// e itself is returned at the full it is stored to, + /// rather than capped at the fifty a one-digit argument would otherwise ask for. + /// + /// + public static PreciseNumber Exp(PreciseNumber x) => + x.IsUnit ? E : Exp(x, DefaultExponentialPrecision(x)); + + /// + /// Returns e raised to a power, to a chosen number of significant digits. + /// + /// The power to raise e to. + /// The number of significant digits to produce. + /// e^x. + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// exp(x) = 10^k · exp(r) with k = round(x / ln 10), so the whole power of ten is a + /// shift of the field and the series only ever sees an r no larger + /// than ln(10) / 2, halved further until it is under + /// and squared back afterwards. + /// + public static PreciseNumber Exp(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return One; + } + + // e itself is stored, so there is no reason to rediscover it a term at a time. + if (x.IsUnit) + { + return ETo(significantDigits); + } + + RequireExponentialArgumentInRange(x); + + int working = significantDigits + ExponentialGuardDigits; + BigInteger powerOfTen = RoundToNearestInteger(Divide(x, Ln10To(working), working)); + int shift = ToExponentShift(powerOfTen); + + // Subtracting k · ln 10 cancels the integer digits of x, so the constant and the series both + // have to be carried that much wider than the answer is wanted. + int consumed = CountDigits(powerOfTen); + PreciseNumber remainder = Subtract(x, Multiply(new(0, powerOfTen), Ln10To(working + consumed))); + PreciseNumber series = ExpReduced(remainder, working + consumed); + return ShiftExponent(series, shift).ReduceSignificance(significantDigits); + } + + /// + /// Returns e raised to a power, minus one. + /// + /// The power to raise e to. + /// e^x - 1. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// Near zero this is the exponential series with its leading one omitted rather than + /// Exp(x) - 1, which would cancel away exactly the digits the function exists to keep: + /// ExpM1(1e-30) is 1e-30, not zero. + /// + public static PreciseNumber ExpM1(PreciseNumber x) => + ExpM1(x, DefaultExponentialPrecision(x)); + + /// + /// Returns e raised to a power, minus one, to a chosen number of significant digits. + /// + /// The power to raise e to. + /// The number of significant digits to produce. + /// e^x - 1. + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + public static PreciseNumber ExpM1(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + // Away from zero the subtraction cancels nothing worth keeping, and the range reduction + // inside Exp is worth more than a series run on an unreduced argument. + if (Abs(x) > DirectSeriesLimit) + { + int wide = significantDigits + ExponentialGuardDigits; + return Subtract(Exp(x, wide), One).ReduceSignificance(significantDigits); + } + + return ExpSeriesWithoutLeadingOne(x, significantDigits + ExponentialGuardDigits) + .ReduceSignificance(significantDigits); + } + + /// + /// Returns two raised to a power. + /// + /// The power to raise two to. + /// 2^x. + /// Thrown when the result needs an exponent outside the range of an . + public static PreciseNumber Exp2(PreciseNumber x) => + Exp2(x, DefaultExponentialPrecision(x)); + + /// + /// Returns two raised to a power, to a chosen number of significant digits. + /// + /// The power to raise two to. + /// The number of significant digits to produce. + /// 2^x. + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// An integer power of two is produced by repeated squaring, which is exact. Anything else goes + /// through exp(x · ln 2). + /// + public static PreciseNumber Exp2(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return One; + } + + if (IsInteger(x)) + { + return Two.Pow(x); + } + + return ExpOfProductWithConstant(x, Ln2To(significantDigits + ExponentialGuardDigits), significantDigits); + } + + /// + /// Returns two raised to a power, minus one. + /// + /// The power to raise two to. + /// 2^x - 1. + /// Thrown when the result needs an exponent outside the range of an . + public static PreciseNumber Exp2M1(PreciseNumber x) => + Exp2M1(x, DefaultExponentialPrecision(x)); + + /// + /// Returns two raised to a power, minus one, to a chosen number of significant digits. + /// + /// The power to raise two to. + /// The number of significant digits to produce. + /// 2^x - 1. + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + /// Routed through , so a small argument keeps its digits. + public static PreciseNumber Exp2M1(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + return ExpM1(Multiply(x, Ln2To(working + IntegerDigitCount(x))), significantDigits); + } + + /// + /// Returns ten raised to a power. + /// + /// The power to raise ten to. + /// 10^x. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// Does not route through the natural logarithm. An integer power of ten is an exponent and + /// nothing else — Exp10(50) is one significand and a shift, with no series run at all — + /// and a fractional power pays only for its fractional part. + /// + public static PreciseNumber Exp10(PreciseNumber x) => + Exp10(x, DefaultExponentialPrecision(x)); + + /// + /// Returns ten raised to a power, to a chosen number of significant digits. + /// + /// The power to raise ten to. + /// The number of significant digits to produce. + /// 10^x. + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + public static PreciseNumber Exp10(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return One; + } + + RequireExponentialArgumentInRange(x); + + // Split off the whole power of ten, which the exponent field carries for free. What is left + // is in [0, 1), so an integer argument never reaches a series at all. + BigInteger wholePart = FloorToInteger(x); + int shift = ToExponentShift(wholePart); + PreciseNumber fraction = Subtract(x, new(0, wholePart)); + + int working = significantDigits + ExponentialGuardDigits; + PreciseNumber series = fraction.Significand.IsZero + ? One + : Exp(Multiply(fraction, Ln10To(working)), working); + + return ShiftExponent(series, shift).ReduceSignificance(significantDigits); + } + + /// + /// Returns ten raised to a power, minus one. + /// + /// The power to raise ten to. + /// 10^x - 1. + /// Thrown when the result needs an exponent outside the range of an . + public static PreciseNumber Exp10M1(PreciseNumber x) => + Exp10M1(x, DefaultExponentialPrecision(x)); + + /// + /// Returns ten raised to a power, minus one, to a chosen number of significant digits. + /// + /// The power to raise ten to. + /// The number of significant digits to produce. + /// 10^x - 1. + /// Thrown when is less than one. + /// Thrown when the result needs an exponent outside the range of an . + /// Routed through , so a small argument keeps its digits. + public static PreciseNumber Exp10M1(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + ExponentialGuardDigits; + return ExpM1(Multiply(x, Ln10To(working + IntegerDigitCount(x))), significantDigits); + } + + /// + /// Returns one value raised to the power of another. + /// + /// The base. + /// The exponent. + /// x^y. + /// Thrown when is negative and is not an integer. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// An integer exponent is exact, by repeated squaring. Anything else is + /// exp(y · ln x), carried wide enough that the digits the exponential's range reduction + /// consumes are digits it was given rather than digits it invents. + /// + public static PreciseNumber Pow(PreciseNumber x, PreciseNumber y) => + x.Pow(y); + + /// + /// Computes exp(y · ln x) for a positive base and a non-integer exponent. + /// + /// The base, which must be positive. + /// The exponent. + /// The number of significant digits to produce. + /// x^y. + /// + /// The exponential factors out 10^k, and forming r cancels every integer digit of + /// y · ln x. Those digits therefore have to be present in the logarithm before the + /// exponential asks for them, which is what the second, wider pass buys. + /// + private static PreciseNumber FractionalPow(PreciseNumber x, PreciseNumber y, int significantDigits) + { + int working = significantDigits + ExponentialGuardDigits; + PreciseNumber product = Multiply(y, Log(x, working)); + int consumed = IntegerDigitCount(product); + + if (consumed > 0) + { + product = Multiply(y, Log(x, working + consumed)); + } + + return Exp(product, significantDigits + consumed).ReduceSignificance(significantDigits); + } + + /// + /// Gets the significant digits an exponential or logarithm produces when the caller does not choose. + /// + /// The value being operated on. + /// The significant digits of , or if that is more. + /// + /// The same rule and the roots follow, so that + /// a constant carrying digits does not silently cap the + /// expression at fifty. + /// + private static int DefaultExponentialPrecision(PreciseNumber value) => + Math.Max(value.SignificantDigits, MinimumDivisionPrecision); + + /// + /// Throws when a caller asks for fewer than one significant digit. + /// + /// The requested significant digits. + /// Thrown when is less than one. + private static void RequireSignificantDigits(int significantDigits) + { + if (significantDigits < 1) + { + throw new ArgumentOutOfRangeException(nameof(significantDigits), significantDigits, "At least one significant digit is required."); + } + } + + /// + /// Throws when an exponential's argument is so large that its result cannot be represented. + /// + /// The argument of the exponential. + /// Thrown when has more integer digits than any representable result could need. + /// + /// The power of ten an exponential factors out has to fit the field, which + /// bounds the argument long before this does. Rejecting the absurd cases up front keeps the + /// intermediate arithmetic — in particular the power of ten a floor divides by — small enough to + /// be worth computing. + /// + private static void RequireExponentialArgumentInRange(PreciseNumber x) + { + if (IntegerDigitCount(x) > ExponentialArgumentDigitLimit) + { + throw new OverflowException( + $"An exponential of a value with {IntegerDigitCount(x).ToString(InvariantCulture)} integer digits needs an exponent outside the range of an int."); + } + } + + /// + /// Narrows a power of ten to the exponent shift that carries it. + /// + /// The power of ten to carry in the exponent. + /// The same value as an . + /// Thrown when the power of ten does not fit an . + private static int ToExponentShift(BigInteger powerOfTen) => + BigInteger.Abs(powerOfTen) > MaximumExponentShift + ? throw new OverflowException( + $"A result scaled by 10^{powerOfTen.ToString(InvariantCulture)} needs an exponent outside the range of an int.") + : (int)powerOfTen; + + /// + /// Multiplies a value by a power of ten by moving its exponent rather than its digits. + /// + /// The value to scale. + /// The power of ten to scale by. + /// multiplied by 10^. + /// Thrown when the shifted exponent does not fit an . + private static PreciseNumber ShiftExponent(PreciseNumber value, int shift) + { + if (shift == 0 || value.Significand.IsZero) + { + return value; + } + + long shifted = (long)value.Exponent + shift; + return shifted is < int.MinValue or > int.MaxValue + ? throw new OverflowException( + $"A result scaled by 10^{shift.ToString(InvariantCulture)} needs an exponent outside the range of an int.") + : new((int)shifted, value.Significand); + } + + /// + /// Counts the digits a value carries ahead of its decimal point. + /// + /// The value to measure. + /// The number of integer digits, or zero when the value is less than one in magnitude. + private static int IntegerDigitCount(PreciseNumber value) + { + if (value.Significand.IsZero) + { + return 0; + } + + long decimalExponent = (long)value.Exponent + value.SignificantDigits - 1; + return decimalExponent < 0 ? 0 : (int)Math.Min(decimalExponent + 1, int.MaxValue); + } + + /// + /// Splits a positive value into a mantissa centred on one and the power of ten it was scaled by. + /// + /// The value to split, which must be positive. + /// A mantissa in [1/√10, √10) and the exponent such that their product is . + /// + /// The representation already holds a mantissa in [1, 10), so the split costs no + /// arithmetic at all. Centring it further halves the worst-case series argument, and the test for + /// it is exact: m > √10 exactly when m² > 10, and squaring is exact. + /// + private static (PreciseNumber Mantissa, long DecimalExponent) SplitAroundRootTen(PreciseNumber value) + { + int digits = value.SignificantDigits; + PreciseNumber mantissa = new(-(digits - 1), value.Significand); + long decimalExponent = (long)value.Exponent + digits - 1; + + if (Multiply(mantissa, mantissa) > Ten) + { + mantissa = new(mantissa.Exponent - 1, mantissa.Significand); + decimalExponent++; + } + + return (mantissa, decimalExponent); + } + + /// + /// Computes the natural logarithm of a mantissa centred on one. + /// + /// The mantissa, in [1/√10, √10). + /// The significant digits to carry through the series. + /// The natural logarithm of . + /// + /// ln m = 2 · atanh((m - 1) / (m + 1)). Over the centred interval the argument never + /// exceeds about 0.52, so each term of the series gains a little over half a digit. + /// + private static PreciseNumber LogMantissa(PreciseNumber mantissa, int workingDigits) + { + if (mantissa.IsUnit) + { + return Zero; + } + + PreciseNumber z = Divide(Subtract(mantissa, One), Add(mantissa, One), workingDigits); + return Multiply(Two, AtanhSeries(z, workingDigits)); + } + + /// + /// Sums the inverse hyperbolic tangent series. + /// + /// The argument, whose magnitude must be below one. + /// The significant digits to carry through the sum. + /// atanh z. + /// Thrown when the series does not converge. + /// + /// atanh z = z + z³/3 + z⁵/5 + …. Every term is positive when z is, and shares its + /// sign otherwise, so nothing cancels and the sum stops as soon as a term falls below the last + /// digit being carried. + /// + private static PreciseNumber AtanhSeries(PreciseNumber z, int workingDigits) + { + if (z.Significand.IsZero) + { + return Zero; + } + + PreciseNumber zSquared = Multiply(z, z).ReduceSignificance(workingDigits); + PreciseNumber term = z; + PreciseNumber sum = z; + + for (int denominator = 3; denominator <= SeriesIterationAllowance(workingDigits); denominator += 2) + { + term = Multiply(term, zSquared).ReduceSignificance(workingDigits); + if (term.Significand.IsZero) + { + return sum; + } + + PreciseNumber next = Add(sum, Divide(term, new(0, denominator), workingDigits)) + .ReduceSignificance(workingDigits); + + if (next == sum) + { + return sum; + } + + sum = next; + } + + throw new ArithmeticException( + $"The logarithm series did not converge to {workingDigits.ToString(InvariantCulture)} significant digits."); + } + + /// + /// Computes the exponential of a value already reduced below ln(10) / 2. + /// + /// The reduced argument. + /// The significant digits to carry through the series. + /// e^r. + /// + /// Halving the argument before the series and squaring the result back afterwards trades a + /// handful of multiplications for most of the terms. Each squaring doubles whatever relative + /// error it is handed, so the sum is carried one digit wider per halving. + /// + private static PreciseNumber ExpReduced(PreciseNumber r, int workingDigits) + { + if (r.Significand.IsZero) + { + return One; + } + + int halvings = 0; + PreciseNumber reduced = r; + + while (halvings < MaximumHalvings && Abs(reduced) > SeriesArgumentLimit) + { + // Halving terminates, so this is exact whatever precision is asked of it. + reduced = Divide(reduced, Two, workingDigits); + halvings++; + } + + int series = workingDigits + halvings; + PreciseNumber result = Add(One, ExpSeriesWithoutLeadingOne(reduced, series)); + + for (int squaring = 0; squaring < halvings; squaring++) + { + result = Multiply(result, result).ReduceSignificance(series); + } + + return result; + } + + /// + /// Sums the exponential series with its leading one omitted. + /// + /// The argument. + /// The significant digits to carry through the sum. + /// e^x - 1. + /// Thrown when the series does not converge. + /// + /// Σ xⁿ/n! from one. Omitting the leading term is what keeps a small argument: the sum is + /// of the order of itself, so carrying it to + /// significant digits keeps them relative to x rather + /// than relative to one. + /// + private static PreciseNumber ExpSeriesWithoutLeadingOne(PreciseNumber x, int workingDigits) + { + if (x.Significand.IsZero) + { + return Zero; + } + + PreciseNumber term = x; + PreciseNumber sum = x; + + for (int n = 2; n <= SeriesIterationAllowance(workingDigits); n++) + { + term = Divide(Multiply(term, x), new(0, n), 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 exponential series did not converge to {workingDigits.ToString(InvariantCulture)} significant digits."); + } + + /// + /// Computes the exponential of a value scaled by a stored constant. + /// + /// The value to scale. + /// The logarithm of the base being raised. + /// The number of significant digits to produce. + /// e^(x · constant). + private static PreciseNumber ExpOfProductWithConstant(PreciseNumber x, PreciseNumber constant, int significantDigits) => + Exp(Multiply(x, constant), significantDigits); + + /// + /// Gets the terms a series is allowed before it is called non-convergent. + /// + /// The significant digits being carried. + /// The largest loop counter the series may reach. + /// + /// The slowest series here is the atanh one at the edge of its centred interval, where the + /// argument squares to about 0.27 and each term is therefore worth a little over half a + /// digit. Four terms per digit leaves that a wide margin, and the constant covers the short sums + /// where a digit count of one would otherwise allow almost no terms at all. + /// + private static int SeriesIterationAllowance(int workingDigits) => + (workingDigits * 4) + 64; + + /// + /// Rounds a value to the nearest integer, half away from negative infinity. + /// + /// The value to round. + /// The nearest integer to . + private static BigInteger RoundToNearestInteger(PreciseNumber value) => + FloorToInteger(Add(value, Half)); + + /// + /// Takes the largest integer no greater than a value. + /// + /// The value to floor. + /// The floor of . + private static BigInteger FloorToInteger(PreciseNumber value) + { + if (value.Significand.IsZero) + { + return BigInteger.Zero; + } + + if (value.Exponent >= 0) + { + return value.Significand * Pow10(value.Exponent); + } + + BigInteger quotient = BigInteger.DivRem(value.Significand, Pow10(-value.Exponent), out BigInteger remainder); + + // BigInteger division truncates towards zero, so a negative value with anything left over + // has been rounded the wrong way for a floor. + return remainder.Sign < 0 ? quotient - BigInteger.One : quotient; + } +} diff --git a/PreciseNumber/PreciseNumber.cs b/PreciseNumber/PreciseNumber.cs index a41c51a..e1741cb 100644 --- a/PreciseNumber/PreciseNumber.cs +++ b/PreciseNumber/PreciseNumber.cs @@ -1703,6 +1703,13 @@ public static PreciseNumber Round(PreciseNumber value, int decimalDigits) => /// /// The power to raise the number to. /// A new instance of that is the result of raising the current instance to the specified power. + /// Thrown when the current instance is negative and is not an integer. + /// Thrown when the result needs an exponent outside the range of an . + /// + /// An integer power is exact, by repeated squaring. A fractional power is + /// exp(power · ln x), produced to the significant digits of the wider operand and never + /// fewer than . + /// public PreciseNumber Pow(PreciseNumber power) { if (power.Significand.IsZero) @@ -1740,28 +1747,18 @@ public PreciseNumber Pow(PreciseNumber power) return power.Significand.Sign < 0 ? One / result : result; } - // Use logarithm and exponential to support decimal powers - double logValue = Math.Log(To()); - return Math.Exp(logValue * power.To()).ToPreciseNumber(); - } - - /// - /// Returns the result of raising e to the specified power. - /// - /// The power to raise e to. - /// A new instance of that is the result of raising e to the specified power. - public static PreciseNumber Exp(PreciseNumber power) - { - if (power.Significand.IsZero) - { - return One; - } - else if (power.IsUnit) + // A fractional power is exp(y · ln x), which has no real value for a negative base. There are + // no complex results here, so this is rejected rather than quietly returned as a NaN. + if (Significand.Sign < 0) { - return E; + throw new ArgumentOutOfRangeException(nameof(power), power, NegativeBaseMessage); } - return Math.Exp(power.To()).ToPreciseNumber(); + int significantDigits = Math.Max( + Math.Max(SignificantDigits, power.SignificantDigits), + MinimumDivisionPrecision); + + return FractionalPow(this, power, significantDigits); } /// diff --git a/README.md b/README.md index 4bb1a83..7601026 100644 --- a/README.md +++ b/README.md @@ -447,11 +447,22 @@ significant digits of the operand, never fewer than `MinimumDivisionPrecision`. through `double`, so a value outside its range — `1e400`, or `1e-400` — roots as accurately as any other. +`Exp`, `Log`, `Pow` and their base-2 and base-10 siblings follow that rule too, and none of them +goes through `double` either. The representation carries most of the work: `ln(m · 10^k)` splits +into `ln m + k · ln 10` against a stored constant, and an exponential factors its whole power of +ten into the exponent field, so only a mantissa centred on one ever reaches a series. `Log10` of a +power of ten, and `Exp10` of an integer, are exact and run no series at all, and an integer `Pow` +is still exact by repeated squaring. + +The `…M1` and `…P1` variants — `ExpM1`, `LogP1` and their siblings — are computed directly rather +than as `Exp(x) - 1` and `Log(1 + x)`, so they keep the digits of a small argument instead of +cancelling them away: `ExpM1(1e-30)` is `1e-30`, not zero. + ## Limitations -- `Exp()`, and `Pow()` with a non-integer power, are computed through `double` and are therefore limited to its precision. Addition, subtraction, multiplication, division and the roots are not +- There is no NaN, so `Sqrt()` of a negative value, `RootN()` of a negative value at an even degree, `Log()` of a value that is not positive, and `Pow()` of a negative value with a fractional exponent, throw `ArgumentOutOfRangeException` where a `double` would return NaN and carry on -- There is no NaN, so `Sqrt()` of a negative value, and `RootN()` of a negative value at an even degree, throw `ArgumentOutOfRangeException` where a `double` would return NaN and carry on +- There is no infinity, so an exponential whose result needs a decimal exponent outside the range of an `int` throws `OverflowException` rather than saturating - A checked conversion to an integer type or `decimal` throws `OverflowException` when the value is out of range. Conversion to `double`, `float`, or `Half` overflows to infinity instead, as it does for every built-in type From d7b24347aaa02f8afc83e496e691346753af8831 Mon Sep 17 00:00:00 2001 From: Matthew Edmondson Date: Mon, 21 Sep 2026 16:54:18 +0000 Subject: [PATCH 2/2] Project the test sweeps instead of remapping the iteration variable [patch] Two CodeQL "missed opportunity to use Select" findings from github-code-quality on the exponential tests: both loops bound a value and then immediately mapped it to the one actually asserted on. TestLogAndExpRoundTripAcrossASweep now iterates the exponents directly, since the value each came from plays no further part. TestExpM1AndLogP1InvertEachOther takes an InversionSweep() of parsed values rather than parsing strings in the body. A .Select(Parse) over the inline array literal would trip CA1861, so this follows the Sweep() already in the file instead, which reads better here anyway. No assertion changed; 330/330 still pass. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01JVVMyN8uWMrqzKf3iCd84T --- .../PreciseNumberExponentialTests.cs | 27 ++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/PreciseNumber.Test/PreciseNumberExponentialTests.cs b/PreciseNumber.Test/PreciseNumberExponentialTests.cs index 625c325..bdc6e13 100644 --- a/PreciseNumber.Test/PreciseNumberExponentialTests.cs +++ b/PreciseNumber.Test/PreciseNumberExponentialTests.cs @@ -3,6 +3,7 @@ namespace ktsu.PreciseNumber.Test; using System.Globalization; +using System.Linq; using System.Numerics; /// @@ -89,6 +90,21 @@ private static PreciseNumber[] Sweep() => Parse("1E300"), ]; + /// + /// Arguments for the M1/P1 inversion, spanning both signs on either side of the + /// magnitude at which those functions stop summing their own series and defer to + /// Exp/Log. + /// + private static PreciseNumber[] InversionSweep() => + [ + Parse("1E-30"), + Parse("-1E-30"), + Parse("0.25"), + Parse("-0.25"), + Parse("3"), + Parse("-3"), + ]; + [TestMethod] public void TestLogMatchesPublishedDigits() { @@ -156,11 +172,11 @@ public void TestExpAndLogRoundTripAcrossASweep() [TestMethod] public void TestLogAndExpRoundTripAcrossASweep() { - foreach (PreciseNumber value in Sweep()) + // Log of an exponential, rather than the other way round, so the range reduction inside Exp + // is the thing being inverted. The sweep is projected to the exponents themselves, since + // the value it came from plays no further part. + foreach (PreciseNumber exponent in Sweep().Select(value => PreciseNumber.Log(value, 50))) { - // Log of an exponential, rather than the other way round, so the range reduction inside - // Exp is the thing being inverted. - PreciseNumber exponent = PreciseNumber.Log(value, 50); AssertAgreesTo(exponent, PreciseNumber.Log(PreciseNumber.Exp(exponent, 50), 50), 48, $"Log(Exp({exponent})) did not return its argument"); } } @@ -230,9 +246,8 @@ public void TestLogP1KeepsTheDigitsOfASmallArgument() [TestMethod] public void TestExpM1AndLogP1InvertEachOther() { - foreach (string text in new[] { "1E-30", "-1E-30", "0.25", "-0.25", "3", "-3" }) + foreach (PreciseNumber value in InversionSweep()) { - PreciseNumber value = Parse(text); AssertAgreesTo(value, PreciseNumber.LogP1(PreciseNumber.ExpM1(value, 55), 55), 45, $"LogP1(ExpM1({value})) did not return its argument"); } }