diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal
--- a/HaskellForMaths.cabal
+++ b/HaskellForMaths.cabal
@@ -1,5 +1,5 @@
    Name:                HaskellForMaths
-   Version:             0.3.1
+   Version:             0.3.2
    Category:            Math
    Description:         A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended for educational purposes, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra
@@ -22,6 +22,8 @@
         Math/Test/TRootSystem.hs,
         Math/Test/TSubquotients.hs,
         Math/Test/TestAll.hs
+        Math/Test/TAlgebras/TVectorSpace.hs
+        Math/Test/TAlgebras/TTensorProduct.hs
         Math/Test/TAlgebras/TStructures.hs
         Math/Test/TAlgebras/TQuaternions.hs
         Math/Test/TAlgebras/TGroupAlgebra.hs
diff --git a/Math/Algebras/AffinePlane.hs b/Math/Algebras/AffinePlane.hs
--- a/Math/Algebras/AffinePlane.hs
+++ b/Math/Algebras/AffinePlane.hs
@@ -41,7 +41,7 @@
     unit 0 = zero -- V []
     unit x = V [(munit,x)] where munit = SL2 (Glex 0 [])
     mult x = x''' where
-        x' = mult $ fmap (\(T (SL2 a) (SL2 b)) -> T a b) x -- perform the multiplication in GlexPoly
+        x' = mult $ fmap ( \(SL2 a, SL2 b) -> (a,b) ) x -- perform the multiplication in GlexPoly
         x'' = x' %% [a*d-b*c-1] -- :: GlexPoly Q ABCD] -- quotient by ad-bc=1 in GlexPoly Q ABCD
         x''' = fmap SL2 x'' -- ie wrap the monomials up as SL2 again
         -- mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs
diff --git a/Math/Algebras/Commutative.hs b/Math/Algebras/Commutative.hs
--- a/Math/Algebras/Commutative.hs
+++ b/Math/Algebras/Commutative.hs
@@ -50,7 +50,7 @@
 instance (Num k, Ord v) => Algebra k (GlexMonomial v) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)] where munit = Glex 0 []
-    mult (V ts) = nf $ fmap (\(T a b) -> a `mmult` b) (V ts)
+    mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts)
         where mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs
 
 {-
diff --git a/Math/Algebras/GroupAlgebra.hs b/Math/Algebras/GroupAlgebra.hs
--- a/Math/Algebras/GroupAlgebra.hs
+++ b/Math/Algebras/GroupAlgebra.hs
@@ -23,14 +23,14 @@
 instance Num k => Algebra k (Permutation Int) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
-    mult = nf . fmap (\(T a b) -> a `mmult` b)
+    mult = nf . fmap (\(a,b) -> a `mmult` b)
 
 -- Set Coalgebra instance
 -- instance SetCoalgebra (Permutation Int) where {}
 
 instance Num k => Coalgebra k (Permutation Int) where
     counit (V ts) = sum [x | (m,x) <- ts] -- trace
-    comult = fmap (\m -> T m m) -- diagonal
+    comult = fmap (\m -> (m,m)) -- diagonal
 
 instance Num k => Bialgebra k (Permutation Int) where {}
 -- should check that the algebra and coalgebra structures are compatible
@@ -44,7 +44,7 @@
 
 
 instance Num k => Module k (Permutation Int) Int where
-    action = nf . fmap (\(T g x) -> x .^ g)
+    action = nf . fmap (\(g,x) -> x .^ g)
 
 -- use *. instead
 -- r *> m = action (r `te` m)
diff --git a/Math/Algebras/LaurentPoly.hs b/Math/Algebras/LaurentPoly.hs
--- a/Math/Algebras/LaurentPoly.hs
+++ b/Math/Algebras/LaurentPoly.hs
@@ -34,7 +34,7 @@
 instance Num k => Algebra k LaurentMonomial where
     unit 0 = zero -- V []
     unit x = V [(munit,x)] 
-    mult (V ts) = nf $ fmap (\(T a b) -> a `mmult` b) (V ts)
+    mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts)
     -- mult (V ts) = nf $ V [(a `mmult` b, x) | (T a b, x) <- ts]
 
 {-
diff --git a/Math/Algebras/Matrix.hs b/Math/Algebras/Matrix.hs
--- a/Math/Algebras/Matrix.hs
+++ b/Math/Algebras/Matrix.hs
@@ -24,7 +24,7 @@
     unit x = x `smultL` V [(E2 i i, 1) | i <- [1..2] ]
     -- mult ab = nf $ ab >>= mult' where
     mult = linear mult' where
-        mult' (T (E2 i j) (E2 k l)) = delta j k `smultL` return (E2 i l)
+        mult' (E2 i j, E2 k l) = delta j k `smultL` return (E2 i l)
 
 -- In other words
 -- unit x = x (1 0)
@@ -35,7 +35,7 @@
 instance Num k => Module k Mat2 EBasis where
     -- action ax = nf $ ax >>= action' where
     action = linear action' where
-        action' (T (E2 i j) (E k)) = delta j k `smultL` return (E i)
+        action' (E2 i j, E k) = delta j k `smultL` return (E i)
 
 -- In other words
 -- action (a b) `te` (x) = (ax+by)
@@ -57,7 +57,7 @@
 instance Num k => Coalgebra k Mat2' where
     counit (V ts) = sum [xij * delta i j | (E2' i j, xij) <- ts]
     -- comult (V ts) = V $ concatMap (\(E2' i j,xij) -> [(T (E2' i k) (E2' k j), xij) | k <- [1..2]]) ts
-    comult = linear (\(E2' i j) -> foldl (<+>) zero [return (T (E2' i k) (E2' k j)) | k <- [1..2]])
+    comult = linear (\(E2' i j) -> foldl (<+>) zero [return (E2' i k, E2' k j) | k <- [1..2]])
 -- In other words
 -- counit (a b) = (1 0)
 --        (c d)   (0 1)
@@ -79,13 +79,13 @@
 instance Num k => Algebra k M3 where
     unit 0 = zero -- V []
     unit x = V [(E3 i i, x) | i <- [1..3] ]
-    mult (V ts) = nf $ V $ map (\(T (E3 i j) (E3 k l), x) -> (E3 i l, delta j k * x)) ts
+    mult (V ts) = nf $ V $ map (\((E3 i j, E3 k l), x) -> (E3 i l, delta j k * x)) ts
 
 {-
 -- Kassel p42
 -- In this coalgebra instance, the E3 i j are to be interpreted as the dual basis, not the original basis
 instance Num k => Coalgebra k M3 where
     counit (V ts) = sum [xij * delta i j | (E3 i j, xij) <- ts]
-    comult (V ts) = V $ concatMap (\(E3 i j,xij) -> [(T (E3 i k) (E3 k j), xij) | k <- [1..3]]) ts
+    comult (V ts) = V $ concatMap (\(E3 i j,xij) -> [((E3 i k, E3 k j), xij) | k <- [1..3]]) ts
 -- (is this order preserving?)
 -}
diff --git a/Math/Algebras/NonCommutative.hs b/Math/Algebras/NonCommutative.hs
--- a/Math/Algebras/NonCommutative.hs
+++ b/Math/Algebras/NonCommutative.hs
@@ -32,14 +32,14 @@
 instance (Num k, Ord v) => Algebra k (NonComMonomial v) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
-    mult = nf . fmap (\(T a b) -> a `mmult` b)
+    mult = nf . fmap (\(a,b) -> a `mmult` b)
 
 {-
 -- This is the monoid algebra for non-commutative monomials (which is the free monoid)
 instance (Num k, Ord v) => Algebra k (NonComMonomial v) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)] where munit = NCM 0 []
-    mult (V ts) = nf $ fmap (\(T a b) -> a `mmult` b) (V ts)
+    mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts)
         where mmult (NCM lu us) (NCM lv vs) = NCM (lu+lv) (us++vs)
     -- mult (V ts) = nf $ V [(a `mmult` b, x) | (T a b, x) <- ts]
 -}
@@ -49,7 +49,7 @@
 -- Also, not used anywhere. Hence commented out
 instance Num k => Coalgebra k (NonComMonomial v) where
     counit (V ts) = sum [x | (m,x) <- ts] -- trace
-    comult = fmap (\m -> T m m)
+    comult = fmap (\m -> (m,m))
 -}
 
 
diff --git a/Math/Algebras/Quaternions.hs b/Math/Algebras/Quaternions.hs
--- a/Math/Algebras/Quaternions.hs
+++ b/Math/Algebras/Quaternions.hs
@@ -28,17 +28,17 @@
     unit x = V [(One,x)]
     -- mult x = nf (x >>= m)
     mult = linear m
-         where m (T One b) = return b
-               m (T b One) = return b
-               m (T I I) = unit (-1)
-               m (T J J) = unit (-1)
-               m (T K K) = unit (-1)
-               m (T I J) = return K
-               m (T J I) = -1 *> return K
-               m (T J K) = return I
-               m (T K J) = -1 *> return I
-               m (T K I) = return J
-               m (T I K) = -1 *> return J
+         where m (One,b) = return b
+               m (b,One) = return b
+               m (I,I) = unit (-1)
+               m (J,J) = unit (-1)
+               m (K,K) = unit (-1)
+               m (I,J) = return K
+               m (J,I) = -1 *> return K
+               m (J,K) = return I
+               m (K,J) = -1 *> return I
+               m (K,I) = return J
+               m (I,K) = -1 *> return J
 
 i,j,k :: Num k => Quaternion k
 i = return I
@@ -55,4 +55,4 @@
 instance Num k => Coalgebra k HBasis where
     counit (V ts) = sum [x | (One,x) <- ts]
     comult = linear cm
-        where cm m = if m == One then return (T m m) else return (T m One) <+> return (T One m)
+        where cm m = if m == One then return (m,m) else return (m,One) <+> return (One,m)
diff --git a/Math/Algebras/Structures.hs b/Math/Algebras/Structures.hs
--- a/Math/Algebras/Structures.hs
+++ b/Math/Algebras/Structures.hs
@@ -4,6 +4,8 @@
 {-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}
 {-# LANGUAGE IncoherentInstances #-}
 
+-- |A module defining various algebraic structures that can be defined on vector spaces
+-- - specifically algebra, coalgebra, bialgebra, Hopf algebra, module, comodule
 module Math.Algebras.Structures where
 
 import Math.Algebras.VectorSpace
@@ -20,12 +22,14 @@
 
 -- ALGEBRAS, COALGEBRAS, BIALGEBRAS, HOPF ALGEBRAS
 
--- |"Vect k b is a k-algebra"
+-- |Caution: If we declare an instance Algebra k b, then we are saying that the vector space Vect k b is a k-algebra.
+-- In other words, we are saying that b is the basis for a k-algebra. So a more accurate name for this class
+-- would have been AlgebraBasis.
 class Algebra k b where
     unit :: k -> Vect k b
     mult :: Vect k (Tensor b b) -> Vect k b
 
--- |"Vect k b is a k-coalgebra"
+-- |An instance declaration for Coalgebra k b is saying that the vector space Vect k b is a k-algebra.
 class Coalgebra k b where
     counit :: Vect k b -> k
     comult :: Vect k b -> Vect k (Tensor b b)
@@ -63,13 +67,15 @@
 instance Num k => Algebra k () where
     unit 0 = zero -- V []
     unit x = V [( (),x)]
-    mult (V [(T () (),x)]) = V [( (),x)]
+    mult (V [( ((),()), x)]) = V [( (),x)]
 
 instance Num k => Coalgebra k () where
     counit (V []) = 0
     counit (V [( (),x)]) = x
-    comult (V [( (),x)]) = V [(T () (),x)]
+    comult (V [( (),x)]) = V [( ((),()), x)]
 
+-- |Trivial k is the field k considered as a k-vector space. In maths, we would not normally make a distinction here,
+-- but in the code, we need this if we want to be able to put k as one side of a tensor product.
 type Trivial k = Vect k ()
 
 unit' :: (Num k, Algebra k b) => Trivial k -> Vect k b
@@ -87,7 +93,7 @@
     unit x = x `smultL` (unit 1 `te` unit 1)
     -- mult x = nf $ x >>= m where
     mult = linear m where
-        m (T (T a b) (T a' b')) = (mult $ return $ T a a') `te` (mult $ return $ T b b')
+        m ((a,b),(a',b')) = (mult $ return (a,a')) `te` (mult $ return (b,b'))
 
 -- Kassel p42
 instance (Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (Tensor a b) where
@@ -101,7 +107,7 @@
 
 instance Num k => Coalgebra k (SetCoalgebra b) where
     counit (V ts) = sum [x | (m,x) <- ts] -- trace
-    comult = fmap (\m -> T m m)           -- diagonal
+    comult = fmap ( \m -> (m,m) )           -- diagonal
 
 
 newtype MonoidCoalgebra m = MC m deriving (Eq,Ord,Show)
@@ -109,7 +115,7 @@
 instance (Num k, Ord m, Mon m) => Coalgebra k (MonoidCoalgebra m) where
     counit (V ts) = sum [if m == MC munit then x else 0 | (m,x) <- ts]
     comult = linear cm
-        where cm m = if m == MC munit then return (T m m) else return (T m (MC munit)) <+> return (T (MC munit) m)
+        where cm m = if m == MC munit then return (m,m) else return (m, MC munit) <+> return (MC munit, m)
 -- Brzezinski and Wisbauer, Corings and Comodules, p5
 
 -- Both of the above can be used to define coalgebra structure on polynomial algebras
@@ -142,13 +148,13 @@
          => Module k (Tensor a a) (Tensor u v) where
     -- action x = nf $ x >>= action'
     action = linear action'
-        where action' (T (T a a') (T u v)) = (action $ return $ T a u) `te` (action $ return $ T a' v)
+        where action' ((a,a'), (u,v)) = (action $ return (a,u)) `te` (action $ return (a',v))
 
 instance (Num k, Ord a, Ord u, Ord v, Bialgebra k a, Module k a u, Module k a v)
          => Module k a (Tensor u v) where
     -- action x = nf $ x >>= action'
     action = linear action'
-        where action' (T a (T u v)) = action $ (comult $ return a) `te` (return $ T u v)
+        where action' (a,(u,v)) = action $ (comult $ return a) `te` (return (u,v))
 -- !! Overlapping instances
 -- If a == Tensor b b, then we have overlapping instance with the previous definition
 -- On the other hand, if a == Tensor u v, then we have overlapping instance with the earlier instance
@@ -157,4 +163,4 @@
 instance (Num k, Ord a, Ord m, Ord n, Bialgebra k a, Comodule k a m, Comodule k a n)
          => Comodule k a (Tensor m n) where
     coaction = (mult `tf` id) . twistm . (coaction `tf` coaction)
-        where twistm x = nf $ fmap (\(T (T h m) (T h' n)) -> T (T h h') (T m n)) x
+        where twistm x = nf $ fmap ( \((h,m), (h',n)) -> ((h,h'), (m,n)) ) x
diff --git a/Math/Algebras/TensorAlgebra.hs b/Math/Algebras/TensorAlgebra.hs
--- a/Math/Algebras/TensorAlgebra.hs
+++ b/Math/Algebras/TensorAlgebra.hs
@@ -23,7 +23,7 @@
 instance (Num k, Ord a) => Algebra k (TensorAlgebra a) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
-    mult = nf . fmap (\(T a b) -> a `mmult` b)
+    mult = nf . fmap (\(a,b) -> a `mmult` b)
 
 
 data SymmetricAlgebra a = Sym Int [a] deriving (Eq,Ord,Show)
@@ -35,7 +35,7 @@
 instance (Num k, Ord a) => Algebra k (SymmetricAlgebra a) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
-    mult = nf . fmap (\(T a b) -> a `mmult` b)
+    mult = nf . fmap (\(a,b) -> a `mmult` b)
 
 
 data ExteriorAlgebra a = Ext Int [a] deriving (Eq,Ord,Show)
@@ -43,7 +43,7 @@
 instance (Num k, Ord a) => Algebra k (ExteriorAlgebra a) where
     unit 0 = zero -- V []
     unit x = V [(Ext 0 [],x)]
-    mult xy = nf $ xy >>= (\(T (Ext i xs) (Ext j ys)) -> signedMerge 1 (0,[]) (i,xs) (j,ys))
+    mult xy = nf $ xy >>= (\(Ext i xs, Ext j ys) -> signedMerge 1 (0,[]) (i,xs) (j,ys))
         where signedMerge s (k,zs) (i,x:xs) (j,y:ys) =
                   case compare x y of
                   EQ -> zero
diff --git a/Math/Algebras/TensorProduct.hs b/Math/Algebras/TensorProduct.hs
--- a/Math/Algebras/TensorProduct.hs
+++ b/Math/Algebras/TensorProduct.hs
@@ -2,26 +2,83 @@
 
 {-# LANGUAGE NoMonomorphismRestriction #-}
 
--- |A module defining tensor products of vector spaces
+-- |A module defining direct sum and tensor product of vector spaces
 module Math.Algebras.TensorProduct where
 
 import Math.Algebras.VectorSpace
 
+-- DIRECT SUM
 
-data Tensor a b = T a b deriving (Eq, Ord, Show)
--- or T !a !b, forcing strictness, but not proven to be better
+-- |A type for constructing a basis for the direct sum of vector spaces.
+-- The direct sum of Vect k a and Vect k b is Vect k (DSum a b)
+type DSum a b = Either a b
 
+-- |Injection of left summand into direct sum
+i1 :: Vect k a -> Vect k (DSum a b)
+i1 = fmap Left
 
--- |Tensor product of two elements
+-- |Injection of right summand into direct sum
+i2 :: Vect k b -> Vect k (DSum a b)
+i2 = fmap Right
+
+-- |The coproduct of two linear functions (with the same target).
+-- Satisfies the universal property that f == coprodf f g . i1 and g == coprodf f g . i2
+coprodf :: (Num k, Ord t) =>
+    (Vect k a -> Vect k t) -> (Vect k b -> Vect k t) -> Vect k (DSum a b) -> Vect k t
+coprodf f g = linear fg' where
+    fg' (Left a) = f (return a)
+    fg' (Right b) = g (return b)
+
+
+-- |Projection onto left summand from direct sum
+p1 :: (Num k, Ord a) => Vect k (DSum a b) -> Vect k a
+p1 = linear p1' where
+    p1' (Left a) = return a
+    p1' (Right b) = zero
+
+-- |Projection onto right summand from direct sum
+p2 :: (Num k, Ord b) => Vect k (DSum a b) -> Vect k b
+p2 = linear p2' where
+    p2' (Left a) = zero
+    p2' (Right b) = return b
+
+-- |The product of two linear functions (with the same source).
+-- Satisfies the universal property that f == p1 . prodf f g and g == p2 . prodf f g
+prodf :: (Num k, Ord a, Ord b) =>
+    (Vect k s -> Vect k a) -> (Vect k s -> Vect k b) -> Vect k s -> Vect k (DSum a b)
+prodf f g = linear fg' where
+    fg' b = fmap Left (f $ return b) <+> fmap Right (g $ return b)
+
+
+-- |The direct sum of two vector space elements
+dsume :: (Num k, Ord a, Ord b) => Vect k a -> Vect k b -> Vect k (DSum a b)
+-- dsume x y = fmap Left x <+> fmap Right y
+dsume x y = i1 x <+> i2 y
+
+-- |The direct sum of two linear functions.
+-- Satisfies the universal property that f == p1 . dsumf f g . i1 and g == p2 . dsumf f g . i2
+dsumf :: (Num k, Ord a, Ord b, Ord a', Ord b') => 
+    (Vect k a -> Vect k a') -> (Vect k b -> Vect k b') -> Vect k (DSum a b) -> Vect k (DSum a' b')
+dsumf f g ab = (i1 . f . p1) ab <+> (i2 . g . p2) ab
+
+
+-- TENSOR PRODUCT
+
+-- |A type for constructing a basis for the tensor product of vector spaces.
+-- The tensor product of Vect k a and Vect k b is Vect k (Tensor a b)
+type Tensor a b = (a,b)
+
+-- |The tensor product of two vector space elements
 te :: Num k => Vect k a -> Vect k b -> Vect k (Tensor a b)
-te (V us) (V vs) = V [(T ei ej, xi*xj) | (ei,xi) <- us, (ej,xj) <- vs]
+te (V us) (V vs) = V [((a,b), x*y) | (a,x) <- us, (b,y) <- vs]
+-- te (V us) (V vs) = V [((ei,ej), xi*xj) | (ei,xi) <- us, (ej,xj) <- vs]
 -- preserves order - that is, if the inputs are correctly ordered, so is the output
 
 -- Implicit assumption - f and g are linear
--- |Tensor product of two (linear) functions
+-- |The tensor product of two linear functions
 tf :: (Num k, Ord a', Ord b') => (Vect k a -> Vect k a') -> (Vect k b -> Vect k b')
    -> Vect k (Tensor a b) -> Vect k (Tensor a' b')
-tf f g (V ts) = sum [te (f $ V [(a, 1)]) (g $ V [(b, x)]) | (T a b, x) <- ts]
+tf f g (V ts) = sum [x *> te (f $ return a) (g $ return b) | ((a,b), x) <- ts]
     where sum = foldl add zero -- (V [])
 
 
@@ -29,19 +86,38 @@
 
 -- in fact, this definition works for any Functor f, not just (Vect k)
 assocL :: Vect k (Tensor u (Tensor v w)) -> Vect k (Tensor (Tensor u v) w)
-assocL = fmap (\(T a (T b c)) -> T (T a b) c)
+assocL = fmap ( \(a,(b,c)) -> ((a,b),c) )
 
 assocR :: Vect k (Tensor (Tensor u v) w) -> Vect k (Tensor u (Tensor v w))
-assocR = fmap (\(T (T a b) c) -> T a (T b c))
+assocR = fmap ( \((a,b),c) -> (a,(b,c)) )
 
-inUnitL = fmap (\a -> T () a)
+unitInL = fmap ( \a -> ((),a) )
 
-inUnitR = fmap (\a -> T a ())
+unitOutL = fmap ( \((),a) -> a )
 
-outUnitL = fmap (\(T () a) -> a)
+unitInR = fmap ( \a -> (a,()) )
 
-outUnitR = fmap (\(T a ()) -> a)
+unitOutR = fmap ( \(a,()) -> a )
 
-twist v = nf $ fmap (\(T a b) -> T b a) v
+twist v = nf $ fmap ( \(a,b) -> (b,a) ) v
 -- note the nf call, as f is not order-preserving
 
+
+distrL :: (Num k, Ord a, Ord b, Ord c)
+    => Vect k (Tensor a (DSum b c)) -> Vect k (DSum (Tensor a b) (Tensor a c))
+distrL v = nf $ fmap (\(a,bc) -> case bc of Left b -> Left (a,b); Right c -> Right (a,c)) v
+
+undistrL :: (Num k, Ord a, Ord b, Ord c)
+    => Vect k (DSum (Tensor a b) (Tensor a c)) -> Vect k (Tensor a (DSum b c))
+undistrL v = nf $ fmap ( \abc -> case abc of Left (a,b) -> (a,Left b); Right (a,c) -> (a,Right c) ) v
+
+distrR :: Vect k (Tensor (DSum a b) c) -> Vect k (DSum (Tensor a c) (Tensor b c))
+distrR v = fmap ( \(ab,c) -> case ab of Left a -> Left (a,c); Right b -> Right (b,c) ) v
+-- order-preserving, so no nf call needed
+
+undistrR :: Vect k (DSum (Tensor a c) (Tensor b c)) -> Vect k (Tensor (DSum a b) c)
+undistrR v = fmap ( \abc -> case abc of Left (a,c) -> (Left a, c); Right (b,c) -> (Right b, c) ) v
+
+-- For example:
+-- > distrL (e1 `te` i1 e2) :: Vect Q (DSum (Tensor EBasis EBasis) (Tensor EBasis EBasis))
+-- Left (e1,e2)
diff --git a/Math/Algebras/VectorSpace.hs b/Math/Algebras/VectorSpace.hs
--- a/Math/Algebras/VectorSpace.hs
+++ b/Math/Algebras/VectorSpace.hs
@@ -17,7 +17,7 @@
 
 -- |Given a field type k (ie a Fractional instance), Vect k b is the type of the free k-vector space over the basis type b.
 -- Elements of Vect k b consist of k-linear combinations of elements of b.
-data Vect k b = V [(b,k)] deriving (Eq,Ord)
+newtype Vect k b = V [(b,k)] deriving (Eq,Ord)
 
 instance (Num k, Show b) => Show (Vect k b) where
     show (V []) = "0"
diff --git a/Math/QuantumAlgebra/OrientedTangle.hs b/Math/QuantumAlgebra/OrientedTangle.hs
--- a/Math/QuantumAlgebra/OrientedTangle.hs
+++ b/Math/QuantumAlgebra/OrientedTangle.hs
@@ -61,30 +61,27 @@
 
 
 
-
-
-
 idV = id
 idV' = id
 
-evalV  = \(T (E i) (E j)) -> if i + j == 0 then return () else zero
-evalV' = \(T (E i) (E j)) -> if i + j == 0 then return () else zero
+evalV  = \(E i, E j) -> if i + j == 0 then return () else zero
+evalV' = \(E i, E j) -> if i + j == 0 then return () else zero
 
 coevalV  m = foldl (<+>) zero [e i `te` e (-i) | i <- [1..m] ]
 coevalV' m = foldl (<+>) zero [e (-i) `te` e i | i <- [1..m] ]
 
 lambda m = q' ^ m -- q^-m
 
-c m (T (E i) (E j)) = case compare i j of
-                      EQ -> (lambda m * q) *> return (T (E i) (E i))
-                      LT -> lambda m *> return (T (E j) (E i))
-                      GT -> lambda m *> (return (T (E j) (E i)) <+> (q - q') *> return (T (E i) (E j)))
+c m (E i, E j) = case compare i j of
+                      EQ -> (lambda m * q) *> return (E i, E i)
+                      LT -> lambda m *> return (E j, E i)
+                      GT -> lambda m *> (return (E j, E i) <+> (q - q') *> return (E i, E j))
 
 -- inverse of c
-c' m (T (E i) (E j)) = case compare i j of
-                       EQ -> (1/(lambda m * q)) *> return (T (E i) (E i))
-                       LT -> (1/lambda m) *> (return (T (E j) (E i)) <+> (q'-q) *> return (T (E i) (E j)))
-                       GT -> (1/lambda m) *> return (T (E j) (E i))
+c' m (E i, E j) = case compare i j of
+                       EQ -> (1/(lambda m * q)) *> return (E i, E i)
+                       LT -> (1/lambda m) *> (return (E j, E i) <+> (q'-q) *> return (E i, E j))
+                       GT -> (1/lambda m) *> return (E j, E i)
 
 testcc' m v = nf $ v >>= c m >>= c' m
 
@@ -97,15 +94,15 @@
 capRL m = coevalV m
 
 capLR m = do
-    T i j <- coevalV' m
+    (i,j) <- coevalV' m
     k <- mu' m j
-    return (T i k)
+    return (i,k)
 
 cupRL m = evalV
 
-cupLR m (T i j) = do
+cupLR m (i,j) = do
     k <- mu m i
-    evalV' (T k j)    
+    evalV' (k,j)    
 -- linear evalV' . (linear (mu' m) `tf` idV)
 
 
@@ -114,81 +111,80 @@
 
 xminus m = c' m
 
-yplus m (T p q) = do
-    T r s <- capRL m
-    T t u <- xplus m (T q r)
-    cupRL m (T p t)
-    return (T u s)
+yplus m (p,q) = do
+    (r,s) <- capRL m
+    (t,u) <- xplus m (q,r)
+    cupRL m (p,t)
+    return (u,s)
 
-yminus m (T p q) = do
-    T r s <- capRL m
-    T t u <- xminus m (T q r)
-    cupRL m (T p t)
-    return (T u s)
+yminus m (p,q) = do
+    (r,s) <- capRL m
+    (t,u) <- xminus m (q,r)
+    cupRL m (p,t)
+    return (u,s)
 
-tplus m (T p q) = do
-    T r s <- capLR m
-    T t u <- xplus m (T s p)
-    cupLR m (T u q)
-    return (T r t)
+tplus m (p,q) = do
+    (r,s) <- capLR m
+    (t,u) <- xplus m (s,p)
+    cupLR m (u,q)
+    return (r,t)
 
-tminus m (T p q) = do
-    T r s <- capLR m
-    T t u <- xminus m (T s p)
-    cupLR m (T u q)
-    return (T r t)
+tminus m (p,q) = do
+    (r,s) <- capLR m
+    (t,u) <- xminus m (s,p)
+    cupLR m (u,q)
+    return (r,t)
 
-zplus m (T p q) = do
-    T r u <- capLR m
-    T s t <- capLR m
-    T v w <- xplus m (T t u)
-    cupLR m (T v q)
-    cupLR m (T w p)
-    return (T r s)
+zplus m (p,q) = do
+    (r,u) <- capLR m
+    (s,t) <- capLR m
+    (v,w) <- xplus m (t,u)
+    cupLR m (v,q)
+    cupLR m (w,p)
+    return (r,s)
 
-zminus m (T p q) = do
-    T r u <- capLR m
-    T s t <- capLR m
-    T v w <- xminus m (T t u)
-    cupLR m (T v q)
-    cupLR m (T w p)
-    return (T r s)
+zminus m (p,q) = do
+    (r,u) <- capLR m
+    (s,t) <- capLR m
+    (v,w) <- xminus m (t,u)
+    cupLR m (v,q)
+    cupLR m (w,p)
+    return (r,s)
 
 {-
 Then we have for example the following:
 > let v = e1 `te` e2 in nf $ v >>= xplus 2 >>= xminus 2
-T e1 e2
+(e1,e2)
 > let v = e (-1) `te` e2 in nf $ v >>= yplus 2 >>= tminus 2
-T e-1 e2
+(e-1,e2)
 > let v = e (-1) `te` e (-2) in nf $ v >>= zplus 2 >>= zminus 2
-T e-1 e-2
-
+(e-1,e-2)
 -}
 
 
 oloop m = nf $ do
-    T a b <- capLR m
-    cupRL m (T a b)
+    (a,b) <- capLR m
+    cupRL m (a,b)
 
 -- oriented trefoil
 otrefoil m = nf $ do
-    T p q <- capLR m
-    T r s <- capLR m
-    T t u <- tminus m (T q r)
-    T v w <- zminus m (T p t)
-    T x y <- xminus m (T u s)
-    cupRL m (T w x)
-    cupRL m (T v y)
+    (p,q) <- capLR m
+    (r,s) <- capLR m
+    (t,u) <- tminus m (q,r)
+    (v,w) <- zminus m (p,t)
+    (x,y) <- xminus m (u,s)
+    cupRL m (w,x)
+    cupRL m (v,y)
 
 -- oriented the other way
 otrefoil' m = nf $ do
-    T p q <- capRL m
-    T r s <- capRL m
-    T t u <- yminus m (T q r)
-    T v w <- xminus m (T p t)
-    T x y <- zminus m (T u s)
-    cupLR m (T w x)
-    cupLR m (T v y)
+    (p,q) <- capRL m
+    (r,s) <- capRL m
+    (t,u) <- yminus m (q,r)
+    (v,w) <- xminus m (p,t)
+    (x,y) <- zminus m (u,s)
+    cupLR m (w,x)
+    cupLR m (v,y)
 
 
 {-
diff --git a/Math/QuantumAlgebra/QuantumPlane.hs b/Math/QuantumAlgebra/QuantumPlane.hs
--- a/Math/QuantumAlgebra/QuantumPlane.hs
+++ b/Math/QuantumAlgebra/QuantumPlane.hs
@@ -56,7 +56,7 @@
     unit 0 = zero -- V []
     unit x = V [(munit,x)] where munit = Aq20 (NCM 0 [])
     mult x = x''' where
-        x' = mult $ fmap (\(T (Aq20 a) (Aq20 b)) -> T a b) x -- unwrap and multiply
+        x' = mult $ fmap ( \(Aq20 a, Aq20 b) -> (a,b) ) x -- unwrap and multiply
         x'' = x' %% aq20 -- quotient by m2q relations while unwrapped
         x''' = fmap Aq20 x'' -- wrap the monomials up as Aq20 again
 
@@ -78,7 +78,7 @@
     unit 0 = zero -- V []
     unit x = V [(munit,x)] where munit = Aq02 (NCM 0 [])
     mult x = x''' where
-        x' = mult $ fmap (\(T (Aq02 a) (Aq02 b)) -> T a b) x -- unwrap and multiply
+        x' = mult $ fmap ( \(Aq02 a, Aq02 b) -> (a,b) ) x -- unwrap and multiply
         x'' = x' %% aq02 -- quotient by m2q relations while unwrapped
         x''' = fmap Aq02 x'' -- wrap the monomials up as Aq02 again
 
@@ -102,7 +102,7 @@
     unit 0 = zero -- V []
     unit x = V [(munit,x)] where munit = M2q (NCM 0 [])
     mult x = x''' where
-        x' = mult $ fmap (\(T (M2q a) (M2q b)) -> T a b) x -- unwrap and multiply
+        x' = mult $ fmap ( \(M2q a, M2q b) -> (a,b) ) x -- unwrap and multiply
         x'' = x' %% m2q -- quotient by m2q relations while unwrapped
         x''' = fmap M2q x'' -- wrap the monomials up as M2q again
 
@@ -191,7 +191,7 @@
     unit 0 = zero -- V []
     unit x = V [(munit,x)] where munit = SL2q (NCM 0 [])
     mult x = x''' where
-        x' = mult $ fmap (\(T (SL2q a) (SL2q b)) -> T a b) x -- unwrap and multiply
+        x' = mult $ fmap ( \(SL2q a, SL2q b) -> (a,b) ) x -- unwrap and multiply
         x'' = x' %% sl2q -- quotient by sl2q relations while unwrapped
         x''' = fmap SL2q x'' -- wrap the monomials up as SL2q again
 
@@ -228,9 +228,9 @@
 -- This is a Yang-Baxter operator, but not the only possible such
 -- Street, p93
 yb x = nf $ x >>= yb' where
-    yb' (T a b) = case compare a b of
-                 GT -> return (T b a)
-                 LT -> return (T b a) + unit (q-q') * return (T a b)
-                 EQ -> unit q * return (T a a)
+    yb' (a,b) = case compare a b of
+                 GT -> return (b,a)
+                 LT -> return (b,a) + unit (q-q') * return (a,b)
+                 EQ -> unit q * return (a,a)
 
 
diff --git a/Math/QuantumAlgebra/Tangle.hs b/Math/QuantumAlgebra/Tangle.hs
--- a/Math/QuantumAlgebra/Tangle.hs
+++ b/Math/QuantumAlgebra/Tangle.hs
@@ -27,7 +27,7 @@
 instance (Num k, Ord a) => Algebra k [a] where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
-    mult = nf . fmap (\(T a b) -> a `mmult` b)
+    mult = nf . fmap (\(a,b) -> a `mmult` b)
 
 -- Could make TensorAlgebra k a into an instance of Category, TensorCategory
     
diff --git a/Math/Test/TAlgebras/TStructures.hs b/Math/Test/TAlgebras/TStructures.hs
--- a/Math/Test/TAlgebras/TStructures.hs
+++ b/Math/Test/TAlgebras/TStructures.hs
@@ -14,19 +14,9 @@
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
 import Math.Algebras.Structures -- what we're testing
--- import MathExperiments.Algebra.MonoidAlgebra
--- import MathExperiments.Algebra.Examples
 
-{-
-prop_VectorSpace (k,l,x,y,z) =
-    smultL k (smultL l x) == smultL (k*l) x &&
-    add x y == add y z &&
-    add x (add y z) == add (add x y) z &&
-    add x zero == x &&
-    add zero x == x
--- !! check definition - have I forgotten anything - yes, additive inverses
--}
 
+
 prop_Linear f (k,l,x,y) =
     f (add (smultL k x) (smultL l y)) == add (smultL k (f x)) (smultL l (f y))
 -- now use this to show algebra and coalgebra ops are linear
@@ -116,14 +106,14 @@
 
 prop_Bialgebra2 (k,xy) =
     (comult . unit') k' + xy == ((unit' `tf` unit') . iso) k' + xy
-    where iso = fmap (\ () -> T () () ) -- the isomorphism k ~= k tensor k
+    where iso = fmap (\ () -> ((),()) ) -- the isomorphism k ~= k tensor k
           k' = unit k :: Trivial Integer -- inject into the trivial algebra
 -- the +xy is just to force the other expression to be of the right type
 
 prop_Bialgebra3 (x,y) =
     (counit' . mult) xy == (iso . (counit' `tf` counit')) xy
     where xy = x `te` y
-          iso = fmap (\(T () ()) -> ())
+          iso = fmap ( \((),()) -> ())
 
 prop_Bialgebra4 (k,x) =
     id k == (counit . (\a -> a+x-x) . unit) k
@@ -161,34 +151,25 @@
 -}
 
 
-
-
-
-
-
-
-
-
-
-
+-- FROBENIUS ALGEBRAS
 
 frobeniusLeft1 = (id `tf` mult) . assocR . (comult `tf` id)
 
 frobeniusLeft2 x = nf $ x >>= fl
-    where fl (T i j) = do
-              T k l <- comultM i
+    where fl (i,j) = do
+              (k,l) <- comultM i
               m <- idM j
               p <- idM k
-              q <- multM (T l m)
-              return (T p q)
+              q <- multM (l,m)
+              return (p,q)
 
 frobeniusMiddle1 = comult . mult
 
 frobeniusMiddle2 x = nf $ x >>= fm
-    where fm (T i j) = do
-              k <- multM (T i j)
-              T l m <- comultM k
-              return (T l m)
+    where fm (i,j) = do
+              k <- multM (i,j)
+              (l,m) <- comultM k
+              return (l,m)
 
 prop_FrobeniusRelation (x,y) =
     let xy = x `te` y
@@ -204,8 +185,3 @@
 -- unit takes k as input, so isn't in the monad
 -- counit gives k as output - what would we do with it
 -- so perhaps we have to use unit' and counit'
-
-
-
-
-
diff --git a/Math/Test/TAlgebras/TTensorProduct.hs b/Math/Test/TAlgebras/TTensorProduct.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TAlgebras/TTensorProduct.hs
@@ -0,0 +1,87 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE EmptyDataDecls, ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies, RankNTypes #-}
+
+module Math.Test.TAlgebras.TTensorProduct where
+
+import Test.QuickCheck
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebra.Field.Base
+import Math.Test.TAlgebras.TVectorSpace hiding (i1, i2)
+
+import Prelude as P
+import Control.Category as C
+import Control.Arrow
+
+type DirectSum k u v =
+    (u ~ Vect k a, v ~ Vect k b) => Vect k (DSum a b)
+
+type TensorProd k u v =
+    (u ~ Vect k a, v ~ Vect k b) => Vect k (Tensor a b)
+
+type En = Vect Q EBasis
+
+{-
+-- But then you need to make sure that you run GHCi with -XTypeFamilies, otherwise:
+
+> e1 `te` e2 :: TensorProd Q En En
+<interactive>:1:0:
+    Illegal equational constraint En ~ Vect Q a
+    (Use -XTypeFamilies to permit this)
+    In an expression type signature: TensorProd Q En En
+    In the expression: e1 `te` e2 :: TensorProd Q En En
+    In the definition of `it': it = e1 `te` e2 :: TensorProd Q En En
+-}
+
+
+-- Now test eg
+-- > quickCheck (\x -> (distrL . undistrL) x == id x)
+-- but need to make x be of interesting type (not just () )
+
+
+data Zero
+-- a type with no inhabitants
+-- so the associated free vector space is the zero space
+
+-- instance Eq Zero where {}
+-- instance Ord Zero where {}
+instance Show Zero where {}
+
+-- > zero :: Vect Q Zero
+-- 0
+
+
+-- ARROW INSTANCE
+-- This isn't currently used anywhere else
+-- It's intended to illustrate the point that tensor product is like doing things in parallel
+
+newtype Linear k a b = Linear (Vect k a -> Vect k b)
+
+instance Category (Linear k) where
+    id = Linear P.id
+    (Linear f) . (Linear g) = Linear (f P.. g)
+
+instance Num k => Arrow (Linear k) where
+    arr f = Linear (fmap f) -- requires nf call afterwards
+    first (Linear f) = Linear $ \(V ts) -> V $
+        concat [let V us = x *> te (f $ return a) (return c) in us | ((a,c),x) <- ts]
+    second (Linear f) = Linear $ \(V ts) -> V $
+        concat [let V us = x *> te (return c) (f $ return a) in us | ((c,a),x) <- ts]
+    Linear f *** Linear g = Linear (f `tf2` g)
+        where tf2 f g (V ts) = V $ concat
+                  [let V us = x *> te (f $ return a) (g $ return b) in us | ((a,b), x) <- ts]
+        -- can't use tf, as it uses add, which assumes Ord instance
+        -- hence we should call nf afterwards
+    -- !! What about &&&
+
+{-
+-- The following are morally correct, but don't work because they require Ord instance
+instance Num k => ArrowChoice (Linear k) where
+    left (Linear f) = Linear (f `dsume` id)
+    right (Linear f) = Linear (id `dsume` f)
+    Linear f +++ Linear g = Linear (f `dsumf` g)
+    Linear f ||| Linear g = Linear (f `coprodf` g)
+-}
+
diff --git a/Math/Test/TAlgebras/TVectorSpace.hs b/Math/Test/TAlgebras/TVectorSpace.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TAlgebras/TVectorSpace.hs
@@ -0,0 +1,208 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE FlexibleInstances, NoMonomorphismRestriction, ScopedTypeVariables, GeneralizedNewtypeDeriving #-}
+
+
+module Math.Test.TAlgebras.TVectorSpace where
+
+import Test.QuickCheck
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebra.Field.Base
+
+-- import Control.Monad -- MonadPlus
+
+
+prop_AddGrp (x,y,z) =
+    x <+> (y <+> z) == (x <+> y) <+> z && -- associativity
+    x <+> y == y <+> x                 && -- commutativity
+    x <+> zero == x                    && -- identity
+    x <+> neg x == zero                   -- inverse
+
+prop_VecSp (a,b,x,y,z) =
+    prop_AddGrp (x,y,z) &&
+    a *> (x <+> y) == a *> x <+> a *> y && -- distributivity through vectors
+    (a+b) *> x == a *> x <+> b *> x     && -- distributivity through scalars
+    (a*b) *> x == a *> (b *> x)         && -- associativity
+    1 *> x == x                            -- unit
+
+instance Arbitrary EBasis where
+    arbitrary = do n <- arbitrary :: Gen Int
+                   return (E n)
+
+instance Arbitrary Q where
+    arbitrary = do n <- arbitrary :: Gen Integer
+                   d <- arbitrary :: Gen Integer
+                   return (if d == 0 then fromInteger n else fromInteger n / fromInteger d)
+
+instance (Num k, Ord b, Arbitrary k, Arbitrary b) => Arbitrary (Vect k b) where
+    arbitrary = do ts <- arbitrary :: Gen [(b, k)] -- ScopedTypeVariables
+                   return $ nf $ V ts
+
+prop_VecSpQn (a,b,x,y,z) = prop_VecSp (a,b,x,y,z)
+    where types = (a,b,x,y,z) :: (Q, Q, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis)
+
+
+prop_Linear f (a,x,y) =
+    f (x <+> y) == f x <+> f y &&
+    f zero == zero             &&
+    f (neg x) == neg (f x)     &&
+    f (a *> x) == a *> f x
+
+prop_LinearQn f (a,x,y) = prop_Linear f (a,x,y)
+    where types = (a,x,y) :: (Q, Vect Q EBasis, Vect Q EBasis)
+
+
+newtype FBasis = F Int deriving (Eq,Ord,Arbitrary)
+
+instance Show FBasis where show (F i) = "f" ++ show i
+
+f i = return (F i) :: Vect Q FBasis
+f1 = f 1
+f2 = f 2
+f3 = f 3
+
+
+-- DIRECT SUM
+
+{-
+instance Num k => MonadPlus (Vect k) where
+    mzero = zero
+    mplus (V xs) (V ys) = V (xs++ys) -- need to call nf afterwards
+-}
+
+
+-- (Alternative versions of prodf and coprodf)
+
+f .*. g = linear fg' where
+    fg' b = fmap Left (f (return b)) <+> fmap Right (g (return b))
+
+f .+. g = linear fg' where
+    fg' (Left a) = f (return a)
+    fg' (Right b) = g (return b)
+
+
+type LinFun k a b = [(a, Vect k b)]
+-- a way of representing a linear function as data
+
+linfun :: (Eq a, Ord b, Num k) => LinFun k a b -> Vect k a -> Vect k b
+linfun avbs = linear f where
+    f a = case lookup a avbs of
+          Just vb -> vb
+          Nothing -> zero
+
+
+prop_Product (f',g',x) =
+    f x == (p1 . fg) x &&
+    g x == (p2 . fg) x
+    where f = linfun f'
+          g = linfun g'
+          fg = prodf f g
+
+prop_Coproduct (f',g',a,b) =
+    f a == (fg . i1) a &&
+    g b == (fg . i2) b
+    where f = linfun f'
+          g = linfun g'
+          fg = coprodf f g
+
+prop_dsumf (f',g',a,b) =
+    f a == (p1 . fg . i1) a &&
+    g b == (p2 . fg . i2) b
+    where f = linfun f'
+          g = linfun g'
+          fg = dsumf f g
+
+
+newtype ABasis = A Int deriving (Eq,Ord,Show,Arbitrary) -- GeneralizedNewtypeDeriving
+newtype BBasis = B Int deriving (Eq,Ord,Show,Arbitrary)
+newtype SBasis = S Int deriving (Eq,Ord,Show,Arbitrary)
+newtype TBasis = T Int deriving (Eq,Ord,Show,Arbitrary)
+
+prop_ProductQn (f,g,x) = prop_Product (f,g,x)
+    where types = (f,g,x) :: (LinFun Q SBasis ABasis, LinFun Q SBasis BBasis, Vect Q SBasis)
+
+prop_CoproductQn (f,g,a,b) = prop_Coproduct (f,g,a,b)
+    where types = (f,g,a,b) :: (LinFun Q ABasis TBasis, LinFun Q BBasis TBasis, Vect Q ABasis, Vect Q BBasis)
+
+prop_dsumfQn (f,g,a,b) = prop_dsumf (f,g,a,b)
+    where types = (f,g,a,b) :: (LinFun Q ABasis SBasis, LinFun Q BBasis TBasis, Vect Q ABasis, Vect Q BBasis)
+
+
+-- TENSOR PRODUCT
+
+dot0 uv = sum [ if a == b then x*y else 0 | (a,x) <- u, (b,y) <- v]
+    where V u = p1 uv
+          V v = p2 uv
+
+dot1 uv = nf $ V [( (), if a == b then x*y else 0) | (a,x) <- u, (b,y) <- v]
+    where V u = p1 uv
+          V v = p2 uv
+
+polymult1 uv = nf $ V [(E (i+j) , x*y) | (E i,x) <- u, (E j,y) <- v]
+    where V u = p1 uv
+          V v = p2 uv
+
+{-
+tensor1 :: (Num k, Ord a, Ord b) => (Vect k a, Vect k b) -> Vect k (a, b)
+tensor1 (V axs, V bys) = nf $ V [((a,b),x*y) | (a,x) <- axs, (b,y) <- bys] 
+
+bilinear1 :: (Num k, Ord a, Ord b, Ord c) =>
+     ((a, b) -> Vect k c) -> (Vect k a, Vect k b) -> Vect k c
+bilinear1 f = linear f . tensor1
+
+prop_Bilinear1 f (a,u1,u2,v1,v2) =
+    prop_Linear (\v -> f (u1,v)) (a,v1,v2) &&
+    prop_Linear (\u -> f (u,v1)) (a,u1,u2)
+
+prop_BilinearQn1 f (a,u1,u2,v1,v2) = prop_Bilinear1 f (a,u1,u2,v1,v2)
+    where types = (a,u1,u2,v1,v2) :: (Q, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis)
+-}
+
+tensor :: (Num k, Ord a, Ord b) => Vect k (Either a b) -> Vect k (a, b)
+tensor uv = nf $ V [( (a,b), x*y) | (a,x) <- u, (b,y) <- v]
+    where V u = p1 uv; V v = p2 uv
+
+bilinear :: (Num k, Ord a, Ord b, Ord c) =>
+    ((a, b) -> Vect k c) -> Vect k (Either a b) -> Vect k c
+bilinear f = linear f . tensor
+
+dot = bilinear (\(a,b) -> if a == b then return () else zero)
+
+polymult = bilinear (\(E i, E j) -> return (E (i+j)))
+
+prop_Bilinear :: (Num k, Ord a, Ord b, Ord t) =>
+     (Vect k (Either a b) -> Vect k t) -> (k, Vect k a, Vect k a, Vect k b, Vect k b) -> Bool
+prop_Bilinear f (a,u1,u2,v1,v2) =
+    prop_Linear (\v -> f (u1 `dsume` v)) (a,v1,v2) &&
+    prop_Linear (\u -> f (u `dsume` v1)) (a,u1,u2)
+
+prop_BilinearQn f (a,u1,u2,v1,v2) = prop_Bilinear f (a,u1,u2,v1,v2)
+    where types = (a,u1,u2,v1,v2) :: (Q, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis)
+
+{-
+> quickCheck (prop_BilinearQn dot1)
++++ OK, passed 100 tests.
+> quickCheck (prop_BilinearQn polymult1)
++++ OK, passed 100 tests.
+*Math.Test.TAlgebras.TVectorSpace> quickCheck (prop_BilinearQn tensor)
++++ OK, passed 100 tests.
+
+> quickCheck (\x -> dot1 x == dot x)
++++ OK, passed 100 tests.
+> quickCheck (\x -> polymult1 x == polymult x)
++++ OK, passed 100 tests.
+
+
+> quickCheck (prop_BilinearQn id)
+*** Failed! Falsifiable (after 2 tests):  
+(1,0,0,e1,0)
+-- fails basically because (0 <+> 0) `dsume` e0 /= (0 `dsume` e0) <+> (0 `dsume` e0)
+
+>  (zero <+> zero) `dsume` e1
+Right e1
+> (zero `dsume` e1) <+> (zero `dsume` e1)
+2Right e1
+
+-}
+