Compare every simd-backed function against Apple's simd and fix the differences found

.github/workflows/ci-macos.yml

@@ -38,5 +38,8 @@ jobs:
       - name: Run Tests (Release, portable path vs Apple simd)
         run: swift test -c release --disable-default-traits --parallel -Xswiftc -warnings-as-errors
 
+      - name: Check Differential Test Coverage
+        run: python3 Scripts/check-simd-differential-coverage.py
+
       - name: Check @inlinable Coverage
         run: Scripts/check-inlinable.sh

Makefile

@@ -11,7 +11,7 @@ DOCS_VERSION_PATH ?= main
 HOSTING_BASE_PATH ?= $(REPO_NAME)/$(DOCS_VERSION_PATH)
 
 # Targets
-.PHONY: setup lint lint-fix check-inlinable test test-coverage clean pre-commit docs docs-preview docs-generate docs-coverage profile-stats benchmark-commits benchmark benchmark-thresholds-check benchmark-thresholds-update
+.PHONY: setup lint lint-fix check-inlinable check-simd-coverage test-portable test test-coverage clean pre-commit docs docs-preview docs-generate docs-coverage profile-stats benchmark-commits benchmark benchmark-thresholds-check benchmark-thresholds-update
 
 setup:
 	@echo "Detected Package Swift Version: $(PACKAGE_SWIFT_VERSION)"
@@ -52,6 +52,13 @@ lint-fix:
 check-inlinable:
 	./Scripts/check-inlinable.sh
 
+check-simd-coverage:
+	python3 Scripts/check-simd-differential-coverage.py
+
+# Portable implementations compared against Apple's simd (Apple platforms only).
+test-portable:
+	swift test $(SWIFT_FLAGS) --disable-default-traits --parallel -Xswiftc -warnings-as-errors
+
 test:
 	swift test $(SWIFT_FLAGS) --parallel --enable-code-coverage -Xswiftc -warnings-as-errors
 

Scripts/check-simd-differential-coverage.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python3
+"""Checks that every FirebladeMath function with a simd branch has a differential test.
+
+Functions whose body contains `#if FRB_MATH_USE_SIMD` have a simd implementation and a
+portable one. The PortableVsSIMD tests compare the two; each test names the functions it
+covers with `// covers: name(ParameterTypes)` comments, e.g.
+
+    // covers: multiply(Mat4x4f, Mat4x4f)
+    // covers: quaternion(angle: Float, axis: SIMD3<Float>)
+
+Fails if a simd-backed function has no `covers:` annotation, or if an annotation names
+a function that does not exist (or no longer has a simd branch).
+"""
+
+import pathlib
+import re
+import sys
+
+ROOT = pathlib.Path(__file__).resolve().parent.parent
+SOURCES = ROOT / "Sources" / "FirebladeMath"
+TESTS = ROOT / "Tests" / "FirebladeMathTests" / "PortableVsSIMD"
+
+DECLARATION = re.compile(
+    r"^(?P<attributes>(?:[ \t]*@[^\n]*\n)*)[ \t]*public\s+(?:static\s+)?func\s+(?P<name>\w+)\s*\((?P<parameters>[^)]*)\)[^{]*\{",
+    re.MULTILINE,
+)
+COVERS = re.compile(r"^[ \t]*//[ \t]*covers:[ \t]*(.+?)[ \t]*$", re.MULTILINE)
+
+
+def split_parameters(parameters):
+    """Splits a parameter list on top-level commas (ignoring commas inside <...>)."""
+    parts, depth, current = [], 0, ""
+    for character in parameters:
+        if character == "<":
+            depth += 1
+        elif character == ">":
+            depth -= 1
+        if character == "," and depth == 0:
+            parts.append(current)
+            current = ""
+        else:
+            current += character
+    if current.strip():
+        parts.append(current)
+    return [part.strip() for part in parts]
+
+
+def key(name, parameters):
+    """`name(Type, label: Type)`: unlabeled parameters by type, labeled ones with their label."""
+    rendered = []
+    for parameter in split_parameters(parameters):
+        names, type_name = parameter.split(":", 1)
+        external = names.split()[0]
+        type_name = " ".join(type_name.split())
+        rendered.append(type_name if external == "_" else f"{external}: {type_name}")
+    return f"{name}({', '.join(rendered)})"
+
+
+def body_end(text, start):
+    depth = 1
+    index = start
+    while depth and index < len(text):
+        depth += {"{": 1, "}": -1}.get(text[index], 0)
+        index += 1
+    return index
+
+
+def simd_backed_functions():
+    functions = {}
+    for path in sorted(SOURCES.rglob("*.swift")):
+        text = path.read_text()
+        for match in DECLARATION.finditer(text):
+            if "deprecated" in match.group("attributes"):
+                continue
+            body = text[match.end():body_end(text, match.end())]
+            if "FRB_MATH_USE_SIMD" in body:
+                functions[key(match.group("name"), match.group("parameters"))] = path.relative_to(ROOT)
+    return functions
+
+
+def covered_functions():
+    covered = {}
+    for path in sorted(TESTS.rglob("*.swift")):
+        for match in COVERS.finditer(path.read_text()):
+            covered[" ".join(match.group(1).split())] = path.relative_to(ROOT)
+    return covered
+
+
+def main():
+    functions = simd_backed_functions()
+    covered = covered_functions()
+
+    missing = sorted(set(functions) - set(covered))
+    unknown = sorted(set(covered) - set(functions))
+
+    if missing:
+        print(f"error: {len(missing)} simd-backed function(s) have no differential test (`// covers:` annotation):")
+        for name in missing:
+            print(f"  {name}    [{functions[name]}]")
+    if unknown:
+        print(f"error: {len(unknown)} `// covers:` annotation(s) name no simd-backed function:")
+        for name in unknown:
+            print(f"  {name}    [{covered[name]}]")
+    if missing or unknown:
+        return 1
+
+    print(f"All {len(functions)} simd-backed functions have differential tests.")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

Sources/FirebladeMath/Functions/clamp.swift

@@ -26,7 +26,8 @@ public func clamp(_ x: Double, _ minVal: Double, _ maxVal: Double) -> Double {
     #if FRB_MATH_USE_SIMD
     return simd_clamp(x, minVal, maxVal)
     #else
-    return Swift.min(Swift.max(x, minVal), maxVal)
+    // IEEE maxNum/minNum like simd_clamp: a NaN `x` yields `minVal`.
+    return Double.minimum(Double.maximum(x, minVal), maxVal)
     #endif
 }
 
@@ -42,7 +43,8 @@ public func clamp(_ x: Float, _ minVal: Float, _ maxVal: Float) -> Float {
     #if FRB_MATH_USE_SIMD
     return simd_clamp(x, minVal, maxVal)
     #else
-    return Swift.min(Swift.max(x, minVal), maxVal)
+    // IEEE maxNum/minNum like simd_clamp: a NaN `x` yields `minVal`.
+    return Float.minimum(Float.maximum(x, minVal), maxVal)
     #endif
 }
 

Sources/FirebladeMath/Functions/fract.swift

@@ -15,7 +15,8 @@ public func fract(_ value: Double) -> Double {
     #if FRB_MATH_USE_SIMD
     return simd_fract(value)
     #else
-    return Swift.min(value - floor(value), Double(0x1.FFFFFFFFFFFFFp-1))
+    // IEEE minNum like simd_fract: NaN and infinities yield the largest value below 1.
+    return Double.minimum(value - floor(value), Double(0x1.FFFFFFFFFFFFFp-1))
     #endif
 }
 
@@ -32,6 +33,7 @@ public func fract(_ value: Float) -> Float {
     #if FRB_MATH_USE_SIMD
     return simd_fract(value)
     #else
-    return Swift.min(value - floor(value), Float(0x1.FFFFFEp-1))
+    // IEEE minNum like simd_fract: NaN and infinities yield the largest value below 1.
+    return Float.minimum(value - floor(value), Float(0x1.FFFFFEp-1))
     #endif
 }

Sources/FirebladeMath/Functions/matrix.swift

@@ -12,16 +12,19 @@ public func matrix3x3(from quat: Quat4f) -> Mat3x3f {
     return Mat3x3f(storage: simd_matrix3x3(quat.storage))
     #else
     let v = Vec4f(quat)
+    // Same diagonal as Apple's simd_matrix3x3/simd_matrix4x4, so non-unit quaternions
+    // produce the same (scaled) matrix on every platform.
+    let xx: Float = v.x * v.x, yy: Float = v.y * v.y, zz: Float = v.z * v.z, ww: Float = v.w * v.w
     return Mat3x3f(
-        Vec3f(1 - 2 * (v.y * v.y + v.z * v.z),
+        Vec3f(xx - yy - zz + ww,
               2 * (v.x * v.y + v.z * v.w),
               2 * (v.x * v.z - v.y * v.w)),
         Vec3f(2 * (v.x * v.y - v.z * v.w),
-              1 - 2 * (v.z * v.z + v.x * v.x),
+              yy - zz + ww - xx,
               2 * (v.y * v.z + v.x * v.w)),
         Vec3f(2 * (v.z * v.x + v.y * v.w),
               2 * (v.y * v.z - v.x * v.w),
-              1 - 2 * (v.y * v.y + v.x * v.x))
+              zz + ww - xx - yy)
     )
     #endif
 }
@@ -35,16 +38,19 @@ public func matrix3x3(from quat: Quat4d) -> Mat3x3d {
     return Mat3x3d(storage: simd_matrix3x3(quat.storage))
     #else
     let v = Vec4d(quat)
+    // Same diagonal as Apple's simd_matrix3x3/simd_matrix4x4, so non-unit quaternions
+    // produce the same (scaled) matrix on every platform.
+    let xx: Double = v.x * v.x, yy: Double = v.y * v.y, zz: Double = v.z * v.z, ww: Double = v.w * v.w
     return Mat3x3d(
-        Vec3d(1 - 2 * (v.y * v.y + v.z * v.z),
+        Vec3d(xx - yy - zz + ww,
               2 * (v.x * v.y + v.z * v.w),
               2 * (v.x * v.z - v.y * v.w)),
         Vec3d(2 * (v.x * v.y - v.z * v.w),
-              1 - 2 * (v.z * v.z + v.x * v.x),
+              yy - zz + ww - xx,
               2 * (v.y * v.z + v.x * v.w)),
         Vec3d(2 * (v.z * v.x + v.y * v.w),
               2 * (v.y * v.z - v.x * v.w),
-              1 - 2 * (v.y * v.y + v.x * v.x))
+              zz + ww - xx - yy)
     )
     #endif
 }
@@ -58,23 +64,26 @@ public func matrix4x4(from quat: Quat4f) -> Mat4x4f {
     return Mat4x4f(storage: simd_matrix4x4(quat.storage))
     #else
     let v = Vec4f(quat)
+    // Same diagonal as Apple's simd_matrix3x3/simd_matrix4x4, so non-unit quaternions
+    // produce the same (scaled) matrix on every platform.
+    let xx: Float = v.x * v.x, yy: Float = v.y * v.y, zz: Float = v.z * v.z, ww: Float = v.w * v.w
     return Mat4x4f(
         Vec4f(
-            1 - 2 * (v.y * v.y + v.z * v.z),
+            xx - yy - zz + ww,
             2 * (v.x * v.y + v.z * v.w),
             2 * (v.x * v.z - v.y * v.w),
             0
         ),
         Vec4f(
             2 * (v.x * v.y - v.z * v.w),
-            1 - 2 * (v.z * v.z + v.x * v.x),
+            yy - zz + ww - xx,
             2 * (v.y * v.z + v.x * v.w),
             0
         ),
         Vec4f(
             2 * (v.z * v.x + v.y * v.w),
             2 * (v.y * v.z - v.x * v.w),
-            1 - 2 * (v.y * v.y + v.x * v.x),
+            zz + ww - xx - yy,
             0
         ),
         Vec4f(0, 0, 0, 1)
@@ -91,23 +100,26 @@ public func matrix4x4(from quat: Quat4d) -> Mat4x4d {
     return Mat4x4d(storage: simd_matrix4x4(quat.storage))
     #else
     let v = Vec4d(quat)
+    // Same diagonal as Apple's simd_matrix3x3/simd_matrix4x4, so non-unit quaternions
+    // produce the same (scaled) matrix on every platform.
+    let xx: Double = v.x * v.x, yy: Double = v.y * v.y, zz: Double = v.z * v.z, ww: Double = v.w * v.w
     return Mat4x4d(
         Vec4d(
-            1 - 2 * (v.y * v.y + v.z * v.z),
+            xx - yy - zz + ww,
             2 * (v.x * v.y + v.z * v.w),
             2 * (v.x * v.z - v.y * v.w),
             0
         ),
         Vec4d(
             2 * (v.x * v.y - v.z * v.w),
-            1 - 2 * (v.z * v.z + v.x * v.x),
+            yy - zz + ww - xx,
             2 * (v.y * v.z + v.x * v.w),
             0
         ),
         Vec4d(
             2 * (v.z * v.x + v.y * v.w),
             2 * (v.y * v.z - v.x * v.w),
-            1 - 2 * (v.y * v.y + v.x * v.x),
+            zz + ww - xx - yy,
             0
         ),
         Vec4d(0, 0, 0, 1)

Sources/FirebladeMath/Functions/normalize.swift

@@ -88,28 +88,40 @@ public func normalize(_ x: SIMD4<Float>) -> SIMD4<Float> {
 
 /// Returns a normalized copy of the quaternion.
 /// - Parameter x: The quaternion to normalize.
-/// - Returns: A quaternion with the same orientation as x but with a length of 1.0.
+/// - Returns: A quaternion with the same orientation as x but with a length of 1.0,
+///   or the identity quaternion if x has zero length.
 @inlinable
 @_disfavoredOverload
 public func normalize(_ x: Quat4f) -> Quat4f {
     #if FRB_MATH_USE_SIMD
     return Quat4f(storage: simd.simd_normalize(x.storage))
     #else
-    let scale: Float = 1.0 / length(x)
+    // Like Apple's simd, a quaternion whose squared length is zero normalizes to the identity.
+    let lengthSquared: Float = x.x * x.x + x.y * x.y + x.z * x.z + x.w * x.w
+    guard lengthSquared != 0 else {
+        return Quat4f()
+    }
+    let scale: Float = 1.0 / lengthSquared.squareRoot()
     return Quat4f(x.x * scale, x.y * scale, x.z * scale, x.w * scale)
     #endif
 }
 
 /// Returns a normalized copy of the quaternion.
 /// - Parameter x: The quaternion to normalize.
-/// - Returns: A quaternion with the same orientation as x but with a length of 1.0.
+/// - Returns: A quaternion with the same orientation as x but with a length of 1.0,
+///   or the identity quaternion if x has zero length.
 @inlinable
 @_disfavoredOverload
 public func normalize(_ x: Quat4d) -> Quat4d {
     #if FRB_MATH_USE_SIMD
     return Quat4d(storage: simd.simd_normalize(x.storage))
     #else
-    let scale = 1.0 / length(x)
+    // Like Apple's simd, a quaternion whose squared length is zero normalizes to the identity.
+    let lengthSquared: Double = x.x * x.x + x.y * x.y + x.z * x.z + x.w * x.w
+    guard lengthSquared != 0 else {
+        return Quat4d()
+    }
+    let scale: Double = 1.0 / lengthSquared.squareRoot()
     return Quat4d(x.x * scale, x.y * scale, x.z * scale, x.w * scale)
     #endif
 }