diff --git a/.travis.yml b/.travis.yml
--- a/.travis.yml
+++ b/.travis.yml
@@ -3,17 +3,18 @@
 # Basic template to test with GHC=[7.8, 7.10, 8.0]. Caches the 'stack' build
 # directory to speed subsequent compilations.
 #
-language: generic
+# https://docs.travis-ci.com/user/reference/overview/
+#
+language: minimal
 dist: trusty
-sudo: false
+sudo: true      # more storage available on sudo-enabled instances
 
 cache:
   directories:
-    - $HOME/.stack
+    - $HOME/.stack/snapshots
     - .stack-work/install
 
 before_cache:
-  - rm -rf $(stack path --local-install-root)/bin
   - rm -rf $(stack path --local-install-root)/doc
 
 addons:
@@ -28,10 +29,6 @@
 
 matrix:
   include:
-    - env: GHC=7.8.4 CABAL=1.24
-      compiler: "GHC 7.8"
-      addons: { apt: { sources: [ *apt_sources ], packages: [ *apt_packages, cabal-install-1.24 ] }}
-
     - env: GHC=7.10.3 CABAL=1.24
       compiler: "GHC 7.10"
       addons: { apt: { sources: [ *apt_sources ], packages: [ *apt_packages, cabal-install-1.24 ] }}
@@ -40,10 +37,14 @@
       compiler: "GHC 8.0"
       addons: { apt: { sources: [ *apt_sources ], packages: [ *apt_packages, cabal-install-1.24 ] }}
 
-    - env: GHC=8.2.0 CABAL=1.24
+    - env: GHC=8.2.2 CABAL=2.0
       compiler: "GHC 8.2"
-      addons: { apt: { sources: [ *apt_sources ], packages: [ *apt_packages, cabal-install-1.24 ] }}
+      addons: { apt: { sources: [ *apt_sources ], packages: [ *apt_packages, cabal-install-2.0 ] }}
 
+    - env: GHC=8.4.1 CABAL=2.0
+      compiler: "GHC 8.4"
+      addons: { apt: { sources: [ *apt_sources ], packages: [ *apt_packages, cabal-install-2.0 ] }}
+
     # - env: GHC=head CABAL=head
     #   compiler: "GHC HEAD"
     #   addons: { apt: { sources: [ *apt_sources ], packages: [ *apt_packages, ghc-head, cabal-install-head ] }}
@@ -61,21 +62,12 @@
   - echo "$(ghc --version) [$(ghc --print-project-git-commit-id 2> /dev/null || echo '?')]"
   - cabal --version
   - stack --version
-  - |
-      if [ ${GHC} == head -o ! -e stack.yaml ]; then
-        travis_retry cabal install --only-dependencies
-      else
-        travis_retry stack build --only-dependencies --no-terminal --no-haddock-deps
-      fi
+  - stack build --no-terminal --dry-run
+  - travis_retry travis_wait 60 stack build --fast --test --only-dependencies --no-terminal --no-copy-bins
 
 script:
-  - |
-      if [ ${GHC} == head -o ! -e stack.yaml ]; then
-        cabal install
-        cabal haddock
-      else
-        stack test --no-terminal --haddock --no-haddock-deps
-      fi
+  - stack build --no-terminal --fast --test --haddock --no-haddock-deps --no-run-tests
+  - stack test
 
 after_success:
   - source .travis/update-accelerate-buildbot.sh
diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -6,39 +6,46 @@
 project adheres to the [Haskell Package Versioning Policy
 (PVP)](https://pvp.haskell.org)
 
-
-## [0.5.0.1]
-
-* Fix [#10]: doctest failure with accelerate-1.1.1.0
-
-## [0.5]
-
-* Fix [#25][acc-llvm#25]: "impossible evaluation" error with accelerate-llvm-native
-
-## [0.4.1]
+## [0.6.0.0] - 2018-04-03
+### Changed
+  * update for accelerate-1.2
 
-* Support `doctest-0.12`.
+### Fixed
+  * Disambiguate uses of `fromInteger` and `fromRational`.
 
-## [0.4]
+## [0.5.0.1] - 2017-09-29
+### Fixed
+  * doctest failure with accelerate-1.1.1.0 ([#10])
 
-* Revamp `Setup.hs` to use `cabal-doctest`. This makes `linear-accelerate`
-  build with `Cabal-2.0`.
-* Fix handedness of `(-^)` and `(^-)`.
-* Re-enable the `doctest`s.
+## [0.5] - 2017-08-06
+### Fixed
+  * "impossible evaluation" error with accelerate-llvm-native ([#25][acc-llvm#25])
 
-## [0.3]
+## [0.4.1] - 2017-07-29
+### Added
+  * Support `doctest-0.12`.
 
-* Expanded API coverage
+## [0.4] - 2017-06-16
+### Changed
+  * Revamp `Setup.hs` to use `cabal-doctest`. This makes `linear-accelerate` build with `Cabal-2.0`.
+  * Re-enable the `doctest`s.
 
-## [0.2]
+### Fixed
+  * Fix handedness of `(-^)` and `(^-)`.
 
-* Fix [#1]: Move `Data.Complex` instances into the `accelerate` package
+## [0.3] - 2017-04-01
+### Changed
+  * Expanded API coverage
 
-## [0.1]
+## [0.2] - 2014-09-15
+### Fixed
+  * Move `Data.Complex` instances into the `accelerate` package ([#1])
 
-* Repository initialized
+## [0.1] - 2014-03-28
+  * Repository initialized
 
 
+[0.6.0.0]:          https://github.com/ekmett/linear-accelerate/compare/v0.5.0.1...v0.6.0.0
 [0.5.0.1]:          https://github.com/ekmett/linear-accelerate/compare/v0.5...v0.5.0.1
 [0.5]:              https://github.com/ekmett/linear-accelerate/compare/v0.4.1...v0.5
 [0.4.1]:            https://github.com/ekmett/linear-accelerate/compare/v0.4...v0.4.1
diff --git a/linear-accelerate.cabal b/linear-accelerate.cabal
--- a/linear-accelerate.cabal
+++ b/linear-accelerate.cabal
@@ -1,6 +1,6 @@
 name:          linear-accelerate
 category:      Math, Algebra, Compilers/Interpreters, Concurrency, Data, Parallelism
-version:       0.5.0.1
+version:       0.6.0.0
 license:       BSD3
 cabal-version: >= 1.10
 license-file:  LICENSE
@@ -21,12 +21,12 @@
   Warning.hs
 
 source-repository head
-  type: git
+  type:     git
   location: git://github.com/ekmett/linear-accelerate.git
 
 source-repository this
-  type: git
-  tag: v0.5.0.1
+  type:     git
+  tag:      v0.6.0.0
   location: git://github.com/ekmett/linear-accelerate.git
 
 custom-setup
@@ -38,15 +38,17 @@
 library
   build-depends:
       base              >= 4.5          && < 5
-    , accelerate        >= 0.16
+    , accelerate        >= 1.2
     , distributive      >= 0.2.2        && < 1
     , lens              >= 4            && < 5
     , linear            >= 1.10         && < 2
 
   exposed-modules:
     Data.Array.Accelerate.Linear
+    Data.Array.Accelerate.Linear.Conjugate
     Data.Array.Accelerate.Linear.Epsilon
     Data.Array.Accelerate.Linear.Matrix
+    Data.Array.Accelerate.Linear.Trace
     Data.Array.Accelerate.Linear.Metric
     Data.Array.Accelerate.Linear.Plucker
     Data.Array.Accelerate.Linear.Quaternion
@@ -61,19 +63,32 @@
   other-modules:
     Data.Array.Accelerate.Linear.Lift
 
-  default-language: Haskell2010
-  ghc-options: -Wall -fwarn-tabs -O2
-  hs-source-dirs: src
+  hs-source-dirs:
+    src
 
+  default-language:
+    Haskell2010
+
+  ghc-options:
+    -Wall
+    -fwarn-tabs
+    -O2
+
 test-suite doctests
-  type:    exitcode-stdio-1.0
-  main-is: doctests.hs
-  default-language: Haskell2010
+  type:               exitcode-stdio-1.0
+  hs-source-dirs:     tests
+  main-is:            doctests.hs
+
   build-depends:
-    base    == 4.*,
-    doctest >= 0.11.1 && < 0.14,
-    linear-accelerate
-  ghc-options: -Wall -threaded
-  hs-source-dirs: tests
+      base              == 4.*
+    , doctest           >= 0.11.1 && < 0.16
+    , linear-accelerate
+
+  default-language:
+    Haskell2010
+
+  ghc-options:
+    -Wall
+    -threaded
 
 -- vim: nospell
diff --git a/src/Data/Array/Accelerate/Linear.hs b/src/Data/Array/Accelerate/Linear.hs
--- a/src/Data/Array/Accelerate/Linear.hs
+++ b/src/Data/Array/Accelerate/Linear.hs
@@ -2,7 +2,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
diff --git a/src/Data/Array/Accelerate/Linear/Conjugate.hs b/src/Data/Array/Accelerate/Linear/Conjugate.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Linear/Conjugate.hs
@@ -0,0 +1,63 @@
+{-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE FlexibleInstances   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE ViewPatterns        #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      : Data.Array.Accelerate.Linear.Conjugate
+-- Copyright   : 2018 Trevor L. McDonell
+-- License     : BSD-style (see the file LICENSE)
+--
+-- Maintainer  : Edward Kmett <ekmett@gmail.com>
+-- Stability   : experimental
+-- Portability : non-portable
+--
+-- Involutive rings
+----------------------------------------------------------------------------
+
+module Data.Array.Accelerate.Linear.Conjugate (
+
+  Conjugate(..),
+  TrivialConjugate,
+
+) where
+
+import Data.Array.Accelerate                                        as A
+import Data.Array.Accelerate.Data.Complex                           ( Complex(..) )
+
+import Linear.Conjugate
+
+instance Conjugate (Exp Int)
+instance Conjugate (Exp Int64)
+instance Conjugate (Exp Int32)
+instance Conjugate (Exp Int16)
+instance Conjugate (Exp Int8)
+instance Conjugate (Exp Word)
+instance Conjugate (Exp Word64)
+instance Conjugate (Exp Word32)
+instance Conjugate (Exp Word16)
+instance Conjugate (Exp Word8)
+instance Conjugate (Exp Double)
+instance Conjugate (Exp Float)
+instance Conjugate (Exp CFloat)
+instance Conjugate (Exp CDouble)
+
+instance (Conjugate (Exp a), A.RealFloat a, Elt (Complex a)) => Conjugate (Exp (Complex a)) where
+  conjugate (unlift -> a :+ b :: Complex (Exp a)) = lift $ conjugate a :+ negate b
+
+instance TrivialConjugate (Exp Int)
+instance TrivialConjugate (Exp Int64)
+instance TrivialConjugate (Exp Int32)
+instance TrivialConjugate (Exp Int16)
+instance TrivialConjugate (Exp Int8)
+instance TrivialConjugate (Exp Word)
+instance TrivialConjugate (Exp Word64)
+instance TrivialConjugate (Exp Word32)
+instance TrivialConjugate (Exp Word16)
+instance TrivialConjugate (Exp Word8)
+instance TrivialConjugate (Exp Double)
+instance TrivialConjugate (Exp Float)
+instance TrivialConjugate (Exp CFloat)
+instance TrivialConjugate (Exp CDouble)
+
diff --git a/src/Data/Array/Accelerate/Linear/Epsilon.hs b/src/Data/Array/Accelerate/Linear/Epsilon.hs
--- a/src/Data/Array/Accelerate/Linear/Epsilon.hs
+++ b/src/Data/Array/Accelerate/Linear/Epsilon.hs
@@ -5,7 +5,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Epsilon
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
diff --git a/src/Data/Array/Accelerate/Linear/Lift.hs b/src/Data/Array/Accelerate/Linear/Lift.hs
--- a/src/Data/Array/Accelerate/Linear/Lift.hs
+++ b/src/Data/Array/Accelerate/Linear/Lift.hs
@@ -5,7 +5,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Lift
--- Copyright   : 2015 Trevor L. McDonell
+-- Copyright   : [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
diff --git a/src/Data/Array/Accelerate/Linear/Matrix.hs b/src/Data/Array/Accelerate/Linear/Matrix.hs
--- a/src/Data/Array/Accelerate/Linear/Matrix.hs
+++ b/src/Data/Array/Accelerate/Linear/Matrix.hs
@@ -6,7 +6,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Matrix
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -18,26 +18,27 @@
 
 module Data.Array.Accelerate.Linear.Matrix (
 
-  M22, M23, M24, M32, M33, M34, M42, M43, M44,
   (!*!), (!+!), (!-!), (!*), (*!), (!!*), (*!!), (!!/),
-  transpose,
+  M22, M23, M24, M32, M33, M34, M42, M43, M44,
+  m33_to_m44, m43_to_m44,
+  det22, det33, det44,
+  inv22, inv33, inv44,
   identity,
+  transpose,
   Trace(..),
+  fromQuaternion,
+  mkTransformation,
+  mkTransformationMat,
 
 ) where
 
 import Data.Array.Accelerate                    as A hiding ( transpose )
-import Data.Array.Accelerate.Data.Complex
 
 import Data.Array.Accelerate.Linear.Lift
-import Data.Array.Accelerate.Linear.Plucker
-import Data.Array.Accelerate.Linear.Quaternion
+import Data.Array.Accelerate.Linear.Trace
 import Data.Array.Accelerate.Linear.Type
-import Data.Array.Accelerate.Linear.V0
-import Data.Array.Accelerate.Linear.V1
-import Data.Array.Accelerate.Linear.V2
+import Data.Array.Accelerate.Linear.Quaternion
 import Data.Array.Accelerate.Linear.V3
-import Data.Array.Accelerate.Linear.V4
 import Data.Array.Accelerate.Linear.Vector
 
 import Linear.Matrix                            ( M22, M23, M24, M32, M33, M34, M42, M43, M44 )
@@ -51,6 +52,9 @@
 
 -- $setup
 -- >>> import Data.Array.Accelerate.Interpreter
+-- >>> import Data.Array.Accelerate.Linear.V2
+-- >>> import Data.Array.Accelerate.Linear.V3
+-- >>> import Data.Array.Accelerate.Linear.V4
 -- >>> :{
 --   let test :: Elt e => Exp e -> e
 --       test e = indexArray (run (unit e)) Z
@@ -185,21 +189,74 @@
 transpose = lift . L.transpose . unlift'
 
 
-class L.Trace m => Trace m where
-  -- | Compute the trace of a matrix
-  trace :: (A.Num a, Box2 m m a) => Exp (m (m a)) -> Exp a
-  trace = lift . L.trace . unlift'
+-- | Build a rotation matrix from a unit 'Quaternion'
+--
+fromQuaternion :: forall a. A.Num a => Exp (Quaternion a) -> Exp (M33 a)
+fromQuaternion = lift1 (L.fromQuaternion :: Quaternion (Exp a) -> M33 (Exp a))
 
-  -- | Compute the diagonal of a matrix
-  diagonal :: Box2 m m a => Exp (m (m a)) -> Exp (m a)
-  diagonal = lift . L.diagonal . unlift'
+-- | Build a transformation matrix from a rotation expressed as a 'Quaternion'
+-- and a translation vector.
+--
+mkTransformation :: forall a. A.Num a => Exp (Quaternion a) -> Exp (V3 a) -> Exp (M44 a)
+mkTransformation = lift2 (L.mkTransformation :: Quaternion (Exp a) -> V3 (Exp a) -> M44 (Exp a))
 
-instance Trace Complex
-instance Trace V0
-instance Trace V1
-instance Trace V2
-instance Trace V3
-instance Trace V4
-instance Trace Plucker
-instance Trace Quaternion
+-- | Build a transformation matrix from a rotation matrix and a translation
+-- vector.
+--
+mkTransformationMat :: A.Num a => Exp (M33 a) -> Exp (V3 a) -> Exp (M44 a)
+mkTransformationMat m v =
+  let r = L.mkTransformationMat (unlift' m) (unlift v)
+  in  lift r
+
+-- | Convert a 4x3 matrix to a 4x4 matrix, extending it with @[ 0 0 0 1 ]@
+-- column vector
+--
+m43_to_m44 :: A.Num a => Exp (M43 a) -> Exp (M44 a)
+m43_to_m44 m43 =
+  let m44 = L.m43_to_m44 (unlift' m43)
+  in  lift m44
+
+-- | Convert a 3x3 matrix to a 4x4 matrix extending it with zeros in the new row
+-- and column.
+--
+m33_to_m44 :: A.Num a => Exp (M33 a) -> Exp (M44 a)
+m33_to_m44 m33 =
+  let m44 = L.m33_to_m44 (unlift' m33)
+  in  lift m44
+
+-- | 2x2 matrix determinant
+--
+det22 :: A.Num a => Exp (M22 a) -> Exp a
+det22 = L.det22 . unlift'
+
+-- | 3x3 matrix determinant
+--
+det33 :: A.Num a => Exp (M33 a) -> Exp a
+det33 = L.det33 . unlift'
+
+-- | 4x4 matrix determinant
+--
+det44 :: A.Num a => Exp (M44 a) -> Exp a
+det44 = L.det44 . unlift'
+
+-- | 2x2 matrix inverse
+--
+inv22 :: A.Fractional a => Exp (M22 a) -> Exp (M22 a)
+inv22 m =
+  let r = L.inv22 (unlift' m)
+  in  lift r
+
+-- | 3x3 matrix inverse
+--
+inv33 :: A.Fractional a => Exp (M33 a) -> Exp (M33 a)
+inv33 m =
+  let r = L.inv33 (unlift' m)
+  in  lift r
+
+-- | 4x4 matrix inverse
+--
+inv44 :: A.Fractional a => Exp (M44 a) -> Exp (M44 a)
+inv44 m =
+  let r = L.inv44 (unlift' m)
+  in  lift r
 
diff --git a/src/Data/Array/Accelerate/Linear/Metric.hs b/src/Data/Array/Accelerate/Linear/Metric.hs
--- a/src/Data/Array/Accelerate/Linear/Metric.hs
+++ b/src/Data/Array/Accelerate/Linear/Metric.hs
@@ -7,7 +7,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Metric
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
diff --git a/src/Data/Array/Accelerate/Linear/Plucker.hs b/src/Data/Array/Accelerate/Linear/Plucker.hs
--- a/src/Data/Array/Accelerate/Linear/Plucker.hs
+++ b/src/Data/Array/Accelerate/Linear/Plucker.hs
@@ -3,16 +3,19 @@
 {-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE FlexibleInstances     #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE RebindableSyntax      #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE StandaloneDeriving    #-}
 {-# LANGUAGE TypeFamilies          #-}
 {-# LANGUAGE UndecidableInstances  #-}
-{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# LANGUAGE ViewPatterns          #-}
+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
+{-# OPTIONS_GHC -fno-warn-orphans        #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Plucker
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -24,25 +27,52 @@
 
 module Data.Array.Accelerate.Linear.Plucker (
 
-  Plucker(..), squaredError, (><), plucker, plucker3D,
+  Plucker(..),
 
+  squaredError,
+  isotropic,
+  (><),
+  plucker,
+  plucker3D,
+
+  -- * operations on lines
+  LinePass(..),
+  parallel,
+  intersects,
+  passes,
+  quadranceToOrigin,
+  closestToOrigin,
+  isLine,
+
+  -- * Basis elements
+       p01, p02, p03,
+  p10,      p12, p13,
+  p20, p21,      p23,
+  p30, p31, p32,
+
 ) where
 
-import Data.Array.Accelerate                    as A
+import Data.Array.Accelerate                    hiding ( fromInteger )
+import Data.Array.Accelerate.Data.Functor
 import Data.Array.Accelerate.Smart
 import Data.Array.Accelerate.Product
 import Data.Array.Accelerate.Array.Sugar
 
+import Data.Array.Accelerate.Linear.Epsilon
 import Data.Array.Accelerate.Linear.Lift
 import Data.Array.Accelerate.Linear.Metric
-import Data.Array.Accelerate.Linear.Vector
+import Data.Array.Accelerate.Linear.V2
 import Data.Array.Accelerate.Linear.V3
 import Data.Array.Accelerate.Linear.V4
+import Data.Array.Accelerate.Linear.Vector
 
-import Data.Typeable
 import Control.Lens
-import Linear.Plucker                           ( Plucker(..) )
+import Data.Function
+import Data.Typeable
+import Linear.Plucker                           ( Plucker(..), LinePass(..) )
+import Prelude                                  ( fromInteger )
 import qualified Linear.Plucker                 as L
+import qualified Prelude                        as P
 
 
 -- | Valid Plücker coordinates @p@ will have @'squaredError' p '==' 0@
@@ -50,30 +80,136 @@
 -- That said, floating point makes a mockery of this claim, so you may want to
 -- use 'nearZero'.
 --
-squaredError :: forall a. A.Num a => Exp (Plucker a) -> Exp a
+squaredError :: forall a. Num a => Exp (Plucker a) -> Exp a
 squaredError = lift1 (L.squaredError :: Plucker (Exp a) -> Exp a)
 
 -- | This isn't the actual metric because this bilinear form gives rise to an
--- isotropic quadratic space
+-- isotropic quadratic space.
 --
 infixl 5 ><
-(><) :: forall a. A.Num a => Exp (Plucker a) -> Exp (Plucker a) -> Exp a
+(><) :: forall a. Num a => Exp (Plucker a) -> Exp (Plucker a) -> Exp a
 (><) = lift2 ((L.><) :: Plucker (Exp a) -> Plucker (Exp a) -> Exp a)
 
+-- | Checks if the line is near-isotropic (isotropic vectors in this quadratic
+-- space represent lines in real 3D space).
+--
+isotropic :: Epsilon a => Exp (Plucker a) -> Exp Bool
+isotropic a = nearZero (a >< a)
+
 -- | Given a pair of points represented by homogeneous coordinates generate
 -- Plücker coordinates for the line through them, directed from the second
 -- towards the first.
 --
-plucker :: forall a. A.Num a => Exp (V4 a) -> Exp (V4 a) -> Exp (Plucker a)
+plucker :: forall a. Num a => Exp (V4 a) -> Exp (V4 a) -> Exp (Plucker a)
 plucker = lift2 (L.plucker :: V4 (Exp a) -> V4 (Exp a) -> Plucker (Exp a))
 
 -- | Given a pair of 3D points, generate Plücker coordinates for the line
 -- through them, directed from the second towards the first.
 --
-plucker3D :: forall a. A.Num a => Exp (V3 a) -> Exp (V3 a) -> Exp (Plucker a)
+plucker3D :: forall a. Num a => Exp (V3 a) -> Exp (V3 a) -> Exp (Plucker a)
 plucker3D = lift2 (L.plucker3D :: V3 (Exp a) -> V3 (Exp a) -> Plucker (Exp a))
 
+-- | Checks if two lines intersect (or nearly intersect).
+--
+intersects :: (Epsilon a, Ord a) => Exp (Plucker a) -> Exp (Plucker a) -> Exp Bool
+intersects a b = not (a `parallel` b) && passes a b == constant Coplanar
 
+-- | Check how two lines pass each other. @passes l1 l2@ describes @l2@ when
+-- looking down @l1@.
+--
+passes :: (Epsilon a, Ord a) => Exp (Plucker a) -> Exp (Plucker a) -> Exp LinePass
+passes a b =
+  if nearZero s then constant Coplanar else
+  if s > 0      then constant Counterclockwise
+                else constant Clockwise
+  where
+    s        = (lift u1 `dot` lift v2) + (lift u2 `dot` lift v1)
+    V2 u1 v1 = toUV a
+    V2 u2 v2 = toUV b
+
+-- | Checks if two lines are parallel.
+--
+parallel :: Epsilon a => Exp (Plucker a) -> Exp (Plucker a) -> Exp Bool
+parallel a b = nearZero $ lift u1 `cross` lift u2
+  where
+    V2 u1 _ = toUV a
+    V2 u2 _ = toUV b
+
+-- | Represent a Plücker coordinate as a pair of 3-tuples, typically denoted
+-- U and V.
+--
+toUV :: Elt a => Exp (Plucker a) -> V2 (V3 (Exp a))
+toUV (unlift -> Plucker a b c d e f) = V2 (V3 a b c) (V3 d e f)
+
+-- | The minimum squared distance of a line from the origin.
+--
+quadranceToOrigin :: Fractional a => Exp (Plucker a) -> Exp a
+quadranceToOrigin p = (lift v `dot` lift v) / (lift u `dot` lift u)
+  where
+    V2 u v = toUV p
+
+-- | The point where a line is closest to the origin.
+--
+closestToOrigin :: Fractional a => Exp (Plucker a) -> Exp (V3 a)
+closestToOrigin p = normalizePoint . lift $ V4 x y z (lift u `dot` lift u)
+  where
+    V2 u v   = toUV p
+    V3 x y z = unlift $ lift v `cross` lift u
+
+-- | Not all 6-dimensional points correspond to a line in 3D. This predicate
+-- tests that a Plücker coordinate lies on the Grassmann manifold, and does
+-- indeed represent a 3D line.
+--
+isLine :: Epsilon a => Exp (Plucker a) -> Exp Bool
+isLine p = nearZero $ lift u `dot` lift v
+  where
+    V2 u v = toUV p
+
+
+-- | These elements form a basis for the Plücker space, or the Grassmanian
+-- manifold @Gr(2,V4)@.
+--
+-- @
+-- 'p01' :: 'Lens'' ('Plucker' a) a
+-- 'p02' :: 'Lens'' ('Plucker' a) a
+-- 'p03' :: 'Lens'' ('Plucker' a) a
+-- 'p23' :: 'Lens'' ('Plucker' a) a
+-- 'p31' :: 'Lens'' ('Plucker' a) a
+-- 'p12' :: 'Lens'' ('Plucker' a) a
+-- @
+--
+p01, p02, p03, p23, p31, p12 :: Elt a => Lens' (Exp (Plucker a)) (Exp a)
+p01 = liftLens (L.p01 :: Lens' (Plucker (Exp a)) (Exp a))
+p02 = liftLens (L.p02 :: Lens' (Plucker (Exp a)) (Exp a))
+p03 = liftLens (L.p03 :: Lens' (Plucker (Exp a)) (Exp a))
+p23 = liftLens (L.p23 :: Lens' (Plucker (Exp a)) (Exp a))
+p31 = liftLens (L.p31 :: Lens' (Plucker (Exp a)) (Exp a))
+p12 = liftLens (L.p12 :: Lens' (Plucker (Exp a)) (Exp a))
+
+-- | These elements form an alternate basis for the Plücker space, or the
+-- Grassmanian manifold @Gr(2,V4)@.
+--
+-- @
+-- 'p10' :: 'Num' a => 'Lens'' ('Plucker' a) a
+-- 'p20' :: 'Num' a => 'Lens'' ('Plucker' a) a
+-- 'p30' :: 'Num' a => 'Lens'' ('Plucker' a) a
+-- 'p32' :: 'Num' a => 'Lens'' ('Plucker' a) a
+-- 'p13' :: 'Num' a => 'Lens'' ('Plucker' a) a
+-- 'p21' :: 'Num' a => 'Lens'' ('Plucker' a) a
+-- @
+--
+p10, p20, p30, p32, p13, p21 :: Num a => Lens' (Exp (Plucker a)) (Exp a)
+p10 = anti p01
+p20 = anti p02
+p30 = anti p03
+p32 = anti p23
+p13 = anti p31
+p21 = anti p21
+
+anti :: (P.Functor f, Num a) => ((Exp a -> f (Exp a)) -> r) -> (Exp a -> f (Exp a)) -> r
+anti k f = k (P.fmap negate . f . negate)
+
+
 -- Instances
 -- ---------
 
@@ -120,4 +256,84 @@
 
 instance (Elt a, Elt b) => Each (Exp (Plucker a)) (Exp (Plucker b)) (Exp a) (Exp b) where
   each = liftLens (each :: Traversal (Plucker (Exp a)) (Plucker (Exp b)) (Exp a) (Exp b))
+
+instance Eq a => Eq (Plucker a) where
+  (==) = (==) `on` t6
+  (/=) = (/=) `on` t6
+
+instance Ord a => Ord (Plucker a) where
+  (<)  = (<) `on` t6
+  (>)  = (>) `on` t6
+  (<=) = (<=) `on` t6
+  (>=) = (>=) `on` t6
+  min  = pl $$ on min t6
+  max  = pl $$ on max t6
+
+t6 :: Elt a => Exp (Plucker a) -> Exp (a,a,a,a,a,a)
+t6 (unlift -> Plucker a b c d e f) = tup6 (a,b,c,d,e,f)
+
+pl :: Elt a => Exp (a,a,a,a,a,a) -> Exp (Plucker a)
+pl (untup6 -> (a,b,c,d,e,f)) = lift (Plucker a b c d e f)
+
+instance Num a => P.Num (Exp (Plucker a)) where
+  (+)           = lift2 ((+) :: Plucker (Exp a) -> Plucker (Exp a) -> Plucker (Exp a))
+  (-)           = lift2 ((-) :: Plucker (Exp a) -> Plucker (Exp a) -> Plucker (Exp a))
+  (*)           = lift2 ((*) :: Plucker (Exp a) -> Plucker (Exp a) -> Plucker (Exp a))
+  negate        = fmap negate
+  abs           = fmap abs
+  signum        = fmap signum
+  fromInteger x = lift (P.fromInteger x :: Plucker (Exp a))
+
+instance Fractional a => P.Fractional (Exp (Plucker a)) where
+  recip          = fmap recip
+  (/)            = lift2 ((/) :: Plucker (Exp a) -> Plucker (Exp a) -> Plucker (Exp a))
+  fromRational x = lift (P.fromRational x :: Plucker (Exp a))
+
+instance Floating a => P.Floating (Exp (Plucker a)) where
+  pi      = lift (pi :: Plucker (Exp a))
+  exp     = fmap exp
+  sqrt    = fmap sqrt
+  log     = fmap log
+  (**)    = lift2 ((**) :: Plucker (Exp a) -> Plucker (Exp a) -> Plucker (Exp a))
+  logBase = lift2 (logBase :: Plucker (Exp a) -> Plucker (Exp a) -> Plucker (Exp a))
+  sin     = fmap sin
+  tan     = fmap tan
+  cos     = fmap cos
+  asin    = fmap asin
+  atan    = fmap atan
+  acos    = fmap acos
+  sinh    = fmap sinh
+  tanh    = fmap tanh
+  cosh    = fmap cosh
+  asinh   = fmap asinh
+  atanh   = fmap atanh
+  acosh   = fmap acosh
+
+instance Epsilon a => Epsilon (Plucker a) where
+  nearZero = nearZero . quadrance
+
+instance Functor Plucker where
+  fmap g (unlift -> Plucker a b c d e f) = lift (Plucker (g a) (g b) (g c) (g d) (g e) (g f))
+  x <$ _                                 = lift (Plucker x x x x x x)
+
+
+type instance EltRepr LinePass = Int8
+
+instance Elt LinePass where
+  eltType _ = eltType (undefined::Int8)
+
+  toElt x = let (==) = (P.==)   -- -XRebindableSyntax hax
+            in  case x of
+                  0 -> Coplanar
+                  1 -> Clockwise
+                  2 -> Counterclockwise
+                  _ -> P.error "LinePass: unhandled constructor"
+
+  fromElt Coplanar         = 0
+  fromElt Clockwise        = 1
+  fromElt Counterclockwise = 2
+
+instance Eq LinePass where
+  x == y = bitcast x == (bitcast y :: Exp Int8)
+  x /= y = bitcast x /= (bitcast y :: Exp Int8)
 
diff --git a/src/Data/Array/Accelerate/Linear/Quaternion.hs b/src/Data/Array/Accelerate/Linear/Quaternion.hs
--- a/src/Data/Array/Accelerate/Linear/Quaternion.hs
+++ b/src/Data/Array/Accelerate/Linear/Quaternion.hs
@@ -1,15 +1,18 @@
 {-# LANGUAGE ConstraintKinds       #-}
+{-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE FlexibleInstances     #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE RebindableSyntax      #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE TypeFamilies          #-}
 {-# LANGUAGE UndecidableInstances  #-}
+{-# LANGUAGE ViewPatterns          #-}
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Quaternion
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -23,22 +26,139 @@
 
   Quaternion(..),
 
+  slerp,
+  asinq,
+  acosq,
+  atanq,
+  asinhq,
+  acoshq,
+  atanhq,
+  absi,
+  pow,
+  rotate,
+  axisAngle,
+
 ) where
 
 import Data.Array.Accelerate
-import Data.Array.Accelerate.Smart
-import Data.Array.Accelerate.Product
 import Data.Array.Accelerate.Array.Sugar
+import Data.Array.Accelerate.Data.Complex       hiding ( conjugate )
+import Data.Array.Accelerate.Data.Functor
+import Data.Array.Accelerate.Product
+import Data.Array.Accelerate.Smart
 
+import Data.Array.Accelerate.Linear.Conjugate
+import Data.Array.Accelerate.Linear.Epsilon
 import Data.Array.Accelerate.Linear.Lift
 import Data.Array.Accelerate.Linear.Metric
-import Data.Array.Accelerate.Linear.Vector
 import Data.Array.Accelerate.Linear.V3
+import Data.Array.Accelerate.Linear.Vector
 
 import Control.Lens
+import Data.Function
 import Linear.Quaternion                        ( Quaternion(..) )
+import qualified Prelude                        as P
 
 
+-- | Spherical linear interpolation between two quaternions
+--
+slerp :: RealFloat a => Exp (Quaternion a) -> Exp (Quaternion a) -> Exp a -> Exp (Quaternion a)
+slerp q p t =
+  if 1.0 - cosphi < 1.0e-8
+    then q
+    else ((sin ((1-t)*phi) *^ q) + sin (t*phi) *^ fp) ^/ sin phi
+  where
+    dqp = dot q p
+    phi = acos cosphi
+    (cosphi, fp) = unlift $ if dqp < 0 then tup2 (-dqp, negate p)
+                                       else tup2 (dqp, p)
+
+-- | 'asin' with a specified branch cut
+--
+asinq :: (RealFloat a, Elt (Complex a)) => Exp (Quaternion a) -> Exp (Quaternion a) -> Exp (Quaternion a)
+asinq q@(unlift -> Quaternion e _) u =
+  if qiq /= 0.0 || e >= -1 && e <= 1
+    then asin q
+    else cutWith (asin (lift $ e :+ sqrt qiq)) u
+  where
+    qiq = qi q
+
+-- | 'acos' with a specified branch cut
+--
+acosq :: (RealFloat a, Elt (Complex a)) => Exp (Quaternion a) -> Exp (Quaternion a) -> Exp (Quaternion a)
+acosq q@(unlift -> Quaternion e _) u =
+  if qiq /= 0.0 || e >= -1 && e <= 1
+    then acos q
+    else cutWith (acos (lift $ e :+ sqrt qiq)) u
+  where
+    qiq = qi q
+
+-- | 'atan' with a specified branch cut
+--
+atanq :: (RealFloat a, Elt (Complex a)) => Exp (Quaternion a) -> Exp (Quaternion a) -> Exp (Quaternion a)
+atanq q@(unlift -> Quaternion e _) u =
+  if e /= 0.0 || qiq >= -1 && qiq <= 1
+    then atan q
+    else cutWith (atan (lift $ e :+ sqrt qiq)) u
+  where
+    qiq = qi q
+
+-- | 'asinh' with a specified branch cut
+--
+asinhq :: (RealFloat a, Elt (Complex a)) => Exp (Quaternion a) -> Exp (Quaternion a) -> Exp (Quaternion a)
+asinhq q@(unlift -> Quaternion e _) u =
+  if e /= 0.0 || qiq >= -1 && qiq <= 1
+    then asinh q
+    else cutWith (asinh (lift $ e :+ sqrt qiq)) u
+  where
+    qiq = qi q
+
+-- | 'acosh' with a specified branch cut
+--
+acoshq :: (RealFloat a, Elt (Complex a)) => Exp (Quaternion a) -> Exp (Quaternion a) -> Exp (Quaternion a)
+acoshq q@(unlift -> Quaternion e _) u =
+  if qiq /= 0.0 || e >= 1
+    then asinh q
+    else cutWith (acosh (lift $ e :+ sqrt qiq)) u
+  where
+    qiq = qi q
+
+-- | 'atanh' with a specified branch cut
+--
+atanhq :: (RealFloat a, Elt (Complex a)) => Exp (Quaternion a) -> Exp (Quaternion a) -> Exp (Quaternion a)
+atanhq q@(unlift -> Quaternion e _) u =
+  if qiq /= 0.0 || e > -1 && e < 1
+    then atanh q
+    else cutWith (atanh (lift $ e :+ sqrt qiq)) u
+  where
+    qiq = qi q
+
+-- | norm of the imaginary component
+--
+absi :: Floating a => Exp (Quaternion a) -> Exp a
+absi = sqrt . qi
+
+-- | raise a 'Quaternion' to a scalar power
+--
+pow :: (RealFloat a, Elt (Complex a)) => Exp (Quaternion a) -> Exp a -> Exp (Quaternion a)
+pow q t = exp (t *^ log q)
+
+-- | Apply a rotation to a vector
+--
+rotate :: forall a. (Conjugate (Exp a), RealFloat a) => Exp (Quaternion a) -> Exp (V3 a) -> Exp (V3 a)
+rotate q v = lift ijk
+  where
+    Quaternion _ ijk = unlift $ q * (lift (Quaternion 0 (unlift v))) * conjugate q :: Quaternion (Exp a)
+
+-- | @'axisAngle' axis theta@ builds a 'Quaternion' representing a rotation of
+-- @theta@ radians about @axis@.
+--
+axisAngle :: (Epsilon a, Floating a) => Exp (V3 a) -> Exp a -> Exp (Quaternion a)
+axisAngle axis theta = lift $ Quaternion (cos half) (unlift (sin half *^ normalize axis))
+  where
+    half = theta / 2
+
+
 -- Instances
 -- ---------
 
@@ -77,4 +197,232 @@
 
 instance (Elt a, Elt b) => Each (Exp (Quaternion a)) (Exp (Quaternion b)) (Exp a) (Exp b) where
   each = liftLens (each :: Traversal (Quaternion (Exp a)) (Quaternion (Exp b)) (Exp a) (Exp b))
+
+instance Eq a => Eq (Quaternion a) where
+  (==) = (==) `on` t4
+  (/=) = (/=) `on` t4
+
+instance Ord a => Ord (Quaternion a) where
+  (<)  = (<) `on` t4
+  (>)  = (>) `on` t4
+  (<=) = (<=) `on` t4
+  (>=) = (>=) `on` t4
+  min  = qu $$ on min t4
+  max  = qu $$ on max t4
+
+t4 :: Elt a => Exp (Quaternion a) -> Exp (a,a,a,a)
+t4 (unlift -> Quaternion x (V3 y z w)) = tup4 (x,y,z,w)
+
+qu :: Elt a => Exp (a,a,a,a) -> Exp (Quaternion a)
+qu (untup4 -> (x,y,z,w)) = lift (Quaternion x (V3 y z w))
+
+instance RealFloat a => P.Num (Exp (Quaternion a)) where
+  (+)           = lift2 ((+) :: Quaternion (Exp a) -> Quaternion (Exp a) -> Quaternion (Exp a))
+  (-)           = lift2 ((-) :: Quaternion (Exp a) -> Quaternion (Exp a) -> Quaternion (Exp a))
+  negate        = fmap negate
+  abs z         = lift (Quaternion (norm z) (V3 0 0 0))
+  fromInteger x = lift (Quaternion (fromInteger x) (V3 0 0 0))
+
+  z1 * z2       = let Quaternion s1 v1' = unlift z1; v1 = lift v1'
+                      Quaternion s2 v2' = unlift z2; v2 = lift v2'
+                  in
+                  lift $ Quaternion (s1*s2 - (v1 `dot` v2))
+                                    (unlift ((v1 `cross` v2) + s1*^v2 + s2*^v1))
+
+  signum q@(unlift -> Quaternion e (V3 i j k)) =
+    if m == 0.0                      then q else
+    if not (isInfinite m || isNaN m) then q ^/ sqrt m else
+    if ne || ni || nj || nk          then qNaN else
+    if not (ii || ij || ik)          then lift $ Quaternion 1 (V3 0 0 0) else
+    if not (ie || ij || ik)          then lift $ Quaternion 0 (V3 1 0 0) else
+    if not (ie || ii || ik)          then lift $ Quaternion 0 (V3 0 1 0) else
+    if not (ie || ii || ij)          then lift $ Quaternion 0 (V3 0 0 1)
+                                     else qNaN
+    where
+      m  = quadrance q
+      ie = isInfinite e
+      ii = isInfinite i
+      ij = isInfinite j
+      ik = isInfinite k
+      ne = isNaN e
+      ni = isNaN i
+      nj = isNaN j
+      nk = isNaN k
+      --
+      qNaN = lift $ Quaternion fNaN (V3 fNaN fNaN fNaN)
+      fNaN = 0/0
+
+instance RealFloat a => P.Fractional (Exp (Quaternion a)) where
+  z1 / z2 =
+    let Quaternion q0 (V3 q1 q2 q3) = unlift z1
+        Quaternion r0 (V3 r1 r2 r3) = unlift z2
+    in
+    lift (Quaternion (r0*q0+r1*q1+r2*q2+r3*q3)
+                     (V3 (r0*q1-r1*q0-r2*q3+r3*q2)
+                         (r0*q2+r1*q3-r2*q0-r3*q1)
+                         (r0*q3-r1*q2+r2*q1-r3*q0)))
+      ^/ (r0*r0 + r1*r1 + r2*r2 + r3*r3)
+
+  recip q = let Quaternion e v = unlift q :: Quaternion (Exp a)
+            in  lift (Quaternion e (P.fmap negate v)) ^/ quadrance q
+
+  fromRational x = lift (Quaternion (fromRational x) (V3 0 0 0))
+
+instance (RealFloat a, Elt (Complex a)) => P.Floating (Exp (Quaternion a)) where
+  pi = lift (Quaternion pi (V3 0 0 0))
+
+  exp q@(unlift -> Quaternion e v) =
+    if qiq == 0
+      then lift (Quaternion exe v)
+      else reimagine (exe * cos ai) (exe * (sin ai / ai)) q
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+      exe = exp e
+
+  log q@(unlift -> Quaternion e v@(V3 _i j k)) =
+    if qiq == 0
+      then if e >= 0
+             then lift $ Quaternion (log e) v
+             else lift $ Quaternion (log (negate e)) (V3 pi j k) -- mmm, pi
+      else reimagine (log m) (atan2 m e / ai) q
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+      m   = sqrte2pqiq e qiq
+
+  x ** y = exp (y * log x)
+
+  sqrt q@(unlift -> Quaternion e v) =
+    if m   == 0 then q else
+    if qiq == 0 then if e > 0
+                        then lift $ Quaternion (sqrt e) (V3 0 0 0)
+                        else lift $ Quaternion 0 (V3 (sqrt (negate e)) 0 0)
+                else lift $ Quaternion (0.5*(m+e)) (unlift (lift v ^* im))
+    where
+      qiq = qi q
+      im  = sqrt (0.5*(m-e)) / sqrt qiq
+      m   = sqrte2pqiq e qiq
+
+  cos q@(unlift -> Quaternion e v) =
+    if qiq == 0 then lift $ Quaternion (cos e) v
+                else reimagine (cos e * cosh ai) (- sin e / ai / sinh ai) q -- 0.15 bits error
+                  -- reimagine (cos e * cosh ai) (- sin e * sinh ai / ai) q -- 13.5 bits worse
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+
+  sin q@(unlift -> Quaternion e v) =
+    if qiq == 0 then lift $ Quaternion (sin e) v
+                else reimagine (sin e * cosh ai) (cos e * sinh ai / ai) q
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+
+  tan q@(unlift -> Quaternion e v) =
+    if qiq == 0 then lift $ Quaternion (tan e) v
+                else reimagine (ce * sin e / d) (tanrhs sai ai d) q
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+      ce  = cos e
+      sai = sinh ai
+      d   = ce*ce + sai*sai
+
+  sinh q@(unlift -> Quaternion e v) =
+    if qiq == 0 then lift $ Quaternion (sinh e) v
+                else reimagine (sinh e * cos ai) (cosh e * sin ai / ai) q
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+
+  cosh q@(unlift -> Quaternion e v) =
+    if qiq == 0 then lift $ Quaternion (cosh e) v
+                else reimagine (cosh e * cos ai) (sin ai * (sinh e / ai)) q
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+
+  tanh q@(unlift -> Quaternion e v) =
+    if qiq == 0 then lift $ Quaternion (tanh e) v
+                else reimagine (cosh e * se / d) (tanhrhs cai ai d) q
+    where
+      qiq = qi q
+      ai  = sqrt qiq
+      se  = sinh e
+      cai = cos ai
+      d   = se*se + cai*cai
+
+  asin = cut asin
+  acos = cut acos
+  atan = cut atan
+
+  asinh = cut asinh
+  acosh = cut acosh
+  atanh = cut atanh
+
+
+reimagine :: RealFloat a => Exp a -> Exp a -> Exp (Quaternion a) -> Exp (Quaternion a)
+reimagine r s (unlift -> Quaternion _ v) =
+  if isNaN s || isInfinite s
+    then let aux x = if x == 0 then 0
+                               else s * x
+         in lift $ Quaternion r (P.fmap aux v)
+    else    lift $ Quaternion r (unlift (lift v ^* s))
+
+-- | Helper for calculating with specific branch cuts
+--
+cut :: (RealFloat a, Elt (Complex a)) => (Exp (Complex a) -> Exp (Complex a)) -> Exp (Quaternion a) -> Exp (Quaternion a)
+cut f q@(unlift -> Quaternion e (V3 _ y z)) =
+  if qiq == 0 then lift $ Quaternion a (V3 b y z)
+              else reimagine a (b / ai) q
+  where
+    qiq    = qi q
+    ai     = sqrt qiq
+    a :+ b = unlift $ f (lift (e :+ ai))
+
+-- | Helper for calculating with specific branch cuts
+--
+cutWith :: (RealFloat a, Elt (Complex a)) => Exp (Complex a) -> Exp (Quaternion a) -> Exp (Quaternion a)
+cutWith (unlift -> r :+ im) q@(unlift -> Quaternion e v) =
+  if e /= 0 || qiq == 0 || isNaN qiq || isInfinite qiq
+    then 0/0 -- error "bad cut"  -- TLM: argh
+    else lift $ Quaternion r (unlift (lift v ^* s))
+  where
+    qiq = qi q
+    s   = im / sqrt qiq
+
+-- | quadrance of the imaginary component
+--
+qi :: Num a => Exp (Quaternion a) -> Exp a
+qi (unlift -> Quaternion _ v :: Quaternion (Exp a)) = quadrance (lift v)
+
+sqrte2pqiq :: (Floating a, Ord a) => Exp a -> Exp a -> Exp a
+sqrte2pqiq e qiq = -- = sqrt (e*e) + qiq
+  if e < - 1.5097698010472593e153 then -(qiq/e) - e else
+  if e < 5.582399551122541e57     then sqrt ((e*e) + qiq) -- direct definition
+                                  else (qiq/e) + e
+
+tanrhs :: (Floating a, Ord a) => Exp a -> Exp a -> Exp a -> Exp a
+tanrhs sai ai d = -- = cosh ai * (sai / ai) / d -- improved from 6.04 bits of error to 0.19 bits
+  if sai < -4.618902267687042e-52 then (sai / d / ai) * cosh ai else
+  if sai < 1.038530535935153e-39  then (cosh ai * sai) / ai / d
+                                  else (sai / d / ai) * cosh ai
+
+tanhrhs :: (Floating a, Ord a) => Exp a -> Exp a -> Exp a -> Exp a
+tanhrhs cai ai d = -- = cai * (sin ai / ai) / d
+  if d >= -4.2173720203427147e-29 && d < 4.446702369113811e64
+    then cai / (d * (ai / sin ai))
+    else cai * (1 / ai / sin ai) / d
+
+instance (RealFloat a, Epsilon a) => Epsilon (Quaternion a) where
+  nearZero = nearZero . quadrance
+
+instance (RealFloat a, Conjugate (Exp a)) => Conjugate (Exp (Quaternion a)) where
+  conjugate (unlift -> Quaternion e v :: Quaternion (Exp a)) =
+    lift (Quaternion (conjugate e) (unlift (negate (lift v :: Exp (V3 a)))))
+
+instance Functor Quaternion where
+  fmap f (unlift -> Quaternion e v) = lift (Quaternion (f e) (P.fmap f v))
+  x <$ _                            = lift (Quaternion x (V3 x x x))
 
diff --git a/src/Data/Array/Accelerate/Linear/Trace.hs b/src/Data/Array/Accelerate/Linear/Trace.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Linear/Trace.hs
@@ -0,0 +1,57 @@
+{-# LANGUAGE ConstraintKinds  #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies     #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      : Data.Array.Accelerate.Linear.Trace
+-- Copyright   : 2014 Edward Kmett
+--               [2015..2018] Trevor L. McDonell
+-- License     : BSD-style (see the file LICENSE)
+--
+-- Maintainer  : Edward Kmett <ekmett@gmail.com>
+-- Stability   : experimental
+-- Portability : non-portable
+--
+-- Simple matrix operations for low-dimensional primitives
+----------------------------------------------------------------------------
+
+module Data.Array.Accelerate.Linear.Trace (
+
+  Trace(..),
+
+) where
+
+import Data.Array.Accelerate                                        as A
+import Data.Array.Accelerate.Data.Complex
+
+import Data.Array.Accelerate.Linear.Lift
+import Data.Array.Accelerate.Linear.Plucker
+import Data.Array.Accelerate.Linear.Quaternion
+import Data.Array.Accelerate.Linear.Type
+import Data.Array.Accelerate.Linear.V0
+import Data.Array.Accelerate.Linear.V1
+import Data.Array.Accelerate.Linear.V2
+import Data.Array.Accelerate.Linear.V3
+import Data.Array.Accelerate.Linear.V4
+
+import qualified Linear.Trace                                       as L
+
+
+class L.Trace m => Trace m where
+  -- | Compute the trace of a matrix
+  trace :: (A.Num a, Box2 m m a) => Exp (m (m a)) -> Exp a
+  trace = lift . L.trace . unlift'
+
+  -- | Compute the diagonal of a matrix
+  diagonal :: Box2 m m a => Exp (m (m a)) -> Exp (m a)
+  diagonal = lift . L.diagonal . unlift'
+
+instance Trace Complex
+instance Trace V0
+instance Trace V1
+instance Trace V2
+instance Trace V3
+instance Trace V4
+instance Trace Plucker
+instance Trace Quaternion
+
diff --git a/src/Data/Array/Accelerate/Linear/Type.hs b/src/Data/Array/Accelerate/Linear/Type.hs
--- a/src/Data/Array/Accelerate/Linear/Type.hs
+++ b/src/Data/Array/Accelerate/Linear/Type.hs
@@ -4,7 +4,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Type
--- Copyright   : 2015 Trevor L. McDonell
+-- Copyright   : [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
diff --git a/src/Data/Array/Accelerate/Linear/V0.hs b/src/Data/Array/Accelerate/Linear/V0.hs
--- a/src/Data/Array/Accelerate/Linear/V0.hs
+++ b/src/Data/Array/Accelerate/Linear/V0.hs
@@ -1,6 +1,5 @@
 {-# LANGUAGE ConstraintKinds       #-}
 {-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE IncoherentInstances   #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE TypeFamilies          #-}
@@ -10,7 +9,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.V0
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -27,6 +26,7 @@
 ) where
 
 import Data.Array.Accelerate                    as A
+import Data.Array.Accelerate.Data.Functor       as A
 import Data.Array.Accelerate.Smart
 import Data.Array.Accelerate.Product
 import Data.Array.Accelerate.Array.Sugar
@@ -66,14 +66,34 @@
 instance Unlift Exp (V0 a) where
   unlift _ = V0
 
+instance (Elt a, Elt b) => Each (Exp (V0 a)) (Exp (V0 b)) (Exp a) (Exp b) where
+  each _ _ = pure (constant V0)
+
 instance A.Num a => P.Num (Exp (V0 a)) where
-  _ + _ = constant V0
-  _ - _ = constant V0
-  _ * _ = constant V0
-  abs _ = constant V0
-  signum _ = constant V0
+  _ + _         = constant V0
+  _ - _         = constant V0
+  _ * _         = constant V0
+  abs _         = constant V0
+  signum _      = constant V0
   fromInteger _ = constant V0
 
-instance (Elt a, Elt b) => Each (Exp (V0 a)) (Exp (V0 b)) (Exp a) (Exp b) where
-  each _ _ = pure (constant V0)
+instance Elt a => A.Eq (V0 a) where
+  _ == _ = constant True
+  _ /= _ = constant False
+
+instance Elt a => A.Ord (V0 a) where
+  _ < _   = constant False
+  _ > _   = constant False
+  _ <= _  = constant True
+  _ >= _  = constant True
+  min _ _ = constant V0
+  max _ _ = constant V0
+
+instance Elt a => P.Bounded (Exp (V0 a)) where
+  minBound = constant V0
+  maxBound = constant V0
+
+instance A.Functor V0 where
+  fmap _ _ = constant V0
+  _ <$ _   = constant V0
 
diff --git a/src/Data/Array/Accelerate/Linear/V1.hs b/src/Data/Array/Accelerate/Linear/V1.hs
--- a/src/Data/Array/Accelerate/Linear/V1.hs
+++ b/src/Data/Array/Accelerate/Linear/V1.hs
@@ -1,7 +1,6 @@
 {-# LANGUAGE ConstraintKinds       #-}
 {-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE IncoherentInstances   #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE TypeFamilies          #-}
@@ -12,7 +11,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.V1
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -29,11 +28,13 @@
 ) where
 
 import Data.Array.Accelerate                    as A
+import Data.Array.Accelerate.Data.Functor       as A
 import Data.Array.Accelerate.Array.Sugar
 import Data.Array.Accelerate.Product
 import Data.Array.Accelerate.Smart
 import Data.Array.Accelerate.Type
 
+import Data.Array.Accelerate.Linear.Epsilon
 import Data.Array.Accelerate.Linear.Lift
 import Data.Array.Accelerate.Linear.Metric
 import Data.Array.Accelerate.Linear.Type
@@ -41,8 +42,8 @@
 
 import Control.Lens
 import Linear.V1                                ( V1(..) )
-import qualified Linear.V1                      as L
 import Prelude                                  as P
+import qualified Linear.V1                      as L
 
 -- $setup
 -- >>> import Data.Array.Accelerate.Interpreter
@@ -80,7 +81,7 @@
 type instance EltRepr (V1 a) = ((), EltRepr a)
 
 instance Elt a => Elt (V1 a) where
-  eltType _ = PairTuple UnitTuple (eltType (undefined :: a))
+  eltType _ = TypeRpair TypeRunit (eltType (undefined :: a))
   toElt ((), x) = V1 (toElt x)
   fromElt (V1 x) = ((), fromElt x)
 
@@ -97,6 +98,25 @@
 instance Elt a => Unlift Exp (V1 (Exp a)) where
   unlift t = V1 $ Exp $ ZeroTupIdx `Prj` t
 
+instance (Elt a, Elt b) => Each (Exp (V1 a)) (Exp (V1 b)) (Exp a) (Exp b) where
+  each = liftLens (each :: Traversal (V1 (Exp a)) (V1 (Exp b)) (Exp a) (Exp b))
+
+instance A.Eq a => A.Eq (V1 a) where
+  (unlift -> V1 x) == (unlift -> V1 y) = x A.== y
+  (unlift -> V1 x) /= (unlift -> V1 y) = x A./= y
+
+instance A.Ord a => A.Ord (V1 a) where
+  (unlift -> V1 x) <  (unlift -> V1 y) = x A.< y
+  (unlift -> V1 x) >  (unlift -> V1 y) = x A.> y
+  (unlift -> V1 x) >= (unlift -> V1 y) = x A.>= y
+  (unlift -> V1 x) <= (unlift -> V1 y) = x A.<= y
+  min (unlift -> V1 x) (unlift -> V1 y) = lift $ V1 (A.min x y)
+  max (unlift -> V1 x) (unlift -> V1 y) = lift $ V1 (A.max x y)
+
+instance A.Bounded a => P.Bounded (Exp (V1 a)) where
+  minBound = lift (V1 (minBound :: Exp a))
+  maxBound = lift (V1 (maxBound :: Exp a))
+
 instance A.Num a => P.Num (Exp (V1 a)) where
   (+)             = lift2 ((+) :: V1 (Exp a) -> V1 (Exp a) -> V1 (Exp a))
   (-)             = lift2 ((-) :: V1 (Exp a) -> V1 (Exp a) -> V1 (Exp a))
@@ -104,12 +124,12 @@
   negate          = lift1 (negate :: V1 (Exp a) -> V1 (Exp a))
   signum          = lift1 (signum :: V1 (Exp a) -> V1 (Exp a))
   abs             = lift1 (signum :: V1 (Exp a) -> V1 (Exp a))
-  fromInteger x   = lift (fromInteger x :: V1 (Exp a))
+  fromInteger x   = lift (P.fromInteger x :: V1 (Exp a))
 
 instance A.Floating a => P.Fractional (Exp (V1 a)) where
   (/)             = lift2 ((/) :: V1 (Exp a) -> V1 (Exp a) -> V1 (Exp a))
   recip           = lift1 (recip :: V1 (Exp a) -> V1 (Exp a))
-  fromRational x  = lift (fromRational x :: V1 (Exp a))
+  fromRational x  = lift (P.fromRational x :: V1 (Exp a))
 
 instance A.Floating a => P.Floating (Exp (V1 a)) where
   pi              = lift (pi :: V1 (Exp a))
@@ -128,6 +148,10 @@
   acosh           = lift1 (acosh :: V1 (Exp a) -> V1 (Exp a))
   atanh           = lift1 (atanh :: V1 (Exp a) -> V1 (Exp a))
 
-instance (Elt a, Elt b) => Each (Exp (V1 a)) (Exp (V1 b)) (Exp a) (Exp b) where
-  each = liftLens (each :: Traversal (V1 (Exp a)) (V1 (Exp b)) (Exp a) (Exp b))
+instance Epsilon a => Epsilon (V1 a) where
+  nearZero (unlift -> V1 x) = nearZero x
+
+instance A.Functor V1 where
+  fmap f (unlift -> V1 x) = lift (V1 (f x))
+  x <$ _                  = lift (V1 x)
 
diff --git a/src/Data/Array/Accelerate/Linear/V2.hs b/src/Data/Array/Accelerate/Linear/V2.hs
--- a/src/Data/Array/Accelerate/Linear/V2.hs
+++ b/src/Data/Array/Accelerate/Linear/V2.hs
@@ -1,7 +1,6 @@
 {-# LANGUAGE ConstraintKinds       #-}
 {-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE IncoherentInstances   #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE TypeFamilies          #-}
@@ -12,7 +11,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.V2
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -32,10 +31,12 @@
 ) where
 
 import Data.Array.Accelerate                    as A
+import Data.Array.Accelerate.Data.Functor       as A
 import Data.Array.Accelerate.Smart
 import Data.Array.Accelerate.Product
 import Data.Array.Accelerate.Array.Sugar
 
+import Data.Array.Accelerate.Linear.Epsilon
 import Data.Array.Accelerate.Linear.Lift
 import Data.Array.Accelerate.Linear.Metric
 import Data.Array.Accelerate.Linear.Type
@@ -43,9 +44,10 @@
 import Data.Array.Accelerate.Linear.Vector
 
 import Control.Lens
+import Data.Function
 import Linear.V2                                ( V2(..) )
-import qualified Linear.V2                      as L
 import Prelude                                  as P
+import qualified Linear.V2                      as L
 
 -- $setup
 -- >>> import Data.Array.Accelerate.Interpreter
@@ -131,6 +133,31 @@
   unlift t = V2 (Exp $ SuccTupIdx ZeroTupIdx `Prj` t)
                 (Exp $ ZeroTupIdx `Prj` t)
 
+instance (Elt a, Elt b) => Each (Exp (V2 a)) (Exp (V2 b)) (Exp a) (Exp b) where
+  each = liftLens (each :: Traversal (V2 (Exp a)) (V2 (Exp b)) (Exp a) (Exp b))
+
+instance A.Eq a => A.Eq (V2 a) where
+  (==) = (A.==) `on` t2
+  (/=) = (A./=) `on` t2
+
+instance A.Ord a => A.Ord (V2 a) where
+  (<)  = (A.<) `on` t2
+  (>)  = (A.>) `on` t2
+  (<=) = (A.<=) `on` t2
+  (>=) = (A.>=) `on` t2
+  min  = v2 $$ on A.min t2
+  max  = v2 $$ on A.max t2
+
+t2 :: Elt a => Exp (V2 a) -> Exp (a,a)
+t2 (unlift -> V2 x y) = tup2 (x,y)
+
+v2 :: Elt a => Exp (a,a) -> Exp (V2 a)
+v2 (untup2 -> (x,y)) = lift (V2 x y)
+
+instance A.Bounded a => P.Bounded (Exp (V2 a)) where
+  minBound = lift (V2 (minBound :: Exp a) (minBound :: Exp a))
+  maxBound = lift (V2 (maxBound :: Exp a) (maxBound :: Exp a))
+
 instance A.Num a => P.Num (Exp (V2 a)) where
   (+)             = lift2 ((+) :: V2 (Exp a) -> V2 (Exp a) -> V2 (Exp a))
   (-)             = lift2 ((-) :: V2 (Exp a) -> V2 (Exp a) -> V2 (Exp a))
@@ -138,12 +165,12 @@
   negate          = lift1 (negate :: V2 (Exp a) -> V2 (Exp a))
   signum          = lift1 (signum :: V2 (Exp a) -> V2 (Exp a))
   abs             = lift1 (signum :: V2 (Exp a) -> V2 (Exp a))
-  fromInteger x   = lift (fromInteger x :: V2 (Exp a))
+  fromInteger x   = lift (P.fromInteger x :: V2 (Exp a))
 
 instance A.Floating a => P.Fractional (Exp (V2 a)) where
   (/)             = lift2 ((/) :: V2 (Exp a) -> V2 (Exp a) -> V2 (Exp a))
   recip           = lift1 (recip :: V2 (Exp a) -> V2 (Exp a))
-  fromRational x  = lift (fromRational x :: V2 (Exp a))
+  fromRational x  = lift (P.fromRational x :: V2 (Exp a))
 
 instance A.Floating a => P.Floating (Exp (V2 a)) where
   pi              = lift (pi :: V2 (Exp a))
@@ -162,6 +189,10 @@
   acosh           = lift1 (acosh :: V2 (Exp a) -> V2 (Exp a))
   atanh           = lift1 (atanh :: V2 (Exp a) -> V2 (Exp a))
 
-instance (Elt a, Elt b) => Each (Exp (V2 a)) (Exp (V2 b)) (Exp a) (Exp b) where
-  each = liftLens (each :: Traversal (V2 (Exp a)) (V2 (Exp b)) (Exp a) (Exp b))
+instance Epsilon a => Epsilon (V2 a) where
+  nearZero = nearZero . quadrance
+
+instance A.Functor V2 where
+  fmap f (unlift -> V2 x y) = lift (V2 (f x) (f y))
+  x <$ _                    = lift (V2 x x)
 
diff --git a/src/Data/Array/Accelerate/Linear/V3.hs b/src/Data/Array/Accelerate/Linear/V3.hs
--- a/src/Data/Array/Accelerate/Linear/V3.hs
+++ b/src/Data/Array/Accelerate/Linear/V3.hs
@@ -1,7 +1,6 @@
 {-# LANGUAGE ConstraintKinds       #-}
 {-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE IncoherentInstances   #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE TypeFamilies          #-}
@@ -12,7 +11,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.V3
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -36,10 +35,12 @@
 ) where
 
 import Data.Array.Accelerate                    as A
+import Data.Array.Accelerate.Data.Functor       as A
 import Data.Array.Accelerate.Smart
 import Data.Array.Accelerate.Product
 import Data.Array.Accelerate.Array.Sugar
 
+import Data.Array.Accelerate.Linear.Epsilon
 import Data.Array.Accelerate.Linear.Lift
 import Data.Array.Accelerate.Linear.Metric
 import Data.Array.Accelerate.Linear.Type
@@ -48,9 +49,10 @@
 import Data.Array.Accelerate.Linear.Vector
 
 import Control.Lens
+import Data.Function
 import Linear.V3                                ( V3(..) )
-import qualified Linear.V3                      as L
 import Prelude                                  as P
+import qualified Linear.V3                      as L
 
 -- $setup
 -- >>> import Data.Array.Accelerate.Interpreter
@@ -140,6 +142,31 @@
                 (Exp $ SuccTupIdx ZeroTupIdx `Prj` t)
                 (Exp $ ZeroTupIdx `Prj` t)
 
+instance (Elt a, Elt b) => Each (Exp (V3 a)) (Exp (V3 b)) (Exp a) (Exp b) where
+  each = liftLens (each :: Traversal (V3 (Exp a)) (V3 (Exp b)) (Exp a) (Exp b))
+
+instance A.Eq a => A.Eq (V3 a) where
+  (==) = (A.==) `on` t3
+  (/=) = (A./=) `on` t3
+
+instance A.Ord a => A.Ord (V3 a) where
+  (<)  = (A.<) `on` t3
+  (>)  = (A.>) `on` t3
+  (<=) = (A.<=) `on` t3
+  (>=) = (A.>=) `on` t3
+  min  = v3 $$ on A.min t3
+  max  = v3 $$ on A.max t3
+
+t3 :: Elt a => Exp (V3 a) -> Exp (a,a,a)
+t3 (unlift -> V3 x y z) = tup3 (x,y,z)
+
+v3 :: Elt a => Exp (a,a,a) -> Exp (V3 a)
+v3 (untup3 -> (x,y,z)) = lift (V3 x y z)
+
+instance A.Bounded a => P.Bounded (Exp (V3 a)) where
+  minBound = lift (V3 (minBound :: Exp a) (minBound :: Exp a) (minBound :: Exp a))
+  maxBound = lift (V3 (maxBound :: Exp a) (maxBound :: Exp a) (maxBound :: Exp a))
+
 instance A.Num a => P.Num (Exp (V3 a)) where
   (+)             = lift2 ((+) :: V3 (Exp a) -> V3 (Exp a) -> V3 (Exp a))
   (-)             = lift2 ((-) :: V3 (Exp a) -> V3 (Exp a) -> V3 (Exp a))
@@ -147,12 +174,12 @@
   negate          = lift1 (negate :: V3 (Exp a) -> V3 (Exp a))
   signum          = lift1 (signum :: V3 (Exp a) -> V3 (Exp a))
   abs             = lift1 (signum :: V3 (Exp a) -> V3 (Exp a))
-  fromInteger x   = lift (fromInteger x :: V3 (Exp a))
+  fromInteger x   = lift (P.fromInteger x :: V3 (Exp a))
 
 instance A.Floating a => P.Fractional (Exp (V3 a)) where
   (/)             = lift2 ((/) :: V3 (Exp a) -> V3 (Exp a) -> V3 (Exp a))
   recip           = lift1 (recip :: V3 (Exp a) -> V3 (Exp a))
-  fromRational x  = lift (fromRational x :: V3 (Exp a))
+  fromRational x  = lift (P.fromRational x :: V3 (Exp a))
 
 instance A.Floating a => P.Floating (Exp (V3 a)) where
   pi              = lift (pi :: V3 (Exp a))
@@ -171,8 +198,12 @@
   acosh           = lift1 (acosh :: V3 (Exp a) -> V3 (Exp a))
   atanh           = lift1 (atanh :: V3 (Exp a) -> V3 (Exp a))
 
-instance (Elt a, Elt b) => Each (Exp (V3 a)) (Exp (V3 b)) (Exp a) (Exp b) where
-  each = liftLens (each :: Traversal (V3 (Exp a)) (V3 (Exp b)) (Exp a) (Exp b))
+instance Epsilon a => Epsilon (V3 a) where
+  nearZero = nearZero . quadrance
+
+instance A.Functor V3 where
+  fmap f (unlift -> V3 x y z) = lift (V3 (f x) (f y) (f z))
+  x <$ _                      = lift (V3 x x x)
 
 -- $liftAcc
 --
diff --git a/src/Data/Array/Accelerate/Linear/V4.hs b/src/Data/Array/Accelerate/Linear/V4.hs
--- a/src/Data/Array/Accelerate/Linear/V4.hs
+++ b/src/Data/Array/Accelerate/Linear/V4.hs
@@ -1,7 +1,6 @@
 {-# LANGUAGE ConstraintKinds       #-}
 {-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE IncoherentInstances   #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE TypeFamilies          #-}
@@ -12,7 +11,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.V4
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -43,10 +42,12 @@
 ) where
 
 import Data.Array.Accelerate                    as A
+import Data.Array.Accelerate.Data.Functor       as A
 import Data.Array.Accelerate.Smart
 import Data.Array.Accelerate.Product
 import Data.Array.Accelerate.Array.Sugar
 
+import Data.Array.Accelerate.Linear.Epsilon
 import Data.Array.Accelerate.Linear.Lift
 import Data.Array.Accelerate.Linear.Metric
 import Data.Array.Accelerate.Linear.Type
@@ -56,9 +57,10 @@
 import Data.Array.Accelerate.Linear.Vector
 
 import Control.Lens
+import Data.Function
 import Linear.V4                                ( V4(..) )
-import qualified Linear.V4                      as L
 import Prelude                                  as P
+import qualified Linear.V4                      as L
 
 -- $setup
 -- >>> import Data.Array.Accelerate.Interpreter
@@ -85,7 +87,7 @@
 -- the projective, homogenous, coordinate @[x:y:z:w]@ is one of many associated
 -- with a single point @(x\/w, y\/w, z\/w)@.
 --
-normalizePoint :: forall a. A.Floating a => Exp (V4 a) -> Exp (V3 a)
+normalizePoint :: forall a. A.Fractional a => Exp (V4 a) -> Exp (V3 a)
 normalizePoint = lift1 (L.normalizePoint :: V4 (Exp a) -> V3 (Exp a))
 
 -- | A space that distinguishes orthogonal basis vectors '_x', '_y', '_z', and '_w'.
@@ -210,6 +212,31 @@
                 (Exp $ SuccTupIdx ZeroTupIdx `Prj` t)
                 (Exp $ ZeroTupIdx `Prj` t)
 
+instance (Elt a, Elt b) => Each (Exp (V4 a)) (Exp (V4 b)) (Exp a) (Exp b) where
+  each = liftLens (each :: Traversal (V4 (Exp a)) (V4 (Exp b)) (Exp a) (Exp b))
+
+instance A.Eq a => A.Eq (V4 a) where
+  (==) = (A.==) `on` t4
+  (/=) = (A./=) `on` t4
+
+instance A.Ord a => A.Ord (V4 a) where
+  (<)  = (A.<) `on` t4
+  (>)  = (A.>) `on` t4
+  (<=) = (A.<=) `on` t4
+  (>=) = (A.>=) `on` t4
+  min  = v4 $$ on A.min t4
+  max  = v4 $$ on A.max t4
+
+t4 :: Elt a => Exp (V4 a) -> Exp (a,a,a,a)
+t4 (unlift -> V4 x y z w) = tup4 (x,y,z,w)
+
+v4 :: Elt a => Exp (a,a,a,a) -> Exp (V4 a)
+v4 (untup4 -> (x,y,z,w)) = lift (V4 x y z w)
+
+instance A.Bounded a => P.Bounded (Exp (V4 a)) where
+  minBound = lift (V4 (minBound :: Exp a) (minBound :: Exp a) (minBound :: Exp a) (minBound :: Exp a))
+  maxBound = lift (V4 (maxBound :: Exp a) (maxBound :: Exp a) (maxBound :: Exp a) (maxBound :: Exp a))
+
 instance A.Num a => P.Num (Exp (V4 a)) where
   (+)             = lift2 ((+) :: V4 (Exp a) -> V4 (Exp a) -> V4 (Exp a))
   (-)             = lift2 ((-) :: V4 (Exp a) -> V4 (Exp a) -> V4 (Exp a))
@@ -217,12 +244,12 @@
   negate          = lift1 (negate :: V4 (Exp a) -> V4 (Exp a))
   signum          = lift1 (signum :: V4 (Exp a) -> V4 (Exp a))
   abs             = lift1 (signum :: V4 (Exp a) -> V4 (Exp a))
-  fromInteger x   = lift (fromInteger x :: V4 (Exp a))
+  fromInteger x   = lift (P.fromInteger x :: V4 (Exp a))
 
 instance A.Floating a => P.Fractional (Exp (V4 a)) where
   (/)             = lift2 ((/) :: V4 (Exp a) -> V4 (Exp a) -> V4 (Exp a))
   recip           = lift1 (recip :: V4 (Exp a) -> V4 (Exp a))
-  fromRational x  = lift (fromRational x :: V4 (Exp a))
+  fromRational x  = lift (P.fromRational x :: V4 (Exp a))
 
 instance A.Floating a => P.Floating (Exp (V4 a)) where
   pi              = lift (pi :: V4 (Exp a))
@@ -241,6 +268,10 @@
   acosh           = lift1 (acosh :: V4 (Exp a) -> V4 (Exp a))
   atanh           = lift1 (atanh :: V4 (Exp a) -> V4 (Exp a))
 
-instance (Elt a, Elt b) => Each (Exp (V4 a)) (Exp (V4 b)) (Exp a) (Exp b) where
-  each = liftLens (each :: Traversal (V4 (Exp a)) (V4 (Exp b)) (Exp a) (Exp b))
+instance Epsilon a => Epsilon (V4 a) where
+  nearZero = nearZero . quadrance
+
+instance A.Functor V4 where
+  fmap f (unlift -> V4 x y z w) = lift (V4 (f x) (f y) (f z) (f w))
+  x <$ _                        = lift (V4 x x x x)
 
diff --git a/src/Data/Array/Accelerate/Linear/Vector.hs b/src/Data/Array/Accelerate/Linear/Vector.hs
--- a/src/Data/Array/Accelerate/Linear/Vector.hs
+++ b/src/Data/Array/Accelerate/Linear/Vector.hs
@@ -8,7 +8,7 @@
 -- |
 -- Module      : Data.Array.Accelerate.Linear.Vector
 -- Copyright   : 2014 Edward Kmett, Charles Durham,
---               2015 Trevor L. McDonell
+--               [2015..2018] Trevor L. McDonell
 -- License     : BSD-style (see the file LICENSE)
 --
 -- Maintainer  : Edward Kmett <ekmett@gmail.com>
@@ -83,7 +83,7 @@
   lerp = lift3 (L.lerp :: Exp a -> f (Exp a) -> f (Exp a) -> f (Exp a))
 
 
-type IsAdditive f a = (Additive f, Box f a)
+-- type IsAdditive f a = (Additive f, Box f a)
 
 
 -- | Basis element