From 09973e78385160737de724969f7cbd54b39f775e Mon Sep 17 00:00:00 2001 From: Christian Treffs Date: Sat, 17 Oct 2026 20:47:02 +0200 Subject: [PATCH 1/2] Make the portable implementations match Apple's simd on edge cases MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The differential tests against Apple's simd found five places where the same call returned different results with and without the simd trait, i.e. on Linux/Windows vs. Apple platforms: - normalize(Quat4f/Quat4d) of a zero-length quaternion: simd returns the identity, portable returned NaN (and so did act(_:_:), which normalizes). Portable now returns the identity when the squared length is zero, like simd. - clamp(NaN, min, max): simd returns min (IEEE maxNum/minNum), portable returned NaN. Now uses Float/Double.maximum/minimum. - fract(NaN/±inf): simd returns the largest value below 1, portable returned NaN. Now uses .minimum like simd_fract. - quaternion(from:to:): portable used only the "reduced" formula, which simd documents as unstable beyond 90 degrees and which returned NaN for opposite vectors. Now a port of simd's algorithm: two-step rotation beyond 90 degrees and the same orthogonal axis for opposite vectors. - matrix3x3/matrix4x4(from:) diagonal: 1 - 2(y^2 + z^2) vs simd's x^2 - y^2 - z^2 + w^2. Equal for unit quaternions, different for non-unit ones. Now uses simd's formula. The simd code path is unchanged. Claude-Session: https://claude.ai/code/session_011WSUtoLptjea3s3iC5wTNC --- Sources/FirebladeMath/Functions/clamp.swift | 6 +- Sources/FirebladeMath/Functions/fract.swift | 6 +- Sources/FirebladeMath/Functions/matrix.swift | 36 ++++++---- .../FirebladeMath/Functions/normalize.swift | 20 ++++-- .../FirebladeMath/Functions/quaternion.swift | 70 ++++++++++++++++--- 5 files changed, 110 insertions(+), 28 deletions(-) diff --git a/Sources/FirebladeMath/Functions/clamp.swift b/Sources/FirebladeMath/Functions/clamp.swift index 7fb8ed2..2980e1d 100644 --- a/Sources/FirebladeMath/Functions/clamp.swift +++ b/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 } diff --git a/Sources/FirebladeMath/Functions/fract.swift b/Sources/FirebladeMath/Functions/fract.swift index 34b4ade..01c410c 100644 --- a/Sources/FirebladeMath/Functions/fract.swift +++ b/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 } diff --git a/Sources/FirebladeMath/Functions/matrix.swift b/Sources/FirebladeMath/Functions/matrix.swift index 9c9e2c6..569f5f5 100644 --- a/Sources/FirebladeMath/Functions/matrix.swift +++ b/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) diff --git a/Sources/FirebladeMath/Functions/normalize.swift b/Sources/FirebladeMath/Functions/normalize.swift index 659e205..67ddb0f 100644 --- a/Sources/FirebladeMath/Functions/normalize.swift +++ b/Sources/FirebladeMath/Functions/normalize.swift @@ -88,28 +88,40 @@ public func normalize(_ x: SIMD4) -> SIMD4 { /// 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 } diff --git a/Sources/FirebladeMath/Functions/quaternion.swift b/Sources/FirebladeMath/Functions/quaternion.swift index 8f8682a..c59f9b8 100644 --- a/Sources/FirebladeMath/Functions/quaternion.swift +++ b/Sources/FirebladeMath/Functions/quaternion.swift @@ -47,10 +47,37 @@ public func quaternion(from: SIMD3, to: SIMD3) -> Quat4f { #if FRB_MATH_USE_SIMD return Quat4f(storage: simd.simd_quaternion(from, to)) #else - let half = normalize(from + to) - let img = cross(from, half) - let real = dot(from, half) - return Quat4f(img.x, img.y, img.z, real) + // Mirrors Apple's simd_quaternion(from:to:), so every platform returns the same rotation. + func reduced(_ from: SIMD3, _ to: SIMD3) -> Quat4f { + // Accurate while `from` and `to` are less than 90 degrees apart. + let half: SIMD3 = normalize(from + to) + let imaginary: SIMD3 = cross(from, half) + return Quat4f(imaginary.x, imaginary.y, imaginary.z, dot(from, half)) + } + + if dot(from, to) >= 0 { + return reduced(from, to) + } + + let sum: SIMD3 = normalize(from) + normalize(to) + if dot(sum, sum) <= 0x1p-46 { + // Nearly opposite: the axis is underspecified, so pick one orthogonal to `from` + // the same way simd does, and rotate by pi. + let absolute = SIMD3(Swift.abs(from.x), Swift.abs(from.y), Swift.abs(from.z)) + let reference: SIMD3 = if absolute.x <= absolute.y, absolute.x <= absolute.z { + SIMD3(1, 0, 0) + } else if absolute.y <= absolute.z { + SIMD3(0, 1, 0) + } else { + SIMD3(0, 0, 1) + } + let axis: SIMD3 = normalize(cross(from, reference)) + return Quat4f(axis.x, axis.y, axis.z, 0) + } + + // More than 90 degrees apart: rotate in two steps through the halfway vector. + let half: SIMD3 = normalize(sum) + return multiply(reduced(from, half), reduced(half, to)) #endif } @@ -69,10 +96,37 @@ public func quaternion(from: SIMD3, to: SIMD3) -> Quat4d { #if FRB_MATH_USE_SIMD return Quat4d(storage: simd.simd_quaternion(from, to)) #else - let half = normalize(from + to) - let img = cross(from, half) - let real = dot(from, half) - return Quat4d(img.x, img.y, img.z, real) + // Mirrors Apple's simd_quaternion(from:to:), so every platform returns the same rotation. + func reduced(_ from: SIMD3, _ to: SIMD3) -> Quat4d { + // Accurate while `from` and `to` are less than 90 degrees apart. + let half: SIMD3 = normalize(from + to) + let imaginary: SIMD3 = cross(from, half) + return Quat4d(imaginary.x, imaginary.y, imaginary.z, dot(from, half)) + } + + if dot(from, to) >= 0 { + return reduced(from, to) + } + + let sum: SIMD3 = normalize(from) + normalize(to) + if dot(sum, sum) <= 0x1p-104 { + // Nearly opposite: the axis is underspecified, so pick one orthogonal to `from` + // the same way simd does, and rotate by pi. + let absolute = SIMD3(Swift.abs(from.x), Swift.abs(from.y), Swift.abs(from.z)) + let reference: SIMD3 = if absolute.x <= absolute.y, absolute.x <= absolute.z { + SIMD3(1, 0, 0) + } else if absolute.y <= absolute.z { + SIMD3(0, 1, 0) + } else { + SIMD3(0, 0, 1) + } + let axis: SIMD3 = normalize(cross(from, reference)) + return Quat4d(axis.x, axis.y, axis.z, 0) + } + + // More than 90 degrees apart: rotate in two steps through the halfway vector. + let half: SIMD3 = normalize(sum) + return multiply(reduced(from, half), reduced(half, to)) #endif } From 041e41ad4008622b48bbe6d3e55b5d29677b444f Mon Sep 17 00:00:00 2001 From: Christian Treffs Date: Sat, 17 Oct 2026 20:47:02 +0200 Subject: [PATCH 2/2] Compare every simd-backed function against Apple's simd Extends the PortableVsSIMD differential tests from the functions changed in #76 to all 108 functions that have both a simd and a portable implementation: matrix multiplication, 3x3 determinant/inverse, transpose, quaternion arithmetic and queries, act, vector dot/length/ normalize/cross, clamp/fract/mix/rsqrt/sign/smoothstep/step, and the quaternion <-> matrix / angle-axis / from-to conversions. Inputs are 256 pseudo-random values per test plus edge cases: identity, diagonal, permutation and singular matrices, small and large scales, zero vectors and quaternions, unit axes, 180-degree rotations, opposite vectors, NaN, infinities and signed zeros. Shared inputs and helpers live in PortableVsSIMDSupport.swift. Every test names the functions it covers with `// covers:` comments. Scripts/check-simd-differential-coverage.py fails if a function with a simd branch is not covered or an annotation is stale; it runs in the macOS CI job (make check-simd-coverage). make test-portable runs the differential tests locally. Claude-Session: https://claude.ai/code/session_011WSUtoLptjea3s3iC5wTNC --- .github/workflows/ci-macos.yml | 3 + Makefile | 9 +- Scripts/check-simd-differential-coverage.py | 112 ++++++ .../PortableVsSIMDCofactorTests.swift} | 68 +--- .../PortableVsSIMDConversionTests.swift | 132 +++++++ .../PortableVsSIMDMatrixTests.swift | 346 ++++++++++++++++++ .../PortableVsSIMDQuaternionTests.swift | 135 +++++++ .../PortableVsSIMDScalarTests.swift | 148 ++++++++ .../PortableVsSIMDSupport.swift | 196 ++++++++++ .../PortableVsSIMDVectorTests.swift | 101 +++++ 10 files changed, 1197 insertions(+), 53 deletions(-) create mode 100755 Scripts/check-simd-differential-coverage.py rename Tests/FirebladeMathTests/{PortableVsSIMDTests.swift => PortableVsSIMD/PortableVsSIMDCofactorTests.swift} (75%) create mode 100644 Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDConversionTests.swift create mode 100644 Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDMatrixTests.swift create mode 100644 Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDQuaternionTests.swift create mode 100644 Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDScalarTests.swift create mode 100644 Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDSupport.swift create mode 100644 Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDVectorTests.swift diff --git a/.github/workflows/ci-macos.yml b/.github/workflows/ci-macos.yml index a23b79c..198f2ab 100644 --- a/.github/workflows/ci-macos.yml +++ b/.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 diff --git a/Makefile b/Makefile index 31de81b..8eaa311 100644 --- a/Makefile +++ b/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 diff --git a/Scripts/check-simd-differential-coverage.py b/Scripts/check-simd-differential-coverage.py new file mode 100755 index 0000000..6668404 --- /dev/null +++ b/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) + +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(?:[ \t]*@[^\n]*\n)*)[ \t]*public\s+(?:static\s+)?func\s+(?P\w+)\s*\((?P[^)]*)\)[^{]*\{", + 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()) diff --git a/Tests/FirebladeMathTests/PortableVsSIMDTests.swift b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDCofactorTests.swift similarity index 75% rename from Tests/FirebladeMathTests/PortableVsSIMDTests.swift rename to Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDCofactorTests.swift index eaeb83f..957e62a 100644 --- a/Tests/FirebladeMathTests/PortableVsSIMDTests.swift +++ b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDCofactorTests.swift @@ -1,16 +1,10 @@ // -// PortableVsSIMDTests.swift +// PortableVsSIMDCofactorTests.swift // FirebladeMathTests // -// Differential tests: FirebladeMath's portable (non-simd) implementations against Apple's -// simd library, which is exactly what FirebladeMath forwards to when the `simd` trait is on. -// -// Only compiled when Apple's simd is available but the `simd` trait is disabled, i.e. -// `swift test --disable-default-traits` on Apple platforms (see the macOS CI workflow). -// -// Covers the functions with portable implementations changed in #76: -// 4x4 inverse, determinant and adjugate, and quaternion multiplication. +// Differential tests for the 4x4 inverse, determinant, adjugate and quaternion multiplication +// against Apple's simd. See PortableVsSIMDSupport.swift. // // Tolerances are about 10x the worst deviation measured over all inputs; both sides differ // only by floating-point rounding. @@ -20,7 +14,7 @@ import FirebladeMath import simd import Testing -struct PortableVsSIMDTests { +struct PortableVsSIMDCofactorTests { // MARK: - Inputs /// Column-major 4x4 element lists: pseudo-random matrices plus structured cases. @@ -62,34 +56,38 @@ struct PortableVsSIMDTests { // MARK: - 4x4 determinant + // covers: determinant(Mat4x4d) @Test(arguments: matrices + singularMatrices) func determinantDouble(_ values: [Double]) { let portable = FirebladeMath.determinant(Mat4x4d(values)) - let reference = simd_determinant(simd_double4x4(values)) + let reference = simd_determinant(simd_double4x4(elements: values)) #expect(isClose(portable, reference, scale: magnitude(values) ^^ 4, tolerance: 1e-14)) } + // covers: determinant(Mat4x4f) @Test(arguments: matrices + singularMatrices) func determinantFloat(_ values: [Double]) { let floats = values.map(Float.init) let portable = FirebladeMath.determinant(Mat4x4f(floats)) - let reference = simd_determinant(simd_float4x4(floats)) + let reference = simd_determinant(simd_float4x4(elements: floats)) #expect(isClose(Double(portable), Double(reference), scale: magnitude(values) ^^ 4, tolerance: 3e-6)) } // MARK: - 4x4 inverse + // covers: inverse(Mat4x4d) @Test(arguments: matrices) func inverseDouble(_ values: [Double]) { - let reference = simd_inverse(simd_double4x4(values)) + let reference = simd_inverse(simd_double4x4(elements: values)) let portable = FirebladeMath.inverse(Mat4x4d(values)) #expect(isClose(portable.elements, reference.elements, tolerance: 5e-13)) } + // covers: inverse(Mat4x4f) @Test(arguments: matrices) func inverseFloat(_ values: [Double]) { let floats = values.map(Float.init) - let reference = simd_inverse(simd_float4x4(floats)) + let reference = simd_inverse(simd_float4x4(elements: floats)) let portable = FirebladeMath.inverse(Mat4x4f(floats)) #expect(isClose(portable.elements.map(Double.init), reference.elements.map(Double.init), tolerance: 2e-4)) } @@ -98,7 +96,7 @@ struct PortableVsSIMDTests { @Test(arguments: matrices) func adjugateDouble(_ values: [Double]) { - let matrix = simd_double4x4(values) + let matrix = simd_double4x4(elements: values) let reference = simd_inverse(matrix) * simd_determinant(matrix) let portable = FirebladeMath.adjugate(Mat4x4d(values)) #expect(isClose(portable.elements, reference.elements, tolerance: 5e-13)) @@ -107,7 +105,7 @@ struct PortableVsSIMDTests { @Test(arguments: matrices) func adjugateFloat(_ values: [Double]) { let floats = values.map(Float.init) - let matrix = simd_float4x4(floats) + let matrix = simd_float4x4(elements: floats) let reference = simd_inverse(matrix) * simd_determinant(matrix) let portable = FirebladeMath.adjugate(Mat4x4f(floats)) #expect(isClose(portable.elements.map(Double.init), reference.elements.map(Double.init), tolerance: 2e-4)) @@ -115,6 +113,7 @@ struct PortableVsSIMDTests { // MARK: - Quaternion multiplication + // covers: multiply(Quat4d, Quat4d) @Test(arguments: quaternionPairs) func quaternionMultiplyDouble(_ pair: ([Double], [Double])) { let (lhs, rhs) = pair @@ -123,6 +122,7 @@ struct PortableVsSIMDTests { #expect(isClose(portable.elements, Array(reference), tolerance: 1e-14)) } + // covers: multiply(Quat4f, Quat4f) @Test(arguments: quaternionPairs) func quaternionMultiplyFloat(_ pair: ([Double], [Double])) { let (lhs, rhs) = pair @@ -136,22 +136,6 @@ struct PortableVsSIMDTests { // MARK: - Helpers -private struct SplitMix { - var state: UInt64 - - init(seed: UInt64) { state = seed } - - mutating func next(in range: ClosedRange) -> Double { - state &+= 0x9E37_79B9_7F4A_7C15 - var z = state - z = (z ^ (z >> 30)) &* 0xBF58_476D_1CE4_E5B9 - z = (z ^ (z >> 27)) &* 0x94D0_49BB_1331_11EB - z ^= z >> 31 - let unit = Double(z >> 11) / Double(1 << 53) - return range.lowerBound + unit * (range.upperBound - range.lowerBound) - } -} - private func rotationAndTranslation(angle: Double, axis: (Double, Double, Double), translation: (Double, Double, Double)) -> [Double] { let rotation = simd_matrix4x4(simd_quatd(angle: angle, axis: simd_normalize(SIMD3(axis.0, axis.1, axis.2)))) var values = rotation.elements @@ -161,26 +145,6 @@ private func rotationAndTranslation(angle: Double, axis: (Double, Double, Double return values } -private extension simd_double4x4 { - init(_ values: [Double]) { - self.init(SIMD4(values[0..<4]), SIMD4(values[4..<8]), SIMD4(values[8..<12]), SIMD4(values[12..<16])) - } - - var elements: [Double] { - [columns.0, columns.1, columns.2, columns.3].flatMap { [$0.x, $0.y, $0.z, $0.w] } - } -} - -private extension simd_float4x4 { - init(_ values: [Float]) { - self.init(SIMD4(values[0..<4]), SIMD4(values[4..<8]), SIMD4(values[8..<12]), SIMD4(values[12..<16])) - } - - var elements: [Float] { - [columns.0, columns.1, columns.2, columns.3].flatMap { [$0.x, $0.y, $0.z, $0.w] } - } -} - private func magnitude(_ values: [Double]) -> Double { Swift.max(1, values.map { Swift.abs($0) }.max() ?? 1) } diff --git a/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDConversionTests.swift b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDConversionTests.swift new file mode 100644 index 0000000..fbedc3a --- /dev/null +++ b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDConversionTests.swift @@ -0,0 +1,132 @@ +// +// PortableVsSIMDConversionTests.swift +// FirebladeMathTests +// + +// Differential tests for quaternion <-> matrix conversions and quaternion construction +// against Apple's simd. See PortableVsSIMDSupport.swift. + +#if canImport(simd) && !FRB_MATH_USE_SIMD +import FirebladeMath +import simd +import Testing + +struct PortableVsSIMDConversionTests { + static let rotations = DifferentialInputs.rotations(seed: 0x5001) + + + // MARK: - Float + + // covers: matrix3x3(from: Quat4f) + // covers: matrix4x4(from: Quat4f) + @Test func matrixFromQuaternionF() { + // Includes non-unit quaternions: both sides must agree on how they are treated. + for (index, values) in (Self.rotations + DifferentialInputs.quaternions(seed: 0x5002)).enumerated() { + let vector = SIMD4(values.floats) + let reference = simd_quatf(vector: vector) + #expect(differentialMatch(FirebladeMath.matrix3x3(from: Quat4f(vector)).elements, simd_matrix3x3(reference).elements, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.matrix4x4(from: Quat4f(vector)).elements, simd_matrix4x4(reference).elements, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: quaternion(matrix: Mat3x3f) + // covers: quaternion(matrix: Mat4x4f) + @Test func quaternionFromMatrixF() { + for (index, values) in Self.rotations.enumerated() { + let rotation = simd_quatf(vector: SIMD4(values.floats)) + let matrix3 = simd_matrix3x3(rotation), matrix4 = simd_matrix4x4(rotation) + let portable3 = FirebladeMath.quaternion(matrix: Mat3x3f(matrix3.elements)) + let portable4 = FirebladeMath.quaternion(matrix: Mat4x4f(matrix4.elements)) + #expect(differentialMatch(portable3.elements, simd_quaternion(matrix3).vector.lanes, tolerance: Tolerance.floatLoose), "input \(index)") + #expect(differentialMatch(portable4.elements, simd_quaternion(matrix4).vector.lanes, tolerance: Tolerance.floatLoose), "input \(index)") + } + } + + // covers: quaternion(angle: Float, axis: SIMD3) + @Test func quaternionFromAngleAxisF() { + let axes = DifferentialInputs.vectors(dimension: 3, seed: 0x5003) + let angles = DifferentialInputs.scalars(seed: 0x5004) + for (index, values) in axes.enumerated() { + let axis = SIMD3(values.floats) + let angle = Float(angles[index % angles.count]) + let portable = FirebladeMath.quaternion(angle: angle, axis: axis) + let reference = simd_quaternion(angle, axis) + #expect(differentialMatch(portable.elements, reference.vector.lanes, tolerance: Tolerance.floatLoose), "input \(index): angle \(angle), axis \(axis)") + } + } + + // covers: quaternion(from: SIMD3, to: SIMD3) + @Test func quaternionFromToF() { + // simd_quaternion(from:to:) expects unit vectors. + let directions = DifferentialInputs.vectors(dimension: 3, seed: 0x5005) + .filter { $0.contains { $0 != 0 } } + .map { simd_normalize(SIMD3($0)) } + for (index, from) in directions.enumerated() { + for to in [directions[(index + 1) % directions.count], from, -from] { + let fromVector = SIMD3(Float(from.x), Float(from.y), Float(from.z)) + let toVector = SIMD3(Float(to.x), Float(to.y), Float(to.z)) + let portable = FirebladeMath.quaternion(from: fromVector, to: toVector) + let reference = simd_quaternion(fromVector, toVector) + #expect(differentialMatch(portable.elements, reference.vector.lanes, tolerance: Tolerance.floatLoose), "input \(index): from \(fromVector) to \(toVector)") + } + } + } + + // MARK: - Double + + // covers: matrix3x3(from: Quat4d) + // covers: matrix4x4(from: Quat4d) + @Test func matrixFromQuaternionD() { + // Includes non-unit quaternions: both sides must agree on how they are treated. + for (index, values) in (Self.rotations + DifferentialInputs.quaternions(seed: 0x5002)).enumerated() { + let vector = SIMD4(values) + let reference = simd_quatd(vector: vector) + #expect(differentialMatch(FirebladeMath.matrix3x3(from: Quat4d(vector)).elements, simd_matrix3x3(reference).elements, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.matrix4x4(from: Quat4d(vector)).elements, simd_matrix4x4(reference).elements, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: quaternion(matrix: Mat3x3d) + // covers: quaternion(matrix: Mat4x4d) + @Test func quaternionFromMatrixD() { + for (index, values) in Self.rotations.enumerated() { + let rotation = simd_quatd(vector: SIMD4(values)) + let matrix3 = simd_matrix3x3(rotation), matrix4 = simd_matrix4x4(rotation) + let portable3 = FirebladeMath.quaternion(matrix: Mat3x3d(matrix3.elements)) + let portable4 = FirebladeMath.quaternion(matrix: Mat4x4d(matrix4.elements)) + #expect(differentialMatch(portable3.elements, simd_quaternion(matrix3).vector.lanes, tolerance: Tolerance.doubleLoose), "input \(index)") + #expect(differentialMatch(portable4.elements, simd_quaternion(matrix4).vector.lanes, tolerance: Tolerance.doubleLoose), "input \(index)") + } + } + + // covers: quaternion(angle: Double, axis: SIMD3) + @Test func quaternionFromAngleAxisD() { + let axes = DifferentialInputs.vectors(dimension: 3, seed: 0x5003) + let angles = DifferentialInputs.scalars(seed: 0x5004) + for (index, values) in axes.enumerated() { + let axis = SIMD3(values) + let angle = Double(angles[index % angles.count]) + let portable = FirebladeMath.quaternion(angle: angle, axis: axis) + let reference = simd_quaternion(angle, axis) + #expect(differentialMatch(portable.elements, reference.vector.lanes, tolerance: Tolerance.doubleLoose), "input \(index): angle \(angle), axis \(axis)") + } + } + + // covers: quaternion(from: SIMD3, to: SIMD3) + @Test func quaternionFromToD() { + // simd_quaternion(from:to:) expects unit vectors. + let directions = DifferentialInputs.vectors(dimension: 3, seed: 0x5005) + .filter { $0.contains { $0 != 0 } } + .map { simd_normalize(SIMD3($0)) } + for (index, from) in directions.enumerated() { + for to in [directions[(index + 1) % directions.count], from, -from] { + let fromVector = SIMD3(Double(from.x), Double(from.y), Double(from.z)) + let toVector = SIMD3(Double(to.x), Double(to.y), Double(to.z)) + let portable = FirebladeMath.quaternion(from: fromVector, to: toVector) + let reference = simd_quaternion(fromVector, toVector) + #expect(differentialMatch(portable.elements, reference.vector.lanes, tolerance: Tolerance.doubleLoose), "input \(index): from \(fromVector) to \(toVector)") + } + } + } +} +#endif diff --git a/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDMatrixTests.swift b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDMatrixTests.swift new file mode 100644 index 0000000..a2ddcda --- /dev/null +++ b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDMatrixTests.swift @@ -0,0 +1,346 @@ +// +// PortableVsSIMDMatrixTests.swift +// FirebladeMathTests +// + +// Differential tests for matrix multiplication, 3x3 determinant/inverse and transpose +// against Apple's simd. See PortableVsSIMDSupport.swift. + +#if canImport(simd) && !FRB_MATH_USE_SIMD +import FirebladeMath +import simd +import Testing + +struct PortableVsSIMDMatrixTests { + + // MARK: - Mat2x2f + + // covers: multiply(Mat2x2f, Mat2x2f) + @Test func multiplyMatrixMatrix2x2f() { + let inputs = DifferentialInputs.matrices(dimension: 2, seed: 0x1001) + for (index, values) in inputs.enumerated() { + let other = inputs[(index + 1) % inputs.count].floats + let lhs = values.floats + let portable = FirebladeMath.multiply(Mat2x2f(lhs), Mat2x2f(other)) + let reference = simd_mul(simd_float2x2(elements: lhs), simd_float2x2(elements: other)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: multiply(Float, Mat2x2f) + @Test func multiplyScalarMatrix2x2f() { + let scalars = DifferentialInputs.scalars(seed: 0x1201) + for (index, values) in DifferentialInputs.matrices(dimension: 2, seed: 0x1001).enumerated() { + let scalar = Float(scalars[index % scalars.count]) + let matrix = values.floats + let portable = FirebladeMath.multiply(scalar, Mat2x2f(matrix)) + let reference = simd_mul(scalar, simd_float2x2(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: multiply(Vec2f, Mat2x2f) + // covers: multiply(Mat2x2f, Vec2f) + @Test func multiplyVectorMatrix2x2f() { + let vectors = DifferentialInputs.vectors(dimension: 2, seed: 0x1101) + for (index, values) in DifferentialInputs.matrices(dimension: 2, seed: 0x1001).enumerated() { + let vector = SIMD2(vectors[index % vectors.count].floats) + let matrix = values.floats + let rowTimesMatrix = FirebladeMath.multiply(vector, Mat2x2f(matrix)) + let matrixTimesColumn = FirebladeMath.multiply(Mat2x2f(matrix), vector) + #expect(differentialMatch(rowTimesMatrix.lanes, simd_mul(vector, simd_float2x2(elements: matrix)).lanes, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(matrixTimesColumn.lanes, simd_mul(simd_float2x2(elements: matrix), vector).lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // MARK: - Mat2x2d + + // covers: multiply(Mat2x2d, Mat2x2d) + @Test func multiplyMatrixMatrix2x2d() { + let inputs = DifferentialInputs.matrices(dimension: 2, seed: 0x1002) + for (index, values) in inputs.enumerated() { + let other = inputs[(index + 1) % inputs.count] + let lhs = values + let portable = FirebladeMath.multiply(Mat2x2d(lhs), Mat2x2d(other)) + let reference = simd_mul(simd_double2x2(elements: lhs), simd_double2x2(elements: other)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: multiply(Double, Mat2x2d) + @Test func multiplyScalarMatrix2x2d() { + let scalars = DifferentialInputs.scalars(seed: 0x1202) + for (index, values) in DifferentialInputs.matrices(dimension: 2, seed: 0x1002).enumerated() { + let scalar = Double(scalars[index % scalars.count]) + let matrix = values + let portable = FirebladeMath.multiply(scalar, Mat2x2d(matrix)) + let reference = simd_mul(scalar, simd_double2x2(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: multiply(Vec2d, Mat2x2d) + // covers: multiply(Mat2x2d, Vec2d) + @Test func multiplyVectorMatrix2x2d() { + let vectors = DifferentialInputs.vectors(dimension: 2, seed: 0x1102) + for (index, values) in DifferentialInputs.matrices(dimension: 2, seed: 0x1002).enumerated() { + let vector = SIMD2(vectors[index % vectors.count]) + let matrix = values + let rowTimesMatrix = FirebladeMath.multiply(vector, Mat2x2d(matrix)) + let matrixTimesColumn = FirebladeMath.multiply(Mat2x2d(matrix), vector) + #expect(differentialMatch(rowTimesMatrix.lanes, simd_mul(vector, simd_double2x2(elements: matrix)).lanes, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(matrixTimesColumn.lanes, simd_mul(simd_double2x2(elements: matrix), vector).lanes, tolerance: Tolerance.double), "input \(index)") + } + } + + // MARK: - Mat3x3f + + // covers: multiply(Mat3x3f, Mat3x3f) + @Test func multiplyMatrixMatrix3x3f() { + let inputs = DifferentialInputs.matrices(dimension: 3, seed: 0x1003) + for (index, values) in inputs.enumerated() { + let other = inputs[(index + 1) % inputs.count].floats + let lhs = values.floats + let portable = FirebladeMath.multiply(Mat3x3f(lhs), Mat3x3f(other)) + let reference = simd_mul(simd_float3x3(elements: lhs), simd_float3x3(elements: other)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: multiply(Float, Mat3x3f) + @Test func multiplyScalarMatrix3x3f() { + let scalars = DifferentialInputs.scalars(seed: 0x1203) + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1003).enumerated() { + let scalar = Float(scalars[index % scalars.count]) + let matrix = values.floats + let portable = FirebladeMath.multiply(scalar, Mat3x3f(matrix)) + let reference = simd_mul(scalar, simd_float3x3(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: multiply(Vec3f, Mat3x3f) + // covers: multiply(Mat3x3f, Vec3f) + @Test func multiplyVectorMatrix3x3f() { + let vectors = DifferentialInputs.vectors(dimension: 3, seed: 0x1103) + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1003).enumerated() { + let vector = SIMD3(vectors[index % vectors.count].floats) + let matrix = values.floats + let rowTimesMatrix = FirebladeMath.multiply(vector, Mat3x3f(matrix)) + let matrixTimesColumn = FirebladeMath.multiply(Mat3x3f(matrix), vector) + #expect(differentialMatch(rowTimesMatrix.lanes, simd_mul(vector, simd_float3x3(elements: matrix)).lanes, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(matrixTimesColumn.lanes, simd_mul(simd_float3x3(elements: matrix), vector).lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: transpose(Mat3x3f) + @Test func transpose3x3f() { + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1003).enumerated() { + let matrix = values.floats + let portable = FirebladeMath.transpose(Mat3x3f(matrix)) + let reference = simd_transpose(simd_float3x3(elements: matrix)) + #expect(portable.elements == reference.elements, "input \(index)") + } + } + + // covers: determinant(Mat3x3f) + @Test func determinant3x3f() { + let singular: [[Double]] = [[Double](repeating: 0, count: 9), [1, 2, 3, 2, 4, 6, 0, 1, 5], [1, 2, 3, 4, 5, 6, 7, 8, 9]] + for (index, values) in (DifferentialInputs.matrices(dimension: 3, seed: 0x1003) + singular).enumerated() { + let matrix = values.floats + let scale = Swift.max(1, values.map { Swift.abs($0) }.max() ?? 1) + let portable = FirebladeMath.determinant(Mat3x3f(matrix)) + let reference = simd_determinant(simd_float3x3(elements: matrix)) + #expect(differentialMatch(Double(portable), Double(reference), tolerance: Tolerance.float, scale: scale * scale * scale), "input \(index)") + } + } + + // covers: inverse(Mat3x3f) + @Test func inverse3x3f() { + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1003).enumerated() { + let matrix = values.floats + let scale = Swift.max(1, values.map { Swift.abs($0) }.max() ?? 1) + // Skip near-singular inputs: both sides would only amplify rounding noise. + guard Swift.abs(Double(simd_determinant(simd_float3x3(elements: matrix)))) > 1e-2 * scale * scale * scale else { continue } + let portable = FirebladeMath.inverse(Mat3x3f(matrix)) + let reference = simd_inverse(simd_float3x3(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.floatLoose), "input \(index)") + } + } + + // MARK: - Mat3x3d + + // covers: multiply(Mat3x3d, Mat3x3d) + @Test func multiplyMatrixMatrix3x3d() { + let inputs = DifferentialInputs.matrices(dimension: 3, seed: 0x1004) + for (index, values) in inputs.enumerated() { + let other = inputs[(index + 1) % inputs.count] + let lhs = values + let portable = FirebladeMath.multiply(Mat3x3d(lhs), Mat3x3d(other)) + let reference = simd_mul(simd_double3x3(elements: lhs), simd_double3x3(elements: other)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: multiply(Double, Mat3x3d) + @Test func multiplyScalarMatrix3x3d() { + let scalars = DifferentialInputs.scalars(seed: 0x1204) + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1004).enumerated() { + let scalar = Double(scalars[index % scalars.count]) + let matrix = values + let portable = FirebladeMath.multiply(scalar, Mat3x3d(matrix)) + let reference = simd_mul(scalar, simd_double3x3(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: multiply(Vec3d, Mat3x3d) + // covers: multiply(Mat3x3d, Vec3d) + @Test func multiplyVectorMatrix3x3d() { + let vectors = DifferentialInputs.vectors(dimension: 3, seed: 0x1104) + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1004).enumerated() { + let vector = SIMD3(vectors[index % vectors.count]) + let matrix = values + let rowTimesMatrix = FirebladeMath.multiply(vector, Mat3x3d(matrix)) + let matrixTimesColumn = FirebladeMath.multiply(Mat3x3d(matrix), vector) + #expect(differentialMatch(rowTimesMatrix.lanes, simd_mul(vector, simd_double3x3(elements: matrix)).lanes, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(matrixTimesColumn.lanes, simd_mul(simd_double3x3(elements: matrix), vector).lanes, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: transpose(Mat3x3d) + @Test func transpose3x3d() { + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1004).enumerated() { + let matrix = values + let portable = FirebladeMath.transpose(Mat3x3d(matrix)) + let reference = simd_transpose(simd_double3x3(elements: matrix)) + #expect(portable.elements == reference.elements, "input \(index)") + } + } + + // covers: determinant(Mat3x3d) + @Test func determinant3x3d() { + let singular: [[Double]] = [[Double](repeating: 0, count: 9), [1, 2, 3, 2, 4, 6, 0, 1, 5], [1, 2, 3, 4, 5, 6, 7, 8, 9]] + for (index, values) in (DifferentialInputs.matrices(dimension: 3, seed: 0x1004) + singular).enumerated() { + let matrix = values + let scale = Swift.max(1, values.map { Swift.abs($0) }.max() ?? 1) + let portable = FirebladeMath.determinant(Mat3x3d(matrix)) + let reference = simd_determinant(simd_double3x3(elements: matrix)) + #expect(differentialMatch(Double(portable), Double(reference), tolerance: Tolerance.double, scale: scale * scale * scale), "input \(index)") + } + } + + // covers: inverse(Mat3x3d) + @Test func inverse3x3d() { + for (index, values) in DifferentialInputs.matrices(dimension: 3, seed: 0x1004).enumerated() { + let matrix = values + let scale = Swift.max(1, values.map { Swift.abs($0) }.max() ?? 1) + // Skip near-singular inputs: both sides would only amplify rounding noise. + guard Swift.abs(Double(simd_determinant(simd_double3x3(elements: matrix)))) > 1e-2 * scale * scale * scale else { continue } + let portable = FirebladeMath.inverse(Mat3x3d(matrix)) + let reference = simd_inverse(simd_double3x3(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.doubleLoose), "input \(index)") + } + } + + // MARK: - Mat4x4f + + // covers: multiply(Mat4x4f, Mat4x4f) + @Test func multiplyMatrixMatrix4x4f() { + let inputs = DifferentialInputs.matrices(dimension: 4, seed: 0x1005) + for (index, values) in inputs.enumerated() { + let other = inputs[(index + 1) % inputs.count].floats + let lhs = values.floats + let portable = FirebladeMath.multiply(Mat4x4f(lhs), Mat4x4f(other)) + let reference = simd_mul(simd_float4x4(elements: lhs), simd_float4x4(elements: other)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: multiply(Float, Mat4x4f) + @Test func multiplyScalarMatrix4x4f() { + let scalars = DifferentialInputs.scalars(seed: 0x1205) + for (index, values) in DifferentialInputs.matrices(dimension: 4, seed: 0x1005).enumerated() { + let scalar = Float(scalars[index % scalars.count]) + let matrix = values.floats + let portable = FirebladeMath.multiply(scalar, Mat4x4f(matrix)) + let reference = simd_mul(scalar, simd_float4x4(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: multiply(Vec4f, Mat4x4f) + // covers: multiply(Mat4x4f, Vec4f) + @Test func multiplyVectorMatrix4x4f() { + let vectors = DifferentialInputs.vectors(dimension: 4, seed: 0x1105) + for (index, values) in DifferentialInputs.matrices(dimension: 4, seed: 0x1005).enumerated() { + let vector = SIMD4(vectors[index % vectors.count].floats) + let matrix = values.floats + let rowTimesMatrix = FirebladeMath.multiply(vector, Mat4x4f(matrix)) + let matrixTimesColumn = FirebladeMath.multiply(Mat4x4f(matrix), vector) + #expect(differentialMatch(rowTimesMatrix.lanes, simd_mul(vector, simd_float4x4(elements: matrix)).lanes, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(matrixTimesColumn.lanes, simd_mul(simd_float4x4(elements: matrix), vector).lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: transpose(Mat4x4f) + @Test func transpose4x4f() { + for (index, values) in DifferentialInputs.matrices(dimension: 4, seed: 0x1005).enumerated() { + let matrix = values.floats + let portable = FirebladeMath.transpose(Mat4x4f(matrix)) + let reference = simd_transpose(simd_float4x4(elements: matrix)) + #expect(portable.elements == reference.elements, "input \(index)") + } + } + + // MARK: - Mat4x4d + + // covers: multiply(Mat4x4d, Mat4x4d) + @Test func multiplyMatrixMatrix4x4d() { + let inputs = DifferentialInputs.matrices(dimension: 4, seed: 0x1006) + for (index, values) in inputs.enumerated() { + let other = inputs[(index + 1) % inputs.count] + let lhs = values + let portable = FirebladeMath.multiply(Mat4x4d(lhs), Mat4x4d(other)) + let reference = simd_mul(simd_double4x4(elements: lhs), simd_double4x4(elements: other)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: multiply(Double, Mat4x4d) + @Test func multiplyScalarMatrix4x4d() { + let scalars = DifferentialInputs.scalars(seed: 0x1206) + for (index, values) in DifferentialInputs.matrices(dimension: 4, seed: 0x1006).enumerated() { + let scalar = Double(scalars[index % scalars.count]) + let matrix = values + let portable = FirebladeMath.multiply(scalar, Mat4x4d(matrix)) + let reference = simd_mul(scalar, simd_double4x4(elements: matrix)) + #expect(differentialMatch(portable.elements, reference.elements, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: multiply(Vec4d, Mat4x4d) + // covers: multiply(Mat4x4d, Vec4d) + @Test func multiplyVectorMatrix4x4d() { + let vectors = DifferentialInputs.vectors(dimension: 4, seed: 0x1106) + for (index, values) in DifferentialInputs.matrices(dimension: 4, seed: 0x1006).enumerated() { + let vector = SIMD4(vectors[index % vectors.count]) + let matrix = values + let rowTimesMatrix = FirebladeMath.multiply(vector, Mat4x4d(matrix)) + let matrixTimesColumn = FirebladeMath.multiply(Mat4x4d(matrix), vector) + #expect(differentialMatch(rowTimesMatrix.lanes, simd_mul(vector, simd_double4x4(elements: matrix)).lanes, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(matrixTimesColumn.lanes, simd_mul(simd_double4x4(elements: matrix), vector).lanes, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: transpose(Mat4x4d) + @Test func transpose4x4d() { + for (index, values) in DifferentialInputs.matrices(dimension: 4, seed: 0x1006).enumerated() { + let matrix = values + let portable = FirebladeMath.transpose(Mat4x4d(matrix)) + let reference = simd_transpose(simd_double4x4(elements: matrix)) + #expect(portable.elements == reference.elements, "input \(index)") + } + } +} +#endif diff --git a/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDQuaternionTests.swift b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDQuaternionTests.swift new file mode 100644 index 0000000..ef5ac94 --- /dev/null +++ b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDQuaternionTests.swift @@ -0,0 +1,135 @@ +// +// PortableVsSIMDQuaternionTests.swift +// FirebladeMathTests +// + +// Differential tests for quaternion arithmetic and queries against Apple's simd. +// See PortableVsSIMDSupport.swift. + +#if canImport(simd) && !FRB_MATH_USE_SIMD +import FirebladeMath +import simd +import Testing + +struct PortableVsSIMDQuaternionTests { + /// Random quaternions, the basis, and the zero quaternion as an edge case. + static let quaternions = DifferentialInputs.quaternions(seed: 0x2001) + [[0, 0, 0, 0]] + + + // MARK: - Quat4f + + // covers: multiply(Float, Quat4f) + // covers: multiply(Quat4f, Float) + // covers: add(Quat4f, Quat4f) + // covers: subtract(Quat4f, Quat4f) + @Test func arithmeticF() { + let inputs = Self.quaternions + let scalars = DifferentialInputs.scalars(seed: 0x2002) + for (index, values) in inputs.enumerated() { + let lhs = SIMD4(values.floats), rhs = SIMD4(inputs[(index + 1) % inputs.count].floats) + let scalar = Float(scalars[index % scalars.count]) + let left = simd_quatf(vector: lhs), right = simd_quatf(vector: rhs) + #expect(differentialMatch(FirebladeMath.multiply(scalar, Quat4f(lhs)).elements, simd_mul(scalar, left).vector.lanes, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.multiply(Quat4f(lhs), scalar).elements, simd_mul(left, scalar).vector.lanes, tolerance: Tolerance.float), "input \(index)") + #expect(FirebladeMath.add(Quat4f(lhs), Quat4f(rhs)).elements == simd_add(left, right).vector.lanes, "input \(index)") + #expect(FirebladeMath.subtract(Quat4f(lhs), Quat4f(rhs)).elements == simd_sub(left, right).vector.lanes, "input \(index)") + } + } + + // covers: conjugate(Quat4f) + // covers: inverse(Quat4f) + // covers: length(Quat4f) + // covers: normalize(Quat4f) + @Test func queriesF() { + for (index, values) in Self.quaternions.enumerated() { + let vector = SIMD4(values.floats) + let reference = simd_quatf(vector: vector) + #expect(FirebladeMath.conjugate(Quat4f(vector)).elements == simd_conjugate(reference).vector.lanes, "input \(index)") + #expect(differentialMatch(FirebladeMath.inverse(Quat4f(vector)).elements, simd_inverse(reference).vector.lanes, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(Quat4f(vector))), Double(simd_length(reference)), tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(Quat4f(vector)).elements, simd_normalize(reference).vector.lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: angle(Quat4f) + // covers: axis(Quat4f) + @Test func angleAndAxisF() { + for (index, values) in Self.quaternions.enumerated() { + let vector = SIMD4(values.floats) + let reference = simd_quatf(vector: vector) + #expect(differentialMatch(Double(FirebladeMath.angle(Quat4f(vector))), Double(simd_angle(reference)), tolerance: Tolerance.floatLoose), "input \(index)") + #expect(differentialMatch(FirebladeMath.axis(Quat4f(vector)).lanes, simd_axis(reference).lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: act(Quat4f, Vec3f) + @Test func actF() { + let vectors = DifferentialInputs.vectors(dimension: 3, seed: 0x2003) + for (index, values) in Self.quaternions.enumerated() { + let quaternion = SIMD4(values.floats) + let vector = Vec3f(vectors[index % vectors.count].floats) + let portable = FirebladeMath.act(Quat4f(quaternion), vector) + let reference = simd_act(simd_normalize(simd_quatf(vector: quaternion)), vector) + #expect(differentialMatch(portable.lanes, reference.lanes, tolerance: Tolerance.floatLoose), "input \(index)") + } + } + + // MARK: - Quat4d + + // covers: multiply(Double, Quat4d) + // covers: multiply(Quat4d, Double) + // covers: add(Quat4d, Quat4d) + // covers: subtract(Quat4d, Quat4d) + @Test func arithmeticD() { + let inputs = Self.quaternions + let scalars = DifferentialInputs.scalars(seed: 0x2002) + for (index, values) in inputs.enumerated() { + let lhs = SIMD4(values), rhs = SIMD4(inputs[(index + 1) % inputs.count]) + let scalar = Double(scalars[index % scalars.count]) + let left = simd_quatd(vector: lhs), right = simd_quatd(vector: rhs) + #expect(differentialMatch(FirebladeMath.multiply(scalar, Quat4d(lhs)).elements, simd_mul(scalar, left).vector.lanes, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.multiply(Quat4d(lhs), scalar).elements, simd_mul(left, scalar).vector.lanes, tolerance: Tolerance.double), "input \(index)") + #expect(FirebladeMath.add(Quat4d(lhs), Quat4d(rhs)).elements == simd_add(left, right).vector.lanes, "input \(index)") + #expect(FirebladeMath.subtract(Quat4d(lhs), Quat4d(rhs)).elements == simd_sub(left, right).vector.lanes, "input \(index)") + } + } + + // covers: conjugate(Quat4d) + // covers: inverse(Quat4d) + // covers: length(Quat4d) + // covers: normalize(Quat4d) + @Test func queriesD() { + for (index, values) in Self.quaternions.enumerated() { + let vector = SIMD4(values) + let reference = simd_quatd(vector: vector) + #expect(FirebladeMath.conjugate(Quat4d(vector)).elements == simd_conjugate(reference).vector.lanes, "input \(index)") + #expect(differentialMatch(FirebladeMath.inverse(Quat4d(vector)).elements, simd_inverse(reference).vector.lanes, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(Quat4d(vector))), Double(simd_length(reference)), tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(Quat4d(vector)).elements, simd_normalize(reference).vector.lanes, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: angle(Quat4d) + // covers: axis(Quat4d) + @Test func angleAndAxisD() { + for (index, values) in Self.quaternions.enumerated() { + let vector = SIMD4(values) + let reference = simd_quatd(vector: vector) + #expect(differentialMatch(Double(FirebladeMath.angle(Quat4d(vector))), Double(simd_angle(reference)), tolerance: Tolerance.doubleLoose), "input \(index)") + #expect(differentialMatch(FirebladeMath.axis(Quat4d(vector)).lanes, simd_axis(reference).lanes, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: act(Quat4d, Vec3d) + @Test func actD() { + let vectors = DifferentialInputs.vectors(dimension: 3, seed: 0x2003) + for (index, values) in Self.quaternions.enumerated() { + let quaternion = SIMD4(values) + let vector = Vec3d(vectors[index % vectors.count]) + let portable = FirebladeMath.act(Quat4d(quaternion), vector) + let reference = simd_act(simd_normalize(simd_quatd(vector: quaternion)), vector) + #expect(differentialMatch(portable.lanes, reference.lanes, tolerance: Tolerance.doubleLoose), "input \(index)") + } + } +} +#endif diff --git a/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDScalarTests.swift b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDScalarTests.swift new file mode 100644 index 0000000..0426e17 --- /dev/null +++ b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDScalarTests.swift @@ -0,0 +1,148 @@ +// +// PortableVsSIMDScalarTests.swift +// FirebladeMathTests +// + +// Differential tests for scalar functions against Apple's simd, including NaN, infinities +// and signed zeros. See PortableVsSIMDSupport.swift. + +#if canImport(simd) && !FRB_MATH_USE_SIMD +import FirebladeMath +import simd +import Testing + +struct PortableVsSIMDScalarTests { + static let values = DifferentialInputs.scalars(seed: 0x4001) + DifferentialInputs.specialScalars + + + // MARK: - Float + + // covers: clamp(Float, Float, Float) + @Test func clampF() { + let values = Self.values + for (index, value) in values.enumerated() { + let a = values[(index + 7) % values.count], b = values[(index + 13) % values.count] + guard a.isFinite, b.isFinite else { continue } + let x = Float(value), low = Float(Swift.min(a, b)), high = Float(Swift.max(a, b)) + #expect(differentialMatch(Double(FirebladeMath.clamp(x, low, high)), Double(simd_clamp(x, low, high)), tolerance: 0, signedZero: true), "input \(index): clamp(\(x), \(low), \(high))") + } + } + + // covers: fract(Float) + // covers: sign(Float) + @Test func fractAndSignF() { + for (index, value) in Self.values.enumerated() { + let x = Float(value) + #expect(differentialMatch(Double(FirebladeMath.fract(x)), Double(simd_fract(x)), tolerance: Tolerance.float, signedZero: true), "input \(index): fract(\(x))") + #expect(differentialMatch(Double(FirebladeMath.sign(x)), Double(simd_sign(x)), tolerance: 0, signedZero: true), "input \(index): sign(\(x))") + } + } + + // covers: mix(Float, Float, Float) + @Test func mixF() { + let values = Self.values + var generator = SplitMix(seed: 0x4002) + for (index, value) in values.enumerated() { + let a = Float(value), b = Float(values[(index + 5) % values.count]) + let t = Float(generator.next(in: -0.5...1.5)) + let scale = Double(Swift.max(Swift.abs(a), Swift.abs(b))) + #expect(differentialMatch(Double(FirebladeMath.mix(a, b, t)), Double(simd_mix(a, b, t)), tolerance: Tolerance.float, scale: scale), "input \(index): mix(\(a), \(b), \(t))") + } + } + + // covers: rsqrt(Float) + @Test func rsqrtF() { + for (index, value) in Self.values.enumerated() where value >= 0 || value.isNaN { + let x = Float(value) + #expect(differentialMatch(Double(FirebladeMath.rsqrt(x)), Double(simd_rsqrt(x)), tolerance: Tolerance.float), "input \(index): rsqrt(\(x))") + } + } + + // covers: smoothstep(Float, Float, Float) + @Test func smoothstepF() { + let values = Self.values + var generator = SplitMix(seed: 0x4003) + for (index, value) in values.enumerated() where value.isFinite { + let edge0 = Float(generator.next(in: -10...0)), edge1 = Float(generator.next(in: 0.5...10)) + let x = Float(value) + #expect(differentialMatch(Double(FirebladeMath.smoothstep(edge0, edge1, x)), Double(simd_smoothstep(edge0, edge1, x)), tolerance: Tolerance.float), "input \(index): smoothstep(\(edge0), \(edge1), \(x))") + } + } + + // covers: step(Float, Float) + @Test func stepF() { + let values = Self.values + for (index, value) in values.enumerated() { + let x = Float(value) + for edge in [Float(values[(index + 3) % values.count]), x] { + #expect(differentialMatch(Double(FirebladeMath.step(x, edge)), Double(simd.step(x, edge: edge)), tolerance: 0, signedZero: true), "input \(index): step(\(x), \(edge))") + } + } + } + + // MARK: - Double + + // covers: clamp(Double, Double, Double) + @Test func clampD() { + let values = Self.values + for (index, value) in values.enumerated() { + let a = values[(index + 7) % values.count], b = values[(index + 13) % values.count] + guard a.isFinite, b.isFinite else { continue } + let x = Double(value), low = Double(Swift.min(a, b)), high = Double(Swift.max(a, b)) + #expect(differentialMatch(Double(FirebladeMath.clamp(x, low, high)), Double(simd_clamp(x, low, high)), tolerance: 0, signedZero: true), "input \(index): clamp(\(x), \(low), \(high))") + } + } + + // covers: fract(Double) + // covers: sign(Double) + @Test func fractAndSignD() { + for (index, value) in Self.values.enumerated() { + let x = Double(value) + #expect(differentialMatch(Double(FirebladeMath.fract(x)), Double(simd_fract(x)), tolerance: Tolerance.double, signedZero: true), "input \(index): fract(\(x))") + #expect(differentialMatch(Double(FirebladeMath.sign(x)), Double(simd_sign(x)), tolerance: 0, signedZero: true), "input \(index): sign(\(x))") + } + } + + // covers: mix(Double, Double, Double) + @Test func mixD() { + let values = Self.values + var generator = SplitMix(seed: 0x4002) + for (index, value) in values.enumerated() { + let a = Double(value), b = Double(values[(index + 5) % values.count]) + let t = Double(generator.next(in: -0.5...1.5)) + let scale = Double(Swift.max(Swift.abs(a), Swift.abs(b))) + #expect(differentialMatch(Double(FirebladeMath.mix(a, b, t)), Double(simd_mix(a, b, t)), tolerance: Tolerance.double, scale: scale), "input \(index): mix(\(a), \(b), \(t))") + } + } + + // covers: rsqrt(Double) + @Test func rsqrtD() { + for (index, value) in Self.values.enumerated() where value >= 0 || value.isNaN { + let x = Double(value) + #expect(differentialMatch(Double(FirebladeMath.rsqrt(x)), Double(simd_rsqrt(x)), tolerance: Tolerance.double), "input \(index): rsqrt(\(x))") + } + } + + // covers: smoothstep(Double, Double, Double) + @Test func smoothstepD() { + let values = Self.values + var generator = SplitMix(seed: 0x4003) + for (index, value) in values.enumerated() where value.isFinite { + let edge0 = Double(generator.next(in: -10...0)), edge1 = Double(generator.next(in: 0.5...10)) + let x = Double(value) + #expect(differentialMatch(Double(FirebladeMath.smoothstep(edge0, edge1, x)), Double(simd_smoothstep(edge0, edge1, x)), tolerance: Tolerance.double), "input \(index): smoothstep(\(edge0), \(edge1), \(x))") + } + } + + // covers: step(Double, Double) + @Test func stepD() { + let values = Self.values + for (index, value) in values.enumerated() { + let x = Double(value) + for edge in [Double(values[(index + 3) % values.count]), x] { + #expect(differentialMatch(Double(FirebladeMath.step(x, edge)), Double(simd.step(x, edge: edge)), tolerance: 0, signedZero: true), "input \(index): step(\(x), \(edge))") + } + } + } +} +#endif diff --git a/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDSupport.swift b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDSupport.swift new file mode 100644 index 0000000..9e0ce0e --- /dev/null +++ b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDSupport.swift @@ -0,0 +1,196 @@ +// +// PortableVsSIMDSupport.swift +// FirebladeMathTests +// + +// Shared inputs and helpers for the PortableVsSIMD differential tests, which compare +// FirebladeMath's portable implementations against Apple's simd library. +// +// Only compiled when Apple's simd is available but the `simd` trait is disabled +// (`swift test --disable-default-traits` on Apple platforms). +// +// Every test is annotated with `// covers: name(ParameterTypes)` lines naming the +// FirebladeMath functions it compares. Scripts/check-simd-differential-coverage.py +// fails CI if a function with a simd branch has no such annotation. + +#if canImport(simd) && !FRB_MATH_USE_SIMD +import FirebladeMath +import simd + +// MARK: - Deterministic inputs + +/// SplitMix64: a small deterministic generator, so failures are reproducible. +struct SplitMix { + var state: UInt64 + + init(seed: UInt64) { state = seed } + + mutating func next(in range: ClosedRange) -> Double { + state &+= 0x9E37_79B9_7F4A_7C15 + var z = state + z = (z ^ (z >> 30)) &* 0xBF58_476D_1CE4_E5B9 + z = (z ^ (z >> 27)) &* 0x94D0_49BB_1331_11EB + z ^= z >> 31 + let unit = Double(z >> 11) / Double(1 << 53) + return range.lowerBound + unit * (range.upperBound - range.lowerBound) + } + + mutating func values(_ count: Int, in range: ClosedRange = -2...2) -> [Double] { + (0.. [[Double]] { + var generator = SplitMix(seed: seed) + var result = (0.. Double in // cyclic permutation + let column = index / n, row = index % n + return row == (column + 1) % n ? 1 : 0 + }) + result.append(identity.map { $0 * 1e-3 }) + result.append(result[0].map { $0 * 1e3 }) + return result + } + + /// Vectors as element lists: pseudo-random plus zero, unit axes, negative, tiny and large. + static func vectors(dimension n: Int, seed: UInt64) -> [[Double]] { + var generator = SplitMix(seed: seed) + var result = (0.. [[Double]] { + var generator = SplitMix(seed: seed) + var result = (0.. [[Double]] { + var generator = SplitMix(seed: seed) + var result = (0.. [Double] in + let q = simd_normalize(simd_quatd(vector: SIMD4(generator.values(4)))) + return Array(q.vector) + } + result += [[0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]] + let half = 0.5.squareRoot() + result += [[half, half, 0, 0], [0, half, half, 0], [half, 0, half, 0]] // 180 degrees about diagonals + return result + } + + /// Scalars: pseudo-random over a wide range plus special values. + static func scalars(seed: UInt64) -> [Double] { + var generator = SplitMix(seed: seed) + return generator.values(randomCount, in: -100...100) + [ + 0, -0.0, 1, -1, 0.5, -0.5, 2.5, -2.5, 1e-30, -1e-30, 1e30, -1e30 + ] + } + + /// Special floating-point values, for functions where both implementations must agree on them. + static let specialScalars: [Double] = [.nan, .infinity, -.infinity, 0, -0.0] +} + +// MARK: - simd conversions + +extension simd_float2x2 { + init(elements values: [Float]) { self.init(SIMD2(values[0..<2]), SIMD2(values[2..<4])) } + var elements: [Float] { [columns.0, columns.1].flatMap { [$0.x, $0.y] } } +} + +extension simd_double2x2 { + init(elements values: [Double]) { self.init(SIMD2(values[0..<2]), SIMD2(values[2..<4])) } + var elements: [Double] { [columns.0, columns.1].flatMap { [$0.x, $0.y] } } +} + +extension simd_float3x3 { + init(elements values: [Float]) { self.init(SIMD3(values[0..<3]), SIMD3(values[3..<6]), SIMD3(values[6..<9])) } + var elements: [Float] { [columns.0, columns.1, columns.2].flatMap { [$0.x, $0.y, $0.z] } } +} + +extension simd_double3x3 { + init(elements values: [Double]) { self.init(SIMD3(values[0..<3]), SIMD3(values[3..<6]), SIMD3(values[6..<9])) } + var elements: [Double] { [columns.0, columns.1, columns.2].flatMap { [$0.x, $0.y, $0.z] } } +} + +extension simd_float4x4 { + init(elements values: [Float]) { + self.init(SIMD4(values[0..<4]), SIMD4(values[4..<8]), SIMD4(values[8..<12]), SIMD4(values[12..<16])) + } + + var elements: [Float] { [columns.0, columns.1, columns.2, columns.3].flatMap { [$0.x, $0.y, $0.z, $0.w] } } +} + +extension simd_double4x4 { + init(elements values: [Double]) { + self.init(SIMD4(values[0..<4]), SIMD4(values[4..<8]), SIMD4(values[8..<12]), SIMD4(values[12..<16])) + } + + var elements: [Double] { [columns.0, columns.1, columns.2, columns.3].flatMap { [$0.x, $0.y, $0.z, $0.w] } } +} + +extension Array where Element == Double { + var floats: [Float] { map(Float.init) } +} + +extension SIMD { + var lanes: [Scalar] { indices.map { self[$0] } } +} + +// MARK: - Comparison + +/// Floating-point comparison for differential tests. +/// +/// Values match if they are bit-for-bit equal in kind (both NaN, same infinity, same signed +/// zero when `signedZero` is set) or differ by at most `tolerance` relative to `scale`. +func differentialMatch( + _ portable: Double, + _ reference: Double, + tolerance: Double, + scale: Double = 1, + signedZero: Bool = false +) -> Bool { + if portable.isNaN || reference.isNaN { return portable.isNaN && reference.isNaN } + if portable.isInfinite || reference.isInfinite { return portable == reference } + if signedZero, portable == 0, reference == 0 { return portable.sign == reference.sign } + return Swift.abs(portable - reference) <= tolerance * Swift.max(1, scale, Swift.abs(reference)) +} + +/// Element-wise `differentialMatch`, relative to the largest finite reference element. +func differentialMatch(_ portable: [Double], _ reference: [Double], tolerance: Double) -> Bool { + let scale = reference.filter(\.isFinite).map { Swift.abs($0) }.max() ?? 1 + return portable.count == reference.count + && zip(portable, reference).allSatisfy { differentialMatch($0, $1, tolerance: tolerance, scale: scale) } +} + +func differentialMatch(_ portable: [Float], _ reference: [Float], tolerance: Double) -> Bool { + differentialMatch(portable.map(Double.init), reference.map(Double.init), tolerance: tolerance) +} + +/// Tolerances relative to the result magnitude; both sides differ only by rounding. +enum Tolerance { + /// Results computed by the same few operations, possibly in a different order or with FMA. + static let float: Double = 2e-6 + static let double: Double = 1e-14 + /// Results with long operation chains (inverses, transcendental functions). + static let floatLoose: Double = 2e-4 + static let doubleLoose: Double = 5e-12 +} +#endif diff --git a/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDVectorTests.swift b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDVectorTests.swift new file mode 100644 index 0000000..0833d47 --- /dev/null +++ b/Tests/FirebladeMathTests/PortableVsSIMD/PortableVsSIMDVectorTests.swift @@ -0,0 +1,101 @@ +// +// PortableVsSIMDVectorTests.swift +// FirebladeMathTests +// + +// Differential tests for vector functions against Apple's simd. See PortableVsSIMDSupport.swift. + +#if canImport(simd) && !FRB_MATH_USE_SIMD +import FirebladeMath +import simd +import Testing + +struct PortableVsSIMDVectorTests { + + // covers: dot(SIMD2, SIMD2) + // covers: length(SIMD2) + // covers: normalize(SIMD2) + // covers: cross(SIMD2, SIMD2) + @Test func vector2f() { + let inputs = DifferentialInputs.vectors(dimension: 2, seed: 0x3021) + for (index, values) in inputs.enumerated() { + let lhs = SIMD2(values.floats), rhs = SIMD2(inputs[(index + 1) % inputs.count].floats) + #expect(differentialMatch(Double(FirebladeMath.dot(lhs, rhs)), Double(simd_dot(lhs, rhs)), tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(lhs)), Double(simd_length(lhs)), tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(lhs).lanes, simd_normalize(lhs).lanes, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.cross(lhs, rhs).lanes, simd_cross(lhs, rhs).lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: dot(SIMD3, SIMD3) + // covers: length(SIMD3) + // covers: normalize(SIMD3) + // covers: cross(SIMD3, SIMD3) + @Test func vector3f() { + let inputs = DifferentialInputs.vectors(dimension: 3, seed: 0x3031) + for (index, values) in inputs.enumerated() { + let lhs = SIMD3(values.floats), rhs = SIMD3(inputs[(index + 1) % inputs.count].floats) + #expect(differentialMatch(Double(FirebladeMath.dot(lhs, rhs)), Double(simd_dot(lhs, rhs)), tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(lhs)), Double(simd_length(lhs)), tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(lhs).lanes, simd_normalize(lhs).lanes, tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.cross(lhs, rhs).lanes, simd_cross(lhs, rhs).lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: dot(SIMD4, SIMD4) + // covers: length(SIMD4) + // covers: normalize(SIMD4) + @Test func vector4f() { + let inputs = DifferentialInputs.vectors(dimension: 4, seed: 0x3041) + for (index, values) in inputs.enumerated() { + let lhs = SIMD4(values.floats), rhs = SIMD4(inputs[(index + 1) % inputs.count].floats) + #expect(differentialMatch(Double(FirebladeMath.dot(lhs, rhs)), Double(simd_dot(lhs, rhs)), tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(lhs)), Double(simd_length(lhs)), tolerance: Tolerance.float), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(lhs).lanes, simd_normalize(lhs).lanes, tolerance: Tolerance.float), "input \(index)") + } + } + + // covers: dot(SIMD2, SIMD2) + // covers: length(SIMD2) + // covers: normalize(SIMD2) + // covers: cross(SIMD2, SIMD2) + @Test func vector2d() { + let inputs = DifferentialInputs.vectors(dimension: 2, seed: 0x3022) + for (index, values) in inputs.enumerated() { + let lhs = SIMD2(values), rhs = SIMD2(inputs[(index + 1) % inputs.count]) + #expect(differentialMatch(Double(FirebladeMath.dot(lhs, rhs)), Double(simd_dot(lhs, rhs)), tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(lhs)), Double(simd_length(lhs)), tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(lhs).lanes, simd_normalize(lhs).lanes, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.cross(lhs, rhs).lanes, simd_cross(lhs, rhs).lanes, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: dot(SIMD3, SIMD3) + // covers: length(SIMD3) + // covers: normalize(SIMD3) + // covers: cross(SIMD3, SIMD3) + @Test func vector3d() { + let inputs = DifferentialInputs.vectors(dimension: 3, seed: 0x3032) + for (index, values) in inputs.enumerated() { + let lhs = SIMD3(values), rhs = SIMD3(inputs[(index + 1) % inputs.count]) + #expect(differentialMatch(Double(FirebladeMath.dot(lhs, rhs)), Double(simd_dot(lhs, rhs)), tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(lhs)), Double(simd_length(lhs)), tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(lhs).lanes, simd_normalize(lhs).lanes, tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.cross(lhs, rhs).lanes, simd_cross(lhs, rhs).lanes, tolerance: Tolerance.double), "input \(index)") + } + } + + // covers: dot(SIMD4, SIMD4) + // covers: length(SIMD4) + // covers: normalize(SIMD4) + @Test func vector4d() { + let inputs = DifferentialInputs.vectors(dimension: 4, seed: 0x3042) + for (index, values) in inputs.enumerated() { + let lhs = SIMD4(values), rhs = SIMD4(inputs[(index + 1) % inputs.count]) + #expect(differentialMatch(Double(FirebladeMath.dot(lhs, rhs)), Double(simd_dot(lhs, rhs)), tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(Double(FirebladeMath.length(lhs)), Double(simd_length(lhs)), tolerance: Tolerance.double), "input \(index)") + #expect(differentialMatch(FirebladeMath.normalize(lhs).lanes, simd_normalize(lhs).lanes, tolerance: Tolerance.double), "input \(index)") + } + } +} +#endif