From cd990ea1e7ed9fb01f03f343a57a71bdaa4c02bf Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 16 Sep 2026 14:37:42 +0000 Subject: [PATCH] Implement IRootFunctions with Sqrt, Cbrt, RootN and Hypot [minor] Every root scales the significand by a power of ten until the degree divides the exponent, then takes an integer Newton root of what is left. Started above the root that iteration is strictly decreasing and stops on the floor of the root exactly, so there is no tolerance to choose and no working precision to carry: a value whose root is exact gets it exactly, however few digits were asked for, the way Divide is exact when a quotient terminates. Nothing routes through double, so 1e400 and 1e-400 root as accurately as 2 does. Precision follows Divide's rule, never below the operand's digits and never below MinimumDivisionPrecision, so rooting a 150 digit constant is not capped at fifty. There is no NaN to return, so Sqrt of a negative value throws rather than answering wrongly. An odd root of a negative value is real and is returned. Fixes #80 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MSvvb3UfRbL65zGQy2fip5 --- CLAUDE.md | 5 +- PreciseNumber.Benchmarks/RootBenchmarks.cs | 72 +++++ PreciseNumber.Test/PreciseNumberRootTests.cs | 302 +++++++++++++++++++ PreciseNumber/PreciseNumber.Roots.cs | 297 ++++++++++++++++++ README.md | 26 +- 5 files changed, 698 insertions(+), 4 deletions(-) create mode 100644 PreciseNumber.Benchmarks/RootBenchmarks.cs create mode 100644 PreciseNumber.Test/PreciseNumberRootTests.cs create mode 100644 PreciseNumber/PreciseNumber.Roots.cs diff --git a/CLAUDE.md b/CLAUDE.md index 71d7425..03dc40e 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -37,6 +37,7 @@ 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` +- 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 - 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 @@ -44,12 +45,12 @@ 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, 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, 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 `PreciseNumber.Benchmarks` is a BenchmarkDotNet suite, one class per area (construction, -comparison, arithmetic, pow, rounding, text, conversion). The library exposes its internals to it +comparison, arithmetic, pow, roots, rounding, text, conversion). The library exposes its internals to it so construction can be measured directly. Most classes are parameterised by `Digits` (8, 30, 200). That axis is the point: digits live in a diff --git a/PreciseNumber.Benchmarks/RootBenchmarks.cs b/PreciseNumber.Benchmarks/RootBenchmarks.cs new file mode 100644 index 0000000..db4728d --- /dev/null +++ b/PreciseNumber.Benchmarks/RootBenchmarks.cs @@ -0,0 +1,72 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Benchmarks; + +using BenchmarkDotNet.Attributes; + +/// +/// Measures the square, cube and n-th roots, and the hypotenuse built on them. +/// +/// +/// A root is the most expensive thing this type does, because it iterates and each iteration +/// divides. Read it against 's division, which is the operation +/// inside the loop: a root that is not several times a division is not iterating enough to be +/// correct, and one that is many times it is iterating more than it needs to. +/// +/// The Digits axis drives both the operand and the digits asked of the answer, since the +/// default precision follows the operand. That makes the first two rows of the axis meet at +/// MinimumDivisionPrecision and cost the same, and it is why the perfect square is the +/// slowest case at 200 digits rather than the fastest: its operand carries twice the digits its +/// root does, so it asks for twice the precision. The degree matters as much as the digit count, +/// so the n-th root is measured at a degree well above the cube root's. +/// +/// +[MemoryDiagnoser] +public class RootBenchmarks +{ + private PreciseNumber value = PreciseNumber.Zero; + private PreciseNumber other = PreciseNumber.Zero; + private PreciseNumber perfectSquare = PreciseNumber.Zero; + + /// + /// Gets or sets the number of significant digits in the operand, and so in the root. + /// + [Params(8, 30, 200)] + public int Digits { get; set; } + + /// + /// Prepares the operands. + /// + [GlobalSetup] + public void Setup() + { + value = Operands.Number(Digits, -10); + other = Operands.Number(Digits, -14, offset: 7); + perfectSquare = value.Squared(); + } + + /// The square root. + /// The root. + [Benchmark(Baseline = true)] + public PreciseNumber Sqrt() => PreciseNumber.Sqrt(value); + + /// The square root of a value that has an exact one, which returns without rounding. + /// The root. + [Benchmark] + public PreciseNumber SqrtOfAPerfectSquare() => PreciseNumber.Sqrt(perfectSquare); + + /// The cube root. + /// The root. + [Benchmark] + public PreciseNumber Cbrt() => PreciseNumber.Cbrt(value); + + /// A root of a degree high enough that raising the estimate dominates the division. + /// The root. + [Benchmark] + public PreciseNumber RootN() => PreciseNumber.RootN(value, 17); + + /// The hypotenuse, which is a square root over an exact sum of squares. + /// The hypotenuse. + [Benchmark] + public PreciseNumber Hypot() => PreciseNumber.Hypot(value, other); +} diff --git a/PreciseNumber.Test/PreciseNumberRootTests.cs b/PreciseNumber.Test/PreciseNumberRootTests.cs new file mode 100644 index 0000000..0c71cad --- /dev/null +++ b/PreciseNumber.Test/PreciseNumberRootTests.cs @@ -0,0 +1,302 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Test; + +using System.Globalization; +using System.Numerics; + +/// +/// Covers on . +/// +/// +/// The digit-for-digit assertions carry published values rather than values this library produced, +/// so a change that makes the roots agree with themselves but not with mathematics still fails. The +/// sweeps then check the property the published values cannot: that the root of an arbitrary value +/// squares back to the value it came from. +/// +[TestClass] +public class PreciseNumberRootTests +{ + /// + /// The first fifty significant digits of the square root of two, three and ten. + /// + private const string Sqrt2Digits = "14142135623730950488016887242096980785696718753769"; + private const string Sqrt3Digits = "17320508075688772935274463415058723669428052538104"; + private const string Sqrt10Digits = "31622776601683793319988935444327185337195551393252"; + + /// + /// The first fifty significant digits of the cube root of two, and of the fifth root of seven. + /// + private const string Cbrt2Digits = "12599210498948731647672106072782283505702514647015"; + private const string Root5Of7Digits = "14757731615945520692769166956322441065440936137402"; + + 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}"); + } + + [TestMethod] + public void TestSqrtMatchesPublishedDigits() + { + Assert.AreEqual(Sqrt2Digits, Digits(PreciseNumber.Sqrt(2.ToPreciseNumber(), 50)), "Sqrt(2) is wrong"); + Assert.AreEqual(Sqrt3Digits, Digits(PreciseNumber.Sqrt(3.ToPreciseNumber(), 50)), "Sqrt(3) is wrong"); + Assert.AreEqual(Sqrt10Digits, Digits(PreciseNumber.Sqrt(10.ToPreciseNumber(), 50)), "Sqrt(10) is wrong"); + } + + [TestMethod] + public void TestSqrtPlacesTheDecimalPoint() + { + Assert.AreEqual("1.4142135623730950488016887242096980785696718753769", PreciseNumber.Sqrt(2.ToPreciseNumber(), 50).ToString()); + Assert.AreEqual("3.1622776601683793319988935444327185337195551393252", PreciseNumber.Sqrt(10.ToPreciseNumber(), 50).ToString()); + Assert.AreEqual("1.414213562", PreciseNumber.Sqrt(2.ToPreciseNumber(), 10).ToString()); + } + + [TestMethod] + public void TestSqrtRoundsTheLastDigitRatherThanTruncatingIt() + { + // Sqrt(3) continues ...0525381038, so the 50th digit is a 3 that rounds up to a 4. A + // truncating implementation leaves it at 3. + Assert.AreEqual('4', Digits(PreciseNumber.Sqrt(3.ToPreciseNumber(), 50))[^1], "Sqrt(3) truncates where it should round"); + } + + [TestMethod] + public void TestSqrtOfAPerfectSquareIsExact() + { + PreciseNumber root = PreciseNumber.Sqrt(144.ToPreciseNumber()); + + Assert.AreEqual(12.ToPreciseNumber(), root, "Sqrt(144) is not 12"); + Assert.AreEqual(2, root.SignificantDigits, "Sqrt(144) carries digits it does not have"); + Assert.AreEqual("12", root.ToString()); + } + + [TestMethod] + public void TestAnExactRootIsExactWhateverPrecisionIsAskedFor() + { + // A perfect square whose root is far longer than the digits requested. An implementation + // that rounds unconditionally answers 1.2e29 for the first of these. + PreciseNumber root = Parse("123456789012345678901234567890"); + PreciseNumber square = root.Squared(); + + Assert.AreEqual(root, PreciseNumber.Sqrt(square, 1), "An exact square root was rounded away"); + Assert.AreEqual(root, PreciseNumber.Sqrt(square, 50), "An exact square root was rounded away"); + } + + [TestMethod] + public void TestSqrtHandlesValuesOutsideTheRangeOfADouble() + { + // 1e400 overflows a double, so an implementation that seeds from one gets these wrong. + PreciseNumber large = PreciseNumber.Sqrt(Parse("2E400"), 50); + PreciseNumber small = PreciseNumber.Sqrt(Parse("2E-400"), 50); + + Assert.AreEqual(Sqrt2Digits, Digits(large), "Sqrt(2e400) is wrong"); + Assert.AreEqual(151, large.Exponent, "Sqrt(2e400) has the wrong exponent"); + + Assert.AreEqual(Sqrt2Digits, Digits(small), "Sqrt(2e-400) is wrong"); + Assert.AreEqual(-249, small.Exponent, "Sqrt(2e-400) has the wrong exponent"); + + // An odd exponent cannot be halved by the exponent alone, so the significand carries it. + PreciseNumber odd = PreciseNumber.Sqrt(Parse("1E401"), 50); + Assert.AreEqual(Sqrt10Digits, Digits(odd), "Sqrt(1e401) is wrong"); + Assert.AreEqual(151, odd.Exponent, "Sqrt(1e401) has the wrong exponent"); + } + + [TestMethod] + public void TestSqrtSquaresBackToItsInput() + { + string[] inputs = ["2", "3", "7", "0.5", "1E-17", "123456.789", "9.87654321E31", "1E400", "6.02214076E-23"]; + + foreach (string input in inputs) + { + PreciseNumber value = Parse(input); + PreciseNumber root = PreciseNumber.Sqrt(value, 50); + + AssertAgreesTo(value, root.Squared(), 48, $"Sqrt({input}) does not square back"); + } + } + + [TestMethod] + public void TestSqrtOfZeroAndOne() + { + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Sqrt(PreciseNumber.Zero)); + Assert.AreEqual(PreciseNumber.One, PreciseNumber.Sqrt(PreciseNumber.One)); + } + + [TestMethod] + public void TestSqrtOfANegativeValueThrows() + { + // There is no NaN to return, so this is a genuine divergence from double and the exception + // type is part of the contract. + Assert.ThrowsExactly(() => PreciseNumber.Sqrt((-1).ToPreciseNumber())); + Assert.ThrowsExactly(() => PreciseNumber.Sqrt((-1).ToPreciseNumber(), 50)); + } + + [TestMethod] + public void TestSqrtDefaultsToTheInputsPrecisionRatherThanCappingIt() + { + // The rule Divide follows: never fewer digits than the operand, and never fewer than + // MinimumDivisionPrecision. Asserted as a floor because a root whose last digit rounds to + // zero has it stripped, as everywhere else in the library: Sqrt(Pi) is one such. + Assert.AreEqual(PreciseNumber.MinimumDivisionPrecision, PreciseNumber.Sqrt(2.ToPreciseNumber()).SignificantDigits); + + Assert.IsGreaterThanOrEqualTo( + PreciseNumber.ConstantPrecision - 1, + PreciseNumber.Sqrt(PreciseNumber.Pi).SignificantDigits, + "Sqrt(Pi) was capped at MinimumDivisionPrecision rather than following its operand"); + } + + [TestMethod] + public void TestCbrtMatchesPublishedDigits() => + Assert.AreEqual(Cbrt2Digits, Digits(PreciseNumber.Cbrt(2.ToPreciseNumber(), 50)), "Cbrt(2) is wrong"); + + [TestMethod] + public void TestCbrtOfANegativeValueIsReal() + { + Assert.AreEqual((-2).ToPreciseNumber(), PreciseNumber.Cbrt((-8).ToPreciseNumber()), "Cbrt(-8) is not -2"); + Assert.AreEqual($"-{Cbrt2Digits}", Digits(PreciseNumber.Cbrt((-2).ToPreciseNumber(), 50)), "Cbrt(-2) has the wrong digits"); + Assert.AreEqual(PreciseNumber.Cbrt(2.ToPreciseNumber(), 50), -PreciseNumber.Cbrt((-2).ToPreciseNumber(), 50), "Cbrt is not odd about zero"); + } + + [TestMethod] + public void TestCbrtOfAPerfectCubeIsExact() + { + Assert.AreEqual(7.ToPreciseNumber(), PreciseNumber.Cbrt(343.ToPreciseNumber()), "Cbrt(343) is not 7"); + Assert.AreEqual(Parse("0.2"), PreciseNumber.Cbrt(Parse("0.008")), "Cbrt(0.008) is not 0.2"); + } + + [TestMethod] + public void TestRootNMatchesPublishedDigits() => + Assert.AreEqual(Root5Of7Digits, Digits(PreciseNumber.RootN(7.ToPreciseNumber(), 5, 50)), "The fifth root of 7 is wrong"); + + [TestMethod] + public void TestRootNRaisesBackToItsInput() + { + int[] degrees = [2, 3, 4, 5, 9, 17]; + + foreach (int degree in degrees) + { + PreciseNumber value = Parse("1234.5678"); + PreciseNumber root = PreciseNumber.RootN(value, degree, 50); + + AssertAgreesTo(value, root.Pow(degree.ToPreciseNumber()), 45, $"The root of degree {degree} does not raise back"); + } + } + + [TestMethod] + public void TestRootNOfDegreeOneIsTheValueItself() => + Assert.AreEqual(Parse("1234.5678"), PreciseNumber.RootN(Parse("1234.5678"), 1)); + + [TestMethod] + public void TestRootNOfANegativeDegreeIsTheReciprocal() + { + PreciseNumber reciprocal = PreciseNumber.RootN(16.ToPreciseNumber(), -2); + + Assert.AreEqual(Parse("0.25"), reciprocal, "The -2 root of 16 is not 1/4"); + AssertAgreesTo(PreciseNumber.One / PreciseNumber.Sqrt(2.ToPreciseNumber()), PreciseNumber.RootN(2.ToPreciseNumber(), -2), 48, "The -2 root of 2 is not 1/Sqrt(2)"); + } + + [TestMethod] + public void TestAnEvenRootOfANegativeValueThrows() + { + Assert.ThrowsExactly(() => PreciseNumber.RootN((-16).ToPreciseNumber(), 4)); + Assert.ThrowsExactly(() => PreciseNumber.RootN((-16).ToPreciseNumber(), -4)); + } + + [TestMethod] + public void TestAnOddRootOfANegativeValueCarriesTheSign() => + Assert.AreEqual((-3).ToPreciseNumber(), PreciseNumber.RootN((-243).ToPreciseNumber(), 5), "The fifth root of -243 is not -3"); + + [TestMethod] + public void TestARootOfDegreeZeroThrows() + { + Assert.ThrowsExactly(() => PreciseNumber.RootN(2.ToPreciseNumber(), 0)); + Assert.ThrowsExactly(() => PreciseNumber.RootN(2.ToPreciseNumber(), int.MinValue)); + } + + [TestMethod] + public void TestARootOfZeroIsZeroUnlessTheDegreeIsNegative() + { + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.RootN(PreciseNumber.Zero, 5)); + Assert.ThrowsExactly(() => PreciseNumber.RootN(PreciseNumber.Zero, -5)); + } + + [TestMethod] + public void TestARootOfFewerThanOneDigitThrows() + { + Assert.ThrowsExactly(() => PreciseNumber.Sqrt(2.ToPreciseNumber(), 0)); + Assert.ThrowsExactly(() => PreciseNumber.Cbrt(2.ToPreciseNumber(), -1)); + Assert.ThrowsExactly(() => PreciseNumber.RootN(2.ToPreciseNumber(), 3, 0)); + Assert.ThrowsExactly(() => PreciseNumber.Hypot(3.ToPreciseNumber(), 4.ToPreciseNumber(), 0)); + } + + [TestMethod] + public void TestHypotOfAPythagoreanTripleIsExact() + { + Assert.AreEqual(5.ToPreciseNumber(), PreciseNumber.Hypot(3.ToPreciseNumber(), 4.ToPreciseNumber()), "Hypot(3, 4) is not 5"); + Assert.AreEqual(13.ToPreciseNumber(), PreciseNumber.Hypot(5.ToPreciseNumber(), 12.ToPreciseNumber()), "Hypot(5, 12) is not 13"); + } + + [TestMethod] + public void TestHypotMatchesTheSquareRootOfTheSumOfSquares() + { + Assert.AreEqual(Sqrt2Digits, Digits(PreciseNumber.Hypot(PreciseNumber.One, PreciseNumber.One, 50)), "Hypot(1, 1) is not Sqrt(2)"); + Assert.AreEqual(5.ToPreciseNumber(), PreciseNumber.Hypot((-3).ToPreciseNumber(), (-4).ToPreciseNumber()), "Hypot does not ignore the signs"); + } + + [TestMethod] + public void TestHypotNeedsNoScalingToAvoidOverflow() + { + // The magnitudes that make a double implementation overflow while squaring. There is no + // exponent range to fall out of here, so the answer is the one the algebra gives. + PreciseNumber large = PreciseNumber.Hypot(Parse("3E400"), Parse("4E400")); + PreciseNumber small = PreciseNumber.Hypot(Parse("3E-400"), Parse("4E-400")); + + Assert.AreEqual(Parse("5E400"), large, "Hypot overflowed at 1e400"); + Assert.AreEqual(Parse("5E-400"), small, "Hypot underflowed at 1e-400"); + } + + [TestMethod] + public void TestRootsAreReachableThroughTheInterface() + { + // The point of implementing IRootFunctions is that generic code can ask for a root without + // knowing which number it holds. + Assert.AreEqual(3.ToPreciseNumber(), RootThroughInterface(9.ToPreciseNumber()), "Sqrt is not reachable generically"); + Assert.AreEqual(PreciseNumber.Pi, PiThroughInterface(), "IFloatingPointConstants is not satisfied"); + } + + private static TNumber RootThroughInterface(TNumber value) + where TNumber : IRootFunctions => + TNumber.Sqrt(value); + + private static TNumber PiThroughInterface() + where TNumber : IRootFunctions => + TNumber.Pi; + + [TestMethod] + public void TestRootsOfLongSignificandsStayCorrect() + { + // 200 digits is the widest of the repository's benchmark digit counts, and wide enough that + // a working precision chosen for the requested digits alone would show up here. + BigInteger significand = BigInteger.Parse(new string('7', 200), NumberStyles.None, CultureInfo.InvariantCulture); + PreciseNumber value = PreciseNumber.CreateFromComponents(-100, significand); + + PreciseNumber root = PreciseNumber.Sqrt(value); + + Assert.AreEqual(200, root.SignificantDigits, "A 200 digit input produced a shorter root"); + AssertAgreesTo(value, root.Squared(), 195, "The root of a 200 digit value does not square back"); + } +} diff --git a/PreciseNumber/PreciseNumber.Roots.cs b/PreciseNumber/PreciseNumber.Roots.cs new file mode 100644 index 0000000..aeb897c --- /dev/null +++ b/PreciseNumber/PreciseNumber.Roots.cs @@ -0,0 +1,297 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber; + +using System; +using System.Numerics; + +/// +/// Square, cube and n-th roots, and the hypotenuse built on them. +/// +/// +/// Every root is taken on the significand as an integer rather than on the number as a whole. A +/// value is significand × 10^exponent, so scaling the significand by a power of ten until the +/// exponent divides by the degree splits the problem in two: the exponent is rooted by division +/// alone, and the significand is rooted by integer Newton, which terminates on an exact answer +/// instead of on a tolerance. A value whose root is exact therefore comes back exact, however many +/// digits were asked for, the same way is exact +/// when a quotient terminates. +/// +public readonly partial record struct PreciseNumber + : IRootFunctions +{ + /// + /// Digits computed past the ones the caller asked for, so that the digit the final rounding + /// decision is made on is itself correct. + /// + private const int RootGuardDigits = 2; + + /// + /// Iterations allowed to a root beyond its degree before it is called non-convergent. + /// + /// + /// The descent below starts at most twice the root and is strictly decreasing, so it cannot + /// loop. The allowance grows with the degree because the first phase closes a fixed fraction of + /// the gap per step, and that fraction is 1/n; measured over degrees up to 2000 the + /// worst case needs four iterations more than the degree. + /// + private const int RootIterationAllowance = 64; + + /// + /// The message carried by the thrown for an even root + /// of a negative value. + /// + private const string NegativeRootMessage = "A negative value has no real root of an even degree."; + + /// + /// Returns the square root of a value. + /// + /// The value to take the square root of. + /// The positive square root of . + /// Thrown when is negative. + /// + /// A perfect square roots exactly. Anything else is produced to the significant digits of + /// , and never fewer than , matching + /// . Use + /// to choose that precision. + /// + /// has no NaN, so a negative value throws where a + /// would quietly return NaN and carry on. An algorithm ported from + /// that relies on that has to test the sign itself. + /// + /// + public static PreciseNumber Sqrt(PreciseNumber x) => + Sqrt(x, DefaultRootPrecision(x)); + + /// + /// Returns the square root of a value, to a chosen number of significant digits. + /// + /// The value to take the square root of. + /// + /// The number of significant digits to produce. A root that is exact is exact regardless of this + /// value. + /// + /// The positive square root of . + /// + /// Thrown when is negative, or when is + /// less than one. + /// + public static PreciseNumber Sqrt(PreciseNumber x, int significantDigits) => + RootN(x, 2, significantDigits); + + /// + /// Returns the cube root of a value. + /// + /// The value to take the cube root of. + /// The cube root of , which carries its sign. + /// + /// A perfect cube roots exactly. Anything else is produced to the significant digits of + /// , and never fewer than . The cube + /// root of a negative value is real, so it is returned rather than rejected. + /// + public static PreciseNumber Cbrt(PreciseNumber x) => + RootN(x, 3); + + /// + /// Returns the cube root of a value, to a chosen number of significant digits. + /// + /// The value to take the cube root of. + /// + /// The number of significant digits to produce. A root that is exact is exact regardless of this + /// value. + /// + /// The cube root of , which carries its sign. + /// Thrown when is less than one. + public static PreciseNumber Cbrt(PreciseNumber x, int significantDigits) => + RootN(x, 3, significantDigits); + + /// + /// Returns the n-th root of a value. + /// + /// The value to take the root of. + /// The degree of the root. + /// The n-th root of . + /// + /// Thrown when is zero or , or when + /// is negative and is even. + /// + /// Thrown when is zero and is negative. + /// + /// A root that is exact is produced exactly. Anything else is produced to the significant digits + /// of , and never fewer than . + /// + public static PreciseNumber RootN(PreciseNumber x, int n) => + RootN(x, n, DefaultRootPrecision(x)); + + /// + /// Returns the n-th root of a value, to a chosen number of significant digits. + /// + /// The value to take the root of. + /// The degree of the root. + /// + /// The number of significant digits to produce. A root that is exact is exact regardless of this + /// value. + /// + /// The n-th root of . + /// + /// Thrown when is zero or , when + /// is negative and is even, or when + /// is less than one. + /// + /// Thrown when is zero and is negative. + /// Thrown when the root does not converge. + public static PreciseNumber RootN(PreciseNumber x, int n, int significantDigits) + { + if (significantDigits < 1) + { + throw new ArgumentOutOfRangeException(nameof(significantDigits), significantDigits, "At least one significant digit is required."); + } + + // A degree of int.MinValue has no negation that fits an int, so the reciprocal below cannot + // express it. Zero asks for a root that is not a number at all. + if (n is 0 or int.MinValue) + { + throw new ArgumentOutOfRangeException(nameof(n), n, "The degree of a root must be a non-zero value other than int.MinValue."); + } + + if (n < 0) + { + return x.Significand.IsZero + ? throw new DivideByZeroException() + : Divide(One, RootN(x, -n, significantDigits), significantDigits); + } + + if (x.Significand.IsZero || n == 1) + { + return x; + } + + if (x.Significand.Sign > 0) + { + return PositiveRootN(x, n, significantDigits); + } + + // An odd root of a negative value is real, so the sign comes out and goes back on. + return int.IsEvenInteger(n) + ? throw new ArgumentOutOfRangeException(nameof(x), x, NegativeRootMessage) + : -PositiveRootN(-x, n, significantDigits); + } + + /// + /// Returns the length of the hypotenuse of a right triangle with the given side lengths. + /// + /// The length of one side. + /// The length of the other side. + /// The square root of x² + y². + /// + /// Computed directly, without the scaling a implementation needs. That + /// scaling exists to keep x² inside a fixed exponent range, and a + /// has no such range: squaring and adding are both exact, so the + /// value handed to the square root is the exact sum whatever the magnitudes involved. + /// + public static PreciseNumber Hypot(PreciseNumber x, PreciseNumber y) => + Hypot(x, y, Math.Max(DefaultRootPrecision(x), DefaultRootPrecision(y))); + + /// + /// Returns the length of the hypotenuse of a right triangle with the given side lengths, to a + /// chosen number of significant digits. + /// + /// The length of one side. + /// The length of the other side. + /// + /// The number of significant digits to produce. A result that is exact is exact regardless of + /// this value. + /// + /// The square root of x² + y². + /// Thrown when is less than one. + public static PreciseNumber Hypot(PreciseNumber x, PreciseNumber y, int significantDigits) => + Sqrt(Add(Multiply(x, x), Multiply(y, y)), significantDigits); + + /// + /// Gets the significant digits a root produces when the caller does not choose. + /// + /// The value being rooted. + /// The significant digits of , or if that is more. + /// + /// The same rule follows, so that rooting a + /// constant carrying digits does not silently cap the expression + /// at fifty. + /// + private static int DefaultRootPrecision(PreciseNumber value) => + Math.Max(value.SignificantDigits, MinimumDivisionPrecision); + + /// + /// Computes 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 two. + /// The number of significant digits to produce. + /// The n-th root of . + /// Thrown when the root needs a scale or an exponent wider than an . + private static PreciseNumber PositiveRootN(PreciseNumber value, int n, int significantDigits) + { + // root(s · 10^(e - k)) == root(s · 10^k') · 10^((e - k) / n) once n divides e - k, so + // scaling the significand until that holds leaves an integer root and an exact exponent. + int digits = CountDigits(value.Significand); + long wanted = (long)n * (significantDigits + RootGuardDigits); + long scale = Math.Max(wanted - digits, 0); + + // Raise the scale to the next value that leaves an exponent the degree divides. + scale += (((value.Exponent - scale) % n) + n) % n; + + long rootExponent = (value.Exponent - scale) / 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 scaled = scale > 0 ? value.Significand * Pow10((int)scale) : value.Significand; + BigInteger root = IntegerRootN(scaled, n); + PreciseNumber result = new((int)rootExponent, root); + + // An exact root has every digit of its own, so rounding it to the requested precision would + // throw away an answer that is already right. + return BigInteger.Pow(root, n) == scaled + ? result + : result.ReduceSignificance(significantDigits); + } + + /// + /// Computes the largest integer whose n-th power does not exceed a value. + /// + /// The value to take the root of, which must be positive. + /// The degree of the root, which must be at least two. + /// The floor of the n-th root of . + /// Thrown when the iteration does not converge. + /// + /// Newton on integers, x ← ((n-1)·x + v / x^(n-1)) / n. Started above the root it is + /// strictly decreasing and stops at the floor of the root exactly, so there is no tolerance to + /// pick and no working precision to carry: every digit it returns is correct. + /// + private static BigInteger IntegerRootN(BigInteger value, int n) + { + if (value.IsZero || value.IsOne) + { + return value; + } + + // value < 2^bits, so 2^ceil(bits / n) is above the root and within a factor of two of it. + int power = (int)((value.GetBitLength() + n - 1) / n); + BigInteger estimate = BigInteger.One << power; + + for (int iteration = 0; iteration <= RootIterationAllowance + n; iteration++) + { + BigInteger next = (((n - 1) * estimate) + (value / BigInteger.Pow(estimate, n - 1))) / n; + if (next >= estimate) + { + return estimate; + } + + estimate = next; + } + + throw new ArithmeticException( + $"The root of degree {n.ToString(InvariantCulture)} of a {CountDigits(value).ToString(InvariantCulture)} digit value did not converge."); + } +} diff --git a/README.md b/README.md index ad8d07c..4bb1a83 100644 --- a/README.md +++ b/README.md @@ -82,7 +82,7 @@ A high-precision numeric type for .NET that provides arbitrary precision arithme - **Value Type**: A `readonly record struct` whose `default` value is zero. Adding, subtracting, multiplying, and comparing allocate nothing when the operands and every intermediate and final significand fit in an `int`. Exponent alignment counts, so `1 + 0.0000000001` allocates because it scales 1 by 10^10, and `99999 * 99999` allocates because its product is 9,999,800,001. -- **Comprehensive Mathematical Support**: Includes advanced mathematical functions like exponential operations (Pow, Exp, Squared, Cubed), constant values (Pi, E, Tau) with high precision, absolute value operations, and specialized numerical checks (isOdd, isEven, etc.)—all with arbitrary precision. +- **Comprehensive Mathematical Support**: Includes advanced mathematical functions like exponential operations (Pow, Exp, Squared, Cubed), roots (Sqrt, Cbrt, RootN, Hypot) through `IRootFunctions`, constant values (Pi, E, Tau) with high precision, absolute value operations, and specialized numerical checks (isOdd, isEven, etc.)—all with arbitrary precision. - **Balanced Performance**: The design prioritizes accuracy and precision while maintaining reasonable performance. For calculations where extreme precision matters more than raw speed, PreciseNumber delivers excellent results, though built-in numeric types remain faster for standard precision needs. @@ -278,6 +278,16 @@ var e = PreciseNumber.E; // Exponential function var expValue = PreciseNumber.Exp(1.ToPreciseNumber()); // e^1 = e +// Roots. A value whose root is exact gets it exactly, whatever precision was asked for +var root = PreciseNumber.Sqrt(2.ToPreciseNumber()); // 1.4142135623730950488016887242096980785696718753769 +var exactRoot = PreciseNumber.Sqrt(144.ToPreciseNumber()); // 12 +var cubeRoot = PreciseNumber.Cbrt((-8).ToPreciseNumber()); // -2, an odd root of a negative value being real +var fifthRoot = PreciseNumber.RootN(7.ToPreciseNumber(), 5); +var hypotenuse = PreciseNumber.Hypot(3.ToPreciseNumber(), 4.ToPreciseNumber()); // 5 + +// Or choose the precision, the same way Divide does +var shortRoot = PreciseNumber.Sqrt(2.ToPreciseNumber(), 10); // 1.414213562 + // Rounding and precision control var roundedValue = number.Round(1); // 2.5 (already at 1 decimal place) var reducedValue = number.ReduceSignificance(1); // 3 (reduced to 1 significant digit) @@ -415,6 +425,8 @@ You can control precision using: - **Divide(left, right, significantDigits)**: Chooses the precision of a quotient +- **Sqrt(value, significantDigits)**, and the same overload on `Cbrt`, `RootN` and `Hypot`: Chooses the precision of a root + Division produces a terminating quotient exactly, however many digits that takes — `1 / 8` is `0.125`, and `1 / 2^64` keeps all 64 decimal places. A repeating quotient is produced to the precision of the wider operand, never fewer than `MinimumDivisionPrecision` (50) significant @@ -429,9 +441,17 @@ Multiplication is exact, so any product involving one of them carries at least 1 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. +`Sqrt`, `Cbrt`, `RootN` and `Hypot` follow the same rule as division: a value whose root is exact +gets that root exactly whatever precision was asked for, and anything else is produced to the +significant digits of the operand, never fewer than `MinimumDivisionPrecision`. None of them goes +through `double`, so a value outside its range — `1e400`, or `1e-400` — roots as accurately as any +other. + ## 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 +- `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, and `RootN()` of a negative value at an even degree, throw `ArgumentOutOfRangeException` where a `double` would return NaN and carry on - 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 @@ -451,6 +471,8 @@ round half away from zero and cache per requested precision. - **Functions**: `Abs()`, `Round()`, `Clamp()`, `Squared()`, `Cubed()`, `Pow()`, `Exp()` +- **Roots**: `Sqrt()`, `Cbrt()`, `RootN()`, `Hypot()`, each with an overload taking the significant digits to produce + - **Utility**: `ToString()`, `Parse()`, `TryParse()`, `To()` - **Generic Conversion**: `TryConvertFromChecked`, `TryConvertFromSaturating`, `TryConvertFromTruncating`, `TryConvertToChecked`, `TryConvertToSaturating`, and `TryConvertToTruncating`, reached through `CreateChecked`, `CreateSaturating`, and `CreateTruncating`