From aa022fd4f951500595107caec3133d213ef67e1c Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 16 Sep 2026 13:34:43 +0000 Subject: [PATCH 1/2] Make StorageMath public, and add Cbrt, RootN and Hypot [minor] MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit StorageMath.Sqrt is the piece of this library that System.Numerics does not offer: a root at the precision of an arbitrary INumber, with no IRootFunctions constraint, which decimal could not have satisfied anyway. It was internal, so an application doing its own vector math over quantities — a norm the generator does not emit — had to reimplement it. The class and Sqrt are now public, with what they guarantee written down rather than left to be discovered: the 256-step cap and the ArithmeticException on a root that does not settle, the double round trip for the primitives, the floor for integers, and OverflowException rather than NaN for a negative value on a type with no NaN. The seeding, the double round trip and the Newton loop stay private. Cbrt, RootN and Hypot ship alongside, on the same seeding and the same loop, since they are what a caller asks for next. Sqrt itself is untouched, so every existing result is unchanged. Cbrt and RootN take the double route only for the binary floating point primitives and refine every integer type in integer arithmetic, so their floor is exact rather than whatever Math.Pow rounded to; each power is taken checked, so a narrow type reports the overflow instead of wrapping into an estimate that goes nowhere. Hypot computes a fractional type from the ratio of its legs, so a pair whose squares leave the type still has its hypotenuse. StorageMathTests pins both halves of the surface — the four methods public, the workings not — so the reason this is public does not erode. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01A9Ke7x6EeGyKHfPHncbNrP --- CLAUDE.md | 12 + Semantics.Quantities/README.md | 1 + Semantics.Quantities/StorageMath.cs | 456 +++++++++++++++++- Semantics.Test/Quantities/StorageMathTests.cs | 214 +++++++- docs/physics-generator.md | 4 + 5 files changed, 680 insertions(+), 7 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 0d53df7a..318ab867 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -305,6 +305,18 @@ back by the matching power of two. A root that does not settle throws `Arithmeti than returning an estimate. The logarithmic scales and the hand-written audio types still compute through `double`. +`StorageMath` is public API, not just the generator's helper (#239): an application doing its own +vector math over quantities would otherwise reimplement the root, and worse. `Cbrt`, `RootN` and +`Hypot` ship alongside `Sqrt` on the same seeding and the same Newton loop. `Sqrt` keeps its round +trip through `double` for every primitive, integers included, because the generated code always +inlined it; `Cbrt` and `RootN` take that route only for the binary floating point primitives and +refine every integer type in integer arithmetic, so their floor is exact rather than whatever +`Math.Pow` rounded to. `Hypot` computes a fractional type from the ratio of its legs, so a pair whose +squares leave the type still has its hypotenuse, and squares an integer type directly, since the +ratio of two integers is not a ratio. The seeding, the double round trip and the Newton loop itself +stay private — `StorageMathTests.TheRootsArePublicAndTheirWorkingsAreNot` pins both halves of that, +since this is a package with a compatibility baseline and the shape is frozen once it ships. + `StorageConversionTests` runs the same conversions, relationships and vector lengths over `double` and `decimal`, exactly where the answer terminates and to a relative tolerance where it does not. Adding a storage type is one derived class. diff --git a/Semantics.Quantities/README.md b/Semantics.Quantities/README.md index fe9575c1..b2c8b103 100644 --- a/Semantics.Quantities/README.md +++ b/Semantics.Quantities/README.md @@ -124,6 +124,7 @@ ForceMagnitude weight = Mass.FromKilogram(70.0) * AccelerationMa | `IVector2` / `IVector3` / `IVector4` | Directional vectors with `X`/`Y`/`Z`/`W`, `Length()`, `LengthSquared()`, `Dot`, `Distance`, `Normalize`; `IVector3` adds `Cross`. | | `Vector0Guards` | `EnsureNonNegative(value, name)` and `EnsurePositive(value, name)`, used by generated `From{Unit}` factories. | | `UnitSystem` | enum classifying units (`SIBase`, `SIDerived`, `Metric`, `Imperial`, ...). | +| `StorageMath` | `Sqrt`, `Cbrt`, `RootN` and `Hypot` over any `T : struct, INumber`, with no `IRootFunctions` constraint, which `decimal` could not meet. Generated `Length()` and `Distance()` use them; an application computing a norm the generator does not emit can too. | ### Generated quantity types diff --git a/Semantics.Quantities/StorageMath.cs b/Semantics.Quantities/StorageMath.cs index 561b8ae6..73adc2b2 100644 --- a/Semantics.Quantities/StorageMath.cs +++ b/Semantics.Quantities/StorageMath.cs @@ -6,9 +6,50 @@ namespace ktsu.Semantics.Quantities; using System.Numerics; /// -/// Numeric operations the generated quantities need that does not declare. +/// Roots at the precision of an arbitrary , which +/// declares only for the types that can satisfy +/// . /// -internal static class StorageMath +/// +/// +/// The generated quantities call these for Length() and Distance(). They are public +/// because an application doing its own vector math over quantities — a physics integrator, an orbit +/// propagator, anything computing a norm the library does not already emit — otherwise has to +/// reimplement them, and no constraint beyond is needed, which is the +/// whole trick: could not meet anyway. +/// +/// +/// What every method here guarantees, so that none of it has to be discovered: +/// +/// +/// +/// A root is refined by Newton steps in T's own arithmetic, from a seed within +/// a factor of two of the answer. At most 256 steps are taken; an estimate still moving after that is +/// cycling rather than settling, and is reported as rather than +/// returned as an answer. +/// +/// +/// The binary floating point primitives round trip through . So +/// Sqrt<double> is exactly and Cbrt<double> +/// exactly , not a refinement of either: a caller expecting digits +/// beyond what holds will not get them from a quantity. +/// +/// +/// An integer type returns the floor of the root. reaches it through +/// for the integer primitives, which is the route the generated code always +/// inlined; and refine every integer type in +/// integer arithmetic, so their floor is exact rather than whatever +/// happened to round to. +/// +/// +/// A root with no real answer — an even root of a negative value — behaves as the +/// round trip always did: where +/// T has one, and where it does not, as +/// does not. +/// +/// +/// +public static class StorageMath { /// /// The most Newton steps taken before the root is reported as not converging. @@ -52,7 +93,7 @@ internal static class StorageMath /// one and an exception where it does not. /// /// - internal static T Sqrt(T value) + public static T Sqrt(T value) where T : struct, INumber { if (IsRoundedThroughDouble() || T.IsNegative(value)) @@ -91,6 +132,141 @@ internal static T Sqrt(T value) throw new ArithmeticException($"The square root did not settle within {MaximumIterations} Newton steps."); } + /// + /// Computes the cube root of a storage value at the precision of its type. + /// + /// The numeric storage type. + /// The value to take the root of. A negative value has a real cube root and gets it. + /// The cube root of , or its floor towards zero for an integer type. + /// + /// The Newton steps do not settle on a root within 256 iterations. + /// + /// + /// The binary floating point primitives take . Every other type, + /// integers included, is refined in its own arithmetic, so a cube root has the + /// digits the type holds and an integer one is the exact floor. + /// + public static T Cbrt(T value) + where T : struct, INumber + => RootN(value, 3); + + /// + /// Computes the th root of a storage value at the precision of its type. + /// + /// The numeric storage type. + /// The value to take the root of. + /// The degree of the root. An odd degree accepts a negative ; an even one does not. + /// + /// The th root of , or its floor towards zero for an + /// integer type. + /// + /// is not positive. + /// + /// is even and is negative, and + /// cannot represent the that results. + /// + /// + /// The Newton steps do not settle on a root within 256 iterations, which a type too narrow to hold + /// the intermediate estimate^(n-1) can cause. + /// + /// + /// + /// A degree of one returns and a degree of two is + /// , contract and all. Any other degree is seeded through + /// and refined by the Newton step + /// x = (((n - 1) * x) + (value / x^(n - 1))) / n in 's arithmetic. + /// + /// + /// The seed for a value outside the range of is taken by scaling by powers of + /// 2^n into [1, 2^n), rooting there, and scaling the root back by the matching power of + /// two — the same trick plays with powers of four, and for the same reason: + /// without it a of 2^2048 spends every step on halving. + /// + /// + public static T RootN(T value, int n) + where T : struct, INumber + { + if (n < 1) + { + throw new ArgumentOutOfRangeException(nameof(n), n, "The degree of a root is a positive integer."); + } + + if (n == 1) + { + return value; + } + + if (n == 2) + { + return Sqrt(value); + } + + if (T.IsNegative(value)) + { + // An even root of a negative value has no real answer, so it keeps the behaviour of the + // double round trip: NaN where the type has one, and an exception where it does not. + return int.IsEvenInteger(n) + ? T.CreateChecked(Math.Pow(double.CreateChecked(value), 1d / n)) + : -RootN(-value, n); + } + + if (T.IsZero(value) || value == T.One) + { + return value; + } + + return IsBinaryFloatingPoint() + ? T.CreateChecked(RootThroughDouble(double.CreateChecked(value), n)) + : RootByNewton(value, n); + } + + /// + /// Computes the length of the hypotenuse of a right triangle at the precision of the storage type. + /// + /// The numeric storage type. + /// One leg of the triangle. + /// The other leg of the triangle. + /// The square root of x² + y², or its floor for an integer type. + /// + /// The Newton steps do not settle on a root within 256 iterations. + /// + /// + /// + /// A type with a fractional part is computed as larger * Sqrt(1 + (smaller / larger)²), so a + /// pair whose squares would leave the range of the type still has its hypotenuse. An integer type + /// squares and sums directly, because the ratio of two integers is not a ratio, and the sum is + /// therefore bounded by the type exactly as the caller's own x * x + y * y would be. + /// + /// + /// The primitives take the round trip, so Hypot<double> is exactly + /// . + /// + /// + public static T Hypot(T x, T y) + where T : struct, INumber + { + if (IsRoundedThroughDouble()) + { + return T.CreateChecked(double.Hypot(double.CreateChecked(x), double.CreateChecked(y))); + } + + T larger = T.Max(T.Abs(x), T.Abs(y)); + T smaller = T.Min(T.Abs(x), T.Abs(y)); + + if (T.IsZero(larger)) + { + return T.Zero; + } + + if (HasFloorDivision()) + { + return Sqrt((x * x) + (y * y)); + } + + T ratio = smaller / larger; + return larger * Sqrt(T.One + (ratio * ratio)); + } + /// /// Reports whether is one of the primitives whose square root has always /// been taken through . @@ -180,7 +356,7 @@ private static T Seed(T value, T two) } /// - /// Takes the root of a value through when the value converts to a normal + /// Takes the square root of a value through when the value converts to a normal /// and the root converts back to something other than zero. /// /// The numeric storage type. @@ -189,13 +365,26 @@ private static T Seed(T value, T two) /// when holds a usable estimate. private static bool TryRootThroughDouble(T value, out T root) where T : struct, INumber + => TryRootThroughDouble(value, 2, out root); + + /// + /// Takes the th root of a value through when the value + /// converts to a normal and the root converts back to something other than zero. + /// + /// The numeric storage type. + /// The positive value. + /// The degree of the root. + /// The root in , or the default when this fails. + /// when holds a usable estimate. + private static bool TryRootThroughDouble(T value, int degree, out T root) + where T : struct, INumber { try { double asDouble = double.CreateChecked(value); if (double.IsNormal(asDouble)) { - root = T.CreateChecked(Math.Sqrt(asDouble)); + root = T.CreateChecked(RootThroughDouble(asDouble, degree)); if (!T.IsZero(root)) { return true; @@ -214,4 +403,261 @@ private static bool TryRootThroughDouble(T value, out T root) root = default; return false; } + + /// + /// Takes a root in , by the most accurate route the framework offers for its degree. + /// + /// The value to take the root of. + /// The degree of the root. + /// The root of . + /// + /// and are correctly rounded where + /// is not, and answers + /// for a negative value however odd the degree, so the sign is taken out + /// first. + /// + private static double RootThroughDouble(double value, int degree) + { + if (degree == 2) + { + return Math.Sqrt(value); + } + + if (degree == 3) + { + return Math.Cbrt(value); + } + + return double.IsNegative(value) && int.IsOddInteger(degree) + ? -Math.Pow(-value, 1d / degree) + : Math.Pow(value, 1d / degree); + } + + /// + /// Reports whether is one of the binary floating point primitives, whose + /// root is already as precise as the type once it has been taken in . + /// + /// The numeric storage type. + /// for , and . + private static bool IsBinaryFloatingPoint() + where T : struct, INumber + => typeof(T) == typeof(double) + || typeof(T) == typeof(float) + || typeof(T) == typeof(Half); + + /// + /// Reports whether division in discards the fractional part, which is what + /// makes a root in it the floor of the root rather than the root. + /// + /// The numeric storage type. + /// when one divided by two is zero. + /// + /// Asked of the arithmetic rather than of a list of types, so a numeric type written outside this + /// library is classified by how it behaves. + /// + private static bool HasFloorDivision() + where T : struct, INumber + => T.IsZero(T.One / (T.One + T.One)); + + /// + /// Refines the th root of a positive value by Newton steps in + /// 's own arithmetic. + /// + /// The numeric storage type. + /// The value to take the root of, greater than zero and not one. + /// The degree of the root, three or more. + /// The root, or its floor for a type whose division floors. + /// The estimate does not settle within steps. + /// + /// + /// A type whose division floors descends: from an estimate at or above the root every step lands no + /// lower than the floor of the root, so the first step that does not descend has passed it and the + /// estimate before it is the answer. The seed is lifted to at or above the root first, since a + /// descent that starts below it would stop on the first step and answer with the seed. + /// + /// + /// Any other type settles instead, on the same terms as : an estimate that + /// stops changing is the root, and one alternating between two neighbours is rounding, so the + /// smaller of the two is taken. + /// + /// + /// Every power is taken checked, because a narrow type can hold a value whose estimate^(n-1) + /// it cannot, and an estimate that wraps sends the steps somewhere arbitrary. A degree so high that + /// the type cannot hold 2^degree is answered before any of that: no value the type holds has + /// a root of two or more, so the answer is one. + /// + /// + private static T RootByNewton(T value, int degree) + where T : struct, INumber + { + T two = T.One + T.One; + T count = T.CreateChecked(degree); + T countLessOne = count - T.One; + bool floors = HasFloorDivision(); + T estimate = SeedForRoot(value, degree, two); + + if (floors) + { + if (!TryPower(two, degree, out _)) + { + // Two raised to this degree leaves the type, so no value the type holds has a root of two + // or more, and the values with a root below one were answered before this was called. + return T.One; + } + + // An estimate whose own power overflows is above the root, since the root's power divides + // a value the type holds, so the lift stops there as well as on a power that exceeds the value. + while (TryPower(estimate, degree, out T raised) && raised < value) + { + estimate *= two; + } + } + + T previous = estimate; + + for (int iteration = 0; iteration < MaximumIterations; iteration++) + { + if (!TryPower(estimate, degree - 1, out T power) || T.IsZero(power)) + { + // The estimate is too large for its own power to be taken in the type. Halving reaches a + // range the type holds, and for a type that floors it stays at or above the root. + estimate = (estimate + T.One) / two; + continue; + } + + T next = ((countLessOne * estimate) + (value / power)) / count; + + if (floors) + { + if (next >= estimate) + { + return estimate; + } + + estimate = next; + continue; + } + + if (next == estimate) + { + return next; + } + + // Rounding can leave the estimate alternating between two neighbours instead of settling. + if (next == previous) + { + return T.Min(estimate, next); + } + + previous = estimate; + estimate = next; + } + + throw new ArithmeticException($"The root of degree {degree} did not settle within {MaximumIterations} Newton steps."); + } + + /// + /// Chooses the first Newton estimate for the th root of a positive value. + /// + /// The numeric storage type. + /// The positive value whose root is wanted. + /// The degree of the root. + /// Two, in . + /// An estimate within a factor of two of the root. + /// + /// The root is used directly when converts to a normal + /// . Otherwise the value is scaled by powers of 2^degree into + /// [1, 2^degree), the root is taken there, and that root is scaled back by the matching power + /// of two — 's trick for a degree other than two. + /// + private static T SeedForRoot(T value, int degree, T two) + where T : struct, INumber + { + if (TryRootThroughDouble(value, degree, out T direct)) + { + return direct; + } + + if (!TryPower(two, degree, out T scale)) + { + // The type cannot hold 2^degree, so it cannot hold a value needing this scaling either. + return T.One; + } + + T scaled = value; + int powerOfTwo = 0; + + while (scaled >= scale) + { + scaled /= scale; + powerOfTwo++; + } + + while (scaled < T.One) + { + scaled *= scale; + powerOfTwo--; + } + + T root = TryRootThroughDouble(scaled, degree, out T scaledRoot) ? scaledRoot : T.One; + + for (; powerOfTwo > 0; powerOfTwo--) + { + root *= two; + } + + for (; powerOfTwo < 0; powerOfTwo++) + { + root /= two; + } + + return root; + } + + /// + /// Raises a value to a non-negative integer power, reporting rather than throwing when the result + /// leaves the range of the type. + /// + /// The numeric storage type. + /// The value to raise. + /// The power to raise it to. + /// The result, or the default when it does not fit. + /// when holds the result. + /// + /// The multiplications are checked, so a narrow type reports the overflow instead of wrapping into a + /// number that would send the Newton steps somewhere arbitrary. + /// + private static bool TryPower(T value, int exponent, out T power) + where T : struct, INumber + { + T result = T.One; + T factor = value; + + try + { + checked + { + for (int remaining = exponent; remaining > 0; remaining >>= 1) + { + if (int.IsOddInteger(remaining)) + { + result *= factor; + } + + if (remaining > 1) + { + factor *= factor; + } + } + } + } + catch (OverflowException) + { + power = default; + return false; + } + + power = result; + return true; + } } diff --git a/Semantics.Test/Quantities/StorageMathTests.cs b/Semantics.Test/Quantities/StorageMathTests.cs index ee0e0a48..7e2292b6 100644 --- a/Semantics.Test/Quantities/StorageMathTests.cs +++ b/Semantics.Test/Quantities/StorageMathTests.cs @@ -4,16 +4,19 @@ namespace ktsu.Semantics.Test.Quantities; using System; using System.Numerics; +using System.Reflection; using ktsu.Semantics.Quantities; using Microsoft.VisualStudio.TestTools.UnitTesting; /// -/// The square root every generated Length() and Distance() takes. +/// The square root every generated Length() and Distance() takes, and the higher roots +/// and hypotenuse that ship beside it. /// /// /// The primitives are asserted against the expression the generated code used to inline, so a change /// to their results fails here. The other types are asserted against the answer at their own -/// precision, which the old round trip through could not reach. +/// precision, which the old round trip through could not reach. The surface +/// itself is asserted too, since these are public API on a package with a compatibility baseline. /// [TestClass] public sealed class StorageMathTests @@ -128,4 +131,211 @@ public void DecimalExtremesHaveTheirRoots() Assert.AreEqual(0.00000000000001m, StorageMath.Sqrt(0.0000000000000000000000000001m)); Assert.IsLessThanOrEqualTo(0.00000000000002m, Math.Abs(StorageMath.Sqrt(decimal.MaxValue) - RootOfMaxValue)); } + + /// + /// The roots are public, because an application computing a norm the library does not emit has to + /// reach them or reimplement them. The seeding and the double round trip are not: they are how the + /// roots are taken, not what they promise. + /// + [TestMethod] + public void TheRootsArePublicAndTheirWorkingsAreNot() + { + Assert.IsTrue(typeof(StorageMath).IsPublic, "StorageMath itself has to be reachable from outside the assembly."); + + string[] contract = [nameof(StorageMath.Sqrt), nameof(StorageMath.Cbrt), nameof(StorageMath.RootN), nameof(StorageMath.Hypot)]; + string[] workings = ["IsRoundedThroughDouble", "IsBinaryFloatingPoint", "HasFloorDivision", "Seed", "SeedForRoot", "TryRootThroughDouble", "RootThroughDouble", "RootByNewton", "TryPower"]; + + foreach (string name in contract) + { + MethodInfo? method = typeof(StorageMath).GetMethod(name, BindingFlags.Public | BindingFlags.Static); + Assert.IsNotNull(method, $"{name} is part of the contract and has to be public."); + } + + foreach (string name in workings) + { + MethodInfo? method = typeof(StorageMath).GetMethod(name, BindingFlags.Public | BindingFlags.Static); + Assert.IsNull(method, $"{name} is implementation and should not be frozen into the public surface."); + } + } + + /// + /// and take , including the + /// negative cube root it defines where answers NaN. + /// + /// The value to take the root of. + [TestMethod] + [DataRow(8d)] + [DataRow(2d)] + [DataRow(0d)] + [DataRow(-27d)] + [DataRow(1e300)] + public void DoubleCubeRootIsUnchanged(double value) => Assert.AreEqual(Math.Cbrt(value), StorageMath.Cbrt(value)); + + /// + /// A cube root is refined in , so cubing it returns + /// further into the type than cubing the root the route gives. + /// + [TestMethod] + public void DecimalCubeRootIsMorePreciseThanTheDoubleRoute() + { + decimal refined = StorageMath.Cbrt(2m); + decimal throughDouble = decimal.CreateChecked(Math.Cbrt(2d)); + + decimal refinedError = Math.Abs((refined * refined * refined) - 2m); + decimal doubleError = Math.Abs((throughDouble * throughDouble * throughDouble) - 2m); + + Assert.IsLessThan(doubleError, refinedError, "The refinement should beat the route it is seeded from, or it is not earning its steps."); + Assert.IsLessThan(0.0000000000000000000000001m, refinedError); + } + + /// + /// A perfect cube comes back exactly, and a negative one keeps its sign, which is the cube root's + /// point of difference from the square root. + /// + /// The value to take the root of, as a literal. + /// Its exact root, as a literal. + [TestMethod] + [DataRow("8", "2")] + [DataRow("-8", "-2")] + [DataRow("1000000000000", "10000")] + [DataRow("0.000000000001", "0.0001")] + [DataRow("-27", "-3")] + public void DecimalPerfectCubeIsExact(string cube, string root) + { + decimal value = decimal.Parse(cube, System.Globalization.CultureInfo.InvariantCulture); + decimal expected = decimal.Parse(root, System.Globalization.CultureInfo.InvariantCulture); + + Assert.AreEqual(expected, StorageMath.Cbrt(value)); + } + + /// + /// An integer type is refined in integer arithmetic for a cube root, so it lands on the exact floor + /// rather than on whatever rounded to. + /// + [TestMethod] + public void IntegerCubeRootIsTheExactFloor() + { + Assert.AreEqual(2, StorageMath.Cbrt(26)); + Assert.AreEqual(3, StorageMath.Cbrt(27)); + Assert.AreEqual(3, StorageMath.Cbrt(28)); + Assert.AreEqual(-3, StorageMath.Cbrt(-27)); + Assert.AreEqual(1290, StorageMath.Cbrt(2147483647)); + } + + /// + /// A beyond the range of gets its exact cube root, + /// through the same scaling the square root uses. + /// + [TestMethod] + public void BigIntegerCubeRootBeyondTheRangeOfDoubleIsExact() + { + BigInteger tenToTheHundred = BigInteger.Pow(10, 100); + + Assert.AreEqual(tenToTheHundred, StorageMath.Cbrt(BigInteger.Pow(10, 300))); + Assert.AreEqual(tenToTheHundred, StorageMath.Cbrt(BigInteger.Pow(10, 300) + BigInteger.One)); + Assert.AreEqual(tenToTheHundred - BigInteger.One, StorageMath.Cbrt(BigInteger.Pow(10, 300) - BigInteger.One)); + Assert.AreEqual(BigInteger.Pow(2, 512), StorageMath.Cbrt(BigInteger.Pow(2, 1536))); + } + + /// + /// The first two degrees are the identity and , so a caller + /// parameterised on the degree does not have to special-case them. + /// + [TestMethod] + public void TheFirstTwoDegreesAreTheIdentityAndTheSquareRoot() + { + Assert.AreEqual(7m, StorageMath.RootN(7m, 1)); + Assert.AreEqual(StorageMath.Sqrt(2m), StorageMath.RootN(2m, 2)); + Assert.AreEqual(StorageMath.Sqrt(26), StorageMath.RootN(26, 2)); + } + + /// + /// A degree that is not a positive integer is a caller error rather than an answer. + /// + /// The degree to ask for. + [TestMethod] + [DataRow(0)] + [DataRow(-1)] + [DataRow(int.MinValue)] + public void ADegreeBelowOneThrows(int degree) + => Assert.ThrowsExactly(() => StorageMath.RootN(8m, degree)); + + /// + /// A higher root is exact where it terminates in the type, for a fractional type and an integer one. + /// + [TestMethod] + public void HigherRootsAreExactWhereTheyTerminate() + { + Assert.AreEqual(2m, StorageMath.RootN(16m, 4)); + Assert.AreEqual(3m, StorageMath.RootN(243m, 5)); + Assert.AreEqual(0.1m, StorageMath.RootN(0.00001m, 5)); + Assert.AreEqual(-2m, StorageMath.RootN(-32m, 5)); + Assert.AreEqual(new BigInteger(10), StorageMath.RootN(BigInteger.Pow(10, 7), 7)); + Assert.AreEqual(BigInteger.Pow(10, 100), StorageMath.RootN(BigInteger.Pow(10, 500), 5)); + } + + /// + /// A degree large enough that the estimate's own power leaves the type still answers, because the + /// power is taken checked and an estimate that overflows it is known to be above the root. + /// + [TestMethod] + public void ADegreeTooLargeForTheTypeToSquareStillAnswers() + { + Assert.AreEqual(1, StorageMath.RootN(100, 20)); + Assert.AreEqual(2, StorageMath.RootN(1048576, 20)); + Assert.AreEqual(1, StorageMath.RootN(int.MaxValue, 40)); + } + + /// + /// An even root of a negative value has no real answer, and says so exactly as the square root does. + /// + [TestMethod] + public void AnEvenRootOfANegativeValueHasNoAnswer() + { + Assert.IsTrue(double.IsNaN(StorageMath.RootN(-16d, 4))); + Assert.ThrowsExactly(static () => StorageMath.RootN(-16m, 4)); + } + + /// + /// The primitives take , which is the answer at their precision. + /// + [TestMethod] + public void PrimitiveHypotenuseIsTheDoubleOne() + { + Assert.AreEqual(double.Hypot(3d, 4d), StorageMath.Hypot(3d, 4d)); + Assert.AreEqual(double.Hypot(-3d, 4d), StorageMath.Hypot(-3d, 4d)); + Assert.AreEqual(5, StorageMath.Hypot(3, 4)); + Assert.AreEqual(5f, StorageMath.Hypot(3f, 4f)); + } + + /// + /// A pair whose squares leave the range of the type still has its hypotenuse, because a type with a + /// fractional part is computed from the ratio of the two legs rather than from their squares. + /// + [TestMethod] + public void DecimalHypotenuseSurvivesLegsWhoseSquaresDoNot() + { + Assert.AreEqual(5m, StorageMath.Hypot(3m, 4m)); + Assert.AreEqual(0m, StorageMath.Hypot(0m, 0m)); + Assert.AreEqual(4m, StorageMath.Hypot(0m, -4m)); + + decimal leg = 1e20m; + Assert.ThrowsExactly(() => leg * leg); + + decimal hypotenuse = StorageMath.Hypot(leg, leg); + decimal expected = leg * StorageMath.Sqrt(2m); + + Assert.IsLessThan(1e6m, Math.Abs(hypotenuse - expected)); + } + + /// + /// An integer type squares and sums instead, since the ratio of two integers is not a ratio. + /// + [TestMethod] + public void BigIntegerHypotenuseIsTheFloorOfTheAnswer() + { + Assert.AreEqual(new BigInteger(5), StorageMath.Hypot(new BigInteger(3), new BigInteger(-4))); + Assert.AreEqual(new BigInteger(2), StorageMath.Hypot(new BigInteger(2), new BigInteger(1))); + Assert.AreEqual(BigInteger.Pow(10, 200), StorageMath.Hypot(BigInteger.Pow(10, 200), BigInteger.Zero)); + } } diff --git a/docs/physics-generator.md b/docs/physics-generator.md index 88446fdc..2e5dcf65 100644 --- a/docs/physics-generator.md +++ b/docs/physics-generator.md @@ -215,6 +215,10 @@ that root with Newton steps in their own arithmetic. A value a `double` cannot h of four into [1, 4) for the seed, and a root that does not settle throws `ArithmeticException` rather than returning an estimate. +`StorageMath` is public, alongside `Cbrt`, `RootN` and `Hypot` on the same seeding and the same loop, +so an application computing a norm the generator does not emit reaches them rather than reimplementing +them. See the type's own documentation for what each one guarantees. + ## Validation, diagnostics, and gotchas - Unknown dimension references in `integrals` / `derivatives` / `dotProducts` / `crossProducts` report **SEM001** and the operator is dropped. From 748c65e07f1e9089c01156c93aea3bafd08383c9 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 16 Sep 2026 13:40:18 +0000 Subject: [PATCH 2/2] Split the root loop into its two disciplines [patch] MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit RootByNewton carried both loops in one method, switching on whether the type's division floors, which Sonar reads as cognitive complexity 19 against a limit of 15 — and it is right that the two have little to say to each other. RootByDescent stops on the first step that does not descend and answers with the estimate before it; RootBySettling stops when the estimate stops changing, or takes the smaller of a pair it alternates between. Each now says only its own rule. No behaviour changes: same seed, same step, same checked powers, and the suite still runs 1274 tests green. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01A9Ke7x6EeGyKHfPHncbNrP --- Semantics.Quantities/StorageMath.cs | 120 ++++++++++++------ Semantics.Test/Quantities/StorageMathTests.cs | 2 +- 2 files changed, 80 insertions(+), 42 deletions(-) diff --git a/Semantics.Quantities/StorageMath.cs b/Semantics.Quantities/StorageMath.cs index 73adc2b2..737bc169 100644 --- a/Semantics.Quantities/StorageMath.cs +++ b/Semantics.Quantities/StorageMath.cs @@ -461,7 +461,8 @@ private static bool HasFloorDivision() /// /// Refines the th root of a positive value by Newton steps in - /// 's own arithmetic. + /// 's own arithmetic, by whichever of the two disciplines the type's + /// division calls for. /// /// The numeric storage type. /// The value to take the root of, greater than zero and not one. @@ -469,75 +470,112 @@ private static bool HasFloorDivision() /// The root, or its floor for a type whose division floors. /// The estimate does not settle within steps. /// - /// - /// A type whose division floors descends: from an estimate at or above the root every step lands no - /// lower than the floor of the root, so the first step that does not descend has passed it and the - /// estimate before it is the answer. The seed is lifted to at or above the root first, since a - /// descent that starts below it would stop on the first step and answer with the seed. - /// - /// - /// Any other type settles instead, on the same terms as : an estimate that - /// stops changing is the root, and one alternating between two neighbours is rounding, so the - /// smaller of the two is taken. - /// - /// - /// Every power is taken checked, because a narrow type can hold a value whose estimate^(n-1) - /// it cannot, and an estimate that wraps sends the steps somewhere arbitrary. A degree so high that - /// the type cannot hold 2^degree is answered before any of that: no value the type holds has - /// a root of two or more, so the answer is one. - /// + /// Every power either loop takes is checked, because a narrow type can hold a value whose + /// estimate^(degree - 1) it cannot, and an estimate that wraps sends the steps somewhere + /// arbitrary. The step itself is x = (((degree - 1) * x) + (value / x^(degree - 1))) / degree + /// in both. /// private static T RootByNewton(T value, int degree) where T : struct, INumber { T two = T.One + T.One; + + return HasFloorDivision() + ? RootByDescent(value, degree, two) + : RootBySettling(value, degree, two); + } + + /// + /// Takes the root in a type whose division floors, by descending to it. + /// + /// The numeric storage type, whose division floors. + /// The value to take the root of, greater than one. + /// The degree of the root, three or more. + /// Two, in . + /// The floor of the root. + /// The descent does not reach the root within steps. + /// + /// From an estimate at or above the root every step lands no lower than the floor of the root, so + /// the first step that does not descend has passed it and the estimate before it is the answer. The + /// seed is lifted to at or above the root first, since a descent that starts below it would stop on + /// the first step and answer with the seed. A degree so high that the type cannot hold + /// 2^degree is answered before any of that: no value the type holds has a root of two or + /// more, and the values with a root below one were answered before this was called. + /// + private static T RootByDescent(T value, int degree, T two) + where T : struct, INumber + { + if (!TryPower(two, degree, out _)) + { + return T.One; + } + T count = T.CreateChecked(degree); T countLessOne = count - T.One; - bool floors = HasFloorDivision(); T estimate = SeedForRoot(value, degree, two); - if (floors) + // An estimate whose own power overflows is above the root, since the root's power divides a + // value the type holds, so the lift stops there as well as on a power that exceeds the value. + while (TryPower(estimate, degree, out T raised) && raised < value) + { + estimate *= two; + } + + for (int iteration = 0; iteration < MaximumIterations; iteration++) { - if (!TryPower(two, degree, out _)) + if (!TryPower(estimate, degree - 1, out T power) || T.IsZero(power)) { - // Two raised to this degree leaves the type, so no value the type holds has a root of two - // or more, and the values with a root below one were answered before this was called. - return T.One; + // The estimate is too large for its own power to be taken in the type. Halving reaches a + // range the type holds, and stays at or above the root. + estimate = (estimate + T.One) / two; + continue; } - // An estimate whose own power overflows is above the root, since the root's power divides - // a value the type holds, so the lift stops there as well as on a power that exceeds the value. - while (TryPower(estimate, degree, out T raised) && raised < value) + T next = ((countLessOne * estimate) + (value / power)) / count; + + if (next >= estimate) { - estimate *= two; + return estimate; } + + estimate = next; } + throw new ArithmeticException($"The root of degree {degree} did not reach its floor within {MaximumIterations} Newton steps."); + } + + /// + /// Takes the root in a type that keeps a fractional part, by settling on it. + /// + /// The numeric storage type. + /// The value to take the root of, greater than zero and not one. + /// The degree of the root, three or more. + /// Two, in . + /// The root, to the precision the type holds. + /// The estimate does not settle within steps. + /// + /// On the same terms as : an estimate that stops changing is the root, and + /// one alternating between two neighbours is rounding, so the smaller of the two is taken. + /// + private static T RootBySettling(T value, int degree, T two) + where T : struct, INumber + { + T count = T.CreateChecked(degree); + T countLessOne = count - T.One; + T estimate = SeedForRoot(value, degree, two); T previous = estimate; for (int iteration = 0; iteration < MaximumIterations; iteration++) { if (!TryPower(estimate, degree - 1, out T power) || T.IsZero(power)) { - // The estimate is too large for its own power to be taken in the type. Halving reaches a - // range the type holds, and for a type that floors it stays at or above the root. + // The estimate is too far from the root for its own power to be taken in the type. estimate = (estimate + T.One) / two; continue; } T next = ((countLessOne * estimate) + (value / power)) / count; - if (floors) - { - if (next >= estimate) - { - return estimate; - } - - estimate = next; - continue; - } - if (next == estimate) { return next; diff --git a/Semantics.Test/Quantities/StorageMathTests.cs b/Semantics.Test/Quantities/StorageMathTests.cs index 7e2292b6..e8df40b8 100644 --- a/Semantics.Test/Quantities/StorageMathTests.cs +++ b/Semantics.Test/Quantities/StorageMathTests.cs @@ -143,7 +143,7 @@ public void TheRootsArePublicAndTheirWorkingsAreNot() Assert.IsTrue(typeof(StorageMath).IsPublic, "StorageMath itself has to be reachable from outside the assembly."); string[] contract = [nameof(StorageMath.Sqrt), nameof(StorageMath.Cbrt), nameof(StorageMath.RootN), nameof(StorageMath.Hypot)]; - string[] workings = ["IsRoundedThroughDouble", "IsBinaryFloatingPoint", "HasFloorDivision", "Seed", "SeedForRoot", "TryRootThroughDouble", "RootThroughDouble", "RootByNewton", "TryPower"]; + string[] workings = ["IsRoundedThroughDouble", "IsBinaryFloatingPoint", "HasFloorDivision", "Seed", "SeedForRoot", "TryRootThroughDouble", "RootThroughDouble", "RootByNewton", "RootByDescent", "RootBySettling", "TryPower"]; foreach (string name in contract) {