diff --git a/Semantics.Quantities/StorageMath.cs b/Semantics.Quantities/StorageMath.cs index 94a752f..662f5d6 100644 --- a/Semantics.Quantities/StorageMath.cs +++ b/Semantics.Quantities/StorageMath.cs @@ -265,11 +265,12 @@ public static T RootN(T value, int n) /// 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 - /// . The 64-bit and 128-bit integers do not: they sum their - /// squares in when they could overflow , and return - /// the exact floor, which throws only when - /// cannot hold it (#520). + /// The floating point primitives take the round trip, so Hypot<double> + /// is exactly . The integer primitives do not: the 64-bit and + /// 128-bit ones sum their squares in when they could overflow + /// (#520), the narrower ones sum theirs in (#537), and + /// both return the exact floor, which throws only when + /// cannot hold it. /// /// public static T Hypot(T x, T y) @@ -280,6 +281,11 @@ public static T Hypot(T x, T y) return HypotOfWideInteger(x, y, T.Zero, T.Zero); } + if (IsNarrowInteger()) + { + return HypotOfNarrowInteger(x, y, T.Zero, T.Zero); + } + if (IsRoundedThroughDouble()) { return T.CreateChecked(double.Hypot(double.CreateChecked(x), double.CreateChecked(y))); @@ -319,6 +325,11 @@ internal static T Hypot(T x, T y, T z) return HypotOfWideInteger(x, y, z, T.Zero); } + if (IsNarrowInteger()) + { + return HypotOfNarrowInteger(x, y, z, T.Zero); + } + if (IsRoundedThroughDouble()) { double dx = double.CreateChecked(x); @@ -361,6 +372,11 @@ internal static T Hypot(T x, T y, T z, T w) return HypotOfWideInteger(x, y, z, w); } + if (IsNarrowInteger()) + { + return HypotOfNarrowInteger(x, y, z, w); + } + if (IsRoundedThroughDouble()) { double dx = double.CreateChecked(x); @@ -433,6 +449,52 @@ private static bool IsWideInteger() || typeof(T) == typeof(UInt128) || ((typeof(T) == typeof(nint) || typeof(T) == typeof(nuint)) && nint.Size == 8); + /// + /// Reports whether is an integer primitive of 32 bits or fewer, whose + /// squares and their sums fit in however large the components are. + /// + /// The numeric storage type. + /// + /// for , , , + /// , , , and and + /// on a 32-bit process. + /// + private static bool IsNarrowInteger() + where T : struct, INumber + => typeof(T) == typeof(byte) + || typeof(T) == typeof(sbyte) + || typeof(T) == typeof(short) + || typeof(T) == typeof(ushort) + || typeof(T) == typeof(int) + || typeof(T) == typeof(uint) + || ((typeof(T) == typeof(nint) || typeof(T) == typeof(nuint)) && nint.Size == 4); + + /// + /// The exact floor of sqrt(x² + y² + z² + w²) for a narrow integer type. + /// + /// A type accepts. + /// The first component. + /// The second component. + /// The third component, zero for fewer components. + /// The fourth component, zero for fewer components. + /// The floor of the length. + /// The length does not fit in . + /// + /// Each component is at most 2^32 in magnitude, so four squares sum within 2^66, which + /// holds exactly. Nesting instead + /// rounds once per call, which can leave an exact length one ulp low, and truncating that back to + /// then answers one too low on a perfect square (#537). + /// + private static T HypotOfNarrowInteger(T x, T y, T z, T w) + where T : struct, INumber + { + Int128 wx = Int128.CreateChecked(x); + Int128 wy = Int128.CreateChecked(y); + Int128 wz = Int128.CreateChecked(z); + Int128 ww = Int128.CreateChecked(w); + return T.CreateChecked(FloorSqrtOfWideInteger((wx * wx) + (wy * wy) + (wz * wz) + (ww * ww))); + } + /// /// The exact floor of the square root of a non-negative wide integer. /// diff --git a/Semantics.Test/Quantities/StorageMathPreciseAndVectorTests.cs b/Semantics.Test/Quantities/StorageMathPreciseAndVectorTests.cs index 32eef76..1c4450b 100644 --- a/Semantics.Test/Quantities/StorageMathPreciseAndVectorTests.cs +++ b/Semantics.Test/Quantities/StorageMathPreciseAndVectorTests.cs @@ -59,6 +59,18 @@ public void IntegerVectorLengthsDoNotWrap() Assert.AreEqual(50000, Position3D.FromMeter(50000, 0, 0).Length()); } + [TestMethod] + public void IntegerVectorLengthsOfPerfectSquaresAreExact() + { + // 36 + 1681 + 1764 = 3481 = 59², and 9 + 9 + 9 + 9 = 36 = 6² (#537). + Position3D position = Position3D.FromMeter(6, 41, 42); + Assert.AreEqual(59, position.Length()); + Assert.AreEqual(59, position.Magnitude().Value); + Assert.AreEqual(59, position.DistanceTo(Position3D.Zero).Value); + Assert.AreEqual((short)59, Position3D.FromMeter(6, 41, 42).Length()); + Assert.AreEqual(6, new Acceleration4D { X = 3, Y = 3, Z = 3, W = 3 }.Length()); + } + [TestMethod] public void DecimalVectorLengthPastTheSquareOfItsRange() { diff --git a/Semantics.Test/Quantities/StorageMathTests.cs b/Semantics.Test/Quantities/StorageMathTests.cs index 31c5b52..67a8907 100644 --- a/Semantics.Test/Quantities/StorageMathTests.cs +++ b/Semantics.Test/Quantities/StorageMathTests.cs @@ -180,7 +180,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", "RootByDescent", "RootBySettling", "TryPower", "IsWideInteger", "FloorSqrtOfWideInteger", "HypotOfWideInteger"]; + string[] workings = ["IsRoundedThroughDouble", "IsBinaryFloatingPoint", "HasFloorDivision", "Seed", "SeedForRoot", "TryRootThroughDouble", "RootThroughDouble", "RootByNewton", "RootByDescent", "RootBySettling", "TryPower", "IsWideInteger", "FloorSqrtOfWideInteger", "HypotOfWideInteger", "IsNarrowInteger", "HypotOfNarrowInteger"]; foreach (string name in contract) { @@ -410,6 +410,77 @@ public void WideIntegerHypotenuseIsTheExactFloor() Assert.AreEqual(2L, StorageMath.Hypot(1L, 1L, 1L, 1L)); } + /// + /// A 32-bit-or-narrower integer vector whose length is a perfect square gets that length, not one + /// less: nesting rounds twice and can land one ulp low, + /// which truncation then turned into an off-by-one (#537). + /// + [TestMethod] + public void NarrowIntegerHypotenuseIsTheExactFloor() + { + Assert.AreEqual(59, StorageMath.Hypot(6, 41, 42)); + Assert.AreEqual(59, StorageMath.Hypot(-6, 41, -42)); + Assert.AreEqual(59, StorageMath.Hypot(14, 39, 42)); + Assert.AreEqual((short)59, StorageMath.Hypot((short)6, (short)41, (short)42)); + Assert.AreEqual(59u, StorageMath.Hypot(6u, 41u, 42u)); + Assert.AreEqual(6, StorageMath.Hypot(3, 3, 3, 3)); + Assert.AreEqual((short)6, StorageMath.Hypot((short)3, (short)3, (short)3, (short)3)); + Assert.AreEqual((byte)15, StorageMath.Hypot((byte)2, (byte)3, (byte)4, (byte)14)); + Assert.AreEqual(15, StorageMath.Hypot(1, 4, 8, 12)); + Assert.AreEqual(5, StorageMath.Hypot(3, -4)); + + // Squares of the extremes sum far past 64 bits, and the length still fits the type. + Assert.AreEqual(uint.MaxValue, StorageMath.Hypot(uint.MaxValue, 0u, 0u, 0u)); + Assert.AreEqual(int.MaxValue, StorageMath.Hypot(int.MinValue + 1, 0, 0)); + Assert.ThrowsExactly(static () => StorageMath.Hypot(int.MinValue, 0, 0)); + Assert.ThrowsExactly(static () => StorageMath.Hypot(uint.MaxValue, uint.MaxValue, uint.MaxValue, uint.MaxValue)); + } + + /// + /// A narrow integer gets the exact floor wherever the components lie: every 3D vector with + /// components up to 60, and every 4D vector up to 20, agrees with an integer square root of the + /// integer sum of squares (#537). + /// + [TestMethod] + public void NarrowIntegerHypotenuseMatchesTheIntegerRootEverywhere() + { + for (int x = 0; x <= 60; x++) + { + for (int y = 0; y <= 60; y++) + { + for (int z = 0; z <= 60; z++) + { + Assert.AreEqual(FloorSqrt((x * x) + (y * y) + (z * z)), StorageMath.Hypot(x, y, z), $"({x}, {y}, {z})"); + } + } + } + + for (int x = 0; x <= 20; x++) + { + for (int y = 0; y <= 20; y++) + { + for (int z = 0; z <= 20; z++) + { + for (int w = 0; w <= 20; w++) + { + Assert.AreEqual(FloorSqrt((x * x) + (y * y) + (z * z) + (w * w)), StorageMath.Hypot(x, y, z, w), $"({x}, {y}, {z}, {w})"); + } + } + } + } + + static int FloorSqrt(int value) + { + int root = 0; + while ((root + 1) * (root + 1) <= value) + { + root++; + } + + return root; + } + } + /// /// A wide-integer hypotenuse the type cannot hold is an overflow, as it was before (#520). ///