diff --git a/Data/Vec.hs b/Data/Vec.hs
new file mode 100644
--- /dev/null
+++ b/Data/Vec.hs
@@ -0,0 +1,137 @@
+{-
+Copyright (c) 2008, Scott E. Dillard
+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.
+
+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.
+-}
+
+{- |
+
+Vec : a library for fixed-length lists and low-dimensional linear algebra
+
+Scott E. Dillard <sedillard@gmail.com>
+
+darcs : <http://graphics.cs.ucdavis.edu/~sdillard/Vec>
+
+/Synopsis/
+
+Vectors are represented by lists with type-encoded lengths. The constructor is
+@:.@, which acts like a cons both at the value and type levels, with @()@
+taking the place of nil. So @x:.y:.z:.()@ is a 3d vector. The library provides
+a set of common list-like functions (map, fold, etc) for working with vectors.
+Built up from these functions are a small but useful set of linear algebra
+operations: matrix multiplication, determinants, solving linear systems,
+inverting matrices.
+
+/Design/
+
+
+* Simplicity : 
+Beyond the initial complexities of type-level lists and
+numbers, I've tried to keep the API simple. There is no vector-space
+class, nor a complicated hierarchy of linear\/affine\/projective
+transformations. These can be added on top of the library easily.
+
+* Purity :
+The library is written in the functional style. For most
+functions this does not hinder performance at all, but some I am still
+working on (Gaussian elimination) so if this library is a bottleneck you
+can easily drop down to C. 
+
+* Low Dimension :
+Although the dimensionality is limited only by what GHC
+will handle, the library is meant for 2,3 and 4 dimensions. For general
+linear algebra, check out the excellent hmatrix library and blas bindings.
+
+To the point of simplicity, vectors and matrices are instances of Num and
+Fractional.  All arithmetic is done component-wise and literals construct
+uniform vectors and matrices. There are many interesting projects aiming to
+overhaul Haskell's number classes, but for now the type of @(*)@ is @a -> a ->
+a@ so that's what we're working with. It is easy to incorporate this library
+into a more mathematically consistent class hierarchy (provided you can design
+one.) 
+
+The rule is simple : 
+  If the method is unary, it's a map. 
+  If it's binary, it's a zipWith.
+
+/Performance/
+
+@(:.)@ is strict in both arguments, but it is also polymorphic, so at runtime
+vectors will be realized as linked lists, albeit with less pattern matching.
+However the library provides packed representations for 2,3 and 4d vectors of
+Ints, Floats and Doubles. @'Vec3F' x y z@ constructs a packed vector of
+unboxed Floats. Functions @'pack'@ and @'unpack'@ convert between packed and
+unpacked types. When vector operations are bracketed by 'pack' and 'unpack',
+GHC can unfold them into very efficient code. The 'Storable' instances for
+vectors also generate fast code.  Without optimizations, the code falls back
+into linked-list mode. The optimizations depend on inlining, so you may need
+to increase your unfolding threshold in certain situations.
+
+/GHC Extensions/
+
+This library makes heavy use of functional dependencies. I have tried to
+tweak things so that they \"just work.\" However, every now and then you will
+get incomprehensible error messages, usually about how this isn't an
+instance of that. These are how type errors typically manifest, so first
+double check to make sure you aren't trying to mix vectors of different
+dimension or component types. If you still get these errors, manual type
+annotations usually make them go away.
+
+
+/Related Work/
+
+See previous work by David Menendez,
+  <http://haskell.org/pipermail/haskell/2005-May/015815.html>
+
+and of course Oleg Kiselyov,
+  <http://okmij.org/ftp/papers/number-parameterized-types.pdf>
+
+Other vector and linear algebra packages :
+
+vector-space, by Conal Elliott : 
+  <http://hackage.haskell.org/cgi-bin/hackage-scripts/package/vector-space>
+
+hmatrix, by Alberto Ruiz :
+  <http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmatrix>
+
+blas bindings, by Patrick Perry :
+  <http://hackage.haskell.org/cgi-bin/hackage-scripts/package/blas>
+
+templatized geometry library (C++), by Oliver Kreylos :
+  <http://graphics.cs.ucdavis.edu/~okreylos/ResDev/Geometry/index.html>
+-}
+
+module Data.Vec 
+  (module Data.Vec.Base
+  ,module Data.Vec.LinAlg
+  ,module Data.Vec.Packed
+  ,module Data.Vec.Nat
+  )
+where
+
+import Data.Vec.Base
+import Data.Vec.LinAlg
+import Data.Vec.Packed
+import Data.Vec.Nat
+import Data.Vec.Instances
+
+
diff --git a/Data/Vec/Base.hs b/Data/Vec/Base.hs
new file mode 100644
--- /dev/null
+++ b/Data/Vec/Base.hs
@@ -0,0 +1,416 @@
+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}
+
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE NoMonomorphismRestriction #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+{-# HADDOCK_OPTIONS prune #-}
+
+module Data.Vec.Base where
+
+import Data.Vec.Nat
+
+import Prelude hiding (map,zipWith,foldl,foldr,reverse,
+                       take,drop,head,tail,sum,last,product,
+                       minimum,maximum)
+import qualified Prelude as P
+
+
+
+-- | The vector constructor. @(:.)@ for vectors is like @(:)@ for lists, and
+-- @()@ takes the place of @[]@. (The list of instances here is not meant to be
+-- readable.)
+
+data a :. b = !a :. !b
+  deriving (Eq,Ord,Read)
+
+infixr :.
+
+--derived show outputs in prefix notation
+instance (Show a, ShowVec v) => Show (a:.v) where
+  show (a:.v) = "(" ++ show a ++ ":." ++ showVec v ++ ")"
+
+
+-- | Helper to keep parentheses at bay. Just use @show@ as usual.
+class ShowVec  v where
+  showVec :: v -> String
+
+instance ShowVec () where
+  showVec = show
+  {-# INLINE showVec #-}
+
+instance (Show a, ShowVec v) => ShowVec (a:.v) where
+  showVec (a:.v) = show a ++ ":." ++ showVec v
+  {-# INLINE showVec #-}
+
+
+-- * Vector Types
+type Vec2  a = a :. a :. ()
+type Vec3  a = a :. (Vec2 a)
+type Vec4  a = a :. (Vec3 a)
+type Vec5  a = a :. (Vec4 a)
+type Vec6  a = a :. (Vec5 a)
+type Vec7  a = a :. (Vec6 a)
+type Vec8  a = a :. (Vec7 a)
+type Vec9  a = a :. (Vec8 a)
+type Vec10 a = a :. (Vec9 a)
+type Vec11 a = a :. (Vec10 a)
+type Vec12 a = a :. (Vec11 a)
+type Vec13 a = a :. (Vec12 a)
+type Vec14 a = a :. (Vec13 a)
+type Vec15 a = a :. (Vec14 a)
+type Vec16 a = a :. (Vec15 a)
+type Vec17 a = a :. (Vec16 a)
+type Vec18 a = a :. (Vec17 a)
+type Vec19 a = a :. (Vec18 a)
+
+
+
+
+-- | The type constraint @Vec n a v@ infers the vector type @v@ from the
+-- length @n@, a type-level natural, and underlying component type @a@.  
+-- So @x :: Vec N4 a v => v@ declares @x@ to be a 4-vector of @a@s.
+
+class Vec n a v | n a -> v, v -> n a where
+  -- | Make a uniform vector of a given length. @n@ is a type-level natural.
+  -- Use `vec` when the length can be inferred.
+  mkVec :: n -> a -> v
+
+  -- | turn a list into a vector of inferred length
+  fromList :: [a] -> v
+
+  -- | get a vector element, which one is determined at runtime
+  getElem :: Int -> v -> a
+
+  -- | set a vector element, which one is determined at runtime
+  setElem :: Int -> a -> v -> v
+
+instance Vec N1 a ( a :. () ) where
+  mkVec _ a = a :. ()
+  fromList (a:_)   = a :. ()
+  fromList []      = error "fromList: list too short"
+  getElem !i (a :. _) 
+    | i == 0    = a
+    | otherwise = error "getElem: index out of bounds"
+  setElem !i a _ 
+    | i == 0    = a :. ()
+    | otherwise = error "setElem: index out of bounds"
+  {-# INLINE setElem #-}
+  {-# INLINE getElem #-}
+  {-# INLINE mkVec #-}
+  {-# INLINE fromList #-}
+
+instance Vec (Succ n) a (a':.v) => Vec (Succ (Succ n)) a (a:.a':.v) where
+  mkVec _ a = a :. (mkVec undefined a)
+  fromList (a:as)  = a :. (fromList as)
+  fromList []      = error "fromList: list too short"
+  getElem !i (a :. v)
+    | i == 0    = a
+    | otherwise = getElem (i-1) v
+  setElem !i a (x :. v)
+    | i == 0    = a :. v
+    | otherwise = x :. (setElem (i-1) a v)
+  {-# INLINE setElem #-}
+  {-# INLINE getElem #-}
+  {-# INLINE mkVec #-}
+  {-# INLINE fromList #-}
+
+
+-- | Make a uniform vector. The length is inferred.
+vec ::  (Vec n a v) => a -> v
+vec = mkVec undefined
+{-# INLINE vec #-}
+
+
+-- | get or set a vector element, known at compile
+--time. Use the Nat types to access vector components. For instance, @get n0@
+--gets the x component, @set n2 44@ sets the z component to 44. 
+
+
+class Access n a v | v -> a where
+  get  :: n -> v -> a
+  set  :: n -> a -> v -> v
+
+instance Access N0 a (a :. v) where
+  get _ (a :. _) = a
+  set _ a (_ :. v) = a :. v
+  {-# INLINE set #-}
+  {-# INLINE get #-}
+
+instance Access n a v => Access (Succ n) a (a :. v) where
+  get _ (_ :. v) = get (undefined::n) v
+  set _ a' (a :. v) = a :. (set (undefined::n) a' v)
+  {-# INLINE set #-}
+  {-# INLINE get #-}
+
+
+-- * List-like functions
+
+-- | The first element.
+
+class Head v a | v -> a  where 
+  head :: v -> a
+
+instance Head (a :. as) a where 
+  head (a :. _) = a
+  {-# INLINE head #-}
+
+
+-- | All but the first element. 
+
+class Tail v v_ | v -> v_ where 
+  tail :: v -> v_
+
+instance Tail (a :. as) as where 
+  tail (_ :. as) = as
+  {-# INLINE tail #-}
+
+
+
+
+-- | Apply a function over each element in a vector. Constraint @Map a b u v@
+-- states that @u@ is a vector of @a@s, @v@ is a vector of @b@s with the same
+-- length as @u@, and the function is of type @a -> b@.
+
+class Map a b u v | u -> a, v -> b, b u -> v, a v -> u where
+  map :: (a -> b) -> u -> v
+
+instance Map a b (a :. ()) (b :. ()) where
+  map f (x :. ()) = (f x) :. ()
+  {-# INLINE map #-}
+
+instance Map a b (a':.u) (b':.v) => Map a b (a:.a':.u) (b:.b':.v) where
+  map f (x:.v) = (f x):.(map f v)
+  {-# INLINE map #-}
+
+
+
+
+-- | Combine two vectors using a binary function. The length of the result is
+-- the min of the lengths of the arguments. The constraint @ZipWith a b c u v
+-- w@ states that @u@ is a vector of @a@s, @v@ is a vector of @b@s, @w@ is a
+-- vector of @c@s, and the binary function is of type @a -> b -> c@.
+
+class ZipWith a b c u v w | u->a, v->b, w->c, u v c -> w where
+  zipWith :: (a -> b -> c) -> u -> v -> w
+
+instance ZipWith a b c (a:.()) (b:.()) (c:.()) where
+  zipWith f (x:._) (y:._) = f x y :.()
+  {-# INLINE zipWith #-}
+
+instance ZipWith a b c (a:.()) (b:.b:.bs) (c:.()) where
+  zipWith f (x:._) (y:._) = f x y :.()
+  {-# INLINE zipWith #-}
+
+instance ZipWith a b c (a:.a:.as) (b:.()) (c:.()) where
+  zipWith f (x:._) (y:._) = f x y :.()
+  {-# INLINE zipWith #-}
+
+instance 
+  ZipWith a b c (a':.u) (b':.v) (c':.w) 
+  => ZipWith a b c (a:.a':.u) (b:.b':.v) (c:.c':.w) 
+    where
+      zipWith f (x:.u) (y:.v) = f x y :. zipWith f u v
+      {-# INLINE zipWith #-}
+
+
+-- | Fold a function over a vector. 
+
+class Fold a v | v -> a where
+  fold  :: (a -> a -> a) -> v -> a
+  foldl :: (b -> a -> b) -> b -> v -> b
+  foldr :: (a -> b -> b) -> b -> v -> b
+
+instance Fold a (a:.()) where
+  fold  f   (a:._) = a 
+  foldl f z (a:._) = (f $! z) $! a
+  foldr f z (a:._) = (f $! a) $! z
+  {-# INLINE fold #-}
+  {-# INLINE foldl #-}
+  {-# INLINE foldr #-}
+
+instance Fold a (a':.u) => Fold a (a:.a':.u) where
+  fold  f   (a:.v) = (f $! a) $! (fold f v)
+  foldl f z (a:.v) = (f $! (foldl f z v)) $! a
+  foldr f z (a:.v) = (f $! a) $! (foldr f z v)
+  {-# INLINE fold #-}
+  {-# INLINE foldl #-}
+  {-# INLINE foldr #-}
+
+-- | Reverse a vector 
+reverse v = reverse' () v
+{-# INLINE reverse #-}
+
+-- Reverse helper function : builds the reversed list as its first argument
+class Reverse' p v v' | p v -> v' where
+  reverse' :: p -> v -> v'
+  
+instance Reverse' p () p where
+  reverse' p () = p
+  {-# INLINE reverse' #-}
+
+instance Reverse' (a:.p) v v' => Reverse' p (a:.v) v' where
+  reverse' p (a:.v) = reverse' (a:.p) v 
+  {-# INLINE reverse' #-}
+
+
+-- | Append two vectors 
+
+class Append v1 v2 v3 | v1 v2 -> v3, v1 v3 -> v2 where 
+  append :: v1 -> v2 -> v3
+
+instance Append () v v where
+  append _ = id
+  {-# INLINE append #-}
+
+instance Append (a:.()) v (a:.v) where
+  append (a:.()) v = a:.v
+  {-# INLINE append #-}
+
+instance (Append (a':.v1) v2 v3) => Append (a:.a':.v1) v2 (a:.v3) where
+  append (a:.u) v  =  a:.(append u v)
+  {-# INLINE append #-}
+
+
+
+-- | @take n v@ constructs a vector from the first @n@ elements of @v@. @n@ is a
+-- type-level natural. For example @take n3 v@ makes a 3-vector of the first
+-- three elements of @v@.
+
+class Take n v v' | n v -> v', n v' -> v where
+  take :: n -> v -> v'
+
+instance Take N0 v () where
+  take _ _ = ()
+  {-# INLINE take #-}
+
+instance Take n v v' => Take (Succ n) (a:.v) (a:.v') where
+  take _ (a:.v) = a:.(take (undefined::n) v)
+  {-# INLINE take #-}
+
+
+-- | @drop n v@ strips the first @n@ elements from @v@. @n@ is a type-level
+-- natural. For example @drop n2 v@ drops the first two elements.
+
+class Drop n v v' | n v -> v', n v' -> v where
+  drop :: n -> v -> v'
+ 
+instance Drop N0 v v where
+  drop _ = id
+  {-# INLINE drop #-}
+
+instance (Tail v' v'', Drop n v v') => Drop (Succ n) v v'' where
+  drop _ = tail . drop (undefined::n)
+  {-# INLINE drop #-}
+
+
+-- | Get the last element, usually significant for some reason (quaternions,
+-- homogenous coordinates, whatever)
+class Last v a | v -> a where
+  last :: v -> a
+
+instance Last (a:.()) a where 
+  last (a:._) = a
+  {-# INLINE last #-}
+
+instance Last (a':.v) a => Last (a:.a':.v) a where
+  last (a:.v) = last v
+  {-# INLINE last #-}
+
+-- | @snoc v a@ appends the element a to the end of v. 
+
+class Snoc v a v' | v a -> v', v' -> v a where 
+  snoc :: v -> a -> v'
+
+instance Snoc () a (a:.()) where
+  snoc _ a = (a:.())
+  {-# INLINE snoc #-}
+
+instance Snoc v a (a:.v) => Snoc (a:.v) a (a:.a:.v) where
+  snoc (b:.v) a = b:.(snoc v a)
+  {-# INLINE snoc #-}
+
+
+
+-- | sum of vector elements
+sum ::  (Fold a v, Num a) => v -> a
+sum x     = fold (+) x
+{-# INLINE sum #-}
+
+-- | product of vector elements
+product ::  (Fold a v, Num a) => v -> a
+product x = fold (*) x
+{-# INLINE product #-}
+
+-- | maximum vector element
+maximum ::  (Fold a v, Ord a) => v -> a
+maximum x = fold max x
+{-# INLINE maximum #-}
+
+-- | minimum vector element
+minimum ::  (Fold a v, Ord a) => v -> a
+minimum x = fold min x
+{-# INLINE minimum #-}
+
+toList ::  (Fold a v) => v -> [a]
+toList = foldr (:) [] 
+{-# INLINE toList #-}
+
+
+
+
+
+
+
+-- * Matrix Types
+
+type Mat22 a = Vec2 (Vec2 a)
+type Mat23 a = Vec2 (Vec3 a)
+type Mat24 a = Vec2 (Vec4 a)
+
+type Mat32 a = Vec3 (Vec2 a)
+type Mat33 a = Vec3 (Vec3 a)
+type Mat34 a = Vec3 (Vec4 a)
+type Mat35 a = Vec3 (Vec5 a)
+type Mat36 a = Vec3 (Vec6 a)
+
+type Mat42 a = Vec4 (Vec2 a)
+type Mat43 a = Vec4 (Vec3 a)
+type Mat44 a = Vec4 (Vec4 a)
+type Mat45 a = Vec4 (Vec5 a)
+type Mat46 a = Vec4 (Vec6 a)
+type Mat47 a = Vec4 (Vec7 a)
+type Mat48 a = Vec4 (Vec8 a)
+
+-- | convert a matrix to a list-of-lists
+matToLists ::  (Fold a v, Fold v m) => m -> [[a]]
+matToLists   = (P.map toList) . toList
+{-# INLINE matToLists   #-}
+
+-- | convert a matrix to a list in row-major order
+matToList  ::  (Fold a v, Fold v m) => m -> [a]
+matToList    = concat . matToLists
+{-# INLINE matToList    #-}
+
+-- | convert a list-of-lists into a matrix
+matFromLists :: (Vec j a v, Vec i v m) => [[a]] -> m
+matFromLists = fromList . (P.map fromList)
+{-# INLINE matFromLists #-}
+
+-- | convert a list into a matrix. (row-major order)
+matFromList :: forall i j v m a. (Vec i v m, Vec j a v, Nat i) => [a] -> m
+matFromList  = matFromLists . groupsOf (nat(undefined::i))
+  where groupsOf n xs = let (a,b) = splitAt n xs in a:(groupsOf n b)
+{-# INLINE matFromList  #-}
+
+
+
diff --git a/Data/Vec/Instances.hs b/Data/Vec/Instances.hs
new file mode 100644
--- /dev/null
+++ b/Data/Vec/Instances.hs
@@ -0,0 +1,117 @@
+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}
+
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE NoMonomorphismRestriction #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+module Data.Vec.Instances where
+
+import Data.Vec.Base as V
+import Data.Vec.Nat
+import Foreign.Storable
+import Foreign.Ptr
+
+-- Storable instances. 
+
+instance Storable a => Storable (a:.()) where
+  sizeOf _ = sizeOf (undefined::a)
+  alignment _ = alignment (undefined::a)
+  peek !p = peek (castPtr p) >>= \a -> return (a:.())
+  peekByteOff !p !o = peek (p`plusPtr`o)
+  peekElemOff !p !i = peek (p`plusPtr`(i*sizeOf(undefined::a)))
+  poke !p (a:._) = poke (castPtr p) a
+  pokeByteOff !p !o !x = poke (p`plusPtr`o) x
+  pokeElemOff !p !i !x = poke (p`plusPtr`(i*sizeOf(undefined::a))) x
+  {-# INLINE sizeOf #-}
+  {-# INLINE alignment #-}
+  {-# INLINE peek #-}
+  {-# INLINE peekByteOff #-}
+  {-# INLINE peekElemOff #-}
+  {-# INLINE poke #-}
+  {-# INLINE pokeByteOff #-}
+  {-# INLINE pokeElemOff #-}
+
+instance (Vec (Succ (Succ n)) a (a:.a:.v), Storable a, Storable (a:.v)) 
+  => Storable (a:.a:.v) 
+  where
+  sizeOf _ = sizeOf (undefined::a) + sizeOf (undefined::(a:.v))
+  alignment _ = alignment (undefined::a)
+  peek !p = 
+    peek (castPtr p) >>= \a -> 
+    peek (castPtr (p`plusPtr`sizeOf(undefined::a))) >>= \v -> 
+    return (a:.v)
+  peekByteOff !p !o = peek (p`plusPtr`o)
+  peekElemOff !p !i = peek (p`plusPtr`(i*sizeOf(undefined::(a:.a:.v))))
+  poke !p (a:.v) = 
+    poke (castPtr p) a >> 
+    poke (castPtr (p`plusPtr`sizeOf(undefined::a))) v
+  pokeByteOff !p !o !x = poke (p`plusPtr`o) x
+  pokeElemOff !p !i !x = poke (p`plusPtr`(i*sizeOf(undefined::(a:.a:.v)))) x
+  {-# INLINE sizeOf #-}
+  {-# INLINE alignment #-}
+  {-# INLINE peek #-}
+  {-# INLINE peekByteOff #-}
+  {-# INLINE peekElemOff #-}
+  {-# INLINE poke #-}
+  {-# INLINE pokeByteOff #-}
+  {-# INLINE pokeElemOff #-}
+
+
+-- Num and Fractional instances : All arithmetic is done component-wise and
+-- literals construct uniform vectors and matrices. 
+--
+-- The rule is simple : 
+--    If the method is unary, it's a map.  
+--    If it's binary, it's a zipWith.
+--
+-- You are free to ignore these instances if the definition of (*) offends you.
+
+instance
+    (Eq (a:.u)
+    ,Show (a:.u)
+    ,Num a
+    ,Map a a (a:.u) (a:.u) 
+    ,ZipWith a a a (a:.u) (a:.u) (a:.u)
+    ,Vec (Succ l) a (a:.u)
+    )
+    => Num (a:.u) 
+  where
+    (+) u v = V.zipWith (+) u v 
+    (-) u v = V.zipWith (-) u v
+    (*) u v = V.zipWith (*) u v
+    abs u = V.map abs u
+    signum u = V.map signum u
+    fromInteger i = vec (fromInteger i)
+    {-# INLINE (+) #-}
+    {-# INLINE (-) #-}
+    {-# INLINE (*) #-}
+    {-# INLINE abs #-}
+    {-# INLINE signum #-}
+    {-# INLINE fromInteger #-}
+
+
+instance 
+    (Fractional a
+    ,Ord (a:.u)
+    ,ZipWith a a a (a:.u) (a:.u) (a:.u)
+    ,Map a a (a:.u) (a:.u)
+    ,Vec (Succ l) a (a:.u)
+    ,Show (a:.u)
+    ) 
+    => Fractional (a:.u) 
+  where
+    (/) u v = V.zipWith (/) u v
+    recip u = V.map recip u
+    fromRational r = vec (fromRational r)
+    {-# INLINE (/) #-}
+    {-# INLINE recip #-}
+    {-# INLINE fromRational #-}
diff --git a/Data/Vec/LinAlg.hs b/Data/Vec/LinAlg.hs
new file mode 100644
--- /dev/null
+++ b/Data/Vec/LinAlg.hs
@@ -0,0 +1,726 @@
+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}
+
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE NoMonomorphismRestriction #-}
+{-# LANGUAGE PatternSignatures #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+{-# OPTIONS_HADDOCK ignore-exports,prune #-}
+
+module Data.Vec.LinAlg 
+  (dot
+  ,normSq
+  ,norm
+  ,normalize
+  ,cross
+  ,homPoint
+  ,homVec
+  ,project
+  ,multvm
+  ,multmv
+  ,multmm
+  ,translate
+  ,column
+  ,row
+  ,Transpose(transpose)
+  ,SetDiagonal(setDiagonal)
+  ,GetDiagonal(getDiagonal)
+  ,scale
+  ,diagonal
+  ,identity
+  ,Det(det)
+  ,cramer'sRule
+  ,NearZero(nearZero)
+  ,GaussElim(gaussElim)
+  ,BackSubstitute(backSubstitute)
+  ,BackSubstitute'(backSubstitute')
+  ,invert
+  ,invertAndDet
+  ,solve
+  ) where
+
+import Prelude hiding (map,zipWith,foldl,foldr,reverse,take,drop,
+                       head,tail,sum,length,last)
+import qualified Prelude as P
+import Data.Vec.Base
+import Data.Vec.Nat
+import Data.Vec.Instances
+
+import Control.Monad
+import Data.Maybe
+
+
+-- | dot / inner / scalar product
+dot ::  (Num a, Num v, Fold a v) => v -> v -> a
+dot u v = sum (u*v)
+{-# INLINE dot #-}
+
+-- | vector norm, squared
+normSq ::  (Num a, Num v, Fold a v) => v -> a
+normSq v = dot v v
+{-# INLINE normSq #-}
+
+-- | vector / L2 / Euclidean norm
+norm ::  (Num v, Floating a, Fold a v) => v -> a
+norm v = sqrt (dot v v)
+{-# INLINE norm #-}
+
+-- | @normalize v@ is a unit vector in the direction of @v@. @v@ is assumed
+-- non-null.
+normalize :: (Floating a, Num v, Fold a v, Map a a v v) => v -> v
+normalize v = map (/(norm v)) v
+{-# INLINE normalize #-}
+
+-- | 3d cross product.
+cross :: Num a => Vec3 a -> Vec3 a -> Vec3 a
+cross (ux:.uy:.uz:.()) (vx:.vy:.vz:.()) =
+  (uy*vz-uz*vy):.(uz*vx-ux*vz):.(ux*vy-uy*vx):.()
+{-# INLINE cross #-}
+
+-- | lift a point into homogenous coordinates
+homPoint ::  (Snoc v a v', Num a) => v -> v'
+homPoint v = snoc v 1
+{-# INLINE homPoint #-}
+
+-- | point-at-infinity in homogenous coordinates
+homVec ::  (Snoc v a v', Num a) => v -> v'
+homVec   v = snoc v 0
+{-# INLINE homVec   #-}
+
+-- | project a vector from homogenous coordinates. Last vector element is
+-- assumed non-zero.
+project :: 
+  ( Reverse' () t1 v'
+  , Fractional t1
+  , Vec a t t1
+  , Reverse' () v (t :. t1)
+  ) => v -> v'
+project  v = case reverse v of (w:.u) -> reverse (u/vec w)
+{-# INLINE project  #-}
+
+
+-- | row vector * matrix
+multvm :: 
+  ( Transpose m mt
+  , Map v a mt v'
+  , Fold a v
+  , Num a
+  , Num v
+  ) => v -> m -> v'
+multvm v m = map (dot v) (transpose m)
+{-# INLINE multvm #-}
+
+-- | matrix * column vector
+multmv :: 
+  ( Map v a m v'
+  , Num v
+  , Fold a v
+  , Num a
+  ) => m -> v -> v'
+multmv m v = map (dot v) m
+{-# INLINE multmv #-}
+
+-- | matrix * matrix 
+multmm :: 
+  (Map v v' m1 m3
+  ,Map v a b v'
+  ,Transpose m2 b
+  ,Fold a v
+  ,Num v
+  ,Num a
+  ) => m1 -> m2 -> m3
+multmm a b = map (\v -> map (dot v) (transpose b)) a
+{-# INLINE multmm #-}
+
+-- | apply a translation to a projective transformation matrix
+translate :: 
+  (Transpose m mt
+  ,Reverse' () mt (v' :. t)
+  ,Reverse' (v' :. ()) t v'1
+  ,Transpose v'1 m
+  ,Num v'
+  ,Num a
+  ,Snoc v a v'
+  ) => v -> m -> m
+translate v m = 
+  case reverse (transpose m) of
+    (h:.t) -> transpose (reverse (((homVec v) + h) :. t))
+{-# INLINE translate #-}
+
+-- | get the @n@-th column as a vector. @n@ is a type-level natural.
+column ::  (Transpose m mt, Access n v mt) => n -> m -> v
+column n = get n . transpose 
+{-# INLINE row #-}
+
+-- | get the @n@-th row as a vector. @n@ is a type-level natural.
+row ::  (Access n a v) => n -> v -> a
+row n = get n
+{-# INLINE column #-}
+
+
+-- Matrix transpose wrapper class: infers type of one argument from the other,
+-- because Transpose` can't do it, the fundeps there can't be bijective
+
+-- | matrix transposition
+class Transpose a b | a -> b, b -> a where 
+  transpose :: a -> b
+
+instance Transpose () () where
+  transpose = id
+
+instance 
+    (Vec (Succ n) s (s:.ra)  --(s:ra) is an n-vector of s'es (row of a)
+    ,Vec (Succ m) (s:.ra) ((s:.ra):.a)  --a is an m-vector of ra's
+    ,Vec (Succ m) s (s:.rb)  --rb is an m-vector of s'es (row of b)
+    ,Vec (Succ n) (s:.rb) ((s:.rb):.b)  --b is an n-vector of rb's
+    ,Transpose' ((s:.ra):.a) ((s:.rb):.b)
+    )
+    => Transpose ((s:.ra):.a) ((s:.rb):.b)
+  where
+    transpose = transpose'
+    {-# INLINE transpose #-}
+
+
+
+class Transpose' a b | a->b
+  where transpose' :: a -> b
+
+instance Transpose' () () where 
+  transpose' = id
+  {-# INLINE transpose' #-}
+
+instance 
+    (Transpose' vs vs') => Transpose' ( () :. vs ) vs'
+  where
+    transpose' (():.vs) = transpose' vs
+    {-# INLINE transpose' #-}
+
+instance Transpose' ((x:.()):.()) ((x:.()):.()) where
+  transpose' = id
+
+instance 
+    (Head xss_h xss_hh
+    ,Map xss_h xss_hh (xss_h:.xss_t) xs'
+    ,Tail xss_h xss_ht
+    ,Map xss_h xss_ht (xss_h:.xss_t) xss_
+    ,Transpose' (xs :. xss_) xss'
+    )
+    => Transpose' ((x:.xs):.(xss_h:.xss_t)) ((x:.xs'):.xss') 
+  where
+    transpose' ((x:.xs):.xss) =
+      (x :. (map head xss)) :. (transpose' (xs :. (map tail xss) :: (xs:.xss_)))
+    {-# INLINE transpose' #-}
+
+
+
+
+
+class SetDiagonal v m | m -> v, v -> m where
+  -- |set the diagonal of an n-by-n matrix to a given n-vector
+  setDiagonal :: v -> m -> m
+
+instance (Vec n a v, Vec n r m, SetDiagonal' N0 v m) => SetDiagonal v m where
+  setDiagonal v m = setDiagonal' (undefined::N0) v m
+  {-# INLINE setDiagonal #-}
+
+class SetDiagonal' n v m  where
+  setDiagonal' :: n -> v -> m -> m
+
+instance SetDiagonal' n () m where
+  setDiagonal' _ _ m = m
+  {-# INLINE setDiagonal' #-}
+
+instance 
+    ( SetDiagonal' (Succ n) v m
+    , Access n a r
+    ) => SetDiagonal' n (a:.v) (r:.m) 
+  where
+    setDiagonal' _ (a:.v) (r:.m) = 
+       (set (undefined::n) a r) :. (setDiagonal' (undefined::Succ n) v m)
+    {-# INLINE setDiagonal' #-}
+
+
+
+class GetDiagonal m v | m -> v, v -> m where
+  -- |get the diagonal of an n-by-n matrix as a vector
+  getDiagonal :: m -> v
+
+instance (Vec n a v, Vec n v m, GetDiagonal' N0 () m v) => GetDiagonal m v where
+  getDiagonal m = getDiagonal' (undefined::N0) () m
+  {-# INLINE getDiagonal #-}
+
+class GetDiagonal' n p m v where
+  getDiagonal' :: n -> p -> m -> v
+
+instance 
+    (Access n a r
+    ,Append p (a:.()) (a:.p)
+    ) => GetDiagonal' n p (r:.()) (a:.p) 
+  where
+    getDiagonal' _ p (r:.()) = append p ((get (undefined::n) r) :. ())
+    {-# INLINE getDiagonal' #-}
+
+instance 
+    (Access n a r
+    ,Append p (a:.()) p'
+    ,GetDiagonal' (Succ n) p' (r:.m) v
+    ) 
+    => GetDiagonal' n p (r:.r:.m) v
+  where
+    getDiagonal' _ p (r:.m) = 
+      getDiagonal' (undefined::Succ n) (append p ((get (undefined::n) r):.())) m
+    {-# INLINE getDiagonal' #-}
+
+
+-- | @scale v m@ multiplies the diagonal of matrix @m@ by the vector @s@, component-wise. So
+-- @scale 5 m@ multiplies the diagonal by 5, whereas @scale 2:.1 m@
+-- only scales the x component.
+scale :: 
+  ( GetDiagonal' N0 () m r
+  , Num r
+  , Vec n a r
+  , Vec n r m
+  , SetDiagonal' N0 r m
+  ) => r -> m -> m
+scale s m = setDiagonal (s * (getDiagonal m)) m
+{-# INLINE scale #-}
+
+
+-- | @diagonal v@ is a square matrix with the vector v as the diagonal, and 0
+-- elsewhere.
+diagonal :: (Vec n a v, Vec n v m, SetDiagonal v m, Num m) => v -> m
+diagonal v = setDiagonal v 0
+{-# INLINE diagonal #-}
+
+
+-- | identity matrix (square)
+identity :: (Vec n a v, Vec n v m, Num v, Num m, SetDiagonal v m) => m
+identity = diagonal 1 
+{-# INLINE identity #-}
+
+
+-- DropConsec: this is a helper function for computing determinants. Given an
+-- n-vector v, drop each element from v and collect the remaning (n-1)-vectors
+-- into an n-vector (ie an n-by-(n-1) matrix)
+class DropConsec v vv | v -> vv where
+  dropConsec :: v -> vv
+
+instance 
+  (Vec n a v
+  ,Pred n n_
+  ,Vec n_ a v_
+  ,Vec n v_ vv
+  ,DropConsec' () v vv
+  ) => DropConsec v vv
+  where
+    dropConsec v = dropConsec' () v 
+    {-# INLINE dropConsec #-}
+
+class DropConsec' p v vv  where
+  dropConsec' :: p -> v -> vv
+    
+instance DropConsec' p (a:.()) (p:.()) where
+  dropConsec' p (a:.()) = (p:.())
+  {-# INLINE dropConsec' #-}
+
+instance 
+    (Append p (a:.v) x
+    ,Append p (a:.()) y
+    ,DropConsec' y (a:.v) z
+    ) 
+    => DropConsec' p (a:.a:.v) (x:.z)
+  where
+    dropConsec' p (a:.v) = 
+      (append p v) :. (dropConsec' (append p (a:.())) v)
+    {-# INLINE dropConsec' #-}
+
+
+
+--Alternating: vector of alternating positive/negative values. This is also a
+--helper for computing determinants
+class Alternating n a v | v -> n a where
+  alternating :: n -> a -> v
+
+instance Alternating N1 a (a:.()) where
+  alternating _ !a = a:.()
+  {-# INLINE alternating #-}
+
+instance (Num a, Alternating n a (a:.v)) => Alternating (Succ n) a (a:.a:.v) where
+  alternating _ !a = a:.(alternating (undefined::n) (negate $! a))
+  {-# INLINE alternating #-}
+
+
+-- The Determinant of a square matrix, by minor expansion. 
+class Det' a m | m -> a where
+  det' :: m -> a
+
+instance Num a => Det' a ((a:.a:.()):.(a:.a:.()):.()) where
+  det' ( (a:.b:.()) :. (c:.d:.()) :. () ) = a*d-b*c
+  {-# INLINE det' #-}
+
+--this instance is particularly ugly in order to avoid overlapping with the one above
+instance
+    (Num a
+    ,Num (a:.a:.a:.v)
+    ,Fold a (a:.a:.a:.v)
+    ,Alternating (Succ (Succ (Succ n))) a (a:.a:.a:.v)
+    ,DropConsec (a:.a:.a:.v) vv
+    ,Map (a:.a:.a:.v) vv ((a:.a:.a:.v):.(a:.a:.a:.v):.m) vmt
+    ,Transpose vmt vm
+    ,Map ((a:.a:.v):.(a:.a:.v):.m_) a vm (a:.a:.a:.v)
+    ,Det' a ((a:.a:.v):.(a:.a:.v):.m_)
+    ,Vec (Succ (Succ (Succ n))) a (a:.a:.a:.v)
+    ,Vec (Succ (Succ (Succ n))) (a:.a:.a:.v) ((a:.a:.a:.v):.(a:.a:.a:.v):.(a:.a:.a:.v):.m)
+    )
+     => 
+    Det' a ((a:.a:.a:.v):.(a:.a:.a:.v):.(a:.a:.a:.v):.m)
+  where
+    det' (mh:.mt) =
+      sum ((alternating undefined 1) * mh *
+          (map det' (transpose (map dropConsec mt :: vmt))))
+    {-# INLINE det' #-}
+
+
+-- For now, use wrapper class to allow type inference. I think maybe the
+-- squareness of the matrix is keeping Det' from inferring properly, so we'll
+-- enforce that here. But really I have no clue.
+
+
+class Det n a m | m -> a where
+  -- | Determinant by minor expansion. Unfolds into a closed form expression.
+  -- This should be the fastest way for 4x4 and smaller, but @snd . gaussElim@
+  -- works too.
+  det :: m -> a
+
+instance (Vec n a r, Vec n r m, Det' a m) => Det n a m where
+  det = det'
+  {-# INLINE det #-}
+
+
+
+--ReplConsec : this is a helper for implementing Cramer's rule.  Given an
+--n-vector v and a value r, replace each consecutive element from v with r,
+--and collect the resulting n-vectors into an n-vector (ie an n-by-n matrix)
+
+class ReplConsec a v vv | v->a, v->vv, vv->v, vv->a where
+  replConsec :: a -> v -> vv
+
+instance 
+  (Vec n a v
+  ,Vec n v vv
+  ,ReplConsec' a () v vv
+  ) => ReplConsec a v vv
+  where
+    replConsec a v = replConsec' a () v :: vv
+    {-# INLINE replConsec #-}
+
+class ReplConsec' a p v vv where
+  replConsec' :: a -> p -> v -> vv
+
+instance ReplConsec' a p () () where
+  replConsec' _ _ () = ()
+  {-# INLINE replConsec' #-}
+
+instance 
+    (Append p (a:.v) x
+    ,Append p (a:.()) y
+    ,ReplConsec' a y v z
+    ) 
+    => ReplConsec' a p (a:.v) (x:.z)
+  where
+    replConsec' r p (a:.v) = 
+      (append p (r:.v)) :. (replConsec' r (append p (a :. ())) v)
+    {-# INLINE replConsec' #-}
+
+
+
+
+-- | @cramer'sRule m v@ computes the solution to @m\`multmv\`x=v@  using the
+-- eponymous method. For larger than 3x3 you will want to use 'solve', which
+-- uses 'gaussElim'. Cramer's rule, however, unfolds into a closed-form
+-- expression, with no branches or allocations (other than the result). You may
+-- need to increase the unfolding threshold to see this.
+
+cramer'sRule :: 
+  (Map a a1 b1 v
+  ,Transpose w b1
+  ,ZipWith a2 b vv v m w
+  ,ReplConsec' a2 () b vv
+  ,Vec n b vv
+  ,Vec n a2 b
+  ,Fractional a1
+  ,Det' a1 m
+  ,Det' a1 a
+  ) => m -> v -> v
+cramer'sRule m b =
+  case map (\m' -> (det' m')/(det' m)) 
+           (transpose (zipWith replConsec b m)) 
+    of b' -> b' `asTypeOf` b 
+{-# INLINE cramer'sRule #-}
+
+
+
+
+
+
+mapFst f (a,b) = (f a,b)
+{-# INLINE mapFst #-}
+
+
+class Num a => NearZero a where
+  -- | @nearZero x@ should be true when x is close enough to 0 to cause
+  -- significant error in division. 
+  nearZero :: a -> Bool
+  nearZero 0 = True
+  nearZero _ = False
+  {-# INLINE nearZero #-}
+
+instance NearZero Float where
+  nearZero x = abs x < 1e-6
+  {-# INLINE nearZero #-}
+
+instance NearZero Double where
+  nearZero x = abs x < 1e-14
+  {-# INLINE nearZero #-}
+
+instance NearZero Rational
+
+
+
+
+-- Pivot1 : find a non-zero pivot column and put a 1 there. Second return
+-- argument tracks value of determinant. Returns nothing if no pivot in the
+-- first row. Does not try to find the 'best' pivot, only an acceptable one:
+-- matrices are assumed small, roundoff error should be negligible. 
+
+class Pivot1 a m | m -> a where
+  pivot1 :: m -> Maybe (m,a)
+
+instance Pivot1 a () where
+  pivot1 _ = Nothing
+
+instance 
+    ( Fractional a, NearZero a
+    ) => Pivot1 a ((a:.()):.()) 
+  where
+    pivot1 ((p:._):._) 
+      | nearZero p = Nothing
+      | otherwise  = Just (1,p)
+    {-# INLINE pivot1 #-}
+
+instance 
+    ( Fractional a, NearZero a 
+    , Map a a (a:.r) (a:.r)
+    ) => Pivot1 a ((a:.(a:.r)):.()) 
+  where
+    pivot1 ((p:.r):._) 
+      | nearZero p = Nothing
+      | otherwise  = Just ((1 :. (map (/p) r)):.(), p)
+    {-# INLINE pivot1 #-}
+
+instance 
+    ( Fractional a, NearZero a
+    , Map a a (a:.r) (a:.r)
+    , ZipWith a a a (a:.r) (a:.r) (a:.r) 
+    , Map (a:.r) (a:.r) ((a:.r):.rs) ((a:.r):.rs)
+    , Pivot1 a ((a:.r):.rs) 
+    ) => Pivot1 a ((a:.r):.(a:.r):.rs) 
+  where
+    pivot1 (row@(p:._):.rows) 
+      | nearZero p = pivot1 rows >>= \(r:.rs,p)-> Just(r:.row:.rs,p)
+      | otherwise  = Just ( first:.(map add rows) , p)
+          where first        = map (/p) row
+                add r@(x:._) = zipWith (-) r . map (*x) $ first 
+    {-# INLINE pivot1 #-}
+
+
+-- Pivot : find a pivot. Second return argument tracks determinant.
+-- Returns Nothing if no pivot anywhere.
+
+class Pivot a m | m -> a where
+  pivot :: m -> Maybe (m,a)
+
+instance Pivot a (():.v) where
+  pivot _ = Nothing
+  {-# INLINE pivot #-}
+
+instance 
+    ( Fractional a
+    , NearZero a
+    , Pivot1 a rs 
+    , Tail (a:.r) r
+    , Map (a:.r) r ((a:.r):.rs) (r:.rs') 
+    , Map r (a:.r) (r:.rs') ((a:.r):.rs)
+    , Pivot1 a ((a:.r):.rs)
+    , Pivot a (r:.rs')
+    ) => Pivot a ((a:.r):.rs) 
+  where
+    pivot m = 
+      mplus (pivot1 m) 
+            (pivot (map tail m) >>= return . mapFst (map (0:.)) )
+    {-# INLINE pivot #-}
+
+
+
+-- | Gaussian elimination, adapted from Mirko Rahn:
+-- <http://www.haskell.org/pipermail/glasgow-haskell-users/2007-May/012648.html>
+--
+-- This is more of a proof of concept. Using a foreign C function will run
+-- slightly faster, and compile much faster. But where is the fun in that?
+-- Set your unfolding threshold as high as possible.
+
+class GaussElim a m | m -> a where
+  -- | @gaussElim m@ returns a pair @(m',d)@ where @m'@ is @m@ in row echelon
+  -- form and @d@ is the determinant of @m@. The determinant of @m'@ is 1 or 0,
+  -- i.e., the leading coefficient of each non-zero row is 1.  
+   
+  gaussElim :: m -> (m,a)
+
+instance (Num a, Pivot a (r:.())) => GaussElim a (r:.())
+  where
+    gaussElim m = fromMaybe (m,1) (pivot m) 
+    {-# INLINE gaussElim #-}
+
+instance 
+    ( Fractional a
+    , Map (a:.r) r ((a:.r):.rs) rs_
+    , Map r (a:.r) rs_ ((a:.r):.rs) 
+    , Pivot a ((a:.r):.(a:.r):.rs)
+    , GaussElim a rs_
+    ) => GaussElim a ((a:.r):.(a:.r):.rs)
+  where
+    gaussElim m =
+      flip (maybe (m,1)) (pivot m) $ \(row:.rows,p) ->
+        case gaussElim (map tail rows)
+          of (rows',p') -> ( row:.(map (0:.) rows') , p*p')
+    {-# INLINE gaussElim #-}
+
+
+
+class BackSubstitute m where
+  -- | backSubstitute takes a full rank matrix from row echelon form to reduced
+  -- row echelon form. Returns @Nothing@ if the matrix is rank deficient. 
+  backSubstitute :: m -> Maybe m 
+
+instance BackSubstitute ((a:.r):.()) where
+  backSubstitute = Just . id
+  {-# INLINE backSubstitute #-}
+
+instance 
+    ( Map (a:.r) r ((a:.r):.rs) rs_ --map tail
+    , Map r (a:.r) rs_ ((a:.r):.rs) --map cons
+    , Fold (a,a:.r) aas
+    , ZipWith a a a (a:.r) (a:.r) (a:.r)
+    , Map a a (a:.r) (a:.r)
+    , ZipWith a (a:.r) (a,a:.r) r ((a:.r):.rs) aas
+    , Num a, NearZero a
+    , BackSubstitute rs_
+    ) => BackSubstitute ((a:.r):.(a:.r):.rs)
+  where
+    backSubstitute (r@(rh:.rt):.rs) 
+      | nearZero (1-rh) = 
+        liftM (map (0:.)) (backSubstitute . map tail $ rs) >>= \rs' -> 
+          return . (:.rs') . foldl (\v (a,w) -> sub v a w) r $ 
+            zipWith (,) rt rs'
+      | otherwise = Nothing -- rank deficient
+          where sub v a = zipWith (-) v . map (*a)
+    {-# INLINE backSubstitute #-}
+
+
+
+
+class BackSubstitute' m where
+  -- | backSubstitute' takes a full rank matrix from row echelon form to reduced
+  -- row echelon form. Returns garbage is matrix is rank deficient.
+  backSubstitute' :: m -> m 
+
+instance BackSubstitute' ((a:.r):.()) where
+  backSubstitute' = id
+  {-# INLINE backSubstitute' #-}
+
+instance 
+    ( Map (a:.r) r ((a:.r):.rs) rs_ --map tail
+    , Map r (a:.r) rs_ ((a:.r):.rs) --map cons
+    , Fold (a,a:.r) aas
+    , ZipWith a a a (a:.r) (a:.r) (a:.r)
+    , Map a a (a:.r) (a:.r)
+    , ZipWith a (a:.r) (a,a:.r) r ((a:.r):.rs) aas
+    , Num a
+    , BackSubstitute' rs_
+    ) => BackSubstitute' ((a:.r):.(a:.r):.rs)
+  where
+    backSubstitute' (r@(_:.rt):.rs) = 
+      case map (0:.) (backSubstitute' . map tail $ rs) 
+        of rs' -> (:.rs') $ foldl (\ v (a,w) -> sub v a w) r 
+                              (zipWith (,) rt rs')
+      where sub v a = zipWith (-) v . map (*a)
+    {-# INLINE backSubstitute' #-}
+
+
+-- | @invert m@ returns @Just@ the inverse of @m@ or @Nothing@ if @m@ is singular.
+invert :: forall n a r m r' m'. 
+  ( Num r, Num m
+  , Vec n a r     -- r is row type
+  , Vec n r m     -- m is matrix type
+  , Append r r r' -- r' is a row of augmented matrix
+  , ZipWith r r r' m m m' -- m' is the augmented matrix
+  , Drop n r' r -- get the right half of an augmented matrix row
+  , Map r' r m' m -- get the right half of the augmented matrix
+  , SetDiagonal r m -- needed to make identity matrix
+  , GaussElim a m'
+  , BackSubstitute m'
+  ) => m -> Maybe m
+invert m = 
+  return i >>= backSubstitute . fst . gaussElim . zipWith append m 
+           >>= return . map dropn
+  where dropn = drop (undefined::n)
+        i = identity :: m
+{-# INLINE invert #-}
+
+-- | inverse and determinant. If det = 0, inverted matrix is garbage.
+invertAndDet :: forall n a r m r' m'. 
+  ( Num r, Num m
+  , Vec n a r     -- r is row type
+  , Vec n r m     -- m is matrix type
+  , Append r r r' -- r' is a row of augmented matrix
+  , ZipWith r r r' m m m' -- m' is the augmented matrix
+  , Drop n r' r -- get the right half of an augmented matrix row
+  , Map r' r m' m -- get the right half of the augmented matrix
+  , SetDiagonal r m -- needed to make identity matrix
+  , GaussElim a m'
+  , BackSubstitute' m'
+  ) => m -> (m,a)
+invertAndDet m = 
+  mapFst ( (map dropn) . backSubstitute') . gaussElim . zipWith append m $ i
+  where dropn = drop (undefined::n)
+        i = identity :: m
+{-# INLINE invertAndDet #-}
+
+
+-- | Solution of linear system by Gaussian elimination. Returns @Nothing@
+-- if no solution. 
+solve :: forall n a v r m r' m'. 
+  ( Num r, Num m
+  , Vec n a r     -- r is row type
+  , Vec n r m     -- m is matrix type
+  , Snoc r a r'   -- a row of the extended matrix is one longer
+  , ZipWith r a r' m r m' -- m' is the augmented matrix
+  , Drop n r' (a:.()) -- get the right part of an augmented matrix row
+  , Map r' a m' r -- get the right part of the augmented matrix
+  , GaussElim a m'
+  , BackSubstitute m'
+  ) => m -> r -> Maybe r
+solve m v = 
+  return v >>= backSubstitute . fst . gaussElim . zipWith snoc m 
+           >>= return . map (head . drop (undefined::n)) 
+{-# INLINE solve #-}
+
diff --git a/Data/Vec/Nat.hs b/Data/Vec/Nat.hs
new file mode 100644
--- /dev/null
+++ b/Data/Vec/Nat.hs
@@ -0,0 +1,69 @@
+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}
+
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE NoMonomorphismRestriction #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-- | Type level naturals. @Ni@ is a type, @ni@ an undefined value of that type,
+-- for @i <- [0..19]@
+module Data.Vec.Nat where
+
+
+data N0
+data Succ a
+
+type N1  = Succ N0
+type N2  = Succ N1
+type N3  = Succ N2
+type N4  = Succ N3
+type N5  = Succ N4
+type N6  = Succ N5
+type N7  = Succ N6
+type N8  = Succ N7
+type N9  = Succ N8
+type N10 = Succ N9
+type N11 = Succ N10
+type N12 = Succ N11
+type N13 = Succ N12
+type N14 = Succ N13
+type N15 = Succ N14
+type N16 = Succ N15
+type N17 = Succ N16
+type N18 = Succ N17
+type N19 = Succ N18
+
+n0  :: N0  ; n0  = undefined
+n1  :: N1  ; n1  = undefined
+n2  :: N2  ; n2  = undefined
+n3  :: N3  ; n3  = undefined
+n4  :: N4  ; n4  = undefined
+n5  :: N5  ; n5  = undefined
+n6  :: N6  ; n6  = undefined
+n7  :: N7  ; n7  = undefined
+n8  :: N8  ; n8  = undefined
+n9  :: N9  ; n9  = undefined
+n10 :: N10 ; n10  = undefined
+n11 :: N11 ; n11  = undefined
+n12 :: N12 ; n12  = undefined
+n13 :: N13 ; n13  = undefined
+n14 :: N14 ; n14  = undefined
+n15 :: N15 ; n15  = undefined
+n16 :: N16 ; n16  = undefined
+n17 :: N17 ; n17  = undefined
+n18 :: N18 ; n18  = undefined
+n19 :: N19 ; n19  = undefined
+
+-- | @nat n@ yields the @Int@ value of the type-level natural @n@.
+class Nat n where nat :: n -> Int
+instance Nat N0 where nat _ = 0
+instance Nat a => Nat (Succ a) where nat _ = 1+(nat (undefined::a))
+
+class Pred x y | x -> y, y -> x
+instance Pred (Succ N0) N0
+instance Pred (Succ n) p => Pred (Succ (Succ n)) (Succ p)
+
diff --git a/Data/Vec/Packed.hs b/Data/Vec/Packed.hs
new file mode 100644
--- /dev/null
+++ b/Data/Vec/Packed.hs
@@ -0,0 +1,152 @@
+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}
+
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE NoMonomorphismRestriction #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-- | Packed vectors : use these whenever possible. The regular vector type is
+-- just a gussied up linked list, but when vector functions are applied to
+-- these types, bracketed by @'pack'@ and @'unpack'@, then things unfold into
+-- perfectly optimized code.
+
+module Data.Vec.Packed where
+
+import Data.Vec.Base as V
+
+-- * Packed Vector Types
+data Vec2I = Vec2I {-#UNPACK#-} !Int 
+                   {-#UNPACK#-} !Int 
+
+data Vec3I = Vec3I {-#UNPACK#-} !Int 
+                   {-#UNPACK#-} !Int 
+                   {-#UNPACK#-} !Int
+
+data Vec4I = Vec4I {-#UNPACK#-} !Int 
+                   {-#UNPACK#-} !Int 
+                   {-#UNPACK#-} !Int
+                   {-#UNPACK#-} !Int
+
+data Vec2F = Vec2F {-#UNPACK#-} !Float 
+                   {-#UNPACK#-} !Float 
+
+data Vec3F = Vec3F {-#UNPACK#-} !Float 
+                   {-#UNPACK#-} !Float 
+                   {-#UNPACK#-} !Float
+
+data Vec4F = Vec4F {-#UNPACK#-} !Float 
+                   {-#UNPACK#-} !Float 
+                   {-#UNPACK#-} !Float
+                   {-#UNPACK#-} !Float
+
+data Vec2D = Vec2D {-#UNPACK#-} !Double 
+                   {-#UNPACK#-} !Double 
+
+data Vec3D = Vec3D {-#UNPACK#-} !Double 
+                   {-#UNPACK#-} !Double 
+                   {-#UNPACK#-} !Double
+
+data Vec4D = Vec4D {-#UNPACK#-} !Double 
+                   {-#UNPACK#-} !Double 
+                   {-#UNPACK#-} !Double
+                   {-#UNPACK#-} !Double
+
+-- * Packed Matrix Types. 
+type Mat22I = Vec2 Vec2I 
+type Mat23I = Vec2 Vec3I 
+type Mat33I = Vec3 Vec3I 
+type Mat34I = Vec3 Vec4I 
+type Mat44I = Vec4 Vec3I 
+
+type Mat22F = Vec2 Vec2F 
+type Mat23F = Vec2 Vec3F 
+type Mat33F = Vec3 Vec3F 
+type Mat34F = Vec3 Vec4F 
+type Mat44F = Vec4 Vec3F 
+
+type Mat22D = Vec2 Vec2D 
+type Mat23D = Vec2 Vec3D 
+type Mat33D = Vec3 Vec3D 
+type Mat34D = Vec3 Vec4D 
+type Mat44D = Vec4 Vec4D 
+
+
+-- | pack a matrix
+packMat ::  (Map v pv m pm, PackedVec pv v) => m -> pm
+packMat = V.map pack 
+{-# INLINE packMat #-}
+
+-- | unpack a matrix
+unpackMat ::  (Map pv v pm m, PackedVec pv v) => pm -> m
+unpackMat = V.map unpack
+{-# INLINE unpackMat #-}
+
+-- | PackedVec class : relates a packed vector type to its unpacked type For
+-- now, the fundep is not bijective -- It may be advantageous to have multiple
+-- packed representations for a canonical vector type. This may change. In the
+-- meantime, you may have to annotate return types.
+class PackedVec pv v | pv -> v  where
+  pack   :: v -> pv
+  unpack :: pv -> v
+
+instance PackedVec Vec2I (Vec2 Int) where
+  pack (x:.y:.()) = Vec2I x y 
+  unpack (Vec2I x y) = x:.y:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+instance PackedVec Vec3I (Vec3 Int) where
+  pack (x:.y:.z:.()) = Vec3I x y z
+  unpack (Vec3I x y z) = x:.y:.z:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+instance PackedVec Vec4I (Vec4 Int) where
+  pack (x:.y:.z:.w:.()) = Vec4I x y z w
+  unpack (Vec4I x y z w) = x:.y:.z:.w:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+
+instance PackedVec Vec2F (Vec2 Float) where
+  pack (x:.y:.()) = Vec2F x y 
+  unpack (Vec2F x y) = x:.y:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+instance PackedVec Vec3F (Vec3 Float) where
+  pack (x:.y:.z:.()) = Vec3F x y z
+  unpack (Vec3F x y z) = x:.y:.z:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+instance PackedVec Vec4F (Vec4 Float) where
+  pack (x:.y:.z:.w:.()) = Vec4F x y z w
+  unpack (Vec4F x y z w) = x:.y:.z:.w:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+
+instance PackedVec Vec2D (Vec2 Double) where
+  pack (x:.y:.()) = Vec2D x y 
+  unpack (Vec2D x y) = x:.y:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+instance PackedVec Vec3D (Vec3 Double) where
+  pack (x:.y:.z:.()) = Vec3D x y z
+  unpack (Vec3D x y z) = x:.y:.z:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+instance PackedVec Vec4D (Vec4 Double) where
+  pack (x:.y:.z:.w:.()) = Vec4D x y z w
+  unpack (Vec4D x y z w) = x:.y:.z:.w:.()
+  {-# INLINE pack #-}
+  {-# INLINE unpack #-}
+
+
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,22 @@
+Copyright (c) 2008, Scott E. Dillard
+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.
+
+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,4 @@
+#! /usr/bin/env runhaskell
+
+import Distribution.Simple
+main = defaultMain
diff --git a/Vec.cabal b/Vec.cabal
new file mode 100644
--- /dev/null
+++ b/Vec.cabal
@@ -0,0 +1,42 @@
+Name:                Vec
+Version:             0.9.0
+License:             BSD3
+License-file:        LICENSE
+Author:              Scott E. Dillard
+Maintainer:          Scott E. Dillard <sedillard@gmail.com>
+Stability:           Experimental
+Synopsis:            Fixed-length lists and low-dimensional linear algebra.
+Description:         
+   Vectors are represented by lists with type-encoded lengths. The constructor
+   is @:.@, which acts like a cons both at the value and type levels, with @()@
+   taking the place of nil. So @x:.y:.z:.()@ is a 3d vector. The library
+   provides a set of common list-like functions (map, fold, etc) for working
+   with vectors. Built up from these functions are a small but useful set of
+   linear algebra operations: matrix multiplication, determinants, solving
+   linear systems, inverting matrices.
+Cabal-version:       >=1.2
+Build-type:          Simple
+Category:            Math,Graphics
+
+library
+    Build-Depends:      base
+
+    Exposed-modules:    Data.Vec 
+                        Data.Vec.Base,
+                        Data.Vec.LinAlg,
+                        Data.Vec.Nat,
+                        Data.Vec.Instances
+                        Data.Vec.Packed
+    Extensions: 
+                        BangPatterns,
+                        EmptyDataDecls,
+                        ExistentialQuantification,
+                        FlexibleInstances, 
+                        FlexibleContexts,
+                        FunctionalDependencies,
+                        MultiParamTypeClasses, 
+                        NoMonomorphismRestriction,
+                        ScopedTypeVariables,
+                        TypeOperators, 
+                        TypeSynonymInstances,
+                        UndecidableInstances