diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal
--- a/HaskellForMaths.cabal
+++ b/HaskellForMaths.cabal
@@ -1,5 +1,5 @@
    Name:                HaskellForMaths
-   Version:             0.2.2
+   Version:             0.3.1
    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
@@ -7,7 +7,7 @@
    License-file:        license.txt
    Author:              David Amos
    Maintainer:          haskellformaths-at-gmail-dot-com
-   Homepage:            http://www.polyomino.f2s.com/haskellformathsv2/HaskellForMathsv2.html
+   Homepage:            http://haskellformaths.blogspot.com/
    Build-Type:          Simple
    Cabal-Version:       >=1.2
 
@@ -22,9 +22,12 @@
         Math/Test/TRootSystem.hs,
         Math/Test/TSubquotients.hs,
         Math/Test/TestAll.hs
+        Math/Test/TAlgebras/TStructures.hs
+        Math/Test/TAlgebras/TQuaternions.hs
+        Math/Test/TAlgebras/TGroupAlgebra.hs
 
    Library
-     Build-Depends:     base >=3 && < 4, containers, array, random, QuickCheck
+     Build-Depends:     base >= 2 && < 5, containers, array, random, QuickCheck
      Exposed-modules:
         Math.Algebra.LinearAlgebra,
         Math.Algebra.Commutative.Monomial,
@@ -34,6 +37,10 @@
         Math.Algebra.Group.RandomSchreierSims, Math.Algebra.Group.Subquotients,
         Math.Algebra.Group.StringRewriting, Math.Algebra.Group.CayleyGraph,
         Math.Algebra.NonCommutative.NCPoly, Math.Algebra.NonCommutative.GSBasis, Math.Algebra.NonCommutative.TensorAlgebra,
+        Math.Algebras.AffinePlane, Math.Algebras.Commutative, Math.Algebras.GroupAlgebra,
+        Math.Algebras.LaurentPoly, Math.Algebras.Matrix, Math.Algebras.NonCommutative,
+        Math.Algebras.Quaternions, Math.Algebras.Structures, Math.Algebras.TensorAlgebra,
+        Math.Algebras.TensorProduct, Math.Algebras.VectorSpace,
         Math.Combinatorics.Graph, Math.Combinatorics.GraphAuts, Math.Combinatorics.StronglyRegularGraph,
         Math.Combinatorics.Design, Math.Combinatorics.FiniteGeometry, Math.Combinatorics.Hypergraph,
         Math.Combinatorics.LatinSquares,
@@ -42,6 +49,8 @@
         Math.Projects.Rubik, Math.Projects.MiniquaternionGeometry,
         Math.Projects.ChevalleyGroup.Classical, Math.Projects.ChevalleyGroup.Exceptional,
         Math.Projects.KnotTheory.Braid,
-        Math.Projects.KnotTheory.LaurentMPoly, Math.Projects.KnotTheory.TemperleyLieb, Math.Projects.KnotTheory.IwahoriHecke
+        Math.Projects.KnotTheory.LaurentMPoly, Math.Projects.KnotTheory.TemperleyLieb, Math.Projects.KnotTheory.IwahoriHecke,
+        Math.QuantumAlgebra.OrientedTangle,
+        Math.QuantumAlgebra.QuantumPlane, Math.QuantumAlgebra.Tangle, Math.QuantumAlgebra.TensorCategory
 
      ghc-options:       -w
diff --git a/Math/Algebra/Commutative/MPoly.hs b/Math/Algebra/Commutative/MPoly.hs
--- a/Math/Algebra/Commutative/MPoly.hs
+++ b/Math/Algebra/Commutative/MPoly.hs
@@ -2,6 +2,7 @@
 
 {-# OPTIONS_GHC -fglasgow-exts #-}
 
+-- |A module providing a type for (commutative) multivariate polynomials, with support for various term orders.
 module Math.Algebra.Commutative.MPoly where
 
 import qualified Data.Map as M
diff --git a/Math/Algebra/Group/PermutationGroup.hs b/Math/Algebra/Group/PermutationGroup.hs
--- a/Math/Algebra/Group/PermutationGroup.hs
+++ b/Math/Algebra/Group/PermutationGroup.hs
@@ -8,12 +8,14 @@
 
 import Math.Common.ListSet (toListSet, union, (\\) ) -- a version of union which assumes the arguments are ascending sets (no repeated elements)
 
+infix 8 ^-, ~^
+
 rotateL (x:xs) = xs ++ [x]
 
 
 -- PERMUTATIONS
 
--- |Type for permutations, considered as group elements.
+-- |A type for permutations, considered as functions or actions which can be performed on an underlying set.
 newtype Permutation a = P (M.Map a a) deriving (Eq,Ord)
 
 fromPairs xys | isValid   = fromPairs' xys
@@ -41,7 +43,7 @@
            Just y  -> y
            Nothing -> x -- if x `notElem` supp (P g), then x is not moved
 
--- |b -^ g returns the image of an edge or block b under the action of the permutation g
+-- |b -^ g returns the image of an edge or block b under the action of the permutation g.
 -- The dash is meant to be a mnemonic for edge or line or block.
 (-^) :: (Ord t) => [t] -> Permutation t -> [t]
 xs -^ g = L.sort [x .^ g | x <- xs]
diff --git a/Math/Algebra/Group/Subquotients.hs b/Math/Algebra/Group/Subquotients.hs
--- a/Math/Algebra/Group/Subquotients.hs
+++ b/Math/Algebra/Group/Subquotients.hs
@@ -58,6 +58,7 @@
 -}
 
 
+isTransitive :: (Ord t) => [Permutation t] -> Bool
 isTransitive gs = length (orbits gs) == 1
 
 
@@ -75,8 +76,12 @@
 -}
 
 -- Seress p81
--- A transitive constituent homomorphism is the restriction of G <= Sym(omega) to an orbit delta <= omega
--- This function returns the kernel and the image
+-- |Given a group gs and a transitive constituent ys, return the kernel and image of the transitive constituent homomorphism.
+-- That is, suppose that gs acts on a set xs, and ys is a subset of xs on which gs acts transitively.
+-- Then the transitive constituent homomorphism is the restriction of the action of gs to an action on the ys.
+transitiveConstituentHomomorphism
+  :: (Ord a, Show a) =>
+     [Permutation a] -> [a] -> ([Permutation a], [Permutation a])
 transitiveConstituentHomomorphism gs delta
     | delta == closure delta [(.^ g) | g <- gs] -- delta is closed under action of gs, hence a union of orbits
         = transitiveConstituentHomomorphism' gs delta
@@ -112,14 +117,17 @@
 -- Because the support of the permutations is not constrained to be [1..n], we have to use a map instead of an array
 -- This probably affects the complexity, but isn't a problem in practice
 
--- Find all block systems
+-- |Given a transitive group gs, find all non-trivial block systems. That is, if gs act on xs,
+-- find all the ways that the xs can be divided into blocks, such that the gs also have a permutation action on the blocks
+blockSystems :: (Ord t) => [Permutation t] -> [[[t]]]
 blockSystems gs
     | isTransitive gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- xs ]
     | otherwise = error "blockSystems: not transitive"
     where x:xs = foldl union [] $ map supp gs
 
 
--- More efficient version if we have an sgs
+-- |A more efficient version of blockSystems, if we have an sgs
+blockSystemsSGS :: (Ord a) => [Permutation a] -> [[[a]]]
 blockSystemsSGS gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- rs ]
     where x:xs = foldl union [] $ map supp gs
           hs = filter (\g -> x < minsupp g) gs -- sgs for stabiliser Gx
@@ -133,14 +141,21 @@
 -- see Holt RandomStab function
 
 
+-- |A permutation group is primitive if it has no non-trivial block systems
+isPrimitive :: (Ord t) => [Permutation t] -> Bool
 isPrimitive gs = null (blockSystems gs)
 
+isPrimitiveSGS :: (Ord a) => [Permutation a] -> Bool
 isPrimitiveSGS gs = null (blockSystemsSGS gs)
 
 -- There are other optimisations we haven't done
 -- see Holt p86
 
-
+-- |Given a transitive group gs, and a block system for gs, return the kernel and image of the block homomorphism
+-- (the homomorphism onto the action of gs on the blocks)
+blockHomomorphism
+  :: (Ord t, Show t) =>
+     [Permutation t] -> [[t]] -> ([Permutation t], [Permutation [t]])
 blockHomomorphism gs bs
     | bs == closure bs [(-^ g) | g <- gs] -- bs is closed under action of gs
         = blockHomomorphism' gs bs
diff --git a/Math/Algebra/NonCommutative/NCPoly.hs b/Math/Algebra/NonCommutative/NCPoly.hs
--- a/Math/Algebra/NonCommutative/NCPoly.hs
+++ b/Math/Algebra/NonCommutative/NCPoly.hs
@@ -1,5 +1,6 @@
 -- Copyright (c) David Amos, 2008. All rights reserved.
 
+-- |A module providing a type for non-commutative polynomials.
 module Math.Algebra.NonCommutative.NCPoly where
 
 import Data.List as L
@@ -14,10 +15,13 @@
     compare (M xs) (M ys) = compare (length xs,xs) (length ys,ys)
 -- Glex ordering
 
-instance Show v => Show (Monomial v) where
+instance (Eq v, Show v) => Show (Monomial v) where
     show (M xs) | null xs = "1"
-                | otherwise = concatMap show xs
--- !! we can do better - we should show "xx" as "x^2"
+                | otherwise = concatMap showPower (L.group xs)
+        where showPower [v] = showVar v
+              showPower vs@(v:_) = showVar v ++ "^" ++ show (length vs)
+              showVar v = filter (/= '"') (show v)
+-- Taken from NonComMonomial - why don't we just use it directly
 
 instance (Eq v, Show v) => Num (Monomial v) where
     M xs * M ys = M (xs ++ ys)
@@ -54,7 +58,6 @@
                                  then "+(" ++ c:cs ++ ")"
                                  else if c == '-' then c:cs else '+':c:cs
 
-
 instance (Ord v, Show v, Num r) => Num (NPoly r v) where
     NP ts + NP us = NP (mergeTerms ts us)
     negate (NP ts) = NP $ map (\(m,c) -> (m,-c)) ts
@@ -97,6 +100,9 @@
     show Y = "y"
     show Z = "z"
 
+-- |Create a non-commutative variable for use in forming non-commutative polynomials.
+-- For example, we could define x = var "x", y = var "y". Then x*y /= y*x.
+var :: (Num k) => v -> NPoly k v
 var v = NP [(M [v], 1)]
 
 x = var X :: NPoly Q Var
diff --git a/Math/Algebras/AffinePlane.hs b/Math/Algebras/AffinePlane.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/AffinePlane.hs
@@ -0,0 +1,90 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+
+-- |A module defining the affine plane and its symmetries
+module Math.Algebras.AffinePlane where
+
+import Math.Algebra.Field.Base hiding (powers)
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+import Math.Algebras.Commutative
+
+
+data XY = X | Y deriving (Eq, Ord)
+
+instance Show XY where show X = "x"; show Y = "y"
+
+x = glexVar X :: GlexPoly Q XY
+y = glexVar Y :: GlexPoly Q XY
+
+
+data ABCD = A | B | C | D deriving (Eq, Ord)
+
+instance Show ABCD where show A = "a"; show B = "b"; show C = "c"; show D = "d"
+
+a,b,c,d :: Monomial m => Vect Q (m ABCD)
+a = var A
+b = var B
+c = var C
+d = var D
+
+
+-- SL2
+
+newtype SL2 v = SL2 (GlexMonomial v) deriving (Eq,Ord)
+
+instance Show v => Show (SL2 v) where show (SL2 m) = show m
+
+instance Algebra Q (SL2 ABCD) where -- to do this for Num k instead of Q we would need a,b,c,d defined for Num k
+    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'' = 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
+
+sl2Var v = V [(SL2 (Glex 1 [(v,1)]), 1)] -- :: Vect Q (SL2 ABCD)
+
+
+-- For example:
+-- > a*d :: Vect Q (SL2 ABCD)
+-- bc+1
+
+instance Monomial SL2 where
+    var = sl2Var
+    powers (SL2 (Glex _ xis)) = xis
+
+
+
+instance Coalgebra Q (SL2 ABCD) where
+    counit x = case x `bind` cu of
+               V [] -> 0
+               V [(SL2 (Glex 0 []), c)] -> c
+        where cu A = 1 :: Vect Q (SL2 ABCD)
+              cu B = 0
+              cu C = 0
+              cu D = 1
+    comult x = x `bind` cm
+        where cm A = a `te` a + b `te` c
+              cm B = a `te` b + b `te` d
+              cm C = c `te` a + d `te` c
+              cm D = c `te` b + d `te` d
+-- In other words
+-- counit (a b) = (1 0)
+--        (c d)   (0 1)
+-- comult (a b) = (a1 b1) `te` (a2 b2)
+--        (c d)   (c1 d1)      (c2 d2)
+
+instance Bialgebra Q (SL2 ABCD) where {}
+
+instance HopfAlgebra Q (SL2 ABCD) where
+    antipode x = x `bind` antipode'
+        where antipode' A = d
+              antipode' B = b
+              antipode' C = c
+              antipode' D = a
+-- in the GL2 case we would need 1/det factor as well
+
diff --git a/Math/Algebras/Commutative.hs b/Math/Algebras/Commutative.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/Commutative.hs
@@ -0,0 +1,160 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+
+-- |A module defining the algebra of commutative polynomials over a field k
+module Math.Algebras.Commutative where
+
+import Math.Algebra.Field.Base hiding (powers)
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+
+
+-- GLEX MONOMIALS
+
+data GlexMonomial v = Glex Int [(v,Int)] deriving (Eq)
+-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons
+
+-- type GlexMonomialS = GlexMonomial String
+
+instance Ord v => Ord (GlexMonomial v) where
+    compare (Glex si xis) (Glex sj yjs) = compare (-si, xis) (-sj, yjs)
+
+instance Show v => Show (GlexMonomial v) where
+    show (Glex _ []) = "1"
+    show (Glex _ xis) = concatMap (\(x,i) -> if i==1 then showVar x else showVar x ++ "^" ++ show i) xis
+        where showVar x = filter ( /= '"' ) (show x) -- in case v == String
+
+{-
+-- GlexMonomial is a functor and a monad
+-- However, this isn't all that much use, and to make proper use of it we'd need a "nf" function
+-- So leaving this commented out
+
+-- map one basis to another
+instance Functor GlexMonomial where
+    fmap f (Glex si xis) = Glex si [(f x, i) | (x,i) <- xis]
+-- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards
+
+-- GlexMonomial is the free commutative monoid, and hence a monad
+instance Monad GlexMonomial where
+    return x = Glex 1 [(x,1)]
+    (Glex _ xis) >>= f = let parts = [(i, sj, yjs) | (x,i) <- xis, let Glex sj yjs = f x]
+                         in Glex (sum [i*sj | (i,sj,_) <- parts])
+                                 (concatMap (\(i,_,yjs)->map (\(y,j)->(y,i*j)) yjs) parts)
+    -- this isn't really much use - it's variable substitution, but we're only allowed to substitute monomials for each var
+-- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards
+-}
+
+-- This is the monoid algebra for commutative monomials (which are the free commutative monoid)
+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)
+        where mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs
+
+{-
+-- This is just the Set Coalgebra, so better to use a generic instance
+-- Also, not used anywhere. Hence commented out
+instance Num k => Coalgebra k (GlexMonomial v) where
+    counit (V ts) = sum [x | (m,x) <- ts] -- trace
+    comult = fmap (\m -> T m m) -- diagonal
+    -- comult (V ts) = V [(T m m, x) | (m, x) <- ts]
+-}
+type GlexPoly k v = Vect k (GlexMonomial v)
+
+
+
+glexVar v = V [(Glex 1 [(v,1)], 1)]
+
+
+class Monomial m where
+    var :: v -> Vect Q (m v)
+    powers :: m v -> [(v,Int)]
+
+-- |In effect, we have (Num k, Monomial m) => Monad (\v -> Vect k (m v)), with return = var, and (>>=) = bind.
+-- However, we can't express this directly in Haskell, firstly because of the Ord b constraint,
+-- secondly because Haskell doesn't support type functions.
+bind :: (Monomial m, Num k, Ord b, Show b, Algebra k b) =>
+     Vect k (m v) -> (v -> Vect k b) -> Vect k b
+V ts `bind` f = sum [c `smultL` product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] 
+-- V ts `bind` f = sum [product [f x ^ i | (x,i) <- powers m] * unit c | (m, c) <- ts] 
+
+
+instance Monomial GlexMonomial where
+    var = glexVar
+    powers (Glex _ xis) = xis
+
+
+-- DIVISION
+
+lt (V (t:ts)) = t
+
+class DivisionBasis b where
+    dividesB :: b -> b -> Bool
+    divB :: b -> b -> b
+
+dividesT (b1,x1) (b2,x2) = dividesB b1 b2
+divT (b1,x1) (b2,x2) = (divB b1 b2, x1/x2)
+
+-- given f, gs, find as, r such that f = sum (zipWith (*) as gs) + r, with r not divisible by any g
+quotRemMP f gs = quotRemMP' f (replicate n 0, 0) where
+    n = length gs
+    quotRemMP' 0 (us,r) = (us,r)
+    quotRemMP' h (us,r) = divisionStep h (gs,[],us,r)
+    divisionStep h (g:gs,us',u:us,r) =
+        if lt g `dividesT` lt h
+        then let t = V [lt h `divT` lt g]
+                 h' = h - t*g
+                 u' = u+t
+             in quotRemMP' h' (reverse us' ++ u':us, r)
+        else divisionStep h (gs,u:us',us,r)
+    divisionStep h ([],us',[],r) =
+        let (lth,h') = splitlt h
+        in quotRemMP' h' (reverse us', r+lth)
+    splitlt (V (t:ts)) = (V [t], V ts)
+
+infixl 7 %%
+
+(%%) :: (Fractional k, Ord b, Show b, Algebra k b, DivisionBasis b)
+     => Vect k b -> [Vect k b] -> Vect k b
+f %% gs = r where (_,r) = quotRemMP f gs
+
+
+instance Ord v => DivisionBasis (GlexMonomial v) where
+    dividesB (Glex si xis) (Glex sj yjs) = si <= sj && dividesB' xis yjs where
+        dividesB' ((x,i):xis) ((y,j):yjs) =
+            case compare x y of
+            LT -> False
+            GT -> dividesB' ((x,i):xis) yjs
+            EQ -> if i<=j then dividesB' xis yjs else False
+        dividesB' [] _ = True
+        dividesB' _ [] = False
+    divB (Glex si xis) (Glex sj yjs) = Glex (si-sj) $ divB' xis yjs where
+        divB' ((x,i):xis) ((y,j):yjs) =
+            case compare x y of
+            LT -> (x,i) : divB' xis ((y,j):yjs)
+            EQ -> if i == j then divB' xis yjs else (x,i-j) : divB' xis yjs -- we don't bother to check i > j
+            GT -> error "divB'" -- (y,-j) : divB' ((x,i):xis) yjs
+        divB' xis [] = xis
+        divB' [] yjs = error "divB'"
+
+{-
+-- Need to thread this through Maybe properly, so perhaps use do notation
+divB2 (Glex si xis) (Glex sj yjs)
+    | si < sj = Nothing
+    | otherwise = case divB' xis yjs of
+                  Nothing -> Nothing
+                  Just zks -> Glex (si-sj) zks
+    where divB' ((x,i):xis) ((y,j):yjs) =
+              case compare x y of
+              LT -> (x,i) : divB' xis ((y,j):yjs)
+              EQ -> case compare i j of
+                    LT -> Nothing
+                    EQ -> divB' xis yjs
+                    GT -> (x,i-j) : divB' xis yjs
+              GT -> Nothing
+-}
+-- !! could change divB to return Maybe, and avoid need for dividesB
+
+
diff --git a/Math/Algebras/GroupAlgebra.hs b/Math/Algebras/GroupAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/GroupAlgebra.hs
@@ -0,0 +1,54 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances #-}
+
+module Math.Algebras.GroupAlgebra where
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+
+import Math.Algebra.Group.PermutationGroup hiding (action)
+
+import Math.Algebra.Field.Base
+
+
+instance Mon (Permutation Int) where
+    munit = 1
+    mmult = (*)
+
+type GroupAlgebra k = Vect k (Permutation Int)
+
+-- Monoid Algebra instance
+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)
+
+-- 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
+
+instance Num k => Bialgebra k (Permutation Int) where {}
+-- should check that the algebra and coalgebra structures are compatible
+
+instance (Num k) => HopfAlgebra k (Permutation Int) where
+    antipode (V ts) = nf $ V [(g^-1,x) | (g,x) <- ts]
+
+-- inject permutation into group algebra
+ip :: [[Int]] -> GroupAlgebra Q
+ip cs = return $ p cs
+
+
+instance Num k => Module k (Permutation Int) Int where
+    action = nf . fmap (\(T g x) -> x .^ g)
+
+-- use *. instead
+-- r *> m = action (r `te` m)
+
+
+
+
diff --git a/Math/Algebras/LaurentPoly.hs b/Math/Algebras/LaurentPoly.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/LaurentPoly.hs
@@ -0,0 +1,81 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances, MultiParamTypeClasses #-}
+
+
+module Math.Algebras.LaurentPoly where
+
+
+import Math.Algebra.Field.Base hiding (powers)
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+import qualified Data.List as L
+
+import Math.Algebras.Commutative -- for DivisionBasis and quotRemMP
+
+
+-- LAURENT MONOMIALS
+
+data LaurentMonomial = LM Int [(String,Int)] deriving (Eq,Ord)
+{-
+instance Ord LaurentMonomial where
+    compare (LM si xis) (LM sj yjs) = compare (-si, xis) (-sj, yjs)
+-}
+instance Show LaurentMonomial where
+    show (LM 0 []) = "1"
+    show (LM _ xis) = concatMap (\(x,i) -> if i==1 then x else x ++ "^" ++ show i) xis
+
+instance Mon LaurentMonomial where
+    munit = LM 0 []
+    mmult (LM si xis) (LM sj yjs) = LM (si+sj) $ addmerge xis yjs
+
+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 $ V [(a `mmult` b, x) | (T a b, x) <- ts]
+
+{-
+-- This is just the Set Coalgebra, so better to use a generic instance
+-- Also, not used anywhere. Hence commented out
+instance Num k => Coalgebra k LaurentMonomial where
+    counit (V ts) = sum [x | (m,x) <- ts] -- trace
+    comult = fmap (\m -> T m m)
+-}
+
+type LaurentPoly k = Vect k LaurentMonomial
+
+lvar v = V [(LM 1 [(v,1)], 1)] :: LaurentPoly Q
+
+instance Fractional k => Fractional (LaurentPoly k) where
+    recip (V [(LM si xis,c)]) = V [(LM (-si) $ map (\(x,i)->(x,-i)) xis, recip c)]
+    recip _ = error "LaurentPoly.recip: only defined for single terms"
+
+q = lvar "q"
+q' = 1/q
+
+
+{-
+-- division doesn't terminate with the derived Ord instance
+-- if we use the graded Ord instance instead, division doesn't continue into negative powers
+-- so we get the negative powers as remained, even if they're divisible
+instance DivisionBasis LaurentMonomial where
+    dividesB (LM si xis) (LM sj yjs) = si <= sj && dividesB' xis yjs where
+        dividesB' ((x,i):xis) ((y,j):yjs) =
+            case compare x y of
+            LT -> False
+            GT -> dividesB' ((x,i):xis) yjs
+            EQ -> if i<=j then dividesB' xis yjs else False
+        dividesB' [] _ = True
+        dividesB' _ [] = False
+    divB (LM si xis) (LM sj yjs) = LM (si-sj) $ divB' xis yjs where
+        divB' ((x,i):xis) ((y,j):yjs) =
+            case compare x y of
+            LT -> (x,i) : divB' xis ((y,j):yjs)
+            EQ -> if i == j then divB' xis yjs else (x,i-j) : divB' xis yjs -- we don't bother to check i > j
+            GT -> error "divB'" -- (y,-j) : divB' ((x,i):xis) yjs
+        divB' xis [] = xis
+        divB' [] yjs = error "divB'"
+-}
diff --git a/Math/Algebras/Matrix.hs b/Math/Algebras/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/Matrix.hs
@@ -0,0 +1,91 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances #-}
+
+
+module Math.Algebras.Matrix where
+
+import Math.Algebra.Field.Base
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+
+
+-- Mat2
+
+delta i j | i == j    = 1
+          | otherwise = 0
+
+data Mat2 = E2 Int Int deriving (Eq,Ord,Show)
+-- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere
+
+instance Num k => Algebra k Mat2 where
+    -- unit 0 = zero -- V []
+    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)
+
+-- In other words
+-- unit x = x (1 0)
+--            (0 1)
+-- mult (a1 b1) `te` (a2 b2) = (a1 b1) * (a2 b2) = (a b)
+--      (c1 d1)      (c2 d2)   (c1 d1)   (c2 d2)   (c d)
+
+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)
+
+-- In other words
+-- action (a b) `te` (x) = (ax+by)
+--        (c d)      (y)   (cx+dy)
+
+toMat2 [[a,b],[c,d]] = sum $ zipWith (\x e -> unit x * return e) [a,b,c,d] [E2 1 1, E2 1 2, E2 2 1, E2 2 2]
+-- fromMat2
+
+toEB2 [x,y] = foldl add zero $ zipWith (\x e -> x `smultL` return e) [x,y] [E 1, E 2]
+
+toEB xs = foldl add zero $ zipWith (\x e -> x `smultL` return e) xs (map E [1..])
+
+
+
+data Mat2' = E2' Int Int deriving (Eq,Ord,Show)
+-- E' i j represents the dual basis element corresponding to E i j
+
+-- Kassel p42
+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]])
+-- In other words
+-- counit (a b) = (1 0)
+--        (c d)   (0 1)
+-- comult (a b) = (a1 b1) `te` (a2 b2)
+--        (c d)   (c1 d1)      (c2 d2)
+-- ??
+-- ?? How does this act on Mat2?
+-- ?? What is the relationship between this and SL2 ABCD, which it seems to resemble
+
+
+
+-- !! Now do the quickchecks
+
+
+
+data M3 = E3 Int Int deriving (Eq,Ord,Show)
+-- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere
+
+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
+
+{-
+-- 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
+-- (is this order preserving?)
+-}
diff --git a/Math/Algebras/NonCommutative.hs b/Math/Algebras/NonCommutative.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/NonCommutative.hs
@@ -0,0 +1,149 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+
+-- |A module defining the algebra of non-commutative polynomials over a field k
+module Math.Algebras.NonCommutative where
+
+import Math.Algebra.Field.Base hiding (powers)
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+import qualified Data.List as L
+
+
+data NonComMonomial v = NCM Int [v] deriving (Eq)
+
+instance Ord v => Ord (NonComMonomial v) where
+    compare (NCM lx xs) (NCM ly ys) = compare (-lx, xs) (-ly, ys)
+-- ie Glex ordering
+
+instance (Eq v, Show v) => Show (NonComMonomial v) where
+    show (NCM _ []) = "1"
+    show (NCM _ vs) = concatMap showPower (L.group vs)
+        where showPower [v] = showVar v
+              showPower vs@(v:_) = showVar v ++ "^" ++ show (length vs)
+              showVar v = filter (/= '"') (show v)
+
+instance Mon (NonComMonomial v) where
+    munit = NCM 0 []
+    mmult (NCM i xs) (NCM j ys) = NCM (i+j) (xs++ys)
+
+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)
+
+{-
+-- 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)
+        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]
+-}
+
+{-
+-- This is just the Set Coalgebra, so better to use a generic instance
+-- 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)
+-}
+
+
+class Monomial m where
+    var :: v -> Vect Q (m v)
+    powers :: Eq v => m v -> [(v,Int)]
+
+V ts `bind` f = sum [c `smultL` product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] 
+-- V ts `bind` f = sum [product [f x ^ i | (x,i) <- powers m] * unit c | (m, c) <- ts] 
+
+instance Monomial NonComMonomial where
+    var v = V [(NCM 1 [v],1)]
+    powers (NCM _ vs) = map power (L.group vs)
+        where power vs@(v:_) = (v,length vs)
+
+
+type NCPoly v = Vect Q (NonComMonomial v)
+
+{-
+x,y,z :: NCPoly String
+x = var "x"
+y = var "y"
+z = var "z"
+-}
+
+
+-- DIVISION
+
+class DivisionBasis m where
+    divM :: m -> m -> Maybe (m,m)
+    -- divM a b tries to find l, r such that a = lbr
+{-
+    findOverlap :: m -> m -> Maybe (m,m,m)
+    -- given two monomials f g, find if possible a,b,c with f=ab g=bc
+-}
+
+instance Eq v => DivisionBasis (NonComMonomial v) where
+    divM (NCM _ a) (NCM _ b) = divM' [] a where
+        divM' ls (r:rs) =
+            if b `L.isPrefixOf` (r:rs)
+            then Just (ncm $ reverse ls, ncm $ drop (length b) (r:rs))
+            else divM' (r:ls) rs
+        divM' _ [] = Nothing
+{-
+    findOverlap (NCM _ xs) (NCM _ ys) = findOverlap' [] xs ys where
+        findOverlap' as [] cs = Nothing -- (reverse as, [], cs)
+        findOverlap' as (b:bs) cs =
+            if (b:bs) `L.isPrefixOf` cs
+            then Just (ncm $ reverse as, ncm $ b:bs, ncm $ drop (length (b:bs)) cs)
+            else findOverlap' (b:as) bs cs
+-}
+ncm xs = NCM (length xs) xs
+
+lm (V ((m,c):ts)) = m
+lc (V ((m,c):ts)) = c
+lt (V (t:ts)) = V [t]
+
+-- given f, gs, find ls, rs, f' such that f = sum (zipWith3 (*) ls gs rs) + f', with f' not divisible by any g
+quotRemNP f gs | all (/=0) gs = quotRemNP' f (replicate n (0,0), 0)
+               | otherwise = error "quotRemNP: division by zero"
+    where
+    n = length gs
+    quotRemNP' 0 (lrs,f') = (lrs,f')
+    quotRemNP' h (lrs,f') = divisionStep h (gs,[],lrs,f')
+    divisionStep h (g:gs, lrs', (l,r):lrs, f') =
+        case lm h `divM` lm g of
+        Just (l',r') -> let l'' = V [(l',lc h / lc g)]
+                            r'' = V [(r',1)]
+                            h' = h - l'' * g * r''
+                        in quotRemNP' h' (reverse lrs' ++ (l+l'',r+r''):lrs, f')
+        Nothing -> divisionStep h (gs,(l,r):lrs',lrs,f')
+    divisionStep h ([],lrs',[],f') =
+        let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
+        in quotRemNP' (h-lth) (reverse lrs', f'+lth)
+
+-- It is only marginally (5-10%) more space/time efficient not to track the (lazily unevaluated) factors
+remNP f gs | all (/=0) gs = remNP' f 0
+           | otherwise = error "remNP: division by zero"
+    where
+    n = length gs
+    remNP' 0 f' = f'
+    remNP' h f' = divisionStep h gs f'
+    divisionStep h (g:gs) f' =
+        case lm h `divM` lm g of
+        Just (l',r') -> let l'' = V [(l',lc h / lc g)]
+                            r'' = V [(r',1)]
+                            h' = h - l'' * g * r''
+                        in remNP' h' f'
+        Nothing -> divisionStep h gs f'
+    divisionStep h [] f' =
+        let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
+        in remNP' (h-lth) (f'+lth)
+
+infixl 7 %%
+-- f %% gs = r where (_,r) = quotRemNP f gs
+f %% gs = remNP f gs
+
diff --git a/Math/Algebras/Quaternions.hs b/Math/Algebras/Quaternions.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/Quaternions.hs
@@ -0,0 +1,58 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction #-}
+
+
+module Math.Algebras.Quaternions where
+
+import Math.Algebra.Field.Base
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+
+
+-- QUATERNIONS
+
+data HBasis = One | I | J | K deriving (Eq,Ord)
+
+type Quaternion k = Vect k HBasis
+
+instance Show HBasis where
+    show One = "1"
+    show I = "i"
+    show J = "j"
+    show K = "k"
+
+instance (Num k) => Algebra k HBasis where
+    unit 0 = zero -- V []
+    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
+
+i,j,k :: Num k => Quaternion k
+i = return I
+j = return J
+k = return K
+
+{-
+-- Set coalgebra instance
+instance Num k => Coalgebra k HBasis where
+    counit (V ts) = sum [x | (m,x) <- ts] -- trace
+    comult = fmap (\m -> T m m)           -- diagonal
+-}
+
+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)
diff --git a/Math/Algebras/Structures.hs b/Math/Algebras/Structures.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/Structures.hs
@@ -0,0 +1,160 @@
+-- Copyright (c) David Amos, 2010. All rights reserved.
+
+{-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction #-}
+{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}
+{-# LANGUAGE IncoherentInstances #-}
+
+module Math.Algebras.Structures where
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+
+
+-- MONOID
+
+-- |Monoid
+class Mon m where
+    munit :: m
+    mmult :: m -> m -> m
+
+
+-- ALGEBRAS, COALGEBRAS, BIALGEBRAS, HOPF ALGEBRAS
+
+-- |"Vect k b is a k-algebra"
+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"
+class Coalgebra k b where
+    counit :: Vect k b -> k
+    comult :: Vect k b -> Vect k (Tensor b b)
+
+
+-- |A bialgebra is an algebra which is also a coalgebra, subject to some compatibility conditions
+class (Algebra k b, Coalgebra k b) => Bialgebra k b where {}
+
+class Bialgebra k b => HopfAlgebra k b where
+    antipode :: Vect k b -> Vect k b
+
+
+instance (Num k, Eq b, Ord b, Show b, Algebra k b) => Num (Vect k b) where
+    x+y = add x y
+    negate (V ts) = V $ map (\(b,x) -> (b, negate x)) ts
+    x*y = mult (x `te` y)
+    fromInteger n = unit (fromInteger n)
+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"
+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"
+
+
+-- This is the Frobenius form, provided some conditions are met
+-- pairing = counit . mult
+
+{-
+-- A class to be used to declare that a type b should be given the set coalgebra structure
+class SetCoalgebra b where {}
+
+instance (Num k, SetCoalgebra b) => Coalgebra k b where
+    counit (V ts) = sum [x | (m,x) <- ts] -- trace
+    comult = fmap (\m -> T m m) -- diagonal
+-}
+
+
+instance Num k => Algebra k () where
+    unit 0 = zero -- V []
+    unit x = V [( (),x)]
+    mult (V [(T () (),x)]) = V [( (),x)]
+
+instance Num k => Coalgebra k () where
+    counit (V []) = 0
+    counit (V [( (),x)]) = x
+    comult (V [( (),x)]) = V [(T () (),x)]
+
+type Trivial k = Vect k ()
+
+unit' :: (Num k, Algebra k b) => Trivial k -> Vect k b
+unit' = unit . unwrap where unwrap = counit :: Num k => Trivial k -> k
+
+counit' :: (Num k, Coalgebra k b) => Vect k b -> Trivial k
+counit' = wrap . counit where wrap = unit :: Num k => k -> Trivial k
+
+-- unit' and counit' enable us to form tensors of these functions
+
+
+-- Kassel p32
+instance (Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (Tensor a b) where
+    unit 0 = V []
+    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')
+
+-- Kassel p42
+instance (Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (Tensor a b) where
+    counit = counit . (counit' `tf` counit')
+    -- counit = counit . linear (\(T x y) -> counit' (return x) * counit' (return y))
+    comult = assocL . (id `tf` assocR) . (id `tf` (twist `tf` id))
+           . (id `tf` assocL) . assocR . (comult `tf` comult)
+
+
+newtype SetCoalgebra b = SC b deriving (Eq,Ord,Show)
+
+instance Num k => Coalgebra k (SetCoalgebra b) where
+    counit (V ts) = sum [x | (m,x) <- ts] -- trace
+    comult = fmap (\m -> T m m)           -- diagonal
+
+
+newtype MonoidCoalgebra m = MC m deriving (Eq,Ord,Show)
+
+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)
+-- Brzezinski and Wisbauer, Corings and Comodules, p5
+
+-- Both of the above can be used to define coalgebra structure on polynomial algebras
+-- by using the definitions above on the generators (ie the indeterminates) and then extending multiplicatively
+-- They are then guaranteed to be algebra morphisms?
+
+
+-- MODULES AND COMODULES
+
+class Algebra k a => Module k a m where
+    action :: Vect k (Tensor a m) -> Vect k m
+
+r *. m = action (r `te` m)
+
+class Coalgebra k c => Comodule k c n where
+    coaction :: Vect k n -> Vect k (Tensor c n)
+
+
+instance Algebra k a => Module k a a where
+    action = mult
+
+instance Coalgebra k c => Comodule k c c where
+    coaction = comult
+
+-- module and comodule instances for tensor products
+
+-- Kassel p57-8
+
+instance (Num k, Ord a, Ord u, Ord v, Algebra k a, Module k a u, Module k a v)
+         => 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)
+
+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)
+-- !! 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
+
+-- Kassel p63
+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
diff --git a/Math/Algebras/TensorAlgebra.hs b/Math/Algebras/TensorAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/TensorAlgebra.hs
@@ -0,0 +1,54 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+
+-- |A module defining the tensor algebra, symmetric algebra, and exterior (or alternating) algebra
+module Math.Algebras.TensorAlgebra where
+
+import qualified Data.List as L
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+
+import Math.Algebra.Field.Base
+
+
+data TensorAlgebra a = TA Int [a] deriving (Eq,Ord,Show)
+
+instance Mon (TensorAlgebra a) where
+    munit = TA 0 []
+    mmult (TA i xs) (TA j ys) = TA (i+j) (xs++ys)
+
+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)
+
+
+data SymmetricAlgebra a = Sym Int [a] deriving (Eq,Ord,Show)
+
+instance Ord a => Mon (SymmetricAlgebra a) where
+    munit = Sym 0 []
+    mmult (Sym i xs) (Sym j ys) = Sym (i+j) $ L.sort (xs++ys)
+
+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)
+
+
+data ExteriorAlgebra a = Ext Int [a] deriving (Eq,Ord,Show)
+
+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))
+        where signedMerge s (k,zs) (i,x:xs) (j,y:ys) =
+                  case compare x y of
+                  EQ -> zero
+                  LT -> signedMerge s (k+1,x:zs) (i-1,xs) (j,y:ys)
+                  GT -> let s' = if even i then s else -s -- we had to commute y past x:xs, with i sign changes
+                        in signedMerge s' (k+1,y:zs) (i,x:xs) (j-1,ys)
+              signedMerge s (k,zs) (i,xs) (0,[]) = s *> (return $ Ext (k+i) $ reverse zs ++ xs)
+              signedMerge s (k,zs) (0,[]) (j,ys) = s *> (return $ Ext (k+j) $ reverse zs ++ ys)
diff --git a/Math/Algebras/TensorProduct.hs b/Math/Algebras/TensorProduct.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/TensorProduct.hs
@@ -0,0 +1,47 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE NoMonomorphismRestriction #-}
+
+-- |A module defining tensor products of vector spaces
+module Math.Algebras.TensorProduct where
+
+import Math.Algebras.VectorSpace
+
+
+data Tensor a b = T a b deriving (Eq, Ord, Show)
+-- or T !a !b, forcing strictness, but not proven to be better
+
+
+-- |Tensor product of two 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]
+-- 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
+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]
+    where sum = foldl add zero -- (V [])
+
+
+-- tensor isomorphisms
+
+-- 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)
+
+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))
+
+inUnitL = fmap (\a -> T () a)
+
+inUnitR = fmap (\a -> T a ())
+
+outUnitL = fmap (\(T () a) -> a)
+
+outUnitR = fmap (\(T a ()) -> a)
+
+twist v = nf $ fmap (\(T a b) -> T b a) v
+-- note the nf call, as f is not order-preserving
+
diff --git a/Math/Algebras/VectorSpace.hs b/Math/Algebras/VectorSpace.hs
new file mode 100644
--- /dev/null
+++ b/Math/Algebras/VectorSpace.hs
@@ -0,0 +1,130 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE NoMonomorphismRestriction #-}
+
+-- |A module defining the type and operations of free k-vector spaces over a basis b (for a field k)
+module Math.Algebras.VectorSpace where
+
+import qualified Data.List as L
+import qualified Data.Set as S -- only needed for toSet
+
+-- toSet = S.toList . S.fromList
+
+infixr 7 *>
+infixl 7 <*
+infixl 6 <+>
+
+
+-- |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)
+
+instance (Num k, Show b) => Show (Vect k b) where
+    show (V []) = "0"
+    show (V ts) = concatWithPlus $ map showTerm ts
+        where showTerm (b,x) | show b == "1" = show x
+                             | show x == "1" = show b
+                             | show x == "-1" = "-" ++ show b
+                             | otherwise = (if isAtomic (show x) then show x else "(" ++ show x ++ ")")
+                                           -- (let (c:cs) = show x in
+                                           -- if any (`elem` "+-") cs then "(" ++ show x ++ ")" else show x)
+                                           ++ show b
+              concatWithPlus (t1:t2:ts) = if head t2 == '-'
+                                          then t1 ++ concatWithPlus (t2:ts)
+                                          else t1 ++ '+' : concatWithPlus (t2:ts)
+              concatWithPlus [t] = t
+              isAtomic (c:cs) = isAtomic' cs
+              isAtomic' ('^':'-':cs) = isAtomic' cs
+              isAtomic' ('+':cs) = False
+              isAtomic' ('-':cs) = False
+              isAtomic' (c:cs) = isAtomic' cs
+              isAtomic' [] = True
+
+-- |The zero vector
+zero :: Vect k b
+zero = V []
+
+-- |Addition of vectors
+add :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
+add (V ts) (V us) = V $ addmerge ts us
+
+-- |Addition of vectors (same as add)
+(<+>) :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
+(<+>) = add
+
+addmerge ((a,x):ts) ((b,y):us) =
+    case compare a b of
+    LT -> (a,x) : addmerge ts ((b,y):us)
+    EQ -> if x+y == 0 then addmerge ts us else (a,x+y) : addmerge ts us
+    GT -> (b,y) : addmerge ((a,x):ts) us
+addmerge ts [] = ts
+addmerge [] us = us
+
+-- |Negation of vector
+neg :: (Num k) => Vect k b -> Vect k b
+neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts
+
+-- |Scalar multiplication (on the left)
+smultL :: (Num k) => k -> Vect k b -> Vect k b
+smultL 0 _ = zero -- V []
+smultL k (V ts) = V [(ei,k*xi) | (ei,xi) <- ts]
+
+-- |Same as smultL. Mnemonic is "multiply through (from the left)"
+(*>) :: (Num k) => k -> Vect k b -> Vect k b
+(*>) = smultL
+
+-- |Scalar multiplication on the right
+smultR :: (Num k) => Vect k b -> k -> Vect k b
+smultR _ 0 = zero -- V []
+smultR (V ts) k = V [(ei,xi*k) | (ei,xi) <- ts]
+
+-- |Same as smultR. Mnemonic is "multiply through (from the right)"
+(<*) :: (Num k) => Vect k b -> k -> Vect k b
+(<*) = smultR
+
+-- same as return
+-- injection of basis elt into vector space
+-- inject b = V [(b,1)]
+
+-- same as fmap
+-- liftFromBasis f (V ts) = V [(f b, x) | (b, x) <- ts]
+-- if f is not order-preserving, then you need to call nf afterwards
+
+-- |Convert an element of Vect k b into normal form. Normal form consists in having the basis elements in ascending order,
+-- with no duplicates, and all coefficients non-zero
+nf :: (Ord b, Num k) => Vect k b -> Vect k b
+nf (V ts) = V $ nf' $ L.sortBy compareFst ts where
+    nf' ((b1,x1):(b2,x2):ts) =
+        case compare b1 b2 of
+        LT -> if x1 == 0 then nf' ((b2,x2):ts) else (b1,x1) : nf' ((b2,x2):ts)
+        EQ -> if x1+x2 == 0 then nf' ts else nf' ((b1,x1+x2):ts)
+        GT -> error "nf': not pre-sorted"
+    nf' [(b,x)] = if x == 0 then [] else [(b,x)]
+    nf' [] = []
+    compareFst (b1,x1) (b2,x2) = compare b1 b2
+    -- compareFst = curry ( uncurry compare . (fst *** fst) )
+
+
+-- lift a function on the basis to a function on the vector space
+instance Functor (Vect k) where
+    fmap f (V ts) = V [(f b, x) | (b,x) <- ts]
+-- Note that if f is not order-preserving, then we need to call "nf" afterwards
+
+instance Num k => Monad (Vect k) where
+    return a = V [(a,1)]
+    V ts >>= f = V $ concat [ [(b,y*x) | let V us = f a, (b,y) <- us] | (a,x) <- ts]
+    -- Note that as we can't assume Ord a in the Monad instance, we need to call "nf" afterwards
+
+linear :: (Ord b, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b
+linear f v = nf $ v >>= f
+
+newtype EBasis = E Int deriving (Eq,Ord)
+
+instance Show EBasis where show (E i) = "e" ++ show i
+
+e i = return (E i)
+e1 = e 1
+e2 = e 2
+e3 = e 3
+
+-- dual (E i) = E (-i)
diff --git a/Math/Combinatorics/Design.hs b/Math/Combinatorics/Design.hs
--- a/Math/Combinatorics/Design.hs
+++ b/Math/Combinatorics/Design.hs
@@ -1,4 +1,4 @@
--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
 module Math.Combinatorics.Design where
 
@@ -12,7 +12,7 @@
 import Math.Algebra.Field.Extension
 import Math.Algebra.Group.PermutationGroup hiding (elts, order, isMember)
 import Math.Algebra.Group.SchreierSims as SS
-import Math.Combinatorics.Graph as G hiding (to1n)
+import Math.Combinatorics.Graph as G hiding (to1n, incidenceMatrix)
 import Math.Combinatorics.GraphAuts (refine', isSingleton, graphAuts, incidenceAuts) -- , removeGens)
 import Math.Combinatorics.FiniteGeometry
 
@@ -89,9 +89,10 @@
 isSquare d@(D xs bs) = is2Design d && length xs == length bs
 
 
--- incidence matrix of a design
--- (rows and columns indexed by blocks and points respectively)
--- (this follows Cameron & van Lint, though elsewhere in the literature it is sometimes the other way round)
+-- (We follow Cameron & van Lint.)
+-- |The incidence matrix of a design, with rows indexed by blocks and columns by points.
+-- (Note that in the literature, the opposite convention is sometimes used instead.)
+incidenceMatrix :: (Eq t) => Design t -> [[Int]]
 incidenceMatrix (D xs bs) = [ [if x `elem` b then 1 else 0 | x <- xs] | b <- bs]
 
 
@@ -109,13 +110,15 @@
     vs = vertices graph
     es = edges graph
 
--- the affine plane AG(2,Fq) - a 2-(q^2,q,1) design
+-- |The affine plane AG(2,Fq), a 2-(q^2,q,1) design
+ag2 :: (FiniteField k, Ord k) => [k] -> Design [k]
 ag2 fq = design (points, lines) where
     points = ptsAG 2 fq
     lines = map line $ tail $ ptsPG 2 fq
     line [a,b,c] = [ [x,y] | [x,y] <- points, a*x+b*y+c==0 ]
 
--- the projective plane PG(2,Fq) - a square 2-(q^2+q+1,q+1,1) design
+-- |The projective plane PG(2,Fq), a square 2-(q^2+q+1,q+1,1) design
+pg2 :: (FiniteField k, Ord k) => [k] -> Design [k]
 pg2 fq = design (points, lines) where
     points = ptsPG 2 fq
     lines = L.sort $ map line points
@@ -168,18 +171,22 @@
 -- NEW DESIGNS FROM OLD
 
 -- Dual of a design. Cameron & van Lint p11
+-- |The dual of a design
+dual :: (Ord t) => Design t -> Design [t]
 dual (D xs bs) = design (bs, map beta xs) where
     beta x = filter (x `elem`) bs
 
 -- Derived design relative to a point. Cameron & van Lint p11
 -- Derived design of a t-(v,k,lambda) is a t-1-(v-1,k-1,lambda) design.
+derivedDesign :: (Ord t) => Design t -> t -> Design t
 derivedDesign (D xs bs) p = design (xs L.\\ [p], [b L.\\ [p] | b <- bs, p `elem` b])
 
 -- Residual design relative to a point. Cameron & van Lint p13
 -- Point-residual of a t-(v,k,lambda) is a t-1-(v-1,k,mu).
+pointResidual :: (Ord t) => Design t -> t -> Design t
 pointResidual (D xs bs) p = design (xs L.\\ [p], [b | b <- bs, p `notElem` b])
 
--- Complementary design relative to a point. Cameron & van Lint p13
+-- Complementary design. Cameron & van Lint p13
 -- Complement of a t-(v,k,lambda) is a t-(v,v-k,mu).
 complementaryDesign (D xs bs) = design (xs, [xs L.\\ b | b <- bs])
 
@@ -187,6 +194,7 @@
 -- This is only a design if (xs,bs) is a square design
 -- It may have repeated blocks - but if so, residuals of the complement will not
 -- Block-residual of a 2-(v,k,lambda) is a 2-(v-k,k-lambda,lambda).
+blockResidual :: (Ord t) => Design t -> [t] -> Design t
 blockResidual d@(D xs bs) b | isSquare d = design (xs L.\\ b, [b' L.\\ b | b' <- bs, b' /= b])
 
 
@@ -195,11 +203,15 @@
 isDesignAut (D xs bs) g | supp g `isSubset` xs = all (`S.member` bs') [b -^ g | b <- bs]
     where bs' = S.fromList bs
 
+-- |The incidence graph of a design
+incidenceGraph :: (Ord a) => Design a -> Graph (Either a [a])
 incidenceGraph (D xs bs) = G vs es where -- graph (vs,es) where
     vs = L.sort $ map Left xs ++ map Right bs
     es = L.sort [ [Left x, Right b] | x <- xs, b <- bs, x `elem` b ]
 
 
+-- |Find a strong generating set for the automorphism group of a design
+designAuts :: (Ord t) => Design t -> [Permutation t]
 designAuts d = incidenceAuts $ incidenceGraph d
 
 -- We find design auts by finding graph auts of the incidence graph of the design
@@ -229,16 +241,21 @@
 deltaM24 = p [[-1],[0],[1,18,4,2,6],[3],[5,21,20,10,7],[8,16,13,9,12],[11,19,22,14,17],[15]]
 -- this is t -> t^3 / 9 (for t a quadratic residue), t -> 9 t^3 (t a non-residue)
 
+-- |Generators for the Mathieu group M24, a finite simple group of order 244823040
+m24 :: [Permutation Integer]
 m24 = [alphaL2_23, betaL2_23, gammaL2_23, deltaM24]
 
+-- |A strong generating set for the Mathieu group M24, a finite simple group of order 244823040
+m24sgs :: [Permutation Integer]
 m24sgs = sgs m24
--- order 244823040
 
+-- |A strong generating set for the Mathieu group M23, a finite simple group of order 10200960
+m23sgs :: [Permutation Integer]
 m23sgs = filter (\g -> (-1).^g == -1) m24sgs
--- order 10200960
 
+-- |A strong generating set for the Mathieu group M22, a finite simple group of order 443520
+m22sgs :: [Permutation Integer]
 m22sgs = filter (\g -> 0.^g == 0) m23sgs
--- order 443520
 
 -- sgs uses the base implied by the Ord instance, which will be [-1,0,..]
 
@@ -248,16 +265,21 @@
 octad = [0,1,2,3,4,7,10,12]
 -- Conway&Sloane p276 - this is a weight 8 codeword from Golay code G24
 
-
+-- |The Steiner system S(5,8,24), with 759 blocks, whose automorphism group is M24
+s_5_8_24 :: Design Integer
 s_5_8_24 = design ([-1..22], octad -^^ l2_23)
 -- S(5,8,24) constructed as the image of a single octad under the action of PSL(2,23)
 -- 759 blocks ( (24 `choose` 5) `div` (8 `choose` 5) )
 -- Automorphism group is M24
 
+-- |The Steiner system S(4,7,23), with 253 blocks, whose automorphism group is M23
+s_4_7_23 :: Design Integer
 s_4_7_23 = derivedDesign s_5_8_24 (-1)
 -- 253 blocks ( (23 `choose` 4) `div` (7 `choose` 4) )
 -- Automorphism group is M23
 
+-- |The Steiner system S(3,6,22), with 77 blocks, whose automorphism group is M22
+s_3_6_22 :: Design Integer
 s_3_6_22 = derivedDesign s_4_7_23 0
 -- 77 blocks
 -- Automorphism group is M22
@@ -293,20 +315,30 @@
 -- the squares (quadratic residues) in F11
 -- http://en.wikipedia.org/wiki/Steiner_system
 
+-- |The Steiner system S(5,6,12), with 132 blocks, whose automorphism group is M12
+s_5_6_12 :: Design Integer
 s_5_6_12 = design ([-1..10], hexad -^^ l2_11)
 -- S(5,6,12) constructed as the image of a single hexad under the action of PSL(2,11)
 -- 132 blocks ( (12 `choose` 5) `div` (6 `choose` 5) )
 -- Automorphism group is M12
 
+-- |The Steiner system S(4,5,11), with 66 blocks, whose automorphism group is M11
+s_4_5_11 :: Design Integer
 s_4_5_11 = derivedDesign s_5_6_12 (-1)
 -- 66 blocks
 -- Automorphism group is M11
 
+-- |Generators for the Mathieu group M12, a finite simple group of order 95040
+m12 :: [Permutation Integer]
 m12 = [alphaL2_11, betaL2_11, gammaL2_11, deltaM12]
 
+-- |A strong generating set for the Mathieu group M12, a finite simple group of order 95040
+m12sgs :: [Permutation Integer]
 m12sgs = sgs m12
 -- order 95040
 
+-- |A strong generating set for the Mathieu group M11, a finite simple group of order 7920
+m11sgs :: [Permutation Integer]
 m11sgs = filter (\g -> (-1).^g == -1) m12sgs
 -- order 7920
 
diff --git a/Math/Combinatorics/Hypergraph.hs b/Math/Combinatorics/Hypergraph.hs
--- a/Math/Combinatorics/Hypergraph.hs
+++ b/Math/Combinatorics/Hypergraph.hs
@@ -1,5 +1,6 @@
--- Copyright (c) David Amos, 2009. All rights reserved.
+-- Copyright (c) 2009, David Amos. All rights reserved.
 
+-- |A module defining a type for hypergraphs.
 module Math.Combinatorics.Hypergraph where
 
 import qualified Data.List as L
@@ -24,7 +25,8 @@
 -- this still doesn't guarantee that all bs are subset of xs
 
 
--- uniform hypergraph - all blocks are same size
+-- |Is this hypergraph uniform - meaning that all blocks are of the same size
+isUniform :: (Ord a) => Hypergraph a -> Bool
 isUniform h@(H xs bs) = isSetSystem xs bs && same (map length bs)
 
 same (x:xs) = all (==x) xs
@@ -46,14 +48,15 @@
 
 -- INCIDENCE GRAPH
 
-data Incidence a = P a | B [a] deriving (Eq, Ord, Show)
+-- data Incidence a = P a | B [a] deriving (Eq, Ord, Show)
 
 -- compare Design, where we just use Left, Right
 
 -- Also called the Levi graph
+incidenceGraph :: (Ord a) => Hypergraph a -> Graph (Either a [a])
 incidenceGraph (H xs bs) = G vs es where
-    vs = map P xs ++ map B bs
-    es = L.sort [ [P x, B b] | b <- bs, x <- b]
+    vs = map Left xs ++ map Right bs
+    es = L.sort [ [Left x, Right b] | b <- bs, x <- b]
 
 
 -- INCIDENCE MATRIX
@@ -81,19 +84,24 @@
 -- We can represent various incidence structures as hypergraphs,
 -- by identifying the lines with the sets of points that they contain
 
+isPartialLinearSpace :: (Ord a) => Hypergraph a -> Bool
 isPartialLinearSpace h@(H ps ls) =
     isSetSystem ps ls &&
     all ( (<=1) . length ) [filter (pair `isSubset`) ls | pair <- combinationsOf 2 ps]
     -- any two points are incident with at most one line
 
 -- Godsil & Royle, p79
+-- |Is this hypergraph a projective plane - meaning that any two lines meet in a unique point,
+-- and any two points lie on a unique line
+isProjectivePlane :: (Ord a) => Hypergraph a -> Bool
 isProjectivePlane h@(H ps ls) =
     isSetSystem ps ls &&
     all ( (==1) . length) [intersect l1 l2 | [l1,l2] <- combinationsOf 2 ls] && -- any two lines meet in a unique point
     all ( (==1) . length) [ filter ([p1,p2] `isSubset`) ls | [p1,p2] <- combinationsOf 2 ps] -- any two points lie in a unique line
 
--- a projective plane with a triangle
--- this is a weak non-degeneracy condition, which eliminates all points on the same line, or all lines through the same point
+-- |Is this hypergraph a projective plane with a triangle.
+-- This is a weak non-degeneracy condition, which eliminates all points on the same line, or all lines through the same point.
+isProjectivePlaneTri :: (Ord a) => Hypergraph a -> Bool
 isProjectivePlaneTri h@(H ps ls) =
     isProjectivePlane h && any triangle (combinationsOf 3 ps)
     where triangle t@[p1,p2,p3] =
@@ -101,8 +109,9 @@
                    (not . null) [l | l <- ls, [p1,p3] `isSubset` l, p2 `notElem` l] &&
                    (not . null) [l | l <- ls, [p2,p3] `isSubset` l, p1 `notElem` l] 
 
--- a projective plane with a quadrangle
--- this is a stronger non-degeneracy condition
+-- |Is this hypergraph a projective plane with a quadrangle.
+-- This is a stronger non-degeneracy condition.
+isProjectivePlaneQuad :: (Ord a) => Hypergraph a -> Bool
 isProjectivePlaneQuad h@(H ps ls) =
     isProjectivePlane h && any quadrangle (combinationsOf 4 ps)
     where quadrangle q = all (not . collinear) (combinationsOf 3 q) -- no three points collinear
@@ -116,6 +125,7 @@
 -- GENERALIZED QUADRANGLES
 
 -- Godsil & Royle p81
+isGeneralizedQuadrangle :: (Ord a) => Hypergraph a -> Bool
 isGeneralizedQuadrangle h@(H ps ls) =
     isPartialLinearSpace h &&
     all (\(l,p) -> unique [p' | p' <- l, collinear (pair p p')]) [(l,p) | l <- ls, p <- ps, p `notElem` l] &&
@@ -143,29 +153,36 @@
 -- CONFIGURATIONS
 
 -- http://en.wikipedia.org/wiki/Projective_configuration
+-- |Is this hypergraph a (projective) configuration.
+isConfiguration :: (Ord a) => Hypergraph a -> Bool
 isConfiguration h@(H ps ls) =
     isUniform h && -- a set system, with each line incident with the same number of points
     same [length (filter (p `elem`) ls) | p <- ps] -- each point is incident with the same number of lines
 
 
-
-
+fanoPlane :: Hypergraph Integer
 fanoPlane = toHypergraph [1..7] [[1,2,4],[2,3,5],[3,4,6],[4,5,7],[5,6,1],[6,7,2],[7,1,3]]
 
+-- |The Heawood graph is the incidence graph of the Fano plane
+heawoodGraph :: Graph (Either Integer [Integer])
 heawoodGraph = incidenceGraph fanoPlane
 
 
+desarguesConfiguration :: Hypergraph [Integer]
 desarguesConfiguration = H xs bs where
     xs = combinationsOf 2 [1..5]
     bs = [ [x | x <- xs, x `isSubset` b] | b <- combinationsOf 3 [1..5] ]
 
+desarguesGraph :: Graph (Either [Integer] [[Integer]])
 desarguesGraph = incidenceGraph desarguesConfiguration
 
 
+pappusConfiguration :: Hypergraph Integer
 pappusConfiguration = H xs bs where
     xs = [1..9]
     bs = L.sort [ [1,2,3], [4,5,6], [7,8,9], [1,5,9], [1,6,8], [2,4,9], [3,4,8], [2,6,7], [3,5,7] ]
 
+pappusGraph :: Graph (Either Integer [Integer])
 pappusGraph = incidenceGraph pappusConfiguration
 
 
@@ -187,6 +204,7 @@
     -- (so these are the projective planes over 7 points)
 -}
 -- Godsil & Royle p69
+coxeterGraph :: Graph [Integer]
 coxeterGraph = G vs es where
     g = p [[1..7]]
     vs = L.sort $ concatMap (orbitB [g]) [[1,2,4],[3,5,7],[3,6,7],[5,6,7]]
@@ -204,7 +222,8 @@
                            d2 <- duads, d2 > d1, disjoint d1 d2,
                            d3 <- duads, d3 > d2, disjoint d1 d3, disjoint d2 d3 ]
 
--- Tutte 8-cage
+-- |The Tutte-Coxeter graph, also called the Tutte 8-cage
+tutteCoxeterGraph :: Graph (Either [Integer] [[Integer]])
 tutteCoxeterGraph = incidenceGraph $ H duads synthemes
 
 
diff --git a/Math/Combinatorics/LatinSquares.hs b/Math/Combinatorics/LatinSquares.hs
--- a/Math/Combinatorics/LatinSquares.hs
+++ b/Math/Combinatorics/LatinSquares.hs
@@ -18,6 +18,7 @@
 
 -- LATIN SQUARES
 
+findLatinSqs :: (Eq a) => [a] -> [[[a]]]
 findLatinSqs xs = findLatinSqs' 1 [xs] where
     n = length xs
     findLatinSqs' i rows
@@ -28,6 +29,7 @@
                                     | r <- rs, r `notElem` col]
     findRows [] ls _ = [reverse ls]
 
+isLatinSq :: (Ord a) => [[a]] -> Bool
 isLatinSq rows = all isOneOfEach rows && all isOneOfEach cols where
     cols = L.transpose rows
 
@@ -54,6 +56,8 @@
 
 -- ORTHOGONAL AND MUTUALLY ORTHOGONAL LATINS SQUARES
 
+-- |Are the two latin squares orthogonal?
+isOrthogonal :: (Ord a, Ord b) => [[a]] -> [[b]] -> Bool
 isOrthogonal greeks latins = isOneOfEach pairs
     where pairs = zip (concat greeks) (concat latins)
 
@@ -65,10 +69,13 @@
         else findMOLS' i ls rs
     findMOLS' _ _ [] = []
 
+-- |Are the latin squares mutually orthogonal (ie each pair is orthogonal)?
+isMOLS :: (Ord a) => [[[a]]] -> Bool
 isMOLS (greek:latins) = all (isOrthogonal greek) latins && isMOLS latins
 isMOLS [] = True
 
--- MOLS from a projective plane
+-- |MOLS from a projective plane
+fromProjectivePlane :: (Ord k, Num k) => Design [k] -> [[[Int]]]
 fromProjectivePlane (D xs bs) = map toLS parallelClasses where
     k = [x | [0,1,x] <- xs] -- the field we're working over
     n = length k            -- the order of the projective plane
diff --git a/Math/Combinatorics/StronglyRegularGraph.hs b/Math/Combinatorics/StronglyRegularGraph.hs
--- a/Math/Combinatorics/StronglyRegularGraph.hs
+++ b/Math/Combinatorics/StronglyRegularGraph.hs
@@ -1,5 +1,6 @@
--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
+-- |A module defining various strongly regular graphs, including the Clebsch, Hoffman-Singleton, Higman-Sims, and McLaughlin graphs
 module Math.Combinatorics.StronglyRegularGraph where
 
 import qualified Data.List as L
diff --git a/Math/Projects/ChevalleyGroup/Classical.hs b/Math/Projects/ChevalleyGroup/Classical.hs
--- a/Math/Projects/ChevalleyGroup/Classical.hs
+++ b/Math/Projects/ChevalleyGroup/Classical.hs
@@ -1,4 +1,4 @@
--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Classical where
 
@@ -20,15 +20,18 @@
 
 -- LINEAR GROUPS
 
--- SL(n,Fq) is generated by elementary transvections
--- sl :: FiniteField k => Int -> k -> [[[k]]]
+-- |The special linear group SL(n,Fq), generated by elementary transvections, returned as matrices
+sl :: FiniteField k => Int -> [k] -> [[[k]]]
 sl n fq = [elemTransvection n (r,c) l | r <- [1..n], c <- [1..n], r /= c, l <- fq']
     where fq' = basisFq undefined -- tail fq
     -- Carter p68 - x_r(t1) x_r(t2) == x_r(t1+t2) - this is true in general, not just in this case
 
 elemTransvection n (r,c) l = fMatrix n (\i j -> if i == j then 1 else if (i,j) == (r,c) then l else 0)
 
--- PSL(n,Fq) == A(n,Fq) == SL(n,Fq)/Z
+-- |The projective special linear group PSL(n,Fq) == A(n,Fq) == SL(n,Fq)/Z,
+-- returned as permutations of the points of PG(n-1,Fq).
+-- This is a finite simple group provided n>2 or q>3.
+l :: (FiniteField k, Ord k) => Int -> [k] -> [Permutation [k]]
 l n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- sl n fq]
     where ps = ptsPG (n-1) fq
 
@@ -39,6 +42,8 @@
 -- SYMPLECTIC GROUPS
 -- Carter p186 and 181-3
 
+-- |The symplectic group Sp(2n,Fq), returned as matrices
+sp2 :: FiniteField k => Int -> [k] -> [[[k]]]
 sp2 n fq =
     [_I <<+>> t *>> (e i j <<->> e (-j) (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
     [_I <<->> t *>> (e (-i) (-j) <<->> e j i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
@@ -52,7 +57,10 @@
         e i j = e' (if i > 0 then i else n-i) (if j > 0 then j else n-j)
         e' i j = fMatrix (2*n) (\k l -> if (k,l) == (i,j) then 1 else 0)
 
--- PSp2n(Fq) == Cn(Fq) == Sp2n(Fq)/Z
+-- |The projective symplectic group PSp(2n,Fq) == Cn(Fq) == Sp(2n,Fq)/Z,
+-- returned as permutations of the points of PG(2n-1,Fq).
+-- This is a finite simple group for n>1, except for PSp(4,F2).
+s2 :: (FiniteField k, Ord k) => Int -> [k] -> [Permutation [k]]
 s2 n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- sp2 n fq]
     where ps = ptsPG (2*n-1) fq
 
diff --git a/Math/Projects/ChevalleyGroup/Exceptional.hs b/Math/Projects/ChevalleyGroup/Exceptional.hs
--- a/Math/Projects/ChevalleyGroup/Exceptional.hs
+++ b/Math/Projects/ChevalleyGroup/Exceptional.hs
@@ -1,4 +1,4 @@
--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Exceptional where
 
@@ -9,7 +9,8 @@
 import Math.Algebra.LinearAlgebra
 
 import Math.Algebra.Group.PermutationGroup hiding (fromList)
-import Math.Algebra.Group.SchreierSims as SS
+-- import Math.Algebra.Group.SchreierSims as SS
+import Math.Algebra.Group.RandomSchreierSims as RSS
 
 import Math.Combinatorics.FiniteGeometry (ptsAG)
 -- import ClassicalChevalleyGroup (ptsAG)
@@ -24,13 +25,14 @@
 
 newtype Octonion k = O [(Int,k)] deriving (Eq, Ord)
 
-i0 = O [(0,1)] :: Octonion Q
-i1 = O [(1,1)] :: Octonion Q
-i2 = O [(2,1)] :: Octonion Q
-i3 = O [(3,1)] :: Octonion Q
-i4 = O [(4,1)] :: Octonion Q
-i5 = O [(5,1)] :: Octonion Q
-i6 = O [(6,1)] :: Octonion Q
+i0, i1, i2, i3, i4, i5, i6 :: Octonion Q
+i0 = O [(0,1)]
+i1 = O [(1,1)]
+i2 = O [(2,1)]
+i3 = O [(3,1)]
+i4 = O [(4,1)]
+i5 = O [(5,1)]
+i6 = O [(6,1)]
 
 fromList as = O $ filter ((/=0) . snd) $ zip [-1..6] as
 
@@ -155,8 +157,11 @@
 beta3' = fromPairs [(x, x %^ beta3) | x <- unitImagOctonions f3]
 gamma3' = fromPairs [(x, x %^ gamma3) | x <- unitImagOctonions f3]
 
+-- |Generators for G2(3), a finite simple group of order 4245696,
+-- as a permutation group on the 702 unit imaginary octonions over F3
+g2_3 :: [Permutation (Octonion F3)]
+g2_3 = [alpha3', beta3', gamma3']
 -- These three together generate a group of order 4245696, which is therefore the whole of G2(3)
--- (But takes nearly 10 minutes to construct the BSGS in interpreter)
 
 -- Unit imaginary octonions form one orbit under the action of G2
 
diff --git a/Math/QuantumAlgebra/OrientedTangle.hs b/Math/QuantumAlgebra/OrientedTangle.hs
new file mode 100644
--- /dev/null
+++ b/Math/QuantumAlgebra/OrientedTangle.hs
@@ -0,0 +1,210 @@
+-- Copyright (c) David Amos, 2010. All rights reserved.
+
+{-# LANGUAGE TypeFamilies, EmptyDataDecls #-}
+
+
+module Math.QuantumAlgebra.OrientedTangle where
+
+import Math.Algebra.Field.Base
+import Math.Algebras.LaurentPoly -- hiding (lvar, q, q')
+
+import Math.QuantumAlgebra.TensorCategory
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+
+-- import MathExperiments.Algebra.TAlgebra
+
+
+-- ORIENTED TANGLE CATEGORY
+
+data Oriented = Plus | Minus deriving (Eq,Ord,Show)
+
+data HorizDir = ToL | ToR deriving (Eq,Ord,Show)
+
+data OrientedTangle
+
+-- In GHCi 6.12.1, we appear to be limited to 8 value constructors within an associated data family
+instance Category OrientedTangle where
+    data Ob OrientedTangle = OT [Oriented] deriving (Eq,Ord,Show)
+    data Ar OrientedTangle = IdT [Oriented]
+                           | CapT HorizDir
+                           | CupT HorizDir
+                           | XPlus | XMinus
+                           | SeqT [Ar OrientedTangle]
+                           | ParT [Ar OrientedTangle]
+                           deriving (Eq,Ord,Show)
+    id_ (OT os) = IdT os
+    source (IdT os) = OT os
+    source (CapT _) = OT []
+    source (CupT toR) = OT [Plus,Minus]
+    source (CupT toL) = OT [Minus,Plus]
+    source XPlus = OT [Plus,Plus]
+    source XMinus = OT [Plus,Plus]
+    source (ParT as) = OT $ concatMap ((\(OT os) -> os) . source) as
+    source (SeqT as) = source (head as)
+    target (IdT os) = OT os
+    target (CapT toR) = OT [Minus,Plus]
+    target (CapT toL) = OT [Plus,Minus]
+    target (CupT _) = OT []
+    target XPlus = OT [Plus,Plus]
+    target XMinus = OT [Plus,Plus]
+    target (ParT as) = OT $ concatMap ((\(OT os) -> os) . target) as
+    target (SeqT as) = target (last as)
+    a >>> b | target a == source b = SeqT [a,b]
+
+instance TensorCategory OrientedTangle where
+    tunit = OT []
+    tob (OT as) (OT bs) = OT (as++bs)
+    tar a b = ParT [a,b]
+
+
+
+
+
+
+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
+
+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)))
+
+-- 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))
+
+testcc' m v = nf $ v >>= c m >>= c' m
+
+mu m (E i) = (1 / (lambda m * q ^ (2*i-1))) *> return (E i)
+
+mu' m (E i) = (lambda m * q ^ (2*i-1)) *> return (E i)
+
+-- The following are modified from Kassel. We compose diagrams downwards, whereas he composes them upwards.
+
+capRL m = coevalV m
+
+capLR m = do
+    T i j <- coevalV' m
+    k <- mu' m j
+    return (T i k)
+
+cupRL m = evalV
+
+cupLR m (T i j) = do
+    k <- mu m i
+    evalV' (T k j)    
+-- linear evalV' . (linear (mu' m) `tf` idV)
+
+
+
+xplus m = c m
+
+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)
+
+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)
+
+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)
+
+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)
+
+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)
+
+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)
+
+{-
+Then we have for example the following:
+> let v = e1 `te` e2 in nf $ v >>= xplus 2 >>= xminus 2
+T e1 e2
+> let v = e (-1) `te` e2 in nf $ v >>= yplus 2 >>= tminus 2
+T e-1 e2
+> let v = e (-1) `te` e (-2) in nf $ v >>= zplus 2 >>= zminus 2
+T e-1 e-2
+
+-}
+
+
+oloop m = nf $ do
+    T a b <- capLR m
+    cupRL m (T 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)
+
+-- 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)
+
+
+{-
+-- REPRESENTATIONS OF THE TANGLE CATEGORY IN VECTOR SPACE CATEGORY
+-- But we need to convert the above code to use TensorAlgebra first
+
+kauffman :: Ar Tangle -> TangleRep [Oriented] -> TangleRep [Oriented]
+kauffman (IdT n) = id -- could be tf of n ids
+kauffman CapT = linear cap
+kauffman CupT = linear cup
+kauffman OverT = linear over
+kauffman UnderT = linear under
+kauffman (SeqT fs) = foldl (>>>) id $ map kauffman fs
+    where g >>> h = h . g
+kauffman (ParT [f]) = kauffman f
+kauffman (ParT (f:fs)) = tf m (kauffman f) (kauffman (ParT fs))
+    where OT m = source f
+          tf m f' fs' = linear (\xs -> let (ls,rs) = splitAt m xs in f' (return ls) * fs' (return rs) )
+-}
diff --git a/Math/QuantumAlgebra/QuantumPlane.hs b/Math/QuantumAlgebra/QuantumPlane.hs
new file mode 100644
--- /dev/null
+++ b/Math/QuantumAlgebra/QuantumPlane.hs
@@ -0,0 +1,236 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances, MultiParamTypeClasses, NoMonomorphismRestriction #-}
+
+-- |A module defining the quantum plane and its symmetries
+module Math.QuantumAlgebra.QuantumPlane where
+
+-- Refs:
+-- Kassel, Quantum Groups
+-- Street, Quantum Groups
+
+import Math.Algebra.Field.Base hiding (powers)
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+import Math.Algebras.LaurentPoly
+import Math.Algebras.NonCommutative
+import qualified Data.List as L
+
+
+
+qvar v = let V [(m,1)] = var v in V [(m,1 :: LaurentPoly Q)]
+
+
+a = qvar "a"
+b = qvar "b"
+c = qvar "c"
+d = qvar "d"
+
+detq = a*d-unit q'*b*c
+
+
+x = qvar "x"
+y = qvar "y"
+-- z = qvar "z"
+
+u = qvar "u"
+v = qvar "v"
+
+
+-- Quantum plane Aq20
+
+aq20 = [y*x-unit q*x*y]
+-- Kassel p72, Street p10
+
+newtype Aq20 v = Aq20 (NonComMonomial v) deriving (Eq,Ord)
+
+instance (Eq v, Show v) => Show (Aq20 v) where show (Aq20 m) = show m
+
+instance Monomial Aq20 where
+    var v = V [(Aq20 (NCM 1 [v]),1)]
+    powers (Aq20 m) = powers m
+
+instance Algebra (LaurentPoly Q) (Aq20 String) where
+    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'' = x' %% aq20 -- quotient by m2q relations while unwrapped
+        x''' = fmap Aq20 x'' -- wrap the monomials up as Aq20 again
+
+
+-- Quantum superplane Aq02
+
+aq02 = [u^2, v^2, u*v+unit q*v*u]
+-- Street p10
+
+newtype Aq02 v = Aq02 (NonComMonomial v) deriving (Eq,Ord)
+
+instance (Eq v, Show v) => Show (Aq02 v) where show (Aq02 m) = show m
+
+instance Monomial Aq02 where
+    var v = V [(Aq02 (NCM 1 [v]),1)]
+    powers (Aq02 m) = powers m
+
+instance Algebra (LaurentPoly Q) (Aq02 String) where
+    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'' = x' %% aq02 -- quotient by m2q relations while unwrapped
+        x''' = fmap Aq02 x'' -- wrap the monomials up as Aq02 again
+
+
+-- M2q
+
+m2q = [a*b-unit q'*b*a, a*c-unit q'*c*a, c*d-unit q'*d*c, b*d-unit q'*d*b,
+       b*c-c*b, a*d-d*a-unit (q'-q)*b*c]
+-- Kassel p78, Street p9
+-- I think this is already a Groebner basis
+
+newtype M2q v = M2q (NonComMonomial v) deriving (Eq,Ord)
+
+instance (Eq v, Show v) => Show (M2q v) where show (M2q m) = show m
+
+instance Monomial M2q where
+    var v = V [(M2q (NCM 1 [v]),1)]
+    powers (M2q m) = powers m
+
+instance Algebra (LaurentPoly Q) (M2q String) where
+    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'' = x' %% m2q -- quotient by m2q relations while unwrapped
+        x''' = fmap M2q x'' -- wrap the monomials up as M2q again
+
+-- Kassel p82-3
+instance Coalgebra (LaurentPoly Q) (M2q String) where
+    counit x = case x `bind` cu of
+               V [] -> 0
+               V [(M2q (NCM 0 []), c)] -> c
+        where cu "a" = 1 :: Vect (LaurentPoly Q) (M2q String)
+              cu "b" = 0
+              cu "c" = 0
+              cu "d" = 1
+    comult x = x `bind` cm
+        where cm "a" = a `te` a + b `te` c
+              cm "b" = a `te` b + b `te` d
+              cm "c" = c `te` a + d `te` c
+              cm "d" = c `te` b + d `te` d
+
+instance Bialgebra (LaurentPoly Q) (M2q String) where {}
+
+{-
+-- The following shows that the M2q relations are *sufficient*
+-- for M2q to be symmetries of Aq20 and Aq02
+
+> let x' = a*x+b*y :: Vect (LaurentPoly Q) (NonComMonomial String)
+> let y' = c*x+d*y :: Vect (LaurentPoly Q) (NonComMonomial String)
+> (y'*x'-unit q*x'*y') %% (m2q ++ aq20 ++ [s*t-t*s | s <- [a,b,c,d], t <- [x,y]])
+0
+
+> let u' = a*u+b*v :: Vect (LaurentPoly Q) (NonComMonomial String)
+> let v' = c*u+d*v :: Vect (LaurentPoly Q) (NonComMonomial String)
+> (u'^2) %% (m2q ++ aq02 ++ [s*t-t*s | s <- [a,b,c,d], t <- [u,v]])
+0
+> (v'^2) %% (m2q ++ aq02 ++ [s*t-t*s | s <- [a,b,c,d], t <- [u,v]])
+0
+> (u'*v'+unit q*v'*u') %% (m2q ++ aq02 ++ [s*t-t*s | s <- [a,b,c,d], t <- [u,v]])
+0
+
+-- To show that the M2q relations are necessary,
+-- set the coefficients of x^2, yx, y^2, and vu == 0 in all of the following
+> (y'*x'-unit q*x'*y') %% (aq20 ++ [p*q-q*p | p <- [a,b,c,d], q <- [x,y]])
+-qx^2ac+x^2ca-yxad-qyxbc+q^-1yxcb+yxda-qy^2bd+y^2db
+> (u'^2) %% (aq02 ++ [p*q-q*p | p <- [a,b,c,d], q <- [u,v]])
+-qvuab+vuba
+> (v'^2) %% (aq02 ++ [p*q-q*p | p <- [a,b,c,d], q <- [u,v]])
+-qvucd+vudc
+> (u'*v'+unit q*v'*u') %% (aq02 ++ [p*q-q*p | p <- [a,b,c,d], q <- [u,v]])
+-qvuad+vubc-q^2vucb+qvuda
+
+-- yx => -ad-qbc+q^-1cb+da == 0
+-- vu => -qad+bc-q^2cb+qda == 0
+-- qyx-vu => -q^2bc+cb-bc+q^2cb == 0 => bc == cb
+-- Now substitute back into yx
+
+-- We could probably have got gb to do this for us
+-}
+
+-- Kassel p85
+instance Comodule (LaurentPoly Q) (M2q String) (Aq20 String) where
+    coaction xy = xy `bind` ca where
+        ca "x" = (a `te` x) + (b `te` y) -- we can use (+) instead of add since Aq20 is an algebra 
+        ca "y" = (c `te` x) + (d `te` y)
+-- coaction (x) = (a b) `te` (x)
+--          (y)   (c d)      (y)
+
+
+-- SL2q
+
+sl2q = [a*b-unit q'*b*a, a*c-unit q'*c*a, c*d-unit q'*d*c, b*d-unit q'*d*b,
+        b*c-c*b, a*d-d*a-unit (q'-q)*b*c,
+        -unit q*c*b + d*a - 1] -- det q, but reduced
+--        a*d-unit q'*b*c-1] -- det_q
+-- We have to hand-reduce detq, or else call gb
+
+
+
+newtype SL2q v = SL2q (NonComMonomial v) deriving (Eq,Ord)
+
+instance (Eq v, Show v) => Show (SL2q v) where show (SL2q m) = show m
+
+instance Monomial SL2q where
+    var v = V [(SL2q (NCM 1 [v]),1)]
+    powers (SL2q m) = powers m
+
+instance Algebra (LaurentPoly Q) (SL2q String) where
+    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'' = x' %% sl2q -- quotient by sl2q relations while unwrapped
+        x''' = fmap SL2q x'' -- wrap the monomials up as SL2q again
+
+instance Coalgebra (LaurentPoly Q) (SL2q String) where
+    counit x = case x `bind` cu of
+               V [] -> 0
+               V [(SL2q (NCM 0 []), c)] -> c
+        where cu "a" = 1 :: Vect (LaurentPoly Q) (SL2q String)
+              cu "b" = 0
+              cu "c" = 0
+              cu "d" = 1
+    comult x = x `bind` cm
+        where cm "a" = a `te` a + b `te` c
+              cm "b" = a `te` b + b `te` d
+              cm "c" = c `te` a + d `te` c
+              cm "d" = c `te` b + d `te` d
+
+instance Bialgebra (LaurentPoly Q) (SL2q String) where {}
+
+-- Kassel p84
+instance HopfAlgebra (LaurentPoly Q) (SL2q String) where
+    antipode x = x `bind` antipode'
+        where antipode' "a" = d
+              antipode' "b" = - unit q * b
+              antipode' "c" = - unit q' * c
+              antipode' "d" = a
+-- in the GL2q case we would need 1/detq factor as well
+
+
+
+-- !! The following probably needs to be rehoused in separate module at some point
+-- YANG-BAXTER OPERATOR
+
+-- 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)
+
+
diff --git a/Math/QuantumAlgebra/Tangle.hs b/Math/QuantumAlgebra/Tangle.hs
new file mode 100644
--- /dev/null
+++ b/Math/QuantumAlgebra/Tangle.hs
@@ -0,0 +1,177 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE MultiParamTypeClasses, TypeFamilies, FlexibleInstances, EmptyDataDecls #-}
+
+-- |A module defining the category of tangles, and representations into the category of vector spaces
+-- (specifically, knot invariants).
+module Math.QuantumAlgebra.Tangle where
+
+-- import qualified Data.List as L
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+
+import Math.Algebra.Field.Base
+import Math.Algebras.LaurentPoly
+
+import Math.QuantumAlgebra.TensorCategory
+
+
+instance Mon [a] where
+    munit = []
+    mmult = (++)
+
+-- type TensorAlgebra k a = Vect k [a]
+
+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)
+
+-- Could make TensorAlgebra k a into an instance of Category, TensorCategory
+    
+
+-- TANGLE CATEGORY
+-- (Unoriented)
+
+data Tangle
+
+instance Category Tangle where
+    data Ob Tangle = OT Int deriving (Eq,Ord,Show)
+    data Ar Tangle = IdT Int
+                   | CapT
+                   | CupT
+                   | OverT
+                   | UnderT
+--                   | SeqT (Ar Tangle) (Ar Tangle)
+                   | SeqT [Ar Tangle]
+--                   | ParT (Ar Tangle) (Ar Tangle)
+                   | ParT [Ar Tangle]
+                   deriving (Eq,Ord,Show)
+    id_ (OT n) = IdT n
+    source (IdT n) = OT n
+    source CapT = OT 0
+    source CupT = OT 2
+    source OverT = OT 2
+    source UnderT = OT 2
+--    source (ParT a b) = OT (sa + sb) where OT sa = source a; OT sb = source b
+    source (ParT as) = OT $ sum [sa | a <- as, let OT sa = source a]
+--    source (SeqT a b) = source a
+    source (SeqT as) = source (head as)
+    target (IdT n) = OT n
+    target CapT = OT 2
+    target CupT = OT 0
+    target OverT = OT 2
+    target UnderT = OT 2
+--    target (ParT a b) = OT (ta + tb) where OT ta = target a; OT tb = target b
+    target (ParT as) = OT $ sum [ta | a <- as, let OT ta = target a]
+--    target (SeqT a b) = target b
+    target (SeqT as) = target (last as)
+--    a >>> b | target a == source b = SeqT a b
+    a >>> b | target a == source b = SeqT [a,b]
+
+instance TensorCategory Tangle where
+    tunit = OT 0
+    tob (OT a) (OT b) = OT (a+b)
+--    tar a b = ParT a b
+    tar a b = ParT [a,b]
+
+
+
+-- KAUFFMAN BRACKET
+
+data Oriented = Plus | Minus deriving (Eq,Ord,Show)
+
+type TangleRep b = Vect (LaurentPoly Q) b
+
+
+-- adapted from http://blog.sigfpe.com/2008/10/untangling-with-continued-fractions.html
+cap :: [Oriented] -> TangleRep [Oriented]
+cap [] = return [Plus, Minus] <+> (-q^2) *> return [Minus, Plus]
+
+cup :: [Oriented] -> TangleRep [Oriented]
+cup [Plus, Minus] = (-q'^2) *> return []
+cup [Minus, Plus] = return []
+cup _ = zero
+
+-- also called xminus
+over :: [Oriented] -> TangleRep [Oriented]
+over [u, v] = q  *> do {[] <- cup [u, v]; cap []}
+          <+> q' *> return [u, v]
+
+{-
+-- if you expand "over" into terms, you find that it equals the following,
+-- which strongly resembles c' below
+over' (T i j) = case compare i j of
+                EQ -> q' *> return (T i i)                                       -- ++ -> q' ++, -- -> q' -- 
+                LT -> q  *> return (T j i)                                       -- +- -> q -+
+                GT -> q  *> (return (T j i) <+> (q'^2 - q^2) *> return (T i j))  -- -+ -> q +- + (q'-q^3) -+
+-}
+-- also called xplus
+under :: [Oriented] -> TangleRep [Oriented]
+under [u, v] = q' *> do {[] <- cup [u, v]; cap []}
+           <+> q  *> return [u, v]
+
+{-
+-- if you expand "under" into terms, you find that it equals the following,
+-- which strongly resembles c below
+under' (T i j) = case compare i j of
+                 EQ -> q  *> return (T i i)                                       -- ++ -> q ++, -- -> q -- 
+                 LT -> q' *> (return (T j i) <+> (q^2 - q'^2) *> return (T i j))  -- +- -> q' -+ + (q-q^-3) -+
+                 GT -> q' *> return (T j i)                                       -- -+ -> q' +-
+-}
+loop = nf $ do {[i, j] <- cap []; cup [i, j]}
+
+trefoil = nf $ do
+    [i, j] <- cap []
+    [k, l] <- cap []
+    [m, n] <- under [j, k]
+    [p, q] <- over [i, m]
+    [r, s] <- over [n, l]
+    cup [p, s]
+    cup [q, r]
+
+
+-- KAUFFMAN BRACKET AS A REPRESENTATION FROM TANGLE TO VECT
+
+-- But this isn't quite the Kauffman bracket - we still need to divide by (-q^2-q^-2)
+kauffman :: Ar Tangle -> TangleRep [Oriented] -> TangleRep [Oriented]
+kauffman (IdT n) = id -- could be tf of n ids
+kauffman CapT = linear cap
+kauffman CupT = linear cup
+kauffman OverT = linear over
+kauffman UnderT = linear under
+kauffman (SeqT fs) = foldl (>>>) id $ map kauffman fs
+    where g >>> h = h . g
+kauffman (ParT [f]) = kauffman f
+kauffman (ParT (f:fs)) = tf m (kauffman f) (kauffman (ParT fs))
+    where OT m = source f
+          tf m f' fs' = linear (\xs -> let (ls,rs) = splitAt m xs in f' (return ls) * fs' (return rs) )
+{-
+kauffman (ParT f g) = tf m n (kauffman f) (kauffman g)
+    where OT m = source f
+          OT n = source g
+          tf m n f' g' = linear (\xs -> let (ls,rs) = splitAt m xs in f' (return ls) * g' (return rs) )
+-}
+
+-- loopT = SeqT CapT CupT
+loopT = SeqT [CapT, CupT]
+
+{-
+trefoilT = (ParT CapT CapT) `SeqT` (ParT (IdT 1) (ParT UnderT (IdT 1)))
+    `SeqT` (ParT OverT OverT) `SeqT` (ParT (IdT 1) (ParT CupT (IdT 1))) `SeqT` CupT
+
+trefoilT = ParT [CapT, CapT]
+    `SeqT` ParT [IdT 1, UnderT, IdT 1]
+    `SeqT` ParT [OverT, OverT]
+    `SeqT` ParT [IdT 1, CupT, IdT 1]
+    `SeqT` CupT
+-}
+trefoilT = SeqT [
+    ParT [CapT, CapT],
+    ParT [IdT 1, UnderT, IdT 1],
+    ParT [OverT, OverT],
+    ParT [IdT 1, CupT, IdT 1],
+    CupT]
+-- eg kauffman (trefoilT) (return [])
diff --git a/Math/QuantumAlgebra/TensorCategory.hs b/Math/QuantumAlgebra/TensorCategory.hs
new file mode 100644
--- /dev/null
+++ b/Math/QuantumAlgebra/TensorCategory.hs
@@ -0,0 +1,128 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE TypeFamilies, EmptyDataDecls #-}
+
+-- |A module defining classes and example instances of categories and tensor categories
+module Math.QuantumAlgebra.TensorCategory where
+
+import Math.Algebra.Group.PermutationGroup
+
+
+class Category c where
+    data Ob c :: *
+    data Ar c :: *
+    id_ :: Ob c -> Ar c
+    source, target :: Ar c -> Ob c
+    (>>>) :: Ar c -> Ar c -> Ar c
+
+-- Kassel p282
+-- The following is actually definition of a strict tensor category
+class Category c => TensorCategory c where
+    tunit :: Ob c
+    tob :: Ob c -> Ob c -> Ob c -- tensor product of objects
+    tar :: Ar c -> Ar c -> Ar c -- tensor product of arrows
+
+class TensorCategory c => StrictTensorCategory c where {}
+-- we want to be able to declare some tensor categories as strict
+
+class TensorCategory c => WeakTensorCategory c where
+    assoc :: Ob c -> Ob c -> Ob c -- (u `tob` v) `tob` w -> u `tob` (v `tob` w)
+    lunit :: Ob c -> Ob c         -- unit `tob` v -> v
+    runit :: Ob c -> Ob c         -- v `tob` unit -> v
+
+{-
+instance (TensorCategory c, Eq (Ar c), Show (Ar c)) => Num (Ar c) where
+    (*) = tar
+-}
+
+-- SYMMETRIC GROUPOID
+
+data SymmetricGroupoid
+
+instance Category SymmetricGroupoid where
+    data Ob SymmetricGroupoid = OS Int deriving (Eq,Ord,Show)
+    data Ar SymmetricGroupoid = AS Int (Permutation Int) deriving (Eq,Ord,Show)
+    id_ (OS n) = AS n 1
+    source (AS n _) = OS n
+    target (AS n _) = OS n
+    AS m g >>> AS n h | m == n = AS m (g*h)
+
+instance TensorCategory SymmetricGroupoid where
+    tunit = OS 0
+    tob (OS m) (OS n) = OS (m+n)
+    tar (AS m g) (AS n h) = AS (m+n) (g * h~^k)
+        where k = p [[1..m+n]] ^ m
+--    tar (AS m g) (AS n h) = AS (m+n) (fromPairs $ toPairs g ++ map (\(x,y)->(x+m,y+m)) (toPairs h))
+
+
+-- BRAID CATEGORY
+
+data Braid
+
+instance Category Braid where
+    data Ob Braid = OB Int deriving (Eq,Ord,Show)
+    data Ar Braid = AB Int [Int] deriving (Eq,Ord,Show)
+    id_ (OB n) = AB n []
+    source (AB n _) = OB n
+    target (AB n _) = OB n
+    AB m is >>> AB n js | m == n = AB m (is ++ js)
+
+s n i | 0 < i && i < n = AB n [i]
+
+instance TensorCategory Braid where
+    tunit = OB 0
+    tob (OB a) (OB b) = OB (a+b)
+    tar (AB m is) (AB n js) = AB (m+n) (is ++ map (+m) js)
+
+
+
+
+data Cob2
+-- works very similar to Tangle category
+
+instance Category Cob2 where
+    data Ob Cob2 = O Int deriving (Eq,Ord,Show)
+    data Ar Cob2 = Id Int
+                 | Unit
+                 | Mult
+                 | Counit
+                 | Comult
+                 | Par (Ar Cob2) (Ar Cob2)
+                 | Seq (Ar Cob2) (Ar Cob2)
+                 deriving (Eq,Ord,Show)
+    id_ (O n) = Id n
+    source (Id n) = O n
+    source Unit = O 0
+    source Mult = O 2
+    source Counit = O 1
+    source Comult = O 1
+    source (Par a b) = O (sa + sb) where O sa = source a; O sb = source b
+    source (Seq a b) = source a
+    target (Id n) = O n
+    target Unit = O 1
+    target Mult = O 1
+    target Counit = O 0
+    target Comult = O 2
+    target (Par a b) = O (ta + tb) where O ta = target a; O tb = target b
+    target (Seq a b) = target b
+    a >>> b | target a == source b = Seq a b
+
+instance TensorCategory Cob2 where
+    tunit = O 0
+    tob (O a) (O b) = O (a+b)
+    tar a b = Par a b
+
+-- rewrite a Cob2 so that it is a Seq of Pars
+-- (this isn't necessarily going to help us towards a normal form - there may not even be one
+rewrite (Par (Seq a1 a2) (Seq b1 b2)) =
+    Seq (Par idSourceA b1')
+        ( (Seq (Par idSourceA b2')
+               (Seq (Par a1' idTargetB)
+                    (Par a2' idTargetB) ) ) )
+    where idSourceA = id_ (source a1)
+          idTargetB = id_ (target b2)
+          a1' = rewrite a1
+          a2' = rewrite a2
+          b1' = rewrite b1
+          b2' = rewrite b2
+rewrite x = x
diff --git a/Math/Test/TAlgebras/TGroupAlgebra.hs b/Math/Test/TAlgebras/TGroupAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TAlgebras/TGroupAlgebra.hs
@@ -0,0 +1,55 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}
+
+module Math.Test.TAlgebras.TGroupAlgebra where
+
+import Test.QuickCheck
+
+import Math.Algebra.Group.PermutationGroup
+import Math.Test.TPermutationGroup -- for instance Arbitrary (Permutation Int)
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Structures
+import Math.Algebras.GroupAlgebra
+
+import Math.Test.TAlgebras.TStructures
+
+{-
+instance Arbitrary (TensorAlgebra Integer) where
+    arbitrary = do ts <- arbitrary :: Gen [([Int], Integer)]
+                   return $ nf $ V ts
+-}
+instance Arbitrary (GroupAlgebra Integer) where
+    arbitrary = do ts <- arbitrary :: Gen [(Permutation Int, Integer)]
+                   return $ nf $ V ts
+
+
+{-
+prop_Algebra_TensorAlgebra (k,x,y,z) = prop_Algebra (k,x,y,z)
+    where types = (k,x,y,z) :: (Integer, TensorAlgebra Integer, TensorAlgebra Integer, TensorAlgebra Integer)
+
+prop_Coalgebra_TensorAlgebra x = prop_Coalgebra x
+    where types = x :: TensorAlgebra Integer
+-}
+
+prop_Algebra_GroupAlgebra (k,x,y,z) = prop_Algebra (k,x,y,z)
+    where types = (k,x,y,z) :: (Integer, GroupAlgebra Integer, GroupAlgebra Integer, GroupAlgebra Integer)
+
+-- have to split the 8-tuple into two 4-tuples to avoid having to write Arbitrary instance
+prop_Algebra_Linear_GroupAlgebra ((k,l,m,n),(x,y,z,w)) = prop_Algebra_Linear (k,l,m,n,x,y,z,w)
+    where types = (k,l,m,n,x,y,z,w) :: (Integer, Integer, Integer, Integer,
+                   GroupAlgebra Integer, GroupAlgebra Integer, GroupAlgebra Integer, GroupAlgebra Integer)
+
+prop_Coalgebra_GroupAlgebra x = prop_Coalgebra x
+    where types = x :: GroupAlgebra Integer
+
+prop_Coalgebra_Linear_GroupAlgebra (k,l,x,y) = prop_Coalgebra_Linear (k,l,x,y)
+    where types = (k,l,x,y) :: (Integer, Integer, GroupAlgebra Integer, GroupAlgebra Integer)
+
+prop_Bialgebra_GroupAlgebra (k,x,y) = prop_Bialgebra (k,x,y)
+    where types = (k,x,y) :: (Integer, GroupAlgebra Integer, GroupAlgebra Integer)
+
+prop_HopfAlgebra_GroupAlgebra x = prop_HopfAlgebra x
+    where types = x :: GroupAlgebra Integer
diff --git a/Math/Test/TAlgebras/TQuaternions.hs b/Math/Test/TAlgebras/TQuaternions.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TAlgebras/TQuaternions.hs
@@ -0,0 +1,39 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}
+
+module Math.Test.TAlgebras.TQuaternions where
+
+import Test.QuickCheck
+
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Quaternions
+
+import Math.Test.TAlgebras.TStructures
+
+instance Arbitrary HBasis where
+    arbitrary = elements [One,I,J,K]
+
+instance Arbitrary (Quaternion Integer) where
+    arbitrary = do ts <- arbitrary :: Gen [(HBasis, Integer)]
+                   return $ nf $ V ts
+
+
+
+prop_Algebra_Quaternion (k,x,y,z) = prop_Algebra (k,x,y,z)
+    where types = (k,x,y,z) :: (Integer, Quaternion Integer, Quaternion Integer, Quaternion Integer)
+
+prop_Coalgebra_Quaternion x = prop_Coalgebra x
+    where types = x :: Quaternion Integer
+
+-- Fails - the algebra and coalgebra structures I've given are not compatible
+prop_Bialgebra_Quaternion (k,x,y) = prop_Bialgebra (k,x,y)
+    where types = (k,x,y) :: (Integer, Quaternion Integer, Quaternion Integer)
+
+{-
+prop_FrobeniusRelation_Quaternion (x,y) = prop_FrobeniusRelation (x,y)
+    where types = (x,y) :: (Quaternion Integer, Quaternion Integer)
+-- !! fails, because the counit we have given is not a Frobenius form
+-}
+
diff --git a/Math/Test/TAlgebras/TStructures.hs b/Math/Test/TAlgebras/TStructures.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TAlgebras/TStructures.hs
@@ -0,0 +1,211 @@
+-- Copyright (c) 2010, David Amos. All rights reserved.
+
+{-# LANGUAGE NoMonomorphismRestriction #-}
+
+
+module Math.Test.TAlgebras.TStructures where
+
+-- import Test.QuickCheck
+-- don't actually need, as we don't define any Arbitrary instances here
+
+import Control.Arrow ( (>>>) ) -- actually you can get this from Category?
+
+
+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
+
+-- in this version we supply z of the intended return type of f,
+-- so that we can make sure we select the correct instance for f polymorphic in return type
+prop_Linear' f (k,l,x,y,z) =
+    f (add (smultL k x) (smultL l y)) `add` z == add (smultL k (f x)) (smultL l (f y)) `add` z
+
+
+-- prop_Bilinear could be defined in terms of prop_Linear over tensor product
+-- if we had a way to convert a bilinear function to a tensor function
+
+prop_Algebra_Linear ::
+    (Num k, Ord b, Algebra k b) =>
+    (k, k, k, k, Vect k b, Vect k b, Vect k b, Vect k b) -> Bool
+prop_Algebra_Linear (k,l,m,n,x,y,z,w) =
+--    (unit (k * m + l * n) :: Vect k b) == (add (smultL k (unit m)) (smultL l (unit n)) :: Vect k b) &&
+    prop_Linear' unit' (k,l,wrap m, wrap n, x) &&
+    prop_Linear mult (k,l, x `te` y, z `te` w)
+    where wrap = (\c -> V [((),c)]) :: k -> Trivial k
+
+prop_Coalgebra_Linear (k,l,x,y) =
+    prop_Linear counit' (k,l,x,y) &&
+    prop_Linear comult (k,l,x,y)
+-- now need instances for GroupAlgebra etc
+
+
+-- ALGEBRAS
+
+prop_Algebra (k,x,y,z) =
+    mult (x `te` mult (y `te` z)) == mult (mult (x `te` y) `te` z)  && -- associativity
+    smultL k x == mult (unit k `te` x)                             && -- left unit
+    -- mult (k' `te` x) == (mult . (unit' `tf` id)) (k' `te` x)         && -- left unit
+    smultR x k == mult (x `te` unit k)                        -- && -- right unit
+    -- mult (x `te` k') == (mult . (id `tf` unit')) (x `te` k')            -- right unit
+    where k' = V [( (),k)]
+-- additionally, unit and mult must be linear
+
+prop_Commutative (x,y) =
+    let xy = x `te` y
+    in (mult . twist) xy == mult xy
+
+
+-- COALGEBRAS
+
+prop_Coalgebra x =
+    ((comult `tf` id) . comult) x == (assocL . (id `tf` comult) . comult) x && -- coassociativity
+    ((counit' `tf` id) . comult) x == V [((),1)] `te` x                     && -- left counit
+    ((id `tf` counit') . comult) x == x `te` V [((),1)]                        -- right counit
+-- additionally, counit and comult must be linear
+
+prop_Cocommutative x =
+    (twist . comult) x == comult x
+
+
+-- MORPHISMS
+
+prop_AlgebraMorphism f (k,l,x,y) =
+    prop_Linear f (k,l,x,y) &&
+    -- (f . unit) k == unit k &&
+    (f . mult) (x `te` y) == (mult . (f `tf` f)) (x `te` y) 
+
+-- in this version we supply z of the intended return type of f,
+-- so that we can make sure we select the correct instance for f polymorphic in return type
+prop_AlgebraMorphism' f (k,l,x,y,z) =
+    prop_Linear f (k,l,x,y) &&
+    (f . unit) k + z == unit k + z &&
+    (f . mult) (x `te` y) == (mult . (f `tf` f)) (x `te` y) 
+
+prop_CoalgebraMorphism f x =
+    -- prop_Linear f (k,l,x,y) &&
+    (counit . f) x == counit x &&
+    ( (f `tf` f) . comult) x == (comult . f) x
+
+
+-- BIALGEBRAS
+
+prop_Bialgebra1 (x,y) =
+    let xy = x `te` y in
+    (comult . mult) xy ==
+    ( (mult `tf` mult) .
+      assocL . (id `tf` assocR) .
+      (id `tf` (twist `tf` id)) .
+      (id `tf` assocL) . assocR .
+      (comult `tf` comult) ) xy
+
+prop_Bialgebra2 (k,xy) =
+    (comult . unit') k' + xy == ((unit' `tf` unit') . iso) k' + xy
+    where iso = fmap (\ () -> T () () ) -- 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 () ()) -> ())
+
+prop_Bialgebra4 (k,x) =
+    id k == (counit . (\a -> a+x-x) . unit) k
+-- so we are using the x just to force the intermediate value to be of the right type
+
+prop_Bialgebra (k,x,y) =
+    prop_Bialgebra1 (x,y) &&
+    prop_Bialgebra2 (k,x `te` y) &&
+    prop_Bialgebra3 (x,y) &&
+    prop_Bialgebra4 (k,x)
+
+
+prop_HopfAlgebra x =
+    (unit . counit) x == (mult . (antipode `tf` id) . comult) x &&
+    (unit . counit) x == (mult . (id `tf` antipode) . comult) x
+
+-- Street p87
+-- we also require that f be invertible
+prop_YangBaxter f (x,y,z) =
+    ( (f `tf` id) >>> assocR >>> (id `tf` f) >>> assocL >>> (f `tf` id) >>> assocR ) xyz == 
+    ( assocR >>> (id `tf` f) >>> assocL >>> (f `tf` id) >>> assocR >>> (id `tf` f) ) xyz
+    where xyz = ( (x `te` y) `te` z) 
+
+
+
+-- MODULES AND COMODULES
+
+prop_Module_Linear (k,l,x,y) = prop_Linear action (k,l,x,y)
+
+prop_Module_Assoc (r,s,m) =
+    (action . (mult `tf` id)) ((r `te` s) `te` m) == (action . (id `tf` action)) (r `te` (s `te` m))
+{-
+prop_Module_Unit (k,m) =
+    (action . (unit' `tf` id)) k' ==  
+-}
+
+
+
+
+
+
+
+
+
+
+
+
+
+frobeniusLeft1 = (id `tf` mult) . assocR . (comult `tf` id)
+
+frobeniusLeft2 x = nf $ x >>= fl
+    where fl (T i j) = do
+              T k l <- comultM i
+              m <- idM j
+              p <- idM k
+              q <- multM (T l m)
+              return (T 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)
+
+prop_FrobeniusRelation (x,y) =
+    let xy = x `te` y
+    in frobeniusLeft1 xy == frobeniusMiddle1 xy
+
+-- (inject == return)
+
+multM = mult . return -- inject
+comultM = comult . return -- inject
+idM = id . return
+
+-- can we do the same with unit, counit?
+-- 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'
+
+
+
+
+