diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,5 +1,10 @@
 # Changelog for geomancy
 
+## 0.2.3.0
+
+- Vec4 moved to ByteArray#.
+- Mat4 !* Vec4 is now SIMD too.
+
 ## 0.2.2.4
 
 + Add Transform trees.
diff --git a/cbits/mat4.c b/cbits/mat4.c
--- a/cbits/mat4.c
+++ b/cbits/mat4.c
@@ -1,22 +1,43 @@
 #include <xmmintrin.h>
 
-void M4x4_SSE(float *A, float *B, float *C) {
-    __m128 row1 = _mm_load_ps(&B[0]);
-    __m128 row2 = _mm_load_ps(&B[4]);
-    __m128 row3 = _mm_load_ps(&B[8]);
-    __m128 row4 = _mm_load_ps(&B[12]);
-    for(int i=0; i<4; i++) {
-        __m128 brod1 = _mm_set1_ps(A[4*i + 0]);
-        __m128 brod2 = _mm_set1_ps(A[4*i + 1]);
-        __m128 brod3 = _mm_set1_ps(A[4*i + 2]);
-        __m128 brod4 = _mm_set1_ps(A[4*i + 3]);
-        __m128 row = _mm_add_ps(
-                    _mm_add_ps(
-                        _mm_mul_ps(brod1, row1),
-                        _mm_mul_ps(brod2, row2)),
-                    _mm_add_ps(
-                        _mm_mul_ps(brod3, row3),
-                        _mm_mul_ps(brod4, row4)));
-        _mm_store_ps(&C[4*i], row);
-    }
+void Mat4xMat4_SSE(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]);
+  __m128 row3 = _mm_load_ps(&B[12]);
+  for(int i=0; i<4; i++) {
+    __m128 brod0 = _mm_set1_ps(A[4*i + 0]);
+    __m128 brod1 = _mm_set1_ps(A[4*i + 1]);
+    __m128 brod2 = _mm_set1_ps(A[4*i + 2]);
+    __m128 brod3 = _mm_set1_ps(A[4*i + 3]);
+    __m128 row = _mm_add_ps(
+      _mm_add_ps(
+        _mm_mul_ps(brod0, row0),
+        _mm_mul_ps(brod1, row1)),
+      _mm_add_ps(
+        _mm_mul_ps(brod2, row2),
+        _mm_mul_ps(brod3, row3)));
+    _mm_store_ps(&O[4*i], row);
+  }
+}
+
+void Mat4xVec4_SSE(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]);
+  __m128 row3 = _mm_load_ps(&M[12]);
+
+  __m128 vec0 = _mm_load_ps(V);
+
+  __m128 vec00 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(0, 0, 0, 0));
+  __m128 vec01 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(1, 1, 1, 1));
+  __m128 vec02 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(2, 2, 2, 2));
+  __m128 vec03 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(3, 3, 3, 3));
+
+  __m128 res = _mm_mul_ps(row0, vec00);
+  res = _mm_add_ps(res, _mm_mul_ps(row1, vec01));
+  res = _mm_add_ps(res, _mm_mul_ps(row2, vec02));
+  res = _mm_add_ps(res, _mm_mul_ps(row3, vec03));
+
+  _mm_storeu_ps(O, res);
 }
diff --git a/geomancy.cabal b/geomancy.cabal
--- a/geomancy.cabal
+++ b/geomancy.cabal
@@ -4,10 +4,10 @@
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: 1480b5bc5e86691e54a70d648214d4951a3685c03bc3eeb8741a8429eaa0b13c
+-- hash: aa93c14151d6ba606892efc64f6795f6f0ec8ebfc5b9f24e104d5b4676cf8bb8
 
 name:           geomancy
-version:        0.2.2.4
+version:        0.2.3.0
 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
@@ -5,6 +5,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE UnboxedTuples #-}
 {-# LANGUAGE UnliftedFFITypes #-}
+{-# LANGUAGE ViewPatterns #-}
 
 -- | General matrix storage and operations.
 
@@ -33,7 +34,7 @@
   ) where
 
 import Prelude hiding (zipWith)
-import GHC.Exts hiding (toList)
+import GHC.Exts hiding (VecCount(..), toList)
 
 import Control.DeepSeq (NFData(rnf))
 import Foreign (Storable(..))
@@ -44,7 +45,7 @@
 import qualified Data.Foldable as Foldable
 import qualified Data.List as List
 
-import Geomancy.Vec4 (Vec4, vec4, withVec4)
+import Geomancy.Vec4 (Vec4(..), unsafeNewVec4)
 
 data Mat4 = Mat4 ByteArray#
 
@@ -401,42 +402,13 @@
 zipWith :: (Float -> Float -> c) -> Mat4 -> Mat4 -> [c]
 zipWith f a b = List.zipWith f (toList a) (toList b)
 
-{-# INLINE scalarMultiply #-}
-scalarMultiply :: Float -> Mat4 -> Mat4
-scalarMultiply x m =
-  withMat4 m
-    \ m00 m01 m02 m03
-      m10 m11 m12 m13
-      m20 m21 m22 m23
-      m30 m31 m32 m33 ->
-      mat4
-        (m00 * x) (m10 * x) (m20 * x) (m30 * x)
-        (m01 * x) (m11 * x) (m21 * x) (m31 * x)
-        (m02 * x) (m12 * x) (m22 * x) (m32 * x)
-        (m03 * x) (m13 * x) (m23 * x) (m33 * x)
-
--- | Matrix - column vector multiplication
-(!*) :: Coercible a Mat4 => a -> Vec4 -> Vec4
-(!*) mat vec =
-  withVec4 vec \v1 v2 v3 v4 ->
-    withColMajor mat
-      \ m11 m12 m13 m14
-        m21 m22 m23 m24
-        m31 m32 m33 m34
-        m41 m42 m43 m44 ->
-          vec4
-            (m11 * v1 + m12 * v2 + m13 * v3 + m14 * v4)
-            (m21 * v1 + m22 * v2 + m23 * v3 + m24 * v4)
-            (m31 * v1 + m32 * v2 + m33 * v3 + m34 * v4)
-            (m41 * v1 + m42 * v2 + m43 * v3 + m44 * v4)
-
-foreign import ccall unsafe "M4x4_SSE" mm4sse :: Addr# -> Addr# -> Addr# -> IO ()
+foreign import ccall unsafe "Mat4xMat4_SSE" m4m4sse :: Addr# -> Addr# -> Addr# -> IO ()
 
 {-# INLINE matrixProduct #-}
 matrixProduct :: Mat4 -> Mat4 -> Mat4
 matrixProduct (Mat4 l) (Mat4 r) = unsafePerformIO do
   result@(Mat4 m) <- unsafeNewMat4
-  mm4sse
+  m4m4sse
     (byteArrayContents# l)
     (byteArrayContents# r)
     (byteArrayContents# m)
@@ -451,6 +423,44 @@
       !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_
     in
       (# world, Mat4 arr #)
+
+{-# INLINE scalarMultiply #-}
+scalarMultiply :: Float -> Mat4 -> Mat4
+scalarMultiply x m =
+  withMat4 m
+    \ m00 m01 m02 m03
+      m10 m11 m12 m13
+      m20 m21 m22 m23
+      m30 m31 m32 m33 ->
+      mat4
+        (m00 * x) (m10 * x) (m20 * x) (m30 * x)
+        (m01 * x) (m11 * x) (m21 * x) (m31 * x)
+        (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 ()
+
+-- | Matrix - column vector multiplication
+(!*) :: Coercible a Mat4 => a -> Vec4 -> Vec4
+(!*) (coerce -> Mat4 m) (Vec4 v) = unsafePerformIO do
+  result@(Vec4 o) <- unsafeNewVec4
+  m4v4sse
+    (byteArrayContents# m)
+    (byteArrayContents# v)
+    (byteArrayContents# o)
+  pure result
+
+  -- withVec4 vec \v1 v2 v3 v4 ->
+  --   withColMajor mat
+  --     \ m11 m12 m13 m14
+  --       m21 m22 m23 m24
+  --       m31 m32 m33 m34
+  --       m41 m42 m43 m44 ->
+  --         vec4
+  --           (m11 * v1 + m12 * v2 + m13 * v3 + m14 * v4)
+  --           (m21 * v1 + m22 * v2 + m23 * v3 + m24 * v4)
+  --           (m31 * v1 + m32 * v2 + m33 * v3 + m34 * v4)
+  --           (m41 * v1 + m42 * v2 + m43 * v3 + m44 * v4)
 
 instance NFData Mat4 where
   rnf Mat4{} = ()
diff --git a/src/Geomancy/Transform.hs b/src/Geomancy/Transform.hs
--- a/src/Geomancy/Transform.hs
+++ b/src/Geomancy/Transform.hs
@@ -29,10 +29,13 @@
 
 import Foreign (Storable(..))
 
-import Geomancy.Mat4 (Mat4, colMajor, withColMajor, inverse)
+import Geomancy.Mat4 (Mat4, colMajor, inverse)
 import Geomancy.Quaternion (Quaternion, withQuaternion)
 import Geomancy.Vec3 (Vec3, vec3, withVec3)
+import Geomancy.Vec4 (fromVec3, withVec4)
 
+import qualified Geomancy.Mat4 as Mat4
+
 newtype Transform = Transform { unTransform :: Mat4 }
   deriving newtype (Show, Semigroup, Monoid, Storable)
 
@@ -43,19 +46,10 @@
 -- | Matrix - column vector multiplication with perspective division
 (!.) :: Transform -> Vec3 -> Vec3
 (!.) mat vec =
-  withVec3 vec \v1 v2 v3 ->
-    withColMajor mat
-      \ m11 m12 m13 m14
-        m21 m22 m23 m24
-        m31 m32 m33 m34
-        m41 m42 m43 m44 ->
-          let
-            px = m11 * v1 + m12 * v2 + m13 * v3 + m14
-            py = m21 * v1 + m22 * v2 + m23 * v3 + m24
-            pz = m31 * v1 + m32 * v2 + m33 * v3 + m34
-            p  = m41 * v1 + m42 * v2 + m43 * v3 + m44
-          in
-            vec3 (px / p) (py / p) (pz / p)
+  withVec4 res \x y z w ->
+    vec3 (x / w) (y / w) (z / w)
+  where
+    res = mat Mat4.!* fromVec3 vec 1.0
 
 -- ** Translation
 
diff --git a/src/Geomancy/Vec4.hs b/src/Geomancy/Vec4.hs
--- a/src/Geomancy/Vec4.hs
+++ b/src/Geomancy/Vec4.hs
@@ -1,12 +1,16 @@
+{-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE BlockArguments #-}
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MagicHash #-}
 {-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE UnboxedTuples #-}
+{-# LANGUAGE UnliftedFFITypes #-}
 {-# LANGUAGE ViewPatterns #-}
 
 -- | Specialized and inlined @V4 Float@.
 
 module Geomancy.Vec4
-  ( Vec4
+  ( Vec4(..)
   , vec4
   , withVec4
   , pattern WithVec4
@@ -21,33 +25,74 @@
 
   , dot
   , normalize
+
+  , unsafeNewVec4
   ) where
 
+import GHC.Exts hiding (VecCount(..), toList)
+
 import Control.DeepSeq (NFData(rnf))
-import Data.Coerce (Coercible, coerce)
-import Foreign (Storable(..), castPtr)
+import Foreign (Storable(..))
+import GHC.IO (IO(..))
+-- import System.IO.Unsafe (unsafePerformIO)
+import Text.Printf (printf)
 
 import Geomancy.Vec2 (Vec2, withVec2)
 import Geomancy.Vec3 (Vec3, withVec3)
 
-data Vec4 = Vec4
-  {-# UNPACK #-} !Float
-  {-# UNPACK #-} !Float
-  {-# UNPACK #-} !Float
-  {-# UNPACK #-} !Float
-  deriving (Eq, Ord, Show)
+data Vec4 = Vec4 ByteArray#
 
 {-# INLINE vec4 #-}
 vec4 :: Float -> Float -> Float -> Float -> Vec4
-vec4 = Vec4
+vec4 (F# v0) (F# v1) (F# v2) (F# v3) =
+  runRW# \world ->
+    let
+      !(# world_, arr #) = newAlignedPinnedByteArray# 16# 16# world
 
+      world0 = writeFloatArray# arr 0x0# v0 world_
+      world1 = writeFloatArray# arr 0x1# v1 world0
+      world2 = writeFloatArray# arr 0x2# v2 world1
+      world3 = writeFloatArray# arr 0x3# v3 world2
+      !(# _world', arr' #) = unsafeFreezeByteArray# arr world3
+    in
+      Vec4 arr'
+
 {-# INLINE withVec4 #-}
 withVec4
   :: Vec4
   -> (Float -> Float -> Float -> Float -> r)
   -> r
-withVec4 (Vec4 a b c d) f = f a b c d
+withVec4 (Vec4 arr) f =
+  f
+    (F# (indexFloatArray# arr 0x0#))
+    (F# (indexFloatArray# arr 0x1#))
+    (F# (indexFloatArray# arr 0x2#))
+    (F# (indexFloatArray# arr 0x3#))
 
+{-# INLINE compareVec4 #-}
+compareVec4 :: Vec4 -> Vec4 -> Ordering
+compareVec4 (Vec4 src1) (Vec4 src2) =
+  compare (I# (compareByteArrays# src1 0# src2 0# 16#)) 0
+
+instance Eq Vec4 where
+  (==) a b =
+    case compareVec4 a b of
+      EQ -> True
+      _  -> False
+
+  (/=) a b =
+    case compareVec4 a b of
+      EQ -> False
+      _  -> True
+
+instance Ord Vec4 where
+  compare = compareVec4
+
+instance Show Vec4 where
+  show v =
+    withVec4 v $
+      printf "Vec4 %.4f %.4f %.4f %.4f"
+
 pattern WithVec4 :: Float -> Float -> Float -> Float -> Vec4
 pattern WithVec4 a b c d <- ((`withVec4` (,,,)) -> (a, b, c, d))
 {-# COMPLETE WithVec4 #-}
@@ -80,53 +125,64 @@
 
 instance Num Vec4 where
   {-# INLINE (+) #-}
-  Vec4 l1 l2 l3 l4 + Vec4 r1 r2 r3 r4 =
-    Vec4
-      (l1 + r1)
-      (l2 + r2)
-      (l3 + r3)
-      (l4 + r4)
+  (+) l r =
+    withVec4 l \l1 l2 l3 l4 ->
+      withVec4 r \r1 r2 r3 r4 ->
+        vec4
+          (l1 + r1)
+          (l2 + r2)
+          (l3 + r3)
+          (l4 + r4)
 
   {-# INLINE (-) #-}
-  Vec4 l1 l2 l3 l4 - Vec4 r1 r2 r3 r4 =
-    Vec4
-      (l1 - r1)
-      (l2 - r2)
-      (l3 - r3)
-      (l4 - r4)
+  (-) l r =
+    withVec4 l \l1 l2 l3 l4 ->
+      withVec4 r \r1 r2 r3 r4 ->
+        vec4
+          (l1 - r1)
+          (l2 - r2)
+          (l3 - r3)
+          (l4 - r4)
 
   {-# INLINE (*) #-}
-  Vec4 l1 l2 l3 l4 * Vec4 r1 r2 r3 r4 =
-    Vec4
-      (l1 * r1)
-      (l2 * r2)
-      (l3 * r3)
-      (l4 * r4)
+  (*) l r =
+    withVec4 l \l1 l2 l3 l4 ->
+      withVec4 r \r1 r2 r3 r4 ->
+        vec4
+          (l1 * r1)
+          (l2 * r2)
+          (l3 * r3)
+          (l4 * r4)
 
   {-# INLINE abs #-}
-  abs (Vec4 a b c d) =
-    Vec4 (abs a) (abs b) (abs c) (abs d)
+  abs v =
+    withVec4 v \a b c d ->
+      vec4 (abs a) (abs b) (abs c) (abs d)
 
   {-# INLINE signum #-}
-  signum (Vec4 a b c d) =
-    Vec4 (signum a) (signum b) (signum c) (signum d)
+  signum v =
+    withVec4 v \a b c d ->
+      vec4 (signum a) (signum b) (signum c) (signum d)
 
   {-# INLINE fromInteger #-}
-  fromInteger x = Vec4 x' x' x' x'
+  fromInteger x = vec4 x' x' x' x'
     where
       x' = fromInteger x
 
 instance Fractional Vec4 where
   {-# INLINE (/) #-}
-  Vec4 l1 l2 l3 l4 / Vec4 r1 r2 r3 r4 =
-    Vec4 (l1 / r1) (l2 / r2) (l3 / r3) (l4 / r4)
+  (/) l r =
+    withVec4 l \l1 l2 l3 l4 ->
+      withVec4 r \r1 r2 r3 r4 ->
+        vec4 (l1 / r1) (l2 / r2) (l3 / r3) (l4 / r4)
 
   {-# INLINE recip #-}
-  recip (Vec4 a b c d) =
-    Vec4 (recip a) (recip b) (recip c) (recip d)
+  recip v =
+    withVec4 v \a b c d ->
+      vec4 (recip a) (recip b) (recip c) (recip d)
 
   {-# INLINE fromRational #-}
-  fromRational x = Vec4 x' x' x' x'
+  fromRational x = vec4 x' x' x' x'
     where
       x' = fromRational x
 
@@ -138,53 +194,50 @@
   alignment _ = 16
 
   {-# INLINE poke #-}
-  poke ptr v4 =
-    withVec4 v4 \a b c d -> do
-      poke ptr' a
-      pokeElemOff ptr' 1 b
-      pokeElemOff ptr' 2 c
-      pokeElemOff ptr' 3 d
-    where
-      ptr' = castPtr ptr
+  poke (Ptr addr) (Vec4 arr) = IO \world ->
+    let
+      world' = copyByteArrayToAddr# arr 0# addr 16# world
+    in
+      (# world', () #)
 
   {-# INLINE peek #-}
-  peek ptr = vec4
-    <$> peek ptr'
-    <*> peekElemOff ptr' 1
-    <*> peekElemOff ptr' 2
-    <*> peekElemOff ptr' 3
-    where
-      ptr' = castPtr ptr
+  peek (Ptr addr) = IO \world ->
+    let
+      !(# world0, arr #)  = newAlignedPinnedByteArray# 16# 16# world
+      world1              = copyAddrToByteArray# addr arr 0# 16# world0
+      !(# world', arr' #) = unsafeFreezeByteArray# arr world1
+    in
+      (# world', Vec4 arr' #)
 
+-- TODO: SIMD
 {-# INLINE (^*) #-}
 (^*) :: Vec4 -> Float -> Vec4
-Vec4 a b c d ^* x =
-  Vec4
-    (a * x)
-    (b * x)
-    (c * x)
-    (d * x)
+(^*) v x =
+  withVec4 v \a b c d ->
+    vec4
+      (a * x)
+      (b * x)
+      (c * x)
+      (d * x)
 
 {-# INLINE (^/) #-}
 (^/) :: Vec4 -> Float -> Vec4
-Vec4 a b c d ^/ x =
-  Vec4
-    (a / x)
-    (b / x)
-    (c / x)
-    (d / x)
+v ^/ x = v ^* recip x
 
 {-# INLINE lerp #-}
 lerp :: Float -> Vec4 -> Vec4 -> Vec4
 lerp alpha u v = u ^* alpha + v ^* (1 - alpha)
 
+-- TODO: SIMD
 {-# INLINE dot #-}
 dot :: Vec4 -> Vec4 -> Float
-dot (Vec4 l1 l2 l3 l4) (Vec4 r1 r2 r3 r4) =
-  l1 * r1 +
-  l2 * r2 +
-  l3 * r3 +
-  l4 * r4
+dot a b =
+  withVec4 a \a1 a2 a3 a4 ->
+    withVec4 b \b1 b2 b3 b4 ->
+      a1 * b1 +
+      a2 * b2 +
+      a3 * b3 +
+      a4 * b4
 
 {-# INLINE normalize #-}
 normalize :: Vec4 -> Vec4
@@ -192,13 +245,20 @@
   if nearZero q || nearZero (1-q) then
     v
   else
-    let
-      Vec4 x y z w = v
-    in
-      Vec4 (x / l) (y / l) (z / l) (w / l)
+    v ^/ l
 
   where
     q = dot v v
     l = sqrt q
 
     nearZero a = abs a <= 1e-6
+
+{-# INLINE unsafeNewVec4 #-}
+unsafeNewVec4 :: IO Vec4
+unsafeNewVec4 =
+  IO \world ->
+    let
+      !(# world_, arr_ #) = newAlignedPinnedByteArray# 16# 16# world
+      !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_
+    in
+      (# world, Vec4 arr #)