diff --git a/CLAUDE.md b/CLAUDE.md index 0d53df7..318ab86 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 fe9575c..b2c8b10 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 561b8ae..737bc16 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,299 @@ 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, 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. + /// 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. + /// + /// 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; + T estimate = SeedForRoot(value, degree, two); + + // 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(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 stays at or above the root. + estimate = (estimate + T.One) / two; + continue; + } + + T next = ((countLessOne * estimate) + (value / power)) / count; + + if (next >= estimate) + { + 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 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 (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 ee0e0a4..e8df40b 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", "RootByDescent", "RootBySettling", "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 88446fd..2e5dcf6 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.