diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal
--- a/HaskellForMaths.cabal
+++ b/HaskellForMaths.cabal
@@ -1,13 +1,14 @@
    Name:                HaskellForMaths
-   Version:             0.4.1
+   Version:             0.4.2
    Category:            Math
-   Description:         A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended for educational purposes, but does have efficient implementations of several fundamental algorithms.
+   Description:         A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended as an educational resource, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra
    License:             BSD3
    License-file:        license.txt
    Author:              David Amos
    Maintainer:          haskellformaths-at-gmail-dot-com
    Homepage:            http://haskellformaths.blogspot.com/
+   Stability:           experimental
    Build-Type:          Simple
    Cabal-Version:       >=1.2
 
@@ -25,10 +26,13 @@
         Math/Test/TAlgebras/TTensorProduct.hs
         Math/Test/TAlgebras/TStructures.hs
         Math/Test/TAlgebras/TQuaternions.hs
+        Math/Test/TAlgebras/TOctonions.hs
         Math/Test/TAlgebras/TMatrix.hs
         Math/Test/TAlgebras/TGroupAlgebra.hs
         Math/Test/TCombinatorics/TPoset.hs
         Math/Test/TCombinatorics/TDigraph.hs
+        Math/Test/TCombinatorics/TFiniteGeometry.hs
+        Math/Test/TCombinatorics/TGraphAuts.hs
         Math/Test/TCombinatorics/TIncidenceAlgebra.hs
         Math/Test/TCombinatorics/TMatroid.hs
         Math/Test/TCommutativeAlgebra/TPolynomial.hs
diff --git a/Math/Algebra/Field/Base.hs b/Math/Algebra/Field/Base.hs
--- a/Math/Algebra/Field/Base.hs
+++ b/Math/Algebra/Field/Base.hs
@@ -1,12 +1,12 @@
 -- Copyright (c) David Amos, 2008. All rights reserved.
 
-{-# OPTIONS_GHC -fglasgow-exts #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving, ScopedTypeVariables #-}
 
 module Math.Algebra.Field.Base where
 
 import Data.Ratio
 import Math.Common.IntegerAsType
-
+import Math.Core.Utils
 
 -- RATIONALS
 
@@ -55,13 +55,17 @@
                    in Fp $ u `mod` p
                    where p = value (undefined :: n)
 
-class Fractional fq => FiniteField fq where
+-- Not sure if Eq fq is required, need to try with ghc >= 7.4.1
+class (Eq fq, Fractional fq) => FiniteField fq where
     eltsFq :: fq -> [fq]  -- return all elts of the field
     basisFq :: fq -> [fq] -- return an additive basis for the field (as Z-module)
 
 instance IntegerAsType p => FiniteField (Fp p) where
     eltsFq _ = map fromInteger [0..p'-1] where p' = value (undefined :: p)
     basisFq _ = [fromInteger 1]
+
+instance IntegerAsType p => FinSet (Fp p) where
+    elts = map fromInteger [0..p'-1] where p' = value (undefined :: p)
 
 primitiveElt fq = head [x | x <- tail fq, length (powers x) == q-1] where
     q = length fq
diff --git a/Math/Algebra/Field/Extension.hs b/Math/Algebra/Field/Extension.hs
--- a/Math/Algebra/Field/Extension.hs
+++ b/Math/Algebra/Field/Extension.hs
@@ -1,6 +1,6 @@
 -- Copyright (c) David Amos, 2008. All rights reserved.
 
-{-# OPTIONS_GHC -fglasgow-exts #-}
+{-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, ScopedTypeVariables, EmptyDataDecls #-}
 
 module Math.Algebra.Field.Extension where
 
@@ -8,6 +8,7 @@
 import Data.List as L (elemIndex)
 
 import Math.Common.IntegerAsType
+import Math.Core.Utils
 import Math.Algebra.Field.Base
 
 
@@ -136,6 +137,12 @@
         d = deg $ pvalue (undefined :: (k,poly))
     basisFq _ = map embed $ take (d-1) $ iterate (*x) 1 where
         d = deg $ pvalue (undefined :: (k,poly))
+
+-- Not sure if Eq fp is required, need to check with ghc >= 7.4.1
+instance (FinSet fp, Eq fp, Num fp, PolynomialAsType fp poly) => FinSet (ExtensionField fp poly) where
+    elts = map Ext (polys (d-1) fp') where
+        fp' = elts
+        d = deg $ pvalue (undefined :: (fp,poly))
 
 embed f = Ext (convert f)
 
diff --git a/Math/Algebra/Group/CayleyGraph.hs b/Math/Algebra/Group/CayleyGraph.hs
--- a/Math/Algebra/Group/CayleyGraph.hs
+++ b/Math/Algebra/Group/CayleyGraph.hs
@@ -4,6 +4,8 @@
 
 module Math.Algebra.Group.CayleyGraph where
 
+import Math.Core.Utils hiding (elts)
+
 import Math.Algebra.Group.StringRewriting as SR
 import Math.Combinatorics.Graph
 -- import Math.Combinatorics.GraphAuts
@@ -11,9 +13,6 @@
 import Math.Algebra.Group.PermutationGroup as P
 
 import qualified Data.List as L
-import qualified Data.Set as S
-
-toSet = S.toList . S.fromList
 
 
 data Digraph a = DG [a] [(a,a)] deriving (Eq,Ord,Show)
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
@@ -1,5 +1,15 @@
--- Copyright (c) David Amos, 2008-2009. All rights reserved.
+-- Copyright (c) David Amos, 2008-2012. All rights reserved.
 
+{-# LANGUAGE NoMonomorphismRestriction #-}
+
+-- |A module for doing arithmetic in permutation groups.
+--
+-- Group elements are represented as permutations of underlying sets, and are entered and displayed
+-- using a Haskell-friendly version of cycle notation. For example, the permutation (1 2 3)(4 5)
+-- would be entered as @p [[1,2,3],[4,5]]@, and displayed as [[1,2,3],[4,5]]. Permutations can be defined
+-- over arbitrary underlying sets (types), not just the integers.
+--
+-- If @g@ and @h@ are group elements, then the expressions @g*h@ and @g^-1@ calculate product and inverse respectively.
 module Math.Algebra.Group.PermutationGroup where
 
 import qualified Data.List as L
@@ -8,8 +18,11 @@
 
 import Math.Common.ListSet (toListSet, union, (\\) ) -- a version of union which assumes the arguments are ascending sets (no repeated elements)
 
-infix 8 ^-, ~^
+import Math.Core.Utils hiding (elts)
+import Math.Algebra.LinearAlgebra hiding (inverse) -- only needed for use in ghci
 
+infix 8 ~^
+
 rotateL (x:xs) = xs ++ [x]
 
 
@@ -18,6 +31,11 @@
 -- |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)
 
+-- |Construct a permutation from a list of cycles.
+-- For example, @p [[1,2,3],[4,5]]@ returns the permutation that sends 1 to 2, 2 to 3, 3 to 1, 4 to 5, 5 to 4.
+p :: (Ord a) => [[a]] -> Permutation a
+p = fromCycles
+
 fromPairs xys | isValid   = fromPairs' xys
               | otherwise = error "Not a permutation"
     where (xs,ys) = unzip xys
@@ -37,27 +55,23 @@
 -- (This is guaranteed not to contain fixed points provided the permutations have been constructed using the supplied constructors)
 
 -- |x .^ g returns the image of a vertex or point x under the action of the permutation g.
+-- For example, @1 .^ p [[1,2,3]]@ returns 2.
 -- The dot is meant to be a mnemonic for point or vertex.
-(.^) :: (Ord k) => k -> Permutation k -> k
+(.^) :: (Ord a) => a -> Permutation a -> a
 x .^ P g = case M.lookup x g of
            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.
+-- For example, @[1,2] -^ p [[1,4],[2,3]]@ returns [3,4].
 -- The dash is meant to be a mnemonic for edge or line or block.
-(-^) :: (Ord t) => [t] -> Permutation t -> [t]
+(-^) :: (Ord a) => [a] -> Permutation a -> [a]
 xs -^ g = L.sort [x .^ g | x <- xs]
 
 -- construct a permutation from cycles
 fromCycles cs = fromPairs $ concatMap fromCycle cs
     where fromCycle xs = zip xs (rotateL xs)
 
--- |Construct a permutation from a list of cycles.
--- For example, p [[1,2,3],[4,5]] returns the permutation that sends 1 to 2, 2 to 3, 3 to 1, 4 to 5, 5 to 4
-p :: (Ord a) => [[a]] -> Permutation a
-p cs = fromCycles cs
--- can't specify in pointfree style because of monomorphism restriction
-
 -- convert a permutation to cycles
 toCycles g = toCycles' $ supp g
     where toCycles' ys@(y:_) = let c = cycleOf g y in c : toCycles' (ys L.\\ c)
@@ -77,23 +91,25 @@
 orderElt g = foldl lcm 1 $ map length $ toCycles g
 -- == order [g]
 
+-- |The Num instance is what enables us to write @g*h@ for the product of group elements and @1@ for the group identity.
+-- Unfortunately we can't of course give sensible definitions for the other functions declared in the Num typeclass.
 instance (Ord a, Show a) => Num (Permutation a) where
     g * h = fromPairs' [(x, x .^ g .^ h) | x <- supp g `union` supp h]
     -- signum = sign -- doesn't work, complains about no (+) instance
     fromInteger 1 = P $ M.empty
+    _ + _ = error "(Permutation a).+: not applicable"
+    negate _ = error "(Permutation a).negate: not applicable"
+    abs _ = error "(Permutation a).abs: not applicable"
+    signum _ = error "(Permutation a).signum: not applicable"
 
-inverse (P g) = P $ M.fromList $ map (\(x,y)->(y,x)) $ M.toList g
+-- |The HasInverses instance is what enables us to write @g^-1@ for the inverse of a group element.
+instance (Ord a, Show a) => HasInverses (Permutation a) where
+    inverse (P g) = P $ M.fromList $ map (\(x,y)->(y,x)) $ M.toList g
 
--- |A trick: g^-1 returns the inverse of g
-(^-) :: (Ord k, Show k) => Permutation k -> Int -> Permutation k
-g ^- n = inverse g ^ n
 
-instance (Ord a, Show a) => Fractional (Permutation a) where
-    recip = inverse
-
--- |g ~^ h returns the conjugate of g by h.
+-- |g ~^ h returns the conjugate of g by h, that is, h^-1*g*h.
 -- The tilde is meant to a mnemonic, because conjugacy is an equivalence relation.
-(~^) :: (Ord t, Show t) => Permutation t -> Permutation t -> Permutation t
+(~^) :: (Ord a, Show a) => Permutation a -> Permutation a -> Permutation a
 g ~^ h = h^-1 * g * h
 
 -- commutator
@@ -122,7 +138,7 @@
 orbitV gs x = orbit (.^) x gs
 
 -- |b -^^ gs returns the orbit of the block or edge b under the action of the gs
-(-^^) :: (Ord t) => [t] -> [Permutation t] -> [[t]]
+(-^^) :: (Ord a) => [a] -> [Permutation a] -> [[a]]
 b -^^ gs = orbit (-^) b gs
 
 orbitB gs b = orbit (-^) b gs
@@ -174,9 +190,9 @@
           t = p [[1,2,3]]
 
 
--- Direct product of groups
--- Given generators for H and K, acting on sets X and Y respectively,
--- return generators for H*K, acting on the disjoint union X+Y (== Either X Y)
+-- |Given generators for groups H and K, acting on sets A and B respectively,
+-- return generators for the direct product H*K, acting on the disjoint union A+B (= Either A B)
+dp :: (Ord a, Ord b) => [Permutation a] -> [Permutation b] -> [Permutation (Either a b)]
 dp hs ks =
     [P $ M.fromList $ map (\(x,x') -> (Left x,Left x')) $ M.toList h' | P h' <- hs] ++
     [P $ M.fromList $ map (\(y,y') -> (Right y,Right y')) $ M.toList k' | P k' <- ks]
@@ -300,9 +316,9 @@
 conjClass gs h = closure [h] [ (~^ g) | g <- gs]
 -- conjClass gs h = h ~^^ gs
 
--- |conjClassReps gs returns a conjugacy class representatives and sizes for the group generated by gs.
+-- |conjClassReps gs returns conjugacy class representatives and sizes for the group generated by gs.
 -- This implementation is only suitable for use with small groups (|G| < 10000).
-conjClassReps :: (Ord t, Show t) => [Permutation t] -> [(Permutation t, Int)]
+conjClassReps :: (Ord a, Show a) => [Permutation a] -> [(Permutation a, Int)]
 conjClassReps gs = conjClassReps' (elts gs) where
     conjClassReps' (h:hs) =
         let cc = conjClass gs h in (h, length cc) : conjClassReps' (hs \\ cc)
@@ -440,4 +456,9 @@
 -- in cube gp, the subgps all appear to correspond to stabilisers of subsets, or of blocks
 
 
+-- right regular permutation representation
+rrpr gs h = rrpr' (elts gs) h
 
+rrpr' gs h = fromPairs [(g, g*h) | g <- gs]
+
+permutationMatrix xs g = [ [if x .^ g == y then 1 else 0 | y <- xs] | x <- xs ]
diff --git a/Math/Algebra/Group/RandomSchreierSims.hs b/Math/Algebra/Group/RandomSchreierSims.hs
--- a/Math/Algebra/Group/RandomSchreierSims.hs
+++ b/Math/Algebra/Group/RandomSchreierSims.hs
@@ -13,6 +13,7 @@
 import System.IO.Unsafe
 
 import Math.Common.ListSet (toListSet)
+import Math.Core.Utils hiding (elts)
 import Math.Algebra.Group.PermutationGroup
 import Math.Algebra.Group.SchreierSims (sift, cosetRepsGx, ss')
 
diff --git a/Math/Algebra/Group/SchreierSims.hs b/Math/Algebra/Group/SchreierSims.hs
--- a/Math/Algebra/Group/SchreierSims.hs
+++ b/Math/Algebra/Group/SchreierSims.hs
@@ -8,6 +8,7 @@
 import qualified Data.Map as M
 import Math.Algebra.Group.PermutationGroup hiding (elts, order, gens, isMember, isSubgp, isNormal, reduceGens, normalClosure, commutatorGp, derivedSubgp)
 import Math.Common.ListSet (toListSet)
+import Math.Core.Utils hiding (elts)
 
 
 -- COSET REPRESENTATIVES FOR STABILISER OF A POINT
diff --git a/Math/Algebra/LinearAlgebra.hs b/Math/Algebra/LinearAlgebra.hs
--- a/Math/Algebra/LinearAlgebra.hs
+++ b/Math/Algebra/LinearAlgebra.hs
@@ -1,6 +1,4 @@
--- Copyright (c) David Amos, 2008. All rights reserved.
-
-{-# OPTIONS_GHC -fglasgow-exts #-}
+-- Copyright (c) 2008-2012, David Amos. All rights reserved.
 
 -- |A module providing elementary operations involving scalars, vectors, and matrices
 -- over a ring or field. Vectors are represented as [a], matrices as [[a]].
@@ -165,7 +163,6 @@
 inverse2 ((1:r):rs) = inverse2' r rs : inverse2 rs where
     inverse2' xs [] = xs
     inverse2' (x:xs) ((1:r):rs) = inverse2' (xs <-> x *> r) rs
--- This is basically reduced row echelon form
 
 xs ! i = xs !! (i-1) -- ie, a 1-based list lookup instead of 0-based
 
@@ -187,6 +184,21 @@
     reduceStep rs@((0:_):_) = zipWith (:) (map head rs) (reduceStep $ map tail rs)
     reduceStep rs = rs
 
+-- Given a matrix m and (column) vector b, either find (column vector) x such that m x == b,
+-- or indicate that there is none
+solveLinearSystem m b =
+    let augmented = zipWith (\r x -> r ++ [x]) m b -- augmented matrix
+        trisystem = inverse1 augmented -- upper triangular form
+        solution = reverse $ solveTriSystem $ reverse $ map reverse trisystem
+    in if length solution == length b then Just solution else Nothing
+    where solveTriSystem ([v,c]:rs) =
+              let x = v/c -- the first row tells us that cx == v
+                  rs' = map (\(v':c':r) -> (v'-c'*x):r) rs
+              in x : solveTriSystem rs'
+          solveTriSystem [] = []
+          solveTriSystem _ = [] -- abnormal termination - m wasn't invertible
+
+
 isZero v = all (==0) v
 
 -- inSpanRE m v returns whether the vector v is in the span of the matrix m, where m is required to be in row echelon form
@@ -213,7 +225,7 @@
     findLeadingCols i (c@(0:_):cs) = findLeadingCols (i+1) cs
     findLeadingCols _ _ = []
 
-m ^- n = recip m ^ n
+-- m ^- n = recip m ^ n
 
 -- t (M m) = M (L.transpose m)
 
diff --git a/Math/Algebras/Commutative.hs b/Math/Algebras/Commutative.hs
--- a/Math/Algebras/Commutative.hs
+++ b/Math/Algebras/Commutative.hs
@@ -47,7 +47,7 @@
 -}
 
 -- This is the monoid algebra for commutative monomials (which are the free commutative monoid)
-instance (Num k, Ord v) => Algebra k (GlexMonomial v) where
+instance (Eq k, Num k, Ord v) => Algebra k (GlexMonomial v) where
     unit x = x *> return munit
         where munit = Glex 0 []
     mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy
@@ -55,7 +55,7 @@
 
 
 -- GlexPoly can be given the set coalgebra structure, which is compatible with the monoid algebra structure
-instance Num k => Coalgebra k (GlexMonomial v) where
+instance (Eq k, Num k) => Coalgebra k (GlexMonomial v) where
     counit = unwrap . nf . fmap (\m -> () )  -- trace
     -- counit (V ts) = sum [x | (m,x) <- ts]  -- trace
     comult = fmap (\m -> (m,m) )             -- diagonal
@@ -78,7 +78,7 @@
 -- |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) =>
+bind :: (Monomial m, Eq k, 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 *> product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] 
 -- flipbind f = linear (\m -> product [f x ^ i | (x,i) <- powers m])
@@ -120,7 +120,7 @@
 infixl 7 %%
 
 -- |(%%) reduces a polynomial with respect to a list of polynomials.
-(%%) :: (Fractional k, Ord b, Show b, Algebra k b, DivisionBasis b)
+(%%) :: (Eq k, 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
 
diff --git a/Math/Algebras/GroupAlgebra.hs b/Math/Algebras/GroupAlgebra.hs
--- a/Math/Algebras/GroupAlgebra.hs
+++ b/Math/Algebras/GroupAlgebra.hs
@@ -1,18 +1,35 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2012, David Amos. All rights reserved.
 
-{-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances #-}
+{-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances #-}
+-- ScopedTypeVariables
 
+-- |A module for doing arithmetic in the group algebra.
+--
+-- Group elements are represented as permutations of the integers, and are entered and displayed
+-- using a Haskell-friendly version of cycle notation. For example, the permutation (1 2 3)(4 5)
+-- would be entered as @p [[1,2,3],[4,5]]@, and displayed as [[1,2,3],[4,5]].
+--
+-- Given a field K and group G, the group algebra KG is the free K-vector space over the elements of G.
+-- Elements of the group algebra consists of arbitrary K-linear combinations of elements of G.
+-- For example, @p [[1,2,3]] + 2 * p [[1,2],[3,4]]@
 module Math.Algebras.GroupAlgebra where
 
+import Math.Core.Field
+import Math.Core.Utils
+
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
 import Math.Algebras.Structures
 
-import Math.Algebra.Group.PermutationGroup hiding (action)
+import Math.Algebra.Group.PermutationGroup hiding (p, action)
+import qualified Math.Algebra.Group.PermutationGroup as P
 
-import Math.Algebra.Field.Base
+import Math.Algebra.LinearAlgebra hiding (inverse, (*>) )
 
+import Math.CommutativeAlgebra.Polynomial
+import Math.CommutativeAlgebra.GroebnerBasis
 
+
 instance Mon (Permutation Int) where
     munit = 1
     mmult = (*)
@@ -20,7 +37,7 @@
 type GroupAlgebra k = Vect k (Permutation Int)
 
 -- Monoid Algebra instance
-instance Num k => Algebra k (Permutation Int) where
+instance (Eq k, Num k) => Algebra k (Permutation Int) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
     mult = nf . fmap (\(a,b) -> a `mmult` b)
@@ -28,27 +45,76 @@
 -- Set Coalgebra instance
 -- instance SetCoalgebra (Permutation Int) where {}
 
-instance Num k => Coalgebra k (Permutation Int) where
+instance (Eq k, Num k) => Coalgebra k (Permutation Int) where
     counit (V ts) = sum [x | (m,x) <- ts] -- trace
     comult = fmap (\m -> (m,m)) -- diagonal
 
-instance Num k => Bialgebra k (Permutation Int) where {}
+instance (Eq k, 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
+instance (Eq k, 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
+-- |Construct a permutation, as an element of the group algebra, from a list of cycles.
+-- For example, @p [[1,2],[3,4,5]]@ constructs the permutation (1 2)(3 4 5), which is displayed
+-- as [[1,2],[3,4,5]].
+p :: [[Int]] -> GroupAlgebra Q
+p cs = return $ P.p cs
 
 
-instance Num k => Module k (Permutation Int) Int where
+instance (Eq k, Num k) => Module k (Permutation Int) Int where
     action = nf . fmap (\(g,x) -> x .^ g)
 
 -- use *. instead
 -- r *> m = action (r `te` m)
 
+newtype X a = X a deriving (Eq,Ord,Show)
 
+-- Find the inverse of a group algebra element using Groebner basis techniques
+-- This is overkill, but it was what I had to hand at first
+inv x@(V ts) =
+    let gs = P.elts $ map fst $ terms x -- all elements in the group generated by the terms
+        cs = map (glexvar . X) gs
+        x' = V $ map (\(g,c) -> (g, unit c)) ts
+        one = x' * (V $ zip gs cs)
+        oneEquations = (coeff 1 one - 1) : [coeff g one - 0 | g <- tail gs]
+        zeroEquations = [coeff g one - 0 | g <- gs]
+        solution = gb oneEquations
+    in if solution == [1]
+       then Left (gb zeroEquations) -- it's a zero divisor
+       else Right solution
+       -- sum [-c *> p g | V [ (Glex (M 1 [(X g, 1)]), 1), (Glex (M 0 []), c) ] <- solution]
+       -- should extract the solution into a group algebra element, but having trouble getting types right
 
+-- The following code can be made to work over an arbitrary field by uncommenting the commented code
+-- However, we should then probably also change the signature of p to p :: Fractional k => [[Int]] -> GroupAlgebra k
+-- instance Fractional k => HasInverses (GroupAlgebra k) where
 
+-- |Note that the inverse of a group algebra element can only be efficiently calculated
+-- if the group generated by the non-zero terms is very small (eg \<100 elements).
+instance HasInverses (GroupAlgebra Q) where
+    inverse x@(V ts) =
+        let gs = P.elts $ map fst $ terms x -- all elements in the group generated by the terms
+            -- cs = map (var . X) gs :: [Vect k (Glex (X (Permutation Int)))]
+            cs = map (glexvar . X) gs
+            x' = V $ map (\(g,c) -> (g, unit c)) ts
+            one = x' * (V $ zip gs cs)
+            m = [ [coeff (mvar (X j)) c | j <- gs] | i <- gs, let c = coeff i one]
+            b = 1 : replicate (length gs - 1) 0
+        in case solveLinearSystem m b of
+            Just v -> nf $ V $ zip gs v
+            Nothing -> error "GroupAlgebra.inverse: not invertible"
+
+maybeInverse x@(V ts) =
+    let gs = P.elts $ map fst $ terms x -- all elements in the group generated by the terms
+        cs = map (glexvar . X) gs
+        x' = V $ map (\(g,c) -> (g, unit c)) ts
+        one = x' * (V $ zip gs cs)
+        m = [ [coeff (mvar (X j)) c | j <- gs] | i <- gs, let c = coeff i one]
+        b = 1 : replicate (length gs - 1) 0
+    in fmap (\v -> nf $ V $ zip gs v) (solveLinearSystem m b)
+{-
+    in case solveLinearSystem m b of
+        Just v -> Just $ nf $ V $ zip gs v
+        Nothing -> Nothing
+-}
diff --git a/Math/Algebras/LaurentPoly.hs b/Math/Algebras/LaurentPoly.hs
--- a/Math/Algebras/LaurentPoly.hs
+++ b/Math/Algebras/LaurentPoly.hs
@@ -31,7 +31,7 @@
     munit = LM 0 []
     mmult (LM si xis) (LM sj yjs) = LM (si+sj) $ addmerge xis yjs
 
-instance Num k => Algebra k LaurentMonomial where
+instance (Eq k, Num k) => Algebra k LaurentMonomial where
     unit 0 = zero -- V []
     unit x = V [(munit,x)] 
     mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts)
@@ -49,7 +49,7 @@
 
 lvar v = V [(LM 1 [(v,1)], 1)] :: LaurentPoly Q
 
-instance Fractional k => Fractional (LaurentPoly k) where
+instance (Eq k, 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"
 
diff --git a/Math/Algebras/Matrix.hs b/Math/Algebras/Matrix.hs
--- a/Math/Algebras/Matrix.hs
+++ b/Math/Algebras/Matrix.hs
@@ -22,7 +22,7 @@
 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
+instance (Eq k, Num k) => Algebra k Mat2 where
     unit x = x *> V [(E2 i i, 1) | i <- [1..2] ]
     mult = linear mult' where
         mult' (E2 i j, E2 k l) = delta j k *> return (E2 i l)
@@ -33,7 +33,7 @@
 -- 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
+instance (Eq k, Num k) => Module k Mat2 EBasis where
     -- action ax = nf $ ax >>= action' where
     action = linear action' where
         action' (E2 i j, E k) = delta j k `smultL` return (E i)
@@ -55,7 +55,7 @@
 -- E2' i j represents the dual basis element corresponding to E i j
 
 -- Kassel p42
-instance Num k => Coalgebra k Mat2' where
+instance (Eq k, 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 (E2' i k, E2' k j) | k <- [1..2]])
@@ -73,7 +73,7 @@
 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
+instance (Eq k, 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 (\((E3 i j, E3 k l), x) -> (E3 i l, delta j k * x)) ts
diff --git a/Math/Algebras/NonCommutative.hs b/Math/Algebras/NonCommutative.hs
--- a/Math/Algebras/NonCommutative.hs
+++ b/Math/Algebras/NonCommutative.hs
@@ -29,7 +29,7 @@
     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
+instance (Eq k, Num k, Ord v) => Algebra k (NonComMonomial v) where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
     mult = nf . fmap (\(a,b) -> a `mmult` b)
diff --git a/Math/Algebras/Octonions.hs b/Math/Algebras/Octonions.hs
--- a/Math/Algebras/Octonions.hs
+++ b/Math/Algebras/Octonions.hs
@@ -5,7 +5,7 @@
 -- |A module defining the (non-associative) algebra of octonions over an arbitrary field.
 --
 -- The octonions are the algebra defined by the basis {1,i0,i1,i2,i3,i4,i5,i6},
--- where each i_n^2 = -1, and i_n+1*i_n+2 = i_n+4 (where the indices are modulo 7).
+-- where each i_n * i_n = -1, and i_n+1 * i_n+2 = i_n+4 (where the indices are modulo 7).
 module Math.Algebras.Octonions where
 
 import Math.Core.Field
@@ -43,7 +43,7 @@
 i_ :: Num k => Int -> Octonion k
 i_ n = return (O n)
 
-instance (Num k) => Algebra k OBasis where
+instance (Eq k, Num k) => Algebra k OBasis where
     unit x = x *> return (O (-1))
     mult = linear m where
         m (O (-1), O n) = return (O n)
@@ -57,7 +57,7 @@
                        5 -> -1 *> i_ ((a+4) `mod` 7) -- i_n+4 * i_n+2 == -i_n+1
                        6 -> -1 *> i_ ((a+2) `mod` 7) -- i_n+2 * i_n+1 == -i_n+4
 
-instance Num k => HasConjugation k OBasis where
+instance (Eq k, Num k) => HasConjugation k OBasis where
     conj = (>>= conj') where
         conj' (O n) = (if n == -1 then 1 else -1) *> return (O n)
     -- ie conj = linear conj', but avoiding unnecessary nf call
diff --git a/Math/Algebras/Quaternions.hs b/Math/Algebras/Quaternions.hs
--- a/Math/Algebras/Quaternions.hs
+++ b/Math/Algebras/Quaternions.hs
@@ -27,7 +27,7 @@
     show J = "j"
     show K = "k"
 
-instance (Num k) => Algebra k HBasis where
+instance (Eq k, Num k) => Algebra k HBasis where
     unit x = x *> return One
     mult = linear mult'
          where mult' (One,b) = return b
@@ -69,7 +69,7 @@
 
 -- |If an algebra has a conjugation operation, then it has multiplicative inverses,
 -- via 1/x = conj x / sqnorm x
-instance (Fractional k, Ord a, Show a, HasConjugation k a) => Fractional (Vect k a) where
+instance (Eq k, Fractional k, Ord a, Show a, HasConjugation k a) => Fractional (Vect k a) where
     recip 0 = error "recip 0"
     recip x = (1 / sqnorm x) *> conj x
     fromRational q = fromRational q *> 1
@@ -79,10 +79,10 @@
 scalarPart = coeff One
 
 -- |The vector part of the quaternion w+xi+yj+zk is xi+yj+zk. Also called the pure part.
-vectorPart :: (Num k) => Quaternion k -> Quaternion k
+vectorPart :: (Eq k, Num k) => Quaternion k -> Quaternion k
 vectorPart q = q - scalarPart q *> 1
 
-instance Num k => HasConjugation k HBasis where
+instance (Eq k, Num k) => HasConjugation k HBasis where
     conj = (>>= conj') where
         conj' One = return One
         conj' imag = -1 *> return imag
@@ -126,7 +126,7 @@
 -- This shows that the multiplicative group of unit quaternions is isomorphic to Spin3, the double cover of SO3.
 --
 -- @reprSO3 q@ returns the 3*3 matrix representing this map.
-reprSO3 :: (Fractional k) => Quaternion k -> [[k]]
+reprSO3 :: (Eq k, Fractional k) => Quaternion k -> [[k]]
 reprSO3 q = reprSO3' q `asMatrix` [i,j,k]
 -- It's clear from the definition that repr3' q leaves scalars invariant
 
@@ -145,7 +145,7 @@
 -- is isomorphic to Spin4, the double cover of SO4.
 --
 -- @reprSO4 (l,r)@ returns the 4*4 matrix representing this map.
-reprSO4 :: (Fractional k) => (Quaternion k, Quaternion k) -> [[k]]
+reprSO4 :: (Eq k, Fractional k) => (Quaternion k, Quaternion k) -> [[k]]
 reprSO4 (l,r) = reprSO4' (l,r) `asMatrix` [1,i,j,k]
 -- could consider checking that l,r are unit length - except that this is hard to achieve working over Q
 
@@ -164,7 +164,7 @@
 
 -- Coalgebra structure on the dual vector space to the quaternions
 -- The comult is the transpose of mult
-instance Num k => Coalgebra k (Dual HBasis) where
+instance (Eq k, Num k) => Coalgebra k (Dual HBasis) where
     counit = unwrap . linear counit'
         where counit' (Dual One) = return ()
               counit' _ = zero
diff --git a/Math/Algebras/Structures.hs b/Math/Algebras/Structures.hs
--- a/Math/Algebras/Structures.hs
+++ b/Math/Algebras/Structures.hs
@@ -43,7 +43,7 @@
     antipode :: Vect k b -> Vect k b
 
 
-instance (Num k, Eq b, Ord b, Show b, Algebra k b) => Num (Vect k b) where
+instance (Eq k, Num k, Eq b, Ord b, Show b, Algebra k b) => Num (Vect k b) where
     x+y = x <+> y
     negate x = neg x
     -- negate (V ts) = V $ map (\(b,x) -> (b, negate x)) ts
@@ -66,7 +66,7 @@
 -}
 
 
-instance Num k => Algebra k () where
+instance (Eq k, Num k) => Algebra k () where
     unit = wrap
     -- unit 0 = zero -- V []
     -- unit x = V [( (),x)]
@@ -74,7 +74,7 @@
     -- mult (V [( ((),()), x)]) = V [( (),x)]
     -- mult (V []) = zerov
 
-instance Num k => Coalgebra k () where
+instance (Eq k, Num k) => Coalgebra k () where
     counit = unwrap
     -- counit (V []) = 0
     -- counit (V [( (),x)]) = x
@@ -82,10 +82,10 @@
     -- comult (V [( (),x)]) = V [( ((),()), x)]
     -- comult (V []) = zerov
 
-unit' :: (Num k, Algebra k b) => Trivial k -> Vect k b
+unit' :: (Eq k, 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' :: (Eq k, 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
@@ -94,7 +94,7 @@
 -- Kassel p4
 -- |The direct sum of k-algebras can itself be given the structure of a k-algebra.
 -- This is the product object in the category of k-algebras.
-instance (Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (DSum a b) where
+instance (Eq k, Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (DSum a b) where
     unit k = i1 (unit k) <+> i2 (unit k)
     -- unit == (i1 . unit) <<+>> (i2 . unit)
     mult = linear mult'
@@ -108,7 +108,7 @@
 
 -- |The direct sum of k-coalgebras can itself be given the structure of a k-coalgebra.
 -- This is the coproduct object in the category of k-coalgebras.
-instance (Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (DSum a b) where
+instance (Eq k, Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (DSum a b) where
     counit = unwrap . linear counit'
         where counit' (Left a) = (wrap . counit) (return a)
               counit' (Right b) = (wrap . counit) (return b)
@@ -123,7 +123,7 @@
 
 -- Kassel p32
 -- |The tensor product of k-algebras can itself be given the structure of a k-algebra
-instance (Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (Tensor a b) where
+instance (Eq k, Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (Tensor a b) where
     -- unit 0 = V []
     unit x = x *> (unit 1 `te` unit 1)
     mult = linear m where
@@ -131,7 +131,7 @@
 
 -- Kassel p42
 -- |The tensor product of k-coalgebras can itself be given the structure of a k-coalgebra
-instance (Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (Tensor a b) where
+instance (Eq k, 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))
@@ -139,20 +139,20 @@
 
 
 -- The set coalgebra - can be defined on any set
-instance Num k => Coalgebra k EBasis where
+instance (Eq k, Num k) => Coalgebra k EBasis where
     counit (V ts) = sum [x | (ei,x) <- ts]  -- trace
     comult = fmap ( \ei -> (ei,ei) )        -- diagonal
 
 newtype SetCoalgebra b = SC b deriving (Eq,Ord,Show)
 
-instance Num k => Coalgebra k (SetCoalgebra b) where
+instance (Eq k, Num k) => Coalgebra k (SetCoalgebra b) where
     counit (V ts) = sum [x | (m,x) <- ts]  -- trace
     comult = fmap ( \m -> (m,m) )          -- diagonal
 
 
 newtype MonoidCoalgebra m = MC m deriving (Eq,Ord,Show)
 
-instance (Num k, Ord m, Mon m) => Coalgebra k (MonoidCoalgebra m) where
+instance (Eq k, 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 (m,m) else return (m, MC munit) <+> return (MC munit, m)
@@ -184,13 +184,13 @@
 
 -- Kassel p57-8
 
-instance (Num k, Ord a, Ord u, Ord v, Algebra k a, Module k a u, Module k a v)
+instance (Eq k, 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' ((a,a'), (u,v)) = (action $ return (a,u)) `te` (action $ return (a',v))
 
-instance (Num k, Ord a, Ord u, Ord v, Bialgebra k a, Module k a u, Module k a v)
+instance (Eq k, 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'
@@ -200,7 +200,7 @@
 -- 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)
+instance (Eq k, 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 ( \((h,m), (h',n)) -> ((h,h'), (m,n)) ) x
diff --git a/Math/Algebras/TensorAlgebra.hs b/Math/Algebras/TensorAlgebra.hs
--- a/Math/Algebras/TensorAlgebra.hs
+++ b/Math/Algebras/TensorAlgebra.hs
@@ -34,7 +34,7 @@
     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
+instance (Eq k, Num k, Ord a) => Algebra k (TensorAlgebra a) where
     unit x = x *> return munit
     mult = nf . fmap (\(a,b) -> a `mmult` b)
 
@@ -50,7 +50,7 @@
 -- The Num k context is not strictly necessary
 
 -- |Inject an element of the set\/type A\/a into the tensor algebra T(A) = Vect k (TensorAlgebra a).
-injectTA' :: Num k => a -> Vect k (TensorAlgebra a)
+injectTA' :: (Eq k, Num k) => a -> Vect k (TensorAlgebra a)
 injectTA' = injectTA . return
 -- injectTA' a = return (TA 1 [a])
 
@@ -59,7 +59,7 @@
 -- where T(A) is the tensor algebra Vect k (TensorAlgebra a).
 -- f' will agree with f on A itself (considered as a subspace of T(A)).
 -- In other words, f = f' . injectTA
-liftTA :: (Num k, Ord b, Show b, Algebra k b) =>
+liftTA :: (Eq k, Num k, Ord b, Show b, Algebra k b) =>
      (Vect k a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k b
 liftTA f = linear (\(TA _ xs) -> product [f (return x) | x <- xs])
 -- The Show b constraint is required because we use product (and Num requires Show)!!
@@ -67,7 +67,7 @@
 -- |Given a set\/type A\/a, and a vector space B = Vect k b, where B is also an algebra,
 -- lift a function f: A -> B to an algebra morphism f': T(A) -> B.
 -- f' will agree with f on A itself. In other words, f = f' . injectTA'
-liftTA' :: (Num k, Ord b, Show b, Algebra k b) =>
+liftTA' :: (Eq k, Num k, Ord b, Show b, Algebra k b) =>
      (a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k b
 liftTA' = liftTA . linear
 -- liftTA' f = linear (\(TA _ xs) -> product [f x | x <- xs])
@@ -77,7 +77,7 @@
 -- |Tensor algebra is a functor from k-Vect to k-Alg.
 -- The action on objects is Vect k a -> Vect k (TensorAlgebra a).
 -- The action on arrows is f -> fmapTA f.
-fmapTA :: (Num k, Ord b, Show b) =>
+fmapTA :: (Eq k, Num k, Ord b, Show b) =>
     (Vect k a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k (TensorAlgebra b)
 fmapTA f = liftTA (injectTA . f)
 -- fmapTA f = linear (\(TA _ xs) -> product [injectTA (f (return x)) | x <- xs])
@@ -86,18 +86,18 @@
 -- we obtain a functor Set -> k-Alg, the free algebra functor.
 -- The action on objects is a -> Vect k (TensorAlgebra a).
 -- The action on arrows is f -> fmapTA' f.
-fmapTA' :: (Num k, Ord b, Show b) =>
+fmapTA' :: (Eq k, Num k, Ord b, Show b) =>
     (a -> b) -> Vect k (TensorAlgebra a) -> Vect k (TensorAlgebra b)
 fmapTA' = fmapTA . fmap
 -- fmapTA' f = liftTA' (injectTA' . f)
 -- fmapTA' f = linear (\(TA _ xs) -> product [injectTA' (f x) | x <- xs])
 
 
-bindTA :: (Num k, Ord b, Show b) =>
+bindTA :: (Eq k, Num k, Ord b, Show b) =>
     Vect k (TensorAlgebra a) -> (Vect k a -> Vect k (TensorAlgebra b)) -> Vect k (TensorAlgebra b)
 bindTA = flip liftTA
 
-bindTA' :: (Num k, Ord b, Show b) =>
+bindTA' :: (Eq k, Num k, Ord b, Show b) =>
     Vect k (TensorAlgebra a) -> (a -> Vect k (TensorAlgebra b)) -> Vect k (TensorAlgebra b)
 bindTA' = flip liftTA'
 -- Another way to think about this is variable substitution
@@ -121,14 +121,14 @@
     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
+instance (Eq k, Num k, Ord a) => Algebra k (SymmetricAlgebra a) where
     unit x = x *> return munit
     mult = nf . fmap (\(a,b) -> a `mmult` b)
 
 -- |Algebra morphism from tensor algebra to symmetric algebra.
 -- The kernel of the morphism is the ideal generated by all
 -- differences of products u&#x2297;v - v&#x2297;u.
-toSym :: (Num k, Ord a) =>
+toSym :: (Eq k, Num k, Ord a) =>
      Vect k (TensorAlgebra a) -> Vect k (SymmetricAlgebra a)
 toSym = linear toSym'
     where toSym' (TA i xs) = return $ Sym i (L.sort xs) 
@@ -147,31 +147,31 @@
 injectSym' = injectSym . return
 -- injectSym' a = return (Sym 1 [a])
 
-liftSym :: (Num k, Ord b, Show b, Algebra k b) =>
+liftSym :: (Eq k, Num k, Ord b, Show b, Algebra k b) =>
      (Vect k a -> Vect k b) -> Vect k (SymmetricAlgebra a) -> Vect k b
 liftSym f = linear (\(Sym _ xs) -> product [f (return x) | x <- xs])
 
-liftSym' :: (Num k, Ord b, Show b, Algebra k b) =>
+liftSym' :: (Eq k, Num k, Ord b, Show b, Algebra k b) =>
      (a -> Vect k b) -> Vect k (SymmetricAlgebra a) -> Vect k b
 liftSym' = liftSym . linear
 -- liftSym' f = linear (\(Sym _ xs) -> product [f x | x <- xs])
 
-fmapSym :: (Num k, Ord b, Show b) =>
+fmapSym :: (Eq k, Num k, Ord b, Show b) =>
     (Vect k a -> Vect k b) -> Vect k (SymmetricAlgebra a) -> Vect k (SymmetricAlgebra b)
 fmapSym f = liftSym (injectSym . f)
 -- fmapSym f = linear (\(Sym _ xs) -> product [injectSym (f (return x)) | x <- xs])
 
-fmapSym' :: (Num k, Ord b, Show b) =>
+fmapSym' :: (Eq k, Num k, Ord b, Show b) =>
     (a -> b) -> Vect k (SymmetricAlgebra a) -> Vect k (SymmetricAlgebra b)
 fmapSym' = fmapSym . fmap
 -- fmapSym' f = liftSym' (injectSym' . f)
 -- fmapSym' f = linear (\(Sym _ xs) -> product [injectSym' (f x) | x <- xs])
 
-bindSym :: (Num k, Ord b, Show b) =>
+bindSym :: (Eq k, Num k, Ord b, Show b) =>
     Vect k (SymmetricAlgebra a) -> (Vect k a -> Vect k (SymmetricAlgebra b)) -> Vect k (SymmetricAlgebra b)
 bindSym = flip liftSym
 
-bindSym' :: (Num k, Ord b, Show b) =>
+bindSym' :: (Eq k, Num k, Ord b, Show b) =>
     Vect k (SymmetricAlgebra a) -> (a -> Vect k (SymmetricAlgebra b)) -> Vect k (SymmetricAlgebra b)
 bindSym' = flip liftSym'
 -- Another way to think about this is variable substitution
@@ -190,7 +190,7 @@
     show (Ext _ xs) = filter (/= '"') $ concat $ L.intersperse "^" $ map show xs
 
 
-instance (Num k, Ord a) => Algebra k (ExteriorAlgebra a) where
+instance (Eq k, Num k, Ord a) => Algebra k (ExteriorAlgebra a) where
     unit x = x *> return (Ext 0 [])
     mult xy = nf $ xy >>= (\(Ext i xs, Ext j ys) -> signedMerge 1 (0,[]) (i,xs) (j,ys))
         where signedMerge s (k,zs) (i,x:xs) (j,y:ys) =
@@ -206,7 +206,7 @@
 -- |Algebra morphism from tensor algebra to exterior algebra.
 -- The kernel of the morphism is the ideal generated by all
 -- self-products u&#x2297;u and sums of products u&#x2297;v + v&#x2297;u
-toExt :: (Num k, Ord a) =>
+toExt :: (Eq k, Num k, Ord a) =>
      Vect k (TensorAlgebra a) -> Vect k (ExteriorAlgebra a)
 toExt = linear toExt'
     where toExt' (TA i xs) = let (sign,xs') = signedSort 1 True [] xs
@@ -230,31 +230,31 @@
 injectExt' = injectExt . return
 -- injectExt' a = return (Ext 1 [a])
 
-liftExt :: (Num k, Ord b, Show b, Algebra k b) =>
+liftExt :: (Eq k, Num k, Ord b, Show b, Algebra k b) =>
      (Vect k a -> Vect k b) -> Vect k (ExteriorAlgebra a) -> Vect k b
 liftExt f = linear (\(Ext _ xs) -> product [f (return x) | x <- xs])
 
-liftExt' :: (Num k, Ord b, Show b, Algebra k b) =>
+liftExt' :: (Eq k, Num k, Ord b, Show b, Algebra k b) =>
      (a -> Vect k b) -> Vect k (ExteriorAlgebra a) -> Vect k b
 liftExt' = liftExt . linear
 -- liftExt' f = linear (\(Ext _ xs) -> product [f x | x <- xs])
 
-fmapExt :: (Num k, Ord b, Show b) =>
+fmapExt :: (Eq k, Num k, Ord b, Show b) =>
     (Vect k a -> Vect k b) -> Vect k (ExteriorAlgebra a) -> Vect k (ExteriorAlgebra b)
 fmapExt f = liftExt (injectExt . f)
 -- fmapExt f = linear (\(Ext _ xs) -> product [injectExt (f (return x)) | x <- xs])
 
-fmapExt' :: (Num k, Ord b, Show b) =>
+fmapExt' :: (Eq k, Num k, Ord b, Show b) =>
     (a -> b) -> Vect k (ExteriorAlgebra a) -> Vect k (ExteriorAlgebra b)
 fmapExt' = fmapExt . fmap
 -- fmapExt' f = liftExt' (injectExt' . f)
 -- fmapExt' f = linear (\(Ext _ xs) -> product [injectExt' (f x) | x <- xs])
 
-bindExt :: (Num k, Ord b, Show b) =>
+bindExt :: (Eq k, Num k, Ord b, Show b) =>
     Vect k (ExteriorAlgebra a) -> (Vect k a -> Vect k (ExteriorAlgebra b)) -> Vect k (ExteriorAlgebra b)
 bindExt = flip liftExt
 
-bindExt' :: (Num k, Ord b, Show b) =>
+bindExt' :: (Eq k, Num k, Ord b, Show b) =>
     Vect k (ExteriorAlgebra a) -> (a -> Vect k (ExteriorAlgebra b)) -> Vect k (ExteriorAlgebra b)
 bindExt' = flip liftExt'
 -- Another way to think about this is variable substitution
@@ -265,7 +265,7 @@
 -- Kassel p67
 data TensorCoalgebra c = TC Int [c] deriving (Eq,Ord,Show)
 
-instance (Num k, Ord c) => Coalgebra k (TensorCoalgebra c) where
+instance (Eq k, Num k, Ord c) => Coalgebra k (TensorCoalgebra c) where
     counit = unwrap . linear counit'
         where counit' (TC 0 []) = return () -- 1
               counit' _ = zerov
@@ -278,14 +278,14 @@
 -- coliftTC f is a coalgebra morphism, and f == projectTC . coliftTC f
 
 -- projection onto the underlying vector space
-projectTC :: (Num k, Ord b) => Vect k (TensorCoalgebra b) -> Vect k b
+projectTC :: (Eq k, Num k, Ord b) => Vect k (TensorCoalgebra b) -> Vect k b
 projectTC = linear projectTC' where projectTC' (TC 1 [b]) = return b; projectTC' _ = zerov 
 -- projectTC t = V [(b,c) | (TC 1 [b], c) <- terms t]
 
 
 -- lift a vector space morphism C -> D to a coalgebra morphism C -> T'(D)
 -- this function returns an approximation, valid only up to second order terms
-coliftTC :: (Num k, Coalgebra k c, Ord d) =>
+coliftTC :: (Eq k, Num k, Coalgebra k c, Ord d) =>
      (Vect k c -> Vect k d) -> Vect k c -> Vect k (TensorCoalgebra d)
 coliftTC f = sumf [coliftTC' i f | i <- [0..2] ]
 
@@ -299,7 +299,7 @@
           fn' c = fmap (\(TC 1 [x], TC _ xs) -> TC n (x:xs)) $ ( (f1' `tf` fn1') . comult) (return c)
 
 
-cobindTC :: (Num k, Ord c, Ord d) =>
+cobindTC :: (Eq k, Num k, Ord c, Ord d) =>
      (Vect k (TensorCoalgebra c) -> Vect k d) -> Vect k (TensorCoalgebra c) -> Vect k (TensorCoalgebra d)
 cobindTC = coliftTC
 
diff --git a/Math/Algebras/TensorProduct.hs b/Math/Algebras/TensorProduct.hs
--- a/Math/Algebras/TensorProduct.hs
+++ b/Math/Algebras/TensorProduct.hs
@@ -27,7 +27,7 @@
 
 -- |The coproduct of two linear functions (with the same target).
 -- Satisfies the universal property that f == coprodf f g . i1 and g == coprodf f g . i2
-coprodf :: (Num k, Ord t) =>
+coprodf :: (Eq k, Num k, Ord t) =>
     (Vect k a -> Vect k t) -> (Vect k b -> Vect k t) -> Vect k (DSum a b) -> Vect k t
 coprodf f g = linear fg' where
     fg' (Left a) = f (return a)
@@ -35,33 +35,33 @@
 
 
 -- |Projection onto left summand from direct sum
-p1 :: (Num k, Ord a) => Vect k (DSum a b) -> Vect k a
+p1 :: (Eq k, Num k, Ord a) => Vect k (DSum a b) -> Vect k a
 p1 = linear p1' where
     p1' (Left a) = return a
     p1' (Right b) = zero
 
 -- |Projection onto right summand from direct sum
-p2 :: (Num k, Ord b) => Vect k (DSum a b) -> Vect k b
+p2 :: (Eq k, Num k, Ord b) => Vect k (DSum a b) -> Vect k b
 p2 = linear p2' where
     p2' (Left a) = zero
     p2' (Right b) = return b
 
 -- |The product of two linear functions (with the same source).
 -- Satisfies the universal property that f == p1 . prodf f g and g == p2 . prodf f g
-prodf :: (Num k, Ord a, Ord b) =>
+prodf :: (Eq k, Num k, Ord a, Ord b) =>
     (Vect k s -> Vect k a) -> (Vect k s -> Vect k b) -> Vect k s -> Vect k (DSum a b)
 prodf f g = linear fg' where
     fg' b = fmap Left (f $ return b) <+> fmap Right (g $ return b)
 
 
 -- |The direct sum of two vector space elements
-dsume :: (Num k, Ord a, Ord b) => Vect k a -> Vect k b -> Vect k (DSum a b)
+dsume :: (Eq k, Num k, Ord a, Ord b) => Vect k a -> Vect k b -> Vect k (DSum a b)
 -- dsume x y = fmap Left x <+> fmap Right y
 dsume x y = i1 x <+> i2 y
 
 -- |The direct sum of two linear functions.
 -- Satisfies the universal property that f == p1 . dsumf f g . i1 and g == p2 . dsumf f g . i2
-dsumf :: (Num k, Ord a, Ord b, Ord a', Ord b') => 
+dsumf :: (Eq k, Num k, Ord a, Ord b, Ord a', Ord b') => 
     (Vect k a -> Vect k a') -> (Vect k b -> Vect k b') -> Vect k (DSum a b) -> Vect k (DSum a' b')
 dsumf f g ab = (i1 . f . p1) ab <+> (i2 . g . p2) ab
 
@@ -80,7 +80,7 @@
 
 -- Implicit assumption - f and g are linear
 -- |The tensor product of two linear functions
-tf :: (Num k, Ord a', Ord b') => (Vect k a -> Vect k a') -> (Vect k b -> Vect k b')
+tf :: (Eq k, 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 [x *> te (f $ return a) (g $ return b) | ((a,b), x) <- ts]
     where sum = foldl add zero -- (V [])
@@ -107,16 +107,16 @@
 unitOutR :: Vect k (Tensor a ()) -> Vect k a
 unitOutR = fmap ( \(a,()) -> a )
 
-twist :: (Num k, Ord a, Ord b) => Vect k (Tensor a b) -> Vect k (Tensor b a)
+twist :: (Eq k, Num k, Ord a, Ord b) => Vect k (Tensor a b) -> Vect k (Tensor b a)
 twist v = nf $ fmap ( \(a,b) -> (b,a) ) v
 -- note the nf call, as f is not order-preserving
 
 
-distrL :: (Num k, Ord a, Ord b, Ord c)
+distrL :: (Eq k, Num k, Ord a, Ord b, Ord c)
     => Vect k (Tensor a (DSum b c)) -> Vect k (DSum (Tensor a b) (Tensor a c))
 distrL v = nf $ fmap (\(a,bc) -> case bc of Left b -> Left (a,b); Right c -> Right (a,c)) v
 
-undistrL :: (Num k, Ord a, Ord b, Ord c)
+undistrL :: (Eq k, Num k, Ord a, Ord b, Ord c)
     => Vect k (DSum (Tensor a b) (Tensor a c)) -> Vect k (Tensor a (DSum b c))
 undistrL v = nf $ fmap ( \abc -> case abc of Left (a,b) -> (a,Left b); Right (a,c) -> (a,Right c) ) v
 
@@ -132,11 +132,11 @@
 -- Left (e1,e2)
 
 
-ev :: (Num k, Ord b) => Vect k (Tensor (Dual b) b) -> k
+ev :: (Eq k, Num k, Ord b) => Vect k (Tensor (Dual b) b) -> k
 ev = unwrap . linear (\(Dual bi, bj) -> delta bi bj *> return ())
 -- slightly cheating, as delta i j is meant to compare indices, not the basis elements themselves
 
 delta i j = if i == j then 1 else 0
 
-reify :: (Num k, Ord b) => Vect k (Dual b) -> (Vect k b -> k)
+reify :: (Eq k, Num k, Ord b) => Vect k (Dual b) -> (Vect k b -> k)
 reify f x = ev (f `te` x)
diff --git a/Math/Algebras/VectorSpace.hs b/Math/Algebras/VectorSpace.hs
--- a/Math/Algebras/VectorSpace.hs
+++ b/Math/Algebras/VectorSpace.hs
@@ -19,7 +19,7 @@
 -- Elements of Vect k b consist of k-linear combinations of elements of b.
 newtype Vect k b = V [(b,k)] deriving (Eq,Ord)
 
-instance (Num k, Show b) => Show (Vect k b) where
+instance (Show k, Eq k, 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
@@ -52,11 +52,11 @@
 zerov = V []
 
 -- |Addition of vectors
-add :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
+add :: (Ord b, Eq k, 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
+(<+>) :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b
 (<+>) = add
 
 addmerge ((a,x):ts) ((b,y):us) =
@@ -68,33 +68,33 @@
 addmerge [] us = us
 
 -- |Sum of a list of vectors
-sumv :: (Ord b, Num k) => [Vect k b] -> Vect k b
+sumv :: (Ord b, Eq k, Num k) => [Vect k b] -> Vect k b
 sumv = foldl (<+>) zerov
 
 -- |Negation of vector
-neg :: (Num k) => Vect k b -> Vect k b
+neg :: (Eq k, Num k) => Vect k b -> Vect k b
 neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts
 
 -- |Subtraction of vectors
-(<->) :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
+(<->) :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b
 (<->) u v = u <+> neg v
 
 -- |Scalar multiplication (on the left)
-smultL :: (Num k) => k -> Vect k b -> Vect k b
+smultL :: (Eq k, 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
+(*>) :: (Eq k, 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 :: (Eq k, 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
+(<*) :: (Eq k, Num k) => Vect k b -> k -> Vect k b
 (<*) = smultR
 
 -- same as return
@@ -107,7 +107,7 @@
 
 -- |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 :: (Ord b, Eq k, 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
@@ -135,7 +135,7 @@
 --
 -- If we have A = Vect k a, B = Vect k b, and f :: a -> Vect k b is a function from the basis elements of A into B,
 -- then @linear f@ is the linear map that this defines by linearity.
-linear :: (Ord b, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b
+linear :: (Ord b, Eq k, 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)
@@ -154,7 +154,7 @@
 -- but in the code, we need this if we want to be able to put k as one side of a tensor product.
 type Trivial k = Vect k ()
 
-wrap :: Num k => k -> Vect k ()
+wrap :: (Eq k, Num k) => k -> Vect k ()
 wrap 0 = zero
 wrap x = V [( (),x)]
 
diff --git a/Math/Combinatorics/Design.hs b/Math/Combinatorics/Design.hs
--- a/Math/Combinatorics/Design.hs
+++ b/Math/Combinatorics/Design.hs
@@ -8,6 +8,8 @@
 import qualified Data.Set as S
 
 import Math.Common.ListSet (intersect, symDiff)
+import Math.Core.Utils (combinationsOf)
+
 import Math.Algebra.Field.Base
 import Math.Algebra.Field.Extension
 import Math.Algebra.Group.PermutationGroup hiding (elts, order, isMember)
diff --git a/Math/Combinatorics/FiniteGeometry.hs b/Math/Combinatorics/FiniteGeometry.hs
--- a/Math/Combinatorics/FiniteGeometry.hs
+++ b/Math/Combinatorics/FiniteGeometry.hs
@@ -1,4 +1,4 @@
--- Copyright (c) David Amos, 2008-2009. All rights reserved.
+-- Copyright (c) David Amos, 2008-2011. All rights reserved.
 
 -- |Constructions of the finite geometries AG(n,Fq) and PG(n,Fq), their points, lines and flats,
 -- together with the incidence graphs between points and lines.
@@ -8,15 +8,15 @@
 import qualified Data.Set as S
 
 import Math.Common.ListSet (toListSet)
+import Math.Core.Utils
 
-import Math.Algebra.Field.Base
-import Math.Algebra.Field.Extension hiding ( (<+>) ) -- , (*>) )
+import Math.Core.Field
 import Math.Algebra.LinearAlgebra -- hiding ( det )
 
 import Math.Combinatorics.Graph
 import Math.Combinatorics.GraphAuts -- for use in GHCi
-import Math.Algebra.Group.PermutationGroup -- for use in GHCi
-import Math.Algebra.Group.SchreierSims as SS -- for use in GHCi
+import Math.Algebra.Group.PermutationGroup hiding (elts) -- for use in GHCi
+import Math.Algebra.Group.SchreierSims as SS hiding (elts) -- for use in GHCi
 
 -- !! The following two functions previously required (FiniteField a) as context
 -- but this has been temporarily removed to enable them to work with Math.Core.Field
@@ -47,13 +47,13 @@
 -- then this is the same as [v <*>> m | v <- vs] == [m' <<*> v | v <- vs]
 
 -- |Given a list of points in AG(n,Fq), return their closure, the smallest flat containing them
-closureAG :: (Ord a, FiniteField a) => [[a]] -> [[a]]
+closureAG :: (Num a, Ord a, FinSet a) => [[a]] -> [[a]]
 closureAG ps =
     let vs = [ (1 - sum xs) : xs | xs <- ptsAG (k-1) fq ] -- k-vectors over fq whose sum is 1
     in toListSet [m' <<*> v | v <- vs]
     where k = length ps -- the dimension of the flat (assuming ps are independent)
           m' = L.transpose ps
-          fq = eltsFq undefined
+          fq = elts
 -- toListSet call sorts the result, and also removes duplicates in case the points weren't independent
 
 {-
@@ -66,21 +66,21 @@
 
 lineAG [p1,p2] = L.sort [ p1 <+> (c *> dp) | c <- fq ] where
     dp = p2 <-> p1
-    fq = eltsFq undefined
+    fq = elts
 
 -- closure of points in PG(n,Fq)
 -- take all linear combinations of the points (ie the subspace generated by the points, considered as points in Fq ^(n+1) )
 -- then discard all which aren't in PNF (thus dropping back into PG(n,Fq))
 
 -- |Given a set of points in PG(n,Fq), return their closure, the smallest flat containing them
-closurePG :: (Ord a, FiniteField a) => [[a]] -> [[a]]
+closurePG :: (Num a, Ord a, FinSet a) => [[a]] -> [[a]]
 closurePG ps = toListSet $ filter ispnf $ map (<*>> ps) $ ptsAG k fq where
     k = length ps
-    fq = eltsFq undefined
+    fq = elts
 -- toListSet call sorts the result, and also removes duplicates in case the points weren't independent
 
 linePG [p1,p2] = toListSet $ filter ispnf [(a *> p1) <+> (b *> p2) | a <- fq, b <- fq]
-    where fq = eltsFq undefined
+    where fq = elts
 
 -- van Lint & Wilson, p325, 332
 qtorial n q | n >= 0 = product [(q^i - 1) `div` (q-1) | i <- [1..n]]
@@ -105,11 +105,9 @@
 -- Cameron, Combinatorics, p126
 
 
-
 -- FLATS VIA REDUCED ROW ECHELON FORMS
 -- Suggested by Cameron p125
 
-
 data ZeroOneStar = Zero | One | Star deriving (Eq)
 
 instance Show ZeroOneStar where
@@ -183,7 +181,7 @@
 -- INCIDENCE GRAPH
 
 -- |Incidence graph of PG(n,fq), considered as an incidence structure between points and lines
-incidenceGraphPG :: (Ord a, FiniteField a) => Int -> [a] -> Graph (Either [a] [[a]])
+incidenceGraphPG :: (Num a, Ord a, FinSet a) => Int -> [a] -> Graph (Either [a] [[a]])
 incidenceGraphPG n fq = G vs es where
     points = ptsPG n fq
     lines = linesPG n fq
@@ -200,7 +198,7 @@
 
 
 -- |Incidence graph of AG(n,fq), considered as an incidence structure between points and lines
-incidenceGraphAG :: (Ord a, FiniteField a) => Int -> [a] -> Graph (Either [a] [[a]])
+incidenceGraphAG :: (Num a, Ord a, FinSet a) => Int -> [a] -> Graph (Either [a] [[a]])
 incidenceGraphAG n fq = G vs es where
     points = ptsAG n fq
     lines = linesAG n fq
diff --git a/Math/Combinatorics/Graph.hs b/Math/Combinatorics/Graph.hs
--- a/Math/Combinatorics/Graph.hs
+++ b/Math/Combinatorics/Graph.hs
@@ -12,6 +12,7 @@
 import Control.Arrow ( (&&&) )
 
 import Math.Common.ListSet as LS
+import Math.Core.Utils
 import Math.Algebra.Group.PermutationGroup hiding (fromDigits, fromBinary)
 import Math.Algebra.Group.SchreierSims as SS
 
@@ -27,14 +28,6 @@
 powerset [] = [[]]
 powerset (x:xs) = let p = powerset xs in p ++ map (x:) p
 
--- |combinationsOf k xs returns the subsets of xs of size k.
--- If xs is in ascending order, then the returned list is in ascending order
-combinationsOf :: (Integral t) => t -> [a] -> [[a]]
-combinationsOf 0 _ = [[]]
-combinationsOf _ [] = []
-combinationsOf k (x:xs) | k > 0 = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs
-
-
 -- GRAPH
 
 -- |Datatype for graphs, represented as a list of vertices and a list of edges.
@@ -338,7 +331,7 @@
 -- kneser v k | v >= 2*k = j v k 0
 -- |kneser n k returns the kneser graph KG n,k -
 -- whose vertices are the k-element subsets of [1..n], with edges joining disjoint subsets
-kneser :: (Integral t) => t -> t -> Graph [t]
+kneser :: Int -> Int -> Graph [Int]
 kneser n k | 2*k <= n = graph (vs,es) where
     vs = combinationsOf k [1..n]
     es = [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, disjoint v1 v2]
diff --git a/Math/Combinatorics/GraphAuts.hs b/Math/Combinatorics/GraphAuts.hs
--- a/Math/Combinatorics/GraphAuts.hs
+++ b/Math/Combinatorics/GraphAuts.hs
@@ -2,12 +2,14 @@
 
 module Math.Combinatorics.GraphAuts where
 
+import Data.Either (lefts)
 import qualified Data.List as L
 import qualified Data.Map as M
 import qualified Data.Set as S
 import Data.Maybe
 
 import Math.Common.ListSet
+import Math.Core.Utils (combinationsOf, pairs)
 import Math.Combinatorics.Graph
 -- import Math.Combinatorics.StronglyRegularGraph
 -- import Math.Combinatorics.Hypergraph -- can't import this, creates circular dependency
@@ -281,13 +283,27 @@
        -- else error (show (src_split, trg_split)) -- for debugging
 
 -- Now, every time we intersect two partitions, refine to an equitable partition
--- |Given a graph g, @graphAuts g@ returns a strong generating set for the automorphism group of g.
---
--- Note that the implementation is currently only valid for connected graphs
+
+-- |Given a graph g, @graphAuts g@ returns generators for the automorphism group of g.
+-- If g is connected, then the generators will be a strong generating set.
 graphAuts :: (Ord a) => Graph a -> [Permutation a]
-graphAuts g@(G vs es)
+graphAuts g = autsWithinComponents ++ isosBetweenComponents
+    where cs = map (inducedSubgraph g) (components g)
+          -- autsWithinComponents = concatMap graphAutsCon cs
+          autsWithinComponents = concatMap graphAuts4 cs
+          isosBetweenComponents = map swapFromIso $ concat [take 1 (graphIsos ci cj) | (ci,cj) <- pairs cs]
+          swapFromIso xys = fromPairs (xys ++ map swap xys)
+          swap (x,y) = (y,x)
+-- Using graphAuts4 instead of graphAutsCon as latter appears to have a bug, eg
+-- > graphAuts4 $ G [1..3] [[1,2],[2,3]]
+-- [[[1,3]]]
+-- > graphAutsCon $ G [1..3] [[1,2],[2,3]]
+-- []
+
+-- Automorphisms of a connected graph
+graphAutsCon g@(G vs es)
     | isConnected g = graphAuts' [] (toEquitable g $ valencyPartition g)
-    | otherwise = error "graphAuts: only implemented for connected graphs"
+    | otherwise = error "graphAutsCon: graph is not connected"
     where graphAuts' us p@((x:ys):pt) =
               let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
               in level us p x ys []
@@ -305,8 +321,6 @@
           dps = M.fromList [(v, distancePartition g v) | v <- vs]
           es' = S.fromList es
           nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
--- To handle disconnected graphs, you not only need to find auts of each component,
--- you also need to find auts that swap components
 
 dfsEquitable (dps,es',nbrs_g) xys p1 p2 = dfs xys p1 p2 where
     dfs xys p1 p2
@@ -333,17 +347,24 @@
 -- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH
 
 -- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
+
 -- |Given the incidence graph of an incidence structure between points and blocks
 -- (for example, a set system),
--- @incidenceAuts g@ returns a strong generating set for the automorphism group of the incidence structure.
+-- @incidenceAuts g@ returns generators for the automorphism group of the incidence structure.
 -- The generators are represented as permutations of the points.
 -- The incidence graph should be represented with the points on the left and the blocks on the right.
---
--- Note that the implementation is currently only valid for connected incidence graphs
+-- If the incidence graph is connected, then the generators will be a strong generating set.
 incidenceAuts :: (Ord p, Ord b) => Graph (Either p b) -> [Permutation p]
-incidenceAuts g@(G vs es) 
+incidenceAuts g = autsWithinComponents ++ isosBetweenComponents
+    where cs = map (inducedSubgraph g) (components g)
+          autsWithinComponents = concatMap incidenceAutsCon cs
+          isosBetweenComponents = map swapFromIso $ concat [take 1 (incidenceIsos ci cj) | (ci,cj) <- pairs cs]
+          swapFromIso xys = fromPairs (xys ++ map swap xys)
+          swap (x,y) = (y,x)
+
+incidenceAutsCon g@(G vs es) 
     | isConnected g = map points (incidenceAuts' [] [vs])
-    | otherwise = error "incidenceAuts: only implemented for connected incidence graphs"
+    | otherwise = error "incidenceAutsCon: graph is not connected"
     where points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
           incidenceAuts' us p@((x@(Left _):ys):pt) =
               -- let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
@@ -351,8 +372,8 @@
               in level us p x ys []
               ++ incidenceAuts' (x:us) p'
           incidenceAuts' us ([]:pt) = incidenceAuts' us pt
-          incidenceAuts' _ (((Right _):_):_) = [] -- if we fix all the points, then the blocks must be fixed too 
-          -- incidenceAuts' _ [] = []
+          incidenceAuts' _ (((Right _):_):_) = [] -- if we fix all the points, then the blocks must be fixed too
+          incidenceAuts' _ [] = []
           level us p@(ph:pt) x (y@(Left _):ys) hs =
               let px = refine' (L.delete x ph : pt) (dps M.! x)
                   py = refine' (L.delete y ph : pt) (dps M.! y)
@@ -370,13 +391,28 @@
 
 -- !! not yet using equitable partitions, so could probably be more efficient
 
+
 -- graphIsos :: (Ord a, Ord b) => Graph a -> Graph b -> [[(a,b)]]
-graphIsos g1 g2 
+graphIsos g1 g2
+    | length cs1 /= length cs2 = []
+    | otherwise = graphIsos' cs1 cs2
+    where cs1 = map (inducedSubgraph g1) (components g1)
+          cs2 = map (inducedSubgraph g2) (components g2)
+          graphIsos' (ci:cis) cjs =
+              [iso ++ iso' | cj <- cjs,
+                             iso <- graphIsosCon ci cj,
+                             let cjs' = L.delete cj cjs,
+                             iso' <- graphIsos' cis cjs']
+          graphIsos' [] [] = [[]]
+
+-- isos between connected graphs
+graphIsosCon g1 g2 
     | isConnected g1 && isConnected g2
-        = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2) | v2 <- vertices g2]
-    | otherwise = error "graphIsos: only implemented for connected graphs"
-    where v1 = head $ vertices g1
-          dfs xys p1 p2
+        = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2)
+                 | v1 <- take 1 (vertices g1), v2 <- vertices g2]
+                 -- the take 1 handles the case where g1 is the null graph
+    | otherwise = error "graphIsosCon: either or both graphs are not connected"
+    where dfs xys p1 p2
               | map length p1 /= map length p2 = []
               | otherwise =
                   let p1' = filter (not . null) p1
@@ -397,10 +433,12 @@
           es2 = S.fromList $ edges g2
 
 
+-- |Are the two graphs isomorphic?
+isGraphIso :: (Ord a, Ord b) => Graph a -> Graph b -> Bool
+isGraphIso g1 g2 = (not . null) (graphIsos g1 g2)
 -- !! If we're only interested in seeing whether or not two graphs are iso,
 -- !! then the cost of calculating distancePartitions may not be warranted
 -- !! (see Math.Combinatorics.Poset: orderIsos01 versus orderIsos)
-isGraphIso g1 g2 = (not . null) (graphIsos g1 g2)
 
 -- !! deprecate
 isIso g1 g2 = (not . null) (graphIsos g1 g2)
@@ -411,12 +449,26 @@
 -- we return only the action on the Lefts, and unLeft it
 -- incidenceIsos :: (Ord p1, Ord b1, Ord p2, Ord b2) =>
 --     Graph (Either p1 b1) -> Graph (Either p2 b2) -> [[(p1,p2)]]
+
 incidenceIsos g1 g2
+    | length cs1 /= length cs2 = []
+    | otherwise = incidenceIsos' cs1 cs2
+    where cs1 = map (inducedSubgraph g1) (filter (not . null . lefts) $ components g1)
+          cs2 = map (inducedSubgraph g2) (filter (not . null . lefts) $ components g2)
+          incidenceIsos' (ci:cis) cjs =
+              [iso ++ iso' | cj <- cjs,
+                             iso <- incidenceIsosCon ci cj,
+                             let cjs' = L.delete cj cjs,
+                             iso' <- incidenceIsos' cis cjs']
+          incidenceIsos' [] [] = [[]]
+
+incidenceIsosCon g1 g2
     | isConnected g1 && isConnected g2
-        = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2) | v2@(Left _) <- vertices g2]
-    | otherwise = error "incidenceIsos: only implemented for connected graphs"
-    where v1@(Left _) = head $ vertices g1
-          dfs xys p1 p2
+        = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2)
+                 | v1@(Left _) <- take 1 (vertices g1), v2@(Left _) <- vertices g2]
+                 -- g1 may have no vertices
+    | otherwise = error "incidenceIsos: one or both graphs not connected"
+    where dfs xys p1 p2
               | map length p1 /= map length p2 = []
               | otherwise =
                   let p1' = filter (not . null) p1
@@ -436,6 +488,9 @@
           es1 = S.fromList $ edges g1
           es2 = S.fromList $ edges g2
 
+-- |Are the two incidence structures represented by these incidence graphs isomorphic?
+isIncidenceIso :: (Ord p1, Ord b1, Ord p2, Ord b2) =>
+     Graph (Either p1 b1) -> Graph (Either p2 b2) -> Bool
 isIncidenceIso g1 g2 = (not . null) (incidenceIsos g1 g2)
 
 {-
diff --git a/Math/Combinatorics/Hypergraph.hs b/Math/Combinatorics/Hypergraph.hs
--- a/Math/Combinatorics/Hypergraph.hs
+++ b/Math/Combinatorics/Hypergraph.hs
@@ -5,6 +5,7 @@
 
 import qualified Data.List as L
 import Math.Common.ListSet
+import Math.Core.Utils (combinationsOf)
 import Math.Combinatorics.Graph hiding (incidenceMatrix)
 import Math.Algebra.Group.PermutationGroup (orbitB, p) -- needed for construction of Coxeter group
 
diff --git a/Math/Combinatorics/IncidenceAlgebra.hs b/Math/Combinatorics/IncidenceAlgebra.hs
--- a/Math/Combinatorics/IncidenceAlgebra.hs
+++ b/Math/Combinatorics/IncidenceAlgebra.hs
@@ -97,7 +97,7 @@
 -- with multiplication defined by concatenation of intervals.
 -- The incidence algebra can also be thought of as the vector space of functions from intervals to k, with multiplication
 -- defined by the convolution (f*g)(x,y) = sum [ f(x,z) g(z,y) | x <= z <= y ].
-instance (Num k, Ord a) => Algebra k (Interval a) where
+instance (Eq k, Num k, Ord a) => Algebra k (Interval a) where
     -- |Note that we are not able to give a generic definition of unit for the incidence algebra,
     -- because it depends on which poset we are working in,
     -- and that information is encoded at the value level rather than the type level. See unitIA.
@@ -111,14 +111,14 @@
 
 
 -- |The unit of the incidence algebra of a poset
-unitIA :: (Num k, Ord t) => Poset t -> Vect k (Interval t)
+unitIA :: (Eq k, Num k, Ord t) => Poset t -> Vect k (Interval t)
 unitIA poset@(Poset (set,_)) = sumv [return (Iv poset (x,x)) | x <- set]
 
 basisIA :: Num k => Poset t -> [Vect k (Interval t)]
 basisIA poset = [return (Iv poset xy) | xy <- intervals poset]
 
 -- |The zeta function of a poset
-zetaIA :: (Num k, Ord t) => Poset t -> Vect k (Interval t)
+zetaIA :: (Eq k, Num k, Ord t) => Poset t -> Vect k (Interval t)
 zetaIA poset = sumv $ basisIA poset
 
 -- Then for example, zeta^2 counts the number of points in each interval
@@ -132,7 +132,7 @@
 
 -- calculate the mobius function of a poset, with memoization
 -- |The Mobius function of a poset
-muIA :: (Num k, Ord t) => Poset t -> Vect k (Interval t)
+muIA :: (Eq k, Num k, Ord t) => Poset t -> Vect k (Interval t)
 muIA poset@(Poset (set,po)) = sumv [mus M.! (x,y) *> return (Iv poset (x,y)) | x <- set, y <- set]
     where mu (x,y) | x == y    = 1
                    | po x y    = negate $ sum [mus M.! (x,z) | z <- set, po x z, po z y, z /= y]
@@ -152,7 +152,7 @@
 -- Stanley, Enumerative Combinatorics I, p144
 -- |The inverse of an element in the incidence algebra of a poset.
 -- This is only defined for elements which are non-zero on all intervals (x,x)
-invIA :: (Fractional k, Ord t) => Vect k (Interval t) -> Maybe (Vect k (Interval t))
+invIA :: (Eq k, Fractional k, Ord t) => Vect k (Interval t) -> Maybe (Vect k (Interval t))
 invIA f | f == zerov = Nothing -- error "invIA 0"
         | any (==0) [f' (x,x) | x <- set] = Nothing -- error "invIA: not invertible"
         | otherwise = Just g
@@ -210,7 +210,7 @@
 -- INCIDENCE COALGEBRA
 -- Schmitt, Incidence Hopf Algebras
 
-instance (Num k, Ord a) => Coalgebra k (Interval a) where
+instance (Eq k, Num k, Ord a) => Coalgebra k (Interval a) where
     counit = unwrap . linear counit'
         where counit' (Iv _ (x,y)) = (if x == y then 1 else 0) *> return ()
     comult = linear comult'
@@ -228,7 +228,7 @@
 -- rather than of intervals themselves.
 -- Note that if this operation is to be performed repeatedly for the same poset,
 -- then it is more efficient to use @toIsoClasses' poset@, which memoizes the isomorphism class lookup table.
-toIsoClasses :: (Num k, Ord a) => Vect k (Interval a) -> Vect k (Interval a)
+toIsoClasses :: (Eq k, Num k, Ord a) => Vect k (Interval a) -> Vect k (Interval a)
 toIsoClasses v
     | v == zerov = zerov
     | otherwise = toIsoClasses' poset v
@@ -236,7 +236,7 @@
 
 -- |Given a poset, @toIsoClasses' poset@ is the linear map from the incidence Hopf algebra of the poset to itself,
 -- in which each interval is mapped to (the minimal representative of) its isomorphism class.
-toIsoClasses' :: (Num k, Ord a) => Poset a -> Vect k (Interval a) -> Vect k (Interval a)
+toIsoClasses' :: (Eq k, Num k, Ord a) => Poset a -> Vect k (Interval a) -> Vect k (Interval a)
 toIsoClasses' poset = linear isoRep
     where isoRep iv = case isoMap M.! iv of
                       Nothing  -> return iv
diff --git a/Math/Combinatorics/LatinSquares.hs b/Math/Combinatorics/LatinSquares.hs
--- a/Math/Combinatorics/LatinSquares.hs
+++ b/Math/Combinatorics/LatinSquares.hs
@@ -14,6 +14,7 @@
 import Math.Combinatorics.Graph
 import Math.Combinatorics.GraphAuts
 import Math.Combinatorics.StronglyRegularGraph
+import Math.Core.Utils (combinationsOf)
 
 
 -- LATIN SQUARES
diff --git a/Math/Combinatorics/Matroid.hs b/Math/Combinatorics/Matroid.hs
--- a/Math/Combinatorics/Matroid.hs
+++ b/Math/Combinatorics/Matroid.hs
@@ -171,13 +171,13 @@
 -- |Given a matrix, represented as a list of rows, number the columns [1..],
 -- and construct the matroid whose independent sets correspond to those sets of columns which are linearly independent
 -- (or in case there are repetitions, those multisets of columns which are sets, and which are linearly independent).
-vectorMatroid :: (Fractional k) => [[k]] -> Matroid Int
+vectorMatroid :: (Eq k, Fractional k) => [[k]] -> Matroid Int
 vectorMatroid = vectorMatroid' . L.transpose
 
 -- |Given a list of vectors (or rows of a matrix), number the vectors (rows) [1..], and construct the matroid whose independent sets
 -- correspond to those sets of vectors (rows) which are linearly independent
 -- (or in case there are repetitions, those multisets which are sets, and which are linearly independent).
-vectorMatroid' :: (Fractional k) => [[k]] -> Matroid Int
+vectorMatroid' :: (Eq k, Fractional k) => [[k]] -> Matroid Int
 vectorMatroid' vs = fromBases (map fst vs') bs
     where vs' = zip [1..] vs
           bs = dfs [] [([],[],vs')]
@@ -247,6 +247,13 @@
           es' = L.sort [ [Left e, Right b] | b <- bs, e <- b ]
 -- incidence graph for the matroid considered as an incidence structure between elements and bases
 
+incidenceGraphC m = G.G vs' es'
+    where es = elements m
+          cs = L.sort $ circuits m
+          vs' = map Left es ++ map Right cs
+          es' = L.sort [ [Left e, Right c] | c <- cs, e <- c ]
+-- incidence graph for the matroid considered as an incidence structure between elements and circuits
+
 incidenceGraphH m = G.G vs' es'
     where es = elements m
           hs = L.sort $ hyperplanes m
@@ -266,12 +273,7 @@
 
 -- |Return the automorphisms of the matroid.
 matroidAuts :: (Ord a) => Matroid a -> [Permutation a]
-matroidAuts m
-    | G.isConnected hgraph = incidenceAuts hgraph
-    | G.isConnected bgraph = incidenceAuts bgraph
-    | otherwise = error "matroidAuts: incidence graph is not connected"
-    where hgraph = incidenceGraphH m
-          bgraph = incidenceGraphB m
+matroidAuts m = incidenceAuts $ incidenceGraphH m
 -- Note that the results aren't always what one intuitively expects from the geometric representation.
 -- This is because geometric representations suggest additional structure beyond matroid structure.
 -- For example, for the Vamos matroid v8,
@@ -613,7 +615,7 @@
 --
 -- A multiset of points in k^n is said to be affinely dependent if it contains two identical points,
 -- or three collinear points, or four coplanar points, or ... - and affinely independent otherwise.
-affineMatroid :: (Fractional k) => [[k]] -> Matroid Int
+affineMatroid :: (Eq k, Fractional k) => [[k]] -> Matroid Int
 affineMatroid vs = vectorMatroid' $ map (1:) vs
 
 -- |fromGeoRep returns a matroid from a geometric representation consisting of dependent flats of various ranks.
@@ -790,11 +792,11 @@
 -- REPRESENTABILITY
 
 -- |@matroidPG n fq@ returns the projective geometry PG(n,Fq), where fq is a list of the elements of Fq
-matroidPG :: (Fractional a) => Int -> [a] -> Matroid Int
+matroidPG :: (Eq a, Fractional a) => Int -> [a] -> Matroid Int
 matroidPG n fq = vectorMatroid' $ ptsPG n fq
 
 -- |@matroidAG n fq@ returns the affine geometry AG(n,Fq), where fq is a list of the elements of Fq
-matroidAG :: (Fractional a) => Int -> [a] -> Matroid Int
+matroidAG :: (Eq a, Fractional a) => Int -> [a] -> Matroid Int
 matroidAG n fq = vectorMatroid' $ ptsAG n fq
 
 
@@ -889,7 +891,7 @@
 
 -- |Find representations of the matroid m over fq. Specifically, this function will find one representative
 -- of each projective equivalence class of representation.
-representations :: (Fractional fq, Ord a) => [fq] -> Matroid a -> [[[fq]]]
+representations :: (Eq fq, Fractional fq, Ord a) => [fq] -> Matroid a -> [[[fq]]]
 representations fq m = map L.transpose $ representations' (reverse $ zip b ir) (zip b' dhash')
     where fq' = tail fq -- fq \ {0}
           b = head $ bases m
@@ -905,7 +907,7 @@
           representations' ls [] = [map snd $ reverse ls]
 
 -- |Is the matroid representable over Fq? For example, to find out whether a matroid m is binary, evaluate @isRepresentable f2 m@.
-isRepresentable :: (Fractional fq, Ord a) => [fq] -> Matroid a -> Bool
+isRepresentable :: (Eq fq, Fractional fq, Ord a) => [fq] -> Matroid a -> Bool
 isRepresentable fq m = (not . null) (representations fq m)
 
 -- |A binary matroid is a matroid which is representable over F2
@@ -1137,7 +1139,7 @@
 
 -- TODO
 
--- 1. Sort out the isomorphism / automorphism code in the case where the incidence graph isn't connected
+-- 1. Sort out the isomorphism code in the case where the incidence graph isn't connected
 
 -- 2. We could generate the geometric representation from a matroid (provided its rank <= 4)
 -- geoRep m = filter (isDependent m) (flats m)
diff --git a/Math/Combinatorics/Poset.hs b/Math/Combinatorics/Poset.hs
--- a/Math/Combinatorics/Poset.hs
+++ b/Math/Combinatorics/Poset.hs
@@ -146,7 +146,7 @@
 -- This is the projective geometry PG(n,q)
 -- |posetL n fq is the lattice of subspaces of the vector space Fq^n, ordered by inclusion.
 -- Subspaces are represented by their reduced row echelon form.
-posetL :: FiniteField fq => Int -> [fq] -> Poset [[fq]]
+posetL :: (Eq fq, FiniteField fq) => Int -> [fq] -> Poset [[fq]]
 posetL n fq = Poset ( subspaces fq n, isSubspace ) 
 
 
diff --git a/Math/Combinatorics/StronglyRegularGraph.hs b/Math/Combinatorics/StronglyRegularGraph.hs
--- a/Math/Combinatorics/StronglyRegularGraph.hs
+++ b/Math/Combinatorics/StronglyRegularGraph.hs
@@ -9,6 +9,7 @@
 import qualified Data.Set as S
 
 import Math.Common.ListSet
+import Math.Core.Utils (combinationsOf)
 import Math.Algebra.Group.PermutationGroup hiding (P)
 import Math.Algebra.Group.SchreierSims as SS
 import Math.Combinatorics.Graph as G hiding (G)
diff --git a/Math/Common/IntegerAsType.hs b/Math/Common/IntegerAsType.hs
--- a/Math/Common/IntegerAsType.hs
+++ b/Math/Common/IntegerAsType.hs
@@ -1,6 +1,6 @@
 -- Copyright (c) David Amos, 2009. All rights reserved.
 
-{-# OPTIONS_GHC -fglasgow-exts #-}
+{-# LANGUAGE EmptyDataDecls, ScopedTypeVariables #-}
 
 module Math.Common.IntegerAsType where
 
diff --git a/Math/CommutativeAlgebra/GroebnerBasis.hs b/Math/CommutativeAlgebra/GroebnerBasis.hs
--- a/Math/CommutativeAlgebra/GroebnerBasis.hs
+++ b/Math/CommutativeAlgebra/GroebnerBasis.hs
@@ -281,10 +281,10 @@
     where h // g = let ([u],_) = quotRemMP h [g] in u
 
 -- |@eliminate vs gs@ returns the elimination ideal obtained from the ideal generated by gs by eliminating the variables vs.
-eliminate :: (Fractional k, Ord k, MonomialConstructor m, Monomial (m v), Ord (m v)) =>
+eliminate :: (Eq k, Fractional k, Ord k, MonomialConstructor m, Monomial (m v), Ord (m v)) =>
     [Vect k (m v)] -> [Vect k (m v)] -> [Vect k (m v)]
 eliminate vs gs = let subs = subFst vs in eliminateFst [g `bind` subs | g <- gs]
-    where subFst :: (Num k, MonomialConstructor m, Eq (m v), Mon (m v)) =>
+    where subFst :: (Eq k, Num k, MonomialConstructor m, Eq (m v), Mon (m v)) =>
               [Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v))
           subFst vs = (\v -> let v' = var v in if v' `elem` vs then toElimFst v' else toElimSnd v')
 
diff --git a/Math/CommutativeAlgebra/Polynomial.hs b/Math/CommutativeAlgebra/Polynomial.hs
--- a/Math/CommutativeAlgebra/Polynomial.hs
+++ b/Math/CommutativeAlgebra/Polynomial.hs
@@ -157,7 +157,7 @@
 lexvar v = return $ Lex $ M 1 [(v,1)]
 -- lexvar = var
 
-instance (Num k, Ord v, Show v) => Algebra k (Lex v) where
+instance (Eq k, Num k, Ord v, Show v) => Algebra k (Lex v) where
     unit x = x *> return munit
     mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy
 
@@ -191,7 +191,7 @@
 glexvar v = return $ Glex $ M 1 [(v,1)]
 -- glexvar = var
 
-instance (Num k, Ord v, Show v) => Algebra k (Glex v) where
+instance (Eq k, Num k, Ord v, Show v) => Algebra k (Glex v) where
     unit x = x *> return munit
     mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy
 
@@ -227,7 +227,7 @@
 grevlexvar v = return $ Grevlex $ M 1 [(v,1)]
 -- grevlexvar = var
 
-instance (Num k, Ord v, Show v) => Algebra k (Grevlex v) where
+instance (Eq k, Num k, Ord v, Show v) => Algebra k (Grevlex v) where
     unit x = x *> return munit
     mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy
 
@@ -258,7 +258,7 @@
     mcoprime (Elim2 a1 b1) (Elim2 a2 b2) = mcoprime a1 a2 && mcoprime b1 b2
     mdeg (Elim2 a b) = mdeg a + mdeg b
 
-instance (Num k, Ord a, Mon a, Ord b, Mon b) => Algebra k (Elim2 a b) where
+instance (Eq k, Num k, Ord a, Mon a, Ord b, Mon b) => Algebra k (Elim2 a b) where
     unit x = x *> return munit
     mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy
 
@@ -286,7 +286,7 @@
 -- This is occasionally useful.
 
 -- |bind performs variable substitution
-bind :: (Num k, MonomialConstructor m, Ord a, Show a, Algebra k a) =>
+bind :: (Eq k, Num k, MonomialConstructor m, Ord a, Show a, Algebra k a) =>
     Vect k (m v) -> (v -> Vect k a) -> Vect k a
 v `bind` f = linear (\m -> product [f x ^ i | (x,i) <- mindices m]) v
 -- V ts `bind` f = sum [c *> product [f x ^ i | (x,i) <- mindices m] | (m, c) <- ts] 
@@ -298,7 +298,7 @@
 
 -- |Evaluate a polynomial at a point.
 -- For example @eval (x^2+y^2) [(x,1),(y,2)]@ evaluates x^2+y^2 at the point (x,y)=(1,2).
-eval :: (Num k, MonomialConstructor m, Eq (m v), Show v) =>
+eval :: (Eq k, Num k, MonomialConstructor m, Eq (m v), Show v) =>
     Vect k (m v) -> [(Vect k (m v), k)] -> k
 eval f vs = unwrap $ f `bind` sub
     where sub x = case lookup (var x) vs of
@@ -307,7 +307,7 @@
 
 -- |Perform variable substitution on a polynomial.
 -- For example @subst (x*z-y^2) [(x,u^2),(y,u*v),(z,v^2)]@ performs the substitution x -> u^2, y -> u*v, z -> v^2.
-subst :: (Num k, MonomialConstructor m, Eq (m u), Show u, Ord (m v), Show (m v), Algebra k (m v)) =>
+subst :: (Eq k, Num k, MonomialConstructor m, Eq (m u), Show u, Ord (m v), Show (m v), Algebra k (m v)) =>
     Vect k (m u) -> [(Vect k (m u), Vect k (m v))] -> Vect k (m v)
 subst f vs = f `bind` sub
     where sub x = case lookup (var x) vs of
@@ -392,7 +392,7 @@
 -- In the case where the gs are a Groebner basis for an ideal I,
 -- then @f %% gs@ is the equivalence class representative of f in R/I,
 -- and is zero if and only if f is in I.
-(%%) :: (Fractional k, Monomial m, Ord m, Algebra k m) =>
+(%%) :: (Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
      Vect k m -> [Vect k m] -> Vect k m
 f %% gs = rewrite f gs
 -- f %% gs = r where (_,r) = quotRemMP f gs
@@ -402,7 +402,7 @@
 -- The instance is well-defined only for scalars, and gives an error if used on other values.
 -- The purpose of this is to allow entry of fractional scalars, in expressions such as @x/2@.
 -- On the other hand, an expression such as @2/x@ will return an error.
-instance (Fractional k, Monomial m, Ord m, Algebra k m) => Fractional (Vect k m) where
+instance (Eq k, Fractional k, Monomial m, Ord m, Algebra k m) => Fractional (Vect k m) where
     recip (V [(m,c)]) | m == munit = V [(m,1/c)]
                       | otherwise = error "Polynomial recip: only defined for scalars"
     fromRational x = V [(munit, fromRational x)]
diff --git a/Math/Core/Field.hs b/Math/Core/Field.hs
--- a/Math/Core/Field.hs
+++ b/Math/Core/Field.hs
@@ -17,6 +17,8 @@
 import Data.Bits
 import Data.List as L
 
+import Math.Core.Utils (FinSet, elts)
+
 -- |Q is just the rationals, but with a better show function than the Prelude version
 newtype Q = Q Rational deriving (Eq,Ord,Num,Fractional)
 
@@ -52,6 +54,8 @@
     recip (F2 1) = F2 1
     fromRational _ = error "F2.fromRational: not well defined"
 
+instance FinSet F2 where elts = f2
+
 -- |f2 is a list of the elements of F2
 f2 :: [F2]
 f2 = map fromInteger [0..1] -- :: [F2]
@@ -77,6 +81,8 @@
     recip (F3 x) = F3 x
     fromRational _ = error "F3.fromRational: not well defined"
 
+instance FinSet F3 where elts = f3
+
 -- |f3 is a list of the elements of F3
 f3 :: [F3]
 f3 = map fromInteger [0..2] -- :: [F3]
@@ -102,6 +108,8 @@
     recip (F5 x) = F5 $ (x^3) `mod` 5
     fromRational _ = error "F5.fromRational: not well defined"
 
+instance FinSet F5 where elts = f5
+
 -- |f5 is a list of the elements of F5
 f5 :: [F5]
 f5 = map fromInteger [0..4]
@@ -127,6 +135,8 @@
     recip (F7 x) = F7 $ (x^5) `mod` 7
     fromRational _ = error "F7.fromRational: not well defined"
 
+instance FinSet F7 where elts = f7
+
 -- |f7 is a list of the elements of F7
 f7 :: [F7]
 f7 = map fromInteger [0..6]
@@ -152,6 +162,8 @@
     recip (F11 x) = F11 $ (x^9) `mod` 11
     fromRational _ = error "F11.fromRational: not well defined"
 
+instance FinSet F11 where elts = f11
+
 -- |f11 is a list of the elements of F11
 f11 :: [F11]
 f11 = map fromInteger [0..10]
@@ -177,6 +189,8 @@
     recip (F13 x) = F13 $ (x5*x5*x) `mod` 13 where x5 = x^5 `mod` 13 -- 12^11 would overflow Int
     fromRational _ = error "F13.fromRational: not well defined"
 
+instance FinSet F13 where elts = f13
+
 -- |f13 is a list of the elements of F13
 f13 :: [F13]
 f13 = map fromInteger [0..12]
@@ -202,6 +216,8 @@
     recip (F17 x) = F17 $ (x5^3) `mod` 17 where x5 = x^5 `mod` 17 -- 16^15 would overflow Int
     fromRational _ = error "F17.fromRational: not well defined"
 
+instance FinSet F17 where elts = f17
+
 -- |f17 is a list of the elements of F17
 f17 :: [F17]
 f17 = map fromInteger [0..16]
@@ -227,6 +243,8 @@
     recip (F19 x) = F19 $ (x4^4*x) `mod` 19 where x4 = x^4 `mod` 19 -- 18^17 would overflow Int
     fromRational _ = error "F19.fromRational: not well defined"
 
+instance FinSet F19 where elts = f19
+
 -- |f19 is a list of the elements of F19
 f19 :: [F19]
 f19 = map fromInteger [0..18]
@@ -252,6 +270,8 @@
     recip (F23 x) = F23 $ (x5^4*x) `mod` 23 where x5 = x^5 `mod` 23 -- 22^21 would overflow Int
     fromRational _ = error "F23.fromRational: not well defined"
 
+instance FinSet F23 where elts = f23
+
 -- |f23 is a list of the elements of F23
 f23 :: [F23]
 f23 = map fromInteger [0..22]
@@ -290,6 +310,8 @@
     recip (F4 x) = F4 (x `xor` 1)
     fromRational _ = error "F4.fromRational: not well defined"
 
+instance FinSet F4 where elts = f4
+
 -- |f4 is a list of the elements of F4
 f4 :: [F4]
 f4 = L.sort $ 0 : powers a4
@@ -330,6 +352,8 @@
     recip x = x^6
     fromRational _ = error "F8.fromRational: not well defined"
 
+instance FinSet F8 where elts = f8
+
 -- |f8 is a list of the elements of F8
 f8 :: [F8]
 f8 = L.sort $ 0 : powers a8
@@ -373,6 +397,8 @@
     recip x = x^7
     fromRational _ = error "F9.fromRational: not well defined"
 
+instance FinSet F9 where elts = f9
+
 -- |f9 is a list of the elements of F9
 f9 :: [F9]
 f9 = L.sort $ 0 : powers a9
@@ -419,6 +445,8 @@
     recip x = x^14
     fromRational _ = error "F16.fromRational: not well defined"
 
+instance FinSet F16 where elts = f16
+
 -- |f16 is a list of the elements of F16
 f16 :: [F16]
 f16 = L.sort $ 0 : powers a16
@@ -459,8 +487,9 @@
     recip x = x^23
     fromRational _ = error "F25.fromRational: not well defined"
 
+instance FinSet F25 where elts = f25
+
 -- |f25 is a list of the elements of F25
 f25 :: [F25]
 f25 = L.sort $ 0 : powers a25
-
 
diff --git a/Math/Core/Utils.hs b/Math/Core/Utils.hs
--- a/Math/Core/Utils.hs
+++ b/Math/Core/Utils.hs
@@ -11,6 +11,14 @@
 
 toSet = S.toList . S.fromList
 
+-- Merge two ordered listsets. Elements appearing in both inputs appear only once in the output
+mergeSet (x:xs) (y:ys) =
+    case compare x y of
+    LT -> x : mergeSet xs (y:ys)
+    EQ -> x : mergeSet xs ys
+    GT -> y : mergeSet (x:xs) ys
+mergeSet xs ys = xs ++ ys
+
 pairs (x:xs) = map (x,) xs ++ pairs xs
 pairs [] = []
 
@@ -62,3 +70,21 @@
 -- |@choose n k@ is the number of ways of choosing k distinct elements from an n-set
 choose :: (Integral a) => a -> a -> a
 choose n k = product [n-k+1..n] `div` product [1..k]
+
+
+-- |The class of finite sets
+class FinSet x where
+    elts :: [x]
+
+-- |A class representing algebraic structures having an inverse operation.
+-- Although strictly speaking the Num precondition means that we are requiring the structure
+-- also to be a ring, we do sometimes bend the rules (eg permutation groups).
+-- Note also that we don't insist that every element has an inverse.
+class Num a => HasInverses a where
+    inverse :: a -> a
+
+infix 8 ^-
+
+-- |A trick: x^-1 returns the inverse of x
+(^-) :: (HasInverses a, Integral b) => a -> b -> a
+x ^- n = inverse x ^ n
diff --git a/Math/NumberTheory/Factor.hs b/Math/NumberTheory/Factor.hs
--- a/Math/NumberTheory/Factor.hs
+++ b/Math/NumberTheory/Factor.hs
@@ -1,6 +1,7 @@
 -- Copyright (c) 2006-2011, David Amos. All rights reserved.
 
-module Math.NumberTheory.Factor where
+-- |A module for finding prime factors.
+module Math.NumberTheory.Factor (pfactors) where
 
 import Math.NumberTheory.Prime
 import Data.Either (lefts)
@@ -105,7 +106,8 @@
 
 
 -- |List the prime factors of n (with multiplicity).
--- The algorithm uses trial division, followed by the elliptic curve method if necessary.
+-- The algorithm uses trial division to find small factors,
+-- followed if necessary by the elliptic curve method to find larger factors.
 -- The running time increases with the size of the second largest prime factor of n.
 -- It can find 10-digit prime factors in seconds, but can struggle with 20-digit prime factors.
 pfactors :: Integer -> [Integer]
diff --git a/Math/NumberTheory/Prime.hs b/Math/NumberTheory/Prime.hs
--- a/Math/NumberTheory/Prime.hs
+++ b/Math/NumberTheory/Prime.hs
@@ -2,7 +2,7 @@
 
 {-# LANGUAGE NoMonomorphismRestriction #-}
 
-
+-- |A module providing functions to test for primality, and find next and previous primes.
 module Math.NumberTheory.Prime where
 
 import System.Random
@@ -13,7 +13,7 @@
     | n > 1 = isNotDivisibleBy primes
     | otherwise = False
     where isNotDivisibleBy (d:ds) | d*d > n         = True
-                                  | n `mod` d == 0  = False
+                                  | n `rem` d == 0  = False
                                   | otherwise       = isNotDivisibleBy ds
 
 -- |A (lazy) list of the primes
@@ -58,7 +58,7 @@
 -- power_mod b t n == b^t mod n
 power_mod b t n = powerMod' b 1 t
     where powerMod' x y 0 = y
-          powerMod' x y t = powerMod' (x*x `mod` n) (if even t then y else x*y `mod` n) (t `div` 2)
+          powerMod' x y t = powerMod' (x*x `rem` n) (if even t then y else x*y `rem` n) (t `div` 2)
 
 isMillerRabinPrime' n
     | n >= 4 =
diff --git a/Math/Projects/KnotTheory/Braid.hs b/Math/Projects/KnotTheory/Braid.hs
--- a/Math/Projects/KnotTheory/Braid.hs
+++ b/Math/Projects/KnotTheory/Braid.hs
@@ -1,6 +1,6 @@
 -- Copyright (c) David Amos, 2008. All rights reserved.
 
-{-# OPTIONS_GHC -XFlexibleInstances -XTypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}
 
 module Math.Projects.KnotTheory.Braid where
 
diff --git a/Math/Projects/KnotTheory/IwahoriHecke.hs b/Math/Projects/KnotTheory/IwahoriHecke.hs
--- a/Math/Projects/KnotTheory/IwahoriHecke.hs
+++ b/Math/Projects/KnotTheory/IwahoriHecke.hs
@@ -1,6 +1,6 @@
 -- Copyright (c) David Amos, 2008. All rights reserved.
 
-{-# OPTIONS_GHC -XFlexibleInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
 
 module Math.Projects.KnotTheory.IwahoriHecke where
 
diff --git a/Math/Projects/KnotTheory/TemperleyLieb.hs b/Math/Projects/KnotTheory/TemperleyLieb.hs
--- a/Math/Projects/KnotTheory/TemperleyLieb.hs
+++ b/Math/Projects/KnotTheory/TemperleyLieb.hs
@@ -1,6 +1,6 @@
 -- Copyright (c) David Amos, 2008. All rights reserved.
 
-{-# OPTIONS_GHC -XFlexibleInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
 
 module Math.Projects.KnotTheory.TemperleyLieb where
 
diff --git a/Math/Projects/MiniquaternionGeometry.hs b/Math/Projects/MiniquaternionGeometry.hs
--- a/Math/Projects/MiniquaternionGeometry.hs
+++ b/Math/Projects/MiniquaternionGeometry.hs
@@ -5,10 +5,11 @@
 import qualified Data.List as L
 
 import Math.Common.ListSet as LS
+import Math.Core.Utils (combinationsOf)
 
 import Math.Algebra.Field.Base
 import Math.Combinatorics.FiniteGeometry (pnf, ispnf, orderPGL)
-import Math.Combinatorics.Graph (combinationsOf)
+-- import Math.Combinatorics.Graph
 import Math.Combinatorics.GraphAuts
 import Math.Algebra.Group.PermutationGroup hiding (order)
 import qualified Math.Algebra.Group.SchreierSims as SS
diff --git a/Math/Projects/RootSystem.hs b/Math/Projects/RootSystem.hs
--- a/Math/Projects/RootSystem.hs
+++ b/Math/Projects/RootSystem.hs
@@ -1,7 +1,5 @@
 -- Copyright (c) David Amos, 2008. All rights reserved.
 
-{-# OPTIONS_GHC -fglasgow-exts #-}
-
 module Math.Projects.RootSystem where
 
 import Data.Ratio
diff --git a/Math/QuantumAlgebra/Tangle.hs b/Math/QuantumAlgebra/Tangle.hs
--- a/Math/QuantumAlgebra/Tangle.hs
+++ b/Math/QuantumAlgebra/Tangle.hs
@@ -24,7 +24,7 @@
 
 -- type TensorAlgebra k a = Vect k [a]
 
-instance (Num k, Ord a) => Algebra k [a] where
+instance (Eq k, Num k, Ord a) => Algebra k [a] where
     unit 0 = zero -- V []
     unit x = V [(munit,x)]
     mult = nf . fmap (\(a,b) -> a `mmult` b)
diff --git a/Math/Test/TAlgebras/TGroupAlgebra.hs b/Math/Test/TAlgebras/TGroupAlgebra.hs
--- a/Math/Test/TAlgebras/TGroupAlgebra.hs
+++ b/Math/Test/TAlgebras/TGroupAlgebra.hs
@@ -4,35 +4,23 @@
 
 module Math.Test.TAlgebras.TGroupAlgebra where
 
+import Test.HUnit
 import Test.QuickCheck
 
-import Math.Algebra.Group.PermutationGroup
+import Math.Algebra.Group.PermutationGroup hiding (p)
 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.Core.Utils
 
 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
--}
+                   return $ nf $ V $ take 10 ts
 
 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)
@@ -53,3 +41,27 @@
 
 prop_HopfAlgebra_GroupAlgebra x = prop_HopfAlgebra x
     where types = x :: GroupAlgebra Integer
+
+
+quickCheckGroupAlgebra = do
+    putStrLn "Checking that group algebra is an algebra, coalgebra, bialgebra, and Hopf algebra..."
+    quickCheck prop_Algebra_GroupAlgebra
+    quickCheck prop_Coalgebra_GroupAlgebra
+    quickCheck prop_Bialgebra_GroupAlgebra
+    quickCheck prop_HopfAlgebra_GroupAlgebra
+
+
+testlistGroupAlgebra = TestList [
+    testlistLeftInverse,
+    testlistRightInverse
+    ]
+
+groupAlgebraElts = [ 1+p[[1,2,3]], 1+p[[1,2,3]]+p[[1,2],[3,4]], 1+2*p[[1,2,3]]+p[[1,2],[3,4]] ]
+
+testcaseLeftInverse x = TestCase $ assertEqual ("inverse " ++ show x) 1 (x^-1 * x)
+
+testlistLeftInverse = TestList $ map testcaseLeftInverse groupAlgebraElts 
+
+testcaseRightInverse x = TestCase $ assertEqual ("inverse " ++ show x) 1 (x * x^-1)
+
+testlistRightInverse = TestList $ map testcaseRightInverse groupAlgebraElts 
diff --git a/Math/Test/TAlgebras/TOctonions.hs b/Math/Test/TAlgebras/TOctonions.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TAlgebras/TOctonions.hs
@@ -0,0 +1,39 @@
+-- Copyright (c) 2011, David Amos. All rights reserved.
+
+module Math.Test.TAlgebras.TOctonions where
+
+import Test.QuickCheck
+
+import Math.Algebra.Field.Base
+import Math.Algebras.VectorSpace
+import Math.Algebras.TensorProduct
+import Math.Algebras.Quaternions
+import Math.Algebras.Octonions
+
+import Math.Test.TAlgebras.TVectorSpace
+import Math.Test.TAlgebras.TStructures
+
+import Math.Algebras.Structures -- not really needed
+
+instance Arbitrary OBasis where
+    arbitrary = elements $ map O [-1..6]
+
+-- TVectorSpace defines an Arbitrary instance for Vect k b, given Arbitrary instances for k and b
+
+-- same as prop_Algebra, but missing associativity axiom
+prop_AlgebraNonAssociative (k,x) =
+    unitOutL (k' `te` x) == (mult . (unit' `tf` id)) (k' `te` x)  && -- left unit
+    unitOutR (x `te` k') == (mult . (id `tf` unit')) (x `te` k')     -- right unit
+    where k' = k *> return ()
+
+prop_AlgebraNonAssociative_Octonions (k,x) = prop_AlgebraNonAssociative (k,x)
+    where types = (k,x) :: (Q, Octonion Q)
+
+prop_InverseLoop (x,y) =
+    x*1 == x && x == 1*x &&
+    (x == 0 ||
+      (x^-1 * (x*y) == y && y == (y*x) * x^-1 &&
+       (x^-1)^-1 == x) )
+
+prop_InverseLoop_Octonions (x,y) = prop_InverseLoop (x,y)
+    where types = (x,y) :: (Octonion Q, Octonion Q)
diff --git a/Math/Test/TAlgebras/TTensorProduct.hs b/Math/Test/TAlgebras/TTensorProduct.hs
--- a/Math/Test/TAlgebras/TTensorProduct.hs
+++ b/Math/Test/TAlgebras/TTensorProduct.hs
@@ -15,6 +15,13 @@
 import Control.Category as C
 import Control.Arrow
 
+
+quickCheckTensorProduct = do
+    putStrLn "Testing that tf is linear, and tensor product is a functor"
+    quickCheck prop_Linear_tf
+    quickCheck prop_TensorFunctor
+
+
 type DirectSum k u v =
     (u ~ Vect k a, v ~ Vect k b) => Vect k (DSum a b)
 
diff --git a/Math/Test/TAlgebras/TVectorSpace.hs b/Math/Test/TAlgebras/TVectorSpace.hs
--- a/Math/Test/TAlgebras/TVectorSpace.hs
+++ b/Math/Test/TAlgebras/TVectorSpace.hs
@@ -28,7 +28,7 @@
 
 instance Arbitrary EBasis where
     arbitrary = do n <- arbitrary :: Gen Int
-                   return (E n)
+                   return (E $ abs n)
 
 instance Arbitrary b => Arbitrary (Dual b) where
     arbitrary = fmap Dual arbitrary
@@ -42,7 +42,8 @@
 
 instance (Num k, Ord b, Arbitrary k, Arbitrary b) => Arbitrary (Vect k b) where
     arbitrary = do ts <- arbitrary :: Gen [(b, k)] -- ScopedTypeVariables
-                   return $ nf $ V ts
+                   return $ nf $ V $ take 10 ts
+                   -- we impose complexity bound of 10 terms, to avoid unbounded running time.
 
 prop_VecSpQn (a,b,x,y,z) = prop_VecSp (a,b,x,y,z)
     where types = (a,b,x,y,z) :: (Q, Q, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis)
diff --git a/Math/Test/TCombinatorics/TFiniteGeometry.hs b/Math/Test/TCombinatorics/TFiniteGeometry.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TCombinatorics/TFiniteGeometry.hs
@@ -0,0 +1,67 @@
+-- Copyright (c) 2011, David Amos. All rights reserved
+
+module Math.Test.TCombinatorics.TFiniteGeometry where
+
+import Test.HUnit
+
+import Math.Combinatorics.FiniteGeometry
+import Math.Core.Field
+import Math.Combinatorics.GraphAuts (incidenceAuts)
+import Math.Algebra.Group.PermutationGroup (orderSGS)
+import Math.NumberTheory.Factor (pfactors)
+
+testlistFiniteGeometry = TestList [
+    testlistFlatsAG,
+    testlistFlatsPG,
+    testlistAutsAG,
+    testlistAutsPG
+    ]
+
+
+-- !! can't make list [f2,f3,f4], because they're different types
+testcaseFlatsAG n fq k = TestCase $
+    assertEqual (show "flatsAG " ++ show n ++ " " ++ show q ++ " " ++ show k)
+                (numFlatsAG n q k) (length (flatsAG n fq k))
+    where q = length fq
+
+testlistFlatsAG = TestList $
+    [testcaseFlatsAG n f2 k | n <- [2,3], k <- [0..n]] ++
+    [testcaseFlatsAG n f3 k | n <- [2,3], k <- [0..n]] ++
+    [testcaseFlatsAG n f4 k | n <- [2,3], k <- [0..n]]
+
+
+testcaseFlatsPG n fq k = TestCase $
+    assertEqual (show "flatsPG " ++ show n ++ " " ++ show q ++ " " ++ show k)
+                (numFlatsPG n q k) (length (flatsPG n fq k))
+    where q = length fq
+
+testlistFlatsPG = TestList $
+    [testcaseFlatsPG n f2 k | n <- [2,3], k <- [0..n]] ++
+    [testcaseFlatsPG n f3 k | n <- [2,3], k <- [0..n]] ++
+    [testcaseFlatsPG n f4 k | n <- [2,3], k <- [0..n]]
+
+
+testcaseAutsAG n fq = TestCase $
+    assertEqual ("autsAG " ++ show n ++ " " ++ show q)
+                (orderAff n q * degree)
+                (orderSGS $ incidenceAuts $ incidenceGraphAG n fq)
+    where q = toInteger $ length fq
+          degree = toInteger $ length $ pfactors $ toInteger q
+
+testlistAutsAG = TestList $ 
+    -- [testcaseAutsAG n f2 | n <- [2,3] ] ++ -- this is the complete graph, so has more auts than expected
+    [testcaseAutsAG n f3 | n <- [2] ] -- ++
+    -- [testcaseAutsAG n f4 | n <- [2,3] ] -- these take too long
+
+
+testcaseAutsPG n fq = TestCase $
+    assertEqual ("autsPG " ++ show n ++ " " ++ show q)
+                (orderPGL (n+1) q * degree)
+                (orderSGS $ incidenceAuts $ incidenceGraphPG n fq)
+    where q = toInteger $ length fq
+          degree = toInteger $ length $ pfactors $ toInteger q
+
+testlistAutsPG = TestList  $ 
+    [testcaseAutsPG n f2 | n <- [2,3] ] ++
+    [testcaseAutsPG n f3 | n <- [2] ] -- ++
+    -- [testcaseAutsPG n f4 | n <- [2,3] ] -- these take too long
diff --git a/Math/Test/TCombinatorics/TGraphAuts.hs b/Math/Test/TCombinatorics/TGraphAuts.hs
new file mode 100644
--- /dev/null
+++ b/Math/Test/TCombinatorics/TGraphAuts.hs
@@ -0,0 +1,127 @@
+-- Copyright (c) 2011, David Amos. All rights reserved.
+
+module Math.Test.TCombinatorics.TGraphAuts where
+
+import Data.List as L
+import Math.Core.Field hiding (f7)
+import Math.Core.Utils (combinationsOf)
+import Math.Algebra.Group.PermutationGroup as P
+import Math.Combinatorics.Graph as G
+import Math.Combinatorics.GraphAuts
+import Math.Combinatorics.Matroid as M
+
+import Test.HUnit
+
+
+testlistGraphAuts = TestList [
+    testlistGraphAutsOrder,
+    testlistGraphAutsGroup,
+    testlistGraphAutsComplement,
+    testlistIncidenceAutsOrder,
+    testlistGraphIsos,
+    testlistIsGraphIso,
+    testlistIncidenceIsos
+    ]
+
+
+-- We know the expected order of the graph automorphism group
+testcaseGraphAutsOrder desc g n = TestCase $
+    assertEqual ("order " ++ desc) n (orderSGS $ graphAuts g)
+
+testlistGraphAutsOrder = TestList [
+    let g = G [1..6] [[1,2],[3,4],[5,6]] in
+        testcaseGraphAutsOrder (show g) g 48, -- 2*2*2*3!
+    testcaseGraphAutsOrder "cube" cube 48,
+    testcaseGraphAutsOrder "dodecahedron" dodecahedron 120
+    ]
+
+
+induced bs g = fromPairs [(b, b -^ g) | b <- bs]
+
+-- We know the expected group of graph automorphisms
+testcaseGraphAutsGroup desc graph group = TestCase $
+    assertEqual ("group " ++ desc) (elts $ graphAuts graph) (elts $ group)
+
+testlistGraphAutsGroup = TestList [
+    testcaseGraphAutsGroup "nullGraph 0" (nullGraph 0) [],
+    testcaseGraphAutsGroup "nullGraph 1" (nullGraph 1) (_S 1),
+    testcaseGraphAutsGroup "nullGraph 2" (nullGraph 2) (_S 2),
+    testcaseGraphAutsGroup "nullGraph 3" (nullGraph 3) (_S 3),
+    testcaseGraphAutsGroup "k 3" (k 3) (_S 3),
+    testcaseGraphAutsGroup "k 4" (k 4) (_S 4),
+    testcaseGraphAutsGroup "k 5" (k 5) (_S 5),
+    testcaseGraphAutsGroup "c 4" (c 4) (_D 8),
+    testcaseGraphAutsGroup "c 5" (c 5) (_D 10),
+    let graph = G [1..3] [[1,2],[2,3]] in
+        testcaseGraphAutsGroup (show graph) graph [p [[1,3]]], -- regression test
+    let graph = G [1..6] [[2,3],[4,5],[5,6]] in
+        testcaseGraphAutsGroup (show graph) graph [p [[2,3]], p [[4,6]]],
+    testcaseGraphAutsGroup "petersen" petersen (map (induced $ combinationsOf 2 [1..5]) $ _S 5)
+    ]
+
+
+-- The automorphisms of the graph should be the same as the auts of its complement
+testcaseGraphAutsComplement desc g = TestCase $
+    assertEqual ("complement " ++ desc) (elts $ graphAuts g) (elts $ graphAuts $ complement g)
+-- the algorithm may not find the same set of generators, so we have to compare the elements
+
+testlistGraphAutsComplement = TestList [
+    testcaseGraphAutsComplement "k 3" (k 3),
+    testcaseGraphAutsComplement "kb 2 3" (kb 2 3), -- complement is not connected
+    testcaseGraphAutsComplement "kb 3 3" (kb 3 3), -- complement is not connected, but components can be swapped
+    testcaseGraphAutsComplement "kt 2 3 3" (kt 2 3 3),
+    testcaseGraphAutsComplement "kt 2 3 4" (kt 2 3 4),
+    testcaseGraphAutsComplement "kt 3 3 3" (kt 3 3 3)
+    ]
+
+kt a b c = graph (vs,es)
+    where vs = [1..a+b+c]
+          es = L.sort $ [[i,j] | i <- [1..a], j <- [a+1..a+b] ]
+                     ++ [[i,k] | i <- [1..a], k <- [a+b+1..a+b+c] ]
+                     ++ [[j,k] | j <- [a+1..a+b], k <- [a+b+1..a+b+c] ]
+
+-- We know the expected order of the incidence structure automorphism group
+testcaseIncidenceAutsOrder desc g n = TestCase $
+    assertEqual ("incidence order " ++ desc) n (P.order $ incidenceAuts g)
+
+-- We use matroids as our incidence structure just because we have a powerful library for constructing them
+testlistIncidenceAutsOrder = TestList [
+    testcaseIncidenceAutsOrder "pg2 f2 (B)" (incidenceGraphB $ matroidPG 2 f2) 168,
+    testcaseIncidenceAutsOrder "pg2 f2 (C)" (incidenceGraphC $ matroidPG 2 f2) 168,
+    testcaseIncidenceAutsOrder "pg2 f2 (H)" (incidenceGraphH $ matroidPG 2 f2) 168,
+    testcaseIncidenceAutsOrder "u 1 3 (B)" (incidenceGraphB $ u 1 3) 6, -- not connected
+    testcaseIncidenceAutsOrder "u 1 3 (C)" (incidenceGraphC $ u 1 3) 6,
+    testcaseIncidenceAutsOrder "u 1 3 (H)" (incidenceGraphH $ u 1 3) 6, -- not connected
+    testcaseIncidenceAutsOrder "u 2 3 `dsum` u 2 3 (H)" (incidenceGraphH $ u 2 3 `dsum` u 2 3) 72, -- 6*6*2
+    testcaseIncidenceAutsOrder "u 2 3 `dsum` u 2 3 (C)" (incidenceGraphC $ u 2 3 `dsum` u 2 3) 72, -- not connected
+    testcaseIncidenceAutsOrder "u 2 3 `dsum` u 3 4 (H)" (incidenceGraphH $ u 2 3 `dsum` u 3 4) 144, -- 6*24
+    testcaseIncidenceAutsOrder "u 2 3 `dsum` u 3 4 (C)" (incidenceGraphC $ u 2 3 `dsum` u 3 4) 144 -- not connected
+    ]
+
+
+testcaseGraphIsos g1 g2 isos = TestCase $
+    assertEqual (show (g1,g2)) isos (graphIsos g1 g2)
+
+testlistGraphIsos = TestList [
+    testcaseGraphIsos (G [1,2] []) (G [3,4] []) [[(1,3),(2,4)],[(1,4),(2,3)]],
+    testcaseGraphIsos (G [1,2,3] [[1,2]]) (G [4,5,6] [[5,6]]) [[(1,5),(2,6),(3,4)],[(1,6),(2,5),(3,4)]]
+    ]
+
+
+testcaseIsGraphIso g1 g2 = TestCase $
+    assertBool (show (g1,g2)) $ isGraphIso g1 g2
+
+testlistIsGraphIso = TestList [
+    testcaseIsGraphIso (nullGraph') (nullGraph')
+    ]
+
+
+testcaseIncidenceIsos g1 g2 isos = TestCase $
+    assertEqual (show (g1,g2)) isos (incidenceIsos g1 g2)
+
+testlistIncidenceIsos = TestList [
+    testcaseIncidenceIsos (G [Left 1, Right 2] []) (G [Left 3, Right 4] []) [[(1,3)]],
+    testcaseIncidenceIsos (G [Left 1, Left 2, Right 1] [[Left 1, Right 1]])
+                          (G [Left 3, Left 4, Right 4] [[Left 4, Right 4]])
+                          [[(1,4),(2,3)]]
+    ]
diff --git a/Math/Test/TCore/TField.hs b/Math/Test/TCore/TField.hs
--- a/Math/Test/TCore/TField.hs
+++ b/Math/Test/TCore/TField.hs
@@ -60,7 +60,8 @@
 
 
 quickCheckField =
-    do putStrLn "Testing F2..."
+    do putStrLn "Testing finite fields"
+       putStrLn "Testing F2..."
        quickCheck (prop_Field :: (F2,F2,F2) -> Bool)
        putStrLn "Testing F3..."
        quickCheck (prop_Field :: (F3,F3,F3) -> Bool)
diff --git a/Math/Test/TFiniteGeometry.hs b/Math/Test/TFiniteGeometry.hs
--- a/Math/Test/TFiniteGeometry.hs
+++ b/Math/Test/TFiniteGeometry.hs
@@ -1,8 +1,9 @@
 module Math.Test.TFiniteGeometry where
 
 import Math.Combinatorics.FiniteGeometry
-import Math.Algebra.Field.Base
-import Math.Algebra.Field.Extension
+import Math.Core.Field
+-- import Math.Algebra.Field.Base
+-- import Math.Algebra.Field.Extension
 import Math.Combinatorics.GraphAuts
 import Math.Algebra.Group.PermutationGroup
 
@@ -23,10 +24,10 @@
     ,numFlatsPG 3 4 1 == length (flatsPG 3 f4 1)
     ,numFlatsPG 3 4 2 == length (flatsPG 3 f4 2)
     ,numFlatsPG 3 4 3 == length (flatsPG 3 f4 3)
-    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f2) == orderAff 2 2 * toInteger (degree f2) 
-    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f3) == orderAff 2 3 * toInteger (degree f3)
-    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f4) == orderAff 2 4 * toInteger (degree f4)
-    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f2) == orderPGL 3 2 * toInteger (degree f2)
-    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f3) == orderPGL 3 3 * toInteger (degree f3)
-    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f4) == orderPGL 3 4 * toInteger (degree f4)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f2) == orderAff 2 2 -- * toInteger (degree f2) 
+    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f3) == orderAff 2 3 -- * toInteger (degree f3)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f4) == orderAff 2 4 * 2 -- * toInteger (degree f4)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f2) == orderPGL 3 2 -- * toInteger (degree f2)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f3) == orderPGL 3 3 -- * toInteger (degree f3)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f4) == orderPGL 3 4 * 2 -- * toInteger (degree f4)
     ]
diff --git a/Math/Test/TPermutationGroup.hs b/Math/Test/TPermutationGroup.hs
--- a/Math/Test/TPermutationGroup.hs
+++ b/Math/Test/TPermutationGroup.hs
@@ -6,6 +6,8 @@
 
 import qualified Data.List as L
 
+import Math.Core.Utils hiding (elts)
+
 import Math.Algebra.Group.PermutationGroup as P
 import Math.Algebra.Group.SchreierSims as SS
 import Math.Algebra.Group.RandomSchreierSims as RSS
diff --git a/Math/Test/TestAll.hs b/Math/Test/TestAll.hs
--- a/Math/Test/TestAll.hs
+++ b/Math/Test/TestAll.hs
@@ -11,7 +11,13 @@
 
 import Math.Test.TCore.TField
 
+import Math.Test.TAlgebras.TGroupAlgebra
+import Math.Test.TAlgebras.TOctonions
+import Math.Test.TAlgebras.TTensorAlgebra
+import Math.Test.TAlgebras.TTensorProduct
 import Math.Test.TCombinatorics.TDigraph
+import Math.Test.TCombinatorics.TFiniteGeometry
+import Math.Test.TCombinatorics.TGraphAuts
 import Math.Test.TCombinatorics.TIncidenceAlgebra
 import Math.Test.TCombinatorics.TMatroid
 import Math.Test.TCombinatorics.TPoset
@@ -21,7 +27,6 @@
 import Math.Test.TProjects.TMiniquaternionGeometry
 
 
-
 import Test.QuickCheck
 import Test.HUnit
 
@@ -37,14 +42,24 @@
 
 quickCheckAll =
     do
-    quickCheck prop_NonCommRingNPoly
+    -- quickCheck prop_NonCommRingNPoly
     quickCheck prop_GroupPerm
     quickCheckField
+    quickCheckTensorProduct
+    quickCheckGroupAlgebra
+    quickCheckTensorAlgebra
+    putStrLn "Testing Octonions..."
+    quickCheck prop_AlgebraNonAssociative_Octonions
+    quickCheck prop_InverseLoop_Octonions
+    putStrLn "Testing miniquaternion geometries..."
     quickCheck prop_NearFieldF9
     quickCheck prop_NearFieldJ9
 
 hunitAll = runTestTT $ TestList [
+    testlistGroupAlgebra,
     testlistDigraph,
+    testlistFiniteGeometry,
+    testlistGraphAuts,
     testlistIncidenceAlgebra,
     testlistMatroid,
     testlistPoset,