diff --git a/Data/Vect/Double/Base.hs b/Data/Vect/Double/Base.hs
--- a/Data/Vect/Double/Base.hs
+++ b/Data/Vect/Double/Base.hs
@@ -1,11 +1,34 @@
 {-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}
 
-module Data.Vect.Flt.Base where
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, GeneralizedNewtypeDeriving #-}
 
+module Data.Vect.Flt.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 Control.Monad
 import System.Random  
 import Foreign
 
+--------------------------------------------------------------------------------
 -- class declarations
 
 class AbelianGroup g where
@@ -20,15 +43,17 @@
 vecSum :: AbelianGroup g => [g] -> g
 vecSum l = foldl (&+) zero l 
 
-class (AbelianGroup r) => 
-      Ring r where
+class MultSemiGroup r where
   (.*.) :: r -> r -> r
   one   :: r
 
+class (AbelianGroup r, MultSemiGroup r) => Ring r 
+
 infixl 7 .*. 
 
-ringProduct :: Ring r => [r] -> r
-ringProduct l = foldl (.*.) one l
+-- 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
@@ -75,7 +100,7 @@
   lensqr = normsqr
   dotprod :: v -> v -> Flt
   normsqr v = (v &. v)  
-  norm = sqrt.lensqr
+  norm = sqrt . lensqr
   dotprod = (&.)
 
 infix 7 &.
@@ -112,26 +137,28 @@
   fromNormalRadius :: Flt -> u -> v
   fromNormalRadius t n = t *& fromNormal n 
 
--- | projects the first vector onto the direction of the second (unit) vector
+-- | Projects the first vector onto the direction of the second (unit) vector
 project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v
 project' what dir = projectUnsafe what (fromNormal dir)
 
--- | direction (second argument) is assumed to be a /unit/ vector!
+-- | Direction (second argument) is assumed to be a /unit/ vector!
 projectUnsafe :: (Vector v, DotProd v) => v -> v -> v
 projectUnsafe what dir = what &- dir &* (what &. dir)
 
 project :: (Vector v, DotProd v) => v -> v -> v
 project what dir = what &- dir &* ((what &. dir) / (dir &. dir))
 
--- | since unit vectors are not a group, we need a separate function.
+-- | Since unit vectors are not a group, we need a separate function.
 flipNormal :: UnitVector v n => n -> n 
 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
@@ -166,13 +193,50 @@
 "inverse is an involution"    forall m. inverse (inverse m) = m
   #-}
   
+class Matrix m => Orthogonal m o | m->o, o->m where  
+  fromOrtho     :: o -> m 
+  toOrthoUnsafe :: m -> o
+  
+class (AbelianGroup m, Matrix m) => MatrixNorms m where
+  frobeniusNorm  :: m -> Flt       -- ^ the frobenius norm (= euclidean norm in the space of matrices)
+  matrixDistance :: m -> m -> Flt  -- ^ euclidean distance in the space of matrices
+  operatorNorm   :: m -> Flt       -- ^ (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 m where
   det :: m -> Flt    
-    
+
+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 v, UnitVector v u, Matrix m, Vector m, Tensor m v) => u -> m
+householder u = idmtx &- (2 *& outer v v) 
+  where v = fromNormal u
+
+householderOrtho :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v, Orthogonal m o) => u -> o
+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 v, Orthogonal n o, Diagonal v n) => Projective 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 -> m
+  toProjectiveUnsafe :: m -> p
+  orthogonal         :: o -> p
+  linear             :: n -> p
+  translation        :: v -> p
+  scaling            :: v -> p
+
+--------------------------------------------------------------------------------
 -- Vec / Mat datatypes
  
 data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt 
@@ -182,16 +246,16 @@
 data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt 
   deriving (Read,Show)
 
--- | these are /row/ vectors 
+-- | The components are /row/ vectors 
 data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show)
 data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show)
 data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)
 
 -- | The assumption when dealing with these is always that they are of unit length.
 -- Also, interpolation works differently.
-newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) 
-newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) 
-newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) 
+newtype Normal2 = Normal2 Vec2 deriving (Read,Show,Storable,DotProd,Dimension) 
+newtype Normal3 = Normal3 Vec3 deriving (Read,Show,Storable,DotProd,Dimension) 
+newtype Normal4 = Normal4 Vec4 deriving (Read,Show,Storable,DotProd,Dimension) 
 
 mkVec2 :: (Flt,Flt) -> Vec2
 mkVec3 :: (Flt,Flt,Flt) -> Vec3
@@ -201,6 +265,19 @@
 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 = Ortho2 Mat2 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+newtype Ortho3 = Ortho3 Mat3 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+newtype Ortho4 = Ortho4 Mat4 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+
+-- | Projective matrices, encoding affine transformations in dimension one less.
+newtype Proj3 = Proj3 Mat3 deriving (Read,Show,Storable,MultSemiGroup)
+newtype Proj4 = Proj4 Mat4 deriving (Read,Show,Storable,MultSemiGroup)
+
+--------------------------------------------------------------------------------
 -- Unit vectors
   
 instance UnitVector Vec2 Normal2 where
@@ -218,47 +295,163 @@
   fromNormal (Normal4 v) = v 
   toNormalUnsafe = Normal4
 
-rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)
-rndUnit g = 
+_rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)
+_rndUnit g = 
   if d > 0.01
     then ( v &* (1.0/d) , h )
-    else rndUnit h
+    else _rndUnit h
   where
     (v,h) = random g
     d = norm v
     
 instance Random Normal2 where
-  random g = let (v,h) = rndUnit g in (Normal2 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal2 v, h)  
   randomR _ = random
 
 instance Random Normal3 where
-  random g = let (v,h) = rndUnit g in (Normal3 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal3 v, h)  
   randomR _ = random
 
 instance Random Normal4 where
-  random g = let (v,h) = rndUnit g in (Normal4 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal4 v, h)  
   randomR _ = random
 
-{-
-instance Storable Normal2 where
-  alignment _ = alignment (undefined::Vec2)
-  sizeOf    _ = sizeOf    (undefined::Vec2)
-  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)
-  poke p (Normal2 v) = poke (castPtr p) v  
- 
-instance Storable Normal3 where
-  alignment _ = alignment (undefined::Vec3)
-  sizeOf    _ = sizeOf    (undefined::Vec3)
-  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)
-  poke p (Normal3 v) = poke (castPtr p) v  
+instance CrossProd Normal3 where
+  crossprod (Normal3 v) (Normal3 w) = mkNormal (crossprod v w)
 
-instance Storable Normal4 where
-  alignment _ = alignment (undefined::Vec4)
-  sizeOf    _ = sizeOf    (undefined::Vec4)
-  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)
-  poke p (Normal4 v) = poke (castPtr p) v  
--}
+--------------------------------------------------------------------------------
+-- Orthogonal matrices
 
+instance Orthogonal Mat2 Ortho2 where
+  fromOrtho (Ortho2 o) = o
+  toOrthoUnsafe = Ortho2
+
+instance Orthogonal Mat3 Ortho3 where
+  fromOrtho (Ortho3 o) = o
+  toOrthoUnsafe = Ortho3 
+
+instance Orthogonal Mat4 Ortho4 where
+  fromOrtho (Ortho4 o) = o
+  toOrthoUnsafe = Ortho4
+
+------
+
+instance Matrix Ortho2 where
+  transpose (Ortho2 o) = Ortho2 (transpose o)
+  idmtx = Ortho2 idmtx
+  inverse = transpose
+
+instance Matrix Ortho3 where
+  transpose (Ortho3 o) = Ortho3 (transpose o)
+  idmtx = Ortho3 idmtx
+  inverse = transpose
+
+instance Matrix Ortho4 where
+  transpose (Ortho4 o) = Ortho4 (transpose o)
+  idmtx = Ortho4 idmtx
+  inverse = transpose
+
+------
+
+instance Random Ortho2 where
+  random g = let (o,h) = _rndOrtho2 g in (toOrthoUnsafe (_flip1stRow2 o), h)
+  randomR _ = random
+
+instance Random Ortho3 where
+  random g = let (o,h) = _rndOrtho3 g in (toOrthoUnsafe (             o), h)
+  randomR _ = random
+
+instance Random Ortho4 where
+  random g = let (o,h) = _rndOrtho4 g in (toOrthoUnsafe (_flip1stRow4 o), h)
+  randomR _ = random
+
+------
+
+-- determinant will be -1
+_rndOrtho2 :: RandomGen g => g -> (Mat2, g)
+_rndOrtho2 g = (h2, g1) where
+  h2 = householder u2 :: Mat2 
+  (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 :: RandomGen g => g -> (Mat3, g) 
+_rndOrtho3 g = ( (h3 .*. m3), g2) where
+  m3 = (extendWith :: Flt -> Mat2 -> Mat3) 1 o2 
+  h3 = householder u3 :: Mat3
+  (u3,g1) = random g
+  (o2,g2) = _rndOrtho2 g1
+
+-- determinant will be -1
+_rndOrtho4 :: RandomGen g => g -> (Mat4, g) 
+_rndOrtho4 g = ( (h4 .*. m4), g2) where
+  m4 = (extendWith :: Flt -> Mat3 -> Mat4) 1 o3 
+  h4 = householder u4 :: Mat4
+  (u4,g1) = random g
+  (o3,g2) = _rndOrtho3 g1
+
+------
+
+_flip1stRow2 :: Mat2 -> Mat2
+_flip1stRow2 (Mat2 a b) = Mat2 (neg a) b
+
+_flip1stRow3 :: Mat3 -> Mat3
+_flip1stRow3 (Mat3 a b c) = Mat3 (neg a) b c
+
+_flip1stRow4 :: Mat4 -> Mat4
+_flip1stRow4 (Mat4 a b c d) = Mat4 (neg a) b c d
+
+--------------------------------------------------------------------------------
+-- projective matrices
+  
+instance Projective Vec2 Mat2 Ortho2 Mat3 Proj3 where
+  fromProjective (Proj3 m) = m
+  toProjectiveUnsafe = Proj3
+  orthogonal = Proj3 . extendWith 1 . fromOrtho
+  linear     = Proj3 . extendWith 1
+  translation v = Proj3 $ Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (extendWith 1 v)
+  scaling     v = Proj3 $ diag (extendWith 1 v)
+  
+instance Projective Vec3 Mat3 Ortho3 Mat4 Proj4 where
+  fromProjective (Proj4 m) = m
+  toProjectiveUnsafe = Proj4
+  orthogonal = Proj4 . extendWith 1 . fromOrtho 
+  linear     = Proj4 . extendWith 1
+  translation v = Proj4 $ Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (extendWith 1 v)
+  scaling     v = Proj4 $ diag (extendWith 1 v)
+
+instance Matrix Proj3 where
+  idmtx = Proj3 idmtx
+  transpose (Proj3 m) = Proj3 (transpose m)
+  inverse = _invertProj3
+
+instance Matrix Proj4 where
+  idmtx = Proj4 idmtx
+  transpose (Proj4 m) = Proj4 (transpose m)
+  inverse = _invertProj4
+
+_invertProj3 :: Proj3 -> Proj3
+_invertProj3 (Proj3 mat@(Mat3 _ _ t)) = 
+  Proj3 $ Mat3 (extendZero a) (extendZero b) (extendWith 1 t') 
+  where
+    t' = neg $ (trim t :: Vec2) .* invm2 
+    invm2@(Mat2 a b) = inverse $ (trim mat :: Mat2)
+
+-- Inverts a projective 4x4 matrix. But you can simply use "inverse" instead.
+-- We assume that the bottom-right corner is 1.
+_invertProj4 :: Proj4 -> Proj4
+_invertProj4 (Proj4 mat@(Mat4 _ _ _ t)) = 
+  Proj4 $ Mat4 (extendZero a) (extendZero b) (extendZero c) (extendWith 1 t') 
+  where
+    t' = neg $ (trim t :: Vec3) .* invm3 
+    invm3@(Mat3 a b c) = inverse $ (trim mat :: Mat3)
+
+--------------------------------------------------------------------------------
 -- Vec2 instances
 
 instance HasCoordinates Vec2 Flt where
@@ -314,7 +507,10 @@
         k = sizeOf (undefined::Flt)
     poke        p   x
     pokeByteOff p k y
-               
+
+instance Dimension Vec2 where dim _ = 2
+
+--------------------------------------------------------------------------------                    
 -- Mat2 instances
 
 instance HasCoordinates Mat2 Vec2 where
@@ -336,19 +532,21 @@
   (&+) (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   -- (zero::Vec2) (zero::Vec2)
-
+  zero = Mat2 zero zero  
+  
 instance Vector 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 Ring Mat2 where
+instance MultSemiGroup Mat2 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 Ring Mat2
+
 instance LeftModule Mat2 Vec2 where
   lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) 
   
@@ -362,11 +560,6 @@
   outer (Vec2 a b) (Vec2 x y) = Mat2
     (Vec2 (a*x) (a*y))
     (Vec2 (b*x) (b*y))
-{-
-  outer v w = 
-    let full = Mat2 (Vec2 1 1) (Vec2 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
 instance Determinant Mat2 where
   det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c 
@@ -393,6 +586,24 @@
     poke        p   r1
     pokeByteOff p k r2
 
+instance Random Mat2 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 Dimension Mat2 where dim _ = 2
+     
+instance MatrixNorms Mat2 where 
+  frobeniusNorm (Mat2 r1 r2) =  
+    sqrt $
+      normsqr r1 + 
+      normsqr r2
+     
+--------------------------------------------------------------------------------     
 -- Vec3 instances
 
 instance HasCoordinates Vec3 Flt where
@@ -454,7 +665,10 @@
     poke        p       x
     pokeByteOff p (k  ) y
     pokeByteOff p (k+k) z
-   
+
+instance Dimension Vec3 where dim _ = 3
+
+--------------------------------------------------------------------------------   
 -- Mat3 instances
 
 instance HasCoordinates Mat3 Vec3 where
@@ -495,13 +709,13 @@
   (&+) (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   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)
+  zero = Mat3 zero zero zero 
 
 instance Vector 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 Ring Mat3 where
+instance MultSemiGroup Mat3 where
   (.*.) (Mat3 r1 r2 r3) n = 
     let (Mat3 c1 c2 c3) = transpose n
     in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))
@@ -509,6 +723,8 @@
             (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))
   one = idmtx 
 
+instance Ring Mat3
+
 instance LeftModule Mat3 Vec3 where
   lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)
   
@@ -523,11 +739,6 @@
     (Vec3 (a*x) (a*y) (a*z))
     (Vec3 (b*x) (b*y) (b*z))
     (Vec3 (c*x) (c*y) (c*z))
-{-
-  outer v w = 
-    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
 instance Determinant Mat3 where
   det (Mat3 r1 r2 r3) = det (r1,r2,r3)
@@ -556,6 +767,26 @@
     pokeByteOff p (k  ) r2
     pokeByteOff p (k+k) r3
 
+instance Random Mat3 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 Dimension Mat3 where dim _ = 3
+  
+instance MatrixNorms Mat3 where 
+  frobeniusNorm (Mat3 r1 r2 r3)  = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 
+    
+--------------------------------------------------------------------------------
 -- Vec4 instances
 
 instance HasCoordinates Vec4 Flt where
@@ -615,6 +846,9 @@
     pokeByteOff p (k+k) z
     pokeByteOff p (3*k) w
 
+instance Dimension Vec4 where dim _ = 4
+
+--------------------------------------------------------------------------------
 -- Mat4 instances
 
 instance HasCoordinates Mat4 Vec4 where
@@ -642,7 +876,7 @@
   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 Ring Mat4 where
+instance MultSemiGroup Mat4 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))
@@ -651,6 +885,8 @@
             (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))
   one = idmtx 
 
+instance Ring Mat4
+
 instance LeftModule Mat4 Vec4 where
   lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)
   
@@ -666,14 +902,10 @@
     (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))
-{-
-  outer v w = 
-    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
---instance Determinant Mat4 where
---  det (Mat4 r1 r2 r3 r4)
+instance Determinant Mat4 where
+  det = error "det/Mat4: not implemented yet" 
+  -- det (Mat4 r1 r2 r3 r4) = 
 
 {-
 instance Show Mat4 where
@@ -701,6 +933,28 @@
     pokeByteOff p (k+k) r3
     pokeByteOff p (3*k) r4
 
+instance Random Mat4 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 Dimension Mat4 where dim _ = 4
+   
+instance MatrixNorms Mat4 where 
+  frobeniusNorm (Mat4 r1 r2 r3 r4) = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 + 
+      normsqr r4  
+    
+--------------------------------------------------------------------------------
 -- Extend instances
 
 instance Extend Vec2 Vec3 where
@@ -719,17 +973,20 @@
   trim (Vec4 x y z _)       = Vec3 x y z
 
 instance Extend Mat2 Mat3 where
-  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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 Extend Mat2 Mat4 where
-  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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 Extend Mat3 Mat4 where
-  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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/Data/Vect/Double/GramSchmidt.hs b/Data/Vect/Double/GramSchmidt.hs
--- a/Data/Vect/Double/GramSchmidt.hs
+++ b/Data/Vect/Double/GramSchmidt.hs
@@ -10,7 +10,7 @@
 
 import Data.Vect.Flt.Base
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
 liftPair :: (a -> b) -> (a,a) -> (b,b)
 liftPair f (x,y) = (f x, f y)
@@ -21,7 +21,7 @@
 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
@@ -33,7 +33,7 @@
 "gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a 
   #-}
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
 instance GramSchmidt (Vec2,Vec2) where
   gramSchmidt = gramSchmidtPair
diff --git a/Data/Vect/Double/OpenGL.hs b/Data/Vect/Double/OpenGL.hs
--- a/Data/Vect/Double/OpenGL.hs
+++ b/Data/Vect/Double/OpenGL.hs
@@ -3,21 +3,92 @@
 -- TODO: the pointer versions of these functions should be really implemented 
 -- via the pointer versions of the original opengl functions...
 
--- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.
+-- | OpenGL support, including 'Vertex', 'TexCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.
  
 module Data.Vect.Flt.OpenGL where
 
 import Control.Monad
 import Data.Vect.Flt.Base
+import Data.Vect.Flt.Util.Projective
 import qualified Graphics.Rendering.OpenGL as GL
 
 import Foreign
 
-import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)
+import Graphics.Rendering.OpenGL hiding 
+  ( Normal3 , rotate , translate , scale
+  , matrix , currentMatrix , withMatrix , multMatrix 
+  )
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
-{-# SPECIALISE radianToDegrees :: Float -> Float #-}
+-- | There should be a big warning here about the different conventions, 
+-- hidden transpositions, and all the confusion this will inevitably cause...
+--
+-- As it stands, 
+--
+-- > glRotate t1 axis1 >> glRotate t2 axis2 >> glRotate t3 axis3
+-- 
+-- has the same result as
+--
+-- > multMatrix (rotMatrixProj4 t3 axis3 .*. rotMatrixProj4 t2 axis2 .*. rotMatrixProj4 t1 axis1)
+--
+-- because at the interface of OpenGL and this library there is a transposition
+-- to compensate for the different conventions. (This transposition is implicit
+-- in the code, because the way the matrices are stored in the memory is also
+-- different: OpenGL stores them column-major, and we store them row-major).
+
+class ToOpenGLMatrix m where
+  makeGLMatrix :: m -> IO (GLmatrix Flt)
+
+class FromOpenGLMatrix m where
+  peekGLMatrix :: GLmatrix Flt -> IO m
+  
+setMatrix :: ToOpenGLMatrix m => Maybe MatrixMode -> m -> IO ()
+setMatrix mode m = makeGLMatrix m >>= \x -> GL.matrix mode $= x
+ 
+getMatrix :: FromOpenGLMatrix m => Maybe MatrixMode -> IO m
+getMatrix mode = get (GL.matrix mode) >>= peekGLMatrix
+
+matrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => Maybe MatrixMode -> StateVar m
+matrix mode = makeStateVar (getMatrix mode) (setMatrix mode)
+
+currentMatrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => StateVar m
+currentMatrix = matrix Nothing
+
+multMatrix :: ToOpenGLMatrix m => m -> IO ()
+multMatrix m = makeGLMatrix m >>= GL.multMatrix
+
+instance ToOpenGLMatrix Mat4 where
+  makeGLMatrix m = GL.withNewMatrix GL.ColumnMajor (flip poke m . castPtr) 
+ 
+instance FromOpenGLMatrix Mat4 where
+  -- huh? GL.withMatrix is strange
+  peekGLMatrix x = GL.withMatrix x $ \_ p -> peek (castPtr p)
+  
+instance ToOpenGLMatrix Mat3 where
+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)
+ 
+instance ToOpenGLMatrix Mat2 where
+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)
+
+instance ToOpenGLMatrix Ortho4 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat4)
+
+instance ToOpenGLMatrix Ortho3 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat3)
+
+instance ToOpenGLMatrix Ortho2 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat2)
+
+instance ToOpenGLMatrix Proj4 where
+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat4)
+
+instance ToOpenGLMatrix Proj3 where
+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat3)
+  
+--------------------------------------------------------------------------------
+
+{-# SPECIALISE radianToDegrees :: Float  -> Float  #-}
 {-# SPECIALISE radianToDegrees :: Double -> Double #-}
 radianToDegrees :: RealFrac a => a -> a
 radianToDegrees x = x * 57.295779513082322
@@ -28,20 +99,61 @@
 degreesToRadian x = x * 1.7453292519943295e-2
 
 -- | The angle is in radians. (WARNING: OpenGL uses degrees!)
-rotate :: Flt -> Vec3 -> IO ()
-rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)
+glRotate :: Flt -> Vec3 -> IO ()
+glRotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)
 
-translate :: Vec3 -> IO ()
-translate (Vec3 x y z) = GL.translate (Vector3 x y z)
+glTranslate :: Vec3 -> IO ()
+glTranslate (Vec3 x y z) = GL.translate (Vector3 x y z)
 
-scale3 :: Vec3 -> IO ()
-scale3 (Vec3 x y z) = GL.scale x y z
+glScale3 :: Vec3 -> IO ()
+glScale3 (Vec3 x y z) = GL.scale x y z
 
-scale :: Flt -> IO ()
-scale x = GL.scale x x x
+glScale :: Flt -> IO ()
+glScale x = GL.scale x x x
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
+ 
+-- | \"Orthogonal projecton\" matrix, a la OpenGL 
+-- (the corresponding functionality is removed in OpenGL 3.1)
+orthoMatrix 
+  :: (Flt,Flt)   -- ^ (left,right)
+  -> (Flt,Flt)   -- ^ (bottom,top)
+  -> (Flt,Flt)   -- ^ (near,far)
+  -> Mat4 
+orthoMatrix (l,r) (b,t) (n,f) = Mat4
+  (Vec4 (2/(r-l)) 0 0 0)
+  (Vec4 0 (2/(t-b)) 0 0)
+  (Vec4 0 0 (-2/(f-n)) 0)
+  (Vec4 (-(r+l)/(r-l)) (-(t+b)/(t-b)) (-(f+n)/(f-n)) 1)
+  
+-- | The same as "orthoMatrix", but with a different parametrization.
+orthoMatrix2 {- ' CPP is sensitive to primes -}
+  :: Vec3     -- ^ (left,top,near)
+  -> Vec3     -- ^ (right,bottom,far)
+  -> Mat4 
+orthoMatrix2 (Vec3 l t n) (Vec3 r b f) = orthoMatrix (l,r) (b,t) (n,f)
 
+-- | \"Perspective projecton\" matrix, a la OpenGL 
+-- (the corresponding functionality is removed in OpenGL 3.1).
+frustumMatrix
+  :: (Flt,Flt)   -- ^ (left,right)
+  -> (Flt,Flt)   -- ^ (bottom,top)
+  -> (Flt,Flt)   -- ^ (near,far)
+  -> Mat4 
+frustumMatrix (l,r) (b,t) (n,f) = Mat4
+  (Vec4 (2*n/(r-l)) 0 0 0)
+  (Vec4 0 (2*n/(t-b)) 0 0)
+  (Vec4 ((r+l)/(r-l)) ((t+b)/(t-b)) (-(f+n)/(f-n)) (-1))
+  (Vec4 0 0 (-2*f*n*(f-n)) 0)
+  
+-- | The same as "frustumMatrix", but with a different parametrization.
+frustumMatrix2 {- ' CPP is sensitive to primes -}
+  :: Vec3     -- ^ (left,top,near)
+  -> Vec3     -- ^ (right,bottom,far)
+  -> Mat4 
+frustumMatrix2 (Vec3 l t n) (Vec3 r b f) = frustumMatrix (l,r) (b,t) (n,f)
+
+--------------------------------------------------------------------------------
 -- Vertex instances
 
 instance GL.Vertex Vec2 where
@@ -56,17 +168,20 @@
   vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)
   vertexv p = peek p >>= vertex   
 
--------------------------------------------------------
-
+--------------------------------------------------------------------------------
 -- the Normal instance
 -- note that there is no Normal2\/Normal4 in the OpenGL binding
 
 instance GL.Normal Normal3 where
-  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)
+  normal u = GL.normal (GL.Normal3 x y z) 
+    where Vec3 x y z = fromNormal u 
   normalv p = peek p >>= normal 
 
--------------------------------------------------------
+instance GL.Normal Vec3 where
+  normal (Vec3 x y z) = GL.normal (GL.Normal3 x y z) 
+  normalv p = peek p >>= normal 
 
+--------------------------------------------------------------------------------
 -- Color instances
   
 instance GL.Color Vec3 where
@@ -88,8 +203,7 @@
   secondaryColorv p = peek p >>= secondaryColor
 -}
 
--------------------------------------------------------
-
+--------------------------------------------------------------------------------
 -- TexCoord instances
 
 instance GL.TexCoord Vec2 where
@@ -110,8 +224,7 @@
   multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)
   multiTexCoordv unit p = peek p >>= multiTexCoord unit
 
--------------------------------------------------------
-    
+--------------------------------------------------------------------------------
 -- Vertex Attributes (experimental)
 
 class VertexAttrib' a where
@@ -130,16 +243,18 @@
   vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w 
 
 instance VertexAttrib' Normal2 where
-  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y
+  vertexAttrib loc u = GL.vertexAttrib2 loc x y
+    where Vec2 x y = fromNormal u 
 
 instance VertexAttrib' Normal3 where
-  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z
+  vertexAttrib loc u = GL.vertexAttrib3 loc x y z
+    where Vec3 x y z = fromNormal u 
 
 instance VertexAttrib' Normal4 where
-  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w
-
--------------------------------------------------------
-
+  vertexAttrib loc u = GL.vertexAttrib4 loc x y z w
+    where Vec4 x y z w = fromNormal u 
+   
+--------------------------------------------------------------------------------
 -- Uniform (again, experimental)
 
 -- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; 
@@ -181,3 +296,5 @@
   uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)
     
 #endif
+
+
diff --git a/Data/Vect/Double/Util/Dim2.hs b/Data/Vect/Double/Util/Dim2.hs
--- a/Data/Vect/Double/Util/Dim2.hs
+++ b/Data/Vect/Double/Util/Dim2.hs
@@ -4,12 +4,13 @@
 
 import Data.Vect.Flt.Base
 
--- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.
+-- | Example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.
 structVec2 :: [Flt] -> [Vec2]
 structVec2 [] = []
 structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls) 
 structVec2 _ = error "structVec2"
 
+-- | The opposite of "structVec2".
 destructVec2 :: [Vec2] -> [Flt]
 destructVec2 [] = []
 destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  
@@ -49,10 +50,13 @@
 angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt
 angle2' = angle2 . fromNormal
 
--- |Rotation matrix by a given angle (in radians), counterclockwise.
+-- | Rotation matrix by a given angle (in radians), counterclockwise.
 rotMatrix2 :: Flt -> Mat2
 rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a
 
+rotMatrixOrtho2 :: Flt -> Ortho2
+rotMatrixOrtho2 = toOrthoUnsafe . rotMatrix2
+
 rotate2 :: Flt -> Vec2 -> Vec2
 rotate2 a v = v .* (rotMatrix2 a) 
 
@@ -63,3 +67,5 @@
 -- |Rotates clockwise by 90 degrees.
 rotateCW :: Vec2 -> Vec2
 rotateCW (Vec2 x y) = Vec2 y (-x)
+
+
diff --git a/Data/Vect/Double/Util/Dim3.hs b/Data/Vect/Double/Util/Dim3.hs
--- a/Data/Vect/Double/Util/Dim3.hs
+++ b/Data/Vect/Double/Util/Dim3.hs
@@ -4,18 +4,26 @@
 
 import Data.Vect.Flt.Base
 
+--------------------------------------------------------------------------------
+
+-- | Example: @structVec3 [1,2,3,4,5,6] = [ Vec3 1 2 3 , Vec3 4 5 6]@.
 structVec3 :: [Flt] -> [Vec3]
 structVec3 [] = []
 structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) 
 structVec3 _ = error "structVec3"
 
+-- | The opposite of "structVec3".
 destructVec3 :: [Vec3] -> [Flt]
 destructVec3 [] = []
 destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)  
 
+--------------------------------------------------------------------------------
+
 det3 :: Vec3 -> Vec3 -> Vec3 -> Flt
 det3 u v w = det (u,v,w)
 
+--------------------------------------------------------------------------------
+
 translate3X :: Flt -> Vec3 -> Vec3
 translate3Y :: Flt -> Vec3 -> Vec3
 translate3Z :: Flt -> Vec3 -> Vec3
@@ -43,26 +51,34 @@
 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 -} Flt       -- ^ angle (in radians)
-         -> Normal3   -- ^ axis (should be a /unit/ vector!) 
-         -> Vec3      -- ^ vector
-         -> Vec3      -- ^ result
+--------------------------------------------------------------------------------
+
+rotate3' {- ' CPP is sensitive to primes -} 
+  :: Flt       -- ^ angle (in radians)
+  -> Normal3   -- ^ axis (should be a /unit/ vector!) 
+  -> Vec3      -- ^ vector
+  -> Vec3      -- ^ result
 rotate3' angle axis v = v .* (rotMatrix3' axis angle)
 
-rotate3 :: Flt    -- ^ angle (in radians)
-        -> Vec3   -- ^ axis (arbitrary nonzero vector)
-        -> Vec3   -- ^ vector
-        -> Vec3   -- ^ result
+rotate3 
+  :: Flt    -- ^ angle (in radians)
+  -> Vec3   -- ^ axis (arbitrary nonzero vector)
+  -> Vec3   -- ^ vector
+  -> Vec3   -- ^ 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/. 
+-- | Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
 rotMatrix3 :: Vec3 -> Flt -> Mat3
 rotMatrix3 v a = rotMatrix3' (mkNormal v) a
 
--- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
+rotMatrixOrtho3 :: Vec3 -> Flt -> Ortho3
+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 -} Normal3 -> Flt -> Mat3
-rotMatrix3' (Normal3 v) a = 
-  let c = cos a
+rotMatrix3' u a = 
+  let v = fromNormal u
+      c = cos a
       s = sin a
       m1 = scalarMul (1-c) (outer v v)
       x = _1 v
@@ -73,4 +89,65 @@
                 (Vec3 ( s*y) (-s*x)   c   )
   in (m1 &+ m2)
 
+rotMatrixOrtho3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Ortho3
+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 :: Normal3 -> Vec3 -> Vec3
+reflect u v = (s *& n) &- v where 
+  n = fromNormal u
+  s = 2 * (n &. v)
+
+reflect' :: Normal3 -> Normal3 -> Normal3
+reflect' u x = toNormalUnsafe $ reflect u (fromNormal x)
+  
+refract :: Flt -> Normal3 -> Vec3 -> Vec3
+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 -} 
+  :: Flt -> Normal3 -> Normal3 -> Normal3
+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 -}  
+  :: Flt -> Normal3 -> Normal3 -> Maybe Normal3
+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/Data/Vect/Double/Util/Dim4.hs b/Data/Vect/Double/Util/Dim4.hs
--- a/Data/Vect/Double/Util/Dim4.hs
+++ b/Data/Vect/Double/Util/Dim4.hs
@@ -40,7 +40,7 @@
 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))@.
@@ -76,8 +76,9 @@
 -- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors.
 -- Intended for multiplication on the /right/!
 rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4
-rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 
+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
diff --git a/Data/Vect/Double/Util/Projective.hs b/Data/Vect/Double/Util/Projective.hs
--- a/Data/Vect/Double/Util/Projective.hs
+++ b/Data/Vect/Double/Util/Projective.hs
@@ -1,6 +1,7 @@
 {-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}
 
--- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on
+-- | 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
 --
 -- >     _____
@@ -13,6 +14,8 @@
 -- 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.Flt.Util.Projective where
 
@@ -21,68 +24,73 @@
 
 import qualified Data.Vect.Flt.Util.Dim4 as Dim4
 
-class ExtendProjective v e | v->e where
-  extendProj     :: v -> e
-  extendProjWith :: Flt -> v -> e
-  extendProj = extendProjWith 1
-  
-instance ExtendProjective Vec2 Vec4 where
-  extendProj       (Vec2 x y) = Vec4 x y 0 1
-  extendProjWith w (Vec2 x y) = Vec4 x y 0 w
-  
-instance ExtendProjective Vec3 Vec4 where
-  extendProj       (Vec3 x y z) = Vec4 x y z 1
-  extendProjWith w (Vec3 x y z) = Vec4 x y z w
-
-instance ExtendProjective Vec4 Vec4 where
-  extendProj = id
-  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w
+--------------------------------------------------------------------------------
 
-instance ExtendProjective Mat2 Mat4 where
-  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)
-  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)
+rotMatrixProj4' :: {- ' CPP is sensitive to primes -}  Flt -> Normal3 -> Proj4
+rotMatrixProj4' angle axis = linear $ rotMatrix3' axis angle
 
-instance ExtendProjective Mat3 Mat4 where
-  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)
-  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)
+rotMatrixProj4 :: Flt -> Vec3 -> Proj4
+rotMatrixProj4 angle axis = linear $ rotMatrix3 axis angle
 
-rotMatrixProj :: Flt -> Normal3 -> Mat4
-rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle
+-- | synonym for "rotateAfterProj4"
+rotateProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateProj4 = rotateAfterProj4
 
-rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4
-rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle
+-- | Synonym for @\m -> m .*. rotMatrixProj4 angle axis@.
+rotateAfterProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateAfterProj4 angle axis m = m .*. (rotMatrixProj4' angle axis) 
 
-translMatrixProj :: Vec3 -> Mat4
-translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)
+-- | Synonym for @\m -> rotMatrixProj4 angle axis .*. m@.
+rotateBeforeProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateBeforeProj4 angle axis m = (rotMatrixProj4' angle axis) .*. m 
 
--- | we assume that the bottom-right corner is 1.
-translWithProj :: Vec3 -> Mat4 -> Mat4
-translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)
+---------------
 
-scaleMatrixProj :: Vec3 -> Mat4
-scaleMatrixProj v = diag $ extendProj v
+--scalingUniformProj3 :: Flt -> Proj3
+--scalingUniformProj3 x = scaling (Vec2 x x)
 
-scaleMatrixUniformProj :: Flt -> Mat4
-scaleMatrixUniformProj s = diag (Vec4 s s s 1)
+scalingUniformProj4 :: Flt -> Proj4
+scalingUniformProj4 x = scaling (Vec3 x x x)
 
-class ProjectiveAction v where
-  actProj :: v -> Mat4 -> v
- 
-instance ProjectiveAction Vec3 where
-  actProj v m = trim $ (extendProj v) .* m 
+-- | Equivalent to @\m -> scaling v .*. m@.
+scaleBeforeProj4 :: Vec3 -> Proj4 -> Proj4
+scaleBeforeProj4 (Vec3 u v w) p4 =  
+  toProjectiveUnsafe $ 
+    Mat4 (u*&a) (v*&b) (w*&c) t
+  where
+    Mat4 a b c t = fromProjective p4
 
-instance ProjectiveAction Vec4 where
-  actProj v m = v .* m 
+-- | Equivalent to @\m -> m .*. scaling v@.
+scaleAfterProj4 :: Vec3 -> Proj4 -> Proj4
+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
+    
+---------------
 
--- | When acting on unit vectors, we ignore the translation part.
-instance ProjectiveAction Normal3 where
-  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))
+-- | Synonym for "translateAfter4"
+translate4 :: Vec3 -> Proj4 -> Proj4
+translate4 = translateAfter4
 
--- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,
--- and the bottom-right corner is 1.
-invertProj :: Mat4 -> Mat4
-invertProj mat@(Mat4 u v w t) = 
-  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)
+-- | Equivalent to @\m -> m .*. translation v@.
+translateAfter4 :: Vec3 -> Proj4 -> Proj4
+translateAfter4 v p4 = 
+  toProjectiveUnsafe $
+    Mat4 r1 r2 r3 (extendWith 0 v &+ r4)
   where
-    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)
-    
+    Mat4 r1 r2 r3 r4 = fromProjective p4 
+
+-- | Equivalent to @\m -> translation v .*. m@.
+translateBefore4 :: Vec3 -> Proj4 -> Proj4
+translateBefore4 v p4 = 
+  toProjectiveUnsafe $ 
+    Mat4 r1 r2 r3 (extendWith 0 u &+ r4) 
+  where 
+   u = v .* (trim mat :: Mat3) 
+   mat@(Mat4 r1 r2 r3 r4) = fromProjective p4
+   
+---------------
+
diff --git a/Data/Vect/Float/Base.hs b/Data/Vect/Float/Base.hs
--- a/Data/Vect/Float/Base.hs
+++ b/Data/Vect/Float/Base.hs
@@ -1,11 +1,34 @@
 {-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}
 
-module Data.Vect.Flt.Base where
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, GeneralizedNewtypeDeriving #-}
 
+module Data.Vect.Flt.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 Control.Monad
 import System.Random  
 import Foreign
 
+--------------------------------------------------------------------------------
 -- class declarations
 
 class AbelianGroup g where
@@ -20,15 +43,17 @@
 vecSum :: AbelianGroup g => [g] -> g
 vecSum l = foldl (&+) zero l 
 
-class (AbelianGroup r) => 
-      Ring r where
+class MultSemiGroup r where
   (.*.) :: r -> r -> r
   one   :: r
 
+class (AbelianGroup r, MultSemiGroup r) => Ring r 
+
 infixl 7 .*. 
 
-ringProduct :: Ring r => [r] -> r
-ringProduct l = foldl (.*.) one l
+-- 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
@@ -75,7 +100,7 @@
   lensqr = normsqr
   dotprod :: v -> v -> Flt
   normsqr v = (v &. v)  
-  norm = sqrt.lensqr
+  norm = sqrt . lensqr
   dotprod = (&.)
 
 infix 7 &.
@@ -112,26 +137,28 @@
   fromNormalRadius :: Flt -> u -> v
   fromNormalRadius t n = t *& fromNormal n 
 
--- | projects the first vector onto the direction of the second (unit) vector
+-- | Projects the first vector onto the direction of the second (unit) vector
 project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v
 project' what dir = projectUnsafe what (fromNormal dir)
 
--- | direction (second argument) is assumed to be a /unit/ vector!
+-- | Direction (second argument) is assumed to be a /unit/ vector!
 projectUnsafe :: (Vector v, DotProd v) => v -> v -> v
 projectUnsafe what dir = what &- dir &* (what &. dir)
 
 project :: (Vector v, DotProd v) => v -> v -> v
 project what dir = what &- dir &* ((what &. dir) / (dir &. dir))
 
--- | since unit vectors are not a group, we need a separate function.
+-- | Since unit vectors are not a group, we need a separate function.
 flipNormal :: UnitVector v n => n -> n 
 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
@@ -166,13 +193,50 @@
 "inverse is an involution"    forall m. inverse (inverse m) = m
   #-}
   
+class Matrix m => Orthogonal m o | m->o, o->m where  
+  fromOrtho     :: o -> m 
+  toOrthoUnsafe :: m -> o
+  
+class (AbelianGroup m, Matrix m) => MatrixNorms m where
+  frobeniusNorm  :: m -> Flt       -- ^ the frobenius norm (= euclidean norm in the space of matrices)
+  matrixDistance :: m -> m -> Flt  -- ^ euclidean distance in the space of matrices
+  operatorNorm   :: m -> Flt       -- ^ (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 m where
   det :: m -> Flt    
-    
+
+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 v, UnitVector v u, Matrix m, Vector m, Tensor m v) => u -> m
+householder u = idmtx &- (2 *& outer v v) 
+  where v = fromNormal u
+
+householderOrtho :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v, Orthogonal m o) => u -> o
+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 v, Orthogonal n o, Diagonal v n) => Projective 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 -> m
+  toProjectiveUnsafe :: m -> p
+  orthogonal         :: o -> p
+  linear             :: n -> p
+  translation        :: v -> p
+  scaling            :: v -> p
+
+--------------------------------------------------------------------------------
 -- Vec / Mat datatypes
  
 data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt 
@@ -182,16 +246,16 @@
 data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt 
   deriving (Read,Show)
 
--- | these are /row/ vectors 
+-- | The components are /row/ vectors 
 data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show)
 data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show)
 data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)
 
 -- | The assumption when dealing with these is always that they are of unit length.
 -- Also, interpolation works differently.
-newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) 
-newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) 
-newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) 
+newtype Normal2 = Normal2 Vec2 deriving (Read,Show,Storable,DotProd,Dimension) 
+newtype Normal3 = Normal3 Vec3 deriving (Read,Show,Storable,DotProd,Dimension) 
+newtype Normal4 = Normal4 Vec4 deriving (Read,Show,Storable,DotProd,Dimension) 
 
 mkVec2 :: (Flt,Flt) -> Vec2
 mkVec3 :: (Flt,Flt,Flt) -> Vec3
@@ -201,6 +265,19 @@
 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 = Ortho2 Mat2 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+newtype Ortho3 = Ortho3 Mat3 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+newtype Ortho4 = Ortho4 Mat4 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+
+-- | Projective matrices, encoding affine transformations in dimension one less.
+newtype Proj3 = Proj3 Mat3 deriving (Read,Show,Storable,MultSemiGroup)
+newtype Proj4 = Proj4 Mat4 deriving (Read,Show,Storable,MultSemiGroup)
+
+--------------------------------------------------------------------------------
 -- Unit vectors
   
 instance UnitVector Vec2 Normal2 where
@@ -218,47 +295,163 @@
   fromNormal (Normal4 v) = v 
   toNormalUnsafe = Normal4
 
-rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)
-rndUnit g = 
+_rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)
+_rndUnit g = 
   if d > 0.01
     then ( v &* (1.0/d) , h )
-    else rndUnit h
+    else _rndUnit h
   where
     (v,h) = random g
     d = norm v
     
 instance Random Normal2 where
-  random g = let (v,h) = rndUnit g in (Normal2 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal2 v, h)  
   randomR _ = random
 
 instance Random Normal3 where
-  random g = let (v,h) = rndUnit g in (Normal3 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal3 v, h)  
   randomR _ = random
 
 instance Random Normal4 where
-  random g = let (v,h) = rndUnit g in (Normal4 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal4 v, h)  
   randomR _ = random
 
-{-
-instance Storable Normal2 where
-  alignment _ = alignment (undefined::Vec2)
-  sizeOf    _ = sizeOf    (undefined::Vec2)
-  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)
-  poke p (Normal2 v) = poke (castPtr p) v  
- 
-instance Storable Normal3 where
-  alignment _ = alignment (undefined::Vec3)
-  sizeOf    _ = sizeOf    (undefined::Vec3)
-  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)
-  poke p (Normal3 v) = poke (castPtr p) v  
+instance CrossProd Normal3 where
+  crossprod (Normal3 v) (Normal3 w) = mkNormal (crossprod v w)
 
-instance Storable Normal4 where
-  alignment _ = alignment (undefined::Vec4)
-  sizeOf    _ = sizeOf    (undefined::Vec4)
-  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)
-  poke p (Normal4 v) = poke (castPtr p) v  
--}
+--------------------------------------------------------------------------------
+-- Orthogonal matrices
 
+instance Orthogonal Mat2 Ortho2 where
+  fromOrtho (Ortho2 o) = o
+  toOrthoUnsafe = Ortho2
+
+instance Orthogonal Mat3 Ortho3 where
+  fromOrtho (Ortho3 o) = o
+  toOrthoUnsafe = Ortho3 
+
+instance Orthogonal Mat4 Ortho4 where
+  fromOrtho (Ortho4 o) = o
+  toOrthoUnsafe = Ortho4
+
+------
+
+instance Matrix Ortho2 where
+  transpose (Ortho2 o) = Ortho2 (transpose o)
+  idmtx = Ortho2 idmtx
+  inverse = transpose
+
+instance Matrix Ortho3 where
+  transpose (Ortho3 o) = Ortho3 (transpose o)
+  idmtx = Ortho3 idmtx
+  inverse = transpose
+
+instance Matrix Ortho4 where
+  transpose (Ortho4 o) = Ortho4 (transpose o)
+  idmtx = Ortho4 idmtx
+  inverse = transpose
+
+------
+
+instance Random Ortho2 where
+  random g = let (o,h) = _rndOrtho2 g in (toOrthoUnsafe (_flip1stRow2 o), h)
+  randomR _ = random
+
+instance Random Ortho3 where
+  random g = let (o,h) = _rndOrtho3 g in (toOrthoUnsafe (             o), h)
+  randomR _ = random
+
+instance Random Ortho4 where
+  random g = let (o,h) = _rndOrtho4 g in (toOrthoUnsafe (_flip1stRow4 o), h)
+  randomR _ = random
+
+------
+
+-- determinant will be -1
+_rndOrtho2 :: RandomGen g => g -> (Mat2, g)
+_rndOrtho2 g = (h2, g1) where
+  h2 = householder u2 :: Mat2 
+  (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 :: RandomGen g => g -> (Mat3, g) 
+_rndOrtho3 g = ( (h3 .*. m3), g2) where
+  m3 = (extendWith :: Flt -> Mat2 -> Mat3) 1 o2 
+  h3 = householder u3 :: Mat3
+  (u3,g1) = random g
+  (o2,g2) = _rndOrtho2 g1
+
+-- determinant will be -1
+_rndOrtho4 :: RandomGen g => g -> (Mat4, g) 
+_rndOrtho4 g = ( (h4 .*. m4), g2) where
+  m4 = (extendWith :: Flt -> Mat3 -> Mat4) 1 o3 
+  h4 = householder u4 :: Mat4
+  (u4,g1) = random g
+  (o3,g2) = _rndOrtho3 g1
+
+------
+
+_flip1stRow2 :: Mat2 -> Mat2
+_flip1stRow2 (Mat2 a b) = Mat2 (neg a) b
+
+_flip1stRow3 :: Mat3 -> Mat3
+_flip1stRow3 (Mat3 a b c) = Mat3 (neg a) b c
+
+_flip1stRow4 :: Mat4 -> Mat4
+_flip1stRow4 (Mat4 a b c d) = Mat4 (neg a) b c d
+
+--------------------------------------------------------------------------------
+-- projective matrices
+  
+instance Projective Vec2 Mat2 Ortho2 Mat3 Proj3 where
+  fromProjective (Proj3 m) = m
+  toProjectiveUnsafe = Proj3
+  orthogonal = Proj3 . extendWith 1 . fromOrtho
+  linear     = Proj3 . extendWith 1
+  translation v = Proj3 $ Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (extendWith 1 v)
+  scaling     v = Proj3 $ diag (extendWith 1 v)
+  
+instance Projective Vec3 Mat3 Ortho3 Mat4 Proj4 where
+  fromProjective (Proj4 m) = m
+  toProjectiveUnsafe = Proj4
+  orthogonal = Proj4 . extendWith 1 . fromOrtho 
+  linear     = Proj4 . extendWith 1
+  translation v = Proj4 $ Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (extendWith 1 v)
+  scaling     v = Proj4 $ diag (extendWith 1 v)
+
+instance Matrix Proj3 where
+  idmtx = Proj3 idmtx
+  transpose (Proj3 m) = Proj3 (transpose m)
+  inverse = _invertProj3
+
+instance Matrix Proj4 where
+  idmtx = Proj4 idmtx
+  transpose (Proj4 m) = Proj4 (transpose m)
+  inverse = _invertProj4
+
+_invertProj3 :: Proj3 -> Proj3
+_invertProj3 (Proj3 mat@(Mat3 _ _ t)) = 
+  Proj3 $ Mat3 (extendZero a) (extendZero b) (extendWith 1 t') 
+  where
+    t' = neg $ (trim t :: Vec2) .* invm2 
+    invm2@(Mat2 a b) = inverse $ (trim mat :: Mat2)
+
+-- Inverts a projective 4x4 matrix. But you can simply use "inverse" instead.
+-- We assume that the bottom-right corner is 1.
+_invertProj4 :: Proj4 -> Proj4
+_invertProj4 (Proj4 mat@(Mat4 _ _ _ t)) = 
+  Proj4 $ Mat4 (extendZero a) (extendZero b) (extendZero c) (extendWith 1 t') 
+  where
+    t' = neg $ (trim t :: Vec3) .* invm3 
+    invm3@(Mat3 a b c) = inverse $ (trim mat :: Mat3)
+
+--------------------------------------------------------------------------------
 -- Vec2 instances
 
 instance HasCoordinates Vec2 Flt where
@@ -314,7 +507,10 @@
         k = sizeOf (undefined::Flt)
     poke        p   x
     pokeByteOff p k y
-               
+
+instance Dimension Vec2 where dim _ = 2
+
+--------------------------------------------------------------------------------                    
 -- Mat2 instances
 
 instance HasCoordinates Mat2 Vec2 where
@@ -336,19 +532,21 @@
   (&+) (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   -- (zero::Vec2) (zero::Vec2)
-
+  zero = Mat2 zero zero  
+  
 instance Vector 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 Ring Mat2 where
+instance MultSemiGroup Mat2 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 Ring Mat2
+
 instance LeftModule Mat2 Vec2 where
   lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) 
   
@@ -362,11 +560,6 @@
   outer (Vec2 a b) (Vec2 x y) = Mat2
     (Vec2 (a*x) (a*y))
     (Vec2 (b*x) (b*y))
-{-
-  outer v w = 
-    let full = Mat2 (Vec2 1 1) (Vec2 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
 instance Determinant Mat2 where
   det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c 
@@ -393,6 +586,24 @@
     poke        p   r1
     pokeByteOff p k r2
 
+instance Random Mat2 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 Dimension Mat2 where dim _ = 2
+     
+instance MatrixNorms Mat2 where 
+  frobeniusNorm (Mat2 r1 r2) =  
+    sqrt $
+      normsqr r1 + 
+      normsqr r2
+     
+--------------------------------------------------------------------------------     
 -- Vec3 instances
 
 instance HasCoordinates Vec3 Flt where
@@ -454,7 +665,10 @@
     poke        p       x
     pokeByteOff p (k  ) y
     pokeByteOff p (k+k) z
-   
+
+instance Dimension Vec3 where dim _ = 3
+
+--------------------------------------------------------------------------------   
 -- Mat3 instances
 
 instance HasCoordinates Mat3 Vec3 where
@@ -495,13 +709,13 @@
   (&+) (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   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)
+  zero = Mat3 zero zero zero 
 
 instance Vector 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 Ring Mat3 where
+instance MultSemiGroup Mat3 where
   (.*.) (Mat3 r1 r2 r3) n = 
     let (Mat3 c1 c2 c3) = transpose n
     in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))
@@ -509,6 +723,8 @@
             (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))
   one = idmtx 
 
+instance Ring Mat3
+
 instance LeftModule Mat3 Vec3 where
   lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)
   
@@ -523,11 +739,6 @@
     (Vec3 (a*x) (a*y) (a*z))
     (Vec3 (b*x) (b*y) (b*z))
     (Vec3 (c*x) (c*y) (c*z))
-{-
-  outer v w = 
-    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
 instance Determinant Mat3 where
   det (Mat3 r1 r2 r3) = det (r1,r2,r3)
@@ -556,6 +767,26 @@
     pokeByteOff p (k  ) r2
     pokeByteOff p (k+k) r3
 
+instance Random Mat3 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 Dimension Mat3 where dim _ = 3
+  
+instance MatrixNorms Mat3 where 
+  frobeniusNorm (Mat3 r1 r2 r3)  = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 
+    
+--------------------------------------------------------------------------------
 -- Vec4 instances
 
 instance HasCoordinates Vec4 Flt where
@@ -615,6 +846,9 @@
     pokeByteOff p (k+k) z
     pokeByteOff p (3*k) w
 
+instance Dimension Vec4 where dim _ = 4
+
+--------------------------------------------------------------------------------
 -- Mat4 instances
 
 instance HasCoordinates Mat4 Vec4 where
@@ -642,7 +876,7 @@
   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 Ring Mat4 where
+instance MultSemiGroup Mat4 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))
@@ -651,6 +885,8 @@
             (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))
   one = idmtx 
 
+instance Ring Mat4
+
 instance LeftModule Mat4 Vec4 where
   lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)
   
@@ -666,14 +902,10 @@
     (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))
-{-
-  outer v w = 
-    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
---instance Determinant Mat4 where
---  det (Mat4 r1 r2 r3 r4)
+instance Determinant Mat4 where
+  det = error "det/Mat4: not implemented yet" 
+  -- det (Mat4 r1 r2 r3 r4) = 
 
 {-
 instance Show Mat4 where
@@ -701,6 +933,28 @@
     pokeByteOff p (k+k) r3
     pokeByteOff p (3*k) r4
 
+instance Random Mat4 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 Dimension Mat4 where dim _ = 4
+   
+instance MatrixNorms Mat4 where 
+  frobeniusNorm (Mat4 r1 r2 r3 r4) = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 + 
+      normsqr r4  
+    
+--------------------------------------------------------------------------------
 -- Extend instances
 
 instance Extend Vec2 Vec3 where
@@ -719,17 +973,20 @@
   trim (Vec4 x y z _)       = Vec3 x y z
 
 instance Extend Mat2 Mat3 where
-  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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 Extend Mat2 Mat4 where
-  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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 Extend Mat3 Mat4 where
-  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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/Data/Vect/Float/GramSchmidt.hs b/Data/Vect/Float/GramSchmidt.hs
--- a/Data/Vect/Float/GramSchmidt.hs
+++ b/Data/Vect/Float/GramSchmidt.hs
@@ -10,7 +10,7 @@
 
 import Data.Vect.Flt.Base
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
 liftPair :: (a -> b) -> (a,a) -> (b,b)
 liftPair f (x,y) = (f x, f y)
@@ -21,7 +21,7 @@
 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
@@ -33,7 +33,7 @@
 "gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a 
   #-}
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
 instance GramSchmidt (Vec2,Vec2) where
   gramSchmidt = gramSchmidtPair
diff --git a/Data/Vect/Float/OpenGL.hs b/Data/Vect/Float/OpenGL.hs
--- a/Data/Vect/Float/OpenGL.hs
+++ b/Data/Vect/Float/OpenGL.hs
@@ -3,21 +3,92 @@
 -- TODO: the pointer versions of these functions should be really implemented 
 -- via the pointer versions of the original opengl functions...
 
--- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.
+-- | OpenGL support, including 'Vertex', 'TexCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.
  
 module Data.Vect.Flt.OpenGL where
 
 import Control.Monad
 import Data.Vect.Flt.Base
+import Data.Vect.Flt.Util.Projective
 import qualified Graphics.Rendering.OpenGL as GL
 
 import Foreign
 
-import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)
+import Graphics.Rendering.OpenGL hiding 
+  ( Normal3 , rotate , translate , scale
+  , matrix , currentMatrix , withMatrix , multMatrix 
+  )
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
-{-# SPECIALISE radianToDegrees :: Float -> Float #-}
+-- | There should be a big warning here about the different conventions, 
+-- hidden transpositions, and all the confusion this will inevitably cause...
+--
+-- As it stands, 
+--
+-- > glRotate t1 axis1 >> glRotate t2 axis2 >> glRotate t3 axis3
+-- 
+-- has the same result as
+--
+-- > multMatrix (rotMatrixProj4 t3 axis3 .*. rotMatrixProj4 t2 axis2 .*. rotMatrixProj4 t1 axis1)
+--
+-- because at the interface of OpenGL and this library there is a transposition
+-- to compensate for the different conventions. (This transposition is implicit
+-- in the code, because the way the matrices are stored in the memory is also
+-- different: OpenGL stores them column-major, and we store them row-major).
+
+class ToOpenGLMatrix m where
+  makeGLMatrix :: m -> IO (GLmatrix Flt)
+
+class FromOpenGLMatrix m where
+  peekGLMatrix :: GLmatrix Flt -> IO m
+  
+setMatrix :: ToOpenGLMatrix m => Maybe MatrixMode -> m -> IO ()
+setMatrix mode m = makeGLMatrix m >>= \x -> GL.matrix mode $= x
+ 
+getMatrix :: FromOpenGLMatrix m => Maybe MatrixMode -> IO m
+getMatrix mode = get (GL.matrix mode) >>= peekGLMatrix
+
+matrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => Maybe MatrixMode -> StateVar m
+matrix mode = makeStateVar (getMatrix mode) (setMatrix mode)
+
+currentMatrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => StateVar m
+currentMatrix = matrix Nothing
+
+multMatrix :: ToOpenGLMatrix m => m -> IO ()
+multMatrix m = makeGLMatrix m >>= GL.multMatrix
+
+instance ToOpenGLMatrix Mat4 where
+  makeGLMatrix m = GL.withNewMatrix GL.ColumnMajor (flip poke m . castPtr) 
+ 
+instance FromOpenGLMatrix Mat4 where
+  -- huh? GL.withMatrix is strange
+  peekGLMatrix x = GL.withMatrix x $ \_ p -> peek (castPtr p)
+  
+instance ToOpenGLMatrix Mat3 where
+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)
+ 
+instance ToOpenGLMatrix Mat2 where
+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)
+
+instance ToOpenGLMatrix Ortho4 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat4)
+
+instance ToOpenGLMatrix Ortho3 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat3)
+
+instance ToOpenGLMatrix Ortho2 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat2)
+
+instance ToOpenGLMatrix Proj4 where
+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat4)
+
+instance ToOpenGLMatrix Proj3 where
+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat3)
+  
+--------------------------------------------------------------------------------
+
+{-# SPECIALISE radianToDegrees :: Float  -> Float  #-}
 {-# SPECIALISE radianToDegrees :: Double -> Double #-}
 radianToDegrees :: RealFrac a => a -> a
 radianToDegrees x = x * 57.295779513082322
@@ -28,20 +99,61 @@
 degreesToRadian x = x * 1.7453292519943295e-2
 
 -- | The angle is in radians. (WARNING: OpenGL uses degrees!)
-rotate :: Flt -> Vec3 -> IO ()
-rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)
+glRotate :: Flt -> Vec3 -> IO ()
+glRotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)
 
-translate :: Vec3 -> IO ()
-translate (Vec3 x y z) = GL.translate (Vector3 x y z)
+glTranslate :: Vec3 -> IO ()
+glTranslate (Vec3 x y z) = GL.translate (Vector3 x y z)
 
-scale3 :: Vec3 -> IO ()
-scale3 (Vec3 x y z) = GL.scale x y z
+glScale3 :: Vec3 -> IO ()
+glScale3 (Vec3 x y z) = GL.scale x y z
 
-scale :: Flt -> IO ()
-scale x = GL.scale x x x
+glScale :: Flt -> IO ()
+glScale x = GL.scale x x x
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
+ 
+-- | \"Orthogonal projecton\" matrix, a la OpenGL 
+-- (the corresponding functionality is removed in OpenGL 3.1)
+orthoMatrix 
+  :: (Flt,Flt)   -- ^ (left,right)
+  -> (Flt,Flt)   -- ^ (bottom,top)
+  -> (Flt,Flt)   -- ^ (near,far)
+  -> Mat4 
+orthoMatrix (l,r) (b,t) (n,f) = Mat4
+  (Vec4 (2/(r-l)) 0 0 0)
+  (Vec4 0 (2/(t-b)) 0 0)
+  (Vec4 0 0 (-2/(f-n)) 0)
+  (Vec4 (-(r+l)/(r-l)) (-(t+b)/(t-b)) (-(f+n)/(f-n)) 1)
+  
+-- | The same as "orthoMatrix", but with a different parametrization.
+orthoMatrix2 {- ' CPP is sensitive to primes -}
+  :: Vec3     -- ^ (left,top,near)
+  -> Vec3     -- ^ (right,bottom,far)
+  -> Mat4 
+orthoMatrix2 (Vec3 l t n) (Vec3 r b f) = orthoMatrix (l,r) (b,t) (n,f)
 
+-- | \"Perspective projecton\" matrix, a la OpenGL 
+-- (the corresponding functionality is removed in OpenGL 3.1).
+frustumMatrix
+  :: (Flt,Flt)   -- ^ (left,right)
+  -> (Flt,Flt)   -- ^ (bottom,top)
+  -> (Flt,Flt)   -- ^ (near,far)
+  -> Mat4 
+frustumMatrix (l,r) (b,t) (n,f) = Mat4
+  (Vec4 (2*n/(r-l)) 0 0 0)
+  (Vec4 0 (2*n/(t-b)) 0 0)
+  (Vec4 ((r+l)/(r-l)) ((t+b)/(t-b)) (-(f+n)/(f-n)) (-1))
+  (Vec4 0 0 (-2*f*n*(f-n)) 0)
+  
+-- | The same as "frustumMatrix", but with a different parametrization.
+frustumMatrix2 {- ' CPP is sensitive to primes -}
+  :: Vec3     -- ^ (left,top,near)
+  -> Vec3     -- ^ (right,bottom,far)
+  -> Mat4 
+frustumMatrix2 (Vec3 l t n) (Vec3 r b f) = frustumMatrix (l,r) (b,t) (n,f)
+
+--------------------------------------------------------------------------------
 -- Vertex instances
 
 instance GL.Vertex Vec2 where
@@ -56,17 +168,20 @@
   vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)
   vertexv p = peek p >>= vertex   
 
--------------------------------------------------------
-
+--------------------------------------------------------------------------------
 -- the Normal instance
 -- note that there is no Normal2\/Normal4 in the OpenGL binding
 
 instance GL.Normal Normal3 where
-  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)
+  normal u = GL.normal (GL.Normal3 x y z) 
+    where Vec3 x y z = fromNormal u 
   normalv p = peek p >>= normal 
 
--------------------------------------------------------
+instance GL.Normal Vec3 where
+  normal (Vec3 x y z) = GL.normal (GL.Normal3 x y z) 
+  normalv p = peek p >>= normal 
 
+--------------------------------------------------------------------------------
 -- Color instances
   
 instance GL.Color Vec3 where
@@ -88,8 +203,7 @@
   secondaryColorv p = peek p >>= secondaryColor
 -}
 
--------------------------------------------------------
-
+--------------------------------------------------------------------------------
 -- TexCoord instances
 
 instance GL.TexCoord Vec2 where
@@ -110,8 +224,7 @@
   multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)
   multiTexCoordv unit p = peek p >>= multiTexCoord unit
 
--------------------------------------------------------
-    
+--------------------------------------------------------------------------------
 -- Vertex Attributes (experimental)
 
 class VertexAttrib' a where
@@ -130,16 +243,18 @@
   vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w 
 
 instance VertexAttrib' Normal2 where
-  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y
+  vertexAttrib loc u = GL.vertexAttrib2 loc x y
+    where Vec2 x y = fromNormal u 
 
 instance VertexAttrib' Normal3 where
-  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z
+  vertexAttrib loc u = GL.vertexAttrib3 loc x y z
+    where Vec3 x y z = fromNormal u 
 
 instance VertexAttrib' Normal4 where
-  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w
-
--------------------------------------------------------
-
+  vertexAttrib loc u = GL.vertexAttrib4 loc x y z w
+    where Vec4 x y z w = fromNormal u 
+   
+--------------------------------------------------------------------------------
 -- Uniform (again, experimental)
 
 -- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; 
@@ -181,3 +296,5 @@
   uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)
     
 #endif
+
+
diff --git a/Data/Vect/Float/Util/Dim2.hs b/Data/Vect/Float/Util/Dim2.hs
--- a/Data/Vect/Float/Util/Dim2.hs
+++ b/Data/Vect/Float/Util/Dim2.hs
@@ -4,12 +4,13 @@
 
 import Data.Vect.Flt.Base
 
--- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.
+-- | Example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.
 structVec2 :: [Flt] -> [Vec2]
 structVec2 [] = []
 structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls) 
 structVec2 _ = error "structVec2"
 
+-- | The opposite of "structVec2".
 destructVec2 :: [Vec2] -> [Flt]
 destructVec2 [] = []
 destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  
@@ -49,10 +50,13 @@
 angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt
 angle2' = angle2 . fromNormal
 
--- |Rotation matrix by a given angle (in radians), counterclockwise.
+-- | Rotation matrix by a given angle (in radians), counterclockwise.
 rotMatrix2 :: Flt -> Mat2
 rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a
 
+rotMatrixOrtho2 :: Flt -> Ortho2
+rotMatrixOrtho2 = toOrthoUnsafe . rotMatrix2
+
 rotate2 :: Flt -> Vec2 -> Vec2
 rotate2 a v = v .* (rotMatrix2 a) 
 
@@ -63,3 +67,5 @@
 -- |Rotates clockwise by 90 degrees.
 rotateCW :: Vec2 -> Vec2
 rotateCW (Vec2 x y) = Vec2 y (-x)
+
+
diff --git a/Data/Vect/Float/Util/Dim3.hs b/Data/Vect/Float/Util/Dim3.hs
--- a/Data/Vect/Float/Util/Dim3.hs
+++ b/Data/Vect/Float/Util/Dim3.hs
@@ -4,18 +4,26 @@
 
 import Data.Vect.Flt.Base
 
+--------------------------------------------------------------------------------
+
+-- | Example: @structVec3 [1,2,3,4,5,6] = [ Vec3 1 2 3 , Vec3 4 5 6]@.
 structVec3 :: [Flt] -> [Vec3]
 structVec3 [] = []
 structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) 
 structVec3 _ = error "structVec3"
 
+-- | The opposite of "structVec3".
 destructVec3 :: [Vec3] -> [Flt]
 destructVec3 [] = []
 destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)  
 
+--------------------------------------------------------------------------------
+
 det3 :: Vec3 -> Vec3 -> Vec3 -> Flt
 det3 u v w = det (u,v,w)
 
+--------------------------------------------------------------------------------
+
 translate3X :: Flt -> Vec3 -> Vec3
 translate3Y :: Flt -> Vec3 -> Vec3
 translate3Z :: Flt -> Vec3 -> Vec3
@@ -43,26 +51,34 @@
 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 -} Flt       -- ^ angle (in radians)
-         -> Normal3   -- ^ axis (should be a /unit/ vector!) 
-         -> Vec3      -- ^ vector
-         -> Vec3      -- ^ result
+--------------------------------------------------------------------------------
+
+rotate3' {- ' CPP is sensitive to primes -} 
+  :: Flt       -- ^ angle (in radians)
+  -> Normal3   -- ^ axis (should be a /unit/ vector!) 
+  -> Vec3      -- ^ vector
+  -> Vec3      -- ^ result
 rotate3' angle axis v = v .* (rotMatrix3' axis angle)
 
-rotate3 :: Flt    -- ^ angle (in radians)
-        -> Vec3   -- ^ axis (arbitrary nonzero vector)
-        -> Vec3   -- ^ vector
-        -> Vec3   -- ^ result
+rotate3 
+  :: Flt    -- ^ angle (in radians)
+  -> Vec3   -- ^ axis (arbitrary nonzero vector)
+  -> Vec3   -- ^ vector
+  -> Vec3   -- ^ 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/. 
+-- | Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
 rotMatrix3 :: Vec3 -> Flt -> Mat3
 rotMatrix3 v a = rotMatrix3' (mkNormal v) a
 
--- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
+rotMatrixOrtho3 :: Vec3 -> Flt -> Ortho3
+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 -} Normal3 -> Flt -> Mat3
-rotMatrix3' (Normal3 v) a = 
-  let c = cos a
+rotMatrix3' u a = 
+  let v = fromNormal u
+      c = cos a
       s = sin a
       m1 = scalarMul (1-c) (outer v v)
       x = _1 v
@@ -73,4 +89,65 @@
                 (Vec3 ( s*y) (-s*x)   c   )
   in (m1 &+ m2)
 
+rotMatrixOrtho3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Ortho3
+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 :: Normal3 -> Vec3 -> Vec3
+reflect u v = (s *& n) &- v where 
+  n = fromNormal u
+  s = 2 * (n &. v)
+
+reflect' :: Normal3 -> Normal3 -> Normal3
+reflect' u x = toNormalUnsafe $ reflect u (fromNormal x)
+  
+refract :: Flt -> Normal3 -> Vec3 -> Vec3
+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 -} 
+  :: Flt -> Normal3 -> Normal3 -> Normal3
+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 -}  
+  :: Flt -> Normal3 -> Normal3 -> Maybe Normal3
+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/Data/Vect/Float/Util/Dim4.hs b/Data/Vect/Float/Util/Dim4.hs
--- a/Data/Vect/Float/Util/Dim4.hs
+++ b/Data/Vect/Float/Util/Dim4.hs
@@ -40,7 +40,7 @@
 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))@.
@@ -76,8 +76,9 @@
 -- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors.
 -- Intended for multiplication on the /right/!
 rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4
-rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 
+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
diff --git a/Data/Vect/Float/Util/Projective.hs b/Data/Vect/Float/Util/Projective.hs
--- a/Data/Vect/Float/Util/Projective.hs
+++ b/Data/Vect/Float/Util/Projective.hs
@@ -1,6 +1,7 @@
 {-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}
 
--- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on
+-- | 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
 --
 -- >     _____
@@ -13,6 +14,8 @@
 -- 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.Flt.Util.Projective where
 
@@ -21,68 +24,73 @@
 
 import qualified Data.Vect.Flt.Util.Dim4 as Dim4
 
-class ExtendProjective v e | v->e where
-  extendProj     :: v -> e
-  extendProjWith :: Flt -> v -> e
-  extendProj = extendProjWith 1
-  
-instance ExtendProjective Vec2 Vec4 where
-  extendProj       (Vec2 x y) = Vec4 x y 0 1
-  extendProjWith w (Vec2 x y) = Vec4 x y 0 w
-  
-instance ExtendProjective Vec3 Vec4 where
-  extendProj       (Vec3 x y z) = Vec4 x y z 1
-  extendProjWith w (Vec3 x y z) = Vec4 x y z w
-
-instance ExtendProjective Vec4 Vec4 where
-  extendProj = id
-  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w
+--------------------------------------------------------------------------------
 
-instance ExtendProjective Mat2 Mat4 where
-  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)
-  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)
+rotMatrixProj4' :: {- ' CPP is sensitive to primes -}  Flt -> Normal3 -> Proj4
+rotMatrixProj4' angle axis = linear $ rotMatrix3' axis angle
 
-instance ExtendProjective Mat3 Mat4 where
-  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)
-  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)
+rotMatrixProj4 :: Flt -> Vec3 -> Proj4
+rotMatrixProj4 angle axis = linear $ rotMatrix3 axis angle
 
-rotMatrixProj :: Flt -> Normal3 -> Mat4
-rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle
+-- | synonym for "rotateAfterProj4"
+rotateProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateProj4 = rotateAfterProj4
 
-rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4
-rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle
+-- | Synonym for @\m -> m .*. rotMatrixProj4 angle axis@.
+rotateAfterProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateAfterProj4 angle axis m = m .*. (rotMatrixProj4' angle axis) 
 
-translMatrixProj :: Vec3 -> Mat4
-translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)
+-- | Synonym for @\m -> rotMatrixProj4 angle axis .*. m@.
+rotateBeforeProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateBeforeProj4 angle axis m = (rotMatrixProj4' angle axis) .*. m 
 
--- | we assume that the bottom-right corner is 1.
-translWithProj :: Vec3 -> Mat4 -> Mat4
-translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)
+---------------
 
-scaleMatrixProj :: Vec3 -> Mat4
-scaleMatrixProj v = diag $ extendProj v
+--scalingUniformProj3 :: Flt -> Proj3
+--scalingUniformProj3 x = scaling (Vec2 x x)
 
-scaleMatrixUniformProj :: Flt -> Mat4
-scaleMatrixUniformProj s = diag (Vec4 s s s 1)
+scalingUniformProj4 :: Flt -> Proj4
+scalingUniformProj4 x = scaling (Vec3 x x x)
 
-class ProjectiveAction v where
-  actProj :: v -> Mat4 -> v
- 
-instance ProjectiveAction Vec3 where
-  actProj v m = trim $ (extendProj v) .* m 
+-- | Equivalent to @\m -> scaling v .*. m@.
+scaleBeforeProj4 :: Vec3 -> Proj4 -> Proj4
+scaleBeforeProj4 (Vec3 u v w) p4 =  
+  toProjectiveUnsafe $ 
+    Mat4 (u*&a) (v*&b) (w*&c) t
+  where
+    Mat4 a b c t = fromProjective p4
 
-instance ProjectiveAction Vec4 where
-  actProj v m = v .* m 
+-- | Equivalent to @\m -> m .*. scaling v@.
+scaleAfterProj4 :: Vec3 -> Proj4 -> Proj4
+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
+    
+---------------
 
--- | When acting on unit vectors, we ignore the translation part.
-instance ProjectiveAction Normal3 where
-  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))
+-- | Synonym for "translateAfter4"
+translate4 :: Vec3 -> Proj4 -> Proj4
+translate4 = translateAfter4
 
--- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,
--- and the bottom-right corner is 1.
-invertProj :: Mat4 -> Mat4
-invertProj mat@(Mat4 u v w t) = 
-  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)
+-- | Equivalent to @\m -> m .*. translation v@.
+translateAfter4 :: Vec3 -> Proj4 -> Proj4
+translateAfter4 v p4 = 
+  toProjectiveUnsafe $
+    Mat4 r1 r2 r3 (extendWith 0 v &+ r4)
   where
-    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)
-    
+    Mat4 r1 r2 r3 r4 = fromProjective p4 
+
+-- | Equivalent to @\m -> translation v .*. m@.
+translateBefore4 :: Vec3 -> Proj4 -> Proj4
+translateBefore4 v p4 = 
+  toProjectiveUnsafe $ 
+    Mat4 r1 r2 r3 (extendWith 0 u &+ r4) 
+  where 
+   u = v .* (trim mat :: Mat3) 
+   mat@(Mat4 r1 r2 r3 r4) = fromProjective p4
+   
+---------------
+
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2008, Balazs Komuves
+Copyright (c) 2008-2009, Balazs Komuves
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -15,15 +15,15 @@
 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
+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
+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/src/flt/Base.hs b/src/flt/Base.hs
--- a/src/flt/Base.hs
+++ b/src/flt/Base.hs
@@ -1,10 +1,33 @@
 
-module Data.Vect.Flt.Base where
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, GeneralizedNewtypeDeriving #-}
 
+module Data.Vect.Flt.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 Control.Monad
 import System.Random  
 import Foreign
 
+--------------------------------------------------------------------------------
 -- class declarations
 
 class AbelianGroup g where
@@ -19,15 +42,17 @@
 vecSum :: AbelianGroup g => [g] -> g
 vecSum l = foldl (&+) zero l 
 
-class (AbelianGroup r) => 
-      Ring r where
+class MultSemiGroup r where
   (.*.) :: r -> r -> r
   one   :: r
 
+class (AbelianGroup r, MultSemiGroup r) => Ring r 
+
 infixl 7 .*. 
 
-ringProduct :: Ring r => [r] -> r
-ringProduct l = foldl (.*.) one l
+-- 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
@@ -74,7 +99,7 @@
   lensqr = normsqr
   dotprod :: v -> v -> Flt
   normsqr v = (v &. v)  
-  norm = sqrt.lensqr
+  norm = sqrt . lensqr
   dotprod = (&.)
 
 infix 7 &.
@@ -111,26 +136,28 @@
   fromNormalRadius :: Flt -> u -> v
   fromNormalRadius t n = t *& fromNormal n 
 
--- | projects the first vector onto the direction of the second (unit) vector
+-- | Projects the first vector onto the direction of the second (unit) vector
 project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v
 project' what dir = projectUnsafe what (fromNormal dir)
 
--- | direction (second argument) is assumed to be a /unit/ vector!
+-- | Direction (second argument) is assumed to be a /unit/ vector!
 projectUnsafe :: (Vector v, DotProd v) => v -> v -> v
 projectUnsafe what dir = what &- dir &* (what &. dir)
 
 project :: (Vector v, DotProd v) => v -> v -> v
 project what dir = what &- dir &* ((what &. dir) / (dir &. dir))
 
--- | since unit vectors are not a group, we need a separate function.
+-- | Since unit vectors are not a group, we need a separate function.
 flipNormal :: UnitVector v n => n -> n 
 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
@@ -165,13 +192,50 @@
 "inverse is an involution"    forall m. inverse (inverse m) = m
   #-}
   
+class Matrix m => Orthogonal m o | m->o, o->m where  
+  fromOrtho     :: o -> m 
+  toOrthoUnsafe :: m -> o
+  
+class (AbelianGroup m, Matrix m) => MatrixNorms m where
+  frobeniusNorm  :: m -> Flt       -- ^ the frobenius norm (= euclidean norm in the space of matrices)
+  matrixDistance :: m -> m -> Flt  -- ^ euclidean distance in the space of matrices
+  operatorNorm   :: m -> Flt       -- ^ (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 m where
   det :: m -> Flt    
-    
+
+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 v, UnitVector v u, Matrix m, Vector m, Tensor m v) => u -> m
+householder u = idmtx &- (2 *& outer v v) 
+  where v = fromNormal u
+
+householderOrtho :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v, Orthogonal m o) => u -> o
+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 v, Orthogonal n o, Diagonal v n) => Projective 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 -> m
+  toProjectiveUnsafe :: m -> p
+  orthogonal         :: o -> p
+  linear             :: n -> p
+  translation        :: v -> p
+  scaling            :: v -> p
+
+--------------------------------------------------------------------------------
 -- Vec / Mat datatypes
  
 data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt 
@@ -181,16 +245,16 @@
 data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt 
   deriving (Read,Show)
 
--- | these are /row/ vectors 
+-- | The components are /row/ vectors 
 data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show)
 data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show)
 data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)
 
 -- | The assumption when dealing with these is always that they are of unit length.
 -- Also, interpolation works differently.
-newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) 
-newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) 
-newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) 
+newtype Normal2 = Normal2 Vec2 deriving (Read,Show,Storable,DotProd,Dimension) 
+newtype Normal3 = Normal3 Vec3 deriving (Read,Show,Storable,DotProd,Dimension) 
+newtype Normal4 = Normal4 Vec4 deriving (Read,Show,Storable,DotProd,Dimension) 
 
 mkVec2 :: (Flt,Flt) -> Vec2
 mkVec3 :: (Flt,Flt,Flt) -> Vec3
@@ -200,6 +264,19 @@
 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 = Ortho2 Mat2 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+newtype Ortho3 = Ortho3 Mat3 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+newtype Ortho4 = Ortho4 Mat4 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)
+
+-- | Projective matrices, encoding affine transformations in dimension one less.
+newtype Proj3 = Proj3 Mat3 deriving (Read,Show,Storable,MultSemiGroup)
+newtype Proj4 = Proj4 Mat4 deriving (Read,Show,Storable,MultSemiGroup)
+
+--------------------------------------------------------------------------------
 -- Unit vectors
   
 instance UnitVector Vec2 Normal2 where
@@ -217,47 +294,163 @@
   fromNormal (Normal4 v) = v 
   toNormalUnsafe = Normal4
 
-rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)
-rndUnit g = 
+_rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)
+_rndUnit g = 
   if d > 0.01
     then ( v &* (1.0/d) , h )
-    else rndUnit h
+    else _rndUnit h
   where
     (v,h) = random g
     d = norm v
     
 instance Random Normal2 where
-  random g = let (v,h) = rndUnit g in (Normal2 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal2 v, h)  
   randomR _ = random
 
 instance Random Normal3 where
-  random g = let (v,h) = rndUnit g in (Normal3 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal3 v, h)  
   randomR _ = random
 
 instance Random Normal4 where
-  random g = let (v,h) = rndUnit g in (Normal4 v, h)  
+  random g = let (v,h) = _rndUnit g in (Normal4 v, h)  
   randomR _ = random
 
-{-
-instance Storable Normal2 where
-  alignment _ = alignment (undefined::Vec2)
-  sizeOf    _ = sizeOf    (undefined::Vec2)
-  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)
-  poke p (Normal2 v) = poke (castPtr p) v  
- 
-instance Storable Normal3 where
-  alignment _ = alignment (undefined::Vec3)
-  sizeOf    _ = sizeOf    (undefined::Vec3)
-  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)
-  poke p (Normal3 v) = poke (castPtr p) v  
+instance CrossProd Normal3 where
+  crossprod (Normal3 v) (Normal3 w) = mkNormal (crossprod v w)
 
-instance Storable Normal4 where
-  alignment _ = alignment (undefined::Vec4)
-  sizeOf    _ = sizeOf    (undefined::Vec4)
-  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)
-  poke p (Normal4 v) = poke (castPtr p) v  
--}
+--------------------------------------------------------------------------------
+-- Orthogonal matrices
 
+instance Orthogonal Mat2 Ortho2 where
+  fromOrtho (Ortho2 o) = o
+  toOrthoUnsafe = Ortho2
+
+instance Orthogonal Mat3 Ortho3 where
+  fromOrtho (Ortho3 o) = o
+  toOrthoUnsafe = Ortho3 
+
+instance Orthogonal Mat4 Ortho4 where
+  fromOrtho (Ortho4 o) = o
+  toOrthoUnsafe = Ortho4
+
+------
+
+instance Matrix Ortho2 where
+  transpose (Ortho2 o) = Ortho2 (transpose o)
+  idmtx = Ortho2 idmtx
+  inverse = transpose
+
+instance Matrix Ortho3 where
+  transpose (Ortho3 o) = Ortho3 (transpose o)
+  idmtx = Ortho3 idmtx
+  inverse = transpose
+
+instance Matrix Ortho4 where
+  transpose (Ortho4 o) = Ortho4 (transpose o)
+  idmtx = Ortho4 idmtx
+  inverse = transpose
+
+------
+
+instance Random Ortho2 where
+  random g = let (o,h) = _rndOrtho2 g in (toOrthoUnsafe (_flip1stRow2 o), h)
+  randomR _ = random
+
+instance Random Ortho3 where
+  random g = let (o,h) = _rndOrtho3 g in (toOrthoUnsafe (             o), h)
+  randomR _ = random
+
+instance Random Ortho4 where
+  random g = let (o,h) = _rndOrtho4 g in (toOrthoUnsafe (_flip1stRow4 o), h)
+  randomR _ = random
+
+------
+
+-- determinant will be -1
+_rndOrtho2 :: RandomGen g => g -> (Mat2, g)
+_rndOrtho2 g = (h2, g1) where
+  h2 = householder u2 :: Mat2 
+  (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 :: RandomGen g => g -> (Mat3, g) 
+_rndOrtho3 g = ( (h3 .*. m3), g2) where
+  m3 = (extendWith :: Flt -> Mat2 -> Mat3) 1 o2 
+  h3 = householder u3 :: Mat3
+  (u3,g1) = random g
+  (o2,g2) = _rndOrtho2 g1
+
+-- determinant will be -1
+_rndOrtho4 :: RandomGen g => g -> (Mat4, g) 
+_rndOrtho4 g = ( (h4 .*. m4), g2) where
+  m4 = (extendWith :: Flt -> Mat3 -> Mat4) 1 o3 
+  h4 = householder u4 :: Mat4
+  (u4,g1) = random g
+  (o3,g2) = _rndOrtho3 g1
+
+------
+
+_flip1stRow2 :: Mat2 -> Mat2
+_flip1stRow2 (Mat2 a b) = Mat2 (neg a) b
+
+_flip1stRow3 :: Mat3 -> Mat3
+_flip1stRow3 (Mat3 a b c) = Mat3 (neg a) b c
+
+_flip1stRow4 :: Mat4 -> Mat4
+_flip1stRow4 (Mat4 a b c d) = Mat4 (neg a) b c d
+
+--------------------------------------------------------------------------------
+-- projective matrices
+  
+instance Projective Vec2 Mat2 Ortho2 Mat3 Proj3 where
+  fromProjective (Proj3 m) = m
+  toProjectiveUnsafe = Proj3
+  orthogonal = Proj3 . extendWith 1 . fromOrtho
+  linear     = Proj3 . extendWith 1
+  translation v = Proj3 $ Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (extendWith 1 v)
+  scaling     v = Proj3 $ diag (extendWith 1 v)
+  
+instance Projective Vec3 Mat3 Ortho3 Mat4 Proj4 where
+  fromProjective (Proj4 m) = m
+  toProjectiveUnsafe = Proj4
+  orthogonal = Proj4 . extendWith 1 . fromOrtho 
+  linear     = Proj4 . extendWith 1
+  translation v = Proj4 $ Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (extendWith 1 v)
+  scaling     v = Proj4 $ diag (extendWith 1 v)
+
+instance Matrix Proj3 where
+  idmtx = Proj3 idmtx
+  transpose (Proj3 m) = Proj3 (transpose m)
+  inverse = _invertProj3
+
+instance Matrix Proj4 where
+  idmtx = Proj4 idmtx
+  transpose (Proj4 m) = Proj4 (transpose m)
+  inverse = _invertProj4
+
+_invertProj3 :: Proj3 -> Proj3
+_invertProj3 (Proj3 mat@(Mat3 _ _ t)) = 
+  Proj3 $ Mat3 (extendZero a) (extendZero b) (extendWith 1 t') 
+  where
+    t' = neg $ (trim t :: Vec2) .* invm2 
+    invm2@(Mat2 a b) = inverse $ (trim mat :: Mat2)
+
+-- Inverts a projective 4x4 matrix. But you can simply use "inverse" instead.
+-- We assume that the bottom-right corner is 1.
+_invertProj4 :: Proj4 -> Proj4
+_invertProj4 (Proj4 mat@(Mat4 _ _ _ t)) = 
+  Proj4 $ Mat4 (extendZero a) (extendZero b) (extendZero c) (extendWith 1 t') 
+  where
+    t' = neg $ (trim t :: Vec3) .* invm3 
+    invm3@(Mat3 a b c) = inverse $ (trim mat :: Mat3)
+
+--------------------------------------------------------------------------------
 -- Vec2 instances
 
 instance HasCoordinates Vec2 Flt where
@@ -313,7 +506,10 @@
         k = sizeOf (undefined::Flt)
     poke        p   x
     pokeByteOff p k y
-               
+
+instance Dimension Vec2 where dim _ = 2
+
+--------------------------------------------------------------------------------                    
 -- Mat2 instances
 
 instance HasCoordinates Mat2 Vec2 where
@@ -335,19 +531,21 @@
   (&+) (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   -- (zero::Vec2) (zero::Vec2)
-
+  zero = Mat2 zero zero  
+  
 instance Vector 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 Ring Mat2 where
+instance MultSemiGroup Mat2 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 Ring Mat2
+
 instance LeftModule Mat2 Vec2 where
   lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) 
   
@@ -361,11 +559,6 @@
   outer (Vec2 a b) (Vec2 x y) = Mat2
     (Vec2 (a*x) (a*y))
     (Vec2 (b*x) (b*y))
-{-
-  outer v w = 
-    let full = Mat2 (Vec2 1 1) (Vec2 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
 instance Determinant Mat2 where
   det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c 
@@ -392,6 +585,24 @@
     poke        p   r1
     pokeByteOff p k r2
 
+instance Random Mat2 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 Dimension Mat2 where dim _ = 2
+     
+instance MatrixNorms Mat2 where 
+  frobeniusNorm (Mat2 r1 r2) =  
+    sqrt $
+      normsqr r1 + 
+      normsqr r2
+     
+--------------------------------------------------------------------------------     
 -- Vec3 instances
 
 instance HasCoordinates Vec3 Flt where
@@ -453,7 +664,10 @@
     poke        p       x
     pokeByteOff p (k  ) y
     pokeByteOff p (k+k) z
-   
+
+instance Dimension Vec3 where dim _ = 3
+
+--------------------------------------------------------------------------------   
 -- Mat3 instances
 
 instance HasCoordinates Mat3 Vec3 where
@@ -494,13 +708,13 @@
   (&+) (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   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)
+  zero = Mat3 zero zero zero 
 
 instance Vector 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 Ring Mat3 where
+instance MultSemiGroup Mat3 where
   (.*.) (Mat3 r1 r2 r3) n = 
     let (Mat3 c1 c2 c3) = transpose n
     in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))
@@ -508,6 +722,8 @@
             (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))
   one = idmtx 
 
+instance Ring Mat3
+
 instance LeftModule Mat3 Vec3 where
   lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)
   
@@ -522,11 +738,6 @@
     (Vec3 (a*x) (a*y) (a*z))
     (Vec3 (b*x) (b*y) (b*z))
     (Vec3 (c*x) (c*y) (c*z))
-{-
-  outer v w = 
-    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
 instance Determinant Mat3 where
   det (Mat3 r1 r2 r3) = det (r1,r2,r3)
@@ -555,6 +766,26 @@
     pokeByteOff p (k  ) r2
     pokeByteOff p (k+k) r3
 
+instance Random Mat3 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 Dimension Mat3 where dim _ = 3
+  
+instance MatrixNorms Mat3 where 
+  frobeniusNorm (Mat3 r1 r2 r3)  = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 
+    
+--------------------------------------------------------------------------------
 -- Vec4 instances
 
 instance HasCoordinates Vec4 Flt where
@@ -614,6 +845,9 @@
     pokeByteOff p (k+k) z
     pokeByteOff p (3*k) w
 
+instance Dimension Vec4 where dim _ = 4
+
+--------------------------------------------------------------------------------
 -- Mat4 instances
 
 instance HasCoordinates Mat4 Vec4 where
@@ -641,7 +875,7 @@
   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 Ring Mat4 where
+instance MultSemiGroup Mat4 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))
@@ -650,6 +884,8 @@
             (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))
   one = idmtx 
 
+instance Ring Mat4
+
 instance LeftModule Mat4 Vec4 where
   lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)
   
@@ -665,14 +901,10 @@
     (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))
-{-
-  outer v w = 
-    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)
-    in  (diag v) .*. full .*. (diag w)
--}
 
---instance Determinant Mat4 where
---  det (Mat4 r1 r2 r3 r4)
+instance Determinant Mat4 where
+  det = error "det/Mat4: not implemented yet" 
+  -- det (Mat4 r1 r2 r3 r4) = 
 
 {-
 instance Show Mat4 where
@@ -700,6 +932,28 @@
     pokeByteOff p (k+k) r3
     pokeByteOff p (3*k) r4
 
+instance Random Mat4 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 Dimension Mat4 where dim _ = 4
+   
+instance MatrixNorms Mat4 where 
+  frobeniusNorm (Mat4 r1 r2 r3 r4) = 
+    sqrt $
+      normsqr r1 + 
+      normsqr r2 + 
+      normsqr r3 + 
+      normsqr r4  
+    
+--------------------------------------------------------------------------------
 -- Extend instances
 
 instance Extend Vec2 Vec3 where
@@ -718,17 +972,20 @@
   trim (Vec4 x y z _)       = Vec3 x y z
 
 instance Extend Mat2 Mat3 where
-  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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 Extend Mat2 Mat4 where
-  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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 Extend Mat3 Mat4 where
-  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero
-  extendWith _ _ = error "extendWith is meaningless for matrices"
+  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/flt/GramSchmidt.hs b/src/flt/GramSchmidt.hs
--- a/src/flt/GramSchmidt.hs
+++ b/src/flt/GramSchmidt.hs
@@ -9,7 +9,7 @@
 
 import Data.Vect.Flt.Base
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
 liftPair :: (a -> b) -> (a,a) -> (b,b)
 liftPair f (x,y) = (f x, f y)
@@ -20,7 +20,7 @@
 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
@@ -32,7 +32,7 @@
 "gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a 
   #-}
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
 instance GramSchmidt (Vec2,Vec2) where
   gramSchmidt = gramSchmidtPair
diff --git a/src/flt/OpenGL.hs b/src/flt/OpenGL.hs
--- a/src/flt/OpenGL.hs
+++ b/src/flt/OpenGL.hs
@@ -2,21 +2,92 @@
 -- TODO: the pointer versions of these functions should be really implemented 
 -- via the pointer versions of the original opengl functions...
 
--- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.
+-- | OpenGL support, including 'Vertex', 'TexCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.
  
 module Data.Vect.Flt.OpenGL where
 
 import Control.Monad
 import Data.Vect.Flt.Base
+import Data.Vect.Flt.Util.Projective
 import qualified Graphics.Rendering.OpenGL as GL
 
 import Foreign
 
-import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)
+import Graphics.Rendering.OpenGL hiding 
+  ( Normal3 , rotate , translate , scale
+  , matrix , currentMatrix , withMatrix , multMatrix 
+  )
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
-{-# SPECIALISE radianToDegrees :: Float -> Float #-}
+-- | There should be a big warning here about the different conventions, 
+-- hidden transpositions, and all the confusion this will inevitably cause...
+--
+-- As it stands, 
+--
+-- > glRotate t1 axis1 >> glRotate t2 axis2 >> glRotate t3 axis3
+-- 
+-- has the same result as
+--
+-- > multMatrix (rotMatrixProj4 t3 axis3 .*. rotMatrixProj4 t2 axis2 .*. rotMatrixProj4 t1 axis1)
+--
+-- because at the interface of OpenGL and this library there is a transposition
+-- to compensate for the different conventions. (This transposition is implicit
+-- in the code, because the way the matrices are stored in the memory is also
+-- different: OpenGL stores them column-major, and we store them row-major).
+
+class ToOpenGLMatrix m where
+  makeGLMatrix :: m -> IO (GLmatrix Flt)
+
+class FromOpenGLMatrix m where
+  peekGLMatrix :: GLmatrix Flt -> IO m
+  
+setMatrix :: ToOpenGLMatrix m => Maybe MatrixMode -> m -> IO ()
+setMatrix mode m = makeGLMatrix m >>= \x -> GL.matrix mode $= x
+ 
+getMatrix :: FromOpenGLMatrix m => Maybe MatrixMode -> IO m
+getMatrix mode = get (GL.matrix mode) >>= peekGLMatrix
+
+matrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => Maybe MatrixMode -> StateVar m
+matrix mode = makeStateVar (getMatrix mode) (setMatrix mode)
+
+currentMatrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => StateVar m
+currentMatrix = matrix Nothing
+
+multMatrix :: ToOpenGLMatrix m => m -> IO ()
+multMatrix m = makeGLMatrix m >>= GL.multMatrix
+
+instance ToOpenGLMatrix Mat4 where
+  makeGLMatrix m = GL.withNewMatrix GL.ColumnMajor (flip poke m . castPtr) 
+ 
+instance FromOpenGLMatrix Mat4 where
+  -- huh? GL.withMatrix is strange
+  peekGLMatrix x = GL.withMatrix x $ \_ p -> peek (castPtr p)
+  
+instance ToOpenGLMatrix Mat3 where
+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)
+ 
+instance ToOpenGLMatrix Mat2 where
+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)
+
+instance ToOpenGLMatrix Ortho4 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat4)
+
+instance ToOpenGLMatrix Ortho3 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat3)
+
+instance ToOpenGLMatrix Ortho2 where
+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat2)
+
+instance ToOpenGLMatrix Proj4 where
+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat4)
+
+instance ToOpenGLMatrix Proj3 where
+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat3)
+  
+--------------------------------------------------------------------------------
+
+{-# SPECIALISE radianToDegrees :: Float  -> Float  #-}
 {-# SPECIALISE radianToDegrees :: Double -> Double #-}
 radianToDegrees :: RealFrac a => a -> a
 radianToDegrees x = x * 57.295779513082322
@@ -27,20 +98,61 @@
 degreesToRadian x = x * 1.7453292519943295e-2
 
 -- | The angle is in radians. (WARNING: OpenGL uses degrees!)
-rotate :: Flt -> Vec3 -> IO ()
-rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)
+glRotate :: Flt -> Vec3 -> IO ()
+glRotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)
 
-translate :: Vec3 -> IO ()
-translate (Vec3 x y z) = GL.translate (Vector3 x y z)
+glTranslate :: Vec3 -> IO ()
+glTranslate (Vec3 x y z) = GL.translate (Vector3 x y z)
 
-scale3 :: Vec3 -> IO ()
-scale3 (Vec3 x y z) = GL.scale x y z
+glScale3 :: Vec3 -> IO ()
+glScale3 (Vec3 x y z) = GL.scale x y z
 
-scale :: Flt -> IO ()
-scale x = GL.scale x x x
+glScale :: Flt -> IO ()
+glScale x = GL.scale x x x
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
+ 
+-- | \"Orthogonal projecton\" matrix, a la OpenGL 
+-- (the corresponding functionality is removed in OpenGL 3.1)
+orthoMatrix 
+  :: (Flt,Flt)   -- ^ (left,right)
+  -> (Flt,Flt)   -- ^ (bottom,top)
+  -> (Flt,Flt)   -- ^ (near,far)
+  -> Mat4 
+orthoMatrix (l,r) (b,t) (n,f) = Mat4
+  (Vec4 (2/(r-l)) 0 0 0)
+  (Vec4 0 (2/(t-b)) 0 0)
+  (Vec4 0 0 (-2/(f-n)) 0)
+  (Vec4 (-(r+l)/(r-l)) (-(t+b)/(t-b)) (-(f+n)/(f-n)) 1)
+  
+-- | The same as "orthoMatrix", but with a different parametrization.
+orthoMatrix2 {- ' CPP is sensitive to primes -}
+  :: Vec3     -- ^ (left,top,near)
+  -> Vec3     -- ^ (right,bottom,far)
+  -> Mat4 
+orthoMatrix2 (Vec3 l t n) (Vec3 r b f) = orthoMatrix (l,r) (b,t) (n,f)
 
+-- | \"Perspective projecton\" matrix, a la OpenGL 
+-- (the corresponding functionality is removed in OpenGL 3.1).
+frustumMatrix
+  :: (Flt,Flt)   -- ^ (left,right)
+  -> (Flt,Flt)   -- ^ (bottom,top)
+  -> (Flt,Flt)   -- ^ (near,far)
+  -> Mat4 
+frustumMatrix (l,r) (b,t) (n,f) = Mat4
+  (Vec4 (2*n/(r-l)) 0 0 0)
+  (Vec4 0 (2*n/(t-b)) 0 0)
+  (Vec4 ((r+l)/(r-l)) ((t+b)/(t-b)) (-(f+n)/(f-n)) (-1))
+  (Vec4 0 0 (-2*f*n*(f-n)) 0)
+  
+-- | The same as "frustumMatrix", but with a different parametrization.
+frustumMatrix2 {- ' CPP is sensitive to primes -}
+  :: Vec3     -- ^ (left,top,near)
+  -> Vec3     -- ^ (right,bottom,far)
+  -> Mat4 
+frustumMatrix2 (Vec3 l t n) (Vec3 r b f) = frustumMatrix (l,r) (b,t) (n,f)
+
+--------------------------------------------------------------------------------
 -- Vertex instances
 
 instance GL.Vertex Vec2 where
@@ -55,17 +167,20 @@
   vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)
   vertexv p = peek p >>= vertex   
 
--------------------------------------------------------
-
+--------------------------------------------------------------------------------
 -- the Normal instance
 -- note that there is no Normal2\/Normal4 in the OpenGL binding
 
 instance GL.Normal Normal3 where
-  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)
+  normal u = GL.normal (GL.Normal3 x y z) 
+    where Vec3 x y z = fromNormal u 
   normalv p = peek p >>= normal 
 
--------------------------------------------------------
+instance GL.Normal Vec3 where
+  normal (Vec3 x y z) = GL.normal (GL.Normal3 x y z) 
+  normalv p = peek p >>= normal 
 
+--------------------------------------------------------------------------------
 -- Color instances
   
 instance GL.Color Vec3 where
@@ -87,8 +202,7 @@
   secondaryColorv p = peek p >>= secondaryColor
 -}
 
--------------------------------------------------------
-
+--------------------------------------------------------------------------------
 -- TexCoord instances
 
 instance GL.TexCoord Vec2 where
@@ -109,8 +223,7 @@
   multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)
   multiTexCoordv unit p = peek p >>= multiTexCoord unit
 
--------------------------------------------------------
-    
+--------------------------------------------------------------------------------
 -- Vertex Attributes (experimental)
 
 class VertexAttrib' a where
@@ -129,16 +242,18 @@
   vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w 
 
 instance VertexAttrib' Normal2 where
-  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y
+  vertexAttrib loc u = GL.vertexAttrib2 loc x y
+    where Vec2 x y = fromNormal u 
 
 instance VertexAttrib' Normal3 where
-  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z
+  vertexAttrib loc u = GL.vertexAttrib3 loc x y z
+    where Vec3 x y z = fromNormal u 
 
 instance VertexAttrib' Normal4 where
-  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w
-
--------------------------------------------------------
-
+  vertexAttrib loc u = GL.vertexAttrib4 loc x y z w
+    where Vec4 x y z w = fromNormal u 
+   
+--------------------------------------------------------------------------------
 -- Uniform (again, experimental)
 
 -- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; 
@@ -180,3 +295,5 @@
   uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)
     
 #endif
+
+
diff --git a/src/flt/Util/Dim2.hs b/src/flt/Util/Dim2.hs
--- a/src/flt/Util/Dim2.hs
+++ b/src/flt/Util/Dim2.hs
@@ -3,12 +3,13 @@
 
 import Data.Vect.Flt.Base
 
--- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.
+-- | Example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.
 structVec2 :: [Flt] -> [Vec2]
 structVec2 [] = []
 structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls) 
 structVec2 _ = error "structVec2"
 
+-- | The opposite of "structVec2".
 destructVec2 :: [Vec2] -> [Flt]
 destructVec2 [] = []
 destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  
@@ -48,10 +49,13 @@
 angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt
 angle2' = angle2 . fromNormal
 
--- |Rotation matrix by a given angle (in radians), counterclockwise.
+-- | Rotation matrix by a given angle (in radians), counterclockwise.
 rotMatrix2 :: Flt -> Mat2
 rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a
 
+rotMatrixOrtho2 :: Flt -> Ortho2
+rotMatrixOrtho2 = toOrthoUnsafe . rotMatrix2
+
 rotate2 :: Flt -> Vec2 -> Vec2
 rotate2 a v = v .* (rotMatrix2 a) 
 
@@ -62,3 +66,5 @@
 -- |Rotates clockwise by 90 degrees.
 rotateCW :: Vec2 -> Vec2
 rotateCW (Vec2 x y) = Vec2 y (-x)
+
+
diff --git a/src/flt/Util/Dim3.hs b/src/flt/Util/Dim3.hs
--- a/src/flt/Util/Dim3.hs
+++ b/src/flt/Util/Dim3.hs
@@ -3,18 +3,26 @@
 
 import Data.Vect.Flt.Base
 
+--------------------------------------------------------------------------------
+
+-- | Example: @structVec3 [1,2,3,4,5,6] = [ Vec3 1 2 3 , Vec3 4 5 6]@.
 structVec3 :: [Flt] -> [Vec3]
 structVec3 [] = []
 structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) 
 structVec3 _ = error "structVec3"
 
+-- | The opposite of "structVec3".
 destructVec3 :: [Vec3] -> [Flt]
 destructVec3 [] = []
 destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)  
 
+--------------------------------------------------------------------------------
+
 det3 :: Vec3 -> Vec3 -> Vec3 -> Flt
 det3 u v w = det (u,v,w)
 
+--------------------------------------------------------------------------------
+
 translate3X :: Flt -> Vec3 -> Vec3
 translate3Y :: Flt -> Vec3 -> Vec3
 translate3Z :: Flt -> Vec3 -> Vec3
@@ -42,26 +50,34 @@
 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 -} Flt       -- ^ angle (in radians)
-         -> Normal3   -- ^ axis (should be a /unit/ vector!) 
-         -> Vec3      -- ^ vector
-         -> Vec3      -- ^ result
+--------------------------------------------------------------------------------
+
+rotate3' {- ' CPP is sensitive to primes -} 
+  :: Flt       -- ^ angle (in radians)
+  -> Normal3   -- ^ axis (should be a /unit/ vector!) 
+  -> Vec3      -- ^ vector
+  -> Vec3      -- ^ result
 rotate3' angle axis v = v .* (rotMatrix3' axis angle)
 
-rotate3 :: Flt    -- ^ angle (in radians)
-        -> Vec3   -- ^ axis (arbitrary nonzero vector)
-        -> Vec3   -- ^ vector
-        -> Vec3   -- ^ result
+rotate3 
+  :: Flt    -- ^ angle (in radians)
+  -> Vec3   -- ^ axis (arbitrary nonzero vector)
+  -> Vec3   -- ^ vector
+  -> Vec3   -- ^ 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/. 
+-- | Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
 rotMatrix3 :: Vec3 -> Flt -> Mat3
 rotMatrix3 v a = rotMatrix3' (mkNormal v) a
 
--- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. 
+rotMatrixOrtho3 :: Vec3 -> Flt -> Ortho3
+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 -} Normal3 -> Flt -> Mat3
-rotMatrix3' (Normal3 v) a = 
-  let c = cos a
+rotMatrix3' u a = 
+  let v = fromNormal u
+      c = cos a
       s = sin a
       m1 = scalarMul (1-c) (outer v v)
       x = _1 v
@@ -72,4 +88,65 @@
                 (Vec3 ( s*y) (-s*x)   c   )
   in (m1 &+ m2)
 
+rotMatrixOrtho3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Ortho3
+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 :: Normal3 -> Vec3 -> Vec3
+reflect u v = (s *& n) &- v where 
+  n = fromNormal u
+  s = 2 * (n &. v)
+
+reflect' :: Normal3 -> Normal3 -> Normal3
+reflect' u x = toNormalUnsafe $ reflect u (fromNormal x)
+  
+refract :: Flt -> Normal3 -> Vec3 -> Vec3
+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 -} 
+  :: Flt -> Normal3 -> Normal3 -> Normal3
+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 -}  
+  :: Flt -> Normal3 -> Normal3 -> Maybe Normal3
+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/flt/Util/Dim4.hs b/src/flt/Util/Dim4.hs
--- a/src/flt/Util/Dim4.hs
+++ b/src/flt/Util/Dim4.hs
@@ -39,7 +39,7 @@
 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))@.
@@ -75,8 +75,9 @@
 -- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors.
 -- Intended for multiplication on the /right/!
 rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4
-rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 
+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
diff --git a/src/flt/Util/Projective.hs b/src/flt/Util/Projective.hs
--- a/src/flt/Util/Projective.hs
+++ b/src/flt/Util/Projective.hs
@@ -1,5 +1,6 @@
 
--- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on
+-- | 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
 --
 -- >     _____
@@ -12,6 +13,8 @@
 -- 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.Flt.Util.Projective where
 
@@ -20,68 +23,73 @@
 
 import qualified Data.Vect.Flt.Util.Dim4 as Dim4
 
-class ExtendProjective v e | v->e where
-  extendProj     :: v -> e
-  extendProjWith :: Flt -> v -> e
-  extendProj = extendProjWith 1
-  
-instance ExtendProjective Vec2 Vec4 where
-  extendProj       (Vec2 x y) = Vec4 x y 0 1
-  extendProjWith w (Vec2 x y) = Vec4 x y 0 w
-  
-instance ExtendProjective Vec3 Vec4 where
-  extendProj       (Vec3 x y z) = Vec4 x y z 1
-  extendProjWith w (Vec3 x y z) = Vec4 x y z w
-
-instance ExtendProjective Vec4 Vec4 where
-  extendProj = id
-  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w
+--------------------------------------------------------------------------------
 
-instance ExtendProjective Mat2 Mat4 where
-  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)
-  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)
+rotMatrixProj4' :: {- ' CPP is sensitive to primes -}  Flt -> Normal3 -> Proj4
+rotMatrixProj4' angle axis = linear $ rotMatrix3' axis angle
 
-instance ExtendProjective Mat3 Mat4 where
-  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)
-  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)
+rotMatrixProj4 :: Flt -> Vec3 -> Proj4
+rotMatrixProj4 angle axis = linear $ rotMatrix3 axis angle
 
-rotMatrixProj :: Flt -> Normal3 -> Mat4
-rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle
+-- | synonym for "rotateAfterProj4"
+rotateProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateProj4 = rotateAfterProj4
 
-rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4
-rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle
+-- | Synonym for @\m -> m .*. rotMatrixProj4 angle axis@.
+rotateAfterProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateAfterProj4 angle axis m = m .*. (rotMatrixProj4' angle axis) 
 
-translMatrixProj :: Vec3 -> Mat4
-translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)
+-- | Synonym for @\m -> rotMatrixProj4 angle axis .*. m@.
+rotateBeforeProj4 :: Flt -> Normal3 -> Proj4 -> Proj4
+rotateBeforeProj4 angle axis m = (rotMatrixProj4' angle axis) .*. m 
 
--- | we assume that the bottom-right corner is 1.
-translWithProj :: Vec3 -> Mat4 -> Mat4
-translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)
+---------------
 
-scaleMatrixProj :: Vec3 -> Mat4
-scaleMatrixProj v = diag $ extendProj v
+--scalingUniformProj3 :: Flt -> Proj3
+--scalingUniformProj3 x = scaling (Vec2 x x)
 
-scaleMatrixUniformProj :: Flt -> Mat4
-scaleMatrixUniformProj s = diag (Vec4 s s s 1)
+scalingUniformProj4 :: Flt -> Proj4
+scalingUniformProj4 x = scaling (Vec3 x x x)
 
-class ProjectiveAction v where
-  actProj :: v -> Mat4 -> v
- 
-instance ProjectiveAction Vec3 where
-  actProj v m = trim $ (extendProj v) .* m 
+-- | Equivalent to @\m -> scaling v .*. m@.
+scaleBeforeProj4 :: Vec3 -> Proj4 -> Proj4
+scaleBeforeProj4 (Vec3 u v w) p4 =  
+  toProjectiveUnsafe $ 
+    Mat4 (u*&a) (v*&b) (w*&c) t
+  where
+    Mat4 a b c t = fromProjective p4
 
-instance ProjectiveAction Vec4 where
-  actProj v m = v .* m 
+-- | Equivalent to @\m -> m .*. scaling v@.
+scaleAfterProj4 :: Vec3 -> Proj4 -> Proj4
+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
+    
+---------------
 
--- | When acting on unit vectors, we ignore the translation part.
-instance ProjectiveAction Normal3 where
-  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))
+-- | Synonym for "translateAfter4"
+translate4 :: Vec3 -> Proj4 -> Proj4
+translate4 = translateAfter4
 
--- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,
--- and the bottom-right corner is 1.
-invertProj :: Mat4 -> Mat4
-invertProj mat@(Mat4 u v w t) = 
-  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)
+-- | Equivalent to @\m -> m .*. translation v@.
+translateAfter4 :: Vec3 -> Proj4 -> Proj4
+translateAfter4 v p4 = 
+  toProjectiveUnsafe $
+    Mat4 r1 r2 r3 (extendWith 0 v &+ r4)
   where
-    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)
-    
+    Mat4 r1 r2 r3 r4 = fromProjective p4 
+
+-- | Equivalent to @\m -> translation v .*. m@.
+translateBefore4 :: Vec3 -> Proj4 -> Proj4
+translateBefore4 v p4 = 
+  toProjectiveUnsafe $ 
+    Mat4 r1 r2 r3 (extendWith 0 u &+ r4) 
+  where 
+   u = v .* (trim mat :: Mat3) 
+   mat@(Mat4 r1 r2 r3 r4) = fromProjective p4
+   
+---------------
+
diff --git a/vect.cabal b/vect.cabal
--- a/vect.cabal
+++ b/vect.cabal
@@ -1,5 +1,5 @@
 Name:                vect
-Version:             0.4.0
+Version:             0.4.5
 Synopsis:            A low-dimensional linear algebra library, tailored to computer graphics.
 Description:         A low-dimensional (2, 3 and 4) linear algebra library, 
                      with lots of useful functions. Intended usage is primarily