Speed up the portable 4x4 inverse, adjugate and quaternion multiply

.github/workflows/ci-macos.yml

@@ -33,5 +33,10 @@ jobs:
       - name: Run Tests (Release)
         run: swift test -c release --parallel -Xswiftc -warnings-as-errors
 
+      # Without the simd trait, FirebladeMath uses its portable implementations while Apple's
+      # simd is still importable, so PortableVsSIMDTests compares the two.
+      - name: Run Tests (Release, portable path vs Apple simd)
+        run: swift test -c release --disable-default-traits --parallel -Xswiftc -warnings-as-errors
+
       - name: Check @inlinable Coverage
         run: Scripts/check-inlinable.sh

Sources/FirebladeMath/Functions/adjugate.swift

@@ -47,124 +47,13 @@ public func adjugate(_ inMat: Mat3x3d) -> Mat3x3d {
 /// - Returns: The adjugate of the input matrix.
 @inlinable
 public func adjugate(_ matrix: Mat4x4f) -> Mat4x4f {
-    // Adjugate is the transpose of the cofactor matrix.
-    // cofactor(i, j) = (-1)^(i+j) * minor(i, j)
-    // Here we use matrix[col, row] indexing.
-    var res = Mat4x4f.identity
-
-    res[0, 0] = matrix[1, 1] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 1] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        + matrix[3, 1] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2])
-    res[1, 0] = -(matrix[1, 0] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 0] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        + matrix[3, 0] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2]))
-    res[2, 0] = matrix[1, 0] * (matrix[2, 1] * matrix[3, 3] - matrix[2, 3] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[1, 1] * matrix[3, 3] - matrix[1, 3] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[1, 1] * matrix[2, 3] - matrix[1, 3] * matrix[2, 1])
-    res[3, 0] = -(matrix[1, 0] * (matrix[2, 1] * matrix[3, 2] - matrix[2, 2] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[1, 1] * matrix[3, 2] - matrix[1, 2] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[1, 1] * matrix[2, 2] - matrix[1, 2] * matrix[2, 1]))
-
-    res[0, 1] = -(matrix[0, 1] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 1] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 1] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2]))
-    res[1, 1] = matrix[0, 0] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 0] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 0] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2])
-    res[2, 1] = -(matrix[0, 0] * (matrix[2, 1] * matrix[3, 3] - matrix[2, 3] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[0, 1] * matrix[3, 3] - matrix[0, 3] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[2, 3] - matrix[0, 3] * matrix[2, 1]))
-    res[3, 1] = matrix[0, 0] * (matrix[2, 1] * matrix[3, 2] - matrix[2, 2] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[0, 1] * matrix[3, 2] - matrix[0, 2] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[2, 2] - matrix[0, 2] * matrix[2, 1])
-
-    res[0, 2] = matrix[0, 1] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        - matrix[1, 1] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 1] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2])
-    res[1, 2] = -(matrix[0, 0] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        - matrix[1, 0] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 0] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2]))
-    res[2, 2] = matrix[0, 0] * (matrix[1, 1] * matrix[3, 3] - matrix[1, 3] * matrix[3, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[3, 3] - matrix[0, 3] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[1, 3] - matrix[0, 3] * matrix[1, 1])
-    res[3, 2] = -(matrix[0, 0] * (matrix[1, 1] * matrix[3, 2] - matrix[1, 2] * matrix[3, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[3, 2] - matrix[0, 2] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[1, 2] - matrix[0, 2] * matrix[1, 1]))
-
-    res[0, 3] = -(matrix[0, 1] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2])
-        - matrix[1, 1] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2])
-        + matrix[2, 1] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2]))
-    res[1, 3] = matrix[0, 0] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2])
-        - matrix[1, 0] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2])
-        + matrix[2, 0] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2])
-    res[2, 3] = -(matrix[0, 0] * (matrix[1, 1] * matrix[2, 3] - matrix[1, 3] * matrix[2, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[2, 3] - matrix[0, 3] * matrix[2, 1])
-        + matrix[2, 0] * (matrix[0, 1] * matrix[1, 3] - matrix[0, 3] * matrix[1, 1]))
-    res[3, 3] = matrix[0, 0] * (matrix[1, 1] * matrix[2, 2] - matrix[1, 2] * matrix[2, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[2, 2] - matrix[0, 2] * matrix[2, 1])
-        + matrix[2, 0] * (matrix[0, 1] * matrix[1, 2] - matrix[0, 2] * matrix[1, 1])
-
-    return res
+    adjugateAndDeterminant(matrix).adjugate
 }
 
 /// Computes the adjugate of a 4x4 matrix.
 /// - Parameter matrix: The input matrix.
 /// - Returns: The adjugate of the input matrix.
 @inlinable
 public func adjugate(_ matrix: Mat4x4d) -> Mat4x4d {
-    var res = Mat4x4d.identity
-
-    res[0, 0] = matrix[1, 1] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 1] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        + matrix[3, 1] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2])
-    res[1, 0] = -(matrix[1, 0] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 0] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        + matrix[3, 0] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2]))
-    res[2, 0] = matrix[1, 0] * (matrix[2, 1] * matrix[3, 3] - matrix[2, 3] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[1, 1] * matrix[3, 3] - matrix[1, 3] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[1, 1] * matrix[2, 3] - matrix[1, 3] * matrix[2, 1])
-    res[3, 0] = -(matrix[1, 0] * (matrix[2, 1] * matrix[3, 2] - matrix[2, 2] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[1, 1] * matrix[3, 2] - matrix[1, 2] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[1, 1] * matrix[2, 2] - matrix[1, 2] * matrix[2, 1]))
-
-    res[0, 1] = -(matrix[0, 1] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 1] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 1] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2]))
-    res[1, 1] = matrix[0, 0] * (matrix[2, 2] * matrix[3, 3] - matrix[2, 3] * matrix[3, 2])
-        - matrix[2, 0] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 0] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2])
-    res[2, 1] = -(matrix[0, 0] * (matrix[2, 1] * matrix[3, 3] - matrix[2, 3] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[0, 1] * matrix[3, 3] - matrix[0, 3] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[2, 3] - matrix[0, 3] * matrix[2, 1]))
-    res[3, 1] = matrix[0, 0] * (matrix[2, 1] * matrix[3, 2] - matrix[2, 2] * matrix[3, 1])
-        - matrix[2, 0] * (matrix[0, 1] * matrix[3, 2] - matrix[0, 2] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[2, 2] - matrix[0, 2] * matrix[2, 1])
-
-    res[0, 2] = matrix[0, 1] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        - matrix[1, 1] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 1] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2])
-    res[1, 2] = -(matrix[0, 0] * (matrix[1, 2] * matrix[3, 3] - matrix[1, 3] * matrix[3, 2])
-        - matrix[1, 0] * (matrix[0, 2] * matrix[3, 3] - matrix[0, 3] * matrix[3, 2])
-        + matrix[3, 0] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2]))
-    res[2, 2] = matrix[0, 0] * (matrix[1, 1] * matrix[3, 3] - matrix[1, 3] * matrix[3, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[3, 3] - matrix[0, 3] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[1, 3] - matrix[0, 3] * matrix[1, 1])
-    res[3, 2] = -(matrix[0, 0] * (matrix[1, 1] * matrix[3, 2] - matrix[1, 2] * matrix[3, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[3, 2] - matrix[0, 2] * matrix[3, 1])
-        + matrix[3, 0] * (matrix[0, 1] * matrix[1, 2] - matrix[0, 2] * matrix[1, 1]))
-
-    res[0, 3] = -(matrix[0, 1] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2])
-        - matrix[1, 1] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2])
-        + matrix[2, 1] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2]))
-    res[1, 3] = matrix[0, 0] * (matrix[1, 2] * matrix[2, 3] - matrix[1, 3] * matrix[2, 2])
-        - matrix[1, 0] * (matrix[0, 2] * matrix[2, 3] - matrix[0, 3] * matrix[2, 2])
-        + matrix[2, 0] * (matrix[0, 2] * matrix[1, 3] - matrix[0, 3] * matrix[1, 2])
-    res[2, 3] = -(matrix[0, 0] * (matrix[1, 1] * matrix[2, 3] - matrix[1, 3] * matrix[2, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[2, 3] - matrix[0, 3] * matrix[2, 1])
-        + matrix[2, 0] * (matrix[0, 1] * matrix[1, 3] - matrix[0, 3] * matrix[1, 1]))
-    res[3, 3] = matrix[0, 0] * (matrix[1, 1] * matrix[2, 2] - matrix[1, 2] * matrix[2, 1])
-        - matrix[1, 0] * (matrix[0, 1] * matrix[2, 2] - matrix[0, 2] * matrix[2, 1])
-        + matrix[2, 0] * (matrix[0, 1] * matrix[1, 2] - matrix[0, 2] * matrix[1, 1])
-
-    return res
+    adjugateAndDeterminant(matrix).adjugate
 }

Sources/FirebladeMath/Functions/determinant.swift

@@ -10,27 +10,7 @@ public func determinant(_ mat: Mat4x4f) -> Float {
     #if FRB_MATH_USE_SIMD
     return simd.simd_determinant(mat.storage)
     #else
-    let subFactor00 = mat[1, 1] * (mat[2, 2] * mat[3, 3] - mat[2, 3] * mat[3, 2])
-    let subFactor01 = mat[1, 2] * (mat[2, 1] * mat[3, 3] - mat[2, 3] * mat[3, 1])
-    let subFactor02 = mat[1, 3] * (mat[2, 1] * mat[3, 2] - mat[2, 2] * mat[3, 1])
-    let det0 = mat[0, 0] * (subFactor00 - subFactor01 + subFactor02)
-
-    let subFactor10 = mat[1, 0] * (mat[2, 2] * mat[3, 3] - mat[2, 3] * mat[3, 2])
-    let subFactor11 = mat[1, 2] * (mat[2, 0] * mat[3, 3] - mat[2, 3] * mat[3, 0])
-    let subFactor12 = mat[1, 3] * (mat[2, 0] * mat[3, 2] - mat[2, 2] * mat[3, 0])
-    let det1 = mat[0, 1] * (subFactor10 - subFactor11 + subFactor12)
-
-    let subFactor20 = mat[1, 0] * (mat[2, 1] * mat[3, 3] - mat[2, 3] * mat[3, 1])
-    let subFactor21 = mat[1, 1] * (mat[2, 0] * mat[3, 3] - mat[2, 3] * mat[3, 0])
-    let subFactor22 = mat[1, 3] * (mat[2, 0] * mat[3, 1] - mat[2, 1] * mat[3, 0])
-    let det2 = mat[0, 2] * (subFactor20 - subFactor21 + subFactor22)
-
-    let subFactor30 = mat[1, 0] * (mat[2, 1] * mat[3, 2] - mat[2, 2] * mat[3, 1])
-    let subFactor31 = mat[1, 1] * (mat[2, 0] * mat[3, 2] - mat[2, 2] * mat[3, 0])
-    let subFactor32 = mat[1, 2] * (mat[2, 0] * mat[3, 1] - mat[2, 1] * mat[3, 0])
-    let det3 = mat[0, 3] * (subFactor30 - subFactor31 + subFactor32)
-
-    return det0 - det1 + det2 - det3
+    return adjugateAndDeterminant(mat).determinant
     #endif
 }
 
@@ -57,27 +37,7 @@ public func determinant(_ mat: Mat4x4d) -> Double {
     #if FRB_MATH_USE_SIMD
     return simd.simd_determinant(mat.storage)
     #else
-    let subFactor00 = mat[1, 1] * (mat[2, 2] * mat[3, 3] - mat[2, 3] * mat[3, 2])
-    let subFactor01 = mat[1, 2] * (mat[2, 1] * mat[3, 3] - mat[2, 3] * mat[3, 1])
-    let subFactor02 = mat[1, 3] * (mat[2, 1] * mat[3, 2] - mat[2, 2] * mat[3, 1])
-    let det0 = mat[0, 0] * (subFactor00 - subFactor01 + subFactor02)
-
-    let subFactor10 = mat[1, 0] * (mat[2, 2] * mat[3, 3] - mat[2, 3] * mat[3, 2])
-    let subFactor11 = mat[1, 2] * (mat[2, 0] * mat[3, 3] - mat[2, 3] * mat[3, 0])
-    let subFactor12 = mat[1, 3] * (mat[2, 0] * mat[3, 2] - mat[2, 2] * mat[3, 0])
-    let det1 = mat[0, 1] * (subFactor10 - subFactor11 + subFactor12)
-
-    let subFactor20 = mat[1, 0] * (mat[2, 1] * mat[3, 3] - mat[2, 3] * mat[3, 1])
-    let subFactor21 = mat[1, 1] * (mat[2, 0] * mat[3, 3] - mat[2, 3] * mat[3, 0])
-    let subFactor22 = mat[1, 3] * (mat[2, 0] * mat[3, 1] - mat[2, 1] * mat[3, 0])
-    let det2 = mat[0, 2] * (subFactor20 - subFactor21 + subFactor22)
-
-    let subFactor30 = mat[1, 0] * (mat[2, 1] * mat[3, 2] - mat[2, 2] * mat[3, 1])
-    let subFactor31 = mat[1, 1] * (mat[2, 0] * mat[3, 2] - mat[2, 2] * mat[3, 0])
-    let subFactor32 = mat[1, 2] * (mat[2, 0] * mat[3, 1] - mat[2, 1] * mat[3, 0])
-    let det3 = mat[0, 3] * (subFactor30 - subFactor31 + subFactor32)
-
-    return det0 - det1 + det2 - det3
+    return adjugateAndDeterminant(mat).determinant
     #endif
 }
 

Sources/FirebladeMath/Functions/inverse.swift

@@ -38,8 +38,7 @@ public func inverse(_ matrix: Mat4x4f) -> Mat4x4f {
     #if FRB_MATH_USE_SIMD
     return Mat4x4f(storage: simd.simd_inverse(matrix.storage))
     #else
-    let adj: Mat4x4f = adjugate(matrix)
-    let det: Float = determinant(matrix)
+    let (adj, det) = adjugateAndDeterminant(matrix)
     return multiply(1.0 / det, adj)
     #endif
 }
@@ -52,8 +51,7 @@ public func inverse(_ matrix: Mat4x4d) -> Mat4x4d {
     #if FRB_MATH_USE_SIMD
     return Mat4x4d(storage: simd.simd_inverse(matrix.storage))
     #else
-    let adj: Mat4x4d = adjugate(matrix)
-    let det: Double = determinant(matrix)
+    let (adj, det) = adjugateAndDeterminant(matrix)
     return multiply(1.0 / det, adj)
     #endif
 }

Sources/FirebladeMath/Matrix/Matrix4x4+Cofactors.swift

@@ -0,0 +1,71 @@
+/// Computes the adjugate and the determinant of a 4x4 matrix in one pass.
+///
+/// Uses the Laplace expansion by 2x2 minors: the six 2x2 sub-determinants of the first two rows
+/// (`s0...s5`) and of the last two rows (`c0...c5`) are computed once and shared by all sixteen
+/// cofactors and the determinant, instead of re-deriving each 3x3 minor from scratch.
+///
+/// - Parameter matrix: The matrix.
+/// - Returns: The adjugate (transposed cofactor matrix) and the determinant of `matrix`.
+@inlinable
+func adjugateAndDeterminant<Storage>(_ matrix: Matrix4x4<Storage>) -> (adjugate: Matrix4x4<Storage>, determinant: Storage.Value) { // swiftlint:disable:this explicit_top_level_acl
+    typealias Value = Storage.Value
+    typealias Vector = Matrix4x4<Storage>.Vector
+
+    // Every intermediate is explicitly typed and kept short: older compilers (Swift 6.1/6.2)
+    // otherwise fail to type-check these generic expressions in reasonable time.
+
+    // Row-major names: aRC is the element in row R, column C.
+    let (col0, col1, col2, col3) = matrix.columns
+    let a00: Value = col0.x, a10: Value = col0.y, a20: Value = col0.z, a30: Value = col0.w
+    let a01: Value = col1.x, a11: Value = col1.y, a21: Value = col1.z, a31: Value = col1.w
+    let a02: Value = col2.x, a12: Value = col2.y, a22: Value = col2.z, a32: Value = col2.w
+    let a03: Value = col3.x, a13: Value = col3.y, a23: Value = col3.z, a33: Value = col3.w
+
+    let s0: Value = a00 * a11 - a10 * a01
+    let s1: Value = a00 * a12 - a10 * a02
+    let s2: Value = a00 * a13 - a10 * a03
+    let s3: Value = a01 * a12 - a11 * a02
+    let s4: Value = a01 * a13 - a11 * a03
+    let s5: Value = a02 * a13 - a12 * a03
+
+    let c5: Value = a22 * a33 - a32 * a23
+    let c4: Value = a21 * a33 - a31 * a23
+    let c3: Value = a21 * a32 - a31 * a22
+    let c2: Value = a20 * a33 - a30 * a23
+    let c1: Value = a20 * a32 - a30 * a22
+    let c0: Value = a20 * a31 - a30 * a21
+
+    let determinant01: Value = s0 * c5 - s1 * c4
+    let determinant23: Value = s2 * c3 + s3 * c2
+    let determinant45: Value = s5 * c0 - s4 * c1
+    let determinant: Value = determinant01 + determinant23 + determinant45
+
+    // bRC is the adjugate element in row R, column C; the matrix is built from its columns.
+    let b00: Value = a11 * c5 - a12 * c4 + a13 * c3
+    let b10: Value = -a10 * c5 + a12 * c2 - a13 * c1
+    let b20: Value = a10 * c4 - a11 * c2 + a13 * c0
+    let b30: Value = -a10 * c3 + a11 * c1 - a12 * c0
+
+    let b01: Value = -a01 * c5 + a02 * c4 - a03 * c3
+    let b11: Value = a00 * c5 - a02 * c2 + a03 * c1
+    let b21: Value = -a00 * c4 + a01 * c2 - a03 * c0
+    let b31: Value = a00 * c3 - a01 * c1 + a02 * c0
+
+    let b02: Value = a31 * s5 - a32 * s4 + a33 * s3
+    let b12: Value = -a30 * s5 + a32 * s2 - a33 * s1
+    let b22: Value = a30 * s4 - a31 * s2 + a33 * s0
+    let b32: Value = -a30 * s3 + a31 * s1 - a32 * s0
+
+    let b03: Value = -a21 * s5 + a22 * s4 - a23 * s3
+    let b13: Value = a20 * s5 - a22 * s2 + a23 * s1
+    let b23: Value = -a20 * s4 + a21 * s2 - a23 * s0
+    let b33: Value = a20 * s3 - a21 * s1 + a22 * s0
+
+    let adjugate = Matrix4x4<Storage>(
+        Vector(b00, b10, b20, b30),
+        Vector(b01, b11, b21, b31),
+        Vector(b02, b12, b22, b32),
+        Vector(b03, b13, b23, b33)
+    )
+    return (adjugate, determinant)
+}

Sources/FirebladeMath/Quat/Quaternion+Multiplication.swift

@@ -33,11 +33,16 @@ public func multiply(_ lhs: Quat4f, _ rhs: Quat4f) -> Quat4f {
     #if FRB_MATH_USE_SIMD
     return Quat4f(storage: simd.simd_mul(lhs.storage, rhs.storage))
     #else
-    let n0 = lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z
-    let n1 = lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y
-    let n2 = lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x
-    let n3 = lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w
-    return Quat4f(n1, n2, n3, n0)
+    // Hamilton product as four 4-wide multiply-adds over the [x, y, z, w] lanes.
+    // Terms are explicitly typed so older compilers type-check them in reasonable time.
+    let lv: Vec4f = lhs.storage.storage
+    let rv: Vec4f = rhs.storage.storage
+    let term0: Vec4f = Vec4f(repeating: lv.w) * rv
+    let term1: Vec4f = Vec4f(lv.x, lv.y, lv.z, lv.x) * Vec4f(rv.w, rv.w, rv.w, -rv.x)
+    let term2: Vec4f = Vec4f(lv.y, lv.z, lv.x, lv.y) * Vec4f(rv.z, rv.x, rv.y, -rv.y)
+    let term3: Vec4f = Vec4f(lv.z, lv.x, lv.y, lv.z) * Vec4f(rv.y, rv.z, rv.x, rv.z)
+    let product: Vec4f = (term0 + term1) + (term2 - term3)
+    return Quat4f(product)
     #endif
 }
 
@@ -135,11 +140,16 @@ public func multiply(_ lhs: Quat4d, _ rhs: Quat4d) -> Quat4d {
     #if FRB_MATH_USE_SIMD
     return Quat4d(storage: simd.simd_mul(lhs.storage, rhs.storage))
     #else
-    let n0 = lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z
-    let n1 = lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y
-    let n2 = lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x
-    let n3 = lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w
-    return Quat4d(n1, n2, n3, n0)
+    // Hamilton product as four 4-wide multiply-adds over the [x, y, z, w] lanes.
+    // Terms are explicitly typed so older compilers type-check them in reasonable time.
+    let lv: Vec4d = lhs.storage.storage
+    let rv: Vec4d = rhs.storage.storage
+    let term0: Vec4d = Vec4d(repeating: lv.w) * rv
+    let term1: Vec4d = Vec4d(lv.x, lv.y, lv.z, lv.x) * Vec4d(rv.w, rv.w, rv.w, -rv.x)
+    let term2: Vec4d = Vec4d(lv.y, lv.z, lv.x, lv.y) * Vec4d(rv.z, rv.x, rv.y, -rv.y)
+    let term3: Vec4d = Vec4d(lv.z, lv.x, lv.y, lv.z) * Vec4d(rv.y, rv.z, rv.x, rv.z)
+    let product: Vec4d = (term0 + term1) + (term2 - term3)
+    return Quat4d(product)
     #endif
 }