From 9725b60d46e8d8d0c1afbbbe2349ba10c34609e1 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Fri, 23 Apr 2021 22:12:57 +0300
Subject: [PATCH 01/16] rename glmm_shuff1x() to glmm_splat()

* mark glmm_shuff1x() as DEPRECATED
---
 include/cglm/simd/sse2/affine.h | 45 ++++++++++++++++-----------------
 include/cglm/simd/sse2/mat4.h   | 17 ++++++-------
 include/cglm/simd/sse2/quat.h   |  8 +++---
 include/cglm/simd/x86.h         |  8 ++++--
 include/cglm/vec4-ext.h         |  6 ++---
 5 files changed, 43 insertions(+), 41 deletions(-)

diff --git a/include/cglm/simd/sse2/affine.h b/include/cglm/simd/sse2/affine.h
index 236408c..b5d64f0 100644
--- a/include/cglm/simd/sse2/affine.h
+++ b/include/cglm/simd/sse2/affine.h
@@ -25,29 +25,28 @@ glm_mul_sse2(mat4 m1, mat4 m2, mat4 dest) {
   
   r = glmm_load(m2[0]);
   glmm_store(dest[0],
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,
-                                   _mm_mul_ps(glmm_shuff1x(r, 2), l2))));
+             glmm_fmadd(glmm_splat(r, 0), l0,
+                        glmm_fmadd(glmm_splat(r, 1), l1,
+                                   _mm_mul_ps(glmm_splat(r, 2), l2))));
   
   r = glmm_load(m2[1]);
   glmm_store(dest[1],
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,
-                                   _mm_mul_ps(glmm_shuff1x(r, 2), l2))));
+             glmm_fmadd(glmm_splat(r, 0), l0,
+                        glmm_fmadd(glmm_splat(r, 1), l1,
+                                   _mm_mul_ps(glmm_splat(r, 2), l2))));
 
   r = glmm_load(m2[2]);
   glmm_store(dest[2],
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,
-                                   _mm_mul_ps(glmm_shuff1x(r, 2), l2))));
+             glmm_fmadd(glmm_splat(r, 0), l0,
+                        glmm_fmadd(glmm_splat(r, 1), l1,
+                                   _mm_mul_ps(glmm_splat(r, 2), l2))));
 
   r = glmm_load(m2[3]);
   glmm_store(dest[3],
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,
-                                   glmm_fmadd(glmm_shuff1x(r, 2), l2,
-                                              _mm_mul_ps(glmm_shuff1x(r, 3),
-                                                         l3)))));
+       glmm_fmadd(glmm_splat(r, 0), l0,
+                  glmm_fmadd(glmm_splat(r, 1), l1,
+                             glmm_fmadd(glmm_splat(r, 2), l2,
+                                        _mm_mul_ps(glmm_splat(r, 3), l3)))));
 }
 
 CGLM_INLINE
@@ -63,22 +62,22 @@ glm_mul_rot_sse2(mat4 m1, mat4 m2, mat4 dest) {
 
   r = glmm_load(m2[0]);
   glmm_store(dest[0],
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,
-                                   _mm_mul_ps(glmm_shuff1x(r, 2), l2))));
+             glmm_fmadd(glmm_splat(r, 0), l0,
+                        glmm_fmadd(glmm_splat(r, 1), l1,
+                                   _mm_mul_ps(glmm_splat(r, 2), l2))));
   
   r = glmm_load(m2[1]);
   glmm_store(dest[1],
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,
-                                   _mm_mul_ps(glmm_shuff1x(r, 2), l2))));
+             glmm_fmadd(glmm_splat(r, 0), l0,
+                        glmm_fmadd(glmm_splat(r, 1), l1,
+                                   _mm_mul_ps(glmm_splat(r, 2), l2))));
   
   
   r = glmm_load(m2[2]);
   glmm_store(dest[2],
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,
-                                   _mm_mul_ps(glmm_shuff1x(r, 2), l2))));
+             glmm_fmadd(glmm_splat(r, 0), l0,
+                        glmm_fmadd(glmm_splat(r, 1), l1,
+                                   _mm_mul_ps(glmm_splat(r, 2), l2))));
 
   glmm_store(dest[3], l3);
 }
diff --git a/include/cglm/simd/sse2/mat4.h b/include/cglm/simd/sse2/mat4.h
index 78fac21..f5d7135 100644
--- a/include/cglm/simd/sse2/mat4.h
+++ b/include/cglm/simd/sse2/mat4.h
@@ -60,11 +60,10 @@ glm_mat4_mul_sse2(mat4 m1, mat4 m2, mat4 dest) {
                                                                               \
   r = glmm_load(m2[C]);                                                       \
   glmm_store(dest[C],                                                         \
-             glmm_fmadd(glmm_shuff1x(r, 0), l0,                               \
-                        glmm_fmadd(glmm_shuff1x(r, 1), l1,                    \
-                                   glmm_fmadd(glmm_shuff1x(r, 2), l2,         \
-                                              _mm_mul_ps(glmm_shuff1x(r, 3),  \
-                                                         l3)))));
+       glmm_fmadd(glmm_splat(r, 0), l0,                                       \
+                   glmm_fmadd(glmm_splat(r, 1), l1,                           \
+                             glmm_fmadd(glmm_splat(r, 2), l2,                 \
+                                         _mm_mul_ps(glmm_splat(r, 3), l3)))));
 
   XX(0);
   XX(1);
@@ -80,11 +79,11 @@ glm_mat4_mulv_sse2(mat4 m, vec4 v, vec4 dest) {
   __m128 x0, x1;
   
   x0 = glmm_load(v);
-  x1 = glmm_fmadd(glmm_load(m[0]), glmm_shuff1x(x0, 0),
-                  glmm_fmadd(glmm_load(m[1]), glmm_shuff1x(x0, 1),
-                             glmm_fmadd(glmm_load(m[2]), glmm_shuff1x(x0, 2),
+  x1 = glmm_fmadd(glmm_load(m[0]), glmm_splat(x0, 0),
+                  glmm_fmadd(glmm_load(m[1]), glmm_splat(x0, 1),
+                             glmm_fmadd(glmm_load(m[2]), glmm_splat(x0, 2),
                                         _mm_mul_ps(glmm_load(m[3]),
-                                                   glmm_shuff1x(x0, 3)))));
+                                                   glmm_splat(x0, 3)))));
   
   glmm_store(dest, x1);
 }
diff --git a/include/cglm/simd/sse2/quat.h b/include/cglm/simd/sse2/quat.h
index 0173f94..894d492 100644
--- a/include/cglm/simd/sse2/quat.h
+++ b/include/cglm/simd/sse2/quat.h
@@ -27,15 +27,15 @@ glm_quat_mul_sse2(versor p, versor q, versor dest) {
   xp = glmm_load(p); /* 3 2 1 0 */
   xq = glmm_load(q);
 
-  r  = _mm_mul_ps(glmm_shuff1x(xp, 3), xq);
+  r  = _mm_mul_ps(glmm_splat(xp, 3), xq);
 
-  x0 = _mm_xor_ps(glmm_shuff1x(xp, 0), _mm_set_ps(-0.f, 0.f, -0.f, 0.f));
+  x0 = _mm_xor_ps(glmm_splat(xp, 0), _mm_set_ps(-0.f, 0.f, -0.f, 0.f));
   r  = _mm_add_ps(r, _mm_mul_ps(x0, glmm_shuff1(xq, 0, 1, 2, 3)));
 
-  x0 = _mm_xor_ps(glmm_shuff1x(xp, 1), _mm_set_ps(-0.f, -0.f, 0.f, 0.f));
+  x0 = _mm_xor_ps(glmm_splat(xp, 1), _mm_set_ps(-0.f, -0.f, 0.f, 0.f));
   r  = _mm_add_ps(r, _mm_mul_ps(x0, glmm_shuff1(xq, 1, 0, 3, 2)));
 
-  x0 = _mm_xor_ps(glmm_shuff1x(xp, 2), _mm_set_ps(-0.f, 0.f, 0.f, -0.f));
+  x0 = _mm_xor_ps(glmm_splat(xp, 2), _mm_set_ps(-0.f, 0.f, 0.f, -0.f));
   r  = _mm_add_ps(r, _mm_mul_ps(x0, glmm_shuff1(xq, 2, 3, 0, 1)));
 
   glmm_store(dest, r);
diff --git a/include/cglm/simd/x86.h b/include/cglm/simd/x86.h
index 5848b7b..43cb413 100644
--- a/include/cglm/simd/x86.h
+++ b/include/cglm/simd/x86.h
@@ -27,7 +27,11 @@
        _mm_shuffle_ps(xmm, xmm, _MM_SHUFFLE(z, y, x, w))
 #endif
 
+#define glmm_splat(x, lane) glmm_shuff1(x, lane, lane, lane, lane)
+
+/* glmm_shuff1x() is DEPRECATED!, use glmm_splat() */
 #define glmm_shuff1x(xmm, x) glmm_shuff1(xmm, x, x, x, x)
+
 #define glmm_shuff2(a, b, z0, y0, x0, w0, z1, y1, x1, w1)                     \
      glmm_shuff1(_mm_shuffle_ps(a, b, _MM_SHUFFLE(z0, y0, x0, w0)),           \
                  z1, y1, x1, w1)
@@ -89,7 +93,7 @@ glmm_vhmin(__m128 v) {
   __m128 x0, x1, x2;
   x0 = _mm_movehl_ps(v, v);     /* [2, 3, 2, 3] */
   x1 = _mm_min_ps(x0, v);       /* [0|2, 1|3, 2|2, 3|3] */
-  x2 = glmm_shuff1x(x1, 1);     /* [1|3, 1|3, 1|3, 1|3] */
+  x2 = glmm_splat(x1, 1);       /* [1|3, 1|3, 1|3, 1|3] */
   return _mm_min_ss(x1, x2);
 }
 
@@ -105,7 +109,7 @@ glmm_vhmax(__m128 v) {
   __m128 x0, x1, x2;
   x0 = _mm_movehl_ps(v, v);     /* [2, 3, 2, 3] */
   x1 = _mm_max_ps(x0, v);       /* [0|2, 1|3, 2|2, 3|3] */
-  x2 = glmm_shuff1x(x1, 1);     /* [1|3, 1|3, 1|3, 1|3] */
+  x2 = glmm_splat(x1, 1);       /* [1|3, 1|3, 1|3, 1|3] */
   return _mm_max_ss(x1, x2);
 }
 
diff --git a/include/cglm/vec4-ext.h b/include/cglm/vec4-ext.h
index dfc6a4e..e4e20cb 100644
--- a/include/cglm/vec4-ext.h
+++ b/include/cglm/vec4-ext.h
@@ -224,10 +224,10 @@ glm_vec4_sign(vec4 v, vec4 dest) {
 
   x0 = glmm_load(v);
   x1 = _mm_set_ps(0.0f, 0.0f, 1.0f, -1.0f);
-  x2 = glmm_shuff1x(x1, 2);
+  x2 = glmm_splat(x1, 2);
 
-  x3 = _mm_and_ps(_mm_cmpgt_ps(x0, x2), glmm_shuff1x(x1, 1));
-  x4 = _mm_and_ps(_mm_cmplt_ps(x0, x2), glmm_shuff1x(x1, 0));
+  x3 = _mm_and_ps(_mm_cmpgt_ps(x0, x2), glmm_splat(x1, 1));
+  x4 = _mm_and_ps(_mm_cmplt_ps(x0, x2), glmm_splat(x1, 0));
 
   glmm_store(dest, _mm_or_ps(x3, x4));
 #else

From 63988a515ce90143892dd5dad77de718f9b0acae Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Fri, 23 Apr 2021 23:32:21 +0300
Subject: [PATCH 02/16] glmm: new defines for splat

---
 include/cglm/simd/arm.h | 8 ++++++++
 include/cglm/simd/x86.h | 8 ++++++++
 2 files changed, 16 insertions(+)

diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index 1153694..4c8a3b4 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -13,6 +13,14 @@
 #define glmm_load(p)      vld1q_f32(p)
 #define glmm_store(p, a)  vst1q_f32(p, a)
 
+#define glmm_set1(x) vdupq_n_f32(x)
+#define glmm_128     float32x4_t
+
+#define glmm_splat_x(x) vdupq_lane_f32(vget_low_f32(x),  0)
+#define glmm_splat_y(x) vdupq_lane_f32(vget_low_f32(x),  1)
+#define glmm_splat_z(x) vdupq_lane_f32(vget_high_f32(x), 0)
+#define glmm_splat_w(x) vdupq_lane_f32(vget_high_f32(x), 1)
+
 static inline
 float32x4_t
 glmm_abs(float32x4_t v) {
diff --git a/include/cglm/simd/x86.h b/include/cglm/simd/x86.h
index 43cb413..df32491 100644
--- a/include/cglm/simd/x86.h
+++ b/include/cglm/simd/x86.h
@@ -18,6 +18,9 @@
 #  define glmm_store(p, a)  _mm_store_ps(p, a)
 #endif
 
+#define glmm_set1(x) _mm_set1_ps(x)
+#define glmm_128     __m128
+
 #ifdef CGLM_USE_INT_DOMAIN
 #  define glmm_shuff1(xmm, z, y, x, w)                                        \
      _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(xmm),                \
@@ -29,6 +32,11 @@
 
 #define glmm_splat(x, lane) glmm_shuff1(x, lane, lane, lane, lane)
 
+#define glmm_splat_x(x) glmm_splat(x, 0)
+#define glmm_splat_y(x) glmm_splat(x, 1)
+#define glmm_splat_z(x) glmm_splat(x, 2)
+#define glmm_splat_w(x) glmm_splat(x, 3)
+
 /* glmm_shuff1x() is DEPRECATED!, use glmm_splat() */
 #define glmm_shuff1x(xmm, x) glmm_shuff1(xmm, x, x, x, x)
 

From 7d5c4da7cf828960dd343abfaaa95eddea9923db Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Fri, 23 Apr 2021 23:32:48 +0300
Subject: [PATCH 03/16] optimize translate functions

---
 include/cglm/affine.h | 57 ++++++++++++++++---------------------------
 1 file changed, 21 insertions(+), 36 deletions(-)

diff --git a/include/cglm/affine.h b/include/cglm/affine.h
index 15f2853..504e13b 100644
--- a/include/cglm/affine.h
+++ b/include/cglm/affine.h
@@ -50,26 +50,23 @@
 CGLM_INLINE
 void
 glm_translate(mat4 m, vec3 v) {
-#if defined( __SSE__ ) || defined( __SSE2__ )
+#if defined(CGLM_SIMD)
+  glmm_128 v0, m0, m1, m2, m3;
+
+  m0 = glmm_load(m[0]);
+  m1 = glmm_load(m[1]);
+  m2 = glmm_load(m[2]);
+  m3 = glmm_load(m[3]);
+  v0 = glmm_load(v);
+
   glmm_store(m[3],
-             _mm_add_ps(_mm_add_ps(_mm_mul_ps(glmm_load(m[0]),
-                                              _mm_set1_ps(v[0])),
-                                   _mm_mul_ps(glmm_load(m[1]),
-                                              _mm_set1_ps(v[1]))),
-                        _mm_add_ps(_mm_mul_ps(glmm_load(m[2]),
-                                              _mm_set1_ps(v[2])),
-                                   glmm_load(m[3]))))
-  ;
+             glmm_fmadd(m0, glmm_splat_x(v0),
+                        glmm_fmadd(m1, glmm_splat_y(v0),
+                                   glmm_fmadd(m2, glmm_splat_z(v0), m3))));
 #else
-  vec4 v1, v2, v3;
-
-  glm_vec4_scale(m[0], v[0], v1);
-  glm_vec4_scale(m[1], v[1], v2);
-  glm_vec4_scale(m[2], v[2], v3);
-
-  glm_vec4_add(v1, m[3], m[3]);
-  glm_vec4_add(v2, m[3], m[3]);
-  glm_vec4_add(v3, m[3], m[3]);
+  glm_vec4_muladds(m[0], v[0], m[3]);
+  glm_vec4_muladds(m[1], v[1], m[3]);
+  glm_vec4_muladds(m[2], v[2], m[3]);
 #endif
 }
 
@@ -99,12 +96,8 @@ glm_translate_to(mat4 m, vec3 v, mat4 dest) {
 CGLM_INLINE
 void
 glm_translate_x(mat4 m, float x) {
-#if defined( __SSE__ ) || defined( __SSE2__ )
-  glmm_store(m[3],
-             _mm_add_ps(_mm_mul_ps(glmm_load(m[0]),
-                                   _mm_set1_ps(x)),
-                        glmm_load(m[3])))
-  ;
+#if defined(CGLM_SIMD)
+  glmm_store(m[3], glmm_fmadd(glmm_load(m[0]), glmm_set1(x), glmm_load(m[3])));
 #else
   vec4 v1;
   glm_vec4_scale(m[0], x, v1);
@@ -121,12 +114,8 @@ glm_translate_x(mat4 m, float x) {
 CGLM_INLINE
 void
 glm_translate_y(mat4 m, float y) {
-#if defined( __SSE__ ) || defined( __SSE2__ )
-  glmm_store(m[3],
-             _mm_add_ps(_mm_mul_ps(glmm_load(m[1]),
-                                   _mm_set1_ps(y)),
-                        glmm_load(m[3])))
-  ;
+#if defined(CGLM_SIMD)
+  glmm_store(m[3], glmm_fmadd(glmm_load(m[1]), glmm_set1(y), glmm_load(m[3])));
 #else
   vec4 v1;
   glm_vec4_scale(m[1], y, v1);
@@ -143,12 +132,8 @@ glm_translate_y(mat4 m, float y) {
 CGLM_INLINE
 void
 glm_translate_z(mat4 m, float z) {
-#if defined( __SSE__ ) || defined( __SSE2__ )
-  glmm_store(m[3],
-             _mm_add_ps(_mm_mul_ps(glmm_load(m[2]),
-                                   _mm_set1_ps(z)),
-                        glmm_load(m[3])))
-  ;
+#if defined(CGLM_SIMD)
+  glmm_store(m[3], glmm_fmadd(glmm_load(m[2]), glmm_set1(z), glmm_load(m[3])));
 #else
   vec4 v1;
   glm_vec4_scale(m[2], z, v1);

From f303984aad6906a823eeda4ab021467996958fa2 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Fri, 23 Apr 2021 23:34:36 +0300
Subject: [PATCH 04/16] use unified glmm api for vec4

---
 include/cglm/vec4.h | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/include/cglm/vec4.h b/include/cglm/vec4.h
index 2453b1b..ce3f017 100644
--- a/include/cglm/vec4.h
+++ b/include/cglm/vec4.h
@@ -590,10 +590,8 @@ glm_vec4_muladd(vec4 a, vec4 b, vec4 dest) {
 CGLM_INLINE
 void
 glm_vec4_muladds(vec4 a, float s, vec4 dest) {
-#if defined( __SSE__ ) || defined( __SSE2__ )
-  glmm_store(dest, glmm_fmadd(glmm_load(a), _mm_set1_ps(s), glmm_load(dest)));
-#elif defined(CGLM_NEON_FP)
-  glmm_store(dest, glmm_fmadd(glmm_load(a), vdupq_n_f32(s), glmm_load(dest)));
+#if defined(CGLM_SIMD)
+  glmm_store(dest, glmm_fmadd(glmm_load(a), glmm_set1(s), glmm_load(dest)));
 #else
   dest[0] += a[0] * s;
   dest[1] += a[1] * s;

From 65292a94a6d3ab7b0d0b2d91668a69a588feb0c1 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sat, 24 Apr 2021 00:00:00 +0300
Subject: [PATCH 05/16] swizzling functions for NEON

---
 CREDITS                 |  5 +++++
 include/cglm/simd/arm.h | 40 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 45 insertions(+)

diff --git a/CREDITS b/CREDITS
index daeb27e..94d9d40 100644
--- a/CREDITS
+++ b/CREDITS
@@ -74,3 +74,8 @@ Link to paper: http://webserver2.tecgraf.puc-rio.br/~mgattass/cg/trbRR/Fast%20Mi
 
 14. ARM NEON: Matrix Vector Multiplication
 https://stackoverflow.com/a/57793352/2676533
+
+15. ARM NEON Vector Swizzling and Permute
+
+https://stackoverflow.com/questions/32536265/how-to-convert-mm-shuffle-ps-sse-intrinsic-to-neon-intrinsic
+http://github.com/microsoft/DirectXMath
diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index 4c8a3b4..1a5adc9 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -21,6 +21,46 @@
 #define glmm_splat_z(x) vdupq_lane_f32(vget_high_f32(x), 0)
 #define glmm_splat_w(x) vdupq_lane_f32(vget_high_f32(x), 1)
 
+#define SWIZZLE(NAME) \
+  static inline float32x4_t NAME(glmm_128 v)
+
+/* TODO:
+ * I'm not sure if glmm_xxxx(), glmm_yyyy()... is better than glmm_0000()...
+ */
+
+SWIZZLE(glmm_0000) { return vdupq_lane_f32(vget_low_f32(v),  0); }
+SWIZZLE(glmm_1111) { return vdupq_lane_f32(vget_low_f32(v),  1); }
+SWIZZLE(glmm_2222) { return vdupq_lane_f32(vget_high_f32(v), 0); }
+SWIZZLE(glmm_3333) { return vdupq_lane_f32(vget_high_f32(v), 1); }
+
+SWIZZLE(glmm_1032) { return vrev64q_f32(v); }
+
+SWIZZLE(glmm_0101) { float32x2_t vt = vget_low_f32(v);              return vcombine_f32(vt, vt); }
+SWIZZLE(glmm_2323) { float32x2_t vt = vget_high_f32(v);             return vcombine_f32(vt, vt); }
+SWIZZLE(glmm_1010) { float32x2_t vt = vrev64_f32(vget_low_f32(v));  return vcombine_f32(vt, vt); }
+SWIZZLE(glmm_3232) { float32x2_t vt = vrev64_f32(vget_high_f32(v)); return vcombine_f32(vt, vt); }
+
+SWIZZLE(glmm_0132) { return vcombine_f32(vget_low_f32(v), vrev64_f32(vget_high_f32(v))); }
+SWIZZLE(glmm_1023) { return vcombine_f32(vrev64_f32(vget_low_f32(v)), vget_high_f32(v)); }
+SWIZZLE(glmm_2310) { return vcombine_f32(vget_high_f32(v), vrev64_f32(vget_low_f32(v))); }
+SWIZZLE(glmm_3201) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vget_low_f32(v)); }
+SWIZZLE(glmm_3210) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vrev64_f32(vget_low_f32(v))); }
+
+SWIZZLE(glmm_0022) { return vtrnq_f32(v, v).val[0]; }
+SWIZZLE(glmm_1133) { return vtrnq_f32(v, v).val[1]; }
+
+SWIZZLE(glmm_0011) { return vzipq_f32(v, v).val[0]; }
+SWIZZLE(glmm_2233) { return vzipq_f32(v, v).val[1]; }
+
+SWIZZLE(glmm_0202) { return vuzpq_f32(v, v).val[0]; }
+SWIZZLE(glmm_1313) { return vuzpq_f32(v, v).val[1]; }
+
+SWIZZLE(glmm_1230) { return vextq_f32(v, v, 1);     }
+SWIZZLE(glmm_2301) { return vextq_f32(v, v, 2);     }
+SWIZZLE(glmm_3012) { return vextq_f32(v, v, 3);     }
+
+#undef SWIZZLE
+
 static inline
 float32x4_t
 glmm_abs(float32x4_t v) {

From 059e5010e6b5ce94e8ca5e9f237063fb8c9a6071 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sat, 24 Apr 2021 00:44:07 +0300
Subject: [PATCH 06/16] arm: define CGLM_ARM64 for identify arm64

---
 include/cglm/simd/arm.h | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index 1a5adc9..17ce15e 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -10,6 +10,10 @@
 #include "intrin.h"
 #ifdef CGLM_SIMD_ARM
 
+#if defined(_M_ARM64) || defined(_M_HYBRID_X86_ARM64) || defined(_M_ARM64EC) || defined(__aarch64__)
+# define CGLM_ARM64 1
+#endif
+
 #define glmm_load(p)      vld1q_f32(p)
 #define glmm_store(p, a)  vst1q_f32(p, a)
 
@@ -70,7 +74,7 @@ glmm_abs(float32x4_t v) {
 static inline
 float
 glmm_hadd(float32x4_t v) {
-#if defined(__aarch64__)
+#if CGLM_ARM64
   return vaddvq_f32(v);
 #else
   v = vaddq_f32(v, vrev64q_f32(v));
@@ -130,7 +134,7 @@ glmm_norm_inf(float32x4_t a) {
 static inline
 float32x4_t
 glmm_fmadd(float32x4_t a, float32x4_t b, float32x4_t c) {
-#if defined(__aarch64__)
+#if CGLM_ARM64
   return vfmaq_f32(c, a, b);
 #else
   return vmlaq_f32(c, a, b);
@@ -140,7 +144,7 @@ glmm_fmadd(float32x4_t a, float32x4_t b, float32x4_t c) {
 static inline
 float32x4_t
 glmm_fnmadd(float32x4_t a, float32x4_t b, float32x4_t c) {
-#if defined(__aarch64__)
+#if CGLM_ARM64
   return vfmsq_f32(c, a, b);
 #else
   return vmlsq_f32(c, a, b);
@@ -150,7 +154,7 @@ glmm_fnmadd(float32x4_t a, float32x4_t b, float32x4_t c) {
 static inline
 float32x4_t
 glmm_fmsub(float32x4_t a, float32x4_t b, float32x4_t c) {
-#if defined(__aarch64__)
+#if CGLM_ARM64
   return vfmsq_f32(c, a, b);
 #else
   return vmlsq_f32(c, a, b);

From d13842e7dece4de1877344b14eee945f91a96c3d Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sat, 24 Apr 2021 00:51:09 +0300
Subject: [PATCH 07/16] arm: optimize vec4 div with NEON

---
 CREDITS                 |  4 ++++
 include/cglm/simd/arm.h | 17 +++++++++++++++++
 include/cglm/simd/x86.h |  6 ++++++
 include/cglm/vec4.h     |  4 ++--
 4 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/CREDITS b/CREDITS
index 94d9d40..96341fc 100644
--- a/CREDITS
+++ b/CREDITS
@@ -79,3 +79,7 @@ https://stackoverflow.com/a/57793352/2676533
 
 https://stackoverflow.com/questions/32536265/how-to-convert-mm-shuffle-ps-sse-intrinsic-to-neon-intrinsic
 http://github.com/microsoft/DirectXMath
+
+16. ARM NEON Div
+
+http://github.com/microsoft/DirectXMath
diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index 17ce15e..c980e54 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -131,6 +131,23 @@ glmm_norm_inf(float32x4_t a) {
   return glmm_hmax(glmm_abs(a));
 }
 
+static inline
+float32x4_t
+glmm_div(float32x4_t a, float32x4_t b) {
+#if CGLM_ARM641
+  return vdivq_f32(a, b);
+#else
+  /* 2 iterations of Newton-Raphson refinement of reciprocal */
+  float32x4_t r0, r1;
+  r0 = vrecpeq_f32(b);
+  r1 = vrecpsq_f32(r0, b);
+  r0 = vmulq_f32(r1, r0);
+  r1 = vrecpsq_f32(r0, b);
+  r0 = vmulq_f32(r1, r0);
+  return vmulq_f32(a, r0);
+#endif
+}
+
 static inline
 float32x4_t
 glmm_fmadd(float32x4_t a, float32x4_t b, float32x4_t c) {
diff --git a/include/cglm/simd/x86.h b/include/cglm/simd/x86.h
index df32491..dbbd0f8 100644
--- a/include/cglm/simd/x86.h
+++ b/include/cglm/simd/x86.h
@@ -209,6 +209,12 @@ glmm_store3(float v[3], __m128 vx) {
   _mm_store_ss(&v[2], glmm_shuff1(vx, 2, 2, 2, 2));
 }
 
+static inline
+__m128
+glmm_div(__m128 a, __m128 b) {
+  return _mm_div_ps(a, b);
+}
+
 /* enable FMA macro for MSVC? */
 #if defined(_MSC_VER) && !defined(__FMA__) && defined(__AVX2__)
 #  define __FMA__ 1
diff --git a/include/cglm/vec4.h b/include/cglm/vec4.h
index ce3f017..8e95ec5 100644
--- a/include/cglm/vec4.h
+++ b/include/cglm/vec4.h
@@ -473,8 +473,8 @@ glm_vec4_scale_as(vec4 v, float s, vec4 dest) {
 CGLM_INLINE
 void
 glm_vec4_div(vec4 a, vec4 b, vec4 dest) {
-#if defined( __SSE__ ) || defined( __SSE2__ )
-  glmm_store(dest, _mm_div_ps(glmm_load(a), glmm_load(b)));
+#if defined(CGLM_SIMD)
+  glmm_store(dest, glmm_div(glmm_load(a), glmm_load(b)));
 #else
   dest[0] = a[0] / b[0];
   dest[1] = a[1] / b[1];

From ce9e5f557530a055f499ff47dbcced800d7ec6ea Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sat, 24 Apr 2021 01:38:04 +0300
Subject: [PATCH 08/16] arm: update glmm swizzling func names

---
 include/cglm/simd/arm.h | 42 ++++++++++++++++++++++-------------------
 1 file changed, 23 insertions(+), 19 deletions(-)

diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index c980e54..791f4cb 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -32,36 +32,40 @@
  * I'm not sure if glmm_xxxx(), glmm_yyyy()... is better than glmm_0000()...
  */
 
+/* Memory layout       Register layout (glmm)
+   0 1 2 3        ->   3 2 1 0
+ */
+
 SWIZZLE(glmm_0000) { return vdupq_lane_f32(vget_low_f32(v),  0); }
 SWIZZLE(glmm_1111) { return vdupq_lane_f32(vget_low_f32(v),  1); }
 SWIZZLE(glmm_2222) { return vdupq_lane_f32(vget_high_f32(v), 0); }
 SWIZZLE(glmm_3333) { return vdupq_lane_f32(vget_high_f32(v), 1); }
 
-SWIZZLE(glmm_1032) { return vrev64q_f32(v); }
+SWIZZLE(glmm_2301) { return vrev64q_f32(v); }
 
-SWIZZLE(glmm_0101) { float32x2_t vt = vget_low_f32(v);              return vcombine_f32(vt, vt); }
-SWIZZLE(glmm_2323) { float32x2_t vt = vget_high_f32(v);             return vcombine_f32(vt, vt); }
-SWIZZLE(glmm_1010) { float32x2_t vt = vrev64_f32(vget_low_f32(v));  return vcombine_f32(vt, vt); }
-SWIZZLE(glmm_3232) { float32x2_t vt = vrev64_f32(vget_high_f32(v)); return vcombine_f32(vt, vt); }
+SWIZZLE(glmm_1010) { float32x2_t vt = vget_low_f32(v);              return vcombine_f32(vt, vt); }
+SWIZZLE(glmm_3232) { float32x2_t vt = vget_high_f32(v);             return vcombine_f32(vt, vt); }
+SWIZZLE(glmm_0101) { float32x2_t vt = vrev64_f32(vget_low_f32(v));  return vcombine_f32(vt, vt); }
+SWIZZLE(glmm_2323) { float32x2_t vt = vrev64_f32(vget_high_f32(v)); return vcombine_f32(vt, vt); }
 
-SWIZZLE(glmm_0132) { return vcombine_f32(vget_low_f32(v), vrev64_f32(vget_high_f32(v))); }
-SWIZZLE(glmm_1023) { return vcombine_f32(vrev64_f32(vget_low_f32(v)), vget_high_f32(v)); }
-SWIZZLE(glmm_2310) { return vcombine_f32(vget_high_f32(v), vrev64_f32(vget_low_f32(v))); }
-SWIZZLE(glmm_3201) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vget_low_f32(v)); }
-SWIZZLE(glmm_3210) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vrev64_f32(vget_low_f32(v))); }
+SWIZZLE(glmm_2310) { return vcombine_f32(vget_low_f32(v), vrev64_f32(vget_high_f32(v))); }
+SWIZZLE(glmm_3201) { return vcombine_f32(vrev64_f32(vget_low_f32(v)), vget_high_f32(v)); }
+SWIZZLE(glmm_0132) { return vcombine_f32(vget_high_f32(v), vrev64_f32(vget_low_f32(v))); }
+SWIZZLE(glmm_1023) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vget_low_f32(v)); }
+SWIZZLE(glmm_0123) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vrev64_f32(vget_low_f32(v))); }
 
-SWIZZLE(glmm_0022) { return vtrnq_f32(v, v).val[0]; }
-SWIZZLE(glmm_1133) { return vtrnq_f32(v, v).val[1]; }
+SWIZZLE(glmm_2200) { return vtrnq_f32(v, v).val[0]; }
+SWIZZLE(glmm_3311) { return vtrnq_f32(v, v).val[1]; }
 
-SWIZZLE(glmm_0011) { return vzipq_f32(v, v).val[0]; }
-SWIZZLE(glmm_2233) { return vzipq_f32(v, v).val[1]; }
+SWIZZLE(glmm_1100) { return vzipq_f32(v, v).val[0]; }
+SWIZZLE(glmm_3322) { return vzipq_f32(v, v).val[1]; }
 
-SWIZZLE(glmm_0202) { return vuzpq_f32(v, v).val[0]; }
-SWIZZLE(glmm_1313) { return vuzpq_f32(v, v).val[1]; }
+SWIZZLE(glmm_2020) { return vuzpq_f32(v, v).val[0]; }
+SWIZZLE(glmm_3131) { return vuzpq_f32(v, v).val[1]; }
 
-SWIZZLE(glmm_1230) { return vextq_f32(v, v, 1);     }
-SWIZZLE(glmm_2301) { return vextq_f32(v, v, 2);     }
-SWIZZLE(glmm_3012) { return vextq_f32(v, v, 3);     }
+SWIZZLE(glmm_0321) { return vextq_f32(v, v, 1);     }
+SWIZZLE(glmm_1032) { return vextq_f32(v, v, 2);     }
+SWIZZLE(glmm_2103) { return vextq_f32(v, v, 3);     }
 
 #undef SWIZZLE
 

From a111693b6b42f71e410d3edc0e420fb9d2bab5cb Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sat, 24 Apr 2021 15:45:36 +0300
Subject: [PATCH 09/16] arm, neon: implement mat4 determinant with neon

---
 include/cglm/mat4.h           |  2 +
 include/cglm/simd/arm.h       | 13 ++++++-
 include/cglm/simd/neon/mat4.h | 70 +++++++++++++++++++++++++++++++++++
 3 files changed, 84 insertions(+), 1 deletion(-)

diff --git a/include/cglm/mat4.h b/include/cglm/mat4.h
index c099574..697f510 100644
--- a/include/cglm/mat4.h
+++ b/include/cglm/mat4.h
@@ -562,6 +562,8 @@ float
 glm_mat4_det(mat4 mat) {
 #if defined( __SSE__ ) || defined( __SSE2__ )
   return glm_mat4_det_sse2(mat);
+#elif defined(CGLM_NEON_FP)
+  return glm_mat4_det_neon(mat);
 #else
   /* [square] det(A) = det(At) */
   float t[6];
diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index 791f4cb..cd4f2b8 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -69,12 +69,23 @@ SWIZZLE(glmm_2103) { return vextq_f32(v, v, 3);     }
 
 #undef SWIZZLE
 
+#define glmm_xor(a, b)                                                        \
+  vreinterpretq_f32_s32(veorq_s32(vreinterpretq_s32_f32(a),                   \
+                                  vreinterpretq_s32_f32(b)))
+
 static inline
 float32x4_t
 glmm_abs(float32x4_t v) {
   return vabsq_f32(v);
 }
 
+static inline
+float32x4_t
+glmm_vhadd(float32x4_t v) {
+  v = vaddq_f32(v, vrev64q_f32(v));
+  return vaddq_f32(v, vcombine_f32(vget_high_f32(v), vget_low_f32(v)));
+}
+
 static inline
 float
 glmm_hadd(float32x4_t v) {
@@ -138,7 +149,7 @@ glmm_norm_inf(float32x4_t a) {
 static inline
 float32x4_t
 glmm_div(float32x4_t a, float32x4_t b) {
-#if CGLM_ARM641
+#if CGLM_ARM64
   return vdivq_f32(a, b);
 #else
   /* 2 iterations of Newton-Raphson refinement of reciprocal */
diff --git a/include/cglm/simd/neon/mat4.h b/include/cglm/simd/neon/mat4.h
index 36d347b..ad32339 100644
--- a/include/cglm/simd/neon/mat4.h
+++ b/include/cglm/simd/neon/mat4.h
@@ -101,5 +101,75 @@ glm_mat4_mulv_neon(mat4 m, vec4 v, vec4 dest) {
   vst1q_f32(dest, l0);
 }
 
+CGLM_INLINE
+float
+glm_mat4_det_neon(mat4 mat) {
+  float32x4_t   r0, r1, r2, r3, x0, x1, x2;
+  float32x2_t   ij, op, mn, kl, nn, mm, jj, ii, gh, ef, t12, t34;
+  float32x4x2_t a1;
+  float32x4_t   x3 = { 0.f, -0.f, 0.f, -0.f };
+
+  /* 127 <- 0, [square] det(A) = det(At) */
+  r0 = glmm_load(mat[0]);              /* d c b a */
+  r1 = vrev64q_f32(glmm_load(mat[1])); /* g h e f */
+  r2 = vrev64q_f32(glmm_load(mat[2])); /* l k i j */
+  r3 = vrev64q_f32(glmm_load(mat[3])); /* o p m n */
+
+  gh = vget_high_f32(r1);
+  ef = vget_low_f32(r1);
+  kl = vget_high_f32(r2);
+  ij = vget_low_f32(r2);
+  op = vget_high_f32(r3);
+  mn = vget_low_f32(r3);
+  mm = vdup_lane_f32(mn, 1);
+  nn = vdup_lane_f32(mn, 0);
+  ii = vdup_lane_f32(ij, 1);
+  jj = vdup_lane_f32(ij, 0);
+  
+  /*
+   t[1] = j * p - n * l;
+   t[2] = j * o - n * k;
+   t[3] = i * p - m * l;
+   t[4] = i * o - m * k;
+   */
+  x0 = glmm_fnmadd(vcombine_f32(kl, kl), vcombine_f32(nn, mm),
+                   vmulq_f32(vcombine_f32(op, op), vcombine_f32(jj, ii)));
+
+  t12 = vget_low_f32(x0);
+  t34 = vget_high_f32(x0);
+  
+  /* 1 3 1 3 2 4 2 4 */
+  a1 = vuzpq_f32(x0, x0);
+  
+  /*
+   t[0] = k * p - o * l;
+   t[0] = k * p - o * l;
+   t[5] = i * n - m * j;
+   t[5] = i * n - m * j;
+   */
+  x1 = glmm_fnmadd(vcombine_f32(vdup_lane_f32(kl, 0), jj),
+                   vcombine_f32(vdup_lane_f32(op, 1), mm),
+                   vmulq_f32(vcombine_f32(vdup_lane_f32(op, 0), nn),
+                             vcombine_f32(vdup_lane_f32(kl, 1), ii)));
+
+  /*
+     a * (f * t[0] - g * t[1] + h * t[2])
+   - b * (e * t[0] - g * t[3] + h * t[4])
+   + c * (e * t[1] - f * t[3] + h * t[5])
+   - d * (e * t[2] - f * t[4] + g * t[5])
+   */
+  x2 = glmm_fnmadd(vcombine_f32(vdup_lane_f32(gh, 1), vdup_lane_f32(ef, 0)),
+                   vcombine_f32(vget_low_f32(a1.val[0]), t34),
+                   vmulq_f32(vcombine_f32(ef, vdup_lane_f32(ef, 1)),
+                             vcombine_f32(vget_low_f32(x1), t12)));
+
+  x2 = glmm_fmadd(vcombine_f32(vdup_lane_f32(gh, 0), gh),
+                  vcombine_f32(vget_low_f32(a1.val[1]), vget_high_f32(x1)), x2);
+
+  x2 = glmm_xor(x2, x3);
+
+  return glmm_hadd(vmulq_f32(x2, r0));
+}
+
 #endif
 #endif /* cglm_mat4_neon_h */

From afac887850dfc03353acc6dbc5f80d74f555284e Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sat, 24 Apr 2021 17:54:01 +0300
Subject: [PATCH 10/16] arm, neon: implement mat4 inv with neon

---
 include/cglm/mat4.h           |   2 +
 include/cglm/simd/neon/mat4.h | 132 ++++++++++++++++++++++++++++++++++
 2 files changed, 134 insertions(+)

diff --git a/include/cglm/mat4.h b/include/cglm/mat4.h
index 697f510..0b64db8 100644
--- a/include/cglm/mat4.h
+++ b/include/cglm/mat4.h
@@ -597,6 +597,8 @@ void
 glm_mat4_inv(mat4 mat, mat4 dest) {
 #if defined( __SSE__ ) || defined( __SSE2__ )
   glm_mat4_inv_sse2(mat, dest);
+#elif defined(CGLM_NEON_FP)
+  glm_mat4_inv_neon(mat, dest);
 #else
   float t[6];
   float det;
diff --git a/include/cglm/simd/neon/mat4.h b/include/cglm/simd/neon/mat4.h
index ad32339..b992acc 100644
--- a/include/cglm/simd/neon/mat4.h
+++ b/include/cglm/simd/neon/mat4.h
@@ -171,5 +171,137 @@ glm_mat4_det_neon(mat4 mat) {
   return glmm_hadd(vmulq_f32(x2, r0));
 }
 
+CGLM_INLINE
+void
+glm_mat4_inv_neon(mat4 mat, mat4 dest) {
+  float32x4_t   r0, r1, r2, r3,
+                v0, v1, v2, v3,
+                t0, t1, t2, t3, t4, t5,
+                x0, x1, x2, x3, x4, x5, x6, x7, x8;
+  float32x4x2_t a1;
+  float32x2_t   lp, ko, hg, jn, im, fe, ae, bf, cg, dh;
+  float32x4_t   x9 = { -0.f, 0.f, -0.f, 0.f };
+
+  x8 = vrev64q_f32(x9);
+
+  /* 127 <- 0 */
+  r0 = glmm_load(mat[0]); /* d c b a */
+  r1 = glmm_load(mat[1]); /* h g f e */
+  r2 = glmm_load(mat[2]); /* l k j i */
+  r3 = glmm_load(mat[3]); /* p o n m */
+  
+  /* l p k o, j n i m */
+  a1  = vzipq_f32(r3, r2);
+  
+  jn  = vget_high_f32(a1.val[0]);
+  im  = vget_low_f32(a1.val[0]);
+  lp  = vget_high_f32(a1.val[1]);
+  ko  = vget_low_f32(a1.val[1]);
+  hg  = vget_high_f32(r1);
+
+  x1  = vcombine_f32(vdup_lane_f32(lp, 0), lp);                   /* l p p p */
+  x2  = vcombine_f32(vdup_lane_f32(ko, 0), ko);                   /* k o o o */
+  x0  = vcombine_f32(vdup_lane_f32(lp, 1), vdup_lane_f32(hg, 1)); /* h h l l */
+  x3  = vcombine_f32(vdup_lane_f32(ko, 1), vdup_lane_f32(hg, 0)); /* g g k k */
+  
+  /* t1[0] = k * p - o * l;
+     t1[0] = k * p - o * l;
+     t2[0] = g * p - o * h;
+     t3[0] = g * l - k * h; */
+  t0 = glmm_fnmadd(x2, x0, vmulq_f32(x3, x1));
+
+  fe = vget_low_f32(r1);
+  x4 = vcombine_f32(vdup_lane_f32(jn, 0), jn);                   /* j n n n */
+  x5 = vcombine_f32(vdup_lane_f32(jn, 1), vdup_lane_f32(fe, 1)); /* f f j j */
+  
+  /* t1[1] = j * p - n * l;
+     t1[1] = j * p - n * l;
+     t2[1] = f * p - n * h;
+     t3[1] = f * l - j * h; */
+   t1 = glmm_fnmadd(x4, x0, vmulq_f32(x5, x1));
+  
+  /* t1[2] = j * o - n * k
+     t1[2] = j * o - n * k;
+     t2[2] = f * o - n * g;
+     t3[2] = f * k - j * g; */
+  t2 = glmm_fnmadd(x4, x3, vmulq_f32(x5, x2));
+  
+  x6 = vcombine_f32(vdup_lane_f32(im, 1), vdup_lane_f32(fe, 0)); /* e e i i */
+  x7 = vcombine_f32(vdup_lane_f32(im, 0), im);                   /* i m m m */
+  
+  /* t1[3] = i * p - m * l;
+     t1[3] = i * p - m * l;
+     t2[3] = e * p - m * h;
+     t3[3] = e * l - i * h; */
+  t3 = glmm_fnmadd(x7, x0, vmulq_f32(x6, x1));
+  
+  /* t1[4] = i * o - m * k;
+     t1[4] = i * o - m * k;
+     t2[4] = e * o - m * g;
+     t3[4] = e * k - i * g; */
+  t4 = glmm_fnmadd(x7, x3, vmulq_f32(x6, x2));
+  
+  /* t1[5] = i * n - m * j;
+     t1[5] = i * n - m * j;
+     t2[5] = e * n - m * f;
+     t3[5] = e * j - i * f; */
+  t5 = glmm_fnmadd(x7, x5, vmulq_f32(x6, x4));
+  
+  /* h d f b, g c e a */
+  a1 = vtrnq_f32(r0, r1);
+  
+  x4 = vrev64q_f32(a1.val[0]); /* c g a e */
+  x5 = vrev64q_f32(a1.val[1]); /* d h b f */
+
+  ae = vget_low_f32(x4);
+  cg = vget_high_f32(x4);
+  bf = vget_low_f32(x5);
+  dh = vget_high_f32(x5);
+  
+  x0 = vcombine_f32(ae, vdup_lane_f32(ae, 1)); /* a a a e */
+  x1 = vcombine_f32(bf, vdup_lane_f32(bf, 1)); /* b b b f */
+  x2 = vcombine_f32(cg, vdup_lane_f32(cg, 1)); /* c c c g */
+  x3 = vcombine_f32(dh, vdup_lane_f32(dh, 1)); /* d d d h */
+  
+  /*
+   dest[0][0] =  f * t1[0] - g * t1[1] + h * t1[2];
+   dest[0][1] =-(b * t1[0] - c * t1[1] + d * t1[2]);
+   dest[0][2] =  b * t2[0] - c * t2[1] + d * t2[2];
+   dest[0][3] =-(b * t3[0] - c * t3[1] + d * t3[2]); */
+  v0 = glmm_xor(glmm_fmadd(x3, t2, glmm_fnmadd(x2, t1, vmulq_f32(x1, t0))), x8);
+  
+  /*
+   dest[2][0] =  e * t1[1] - f * t1[3] + h * t1[5];
+   dest[2][1] =-(a * t1[1] - b * t1[3] + d * t1[5]);
+   dest[2][2] =  a * t2[1] - b * t2[3] + d * t2[5];
+   dest[2][3] =-(a * t3[1] - b * t3[3] + d * t3[5]);*/
+  v2 = glmm_xor(glmm_fmadd(x3, t5, glmm_fnmadd(x1, t3, vmulq_f32(x0, t1))), x8);
+
+  /*
+   dest[1][0] =-(e * t1[0] - g * t1[3] + h * t1[4]);
+   dest[1][1] =  a * t1[0] - c * t1[3] + d * t1[4];
+   dest[1][2] =-(a * t2[0] - c * t2[3] + d * t2[4]);
+   dest[1][3] =  a * t3[0] - c * t3[3] + d * t3[4]; */
+  v1 = glmm_xor(glmm_fmadd(x3, t4, glmm_fnmadd(x2, t3, vmulq_f32(x0, t0))), x9);
+  
+  /*
+   dest[3][0] =-(e * t1[2] - f * t1[4] + g * t1[5]);
+   dest[3][1] =  a * t1[2] - b * t1[4] + c * t1[5];
+   dest[3][2] =-(a * t2[2] - b * t2[4] + c * t2[5]);
+   dest[3][3] =  a * t3[2] - b * t3[4] + c * t3[5]; */
+  v3 = glmm_xor(glmm_fmadd(x2, t5, glmm_fnmadd(x1, t4, vmulq_f32(x0, t2))), x9);
+
+  /* determinant */
+  x0 = vcombine_f32(vget_low_f32(vzipq_f32(v0, v1).val[0]),
+                    vget_low_f32(vzipq_f32(v2, v3).val[0]));
+
+  x0 = glmm_div(glmm_set1(1.0f), glmm_vhadd(vmulq_f32(x0, r0)));
+
+  glmm_store(dest[0], vmulq_f32(v0, x0));
+  glmm_store(dest[1], vmulq_f32(v1, x0));
+  glmm_store(dest[2], vmulq_f32(v2, x0));
+  glmm_store(dest[3], vmulq_f32(v3, x0));
+}
+
 #endif
 #endif /* cglm_mat4_neon_h */

From 92151c63287d2e2893ad5428ddde61d4e460cd74 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sat, 24 Apr 2021 18:02:47 +0300
Subject: [PATCH 11/16] arm, neon: use div instead of mul by  1 / det for
 mat4_inv

---
 include/cglm/simd/neon/mat4.h | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/include/cglm/simd/neon/mat4.h b/include/cglm/simd/neon/mat4.h
index b992acc..2fecebe 100644
--- a/include/cglm/simd/neon/mat4.h
+++ b/include/cglm/simd/neon/mat4.h
@@ -295,12 +295,21 @@ glm_mat4_inv_neon(mat4 mat, mat4 dest) {
   x0 = vcombine_f32(vget_low_f32(vzipq_f32(v0, v1).val[0]),
                     vget_low_f32(vzipq_f32(v2, v3).val[0]));
 
+  /*
   x0 = glmm_div(glmm_set1(1.0f), glmm_vhadd(vmulq_f32(x0, r0)));
 
   glmm_store(dest[0], vmulq_f32(v0, x0));
   glmm_store(dest[1], vmulq_f32(v1, x0));
   glmm_store(dest[2], vmulq_f32(v2, x0));
   glmm_store(dest[3], vmulq_f32(v3, x0));
+  */
+
+  x0 = glmm_vhadd(vmulq_f32(x0, r0));
+
+  glmm_store(dest[0], glmm_div(v0, x0));
+  glmm_store(dest[1], glmm_div(v1, x0));
+  glmm_store(dest[2], glmm_div(v2, x0));
+  glmm_store(dest[3], glmm_div(v3, x0));
 }
 
 #endif

From 1fb941a41b2a6e0627633185132afd590eacf6c9 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sun, 25 Apr 2021 02:35:55 +0300
Subject: [PATCH 12/16] drop swizzling helpers fro now for simplicity

---
 CREDITS                 |  5 -----
 include/cglm/simd/arm.h | 44 -----------------------------------------
 2 files changed, 49 deletions(-)

diff --git a/CREDITS b/CREDITS
index 96341fc..c388348 100644
--- a/CREDITS
+++ b/CREDITS
@@ -75,11 +75,6 @@ Link to paper: http://webserver2.tecgraf.puc-rio.br/~mgattass/cg/trbRR/Fast%20Mi
 14. ARM NEON: Matrix Vector Multiplication
 https://stackoverflow.com/a/57793352/2676533
 
-15. ARM NEON Vector Swizzling and Permute
-
-https://stackoverflow.com/questions/32536265/how-to-convert-mm-shuffle-ps-sse-intrinsic-to-neon-intrinsic
-http://github.com/microsoft/DirectXMath
-
 16. ARM NEON Div
 
 http://github.com/microsoft/DirectXMath
diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index cd4f2b8..a7009e5 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -25,50 +25,6 @@
 #define glmm_splat_z(x) vdupq_lane_f32(vget_high_f32(x), 0)
 #define glmm_splat_w(x) vdupq_lane_f32(vget_high_f32(x), 1)
 
-#define SWIZZLE(NAME) \
-  static inline float32x4_t NAME(glmm_128 v)
-
-/* TODO:
- * I'm not sure if glmm_xxxx(), glmm_yyyy()... is better than glmm_0000()...
- */
-
-/* Memory layout       Register layout (glmm)
-   0 1 2 3        ->   3 2 1 0
- */
-
-SWIZZLE(glmm_0000) { return vdupq_lane_f32(vget_low_f32(v),  0); }
-SWIZZLE(glmm_1111) { return vdupq_lane_f32(vget_low_f32(v),  1); }
-SWIZZLE(glmm_2222) { return vdupq_lane_f32(vget_high_f32(v), 0); }
-SWIZZLE(glmm_3333) { return vdupq_lane_f32(vget_high_f32(v), 1); }
-
-SWIZZLE(glmm_2301) { return vrev64q_f32(v); }
-
-SWIZZLE(glmm_1010) { float32x2_t vt = vget_low_f32(v);              return vcombine_f32(vt, vt); }
-SWIZZLE(glmm_3232) { float32x2_t vt = vget_high_f32(v);             return vcombine_f32(vt, vt); }
-SWIZZLE(glmm_0101) { float32x2_t vt = vrev64_f32(vget_low_f32(v));  return vcombine_f32(vt, vt); }
-SWIZZLE(glmm_2323) { float32x2_t vt = vrev64_f32(vget_high_f32(v)); return vcombine_f32(vt, vt); }
-
-SWIZZLE(glmm_2310) { return vcombine_f32(vget_low_f32(v), vrev64_f32(vget_high_f32(v))); }
-SWIZZLE(glmm_3201) { return vcombine_f32(vrev64_f32(vget_low_f32(v)), vget_high_f32(v)); }
-SWIZZLE(glmm_0132) { return vcombine_f32(vget_high_f32(v), vrev64_f32(vget_low_f32(v))); }
-SWIZZLE(glmm_1023) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vget_low_f32(v)); }
-SWIZZLE(glmm_0123) { return vcombine_f32(vrev64_f32(vget_high_f32(v)), vrev64_f32(vget_low_f32(v))); }
-
-SWIZZLE(glmm_2200) { return vtrnq_f32(v, v).val[0]; }
-SWIZZLE(glmm_3311) { return vtrnq_f32(v, v).val[1]; }
-
-SWIZZLE(glmm_1100) { return vzipq_f32(v, v).val[0]; }
-SWIZZLE(glmm_3322) { return vzipq_f32(v, v).val[1]; }
-
-SWIZZLE(glmm_2020) { return vuzpq_f32(v, v).val[0]; }
-SWIZZLE(glmm_3131) { return vuzpq_f32(v, v).val[1]; }
-
-SWIZZLE(glmm_0321) { return vextq_f32(v, v, 1);     }
-SWIZZLE(glmm_1032) { return vextq_f32(v, v, 2);     }
-SWIZZLE(glmm_2103) { return vextq_f32(v, v, 3);     }
-
-#undef SWIZZLE
-
 #define glmm_xor(a, b)                                                        \
   vreinterpretq_f32_s32(veorq_s32(vreinterpretq_s32_f32(a),                   \
                                   vreinterpretq_s32_f32(b)))

From 701e015bfd54ffd2392b3b9e53e646d0e2fc1a61 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sun, 25 Apr 2021 02:36:06 +0300
Subject: [PATCH 13/16] avoid loading vec3 by glmm_load()

* use glmm_set1() for each for now
---
 include/cglm/affine.h | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/include/cglm/affine.h b/include/cglm/affine.h
index 504e13b..d0e5bc9 100644
--- a/include/cglm/affine.h
+++ b/include/cglm/affine.h
@@ -51,18 +51,17 @@ CGLM_INLINE
 void
 glm_translate(mat4 m, vec3 v) {
 #if defined(CGLM_SIMD)
-  glmm_128 v0, m0, m1, m2, m3;
+  glmm_128 m0, m1, m2, m3;
 
   m0 = glmm_load(m[0]);
   m1 = glmm_load(m[1]);
   m2 = glmm_load(m[2]);
   m3 = glmm_load(m[3]);
-  v0 = glmm_load(v);
 
   glmm_store(m[3],
-             glmm_fmadd(m0, glmm_splat_x(v0),
-                        glmm_fmadd(m1, glmm_splat_y(v0),
-                                   glmm_fmadd(m2, glmm_splat_z(v0), m3))));
+             glmm_fmadd(m0, glmm_set1(v[0]),
+                        glmm_fmadd(m1, glmm_set1(v[1]),
+                                   glmm_fmadd(m2, glmm_set1(v[2]), m3))));
 #else
   glm_vec4_muladds(m[0], v[0], m[3]);
   glm_vec4_muladds(m[1], v[1], m[3]);

From 0ab50f720809e91c1d4f1edc5849ed0fc920d871 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sun, 25 Apr 2021 03:41:39 +0300
Subject: [PATCH 14/16] arm, neon: update mat4_mul to use FMA

---
 include/cglm/simd/arm.h       |  2 +-
 include/cglm/simd/neon/mat4.h | 57 +++++++++++++++++------------------
 2 files changed, 28 insertions(+), 31 deletions(-)

diff --git a/include/cglm/simd/arm.h b/include/cglm/simd/arm.h
index a7009e5..e943ce3 100644
--- a/include/cglm/simd/arm.h
+++ b/include/cglm/simd/arm.h
@@ -123,7 +123,7 @@ static inline
 float32x4_t
 glmm_fmadd(float32x4_t a, float32x4_t b, float32x4_t c) {
 #if CGLM_ARM64
-  return vfmaq_f32(c, a, b);
+  return vfmaq_f32(c, a, b); /* why vfmaq_f32 is slower than vmlaq_f32 ??? */
 #else
   return vmlaq_f32(c, a, b);
 #endif
diff --git a/include/cglm/simd/neon/mat4.h b/include/cglm/simd/neon/mat4.h
index 2fecebe..f1b3b60 100644
--- a/include/cglm/simd/neon/mat4.h
+++ b/include/cglm/simd/neon/mat4.h
@@ -41,42 +41,39 @@ glm_mat4_transp_neon(mat4 m, mat4 dest) {
 CGLM_INLINE
 void
 glm_mat4_mul_neon(mat4 m1, mat4 m2, mat4 dest) {
-  /* D = R * L (Column-Major) */
-  float32x4_t l0, l1, l2, l3, r, d0, d1, d2, d3;
+   /* D = R * L (Column-Major) */
 
-  l0 = vld1q_f32(m2[0]);
-  l1 = vld1q_f32(m2[1]);
-  l2 = vld1q_f32(m2[2]);
-  l3 = vld1q_f32(m2[3]);
+  glmm_128 l0, l1, l2, l3, r0, r1, r2, r3, v0, v1, v2, v3;
 
-  r  = vld1q_f32(m1[0]);
-  d0 = vmulq_lane_f32(r, vget_low_f32(l0), 0);
-  d1 = vmulq_lane_f32(r, vget_low_f32(l1), 0);
-  d2 = vmulq_lane_f32(r, vget_low_f32(l2), 0);
-  d3 = vmulq_lane_f32(r, vget_low_f32(l3), 0);
+  l0 = glmm_load(m1[0]);    r0 = glmm_load(m2[0]);
+  l1 = glmm_load(m1[1]);    r1 = glmm_load(m2[1]);
+  l2 = glmm_load(m1[2]);    r2 = glmm_load(m2[2]);
+  l3 = glmm_load(m1[3]);    r3 = glmm_load(m2[3]);
 
-  r  = vld1q_f32(m1[1]);
-  d0 = vmlaq_lane_f32(d0, r, vget_low_f32(l0), 1);
-  d1 = vmlaq_lane_f32(d1, r, vget_low_f32(l1), 1);
-  d2 = vmlaq_lane_f32(d2, r, vget_low_f32(l2), 1);
-  d3 = vmlaq_lane_f32(d3, r, vget_low_f32(l3), 1);
+  v0 = vmulq_f32(glmm_splat_x(r0), l0);
+  v1 = vmulq_f32(glmm_splat_x(r1), l0);
+  v2 = vmulq_f32(glmm_splat_x(r2), l0);
+  v3 = vmulq_f32(glmm_splat_x(r3), l0);
 
-  r  = vld1q_f32(m1[2]);
-  d0 = vmlaq_lane_f32(d0, r, vget_high_f32(l0), 0);
-  d1 = vmlaq_lane_f32(d1, r, vget_high_f32(l1), 0);
-  d2 = vmlaq_lane_f32(d2, r, vget_high_f32(l2), 0);
-  d3 = vmlaq_lane_f32(d3, r, vget_high_f32(l3), 0);
+  v0 = glmm_fmadd(glmm_splat_y(r0), l1, v0);
+  v1 = glmm_fmadd(glmm_splat_y(r1), l1, v1);
+  v2 = glmm_fmadd(glmm_splat_y(r2), l1, v2);
+  v3 = glmm_fmadd(glmm_splat_y(r3), l1, v3);
 
-  r  = vld1q_f32(m1[3]);
-  d0 = vmlaq_lane_f32(d0, r, vget_high_f32(l0), 1);
-  d1 = vmlaq_lane_f32(d1, r, vget_high_f32(l1), 1);
-  d2 = vmlaq_lane_f32(d2, r, vget_high_f32(l2), 1);
-  d3 = vmlaq_lane_f32(d3, r, vget_high_f32(l3), 1);
+  v0 = glmm_fmadd(glmm_splat_z(r0), l2, v0);
+  v1 = glmm_fmadd(glmm_splat_z(r1), l2, v1);
+  v2 = glmm_fmadd(glmm_splat_z(r2), l2, v2);
+  v3 = glmm_fmadd(glmm_splat_z(r3), l2, v3);
 
-  vst1q_f32(dest[0], d0);
-  vst1q_f32(dest[1], d1);
-  vst1q_f32(dest[2], d2);
-  vst1q_f32(dest[3], d3);
+  v0 = glmm_fmadd(glmm_splat_w(r0), l3, v0);
+  v1 = glmm_fmadd(glmm_splat_w(r1), l3, v1);
+  v2 = glmm_fmadd(glmm_splat_w(r2), l3, v2);
+  v3 = glmm_fmadd(glmm_splat_w(r3), l3, v3);
+
+  glmm_store(dest[0], v0);
+  glmm_store(dest[1], v1);
+  glmm_store(dest[2], v2);
+  glmm_store(dest[3], v3);
 }
 
 CGLM_INLINE

From 2903813765c7f08c8545e4461e334b5d752811e9 Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sun, 25 Apr 2021 03:45:00 +0300
Subject: [PATCH 15/16] arm, neon: neon/fma support for glm_mul()

---
 include/cglm/affine-mat.h       |  6 ++++
 include/cglm/simd/neon/affine.h | 51 +++++++++++++++++++++++++++++++++
 include/cglm/simd/neon/mat4.h   |  2 +-
 3 files changed, 58 insertions(+), 1 deletion(-)
 create mode 100644 include/cglm/simd/neon/affine.h

diff --git a/include/cglm/affine-mat.h b/include/cglm/affine-mat.h
index f58749b..8383283 100644
--- a/include/cglm/affine-mat.h
+++ b/include/cglm/affine-mat.h
@@ -26,6 +26,10 @@
 #  include "simd/avx/affine.h"
 #endif
 
+#ifdef CGLM_NEON_FP
+#  include "simd/neon/affine.h"
+#endif
+
 /*!
  * @brief this is similar to glm_mat4_mul but specialized to affine transform
  *
@@ -49,6 +53,8 @@ glm_mul(mat4 m1, mat4 m2, mat4 dest) {
   glm_mul_avx(m1, m2, dest);
 #elif defined( __SSE__ ) || defined( __SSE2__ )
   glm_mul_sse2(m1, m2, dest);
+#elif defined(CGLM_NEON_FP)
+  glm_mul_neon(m1, m2, dest);
 #else
   float a00 = m1[0][0], a01 = m1[0][1], a02 = m1[0][2], a03 = m1[0][3],
         a10 = m1[1][0], a11 = m1[1][1], a12 = m1[1][2], a13 = m1[1][3],
diff --git a/include/cglm/simd/neon/affine.h b/include/cglm/simd/neon/affine.h
new file mode 100644
index 0000000..90c9ab4
--- /dev/null
+++ b/include/cglm/simd/neon/affine.h
@@ -0,0 +1,51 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglm_affine_neon_h
+#define cglm_affine_neon_h
+#if defined(__ARM_NEON_FP)
+
+#include "../../common.h"
+#include "../intrin.h"
+
+CGLM_INLINE
+void
+glm_mul_neon(mat4 m1, mat4 m2, mat4 dest) {
+  /* D = R * L (Column-Major) */
+
+  glmm_128 l0, l1, l2, l3, r0, r1, r2, r3, v0, v1, v2, v3;
+
+  l0 = glmm_load(m1[0]);    r0 = glmm_load(m2[0]);
+  l1 = glmm_load(m1[1]);    r1 = glmm_load(m2[1]);
+  l2 = glmm_load(m1[2]);    r2 = glmm_load(m2[2]);
+  l3 = glmm_load(m1[3]);    r3 = glmm_load(m2[3]);
+
+  v0 = vmulq_f32(glmm_splat_x(r0), l0);
+  v1 = vmulq_f32(glmm_splat_x(r1), l0);
+  v2 = vmulq_f32(glmm_splat_x(r2), l0);
+  v3 = vmulq_f32(glmm_splat_x(r3), l0);
+
+  v0 = glmm_fmadd(glmm_splat_y(r0), l1, v0);
+  v1 = glmm_fmadd(glmm_splat_y(r1), l1, v1);
+  v2 = glmm_fmadd(glmm_splat_y(r2), l1, v2);
+  v3 = glmm_fmadd(glmm_splat_y(r3), l1, v3);
+
+  v0 = glmm_fmadd(glmm_splat_z(r0), l2, v0);
+  v1 = glmm_fmadd(glmm_splat_z(r1), l2, v1);
+  v2 = glmm_fmadd(glmm_splat_z(r2), l2, v2);
+  v3 = glmm_fmadd(glmm_splat_z(r3), l2, v3);
+
+  v3 = glmm_fmadd(glmm_splat_w(r3), l3, v3);
+
+  glmm_store(dest[0], v0);
+  glmm_store(dest[1], v1);
+  glmm_store(dest[2], v2);
+  glmm_store(dest[3], v3);
+}
+
+#endif
+#endif /* cglm_affine_neon_h */
diff --git a/include/cglm/simd/neon/mat4.h b/include/cglm/simd/neon/mat4.h
index f1b3b60..bb00a41 100644
--- a/include/cglm/simd/neon/mat4.h
+++ b/include/cglm/simd/neon/mat4.h
@@ -41,7 +41,7 @@ glm_mat4_transp_neon(mat4 m, mat4 dest) {
 CGLM_INLINE
 void
 glm_mat4_mul_neon(mat4 m1, mat4 m2, mat4 dest) {
-   /* D = R * L (Column-Major) */
+  /* D = R * L (Column-Major) */
 
   glmm_128 l0, l1, l2, l3, r0, r1, r2, r3, v0, v1, v2, v3;
 

From 155eb109a8bb21e82992b628e02304448b387a7f Mon Sep 17 00:00:00 2001
From: Recep Aslantas <info@recp.me>
Date: Sun, 25 Apr 2021 03:49:35 +0300
Subject: [PATCH 16/16] arm, neon: neon/fma support for glm_mul_rot()

---
 include/cglm/affine-mat.h       |  2 ++
 include/cglm/simd/neon/affine.h | 29 +++++++++++++++++++++++++++++
 2 files changed, 31 insertions(+)

diff --git a/include/cglm/affine-mat.h b/include/cglm/affine-mat.h
index 8383283..c555eae 100644
--- a/include/cglm/affine-mat.h
+++ b/include/cglm/affine-mat.h
@@ -109,6 +109,8 @@ void
 glm_mul_rot(mat4 m1, mat4 m2, mat4 dest) {
 #if defined( __SSE__ ) || defined( __SSE2__ )
   glm_mul_rot_sse2(m1, m2, dest);
+#elif defined(CGLM_NEON_FP)
+  glm_mul_rot_neon(m1, m2, dest);
 #else
   float a00 = m1[0][0], a01 = m1[0][1], a02 = m1[0][2], a03 = m1[0][3],
         a10 = m1[1][0], a11 = m1[1][1], a12 = m1[1][2], a13 = m1[1][3],
diff --git a/include/cglm/simd/neon/affine.h b/include/cglm/simd/neon/affine.h
index 90c9ab4..b511083 100644
--- a/include/cglm/simd/neon/affine.h
+++ b/include/cglm/simd/neon/affine.h
@@ -47,5 +47,34 @@ glm_mul_neon(mat4 m1, mat4 m2, mat4 dest) {
   glmm_store(dest[3], v3);
 }
 
+CGLM_INLINE
+void
+glm_mul_rot_neon(mat4 m1, mat4 m2, mat4 dest) {
+  /* D = R * L (Column-Major) */
+
+  glmm_128 l0, l1, l2, r0, r1, r2, v0, v1, v2;
+
+  l0 = glmm_load(m1[0]);    r0 = glmm_load(m2[0]);
+  l1 = glmm_load(m1[1]);    r1 = glmm_load(m2[1]);
+  l2 = glmm_load(m1[2]);    r2 = glmm_load(m2[2]);
+
+  v0 = vmulq_f32(glmm_splat_x(r0), l0);
+  v1 = vmulq_f32(glmm_splat_x(r1), l0);
+  v2 = vmulq_f32(glmm_splat_x(r2), l0);
+
+  v0 = glmm_fmadd(glmm_splat_y(r0), l1, v0);
+  v1 = glmm_fmadd(glmm_splat_y(r1), l1, v1);
+  v2 = glmm_fmadd(glmm_splat_y(r2), l1, v2);
+
+  v0 = glmm_fmadd(glmm_splat_z(r0), l2, v0);
+  v1 = glmm_fmadd(glmm_splat_z(r1), l2, v1);
+  v2 = glmm_fmadd(glmm_splat_z(r2), l2, v2);
+
+  glmm_store(dest[0], v0);
+  glmm_store(dest[1], v1);
+  glmm_store(dest[2], v2);
+  glmm_store(dest[3], glmm_load(m1[3]));
+}
+
 #endif
 #endif /* cglm_affine_neon_h */