|
| 1 | +// Copyright (c) 2023-2026 ktsu-dev contributors |
| 2 | + |
| 3 | +namespace ktsu.PreciseNumber.Test; |
| 4 | + |
| 5 | +using System.Globalization; |
| 6 | +using System.Numerics; |
| 7 | + |
| 8 | +/// <summary> |
| 9 | +/// Pins the mathematical constants against independent computations of the same values. |
| 10 | +/// </summary> |
| 11 | +/// <remarks> |
| 12 | +/// Every constant is checked digit for digit against a series that shares nothing with the literal |
| 13 | +/// in the source, so an edit that corrupts a literal cannot pass. Each series is evaluated in |
| 14 | +/// <see cref="BigInteger"/> fixed point with guard digits, then rounded the same way the literal |
| 15 | +/// was. |
| 16 | +/// </remarks> |
| 17 | +[TestClass] |
| 18 | +public class PreciseNumberConstantTests |
| 19 | +{ |
| 20 | + /// <summary> |
| 21 | + /// Digits carried beyond <see cref="PreciseNumber.ConstantPrecision"/> while a series is summed, |
| 22 | + /// so that truncation in the guard region cannot reach the digits under test. |
| 23 | + /// </summary> |
| 24 | + private const int GuardDigits = 40; |
| 25 | + |
| 26 | + /// <summary> |
| 27 | + /// The fixed point scale every series below is evaluated in. A value <c>v</c> is held as |
| 28 | + /// <c>round(v * Scale)</c>. |
| 29 | + /// </summary> |
| 30 | + private static readonly BigInteger Scale = BigInteger.Pow(10, PreciseNumber.ConstantPrecision + GuardDigits); |
| 31 | + |
| 32 | + [TestMethod] |
| 33 | + public void TestPiMatchesMachinsFormula() |
| 34 | + { |
| 35 | + // pi = 16*arctan(1/5) - 4*arctan(1/239) |
| 36 | + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); |
| 37 | + PreciseNumber expected = RoundToSignificantDigits(pi, PreciseNumber.ConstantPrecision); |
| 38 | + |
| 39 | + Assert.AreEqual(expected, PreciseNumber.Pi, "Pi does not match Machin's formula"); |
| 40 | + Assert.AreEqual(expected.Significand, PreciseNumber.Pi.Significand); |
| 41 | + Assert.AreEqual(expected.Exponent, PreciseNumber.Pi.Exponent); |
| 42 | + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.Pi.SignificantDigits); |
| 43 | + } |
| 44 | + |
| 45 | + [TestMethod] |
| 46 | + public void TestPiIsRoundedRatherThanTruncated() |
| 47 | + { |
| 48 | + // pi = 3.14159265358979323846264338327950288..., so a 26 digit pi ends 434 when it is |
| 49 | + // rounded and 433 when it is truncated. The literal used to end 433. |
| 50 | + PreciseNumber pi26 = PreciseNumber.PiTo(26); |
| 51 | + |
| 52 | + Assert.AreEqual( |
| 53 | + "31415926535897932384626434", |
| 54 | + pi26.Significand.ToString(CultureInfo.InvariantCulture), |
| 55 | + "Pi truncates where it should round"); |
| 56 | + |
| 57 | + // The same at full precision: the 151st significant digit of pi is 8, so the 150th rounds |
| 58 | + // up from 2 to 3. |
| 59 | + string digits = PreciseNumber.Pi.Significand.ToString(CultureInfo.InvariantCulture); |
| 60 | + Assert.AreEqual('3', digits[^1], "The last digit of Pi is not rounded up"); |
| 61 | + } |
| 62 | + |
| 63 | + [TestMethod] |
| 64 | + public void TestTauMatchesMachinsFormulaDoubled() |
| 65 | + { |
| 66 | + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); |
| 67 | + PreciseNumber expected = RoundToSignificantDigits(2 * pi, PreciseNumber.ConstantPrecision); |
| 68 | + |
| 69 | + Assert.AreEqual(expected, PreciseNumber.Tau, "Tau does not match Machin's formula doubled"); |
| 70 | + Assert.AreEqual(expected.Significand, PreciseNumber.Tau.Significand); |
| 71 | + Assert.AreEqual(expected.Exponent, PreciseNumber.Tau.Exponent); |
| 72 | + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.Tau.SignificantDigits); |
| 73 | + } |
| 74 | + |
| 75 | + [TestMethod] |
| 76 | + public void TestTauIsExactlyPiDoubled() |
| 77 | + { |
| 78 | + // The two are independent literals. Nothing but this assertion stops them drifting apart. |
| 79 | + Assert.AreEqual(PreciseNumber.Pi * 2.ToPreciseNumber(), PreciseNumber.Tau); |
| 80 | + } |
| 81 | + |
| 82 | + [TestMethod] |
| 83 | + public void TestEMatchesItsSeries() |
| 84 | + { |
| 85 | + // e = sum 1/k! |
| 86 | + BigInteger e = Scale; |
| 87 | + BigInteger term = Scale; |
| 88 | + for (int k = 1; term != BigInteger.Zero; k++) |
| 89 | + { |
| 90 | + term /= k; |
| 91 | + e += term; |
| 92 | + } |
| 93 | + |
| 94 | + PreciseNumber expected = RoundToSignificantDigits(e, PreciseNumber.ConstantPrecision); |
| 95 | + |
| 96 | + Assert.AreEqual(expected, PreciseNumber.E, "E does not match its series"); |
| 97 | + Assert.AreEqual(expected.Significand, PreciseNumber.E.Significand); |
| 98 | + Assert.AreEqual(expected.Exponent, PreciseNumber.E.Exponent); |
| 99 | + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.E.SignificantDigits); |
| 100 | + } |
| 101 | + |
| 102 | + [TestMethod] |
| 103 | + public void TestLn2MatchesItsSeries() |
| 104 | + { |
| 105 | + // ln(2) = 2*artanh(1/3) |
| 106 | + PreciseNumber expected = RoundToSignificantDigits(2 * ArtanhReciprocal(3), PreciseNumber.ConstantPrecision); |
| 107 | + |
| 108 | + Assert.AreEqual(expected, PreciseNumber.Ln2, "Ln2 does not match its series"); |
| 109 | + Assert.AreEqual(expected.Significand, PreciseNumber.Ln2.Significand); |
| 110 | + Assert.AreEqual(expected.Exponent, PreciseNumber.Ln2.Exponent); |
| 111 | + Assert.AreEqual(PreciseNumber.ConstantPrecision, PreciseNumber.Ln2.SignificantDigits); |
| 112 | + } |
| 113 | + |
| 114 | + [TestMethod] |
| 115 | + public void TestLn10MatchesItsSeries() |
| 116 | + { |
| 117 | + // ln(10) = ln(8) + ln(10/8) = 6*artanh(1/3) + 2*artanh(1/9) |
| 118 | + PreciseNumber expected = RoundToSignificantDigits( |
| 119 | + (6 * ArtanhReciprocal(3)) + (2 * ArtanhReciprocal(9)), |
| 120 | + PreciseNumber.ConstantPrecision); |
| 121 | + |
| 122 | + Assert.AreEqual(expected, PreciseNumber.Ln10, "Ln10 does not match its series"); |
| 123 | + Assert.AreEqual(expected.Significand, PreciseNumber.Ln10.Significand); |
| 124 | + Assert.AreEqual(expected.Exponent, PreciseNumber.Ln10.Exponent); |
| 125 | + |
| 126 | + // The 150th significant digit of ln(10) is a zero, which the constructor strips along with |
| 127 | + // any other trailing zero. The value is still correct to 150 digits. |
| 128 | + Assert.AreEqual(PreciseNumber.ConstantPrecision - 1, PreciseNumber.Ln10.SignificantDigits); |
| 129 | + } |
| 130 | + |
| 131 | + [TestMethod] |
| 132 | + public void TestConstantsCarryAtLeastMinimumDivisionPrecision() |
| 133 | + { |
| 134 | + // A constant shorter than this caps any expression that mixes it with a quotient, silently. |
| 135 | + foreach (PreciseNumber constant in new[] |
| 136 | + { |
| 137 | + PreciseNumber.E, |
| 138 | + PreciseNumber.Pi, |
| 139 | + PreciseNumber.Tau, |
| 140 | + PreciseNumber.Ln2, |
| 141 | + PreciseNumber.Ln10, |
| 142 | + }) |
| 143 | + { |
| 144 | + Assert.IsGreaterThanOrEqualTo( |
| 145 | + PreciseNumber.MinimumDivisionPrecision, |
| 146 | + constant.SignificantDigits, |
| 147 | + $"a constant carries only {constant.SignificantDigits} significant digits"); |
| 148 | + } |
| 149 | + } |
| 150 | + |
| 151 | + [TestMethod] |
| 152 | + public void TestQuotientOfPiIsCorrectPastTheOldPrecision() |
| 153 | + { |
| 154 | + // pi/3 to 50 digits needs a pi of at least 50 digits. The 26 digit literal could not do it. |
| 155 | + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); |
| 156 | + PreciseNumber expected = RoundToSignificantDigits(pi / 3, PreciseNumber.MinimumDivisionPrecision); |
| 157 | + |
| 158 | + PreciseNumber actual = (PreciseNumber.Pi / 3.ToPreciseNumber()) |
| 159 | + .ReduceSignificance(PreciseNumber.MinimumDivisionPrecision); |
| 160 | + |
| 161 | + Assert.AreEqual(expected, actual); |
| 162 | + } |
| 163 | + |
| 164 | + [TestMethod] |
| 165 | + public void TestPiToReducesToTheRequestedPrecision() |
| 166 | + { |
| 167 | + BigInteger pi = (16 * ArctanReciprocal(5)) - (4 * ArctanReciprocal(239)); |
| 168 | + |
| 169 | + for (int digits = 1; digits <= PreciseNumber.ConstantPrecision; digits++) |
| 170 | + { |
| 171 | + PreciseNumber expected = RoundToSignificantDigits(pi, digits); |
| 172 | + Assert.AreEqual(expected, PreciseNumber.PiTo(digits), $"PiTo({digits}) is not correctly rounded"); |
| 173 | + } |
| 174 | + } |
| 175 | + |
| 176 | + [TestMethod] |
| 177 | + public void TestConstantAccessorsReduceToTheRequestedPrecision() |
| 178 | + { |
| 179 | + Assert.AreEqual(PreciseNumber.E.ReduceSignificance(20), PreciseNumber.ETo(20)); |
| 180 | + Assert.AreEqual(PreciseNumber.Pi.ReduceSignificance(20), PreciseNumber.PiTo(20)); |
| 181 | + Assert.AreEqual(PreciseNumber.Tau.ReduceSignificance(20), PreciseNumber.TauTo(20)); |
| 182 | + Assert.AreEqual(PreciseNumber.Ln2.ReduceSignificance(20), PreciseNumber.Ln2To(20)); |
| 183 | + Assert.AreEqual(PreciseNumber.Ln10.ReduceSignificance(20), PreciseNumber.Ln10To(20)); |
| 184 | + |
| 185 | + Assert.AreEqual(20, PreciseNumber.PiTo(20).SignificantDigits); |
| 186 | + } |
| 187 | + |
| 188 | + [TestMethod] |
| 189 | + public void TestConstantAccessorsServeRepeatedRequestsIdentically() |
| 190 | + { |
| 191 | + // The second call comes from the cache; it has to be the same number as the first. |
| 192 | + Assert.AreEqual(PreciseNumber.PiTo(30), PreciseNumber.PiTo(30)); |
| 193 | + Assert.AreEqual(PreciseNumber.PiTo(30).Significand, PreciseNumber.PiTo(30).Significand); |
| 194 | + Assert.AreEqual(PreciseNumber.PiTo(30).Exponent, PreciseNumber.PiTo(30).Exponent); |
| 195 | + |
| 196 | + // Two precisions of the same constant must not collide in that cache. |
| 197 | + Assert.AreNotEqual(PreciseNumber.PiTo(30), PreciseNumber.PiTo(20)); |
| 198 | + |
| 199 | + // Nor may two constants asked for the same precision. |
| 200 | + Assert.AreNotEqual(PreciseNumber.PiTo(30), PreciseNumber.TauTo(30)); |
| 201 | + } |
| 202 | + |
| 203 | + [TestMethod] |
| 204 | + public void TestConstantAccessorsReturnTheWholeConstantWhenAskedForMore() |
| 205 | + { |
| 206 | + Assert.AreEqual(PreciseNumber.Pi, PreciseNumber.PiTo(PreciseNumber.ConstantPrecision)); |
| 207 | + Assert.AreEqual(PreciseNumber.Pi, PreciseNumber.PiTo(PreciseNumber.ConstantPrecision + 100)); |
| 208 | + Assert.AreEqual(PreciseNumber.E, PreciseNumber.ETo(int.MaxValue)); |
| 209 | + Assert.AreEqual(PreciseNumber.Tau, PreciseNumber.TauTo(int.MaxValue)); |
| 210 | + Assert.AreEqual(PreciseNumber.Ln2, PreciseNumber.Ln2To(int.MaxValue)); |
| 211 | + Assert.AreEqual(PreciseNumber.Ln10, PreciseNumber.Ln10To(int.MaxValue)); |
| 212 | + } |
| 213 | + |
| 214 | + [TestMethod] |
| 215 | + public void TestConstantAccessorsRejectFewerThanOneDigit() |
| 216 | + { |
| 217 | + Assert.ThrowsExactly<ArgumentOutOfRangeException>(() => PreciseNumber.ETo(0)); |
| 218 | + Assert.ThrowsExactly<ArgumentOutOfRangeException>(() => PreciseNumber.PiTo(0)); |
| 219 | + Assert.ThrowsExactly<ArgumentOutOfRangeException>(() => PreciseNumber.TauTo(-1)); |
| 220 | + Assert.ThrowsExactly<ArgumentOutOfRangeException>(() => PreciseNumber.Ln2To(-1)); |
| 221 | + Assert.ThrowsExactly<ArgumentOutOfRangeException>(() => PreciseNumber.Ln10To(int.MinValue)); |
| 222 | + } |
| 223 | + |
| 224 | + /// <summary> |
| 225 | + /// Sums <c>arctan(1/n) = sum (-1)^k / ((2k+1) n^(2k+1))</c> in fixed point. |
| 226 | + /// </summary> |
| 227 | + /// <param name="n">The reciprocal of the argument.</param> |
| 228 | + /// <returns><c>arctan(1/n) * Scale</c>.</returns> |
| 229 | + private static BigInteger ArctanReciprocal(int n) |
| 230 | + { |
| 231 | + BigInteger total = Scale / n; |
| 232 | + BigInteger term = total; |
| 233 | + BigInteger squared = (BigInteger)n * n; |
| 234 | + int k = 1; |
| 235 | + |
| 236 | + while (term != BigInteger.Zero) |
| 237 | + { |
| 238 | + term /= squared; |
| 239 | + k += 2; |
| 240 | + total += k % 4 == 3 ? -(term / k) : term / k; |
| 241 | + } |
| 242 | + |
| 243 | + return total; |
| 244 | + } |
| 245 | + |
| 246 | + /// <summary> |
| 247 | + /// Sums <c>artanh(1/n) = sum 1 / ((2k+1) n^(2k+1))</c> in fixed point. |
| 248 | + /// </summary> |
| 249 | + /// <param name="n">The reciprocal of the argument.</param> |
| 250 | + /// <returns><c>artanh(1/n) * Scale</c>.</returns> |
| 251 | + private static BigInteger ArtanhReciprocal(int n) |
| 252 | + { |
| 253 | + BigInteger total = Scale / n; |
| 254 | + BigInteger term = total; |
| 255 | + BigInteger squared = (BigInteger)n * n; |
| 256 | + int k = 1; |
| 257 | + |
| 258 | + while (term != BigInteger.Zero) |
| 259 | + { |
| 260 | + term /= squared; |
| 261 | + k += 2; |
| 262 | + total += term / k; |
| 263 | + } |
| 264 | + |
| 265 | + return total; |
| 266 | + } |
| 267 | + |
| 268 | + /// <summary> |
| 269 | + /// Rounds a fixed point value to a number of significant digits, half away from zero. |
| 270 | + /// </summary> |
| 271 | + /// <param name="value">The value, scaled by <see cref="Scale"/>.</param> |
| 272 | + /// <param name="significantDigits">The number of significant digits to keep.</param> |
| 273 | + /// <returns>The rounded value.</returns> |
| 274 | + private static PreciseNumber RoundToSignificantDigits(BigInteger value, int significantDigits) |
| 275 | + { |
| 276 | + int integerDigits = DigitCount(value) - DigitCount(Scale) + 1; |
| 277 | + int shift = significantDigits - integerDigits; |
| 278 | + |
| 279 | + // One digit beyond the ones being kept, to round on. |
| 280 | + BigInteger scaled = value * BigInteger.Pow(10, shift + 1) / Scale; |
| 281 | + BigInteger rounded = BigInteger.DivRem(scaled, 10, out BigInteger remainder); |
| 282 | + if (BigInteger.Abs(remainder) >= 5) |
| 283 | + { |
| 284 | + rounded += value.Sign; |
| 285 | + } |
| 286 | + |
| 287 | + int exponent = -shift; |
| 288 | + if (DigitCount(rounded) > significantDigits) |
| 289 | + { |
| 290 | + rounded /= 10; |
| 291 | + exponent++; |
| 292 | + } |
| 293 | + |
| 294 | + return PreciseNumber.CreateFromComponents(exponent, rounded); |
| 295 | + } |
| 296 | + |
| 297 | + /// <summary> |
| 298 | + /// Counts the decimal digits of a <see cref="BigInteger"/>. |
| 299 | + /// </summary> |
| 300 | + /// <param name="value">The value to count.</param> |
| 301 | + /// <returns>The number of decimal digits, ignoring any sign.</returns> |
| 302 | + private static int DigitCount(BigInteger value) => |
| 303 | + BigInteger.Abs(value).ToString(CultureInfo.InvariantCulture).Length; |
| 304 | +} |
0 commit comments