diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2012, Jason Dagit
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Jason Dagit nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/lin-alg.cabal b/lin-alg.cabal
new file mode 100644
--- /dev/null
+++ b/lin-alg.cabal
@@ -0,0 +1,25 @@
+name:                lin-alg
+version:             0.1.0.0
+synopsis:            Low-dimensional matrices and vectors for graphics and physics
+-- description:         
+license:             BSD3
+license-file:        LICENSE
+author:              Jason Dagit
+maintainer:          dagitj@gmail.com
+copyright:           2012
+category:            Math
+build-type:          Simple
+cabal-version:       >=1.8
+
+library
+  exposed-modules:     Numeric.LinearAlgebra.Matrix
+                     , Numeric.LinearAlgebra.Matrix.Class
+                     , Numeric.LinearAlgebra.Matrix.Mat44
+                     , Numeric.LinearAlgebra.Matrix.Mat33
+                     , Numeric.LinearAlgebra.Vector
+                     , Numeric.LinearAlgebra.OrthoNormalBasis
+  -- other-modules:       
+  build-depends:       base ==4.5.*
+                     , NumInstances == 1.*
+  hs-source-dirs:      src
+  ghc-options: -Wall -O2 -funbox-strict-fields
diff --git a/src/Numeric/LinearAlgebra/Matrix.hs b/src/Numeric/LinearAlgebra/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearAlgebra/Matrix.hs
@@ -0,0 +1,14 @@
+module Numeric.LinearAlgebra.Matrix
+( module C
+, module M44
+, module M33
+, Mat44
+, Mat33
+) where
+
+import Numeric.LinearAlgebra.Matrix.Class as C
+import Numeric.LinearAlgebra.Matrix.Mat44 as M44 hiding (Mat44(..))
+import Numeric.LinearAlgebra.Matrix.Mat44 (Mat44)
+import Numeric.LinearAlgebra.Matrix.Mat33 as M33 hiding (Mat33(..))
+import Numeric.LinearAlgebra.Matrix.Mat33 (Mat33)
+
diff --git a/src/Numeric/LinearAlgebra/Matrix/Class.hs b/src/Numeric/LinearAlgebra/Matrix/Class.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearAlgebra/Matrix/Class.hs
@@ -0,0 +1,39 @@
+module Numeric.LinearAlgebra.Matrix.Class where
+
+-- import Numeric.LinearAlgebra.Vector
+
+class Functor m => Matrix m where
+  mDim     :: m a -> Int
+  mElement :: m a -> Int -> Int -> a
+  mIndexOf :: (Ord a) => (a -> a -> Bool) -> m a -> (Int, Int)
+  mZip     :: (a -> b -> c) -> m a -> m b -> m c
+  -- | mFold is foldl1'
+  mFold    :: (a -> a -> a) -> m a -> a
+  det      :: Num a => m a -> a
+
+{-# INLINE mApply #-}
+mApply :: Functor f => f (a -> b) -> a -> f b
+mApply f m = fmap ($ m) f
+
+(.+.) :: (Num k, Matrix m) => m k -> m k -> m k
+(.+.) = mZip (+)
+(.-.) :: (Num k, Matrix m) => m k -> m k -> m k
+(.-.) = mZip (-)
+
+{-
+(.*.) :: (k ~ Element m, Num k, Matrix m) => m -> m -> m
+m .*. n = mIdxMap m $ \i j -> sum [ mElement m i k * mElement n k j | k <- [ 0 .. 3 ] ]
+-}
+{-
+(.*>) :: (k ~ Element m, k ~ Scalar v, Num k, Matrix m, Vector v)
+      => m -> v -> v
+m .*> v | vDim v == mDim m = flip vIdxMap v $ \k -> sum [ mElement m i k * vElement v k | i <- [ 0 .. mDim m ] ]
+        | otherwise        = error "Dimensions do not match"
+-}
+-- (*.) :: (k ~ Element m, Num k, Matrix m) => m -> k -> m
+-- m *. k = mMap (k*) m
+
+{-
+transpose :: Matrix m => m -> m
+transpose m = flip mIdxMap m $ \i j -> mElement j i m
+-}
diff --git a/src/Numeric/LinearAlgebra/Matrix/Mat33.hs b/src/Numeric/LinearAlgebra/Matrix/Mat33.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearAlgebra/Matrix/Mat33.hs
@@ -0,0 +1,122 @@
+module Numeric.LinearAlgebra.Matrix.Mat33 where
+
+import Numeric.LinearAlgebra.Matrix.Class
+
+import Data.List ( foldl1' )
+
+data Mat33 a = Mat33
+  { m00 :: !a, m01 :: !a, m02 :: !a
+  , m10 :: !a, m11 :: !a, m12 :: !a
+  , m20 :: !a, m21 :: !a, m22 :: !a
+  }
+  deriving (Read, Show, Eq, Ord)
+
+instance Matrix Mat33 where
+  {-# INLINE mDim #-}
+  mDim _ = 3
+  {-# INLINE mElement #-}
+  mElement m 0 0 = m00 m
+  mElement m 0 1 = m01 m
+  mElement m 0 2 = m02 m
+  mElement m 1 0 = m10 m
+  mElement m 1 1 = m11 m
+  mElement m 1 2 = m12 m
+  mElement m 2 0 = m20 m
+  mElement m 2 1 = m21 m
+  mElement m 2 2 = m22 m
+  mElement _ i j = error ("Index " ++ show i ++ ", " ++ show j ++ ": out of range, must be 0,0 to 2,2")
+  {-# INLINE mZip #-}
+  mZip f m n = Mat33 (f (m00 m) (m00 n)) (f (m01 m) (m01 n)) (f (m02 m) (m02 n))
+                     (f (m10 m) (m10 n)) (f (m11 m) (m11 n)) (f (m12 m) (m12 n))
+                     (f (m20 m) (m20 n)) (f (m21 m) (m21 n)) (f (m22 m) (m22 n))
+  {-# INLINE mFold #-}
+  mFold f m = foldl1' f [ mElement m i j | i <- [ 0 .. 2 ], j <- [ 0 .. 2 ] ]
+  {-# INLINE mIndexOf #-}
+  mIndexOf p m = fst (foldl1' p' [ ((i,j), mElement m i j) | j <- [ 2, 1 .. 0 ], i <- [ 2, 1 .. 0 ] ])
+    where
+    p' x@(_, a) y@(_, a') | a `p` a'  = x
+                          | otherwise = y
+  {-# INLINE det #-}
+  det (Mat33 a b c d e f g h i) = det33 a b c d e f g h i
+
+instance Functor Mat33 where
+  {-# INLINE fmap #-}
+  fmap f m = Mat33 (f (m00 m)) (f (m01 m)) (f (m02 m))
+                   (f (m10 m)) (f (m11 m)) (f (m12 m))
+                   (f (m20 m)) (f (m21 m)) (f (m22 m))
+
+{-# SPECIALIZE det33 :: Double -> Double -> Double
+                     -> Double -> Double -> Double
+                     -> Double -> Double -> Double
+                     -> Double #-}
+{-# SPECIALIZE det33 :: Float -> Float -> Float
+                     -> Float -> Float -> Float
+                     -> Float -> Float -> Float
+                     -> Float #-}
+det33 :: Num a
+      => a -> a -> a
+      -> a -> a -> a
+      -> a -> a -> a
+      -> a
+det33 a b c d e f g h i = a*e*i + d*h*c + g*b*f - g*e*c - d*b*i - a*h*f
+
+{- SPECIALIZE inv33 :: Mat33 Double -> Mat33 Double -}
+{- SPECIALIZE inv33 :: Mat33 Float  -> Mat33 Float -}
+{-
+inv33 :: Fractional a => Mat33 a -> Mat33 a
+inv33 m = m'
+  where
+  d  = det m
+  m' = Mat33 (m11 * m22 - m21 * m21)
+  m' = Mat44 (det3 (m11 m) (m12 m) (m13 m)
+                   (m21 m) (m22 m) (m23 m)
+                   (m31 m) (m32 m) (m33 m) / d)
+           (-(det3 (m01 m) (m02 m) (m03 m)
+                   (m21 m) (m22 m) (m23 m)
+                   (m31 m) (m32 m) (m33 m) / d))
+             (det3 (m01 m) (m02 m) (m03 m)
+                   (m11 m) (m12 m) (m13 m)
+                   (m31 m) (m32 m) (m33 m) / d)
+           (-(det3 (m01 m) (m02 m) (m03 m)
+                   (m11 m) (m12 m) (m13 m)
+                   (m21 m) (m22 m) (m23 m) / d))
+             
+           (-(det3 (m10 m) (m12 m) (m13 m)
+                   (m20 m) (m22 m) (m23 m)
+                   (m30 m) (m32 m) (m33 m) / d))
+             (det3 (m00 m) (m02 m) (m03 m)
+                   (m20 m) (m22 m) (m23 m)
+                   (m30 m) (m32 m) (m33 m) / d)
+           (-(det3 (m00 m) (m02 m) (m03 m)
+                   (m10 m) (m12 m) (m13 m)
+                   (m30 m) (m32 m) (m33 m) / d))
+             (det3 (m00 m) (m02 m) (m03 m)
+                   (m10 m) (m12 m) (m13 m)
+                   (m20 m) (m22 m) (m23 m) / d)
+
+             (det3 (m10 m) (m11 m) (m13 m)
+                   (m20 m) (m21 m) (m23 m)
+                   (m30 m) (m31 m) (m33 m) / d)
+           (-(det3 (m00 m) (m01 m) (m03 m)
+                   (m20 m) (m21 m) (m23 m)
+                   (m30 m) (m31 m) (m33 m) / d))
+             (det3 (m00 m) (m01 m) (m03 m)
+                   (m10 m) (m11 m) (m13 m)
+                   (m30 m) (m31 m) (m33 m) / d)
+           (-(det3 (m00 m) (m01 m) (m03 m)
+                   (m10 m) (m11 m) (m13 m)
+                   (m20 m) (m21 m) (m23 m) / d))
+
+           (-(det3 (m10 m) (m11 m) (m12 m)
+                   (m20 m) (m21 m) (m22 m)
+                   (m30 m) (m31 m) (m32 m) / d))
+             (det3 (m00 m) (m01 m) (m02 m)
+                   (m20 m) (m21 m) (m22 m)
+                   (m30 m) (m31 m) (m32 m) / d)
+           (-(det3 (m00 m) (m01 m) (m02 m)
+                   (m10 m) (m11 m) (m12 m)
+                   (m30 m) (m31 m) (m32 m) / d))
+             (det3 (m00 m) (m01 m) (m02 m)
+                   (m10 m) (m11 m) (m12 m)
+                   (m20 m) (m21 m) (m22 m) / d)
+-}
diff --git a/src/Numeric/LinearAlgebra/Matrix/Mat44.hs b/src/Numeric/LinearAlgebra/Matrix/Mat44.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearAlgebra/Matrix/Mat44.hs
@@ -0,0 +1,132 @@
+module Numeric.LinearAlgebra.Matrix.Mat44 where
+
+import Numeric.LinearAlgebra.Matrix.Class
+import Numeric.LinearAlgebra.Matrix.Mat33 ( det33 )
+
+import Data.List ( foldl1' )
+import Data.NumInstances() -- Add function instances to Num
+
+data Mat44 a = Mat44
+  { m00 :: !a, m01 :: !a, m02 :: !a, m03 :: !a
+  , m10 :: !a, m11 :: !a, m12 :: !a, m13 :: !a
+  , m20 :: !a, m21 :: !a, m22 :: !a, m23 :: !a
+  , m30 :: !a, m31 :: !a, m32 :: !a, m33 :: !a
+  }
+  deriving (Read, Show, Eq, Ord)
+
+instance Matrix Mat44 where
+  {-# INLINE mDim #-}
+  mDim _ = 4
+  {-# INLINE mElement #-}
+  mElement m 0 0 = m00 m
+  mElement m 0 1 = m01 m
+  mElement m 0 2 = m02 m
+  mElement m 0 3 = m03 m
+  mElement m 1 0 = m10 m
+  mElement m 1 1 = m11 m
+  mElement m 1 2 = m12 m
+  mElement m 1 3 = m13 m
+  mElement m 2 0 = m20 m
+  mElement m 2 1 = m21 m
+  mElement m 2 2 = m22 m
+  mElement m 2 3 = m23 m
+  mElement m 3 0 = m30 m
+  mElement m 3 1 = m31 m
+  mElement m 3 2 = m32 m
+  mElement m 3 3 = m33 m
+  mElement _ i j = error ("Index " ++ show i ++ ", " ++ show j ++ ": out of range, must be 0,0 to 3,3")
+  {-# INLINE mZip #-}
+  mZip f m n = Mat44 (f (m00 m) (m00 n)) (f (m01 m) (m01 n)) (f (m02 m) (m02 n)) (f (m03 m) (m03 n))
+                     (f (m10 m) (m10 n)) (f (m11 m) (m11 n)) (f (m12 m) (m12 n)) (f (m13 m) (m13 n))
+                     (f (m20 m) (m20 n)) (f (m21 m) (m21 n)) (f (m22 m) (m22 n)) (f (m23 m) (m23 n))
+                     (f (m30 m) (m30 n)) (f (m31 m) (m31 n)) (f (m32 m) (m32 n)) (f (m33 m) (m33 n))
+  {-# INLINE mFold #-}
+  mFold f m = foldl1' f [ mElement m i j | i <- [ 0 .. 3 ], j <- [ 0 .. 3 ] ]
+  {-# INLINE mIndexOf #-}
+  mIndexOf p m = fst (foldl1' p' [ ((i,j), mElement m i j) | j <- [ 3, 2 .. 0 ], i <- [ 3, 2 .. 0 ] ])
+    where
+    p' x@(_, a) y@(_, a') | a `p` a'  = x
+                          | otherwise = y
+  {-# INLINE det #-}
+  det = det44
+
+instance Functor Mat44 where
+  {-# INLINE fmap #-}
+  fmap f m = Mat44 (f (m00 m)) (f (m01 m)) (f (m02 m)) (f (m03 m))
+                   (f (m10 m)) (f (m11 m)) (f (m12 m)) (f (m13 m))
+                   (f (m20 m)) (f (m21 m)) (f (m22 m)) (f (m23 m))
+                   (f (m30 m)) (f (m31 m)) (f (m32 m)) (f (m33 m))
+
+{-# SPECIALIZE det44 :: Mat44 Double -> Double #-}
+{-# SPECIALIZE det44 :: Mat44 Float  -> Float #-}
+det44 :: Num a => Mat44 a -> a
+det44 = m00 * det33 m11 m12 m13
+                    m21 m22 m23
+                    m31 m32 m33
+      - m01 * det33 m10 m12 m13
+                    m20 m22 m23
+                    m30 m32 m33
+      + m02 * det33 m10 m11 m13
+                    m20 m21 m23
+                    m30 m31 m33
+      - m03 * det33 m10 m11 m12
+                    m20 m21 m22
+                    m30 m31 m32
+
+{-# SPECIALIZE inv44 :: Mat44 Double -> Mat44 Double #-}
+{-# SPECIALIZE inv44 :: Mat44 Float  -> Mat44 Float  #-}
+inv44 :: Fractional a => Mat44 a -> Mat44 a
+inv44 m = mApply m' m
+  where
+  d  = const $ det44 m
+  m' = Mat44 (det33 m11 m12 m13
+                    m21 m22 m23
+                    m31 m32 m33 / d)
+           (-(det33 m01 m02 m03
+                    m21 m22 m23
+                    m31 m32 m33 / d))
+             (det33 m01 m02 m03
+                    m11 m12 m13
+                    m31 m32 m33 / d)
+           (-(det33 m01 m02 m03
+                    m11 m12 m13
+                    m21 m22 m23 / d))
+               
+           (-(det33 m10 m12 m13
+                    m20 m22 m23
+                    m30 m32 m33 / d))
+             (det33 m00 m02 m03
+                    m20 m22 m23
+                    m30 m32 m33 / d)
+           (-(det33 m00 m02 m03
+                    m10 m12 m13
+                    m30 m32 m33 / d))
+             (det33 m00 m02 m03
+                    m10 m12 m13
+                    m20 m22 m23 / d)
+
+             (det33 m10 m11 m13
+                    m20 m21 m23
+                    m30 m31 m33 / d)
+           (-(det33 m00 m01 m03
+                    m20 m21 m23
+                    m30 m31 m33 / d))
+             (det33 m00 m01 m03
+                    m10 m11 m13
+                    m30 m31 m33 / d)
+           (-(det33 m00 m01 m03
+                    m10 m11 m13
+                    m20 m21 m23 / d))
+
+           (-(det33 m10 m11 m12
+                    m20 m21 m22
+                    m30 m31 m32 / d))
+             (det33 m00 m01 m02
+                    m20 m21 m22
+                    m30 m31 m32 / d)
+           (-(det33 m00 m01 m02
+                    m10 m11 m12
+                    m30 m31 m32 / d))
+             (det33 m00 m01 m02
+                    m10 m11 m12
+                    m20 m21 m22 / d)
diff --git a/src/Numeric/LinearAlgebra/OrthoNormalBasis.hs b/src/Numeric/LinearAlgebra/OrthoNormalBasis.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearAlgebra/OrthoNormalBasis.hs
@@ -0,0 +1,87 @@
+module Numeric.LinearAlgebra.OrthoNormalBasis where
+
+import Numeric.LinearAlgebra.Vector
+
+epsilon :: (Ord a, Floating a) => a
+epsilon = 0.01
+
+data ONB a = ONB
+  { onbU :: Vec3 a
+  , onbV :: Vec3 a
+  , onbW :: Vec3 a
+  } deriving (Read, Show, Eq, Ord)
+
+mkFromU :: (Ord a, Floating a) => Vec3 a -> ONB a
+mkFromU u' = ONB u v w
+  where
+  n = Vec3 1 0 0
+  m = Vec3 0 1 0
+  u = unitVector u'
+  v = if len (u <%> n) < epsilon
+        then u <%> m
+        else u <%> n
+  w = u <%> v
+
+mkFromV :: (Ord a, Floating a) => Vec3 a -> ONB a
+mkFromV v' = ONB u v w
+  where
+  n = Vec3 1 0 0
+  m = Vec3 0 1 0
+  v = unitVector v'
+  u = if lenSquared (v <%> n) < epsilon
+        then v <%> m
+        else v <%> n
+  w = u <%> v
+
+mkFromW :: (Ord a, Floating a) => Vec3 a -> ONB a
+mkFromW w' = ONB u v w
+  where
+  n = Vec3 1 0 0
+  m = Vec3 0 1 0
+  w = unitVector w'
+  u = if len (w <%> n) < epsilon
+        then w <%> m
+        else w <%> n
+  v = w <%> u
+
+mkFromUV :: Floating a => Vec3 a -> Vec3 a -> ONB a
+mkFromUV u' v' = ONB u v w
+  where
+  u = unitVector u'
+  w = unitVector (u' <%> v')
+  v = w <%> u
+
+mkFromVU :: Floating a => Vec3 a -> Vec3 a -> ONB a
+mkFromVU v' u' = ONB u v w
+  where
+  v = unitVector v'
+  w = unitVector (u' <%> v')
+  u = v <%> w
+
+mkFromUW :: Floating a => Vec3 a -> Vec3 a -> ONB a
+mkFromUW u' w' = ONB u v w
+  where
+  u = unitVector u'
+  v = unitVector (w' <%> u')
+  w = u <%> v
+
+mkFromWU :: Floating a => Vec3 a -> Vec3 a -> ONB a
+mkFromWU w' u' = ONB u v w
+  where
+  w = unitVector w'
+  v = unitVector (w' <%> u')
+  u = v <%> w
+
+mkFromVW :: Floating a => Vec3 a -> Vec3 a -> ONB a
+mkFromVW v' w' = ONB u v w
+  where
+  v = unitVector v'
+  u = unitVector (v' <%> w')
+  w = u <%> v
+
+mkFromWV :: Floating a => Vec3 a -> Vec3 a -> ONB a
+mkFromWV w' v' = ONB u v w
+  where
+  w = unitVector w'
+  u = unitVector (v' <%> w')
+  v = w <%> u 
diff --git a/src/Numeric/LinearAlgebra/Vector.hs b/src/Numeric/LinearAlgebra/Vector.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearAlgebra/Vector.hs
@@ -0,0 +1,147 @@
+module Numeric.LinearAlgebra.Vector where
+
+data Vec4 a = Vec4 !a !a !a !a
+  deriving (Read, Show, Eq, Ord)
+data Vec3 a = Vec3 !a !a !a
+  deriving (Read, Show, Eq, Ord)
+data Vec2 a = Vec2 !a !a
+  deriving (Read, Show, Eq, Ord)
+
+class Functor v => Vector v where
+  vDim     :: v a -> Int
+  vElement :: v a -> Int -> a
+  vIndexOf :: (Ord a) => (a -> a -> Bool) -> v a -> Int
+  vZip     :: (a -> b -> c) -> v a -> v b -> v c
+  -- | vFold is foldl1'
+  vFold    :: (a -> a -> a) -> v a -> a
+
+instance Vector Vec4 where
+  {-# INLINE vDim #-}
+  vDim _ = 4
+  {-# INLINE vElement #-}
+  vElement (Vec4 x _ _ _) 0 = x
+  vElement (Vec4 _ y _ _) 1 = y
+  vElement (Vec4 _ _ z _) 2 = z
+  vElement (Vec4 _ _ _ w) 3 = w
+  vElement _ i = error ("Index " ++ show i ++ ": out of range, must be 0 to 3")
+  {-# INLINE vZip #-}
+  vZip f (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (f x1 x2) (f y1 y2) (f z1 z2) (f w1 w2)
+  {-# INLINE vFold #-}
+  vFold f (Vec4 x y z w) = f (f (f x y) z) w
+  {-# INLINE vIndexOf #-}
+  vIndexOf p (Vec4 x y z w) | w `p` x && w `p` y && w `p` z = 3
+                            | z `p` x && z `p` y && z `p` w = 2
+                            | y `p` x && y `p` z && y `p` w = 1
+                            | otherwise = 0
+
+instance Functor Vec4 where
+  {-# INLINE fmap #-}
+  fmap f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w)
+
+instance Vector Vec3 where
+  {-# INLINE vDim #-}
+  vDim _ = 3
+  {-# INLINE vElement #-}
+  vElement (Vec3 x _ _) 0 = x
+  vElement (Vec3 _ y _) 1 = y
+  vElement (Vec3 _ _ z) 2 = z
+  vElement _ i = error ("Index " ++ show i ++ ": out of range, must be 0 to 2")
+  {-# INLINE vZip #-}
+  vZip f (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (f x1 x2) (f y1 y2) (f z1 z2)
+  {-# INLINE vFold #-}
+  vFold f (Vec3 x y z) = f (f x y) z
+  {-# INLINE vIndexOf #-}
+  vIndexOf p (Vec3 x y z) | z `p` x && z `p` y = 2
+                          | y `p` x && y `p` z = 1
+                          | otherwise = 0
+
+instance Functor Vec3 where
+  {-# INLINE fmap #-}
+  fmap f (Vec3 x y z) = Vec3 (f x) (f y) (f z)
+
+instance Vector Vec2 where
+  {-# INLINE vDim #-}
+  vDim _ = 2
+  {-# INLINE vElement #-}
+  vElement (Vec2 x _) 0 = x
+  vElement (Vec2 _ y) 1 = y
+  vElement _ i = error ("Index " ++ show i ++ ": out of range, must be 0 or 1")
+  {-# INLINE vZip #-}
+  vZip f (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (f x1 x2) (f y1 y2)
+  {-# INLINE vFold #-}
+  vFold f (Vec2 x y) = f x y
+  {-# INLINE vIndexOf #-}
+  vIndexOf p (Vec2 x y) | y `p` x   = 1
+                        | otherwise = 0
+
+instance Functor Vec2 where
+  {-# INLINE fmap #-}
+  fmap f (Vec2 x y) = Vec2 (f x) (f y)
+
+{-# INLINE vNegate #-}
+vNegate :: (Num a, Vector v) => v a -> v a
+vNegate = fmap negate
+{-# INLINE (<+>) #-}
+(<+>) :: (Num a, Vector v) => v a -> v a -> v a
+(<+>) = vZip (+)
+{-# INLINE (<*>) #-}
+(<*>) :: (Num a, Vector v) => v a -> v a -> v a
+(<*>) = vZip (*)
+{-# INLINE (<->) #-}
+(<->) :: (Num a, Vector v) => v a -> v a -> v a
+(<->) = vZip (-)
+{-# INLINE (</>) #-}
+(</>) :: (Fractional a, Vector v) => v a -> v a -> v a
+(</>) = vZip (/)
+{-# INLINE (<.>) #-}
+(<.>) :: (Num a, Vector v) => v a -> v a -> a -- dot product
+v1 <.> v2 = vFold (+) (v1 <*> v2)
+{-# INLINE (*>) #-}
+(*>) :: (Num a, Vector v) => a -> v a -> v a
+k *> v = fmap (k*) v
+{-# INLINE (</) #-}
+(</) :: (Fractional a, Vector v) => v a -> a -> v a
+v </ k = fmap (/k) v
+{-# INLINE len #-}
+len :: (Floating a, Vector v) => v a -> a
+len v = sqrt $ lenSquared v
+{-# INLINE lenSquared #-}
+lenSquared :: (Num a, Vector v) => v a -> a
+lenSquared v = v <.> v
+
+maxVec :: (Ord k, Vector v) => v k -> v k -> v k
+maxVec = vZip max
+minVec :: (Ord k, Vector v) => v k -> v k -> v k
+minVec = vZip min
+minComponent :: (Ord k, Vector v) => v k -> k
+minComponent = vFold min
+maxComponent :: (Ord k, Vector v) => v k -> k
+maxComponent = vFold max
+minAbsComponent :: (Num k, Ord k, Vector v) => v k -> k
+minAbsComponent = vFold (\x y -> min (abs x) (abs y))
+maxAbsComponent :: (Num k, Ord k, Vector v) => v k -> k
+maxAbsComponent = vFold (\x y -> max (abs x) (abs y))
+vIndexOfMinComponent :: (Ord k, Vector v) => v k -> Int
+vIndexOfMinComponent = vIndexOf (<)
+vIndexOfMaxComponent :: (Ord k, Vector v) => v k -> Int
+vIndexOfMaxComponent = vIndexOf (>)
+vIndexOfMinAbsComponent :: (Ord k, Num k, Vector v)
+                        => v k -> Int
+vIndexOfMinAbsComponent = vIndexOf (\x y -> abs x < abs y)
+vIndexOfMaxAbsComponent :: (Ord k, Num k, Vector v)
+                        => v k -> Int
+vIndexOfMaxAbsComponent = vIndexOf (\x y -> abs x > abs y)
+
+-- |Cross product
+{-# INLINE (<%>) #-}
+(<%>) :: Floating a => Vec3 a -> Vec3 a -> Vec3 a
+(Vec3 x1 y1 z1) <%> (Vec3 x2 y2 z2) = Vec3 (y1*z2 - z1*y2)
+                                           (z1*x2 - x1*z2)
+                                           (x1*y2 - y1*x2)
+
+{-# INLINE unitVector #-}
+unitVector :: (Floating a, Vector v) => v a -> v a
+unitVector v = v </ len v
+
+tripleProduct :: Floating a => Vec3 a -> Vec3 a -> Vec3 a -> a
+tripleProduct v1 v2 v3 = (v1 <%> v2) <.> v3