diff --git a/.ghci b/.ghci
new file mode 100644
--- /dev/null
+++ b/.ghci
@@ -0,0 +1,1 @@
+:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
diff --git a/.gitignore b/.gitignore
new file mode 100644
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+dist
+docs
+wiki
+TAGS
+tags
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,8 @@
+language: haskell
+# Uncomment the next 4 lines whenever hackage is down.
+# before_install:
+#   - mkdir -p ~/.cabal
+#   - cp config ~/.cabal/config
+#   - cabal update
+notifications:
+  irc: "irc.freenode.org#haskell-lens"
diff --git a/.vim.custom b/.vim.custom
new file mode 100644
--- /dev/null
+++ b/.vim.custom
@@ -0,0 +1,21 @@
+" Add the following to your .vimrc to automatically load this on startup
+" if filereadable(".vim.custom")
+"     so .vim.custom
+" endif
+
+function StripTrailingWhitespace()
+  let myline=line(".")
+  let mycolumn = col(".")
+  silent %s/  *$//
+  call cursor(myline, mycolumn)
+endfunction
+
+syntax on
+set tags=TAGS;/
+set listchars=tab:‗‗,trail:‗
+set list
+
+map <F2> :exec ":!hasktags -x -c --ignore src"<CR><CR>
+
+au BufWritePre *.hs,*.markdown silent! cal StripTrailingWhitespace()
+au BufWritePost *.hs silent! :exec ":!hasktags -x -c --ignore src"
diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.markdown
@@ -0,0 +1,3 @@
+0.2
+---
+* Initial hackage release
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright 2011-12 Edward Kmett
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.markdown b/README.markdown
new file mode 100644
--- /dev/null
+++ b/README.markdown
@@ -0,0 +1,13 @@
+linear
+======
+
+[![Build Status](https://secure.travis-ci.org/ekmett/linear.png?branch=master)](http://travis-ci.org/ekmett/linear)
+
+Contact Information
+-------------------
+
+Contributions and bug reports are welcome!
+
+Please feel free to contact me through github or on the #haskell IRC channel on irc.freenode.net.
+
+-Edward Kmett
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,7 @@
+#!/usr/bin/runhaskell
+> module Main (main) where
+
+> import Distribution.Simple
+
+> main :: IO ()
+> main = defaultMain
diff --git a/config b/config
new file mode 100644
--- /dev/null
+++ b/config
@@ -0,0 +1,16 @@
+-- This provides a custom ~/.cabal/config file for use when hackage is down that should work on unix
+--
+-- This is particularly useful for travis-ci to get it to stop complaining
+-- about a broken build when everything is still correct on our end.
+--
+-- This uses Luite Stegeman's mirror of hackage provided by his 'hdiff' site instead
+--
+-- To enable this, uncomment the before_script in .travis.yml
+
+remote-repo: hdiff.luite.com:http://hdiff.luite.com/packages/archive
+remote-repo-cache: ~/.cabal/packages
+world-file: ~/.cabal/world
+build-summary: ~/.cabal/logs/build.log
+remote-build-reporting: anonymous
+install-dirs user
+install-dirs global
diff --git a/linear.cabal b/linear.cabal
new file mode 100644
--- /dev/null
+++ b/linear.cabal
@@ -0,0 +1,62 @@
+name:          linear
+category:      Math, Algebra
+version:       0.2
+license:       BSD3
+cabal-version: >= 1.8
+license-file:  LICENSE
+author:        Edward A. Kmett
+maintainer:    Edward A. Kmett <ekmett@gmail.com>
+stability:     provisional
+homepage:      http://github.com/ekmett/linear/
+bug-reports:   http://github.com/ekmett/linear/issues
+copyright:     Copyright (C) 2012 Edward A. Kmett
+synopsis:      Linear Algebra
+description:   Types and combinators for low-dimension-count linear algebra on free vectors spaces
+build-type:    Simple
+tested-with:   GHC == 7.4.1
+extra-source-files:
+  .travis.yml
+  .ghci
+  .gitignore
+  .vim.custom
+  config
+  README.markdown
+  CHANGELOG.markdown
+
+source-repository head
+  type: git
+  location: git://github.com/ekmett/linear.git
+
+library
+  build-depends:
+    base == 4.*,
+    distributive >= 0.2.2 && < 0.3,
+    lens == 2.9.*
+
+  exposed-modules:
+    Linear.Conjugate
+    Linear.Epsilon
+    Linear.Matrix
+    Linear.Metric
+    Linear.Plucker
+    Linear.Quaternion
+    Linear.V2
+    Linear.V3
+    Linear.V4
+    Linear.Vector
+
+  ghc-options: -Wall -fwarn-tabs -O2 -fdicts-cheap -funbox-strict-fields
+  hs-source-dirs: src
+
+-- Verify the results of the examples
+test-suite doctests
+  type:    exitcode-stdio-1.0
+  main-is: doctests.hs
+  build-depends:
+    base == 4.*,
+    directory >= 1.0 && < 1.2,
+    doctest >= 0.8 && <= 0.9,
+    filepath >= 1.3 && < 1.4
+  ghc-options: -Wall -Werror -threaded
+  hs-source-dirs: tests
+
diff --git a/src/Linear/Conjugate.hs b/src/Linear/Conjugate.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/Conjugate.hs
@@ -0,0 +1,18 @@
+module Linear.Conjugate
+  ( Conjugate(..)
+  ) where
+
+import Data.Complex hiding (conjugate)
+
+-- | An involutive ring
+class Num a => Conjugate a where
+  -- | Conjugate a value. This defaults to the trivial involution.
+  conjugate :: a -> a
+  conjugate = id
+
+instance Conjugate Double
+instance Conjugate Float
+instance (Conjugate a, RealFloat a) => Conjugate (Complex a) where
+  {-# SPECIALIZE instance Conjugate (Complex Float) #-}
+  {-# SPECIALIZE instance Conjugate (Complex Double) #-}
+  conjugate (a :+ b) = conjugate a :+ negate b
diff --git a/src/Linear/Epsilon.hs b/src/Linear/Epsilon.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/Epsilon.hs
@@ -0,0 +1,26 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Linear.Epsilon
+-- Copyright   :  (C) 2012 Edward Kmett
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-----------------------------------------------------------------------------
+module Linear.Epsilon
+  ( Epsilon(..)
+  ) where
+
+-- | Provides a fairly subjective test to see if a quantity is near zero.
+class Num a => Epsilon a where
+  -- | Determine if a quantity is near zero.
+  nearZero :: a -> Bool
+
+-- | @'abs' a '<=' 1e-6@
+instance Epsilon Float where
+  nearZero a = abs a <= 1e-6
+
+-- | @'abs' a '<=' 1e-12@
+instance Epsilon Double where
+  nearZero a = abs a <= 1e-12
diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/Matrix.hs
@@ -0,0 +1,95 @@
+module Linear.Matrix
+  ( (!*!), (!*) , (*!)
+  , adjoint
+  , M33, M44, M43, m33_to_m44, m43_to_m44
+  , eye3, eye4
+  , trace
+  , translation
+  , fromQuaternion
+  , mkTransformation
+  ) where
+
+import Control.Applicative
+import Control.Lens
+import Data.Distributive
+import Data.Foldable as Foldable
+import Linear.Quaternion
+import Linear.V3
+import Linear.V4
+import Linear.Metric
+import Linear.Conjugate
+
+infixl 7 !*!
+-- | matrix product
+(!*!) :: (Functor m, Foldable r, Applicative r, Distributive n, Num a) => m (r a) -> r (n a) -> m (n a)
+f !*! g = fmap (\r -> Foldable.foldr (+) 0 . liftA2 (*) r <$> g') f
+  where g' = distribute g
+
+-- | matrix * column vector
+infixl 7 *!
+(!*) :: (Functor m, Metric r, Num a) => m (r a) -> r a -> m a
+m !* v = dot v <$> m
+
+infixl 7 !*
+-- | row vector * matrix
+(*!) :: (Metric r, Distributive n, Num a) => r a -> r (n a) -> n a
+f *! g = dot f <$> distribute g
+
+-- | hermitian conjugate or conjugate transpose
+adjoint :: (Functor m, Distributive n, Conjugate a) => m (n a) -> n (m a)
+adjoint = collect (fmap conjugate)
+{-# INLINE adjoint #-}
+
+-- | Compute the trace of a matrix
+trace :: (Monad f, Foldable f, Num a) => f (f a) -> a
+trace m = Foldable.sum (m >>= id)
+{-# INLINE trace #-}
+
+-- | Matrices use a row-major representation.
+type M33 a = V3 (V3 a)
+type M44 a = V4 (V4 a)
+type M43 a = V4 (V3 a)
+
+-- | Build a rotation matrix from a 'Quaternion'.
+fromQuaternion :: Num a => Quaternion a -> M33 a
+fromQuaternion (Quaternion a (V3 b c d)) =
+  V3 (V3 (a*a+b*b-c*c-d*d) (2*b*c-2*a*d) (2*b*d+2*a*c))
+     (V3 (2*b*c+2*a*d) (a*a-b*b+c*c-d*d) (2*c*d-2*a*b))
+     (V3 (2*b*d-2*a*c) (2*c*d+2*a*b) (a*a-b*b-c*c+d*d))
+
+mkTransformationMat :: Num a => M33 a -> V3 a -> M44 a
+mkTransformationMat (V3 r1 r2 r3) (V3 tx ty tz) =
+  V4 (snoc3 r1 tx) (snoc3 r2 ty) (snoc3 r3 tz) (set _w 1 0)
+  where snoc3 (V3 x y z) w = V4 x y z w
+
+-- |Build a transformation matrix from a rotation expressed as a
+-- 'Quaternion' and a translation vector.
+mkTransformation :: Num a => Quaternion a -> V3 a -> M44 a
+mkTransformation = mkTransformationMat . fromQuaternion
+
+m43_to_m44 :: Num a => M43 a -> M44 a
+m43_to_m44
+  (V4 (V3 a b c)
+      (V3 d e f)
+      (V3 g h i)
+      (V3 j k l)) =
+  (V4 (V4 a b c 0)
+      (V4 d e f 0)
+      (V4 g h i 0)
+      (V4 j k l 1))
+
+m33_to_m44 :: Num a => M33 a -> M44 a
+m33_to_m44 (V3 r1 r2 r3) = V4 (vector r1) (vector r2) (vector r3) (point 0)
+
+-- |3x3 identity matrix.
+eye3 :: Num a => M33 a
+eye3 = V3 (set _x 1 0) (set _y 1 0) (set _z 1 0)
+
+-- |4x4 identity matrix.
+eye4 :: Num a => M44 a
+eye4 = V4 (set _x 1 0) (set _y 1 0) (set _z 1 0) (set _w 1 0)
+
+-- |Extract the translation vector (first three entries of the last
+-- column) from a 3x4 or 4x4 matrix
+translation :: (R3 t, R4 v, Functor f, Functor t) => (V3 a -> f (V3 a)) -> t (v a) -> f (t a)
+translation = (. fmap (^._w)) . _xyz
diff --git a/src/Linear/Metric.hs b/src/Linear/Metric.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/Metric.hs
@@ -0,0 +1,40 @@
+module Linear.Metric
+  ( Metric(..), normalize
+  ) where
+
+import Control.Applicative
+import Linear.Epsilon
+
+-- | A free inner product/metric space
+class Applicative f => Metric f where
+  -- | Compute the inner product of two vectors or (equivalently)
+  -- convert a vector @f a@ into a covector @f a -> a@.
+  dot :: Num a => f a -> f a -> a
+
+  -- | Compute the squared norm. The name quadrance arises from
+  -- Norman J. Wildberger's rational trigonometry.
+  quadrance :: Num a => f a -> a
+  quadrance v = dot v v
+
+  -- | Compute the quadrance of the difference
+  qd :: Num a => f a -> f a -> a
+  qd f g = quadrance (liftA2 (-) f g)
+
+  -- | Compute the distance between two vectors in a metric space
+  distance :: Floating a => f a -> f a -> a
+  distance f g = norm (liftA2 (-) f g)
+
+  -- | Compute the norm of a vector in a metric space
+  norm :: Floating a => f a -> a
+  norm v = sqrt (dot v v)
+
+  -- | Convert a non-zero vector to unit vector.
+  signorm :: Floating a => f a -> f a
+  signorm v = fmap (/m) v where
+    m = norm v
+
+-- | Normalize a 'Metric' functor to have unit 'norm'. This function
+-- does not change the functor if its 'norm' is 0 or 1.
+normalize :: (Floating a, Metric f, Epsilon a) => f a -> f a
+normalize v = if nearZero l || nearZero (1-l) then v else fmap (/sqrt l) v
+  where l = quadrance v
diff --git a/src/Linear/Plucker.hs b/src/Linear/Plucker.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/Plucker.hs
@@ -0,0 +1,114 @@
+module Linear.Plucker
+  ( Plucker(..)
+  , squaredError
+  , isotropic
+  , (><)
+  , plucker
+  , intersects
+  ) where
+
+import Control.Applicative
+import Data.Distributive
+import Data.Foldable as Foldable
+import Data.Monoid
+import Data.Traversable
+import Linear.Epsilon
+import Linear.Metric
+import Control.Lens
+import Linear.V4
+
+-- | Plücker coordinates for lines in a 3-dimensional space.
+data Plucker a = Plucker a a a a a a deriving (Eq,Ord,Show,Read)
+
+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)
+
+instance Applicative Plucker where
+  pure a = Plucker a a a a a a
+  Plucker a b c d e f <*> Plucker g h i j k l =
+    Plucker (a g) (b h) (c i) (d j) (e k) (f l)
+
+instance Monad Plucker where
+  return a = Plucker a a a a a a
+  (>>=) = bindRep
+
+instance Distributive Plucker where
+  distribute = distributeRep
+
+instance Representable Plucker where
+  rep f = Plucker (f p01) (f p02) (f p03) (f p23) (f p31) (f p12)
+
+instance Foldable Plucker where
+  foldMap g (Plucker a b c d e f) =
+    g a `mappend` g b `mappend` g c `mappend` g d `mappend` g e `mappend` g f
+
+instance Traversable Plucker where
+  traverse g (Plucker a b c d e f) =
+    Plucker <$> g a <*> g b <*> g c <*> g d <*> g e <*> g f
+
+instance Num a => Num (Plucker a) where
+  (+) = liftA2 (+)
+  (*) = liftA2 (*)
+  negate = fmap negate
+  abs = fmap abs
+  signum = fmap signum
+  fromInteger = pure . fromInteger
+
+instance Fractional a => Fractional (Plucker a) where
+  recip = fmap recip
+  (/) = liftA2 (/)
+  fromRational = pure . fromRational
+
+-- | Given a pair of points represented by homogeneous coordinates generate Plücker coordinates
+-- for the line through them.
+plucker :: Num a => V4 a -> V4 a -> Plucker a
+plucker (V4 a b c d)
+        (V4 e f g h) =
+  Plucker (a*f-b*e)
+          (a*g-c*e)
+          (a*d-h*e)
+          (c*h-d*g)
+          (d*f-b*h)
+          (b*g-c*f)
+
+-- | These elements form a basis for the Plücker space, or the Grassmanian manifold @Gr(2,V4)@.
+p01, p02, p03, p23, p31, p12 :: Functor f => (a -> f a) -> Plucker a -> f (Plucker a)
+p01 g (Plucker a b c d e f) = (\a' -> Plucker a' b c d e f) <$> g a
+p02 g (Plucker a b c d e f) = (\b' -> Plucker a b' c d e f) <$> g b
+p03 g (Plucker a b c d e f) = (\c' -> Plucker a b c' d e f) <$> g c
+p23 g (Plucker a b c d e f) = (\d' -> Plucker a b c d' e f) <$> g d
+p31 g (Plucker a b c d e f) = (\e' -> Plucker a b c d e' f) <$> g e
+p12 g (Plucker a b c d e f) = Plucker a b c d e <$> g f
+{-# INLINE p01 #-}
+{-# INLINE p02 #-}
+{-# INLINE p03 #-}
+{-# INLINE p23 #-}
+{-# INLINE p31 #-}
+{-# INLINE p12 #-}
+
+-- | Valid Plücker coordinates @p@ will have @'squaredError' p '==' 0@
+--
+-- That said, floating point makes a mockery of this claim, so you may want to use 'nearZero'.
+squaredError :: (Eq a, Num a) => Plucker a -> a
+squaredError v = v >< v
+
+-- | This isn't th actual metric because this bilinear form gives rise to an isotropic quadratic space
+infixl 5 ><
+(><) :: Num a => Plucker a -> Plucker a -> a
+Plucker a b c d e f >< Plucker g h i j k l = a*g+b*h+c*i-d*j-e*k-f*l
+
+-- | Checks if the line is near-isotropic (isotropic vectors in this quadratic space represent lines in real 3d space)
+isotropic :: Epsilon a => Plucker a -> Bool
+isotropic a = nearZero (a >< a)
+
+-- | Checks if the two vectors intersect (or nearly intersect)
+intersects :: Epsilon a => Plucker a -> Plucker a -> Bool
+intersects a b = nearZero (a >< b)
+
+instance Metric Plucker where
+  dot (Plucker a b c d e f) (Plucker g h i j k l) = a*g+b*h+c*i+d*j+e*k+f*l
+
+instance Epsilon a => Epsilon (Plucker a) where
+  nearZero = nearZero . quadrance
+
+-- TODO: drag some more stuff out of my thesis
diff --git a/src/Linear/Quaternion.hs b/src/Linear/Quaternion.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/Quaternion.hs
@@ -0,0 +1,337 @@
+{-# LANGUAGE DeriveDataTypeable, PatternGuards, ScopedTypeVariables #-}
+module Linear.Quaternion
+  ( Quaternion(..)
+  , Complicated(..)
+  , Hamiltonian(..)
+  , slerp
+  , asinq
+  , acosq
+  , atanq
+  , asinhq
+  , acoshq
+  , atanhq
+  , absi
+  , pow
+  , rotate
+  , axisAngle
+  ) where
+import Control.Applicative
+import Control.Lens
+import Data.Complex (Complex((:+)))
+import Data.Data
+import Data.Distributive
+import Data.Foldable
+import qualified Data.Foldable as F
+import Data.Monoid
+import Foreign.Ptr (castPtr, plusPtr)
+import Foreign.Storable (Storable(..))
+import Linear.Epsilon
+import Linear.Conjugate
+import Linear.Metric
+import Linear.V3
+import Linear.Vector
+import Prelude hiding (any)
+
+data Quaternion a = Quaternion a {-# UNPACK #-}!(V3 a)
+                    deriving (Eq,Ord,Read,Show,Data,Typeable)
+
+instance Functor Quaternion where
+  fmap f (Quaternion e v) = Quaternion (f e) (fmap f v)
+  a <$ _ = Quaternion a (V3 a a a)
+
+instance Applicative Quaternion where
+  pure a = Quaternion a (pure a)
+  Quaternion f fv <*> Quaternion a v = Quaternion (f a) (fv <*> v)
+
+instance Monad Quaternion where
+  return = pure
+  (>>=) = bindRep -- the diagonal of a sedenion is super useful!
+
+instance Representable Quaternion where
+  rep f = Quaternion (f _e) (V3 (f _i) (f _j) (f _k))
+
+instance Foldable Quaternion where
+  foldMap f (Quaternion e v) = f e `mappend` foldMap f v
+  foldr f z (Quaternion e v) = f e (F.foldr f z v)
+
+instance Traversable Quaternion where
+  traverse f (Quaternion e v) = Quaternion <$> f e <*> traverse f v
+
+instance forall a. Storable a => Storable (Quaternion a) where
+  sizeOf _ = 4 * sizeOf (undefined::a)
+  alignment _ = alignment (undefined::a)
+  poke ptr (Quaternion e v) = poke (castPtr ptr) e >>
+                              poke (castPtr (ptr `plusPtr` sz)) v
+    where sz = sizeOf (undefined::a)
+  peek ptr = Quaternion <$> peek (castPtr ptr)
+                        <*> peek (castPtr (ptr `plusPtr` sz))
+    where sz = sizeOf (undefined::a)
+
+instance RealFloat a => Num (Quaternion a) where
+  {-# SPECIALIZE instance Num (Quaternion Float) #-}
+  {-# SPECIALIZE instance Num (Quaternion Double) #-}
+  (+) = liftA2 (+)
+  (-) = liftA2 (-)
+  negate = fmap negate
+  Quaternion s1 v1 * Quaternion s2 v2 = Quaternion (s1*s2 - (v1 `dot` v2)) $
+                                        (v1 `cross` v2) + s1*^v2 + s2*^v1
+  fromInteger x = Quaternion (fromInteger x) 0
+  abs z = Quaternion (norm z) 0
+  signum q@(Quaternion e (V3 i j k))
+    | m == 0.0 = q
+    | not (isInfinite m || isNaN m) = q ^/ sqrt m
+    | any isNaN q = qNaN
+    | not (ii || ij || ik) = Quaternion 1 (V3 0 0 0)
+    | not (ie || ij || ik) = Quaternion 0 (V3 1 0 0)
+    | not (ie || ii || ik) = Quaternion 0 (V3 0 1 0)
+    | not (ie || ii || ij) = Quaternion 0 (V3 0 0 1)
+    | otherwise = qNaN
+    where
+      m = quadrance q
+      ie = isInfinite e
+      ii = isInfinite i
+      ij = isInfinite j
+      ik = isInfinite k
+
+  -- abs    = error "Quaternion.abs: use norm"
+  -- signum = error "Quaternion.signum: use signorm"
+
+qNaN :: RealFloat a => Quaternion a
+qNaN = Quaternion fNaN (V3 fNaN fNaN fNaN) where fNaN = 0/0
+
+-- {-# RULES "abs/norm" abs x = Quaternion (norm x) 0 #-}
+-- {-# RULES "signum/signorm" signum = signorm #-}
+
+-- this will attempt to rewrite calls to abs to use norm intead when it is available.
+
+instance RealFloat a => Fractional (Quaternion a) where
+  {-# SPECIALIZE instance Fractional (Quaternion Float) #-}
+  {-# SPECIALIZE instance Fractional (Quaternion Double) #-}
+  Quaternion q0 (V3 q1 q2 q3) / Quaternion r0 (V3 r1 r2 r3) =
+    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 = q ^/ quadrance q
+  fromRational x = Quaternion (fromRational x) 0
+
+instance Metric Quaternion where
+  Quaternion e v `dot` Quaternion e' v' = e*e' + (v `dot` v')
+
+class Complicated t where
+  _e :: Functor f => (a -> f a) -> t a -> f (t a)
+  _i :: Functor f => (a -> f a) -> t a -> f (t a)
+
+instance Complicated Complex where
+  _e f (a :+ b) = (:+ b) <$> f a
+  _i f (a :+ b) = (a :+) <$> f b
+
+instance Complicated Quaternion where
+  _e f (Quaternion a v) = (\a' -> Quaternion a' v) <$> f a
+  _i f (Quaternion a v) = Quaternion a <$> traverseOf _x f v
+  --_i f (Quaternion a (V3 b c d)) = (\b' -> Quaternion a (V3 b' c d)) <$> f b
+
+class Complicated t => Hamiltonian t where
+  _j :: Functor f => (a -> f a) -> t a -> f (t a)
+  _k :: Functor f => (a -> f a) -> t a -> f (t a)
+  _ijk :: Functor f => (V3 a -> f (V3 a)) -> t a -> f (t a)
+
+instance Hamiltonian Quaternion where
+  _j f (Quaternion a v) = Quaternion a <$> traverseOf _y f v
+  _k f (Quaternion a v) = Quaternion a <$> traverseOf _z f v
+  -- _j f (Quaternion a (V3 b c d)) = (\c' -> Quaternion a (V3 b c' d)) <$> f c
+  -- _k f (Quaternion a (V3 b c d)) = Quaternion a . V3 b c <$> f d
+
+  _ijk f (Quaternion a v) = Quaternion a <$> f v
+
+instance Distributive Quaternion where
+  distribute = distributeRep
+
+instance (Conjugate a, RealFloat a) => Conjugate (Quaternion a) where
+  conjugate (Quaternion e v) = Quaternion (conjugate e) (negate v)
+
+reimagine :: RealFloat a => a -> a -> Quaternion a -> Quaternion a
+reimagine r s (Quaternion _ v)
+  | isNaN s || isInfinite s = let aux 0 = 0
+                                  aux x = s * x
+                              in Quaternion r (aux <$> v)
+  | otherwise = Quaternion r (v^*s)
+
+-- | quadrance of the imaginary component
+qi :: Num a => Quaternion a -> a
+qi (Quaternion _ v) = quadrance v
+
+-- | norm of the imaginary component
+absi :: Floating a => Quaternion a -> a
+absi = sqrt . qi
+
+-- | raise a 'Quaternion' to a scalar power
+pow :: RealFloat a => Quaternion a -> a -> Quaternion a
+pow q t = exp (t *^ log q)
+
+-- ehh..
+instance RealFloat a => Floating (Quaternion a) where
+  {-# SPECIALIZE instance Floating (Quaternion Float) #-}
+  {-# SPECIALIZE instance Floating (Quaternion Double) #-}
+  pi = Quaternion pi 0
+  exp q@(Quaternion e v)
+    | qiq == 0 = Quaternion (exp e) v
+    | ai <- sqrt qiq, ee <- exp e = reimagine (ee * cos ai) (ee * (sin ai / ai)) q
+    where qiq = qi q
+  log q@(Quaternion e v@(V3 _i j k))
+    | qiq == 0 = if e >= 0
+                 then Quaternion (log e) v
+                 else Quaternion (log (negate e)) (V3 pi j k) -- mmm, pi
+    | ai <- sqrt qiq, m <- sqrt (e*e + qiq) = reimagine (log m) (atan2 m e / ai) q
+    where qiq = qi q
+  x ** y = exp (y * log x)
+  sqrt q@(Quaternion e v)
+    | m   == 0 = q
+    | qiq == 0 = if e > 0
+                 then Quaternion (sqrt e) 0
+                 else Quaternion 0 (V3 (sqrt (negate e)) 0 0)
+    | im <- sqrt (0.5*(m-e)) / sqrt qiq = Quaternion (0.5*(m+e)) (v^*im)
+    where qiq = qi q
+          m = sqrt (e*e + qiq)
+  cos q@(Quaternion e v)
+    | qiq == 0 = Quaternion (cos e) v
+    | ai <- sqrt qiq = reimagine (cos e * cosh ai) (- sin e * (sinh ai / ai)) q
+    where qiq = qi q
+  sin q@(Quaternion e v)
+    | qiq == 0 = Quaternion (sin e) v
+    | ai <- sqrt qiq = reimagine (sin e * cosh ai) (cos e * (sinh ai / ai)) q
+    where qiq = qi q
+  tan q@(Quaternion e v)
+    | qiq == 0 = Quaternion (tan e) v
+    | ai <- sqrt qiq, ce <- cos e, sai <- sinh ai, d <- ce*ce + sai*sai =
+      reimagine (ce * sin e / d) (cosh ai * (sai / ai) / d) q
+    where qiq = qi q
+  sinh q@(Quaternion e v)
+    | qiq == 0 = Quaternion (sinh e) v
+    | ai <- sqrt qiq = reimagine (sinh e * cos ai) (cosh e * (sin ai / ai)) q
+    where qiq = qi q
+  cosh q@(Quaternion e v)
+    | qiq == 0 = Quaternion (cosh e) v
+    | ai <- sqrt qiq = reimagine (cosh e * cos ai) ((sinh e * sin ai) / ai) q
+    where qiq = qi q
+  tanh q@(Quaternion e v)
+    | qiq == 0 = Quaternion (tanh e) v
+    | ai <- sqrt qiq, se <- sinh e, cai <- cos ai, d <- se*se + cai*cai =
+      reimagine ((cosh e * se) / d) ((cai * (sin ai / ai)) / d) q
+    where qiq = qi q
+
+  asin q = cut asin q
+  acos q = cut acos q
+  atan q = cut atan q
+
+  asinh q = cut asinh q
+  acosh q = cut acosh q
+  atanh q = cut atanh q
+
+-- | Helper for calculating with specific branch cuts
+cut :: RealFloat a => (Complex a -> Complex a) -> Quaternion a -> Quaternion a
+cut f q@(Quaternion e v)
+  | qiq == 0 = Quaternion a (_x.~b$v)
+  | otherwise = reimagine a (b / ai) q
+  where qiq = qi q
+        ai = sqrt qiq
+        a :+ b = f (e :+ ai)
+
+-- | Helper for calculating with specific branch cuts
+cutWith :: RealFloat a => Complex a -> Quaternion a -> Quaternion a
+cutWith (r :+ im) q@(Quaternion e v)
+  | e /= 0 || qiq == 0 || isNaN qiq || isInfinite qiq = error "bad cut"
+  | s <- im / sqrt qiq = Quaternion r (v^*s)
+  where qiq = qi q
+
+-- | 'asin' with a specified branch cut.
+asinq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
+asinq q@(Quaternion e _) u
+  | qiq /= 0.0 || e >= -1 && e <= 1 = asin q
+  | otherwise = cutWith (asin (e :+ sqrt qiq)) u
+  where qiq = qi q
+
+-- | 'acos' with a specified branch cut.
+acosq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
+acosq q@(Quaternion e _) u
+  | qiq /= 0.0 || e >= -1 && e <= 1 = acos q
+  | otherwise = cutWith (acos (e :+ sqrt qiq)) u
+  where qiq = qi q
+
+-- | 'atan' with a specified branch cut.
+atanq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
+atanq q@(Quaternion e _) u
+  | e /= 0.0 || qiq >= -1 && qiq <= 1 = atan q
+  | otherwise = cutWith (atan (e :+ sqrt qiq)) u
+  where qiq = qi q
+
+-- | 'asinh' with a specified branch cut.
+asinhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
+asinhq q@(Quaternion e _) u
+  | e /= 0.0 || qiq >= -1 && qiq <= 1 = asinh q
+  | otherwise = cutWith (asinh (e :+ sqrt qiq)) u
+  where qiq = qi q
+
+-- | 'acosh' with a specified branch cut.
+acoshq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
+acoshq q@(Quaternion e _) u
+  | qiq /= 0.0 || e >= 1 = asinh q
+  | otherwise = cutWith (acosh (e :+ sqrt qiq)) u
+  where qiq = qi q
+
+-- | 'atanh' with a specified branch cut.
+atanhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
+atanhq q@(Quaternion e _) u
+  | qiq /= 0.0 || e > -1 && e < 1 = atanh q
+  | otherwise = cutWith (atanh (e :+ sqrt qiq)) u
+  where qiq = qi q
+
+-- | Spherical linear interpolation between two quaternions.
+slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a
+slerp q p t
+  | 1.0 - cosphi < 1e-8 = q
+  | phi <- acos cosphi, r <- recip (sin phi)
+  = (sin ((1-t)*phi)*r *^ q ^+^ f (sin (t*phi)*r) *^ p) ^/ sin phi
+  where
+   dqp = dot q p
+   (cosphi, f) = if dqp < 0 then (-dqp, negate) else (dqp, id)
+{-# SPECIALIZE slerp :: Quaternion Float -> Quaternion Float -> Float -> Quaternion Float #-}
+{-# SPECIALIZE slerp :: Quaternion Double -> Quaternion Double -> Double -> Quaternion Double #-}
+
+--slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a
+--slerp q0 q1 = let q10 = q1 / q0 in \t -> pow q10 t * q0
+
+-- | Apply a rotation to a vector.
+rotate :: (Conjugate a, RealFloat a) => Quaternion a -> V3 a -> V3 a
+rotate q v = (q * Quaternion 0 v * conjugate q)^._ijk
+
+{-
+rotate :: Num a => Quaternion a -> V3 a -> V3 a
+rotate (Quaternion a' b c d) (V3 x y z) = V3
+  (2*((t8+t10)*x+(t6- t4)*y+(t3+t7)*z)+x)
+  (2*((t4+ t6)*y+(t5+t10)*y+(t9-t2)*z)+y)
+  (2*((t7- t3)*z+(t2+ t9)*z+(t5+t8)*z)+z)
+  where
+    a = -a'
+    t2 = a*b
+    t3 = a*c
+    t4 = a*d
+    t5 = -b*b
+    t6 = b*c
+    t7 = b*d
+    t8 = -c*c
+    t9 = c*d
+    t10 = -d*d
+-}
+{-# SPECIALIZE rotate :: Quaternion Float -> V3 Float -> V3 Float #-}
+{-# SPECIALIZE rotate :: Quaternion Double -> V3 Double -> V3 Double #-}
+
+instance (RealFloat a, Epsilon a) => Epsilon (Quaternion a) where
+  nearZero = nearZero . quadrance
+
+-- | @'axisAngle' axis theta@ builds a 'Quaternion' representing a
+-- rotation of @theta@ radians about @axis@.
+axisAngle :: (Epsilon a, Floating a) => V3 a -> a -> Quaternion a
+axisAngle axis theta = normalize $ Quaternion (cos half) $ (sin half) *^ axis
+  where half = theta / 2
diff --git a/src/Linear/V2.hs b/src/Linear/V2.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/V2.hs
@@ -0,0 +1,92 @@
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+-- {-# OPTIONS_GHC -fno-warn-name-shadowing #-}
+module Linear.V2
+  ( V2(..)
+  , R2(..)
+  , perp
+  ) where
+
+import Control.Applicative
+import Control.Lens
+import Data.Data
+import Data.Distributive
+import Data.Foldable
+import Data.Monoid
+import Foreign.Ptr (castPtr)
+import Foreign.Storable (Storable(..))
+import Linear.Metric
+import Linear.Epsilon
+
+-- | A 2-dimensional vector
+data V2 a = V2 a a deriving (Eq,Ord,Show,Read,Data,Typeable)
+
+instance Functor V2 where
+  fmap f (V2 a b) = V2 (f a) (f b)
+
+instance Foldable V2 where
+  foldMap f (V2 a b) = f a `mappend` f b
+
+instance Traversable V2 where
+  traverse f (V2 a b) = V2 <$> f a <*> f b
+
+instance Applicative V2 where
+  pure a = V2 a a
+  V2 a b <*> V2 d e = V2 (a d) (b e)
+
+instance Monad V2 where
+  return a = V2 a a
+  (>>=) = bindRep
+
+instance Num a => Num (V2 a) where
+  (+) = liftA2 (+)
+  (*) = liftA2 (*)
+  negate = fmap negate
+  abs = fmap abs
+  signum = fmap signum
+  fromInteger = pure . fromInteger
+
+instance Fractional a => Fractional (V2 a) where
+  recip = fmap recip
+  (/) = liftA2 (/)
+  fromRational = pure . fromRational
+
+instance Metric V2 where
+  dot (V2 a b) (V2 c d) = a * c + b * d
+
+-- | A space that distinguishes 2 orthogonal basis vectors '_x' and '_y', but may have more.
+class R2 t where
+  _x :: Functor f => (a -> f a) -> t a -> f (t a)
+  _x = _xy._x
+
+  _y :: Functor f => (a -> f a) -> t a -> f (t a)
+  _y = _xy._y
+
+  _xy :: Functor f => (V2 a -> f (V2 a)) -> t a -> f (t a)
+
+instance R2 V2 where
+  _x f (V2 a b) = (`V2` b) <$> f a
+  _y f (V2 a b) = (V2 a) <$> f b
+  _xy = id
+
+instance Representable V2 where
+  rep f = V2 (f _x) (f _y)
+
+instance Distributive V2 where
+  distribute f = V2 (fmap (^._x) f) (fmap (^._y) f)
+
+-- | the counter-clockwise perpendicular vector
+perp :: Num a => V2 a -> V2 a
+perp (V2 a b) = V2 (negate b) a
+
+instance Epsilon a => Epsilon (V2 a) where
+  nearZero = nearZero . quadrance
+
+instance forall a. Storable a => Storable (V2 a) where
+  sizeOf _ = 2 * sizeOf (undefined::a)
+  alignment _ = alignment (undefined::a)
+  poke ptr (V2 x y) = poke ptr' x >> pokeElemOff ptr' 1 y
+    where ptr' = castPtr ptr
+  peek ptr = V2 <$> peek ptr' <*> peekElemOff ptr' 1
+    where ptr' = castPtr ptr
diff --git a/src/Linear/V3.hs b/src/Linear/V3.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/V3.hs
@@ -0,0 +1,98 @@
+{-# LANGUAGE DeriveDataTypeable, ScopedTypeVariables #-}
+module Linear.V3
+  ( V3(..)
+  , cross, triple
+  , R2(..)
+  , R3(..)
+  ) where
+
+import Control.Applicative
+import Control.Lens
+import Data.Data
+import Data.Distributive
+import Data.Foldable
+import Data.Monoid
+import Foreign.Ptr (castPtr)
+import Foreign.Storable (Storable(..))
+import Linear.Epsilon
+import Linear.Metric
+import Linear.V2
+
+-- | A 3-dimensional vector
+data V3 a = V3 a a a deriving (Eq,Ord,Show,Read,Data,Typeable)
+
+instance Functor V3 where
+  fmap f (V3 a b c) = V3 (f a) (f b) (f c)
+
+instance Foldable V3 where
+  foldMap f (V3 a b c) = f a `mappend` f b `mappend` f c
+
+instance Traversable V3 where
+  traverse f (V3 a b c) = V3 <$> f a <*> f b <*> f c
+
+instance Applicative V3 where
+  pure a = V3 a a a
+  V3 a b c <*> V3 d e f = V3 (a d) (b e) (c f)
+
+instance Monad V3 where
+  return a = V3 a a a
+  (>>=) = bindRep
+
+instance Num a => Num (V3 a) where
+  (+) = liftA2 (+)
+  (*) = liftA2 (*)
+  negate = fmap negate
+  abs = fmap abs
+  signum = fmap signum
+  fromInteger = pure . fromInteger
+
+instance Fractional a => Fractional (V3 a) where
+  recip = fmap recip
+  (/) = liftA2 (/)
+  fromRational = pure . fromRational
+
+instance Metric V3 where
+  dot (V3 a b c) (V3 d e f) = a * d + b * e + c * f
+  {-# INLINABLE dot #-}
+
+instance Distributive V3 where
+  distribute f = V3 (fmap (^._x) f) (fmap (^._y) f) (fmap (^._z) f)
+
+-- | A space that distinguishes 3 orthogonal basis vectors: '_x', '_y', and '_z'. (It may have more)
+class R2 t => R3 t where
+  _z :: Functor f => (a -> f a) -> t a -> f (t a)
+  _xyz :: Functor f => (V3 a -> f (V3 a)) -> t a -> f (t a)
+
+instance R2 V3 where
+  _x f (V3 a b c) = (\a' -> V3 a' b c) <$> f a
+  _y f (V3 a b c) = (\b' -> V3 a b' c) <$> f b
+  _xy f (V3 a b c) = (\(V2 a' b') -> V3 a' b' c) <$> f (V2 a b)
+
+instance R3 V3 where
+  _z f (V3 a b c) = V3 a b <$> f c
+  _xyz = id
+
+instance Representable V3 where
+  rep f = V3 (f _x) (f _y) (f _z)
+
+instance forall a. Storable a => Storable (V3 a) where
+  sizeOf _ = 3 * sizeOf (undefined::a)
+  alignment _ = alignment (undefined::a)
+  poke ptr (V3 x y z) = do poke ptr' x
+                           pokeElemOff ptr' 1 y
+                           pokeElemOff ptr' 2 z
+    where ptr' = castPtr ptr
+  peek ptr = V3 <$> peek ptr' <*> peekElemOff ptr' 1 <*> peekElemOff ptr' 2
+    where ptr' = castPtr ptr
+
+-- | cross product
+cross :: Num a => V3 a -> V3 a -> V3 a
+cross (V3 a b c) (V3 d e f) = V3 (b*f-c*e) (c*d-a*f) (a*e-b*d)
+{-# INLINABLE cross #-}
+
+-- | scalar triple product
+triple :: Num a => V3 a -> V3 a -> V3 a -> a
+triple a b c = dot a (cross b c)
+
+instance Epsilon a => Epsilon (V3 a) where
+  nearZero = nearZero . quadrance
diff --git a/src/Linear/V4.hs b/src/Linear/V4.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/V4.hs
@@ -0,0 +1,104 @@
+{-# LANGUAGE DeriveDataTypeable, ScopedTypeVariables #-}
+module Linear.V4
+  ( V4(..)
+  , vector, point
+  , R2(..)
+  , R3(..)
+  , R4(..)
+  ) where
+
+import Control.Applicative
+import Control.Lens
+import Data.Data
+import Data.Distributive
+import Data.Foldable
+import Data.Monoid
+import Foreign.Ptr (castPtr)
+import Foreign.Storable (Storable(..))
+import Linear.Epsilon
+import Linear.Metric
+import Linear.V2
+import Linear.V3
+
+-- | A 4-dimensional vector.
+data V4 a = V4 a a a a deriving (Eq,Ord,Show,Read,Data,Typeable)
+
+instance Functor V4 where
+  fmap f (V4 a b c d) = V4 (f a) (f b) (f c) (f d)
+
+instance Foldable V4 where
+  foldMap f (V4 a b c d) = f a `mappend` f b `mappend` f c `mappend` f d
+
+instance Traversable V4 where
+  traverse f (V4 a b c d) = V4 <$> f a <*> f b <*> f c <*> f d
+
+instance Applicative V4 where
+  pure a = V4 a a a a
+  V4 a b c d <*> V4 e f g h = V4 (a e) (b f) (c g) (d h)
+
+instance Monad V4 where
+  return a = V4 a a a a
+  (>>=) = bindRep
+
+instance Num a => Num (V4 a) where
+  (+) = liftA2 (+)
+  (*) = liftA2 (*)
+  negate = fmap negate
+  abs = fmap abs
+  signum = fmap signum
+  fromInteger = pure . fromInteger
+
+instance Fractional a => Fractional (V4 a) where
+  recip = fmap recip
+  (/) = liftA2 (/)
+  fromRational = pure . fromRational
+
+instance Metric V4 where
+  dot (V4 a b c d) (V4 e f g h) = a * e + b * f + c * g + d * h
+
+instance Distributive V4 where
+  distribute f = V4 (fmap (^._x) f) (fmap (^._y) f) (fmap (^._z) f) (fmap (^._w) f)
+
+-- | A space that distinguishes orthogonal basis vectors '_x', '_y', '_z', '_w'. (It may have more.)
+class R3 t => R4 t where
+  _w :: Functor f => (a -> f a) -> t a -> f (t a)
+  _xyzw :: Functor f => (V4 a -> f (V4 a)) -> t a -> f (t a)
+
+instance R2 V4 where
+  _x f (V4 a b c d) = (\a' -> V4 a' b c d) <$> f a
+  _y f (V4 a b c d) = (\b' -> V4 a b' c d) <$> f b
+  _xy f (V4 a b c d) = (\(V2 a' b') -> V4 a' b' c d) <$> f (V2 a b)
+
+instance R3 V4 where
+  _z f (V4 a b c d) = (\c' -> V4 a b c' d) <$> f c
+  _xyz f (V4 a b c d) = (\(V3 a' b' c') -> V4 a' b' c' d) <$> f (V3 a b c)
+
+instance R4 V4 where
+  _w f (V4 a b c d) = V4 a b c <$> f d
+  _xyzw = id
+
+instance Representable V4 where
+  rep f = V4 (f _x) (f _y) (f _z) (f _w)
+
+instance forall a. Storable a => Storable (V4 a) where
+  sizeOf _ = 4 * sizeOf (undefined::a)
+  alignment _ = alignment (undefined::a)
+  poke ptr (V4 x y z w) = do poke ptr' x
+                             pokeElemOff ptr' 1 y
+                             pokeElemOff ptr' 2 z
+                             pokeElemOff ptr' 3 w
+    where ptr' = castPtr ptr
+  peek ptr = V4 <$> peek ptr' <*> peekElemOff ptr' 1
+                <*> peekElemOff ptr' 2 <*> peekElemOff ptr' 3
+    where ptr' = castPtr ptr
+
+-- | Convert a 3-dimensional affine vector into a 4-dimensional homogeneous vector.
+vector :: Num a => V3 a -> V4 a
+vector (V3 a b c) = V4 a b c 0
+
+-- | Convert a 3-dimensional affine point into a 4-dimensional homogeneous vector.
+point :: Num a => V3 a -> V4 a
+point (V3 a b c) = V4 a b c 1
+
+instance Epsilon a => Epsilon (V4 a) where
+  nearZero = nearZero . quadrance
diff --git a/src/Linear/Vector.hs b/src/Linear/Vector.hs
new file mode 100644
--- /dev/null
+++ b/src/Linear/Vector.hs
@@ -0,0 +1,60 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Linear.Epsilon
+-- Copyright   :  (C) 2012 Edward Kmett
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Operations on free vector spaces.
+-----------------------------------------------------------------------------
+module Linear.Vector
+  ( (^+^)
+  , gnegate
+  , (^-^)
+  , (^*)
+  , (*^)
+  , (^/)
+  , lerp
+  ) where
+
+import Control.Applicative
+
+infixl 6 ^+^, ^-^
+infixl 7 ^*, *^, ^/
+
+-- | Compute the sum of two vectors
+(^+^) :: (Applicative f, Num a) => f a -> f a -> f a
+(^+^) = liftA2 (+)
+{-# INLINE (^+^) #-}
+
+-- | Compute the negation of a vector
+gnegate :: (Functor f, Num a) => f a -> f a
+gnegate = fmap negate
+{-# INLINE gnegate #-}
+
+-- | Compute the difference between two vectors
+(^-^) :: (Applicative f, Num a) => f a -> f a -> f a
+(^-^) = liftA2 (-)
+{-# INLINE (^-^) #-}
+
+-- | Compute the left scalar product
+(*^) :: (Functor f, Num a) => a -> f a -> f a
+(*^) a = fmap (a*)
+{-# INLINE (*^) #-}
+
+-- | Compute the right scalar product
+(^*) :: (Functor f, Num a) => f a -> a -> f a
+f ^* a = fmap (*a) f
+{-# INLINE (^*) #-}
+
+-- | Compute division by a scalar on the right.
+(^/) :: (Functor f, Fractional a) => f a -> a -> f a
+f ^/ a = fmap (/a) f
+{-# INLINE (^/) #-}
+
+-- | Linearly interpolate between two vectors.
+lerp :: (Applicative f, Num a) => a -> f a -> f a -> f a
+lerp alpha u v = alpha *^ u ^+^ (1 - alpha) *^ v
+{-# INLINE lerp #-}
diff --git a/tests/doctests.hs b/tests/doctests.hs
new file mode 100644
--- /dev/null
+++ b/tests/doctests.hs
@@ -0,0 +1,28 @@
+module Main where
+
+import Test.DocTest
+import System.Directory
+import System.FilePath
+import Control.Applicative
+import Control.Monad
+import Data.List
+
+main :: IO ()
+main = getSources >>= \sources -> doctest $
+    "-isrc"
+  : "-idist/build/autogen"
+  : "-optP-include"
+  : "-optPdist/build/autogen/cabal_macros.h"
+  : sources
+
+getSources :: IO [FilePath]
+getSources = filter (isSuffixOf ".hs") <$> go "src"
+  where
+    go dir = do
+      (dirs, files) <- getFilesAndDirectories dir
+      (files ++) . concat <$> mapM go dirs
+
+getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath])
+getFilesAndDirectories dir = do
+  c <- map (dir </>) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir
+  (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c