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..bdc6e13 --- /dev/null +++ b/PreciseNumber.Test/PreciseNumberExponentialTests.cs @@ -0,0 +1,355 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Test; + +using System.Globalization; +using System.Linq; +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"), + ]; + + /// + /// 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() + { + 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() + { + // 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))) + { + 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 (PreciseNumber value in InversionSweep()) + { + 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