diff --git a/Algebra/Classes.hs b/Algebra/Classes.hs
--- a/Algebra/Classes.hs
+++ b/Algebra/Classes.hs
@@ -1,8 +1,15 @@
+{-# LANGUAGE TupleSections #-}
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE GADTs #-}
 {-# LANGUAGE MultiParamTypeClasses, ConstraintKinds, FlexibleContexts, FlexibleInstances, DeriveGeneric #-}
 module Algebra.Classes where
 
-import Prelude as Algebra.Classes (Int,Integer,Float,Double, Foldable (..), (==), Monoid(..), Ord(..)
-                                  ,Real(..), Enum(..), snd, Rational, Functor(..), Eq(..), Bool(..), Semigroup(..))
+import Prelude (Int,Integer,Float,Double, (==), Monoid(..), Ord(..), Foldable,
+                foldMap, Char,
+                 Real(..), Enum(..), snd, Rational, Functor(..), Eq(..), Bool(..), Semigroup(..), Show(..), uncurry)
+
 import qualified Prelude
 import qualified Data.Ratio
 import qualified Data.Map.Strict as M
@@ -10,8 +17,13 @@
 import Foreign.C
 import Data.Word
 import Data.Binary
+import Data.Complex
 import GHC.Generics
+import Test.QuickCheck
+import Control.Applicative
 
+-- import Data.Functor.Utils ((#.))
+
 infixl 6 -
 infixl 6 +
 
@@ -20,6 +32,8 @@
 infixl 7 /
 infixl 7 `mod`
 infixl 7 `div`
+infixr 8 ^
+infixr 8 ^+
 
 type Natural = Integer
 
@@ -48,25 +62,62 @@
 instance Additive a => Multiplicative (Exponential a) where
   Exponential a * Exponential b = Exponential (a + b)
   one = Exponential zero
-  Exponential a ^ n = Exponential (times n a)
+  Exponential a ^+ n = Exponential (times n a)
 
 instance Group a => Division (Exponential a) where
   recip (Exponential a) = Exponential (negate a)
   Exponential a / Exponential b = Exponential (a - b)
 
+timesDefault :: (Additive a2, Prelude.Integral a1) => a1 -> a2 -> a2
+timesDefault n0 = if n0 < 0 then Prelude.error "Algebra.Classes.times: negative number of times" else go n0
+    where go 0 _ = zero
+          go n x = if r == 0 then y + y else x + y + y
+            where (m,r) = n `Prelude.divMod` 2
+                  y = go m x
+
 -- | Additive monoid
 class Additive a where
   (+) :: a -> a -> a
   zero :: a
   times :: Natural -> a -> a
-  times 0 _ = zero
-  times n x = if r == 0 then y + y else x + y + y
-    where (m,r) = n `Prelude.divMod` 2
-          y = times m x
+  times = timesDefault
 
-add :: (Foldable t, Additive a) => t a -> a
-add xs = fromSum (foldMap Sum xs)
+class (Arbitrary a, Show a) => TestEqual a where
+  (=.=) :: a -> a -> Property
 
+infix 0 =.=
+
+instance Multiplicative Property where
+  one = property True
+  (*) = (.&&.)
+
+nameLaw :: Testable prop => Prelude.String -> prop -> Property
+nameLaw x p = label x (counterexample x p)
+
+law_zero_plus :: forall a. (Additive a, TestEqual a) => a -> Property
+law_zero_plus n = nameLaw "zero/plus" (zero + n =.= n)
+
+law_plus_zero :: (Additive a, TestEqual a) => a -> Property
+law_plus_zero n = nameLaw "plus/zero" (n + zero =.= n)
+
+law_plus_assoc :: (Additive a, TestEqual a) => a -> a -> a -> Property
+law_plus_assoc m n o = nameLaw "plus/assoc" (n + (m + o) =.= (n + m) + o)
+
+law_times :: (TestEqual a, Additive a) => Positive Integer -> a -> Property
+law_times (Positive m) n = nameLaw "times" (times m n =.= timesDefault m n)
+
+laws_additive :: forall a. (Additive a, TestEqual a) => Property
+laws_additive = product [property (law_zero_plus @a)
+                        ,property (law_plus_zero @a)
+                        ,property (law_plus_assoc @a)
+                        ,property (law_times @a)]
+
+instance TestEqual Int where (=.=) = (===)
+
+
+sum :: (Foldable t, Additive a) => t a -> a
+sum xs = fromSum (foldMap Sum xs)
+
 instance Additive Integer where
   (+) = (Prelude.+)
   zero = 0
@@ -115,6 +166,10 @@
 class Additive r => DecidableZero r where
   isZero :: r -> Bool
 
+law_decidable_zero :: forall a. (DecidableZero a, TestEqual a) => Property
+law_decidable_zero = property (isZero (zero @a))
+
+
 instance DecidableZero Integer where
   isZero = (== 0)
 instance DecidableZero CInt where
@@ -136,6 +191,13 @@
 
 class Additive a => AbelianAdditive a
   -- just a law.
+
+law_plus_comm :: (TestEqual a, Additive a) => a -> a -> Property
+law_plus_comm m n = nameLaw "plus/comm" (m + n =.= n + m)
+
+laws_abelian_additive :: forall a. (Group a, TestEqual a) => Property
+laws_abelian_additive = laws_additive @a .&&. product [property (law_plus_comm @a)]
+
 instance AbelianAdditive Integer
 instance AbelianAdditive CInt
 instance AbelianAdditive Int
@@ -143,6 +205,9 @@
 instance AbelianAdditive Float
 instance (Ord k,AbelianAdditive v) => AbelianAdditive (Map k v)
 
+multDefault :: Group a => Natural -> a -> a
+multDefault n x = if n < 0 then negate (times (negate n) x) else times n x
+
 class Additive a => Group a where
   {-# MINIMAL (negate | (-)) #-}
   (-) :: a -> a -> a
@@ -150,8 +215,23 @@
   negate :: a -> a
   negate b = zero - b
   mult :: Integer -> a -> a
-  mult n x = if n < 0 then negate (times (negate n) x) else times n x
+  mult = multDefault
 
+law_negate_minus :: (TestEqual a, Group a) => a -> a -> Property
+law_negate_minus m n = nameLaw "minus/negate" (m + negate n =.= m - n)
+
+law_mult :: (TestEqual a, Group a) => Integer -> a -> Property
+law_mult m n = nameLaw "mult" (mult m n =.= multDefault m n)
+
+
+laws_group :: forall a. (Group a, TestEqual a) => Property
+laws_group = laws_additive @a .&&. product [property (law_negate_minus @a)
+                                           ,property (law_mult @a)]
+
+laws_abelian_group :: forall a. (Group a, TestEqual a) => Property
+laws_abelian_group = laws_group @a .&&. product [property (law_plus_comm @a)]
+
+
 instance Group Integer where
   (-) = (Prelude.-)
   negate = Prelude.negate
@@ -194,6 +274,36 @@
 class (AbelianAdditive a, PreRing scalar) => Module scalar a where
   (*^) :: scalar -> a -> a
 
+law_module_zero :: forall s a. (Module s a, TestEqual a) => s -> Property
+law_module_zero s = nameLaw "module/zero" (s *^ zero =.= zero @a)
+
+law_module_one :: forall s a. (Module s a, TestEqual a) => a -> Property
+law_module_one x = nameLaw "module/one" ((one @s) *^ x =.= x)
+
+law_module_sum :: forall s a. (Module s a, TestEqual a) => s -> a -> a -> Property
+law_module_sum s x y = nameLaw "module/distr/left" (s *^ (x + y) =.= s*^x + s *^ y)
+
+law_module_sum_left :: forall s a. (Module s a, TestEqual a) => s -> s -> a -> Property
+law_module_sum_left s t x = nameLaw "module/distr/right" ((s + t) *^ x =.= s*^x + t *^ x)
+
+law_module_mul :: forall s a. (Module s a, TestEqual a) => s -> s -> a -> Property
+law_module_mul s t x = nameLaw "module/mul/assoc" ((s * t) *^ x =.= s *^ t *^ x)
+
+laws_module :: forall s a. (Module s a, TestEqual a, Arbitrary s, Show s) => Property
+laws_module = laws_additive @a .&&. product [property (law_module_zero @s @a)
+                                            ,property (law_module_one @s @a)
+                                            ,property (law_module_sum @s @a)
+                                            ,property (law_module_sum_left @s @a)
+                                            ,property (law_module_mul @s @a)
+                                            ]
+
+-- Comparision of maps with absence of a key equivalent to zero value.
+instance (Ord x, Show x, Arbitrary x,TestEqual a,Additive a) => TestEqual (Map x a) where
+  x =.= y = product (uncurry (=.=) <$> M.unionWith collapse ((,zero) <$> x) ((zero,) <$> y))
+    where collapse :: (a,b) -> (c,d) -> (a,d)
+          collapse (a,_) (_,b) = (a,b)
+
+
 instance Module Integer Integer where
   (*^) = (*)
 
@@ -209,64 +319,67 @@
 instance Module Float Float where
   (*^) = (*)
 
-instance (Ord k, Module v v) => Module v (Map k v) where
-  s *^ m = fmap (s *) m
+instance (Ord k, Module a b) => Module a (Map k b) where
+  s *^ m = fmap (s *^) m
 
 -- | Multiplicative monoid
 class Multiplicative a where
   (*) :: a -> a -> a
   one :: a
-  (^) :: a -> Natural -> a
+  (^+) :: a -> Natural -> a
 
-  (^) _ 0 = one
-  (^) x n = if r == 0 then y * y else x * y * y
-    where (m,r) = n `Prelude.divMod` 2
-          y = (^) y m
+  x0 ^+ n0 = if n0 < 0 then Prelude.error "Algebra.Classes.^: negative exponent" else go x0 n0
+    where go _ 0 = one
+          go x n = if r == 0 then y * y else x * y * y
+            where (m,r) = n `Prelude.divMod` 2
+                  y = go x m
 
-multiply :: (Multiplicative a, Foldable f) => f a -> a
-multiply xs = fromProduct (foldMap Product xs)
 
+product :: (Multiplicative a, Foldable f) => f a -> a
+product xs = fromProduct (foldMap Product xs)
+
 instance Multiplicative Integer where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 instance Multiplicative CInt where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 instance Multiplicative Word32 where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 instance Multiplicative Word16 where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 instance Multiplicative Word8 where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 instance Multiplicative Int where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 instance Multiplicative Double where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 instance Multiplicative Float where
   (*) = (Prelude.*)
   one = 1
-  (^) = (Prelude.^)
+  (^+) = (Prelude.^)
 
 
+
 type SemiRing a = (Multiplicative a, AbelianAdditive a)
 type PreRing a = (SemiRing a, Group a)
 
@@ -300,11 +413,19 @@
   (/) :: a -> a -> a
   x / y           =  x * recip y
 
+  (^) :: a -> Integer -> a
+  b ^ n | n < 0 = recip b ^+ negate n
+        | True  = b ^+ n
+
 instance Division Double where
   (/) = (Prelude./)
+  recip = Prelude.recip
+  (^) = (Prelude.^^)
 
 instance Division Float where
   (/) = (Prelude./)
+  recip = Prelude.recip
+  (^) = (Prelude.^^)
 
 class (Ring a, Division a) => Field a where
   fromRational :: Rational -> a
@@ -317,7 +438,6 @@
 instance Field Float where
   fromRational = Prelude.fromRational
 
-type VectorSpace scalar a = (Field scalar, Module scalar a)
 
 class Ring a => EuclideanDomain a where
     {-# MINIMAL (stdUnit | normalize) , (divMod | (div , mod)) #-}
@@ -370,19 +490,7 @@
     rem       =  Prelude.rem
     toInteger = Prelude.toInteger
 
-{-
-Note: the following is not quite what we intuitively want, because
 
-class Field a => AlgebraicallyClosed  a where
-  sqrt :: a -> (a,a)
-
-AlgebraicallyClosed numbers have two square roots.
-
--}
-
-data Ratio a = !a :% !a  deriving (Eq)
-type MyRational = Ratio Integer
-
 gcd             :: (Integral a) => a -> a -> a
 {-# NOINLINE [1] gcd #-}
 gcd x y         =  gcd' (stdAssociate x) (stdAssociate y)
@@ -391,50 +499,12 @@
    gcd' a 0  =  a
    gcd' a b  =  gcd' b (a `rem` b)
 
-{-
--- | 'reduce' is a subsidiary function used only in this module.
--- It normalises a ratio by dividing both numerator and denominator by
--- their greatest common divisor.
-reduce :: (Eq a, Integral a) => a -> a -> Ratio a
-{-# SPECIALISE reduce :: Integer -> Integer -> MyRational #-}
-reduce _ 0              =  error "reduce: division by zero"
-reduce x y              =  (x `quot` d) :% (y `quot` d)
-                           where d = gcd x y
+{- -}
 
-(%) :: Integral a => a -> a -> Ratio a
-x % y                   =  reduce (x * stdUnit y) (stdAssociate y)
 
 
-instance Integral a => AbelianAdditive (Ratio a) where
-
-instance Integral a => Additive (Ratio a) where
-    (x:%y) + (x':%y')   =  reduce (x*y' + x'*y) (y*y')
-    zero = 0
-    times n (x :% y) = reduce (times n x) y
-
-instance Integral a => Multiplicative (Ratio a) where
-    (x:%y) * (x':%y')   =  reduce (x * x') (y * y')
-    one = 1 :% 1
-
-instance Integral a => Group (Ratio a) where
-    (x:%y) - (x':%y')   =  reduce (x*y' - x'*y) (y*y')
-    negate (x:%y)       =  (-x) :% y
-
-instance Integral a => EuclideanDomain (Ratio a) where
-    stdAssociate (x:%y)          =  stdAssociate x :% y
-    stdUnit (x:%_)       =  stdUnit x :% 1
-
-instance Integral a => Ring (Ratio a) where
-    fromInteger x       =  fromInteger x :% 1
-
-instance Integral a => Division (Ratio a) where
-  recip (x:%y) = y:%x
-
-instance Integral a => Field (Ratio a) where
-
-
--}
-
+--------------------------
+-- Ratio instances
 instance Prelude.Integral a => Additive (Data.Ratio.Ratio a) where
   zero = Prelude.fromInteger 0
   (+) = (Prelude.+)
@@ -448,31 +518,94 @@
 instance Prelude.Integral a => Multiplicative (Data.Ratio.Ratio a) where
   one = Prelude.fromInteger 1
   (*) = (Prelude.*)
+  (^+) = (Prelude.^)
 
 instance Prelude.Integral a => Division (Data.Ratio.Ratio a) where
   recip = Prelude.recip
   (/) = (Prelude./)
-
+  (^) = (Prelude.^^)
 instance Prelude.Integral a => Module (Data.Ratio.Ratio a) (Data.Ratio.Ratio a) where
   (*^) = (*)
-
 instance Prelude.Integral a => Ring (Data.Ratio.Ratio a) where
   fromInteger = Prelude.fromInteger
-
 instance Prelude.Integral a => Field (Data.Ratio.Ratio a) where
   fromRational = Prelude.fromRational
 
+
+----------------------
+-- Complex instances
+instance Module Rational Double where
+    r *^ d = fromRational r * d
+instance Additive a => Additive (Complex a) where
+    (x:+y) + (x':+y')   =  (x+x') :+ (y+y')
+    zero = zero :+ zero
+instance Ring a => Multiplicative (Complex a) where
+    (x:+y) * (x':+y')   =  (x*x'-y*y') :+ (x*y'+y*x')
+    one = one :+ zero
+instance Group a => Group  (Complex a) where
+    (x:+y) - (x':+y')   =  (x-x') :+ (y-y')
+    negate (x:+y)       =  negate x :+ negate y
+instance AbelianAdditive a => AbelianAdditive (Complex a)
+instance Ring a => Module (Complex a) (Complex a) where
+  (*^) = (*)
+instance Ring a => Module a (Complex a) where
+  s *^ (x :+ y) =  (s *^ x :+ s *^ y)
+instance Ring a => Ring (Complex a) where
+    fromInteger n  =  fromInteger n :+ zero
+
+instance  Field a => Division (Complex a)  where
+    {-# SPECIALISE instance Division (Complex Double) #-}
+    (x:+y) / (x':+y')   =  (x*x'+y*y') / d :+ (y*x'-x*y') / d
+      where d   = x'*x' + y'*y'
+
+instance Field a => Field (Complex a) where
+    fromRational a =  fromRational a :+ zero
+
+{-data Expr a where
+  Embed :: a -> Expr a
+  Add :: Expr a -> Expr a -> Expr a
+  Mul :: Expr a -> Expr a -> Expr a
+  Zero :: Expr a
+  One :: Expr a
+  deriving (Prelude.Show)
+
+
+instance Additive (Expr a) where
+  zero = Zero
+  Zero + x = x
+  x + Zero = x
+  x + y = Add x y
+
+instance Multiplicative (Expr a) where
+  one = One
+  One * x = x
+  x * One = x
+  x * y = Mul x y
+-}
+
+-- Syntax
+
 ifThenElse :: Bool -> t -> t -> t
 ifThenElse True a _ = a
 ifThenElse False _ a = a
 
 
-data InitialAdditive = InitialAdditive :+ InitialAdditive | Zero
-  deriving (Prelude.Show)
 
-instance Additive InitialAdditive where
-  zero = Zero
-  (+) = (:+)
+-- >>> times 5 (Embed "x")
+-- Add (Add (Embed "x") (Add (Embed "x") (Embed "x"))) (Add (Embed "x") (Embed "x"))
 
-instance Module Rational Double where
-  r *^ d = fromRational r * d
+
+-- >>> (Embed "x")
+-- Zero
+
+
+{-
+Note: the following is not quite what we intuitively want, because
+
+class Field a => AlgebraicallyClosed  a where
+  sqrt :: a -> (a,a)
+
+AlgebraicallyClosed numbers have two square roots.
+
+-}
+
diff --git a/Algebra/Linear.hs b/Algebra/Linear.hs
--- a/Algebra/Linear.hs
+++ b/Algebra/Linear.hs
@@ -21,16 +21,36 @@
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE RebindableSyntax #-}
 
 module Algebra.Linear where
 
-import Algebra.Classes
-import Prelude (cos,sin,Floating(..),Functor(..),Show(..),Eq(..),Int,fst,flip,($),Ord)
+import Algebra.Classes hiding ((*<))
+import Prelude (cos,sin,Floating(..),Functor(..),Show(..),Eq(..),Int,fst,($),Ord,Double)
 import Control.Applicative
-import Data.Foldable
+import Data.Foldable hiding (sum,product)
 import Data.Traversable
 import Control.Monad.State
 import Algebra.Category
+
+infixr 7 *<
+
+type VectorSpace scalar a = (Field scalar, Module scalar a, Group a)
+-- Because of the existence of bases, vector spaces can always be made representable (Traversable, Applicative) functors.
+-- So we'd be better off using the following definition:
+
+-- | Representation of vector as traversable functor
+type VectorR v = (Applicative v,Traversable v)
+-- ... but this is missing the link with *^ for module.  We should be
+-- able to add forall s. PreRing s => Module s (v s), but GHC does not
+-- like it. (In fact, QuantifiedConstraints is very buggy in ghc 8.6)
+
+class VectorR v => InnerProdSpace v where
+  inner :: Field s => v s -> v s -> s
+
+--------------------------------------------------------------
+-- Construction of finite vectors
+
 data VZero a = VZero deriving (Functor,Foldable,Traversable,Show,Eq,Ord)
 instance Applicative VZero where
   pure _ = VZero
@@ -53,7 +73,10 @@
 pattern V3' :: forall a. a -> a -> a -> V3' a
 pattern V3' x y z = VNext (V2' x y) z
 
--- | Make a Euclidean vector out of a traversable functor
+--------------------------------------------
+-- Euclidean spaces with a (inner product)
+
+-- | Make a Euclidean vector out of a traversable functor. (The p)
 newtype Euclid f a = Euclid {fromEuclid :: f a} deriving (Functor,Foldable,Traversable,Show,Eq,Ord,Applicative)
 
 type V3 = Euclid V3'
@@ -78,49 +101,38 @@
 pureMat :: (Applicative v, Applicative w) => s -> Mat s v w
 pureMat x = Mat (pure (pure x))
 
-(>*<) :: (Applicative v, Applicative w) =>
-               Mat (a -> s) v w -> Mat a v w -> Mat s v w
-Mat f >*< Mat x = Mat (((<*>) <$> f) <*> x)
-
-(>$<) :: (Applicative v, Applicative w) =>
-               (a -> s) -> Mat a v w -> Mat s v w
-f >$< x = pureMat f >*< x
-
 instance (Applicative f,Applicative g,Additive a) => Additive (Mat a f g) where
   zero = pureMat zero
-  x + y =  (+) >$< x >*< y
+  x + y = matFlat ((+) <$> flatMat x <*> flatMat y)
 instance (Applicative f,Applicative g,AbelianAdditive a) => AbelianAdditive (Mat a f g) where
 instance (Applicative f,Applicative g,Group a) => Group (Mat a f g) where
-  negate x = negate >$< x
-  x - y = (-) >$< x >*< y
+  negate x = matFlat (negate <$> flatMat x)
+  x - y = matFlat ((-) <$> flatMat x <*> flatMat y)
 
 instance (Applicative f, Applicative g,Module s a) => Module s (Mat a f g) where
   s *^ Mat t = Mat (((s*^) <$>) <$> t)
 
-class VectorSpace (Scalar v) v => InnerProdSpace v where
-  type Scalar v
-  dotProd :: v -> v -> Scalar v
 
+-- | Hadamard product
 (⊙) :: Applicative v => Multiplicative s => v s -> v s -> v s
 x ⊙ y = (*) <$> x <*> y
 
-instance (Ring a, Field a, Applicative f, Foldable f) => InnerProdSpace (Euclid f a) where
-  type Scalar (Euclid f a) = a
-  dotProd x y = add (x ⊙ y)
+instance (VectorR f) => InnerProdSpace (Euclid f) where
+  inner x y = sum (x ⊙ y) -- fixme
 
-(·) :: InnerProdSpace v => v -> v -> Scalar v
-(·) = dotProd
+(·) :: Field s => InnerProdSpace v => v s -> v s -> s
+(·) = inner
 
-sqNorm :: InnerProdSpace v => v -> Scalar v
-sqNorm x = dotProd x x
+sqNorm :: Field s => InnerProdSpace v => v s -> s
+sqNorm x = inner x x
 
-norm :: InnerProdSpace v => Floating (Scalar v) => v  -> Scalar v
+norm :: Field s => InnerProdSpace v => Floating s => v s  -> s
 norm = sqrt . sqNorm
 
-normalize :: Floating (Scalar v) => InnerProdSpace v => v -> v
+normalize :: (VectorSpace s (v s)) => Floating s => InnerProdSpace v => v s -> v s
 normalize v = recip (norm v) *^ v
 
--- | Cross product https://en.wikipedia.org/wiki/Cross_product
+-- | Cross product in 3 dimensions https://en.wikipedia.org/wiki/Cross_product
 (×) :: Ring a => V3 a -> V3 a -> V3 a
 (V3 a1 a2 a3) × (V3 b1 b2 b3) = V3 (a2*b3 - a3*b2)  (negate (a1*b3 - a3*b1)) (a1*b2 - a2*b1)
 
@@ -129,52 +141,77 @@
   where increment = do x <- get; put (x+1); return x
 
 type SqMat v s = Mat s v v
-newtype Mat s w v = Mat {fromMat :: v (w s)} deriving Show
 
+-- | Matrix type. (w s) is a column. (v s) is a row.
+newtype Mat s w v = Mat {fromMat :: w (v s)} deriving Show
+
+-- | View of the matrix as a composition of functors.
+newtype Flat w v s = Flat {fromFlat :: w (v s)} deriving (Show,Functor,Foldable)
+flatMat :: Mat s w v -> Flat w v s
+flatMat (Mat x) = (Flat x)
+matFlat :: Flat w v s -> Mat s w v
+matFlat (Flat x) = (Mat x)
+
+instance (Applicative w, Applicative v) => Applicative (Flat w v) where
+  pure x = Flat (pure (pure x))
+  Flat f <*> Flat a = Flat (((<*>) <$> f) <*> a)
+
+
 instance Ring s => Category (Mat s) where
-  type Con v = (Applicative v, Traversable v)
+  type Con v = VectorR v
   (.) = matMul
   id = identity
 
-type Mat3x3 s = SqMat V3' s
-type Mat2x2 s = SqMat V2' s
+type Mat3x3 s = SqMat V3 s
+type Mat2x2 s = SqMat V2 s
 
-pattern Mat2x2 :: forall s. s -> s -> s -> s -> Mat s V2' V2'
-pattern Mat2x2 a b c d = Mat (V2' (V2' a b) (V2' c d))
+pattern Mat2x2 :: forall s. s -> s -> s -> s -> Mat s V2 V2
+pattern Mat2x2 a b c d = Mat (V2 (V2 a c)
+                                 (V2 b d))
 
-pattern Mat3x3 :: forall s. s -> s -> s -> s -> s -> s -> s -> s -> s -> Mat s V3' V3'
-pattern Mat3x3 a b c d e f g h i = Mat (V3' (V3' a b c) (V3' d e f) (V3' g h i))
+pattern Mat3x3 :: forall s. s -> s -> s -> s -> s -> s -> s -> s -> s -> Mat s V3 V3
+pattern Mat3x3 a b c d e f g h i = Mat (V3 (V3 a d g)
+                                           (V3 b e h)
+                                           (V3 c f i))
 
-matVecMul :: (Foldable f1, Ring b, Applicative f1, Functor f2) => Mat b f1 f2 -> Euclid f1 b -> Euclid f2 b
-matVecMul (Mat m) v = Euclid (euclideanDotProd v <$> (Euclid <$> m))
-   where euclideanDotProd x y = add (Euclid x ⊙ Euclid y)
+-- | Vector scaling. If Module a (f a), then (*^) must be the same as (*<).
+(*<) :: (Functor f, Multiplicative b) => b -> f b -> f b
+s *< v = (s*) <$> v
 
-rotation2d :: Floating a => a -> Mat2x2 a
-rotation2d θ = Mat $ V2' (V2' (cos θ) (-sin θ))
-                         (V2' (sin θ)  (cos θ))
 
+(<+>) :: (Applicative f, Additive b) => f b -> f b -> f b
+u <+> v = (+) <$> u <*> v
+
+
+matVecMul :: forall s v w. (Ring s, Foldable v,Applicative v,Applicative w) => Mat s v w -> v s -> w s
+matVecMul (Mat m) x = foldr (<+>) (pure zero) ((*<) <$> x <*> m) -- If GHC gets fixed: use VectorR constraint instead of Applicative, and add instead of foldr.
+
+rotation2d :: (Group a,Floating a) => a -> Mat2x2 a
+rotation2d θ = transpose $ Mat $ V2 (V2 (cos θ) (-sin θ))
+                                    (V2 (sin θ)  (cos θ))
+
 -- >>> rotation2d (pi/2)
 -- Mat {fromMat = V2' (V2' 6.123233995736766e-17 (-1.0)) (V2' 1.0 6.123233995736766e-17)}
 
 crossProductMatrix :: Group a => V3 a -> Mat3x3 a
-crossProductMatrix (V3 a1 a2 a3) = Mat (V3'  (V3' zero (negate a3) a2)
-                                             (V3' a3 zero (negate a1))
-                                             (V3' (negate a2) a1 zero))
+crossProductMatrix (V3 a1 a2 a3) = Mat3x3 zero  (-a3) a2
+                                          a3    zero  (-a1)
+                                          (-a2) a1    zero
 
 -- | Tensor product
 (⊗) :: (Applicative v, Applicative w, Multiplicative s)
-    => Euclid w s -> Euclid v s -> Mat s w v
-Euclid v1 ⊗ Euclid v2 = flip (tensorWith (*)) v1 v2
+    => w s -> v s -> Mat s w v
+v1 ⊗ v2 = tensorWith (*) v2 v1
 
 tensorWith :: (Applicative v, Applicative w)
            => (s -> t -> u) -> w s -> v t -> Mat u v w
-tensorWith f v1 v2 = flip f >$< Mat (pure v2) >*< Mat (pure <$> v1)
+tensorWith f v1 v2 = matFlat (f <$> Flat (pure v1) <*> Flat (pure <$> v2))
 
 identity :: Traversable v => Ring s => Applicative v => SqMat v s
 identity = tensorWith (\x y -> if x == y then one else zero) index index
 
-diagonal :: Traversable v => Ring s => Applicative v => Euclid v s -> SqMat v s
-diagonal (Euclid v) = tensorWith (\x (y,a) -> if x == y then a else zero) index ((,) <$> index <*> v)
+diagonal :: Traversable v => Ring s => Applicative v => v s -> SqMat v s
+diagonal v = tensorWith (\x (y,a) -> if x == y then a else zero) index ((,) <$> index <*> v)
 
 -- | 3d rotation around given axis
 rotation3d :: Ring a => Floating a => a -> V3 a -> Mat3x3 a
@@ -182,27 +219,30 @@
                  sin θ *^ crossProductMatrix u +
                  (1 - cos θ) *^ (u ⊗ u)
 
+-- | 3d rotation mapping the direction of 'from' to that of 'to'
 rotationFromTo :: (Floating a, Module a a,Field a)
                => V3 a -> V3 a -> Mat3x3 a
 rotationFromTo from to = c *^ identity + s *^ crossProductMatrix v + (1-c) *^ (v ⊗ v)
   where y = to
         x = from
-        v = x × y
-        c = dotProd x y
-        s = norm v
+        v = x × y -- axis of rotation
+        c = inner x y -- cos of angle
+        s = norm v -- sin of angle
 
-transpose :: Applicative f => Traversable g => Mat a f g -> Mat a g f
-transpose = Mat . sequenceA . fromMat
+-- >>> let u = (V3 (1::Double) 0 0); v = (V3 0 1 1); in (rotationFromTo u v) `matVecMul` u
+-- Euclid {fromEuclid = VNext (VNext (VNext VZero 0.0) 1.4142135623730951) 1.4142135623730951}
 
-matMul' :: (Traversable u, Ring s, Applicative w, Applicative v, Applicative u) => Mat s v u -> Mat s u w -> Mat s v w
-matMul' (transpose -> Mat y) (Mat x)  = tensorWith (\a b -> add (a ⊙ b)) x y
+transpose :: Applicative g => Traversable f => Mat a f g -> Mat a g f
+transpose = Mat . sequenceA . fromMat
 
 matMul :: (Traversable u, Ring s, Applicative w, Applicative v, Applicative u) => Mat s u w -> Mat s v u -> Mat s v w
-matMul = flip matMul'
+matMul a (Mat b) = Mat (matVecMul a <$> b)
 
 
--- >>> let t1 = rotation2d (1::Double) in matMul t1 (transpose t1)
--- Mat {fromMat = V2' (V2' 1.0 0.0) (V2' 0.0 1.0)}
+-- >>> let t1 = rotation2d (1::Double) in matMul (transpose t1) t1
+-- Mat {fromMat = VNext (VNext VZero (VNext (VNext VZero 1.0) 0.0)) (VNext (VNext VZero 0.0) 1.0)}
+
+
 
 -- The group of Orthogonal matrices, using "Multiplicative" for respecting conventions a bit better
 newtype OrthoMat v s = OrthoMat (SqMat v s)
diff --git a/gasp.cabal b/gasp.cabal
--- a/gasp.cabal
+++ b/gasp.cabal
@@ -1,5 +1,5 @@
 name:           gasp
-version:        1.2.0.0
+version:        1.3.0.0
 category:       Algebra
 synopsis:       A framework of algebraic classes
 description:
@@ -22,6 +22,7 @@
   build-depends: containers
   build-depends: binary
   build-depends: mtl
+  build-depends: QuickCheck
 
   exposed-modules:
        Algebra.Classes