diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,5 +1,30 @@
 # Changelog for geomancy
 
+## 0.3.0.0
+
+> `[0]` releases without row/col-major, handedness, and counter/clockwise confusion.
+
+This release brings sanity and consistency around the transform stack.
+- `Transform`s have proper *col-major* definitions.
+- Operation ordering follows the GLSL convention.
+  * The composition order is `(p <> v <> m <> ...)` (from global to local).
+    Reverse your transforms to adjust.
+  * The application order is `m !* v` (with the operators becoming infixr 5, just under the `<>`).
+    It will do `Mᵀ * v` inside for SIMD reasons and to match what GLSL does for its `M * v` ops.
+- `Geomancy.Vulkan.Projection` is right-handed, BUT produces *reverse-depth* range ([1; 0], In 3D infinite-Z converges to 0).
+  * Replace your depth tests with `OP_GREATER` and clear to `0.0` - get better precision for your migration troubles.
+- `Geomancy.Vulkan.View` is *right-handed*, with +Z being forward.
+  * The intended up vector is still `vec3 0 (-1) 0` -- +Y down.
+    Silly as it sounds, this matches the XY plane of the window with XY plane in front of a "first person" camera.
+- Axis rotations (using `rotateQ`) will appear clockwise when looking along the axis.
+- Angle rotations follow Tait-Bryan angles (heading/elevation/bank or yaw/pitch/roll) in the y-x-z frame.
+  * `rotateZ (time * rate)` will follow the clock hands in 2D scenes and roll in 3D.
+  * `rotateX` will follow the sun from sunrise to sunset, pitching UP / increasing elevation.
+  * `rotateY` will turn you right, increasing yaw / heading eastwards.
+- You're of course free to define your own transforms, just copy the modules and tune to your liking.
+  Just make sure that you use matching row/column constructors and the math layer will do the rest, fast.
+- Added `webcolor-labels` instances for UVec3/Vec3/Vec4.
+
 ## 0.2.6.0
 
 * `Geomancy.Gl.Block` extract to `gl-block` package as `Graphics.Gl.Block`.
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -6,12 +6,6 @@
 * `Mat4` and `Vec4` are `ByteArray#`.
 * `Mat4`x`Mat4` and `Mat4`x`Vec4` is done with SIMD.
 
-## 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):
@@ -64,3 +58,45 @@
 Add that time to the poking that'll follow.
 
 Sure, it is in the lower microseconds range, but this budget can be used elsewhere.
+
+## Conventions
+
+### Matrix layout
+
+Transforms produced, composed, and applied to mimic the GLSL order (col-major):
+
+- `vec4 vPosOut = P * V * M * vPosIn;`
+- `vPosOut = (p <> v <> m) !* vPosIn`
+
+This way you don't have to transpose your transforms or fiddle with layout annotations.
+
+### Projections / Views
+
+`Geomancy.Vulkan.Projection` is using the "reverse-depth" trick that remaps the vulkan default `[0; 1]` range to `[1; 0]`.
+This grants extra precision with one less parameter to specify (you only need "near" now), but makes handedness reasoning tricky.
+The default depth range the coordinate is left-handed (+X right, +Y down, +Z forward).
+But after reversing the depth it has to be paired with a right-handed view function like `Geomancy.Vulkan.View.lookAtRH`.
+
+The intended up vector is still `vec3 0 (-1) 0` -- +Y down.
+Silly as it sounds, this matches the XY plane of the window with XY plane in front of a "first person" camera.
+
+### Rotations
+
+Axis rotations (using `rotateQ`) will appear clockwise when looking along the axis.
+
+Angle rotations follow Tait-Bryan angles (heading/elevation/bank or yaw/pitch/roll) in the y-x-z frame.
+- `rotateZ (time * rate)` will follow the clock hands in 2D scenes and roll in 3D.
+- `rotateX` will follow the sun from sunrise to sunset, increasing elevation / pitching UP.
+- `rotateY` will turn you right, increasing yaw/heading eastwards.
+
+Using `Geomancy.Quaternion.intrinsic roll pitch yaw` will make a rotation from the 3 angles in one go.
+You can use it to `rotate` a point directly (e.g `vec3 0 0 1` to get a direction vector from Quaternion) or commit to a matrix using `Transform.rotateQ`.
+
+![yaw-pitch-roll](./yaw-pitch-roll.svg)
+
+You're of course free to define your own transforms, just copy the modules and tune to your liking.
+Just make sure that you use matching row/column constructors and the math layer will do the rest, fast.
+
+## GLSL-like functions
+
+To further facilitate conversion between the host and shader code `Geomancy.Gl.Funs` provides common functions like `glFract` and `smoothstep`.
diff --git a/bench/Main.hs b/bench/Main.hs
--- a/bench/Main.hs
+++ b/bench/Main.hs
@@ -10,6 +10,7 @@
 
 import qualified Geomancy
 import qualified Geomancy.Mat4
+import qualified Geomancy.Quaternion
 
 import qualified Linear
 import qualified Linear.Matrix
@@ -26,6 +27,10 @@
   , bgroup "4x4 poke" $ flip map [10, 100, 1000, 10000] \size ->
       env (clones size mempty Linear.identity) $
         mat4poke size
+
+  , bgroup "quaternion"
+      [ bgroup "euler" quatFromEuler
+      ]
   ]
 
 mat4transpose :: [Benchmark]
@@ -73,3 +78,9 @@
     pokeLinearTranspose =
       Foreign.withArray (map Linear.transpose m44s) \_ptr ->
         pure ()
+
+quatFromEuler :: [Benchmark]
+quatFromEuler =
+  [ bench "fused" $ whnf (\a -> Geomancy.Quaternion.intrinsic a a a) (pi/3)
+  , bench "composed" $ whnf (\a -> Geomancy.Quaternion.extrinsic a a a) (pi/3)
+  ]
diff --git a/cbits/mat4.c b/cbits/mat4.c
--- a/cbits/mat4.c
+++ b/cbits/mat4.c
@@ -3,15 +3,15 @@
 #include <xmmintrin.h>
 
 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]);
-  __m128 row3 = _mm_load_ps(&B[12]);
+  __m128 row0 = _mm_load_ps(&A[0]);
+  __m128 row1 = _mm_load_ps(&A[4]);
+  __m128 row2 = _mm_load_ps(&A[8]);
+  __m128 row3 = _mm_load_ps(&A[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 brod0 = _mm_set1_ps(B[4*i + 0]);
+    __m128 brod1 = _mm_set1_ps(B[4*i + 1]);
+    __m128 brod2 = _mm_set1_ps(B[4*i + 2]);
+    __m128 brod3 = _mm_set1_ps(B[4*i + 3]);
     __m128 row = _mm_add_ps(
       _mm_add_ps(
         _mm_mul_ps(brod0, row0),
@@ -86,28 +86,28 @@
     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);
+    C0 = vfmaq_laneq_f32(C0, A0, B0, 0);
+    C0 = vfmaq_laneq_f32(C0, A1, B0, 1);
+    C0 = vfmaq_laneq_f32(C0, A2, B0, 2);
+    C0 = vfmaq_laneq_f32(C0, A3, B0, 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);
+    C1 = vfmaq_laneq_f32(C1, A0, B1, 0);
+    C1 = vfmaq_laneq_f32(C1, A1, B1, 1);
+    C1 = vfmaq_laneq_f32(C1, A2, B1, 2);
+    C1 = vfmaq_laneq_f32(C1, A3, B1, 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);
+    C2 = vfmaq_laneq_f32(C2, A0, B2, 0);
+    C2 = vfmaq_laneq_f32(C2, A1, B2, 1);
+    C2 = vfmaq_laneq_f32(C2, A2, B2, 2);
+    C2 = vfmaq_laneq_f32(C2, A3, B2, 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);
+    C3 = vfmaq_laneq_f32(C3, A0, B3, 0);
+    C3 = vfmaq_laneq_f32(C3, A1, B3, 1);
+    C3 = vfmaq_laneq_f32(C3, A2, B3, 2);
+    C3 = vfmaq_laneq_f32(C3, A3, B3, 3);
     vst1q_f32(&O[12], C3);
 
 }
diff --git a/geomancy.cabal b/geomancy.cabal
--- a/geomancy.cabal
+++ b/geomancy.cabal
@@ -1,18 +1,18 @@
 cabal-version: 1.12
 
--- This file has been generated from package.yaml by hpack version 0.35.2.
+-- This file has been generated from package.yaml by hpack version 0.38.1.
 --
 -- see: https://github.com/sol/hpack
 
 name:           geomancy
-version:        0.2.6.0
-synopsis:       Geometry and matrix manipulation
+version:        0.3.0.0
+synopsis:       Vectors 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
+author:         IC Rainbow
 maintainer:     aenor.realm@gmail.com
-copyright:      2021 Alexander Bondarenko
+copyright:      2021-2025 IC Rainbow
 license:        BSD3
 license-file:   LICENSE
 build-type:     Simple
@@ -34,6 +34,7 @@
       Geomancy.IVec3
       Geomancy.IVec4
       Geomancy.Mat4
+      Geomancy.Mat4.Internal
       Geomancy.Point
       Geomancy.Quaternion
       Geomancy.Swizzle
@@ -64,6 +65,7 @@
     , mono-traversable
     , ptrdiff
     , simple-affine-space
+    , webcolor-labels
   default-language: Haskell2010
 
 test-suite geomancy-test
diff --git a/src/Geomancy/Gl/Funs.hs b/src/Geomancy/Gl/Funs.hs
--- a/src/Geomancy/Gl/Funs.hs
+++ b/src/Geomancy/Gl/Funs.hs
@@ -6,6 +6,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
 
 module Geomancy.Gl.Funs
   ( GlClamp(..)
diff --git a/src/Geomancy/Interpolate.hs b/src/Geomancy/Interpolate.hs
--- a/src/Geomancy/Interpolate.hs
+++ b/src/Geomancy/Interpolate.hs
@@ -1,8 +1,9 @@
 {-# LANGUAGE BlockArguments #-}
-{-# LANGUAGE RankNTypes #-}
-{-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE PolyKinds #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
 
 module Geomancy.Interpolate
   ( linear
diff --git a/src/Geomancy/Mat4.hs b/src/Geomancy/Mat4.hs
--- a/src/Geomancy/Mat4.hs
+++ b/src/Geomancy/Mat4.hs
@@ -1,58 +1,65 @@
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE BlockArguments #-}
-{-# LANGUAGE DataKinds #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MagicHash #-}
 {-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE UnboxedTuples #-}
 {-# LANGUAGE UnliftedFFITypes #-}
 {-# LANGUAGE ViewPatterns #-}
 
+{-# OPTIONS_GHC -Wno-orphans #-} -- XXX: until there are order-specific types
+
 -- | General matrix storage and operations.
 
 module Geomancy.Mat4
   ( Mat4
-
-  , rowMajor
-  , withRowMajor
-  , toListRowMajor
-  , toListRowMajor2d
-  , fromRowMajor2d
-
-  , colMajor
-  , withColMajor
-  , toListColMajor
-  , toListColMajor2d
-
+    -- * Operations
   , identity
   , transpose
   , inverse
+  , det
   , pointwise
   , zipWith
   , matrixProduct
   , scalarMultiply
   , (!*)
+
+    -- * Construction
+
+    -- ** Col-major
+  , colMajor
+  , withColMajor
+  , toListColMajor
+  , toListColMajor2d
+  , showColumns
+
+    -- ** Row-major
+  , rowMajor
+  , withRowMajor
+  , toListRowMajor
+  , toListRowMajor2d
+  , fromRowMajor2d
+  , showRows
   ) where
 
 import Prelude hiding (zipWith)
 import GHC.Exts hiding (VecCount(..), toList)
+import Geomancy.Mat4.Internal
 
-import Control.DeepSeq (NFData(rnf))
-import Foreign (Storable(..))
-import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)
-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 Graphics.Gl.Block (Block(..))
-import Geomancy.Vec4 (Vec4(..), unsafeNewVec4)
+import Geomancy.Vec4 (Vec4(..), newVec4)
 
-data Mat4 = Mat4 ByteArray#
+toListColMajor :: Coercible a Mat4 => a -> [Float]
+toListColMajor = toListTrans . coerce
 
+toListColMajor2d :: Coercible a Mat4 => a -> [[Float]]
+toListColMajor2d = toList2dTrans . coerce
+
 {- | Construct 'Mat4' from @row@ notation.
 -}
 rowMajor
@@ -62,28 +69,17 @@
   -> 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
+rowMajor
+  e0 e1 e2 e3
+  e4 e5 e6 e7
+  e8 e9 eA eB
+  eC eD eE eF =
+    -- just store
+    coerce $ fromMemory
+      e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF
 
 {- |
   Build a Mat4 from a list-of-lists kind of container
@@ -103,15 +99,15 @@
 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
+      withRow r0 \e0 e1 e2 e3 ->
+      withRow r1 \e4 e5 e6 e7 ->
+      withRow r2 \e8 e9 eA eB ->
+      withRow r3 \eC eD eE eF ->
+        Just . coerce @a $ rowMajor
+          e0 e1 e2 e3
+          e4 e5 e6 e7
+          e8 e9 eA eB
+          eC eD eE eF
     _ ->
       Nothing
   where
@@ -122,31 +118,10 @@
         _ ->
           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.
+{- | Reduce 'Mat4' with a function with @row@ notation for arguments.
 -}
-{-# INLINE withColMajor #-}
-withColMajor
+{-# INLINE withRowMajor #-}
+withRowMajor
   :: Coercible a Mat4
   => a
   ->
@@ -157,70 +132,44 @@
       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
+withRowMajor m f = withMemory (coerce m)
+  \ e0 e1 e2 e3
+    e4 e5 e6 e7
+    e8 e9 eA eB
+    eC eD eE eF ->
+  f -- preserved element order
+    e0 e1 e2 e3
+    e4 e5 e6 e7
+    e8 e9 eA eB
+    eC eD eE eF
 
-{- | Construct 'Mat4' from elements in memory order.
+{- | Construct a 'Mat4' from @column@ notation.
 -}
-{-# INLINE mat4 #-}
-mat4
-  :: Float -> Float -> Float -> Float
+{-# INLINE colMajor #-}
+colMajor
+  :: forall a . Coercible Mat4 a
+  => 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'
+  -> a
+colMajor
+  e0 e4 e8 eC
+  e1 e5 e9 eD
+  e2 e6 eA eE
+  e3 e7 eB eF =
+    coerce $ fromMemory
+      e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF
 
-{- | Reduce 'Mat4' with a function with @memory@ notation for arguments.
+{- | Reduce 'Mat4' with a function with @column@ notation for arguments.
 -}
-{-# INLINE withMat4 #-}
-withMat4
-  :: Mat4
+{-# INLINE withColMajor #-}
+withColMajor
+  :: Coercible a Mat4
+  => a
   ->
     ( Float -> Float -> Float -> Float ->
       Float -> Float -> Float -> Float ->
@@ -229,68 +178,64 @@
       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#))
+withColMajor m f = withMemory (coerce m)
+  \ e0 e4 e8 eC
+    e1 e5 e9 eD
+    e2 e6 eA eE
+    e3 e7 eB eF ->
+  f -- transposed element order
+    e0 e1 e2 e3
+    e4 e5 e6 e7
+    e8 e9 eA eB
+    eC eD eE eF
 
-{- | @I@, the identity matrix.
+toListRowMajor :: Coercible a Mat4 => a -> [Float]
+toListRowMajor = toListMemory . coerce
 
-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
+toListRowMajor2d :: Coercible a Mat4 => a -> [[Float]]
+toListRowMajor2d = toList2dMemory . coerce
 
-{-# INLINE transpose #-}
-transpose :: Mat4 -> Mat4
-transpose =
-  flip withMat4
+det :: forall a . (Coercible Mat4 a, Coercible Mat4 a) => a -> Float
+det m =
+  withRowMajor (coerce @_ @a m)
     \ 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
+        let
+          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
+        in
+            s0 * c5
+          - s1 * c4
+          + s2 * c3
+          + s3 * c2
+          - s4 * c1
+          + s5 * c0
+
 {- | Compute an inverse matrix, slowly.
 -}
-inverse :: (Coercible Mat4 a, Coercible Mat4 a) => a -> a
+inverse :: forall a . (Coercible Mat4 a, Coercible Mat4 a) => a -> a
 inverse m =
-  coerce $ withMat4 (coerce m)
+  coerce @a $ withRowMajor (coerce @_ @a m)
     \ m00 m01 m02 m03
       m10 m11 m12 m13
       m20 m21 m22 m23
       m30 m31 m32 m33 ->
         let
-          invDet = recip det
-
-          det
-            = s0 * c5
+          invDet = recip $
+              s0 * c5
             - s1 * c4
             + s2 * c3
             + s3 * c2
@@ -331,143 +276,36 @@
           i32 = (-m30 * s3 + m31 * s1 - m32 * s0) * invDet
           i33 = ( m20 * s3 - m21 * s1 + m22 * s0) * invDet
         in
-          mat4
+          rowMajor
             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)
-
-foreign import ccall unsafe "Mat4xMat4_SIMD" m4m4simd :: Addr# -> Addr# -> Addr# -> IO ()
+foreign import ccall unsafe "Mat4xMat4_SIMD" m4m4simd :: ByteArray# -> ByteArray# -> ByteArray# -> IO ()
 
 {-# INLINE matrixProduct #-}
 matrixProduct :: Mat4 -> Mat4 -> Mat4
 matrixProduct (Mat4 l) (Mat4 r) = unsafePerformIO do
-  result@(Mat4 m) <- unsafeNewMat4
-  m4m4simd
-    (byteArrayContents# l)
-    (byteArrayContents# r)
-    (byteArrayContents# m)
+  result@(Mat4 m) <- newMat4
+  m4m4simd l r 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 #)
-
-{-# 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_SIMD" m4v4simd :: ByteArray# -> ByteArray# -> ByteArray# -> IO ()
 
-foreign import ccall unsafe "Mat4xVec4_SIMD" m4v4simd :: Addr# -> Addr# -> Addr# -> IO ()
+{- | Matrix - column vector multiplication (post)
 
--- | Matrix - column vector multiplication
+@
+vOut = p <> v <> m !* vIn
+@
+-}
 (!*) :: Coercible a Mat4 => a -> Vec4 -> Vec4
 (!*) (coerce -> Mat4 m) (Vec4 v) = unsafePerformIO do
-  result@(Vec4 o) <- unsafeNewVec4
-  m4v4simd
-    (byteArrayContents# m)
-    (byteArrayContents# v)
-    (byteArrayContents# o)
+  result@(Vec4 o) <- newVec4
+  !() <- m4v4simd m v 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{} = ()
+infixr 5 !*
 
 instance Semigroup Mat4 where
   {-# INLINE (<>) #-}
@@ -478,60 +316,30 @@
   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', () #)
+  show = showColumns
 
-  {-# 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' #)
+showColumns :: Mat4 -> String
+showColumns cm = withColMajor cm
+  \ e0 e4 e8 eC
+    e1 e5 e9 eD
+    e2 e6 eA eE
+    e3 e7 eB eF ->
+  unlines
+    [ printf  "/ %.4f %.4f %.4f %.4f \\" e0 e4 e8 eC
+    , printf  "| %.4f %.4f %.4f %.4f |"  e1 e5 e9 eD
+    , printf  "| %.4f %.4f %.4f %.4f |"  e2 e6 eA eE
+    , printf "\\ %.4f %.4f %.4f %.4f /"  e3 e7 eB eF
+    ]
 
-instance Block Mat4 where
-  type PackedSize Mat4 = 64
-  alignment140 _  = 16
-  sizeOf140       = sizeOfPacked
-  alignment430    = alignment140
-  sizeOf430       = sizeOf140
-  isStruct _      = False
-  read140     = peekDiffOff
-  write140    = pokeDiffOff
-  read430     = read140
-  write430    = write140
-  readPacked  = read140
-  writePacked = write140
-  {-# INLINE alignment140 #-}
-  {-# INLINE sizeOf140 #-}
-  {-# INLINE alignment430 #-}
-  {-# INLINE sizeOf430 #-}
-  {-# INLINE isStruct #-}
-  {-# INLINE read140 #-}
-  {-# INLINE write140 #-}
-  {-# INLINE read430 #-}
-  {-# INLINE write430 #-}
-  {-# INLINE readPacked #-}
-  {-# INLINE writePacked #-}
+showRows :: Mat4 -> String
+showRows rm = withRowMajor rm
+  \ e0 e1 e2 e3
+    e4 e5 e6 e7
+    e8 e9 eA eB
+    eC eD eE eF ->
+  unlines
+    [ printf "[ %.4f %.4f %.4f %.4f |" e0 e1 e2 e3
+    , printf "| %.4f %.4f %.4f %.4f |" e4 e5 e6 e7
+    , printf "| %.4f %.4f %.4f %.4f |" e8 e9 eA eB
+    , printf "| %.4f %.4f %.4f %.4f ]" eC eD eE eF
+    ]
diff --git a/src/Geomancy/Mat4/Internal.hs b/src/Geomancy/Mat4/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Geomancy/Mat4/Internal.hs
@@ -0,0 +1,286 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE MagicHash #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE UnboxedTuples #-}
+
+-- | General matrix storage and operations.
+
+module Geomancy.Mat4.Internal
+  ( Mat4(..)
+  , fromMemory
+  , withMemory
+  , newMat4
+    -- * Order-independent operations
+  , identity
+  , transpose
+  , pointwise
+  , scalarMultiply
+    -- * List operations in memory order
+  , toListMemory
+  , toList2dMemory
+  , toListTrans
+  , toList2dTrans
+  , zipWith
+  ) where
+
+import Prelude hiding (zipWith)
+import GHC.Exts hiding (VecCount(..), toList)
+
+import Control.DeepSeq (NFData(rnf))
+import Foreign (Storable(..))
+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)
+import GHC.IO (IO(..))
+
+import qualified Data.List as List
+
+import Graphics.Gl.Block (Block(..))
+
+data Mat4 = Mat4 ByteArray#
+
+instance NFData Mat4 where
+  rnf Mat4{} = ()
+
+{- | Construct 'Mat4' from elements in memory order.
+-}
+{-# INLINE fromMemory #-}
+fromMemory
+  :: Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Float -> Float -> Float -> Float
+  -> Mat4
+fromMemory
+  (F# e0) (F# e1) (F# e2) (F# e3)
+  (F# e4) (F# e5) (F# e6) (F# e7)
+  (F# e8) (F# e9) (F# ea) (F# eb)
+  (F# ec) (F# ed) (F# ee) (F# ef) =
+  runRW# \world ->
+    let
+      !(# world_, arr #) = newAlignedPinnedByteArray# 64# 16# world
+
+      world0 = writeFloatArray# arr 0x0# e0 world_
+      world1 = writeFloatArray# arr 0x1# e1 world0
+      world2 = writeFloatArray# arr 0x2# e2 world1
+      world3 = writeFloatArray# arr 0x3# e3 world2
+
+      world4 = writeFloatArray# arr 0x4# e4 world3
+      world5 = writeFloatArray# arr 0x5# e5 world4
+      world6 = writeFloatArray# arr 0x6# e6 world5
+      world7 = writeFloatArray# arr 0x7# e7 world6
+
+      world8 = writeFloatArray# arr 0x8# e8 world7
+      world9 = writeFloatArray# arr 0x9# e9 world8
+      worldA = writeFloatArray# arr 0xA# ea world9
+      worldB = writeFloatArray# arr 0xB# eb worldA
+
+      worldC = writeFloatArray# arr 0xC# ec worldB
+      worldD = writeFloatArray# arr 0xD# ed worldC
+      worldE = writeFloatArray# arr 0xE# ee worldD
+      worldF = writeFloatArray# arr 0xF# ef worldE
+
+      !(# _world', arr' #) = unsafeFreezeByteArray# arr worldF
+    in
+      Mat4 arr'
+
+{- | Reduce 'Mat4' with a function with @memory@ notation for arguments.
+-}
+{-# INLINE withMemory #-}
+withMemory
+  :: Mat4
+  ->
+    ( Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      Float -> Float -> Float -> Float ->
+      r
+    )
+  -> r
+withMemory (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 = fromMemory
+  1 0 0 0
+  0 1 0 0
+  0 0 1 0
+  0 0 0 1
+
+-- TODO: simdify
+{-# INLINE transpose #-}
+transpose :: Mat4 -> Mat4
+transpose =
+  flip withMemory
+    \ e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF ->
+    fromMemory
+      e0 e4 e8 eC
+      e1 e5 e9 eD
+      e2 e6 eA eE
+      e3 e7 eB eF
+
+pointwise :: Mat4 -> Mat4 -> (Float -> Float -> Float) -> Mat4
+pointwise a b f =
+  withMemory a
+    \ a0 a1 a2 a3
+      a4 a5 a6 a7
+      a8 a9 aA aB
+      aC aD aE aF ->
+  withMemory b
+    \ b0 b1 b2 b3
+      b4 b5 b6 b7
+      b8 b9 bA bB
+      bC bD bE bF ->
+  fromMemory
+    (f a0 b0) (f a1 b1) (f a2 b2) (f a3 b3)
+    (f a4 b4) (f a5 b5) (f a6 b6) (f a7 b7)
+    (f a8 b8) (f a9 b9) (f aA bA) (f aB bB)
+    (f aC bC) (f aD bD) (f aE bE) (f aF bF)
+
+{-# INLINE newMat4 #-}
+newMat4 :: IO Mat4
+newMat4 =
+  IO \world ->
+    let
+      !(# world_, arr_ #) = newAlignedPinnedByteArray# 64# 16# world
+      !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_
+    in
+      (# world, Mat4 arr #)
+
+{-# INLINE scalarMultiply #-}
+scalarMultiply :: Float -> Mat4 -> Mat4
+scalarMultiply x m =
+  withMemory m
+    \ e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF ->
+      fromMemory
+      (e0 * x) (e1 * x) (e2 * x) (e3 * x)
+      (e4 * x) (e5 * x) (e6 * x) (e7 * x)
+      (e8 * x) (e9 * x) (eA * x) (eB * x)
+      (eC * x) (eD * x) (eE * x) (eF * x)
+
+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' #)
+
+instance Block Mat4 where
+  type PackedSize Mat4 = 64
+  alignment140 _  = 16
+  sizeOf140       = sizeOfPacked
+  alignment430    = alignment140
+  sizeOf430       = sizeOf140
+  isStruct _      = False
+  read140     = peekDiffOff
+  write140    = pokeDiffOff
+  read430     = read140
+  write430    = write140
+  readPacked  = read140
+  writePacked = write140
+  {-# INLINE alignment140 #-}
+  {-# INLINE sizeOf140 #-}
+  {-# INLINE alignment430 #-}
+  {-# INLINE sizeOf430 #-}
+  {-# INLINE isStruct #-}
+  {-# INLINE read140 #-}
+  {-# INLINE write140 #-}
+  {-# INLINE read430 #-}
+  {-# INLINE write430 #-}
+  {-# INLINE readPacked #-}
+  {-# INLINE writePacked #-}
+
+toListMemory :: Mat4 -> [Float]
+toListMemory = flip withMemory
+    \ e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF ->
+    [ e0, e1, e2, e3
+    , e4, e5, e6, e7
+    , e8, e9, eA, eB
+    , eC, eD, eE, eF
+    ]
+
+toList2dMemory :: Mat4 -> [[Float]]
+toList2dMemory = flip withMemory
+    \ e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF ->
+    [ [e0, e1, e2, e3]
+    , [e4, e5, e6, e7]
+    , [e8, e9, eA, eB]
+    , [eC, eD, eE, eF]
+    ]
+
+toListTrans :: Mat4 -> [Float]
+toListTrans = flip withMemory
+    \ e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF ->
+    [ e0, e4, e8, eC
+    , e1, e5, e9, eD
+    , e2, e6, eA, eE
+    , e3, e7, eB, eF
+    ]
+
+toList2dTrans :: Mat4 -> [[Float]]
+toList2dTrans = flip withMemory
+    \ e0 e1 e2 e3
+      e4 e5 e6 e7
+      e8 e9 eA eB
+      eC eD eE eF ->
+    [ [e0, e4, e8, eC]
+    , [e1, e5, e9, eD]
+    , [e2, e6, eA, eE]
+    , [e3, e7, eB, eF]
+    ]
+
+zipWith :: (Float -> Float -> c) -> Mat4 -> Mat4 -> [c]
+zipWith f a b = List.zipWith f (toListMemory a) (toListMemory b)
diff --git a/src/Geomancy/Quaternion.hs b/src/Geomancy/Quaternion.hs
--- a/src/Geomancy/Quaternion.hs
+++ b/src/Geomancy/Quaternion.hs
@@ -10,6 +10,8 @@
   , withQuaternion
 
   , axisAngle
+  , intrinsic
+  , extrinsic
   , rotate
   , rotatePoint
   , rotationBetween
@@ -149,14 +151,14 @@
       (d - h)
 
   {-# INLINE (*) #-}
-  Quaternion a b c d * Quaternion e f g h =
-    withVec3 v \y z w ->
-      Quaternion x y z w
+  Quaternion as ax ay az * Quaternion bs bx by bz =
+    withVec3 v \x y z ->
+      Quaternion s x y z
     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
+      s = as * bs - Vec3.dot v1 v2
+      v = Vec3.cross v1 v2 + v2 Vec3.^* as + v1 Vec3.^* bs
+      v1 = vec3 ax ay az
+      v2 = vec3 bx by bz
 
   {-# INLINE fromInteger #-}
   fromInteger x = Quaternion (fromInteger x) 0 0 0
@@ -217,7 +219,10 @@
     where
       ptr' = castPtr ptr
 
--- | Quaternion construction from axis and angle.
+{- | Quaternion construction from axis and angle.
+
+In a right-handed system, the rotation would appear clockwise when looking in the axis direction.
+-}
 {-# INLINE axisAngle #-}
 axisAngle :: Vec3 -> Float -> Quaternion
 axisAngle axis rads =
@@ -226,6 +231,55 @@
   where
     half = rads / 2
 
+{- A composition of roll-pitch-yaw (@Z-X'-Y''@) rotations using local/object's own frame.
+
+Useful for airplane-like controls and first-person camera directions.
+
+The order of application is:
+- Roll is applied first, turning around forward and adjusting the "up" and the "right".
+- Pitch is applied second, turning around the *new* "right" and adjusting the "up" again.
+- Yaw is applied last, turning around the twice-adjusted "up".
+
+Apply the resulting `Quaternion` to the previous local basis to get an updated local basis.
+-}
+intrinsic :: Float -> Float -> Float -> Quaternion
+intrinsic roll pitch yaw =
+  -- XXX: fused construction is a bit faster than chaining axisAngles, but less legible.
+  -- see Spec.hs for reference
+  quaternion
+    -- scalar
+    (cx * cy * cz + sx * sy * sz)
+    -- vector
+    (sx * cy * cz + cx * sy * sz)
+    (cx * sy * cz - sx * cy * sz)
+    (cx * cy * sz - sx * sy * cz)
+  where
+    cx = cos (pitch * 0.5)
+    sx = sin (pitch * 0.5)
+    cy = cos (yaw * 0.5)
+    sy = sin (yaw * 0.5)
+    cz = cos (roll * 0.5)
+    sz = sin (roll * 0.5)
+
+{- A composition of rotations using global/parent/world frame.
+
+Useful for kinematics problems. Use intrinsic if unsure.
+
+The order of application is:
+- Heading is applied first, turning around the "gravity" axis (west-east rotation).
+- Elevation is applied second, adding some down/up rotation wrt. the original right direction.
+- Tilt is applied last, adding some down/up rotation wrt. the original forward direction.
+
+Apply the resulting `Quaternion` to the canonical/parent basis to get a local basis.
+-}
+extrinsic :: Float -> Float -> Float -> Quaternion
+extrinsic heading elevation tilt =
+  axisAngle (vec3 0 0 1) tilt *
+  axisAngle (vec3 1 0 0) elevation *
+  axisAngle (vec3 0 1 0) heading
+
+{- | Rotate point with a unit quaternion.
+-}
 {-# INLINE rotate #-}
 rotate :: Quaternion -> Vec3 -> Vec3
 rotate q v = withQuaternion q' \_a b c d -> vec3 b c d
@@ -233,6 +287,8 @@
     q' = withVec3 v \x y z ->
       q * quaternion 0 x y z * conjugate q
 
+{- | Rotate point around another point with a unit quaternion.
+-}
 {-# INLINE rotatePoint #-}
 rotatePoint :: Quaternion -> Vec3 -> Vec3 -> Vec3
 rotatePoint q origin point =
diff --git a/src/Geomancy/Transform.hs b/src/Geomancy/Transform.hs
--- a/src/Geomancy/Transform.hs
+++ b/src/Geomancy/Transform.hs
@@ -2,11 +2,12 @@
 {-# LANGUAGE DataKinds #-}
 {-# LANGUAGE DerivingStrategies #-}
 {-# LANGUAGE GeneralisedNewtypeDeriving #-}
+{-# LANGUAGE ImportQualifiedPost #-}
 {-# LANGUAGE TypeFamilies #-}
 
 module Geomancy.Transform
   ( Transform(..)
-  , inverse
+  , Mat4.inverse
 
   , apply
   , (!.)
@@ -25,19 +26,21 @@
   , scaleZ
   , scaleXY
   , scale3
+  , scaleV
 
   , dirPos
+  , node
   ) where
 
 import Foreign (Storable(..))
 import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)
 
-import Geomancy.Mat4 (Mat4, colMajor, inverse)
+import Geomancy.Mat4 (Mat4, colMajor)
 import Geomancy.Quaternion (Quaternion, withQuaternion)
 import Geomancy.Vec3 (Vec3, vec3, withVec3)
 import Geomancy.Vec4 (fromVec3, withVec4)
+import Geomancy.Mat4 qualified as Mat4
 
-import qualified Geomancy.Mat4 as Mat4
 import Graphics.Gl.Block (Block(..))
 
 newtype Transform = Transform { unTransform :: Mat4 }
@@ -72,7 +75,12 @@
 apply :: Vec3 -> Transform -> Vec3
 apply = flip (!.)
 
--- | Matrix - column vector multiplication with perspective division
+{- | Matrix - row vector multiplication with perspective division
+
+@
+vOut = pv <> translate !. vIn
+@
+-}
 (!.) :: Transform -> Vec3 -> Vec3
 (!.) mat vec =
   withVec4 res \x y z w ->
@@ -80,6 +88,8 @@
   where
     res = mat Mat4.!* fromVec3 vec 1.0
 
+infixr 5 !.
+
 -- ** Translation
 
 {-# INLINE translate #-}
@@ -124,55 +134,62 @@
 scaleXY :: Float -> Float -> Transform
 scaleXY x y = scale3 x y 1
 
--- ** Euler angle rotations
+{-# INLINE scaleV #-}
+scaleV :: Vec3 -> Transform
+scaleV s = withVec3 s scale3
 
+-- ** Rotation
+
+{- | Clockwise rotation around positive X axis.
+
+Matches @\a -> rotateQ (axisAngle (vec3 1 0 0) a)@.
+-}
 {-# INLINE rotateX #-}
 rotateX :: Float -> Transform
-rotateX rads = colMajor
-  1 0   0   0
-  0 t11 t21 0
-  0 t12 t22 0
-  0 0   0   1
+rotateX rads =
+  colMajor
+    1 0   0  0
+    0 c (-s) 0
+    0 s   c  0
+    0 0   0  1
   where
-    t11 = cost
-    t12 = -sint
-    t21 = sint
-    t22 = cost
+    c = cos rads
+    s = sin rads
 
-    cost = cos rads
-    sint = sin rads
+{- | Clockwise rotation around positive Y axis.
 
+Matches @\a -> rotateQ (axisAngle (vec3 0 1 0) a)@.
+-}
 {-# INLINE rotateY #-}
 rotateY :: Float -> Transform
-rotateY rads = colMajor
-  t00 0 t20 0
-  0   1 0   0
-  t02 0 t22 0
-  0   0 0   1
+rotateY rads =
+  colMajor
+    c  0 s 0
+    0  1 0 0
+  (-s) 0 c 0
+    0  0 0 1
   where
-    cost = cos rads
-    sint = sin rads
+    c = cos rads
+    s = sin rads
 
-    t00 = cost
-    t02 = sint
-    t20 = -sint
-    t22 = cost
+{- | Clockwise rotation around positive Z axis.
 
+Can be used for 2D rotation in the XY plane.
+Matches @\a -> rotateQ (axisAngle (vec3 0 0 1) a)@.
+
+In the right-handed "window coordinates" (e.g. top-left corner is 0,0) "right" becomes "down" after 90deg turn.
+-}
 {-# INLINE rotateZ #-}
 rotateZ :: Float -> Transform
-rotateZ rads = colMajor
-  t00 t10 0 0
-  t01 t11 0 0
-  0   0   1 0
-  0   0   0 1
+rotateZ rads =
+  colMajor
+    c (-s) 0 0
+    s   c  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
+   c = cos rads
+   s = sin rads
 
 {-# INLINE rotateQ #-}
 rotateQ :: Quaternion -> Transform
@@ -198,4 +215,9 @@
         (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
+        0                   0                   0                    1
+
+node :: Vec3 -> Quaternion -> Vec3 -> Transform
+node t r s
+ | s == 1.0 = dirPos r t
+ | otherwise = dirPos r t <> scaleV s -- TODO: finish derivation
diff --git a/src/Geomancy/UVec3.hs b/src/Geomancy/UVec3.hs
--- a/src/Geomancy/UVec3.hs
+++ b/src/Geomancy/UVec3.hs
@@ -1,8 +1,12 @@
 {-# LANGUAGE BlockArguments #-}
 {-# LANGUAGE DataKinds #-}
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE ViewPatterns #-}
 
@@ -28,6 +32,8 @@
 import Foreign (Storable(..))
 import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)
 import GHC.Ix (Ix(..))
+import GHC.OverloadedLabels (IsLabel(..))
+import WebColor.Labels (IsWebColor(..))
 
 import Geomancy.Elementwise (Elementwise(..))
 import Graphics.Gl.Block (Block(..))
@@ -67,6 +73,12 @@
 dot :: UVec3 -> UVec3 -> Word32
 dot (UVec3 l1 l2 l3) (UVec3 r1 r2 r3) =
   l1 * r1 + l2 * r2 + l3 * r3
+
+instance IsWebColor s => IsLabel s UVec3 where
+  {-# INLINE fromLabel #-}
+  fromLabel =
+    webColor @s \r g b ->
+      uvec3 (fromIntegral r) (fromIntegral g) (fromIntegral b)
 
 instance NFData UVec3 where
   rnf UVec3{} = ()
diff --git a/src/Geomancy/Vec3.hs b/src/Geomancy/Vec3.hs
--- a/src/Geomancy/Vec3.hs
+++ b/src/Geomancy/Vec3.hs
@@ -1,13 +1,16 @@
 {-# LANGUAGE BlockArguments #-}
 {-# LANGUAGE DataKinds #-}
+{-# LANGUAGE DerivingStrategies #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE ImportQualifiedPost #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE ViewPatterns #-}
-{-# LANGUAGE DerivingStrategies #-}
 
 -- | Specialized and inlined @V3 Float@.
 
@@ -40,9 +43,11 @@
 import Data.Coerce (Coercible, coerce)
 import Data.MonoTraversable (Element, MonoFunctor(..), MonoPointed(..))
 import Data.VectorSpace (VectorSpace)
+import Data.VectorSpace qualified as VectorSpace
 import Foreign (Storable(..), castPtr)
 import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)
-import qualified Data.VectorSpace as VectorSpace
+import GHC.OverloadedLabels (IsLabel(..))
+import WebColor.Labels (IsWebColor(..))
 
 import Geomancy.Elementwise (Elementwise(..))
 import Graphics.Gl.Block (Block(..))
@@ -85,6 +90,15 @@
 {-# INLINE fromTuple #-}
 fromTuple :: Coercible Vec3 a => (Float, Float, Float) -> a
 fromTuple (x, y, z) = coerce (vec3 x y z)
+
+instance IsWebColor s => IsLabel s Vec3 where
+  {-# INLINE fromLabel #-}
+  fromLabel =
+    webColor @s \r g b ->
+      vec3
+        (fromIntegral r / 255)
+        (fromIntegral g / 255)
+        (fromIntegral b / 255)
 
 instance NFData Vec3 where
   rnf Vec3{} = ()
diff --git a/src/Geomancy/Vec4.hs b/src/Geomancy/Vec4.hs
--- a/src/Geomancy/Vec4.hs
+++ b/src/Geomancy/Vec4.hs
@@ -2,9 +2,13 @@
 {-# LANGUAGE BlockArguments #-}
 {-# LANGUAGE DataKinds #-}
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE ImportQualifiedPost #-}
 {-# LANGUAGE MagicHash #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE UnboxedTuples #-}
 {-# LANGUAGE UnliftedFFITypes #-}
@@ -30,7 +34,7 @@
   , dot
   , normalize
 
-  , unsafeNewVec4
+  , newVec4
   ) where
 
 import GHC.Exts hiding (VecCount(..), toList)
@@ -38,11 +42,13 @@
 import Control.DeepSeq (NFData(rnf))
 import Data.MonoTraversable (Element, MonoFunctor(..), MonoPointed(..))
 import Data.VectorSpace (VectorSpace)
+import Data.VectorSpace qualified as VectorSpace
 import Foreign (Storable(..))
 import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)
 import GHC.IO (IO(..))
+import GHC.OverloadedLabels (IsLabel(..))
 import Text.Printf (printf)
-import qualified Data.VectorSpace as VectorSpace
+import WebColor.Labels (IsWebColorAlpha(..))
 
 import Geomancy.Elementwise (Elementwise(..))
 import Graphics.Gl.Block (Block(..))
@@ -136,6 +142,16 @@
 fromTuple :: (Float, Float, Float, Float) -> Vec4
 fromTuple (x, y, z, w) = vec4 x y z w
 
+instance IsWebColorAlpha s => IsLabel s Vec4 where
+  {-# INLINE fromLabel #-}
+  fromLabel =
+    webColorAlpha @s \r g b a ->
+      vec4
+        (fromIntegral r / 255)
+        (fromIntegral g / 255)
+        (fromIntegral b / 255)
+        (fromIntegral a / 255)
+
 instance NFData Vec4 where
   rnf Vec4{} = ()
 
@@ -402,9 +418,9 @@
 
     nearZero a = abs a <= 1e-6
 
-{-# INLINE unsafeNewVec4 #-}
-unsafeNewVec4 :: IO Vec4
-unsafeNewVec4 =
+{-# INLINE newVec4 #-}
+newVec4 :: IO Vec4
+newVec4 =
   IO \world ->
     let
       !(# world_, arr_ #) = newAlignedPinnedByteArray# 16# 16# world
diff --git a/src/Geomancy/Vulkan/Projection.hs b/src/Geomancy/Vulkan/Projection.hs
--- a/src/Geomancy/Vulkan/Projection.hs
+++ b/src/Geomancy/Vulkan/Projection.hs
@@ -1,72 +1,79 @@
 module Geomancy.Vulkan.Projection
-  ( perspective
-  , infinitePerspective
-  , orthoOffCenter
+  ( reverseDepthRH
+  , reverseDepthOrthoRH
+  , orthoRH
   ) where
 
 import Geomancy.Mat4 (colMajor)
 import Geomancy.Transform (Transform(..))
 
-perspective
-  :: Integral side
-  => Float
+-- | Construct a view-to-NDC transformation.
+--
+-- This will shove a camera frustum (an expanding pyramid) into a Vulkan "clip space" box
+-- with the dimensions [-1; 1] left-to-right, [-1; 1] top-to-bottom, and [1; 0] into the screen.
+--
+-- That is, things further away in the view will be pulled into [0; 0; 0] point on the *back* of the box.
+-- And things that on the near plane set by the argument to this function will be scaled to match
+-- the aspect ratio and the field of view.
+--
+-- When using FoV @pi/2@ and @width@=@height the points with Z=near will keep their positions in the XY plane.
+--
+-- NOTE: To update your code using a vanilla perspective:
+--
+-- Change your depth buffer's clear value.
+--    Instead of clearing to 1.0 (farthest), you now clear to 0.0.
+-- Change your depth comparison function.
+--    Instead of VK_COMPARE_OP_LESS, you must use VK_COMPARE_OP_GREATER.
+reverseDepthRH
+  :: Float
+  -> 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
-
+reverseDepthRH fovRads zn width height =
+  colMajor
+    sx  0 0  0
+    0  sy 0  0
+    0   0 0 zn
+    0   0 1  0
   where
-    x = cotFoV * aspectX
-    y = -cotFoV
-    z = far / (near - far)
-    w23 = near * far / (near - far)
-    w32 = -1
-
-    cotFoV = recip . tan $ 0.5 * fovRads
+    sx = sy * recipAspect
+      where
+        recipAspect = height / width
+    sy = 1 / tan (fovRads / 2)
 
-    aspectX = fromIntegral height / fromIntegral width
+{- | Vanilla orthographic projection centered on @0,0@
 
-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
+@orthoRH 0 1 2 2@ gives identity transform.
+@orthoRH 0 1 800 600@ will map @vec3 400 300 0@ to @vec3 0 0 0@.
+-}
+orthoRH :: Float -> Float -> Float -> Float -> Transform
+orthoRH near far width height =
+  colMajor
+    sx 0 0 0
+    0 sy 0 0
+    0  0 z w
+    0  0 0 1
   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
+    sx = 2 / width
+    sy = 2 / height
+    z = 1 / (far - near)
+    w = near / (near - far)
 
-    w = -2 * near
-    near = 1/128 -- 2048
+{- | Reverse-depth orthographic projection centered on @0,0@
 
-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
+Can be used in the same render pass with reverseDepthRH/VK_COMPARE_OP_GREATER pipelines.
 
+@reverseDepthOrthoRH 0 1 800 600@ will map @vec3 400 300 0@ to @vec3 0 0 1@.
+-}
+reverseDepthOrthoRH :: Float -> Float -> Float -> Float -> Transform
+reverseDepthOrthoRH near far width height =
+  colMajor
+    sx 0 0 0
+    0 sy 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)
+    sx = 2 / width
+    sy = 2 / height
+    z = 1 / (near - far)
+    w = far / (far - near)
diff --git a/src/Geomancy/Vulkan/View.hs b/src/Geomancy/Vulkan/View.hs
--- a/src/Geomancy/Vulkan/View.hs
+++ b/src/Geomancy/Vulkan/View.hs
@@ -1,41 +1,49 @@
 {-# LANGUAGE BlockArguments #-}
 
 module Geomancy.Vulkan.View
-  ( orthoFitScreen
-  , lookAt
+  ( lookAtRH
+  , lookAtRH_
   ) where
 
-import Geomancy.Mat4 (Mat4, rowMajor)
-import Geomancy.Vec3 (Vec3, withVec3)
+import Geomancy.Mat4 (colMajor)
 import Geomancy.Transform (Transform(..))
+import Geomancy.Vec3 (Vec3, cross, dot, normalize, vec3, withVec3)
 
-import qualified Geomancy.Vec3 as Vec3
+{- | Construct a right-handed world-to-view transformation.
 
-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
+This will rotate and translate the world around a static camera.
+This would produce an identity when looking from @vec3 0 0 0@ to @vec3 0 0 1@ using (-Y up).
+
+The fallback "up" should be used when you're looking close to your UP in either direction.
+-}
+lookAtRH
+  :: Vec3 -- ^ Eye position
+  -> Vec3 -- ^ Target
+  -> Vec3 -- ^ Up direction
+  -> Transform
+lookAtRH eye target up =
+  withVec3 rgt \rx ry rz ->
+  withVec3 up' \ux uy uz ->
+  withVec3 fwd \fx fy fz ->
+  colMajor
+    rx ry rz (-er)
+    ux uy uz (-eu)
+    fx fy fz (-ef)
+     0  0  0    1
   where
-    xa = Vec3.normalize $ Vec3.cross za up
-    ya = Vec3.cross xa za
-    za = Vec3.normalize $ center - eye
+    fwd = normalize (target - eye)
+    rgt = normalize (cross fwd up)
+    up' = normalize (cross fwd rgt)
+    er = dot eye rgt
+    eu = dot eye up'
+    ef = dot eye fwd
 
-    xd = - Vec3.dot xa eye
-    yd = - Vec3.dot ya eye
-    zd =   Vec3.dot za eye
+{- | A shortcut for using the -Y as the up axis.
 
-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
+Looking forward from the origin does nothing.
+-}
+{-# INLINE lookAtRH_ #-}
+lookAtRH_ :: Vec3 -> Vec3 -> Transform
+lookAtRH_ eye target = lookAtRH eye target yNeg
   where
-    s = min (screenWidth / targetWidth) (screenHeight / targetHeight)
+    yNeg = vec3 0 (-1) 0
diff --git a/test/Spec.hs b/test/Spec.hs
--- a/test/Spec.hs
+++ b/test/Spec.hs
@@ -1,27 +1,31 @@
 {-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE ImportQualifiedPost #-}
 {-# LANGUAGE NumericUnderscores #-}
 {-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE RecordWildCards #-}
 {-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE TypeApplications #-}
 
+import Geomancy
 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 Control.Monad (when)
 import Data.Foldable (toList)
-import Data.Maybe (catMaybes)
 import GHC.Stack (withFrozenCallStack)
+import Hedgehog.Gen qualified as Gen
+import Hedgehog.Range qualified as Range
+import Linear qualified
+import Linear.Quaternion qualified
 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 Geomancy.Mat4 (withRowMajor, (!*))
+import Geomancy.Mat4 qualified
+import Geomancy.Quaternion qualified
+import Geomancy.Transform ((!.))
+import Geomancy.Transform qualified as Transform
+import Geomancy.Vulkan.Projection qualified as Projection
+import Geomancy.Vulkan.View qualified as View
 
 import Linear ((!*!))
 
@@ -41,43 +45,286 @@
 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
+prop_mat4_assoc :: Property
+prop_mat4_assoc = withTests PROP_TESTS $ property do
+  (p, p_) <- forAllTransform
+  (v, v_) <- forAllTransform
+  (m, m_) <- forAllTransform
 
   let
-    ab'c = (a <> b) <> c
-    a'bc = a <> (b <> c)
-    delta' = nearlyEqualMat4 ab'c a'bc
-  annotateShow delta'
+    pv'm = (p <> v) <> m
+    p'vm = p <> (v <> m)
+    delta' = nearlyEqualMat4 pv'm p'vm
+  replicate 16 Nothing === 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_
+    mv_p = (m_ !*! v_) !*! p_
+    m_vp = m_ !*! (v_ !*! p_)
+    delta_ = nearlyEqualM44 mv_p m_vp
+  replicate 16 Nothing === delta_
 
   -- Inter-library calculated values nearlyEqual
-  ab'c_ <- toM44 ab'c
-  catMaybes (nearlyEqualM44 ab_c ab'c_) === []
+  replicate 16 Nothing === nearlyEqualM44 m_vp (toM44 p'vm)
+  replicate 16 Nothing === nearlyEqualM44 mv_p (toM44 pv'm)
 
-  a'bc_ <- toM44 a'bc
-  catMaybes (nearlyEqualM44 a_bc a'bc_) === []
+prop_mat4_order :: Property
+prop_mat4_order = withTests 1 $ property do
+  let
+    cm = Geomancy.Mat4.colMajor @Geomancy.Mat4.Mat4
+      0 1 2 3
+      4 5 6 7
+      8 9 0 1
+      2 3 4 5
+    rm = Geomancy.Mat4.rowMajor @Geomancy.Mat4.Mat4
+      0 4 8 2
+      1 5 9 3
+      2 6 0 4
+      3 7 1 5
+  Geomancy.Mat4.toListColMajor2d cm === Geomancy.Mat4.toListColMajor2d rm
 
+  let
+    cmT = Geomancy.Mat4.transpose cm
+    rmT = Geomancy.Mat4.transpose rm
+  Geomancy.Mat4.toListColMajor2d cmT === Geomancy.Mat4.toListColMajor2d rmT
+
+prop_transform_det :: Property
+prop_transform_det = withTests 1 $ property do
+  1     === Geomancy.Mat4.det Geomancy.Mat4.identity
+  1     === Geomancy.Mat4.det (Transform.rotateY (pi/6))  -- area-preserving
+  1     === Geomancy.Mat4.det (Transform.translate 1 2 3) -- area-preserving too
+  0     === Geomancy.Mat4.det (Transform.scaleX 0)        -- space collapse, no inverse
+  -1    === Geomancy.Mat4.det (Transform.scaleZ (-1))     -- flips the space
+  8     === Geomancy.Mat4.det (Transform.scale 2)         -- stretches a little
+  (1/8) === Geomancy.Mat4.det (Transform.scale (1/2))     -- stretches backwards
+
+prop_projection_reverseDepthRH :: Property
+prop_projection_reverseDepthRH = withTests 1 $ property do
+  let fov90 = pi / 2
+  let zNear = 0.1
+  let Transform p = Projection.reverseDepthRH fov90 zNear 800 600
+  annotateShow p
+  let
+    -- stating the expected in transposed order for testing
+    Transform expected = Geomancy.Mat4.rowMajor
+      0.7500 0.0000 0.0000 0.0000
+      0.0000 1.0000 0.0000 0.0000
+      0.0000 0.0000 0.0000 1.0000
+      0.0000 0.0000  zNear 0.0000
+  annotateShow expected
+  let delta_ = nearlyEqualMat4 p expected
+  replicate 16 Nothing === delta_
+
+  let
+    centerNearIn = vec4 0 0 zNear 1
+    centerNearOut = p !* centerNearIn
+    centerNearNDC = toNDC centerNearOut
+  annotateShow (centerNearIn, centerNearOut)
+  vec3 0 0 1 === centerNearNDC
+
+prop_projection_orthoRH_identity :: Property
+prop_projection_orthoRH_identity = withTests 1 $ property do
+  let Transform pNative = Projection.orthoRH 0 1 2 2
+  annotateShow pNative
+  let delta_ = nearlyEqualMat4 pNative mempty
+  replicate 16 Nothing === delta_
+
+prop_projection_reverseDepthOrthoRH :: Property
+prop_projection_reverseDepthOrthoRH = withTests 1 $ property do
+  let t = Projection.reverseDepthOrthoRH 1 10 800 600
+  annotateShow t
+
+  vec4 0 0 1 1 === t !* vec4 0 0 1 1 -- center/near
+  vec4 0 0 0 1 === t !* vec4 0 0 10 1 -- center/far
+  vec4 1 1 0 1 === t !* vec4 400 300 10 1 -- edge/far
+  vec4 (-1/400) (-1/300) 1 1 === t !* vec4 (-1) (-1) 1 1 -- step/near
+  vec4 (1/400) (1/300) 0 1 === t !* vec4 1 1 10 1 -- step/far
+
+toNDC :: Vec4 -> Vec3
+toNDC v = withVec4 v \x y z w -> vec3 (x/w) (y/w) (z/w)
+
+prop_view_lookAtRH_identity :: Property
+prop_view_lookAtRH_identity = property do
+  let alreadyThere = View.lookAtRH_ (vec3 0 0 0) (vec3 0 0 1)
+  annotateShow alreadyThere
+
+  p <- forAll genVec3
+
+  p === Transform.apply p alreadyThere
+
+xPos :: Vec3
+xPos = vec3 1 0 0
+
+xNeg :: Vec3
+xNeg = vec3 (-1) 0 0
+
+yPos :: Vec3
+yPos = vec3 0 1 0
+
+yNeg :: Vec3
+yNeg = vec3 0 (-1) 0
+
+zPos :: Vec3
+zPos = vec3 0 0 1
+
+zNeg :: Vec3
+zNeg = vec3 0 0 (-1)
+
+prop_transform_rotateX :: Property
+prop_transform_rotateX = property do
+  let angle = pi / 2
+  let t = Transform.rotateX angle
+  let q = Transform.rotateQ (Geomancy.Quaternion.axisAngle (vec3 1 0 0) angle)
+  let yNegT = t !. zPos
+  annotateShow yNegT
+  let yNegQ = q !. zPos
+  annotateShow yNegQ
+  replicate 3 Nothing === nearlyEqualVec3 yNegT yNeg
+  replicate 3 Nothing === nearlyEqualVec3 yNegQ yNeg
+
+prop_transform_rotateY :: Property
+prop_transform_rotateY = property do
+  let angle = pi / 2
+  let t = Transform.rotateY angle
+  let q = Transform.rotateQ (Geomancy.Quaternion.axisAngle (vec3 0 1 0) angle)
+  let xPosT = t !. zPos
+  annotateShow xPosT
+  let xPosQ = q !. zPos
+  annotateShow xPosQ
+  replicate 3 Nothing === nearlyEqualVec3 xPosT xPos
+  replicate 3 Nothing === nearlyEqualVec3 xPosQ xPos
+
+prop_transform_rotateZ :: Property
+prop_transform_rotateZ = property do
+  let angle = pi / 2
+  let t = Transform.rotateZ angle
+  let q = Transform.rotateQ (Geomancy.Quaternion.axisAngle (vec3 0 0 1) angle)
+  let yPosT = t !. xPos
+  annotateShow yPosT
+  let yPosQ = q !. xPos
+  annotateShow yPosQ
+  replicate 3 Nothing === nearlyEqualVec3 yPosT yPos
+  replicate 3 Nothing === nearlyEqualVec3 yPosQ yPos
+
+prop_quaternion_ref :: Property
+prop_quaternion_ref = property do
+  annotate "Construction"
+  (axis, angle) <- forAll $ (,) <$> genVec3 <*> genAngle
+  let qg = Geomancy.Quaternion.axisAngle axis angle
+  let ql = Linear.Quaternion.axisAngle (withVec3 axis Linear.V3) angle
+  show qg === show (fromLQ ql)
+
+  annotate "Rotation"
+  p <- forAll genVec3
+  let rpg = Geomancy.Quaternion.rotate qg p
+  let rpl = Linear.Quaternion.rotate ql (withVec3 p Linear.V3)
+  show rpg === show (fromV3 rpl)
+
+  annotate "Multiplication"
+  (axisB, angleB) <- forAll $ (,) <$> genVec3 <*> genAngle
+  let qgB = Geomancy.Quaternion.axisAngle axisB angleB
+  let qlB = Linear.Quaternion.axisAngle (withVec3 axisB Linear.V3) angleB
+  let prodg = qg * qgB
+  let prodl = ql * qlB
+  show prodg === show (fromLQ prodl)
+
+prop_quaternion_associativity :: Property
+prop_quaternion_associativity = property do
+  (a, b, c) <- forAll $ (,,)
+    <$> (Geomancy.Quaternion.axisAngle <$> genVec3 <*> genAngle)
+    <*> (Geomancy.Quaternion.axisAngle <$> genVec3 <*> genAngle)
+    <*> (Geomancy.Quaternion.axisAngle <$> genVec3 <*> genAngle)
+  let ab'c = (a * b) * c
+  annotateShow ab'c
+  let a'bc = a * (b * c)
+  annotateShow a'bc
+  replicate 4 Nothing === nearlyEqualQ ab'c a'bc
+
+prop_quaternion_intrinsic :: Property
+prop_quaternion_intrinsic = withTests 1 $ property do
+  let yaw = pi / 2
+  let pitch = pi / 2
+  let roll = 0
+  let q = Geomancy.Quaternion.intrinsic roll pitch yaw
+  let forward' = Geomancy.Quaternion.rotate q zPos
+  let right' = Geomancy.Quaternion.rotate q xPos
+  let up' = Geomancy.Quaternion.rotate q yNeg
+  replicate 3 Nothing === nearlyEqualVec3 yNeg forward'
+  replicate 3 Nothing === nearlyEqualVec3 zNeg right'
+  replicate 3 Nothing === nearlyEqualVec3 xNeg up'
+
+prop_quaternion_extrinsic :: Property
+prop_quaternion_extrinsic = withTests 1 $ property do
+  let heading = pi / 2
+  let elevation = pi / 2
+  let tilt = 0
+  let q = Geomancy.Quaternion.extrinsic heading elevation tilt
+  let forward' = Geomancy.Quaternion.rotate q zPos
+  let right' = Geomancy.Quaternion.rotate q xPos
+  let up' = Geomancy.Quaternion.rotate q yNeg
+  replicate 3 Nothing === nearlyEqualVec3 xPos forward'
+  replicate 3 Nothing === nearlyEqualVec3 yPos right'
+  replicate 3 Nothing === nearlyEqualVec3 zNeg up'
+
+prop_quaternion_intrinsic_ref :: Property
+prop_quaternion_intrinsic_ref = property do
+  (roll, pitch, yaw) <- forAll $ (,,) <$> genAngle <*> genAngle <*> genAngle
+  let fused = Geomancy.Quaternion.intrinsic roll pitch yaw
+  annotateShow fused
+  let composed = intrinsicComposed roll pitch yaw
+  annotateShow composed
+  let viaFused = Geomancy.Quaternion.rotate fused zPos
+  annotateShow viaFused
+  let viaComposed = Geomancy.Quaternion.rotate composed zPos
+  replicate 3 Nothing === nearlyEqualVec3 viaFused viaComposed
+
+intrinsicComposed :: Float -> Float -> Float -> Quaternion
+intrinsicComposed roll pitch yaw =
+  Geomancy.Quaternion.axisAngle (vec3 0 1 0) yaw *
+  Geomancy.Quaternion.axisAngle (vec3 1 0 0) pitch *
+  Geomancy.Quaternion.axisAngle (vec3 0 0 1) roll
+
+prop_transform_node :: Property
+prop_transform_node = withTests 1 $ property do
+  let t = vec3 2 3 5
+  let r = intrinsicComposed 0 0 (pi/7)
+  let s = 11
+  let dps = Transform.node t r s
+  let node = trs t r s
+  let v = vec3 0 0 1
+  replicate 3 Nothing === nearlyEqualVec3 (dps !. v) (node !. v)
+
+trs :: Vec3 -> Quaternion -> Vec3 -> Transform
+trs t r s = mconcat
+  [ Transform.translateV t
+  , Transform.rotateQ r
+  , Transform.scaleV s
+  ]
+
+prop_transform_order :: Property
+prop_transform_order = withTests 1 $ property do
+  let parent = Transform.node (vec3 4 5 6) (intrinsicComposed 0 0 (pi/6)) 2
+  let local = Transform.node (vec3 7 8 9) (intrinsicComposed 0 (pi/6) 0) (1/3)
+  let precomp = parent <> local !. vec3 0 0 1
+  let nocomp = parent !. local !. vec3 0 0 1
+
+  replicate 3 Nothing === nearlyEqualVec3 precomp nocomp
+
 forAllTransform :: PropertyT IO (Geomancy.Mat4.Mat4, Linear.M44 Float)
 forAllTransform = withFrozenCallStack do
   (_name, Transform g) <- forAllWith fst genTransform
-  l <- toM44 g
-  pure (g, l)
+  pure (g, toM44 g)
 
-genTransform :: Gen ([Char], Transform)
+genVec3 :: Gen Vec3
+genVec3 =
+  vec3
+    <$> Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)
+    <*> Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)
+    <*> Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)
+
+genAngle :: Gen Float
+genAngle = Gen.float (Range.linearFracFrom 0.0 (-8*pi) (8*pi))
+
+genTransform :: Gen ([Char], Transform) -- TODO: cross-check with Linear
 genTransform = Gen.choice
   [ genIdentity
   , genTranslate
@@ -91,6 +338,7 @@
       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)
+      when (abs x == 0 && abs y == 0 && abs z == 0) Gen.discard
       pure
         ( printf "translate %0.4f %0.4f %0.4f" x y z
         , Transform.translate x y z
@@ -102,7 +350,8 @@
         , ("rotate/y", Transform.rotateY)
         , ("rotate/z", Transform.rotateZ)
         ]
-      angle <- Gen.float (Range.linearFracFrom 0.0 (-4 * pi) (4 * pi))
+      angle <- genAngle
+      when (abs angle == 0) Gen.discard
       pure
         ( printf "%s %0.4f" name angle
         , axis angle
@@ -112,40 +361,70 @@
       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)
+      when (abs x == 1 && abs y == 1 && abs z == 1) Gen.discard
       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)
+fromV3 :: Linear.V3 Float -> Vec3
+fromV3 (Linear.V3 x y z) = Geomancy.vec3 x y z
 
-toM44 :: MonadIO io => Geomancy.Mat4 -> io (Linear.M44 Float)
-toM44 mat4 =
-  liftIO $
-    Foreign.with mat4 $
-      Foreign.peek . Foreign.castPtr
+fromLQ :: Linear.Quaternion.Quaternion Float -> Geomancy.Quaternion
+fromLQ (Linear.Quaternion.Quaternion s (Linear.V3 x y z)) = Geomancy.quaternion s x y z
 
-type NearlyEqual = Maybe (Float, Float, Float)
+toM44 :: Geomancy.Mat4 -> Linear.M44 Float
+toM44 mat4 =
+  withRowMajor mat4
+    \ m00 m10 m20 m30
+      m01 m11 m21 m31
+      m02 m12 m22 m32
+      m03 m13 m23 m33 ->
+        Linear.V4
+          (Linear.V4 m00 m10 m20 m30)
+          (Linear.V4 m01 m11 m21 m31)
+          (Linear.V4 m02 m12 m22 m32)
+          (Linear.V4 m03 m13 m23 m33)
 
-nearlyEqualMat4 :: Geomancy.Mat4.Mat4 -> Geomancy.Mat4.Mat4 -> [NearlyEqual]
+nearlyEqualMat4 :: Geomancy.Mat4.Mat4 -> Geomancy.Mat4.Mat4 -> [Maybe Oops]
 nearlyEqualMat4 a b = Geomancy.Mat4.zipWith nearlyEqual a b
 
-nearlyEqualM44 :: Linear.M44 Float -> Linear.M44 Float -> [NearlyEqual]
+nearlyEqualM44 :: Linear.M44 Float -> Linear.M44 Float -> [Maybe Oops]
 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
+nearlyEqualVec3 :: Vec3 -> Vec3 -> [Maybe Oops]
+nearlyEqualVec3 a b =
+  withVec3 a \ax ay az ->
+  withVec3 b \bx by bz ->
+    [ nearlyEqual ax bx
+    , nearlyEqual ay by
+    , nearlyEqual az bz
+    ]
+
+nearlyEqualQ :: Geomancy.Quaternion -> Geomancy.Quaternion -> [Maybe Oops]
+nearlyEqualQ a b =
+  Geomancy.withQuaternion a \as ax ay az ->
+    Geomancy.withQuaternion b \bs bx by bz ->
+      [ nearlyEqual as bs
+      , nearlyEqual ax bx
+      , nearlyEqual ay by
+      , nearlyEqual az bz
+      ]
+
+nearlyEqual :: Float -> Float -> Maybe Oops
+nearlyEqual lhs rhs =
+  if lhs == rhs || absDiff < 1e-4 || relDiff < 1e-4 then
     Nothing
   else
-    Just (x, y, abs (1 - x / y))
+    Just Oops{..}
+  where
+    absDiff = abs $ lhs - rhs
+    relDiff = absDiff / (abs lhs + abs rhs)
+    -- fltEps = 1.19209290e-07
+    -- fltMin = 1.175494e-38
+
+data Oops = Oops {lhs :: Float, rhs :: Float, absDiff :: Float, relDiff :: Float}
+  deriving (Eq, Show)
 
 discovered :: Group
 discovered = $$(discover)