diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,5 +1,9 @@
 # Changelog for geomancy
 
+## 0.2.4.2
+
+* Support ARM/aarch64 SIMD.
+
 ## 0.2.4.1
 
 * Support simple-affine-space-0.2
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,14 +1,18 @@
 # Geomancy
 
-Linear is nice, but slow. Those are naughty, but a bit faster.
+> Linear is nice, but slow. Those are naughty, but a bit faster.
 
-* All data types are monomorphic, unpacked and specialized to `Float`.
-* `Mat4` is a `ByteArray#`.
-* `Mat4` multiplication with SSE.
-* Matrix construction states their element order.
-* Transforms don't require transposition for GLSL.
+* All data types are monomorphic, unpacked and specialized.
+* `Mat4` and `Vec4` are `ByteArray#`.
+* `Mat4`x`Mat4` and `Mat4`x`Vec4` is done with SIMD.
 
-### The Numbers
+## Matrix layout
+
+CPU-side matrices compose in MVP order, optimized for `mconcat (local1 : local2 : ... : root)` operation.
+
+GPU-side, in GLSL, it is `PVM * v`.
+
+## The Numbers
 
 Storing a list of 1000 transformations (e.g. rendering instance data):
 
diff --git a/cbits/mat4.c b/cbits/mat4.c
--- a/cbits/mat4.c
+++ b/cbits/mat4.c
@@ -1,6 +1,8 @@
+#if defined(__i386__) || defined(__x86_64__)
+
 #include <xmmintrin.h>
 
-void Mat4xMat4_SSE(float *A, float *B, float *O) {
+void Mat4xMat4_SIMD(float *A, float *B, float *O) {
   __m128 row0 = _mm_load_ps(&B[0]);
   __m128 row1 = _mm_load_ps(&B[4]);
   __m128 row2 = _mm_load_ps(&B[8]);
@@ -21,7 +23,7 @@
   }
 }
 
-void Mat4xVec4_SSE(float *M, float *V, float *O) {
+void Mat4xVec4_SIMD(float *M, float *V, float *O) {
   __m128 row0 = _mm_load_ps(&M[0]);
   __m128 row1 = _mm_load_ps(&M[4]);
   __m128 row2 = _mm_load_ps(&M[8]);
@@ -41,3 +43,99 @@
 
   _mm_storeu_ps(O, res);
 }
+
+#elif defined(__arm__) || defined(__aarch64__)
+
+#include <arm_neon.h>
+// Assumes 32bit floats
+
+void Mat4xMat4_SIMD(float *A, float *B, float *O) {
+
+    // Based on the column major ARM example:
+    // https://developer.arm.com/documentation/102467/0100/Example---matrix-multiplication
+    //
+    // NEON intrinsics mapping:
+    // https://arm-software.github.io/acle/neon_intrinsics/advsimd.html
+
+    // vfmaq_laneq_f32 -> FMLA
+    // Floating-point fused Multiply-Add to accumulator (by element). This
+    // instruction multiplies the vector elements in the first source SIMD and
+    // FP register by the specified value in the second source SIMD and FP
+    // register, and accumulates the results in the vector elements of the
+    // destination SIMD and FP register. All the values in this instruction are
+    // floating-point values.
+
+    // Contrary to the SSE function above we unroll the loop manually
+
+    // these are the rows of A
+    float32x4_t A0 = vld1q_f32(A);
+    float32x4_t A1 = vld1q_f32(&A[4]);
+    float32x4_t A2 = vld1q_f32(&A[8]);
+    float32x4_t A3 = vld1q_f32(&A[12]);
+
+    // these are the rows B
+    float32x4_t B0 = vld1q_f32(B);
+    float32x4_t B1 = vld1q_f32(&B[4]);
+    float32x4_t B2 = vld1q_f32(&B[8]);
+    float32x4_t B3 = vld1q_f32(&B[12]);
+
+    // Zero accumulators for output rows
+    float32x4_t C0 = vmovq_n_f32(0);
+    float32x4_t C1 = vmovq_n_f32(0);
+    float32x4_t C2 = vmovq_n_f32(0);
+    float32x4_t C3 = vmovq_n_f32(0);
+
+    // Multiply accumulate in 4x1 blocks to output row
+    C0 = vfmaq_laneq_f32(C0, B0, A0, 0);
+    C0 = vfmaq_laneq_f32(C0, B1, A0, 1);
+    C0 = vfmaq_laneq_f32(C0, B2, A0, 2);
+    C0 = vfmaq_laneq_f32(C0, B3, A0, 3);
+    vst1q_f32(O, C0);
+
+    C1 = vfmaq_laneq_f32(C1, B0, A1, 0);
+    C1 = vfmaq_laneq_f32(C1, B1, A1, 1);
+    C1 = vfmaq_laneq_f32(C1, B2, A1, 2);
+    C1 = vfmaq_laneq_f32(C1, B3, A1, 3);
+    vst1q_f32(&O[4], C1);
+
+    C2 = vfmaq_laneq_f32(C2, B0, A2, 0);
+    C2 = vfmaq_laneq_f32(C2, B1, A2, 1);
+    C2 = vfmaq_laneq_f32(C2, B2, A2, 2);
+    C2 = vfmaq_laneq_f32(C2, B3, A2, 3);
+    vst1q_f32(&O[8], C2);
+
+    C3 = vfmaq_laneq_f32(C3, B0, A3, 0);
+    C3 = vfmaq_laneq_f32(C3, B1, A3, 1);
+    C3 = vfmaq_laneq_f32(C3, B2, A3, 2);
+    C3 = vfmaq_laneq_f32(C3, B3, A3, 3);
+    vst1q_f32(&O[12], C3);
+
+}
+
+void Mat4xVec4_SIMD(float *M, float *V, float *O) {
+
+    // Rows
+    float32x4_t M0 = vld1q_f32(M);
+    float32x4_t M1 = vld1q_f32(&M[4]);
+    float32x4_t M2 = vld1q_f32(&M[8]);
+    float32x4_t M3 = vld1q_f32(&M[12]);
+
+    // Col
+    float32x4_t VC = vld1q_f32(V);
+
+    // Local output vector
+    float32x4_t L0 = vmovq_n_f32(0);
+
+    L0 = vfmaq_laneq_f32(L0, M0, VC, 0);
+    L0 = vfmaq_laneq_f32(L0, M1, VC, 1);
+    L0 = vfmaq_laneq_f32(L0, M2, VC, 2);
+    L0 = vfmaq_laneq_f32(L0, M3, VC, 3);
+
+    vst1q_f32(O, L0);
+}
+
+#else
+
+#error "Architecture not supported! Please submit an issue."
+
+#endif
diff --git a/geomancy.cabal b/geomancy.cabal
--- a/geomancy.cabal
+++ b/geomancy.cabal
@@ -1,13 +1,13 @@
 cabal-version: 1.12
 
--- This file has been generated from package.yaml by hpack version 0.35.0.
+-- This file has been generated from package.yaml by hpack version 0.35.1.
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: c412ee063c335404b0a3d322adb6e19f362c3bfdedf86b7720b8fe6db4a7a741
+-- hash: 13d331028f6ae562f3d1e96d1909faf808cab9b792c467ef78c9d3fa8285c72e
 
 name:           geomancy
-version:        0.2.4.1
+version:        0.2.4.2
 synopsis:       Geometry and matrix manipulation
 description:    Sometimes it is unavoidable you have to do stuff on CPU.
                 Let's at least do it faster.
diff --git a/src/Geomancy/Mat4.hs b/src/Geomancy/Mat4.hs
--- a/src/Geomancy/Mat4.hs
+++ b/src/Geomancy/Mat4.hs
@@ -402,13 +402,13 @@
 zipWith :: (Float -> Float -> c) -> Mat4 -> Mat4 -> [c]
 zipWith f a b = List.zipWith f (toList a) (toList b)
 
-foreign import ccall unsafe "Mat4xMat4_SSE" m4m4sse :: Addr# -> Addr# -> Addr# -> IO ()
+foreign import ccall unsafe "Mat4xMat4_SIMD" m4m4simd :: Addr# -> Addr# -> Addr# -> IO ()
 
 {-# INLINE matrixProduct #-}
 matrixProduct :: Mat4 -> Mat4 -> Mat4
 matrixProduct (Mat4 l) (Mat4 r) = unsafePerformIO do
   result@(Mat4 m) <- unsafeNewMat4
-  m4m4sse
+  m4m4simd
     (byteArrayContents# l)
     (byteArrayContents# r)
     (byteArrayContents# m)
@@ -438,13 +438,13 @@
         (m02 * x) (m12 * x) (m22 * x) (m32 * x)
         (m03 * x) (m13 * x) (m23 * x) (m33 * x)
 
-foreign import ccall unsafe "Mat4xVec4_SSE" m4v4sse :: Addr# -> Addr# -> Addr# -> IO ()
+foreign import ccall unsafe "Mat4xVec4_SIMD" m4v4simd :: Addr# -> Addr# -> Addr# -> IO ()
 
 -- | Matrix - column vector multiplication
 (!*) :: Coercible a Mat4 => a -> Vec4 -> Vec4
 (!*) (coerce -> Mat4 m) (Vec4 v) = unsafePerformIO do
   result@(Vec4 o) <- unsafeNewVec4
-  m4v4sse
+  m4v4simd
     (byteArrayContents# m)
     (byteArrayContents# v)
     (byteArrayContents# o)
diff --git a/src/Geomancy/Tree.hs b/src/Geomancy/Tree.hs
--- a/src/Geomancy/Tree.hs
+++ b/src/Geomancy/Tree.hs
@@ -39,7 +39,6 @@
   -> Tree (acc, a)
 applyWith f = mapAccum next
   where
-    -- nextAcc = f ann acc
     next acc (ann, item) =
       let
         acc' = f ann acc