Add macro-controlled SIMD support (SSE/NEON) to math library

.gitignore

@@ -15,3 +15,4 @@
 /active_project
 /test/assets/assets.db
 /test/products/
+_codeql_detected_source_root

cmake/configuration.cmake

@@ -27,3 +27,7 @@ if (SKY_EDITOR OR SKY_BUILD_TOOL)
 endif ()
 
 add_compile_definitions("$<$<CONFIG:Debug>:_DEBUG;DEBUG>")
+
+if (SKY_MATH_SIMD)
+    add_compile_definitions(SKY_MATH_SIMD=1)
+endif ()

cmake/options.cmake

@@ -13,4 +13,5 @@ option(SKY_USE_TRACY "use tracy profiler" OFF)
 
 # todo: controlled by project config json.
 option(SKY_BUILD_XR "xr plugin" OFF)
-option(SKY_BUILD_PYTHON "python plugin" OFF)
+option(SKY_BUILD_PYTHON "python plugin" OFF)
+option(SKY_MATH_SIMD "enable simd math" ON)

runtime/core/include/core/math/MathSimd.h

@@ -0,0 +1,31 @@
+//
+// Created for SIMD math support.
+//
+
+#pragma once
+
+// SIMD platform detection
+#if defined(SKY_MATH_SIMD) && SKY_MATH_SIMD
+
+    #if defined(__SSE__) || defined(_M_X64) || (defined(_M_IX86_FP) && _M_IX86_FP >= 1)
+        #define SKY_SIMD_SSE 1
+        #include <xmmintrin.h> // SSE
+        #include <emmintrin.h> // SSE2
+    #endif
+
+    #if defined(__ARM_NEON) || defined(__ARM_NEON__)
+        #define SKY_SIMD_NEON 1
+        #include <arm_neon.h>
+    #endif
+
+#endif
+
+#ifndef SKY_SIMD_SSE
+    #define SKY_SIMD_SSE 0
+#endif
+
+#ifndef SKY_SIMD_NEON
+    #define SKY_SIMD_NEON 0
+#endif
+
+#define SKY_SIMD_ENABLED (SKY_SIMD_SSE || SKY_SIMD_NEON)

runtime/core/include/core/math/Matrix4.inl

@@ -32,13 +32,37 @@ namespace sky {
 
     inline Matrix4 Matrix4::operator*(const Matrix4& rhs) const
     {
+#if SKY_SIMD_SSE
+        Matrix4 ret;
+        for (int i = 0; i < 4; ++i) {
+            __m128 e0 = _mm_set1_ps(rhs.m[i].v[0]);
+            __m128 e1 = _mm_set1_ps(rhs.m[i].v[1]);
+            __m128 e2 = _mm_set1_ps(rhs.m[i].v[2]);
+            __m128 e3 = _mm_set1_ps(rhs.m[i].v[3]);
+            ret.m[i].simd = _mm_add_ps(
+                _mm_add_ps(_mm_mul_ps(m[0].simd, e0), _mm_mul_ps(m[1].simd, e1)),
+                _mm_add_ps(_mm_mul_ps(m[2].simd, e2), _mm_mul_ps(m[3].simd, e3)));
+        }
+        return ret;
+#elif SKY_SIMD_NEON
+        Matrix4 ret;
+        for (int i = 0; i < 4; ++i) {
+            float32x4_t r = vmulq_n_f32(m[0].simd, rhs.m[i].v[0]);
+            r = vmlaq_n_f32(r, m[1].simd, rhs.m[i].v[1]);
+            r = vmlaq_n_f32(r, m[2].simd, rhs.m[i].v[2]);
+            r = vmlaq_n_f32(r, m[3].simd, rhs.m[i].v[3]);
+            ret.m[i].simd = r;
+        }
+        return ret;
+#else
         const auto& rw = rhs.m;
         Matrix4 ret;
         ret[0] = m[0] * rw[0][0] + m[1] * rw[0][1] + m[2] * rw[0][2] + m[3] * rw[0][3];
         ret[1] = m[0] * rw[1][0] + m[1] * rw[1][1] + m[2] * rw[1][2] + m[3] * rw[1][3];
         ret[2] = m[0] * rw[2][0] + m[1] * rw[2][1] + m[2] * rw[2][2] + m[3] * rw[2][3];
         ret[3] = m[0] * rw[3][0] + m[1] * rw[3][1] + m[2] * rw[3][2] + m[3] * rw[3][3];
         return ret;
+#endif
     }
 
     inline Matrix4 Matrix4::operator*(float multiplier) const
@@ -95,11 +119,31 @@ namespace sky {
 
     inline Vector4 Matrix4::operator*(const Vector4& rhs) const
     {
+#if SKY_SIMD_SSE
+        __m128 e0 = _mm_set1_ps(rhs.v[0]);
+        __m128 e1 = _mm_set1_ps(rhs.v[1]);
+        __m128 e2 = _mm_set1_ps(rhs.v[2]);
+        __m128 e3 = _mm_set1_ps(rhs.v[3]);
+        Vector4 result;
+        result.simd = _mm_add_ps(
+            _mm_add_ps(_mm_mul_ps(m[0].simd, e0), _mm_mul_ps(m[1].simd, e1)),
+            _mm_add_ps(_mm_mul_ps(m[2].simd, e2), _mm_mul_ps(m[3].simd, e3)));
+        return result;
+#elif SKY_SIMD_NEON
+        float32x4_t r = vmulq_n_f32(m[0].simd, rhs.v[0]);
+        r = vmlaq_n_f32(r, m[1].simd, rhs.v[1]);
+        r = vmlaq_n_f32(r, m[2].simd, rhs.v[2]);
+        r = vmlaq_n_f32(r, m[3].simd, rhs.v[3]);
+        Vector4 result;
+        result.simd = r;
+        return result;
+#else
         Vector4 v0 = m[0] * rhs[0];
         Vector4 v1 = m[1] * rhs[1];
         Vector4 v2 = m[2] * rhs[2];
         Vector4 v3 = m[3] * rhs[3];
         return (v0 + v1) + (v2 + v3);
+#endif
     }
 
     inline Vector4 &Matrix4::operator[](uint32_t i)

runtime/core/include/core/math/Vector4.h

@@ -6,10 +6,11 @@
 #pragma once
 
 #include "core/math/Math.h"
+#include "core/math/MathSimd.h"
 
 namespace sky {
 
-    struct Vector4 {
+    struct alignas(16) Vector4 {
         union {
             float v[4];
             struct {
@@ -18,6 +19,11 @@ namespace sky {
                 float z;
                 float w;
             };
+#if SKY_SIMD_SSE
+            __m128 simd;
+#elif SKY_SIMD_NEON
+            float32x4_t simd;
+#endif
         };
 
         inline constexpr Vector4();

runtime/core/include/core/math/Vector4.inl

@@ -53,55 +53,91 @@ namespace sky {
 
     inline Vector4& Vector4::operator+=(const Vector4& rhs)
     {
+#if SKY_SIMD_SSE
+        simd = _mm_add_ps(simd, rhs.simd);
+#elif SKY_SIMD_NEON
+        simd = vaddq_f32(simd, rhs.simd);
+#else
         x += rhs.x;
         y += rhs.y;
         z += rhs.z;
         w += rhs.w;
+#endif
         return *this;
     }
 
     inline Vector4& Vector4::operator-=(const Vector4& rhs)
     {
+#if SKY_SIMD_SSE
+        simd = _mm_sub_ps(simd, rhs.simd);
+#elif SKY_SIMD_NEON
+        simd = vsubq_f32(simd, rhs.simd);
+#else
         x -= rhs.x;
         y -= rhs.y;
         z -= rhs.z;
         w -= rhs.w;
+#endif
         return *this;
     }
 
     inline Vector4& Vector4::operator*=(const Vector4& rhs)
     {
+#if SKY_SIMD_SSE
+        simd = _mm_mul_ps(simd, rhs.simd);
+#elif SKY_SIMD_NEON
+        simd = vmulq_f32(simd, rhs.simd);
+#else
         x *= rhs.x;
         y *= rhs.y;
         z *= rhs.z;
         w *= rhs.w;
+#endif
         return *this;
     }
 
     inline Vector4& Vector4::operator/=(const Vector4& rhs)
     {
+#if SKY_SIMD_SSE
+        simd = _mm_div_ps(simd, rhs.simd);
+#elif SKY_SIMD_NEON
+        simd = vdivq_f32(simd, rhs.simd);
+#else
         x /= rhs.x;
         y /= rhs.y;
         z /= rhs.z;
         w /= rhs.w;
+#endif
         return *this;
     }
 
     inline Vector4& Vector4::operator*=(float m)
     {
+#if SKY_SIMD_SSE
+        simd = _mm_mul_ps(simd, _mm_set1_ps(m));
+#elif SKY_SIMD_NEON
+        simd = vmulq_n_f32(simd, m);
+#else
         x *= m;
         y *= m;
         z *= m;
         w *= m;
+#endif
         return *this;
     }
 
     inline Vector4& Vector4::operator/=(float d)
     {
+#if SKY_SIMD_SSE
+        simd = _mm_div_ps(simd, _mm_set1_ps(d));
+#elif SKY_SIMD_NEON
+        simd = vdivq_f32(simd, vdupq_n_f32(d));
+#else
         x /= d;
         y /= d;
         z /= d;
         w /= d;
+#endif
         return *this;
     }
 
@@ -117,13 +153,51 @@ namespace sky {
 
     inline float Vector4::Dot(const Vector4 &rhs) const
     {
+#if SKY_SIMD_SSE
+        __m128 mul = _mm_mul_ps(simd, rhs.simd);
+        __m128 shuf1 = _mm_shuffle_ps(mul, mul, _MM_SHUFFLE(2, 3, 0, 1));
+        __m128 sums1 = _mm_add_ps(mul, shuf1);
+        __m128 shuf2 = _mm_shuffle_ps(sums1, sums1, _MM_SHUFFLE(0, 1, 2, 3));
+        __m128 sums2 = _mm_add_ps(sums1, shuf2);
+        return _mm_cvtss_f32(sums2);
+#elif SKY_SIMD_NEON
+        float32x4_t mul = vmulq_f32(simd, rhs.simd);
+        float32x2_t sum = vadd_f32(vget_low_f32(mul), vget_high_f32(mul));
+        sum = vpadd_f32(sum, sum);
+        return vget_lane_f32(sum, 0);
+#else
         Vector4 ret = (*this) * rhs;
         return (ret.x + ret.y) + (ret.z + ret.w);
+#endif
     }
 
     inline void Vector4::Normalize()
     {
+#if SKY_SIMD_SSE
+        __m128 dot = _mm_mul_ps(simd, simd);
+        __m128 shuf1 = _mm_shuffle_ps(dot, dot, _MM_SHUFFLE(2, 3, 0, 1));
+        __m128 sums1 = _mm_add_ps(dot, shuf1);
+        __m128 shuf2 = _mm_shuffle_ps(sums1, sums1, _MM_SHUFFLE(0, 1, 2, 3));
+        __m128 sums2 = _mm_add_ps(sums1, shuf2);
+        __m128 inv = _mm_rsqrt_ps(sums2);
+        // Newton-Raphson refinement: inv = inv * (1.5 - 0.5 * sums2 * inv * inv)
+        __m128 half = _mm_set1_ps(0.5f);
+        __m128 three_half = _mm_set1_ps(1.5f);
+        __m128 muls = _mm_mul_ps(_mm_mul_ps(half, sums2), _mm_mul_ps(inv, inv));
+        inv = _mm_mul_ps(inv, _mm_sub_ps(three_half, muls));
+        simd = _mm_mul_ps(simd, inv);
+#elif SKY_SIMD_NEON
+        float32x4_t mul = vmulq_f32(simd, simd);
+        float32x2_t sum = vadd_f32(vget_low_f32(mul), vget_high_f32(mul));
+        sum = vpadd_f32(sum, sum);
+        float32x4_t dotVal = vdupq_lane_f32(sum, 0);
+        float32x4_t inv = vrsqrteq_f32(dotVal);
+        // Newton-Raphson refinement
+        inv = vmulq_f32(inv, vrsqrtsq_f32(vmulq_f32(dotVal, inv), inv));
+        simd = vmulq_f32(simd, inv);
+#else
         float inverseSqrt = 1 / sqrt(Dot(*this));
         Vector4::operator*=(inverseSqrt);
+#endif
     }
 }