Account for over and underflow in quaternion length

Author: markuswntr

Sources/QuaternionModule/Norms.swift

@@ -77,8 +77,19 @@ extension Quaternion {
   /// - `.lengthSquared`
   @_transparent
   public var length: RealType {
+    let naive = lengthSquared
+    guard naive.isNormal else { return carefulLength }
+    return .sqrt(naive)
+  }
+
+  //  Internal implementation detail of `length`, moving slow path off
+  //  of the inline function.
+  @usableFromInline
+  internal var carefulLength: RealType {
     guard isFinite else { return .infinity }
-    return .sqrt(lengthSquared)
+    guard !magnitude.isZero else { return .zero }
+    // Unscale the quaternion, calculate its length and rescale the result
+    return divided(by: magnitude).length * magnitude
   }
 
   /// The squared length `(r*r + x*x + y*y + z*z)`.

Tests/QuaternionTests/PropertyTests.swift

@@ -34,6 +34,11 @@ final class PropertyTests: XCTestCase {
     XCTAssertEqual(Quaternion<T>.zero.length, .zero)
     XCTAssertEqual(Quaternion<T>(real: .zero, imaginary: -.zero).length, .zero)
     XCTAssertEqual(Quaternion<T>(real: -.zero, imaginary: -.zero).length, .zero)
+    // Check for overflow and underflow when calculating the length
+    XCTAssertEqual(Quaternion<T>(real: .greatestFiniteMagnitude, imaginary: 0, 0, 0).length, .greatestFiniteMagnitude)
+    XCTAssertEqual(Quaternion<T>(real: 0, imaginary: -.greatestFiniteMagnitude, 0, 0).length, .greatestFiniteMagnitude)
+    XCTAssertEqual(Quaternion<T>(real: 0, imaginary: 0, .leastNormalMagnitude, 0).length, .leastNormalMagnitude)
+    XCTAssertEqual(Quaternion<T>(real: 0, imaginary: 0, 0, -.leastNormalMagnitude).length, .leastNormalMagnitude)
   }
 
   func testProperties() {