diff --git a/PreciseNumber.Test/PreciseNumberRootTests.cs b/PreciseNumber.Test/PreciseNumberRootTests.cs index 1b4c002..43e7a6c 100644 --- a/PreciseNumber.Test/PreciseNumberRootTests.cs +++ b/PreciseNumber.Test/PreciseNumberRootTests.cs @@ -307,4 +307,39 @@ public void TestNegativeDegreeRootRoundsTheReciprocalOfAnUnroundedRoot() // exactly 0.9765625. The true reciprocal is 0.976551…, so it must come from a wider root. Assert.AreEqual(Parse("0.9766"), PreciseNumber.RootN(Parse("1.0486"), -2, 4), "RootN(1.0486, -2) is wrong"); } + + [TestMethod] + public void TestNegativeDegreeRootRoundsOnceFromTheTrueValue() + { + // Rooting and then taking the reciprocal rounded twice before the final rounding, and two + // guard digits could not absorb that near a rounding boundary: 3^(-1/3) is 0.69336127435… + Assert.AreEqual(Parse("0.6933612744"), PreciseNumber.RootN(3.ToPreciseNumber(), -3, 10), "RootN(3, -3, 10) is wrong"); + Assert.AreEqual(Parse("0.2236067977"), PreciseNumber.RootN(20.ToPreciseNumber(), -2, 10), "RootN(20, -2, 10) is wrong"); + Assert.AreEqual(Parse("0.14586"), PreciseNumber.RootN(47.ToPreciseNumber(), -2, 5), "RootN(47, -2, 5) is wrong"); + Assert.AreEqual(Parse("0.29"), PreciseNumber.RootN(41.ToPreciseNumber(), -3, 5), "RootN(41, -3, 5) is wrong"); + Assert.AreEqual(Parse("1.2599210"), PreciseNumber.RootN(Parse("0.5"), -3, 8), "RootN(0.5, -3, 8) is wrong"); + Assert.AreEqual(Parse("0.87358"), PreciseNumber.RootN(Parse("1.5"), -3, 6), "RootN(1.5, -3, 6) is wrong"); + Assert.AreEqual(Parse("-0.6933612744"), PreciseNumber.RootN((-3).ToPreciseNumber(), -3, 10), "RootN(-3, -3, 10) is wrong"); + } + + [TestMethod] + public void TestNegativeDegreeRootsAgreeWithAWideReciprocalAcrossASweep() + { + // The reference goes through the positive-degree root and a division, forty digits wider, so + // it shares no code with the path under test and its own rounding is far from the last digit. + foreach (int degree in new[] { 2, 3, 5 }) + { + foreach (int digits in new[] { 5, 10, 20, 50 }) + { + for (int k = 2; k <= 300; k++) + { + PreciseNumber value = k.ToPreciseNumber(); + int wide = digits + 40; + PreciseNumber expected = PreciseNumber.Divide(PreciseNumber.One, PreciseNumber.RootN(value, degree, wide), wide) + .ReduceSignificance(digits); + Assert.AreEqual(expected, PreciseNumber.RootN(value, -degree, digits), $"RootN({k}, -{degree}, {digits}) is wrong"); + } + } + } + } } diff --git a/PreciseNumber/PreciseNumber.Roots.cs b/PreciseNumber/PreciseNumber.Roots.cs index 465b923..1ad4db6 100644 --- a/PreciseNumber/PreciseNumber.Roots.cs +++ b/PreciseNumber/PreciseNumber.Roots.cs @@ -156,9 +156,14 @@ public static PreciseNumber RootN(PreciseNumber x, int n, int significantDigits) if (n < 0) { - return x.Significand.IsZero - ? throw new DivideByZeroException() - : DivideApproximation(One, RootN(x, -n, significantDigits + RootGuardDigits), significantDigits + RootGuardDigits, significantDigits); + return x.Significand.Sign switch + { + 0 => throw new DivideByZeroException(), + > 0 => PositiveReciprocalRootN(x, -n, significantDigits), + _ => int.IsEvenInteger(n) + ? throw new ArgumentOutOfRangeException(nameof(x), x, NegativeRootMessage) + : -PositiveReciprocalRootN(-x, -n, significantDigits), + }; } if (x.Significand.IsZero || n == 1) @@ -257,6 +262,49 @@ private static PreciseNumber PositiveRootN(PreciseNumber value, int n, int signi : result.ReduceSignificance(significantDigits); } + /// + /// Computes the reciprocal of the n-th root of a positive value. + /// + /// The value to take the root of, which must be positive. + /// The degree of the root, which must be at least one. + /// The number of significant digits to produce. + /// The n-th root of the reciprocal of . + /// Thrown when the root needs a scale or an exponent wider than an . + /// + /// Taking the root and then its reciprocal rounds twice before the final rounding, and no fixed + /// number of guard digits stops that from landing on the wrong side of a rounding boundary. This + /// instead roots the reciprocal in integers: 1 / (s · 10^e) is (10^k / s) · 10^(-k - e), + /// and q = floor(10^k / s) loses nothing a root can see, because the floor of the n-th root + /// of floor(y) is the floor of the n-th root of y. The integer root of q is + /// therefore the true value truncated, and the one rounding to the requested digits is correct. + /// + private static PreciseNumber PositiveReciprocalRootN(PreciseNumber value, int n, int significantDigits) + { + // q has at least digits(10^k) - digits(s) digits, so this leaves n · (d + guard) of them for + // the root to work on, as in PositiveRootN. + int digits = CountDigits(value.Significand); + long scale = ((long)n * (significantDigits + RootGuardDigits)) + digits; + + // Raise the scale until the degree divides the exponent -k - e. + scale += (((-scale - value.Exponent) % n) + n) % n; + + long rootExponent = (-scale - value.Exponent) / n; + if (scale > int.MaxValue || rootExponent is < int.MinValue or > int.MaxValue) + { + throw new OverflowException( + $"A root of degree {(-n).ToString(InvariantCulture)} to {significantDigits.ToString(InvariantCulture)} significant digits needs an exponent outside the range of an int."); + } + + BigInteger quotient = BigInteger.DivRem(Pow10((int)scale), value.Significand, out BigInteger remainder); + BigInteger root = n == 1 ? quotient : IntegerRootN(quotient, n); + PreciseNumber result = new((int)rootExponent, root); + + // Exact only when both the reciprocal and its root were. + return remainder.IsZero && BigInteger.Pow(root, n) == quotient + ? result + : result.ReduceSignificance(significantDigits); + } + /// /// Computes the largest integer whose n-th power does not exceed a value. ///