diff --git a/PreciseNumber.Test/PreciseNumberConstantTests.cs b/PreciseNumber.Test/PreciseNumberConstantTests.cs new file mode 100644 index 0000000..330b46a --- /dev/null +++ b/PreciseNumber.Test/PreciseNumberConstantTests.cs @@ -0,0 +1,304 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Test; + +using System.Globalization; +using System.Numerics; + +/// +/// Pins the mathematical constants against independent computations of the same values. +/// +/// +/// Every constant is checked digit for digit against a series that shares nothing with the literal +/// in the source, so an edit that corrupts a literal cannot pass. Each series is evaluated in +/// fixed point with guard digits, then rounded the same way the literal +/// was. +/// +[TestClass] +public class PreciseNumberConstantTests +{ + /// + /// Digits carried beyond while a series is summed, + /// so that truncation in the guard region cannot reach the digits under test. + /// + private const int GuardDigits = 40; + + /// + /// The fixed point scale every series below is evaluated in. A value v is held as + /// round(v * Scale). + /// + private static readonly BigInteger Scale = BigInteger.Pow(10, PreciseNumber.ConstantPrecision + GuardDigits); + + [TestMethod] + public void TestPiMatchesMachinsFormula() + { + // pi = 16*arctan(1/5) - 4*arctan(1/239) + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); + PreciseNumber expected = RoundToSignificantDigits(pi, PreciseNumber.ConstantPrecision); + + Assert.AreEqual(expected, PreciseNumber.Pi, "Pi does not match Machin's formula"); + Assert.AreEqual(expected.Significand, PreciseNumber.Pi.Significand); + Assert.AreEqual(expected.Exponent, PreciseNumber.Pi.Exponent); + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.Pi.SignificantDigits); + } + + [TestMethod] + public void TestPiIsRoundedRatherThanTruncated() + { + // pi = 3.14159265358979323846264338327950288..., so a 26 digit pi ends 434 when it is + // rounded and 433 when it is truncated. The literal used to end 433. + PreciseNumber pi26 = PreciseNumber.PiTo(26); + + Assert.AreEqual( + "31415926535897932384626434", + pi26.Significand.ToString(CultureInfo.InvariantCulture), + "Pi truncates where it should round"); + + // The same at full precision: the 151st significant digit of pi is 8, so the 150th rounds + // up from 2 to 3. + string digits = PreciseNumber.Pi.Significand.ToString(CultureInfo.InvariantCulture); + Assert.AreEqual('3', digits[^1], "The last digit of Pi is not rounded up"); + } + + [TestMethod] + public void TestTauMatchesMachinsFormulaDoubled() + { + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); + PreciseNumber expected = RoundToSignificantDigits(2 * pi, PreciseNumber.ConstantPrecision); + + Assert.AreEqual(expected, PreciseNumber.Tau, "Tau does not match Machin's formula doubled"); + Assert.AreEqual(expected.Significand, PreciseNumber.Tau.Significand); + Assert.AreEqual(expected.Exponent, PreciseNumber.Tau.Exponent); + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.Tau.SignificantDigits); + } + + [TestMethod] + public void TestTauIsExactlyPiDoubled() + { + // The two are independent literals. Nothing but this assertion stops them drifting apart. + Assert.AreEqual(PreciseNumber.Pi * 2.ToPreciseNumber(), PreciseNumber.Tau); + } + + [TestMethod] + public void TestEMatchesItsSeries() + { + // e = sum 1/k! + BigInteger e = Scale; + BigInteger term = Scale; + for (int k = 1; term != BigInteger.Zero; k++) + { + term /= k; + e += term; + } + + PreciseNumber expected = RoundToSignificantDigits(e, PreciseNumber.ConstantPrecision); + + Assert.AreEqual(expected, PreciseNumber.E, "E does not match its series"); + Assert.AreEqual(expected.Significand, PreciseNumber.E.Significand); + Assert.AreEqual(expected.Exponent, PreciseNumber.E.Exponent); + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.E.SignificantDigits); + } + + [TestMethod] + public void TestLn2MatchesItsSeries() + { + // ln(2) = 2*artanh(1/3) + PreciseNumber expected = RoundToSignificantDigits(2 * ArtanhReciprocal(3), PreciseNumber.ConstantPrecision); + + Assert.AreEqual(expected, PreciseNumber.Ln2, "Ln2 does not match its series"); + Assert.AreEqual(expected.Significand, PreciseNumber.Ln2.Significand); + Assert.AreEqual(expected.Exponent, PreciseNumber.Ln2.Exponent); + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.Ln2.SignificantDigits); + } + + [TestMethod] + public void TestLn10MatchesItsSeries() + { + // ln(10) = ln(8) + ln(10/8) = 6*artanh(1/3) + 2*artanh(1/9) + PreciseNumber expected = RoundToSignificantDigits( + (6 * ArtanhReciprocal(3)) + (2 * ArtanhReciprocal(9)), + PreciseNumber.ConstantPrecision); + + Assert.AreEqual(expected, PreciseNumber.Ln10, "Ln10 does not match its series"); + Assert.AreEqual(expected.Significand, PreciseNumber.Ln10.Significand); + Assert.AreEqual(expected.Exponent, PreciseNumber.Ln10.Exponent); + + // The 150th significant digit of ln(10) is a zero, which the constructor strips along with + // any other trailing zero. The value is still correct to 150 digits. + Assert.AreEqual(PreciseNumber.ConstantPrecision - 1, PreciseNumber.Ln10.SignificantDigits); + } + + [TestMethod] + public void TestConstantsCarryAtLeastMinimumDivisionPrecision() + { + // A constant shorter than this caps any expression that mixes it with a quotient, silently. + foreach (PreciseNumber constant in new[] + { + PreciseNumber.E, + PreciseNumber.Pi, + PreciseNumber.Tau, + PreciseNumber.Ln2, + PreciseNumber.Ln10, + }) + { + Assert.IsGreaterThanOrEqualTo( + PreciseNumber.MinimumDivisionPrecision, + constant.SignificantDigits, + $"a constant carries only {constant.SignificantDigits} significant digits"); + } + } + + [TestMethod] + public void TestQuotientOfPiIsCorrectPastTheOldPrecision() + { + // pi/3 to 50 digits needs a pi of at least 50 digits. The 26 digit literal could not do it. + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); + PreciseNumber expected = RoundToSignificantDigits(pi / 3, PreciseNumber.MinimumDivisionPrecision); + + PreciseNumber actual = (PreciseNumber.Pi / 3.ToPreciseNumber()) + .ReduceSignificance(PreciseNumber.MinimumDivisionPrecision); + + Assert.AreEqual(expected, actual); + } + + [TestMethod] + public void TestPiToReducesToTheRequestedPrecision() + { + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); + + for (int digits = 1; digits <= PreciseNumber.ConstantPrecision; digits++) + { + PreciseNumber expected = RoundToSignificantDigits(pi, digits); + Assert.AreEqual(expected, PreciseNumber.PiTo(digits), $"PiTo({digits}) is not correctly rounded"); + } + } + + [TestMethod] + public void TestConstantAccessorsReduceToTheRequestedPrecision() + { + Assert.AreEqual(PreciseNumber.E.ReduceSignificance(20), PreciseNumber.ETo(20)); + Assert.AreEqual(PreciseNumber.Pi.ReduceSignificance(20), PreciseNumber.PiTo(20)); + Assert.AreEqual(PreciseNumber.Tau.ReduceSignificance(20), PreciseNumber.TauTo(20)); + Assert.AreEqual(PreciseNumber.Ln2.ReduceSignificance(20), PreciseNumber.Ln2To(20)); + Assert.AreEqual(PreciseNumber.Ln10.ReduceSignificance(20), PreciseNumber.Ln10To(20)); + + Assert.AreEqual(20, PreciseNumber.PiTo(20).SignificantDigits); + } + + [TestMethod] + public void TestConstantAccessorsServeRepeatedRequestsIdentically() + { + // The second call comes from the cache; it has to be the same number as the first. + Assert.AreEqual(PreciseNumber.PiTo(30), PreciseNumber.PiTo(30)); + Assert.AreEqual(PreciseNumber.PiTo(30).Significand, PreciseNumber.PiTo(30).Significand); + Assert.AreEqual(PreciseNumber.PiTo(30).Exponent, PreciseNumber.PiTo(30).Exponent); + + // Two precisions of the same constant must not collide in that cache. + Assert.AreNotEqual(PreciseNumber.PiTo(30), PreciseNumber.PiTo(20)); + + // Nor may two constants asked for the same precision. + Assert.AreNotEqual(PreciseNumber.PiTo(30), PreciseNumber.TauTo(30)); + } + + [TestMethod] + public void TestConstantAccessorsReturnTheWholeConstantWhenAskedForMore() + { + Assert.AreEqual(PreciseNumber.Pi, PreciseNumber.PiTo(PreciseNumber.ConstantPrecision)); + Assert.AreEqual(PreciseNumber.Pi, PreciseNumber.PiTo(PreciseNumber.ConstantPrecision + 100)); + Assert.AreEqual(PreciseNumber.E, PreciseNumber.ETo(int.MaxValue)); + Assert.AreEqual(PreciseNumber.Tau, PreciseNumber.TauTo(int.MaxValue)); + Assert.AreEqual(PreciseNumber.Ln2, PreciseNumber.Ln2To(int.MaxValue)); + Assert.AreEqual(PreciseNumber.Ln10, PreciseNumber.Ln10To(int.MaxValue)); + } + + [TestMethod] + public void TestConstantAccessorsRejectFewerThanOneDigit() + { + Assert.ThrowsExactly(() => PreciseNumber.ETo(0)); + Assert.ThrowsExactly(() => PreciseNumber.PiTo(0)); + Assert.ThrowsExactly(() => PreciseNumber.TauTo(-1)); + Assert.ThrowsExactly(() => PreciseNumber.Ln2To(-1)); + Assert.ThrowsExactly(() => PreciseNumber.Ln10To(int.MinValue)); + } + + /// + /// Sums arctan(1/n) = sum (-1)^k / ((2k+1) n^(2k+1)) in fixed point. + /// + /// The reciprocal of the argument. + /// arctan(1/n) * Scale. + private static BigInteger ArctanReciprocal(int n) + { + BigInteger total = Scale / n; + BigInteger term = total; + BigInteger squared = (BigInteger)n * n; + int k = 1; + + while (term != BigInteger.Zero) + { + term /= squared; + k += 2; + total += k % 4 == 3 ? -(term / k) : term / k; + } + + return total; + } + + /// + /// Sums artanh(1/n) = sum 1 / ((2k+1) n^(2k+1)) in fixed point. + /// + /// The reciprocal of the argument. + /// artanh(1/n) * Scale. + private static BigInteger ArtanhReciprocal(int n) + { + BigInteger total = Scale / n; + BigInteger term = total; + BigInteger squared = (BigInteger)n * n; + int k = 1; + + while (term != BigInteger.Zero) + { + term /= squared; + k += 2; + total += term / k; + } + + return total; + } + + /// + /// Rounds a fixed point value to a number of significant digits, half away from zero. + /// + /// The value, scaled by . + /// The number of significant digits to keep. + /// The rounded value. + private static PreciseNumber RoundToSignificantDigits(BigInteger value, int significantDigits) + { + int integerDigits = DigitCount(value) - DigitCount(Scale) + 1; + int shift = significantDigits - integerDigits; + + // One digit beyond the ones being kept, to round on. + BigInteger scaled = value * BigInteger.Pow(10, shift + 1) / Scale; + BigInteger rounded = BigInteger.DivRem(scaled, 10, out BigInteger remainder); + if (BigInteger.Abs(remainder) >= 5) + { + rounded += value.Sign; + } + + int exponent = -shift; + if (DigitCount(rounded) > significantDigits) + { + rounded /= 10; + exponent++; + } + + return PreciseNumber.CreateFromComponents(exponent, rounded); + } + + /// + /// Counts the decimal digits of a . + /// + /// The value to count. + /// The number of decimal digits, ignoring any sign. + private static int DigitCount(BigInteger value) => + BigInteger.Abs(value).ToString(CultureInfo.InvariantCulture).Length; +} diff --git a/PreciseNumber.Test/PreciseNumberTests.cs b/PreciseNumber.Test/PreciseNumberTests.cs index f5cf994..e53650f 100644 --- a/PreciseNumber.Test/PreciseNumberTests.cs +++ b/PreciseNumber.Test/PreciseNumberTests.cs @@ -1881,37 +1881,8 @@ public void TestParseRejectsExponentsOutsideIntRange() Assert.AreEqual(int.MaxValue, PreciseNumber.Parse("1E2147483647", CultureInfo.InvariantCulture).Exponent); } - [TestMethod] - public void TestEValue() - { - BigInteger expectedSignificand = BigInteger.Parse("27182818284590452353602874713526624977572", CultureInfo.InvariantCulture); - int expectedExponent = -40; - PreciseNumber eValue = PreciseNumber.E; - - Assert.AreEqual(expectedSignificand, eValue.Significand); - Assert.AreEqual(expectedExponent, eValue.Exponent); - } - - [TestMethod] - public void TestTauValue() - { - BigInteger expectedSignificand = BigInteger.Parse("6283185307179586476925287", CultureInfo.InvariantCulture); - int expectedExponent = -24; - PreciseNumber tauValue = PreciseNumber.Tau; - - Assert.AreEqual(expectedSignificand, tauValue.Significand); - Assert.AreEqual(expectedExponent, tauValue.Exponent); - } - - [TestMethod] - public void TestPiValue() - { - BigInteger expectedSignificand = BigInteger.Parse("31415926535897932384626433", CultureInfo.InvariantCulture); - int expectedExponent = -25; - PreciseNumber piValue = PreciseNumber.Pi; - Assert.AreEqual(expectedSignificand, piValue.Significand); - Assert.AreEqual(expectedExponent, piValue.Exponent); - } + // E, Tau and Pi are pinned by PreciseNumberConstantTests, which checks each of them digit for + // digit against an independent computation rather than against a second copy of the literal. [TestMethod] public void TestNotEqual() diff --git a/PreciseNumber/PreciseNumber.cs b/PreciseNumber/PreciseNumber.cs index 78c9f9e..a41c51a 100644 --- a/PreciseNumber/PreciseNumber.cs +++ b/PreciseNumber/PreciseNumber.cs @@ -4,6 +4,7 @@ namespace ktsu.PreciseNumber; using System; using System.Buffers; +using System.Collections.Concurrent; using System.Diagnostics; using System.Diagnostics.CodeAnalysis; using System.Globalization; @@ -250,26 +251,168 @@ internal PreciseNumber(int exponent, BigInteger significand, bool sanitize) /// public static PreciseNumber Zero { get; } = new(0, 0); - private const int EExponent = -40; + /// + /// The number of significant digits carried by , , + /// , , and . + /// + /// + /// Chosen to match the 150 digits ktsu.Semantics standardizes on for the conversion + /// factors it derives from pi, so that the two libraries cannot disagree about pi. It is also + /// well above , so an expression mixing a constant with a + /// quotient is no longer capped at the constant's precision. + /// + /// Argument reduction cannot be more accurate than the constant it reduces by: reducing an angle + /// of magnitude 10^d modulo tau to n correct digits spends roughly d of the constant's digits + /// before it starts on the answer. 150 leaves room for that. + /// + /// + public const int ConstantPrecision = 150; + + private const int EExponent = -149; + + /// + /// Gets the value of e for the type, correctly rounded to + /// significant digits. + /// + /// Its own literal, never computed from another constant, so that an error in one cannot reach the others. + public static PreciseNumber E { get; } = new(EExponent, BigInteger.Parse("271828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526", InvariantCulture)); + + private const int PiExponent = -149; /// - /// Gets the value of e for the type. + /// Gets the value of pi for the type, correctly rounded to + /// significant digits. /// - public static PreciseNumber E { get; } = new(EExponent, BigInteger.Parse("27182818284590452353602874713526624977572", InvariantCulture)); + /// Its own literal, never computed from another constant, so that an error in one cannot reach the others. + public static PreciseNumber Pi { get; } = new(PiExponent, BigInteger.Parse("314159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813", InvariantCulture)); - private const int PiExponent = -25; + private const int TauExponent = -149; /// - /// Gets the value of pi for the type. + /// Gets the value of tau for the type, correctly rounded to + /// significant digits. /// - public static PreciseNumber Pi { get; } = new(PiExponent, BigInteger.Parse("31415926535897932384626433", InvariantCulture)); + /// + /// Its own literal, never computed from another constant, so that an error in one cannot reach + /// the others. It nonetheless agrees exactly with doubled, which + /// TestTauIsExactlyPiDoubled pins so that the two cannot drift apart. + /// + public static PreciseNumber Tau { get; } = new(TauExponent, BigInteger.Parse("628318530717958647692528676655900576839433879875021164194988918461563281257241799725606965068423413596429617302656461329418768921910116446345071881626", InvariantCulture)); - private const int TauExponent = -24; + private const int Ln2Exponent = -150; /// - /// Gets the value of tau for the type. + /// Gets the natural logarithm of two, correctly rounded to + /// significant digits. /// - public static PreciseNumber Tau { get; } = new(TauExponent, BigInteger.Parse("6283185307179586476925287", InvariantCulture)); + /// Its own literal, never computed from another constant, so that an error in one cannot reach the others. + public static PreciseNumber Ln2 { get; } = new(Ln2Exponent, BigInteger.Parse("693147180559945309417232121458176568075500134360255254120680009493393621969694715605863326996418687542001481020570685733685520235758130557032670751635", InvariantCulture)); + + private const int Ln10Exponent = -149; + + /// + /// Gets the natural logarithm of ten, correctly rounded to + /// significant digits. + /// + /// + /// Its own literal, never computed from another constant, so that an error in one cannot reach + /// the others. Its 150th significant digit is a zero, which the constructor removes along with + /// any other trailing zero, so it stores 149 digits for the same value. + /// + public static PreciseNumber Ln10 { get; } = new(Ln10Exponent, BigInteger.Parse("230258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834196778404228624863340952546508280675666628736910", InvariantCulture)); + + /// + /// Reduced forms of the constants above, keyed by the constant and the number of significant + /// digits asked of it. + /// + /// + /// Multiplication is exact, so any product involving a digit + /// constant carries at least that many digits. A caller that only wants fifteen pays for all of + /// them unless it asks for fifteen, and asking repeatedly should not re-round every time. + /// + private static readonly ConcurrentDictionary<(PreciseNumber Constant, int SignificantDigits), PreciseNumber> constantCache = new(); + + /// + /// Gets the value of e reduced to the specified number of significant digits. + /// + /// The number of significant digits to produce. + /// + /// rounded to significant digits, or + /// itself when that is no fewer digits than it carries. + /// + /// Thrown when is less than one. + public static PreciseNumber ETo(int significantDigits) => + ConstantTo(E, significantDigits); + + /// + /// Gets the value of pi reduced to the specified number of significant digits. + /// + /// The number of significant digits to produce. + /// + /// rounded to significant digits, or + /// itself when that is no fewer digits than it carries. + /// + /// Thrown when is less than one. + public static PreciseNumber PiTo(int significantDigits) => + ConstantTo(Pi, significantDigits); + + /// + /// Gets the value of tau reduced to the specified number of significant digits. + /// + /// The number of significant digits to produce. + /// + /// rounded to significant digits, or + /// itself when that is no fewer digits than it carries. + /// + /// Thrown when is less than one. + public static PreciseNumber TauTo(int significantDigits) => + ConstantTo(Tau, significantDigits); + + /// + /// Gets the natural logarithm of two reduced to the specified number of significant digits. + /// + /// The number of significant digits to produce. + /// + /// rounded to significant digits, or + /// itself when that is no fewer digits than it carries. + /// + /// Thrown when is less than one. + public static PreciseNumber Ln2To(int significantDigits) => + ConstantTo(Ln2, significantDigits); + + /// + /// Gets the natural logarithm of ten reduced to the specified number of significant digits. + /// + /// The number of significant digits to produce. + /// + /// rounded to significant digits, or + /// itself when that is no fewer digits than it carries. + /// + /// Thrown when is less than one. + public static PreciseNumber Ln10To(int significantDigits) => + ConstantTo(Ln10, significantDigits); + + /// + /// Reduces one of the constants to the specified number of significant digits, serving it from + /// a cache. + /// + /// The full precision constant. + /// The number of significant digits to produce. + /// The constant at the requested precision, rounded half away from zero. + /// Thrown when is less than one. + private static PreciseNumber ConstantTo(PreciseNumber constant, int significantDigits) + { + if (significantDigits < 1) + { + throw new ArgumentOutOfRangeException(nameof(significantDigits), significantDigits, "At least one significant digit is required."); + } + + return significantDigits >= constant.SignificantDigits + ? constant + : constantCache.GetOrAdd( + (constant, significantDigits), + static key => key.Constant.ReduceSignificance(key.SignificantDigits)); + } /// /// Gets the exponent of the number. diff --git a/README.md b/README.md index fe4b64f..ad8d07c 100644 --- a/README.md +++ b/README.md @@ -421,6 +421,14 @@ precision of the wider operand, never fewer than `MinimumDivisionPrecision` (50) digits, with the last digit rounded half away from zero. Pass an explicit precision to the three-argument overload when you want something other than that. +`Pi`, `Tau`, `E`, `Ln2` and `Ln10` are each carried to `ConstantPrecision` (150) significant +digits, correctly rounded, and each is its own literal rather than being computed from a sibling. +150 matches what `ktsu.Semantics` standardises on for the factors it derives from pi, and leaves +room for argument reduction, which cannot be more accurate than the constant it reduces by. +Multiplication is exact, so any product involving one of them carries at least 150 digits; a +caller that only needs fifteen should ask for fifteen with `PiTo(15)` and its siblings, which +round half away from zero and cache per requested precision. + ## Limitations - `Exp()`, and `Pow()` with a non-integer power, are computed through `double` and are therefore limited to its precision. Addition, subtraction, multiplication and division are not @@ -433,7 +441,9 @@ three-argument overload when you want something other than that. ### PreciseNumber Class -- **Constants**: `Zero`, `One`, `NegativeOne`, `Pi`, `E`, `Tau` +- **Constants**: `Zero`, `One`, `NegativeOne`, `Pi`, `E`, `Tau`, `Ln2`, `Ln10` + +- **Constants at a chosen precision**: `PiTo()`, `ETo()`, `TauTo()`, `Ln2To()`, `Ln10To()` - **Arithmetic**: `+`, `-`, `*`, `/`, `%`, `++`, `--`