diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,59 @@
+Copyright (c) 2014, Charles Durham
+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 names of the copyright holders nor the names of the 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.
+
+
+Copyright (c) 2008-2010, Balazs Komuves
+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 names of the copyright holders nor the names of the 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/src/Data/Vect/Floating.hs b/src/Data/Vect/Floating.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating.hs
@@ -0,0 +1,13 @@
+module Data.Vect.Floating
+       ( module Data.Vect.Floating.Base
+       , module Data.Vect.Floating.Interpolate
+       , module Data.Vect.Floating.Util.Dim2
+       , module Data.Vect.Floating.Util.Dim3
+       , module Data.Vect.Floating.Util.Projective
+       ) where
+
+import Data.Vect.Floating.Base
+import Data.Vect.Floating.Interpolate
+import Data.Vect.Floating.Util.Dim2
+import Data.Vect.Floating.Util.Dim3
+import Data.Vect.Floating.Util.Projective
diff --git a/src/Data/Vect/Floating/Base.hs b/src/Data/Vect/Floating/Base.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/Base.hs
@@ -0,0 +1,1005 @@
+{-# LANGUAGE DeriveDataTypeable, MultiParamTypeClasses, FunctionalDependencies, GeneralizedNewtypeDeriving, StandaloneDeriving, GADTs, FlexibleInstances, FlexibleContexts, ScopedTypeVariables #-}
+
+module Data.Vect.Floating.Base
+  ( AbelianGroup(..) , vecSum
+  , MultSemiGroup(..) , Ring , semigroupProduct
+  , LeftModule(..) , RightModule(..)
+  , Vector(..) , DotProd(..) , CrossProd(..)
+  , normalize , distance , angle , angle'
+  , UnitVector(..)
+  , Pointwise(..)
+  , Extend(..) , HasCoordinates(..) , Dimension(..)
+  , Matrix(..) , Tensor(..) , Diagonal (..) , Determinant(..)
+  , Orthogonal(..) , Projective(..) , MatrixNorms(..)
+  , Vec2(..) , Vec3(..) , Vec4(..)
+  , Mat2(..) , Mat3(..) , Mat4(..)
+  , Ortho2 , Ortho3 , Ortho4
+  , Normal2 , Normal3 , Normal4
+  , Proj3 , Proj4
+  , mkVec2 , mkVec3 , mkVec4
+  , project , project' , projectUnsafe , flipNormal
+  , householder, householderOrtho
+  )
+  where
+
+import Data.Typeable
+import Control.Monad
+import System.Random  
+import Foreign
+
+--------------------------------------------------------------------------------
+-- class declarations
+
+class AbelianGroup g where
+  (&+) :: g -> g -> g
+  (&-) :: g -> g -> g
+  neg  :: g -> g
+  zero :: g
+
+infixl 6 &+
+infixl 6 &- 
+
+vecSum :: AbelianGroup g => [g] -> g
+vecSum l = foldl (&+) zero l 
+
+class MultSemiGroup r where
+  (.*.) :: r -> r -> r
+  one   :: r
+
+class (AbelianGroup r, MultSemiGroup r) => Ring r 
+
+infixl 7 .*. 
+
+-- was: ringProduct :: Ring r => [r] -> r
+semigroupProduct :: MultSemiGroup r => [r] -> r 
+semigroupProduct l = foldl (.*.) one l
+
+class LeftModule r m where
+  lmul :: r -> m -> m
+  (*.) :: r -> m -> m
+  (*.) = lmul
+
+class RightModule m r | m -> r, r -> m where
+  rmul :: m -> r -> m
+  (.*) :: m -> r -> m
+  (.*) = rmul
+
+-- I'm not really sure about this.. may actually degrade the performance in some cases?  
+{- RULES
+"matrix multiplication left"   forall m n x.  (n .*. m) *. x = n *. (m *. x)  
+"matrix multiplication right"  forall m n x.  x .* (m .*. n) = (x .* m) .* n
+  -}
+
+infixr 7 *.
+infixl 7 .*
+
+class AbelianGroup (v a) => Vector a v where
+  mapVec    :: (a -> a) -> v a -> v a
+  scalarMul :: a -> v a -> v a
+  (*&) ::      a -> v a -> v a
+  (&*) ::      v a -> a -> v a 
+  (*&) s v = scalarMul s v
+  (&*) v s = scalarMul s v
+
+infixr 7 *&
+infixl 7 &*
+
+{-# RULES
+"scalar multiplication left"   forall (s :: Num s => s) (t :: Num t => t) x. t *& (s *& x) = (t*s) *& x 
+"scalar multiplication right"  forall (s :: Num s => s) (t :: Num t => t) x.  (x &* s) &* t = x &* (s*t)  
+  #-}
+
+class Floating a => DotProd a v where
+  (&.) :: v a -> v a -> a
+  norm    :: v a -> a
+  normsqr :: v a -> a
+  len     :: v a -> a
+  lensqr  :: v a -> a
+  len = norm
+  lensqr = normsqr
+  dotprod :: v a -> v a -> a
+  normsqr v = (v &. v)
+  norm = sqrt . lensqr
+  dotprod = (&.)
+
+infix 7 &.
+
+{-# RULES
+"len/square 1"   forall x.  (len x)*(len x) = lensqr x
+"len/square 2"   forall x.  (len x)^2 = lensqr x
+"norm/square 1"  forall x.  (norm x)*(norm x) = normsqr x
+"norm/square 2"  forall x.  (norm x)^2 = normsqr x
+  #-}
+
+
+normalize :: (Vector a v, DotProd a v) => v a -> v a
+normalize v = scalarMul (recip (len v)) v
+
+distance :: (Vector a v, DotProd a v) => v a -> v a -> a
+distance x y = norm (x &- y)
+
+-- | the angle between two vectors
+angle :: (Vector a v, DotProd a v) => v a -> v a -> a 
+angle x y = acos $ (x &. y) / (norm x * norm y)
+
+-- | the angle between two unit vectors
+angle' {- ' CPP is sensitive to primes -} :: (Vector a v, UnitVector a v u, DotProd a v) => u a -> u a -> a
+angle' x y = acos (fromNormal x &. fromNormal y)
+
+{-# RULES
+"normalize is idempotent"  forall x. normalize (normalize x) = normalize x
+  #-}
+
+class (Vector a v, DotProd a v) => UnitVector a v u | u -> v, v -> u where
+  mkNormal         :: v a -> u a       -- ^ normalizes the input
+  toNormalUnsafe   :: v a -> u a       -- ^ does not normalize the input!
+  fromNormal       :: u a -> v a
+  fromNormalRadius :: a -> u a -> v a
+  fromNormalRadius t n = t *& fromNormal n 
+
+-- | Projects the first vector down to the hyperplane orthogonal to the second (unit) vector
+project' :: (Vector a v, UnitVector a v u, DotProd a v) => v a -> u a -> v a
+project' what dir = projectUnsafe what (fromNormal dir)
+
+-- | Direction (second argument) is assumed to be a /unit/ vector!
+projectUnsafe :: (Vector a v, DotProd a v) => v a -> v a -> v a
+projectUnsafe what dir = what &- dir &* (what &. dir)
+
+project :: (Vector a v, DotProd a v) => v a -> v a -> v a
+project what dir = what &- dir &* ((what &. dir) / (dir &. dir))
+
+-- | Since unit vectors are not a group, we need a separate function.
+flipNormal :: UnitVector a v n => n a -> n a
+flipNormal = toNormalUnsafe . neg . fromNormal 
+
+-- | Cross product
+class CrossProd v where
+  crossprod :: v -> v -> v
+  (&^)      :: v -> v -> v
+  (&^) = crossprod
+ 
+-- | Pointwise multiplication 
+class Pointwise v where
+  pointwise :: v -> v -> v
+  (&!)      :: v -> v -> v
+  (&!) = pointwise 
+
+infix 7 &^
+infix 7 &!
+
+class HasCoordinates v x | v->x where
+  _1 :: v -> x
+  _2 :: v -> x
+  _3 :: v -> x
+  _4 :: v -> x
+
+-- | conversion between vectors (and matrices) of different dimensions
+class Extend a u v where
+  extendZero :: u a -> v a          -- ^ example: @extendZero (Vec2 5 6) = Vec4 5 6 0 0@
+  extendWith :: a -> u a -> v a   -- ^ example: @extendWith 1 (Vec2 5 6) = Vec4 5 6 1 1@
+  trim :: v a -> u a                -- ^ example: @trim (Vec4 5 6 7 8) = Vec2 5 6@
+
+-- | makes a diagonal matrix from a vector
+class Diagonal s t | t->s where
+  diag :: s -> t
+
+class Matrix m where
+  transpose :: m -> m 
+  inverse :: m -> m
+  idmtx :: m
+
+{-# RULES
+"transpose is an involution"  forall m. transpose (transpose m) = m
+"inverse is an involution"    forall m. inverse (inverse m) = m
+  #-}
+  
+class Matrix (m a) => Orthogonal a m o | m -> o, o -> m where
+  fromOrtho     :: o a -> m a
+  toOrthoUnsafe :: m a -> o a
+  
+class (AbelianGroup m, Matrix m) => MatrixNorms a m where
+  frobeniusNorm  :: m -> a       -- ^ the frobenius norm (= euclidean norm in the space of matrices)
+  matrixDistance :: m -> m -> a  -- ^ euclidean distance in the space of matrices
+  operatorNorm   :: m -> a      -- ^ (euclidean) operator norm (not implemented yet)
+  matrixDistance m n = frobeniusNorm (n &- m)
+  operatorNorm = error "operatorNorm: not implemented yet"
+  
+-- | Outer product (could be unified with Diagonal?)
+class Tensor t v | t->v where
+  outer :: v -> v -> t
+    
+class Determinant a m where
+  det :: m -> a
+
+class Dimension a where
+  dim :: a -> Int
+     
+-- | Householder matrix, see <http://en.wikipedia.org/wiki/Householder_transformation>.  
+-- In plain words, it is the reflection to the hyperplane orthogonal to the input vector.
+householder :: (Vector a v, UnitVector a v u, Matrix (m a), Vector a m, Tensor (m a) (v a)) => u a -> m a
+householder u = idmtx &- (2 *& outer v v) 
+  where v = fromNormal u
+
+householderOrtho :: (Vector a v, UnitVector a v u, Matrix (m a), Vector a m, Tensor (m a) (v a), Orthogonal a m o) => u a -> o a
+householderOrtho = toOrthoUnsafe . householder
+
+-- | \"Projective\" matrices have the following form: the top left corner
+-- is an any matrix, the bottom right corner is 1, and the top-right
+-- column is zero. These describe the affine orthogonal transformation of
+-- the space one dimension less.
+class (Vector a v, Orthogonal a n o, Diagonal (v a) (n a)) => Projective a v n o m p | m -> p, p -> m, p -> o, o -> p, p -> n, n -> p, p -> v, v -> p, n -> o, n -> v, v -> n where
+  fromProjective     :: p a -> m a
+  toProjectiveUnsafe :: m a -> p a
+  orthogonal         :: o a -> p a
+  linear             :: n a -> p a
+  translation        :: v a -> p a
+  scaling            :: v a -> p a
+
+--------------------------------------------------------------------------------
+-- Vec / Mat datatypes
+  
+data Vec2 a = Vec2 !a !a 
+  deriving (Read,Show,Typeable)
+data Vec3 a = Vec3 !a !a !a
+  deriving (Read,Show,Typeable)
+data Vec4 a = Vec4 !a !a !a !a
+  deriving (Read,Show,Typeable)
+
+-- | The components are /row/ vectors 
+
+data Mat2 a = Mat2 !(Vec2 a) !(Vec2 a) deriving (Read,Show)
+data Mat3 a = Mat3 !(Vec3 a) !(Vec3 a) !(Vec3 a) deriving (Read,Show)
+data Mat4 a = Mat4 !(Vec4 a) !(Vec4 a) !(Vec4 a) !(Vec4 a) deriving (Read,Show)
+
+-- | The assumption when dealing with these is always that they are of unit length.
+-- Also, interpolation works differently.
+newtype Normal2 a = Normal2 (Vec2 a) deriving (Read,Show,Storable,Dimension,Typeable) 
+newtype Normal3 a = Normal3 (Vec3 a) deriving (Read,Show,Storable,Dimension,Typeable)
+newtype Normal4 a = Normal4 (Vec4 a) deriving (Read,Show,Storable,Dimension,Typeable)
+
+deriving instance Floating a => DotProd a Normal2
+deriving instance Floating a => DotProd a Normal3
+deriving instance Floating a => DotProd a Normal4
+
+mkVec2 :: (a,a) -> Vec2 a
+mkVec3 :: (a,a,a) -> Vec3 a
+mkVec4 :: (a,a,a,a) -> Vec4 a
+
+mkVec2 (x,y)     = Vec2 x y
+mkVec3 (x,y,z)   = Vec3 x y z
+mkVec4 (x,y,z,w) = Vec4 x y z w
+
+-- | Orthogonal matrices.
+--
+-- Note: the "Random" instances generates orthogonal matrices with determinant 1
+-- (that is, orientation-preserving orthogonal transformations)!
+newtype Ortho2 a = Ortho2 (Mat2 a) deriving (Read,Show,Storable,MultSemiGroup,Determinant a,Dimension)
+newtype Ortho3 a = Ortho3 (Mat3 a) deriving (Read,Show,Storable,MultSemiGroup,Determinant a,Dimension)
+newtype Ortho4 a = Ortho4 (Mat4 a) deriving (Read,Show,Storable,MultSemiGroup,Determinant a,Dimension)
+
+-- | Projective matrices, encoding affine transformations in dimension one less.
+newtype Proj3 a = Proj3 (Mat3 a) deriving (Read,Show,Storable,MultSemiGroup)
+newtype Proj4 a = Proj4 (Mat4 a) deriving (Read,Show,Storable,MultSemiGroup)
+
+--------------------------------------------------------------------------------
+-- Unit vectors
+  
+instance Floating a => UnitVector a Vec2 Normal2 where
+  mkNormal v = Normal2 (normalize v)
+  fromNormal (Normal2 v) = v 
+  toNormalUnsafe = Normal2
+
+instance Floating a => UnitVector a Vec3 Normal3 where
+  mkNormal v = Normal3 (normalize v)
+  fromNormal (Normal3 v) = v 
+  toNormalUnsafe = Normal3
+
+instance Floating a => UnitVector a Vec4 Normal4 where
+  mkNormal v = Normal4 (normalize v)
+  fromNormal (Normal4 v) = v 
+  toNormalUnsafe = Normal4
+
+_rndUnit :: (Ord a, RandomGen g, Random (v a), Vector a v, DotProd a v) => g -> (v a,g)
+_rndUnit g = 
+  if d > 0.01
+    then ( v &* (1.0/d) , h )
+    else _rndUnit h
+  where
+    (v,h) = random g
+    d = norm v
+    
+instance (Floating a, Random a, Ord a) => Random (Normal2 a) where
+  random g = let (v,h) = _rndUnit g in (Normal2 v, h)  
+  randomR _ = random
+
+instance (Floating a, Random a, Ord a) => Random (Normal3 a) where
+  random g = let (v,h) = _rndUnit g in (Normal3 v, h)  
+  randomR _ = random
+
+instance (Floating a, Random a, Ord a) => Random (Normal4 a) where
+  random g = let (v,h) = _rndUnit g in (Normal4 v, h)  
+  randomR _ = random
+
+instance Floating a => CrossProd (Normal3 a) where
+  crossprod (Normal3 v) (Normal3 w) = mkNormal (crossprod v w)
+
+--------------------------------------------------------------------------------
+-- Orthogonal matrices
+
+instance Floating a =>  Orthogonal a Mat2 Ortho2 where
+  fromOrtho (Ortho2 o) = o
+  toOrthoUnsafe = Ortho2
+
+instance Floating a => Orthogonal a Mat3 Ortho3 where
+  fromOrtho (Ortho3 o) = o
+  toOrthoUnsafe = Ortho3 
+
+instance Floating a => Orthogonal a Mat4 Ortho4 where
+  fromOrtho (Ortho4 o) = o
+  toOrthoUnsafe = Ortho4
+
+------
+
+instance Floating a => Matrix (Ortho2 a) where
+  transpose (Ortho2 o) = Ortho2 (transpose o)
+  idmtx = Ortho2 idmtx
+  inverse = transpose
+
+instance Floating a => Matrix (Ortho3 a) where
+  transpose (Ortho3 o) = Ortho3 (transpose o)
+  idmtx = Ortho3 idmtx
+  inverse = transpose
+
+instance Floating a => Matrix (Ortho4 a) where
+  transpose (Ortho4 o) = Ortho4 (transpose o)
+  idmtx = Ortho4 idmtx
+  inverse = transpose
+
+------
+
+instance (Floating a, Ord a, Random a) => Random (Ortho2 a) where
+  random g = let (o,h) = _rndOrtho2 g in (toOrthoUnsafe (_flip1stRow2 o), h)
+  randomR _ = random
+
+instance (Floating a, Ord a, Random a) => Random (Ortho3 a) where
+  random g = let (o,h) = _rndOrtho3 g in (toOrthoUnsafe (             o), h)
+  randomR _ = random
+
+instance (Floating a, Ord a, Random a) => Random (Ortho4 a) where
+  random g = let (o,h) = _rndOrtho4 g in (toOrthoUnsafe (_flip1stRow4 o), h)
+  randomR _ = random
+
+------
+
+-- determinant will be -1
+_rndOrtho2 :: (Floating a, Random a, Ord a, RandomGen g) => g -> (Mat2 a, g)
+_rndOrtho2 g = (h2, g1) where
+  h2 = householder u2
+  (u2,g1) = random g   
+
+-- generates a uniformly random orthogonal 3x3 matrix 
+-- /with determinant +1/, with respect to the Haar measure of SO3.
+--
+-- see Theorem 4 in:
+-- Francesco Mezzadri: How to Generate Random Matrices from the Classical Compact Groups 
+-- Notices of the AMS, May 2007 issue
+-- <http://www.ams.org/notices/200705/fea-mezzadri-web.ps>
+_rndOrtho3 :: (Floating a, Random a, Ord a, RandomGen g) => g -> (Mat3 a, g) 
+_rndOrtho3 g = ( (h3 .*. m3), g2) where
+  m3 = (extendWith :: Floating a => a -> Mat2 a -> Mat3 a) 1 o2 
+  h3 = householder u3
+  (u3,g1) = random g
+  (o2,g2) = _rndOrtho2 g1
+
+-- determinant will be -1
+_rndOrtho4 :: (Floating a, Random a, Ord a, RandomGen g) => g -> (Mat4 a, g) 
+_rndOrtho4 g = ( (h4 .*. m4), g2) where
+  m4 = (extendWith :: Floating a => a -> Mat3 a -> Mat4 a) 1 o3 
+  h4 = householder u4
+  (u4,g1) = random g
+  (o3,g2) = _rndOrtho3 g1
+
+------
+
+_flip1stRow2 :: Floating a => Mat2 a -> Mat2 a
+_flip1stRow2 (Mat2 a b) = Mat2 (neg a) b
+
+_flip1stRow3 :: Floating a => Mat3 a -> Mat3 a
+_flip1stRow3 (Mat3 a b c) = Mat3 (neg a) b c
+
+_flip1stRow4 :: Floating a => Mat4 a -> Mat4 a
+_flip1stRow4 (Mat4 a b c d) = Mat4 (neg a) b c d
+
+--------------------------------------------------------------------------------
+-- projective matrices
+  
+instance Floating a => Projective a Vec2 Mat2 Ortho2 Mat3 Proj3 where
+  fromProjective (Proj3 m) = m
+  toProjectiveUnsafe = Proj3
+  orthogonal = Proj3 . extendWith (1 :: a) . fromOrtho
+  linear     = Proj3 . extendWith (1 :: a)
+  translation v = Proj3 $ Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (extendWith (1 :: a) v)
+  scaling     v = Proj3 $ diag (extendWith (1 :: a) v)
+  
+instance Floating a => Projective a Vec3 Mat3 Ortho3 Mat4 Proj4 where
+  fromProjective (Proj4 m) = m
+  toProjectiveUnsafe = Proj4
+  orthogonal = Proj4 . extendWith (1 :: a) . fromOrtho 
+  linear     = Proj4 . extendWith (1 :: a)
+  translation v = Proj4 $ Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (extendWith (1 :: a) v)
+  scaling     v = Proj4 $ diag (extendWith (1 :: a) v)
+
+instance Floating a => Matrix (Proj3 a) where
+  idmtx = Proj3 idmtx
+  transpose (Proj3 m) = Proj3 (transpose m)
+  inverse = _invertProj3
+
+instance Floating a => Matrix (Proj4 a) where
+  idmtx = Proj4 idmtx
+  transpose (Proj4 m) = Proj4 (transpose m)
+  inverse = _invertProj4
+
+_invertProj3 :: (Floating a, Extend a Vec2 Vec3) => Proj3 a -> Proj3 a
+_invertProj3 (Proj3 mat@(Mat3 _ _ t)) = 
+  Proj3 $ Mat3 (extendZero a) (extendZero b) (extendWith 1 t')
+  where
+    t' = neg $ (trim :: Extend a Vec2 Vec3 => Vec3 a -> Vec2 a) t .* invm2
+    invm2@(Mat2 a b) = inverse $ trim mat
+
+-- Inverts a projective 4x4 matrix. But you can simply use "inverse" instead.
+-- We assume that the bottom-right corner is 1.
+_invertProj4 :: Floating a => Proj4 a -> Proj4 a
+_invertProj4 (Proj4 mat@(Mat4 _ _ _ t)) = 
+  Proj4 $ Mat4 (extendZero a) (extendZero b) (extendZero c) (extendWith 1 t') 
+  where
+    t' = neg $ (trim :: Extend a Vec3 Vec4 => Vec4 a -> Vec3 a) t .* invm3 
+    invm3@(Mat3 a b c) = inverse $ trim mat
+
+--------------------------------------------------------------------------------
+-- Vec2 instances
+
+instance Floating a => HasCoordinates (Vec2 a) a where
+  _1 (Vec2 x _) = x
+  _2 (Vec2 _ y) = y
+  _3 _ = error "has only 2 coordinates"
+  _4 _ = error "has only 2 coordinates"
+
+instance Floating a => AbelianGroup (Vec2 a) where
+  (&+) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1+x2) (y1+y2) 
+  (&-) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1-x2) (y1-y2)
+  neg  (Vec2 x y)                = Vec2 (-x) (-y)
+  zero = Vec2 0 0
+  
+instance Floating a => Vector a Vec2 where
+  scalarMul s (Vec2 x y) = Vec2 (s*x) (s*y)
+  mapVec    f (Vec2 x y) = Vec2 (f x) (f y)
+  
+instance Floating a => DotProd a Vec2 where
+  (&.) (Vec2 x1 y1) (Vec2 x2 y2) = x1*x2 + y1*y2
+
+instance Floating a => Pointwise (Vec2 a) where
+  pointwise (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1*x2) (y1*y2)
+
+instance Floating a => Determinant a (Vec2 a, Vec2 a) where
+  det (Vec2 x1 y1 , Vec2 x2 y2) = x1*y2 - x2*y1
+
+{-     
+instance Show Vec2 where
+  show (Vec2 x y) = "( " ++ show x ++ " , " ++ show y ++ " )"
+-}
+
+instance (Floating a, Random a) => Random (Vec2 a) where
+  random = randomR (Vec2 (-1) (-1),Vec2 1 1)
+  randomR (Vec2 a b, Vec2 c d) gen = 
+    let (x,gen1) = randomR (a,c) gen
+        (y,gen2) = randomR (b,d) gen1
+    in (Vec2 x y, gen2)
+     
+instance (Floating a, Storable a) => Storable (Vec2 a) where
+  sizeOf    _ = 2 * sizeOf (undefined :: a)
+  alignment _ = sizeOf (undefined :: a)
+  
+  peek q = do
+    let p = castPtr q :: Ptr a
+        k = sizeOf (undefined :: a)
+    x <- peek        p 
+    y <- peekByteOff p k
+    return (Vec2 x y)
+    
+  poke q (Vec2 x y) = do
+    let p = castPtr q :: Ptr a
+        k = sizeOf (undefined :: a)
+    poke        p   x
+    pokeByteOff p k y
+
+instance Floating a => Dimension (Vec2 a) where dim _ = 2
+
+--------------------------------------------------------------------------------                    
+-- Mat2 instances
+
+instance Floating a => HasCoordinates (Mat2 a) (Vec2 a) where
+  _1 (Mat2 x _) = x
+  _2 (Mat2 _ y) = y
+  _3 _ = error "has only 2 coordinates"
+  _4 _ = error "has only 2 coordinates"
+
+instance Floating a => Matrix (Mat2 a) where
+  transpose (Mat2 row1 row2) = 
+    Mat2 (Vec2 (_1 row1) (_1 row2)) 
+         (Vec2 (_2 row1) (_2 row2)) 
+  idmtx = Mat2 (Vec2 1 0) (Vec2 0 1)
+  inverse (Mat2 (Vec2 a b) (Vec2 c d)) = 
+    Mat2 (Vec2 (d*r) (-b*r)) (Vec2 (-c*r) (a*r)) 
+    where r = 1.0 / (a*d - b*c)
+
+instance Floating a => AbelianGroup (Mat2 a) where
+  (&+) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &+ s1) (r2 &+ s2)
+  (&-) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &- s1) (r2 &- s2)
+  neg  (Mat2 r1 r2)              = Mat2 (neg r1) (neg r2)  
+  zero = Mat2 zero zero  
+  
+instance Floating a => Vector a Mat2 where
+  scalarMul s (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = scalarMul s
+  mapVec    f (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = mapVec f
+
+instance Floating a => MultSemiGroup (Mat2 a) where
+  (.*.) (Mat2 r1 r2) n = 
+    let (Mat2 c1 c2) = transpose n
+    in Mat2 (Vec2 (r1 &. c1) (r1 &. c2))
+            (Vec2 (r2 &. c1) (r2 &. c2))
+  one = idmtx 
+
+instance Floating a => Ring (Mat2 a)
+
+instance Floating a => LeftModule (Mat2 a) (Vec2 a) where
+  lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) 
+  
+instance Floating a => RightModule (Vec2 a) (Mat2 a) where
+  rmul v mt = lmul (transpose mt) v
+
+instance Floating a => Diagonal (Vec2 a) (Mat2 a) where
+  diag (Vec2 x y) = Mat2 (Vec2 x 0) (Vec2 0 y)
+
+instance Floating a => Tensor (Mat2 a) (Vec2 a) where
+  outer (Vec2 a b) (Vec2 x y) = Mat2
+    (Vec2 (a*x) (a*y))
+    (Vec2 (b*x) (b*y))
+
+instance Floating a => Determinant a (Mat2 a) where
+  det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c 
+
+{-
+instance Show Mat2 where
+  show (Mat2 r1 r2) = show r1 ++ "\n" ++ show r2
+-}
+
+instance (Floating a, Storable a) => Storable (Mat2 a) where
+  sizeOf    _ = 2 * sizeOf (undefined :: Vec2 a)
+  alignment _ = alignment  (undefined :: Vec2 a)
+  
+  peek q = do
+    let p = castPtr q :: Ptr (Vec2 a)
+        k = sizeOf (undefined :: Vec2 a)
+    r1 <- peek        p 
+    r2 <- peekByteOff p k
+    return (Mat2 r1 r2)
+    
+  poke q (Mat2 r1 r2) = do
+    let p = castPtr q :: Ptr (Vec2 a)
+        k = sizeOf (undefined :: Vec2 a)
+    poke        p   r1
+    pokeByteOff p k r2
+
+instance (Floating a, Random a) => Random (Mat2 a) where
+  random = randomR (Mat2 v1 v1 , Mat2 v2 v2) where 
+    v1 = Vec2 (-1) (-1) 
+    v2 = Vec2   1    1
+  randomR (Mat2 a b, Mat2 c d) gen = 
+    let (x,gen1) = randomR (a,c) gen
+        (y,gen2) = randomR (b,d) gen1
+    in (Mat2 x y, gen2)
+          
+instance Floating a => Dimension (Mat2 a) where dim _ = 2
+     
+instance (Floating a) => MatrixNorms a (Mat2 a) where 
+  frobeniusNorm (Mat2 r1 r2) =  
+    sqrt $
+      normsqr r1 + 
+      normsqr r2
+     
+instance Floating a => Pointwise (Mat2 a) where
+  pointwise (Mat2 x1 y1) (Mat2 x2 y2) = Mat2 (x1 &! x2) (y1 &! y2)
+       
+--------------------------------------------------------------------------------     
+-- Vec3 instances
+
+instance Floating a => HasCoordinates (Vec3 a) a where
+  _1 (Vec3 x _ _) = x
+  _2 (Vec3 _ y _) = y
+  _3 (Vec3 _ _ z) = z
+  _4 _ = error "has only 3 coordinates"
+
+instance Floating a => AbelianGroup (Vec3 a) where
+  (&+) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1+x2) (y1+y2) (z1+z2) 
+  (&-) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1-x2) (y1-y2) (z1-z2) 
+  neg  (Vec3 x y z)                    = Vec3 (-x) (-y) (-z)
+  zero = Vec3 0 0 0
+  
+instance Floating a => Vector a Vec3 where
+  scalarMul s (Vec3 x y z) = Vec3 (s*x) (s*y) (s*z)
+  mapVec    f (Vec3 x y z) = Vec3 (f x) (f y) (f z)
+
+instance Floating a => DotProd a Vec3 where
+  (&.) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = x1*x2 + y1*y2 + z1*z2
+
+instance Num a => Pointwise (Vec3 a) where
+  pointwise (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1*x2) (y1*y2) (z1*z2)
+
+{-
+instance Show Vec3 where
+  show (Vec3 x y z) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " )"
+-}
+
+instance (Floating a, Random a) => Random (Vec3 a) where
+  random = randomR (Vec3 (-1) (-1) (-1),Vec3 1 1 1)
+  randomR (Vec3 a b c, Vec3 d e f) gen = 
+    let (x,gen1) = randomR (a,d) gen
+        (y,gen2) = randomR (b,e) gen1
+        (z,gen3) = randomR (c,f) gen2  
+    in (Vec3 x y z, gen3)
+      
+instance Floating a => CrossProd (Vec3 a) where
+  crossprod (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (y1*z2-y2*z1) (z1*x2-z2*x1) (x1*y2-x2*y1) 
+
+instance Floating a => Determinant a (Vec3 a, Vec3 a, Vec3 a) where
+  det (u,v,w) = u &. (v &^ w)  
+ 
+instance (Floating a, Storable a) => Storable (Vec3 a) where
+  sizeOf    _ = 3 * sizeOf (undefined :: a)
+  alignment _ = sizeOf (undefined :: a)
+  
+  peek q = do
+    let p = castPtr q :: Ptr a
+        k = sizeOf (undefined :: a)
+    x <- peek        p 
+    y <- peekByteOff p (k  )
+    z <- peekByteOff p (k+k)
+    return (Vec3 x y z)
+    
+  poke q (Vec3 x y z) = do
+    let p = castPtr q :: Ptr a
+        k = sizeOf (undefined :: a)
+    poke        p       x
+    pokeByteOff p (k  ) y
+    pokeByteOff p (k+k) z
+
+instance Floating a => Dimension (Vec3 a) where dim _ = 3
+
+--------------------------------------------------------------------------------   
+-- Mat3 instances
+
+instance Floating a => HasCoordinates (Mat3 a) (Vec3 a) where
+  _1 (Mat3 x _ _) = x
+  _2 (Mat3 _ y _) = y
+  _3 (Mat3 _ _ z) = z
+  _4 _ = error "has only 3 coordinates"  
+
+instance Floating a => Matrix (Mat3 a) where
+
+  transpose (Mat3 row1 row2 row3) = 
+    Mat3 (Vec3 (_1 row1) (_1 row2) (_1 row3)) 
+         (Vec3 (_2 row1) (_2 row2) (_2 row3)) 
+         (Vec3 (_3 row1) (_3 row2) (_3 row3)) 
+         
+  idmtx = Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (Vec3 0 0 1)
+  
+  inverse (Mat3 (Vec3 a b c) (Vec3 e f g) (Vec3 i j k)) = 
+    Mat3 (Vec3 (d11*r) (d21*r) (d31*r))  
+         (Vec3 (d12*r) (d22*r) (d32*r))  
+         (Vec3 (d13*r) (d23*r) (d33*r))  
+    where
+      r = 1.0 / ( a*d11 + b*d12 + c*d13 )
+
+      d11 = f*k - g*j
+      d12 = g*i - e*k
+      d13 = e*j - f*i
+
+      d31 = b*g - c*f
+      d32 = c*e - a*g
+      d33 = a*f - b*e
+
+      d21 = c*j - b*k 
+      d22 = a*k - c*i 
+      d23 = b*i - a*j 
+
+instance Floating a => AbelianGroup (Mat3 a) where
+  (&+) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3)
+  (&-) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &- s1) (r2 &- s2) (r3 &- s3)
+  neg  (Mat3 r1 r2 r3)                 = Mat3 (neg r1) (neg r2) (neg r3) 
+  zero = Mat3 zero zero zero 
+
+instance Floating a => Vector a Mat3 where
+  scalarMul s (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = scalarMul s
+  mapVec    f (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = mapVec f
+
+instance Floating a => MultSemiGroup (Mat3 a) where
+  (.*.) (Mat3 r1 r2 r3) n = 
+    let (Mat3 c1 c2 c3) = transpose n
+    in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))
+            (Vec3 (r2 &. c1) (r2 &. c2) (r2 &. c3))
+            (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))
+  one = idmtx 
+
+instance Floating a => Ring (Mat3 a)
+
+instance Floating a => LeftModule (Mat3 a) (Vec3 a) where
+  lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)
+  
+instance Floating a => RightModule (Vec3 a) (Mat3 a) where
+  rmul v mt = lmul (transpose mt) v
+
+instance Floating a => Diagonal (Vec3 a) (Mat3 a) where
+  diag (Vec3 x y z) = Mat3 (Vec3 x 0 0) (Vec3 0 y 0) (Vec3 0 0 z)
+
+instance Floating a => Tensor (Mat3 a) (Vec3 a) where
+  outer (Vec3 a b c) (Vec3 x y z) = Mat3
+    (Vec3 (a*x) (a*y) (a*z))
+    (Vec3 (b*x) (b*y) (b*z))
+    (Vec3 (c*x) (c*y) (c*z))
+
+instance Floating a => Determinant a (Mat3 a) where
+  det (Mat3 r1 r2 r3) = det (r1,r2,r3)
+
+{-
+instance Show Mat3 where
+  show (Mat3 r1 r2 r3) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3
+-}
+
+instance (Floating a, Storable a) => Storable (Mat3 a) where
+  sizeOf    _ = 3 * sizeOf (undefined::Vec3 a)
+  alignment _ = alignment  (undefined::Vec3 a)
+  
+  peek q = do
+    let p = castPtr q :: Ptr (Vec3 a)
+        k = sizeOf (undefined::Vec3 a)
+    r1 <- peek        p 
+    r2 <- peekByteOff p (k  )
+    r3 <- peekByteOff p (k+k)
+    return (Mat3 r1 r2 r3)
+    
+  poke q (Mat3 r1 r2 r3) = do
+    let p = castPtr q :: Ptr (Vec3 a)
+        k = sizeOf (undefined::Vec3 a)
+    poke        p       r1
+    pokeByteOff p (k  ) r2
+    pokeByteOff p (k+k) r3
+
+instance (Floating a, Random a) => Random (Mat3 a) where
+  random = randomR (Mat3 v1 v1 v1 , Mat3 v2 v2 v2) where
+    v1 = Vec3 (-1) (-1) (-1)
+    v2 = Vec3   1    1    1
+  randomR (Mat3 a b c, Mat3 d e f) gen = 
+    let (x,gen1) = randomR (a,d) gen
+        (y,gen2) = randomR (b,e) gen1
+        (z,gen3) = randomR (c,f) gen2  
+    in (Mat3 x y z, gen3)
+   
+instance Floating a => Dimension (Mat3 a) where dim _ = 3
+  
+instance Floating a => MatrixNorms a (Mat3 a) where 
+  frobeniusNorm (Mat3 r1 r2 r3)  = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 
+
+instance Floating a => Pointwise (Mat3 a) where
+  pointwise (Mat3 x1 y1 z1) (Mat3 x2 y2 z2) = Mat3 (x1 &! x2) (y1 &! y2) (z1 &! z2)
+    
+--------------------------------------------------------------------------------
+-- Vec4 instances
+
+instance Floating a => HasCoordinates (Vec4 a) a where
+  _1 (Vec4 x _ _ _) = x
+  _2 (Vec4 _ y _ _) = y
+  _3 (Vec4 _ _ z _) = z
+  _4 (Vec4 _ _ _ w) = w
+
+instance Floating a => AbelianGroup (Vec4 a) where
+  (&+) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1+x2) (y1+y2) (z1+z2) (w1+w2)
+  (&-) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1-x2) (y1-y2) (z1-z2) (w1-w2)
+  neg  (Vec4 x y z w)                        = Vec4 (-x) (-y) (-z) (-w)
+  zero = Vec4 0 0 0 0
+  
+instance Floating a => Vector a Vec4 where
+  scalarMul s (Vec4 x y z w) = Vec4 (s*x) (s*y) (s*z) (s*w)
+  mapVec    f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w)
+
+instance Floating a => DotProd a Vec4 where
+  (&.) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = x1*x2 + y1*y2 + z1*z2 + w1*w2
+
+instance Num a => Pointwise (Vec4 a) where
+  pointwise (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1*x2) (y1*y2) (z1*z2) (w1*w2)
+
+{-
+instance Show Vec4 where
+  show (Vec4 x y z w) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " , " ++ show w ++ " )"
+-}
+
+instance (Floating a, Random a) => Random (Vec4 a) where
+  random = randomR (Vec4 (-1) (-1) (-1) (-1),Vec4 1 1 1 1)
+  randomR (Vec4 a b c d, Vec4 e f g h) gen = 
+    let (x,gen1) = randomR (a,e) gen
+        (y,gen2) = randomR (b,f) gen1
+        (z,gen3) = randomR (c,g) gen2  
+        (w,gen4) = randomR (d,h) gen3  
+    in (Vec4 x y z w, gen4)
+           
+instance (Floating a, Storable a) => Storable (Vec4 a) where
+  sizeOf    _ = 4 * sizeOf (undefined :: a)
+  alignment _ = sizeOf (undefined :: a)
+  
+  peek q = do
+    let p = castPtr q :: Ptr a
+        k = sizeOf (undefined :: a)
+    x <- peek        p 
+    y <- peekByteOff p (k  )
+    z <- peekByteOff p (k+k)
+    w <- peekByteOff p (3*k)
+    return (Vec4 x y z w)
+    
+  poke q (Vec4 x y z w) = do
+    let p = castPtr q :: Ptr a
+        k = sizeOf (undefined :: a)
+    poke        p       x
+    pokeByteOff p (k  ) y
+    pokeByteOff p (k+k) z
+    pokeByteOff p (3*k) w
+
+instance Floating a => Dimension (Vec4 a) where dim _ = 4
+
+--------------------------------------------------------------------------------
+-- Mat4 instances
+
+instance Floating a => HasCoordinates (Mat4 a) (Vec4 a) where
+  _1 (Mat4 x _ _ _) = x
+  _2 (Mat4 _ y _ _) = y
+  _3 (Mat4 _ _ z _) = z
+  _4 (Mat4 _ _ _ w) = w
+
+instance Floating a => Matrix (Mat4 a) where
+  transpose (Mat4 row1 row2 row3 row4) = 
+    Mat4 (Vec4 (_1 row1) (_1 row2) (_1 row3) (_1 row4)) 
+         (Vec4 (_2 row1) (_2 row2) (_2 row3) (_2 row4)) 
+         (Vec4 (_3 row1) (_3 row2) (_3 row3) (_3 row4)) 
+         (Vec4 (_4 row1) (_4 row2) (_4 row3) (_4 row4)) 
+  idmtx = Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (Vec4 0 0 0 1)
+  inverse = error "inverse/Mat4: not implemented yet"
+
+instance Floating a => AbelianGroup (Mat4 a) where
+  (&+) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3) (r4 &+ s4)
+  (&-) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &- s1) (r2 &- s2) (r3 &- s3) (r4 &- s4)
+  neg  (Mat4 r1 r2 r3 r4)                    = Mat4 (neg r1) (neg r2) (neg r3) (neg r4) 
+  zero = Mat4 zero zero zero zero
+  
+instance Floating a => Vector a Mat4 where
+  scalarMul s (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = scalarMul s
+  mapVec    f (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = mapVec f
+
+instance Floating a => MultSemiGroup (Mat4 a) where
+  (.*.) (Mat4 r1 r2 r3 r4) n = 
+    let (Mat4 c1 c2 c3 c4) = transpose n
+    in Mat4 (Vec4 (r1 &. c1) (r1 &. c2) (r1 &. c3) (r1 &. c4))
+            (Vec4 (r2 &. c1) (r2 &. c2) (r2 &. c3) (r2 &. c4))
+            (Vec4 (r3 &. c1) (r3 &. c2) (r3 &. c3) (r3 &. c4))
+            (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))
+  one = idmtx 
+
+instance Floating a => Ring (Mat4 a)
+
+instance Floating a => LeftModule (Mat4 a) (Vec4 a) where
+  lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)
+  
+instance Floating a => RightModule (Vec4 a) (Mat4 a) where
+  rmul v mt = lmul (transpose mt) v
+
+instance Floating a => Diagonal (Vec4 a) (Mat4 a) where
+  diag (Vec4 x y z w) = Mat4 (Vec4 x 0 0 0) (Vec4 0 y 0 0) (Vec4 0 0 z 0) (Vec4 0 0 0 w)
+
+instance Floating a => Tensor (Mat4 a) (Vec4 a) where
+  outer (Vec4 a b c d) (Vec4 x y z w) = Mat4
+    (Vec4 (a*x) (a*y) (a*z) (a*w))
+    (Vec4 (b*x) (b*y) (b*z) (b*w))
+    (Vec4 (c*x) (c*y) (c*z) (c*w))
+    (Vec4 (d*x) (d*y) (d*z) (d*w))
+
+instance Floating a => Determinant a (Mat4 a) where
+  det = error "det/Mat4: not implemented yet" 
+  -- det (Mat4 r1 r2 r3 r4) = 
+
+{-
+instance Show Mat4 where
+  show (Mat4 r1 r2 r3 r4) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 ++ "\n" ++ show r4
+-}
+
+instance (Floating a, Storable a) => Storable (Mat4 a) where
+  sizeOf    _ = 4 * sizeOf (undefined::Vec4 a)
+  alignment _ = alignment  (undefined::Vec4 a)
+  
+  peek q = do
+    let p = castPtr q :: Ptr (Vec4 a)
+        k = sizeOf (undefined :: Vec4 a)
+    r1 <- peek        p 
+    r2 <- peekByteOff p (k  )
+    r3 <- peekByteOff p (k+k)
+    r4 <- peekByteOff p (3*k)
+    return (Mat4 r1 r2 r3 r4)
+    
+  poke q (Mat4 r1 r2 r3 r4) = do
+    let p = castPtr q :: Ptr (Vec4 a)
+        k = sizeOf (undefined :: Vec4 a)
+    poke        p       r1
+    pokeByteOff p (k  ) r2
+    pokeByteOff p (k+k) r3
+    pokeByteOff p (3*k) r4
+
+instance (Floating a, Random a) => Random (Mat4 a) where
+  random = randomR (Mat4 v1 v1 v1 v1, Mat4 v2 v2 v2 v2) where
+    v1 = Vec4 (-1) (-1) (-1) (-1)
+    v2 = Vec4   1    1    1    1
+  randomR (Mat4 a b c d, Mat4 e f g h) gen = 
+    let (x,gen1) = randomR (a,e) gen
+        (y,gen2) = randomR (b,f) gen1
+        (z,gen3) = randomR (c,g) gen2  
+        (w,gen4) = randomR (d,h) gen3  
+    in (Mat4 x y z w, gen4)
+    
+instance Floating a => Dimension (Mat4 a) where dim _ = 4
+   
+instance Floating a => MatrixNorms a (Mat4 a) where 
+  frobeniusNorm (Mat4 r1 r2 r3 r4) = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 + 
+      normsqr r4  
+    
+instance Floating a => Pointwise (Mat4 a) where
+  pointwise (Mat4 x1 y1 z1 w1) (Mat4 x2 y2 z2 w2) = Mat4 (x1 &! x2) (y1 &! y2) (z1 &! z2) (w1 &! w2)
+    
+--------------------------------------------------------------------------------
+-- Extend instances
+
+instance Floating a => Extend a Vec2 Vec3 where
+  extendZero   (Vec2 x y) = Vec3 x y 0
+  extendWith t (Vec2 x y) = Vec3 x y t
+  trim (Vec3 x y _)       = Vec2 x y
+
+instance Floating a => Extend a Vec2 Vec4 where
+  extendZero   (Vec2 x y) = Vec4 x y 0 0
+  extendWith t (Vec2 x y) = Vec4 x y t t
+  trim (Vec4 x y _ _)     = Vec2 x y 
+
+instance Floating a => Extend a Vec3 Vec4 where
+  extendZero   (Vec3 x y z) = Vec4 x y z 0
+  extendWith t (Vec3 x y z) = Vec4 x y z t
+  trim (Vec4 x y z _)       = Vec3 x y z
+
+instance Floating a => Extend a Mat2 Mat3 where
+  extendZero   (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero
+  extendWith w (Mat2 p q) = Mat3 (extendZero p) (extendZero q) (Vec3 0 0 w)
+  trim (Mat3 p q _) = Mat2 (trim p) (trim q)
+
+instance Floating a => Extend a Mat2 Mat4 where
+  extendZero   (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero
+  extendWith w (Mat2 p q) = Mat4 (extendZero p) (extendZero q) (Vec4 0 0 w 0) (Vec4 0 0 0 w)
+  trim (Mat4 p q _ _) = Mat2 (trim p) (trim q)
+
+instance Floating a => Extend a Mat3 Mat4 where
+  extendZero   (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero
+  extendWith w (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) (Vec4 0 0 0 w)
+  trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r)
+  
+--------------------------------------------------------------------------------
+
+  
diff --git a/src/Data/Vect/Floating/GramSchmidt.hs b/src/Data/Vect/Floating/GramSchmidt.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/GramSchmidt.hs
@@ -0,0 +1,135 @@
+{-# LANGUAGE FlexibleInstances #-}
+
+-- | Gram-Schmidt orthogonalization.
+-- This module is not re-exported by "Data.Vect".
+
+module Data.Vect.Floating.GramSchmidt 
+  ( GramSchmidt(..)
+  )
+  where
+
+import Data.Vect.Floating.Base
+
+--------------------------------------------------------------------------------
+
+liftPair :: (a -> b) -> (a,a) -> (b,b)
+liftPair f (x,y) = (f x, f y)
+
+liftTriple :: (a -> b) -> (a,a,a) -> (b,b,b)
+liftTriple f (x,y,z) = (f x, f y, f z)
+
+liftQuadruple :: (a -> b) -> (a,a,a,a) -> (b,b,b,b)
+liftQuadruple f (x,y,z,w) = (f x, f y, f z, f w)
+
+--------------------------------------------------------------------------------
+    
+-- | produces orthogonal\/orthonormal vectors from a set of vectors    
+class GramSchmidt a where
+  gramSchmidt          :: a -> a   -- ^ does not normalize the vectors!
+  gramSchmidtNormalize :: a -> a   -- ^ normalizes the vectors.
+
+{-# RULES
+"gramSchmidt is idempotent"  forall a. gramSchmidt (gramSchmidt a) = gramSchmidt a 
+"gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a 
+  #-}
+
+--------------------------------------------------------------------------------
+
+instance Floating a => GramSchmidt (Vec2 a,Vec2 a) where
+  gramSchmidt = gramSchmidtPair
+  gramSchmidtNormalize = gramSchmidtNormalizePair
+  
+instance Floating a => GramSchmidt (Vec3 a,Vec3 a) where
+  gramSchmidt = gramSchmidtPair
+  gramSchmidtNormalize = gramSchmidtNormalizePair
+  
+instance Floating a => GramSchmidt (Vec4 a,Vec4 a) where
+  gramSchmidt = gramSchmidtPair
+  gramSchmidtNormalize = gramSchmidtNormalizePair
+
+----------
+
+instance Floating a => GramSchmidt (Normal2 a,Normal2 a) where
+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal2!"
+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal
+
+instance Floating a => GramSchmidt (Normal3 a,Normal3 a) where
+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"
+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal
+
+instance Floating a => GramSchmidt (Normal4 a,Normal4 a) where
+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"
+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal
+
+----------
+  
+gramSchmidtPair :: (Vector a v, DotProd a v) => (v a,v a) -> (v a,v a)
+gramSchmidtPair (u,v) = (u',v') where 
+  u' = u
+  v' = project v u'     
+  
+gramSchmidtNormalizePair :: (Vector a v, DotProd a v) => (v a,v a) -> (v a,v a)
+gramSchmidtNormalizePair (u,v) = (u',v') where
+  u' = normalize u 
+  v' = normalize $ projectUnsafe v u'     
+
+----------
+
+instance Floating a => GramSchmidt (Vec3 a,Vec3 a,Vec3 a) where
+  gramSchmidt = gramSchmidtTriple
+  gramSchmidtNormalize = gramSchmidtNormalizeTriple
+     
+instance Floating a => GramSchmidt (Vec4 a,Vec4 a,Vec4 a) where
+  gramSchmidt = gramSchmidtTriple
+  gramSchmidtNormalize = gramSchmidtNormalizeTriple
+
+instance Floating a => GramSchmidt (Normal3 a,Normal3 a,Normal3 a) where
+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"
+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal
+
+instance Floating a => GramSchmidt (Normal4 a,Normal4 a,Normal4 a) where
+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"
+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal
+
+----------
+
+gramSchmidtTriple :: (Vector a v, DotProd a v) => (v a,v a,v a) -> (v a,v a,v a)
+gramSchmidtTriple (u,v,w) = (u',v',w') where 
+  u' = u
+  v' = project v u'     
+  w' = project (project w u') v' 
+  
+gramSchmidtNormalizeTriple :: (Vector a v, DotProd a v) => (v a,v a,v a) -> (v a,v a,v a)
+gramSchmidtNormalizeTriple (u,v,w) = (u',v',w') where
+  u' = normalize $ u 
+  v' = normalize $ projectUnsafe v u'     
+  w' = normalize $ projectUnsafe (projectUnsafe w u') v'     
+
+----------
+
+instance Floating a => GramSchmidt (Vec4 a,Vec4 a,Vec4 a,Vec4 a) where
+  gramSchmidt          = gramSchmidtQuadruple
+  gramSchmidtNormalize = gramSchmidtNormalizeQuadruple 
+
+instance Floating a => GramSchmidt (Normal4 a,Normal4 a,Normal4 a,Normal4 a) where
+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"
+  gramSchmidtNormalize = liftQuadruple toNormalUnsafe . gramSchmidtNormalizeQuadruple . liftQuadruple fromNormal
+
+----------
+  
+gramSchmidtQuadruple :: (Vector a v, DotProd a v) => (v a,v a,v a,v a) -> (v a,v a,v a,v a)
+gramSchmidtQuadruple (u,v,w,z) = (u',v',w',z') where 
+  u' = u
+  v' = project v u'     
+  w' = project (project w u') v' 
+  z' = project (project (project z u') v') w'
+
+gramSchmidtNormalizeQuadruple :: (Vector a v, DotProd a v) => (v a,v a,v a,v a) -> (v a,v a,v a,v a)
+gramSchmidtNormalizeQuadruple (u,v,w,z) = (u',v',w',z') where
+  u' = normalize $ u
+  v' = normalize $ projectUnsafe v u'     
+  w' = normalize $ projectUnsafe (projectUnsafe w u') v' 
+  z' = normalize $ projectUnsafe (projectUnsafe (projectUnsafe z u') v') w'
+  
+----------
+  
diff --git a/src/Data/Vect/Floating/Interpolate.hs b/src/Data/Vect/Floating/Interpolate.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/Interpolate.hs
@@ -0,0 +1,83 @@
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+
+-- TODO: interpolation for Ortho3 matrices using the (short) quaternion 'slerpU'
+
+-- | Interpolation of vectors. 
+-- Note: we interpolate unit vectors differently from ordinary vectors.
+
+module Data.Vect.Floating.Interpolate where
+
+--------------------------------------------------------------------------------
+
+import Data.Vect.Floating.Base
+import Data.Vect.Floating.Util.Dim2 (sinCos',angle2')
+import Data.Vect.Floating.Util.Dim3 (rotate3')
+
+--------------------------------------------------------------------------------
+
+class Interpolate a v where
+  interpolate :: a -> v -> v -> v
+  
+instance Num a => Interpolate a a where
+  interpolate t x y = x + t*(y-x)
+
+--------------------------------------------------------------------------------
+
+instance Floating a => Interpolate a (Vec2 a) where interpolate t x y = x &+ t *& (y &- x)
+instance Floating a => Interpolate a (Vec3 a) where interpolate t x y = x &+ t *& (y &- x)
+instance Floating a => Interpolate a (Vec4 a) where interpolate t x y = x &+ t *& (y &- x)
+
+--------------------------------------------------------------------------------
+
+{-
+instance Interpolate Normal2 where
+  interpolate t nx ny = sinCos' $ ax + t*adiff where
+    ax = angle2' nx
+    ay = angle2' ny
+    adiff = helper (ay - ax)
+    helper d 
+      | d < -pi   = d + twopi
+      | d >  pi   = d - twopi
+      | otherwise = d
+    twopi = 2*pi
+    
+instance Interpolate Normal3 where 
+  interpolate t nx ny = 
+    if maxAngle < 0.001  -- more or less ad-hoc critical angle
+      then mkNormal $ interpolate t x y
+      else toNormalUnsafe $ rotate3' (t*maxAngle) (mkNormal axis) x where
+    x = fromNormal nx
+    y = fromNormal ny
+    axis = (x &^ y)
+    maxAngle = acos (x &. y)
+-}        
+
+instance Floating a => Interpolate a (Normal2 a) where interpolate = slerp
+instance Floating a => Interpolate a (Normal3 a) where interpolate = slerp
+instance Floating a => Interpolate a (Normal4 a) where interpolate = slerp
+        
+--------------------------------------------------------------------------------
+  
+{-# SPECIALIZE slerp :: Float -> Normal2 Float -> Normal2 Float -> Normal2 Float #-}
+{-# SPECIALIZE slerp :: Float -> Normal3 Float -> Normal3 Float -> Normal3 Float #-}
+{-# SPECIALIZE slerp :: Float -> Normal4 Float -> Normal4 Float -> Normal4 Float #-}
+    
+{-# SPECIALIZE slerp :: Double -> Normal2 Double -> Normal2 Double -> Normal2 Double #-}
+{-# SPECIALIZE slerp :: Double -> Normal3 Double -> Normal3 Double -> Normal3 Double #-}
+{-# SPECIALIZE slerp :: Double -> Normal4 Double -> Normal4 Double -> Normal4 Double #-}
+
+-- | Spherical linear interpolation.
+-- See <http://en.wikipedia.org/wiki/Slerp>    
+slerp :: UnitVector a v u => a -> u a -> u a -> u a
+slerp t n0 n1 = toNormalUnsafe v where
+  v = (p0 &* y0) &+ (p1 &* y1) 
+  p0 = fromNormal n0
+  p1 = fromNormal n1
+  omega = acos (p0 &. p1)
+  s = sin omega
+  y0 = sin (omega*(1-t)) / s 
+  y1 = sin (omega*   t ) / s
+  
+--------------------------------------------------------------------------------
+
+  
diff --git a/src/Data/Vect/Floating/Util/Dim2.hs b/src/Data/Vect/Floating/Util/Dim2.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/Util/Dim2.hs
@@ -0,0 +1,69 @@
+module Data.Vect.Floating.Util.Dim2 where
+
+import Data.Vect.Floating.Base
+
+-- | Example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.
+structVec2 :: [a] -> [Vec2 a]
+structVec2 [] = []
+structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls) 
+structVec2 _ = error "structVec2"
+
+-- | The opposite of "structVec2".
+destructVec2 :: [Vec2 a] -> [a]
+destructVec2 [] = []
+destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  
+
+det2 :: Floating a => Vec2 a -> Vec2 a -> a
+det2 u v = det (u,v)
+
+vec2X :: Num a => Vec2 a
+vec2Y :: Num a => Vec2 a
+
+vec2X = Vec2 1 0 
+vec2Y = Vec2 0 1 
+
+translate2X :: Num a => a -> Vec2 a -> Vec2 a
+translate2Y :: Num a => a -> Vec2 a -> Vec2 a
+
+translate2X t (Vec2 x y) = Vec2 (x+t) y 
+translate2Y t (Vec2 x y) = Vec2 x (y+t) 
+
+-- | unit vector with given angle relative to the positive X axis (in the positive direction, that is, CCW).
+-- A more precise name would be @cosSin@, but that sounds bad :)
+sinCos :: Floating a => a -> Vec2 a
+sinCos a = Vec2 (cos a) (sin a)
+
+sinCos' {- ' CPP is sensitive to primes -} :: Floating a => a -> Normal2 a
+sinCos' = toNormalUnsafe . sinCos
+
+sinCosRadius :: Floating a => a    -- ^ angle (in radians)
+             -> a    -- ^ radius
+             -> Vec2 a
+sinCosRadius a r = Vec2 (r * cos a) (r * sin a)
+
+-- | The angle relative to the positive X axis
+angle2 :: RealFloat a => Vec2 a -> a
+angle2 (Vec2 x y) = atan2 y x
+
+angle2' {- ' CPP is sensitive to primes -} :: RealFloat a => Normal2 a -> a
+angle2' = angle2 . fromNormal
+
+-- | Rotation matrix by a given angle (in radians), counterclockwise.
+rotMatrix2 :: Floating a => a -> Mat2 a
+rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a
+
+rotMatrixOrtho2 :: Floating a => a -> Ortho2 a
+rotMatrixOrtho2 = toOrthoUnsafe . rotMatrix2
+
+rotate2 :: Floating a => a -> Vec2 a -> Vec2 a
+rotate2 a v = v .* (rotMatrix2 a) 
+
+-- |Rotates counterclockwise by 90 degrees.
+rotateCCW :: Floating a => Vec2 a -> Vec2 a
+rotateCCW (Vec2 x y) = Vec2 (-y) x
+
+-- |Rotates clockwise by 90 degrees.
+rotateCW :: Floating a => Vec2 a -> Vec2 a
+rotateCW (Vec2 x y) = Vec2 y (-x)
+
+
diff --git a/src/Data/Vect/Floating/Util/Dim3.hs b/src/Data/Vect/Floating/Util/Dim3.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/Util/Dim3.hs
@@ -0,0 +1,151 @@
+module Data.Vect.Floating.Util.Dim3 where
+
+import Data.Vect.Floating.Base
+
+--------------------------------------------------------------------------------
+
+-- | Example: @structVec3 [1,2,3,4,5,6] = [ Vec3 1 2 3 , Vec3 4 5 6]@.
+structVec3 :: [a] -> [Vec3 a]
+structVec3 [] = []
+structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) 
+structVec3 _ = error "structVec3"
+
+-- | The opposite of "structVec3".
+destructVec3 :: [Vec3 a] -> [a]
+destructVec3 [] = []
+destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)  
+
+--------------------------------------------------------------------------------
+
+det3 :: Floating a => Vec3 a -> Vec3 a -> Vec3 a -> a
+det3 u v w = det (u,v,w)
+
+--------------------------------------------------------------------------------
+
+translate3X :: Num a => a -> Vec3 a -> Vec3 a
+translate3Y :: Num a => a -> Vec3 a -> Vec3 a
+translate3Z :: Num a => a -> Vec3 a -> Vec3 a
+
+translate3X t (Vec3 x y z) = Vec3 (x+t) y z 
+translate3Y t (Vec3 x y z) = Vec3 x (y+t) z 
+translate3Z t (Vec3 x y z) = Vec3 x y (z+t) 
+
+vec3X :: Num a => Vec3 a
+vec3Y :: Num a => Vec3 a
+vec3Z :: Num a => Vec3 a
+
+vec3X = Vec3 1 0 0
+vec3Y = Vec3 0 1 0
+vec3Z = Vec3 0 0 1
+
+rotMatrixZ :: Floating a => a -> Mat3 a
+rotMatrixY :: Floating a => a -> Mat3 a
+rotMatrixX :: Floating a => a -> Mat3 a
+
+-- These are intended for multiplication on the /right/.
+-- Should be consistent with the rotation around an arbitrary axis 
+-- (eg, @rotMatrixY a == rotate3 a vec3Y@)
+rotMatrixZ a = Mat3 (Vec3 c s 0) (Vec3 (-s) c 0) (Vec3 0 0 1) where c = cos a; s = sin a
+rotMatrixY a = Mat3 (Vec3 c 0 (-s)) (Vec3 0 1 0) (Vec3 s 0 c) where c = cos a; s = sin a
+rotMatrixX a = Mat3 (Vec3 1 0 0) (Vec3 0 c s) (Vec3 0 (-s) c) where c = cos a; s = sin a
+
+--------------------------------------------------------------------------------
+
+rotate3' {- ' CPP is sensitive to primes -} 
+  :: Floating a => a       -- ^ angle (in radians)
+  -> Normal3 a   -- ^ axis (should be a /unit/ vector!) 
+  -> Vec3 a      -- ^ vector
+  -> Vec3 a      -- ^ result
+rotate3' angle axis v = v .* (rotMatrix3' axis angle)
+
+rotate3 
+  :: Floating a => a    -- ^ angle (in radians)
+  -> Vec3 a   -- ^ axis (arbitrary nonzero vector)
+  -> Vec3 a   -- ^ vector
+  -> Vec3 a  -- ^ result
+rotate3 angle axis v = v .* (rotMatrix3 axis angle)
+      
+-- | Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
+rotMatrix3 :: Floating a => Vec3 a -> a -> Mat3 a
+rotMatrix3 v a = rotMatrix3' (mkNormal v) a
+
+rotMatrixOrtho3 :: Floating a => Vec3 a -> a -> Ortho3 a
+rotMatrixOrtho3 v a = toOrthoUnsafe $ rotMatrix3 v a
+
+-- | Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
+rotMatrix3' :: {- ' CPP is sensitive to primes -} Floating a => Normal3 a -> a -> Mat3 a
+rotMatrix3' u a = 
+  let v = fromNormal u
+      c = cos a
+      s = sin a
+      m1 = scalarMul (1-c) (outer v v)
+      x = _1 v
+      y = _2 v
+      z = _3 v
+      m2 = Mat3 (Vec3   c    ( s*z) (-s*y))
+                (Vec3 (-s*z)   c    ( s*x))
+                (Vec3 ( s*y) (-s*x)   c   )
+  in (m1 &+ m2)
+
+rotMatrixOrtho3' :: {- ' CPP is sensitive to primes -} Floating a => Normal3 a -> a -> Ortho3 a
+rotMatrixOrtho3' u a = toOrthoUnsafe $ rotMatrix3' u a
+
+--------------------------------------------------------------------------------
+
+-- | Reflects a vector to an axis: that is, the result of @reflect n v@ is
+-- 2\<n,v\>n - v
+reflect :: Floating a => Normal3 a -> Vec3 a -> Vec3 a
+reflect u v = (s *& n) &- v where 
+  n = fromNormal u
+  s = 2 * (n &. v)
+
+reflect' :: Floating a => Normal3 a -> Normal3 a -> Normal3 a
+reflect' u x = toNormalUnsafe $ reflect u (fromNormal x)
+  
+refract :: (Floating a, Ord a) => a -> Normal3 a -> Vec3 a -> Vec3 a
+refract eta u v = s *& fromNormal w where
+  s = norm v 
+  w = refract' eta u (toNormalUnsafe $ v &* (1.0/s))
+  
+-- | Refraction.
+-- First parameter (@eta@) is the relative refraction index 
+--
+-- >        refl_inside
+-- > eta = --------------
+-- >        refl_outside
+--
+-- where \"inside\" is the direction of the second argument 
+-- (to vector normal to plane which models the boundary 
+-- between the two materials). That is, total internal reflection
+-- can occur when @eta>1@.
+--
+-- The convention is that the origin is the point of intersection
+-- of the ray and the surface, and all the vectors \"point away\"
+-- from here (unlike, say, GLSL's @refract@, where the incident
+-- vector \"points towards\" the material)
+refract' {- ' CPP is sensitive to primes -} 
+  :: (Floating a, Ord a) => a -> Normal3 a -> Normal3 a -> Normal3 a
+refract' eta u i = 
+  if k<0
+    then reflect' u i 
+    else toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n 
+  where
+    n = fromNormal u
+    v = fromNormal i
+    c = n &. v
+    k = 1 - eta*eta*(1-c*c)
+
+-- | When total internal reflection would occur, we return "Nothing".
+refractOnly' {- ' CPP is sensitive to primes -}  
+  :: (Floating a, Ord a) => a -> Normal3 a -> Normal3 a -> Maybe (Normal3 a)
+refractOnly' eta u i = 
+  if k<0
+    then Nothing 
+    else Just $ toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n 
+  where
+    n = fromNormal u
+    v = fromNormal i
+    c = n &. v
+    k = 1 - eta*eta*(1-c*c)
+
+--------------------------------------------------------------------------------
diff --git a/src/Data/Vect/Floating/Util/Dim4.hs b/src/Data/Vect/Floating/Util/Dim4.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/Util/Dim4.hs
@@ -0,0 +1,92 @@
+-- | Rotation around an arbitrary plane in four dimensions, and other miscellanea.
+-- Not very useful for most people, and not re-exported by "Data.Vect".
+
+module Data.Vect.Floating.Util.Dim4 where
+
+import Data.Vect.Floating.Base
+import Data.Vect.Floating.GramSchmidt
+
+structVec4 :: [a] -> [Vec4 a]
+structVec4 [] = []
+structVec4 (x:y:z:w:ls) = (Vec4 x y z w):(structVec4 ls) 
+structVec4 _ = error "structVec4"
+
+destructVec4 :: [Vec4 a] -> [a]
+destructVec4 [] = []
+destructVec4 ((Vec4 x y z w):ls) = x:y:z:w:(destructVec4 ls)  
+
+--det4 :: Vec4 -> Vec4 -> Vec4 -> Vec4 -> Flt
+--det4 u v w z = det (u,v,w,z)
+
+translate4X :: Num a => a -> Vec4 a -> Vec4 a
+translate4Y :: Num a => a -> Vec4 a -> Vec4 a
+translate4Z :: Num a => a -> Vec4 a -> Vec4 a
+translate4W :: Num a => a -> Vec4 a -> Vec4 a
+
+translate4X t (Vec4 x y z w) = Vec4 (x+t) y z w 
+translate4Y t (Vec4 x y z w) = Vec4 x (y+t) z w 
+translate4Z t (Vec4 x y z w) = Vec4 x y (z+t) w
+translate4W t (Vec4 x y z w) = Vec4 x y z (w+t) 
+
+vec4X :: Num a => Vec4 a
+vec4Y :: Num a => Vec4 a
+vec4Z :: Num a => Vec4 a
+vec4W :: Num a => Vec4 a
+
+vec4X = Vec4 1 0 0 0
+vec4Y = Vec4 0 1 0 0
+vec4Z = Vec4 0 0 1 0
+vec4W = Vec4 0 0 0 1
+
+--------------------------------------------------------------------------------
+
+-- |If @(x,y,u,v)@ is an orthonormal system, then (written in pseudo-code)
+-- @biVector4 (x,y) = plusMinus (reverse $ biVector4 (u,v))@.
+-- This is a helper function for the 4 dimensional rotation code.
+-- If @(x,y,z,p,q,r) = biVector4 a b@, then the corresponding antisymmetric tensor is
+--
+-- > [  0  r  q  p ]
+-- > [ -r  0  z -y ]
+-- > [ -q -z  0  x ]
+-- > [ -p  y -x  0 ]
+biVector4 :: Num a => Vec4 a -> Vec4 a -> (a,a,a,a,a,a)
+biVector4 (Vec4 x y z w) (Vec4 a b c d) = 
+  ( x*b-y*a , x*c-z*a , x*d-w*a , y*c-z*b , -y*d+w*b , z*d-w*c )
+
+-- | the corresponding antisymmetric tensor
+biVector4AsTensor :: Num a => Vec4 a -> Vec4 a -> Mat4 a
+biVector4AsTensor v w = 
+  Mat4 ( Vec4   0  ( r) ( q) ( p) )
+       ( Vec4 (-r)   0  ( z) (-y) )
+       ( Vec4 (-q) (-z)   0  ( x) )
+       ( Vec4 (-p) ( y) (-x)   0  )
+  where 
+    (x,y,z,p,q,r) = biVector4 v w
+
+-- | We assume that the axes are normalized and /orthogonal/ to each other!
+rotate4' :: {- ' CPP is sensitive to primes -} Floating a => a -> (Normal4 a,Normal4 a) -> Vec4 a -> Vec4 a
+rotate4' angle axes v = v .* (rotMatrix4' angle axes)
+
+-- | We assume only that the axes are independent vectors.
+rotate4 :: Floating a => a -> (Vec4 a,Vec4 a) -> Vec4 a -> Vec4 a
+rotate4 angle axes v = v .* (rotMatrix4 angle axes)
+
+-- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors.
+-- Intended for multiplication on the /right/!
+rotMatrix4' :: {- ' CPP is sensitive to primes -} Floating a => a -> (Normal4 a,Normal4 a) -> Mat4 a
+rotMatrix4' angle (u1,u2) = m1 &+ (s *& m2) &+ m3 
+  where
+    v = fromNormal u1 ; w = fromNormal u2
+    c = cos angle ; s = sin angle
+    m1 = scalarMul (1-c) ( outer v v  &+  outer w w )
+    m2 = biVector4AsTensor v w
+    m3 = diag (Vec4 c c c c)
+
+-- | We assume only that the axes are independent vectors.
+rotMatrix4 :: Floating a => a -> (Vec4 a,Vec4 a) -> Mat4 a  
+rotMatrix4 angle axes = 
+  rotMatrix4' angle $ liftPair toNormalUnsafe $ gramSchmidtNormalize axes 
+  where 
+    liftPair f (x,y) = (f x, f y)
+    
+    
diff --git a/src/Data/Vect/Floating/Util/Projective.hs b/src/Data/Vect/Floating/Util/Projective.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/Util/Projective.hs
@@ -0,0 +1,95 @@
+{-# LANGUAGE FlexibleContexts, ScopedTypeVariables #-}
+-- | Classic 4x4 projective matrices, encoding the affine transformations of R^3.
+-- Our convention is that they are intended for multiplication on
+-- the /right/, that is, they are of the form
+--
+-- >     _____
+-- > [  |     |  0  ]
+-- > [  | 3x3 |  0  ]
+-- > [  |_____|  0  ]
+-- > [  p  q  r  1  ]
+--
+-- Please note that by default, OpenGL stores the matrices (in memory) by columns, while we 
+-- store them by rows; but OpenGL also use the opposite convention (so the OpenGL projective matrices 
+-- are intended for multiplication on the /left/). So in effect, they are the same when stored in the memory,
+-- say with @poke :: Ptr Mat4 -> Mat4 -> IO ()@.
+--
+-- Warning: The naming conventions will probably change in the future.
+
+module Data.Vect.Floating.Util.Projective where
+
+import Data.Vect.Floating.Base
+import Data.Vect.Floating.Util.Dim3
+
+import qualified Data.Vect.Floating.Util.Dim4 as Dim4
+
+--------------------------------------------------------------------------------
+
+rotMatrixProj4' :: {- ' CPP is sensitive to primes -} (Floating a, Projective a Vec3 Mat3 Ortho3 b Proj4) => a -> Normal3 a -> Proj4 a
+rotMatrixProj4' angle axis = linear $ rotMatrix3' axis angle
+
+rotMatrixProj4 :: Floating a => a -> Vec3 a -> Proj4 a
+rotMatrixProj4 angle axis = linear $ rotMatrix3 axis angle
+
+-- | synonym for "rotateAfterProj4"
+rotateProj4 :: Floating a => a -> Normal3 a -> Proj4 a -> Proj4 a
+rotateProj4 = rotateAfterProj4
+
+-- | Synonym for @\m -> m .*. rotMatrixProj4 angle axis@.
+rotateAfterProj4 :: Floating a => a -> Normal3 a -> Proj4 a -> Proj4 a
+rotateAfterProj4 angle axis m = m .*. (rotMatrixProj4' angle axis) 
+
+-- | Synonym for @\m -> rotMatrixProj4 angle axis .*. m@.
+rotateBeforeProj4 :: Floating a => a -> Normal3 a -> Proj4 a -> Proj4 a
+rotateBeforeProj4 angle axis m = (rotMatrixProj4' angle axis) .*. m 
+
+---------------
+
+--scalingUniformProj3 :: Flt -> Proj3
+--scalingUniformProj3 x = scaling (Vec2 x x)
+
+scalingUniformProj4 :: Floating a => a -> Proj4 a
+scalingUniformProj4 x = scaling (Vec3 x x x)
+
+-- | Equivalent to @\m -> scaling v .*. m@.
+scaleBeforeProj4 :: Floating a => Vec3 a -> Proj4 a -> Proj4 a
+scaleBeforeProj4 (Vec3 u v w) p4 =  
+  toProjectiveUnsafe $ 
+    Mat4 (u*&a) (v*&b) (w*&c) t
+  where
+    Mat4 a b c t = fromProjective p4
+
+-- | Equivalent to @\m -> m .*. scaling v@.
+scaleAfterProj4 :: Floating a => Vec3 a -> Proj4 a -> Proj4 a
+scaleAfterProj4 v p4 =
+  toProjectiveUnsafe $ 
+    Mat4 (a&!w) (b&!w) (c&!w) (t&!w)
+  where
+    w = extendWith 1 v
+    Mat4 a b c t = fromProjective p4
+    
+---------------
+
+-- | Synonym for "translateAfter4"
+translate4 :: Floating a => Vec3 a -> Proj4 a -> Proj4 a
+translate4 = translateAfter4
+
+-- | Equivalent to @\m -> m .*. translation v@.
+translateAfter4 :: Floating a => Vec3 a -> Proj4 a -> Proj4 a
+translateAfter4 v p4 = 
+  toProjectiveUnsafe $
+    Mat4 r1 r2 r3 (extendWith 0 v &+ r4)
+  where
+    Mat4 r1 r2 r3 r4 = fromProjective p4 
+
+-- | Equivalent to @\m -> translation v .*. m@.
+translateBefore4 :: (Floating a, Extend a Mat3 Mat4, RightModule (Vec3 a) (Mat3 a)) => Vec3 a -> Proj4 a -> Proj4 a
+translateBefore4 v p4 = 
+  toProjectiveUnsafe $ 
+    Mat4 r1 r2 r3 (extendWith 0 u &+ r4) 
+  where 
+   u = v .*  trim mat 
+   mat@(Mat4 r1 r2 r3 r4) = fromProjective p4
+   
+---------------
+
diff --git a/src/Data/Vect/Floating/Util/Quaternion.hs b/src/Data/Vect/Floating/Util/Quaternion.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vect/Floating/Util/Quaternion.hs
@@ -0,0 +1,269 @@
+-- | The unit sphere in the space of quaternions has the group structure
+-- SU(2) coming from the quaternion multiplication, which is the double
+-- cover of the group SO(3) of rotations in R^3. Thus, unit quaternions can
+-- be used to encode rotations in 3D, which is a more compact encoding 
+-- (4 floats) than a 3x3 matrix; however, there are /two/ quaternions
+-- corresponding to each rotation.
+--
+-- See <http://en.wikipedia.org/wiki/Quaternion> and 
+-- <http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation>
+-- for more information.
+
+{-# LANGUAGE DeriveDataTypeable, FlexibleInstances, MultiParamTypeClasses, StandaloneDeriving, GeneralizedNewtypeDeriving #-}
+module Data.Vect.Floating.Util.Quaternion where
+
+--------------------------------------------------------------------------------
+
+import Data.Typeable
+
+import Data.Vect.Floating.Base
+import Data.Vect.Floating.Interpolate
+
+import Foreign.Storable
+import System.Random
+
+--------------------------------------------------------------------------------
+-- * types
+
+-- | The type for quaternions. 
+newtype Quaternion a = Q (Vec4 a)
+  deriving (Read,Show,Storable,AbelianGroup,Random,Interpolate a, Typeable)
+           
+deriving instance Floating a => Vector a Quaternion
+deriving instance Floating a => DotProd a Quaternion
+
+-- | The type for unit quaternions. 
+newtype UnitQuaternion a = U (Vec4 a)
+  deriving (Read,Show,Storable,Typeable)
+           
+deriving instance Floating a => DotProd a UnitQuaternion
+
+-- | An abbreviated type synonym for quaternions
+type Q = Quaternion
+
+-- | An abbreviated type synonym for unit quaternions
+type U = UnitQuaternion
+   
+--------------------------------------------------------------------------------
+
+instance Floating a => UnitVector a Quaternion UnitQuaternion where
+  mkNormal (Q v) = U (normalize v)
+  toNormalUnsafe (Q v) = U v
+  fromNormal (U v) = Q v
+  fromNormalRadius r (U v) = Q (v &* r)
+
+--------------------------------------------------------------------------------
+-- * general quaternions
+
+unitQ :: Num a => Q a
+unitQ = Q (Vec4 1 0 0 0)
+
+zeroQ :: Num a => Q a
+zeroQ = Q (Vec4 0 0 0 0)
+
+multQ :: Num a => Q a -> Q a -> Q a
+multQ (Q (Vec4 a1 b1 c1 d1)) (Q (Vec4 a2 b2 c2 d2)) = Q $ Vec4 
+  (a1*a2 - b1*b2 - c1*c2 - d1*d2)
+  (a1*b2 + b1*a2 + c1*d2 - d1*c2)
+  (a1*c2 - b1*d2 + c1*a2 + d1*b2)
+  (a1*d2 + b1*c2 - c1*b2 + d1*a2)
+
+negQ :: Floating a => Q a -> Q a
+negQ (Q v) = Q (neg v)
+
+normalizeQ :: Floating a => Q a -> Q a
+normalizeQ (Q v) = Q (normalize v)
+
+-- | The inverse quaternion
+invQ :: Floating a => Q a -> Q a
+invQ (Q (Vec4 a b c d)) = Q (v &* (1 / normsqr v)) where 
+  v = Vec4 a (-b) (-c) (-d)
+  
+fromQ :: Q a -> Vec4 a
+fromQ (Q v) = v 
+
+toQ :: Vec4 a -> Q a
+toQ = Q 
+  
+--------------------------------------------------------------------------------
+
+{- 
+-- we use newtype deriving instead
+
+instance AbelianGroup Quaternion where
+  (Q v1) &+ (Q v2) = Q (v1 &+ v2)   
+  (Q v1) &- (Q v2) = Q (v1 &+ v2)   
+  neg (Q v) = Q (neg v)
+  zero = zeroQ
+
+instance DotProd Quaternion where
+  dotprod (Q v1) (Q v2) = dotprod v1 v2
+
+-}
+
+instance Num a => MultSemiGroup (Quaternion a) where
+  one   = unitQ  
+  (.*.) = multQ
+
+
+--------------------------------------------------------------------------------
+-- * unit quaternions
+
+unitU :: Num a => U a
+unitU = U (Vec4 1 0 0 0)
+
+multU :: Num a => U a -> U a -> U a
+multU (U (Vec4 a1 b1 c1 d1)) (U (Vec4 a2 b2 c2 d2)) = U $ Vec4 
+  (a1*a2 - b1*b2 - c1*c2 - d1*d2)
+  (a1*b2 + b1*a2 + c1*d2 - d1*c2)
+  (a1*c2 - b1*d2 + c1*a2 + d1*b2)
+  (a1*d2 + b1*c2 - c1*b2 + d1*a2)
+  
+-- | The opposite quaternion (which encodes the same rotation)
+negU :: Floating a => U a -> U a
+negU (U v) = U (neg v)
+
+-- | This is no-op, up to numerical imprecision.
+-- However, if you multiply together a large number of unit quaternions, 
+-- it may be a good idea to normalize the end result.
+normalizeU :: Floating a => U a -> U a
+normalizeU (U v) = U (normalize v)
+
+-- | The inverse of a unit quaternion
+invU :: Num a => U a -> U a
+invU (U (Vec4 a b c d)) = U $ Vec4 a (-b) (-c) (-d)
+
+--------------------------------------------------------------------------------
+ 
+fromU :: Num a => U a -> Vec4 a
+fromU (U v) = v 
+
+fromU' :: Floating a => U a -> Normal4 a
+fromU' (U v) = toNormalUnsafe v
+
+mkU :: Floating a => Vec4 a -> U a
+mkU = U . normalize
+
+toU :: Floating a => Normal4 a -> U a
+toU = U . fromNormal
+
+unsafeToU :: Num a => Vec4 a -> U a
+unsafeToU = U 
+
+--------------------------------------------------------------------------------
+
+{-  
+-- we use newtype deriving instead
+
+instance DotProd UnitQuaternion where
+  dotprod (Q v1) (Q v2) = dotprod v1 v2
+  
+-}
+   
+instance Floating a => MultSemiGroup (UnitQuaternion a) where
+  one   = unitU  
+  (.*.) = multU
+ 
+instance Num a => LeftModule (UnitQuaternion a) (Vec3 a) where
+  lmul u v = actU u v
+
+instance (Floating a, Ord a, Random a) => Random (UnitQuaternion a) where
+  random g = let (n, h) = random g 
+                 v = fromNormal n
+             in  (U v, h)  
+  randomR _ = random     
+  
+--------------------------------------------------------------------------------
+-- * unit quaternions as rotations
+
+-- | The /left/ action of unit quaternions on 3D vectors.
+-- That is,
+-- 
+-- > actU q1 $ actU q2 v  ==  actU (q1 `multU` q2) v 
+actU :: Num a => U a -> Vec3 a -> Vec3 a
+actU (U (Vec4 a b c d)) (Vec3 x y z) = Vec3 x' y' z' where
+  x' =  x*(aa + bb - cc - dd)  +  y*(  2 * (bc - ad)  )  +  z*(  2 * (bd + ac)  )
+  y' =  x*(  2 * (bc + ad)  )  +  y*(aa - bb + cc - dd)  +  z*(  2 * (cd - ab)  )
+  z' =  x*(  2 * (bd - ac)  )  +  y*(  2 * (cd + ab)  )  +  z*(aa - bb - cc + dd)
+  --
+  aa = a*a ; bb = b*b ; cc = c*c ; dd = d*d
+  ab = a*b ; ac = a*c ; ad = a*d
+  bc = b*c ; bd = b*d ; cd = c*d
+  
+actU' :: Floating a => U a -> Normal3 a -> Normal3 a
+actU' u n = toNormalUnsafe $ actU u (fromNormal n)
+
+-- | The quaternion to encode rotation around an axis. Please note
+-- that quaternions act on the /left/, that is
+--
+-- > rotU axis1 angl1 *. rotU axis2 angl2 *. v  ==  (rotU axis1 angl1 .*. rotU axis2 angl2) *. v 
+--
+rotU :: Floating a => Vec3 a -> a -> U a
+rotU axis angle = rotU' (mkNormal axis) angle
+
+rotU' {- ' CPP is sensitive to primes -} :: Floating a => Normal3 a -> a -> U a
+rotU' axis angle = U (Vec4 c (x*s) (y*s) (z*s)) where
+  Vec3 x y z = fromNormal axis 
+  half = 0.5 * angle
+  c = cos half
+  s = sin half
+
+-- | Interpolation of unit quaternions. Note that when applied to rotations,
+-- this may be not what you want, since it is possible that the shortest path
+-- in the space of unit quaternions is not the shortest path in the space of
+-- rotations; see 'slerpU'!
+longSlerpU :: Floating a => a -> U a -> U a -> U a
+longSlerpU t (U p0) (U p1) = U v where
+  v = (p0 &* y0) &+ (p1 &* y1) 
+  omega = acos (p0 &. p1)
+  s = sin omega
+  y0 = sin (omega*(1-t)) / s 
+  y1 = sin (omega*   t ) / s
+
+-- | This is shortest path interpolation in the space of rotations; however
+-- this is achieved by possibly flipping the first endpoint in the space of
+-- quaternions. Thus @slerpU 0.001 q1 q2@ may be very far from @q1@ (and very
+-- close to @negU q1@) in the space of quaternions (but they are very close
+-- in the space of rotations). 
+slerpU :: (Floating a, Ord a) => a -> U a -> U a -> U a
+slerpU t (U p0') (U p1) = U v where
+  v = (p0 &* y0) &+ (p1 &* y1) 
+  
+  d' = p0' &. p1  
+  (d,p0) = if d' >= 0 
+    then ( d',     p0')
+    else (-d', neg p0')
+    
+  omega = acos d
+  s = sin omega
+  y0 = sin (omega*(1-t)) / s 
+  y1 = sin (omega*   t ) / s
+  
+-- | Makes a rotation matrix (to be multiplied with on the /right/) out of a unit quaternion:
+--
+-- > v .* rightOrthoU (rotU axis angl)  ==  v .* rotMatrix3 axis angl
+-- 
+-- Please note that while these matrices act on the /right/, quaternions act on the /left/; thus
+-- 
+-- > rightOrthoU q1 .*. rightOrthoU q2  ==  rightOrthoU (q2 .*. q1)
+--
+rightOrthoU :: Floating a => U a -> Ortho3 a
+rightOrthoU = toOrthoUnsafe . transpose . fromOrtho . leftOrthoU
+
+-- | Makes a rotation matrix (to be multiplied with on the /left/) out of a unit quaternion.
+--
+-- > leftOrthoU (rotU axis angl) *. v  ==  v .* rotMatrix3 axis angl
+-- 
+leftOrthoU :: Floating a => U a -> Ortho3 a
+leftOrthoU (U (Vec4 a b c d)) = toOrthoUnsafe $ Mat3 row1 row2 row3 where
+  row1 = Vec3  (aa + bb - cc - dd)  (  2 * (bc - ad)  )  (  2 * (bd + ac)  )
+  row2 = Vec3  (  2 * (bc + ad)  )  (aa - bb + cc - dd)  (  2 * (cd - ab)  )
+  row3 = Vec3  (  2 * (bd - ac)  )  (  2 * (cd + ab)  )  (aa - bb - cc + dd)
+  --
+  aa = a*a ; bb = b*b ; cc = c*c ; dd = d*d
+  ab = a*b ; ac = a*c ; ad = a*d
+  bc = b*c ; bd = b*d ; cd = c*d
+  
+--------------------------------------------------------------------------------
+  
+  
diff --git a/vect-floating.cabal b/vect-floating.cabal
new file mode 100644
--- /dev/null
+++ b/vect-floating.cabal
@@ -0,0 +1,27 @@
+-- Initial vect-floating.cabal generated by cabal init.  For further 
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name:                vect-floating
+version:             0.1.0.0
+synopsis:            A low-dimensional linear algebra library, operating on the Floating typeclass
+description:         A low-dimensional (2, 3, 4 and Quaternions) linear algebra library, 
+		     with lots of useful functions. Forked from vect library at <http://hackage.haskell.org/package/vect>
+		     by Balazs Komuves in order to use Floating typeclass instead of Float, Double. 
+		     Will be compatible with accelerate in a future package. 
+		     Many thanks to Balazs Komuves.
+homepage:            http://github.com/cdurham/vect-floating
+license:             BSD3
+license-file:        LICENSE
+author:              Balazs Komuves, Charles Durham
+maintainer:          Charles Durham <cpdurham@gmail.com>
+copyright:           (c) 2008-2011 Balazs Komuves (c) 2014 Charles Durham
+category:            Math
+build-type:          Simple
+cabal-version:       >=1.8
+
+library
+  exposed-modules: Data.Vect.Floating, Data.Vect.Floating.Base, Data.Vect.Floating.Interpolate, Data.Vect.Floating.GramSchmidt, Data.Vect.Floating.Util.Projective, Data.Vect.Floating.Util.Dim2, Data.Vect.Floating.Util.Dim3, Data.Vect.Floating.Util.Dim4, Data.Vect.Floating.Util.Quaternion
+		       
+  -- other-modules:       
+  build-depends:       base ==4.6.*, random ==1.0.*
+  hs-source-dirs:      src