diff --git a/.travis.yml b/.travis.yml
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@
 #
 # For more information, see https://github.com/haskell-CI/haskell-ci
 #
-# version: 0.3.20190425
+# version: 0.5.20180830
 #
 language: c
 dist: xenial
@@ -35,38 +35,37 @@
   include:
     - compiler: ghc-8.8.1
       addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.8.1","cabal-install-3.0"]}}
-      env: GHCHEAD=true
     - compiler: ghc-8.6.5
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.6.5","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.6.5","cabal-install-3.0"]}}
     - compiler: ghc-8.4.4
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.4.4","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.4.4","cabal-install-3.0"]}}
     - compiler: ghc-8.2.2
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.2.2","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.2.2","cabal-install-3.0"]}}
     - compiler: ghc-8.0.2
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.0.2","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.0.2","cabal-install-3.0"]}}
     - compiler: ghc-7.10.3
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.10.3","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.10.3","cabal-install-3.0"]}}
     - compiler: ghc-7.8.4
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.8.4","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.8.4","cabal-install-3.0"]}}
     - compiler: ghc-7.6.3
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.6.3","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.6.3","cabal-install-3.0"]}}
     - compiler: ghc-7.4.2
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.4.2","cabal-install-2.4"]}}
+      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.4.2","cabal-install-3.0"]}}
     - compiler: ghc-head
       addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-head","cabal-install-head"]}}
-      env: GHCHEAD=true
   allow_failures:
     - compiler: ghc-head
-    - compiler: ghc-8.8.1
 before_install:
   - HC=$(echo "/opt/$CC/bin/ghc" | sed 's/-/\//')
+  - WITHCOMPILER="-w $HC"
+  - HADDOCK=$(echo "/opt/$CC/bin/haddock" | sed 's/-/\//')
   - HCPKG="$HC-pkg"
   - unset CC
   - CABAL=/opt/ghc/bin/cabal
   - CABALHOME=$HOME/.cabal
   - export PATH="$CABALHOME/bin:$PATH"
   - TOP=$(pwd)
-  - HCNUMVER=$(( $(${HC} --numeric-version|sed -E 's/([0-9]+)\.([0-9]+)\.([0-9]+).*/\1 * 10000 + \2 * 100 + \3/') ))
+  - "HCNUMVER=$(${HC} --numeric-version|perl -ne '/^(\\d+)\\.(\\d+)\\.(\\d+)(\\.(\\d+))?$/; print(10000 * $1 + 100 * $2 + ($3 == 0 ? $5 != 1 : $3))')"
   - echo $HCNUMVER
   - CABAL="$CABAL -vnormal+nowrap+markoutput"
   - set -o pipefail
@@ -93,11 +92,13 @@
   - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"
   - TEST=--enable-tests
   - BENCH=--enable-benchmarks
-  - GHCHEAD=${GHCHEAD-false}
+  - HEADHACKAGE=false
+  - if [ $HCNUMVER -gt 80801 ] ; then HEADHACKAGE=true ; fi
   - rm -f $CABALHOME/config
   - |
     echo "verbose: normal +nowrap +markoutput"          >> $CABALHOME/config
     echo "remote-build-reporting: anonymous"            >> $CABALHOME/config
+    echo "write-ghc-environment-files: always"          >> $CABALHOME/config
     echo "remote-repo-cache: $CABALHOME/packages"       >> $CABALHOME/config
     echo "logs-dir:          $CABALHOME/logs"           >> $CABALHOME/config
     echo "world-file:        $CABALHOME/world"          >> $CABALHOME/config
@@ -111,15 +112,14 @@
     echo "repository hackage.haskell.org"               >> $CABALHOME/config
     echo "  url: http://hackage.haskell.org/"           >> $CABALHOME/config
   - |
-    if $GHCHEAD; then
+    if $HEADHACKAGE; then
     echo "allow-newer: $($HCPKG list --simple-output | sed -E 's/([a-zA-Z-]+)-[0-9.]+/*:\1/g')" >> $CABALHOME/config
-    
-    echo "repository head.hackage"                                                        >> $CABALHOME/config
-    echo "   url: http://head.hackage.haskell.org/"                                       >> $CABALHOME/config
+    echo "repository head.hackage.ghc.haskell.org"                                        >> $CABALHOME/config
+    echo "   url: https://ghc.gitlab.haskell.org/head.hackage/"                           >> $CABALHOME/config
     echo "   secure: True"                                                                >> $CABALHOME/config
-    echo "   root-keys: 07c59cb65787dedfaef5bd5f987ceb5f7e5ebf88b904bbd4c5cbdeb2ff71b740" >> $CABALHOME/config
-    echo "              2e8555dde16ebd8df076f1a8ef13b8f14c66bad8eafefd7d9e37d0ed711821fb" >> $CABALHOME/config
-    echo "              8f79fd2389ab2967354407ec852cbe73f2e8635793ac446d09461ffb99527f6e" >> $CABALHOME/config
+    echo "   root-keys: 7541f32a4ccca4f97aea3b22f5e593ba2c0267546016b992dfadcd2fe944e55d" >> $CABALHOME/config
+    echo "              26021a13b401500c8eb2761ca95c61f2d625bfef951b939a8124ed12ecf07329" >> $CABALHOME/config
+    echo "              f76d08be13e9a61a377a85e2fb63f4c5435d40f8feb3e12eb05905edb8cdea89" >> $CABALHOME/config
     echo "   key-threshold: 3"                                                            >> $CABALHOME/config
     fi
   - cat $CABALHOME/config
@@ -129,17 +129,16 @@
   - rm -rf cabal.project cabal.project.local cabal.project.freeze
   - touch cabal.project
   - |
-    echo 'packages: "."' >> cabal.project
+    echo "packages: ." >> cabal.project
   - |
-    echo "write-ghc-environment-files: always" >> cabal.project
-  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | grep -vE -- '^(linear)$' | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
+  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(linear)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
   - cat cabal.project || true
   - cat cabal.project.local || true
   - if [ -f "./configure.ac" ]; then (cd "." && autoreconf -i); fi
-  - ${CABAL} v2-freeze -w ${HC} ${TEST} ${BENCH} | color_cabal_output
+  - ${CABAL} v2-freeze $WITHCOMPILER ${TEST} ${BENCH} | color_cabal_output
   - "cat cabal.project.freeze | sed -E 's/^(constraints: *| *)//' | sed 's/any.//'"
   - rm  cabal.project.freeze
-  - ${CABAL} v2-build -w ${HC} ${TEST} ${BENCH} --dep -j2 all | color_cabal_output
+  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} --dep -j2 all | color_cabal_output
 script:
   - DISTDIR=$(mktemp -d /tmp/dist-test.XXXX)
   # Packaging...
@@ -147,26 +146,27 @@
   # Unpacking...
   - mv dist-newstyle/sdist/*.tar.gz ${DISTDIR}/
   - cd ${DISTDIR} || false
-  - find . -maxdepth 1 -name '*.tar.gz' -exec tar -xvf '{}' \;
+  - find . -maxdepth 1 -type f -name '*.tar.gz' -exec tar -xvf '{}' \;
+  - find . -maxdepth 1 -type f -name '*.tar.gz' -exec rm       '{}' \;
+  - PKGDIR_linear="$(find . -maxdepth 1 -type d -regex '.*/linear-[0-9.]*')"
   # Generate cabal.project
   - rm -rf cabal.project cabal.project.local cabal.project.freeze
   - touch cabal.project
   - |
-    echo 'packages: "linear-*/*.cabal"' >> cabal.project
+    echo "packages: ${PKGDIR_linear}" >> cabal.project
   - |
-    echo "write-ghc-environment-files: always" >> cabal.project
-  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | grep -vE -- '^(linear)$' | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
+  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(linear)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
   - cat cabal.project || true
   - cat cabal.project.local || true
   # Building with tests and benchmarks...
   # build & run tests, build benchmarks
-  - ${CABAL} v2-build -w ${HC} ${TEST} ${BENCH} all | color_cabal_output
+  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} all | color_cabal_output
   # Testing...
-  - ${CABAL} v2-test -w ${HC} ${TEST} ${BENCH} all | color_cabal_output
+  - ${CABAL} v2-test $WITHCOMPILER ${TEST} ${BENCH} all | color_cabal_output
   # cabal check...
-  - (cd linear-* && ${CABAL} -vnormal check)
+  - (cd ${PKGDIR_linear} && ${CABAL} -vnormal check)
   # haddock...
-  - ${CABAL} v2-haddock -w ${HC} ${TEST} ${BENCH} all | color_cabal_output
+  - ${CABAL} v2-haddock $WITHCOMPILER --with-haddock $HADDOCK ${TEST} ${BENCH} all | color_cabal_output
 
 # REGENDATA ["--output=.travis.yml","--config=cabal.haskell-ci","cabal.project"]
 # EOF
diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,19 @@
+1.21 [2020.02.03]
+-----------------
+* Add instances for direct sums (`Product`) and tensor products (`Compose`) of
+  other vector spaces. This makes is much more convenient to do things like treat
+  a matrix temporarily as a vector through Compose, or to consider things like
+  Gauss-Jordan elimination, which wants augmented structures.
+* Add `frobenius` for computing the Frobenius norm of a matrix.
+* Added `Random` instances for `System.Random`. We had an indirect dependency
+  through `vector` anyways.
+* Add "obvious" zipping `Semigroup` and `Monoid` instances to all the
+  representable vector spaces.
+* Add `R1`..`R4` instances to `Quaternion`. `_w` is the scalar component so that
+  `_x`,`_y`,`_z` can be directional.
+* Add more solvers to `Linear.Matrix`, available with `base-4.8` or later.
+* Add `unangle` function to `Linear.V2`.
+
 1.20.9 [2019.05.02]
 -------------------
 * Derive `Lift` instances for `Plucker`, `Quaternion`, and `V{0,1,2,3,4}`.
diff --git a/linear.cabal b/linear.cabal
--- a/linear.cabal
+++ b/linear.cabal
@@ -1,6 +1,6 @@
 name:          linear
 category:      Math, Algebra
-version:       1.20.9
+version:       1.21
 license:       BSD3
 cabal-version: >= 1.8
 license-file:  LICENSE
@@ -67,8 +67,9 @@
     deepseq              >= 1.1   && < 1.5,
     distributive         >= 0.2.2 && < 1,
     ghc-prim,
-    hashable             >= 1.1   && < 1.3,
+    hashable             >= 1.1   && < 1.4,
     lens                 >= 4.15.2 && < 5,
+    random               >= 1.0   && < 1.2,
     reflection           >= 1.3.2 && < 3,
     semigroups           >= 0.9   && < 1,
     semigroupoids        >= 5.2.1 && < 6,
diff --git a/src/Linear/Affine.hs b/src/Linear/Affine.hs
--- a/src/Linear/Affine.hs
+++ b/src/Linear/Affine.hs
@@ -47,6 +47,7 @@
 import Data.Foldable as Foldable
 import Data.Functor.Bind
 import Data.Functor.Classes
+import Data.Functor.Product
 import Data.Functor.Rep as Rep
 import Data.HashMap.Lazy (HashMap)
 import Data.Hashable
@@ -56,6 +57,12 @@
 import Data.IntMap (IntMap)
 import Data.Ix
 import Data.Map (Map)
+#if !(MIN_VERSION_base(4,8,0))
+import Data.Monoid (Monoid(..))
+#endif
+#if !(MIN_VERSION_base(4,11,0))
+import Data.Semigroup (Semigroup)
+#endif
 import Data.Serialize as Cereal
 import Data.Vector (Vector)
 import qualified Data.Vector.Generic.Mutable as M
@@ -79,6 +86,7 @@
 import Linear.V3
 import Linear.V4
 import Linear.Vector
+import System.Random
 
 #ifdef HLINT
 {-# ANN module "HLint: ignore Unused LANGUAGE pragma" #-}
@@ -107,6 +115,12 @@
   p .-^ v = p .+^ negated v
   {-# INLINE (.-^) #-}
 
+instance (Affine f, Affine g) => Affine (Product f g) where
+  type Diff (Product f g) = Product (Diff f) (Diff g)
+  Pair a b .-. Pair c d = Pair (a .-. c) (b .-. d)
+  Pair a b .+^ Pair c d = Pair (a .+^ c) (b .+^ d)
+  Pair a b .-^ Pair c d = Pair (a .+^ c) (b .+^ d)
+
 -- | Compute the quadrance of the difference (the square of the distance)
 qdA :: (Affine p, Foldable (Diff p), Num a) => p a -> p a -> a
 qdA a b = Foldable.sum (fmap (join (*)) (a .-. b))
@@ -148,7 +162,8 @@
            , Eq1, Ord1, Show1, Read1
            , Traversable, Apply, Additive, Metric
            , Fractional , Num, Ix, Storable, Epsilon
-           , Hashable
+           , Semigroup, Monoid
+           , Random, Hashable
 #if __GLASGOW_HASKELL__ >= 702
            , Generic
 #endif
diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -1,5 +1,7 @@
+{-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE CPP #-}
 {-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 #if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
 {-# LANGUAGE Trustworthy #-}
@@ -33,6 +35,20 @@
   , _m22, _m23, _m24
   , _m32, _m33, _m34
   , _m42, _m43, _m44
+#if MIN_VERSION_base(4,8,0)
+  , lu
+  , luFinite
+  , forwardSub
+  , forwardSubFinite
+  , backwardSub
+  , backwardSubFinite
+  , luSolve
+  , luSolveFinite
+  , luInv
+  , luInvFinite
+  , luDet
+  , luDetFinite
+#endif
   ) where
 
 #if __GLASGOW_HASKELL__ < 710
@@ -51,6 +67,11 @@
 import Linear.Conjugate
 import Linear.Trace
 
+#if MIN_VERSION_base(4,8,0)
+import GHC.TypeLits
+import Linear.V
+#endif
+
 #ifdef HLINT
 {-# ANN module "HLint: ignore Reduce duplication" #-}
 #endif
@@ -418,3 +439,307 @@
                        (-i30 * s3 + i31 * s1 - i32 * s0)
                        (i20 * s3 - i21 * s1 + i22 * s0))
 {-# INLINE inv44 #-}
+
+#if MIN_VERSION_base(4,8,0)
+-- | Compute the (L, U) decomposition of a square matrix using Crout's
+--   algorithm. The 'Index' of the vectors must be 'Integral'.
+lu :: ( Num a
+      , Fractional a
+      , Foldable m
+      , Traversable m
+      , Applicative m
+      , Additive m
+      , Ixed (m a)
+      , Ixed (m (m a))
+      , i ~ Index (m a)
+      , i ~ Index (m (m a))
+      , Eq i
+      , Integral i
+      , a ~ IxValue (m a)
+      , m a ~ IxValue (m (m a))
+      , Num (m a)
+      )
+   => m (m a)
+   -> (m (m a), m (m a))
+lu a =
+    let n = fromIntegral (length a)
+        initU = identity
+        initL = zero
+        buildLVal !i !j !l !u =
+            let go !k !s
+                    | k == j = s
+                    | otherwise = go (k+1)
+                                     ( s
+                                      + ( (l ^?! ix i ^?! ix k)
+                                        * (u ^?! ix k ^?! ix j)
+                                        )
+                                      )
+                s' = go 0 0
+            in l & (ix i . ix j) .~ ((a ^?! ix i ^?! ix j) - s')
+        buildL !i !j !l !u
+            | i == n = l
+            | otherwise = buildL (i+1) j (buildLVal i j l u) u
+        buildUVal !i !j !l !u =
+            let go !k !s
+                    | k == j = s
+                    | otherwise = go (k+1)
+                                     ( s
+                                     + ( (l ^?! ix j ^?! ix k)
+                                       * (u ^?! ix k ^?! ix i)
+                                       )
+                                     )
+                s' = go 0 0
+            in u & (ix j . ix i) .~ ( ((a ^?! ix j ^?! ix i) - s')
+                                    / (l ^?! ix j ^?! ix j)
+                                    )
+        buildU !i !j !l !u
+            | i == n = u
+            | otherwise = buildU (i+1) j l (buildUVal i j l u)
+        buildLU !j !l !u
+            | j == n = (l, u)
+            | otherwise =
+                let l' = buildL j j l u
+                    u' = buildU j j l' u
+                in buildLU (j+1) l' u'
+    in buildLU 0 initL initU
+
+-- | Compute the (L, U) decomposition of a square matrix using Crout's
+--   algorithm, using the vector's 'Finite' instance to provide an index.
+luFinite :: ( Num a
+            , Fractional a
+            , Functor m
+            , Finite m
+            , n ~ Size m
+            , KnownNat n
+            , Num (m a)
+            )
+         => m (m a)
+         -> (m (m a), m (m a))
+luFinite a =
+    bimap (fmap fromV . fromV)
+          (fmap fromV . fromV)
+          (lu (fmap toV (toV a)))
+
+-- | Solve a linear system with a lower-triangular matrix of coefficients with
+--   forwards substitution.
+forwardSub :: ( Num a
+              , Fractional a
+              , Foldable m
+              , Additive m
+              , Ixed (m a)
+              , Ixed (m (m a))
+              , i ~ Index (m a)
+              , i ~ Index (m (m a))
+              , Eq i
+              , Ord i
+              , Integral i
+              , a ~ IxValue (m a)
+              , m a ~ IxValue (m (m a))
+              )
+           => m (m a)
+           -> m a
+           -> m a
+forwardSub a b =
+    let n = fromIntegral (length b)
+        initX = zero
+        coeff !i !j !s !x
+            | j == i = s
+            | otherwise = coeff i (j+1) (s + ((a ^?! ix i ^?! ix j) * (x ^?! ix j))) x
+        go !i !x
+            | i == n = x
+            | otherwise = go (i + 1) (x & ix i .~ ( ((b ^?! ix i) - coeff i 0 0 x)
+                                                  / (a ^?! ix i ^?! ix i)
+                                                  ))
+    in go 0 initX
+
+-- | Solve a linear system with a lower-triangular matrix of coefficients with
+--   forwards substitution, using the vector's 'Finite' instance to provide an
+--   index.
+forwardSubFinite :: ( Num a
+                    , Fractional a
+                    , Foldable m
+                    , n ~ Size m
+                    , KnownNat n
+                    , Additive m
+                    , Finite m
+                    )
+                 => m (m a)
+                 -> m a
+                 -> m a
+forwardSubFinite a b = fromV (forwardSub (fmap toV (toV a)) (toV b))
+
+-- | Solve a linear system with an upper-triangular matrix of coefficients with
+--   backwards substitution.
+backwardSub :: ( Num a
+               , Fractional a
+               , Foldable m
+               , Additive m
+               , Ixed (m a)
+               , Ixed (m (m a))
+               , i ~ Index (m a)
+               , i ~ Index (m (m a))
+               , Eq i
+               , Ord i
+               , Integral i
+               , a ~ IxValue (m a)
+               , m a ~ IxValue (m (m a))
+               )
+            => m (m a)
+            -> m a
+            -> m a
+backwardSub a b =
+    let n = fromIntegral (length b)
+        initX = zero
+        coeff !i !j !s !x
+            | j == n = s
+            | otherwise = coeff i
+                                (j+1)
+                                (s + ((a ^?! ix i ^?! ix j) * (x ^?! ix j)))
+                                x
+        go !i !x
+            | i < 0 = x
+            | otherwise = go (i-1)
+                             (x & ix i .~ ( ((b ^?! ix i) - coeff i (i+1) 0 x)
+                                          / (a ^?! ix i ^?! ix i)
+                                          ))
+    in go (n-1) initX
+
+-- | Solve a linear system with an upper-triangular matrix of coefficients with
+--   backwards substitution, using the vector's 'Finite' instance to provide an
+--   index.
+backwardSubFinite :: ( Num a
+                     , Fractional a
+                     , Foldable m
+                     , n ~ Size m
+                     , KnownNat n
+                     , Additive m
+                     , Finite m
+                     )
+                  => m (m a)
+                  -> m a
+                  -> m a
+backwardSubFinite a b = fromV (backwardSub (fmap toV (toV a)) (toV b))
+
+-- | Solve a linear system with LU decomposition.
+luSolve :: ( Num a
+           , Fractional a
+           , Foldable m
+           , Traversable m
+           , Applicative m
+           , Additive m
+           , Ixed (m a)
+           , Ixed (m (m a))
+           , i ~ Index (m a)
+           , i ~ Index (m (m a))
+           , Eq i
+           , Integral i
+           , a ~ IxValue (m a)
+           , m a ~ IxValue (m (m a))
+           , Num (m a)
+           )
+        => m (m a)
+        -> m a
+        -> m a
+luSolve a b =
+    let (l, u) = lu a
+    in backwardSub u (forwardSub l b)
+
+-- | Solve a linear system with LU decomposition, using the vector's 'Finite'
+--   instance to provide an index.
+luSolveFinite :: ( Num a
+                 , Fractional a
+                 , Functor m
+                 , Finite m
+                 , n ~ Size m
+                 , KnownNat n
+                 , Num (m a)
+                 )
+              => m (m a)
+              -> m a
+              -> m a
+luSolveFinite a b = fromV (luSolve (fmap toV (toV a)) (toV b))
+
+-- | Invert a matrix with LU decomposition.
+luInv :: ( Num a
+         , Fractional a
+         , Foldable m
+         , Traversable m
+         , Applicative m
+         , Additive m
+         , Distributive m
+         , Ixed (m a)
+         , Ixed (m (m a))
+         , i ~ Index (m a)
+         , i ~ Index (m (m a))
+         , Eq i
+         , Integral i
+         , a ~ IxValue (m a)
+         , m a ~ IxValue (m (m a))
+         , Num (m a)
+         )
+      => m (m a)
+      -> m (m a)
+luInv a =
+    let n = fromIntegral (length a)
+        initA' = zero
+        (l, u) = lu a
+        go !i !a'
+            | i == n = a'
+            | otherwise = let e   = zero & ix i .~ 1
+                              a'r = backwardSub u (forwardSub l e)
+                          in go (i+1) (a' & ix i .~ a'r)
+    in transpose (go 0 initA')
+
+-- | Invert a matrix with LU decomposition, using the vector's 'Finite' instance
+--   to provide an index.
+luInvFinite :: ( Num a
+               , Fractional a
+               , Functor m
+               , Finite m
+               , n ~ Size m
+               , KnownNat n
+               , Num (m a)
+               )
+            => m (m a)
+            -> m (m a)
+luInvFinite a = fmap fromV (fromV (luInv (fmap toV (toV a))))
+
+-- | Compute the determinant of a matrix using LU decomposition.
+luDet :: ( Num a
+         , Fractional a
+         , Foldable m
+         , Traversable m
+         , Applicative m
+         , Additive m
+         , Trace m
+         , Ixed (m a)
+         , Ixed (m (m a))
+         , i ~ Index (m a)
+         , i ~ Index (m (m a))
+         , Eq i
+         , Integral i
+         , a ~ IxValue (m a)
+         , m a ~ IxValue (m (m a))
+         , Num (m a)
+         )
+      => m (m a)
+      -> a
+luDet a =
+    let (l, u) = lu a
+        p      = Foldable.foldl (*) 1
+    in (p (diagonal l)) * (p (diagonal u))
+
+-- | Compute the determinant of a matrix using LU decomposition, using the
+--   vector's 'Finite' instance to provide an index.
+luDetFinite :: ( Num a
+               , Fractional a
+               , Functor m
+               , Finite m
+               , n ~ Size m
+               , KnownNat n
+               , Num (m a)
+               )
+            => m (m a)
+            -> a
+luDetFinite = luDet . fmap toV . toV
+#endif
diff --git a/src/Linear/Metric.hs b/src/Linear/Metric.hs
--- a/src/Linear/Metric.hs
+++ b/src/Linear/Metric.hs
@@ -20,7 +20,9 @@
 
 import Control.Applicative
 import Data.Foldable as Foldable
+import Data.Functor.Compose
 import Data.Functor.Identity
+import Data.Functor.Product
 import Data.Vector (Vector)
 import Data.IntMap (IntMap)
 import Data.Map (Map)
@@ -31,6 +33,7 @@
 
 -- $setup
 -- >>> import Linear
+--
 
 -- | Free and sparse inner product/metric spaces.
 class Additive f => Metric f where
@@ -67,6 +70,18 @@
   signorm v = fmap (/m) v where
     m = norm v
 
+instance (Metric f, Metric g) => Metric (Product f g) where
+  dot (Pair a b) (Pair c d) = dot a c + dot b d
+  quadrance (Pair a b) = quadrance a + quadrance b
+  qd (Pair a b) (Pair c d) = qd a c + qd b d
+  distance p q = sqrt (qd p q)
+
+instance (Metric f, Metric g) => Metric (Compose f g) where
+  dot (Compose a) (Compose b) = quadrance (liftI2 dot a b)
+  quadrance = quadrance . fmap quadrance . getCompose
+  qd (Compose a) (Compose b) = quadrance (liftI2 qd a b)
+  distance (Compose a) (Compose b) = norm (liftI2 qd a b)
+
 instance Metric Identity where
   dot (Identity x) (Identity y) = x * y
 
@@ -95,4 +110,3 @@
 -- | @project u v@ computes the projection of @v@ onto @u@.
 project :: (Metric v, Fractional a) => v a -> v a -> v a
 project u v = ((v `dot` u) / quadrance u) *^ u
-
diff --git a/src/Linear/Plucker.hs b/src/Linear/Plucker.hs
--- a/src/Linear/Plucker.hs
+++ b/src/Linear/Plucker.hs
@@ -104,6 +104,7 @@
 import Linear.V3
 import Linear.V4
 import Linear.Vector
+import System.Random
 
 {-# ANN module "HLint: ignore Reduce duplication" #-}
 
@@ -127,6 +128,22 @@
   fromV (V v) = Plucker (v V.! 0) (v V.! 1) (v V.! 2) (v V.! 3) (v V.! 4) (v V.! 5)
 #endif
 
+instance Random a => Random (Plucker a) where
+  random g = case random g of
+    (a, g1) -> case random g1 of
+      (b, g2) -> case random g2 of
+        (c, g3) -> case random g3 of
+          (d, g4) -> case random g4 of
+            (e, g5) -> case random g5 of
+              (f, g6) -> (Plucker a b c d e f, g6)
+  randomR (Plucker a b c d e f, Plucker a' b' c' d' e' f') g = case randomR (a,a') g of
+    (a'', g1) -> case randomR (b,b') g1 of
+      (b'', g2) -> case randomR (c,c') g2 of
+        (c'', g3) -> case randomR (d,d') g3 of
+          (d'', g4) -> case randomR (e,e') g4 of
+            (e'', g5) -> case randomR (f,f') g5 of
+              (f'', g6) -> (Plucker a'' b'' c'' d'' e'' f'', g6)
+
 instance Functor Plucker where
   fmap g (Plucker a b c d e f) = Plucker (g a) (g b) (g c) (g d) (g e) (g f)
   {-# INLINE fmap #-}
@@ -671,3 +688,12 @@
 
 instance Field6 (Plucker a) (Plucker a) a a where
   _6 f (Plucker x y z u v w) = f w <&> \w' -> Plucker x y z u v w'
+
+instance Semigroup a => Semigroup (Plucker a) where
+ (<>) = liftA2 (<>)
+
+instance Monoid a => Monoid (Plucker a) where
+  mempty = pure mempty
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = liftA2 mappend
+#endif
diff --git a/src/Linear/Quaternion.hs b/src/Linear/Quaternion.hs
--- a/src/Linear/Quaternion.hs
+++ b/src/Linear/Quaternion.hs
@@ -24,6 +24,11 @@
 #ifndef MIN_VERSION_vector
 #define MIN_VERSION_vector(x,y,z) 1
 #endif
+
+#ifndef MIN_VERSION_base
+#define MIN_VERSION_base(x,y,z) 1
+#endif
+
 -----------------------------------------------------------------------------
 -- |
 -- Copyright   :  (C) 2012-2015 Edward Kmett
@@ -72,10 +77,15 @@
 #if (MIN_VERSION_hashable(1,2,5))
 import Data.Hashable.Lifted
 #endif
+#if !(MIN_VERSION_base(4,11,0))
+import Data.Semigroup (Semigroup(..))
+#endif
 import Data.Serialize as Cereal
 import GHC.Arr (Ix(..))
 import qualified Data.Foldable as F
-import Data.Monoid
+#if !(MIN_VERSION_base(4,8,0))
+import Data.Monoid (Monoid(..))
+#endif
 #if __GLASGOW_HASKELL__ >= 707
 import qualified Data.Vector as V
 #endif
@@ -99,9 +109,12 @@
 #if __GLASGOW_HASKELL__ >= 707
 import Linear.V
 #endif
+import Linear.V2
 import Linear.V3
+import Linear.V4
 import Linear.Vector
 import Prelude hiding (any)
+import System.Random
 
 {-# ANN module "HLint: ignore Reduce duplication" #-}
 
@@ -126,6 +139,14 @@
   fromV (V v) = Quaternion (v V.! 0) (V3 (v V.! 1) (v V.! 2) (v V.! 3))
 #endif
 
+instance Random a => Random (Quaternion a) where
+  random g = case random g of
+    (a, g') -> case random g' of
+      (b, g'') -> (Quaternion a b, g'')
+  randomR (Quaternion a b, Quaternion c d) g = case randomR (a,c) g of
+    (e, g') -> case randomR (b,d) g' of
+      (f, g'') -> (Quaternion e f, g'')
+
 instance Functor Quaternion where
   fmap f (Quaternion e v) = Quaternion (f e) (fmap f v)
   {-# INLINE fmap #-}
@@ -693,3 +714,28 @@
 
 instance Field4 (Quaternion a) (Quaternion a) a a where
   _4 f (Quaternion w (V3 x y z)) = f z <&> \z' -> Quaternion w (V3 x y z')
+
+instance Semigroup a => Semigroup (Quaternion a) where
+ (<>) = liftA2 (<>)
+
+instance Monoid a => Monoid (Quaternion a) where
+  mempty = pure mempty
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = liftA2 mappend
+#endif
+
+instance R1 Quaternion where
+  _x f (Quaternion w (V3 x y z)) = f x <&> \x' -> Quaternion w (V3 x' y z)
+
+instance R2 Quaternion where
+  _y f (Quaternion w (V3 x y z)) = f y <&> \y' -> Quaternion w (V3 x y' z)
+  _xy f (Quaternion w (V3 x y z)) = f (V2 x y) <&> \(V2 x' y') -> Quaternion w (V3 x' y' z)
+
+instance R3 Quaternion where
+  _z f (Quaternion w (V3 x y z)) = f z <&> \z' -> Quaternion w (V3 x y z')
+  _xyz f (Quaternion w xyz) = Quaternion w <$> f xyz
+
+instance R4 Quaternion where
+  _w f (Quaternion w xyz) = f w <&> \w' -> Quaternion w' xyz
+  _xyzw f (Quaternion w (V3 x y z)) = f (V4 x y z w) <&> \(V4 x' y' z' w') -> Quaternion w' (V3 x' y' z')
+
diff --git a/src/Linear/Trace.hs b/src/Linear/Trace.hs
--- a/src/Linear/Trace.hs
+++ b/src/Linear/Trace.hs
@@ -20,6 +20,7 @@
 ---------------------------------------------------------------------------
 module Linear.Trace
   ( Trace(..)
+  , frobenius
   ) where
 
 import Control.Monad as Monad
@@ -31,6 +32,7 @@
 import Linear.Plucker
 import Linear.Quaternion
 import Linear.V
+import Linear.Vector
 #if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ > 704
 import Data.Complex
 #endif
@@ -116,3 +118,7 @@
   {-# INLINE trace #-}
   diagonal = Compose . fmap diagonal . diagonal . fmap distribute . getCompose . fmap getCompose
   {-# INLINE diagonal #-}
+
+-- | Compute the <http://mathworld.wolfram.com/FrobeniusNorm.html Frobenius norm> of a matrix.
+frobenius :: (Num a, Foldable f, Additive f, Additive g, Distributive g, Trace g) => f (g a) -> a
+frobenius m = trace $ fmap (\ f' -> Foldable.foldl' (^+^) zero $ liftI2 (*^) f' m) (distribute m)
diff --git a/src/Linear/V.hs b/src/Linear/V.hs
--- a/src/Linear/V.hs
+++ b/src/Linear/V.hs
@@ -32,6 +32,10 @@
 #define MIN_VERSION_transformers(x,y,z) 1
 #endif
 
+#ifndef MIN_VERSION_base
+#define MIN_VERSION_base(x,y,z) 1
+#endif
+
 -----------------------------------------------------------------------------
 -- |
 -- Copyright   :  (C) 2012-2015 Edward Kmett
@@ -64,12 +68,11 @@
 #endif
   ) where
 
-#if __GLASGOW_HASKELL__ < 710
 import Control.Applicative
-#endif
 import Control.DeepSeq (NFData)
 import Control.Monad
 import Control.Monad.Fix
+import Control.Monad.Trans.State
 import Control.Monad.Zip
 import Control.Lens as Lens
 import Data.Binary as Binary
@@ -118,10 +121,11 @@
 import Linear.Vector
 #if (MIN_VERSION_transformers(0,5,0)) || !(MIN_VERSION_transformers(0,4,0))
 import Prelude as P
-#if __GLASGOW_HASKELL__ < 710
-import Data.Monoid
 #endif
+#if !(MIN_VERSION_base(4,11,0))
+import Data.Semigroup
 #endif
+import System.Random
 
 
 #ifdef HLINT
@@ -176,6 +180,10 @@
   {-# INLINE reflectDim #-}
 #endif
 
+instance (Dim n, Random a) => Random (V n a) where
+  random = runState (V <$> V.replicateM (reflectDim (Proxy :: Proxy n)) (state random))
+  randomR (V ls,V hs) = runState (V <$> V.zipWithM (\l h -> state $ randomR (l,h)) ls hs)
+
 data ReifiedDim (s :: *)
 
 retagDim :: (Proxy s -> a) -> proxy (ReifiedDim s) -> a
@@ -211,6 +219,15 @@
 instance Dim n => Dim (V n a) where
   reflectDim _ = reflectDim (Proxy :: Proxy n)
   {-# INLINE reflectDim #-}
+
+instance (Dim n, Semigroup a) => Semigroup (V n a) where
+ (<>) = liftA2 (<>)
+
+instance (Dim n, Monoid a) => Monoid (V n a) where
+  mempty = pure mempty
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = liftA2 mappend
+#endif
 
 instance Functor (V n) where
   fmap f (V as) = V (fmap f as)
diff --git a/src/Linear/V0.hs b/src/Linear/V0.hs
--- a/src/Linear/V0.hs
+++ b/src/Linear/V0.hs
@@ -30,6 +30,9 @@
 #define MIN_VERSION_transformers(x,y,z) 1
 #endif
 
+#ifndef MIN_VERSION_base
+#define MIN_VERSION_base(x,y,z) 1
+#endif
 
 -----------------------------------------------------------------------------
 -- |
@@ -90,6 +93,7 @@
 #if __GLASGOW_HASKELL__ >= 707
 import Linear.V
 #endif
+import System.Random
 import Prelude hiding (sum)
 
 -- $setup
@@ -122,6 +126,14 @@
   fromV _ = V0
 #endif
 
+instance Random (V0 a) where
+  random g = (V0, g)
+  randomR _ g = (V0, g)
+  randomRs _ _ = repeat V0
+  randoms _ = repeat V0
+  randomRIO _ = pure V0
+  randomIO = pure V0
+
 instance Serial1 V0 where
   serializeWith _ = serialize
   deserializeWith _ = deserialize
@@ -165,6 +177,15 @@
   {-# INLINE pure #-}
   V0 <*> V0 = V0
   {-# INLINE (<*>) #-}
+
+instance Semigroup (V0 a) where
+  _ <> _ = V0
+
+instance Monoid (V0 a) where
+  mempty = V0
+#if !(MIN_VERSION_base(4,11,0))
+  mappend _ _ = V0
+#endif
 
 instance Additive V0 where
   zero = V0
diff --git a/src/Linear/V1.hs b/src/Linear/V1.hs
--- a/src/Linear/V1.hs
+++ b/src/Linear/V1.hs
@@ -32,6 +32,10 @@
 #define MIN_VERSION_transformers(x,y,z) 1
 #endif
 
+#ifndef MIN_VERSION_base
+#define MIN_VERSION_base(x,y,z) 1
+#endif
+
 -----------------------------------------------------------------------------
 -- |
 -- Copyright   :  (C) 2012-2015 Edward Kmett
@@ -88,6 +92,10 @@
 import Linear.Epsilon
 import Linear.Vector
 import Prelude hiding (sum)
+import System.Random
+#if !(MIN_VERSION_base(4,11,0))
+import Data.Semigroup
+#endif
 
 import qualified Data.Vector.Generic.Mutable as M
 import qualified Data.Vector.Generic as G
@@ -392,6 +400,13 @@
   put = serializeWith Cereal.put
   get = deserializeWith Cereal.get
 
+instance Random a => Random (V1 a) where
+  random g = case random g of (a, g') -> (V1 a, g')
+  randoms g = V1 <$> randoms g
+  randomR (V1 a, V1 b) g = case randomR (a, b) g of (a', g') -> (V1 a', g')
+  randomRs (V1 a, V1 b) g = V1 <$> randomRs (a, b) g
+  randomIO = V1 <$> randomIO
+
 #if (MIN_VERSION_transformers(0,5,0)) || !(MIN_VERSION_transformers(0,4,0))
 instance Eq1 V1 where
   liftEq f (V1 a) (V1 b) = f a b
@@ -410,3 +425,13 @@
 
 instance Field1 (V1 a) (V1 b) a b where
   _1 f (V1 x) = V1 <$> f x
+
+instance Semigroup a => Semigroup (V1 a) where
+ (<>) = liftA2 (<>)
+
+instance Monoid a => Monoid (V1 a) where
+  mempty = pure mempty
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = liftA2 mappend
+#endif
+
diff --git a/src/Linear/V2.hs b/src/Linear/V2.hs
--- a/src/Linear/V2.hs
+++ b/src/Linear/V2.hs
@@ -29,6 +29,9 @@
 #define MIN_VERSION_transformers(x,y,z) 1
 #endif
 
+#ifndef MIN_VERSION_base
+#define MIN_VERSION_base(x,y,z) 1
+#endif
 
 -----------------------------------------------------------------------------
 -- |
@@ -49,6 +52,7 @@
   , ex, ey
   , perp
   , angle
+  , unangle
   , crossZ
   ) where
 
@@ -99,6 +103,7 @@
 import Linear.Vector
 import Linear.V1 (R1(..),ex)
 import Prelude hiding (sum)
+import System.Random
 
 -- $setup
 -- >>> import Control.Lens
@@ -137,6 +142,16 @@
   fromV (V v) = V2 (v V.! 0) (v V.! 1)
 #endif
 
+instance Random a => Random (V2 a) where
+  random g = case random g of
+   (a, g') -> case random g' of
+     (b, g'') -> (V2 a b, g'')
+  {-# inline random #-}
+  randomR (V2 a b, V2 c d) g = case randomR (a, c) g of
+    (x, g') -> case randomR (b, d) g' of
+      (y, g'') -> (V2 x y, g'')
+  {-# inline randomR #-}
+
 instance Functor V2 where
   fmap f (V2 a b) = V2 (f a) (f b)
   {-# INLINE fmap #-}
@@ -434,6 +449,13 @@
 angle :: Floating a => a -> V2 a
 angle a = V2 (cos a) (sin a)
 
+unangle :: (Floating a, Ord a) => V2 a -> a
+unangle a@(V2 ax ay) =
+  let alpha = asin $ ay / norm a
+  in if ax < 0
+       then pi - alpha
+       else alpha
+
 -- | The Z-component of the cross product of two vectors in the XY-plane.
 --
 -- >>> crossZ (V2 1 0) (V2 0 1)
@@ -488,3 +510,12 @@
 
 instance Field2 (V2 a) (V2 a) a a where
   _2 f (V2 x y) = f y <&> \y' -> V2 x y'
+
+instance Semigroup a => Semigroup (V2 a) where
+ (<>) = liftA2 (<>)
+
+instance Monoid a => Monoid (V2 a) where
+  mempty = pure mempty
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = liftA2 mappend
+#endif
diff --git a/src/Linear/V3.hs b/src/Linear/V3.hs
--- a/src/Linear/V3.hs
+++ b/src/Linear/V3.hs
@@ -100,6 +100,7 @@
 #endif
 import Linear.V2
 import Linear.Vector
+import System.Random
 
 {-# ANN module "HLint: ignore Reduce duplication" #-}
 
@@ -137,6 +138,16 @@
   length _ = 3
 #endif
 
+instance Random a => Random (V3 a) where
+  random g = case random g of
+    (a, g') -> case random g' of
+      (b, g'') -> case random g'' of
+        (c, g''') -> (V3 a b c, g''')
+  randomR (V3 a b c, V3 a' b' c') g = case randomR (a,a') g of
+    (a'', g') -> case randomR (b,b') g' of
+      (b'', g'') -> case randomR (c,c') g'' of
+        (c'', g''') -> (V3 a'' b'' c'', g''')
+
 instance Traversable V3 where
   traverse f (V3 a b c) = V3 <$> f a <*> f b <*> f c
   {-# INLINE traverse #-}
@@ -513,3 +524,13 @@
 
 instance Field3 (V3 a) (V3 a) a a where
   _3 f (V3 x y z) = f z <&> \z' -> V3 x y z'
+
+instance Semigroup a => Semigroup (V3 a) where
+ (<>) = liftA2 (<>)
+
+instance Monoid a => Monoid (V3 a) where
+  mempty = pure mempty
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = liftA2 mappend
+#endif
+
diff --git a/src/Linear/V4.hs b/src/Linear/V4.hs
--- a/src/Linear/V4.hs
+++ b/src/Linear/V4.hs
@@ -107,6 +107,7 @@
 import Linear.V2
 import Linear.V3
 import Linear.Vector
+import System.Random
 
 {-# ANN module "HLint: ignore Reduce duplication" #-}
 
@@ -144,6 +145,18 @@
   length _ = 4
 #endif
 
+instance Random a => Random (V4 a) where
+  random g = case random g of
+    (a, g') -> case random g' of
+      (b, g'') -> case random g'' of
+        (c, g''') -> case random g''' of
+          (d, g'''') -> (V4 a b c d, g'''')
+  randomR (V4 a b c d, V4 a' b' c' d') g = case randomR (a,a') g of
+    (a'', g') -> case randomR (b,b') g' of
+      (b'', g'') -> case randomR (c,c') g'' of
+        (c'', g''') -> case randomR (d,d') g''' of
+          (d'', g'''') -> (V4 a'' b'' c'' d'', g'''')
+
 instance Traversable V4 where
   traverse f (V4 a b c d) = V4 <$> f a <*> f b <*> f c <*> f d
   {-# INLINE traverse #-}
@@ -654,3 +667,13 @@
 
 instance Field4 (V4 a) (V4 a) a a where
   _4 f (V4 x y z w) = f w <&> \w' -> V4 x y z w'
+
+instance Semigroup a => Semigroup (V4 a) where
+ (<>) = liftA2 (<>)
+
+instance Monoid a => Monoid (V4 a) where
+  mempty = pure mempty
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = liftA2 mappend
+#endif
+
diff --git a/src/Linear/Vector.hs b/src/Linear/Vector.hs
--- a/src/Linear/Vector.hs
+++ b/src/Linear/Vector.hs
@@ -2,6 +2,7 @@
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 #if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
 {-# LANGUAGE Trustworthy #-}
@@ -41,6 +42,8 @@
 #else
 import Data.Foldable as Foldable (forM_, foldl')
 #endif
+import Data.Functor.Compose
+import Data.Functor.Product
 import Data.HashMap.Lazy as HashMap
 import Data.Hashable
 import Data.IntMap as IntMap
@@ -88,6 +91,14 @@
   gliftI2 f (a :*: b) (c :*: d) = gliftI2 f a c :*: gliftI2 f b d
   {-# INLINE gliftI2 #-}
 
+instance (Additive f, GAdditive g) => GAdditive (f :.: g) where
+  gzero = Comp1 $ gzero <$ (zero :: f Int)
+  {-# INLINE gzero #-}
+  gliftU2 f (Comp1 a) (Comp1 b) = Comp1 $ liftU2 (gliftU2 f) a b
+  {-# INLINE gliftU2 #-}
+  gliftI2 f (Comp1 a) (Comp1 b) = Comp1 $ liftI2 (gliftI2 f) a b
+  {-# INLINE gliftI2 #-}
+
 instance Additive f => GAdditive (Rec1 f) where
   gzero = Rec1 zero
   {-# INLINE gzero #-}
@@ -112,7 +123,6 @@
   {-# INLINE gliftI2 #-}
 #endif
 
-
 -- | A vector is an additive group with additional structure.
 class Functor f => Additive f where
   -- | The zero vector
@@ -171,6 +181,26 @@
   {-# INLINE liftI2 #-}
 #endif
 #endif
+
+instance (Additive f, Additive g) => Additive (Product f g) where
+  zero = Pair zero zero
+  liftU2 f (Pair a b) (Pair c d) = Pair (liftU2 f a c) (liftU2 f b d)
+  liftI2 f (Pair a b) (Pair c d) = Pair (liftI2 f a c) (liftI2 f b d)
+  Pair a b ^+^ Pair c d = Pair (a ^+^ c) (b ^+^ d)
+  Pair a b ^-^ Pair c d = Pair (a ^-^ c) (b ^-^ d)
+  lerp alpha (Pair a b) (Pair c d) = Pair (lerp alpha a c) (lerp alpha b d)
+
+instance (Additive f, Additive g) => Additive (Compose f g) where
+  zero = Compose $ zero <$ (zero :: f Int)
+  {-# INLINE zero #-}
+  Compose a ^+^ Compose b = Compose $ liftU2 (^+^) a b
+  {-# INLINE (^+^) #-}
+  Compose a ^-^ Compose b = Compose $ liftU2 (^-^) a b
+  {-# INLINE (^-^) #-}
+  liftU2 f (Compose a) (Compose b) = Compose $ liftU2 (liftU2 f) a b
+  {-# INLINE liftU2 #-}
+  liftI2 f (Compose a) (Compose b) = Compose $ liftI2 (liftI2 f) a b
+  {-# INLINE liftI2 #-}
 
 instance Additive ZipList where
   zero = ZipList []