diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,76 @@
+# Changelog for geomancy
+
+## 0.2.2.3
+
++ Add IVec and UVec 32-bit integer vectors.
+
+## 0.2.2.2
+
++ Fixed using Transform.apply with projection inverses.
+
+## 0.2.2.1
+
++ Add Vec2.
++ Add pattern synonyms to `Geomancy` re-exports.
++ Add Fractional instances.
++ Add conversions from tuples and lower-dimension vecs.
+
+## 0.2.2.0
+
++ Add pattern synonym alternatives to `withVecN`.
++ Add Mat4 converstion from `Linear.M44`.
+- Hide `toList` and `toListTrans`.
++ Add `toList2d` and publish its element-order wrappers.
+
+## 0.2.1.0
+
+Transform rewrite
+
+* Change perspective FoV to radians.
++ Add `infinitePerspective`.
++ Add `Transform.inverse`.
+- Hide `mat4`, `withMat4`.
+
+## 0.2.0.0
+
+Mat4 rewrite
+
++ Add `rowMajor`, `withRowMajor`, `toListRowMajor`.
++ Add `colMajor`, `withColMajor`, `toListColMajor`.
++ Add `Mat4.inverse`.
+* Expose `matrixProduct`.
++ Add `scalarMultiply`.
++ Add `Mat4.(!*)` to use with `Vec4`.
+* Extract transformations to `Geomancy.Transform` and use column notation.
++ Add `Transform.(!.)` and `apply` to use with `Vec3`.
+* Rename `scale` to `scale3`.
++ Add uniform `scale`.
++ Add `scaleXY` for flat meshes.
+* Rename `mkTransformation` to `dirPos`.
++ Add `rotateQ` via `dirPos` with empty translation.
++ Add `Vec3.Packed` newtype without extra padding.
+
+## 0.1.3.0
+
+* Update tests
++ Add `zipWith`
+* Rename `elementwise` to `pointwise`
+* Rename `colMajor` to `toList`
+* Rename `rowMajor` to `toListTrans`
+
+## 0.1.2.1
+
++ Add Mat4 multiplication test via `linear`
++ Add `elementwise`, `colMajor`, `rowMajor`
+
+## 0.1.2.0
+
+* Move projections and views to Vulkan namespace.
+
+## 0.1.1.2
+
+* Fix bug in quaternion rotationBetween.
+
+## 0.1.1.1
+
+* Add lookAtUp and rotationBetween.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Alexander Bondarenko (c) 2020
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Alexander Bondarenko nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,62 @@
+# Geomancy
+
+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.
+
+### The Numbers
+
+Storing a list of 1000 transformations (e.g. rendering instance data):
+
+```
+benchmarking 4x4 poke/1000/geomancy
+time                 11.76 μs   (11.66 μs .. 11.92 μs)
+                     0.999 R²   (0.998 R² .. 1.000 R²)
+mean                 11.75 μs   (11.69 μs .. 11.86 μs)
+std dev              283.4 ns   (199.0 ns .. 399.0 ns)
+variance introduced by outliers: 26% (moderately inflated)
+```
+
+If you're willing to adjust your shaders, it's only 2.4 times slower.
+
+```
+benchmarking 4x4 poke/1000/linear
+time                 28.29 μs   (28.21 μs .. 28.38 μs)
+                     1.000 R²   (1.000 R² .. 1.000 R²)
+mean                 28.40 μs   (28.34 μs .. 28.50 μs)
+std dev              267.4 ns   (145.5 ns .. 419.9 ns)
+```
+
+Keeping your shaders straight make the affair 6.1x slower.
+
+```
+benchmarking 4x4 poke/1000/linear/T
+time                 73.70 μs   (73.06 μs .. 74.49 μs)
+                     1.000 R²   (0.999 R² .. 1.000 R²)
+mean                 72.77 μs   (72.50 μs .. 73.22 μs)
+std dev              1.129 μs   (793.5 ns .. 1.580 μs)
+```
+
+Folding down a `gloss`-style scene graph is where it is all started:
+
+```
+benchmarking 4x4 multiply/1000/geomancy
+time                 20.79 μs   (20.77 μs .. 20.83 μs)
+                     1.000 R²   (1.000 R² .. 1.000 R²)
+mean                 20.80 μs   (20.78 μs .. 20.83 μs)
+std dev              76.71 ns   (60.01 ns .. 99.06 ns)
+
+benchmarking 4x4 multiply/1000/linear
+time                 173.9 μs   (173.6 μs .. 174.4 μs)
+                     1.000 R²   (1.000 R² .. 1.000 R²)
+mean                 173.5 μs   (173.2 μs .. 174.4 μs)
+std dev              1.733 μs   (727.8 ns .. 3.422 μs)
+```
+
+Add that time to the poking that'll follow.
+
+Sure, it is in the lower microseconds range, but this budget can be used elsewhere.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/bench/Main.hs b/bench/Main.hs
new file mode 100644
--- /dev/null
+++ b/bench/Main.hs
@@ -0,0 +1,75 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE TypeApplications #-}
+
+module Main where
+
+import Criterion.Main
+
+import Data.List (foldl1')
+import qualified Foreign
+
+import qualified Geomancy
+import qualified Geomancy.Mat4
+
+import qualified Linear
+import qualified Linear.Matrix
+
+main :: IO ()
+main = defaultMain
+  [ bgroup "4x4 transpose"
+      mat4transpose
+
+  , bgroup "4x4 multiply" $ flip map [2, 3, 10, 100, 1000] \size ->
+      env (clones size mempty Linear.identity) $
+        mat4multiply size
+
+  , bgroup "4x4 poke" $ flip map [10, 100, 1000, 10000] \size ->
+      env (clones size mempty Linear.identity) $
+        mat4poke size
+  ]
+
+mat4transpose :: [Benchmark]
+mat4transpose =
+  [ bench "geomancy" $ whnf Geomancy.Mat4.transpose mempty
+  , bench "linear"   $ whnf Linear.Matrix.transpose (Linear.identity :: Linear.M44 Float)
+  ]
+
+mat4multiply
+  :: Int
+  -> ([Geomancy.Mat4], [Linear.M44 Float])
+  -> Benchmark
+mat4multiply size ~(mat4s, m44s) = bgroup (show size)
+  [ bench "geomancy" $ whnf multGeomancy mat4s
+  , bench "linear"   $ whnf multLinear m44s
+  ]
+  where
+    multGeomancy = foldl1' (<>)
+    multLinear = foldl1' (Linear.!*!)
+
+clones :: Applicative f => Int -> a -> b -> f ([a], [b])
+clones n geo lin = (,) <$> genMat4s <*> genM44s
+  where
+    genMat4s = pure $! replicate n geo
+    genM44s = pure $! replicate n lin
+
+mat4poke
+  :: Int
+  -> ([Geomancy.Mat4], [Linear.M44 Float])
+  -> Benchmark
+mat4poke size ~(mat4s, m44s) = bgroup (show size)
+  [ bench "geomancy" $ nfIO pokeGeomancy --  mat4s
+  , bench "linear"   $ nfIO pokeLinear --  m44s
+  , bench "linear/T" $ nfIO pokeLinearTranspose
+  ]
+  where
+    pokeGeomancy =
+      Foreign.withArray mat4s \_ptr ->
+        pure ()
+
+    pokeLinear =
+      Foreign.withArray m44s \_ptr ->
+        pure ()
+
+    pokeLinearTranspose =
+      Foreign.withArray (map Linear.transpose m44s) \_ptr ->
+        pure ()
diff --git a/cbits/mat4.c b/cbits/mat4.c
new file mode 100644
--- /dev/null
+++ b/cbits/mat4.c
@@ -0,0 +1,22 @@
+#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);
+    }
+}
diff --git a/geomancy.cabal b/geomancy.cabal
new file mode 100644
--- /dev/null
+++ b/geomancy.cabal
@@ -0,0 +1,88 @@
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.33.0.
+--
+-- see: https://github.com/sol/hpack
+--
+-- hash: 0bbbcdde42a68176d0c54cb31c5798701541c1ce73877ad066e6aec522450641
+
+name:           geomancy
+version:        0.2.2.3
+synopsis:       Geometry and matrix manipulation
+description:    Sometimes it is unavoidable you have to do stuff on CPU.
+                Let's at least do it faster.
+category:       Graphics
+author:         Alexander Bondarenko
+maintainer:     aenor.realm@gmail.com
+copyright:      2021 Alexander Bondarenko
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+extra-source-files:
+    README.md
+    ChangeLog.md
+
+source-repository head
+  type: git
+  location: https://gitlab.com/dpwiz/geomancy
+
+library
+  exposed-modules:
+      Geomancy
+      Geomancy.IVec2
+      Geomancy.IVec3
+      Geomancy.IVec4
+      Geomancy.Mat4
+      Geomancy.Quaternion
+      Geomancy.Transform
+      Geomancy.UVec2
+      Geomancy.UVec3
+      Geomancy.UVec4
+      Geomancy.Vec2
+      Geomancy.Vec3
+      Geomancy.Vec4
+      Geomancy.Vulkan.Projection
+      Geomancy.Vulkan.View
+  other-modules:
+      Paths_geomancy
+  hs-source-dirs:
+      src
+  ghc-options: -Wall
+  c-sources:
+      cbits/mat4.c
+  build-depends:
+      base >=4.7 && <5
+    , deepseq
+  default-language: Haskell2010
+
+test-suite geomancy-test
+  type: exitcode-stdio-1.0
+  main-is: Spec.hs
+  other-modules:
+      Paths_geomancy
+  hs-source-dirs:
+      test
+  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N -with-rtsopts=-qn1 -with-rtsopts=-A128m
+  build-depends:
+      base >=4.7 && <5
+    , deepseq
+    , geomancy
+    , hedgehog
+    , linear
+  default-language: Haskell2010
+
+benchmark geomancy-bench
+  type: exitcode-stdio-1.0
+  main-is: Main.hs
+  other-modules:
+      Paths_geomancy
+  hs-source-dirs:
+      bench
+  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N -with-rtsopts=-qn1 -with-rtsopts=-A16m
+  build-depends:
+      base >=4.7 && <5
+    , criterion
+    , deepseq
+    , geomancy
+    , linear
+  default-language: Haskell2010
diff --git a/src/Geomancy.hs b/src/Geomancy.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy.hs
@@ -0,0 +1,85 @@
+{-# LANGUAGE PatternSynonyms #-}
+
+module Geomancy
+  (
+  -- * Vectors
+
+  -- ** Single-precision / float32s
+
+    Vec2
+  , vec2
+  , withVec2
+  , pattern WithVec2
+
+  , Vec3
+  , vec3
+  , withVec3
+  , pattern WithVec3
+
+  , Vec4
+  , vec4
+  , withVec4
+  , pattern WithVec4
+
+  -- ** Signed / int32s
+
+  , IVec2
+  , ivec2
+  , withIVec2
+  , pattern WithIVec2
+
+  , IVec3
+  , ivec3
+  , withIVec3
+  , pattern WithIVec3
+
+  , IVec4
+  , ivec4
+  , withIVec4
+  , pattern WithIVec4
+
+  -- ** Unsigned / word32s
+
+  , UVec2
+  , uvec2
+  , withUVec2
+  , pattern WithUVec2
+
+  , UVec3
+  , uvec3
+  , withUVec3
+  , pattern WithUVec3
+
+  , UVec4
+  , uvec4
+  , withUVec4
+  , pattern WithUVec4
+
+  -- * Matrices
+
+  , Mat4
+  , Transform(..)
+
+  -- * Other beasts
+
+  , Quaternion
+  , quaternion
+  , withQuaternion
+  ) where
+
+import Geomancy.Vec2 (Vec2, vec2, withVec2, pattern WithVec2)
+import Geomancy.Vec3 (Vec3, vec3, withVec3, pattern WithVec3)
+import Geomancy.Vec4 (Vec4, vec4, withVec4, pattern WithVec4)
+
+import Geomancy.IVec2 (IVec2, ivec2, withIVec2, pattern WithIVec2)
+import Geomancy.IVec3 (IVec3, ivec3, withIVec3, pattern WithIVec3)
+import Geomancy.IVec4 (IVec4, ivec4, withIVec4, pattern WithIVec4)
+
+import Geomancy.UVec2 (UVec2, uvec2, withUVec2, pattern WithUVec2)
+import Geomancy.UVec3 (UVec3, uvec3, withUVec3, pattern WithUVec3)
+import Geomancy.UVec4 (UVec4, uvec4, withUVec4, pattern WithUVec4)
+
+import Geomancy.Mat4 (Mat4)
+import Geomancy.Transform (Transform(..))
+
+import Geomancy.Quaternion (Quaternion, quaternion, withQuaternion)
diff --git a/src/Geomancy/IVec2.hs b/src/Geomancy/IVec2.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/IVec2.hs
@@ -0,0 +1,95 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V2 Int32@.
+
+module Geomancy.IVec2
+  ( IVec2
+  , ivec2
+  , withIVec2
+  , pattern WithIVec2
+  , fromTuple
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Int (Int32)
+import Foreign (Storable(..))
+
+data IVec2 = IVec2
+  {-# UNPACK #-} !Int32
+  {-# UNPACK #-} !Int32
+  deriving (Eq, Ord, Show)
+
+{-# INLINE ivec2 #-}
+ivec2 :: Int32 -> Int32 -> IVec2
+ivec2 = IVec2
+
+{-# INLINE withIVec2 #-}
+withIVec2
+  :: IVec2
+  -> (Int32 -> Int32 -> r)
+  -> r
+withIVec2 (IVec2 a b) f = f a b
+
+pattern WithIVec2 :: Int32 -> Int32 -> IVec2
+pattern WithIVec2 a b <- ((`withIVec2` (,)) -> (a, b))
+{-# COMPLETE WithIVec2 #-}
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Int32, Int32) -> IVec2
+fromTuple (x, y) = ivec2 x y
+
+instance NFData IVec2 where
+  rnf IVec2{} = ()
+
+-- XXX: That's one nasty instance...
+instance Num IVec2 where
+  {-# INLINE (+) #-}
+  IVec2 l1 l2 + IVec2 r1 r2 =
+    IVec2
+      (l1 + r1)
+      (l2 + r2)
+
+  {-# INLINE (-) #-}
+  IVec2 l1 l2 - IVec2 r1 r2 =
+    IVec2
+      (l1 - r1)
+      (l2 - r2)
+
+  {-# INLINE (*) #-}
+  IVec2 l1 l2 * IVec2 r1 r2 =
+    IVec2
+      (l1 * r1)
+      (l2 * r2)
+
+  {-# INLINE abs #-}
+  abs (IVec2 a b) =
+    IVec2 (abs a) (abs b)
+
+  {-# INLINE signum #-}
+  signum v2 = withIVec2 v2 \a b ->
+    ivec2 (signum a) (signum b)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = IVec2 x' x'
+    where
+      x' = fromInteger x
+
+instance Storable IVec2 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 8
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr v4 =
+    withIVec2 v4 \a b -> do
+      pokeByteOff ptr 0 a
+      pokeByteOff ptr 4 b
+
+  {-# INLINE peek #-}
+  peek ptr = ivec2
+    <$> peekByteOff ptr 0
+    <*> peekByteOff ptr 4
diff --git a/src/Geomancy/IVec3.hs b/src/Geomancy/IVec3.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/IVec3.hs
@@ -0,0 +1,131 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V2 Int32@.
+
+module Geomancy.IVec3
+  ( IVec3
+  , ivec3
+  , withIVec3
+  , pattern WithIVec3
+  , fromTuple
+
+  , Packed(..)
+  , packed
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Int (Int32)
+import Foreign (Storable(..))
+
+data IVec3 = IVec3
+  {-# UNPACK #-} !Int32
+  {-# UNPACK #-} !Int32
+  {-# UNPACK #-} !Int32
+  deriving (Eq, Ord, Show)
+
+{-# INLINE ivec3 #-}
+ivec3 :: Int32 -> Int32 -> Int32 -> IVec3
+ivec3 = IVec3
+
+{-# INLINE withIVec3 #-}
+withIVec3
+  :: IVec3
+  -> (Int32 -> Int32 -> Int32 -> r)
+  -> r
+withIVec3 (IVec3 a b c) f = f a b c
+
+pattern WithIVec3 :: Int32 -> Int32 -> Int32 -> IVec3
+pattern WithIVec3 a b c <- ((`withIVec3` (,,)) -> (a, b, c))
+{-# COMPLETE WithIVec3 #-}
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Int32, Int32, Int32) -> IVec3
+fromTuple (a, b, c) = ivec3 a b c
+
+instance NFData IVec3 where
+  rnf IVec3{} = ()
+
+-- XXX: That's another nasty instance...
+instance Num IVec3 where
+  {-# INLINE (+) #-}
+  IVec3 l1 l2 l3 + IVec3 r1 r2 r3 =
+    IVec3
+      (l1 + r1)
+      (l2 + r2)
+      (l3 + r3)
+
+  {-# INLINE (-) #-}
+  IVec3 l1 l2 l3 - IVec3 r1 r2 r3 =
+    IVec3
+      (l1 - r1)
+      (l2 - r2)
+      (l3 - r3)
+
+  {-# INLINE (*) #-}
+  IVec3 l1 l2 l3 * IVec3 r1 r2 r3 =
+    IVec3
+      (l1 * r1)
+      (l2 * r2)
+      (l3 * r3)
+
+  {-# INLINE abs #-}
+  abs x = x
+
+  {-# INLINE signum #-}
+  signum v3 = withIVec3 v3 \a b c ->
+    ivec3 (signum a) (signum b) (signum c)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = IVec3 x' x' x'
+    where
+      x' = fromInteger x
+
+instance Storable IVec3 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 16
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr v3 =
+    withIVec3 v3 \a b c -> do
+      pokeByteOff ptr  0 a
+      pokeByteOff ptr  4 b
+      pokeByteOff ptr  8 c
+
+  {-# INLINE peek #-}
+  peek ptr = ivec3
+    <$> peekByteOff ptr  0
+    <*> peekByteOff ptr  4
+    <*> peekByteOff ptr  8
+
+newtype Packed = Packed { unPacked :: IVec3 }
+  deriving (Eq, Ord, Show, NFData, Num)
+
+{-# INLINE packed #-}
+packed :: Int32 -> Int32 -> Int32 -> Packed
+packed a b c = Packed (ivec3 a b c)
+
+instance Storable Packed where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 12
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr (Packed v3) =
+    withIVec3 v3 \a b c -> do
+      pokeByteOff ptr 0 a
+      pokeByteOff ptr 4 b
+      pokeByteOff ptr 8 c
+
+  {-# INLINE peek #-}
+  peek ptr = packed
+    <$> peekByteOff ptr 0
+    <*> peekByteOff ptr 4
+    <*> peekByteOff ptr 8
diff --git a/src/Geomancy/IVec4.hs b/src/Geomancy/IVec4.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/IVec4.hs
@@ -0,0 +1,106 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V2 Int32@.
+
+module Geomancy.IVec4
+  ( IVec4
+  , ivec4
+  , withIVec4
+  , pattern WithIVec4
+  , fromTuple
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Int (Int32)
+import Foreign (Storable(..))
+
+data IVec4 = IVec4
+  {-# UNPACK #-} !Int32
+  {-# UNPACK #-} !Int32
+  {-# UNPACK #-} !Int32
+  {-# UNPACK #-} !Int32
+  deriving (Eq, Ord, Show)
+
+{-# INLINE ivec4 #-}
+ivec4 :: Int32 -> Int32 -> Int32 -> Int32 -> IVec4
+ivec4 = IVec4
+
+{-# INLINE withIVec4 #-}
+withIVec4
+  :: IVec4
+  -> (Int32 -> Int32 -> Int32 -> Int32 -> r)
+  -> r
+withIVec4 (IVec4 a b c d) f = f a b c d
+
+pattern WithIVec4 :: Int32 -> Int32 -> Int32 -> Int32 -> IVec4
+pattern WithIVec4 a b c d <- ((`withIVec4` (,,,)) -> (a, b, c, d))
+{-# COMPLETE WithIVec4 #-}
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Int32, Int32, Int32, Int32) -> IVec4
+fromTuple (x, y, z, w) = ivec4 x y z w
+
+instance NFData IVec4 where
+  rnf IVec4{} = ()
+
+-- XXX: That's another nasty instance...
+instance Num IVec4 where
+  {-# INLINE (+) #-}
+  IVec4 l1 l2 l3 l4 + IVec4 r1 r2 r3 r4 =
+    IVec4
+      (l1 + r1)
+      (l2 + r2)
+      (l3 + r3)
+      (l4 + r4)
+
+  {-# INLINE (-) #-}
+  IVec4 l1 l2 l3 l4 - IVec4 r1 r2 r3 r4 =
+    IVec4
+      (l1 - r1)
+      (l2 - r2)
+      (l3 - r3)
+      (l4 - r4)
+
+  {-# INLINE (*) #-}
+  IVec4 l1 l2 l3 l4 * IVec4 r1 r2 r3 r4 =
+    IVec4
+      (l1 * r1)
+      (l2 * r2)
+      (l3 * r3)
+      (l4 * r4)
+
+  {-# INLINE abs #-}
+  abs x = x
+
+  {-# INLINE signum #-}
+  signum v4 = withIVec4 v4 \a b c d ->
+    ivec4 (signum a) (signum b) (signum c) (signum d)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = IVec4 x' x' x' x'
+    where
+      x' = fromInteger x
+
+instance Storable IVec4 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 16
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr v4 =
+    withIVec4 v4 \a b c d -> do
+      pokeByteOff ptr  0 a
+      pokeByteOff ptr  4 b
+      pokeByteOff ptr  8 c
+      pokeByteOff ptr 12 d
+
+  {-# INLINE peek #-}
+  peek ptr = ivec4
+    <$> peekByteOff ptr  0
+    <*> peekByteOff ptr  4
+    <*> peekByteOff ptr  8
+    <*> peekByteOff ptr 12
diff --git a/src/Geomancy/Mat4.hs b/src/Geomancy/Mat4.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Mat4.hs
@@ -0,0 +1,498 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MagicHash #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE UnboxedTuples #-}
+{-# LANGUAGE UnliftedFFITypes #-}
+
+-- | General matrix storage and operations.
+
+module Geomancy.Mat4
+  ( Mat4
+
+  , rowMajor
+  , withRowMajor
+  , toListRowMajor
+  , toListRowMajor2d
+  , fromRowMajor2d
+
+  , colMajor
+  , withColMajor
+  , toListColMajor
+  , toListColMajor2d
+
+  , identity
+  , transpose
+  , inverse
+  , pointwise
+  , zipWith
+  , matrixProduct
+  , scalarMultiply
+  , (!*)
+  ) where
+
+import Prelude hiding (zipWith)
+import GHC.Exts hiding (toList)
+
+import Control.DeepSeq (NFData(rnf))
+import Foreign (Storable(..))
+import GHC.IO (IO(..))
+import System.IO.Unsafe (unsafePerformIO)
+import Text.Printf (printf)
+
+import qualified Data.Foldable as Foldable
+import qualified Data.List as List
+
+import Geomancy.Vec4 (Vec4, vec4, withVec4)
+
+data Mat4 = Mat4 ByteArray#
+
+{- | Construct 'Mat4' from @row@ notation.
+-}
+rowMajor
+  :: Coercible Mat4 a
+  => Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> a
+rowMajor = coerce mat4
+
+{- | Reduce 'Mat4' with a function with @row@ notation of arguments.
+-}
+withRowMajor
+  :: Coercible a Mat4
+  => a
+  ->
+    ( Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      r
+    )
+  -> r
+withRowMajor m = withMat4 (coerce m)
+
+toListRowMajor :: Coercible a Mat4 => a -> [Float]
+toListRowMajor = toList . coerce
+
+toListRowMajor2d :: Coercible a Mat4 => a -> [[Float]]
+toListRowMajor2d = toList2d . coerce
+
+{- |
+  Build a Mat4 from a list-of-lists kind of container
+  with row-major ordering of elements.
+
+@
+  fromRowMajor2d (Linear.mkTransformation dir pos) :: Transform
+@
+-}
+fromRowMajor2d
+  :: forall t a
+  .  ( Foldable t
+     , Coercible Mat4 a
+     )
+  => t (t Float)
+  -> Maybe a
+fromRowMajor2d rows =
+  case Foldable.toList rows of
+    [r0, r1, r2, r3] ->
+      withRow r0 \m00 m01 m02 m03 ->
+      withRow r1 \m10 m11 m12 m13 ->
+      withRow r2 \m20 m21 m22 m23 ->
+      withRow r3 \m30 m31 m32 m33 ->
+        Just . coerce $ mat4
+          m00 m01 m02 m03
+          m10 m11 m12 m13
+          m20 m21 m22 m23
+          m30 m31 m32 m33
+    _ ->
+      Nothing
+  where
+    withRow row f =
+      case Foldable.toList row of
+        [c0, c1, c2, c3] ->
+          f c0 c1 c2 c3
+        _ ->
+          Nothing
+
+{- | Construct a 'Mat4' from @column@ notation.
+-}
+{-# INLINE colMajor #-}
+colMajor
+  :: Coercible Mat4 a
+  => Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> a
+colMajor
+  m00 m01 m02 m03
+  m10 m11 m12 m13
+  m20 m21 m22 m23
+  m30 m31 m32 m33 =
+    coerce $ mat4
+      m00 m10 m20 m30
+      m01 m11 m21 m31
+      m02 m12 m22 m32
+      m03 m13 m23 m33
+
+{- | Reduce 'Mat4' with a function with @column@ notation for arguments.
+-}
+{-# INLINE withColMajor #-}
+withColMajor
+  :: Coercible a Mat4
+  => a
+  ->
+    ( Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      r
+    )
+  -> r
+withColMajor m f = withMat4 (coerce m)
+  \ m00 m01 m02 m03
+    m10 m11 m12 m13
+    m20 m21 m22 m23
+    m30 m31 m32 m33 ->
+  f
+    m00 m10 m20 m30
+    m01 m11 m21 m31
+    m02 m12 m22 m32
+    m03 m13 m23 m33
+
+toListColMajor :: Coercible a Mat4 => a -> [Float]
+toListColMajor = toListTrans . coerce
+
+toListColMajor2d :: Coercible a Mat4 => a -> [[Float]]
+toListColMajor2d = toList2dTrans . coerce
+
+{- | Construct 'Mat4' from elements in memory order.
+-}
+{-# INLINE mat4 #-}
+mat4
+  :: Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Mat4
+mat4
+  (F# m00) (F# m01) (F# m02) (F# m03)
+  (F# m10) (F# m11) (F# m12) (F# m13)
+  (F# m20) (F# m21) (F# m22) (F# m23)
+  (F# m30) (F# m31) (F# m32) (F# m33) =
+  runRW# \world ->
+    let
+      !(# world_, arr #) = newAlignedPinnedByteArray# 64# 16# world
+
+      world00 = writeFloatArray# arr 0x0# m00 world_
+      world01 = writeFloatArray# arr 0x1# m01 world00
+      world02 = writeFloatArray# arr 0x2# m02 world01
+      world03 = writeFloatArray# arr 0x3# m03 world02
+
+      world10 = writeFloatArray# arr 0x4# m10 world03
+      world11 = writeFloatArray# arr 0x5# m11 world10
+      world12 = writeFloatArray# arr 0x6# m12 world11
+      world13 = writeFloatArray# arr 0x7# m13 world12
+
+      world20 = writeFloatArray# arr 0x8# m20 world13
+      world21 = writeFloatArray# arr 0x9# m21 world20
+      world22 = writeFloatArray# arr 0xA# m22 world21
+      world23 = writeFloatArray# arr 0xB# m23 world22
+
+      world30 = writeFloatArray# arr 0xC# m30 world23
+      world31 = writeFloatArray# arr 0xD# m31 world30
+      world32 = writeFloatArray# arr 0xE# m32 world31
+      world33 = writeFloatArray# arr 0xF# m33 world32
+
+      !(# _world', arr' #) = unsafeFreezeByteArray# arr world33
+    in
+      Mat4 arr'
+
+{- | Reduce 'Mat4' with a function with @memory@ notation for arguments.
+-}
+{-# INLINE withMat4 #-}
+withMat4
+  :: Mat4
+  ->
+    ( Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      r
+    )
+  -> r
+withMat4 (Mat4 arr) f =
+  f
+    (F# (indexFloatArray# arr 0x0#))
+    (F# (indexFloatArray# arr 0x1#))
+    (F# (indexFloatArray# arr 0x2#))
+    (F# (indexFloatArray# arr 0x3#))
+
+    (F# (indexFloatArray# arr 0x4#))
+    (F# (indexFloatArray# arr 0x5#))
+    (F# (indexFloatArray# arr 0x6#))
+    (F# (indexFloatArray# arr 0x7#))
+
+    (F# (indexFloatArray# arr 0x8#))
+    (F# (indexFloatArray# arr 0x9#))
+    (F# (indexFloatArray# arr 0xA#))
+    (F# (indexFloatArray# arr 0xB#))
+
+    (F# (indexFloatArray# arr 0xC#))
+    (F# (indexFloatArray# arr 0xD#))
+    (F# (indexFloatArray# arr 0xE#))
+    (F# (indexFloatArray# arr 0xF#))
+
+{- | @I@, the identity matrix.
+
+Neutral element of its monoid, so you can use 'mempty'.
+-}
+{-# INLINE identity #-}
+identity :: Mat4
+identity = mat4
+  1 0 0 0
+  0 1 0 0
+  0 0 1 0
+  0 0 0 1
+
+{-# INLINE transpose #-}
+transpose :: Mat4 -> Mat4
+transpose =
+  flip withMat4
+    \ m00 m01 m02 m03
+      m10 m11 m12 m13
+      m20 m21 m22 m23
+      m30 m31 m32 m33 ->
+    mat4
+      m00 m10 m20 m30
+      m01 m11 m21 m31
+      m02 m12 m22 m32
+      m03 m13 m23 m33
+
+{- | Compute an inverse matrix, slowly.
+-}
+inverse :: (Coercible Mat4 a, Coercible Mat4 a) => a -> a
+inverse m =
+  coerce $ withMat4 (coerce m)
+    \ m00 m01 m02 m03
+      m10 m11 m12 m13
+      m20 m21 m22 m23
+      m30 m31 m32 m33 ->
+        let
+          invDet = recip det
+
+          det
+            = s0 * c5
+            - s1 * c4
+            + s2 * c3
+            + s3 * c2
+            - s4 * c1
+            + s5 * c0
+
+          s0 = m00 * m11 - m10 * m01
+          s1 = m00 * m12 - m10 * m02
+          s2 = m00 * m13 - m10 * m03
+          s3 = m01 * m12 - m11 * m02
+          s4 = m01 * m13 - m11 * m03
+          s5 = m02 * m13 - m12 * m03
+
+          c5 = m22 * m33 - m32 * m23
+          c4 = m21 * m33 - m31 * m23
+          c3 = m21 * m32 - m31 * m22
+          c2 = m20 * m33 - m30 * m23
+          c1 = m20 * m32 - m30 * m22
+          c0 = m20 * m31 - m30 * m21
+
+          i00 = ( m11 * c5 - m12 * c4 + m13 * c3) * invDet
+          i01 = (-m01 * c5 + m02 * c4 - m03 * c3) * invDet
+          i02 = ( m31 * s5 - m32 * s4 + m33 * s3) * invDet
+          i03 = (-m21 * s5 + m22 * s4 - m23 * s3) * invDet
+
+          i10 = (-m10 * c5 + m12 * c2 - m13 * c1) * invDet
+          i11 = ( m00 * c5 - m02 * c2 + m03 * c1) * invDet
+          i12 = (-m30 * s5 + m32 * s2 - m33 * s1) * invDet
+          i13 = ( m20 * s5 - m22 * s2 + m23 * s1) * invDet
+
+          i20 = ( m10 * c4 - m11 * c2 + m13 * c0) * invDet
+          i21 = (-m00 * c4 + m01 * c2 - m03 * c0) * invDet
+          i22 = ( m30 * s4 - m31 * s2 + m33 * s0) * invDet
+          i23 = (-m20 * s4 + m21 * s2 - m23 * s0) * invDet
+
+          i30 = (-m10 * c3 + m11 * c1 - m12 * c0) * invDet
+          i31 = ( m00 * c3 - m01 * c1 + m02 * c0) * invDet
+          i32 = (-m30 * s3 + m31 * s1 - m32 * s0) * invDet
+          i33 = ( m20 * s3 - m21 * s1 + m22 * s0) * invDet
+        in
+          mat4
+            i00 i01 i02 i03
+            i10 i11 i12 i13
+            i20 i21 i22 i23
+            i30 i31 i32 i33
+
+pointwise :: Mat4 -> Mat4 -> (Float -> Float -> Float) -> Mat4
+pointwise a b f =
+  withMat4 a
+    \ a00 a01 a02 a03
+      a10 a11 a12 a13
+      a20 a21 a22 a23
+      a30 a31 a32 a33 ->
+  withMat4 b
+    \ b00 b01 b02 b03
+      b10 b11 b12 b13
+      b20 b21 b22 b23
+      b30 b31 b32 b33 ->
+  mat4
+    (f a00 b00) (f a01 b01) (f a02 b02) (f a03 b03)
+    (f a10 b10) (f a11 b11) (f a12 b12) (f a13 b13)
+    (f a20 b20) (f a21 b21) (f a22 b22) (f a23 b23)
+    (f a30 b30) (f a31 b31) (f a32 b32) (f a33 b33)
+
+toList :: Mat4 -> [Float]
+toList = flip withMat4
+    \ a00 a01 a02 a03
+      a10 a11 a12 a13
+      a20 a21 a22 a23
+      a30 a31 a32 a33 ->
+    [ a00, a01, a02, a03
+    , a10, a11, a12, a13
+    , a20, a21, a22, a23
+    , a30, a31, a32, a33
+    ]
+
+toList2d :: Mat4 -> [[Float]]
+toList2d = flip withMat4
+    \ a00 a01 a02 a03
+      a10 a11 a12 a13
+      a20 a21 a22 a23
+      a30 a31 a32 a33 ->
+    [ [a00, a01, a02, a03]
+    , [a10, a11, a12, a13]
+    , [a20, a21, a22, a23]
+    , [a30, a31, a32, a33]
+    ]
+
+toListTrans :: Mat4 -> [Float]
+toListTrans = flip withMat4
+    \ a00 a01 a02 a03
+      a10 a11 a12 a13
+      a20 a21 a22 a23
+      a30 a31 a32 a33 ->
+    [ a00, a10, a20, a30
+    , a01, a11, a21, a31
+    , a02, a12, a22, a32
+    , a03, a13, a23, a33
+    ]
+
+toList2dTrans :: Mat4 -> [[Float]]
+toList2dTrans = flip withMat4
+    \ a00 a01 a02 a03
+      a10 a11 a12 a13
+      a20 a21 a22 a23
+      a30 a31 a32 a33 ->
+    [ [a00, a10, a20, a30]
+    , [a01, a11, a21, a31]
+    , [a02, a12, a22, a32]
+    , [a03, a13, a23, a33]
+    ]
+
+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 ()
+
+{-# INLINE matrixProduct #-}
+matrixProduct :: Mat4 -> Mat4 -> Mat4
+matrixProduct (Mat4 l) (Mat4 r) = unsafePerformIO do
+  result@(Mat4 m) <- unsafeNewMat4
+  mm4sse
+    (byteArrayContents# l)
+    (byteArrayContents# r)
+    (byteArrayContents# m)
+  pure result
+
+{-# INLINE unsafeNewMat4 #-}
+unsafeNewMat4 :: IO Mat4
+unsafeNewMat4 =
+  IO \world ->
+    let
+      !(# world_, arr_ #) = newAlignedPinnedByteArray# 64# 16# world
+      !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_
+    in
+      (# world, Mat4 arr #)
+
+instance NFData Mat4 where
+  rnf Mat4{} = ()
+
+instance Semigroup Mat4 where
+  {-# INLINE (<>) #-}
+  (<>) = matrixProduct
+
+instance Monoid Mat4 where
+  {-# INLINE mempty #-}
+  mempty = identity
+
+instance Show Mat4 where
+  show = flip withMat4
+    \ m00 m01 m02 m03
+      m10 m11 m12 m13
+      m20 m21 m22 m23
+      m30 m31 m32 m33 ->
+    unlines
+      [ printf "| %.4f %.4f %.4f %.4f |" m00 m01 m02 m03
+      , printf "| %.4f %.4f %.4f %.4f |" m10 m11 m12 m13
+      , printf "| %.4f %.4f %.4f %.4f |" m20 m21 m22 m23
+      , printf "| %.4f %.4f %.4f %.4f |" m30 m31 m32 m33
+      ]
+
+instance Storable Mat4 where
+  sizeOf _mat4 = 64
+
+  alignment _mat4 = 16
+
+  {-# INLINE poke #-}
+  poke (Ptr addr) (Mat4 arr) = IO \world ->
+    let
+      world' = copyByteArrayToAddr# arr 0# addr 64# world
+    in
+      (# world', () #)
+
+  {-# INLINE peek #-}
+  peek (Ptr addr) = IO \world ->
+    let
+      !(# world0, arr #)  = newAlignedPinnedByteArray# 64# 16# world
+      world1              = copyAddrToByteArray# addr arr 0# 64# world0
+      !(# world', arr' #) = unsafeFreezeByteArray# arr world1
+    in
+      (# world', Mat4 arr' #)
diff --git a/src/Geomancy/Quaternion.hs b/src/Geomancy/Quaternion.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Quaternion.hs
@@ -0,0 +1,265 @@
+{-# LANGUAGE BlockArguments #-}
+
+-- | Specialized and inlined @Quaternion Float@.
+
+module Geomancy.Quaternion
+  ( Quaternion
+  , quaternion
+  , withQuaternion
+
+  , axisAngle
+  , rotate
+  , rotatePoint
+  , rotationBetween
+  , lookAtUp
+
+  , (^*)
+  , (^/)
+  , slerp
+
+  , conjugate
+  , norm
+  , quadrance
+  , dot
+  , normalize
+  , qNaN
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Foreign (Storable(..), castPtr)
+
+import Geomancy.Vec3 (Vec3, vec3, withVec3)
+
+import qualified Geomancy.Vec3 as Vec3
+
+data Quaternion = Quaternion
+  {-# UNPACK #-} !Float
+  {-# UNPACK #-} !Float
+  {-# UNPACK #-} !Float
+  {-# UNPACK #-} !Float
+  deriving (Eq, Ord, Show)
+
+{-# INLINE quaternion #-}
+quaternion :: Float -> Float -> Float -> Float -> Quaternion
+quaternion = Quaternion
+
+{-# INLINE withQuaternion #-}
+withQuaternion
+  :: Quaternion
+  -> (Float -> Float -> Float -> Float -> r)
+  -> r
+withQuaternion (Quaternion a b c d) f = f a b c d
+
+{-# INLINE (^*) #-}
+(^*) :: Quaternion -> Float -> Quaternion
+Quaternion a b c d ^* x =
+  Quaternion
+    (a * x)
+    (b * x)
+    (c * x)
+    (d * x)
+
+{-# INLINE (^/) #-}
+(^/) :: Quaternion -> Float -> Quaternion
+Quaternion a b c d ^/ x =
+  Quaternion
+    (a / x)
+    (b / x)
+    (c / x)
+    (d / x)
+
+slerp :: Quaternion -> Quaternion -> Float -> Quaternion
+slerp q p t
+  | 1.0 - cosphi < 1e-8 =
+      q
+  | otherwise =
+      ( (q   ^* sin ((1 - t) * phi)) +
+         f p ^* sin (t * phi)
+      ) ^/ sin phi
+  where
+    phi = acos cosphi
+
+    (cosphi, f) =
+      if dqp < 0 then
+        (-dqp, negate)
+      else
+        (dqp, id)
+
+    dqp = dot q p
+
+{-# INLINE conjugate #-}
+conjugate :: Quaternion -> Quaternion
+conjugate (Quaternion e x y z) = Quaternion e (-x) (-y) (-z)
+
+{-# INLINE norm #-}
+norm :: Quaternion -> Float
+norm = sqrt . quadrance
+
+{-# INLINE quadrance #-}
+quadrance :: Quaternion -> Float
+quadrance q = dot q q
+
+{-# INLINE dot #-}
+dot :: Quaternion -> Quaternion -> Float
+dot (Quaternion a b c d) (Quaternion e f g h) =
+  a * e +
+  b * f +
+  c * g +
+  d * h -- XXX: SIMD time!
+
+{-# INLINE normalize #-}
+normalize :: Quaternion -> Quaternion
+normalize v =
+  if nearZero q || nearZero (1-q) then
+    v
+  else
+    let
+      Quaternion e i j k = v
+    in
+      Quaternion (e / l) (i / l) (j / l) (k / l)
+
+  where
+    q = dot v v
+    l = sqrt q
+
+    nearZero a = abs a <= 1e-6
+
+instance NFData Quaternion where
+  rnf Quaternion{} = ()
+
+instance Num Quaternion where
+  {-# INLINE (+) #-}
+  Quaternion a b c d + Quaternion e f g h =
+    Quaternion
+      (a + e)
+      (b + f)
+      (c + g)
+      (d + h)
+
+  {-# INLINE (-) #-}
+  Quaternion a b c d - Quaternion e f g h =
+    Quaternion
+      (a - e)
+      (b - f)
+      (c - g)
+      (d - h)
+
+  {-# INLINE (*) #-}
+  Quaternion a b c d * Quaternion e f g h =
+    withVec3 v \y z w ->
+      Quaternion x y z w
+    where
+      x = a * e - Vec3.dot v1 v2
+      v = Vec3.cross v1 v2 + v2 Vec3.^* a + v1 Vec3.^* e
+      v1 = vec3 b c d
+      v2 = vec3 f g h
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = Quaternion (fromInteger x) 0 0 0
+
+  {-# INLINE abs #-}
+  abs z = Quaternion (norm z) 0 0 0
+
+  {-# INLINE signum #-}
+  signum q@(Quaternion e i j k)
+    | m == 0 =
+        q
+    | not (isInfinite m || isNaN m) =
+        Quaternion (e * misqrt) (i * misqrt) (j * misqrt) (k * misqrt)
+    | any isNaN [e, i, j, k] = qNaN
+    | not (ii || ij || ik) = Quaternion 1 0 0 0
+    | not (ie || ij || ik) = Quaternion 0 1 0 0
+    | not (ie || ii || ik) = Quaternion 0 0 1 0
+    | not (ie || ii || ij) = Quaternion 0 0 0 1
+    | otherwise = qNaN
+    where
+      m = quadrance q
+      misqrt = recip (sqrt m)
+
+      ie = isInfinite e
+      ii = isInfinite i
+      ij = isInfinite j
+      ik = isInfinite k
+
+{-# INLINE qNaN #-}
+qNaN :: Quaternion
+qNaN = Quaternion fNaN fNaN fNaN fNaN
+  where
+    fNaN = 0/0
+
+-- XXX: GPU layouts call for some padding.
+instance Storable Quaternion where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 16
+
+  {-# INLINE alignment #-}
+  alignment _ = 16
+
+  {-# INLINE poke #-}
+  poke ptr (Quaternion a b c d) = do
+    poke ptr' a
+    pokeElemOff ptr' 1 b
+    pokeElemOff ptr' 2 c
+    pokeElemOff ptr' 3 d
+    where
+      ptr' = castPtr ptr
+
+  {-# INLINE peek #-}
+  peek ptr = Quaternion
+    <$> peek ptr'
+    <*> peekElemOff ptr' 1
+    <*> peekElemOff ptr' 2
+    <*> peekElemOff ptr' 3
+    where
+      ptr' = castPtr ptr
+
+-- | Quaternion construction from axis and angle.
+{-# INLINE axisAngle #-}
+axisAngle :: Vec3 -> Float -> Quaternion
+axisAngle axis rads =
+  withVec3 (Vec3.normalize axis Vec3.^* sin half) $
+    quaternion (cos half)
+  where
+    half = rads / 2
+
+{-# INLINE rotate #-}
+rotate :: Quaternion -> Vec3 -> Vec3
+rotate q v = withQuaternion q' \_a b c d -> vec3 b c d
+  where
+    q' = withVec3 v \x y z ->
+      q * quaternion 0 x y z * conjugate q
+
+{-# INLINE rotatePoint #-}
+rotatePoint :: Quaternion -> Vec3 -> Vec3 -> Vec3
+rotatePoint q origin point =
+  origin + rotate q (point - origin)
+
+{- | Rotation between vectors.
+
+(in other words: the quaternion needed to rotate @v1@ so that it matches @v2@)
+-}
+rotationBetween :: Vec3 -> Vec3 -> Quaternion
+rotationBetween v1 v2 = axisAngle axis angle
+  where
+    axis = Vec3.cross v1 v2
+    angle = acos cosAngle
+    cosAngle =
+      max (-1) . min 1 $
+        Vec3.dot (Vec3.normalize v1) (Vec3.normalize v2)
+
+{- | Orient towards a point.
+
+Use "rotationBetween" if you don't need to keep the object upright.
+-}
+lookAtUp :: Vec3 -> Vec3 -> Vec3 -> Quaternion
+lookAtUp src dst up = rot2 * rot1
+  where
+    dir3 = dst - src
+
+    -- XXX: turn "eye"
+    rot1 = rotationBetween (vec3 0 0 1) dir3
+
+    rot2 = rotationBetween newUp fixedUp
+
+    newUp = rotate rot1 up
+    fixedUp = Vec3.cross (Vec3.cross dir3 up) dir3
diff --git a/src/Geomancy/Transform.hs b/src/Geomancy/Transform.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Transform.hs
@@ -0,0 +1,178 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralisedNewtypeDeriving #-}
+
+module Geomancy.Transform
+  ( Transform(..)
+  , inverse
+
+  , apply
+  , (!.)
+
+  , translate
+  , translateV
+
+  , rotateX
+  , rotateY
+  , rotateZ
+  , rotateQ
+
+  , scale
+  , scaleX
+  , scaleY
+  , scaleZ
+  , scaleXY
+  , scale3
+
+  , dirPos
+  ) where
+
+import Foreign (Storable(..))
+
+import Geomancy.Mat4 (Mat4, colMajor, withColMajor, inverse)
+import Geomancy.Quaternion (Quaternion, withQuaternion)
+import Geomancy.Vec3 (Vec3, vec3, withVec3)
+
+newtype Transform = Transform { unTransform :: Mat4 }
+  deriving newtype (Show, Semigroup, Monoid, Storable)
+
+-- | Apply transformation to a vector, then normalize with perspective division
+apply :: Vec3 -> Transform -> Vec3
+apply = flip (!.)
+
+-- | 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)
+
+-- ** Translation
+
+{-# INLINE translate #-}
+translate :: Float -> Float -> Float -> Transform
+translate x y z = colMajor
+  1 0 0 x
+  0 1 0 y
+  0 0 1 z
+  0 0 0 1
+
+{-# INLINE translateV #-}
+translateV :: Vec3 -> Transform
+translateV vec = withVec3 vec translate
+
+-- ** Scaling
+
+{-# INLINE scale3 #-}
+scale3 :: Float -> Float -> Float -> Transform
+scale3 x y z = colMajor
+  x 0 0 0
+  0 y 0 0
+  0 0 z 0
+  0 0 0 1
+
+{-# INLINE scale #-}
+scale :: Float -> Transform
+scale s = scale3 s s s
+
+{-# INLINE scaleX #-}
+scaleX :: Float -> Transform
+scaleX x = scale3 x 1 1
+
+{-# INLINE scaleY #-}
+scaleY :: Float -> Transform
+scaleY y = scale3 1 y 1
+
+{-# INLINE scaleZ #-}
+scaleZ :: Float -> Transform
+scaleZ z = scale3 1 1 z
+
+{-# INLINE scaleXY #-}
+scaleXY :: Float -> Float -> Transform
+scaleXY x y = scale3 x y 1
+
+-- ** Euler angle rotations
+
+{-# INLINE rotateX #-}
+rotateX :: Float -> Transform
+rotateX rads = colMajor
+  1 0   0   0
+  0 t11 t21 0
+  0 t12 t22 0
+  0 0   0   1
+  where
+    t11 = cost
+    t12 = -sint
+    t21 = sint
+    t22 = cost
+
+    cost = cos rads
+    sint = sin rads
+
+{-# INLINE rotateY #-}
+rotateY :: Float -> Transform
+rotateY rads = colMajor
+  t00 0 t20 0
+  0   1 0   0
+  t02 0 t22 0
+  0   0 0   1
+  where
+    cost = cos rads
+    sint = sin rads
+
+    t00 = cost
+    t02 = sint
+    t20 = -sint
+    t22 = cost
+
+{-# INLINE rotateZ #-}
+rotateZ :: Float -> Transform
+rotateZ rads = colMajor
+  t00 t10 0 0
+  t01 t11 0 0
+  0   0   1 0
+  0   0   0 1
+  where
+   t00 = cost
+   t01 = -sint
+   t10 = sint
+   t11 = cost
+
+   cost = cos rads
+   sint = sin rads
+
+{-# INLINE rotateQ #-}
+rotateQ :: Quaternion -> Transform
+rotateQ dir = dirPos dir 0
+
+{-# INLINE dirPos #-}
+dirPos :: Quaternion -> Vec3 -> Transform
+dirPos rs t =
+  withQuaternion rs \w x y z ->
+  withVec3 t \tx ty tz ->
+    let
+      x2 = x * x
+      y2 = y * y
+      z2 = z * z
+      xy = x * y
+      xz = x * z
+      xw = x * w
+      yz = y * z
+      yw = y * w
+      zw = z * w
+    in
+      colMajor
+        (1 - 2 * (y2 + z2)) (    2 * (xy - zw)) (    2 * (xz + yw)) tx
+        (    2 * (xy + zw)) (1 - 2 * (x2 + z2)) (    2 * (yz - xw)) ty
+        (    2 * (xz - yw)) (    2 * (yz + xw)) (1 - 2 * (x2 + y2)) tz
+         0                   0                   0                  1
diff --git a/src/Geomancy/UVec2.hs b/src/Geomancy/UVec2.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/UVec2.hs
@@ -0,0 +1,95 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V2 Word32@.
+
+module Geomancy.UVec2
+  ( UVec2
+  , uvec2
+  , withUVec2
+  , pattern WithUVec2
+  , fromTuple
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Word (Word32)
+import Foreign (Storable(..))
+
+data UVec2 = UVec2
+  {-# UNPACK #-} !Word32
+  {-# UNPACK #-} !Word32
+  deriving (Eq, Ord, Show)
+
+{-# INLINE uvec2 #-}
+uvec2 :: Word32 -> Word32 -> UVec2
+uvec2 = UVec2
+
+{-# INLINE withUVec2 #-}
+withUVec2
+  :: UVec2
+  -> (Word32 -> Word32 -> r)
+  -> r
+withUVec2 (UVec2 a b) f = f a b
+
+pattern WithUVec2 :: Word32 -> Word32 -> UVec2
+pattern WithUVec2 a b <- ((`withUVec2` (,)) -> (a, b))
+{-# COMPLETE WithUVec2 #-}
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Word32, Word32) -> UVec2
+fromTuple (x, y) = uvec2 x y
+
+instance NFData UVec2 where
+  rnf UVec2{} = ()
+
+-- XXX: That's one nasty instance...
+instance Num UVec2 where
+  {-# INLINE (+) #-}
+  UVec2 l1 l2 + UVec2 r1 r2 =
+    UVec2
+      (l1 + r1)
+      (l2 + r2)
+
+  {-# INLINE (-) #-}
+  UVec2 l1 l2 - UVec2 r1 r2 =
+    UVec2
+      (l1 - r1)
+      (l2 - r2)
+
+  {-# INLINE (*) #-}
+  UVec2 l1 l2 * UVec2 r1 r2 =
+    UVec2
+      (l1 * r1)
+      (l2 * r2)
+
+  {-# INLINE abs #-}
+  abs (UVec2 a b) =
+    UVec2 (abs a) (abs b)
+
+  {-# INLINE signum #-}
+  signum v2 = withUVec2 v2 \a b ->
+    uvec2 (signum a) (signum b)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = UVec2 x' x'
+    where
+      x' = fromInteger x
+
+instance Storable UVec2 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 8
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr v4 =
+    withUVec2 v4 \a b -> do
+      pokeByteOff ptr 0 a
+      pokeByteOff ptr 4 b
+
+  {-# INLINE peek #-}
+  peek ptr = uvec2
+    <$> peekByteOff ptr 0
+    <*> peekByteOff ptr 4
diff --git a/src/Geomancy/UVec3.hs b/src/Geomancy/UVec3.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/UVec3.hs
@@ -0,0 +1,131 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V2 Word32@.
+
+module Geomancy.UVec3
+  ( UVec3
+  , uvec3
+  , withUVec3
+  , pattern WithUVec3
+  , fromTuple
+
+  , Packed(..)
+  , packed
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Word (Word32)
+import Foreign (Storable(..))
+
+data UVec3 = UVec3
+  {-# UNPACK #-} !Word32
+  {-# UNPACK #-} !Word32
+  {-# UNPACK #-} !Word32
+  deriving (Eq, Ord, Show)
+
+{-# INLINE uvec3 #-}
+uvec3 :: Word32 -> Word32 -> Word32 -> UVec3
+uvec3 = UVec3
+
+{-# INLINE withUVec3 #-}
+withUVec3
+  :: UVec3
+  -> (Word32 -> Word32 -> Word32 -> r)
+  -> r
+withUVec3 (UVec3 a b c) f = f a b c
+
+pattern WithUVec3 :: Word32 -> Word32 -> Word32 -> UVec3
+pattern WithUVec3 a b c <- ((`withUVec3` (,,)) -> (a, b, c))
+{-# COMPLETE WithUVec3 #-}
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Word32, Word32, Word32) -> UVec3
+fromTuple (a, b, c) = uvec3 a b c
+
+instance NFData UVec3 where
+  rnf UVec3{} = ()
+
+-- XXX: That's another nasty instance...
+instance Num UVec3 where
+  {-# INLINE (+) #-}
+  UVec3 l1 l2 l3 + UVec3 r1 r2 r3 =
+    UVec3
+      (l1 + r1)
+      (l2 + r2)
+      (l3 + r3)
+
+  {-# INLINE (-) #-}
+  UVec3 l1 l2 l3 - UVec3 r1 r2 r3 =
+    UVec3
+      (l1 - r1)
+      (l2 - r2)
+      (l3 - r3)
+
+  {-# INLINE (*) #-}
+  UVec3 l1 l2 l3 * UVec3 r1 r2 r3 =
+    UVec3
+      (l1 * r1)
+      (l2 * r2)
+      (l3 * r3)
+
+  {-# INLINE abs #-}
+  abs x = x
+
+  {-# INLINE signum #-}
+  signum v3 = withUVec3 v3 \a b c ->
+    uvec3 (signum a) (signum b) (signum c)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = UVec3 x' x' x'
+    where
+      x' = fromInteger x
+
+instance Storable UVec3 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 16
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr v3 =
+    withUVec3 v3 \a b c -> do
+      pokeByteOff ptr  0 a
+      pokeByteOff ptr  4 b
+      pokeByteOff ptr  8 c
+
+  {-# INLINE peek #-}
+  peek ptr = uvec3
+    <$> peekByteOff ptr  0
+    <*> peekByteOff ptr  4
+    <*> peekByteOff ptr  8
+
+newtype Packed = Packed { unPacked :: UVec3 }
+  deriving (Eq, Ord, Show, NFData, Num)
+
+{-# INLINE packed #-}
+packed :: Word32 -> Word32 -> Word32 -> Packed
+packed a b c = Packed (uvec3 a b c)
+
+instance Storable Packed where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 12
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr (Packed v3) =
+    withUVec3 v3 \a b c -> do
+      pokeByteOff ptr 0 a
+      pokeByteOff ptr 4 b
+      pokeByteOff ptr 8 c
+
+  {-# INLINE peek #-}
+  peek ptr = packed
+    <$> peekByteOff ptr 0
+    <*> peekByteOff ptr 4
+    <*> peekByteOff ptr 8
diff --git a/src/Geomancy/UVec4.hs b/src/Geomancy/UVec4.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/UVec4.hs
@@ -0,0 +1,106 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V2 Word32@.
+
+module Geomancy.UVec4
+  ( UVec4
+  , uvec4
+  , withUVec4
+  , pattern WithUVec4
+  , fromTuple
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Word (Word32)
+import Foreign (Storable(..))
+
+data UVec4 = UVec4
+  {-# UNPACK #-} !Word32
+  {-# UNPACK #-} !Word32
+  {-# UNPACK #-} !Word32
+  {-# UNPACK #-} !Word32
+  deriving (Eq, Ord, Show)
+
+{-# INLINE uvec4 #-}
+uvec4 :: Word32 -> Word32 -> Word32 -> Word32 -> UVec4
+uvec4 = UVec4
+
+{-# INLINE withUVec4 #-}
+withUVec4
+  :: UVec4
+  -> (Word32 -> Word32 -> Word32 -> Word32 -> r)
+  -> r
+withUVec4 (UVec4 a b c d) f = f a b c d
+
+pattern WithUVec4 :: Word32 -> Word32 -> Word32 -> Word32 -> UVec4
+pattern WithUVec4 a b c d <- ((`withUVec4` (,,,)) -> (a, b, c, d))
+{-# COMPLETE WithUVec4 #-}
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Word32, Word32, Word32, Word32) -> UVec4
+fromTuple (x, y, z, w) = uvec4 x y z w
+
+instance NFData UVec4 where
+  rnf UVec4{} = ()
+
+-- XXX: That's another nasty instance...
+instance Num UVec4 where
+  {-# INLINE (+) #-}
+  UVec4 l1 l2 l3 l4 + UVec4 r1 r2 r3 r4 =
+    UVec4
+      (l1 + r1)
+      (l2 + r2)
+      (l3 + r3)
+      (l4 + r4)
+
+  {-# INLINE (-) #-}
+  UVec4 l1 l2 l3 l4 - UVec4 r1 r2 r3 r4 =
+    UVec4
+      (l1 - r1)
+      (l2 - r2)
+      (l3 - r3)
+      (l4 - r4)
+
+  {-# INLINE (*) #-}
+  UVec4 l1 l2 l3 l4 * UVec4 r1 r2 r3 r4 =
+    UVec4
+      (l1 * r1)
+      (l2 * r2)
+      (l3 * r3)
+      (l4 * r4)
+
+  {-# INLINE abs #-}
+  abs x = x
+
+  {-# INLINE signum #-}
+  signum v4 = withUVec4 v4 \a b c d ->
+    uvec4 (signum a) (signum b) (signum c) (signum d)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = UVec4 x' x' x' x'
+    where
+      x' = fromInteger x
+
+instance Storable UVec4 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 16
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr v4 =
+    withUVec4 v4 \a b c d -> do
+      pokeByteOff ptr  0 a
+      pokeByteOff ptr  4 b
+      pokeByteOff ptr  8 c
+      pokeByteOff ptr 12 d
+
+  {-# INLINE peek #-}
+  peek ptr = uvec4
+    <$> peekByteOff ptr  0
+    <*> peekByteOff ptr  4
+    <*> peekByteOff ptr  8
+    <*> peekByteOff ptr 12
diff --git a/src/Geomancy/Vec2.hs b/src/Geomancy/Vec2.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Vec2.hs
@@ -0,0 +1,154 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V2 Float@.
+
+module Geomancy.Vec2
+  ( Vec2
+  , vec2
+  , withVec2
+  , pattern WithVec2
+  , fromTuple
+
+  , (^*)
+  , (^/)
+  , lerp
+
+  , dot
+  , normalize
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Foreign (Storable(..))
+
+data Vec2 = Vec2
+  {-# UNPACK #-} !Float
+  {-# UNPACK #-} !Float
+  deriving (Eq, Ord, Show)
+
+{-# INLINE vec2 #-}
+vec2 :: Float -> Float -> Vec2
+vec2 = Vec2
+
+{-# INLINE withVec2 #-}
+withVec2
+  :: Vec2
+  -> (Float -> Float -> r)
+  -> r
+withVec2 (Vec2 a b) f = f a b
+
+pattern WithVec2 :: Float -> Float -> Vec2
+pattern WithVec2 a b <- ((`withVec2` (,)) -> (a, b))
+{-# COMPLETE WithVec2 #-}
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Float, Float) -> Vec2
+fromTuple (x, y) = vec2 x y
+
+instance NFData Vec2 where
+  rnf Vec2{} = ()
+
+instance Num Vec2 where
+  {-# INLINE (+) #-}
+  Vec2 l1 l2 + Vec2 r1 r2 =
+    Vec2
+      (l1 + r1)
+      (l2 + r2)
+
+  {-# INLINE (-) #-}
+  Vec2 l1 l2 - Vec2 r1 r2 =
+    Vec2
+      (l1 - r1)
+      (l2 - r2)
+
+  {-# INLINE (*) #-}
+  Vec2 l1 l2 * Vec2 r1 r2 =
+    Vec2
+      (l1 * r1)
+      (l2 * r2)
+
+  {-# INLINE abs #-}
+  abs (Vec2 a b) =
+    Vec2 (abs a) (abs b)
+
+  {-# INLINE signum #-}
+  signum (Vec2 a b) =
+    Vec2 (signum a) (signum b)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = Vec2 x' x'
+    where
+      x' = fromInteger x
+
+instance Fractional Vec2 where
+  {-# INLINE (/) #-}
+  Vec2 l1 l2 / Vec2 r1 r2 =
+    Vec2 (l1 / r1) (l2 / r2)
+
+  {-# INLINE recip #-}
+  recip (Vec2 a b) =
+    Vec2 (recip a) (recip b)
+
+  {-# INLINE fromRational #-}
+  fromRational x = Vec2 x' x'
+    where
+      x' = fromRational x
+
+{-# INLINE (^*) #-}
+(^*) :: Vec2 -> Float -> Vec2
+Vec2 a b ^* x =
+  Vec2
+    (a * x)
+    (b * x)
+
+{-# INLINE (^/) #-}
+(^/) :: Vec2 -> Float -> Vec2
+Vec2 a b ^/ x =
+  Vec2
+    (a / x)
+    (b / x)
+
+{-# INLINE lerp #-}
+lerp :: Float -> Vec2 -> Vec2 -> Vec2
+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)
+
+{-# INLINE dot #-}
+dot :: Vec2 -> Vec2 -> Float
+dot (Vec2 l1 l2) (Vec2 r1 r2) =
+  l1 * r1 + l2 * r2
+
+{-# INLINE normalize #-}
+normalize :: Vec2 -> Vec2
+normalize v =
+  if nearZero q || nearZero (1 - q) then
+    v
+  else
+    let
+      Vec2 x y = v
+    in
+      Vec2 (x / l) (y / l)
+
+  where
+    q = dot v v
+    l = sqrt q
+
+    nearZero a = abs a <= 1e-6
+
+instance Storable Vec2 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 8
+
+  {-# INLINE alignment #-}
+  alignment _ = 8
+
+  {-# INLINE poke #-}
+  poke ptr v4 =
+    withVec2 v4 \a b -> do
+      pokeByteOff ptr 0 a
+      pokeByteOff ptr 4 b
+
+  {-# INLINE peek #-}
+  peek ptr = vec2
+    <$> peekByteOff ptr 0
+    <*> peekByteOff ptr 4
diff --git a/src/Geomancy/Vec3.hs b/src/Geomancy/Vec3.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Vec3.hs
@@ -0,0 +1,229 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V3 Float@.
+
+module Geomancy.Vec3
+  ( Vec3
+  , vec3
+  , withVec3
+  , pattern WithVec3
+  , fromVec2
+  , fromTuple
+
+  , (^*)
+  , (^/)
+  , lerp
+
+  , cross
+  , dot
+  , normalize
+
+  , Packed(..)
+  , packed
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Coerce (Coercible, coerce)
+import Foreign (Storable(..), castPtr)
+
+import Geomancy.Vec2 (Vec2, withVec2)
+
+data Vec3 = Vec3
+  {-# UNPACK #-} !Float
+  {-# UNPACK #-} !Float
+  {-# UNPACK #-} !Float
+  deriving (Eq, Ord, Show)
+
+{-# INLINE vec3 #-}
+vec3 :: Float -> Float -> Float -> Vec3
+vec3 = Vec3
+
+{-# INLINE withVec3 #-}
+withVec3
+  :: Vec3
+  -> (Float -> Float -> Float -> r)
+  -> r
+withVec3 (Vec3 a b c) f = f a b c
+
+pattern WithVec3 :: Float -> Float -> Float -> Vec3
+pattern WithVec3 a b c <- ((`withVec3` (,,)) -> (a, b, c))
+{-# COMPLETE WithVec3 #-}
+
+{-# INLINE fromVec2 #-}
+fromVec2 :: Coercible Vec3 a => Vec2 -> Float -> a
+fromVec2 xy z =
+  withVec2 xy \x y ->
+    coerce (vec3 x y z)
+
+{-# INLINE fromTuple #-}
+fromTuple :: Coercible Vec3 a => (Float, Float, Float) -> a
+fromTuple (x, y, z) = coerce (vec3 x y z)
+
+instance NFData Vec3 where
+  rnf Vec3{} = ()
+
+instance Num Vec3 where
+  {-# INLINE (+) #-}
+  Vec3 a b c + Vec3 d e f =
+    Vec3
+      (a + d)
+      (b + e)
+      (c + f)
+
+  {-# INLINE (-) #-}
+  Vec3 a b c - Vec3 d e f =
+    Vec3
+      (a - d)
+      (b - e)
+      (c - f)
+
+  {-# INLINE (*) #-}
+  Vec3 a b c * Vec3 d e f =
+    Vec3
+      (a * d)
+      (b * e)
+      (c * f)
+
+  {-# INLINE abs #-}
+  abs (Vec3 a b c) =
+    Vec3 (abs a) (abs b) (abs c)
+
+  {-# INLINE signum #-}
+  signum (Vec3 a b c) =
+    Vec3 (signum a) (signum b) (signum c)
+
+  {-# INLINE fromInteger #-}
+  fromInteger x = Vec3 x' x' x'
+    where
+      x' = fromInteger x
+
+instance Fractional Vec3 where
+  {-# INLINE (/) #-}
+  Vec3 l1 l2 l3 / Vec3 r1 r2 r3 =
+    Vec3 (l1 / r1) (l2 / r2) (l3 / r3)
+
+  {-# INLINE recip #-}
+  recip (Vec3 a b c) =
+    Vec3 (recip a) (recip b) (recip c)
+
+  {-# INLINE fromRational #-}
+  fromRational x = Vec3 x' x' x'
+    where
+      x' = fromRational x
+
+{-
+  XXX: GPU layouts call for some padding.
+
+  Maybe it would be worth it to flip the sizeOf-s.
+-}
+instance Storable Vec3 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 16
+
+  {-# INLINE alignment #-}
+  alignment _ = 4
+
+  {-# INLINE poke #-}
+  poke ptr v3 =
+    withVec3 v3 \a b c -> do
+      poke ptr' a
+      pokeElemOff ptr' 1 b
+      pokeElemOff ptr' 2 c
+      pokeElemOff ptr' 3 (1.0 :: Float)
+    where
+      ptr' = castPtr ptr
+
+  {-# INLINE peek #-}
+  peek ptr =
+    vec3 <$> peek ptr' <*> peekElemOff ptr' 1 <*> peekElemOff ptr' 2
+    where
+      ptr' = castPtr ptr
+
+{-# INLINE (^*) #-}
+(^*) :: Vec3 -> Float -> Vec3
+Vec3 a b c ^* x =
+  Vec3
+    (a * x)
+    (b * x)
+    (c * x)
+
+{-# INLINE (^/) #-}
+(^/) :: Vec3 -> Float -> Vec3
+Vec3 a b c ^/ x =
+  Vec3
+    (a / x)
+    (b / x)
+    (c / x)
+
+{-# INLINE lerp #-}
+lerp :: Float -> Vec3 -> Vec3 -> Vec3
+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)
+
+{-# INLINE cross #-}
+cross :: Vec3 -> Vec3 -> Vec3
+cross (Vec3 a b c) (Vec3 d e f) =
+  Vec3
+    (b * f - c * e)
+    (c * d - a * f)
+    (a * e - b * d)
+
+{-# INLINE dot #-}
+dot :: Vec3 -> Vec3 -> Float
+dot (Vec3 a b c) (Vec3 d e f) =
+  a * d +
+  b * e +
+  c * f
+
+{-# INLINE normalize #-}
+normalize :: Vec3 -> Vec3
+normalize v =
+  if nearZero q || nearZero (1-q) then
+    v
+  else
+    let
+      Vec3 x y z = v
+    in
+      Vec3 (x / l) (y / l) (z / l)
+
+  where
+    q = dot v v
+    l = sqrt q
+
+    nearZero a = abs a <= 1e-6
+
+-- * Unpadded
+
+newtype Packed = Packed { unPacked :: Vec3 }
+  deriving (Eq, Ord, Show, NFData, Num, Fractional)
+
+{-# INLINE packed #-}
+packed :: Float -> Float -> Float -> Packed
+packed x y z = Packed (vec3 x y z)
+
+instance Storable Packed where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 12
+
+  {-# INLINE alignment #-}
+  alignment _ = 4
+
+  {-# INLINE poke #-}
+  poke ptr (Packed v3) =
+    withVec3 v3 \a b c -> do
+      poke ptr' a
+      pokeElemOff ptr' 1 b
+      pokeElemOff ptr' 2 c
+    where
+      ptr' = castPtr ptr
+
+  {-# INLINE peek #-}
+  peek ptr = packed
+    <$> peek ptr'
+    <*> peekElemOff ptr' 1
+    <*> peekElemOff ptr' 2
+    where
+      ptr' = castPtr ptr
diff --git a/src/Geomancy/Vec4.hs b/src/Geomancy/Vec4.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Vec4.hs
@@ -0,0 +1,204 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- | Specialized and inlined @V4 Float@.
+
+module Geomancy.Vec4
+  ( Vec4
+  , vec4
+  , withVec4
+  , pattern WithVec4
+  , fromVec2
+  , fromVec22
+  , fromVec3
+  , fromTuple
+
+  , (^*)
+  , (^/)
+  , lerp
+
+  , dot
+  , normalize
+  ) where
+
+import Control.DeepSeq (NFData(rnf))
+import Data.Coerce (Coercible, coerce)
+import Foreign (Storable(..), castPtr)
+
+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)
+
+{-# INLINE vec4 #-}
+vec4 :: Float -> Float -> Float -> Float -> Vec4
+vec4 = Vec4
+
+{-# INLINE withVec4 #-}
+withVec4
+  :: Vec4
+  -> (Float -> Float -> Float -> Float -> r)
+  -> r
+withVec4 (Vec4 a b c d) f = f a b c d
+
+pattern WithVec4 :: Float -> Float -> Float -> Float -> Vec4
+pattern WithVec4 a b c d <- ((`withVec4` (,,,)) -> (a, b, c, d))
+{-# COMPLETE WithVec4 #-}
+
+{-# INLINE fromVec2 #-}
+fromVec2 :: Vec2 -> Float -> Float -> Vec4
+fromVec2 xy z w =
+  withVec2 xy \x y ->
+    vec4 x y z w
+
+{-# INLINE fromVec22 #-}
+fromVec22 :: Vec2 -> Vec2 -> Vec4
+fromVec22 xy zw =
+  withVec2 xy \x y ->
+  withVec2 zw \z w ->
+    vec4 x y z w
+
+{-# INLINE fromVec3 #-}
+fromVec3 :: Coercible a Vec3 => a -> Float -> Vec4
+fromVec3 xyz w =
+  withVec3 (coerce xyz) \x y z ->
+    vec4 x y z w
+
+{-# INLINE fromTuple #-}
+fromTuple :: (Float, Float, Float, Float) -> Vec4
+fromTuple (x, y, z, w) = vec4 x y z w
+
+instance NFData Vec4 where
+  rnf Vec4{} = ()
+
+instance Num Vec4 where
+  {-# INLINE (+) #-}
+  Vec4 l1 l2 l3 l4 + Vec4 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)
+
+  {-# INLINE (*) #-}
+  Vec4 l1 l2 l3 l4 * Vec4 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)
+
+  {-# INLINE signum #-}
+  signum (Vec4 a b c d) =
+    Vec4 (signum a) (signum b) (signum c) (signum d)
+
+  {-# INLINE fromInteger #-}
+  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)
+
+  {-# INLINE recip #-}
+  recip (Vec4 a b c d) =
+    Vec4 (recip a) (recip b) (recip c) (recip d)
+
+  {-# INLINE fromRational #-}
+  fromRational x = Vec4 x' x' x' x'
+    where
+      x' = fromRational x
+
+instance Storable Vec4 where
+  {-# INLINE sizeOf #-}
+  sizeOf _ = 16
+
+  {-# INLINE alignment #-}
+  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
+
+  {-# INLINE peek #-}
+  peek ptr = vec4
+    <$> peek ptr'
+    <*> peekElemOff ptr' 1
+    <*> peekElemOff ptr' 2
+    <*> peekElemOff ptr' 3
+    where
+      ptr' = castPtr ptr
+
+{-# INLINE (^*) #-}
+(^*) :: Vec4 -> Float -> Vec4
+Vec4 a b c d ^* x =
+  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)
+
+{-# INLINE lerp #-}
+lerp :: Float -> Vec4 -> Vec4 -> Vec4
+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)
+
+{-# 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
+
+{-# INLINE normalize #-}
+normalize :: Vec4 -> Vec4
+normalize v =
+  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)
+
+  where
+    q = dot v v
+    l = sqrt q
+
+    nearZero a = abs a <= 1e-6
diff --git a/src/Geomancy/Vulkan/Projection.hs b/src/Geomancy/Vulkan/Projection.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Vulkan/Projection.hs
@@ -0,0 +1,72 @@
+module Geomancy.Vulkan.Projection
+  ( perspective
+  , infinitePerspective
+  , orthoOffCenter
+  ) where
+
+import Geomancy.Mat4 (colMajor)
+import Geomancy.Transform (Transform(..))
+
+perspective
+  :: Integral side
+  => Float
+  -> Float -> Float
+  -> side -> side
+  -> Transform
+perspective fovRads near far width height = colMajor
+  x 0   0   0
+  0 y   0   0
+  0 0   z w23
+  0 0 w32   1
+
+  where
+    x = cotFoV * aspectX
+    y = -cotFoV
+    z = far / (near - far)
+    w23 = near * far / (near - far)
+    w32 = -1
+
+    cotFoV = recip . tan $ 0.5 * fovRads
+
+    aspectX = fromIntegral height / fromIntegral width
+
+infinitePerspective
+  :: Integral side
+  => Float
+  -> side
+  -> side
+  -> Transform
+infinitePerspective fovRads width height = colMajor
+  x   0   0  0
+  0 (-y)  0  0
+  0   0 (-1) w
+  0   0 (-1) 0
+  where
+    (x, y) =
+      if width > height then
+        ( cotFoV / camAspect
+        , cotFoV
+        )
+      else
+        ( cotFoV
+        , cotFoV * camAspect
+        )
+    camAspect = fromIntegral width / fromIntegral height
+    cotFoV = recip . tan $ 0.5 * fovRads
+
+    w = -2 * near
+    near = 1/128 -- 2048
+
+orthoOffCenter :: Integral side => Float -> Float -> side -> side -> Transform
+orthoOffCenter near far width height = colMajor
+  x 0 0 0
+  0 y 0 0
+  0 0 z w
+  0 0 0 1
+
+  where
+    x = 2 / fromIntegral width
+    y = 2 / fromIntegral height
+    z = 1 / (far - near)
+
+    w = near * (near - far)
diff --git a/src/Geomancy/Vulkan/View.hs b/src/Geomancy/Vulkan/View.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Vulkan/View.hs
@@ -0,0 +1,41 @@
+{-# LANGUAGE BlockArguments #-}
+
+module Geomancy.Vulkan.View
+  ( orthoFitScreen
+  , lookAt
+  ) where
+
+import Geomancy.Mat4 (Mat4, rowMajor)
+import Geomancy.Vec3 (Vec3, withVec3)
+import Geomancy.Transform (Transform(..))
+
+import qualified Geomancy.Vec3 as Vec3
+
+lookAt :: Vec3 -> Vec3 -> Vec3 -> Transform
+lookAt eye center up =
+  withVec3 xa \xaX xaY xaZ ->
+  withVec3 ya \yaX yaY yaZ ->
+  withVec3 za \zaX zaY zaZ ->
+  rowMajor
+    xaX yaX (-zaX) 0
+    xaY yaY (-zaY) 0
+    xaZ yaZ (-zaZ) 0
+    xd  yd    zd   1
+  where
+    xa = Vec3.normalize $ Vec3.cross za up
+    ya = Vec3.cross xa za
+    za = Vec3.normalize $ center - eye
+
+    xd = - Vec3.dot xa eye
+    yd = - Vec3.dot ya eye
+    zd =   Vec3.dot za eye
+
+orthoFitScreen :: Float -> Float -> Float -> Float -> Mat4
+orthoFitScreen screenWidth screenHeight targetWidth targetHeight =
+  rowMajor
+    s 0 0 0
+    0 s 0 0
+    0 0 1 0
+    0 0 0 1
+  where
+    s = min (screenWidth / targetWidth) (screenHeight / targetHeight)
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,151 @@
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE NumericUnderscores #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE TemplateHaskell #-}
+
+import Hedgehog
+import qualified Hedgehog.Gen as Gen
+import qualified Hedgehog.Range as Range
+
+import Control.Monad (unless)
+import Control.Monad.IO.Class (MonadIO(..))
+import Data.Foldable (toList)
+import Data.Maybe (catMaybes)
+import GHC.Stack (withFrozenCallStack)
+import System.Exit (exitFailure, exitSuccess)
+import Text.Printf (printf)
+
+import Geomancy.Transform (Transform(..))
+
+import qualified Foreign
+import qualified Linear
+import qualified Geomancy
+import qualified Geomancy.Mat4
+import qualified Geomancy.Transform as Transform
+
+import Linear ((!*!))
+
+main :: IO ()
+main = do
+  passed <- checkParallel discovered
+  if passed then
+    exitSuccess
+  else
+    exitFailure
+
+-- | Enough tests to cover the principles
+pattern PROP_TESTS :: TestLimit
+pattern PROP_TESTS = 10_000
+
+-- | Try harder to catch FP precision errors
+pattern PROP_TESTS_BRUTAL :: TestLimit
+pattern PROP_TESTS_BRUTAL = 10_000_000
+
+prop_assoc_multiply :: Property
+prop_assoc_multiply = withTests PROP_TESTS $ property do
+  (a, a_) <- forAllTransform
+  (b, b_) <- forAllTransform
+  (c, c_) <- forAllTransform
+
+  let
+    ab'c = (a <> b) <> c
+    a'bc = a <> (b <> c)
+    delta' = nearlyEqualMat4 ab'c a'bc
+  annotateShow delta'
+
+  let
+    ab_c = (a_ !*! b_) !*! c_
+    a_bc = a_ !*! (b_ !*! c_)
+    delta_ = nearlyEqualM44 ab_c a_bc
+  annotateShow delta_
+
+  -- Intra-library transitivity
+  unless (null $ catMaybes delta') do
+    -- XXX: check only if there is some outstanding error
+    delta' === delta_
+
+  -- Inter-library calculated values nearlyEqual
+  ab'c_ <- toM44 ab'c
+  catMaybes (nearlyEqualM44 ab_c ab'c_) === []
+
+  a'bc_ <- toM44 a'bc
+  catMaybes (nearlyEqualM44 a_bc a'bc_) === []
+
+forAllTransform :: PropertyT IO (Geomancy.Mat4.Mat4, Linear.M44 Float)
+forAllTransform = withFrozenCallStack do
+  (_name, Transform g) <- forAllWith fst genTransform
+  l <- toM44 g
+  pure (g, l)
+
+genTransform :: Gen ([Char], Transform)
+genTransform = Gen.choice
+  [ genIdentity
+  , genTranslate
+  , genRotate
+  , genScale
+  ]
+  where
+    genIdentity = pure ("identity", mempty)
+
+    genTranslate = do
+      x <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)
+      y <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)
+      z <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)
+      pure
+        ( printf "translate %0.4f %0.4f %0.4f" x y z
+        , Transform.translate x y z
+        )
+
+    genRotate = do
+      (name, axis) <- Gen.element
+        [ ("rotate/x", Transform.rotateX)
+        , ("rotate/y", Transform.rotateY)
+        , ("rotate/z", Transform.rotateZ)
+        ]
+      angle <- Gen.float (Range.linearFracFrom 0.0 (-4 * pi) (4 * pi))
+      pure
+        ( printf "%s %0.4f" name angle
+        , axis angle
+        )
+
+    genScale = do
+      x <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)
+      y <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)
+      z <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)
+      pure
+        ( printf "scale %0.4f %0.4f %0.4f" x y z
+        , Transform.scale3 x y z
+        )
+
+-- toVulkan :: Linear.M44 Float -> Linear.M44 Float
+-- toVulkan = {- Linear.transpose . -} (correction Linear.!*!)
+--   where
+--     correction = Linear.V4
+--       (Linear.V4 1   0  0   0)
+--       (Linear.V4 0 (-1) 0   0)
+--       (Linear.V4 0   0  0.5 0.5)
+--       (Linear.V4 0   0  0   1)
+
+toM44 :: MonadIO io => Geomancy.Mat4 -> io (Linear.M44 Float)
+toM44 mat4 =
+  liftIO $
+    Foreign.with mat4 $
+      Foreign.peek . Foreign.castPtr
+
+type NearlyEqual = Maybe (Float, Float, Float)
+
+nearlyEqualMat4 :: Geomancy.Mat4.Mat4 -> Geomancy.Mat4.Mat4 -> [NearlyEqual]
+nearlyEqualMat4 a b = Geomancy.Mat4.zipWith nearlyEqual a b
+
+nearlyEqualM44 :: Linear.M44 Float -> Linear.M44 Float -> [NearlyEqual]
+nearlyEqualM44 a b = zipWith nearlyEqual (concatMap toList a) (concatMap toList b)
+
+nearlyEqual :: Float -> Float -> NearlyEqual
+nearlyEqual x y =
+  if x == y || abs (1 - x / y) < 0.001 then
+    Nothing
+  else
+    Just (x, y, abs (1 - x / y))
+
+discovered :: Group
+discovered = $$(discover)