diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal
--- a/HaskellForMaths.cabal
+++ b/HaskellForMaths.cabal
@@ -1,5 +1,5 @@
    Name:                HaskellForMaths
-   Version:             0.4.5
+   Version:             0.4.6
    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 as an educational resource, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra
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,9 +1,11 @@
--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) David Amos, 2008-2015. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, ScopedTypeVariables, EmptyDataDecls, FlexibleInstances #-}
 
 module Math.Algebra.Field.Extension where
 
+import Prelude hiding ( (<*>) )
+
 import Data.Ratio
 import Data.List as L (elemIndex)
 
@@ -41,9 +43,9 @@
                         | i > 1  = v ++ "^" ++ show i -- "x^" ++ show i
 
 instance (Eq a, Num a) => Num (UPoly a) where
-    UP as + UP bs = toUPoly $ as <+> bs
+    UP as + UP bs = UP $ as <+> bs
     negate (UP as) = UP $ map negate as
-    UP as * UP bs = toUPoly $ as <*> bs
+    UP as * UP bs = UP $ as <*> bs
     fromInteger 0 = UP []
     fromInteger a = UP [fromInteger a]
     abs _ = error "Prelude.Num.abs: inappropriate abstraction"
@@ -51,14 +53,30 @@
 
 toUPoly as = UP (reverse (dropWhile (== 0) (reverse as)))
 
-(a:as) <+> (b:bs) = (a+b) : (as <+> bs)
+-- The fussiness of the code is to avoid adding trailing zeroes, eg [3] <+> [-3]
+-- Otherwise we would have to normalise after every addition
 as <+> [] = as
 [] <+> bs = bs
+-- (a:as) <+> (b:bs) = (a+b) : (as <+> bs)
+(a:as) <+> (b:bs) = let c = a+b
+                        cs = as <+> bs
+                    in if c == 0 && null cs then [] else c:cs
 
+-- The fussiness of the code is to avoid adding trailing zeroes.
+-- Note that since we call <+>, we rely on it having similar properties.
 [] <*> _ = []
+_ <*> [] = [] -- to avoid [0,1] <*> [] -> [0]
+(a:as) <*> bs = if null as then map (a*) bs else map (a*) bs <+> (0 : as <*> bs)
+
+-- > let valid xs = null xs || last xs /= 0
+-- > quickCheck (\as bs -> not (valid as) || not (valid bs) || valid (as <*> bs))
+
+{-
+-- The following definition introduces unnecessary trailing zeroes, eg [3] <*> [2] -> [6,0]
+[] <*> _ = []
 _ <*> [] = []
 (a:as) <*> (b:bs) = [a*b] <+> (0 : map (a*) bs) <+> (0 : map (*b) as) <+> (0 : 0 : as <*> bs)
-
+-}
 
 convert (UP as) = toUPoly $ map fromInteger as
 -- Can be used with type annotations to construct polynomials over other types, eg
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
@@ -31,6 +31,9 @@
 -- |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)
 
+-- Can't make a Functor instance because we need an Ord instance
+fmapP f = fromPairs . map (\(x,y) -> (f x, f y)) . toPairs
+
 -- |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
@@ -69,7 +72,8 @@
 xs -^ g = L.sort [x .^ g | x <- xs]
 
 -- construct a permutation from cycles
-fromCycles cs = fromPairs $ concatMap fromCycle cs
+-- fromCycles cs = fromPairs $ concatMap fromCycle cs
+fromCycles cs = product $ map (fromPairs . fromCycle) cs
     where fromCycle xs = zip xs (rotateL xs)
 
 -- convert a permutation to cycles
@@ -94,7 +98,7 @@
 
 -- |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
+instance Ord 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
@@ -104,13 +108,13 @@
     signum _ = error "(Permutation a).signum: not applicable"
 
 -- |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
+instance Ord a => HasInverses (Permutation a) where
     inverse (P g) = P $ M.fromList $ map (\(x,y)->(y,x)) $ M.toList g
 
 
 -- |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 a, Show a) => Permutation a -> Permutation a -> Permutation a
+(~^) :: Ord a => Permutation a -> Permutation a -> Permutation a
 g ~^ h = h^-1 * g * h
 
 -- commutator
@@ -119,6 +123,7 @@
 
 -- ORBITS
 
+{-
 closureS xs fs = closure' S.empty (S.fromList xs) where
     closure' interior boundary
         | S.null boundary = interior
@@ -126,6 +131,12 @@
             let interior' = S.union interior boundary
                 boundary' = S.fromList [f x | x <- S.toList boundary, f <- fs] S.\\ interior'
             in closure' interior' boundary'
+-}
+closureS xs fs = closure' S.empty xs where
+    closure' interior (x:xs)
+        | S.member x interior = closure' interior xs
+        | otherwise = closure' (S.insert x interior) ([f x | f <- fs] ++ xs)
+    closure' interior [] = interior
 
 closure xs fs = S.toList $ closureS xs fs
 
@@ -293,11 +304,19 @@
 -- For example, sgs (_A 5) == [[[1,2,3]],[[2,4,5]],[[3,4,5]]]
 -- So we need all three to generate the first transversal, then the last two to generate the second transversal, etc
 
--- |Given a strong generating set, return the order of the group it generates
+-- |Given a strong generating set, return the order of the group it generates.
+-- Note that the SGS is assumed to be relative to the natural order of the points on which the group acts.
 orderSGS :: (Ord a) => [Permutation a] -> Integer
 orderSGS sgs = product $ map (L.genericLength . fundamentalOrbit) bs where
     bs = toListSet $ map minsupp sgs
     fundamentalOrbit b = b .^^ filter ( (b <=) . minsupp ) sgs
+
+-- !! Needs more testing
+-- |Given a base and strong generating set, return the order of the group it generates.
+orderBSGS :: (Ord a) => ([a],[Permutation a]) -> Integer
+orderBSGS (bs,sgs) = go 1 bs sgs where
+    go n [] _ = n
+    go n (b:bs) gs = go (n * L.genericLength (b .^^ gs)) bs (filter (\g -> b .^ g == b) gs)
 
 
 -- MORE INVESTIGATIONS
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
@@ -6,7 +6,7 @@
 import Data.Maybe (isNothing, isJust)
 import qualified Data.Set as S
 import qualified Data.Map as M
-import Math.Algebra.Group.PermutationGroup hiding (elts, order, gens, isMember, isSubgp, isNormal, reduceGens, normalClosure, commutatorGp, derivedSubgp)
+import Math.Algebra.Group.PermutationGroup hiding (elts, order, orderBSGS, gens, isMember, isSubgp, isNormal, reduceGens, normalClosure, commutatorGp, derivedSubgp)
 import Math.Common.ListSet (toListSet)
 import Math.Core.Utils hiding (elts)
 
diff --git a/Math/Algebra/LinearAlgebra.hs b/Math/Algebra/LinearAlgebra.hs
--- a/Math/Algebra/LinearAlgebra.hs
+++ b/Math/Algebra/LinearAlgebra.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2008-2012, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, 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]].
@@ -9,6 +9,8 @@
 -- on each side indicates the dimension of the argument on that side. For example,
 -- v \<*\>\> m is multiplication of a vector on the left by a matrix on the right.
 module Math.Algebra.LinearAlgebra where
+
+import Prelude hiding ( (*>), (<*>) )
 
 import qualified Data.List as L
 import Math.Core.Field -- not actually used in this module
diff --git a/Math/Algebras/AffinePlane.hs b/Math/Algebras/AffinePlane.hs
--- a/Math/Algebras/AffinePlane.hs
+++ b/Math/Algebras/AffinePlane.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
 
@@ -38,7 +38,7 @@
 instance Show v => Show (SL2 v) where show (SL2 m) = show m
 
 instance Algebra Q (SL2 ABCD) where -- to do this for Num k instead of Q we would need a,b,c,d defined for Num k
-    unit 0 = zero -- V []
+    unit 0 = zerov -- V []
     unit x = V [(munit,x)] where munit = SL2 (Glex 0 [])
     mult x = x''' where
         x' = mult $ fmap ( \(SL2 a, SL2 b) -> (a,b) ) x -- perform the multiplication in GlexPoly
diff --git a/Math/Algebras/Commutative.hs b/Math/Algebras/Commutative.hs
--- a/Math/Algebras/Commutative.hs
+++ b/Math/Algebras/Commutative.hs
@@ -1,10 +1,16 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
 
--- |A module defining the algebra of commutative polynomials over a field k
+-- |A module defining the algebra of commutative polynomials over a field k.
+--
+-- Most users should probably use Math.CommutativeAlgebra.Polynomial instead, which is basically the same thing
+-- but more fully-featured. This module will probably be deprecated at some point, but remains for now because
+-- it has a simpler implementation which may be more helpful for people wanting to understand the code.
 module Math.Algebras.Commutative where
 
+import Prelude hiding ( (*>) )
+
 import Math.Algebra.Field.Base hiding (powers)
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
@@ -35,6 +41,10 @@
 instance Functor GlexMonomial where
     fmap f (Glex si xis) = Glex si [(f x, i) | (x,i) <- xis]
 -- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards
+
+instance Applicative GlexMonomial where
+    pure = return
+    (<*>) = ap
 
 -- GlexMonomial is the free commutative monoid, and hence a monad
 instance Monad GlexMonomial where
diff --git a/Math/Algebras/GroupAlgebra.hs b/Math/Algebras/GroupAlgebra.hs
--- a/Math/Algebras/GroupAlgebra.hs
+++ b/Math/Algebras/GroupAlgebra.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2010-2012, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances #-}
 -- ScopedTypeVariables
@@ -13,6 +13,8 @@
 -- Elements of the group algebra consist 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 (GroupAlgebra, p) where
+
+import Prelude hiding ( (*>) )
 
 import Math.Core.Field
 import Math.Core.Utils hiding (elts)
diff --git a/Math/Algebras/LaurentPoly.hs b/Math/Algebras/LaurentPoly.hs
--- a/Math/Algebras/LaurentPoly.hs
+++ b/Math/Algebras/LaurentPoly.hs
@@ -32,7 +32,7 @@
     mmult (LM si xis) (LM sj yjs) = LM (si+sj) $ addmerge xis yjs
 
 instance (Eq k, Num k) => Algebra k LaurentMonomial where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- V []
     unit x = V [(munit,x)] 
     mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts)
     -- mult (V ts) = nf $ V [(a `mmult` b, x) | (T a b, x) <- ts]
diff --git a/Math/Algebras/Matrix.hs b/Math/Algebras/Matrix.hs
--- a/Math/Algebras/Matrix.hs
+++ b/Math/Algebras/Matrix.hs
@@ -1,10 +1,12 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances #-}
 
 
 module Math.Algebras.Matrix where
 
+import Prelude hiding ( (*>) )
+
 import Math.Algebra.Field.Base
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
@@ -45,9 +47,9 @@
 toMat2 [[a,b],[c,d]] = sum $ zipWith (\x e -> unit x * return e) [a,b,c,d] [E2 1 1, E2 1 2, E2 2 1, E2 2 2]
 -- fromMat2
 
-toEB2 [x,y] = foldl add zero $ zipWith (\x e -> x `smultL` return e) [x,y] [E 1, E 2]
+toEB2 [x,y] = foldl add zerov $ zipWith (\x e -> x `smultL` return e) [x,y] [E 1, E 2]
 
-toEB xs = foldl add zero $ zipWith (\x e -> x `smultL` return e) xs (map E [1..])
+toEB xs = foldl add zerov $ zipWith (\x e -> x `smultL` return e) xs (map E [1..])
 
 
 
@@ -58,7 +60,7 @@
 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]])
+    comult = linear (\(E2' i j) -> foldl (<+>) zerov [return (E2' i k, E2' k j) | k <- [1..2]])
 -- In other words
 -- counit (a b) = (1 0)
 --        (c d)   (0 1)
@@ -74,7 +76,7 @@
 -- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere
 
 instance (Eq k, Num k) => Algebra k M3 where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- 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
     mult = linear mult' where
diff --git a/Math/Algebras/NonCommutative.hs b/Math/Algebras/NonCommutative.hs
--- a/Math/Algebras/NonCommutative.hs
+++ b/Math/Algebras/NonCommutative.hs
@@ -1,10 +1,12 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
 
 -- |A module defining the algebra of non-commutative polynomials over a field k
 module Math.Algebras.NonCommutative where
 
+import Prelude hiding ( (*>) )
+
 import Math.Algebra.Field.Base hiding (powers)
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
@@ -30,7 +32,7 @@
     mmult (NCM i xs) (NCM j ys) = NCM (i+j) (xs++ys)
 
 instance (Eq k, Num k, Ord v) => Algebra k (NonComMonomial v) where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- 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
@@ -1,4 +1,4 @@
--- Copyright (c) 2011, David Amos. All rights reserved.
+-- Copyright (c) 2011-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances, NoMonomorphismRestriction #-}
 
@@ -7,6 +7,8 @@
 -- The octonions are the algebra defined by the basis {1,i0,i1,i2,i3,i4,i5,i6},
 -- 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 Prelude hiding ( (*>) )
 
 import Math.Core.Field
 import Math.Algebras.VectorSpace
diff --git a/Math/Algebras/Quaternions.hs b/Math/Algebras/Quaternions.hs
--- a/Math/Algebras/Quaternions.hs
+++ b/Math/Algebras/Quaternions.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances, NoMonomorphismRestriction #-}
 
@@ -7,6 +7,8 @@
 -- The quaternions are the algebra defined by the basis {1,i,j,k}, where i^2 = j^2 = k^2 = ijk = -1
 module Math.Algebras.Quaternions where
 
+import Prelude hiding ( (*>) )
+
 import Math.Core.Field
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
@@ -68,7 +70,7 @@
     sqnorm :: Vect k a -> k
 
 -- |If an algebra has a conjugation operation, then it has multiplicative inverses,
--- via 1/x = conj x / sqnorm x
+-- via 1\/x = conj x \/ sqnorm x
 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
@@ -167,7 +169,7 @@
 instance (Eq k, Num k) => Coalgebra k (Dual HBasis) where
     counit = unwrap . linear counit'
         where counit' (Dual One) = return ()
-              counit' _ = zero
+              counit' _ = zerov
     comult = linear comult'
         where comult' (Dual One) = return (Dual One, Dual One) <+>
                   (-1) *> ( return (Dual I, Dual I) <+> return (Dual J, Dual J) <+> return (Dual K, Dual K) )
diff --git a/Math/Algebras/Structures.hs b/Math/Algebras/Structures.hs
--- a/Math/Algebras/Structures.hs
+++ b/Math/Algebras/Structures.hs
@@ -1,4 +1,4 @@
--- Copyright (c) David Amos, 2010-2012. All rights reserved.
+-- Copyright (c) David Amos, 2010-2015. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction #-}
 {-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}
@@ -8,6 +8,8 @@
 -- - specifically algebra, coalgebra, bialgebra, Hopf algebra, module, comodule
 module Math.Algebras.Structures where
 
+import Prelude hiding ( (*>) )
+
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
 
@@ -103,7 +105,7 @@
     mult = linear mult'
         where mult' (Left a1, Left a2) = i1 $ mult $ return (a1,a2)
               mult' (Right b1, Right b2) = i2 $ mult $ return (b1,b2)
-              mult' _ = zero
+              mult' _ = zerov
 -- This is the product algebra, which is the product in the category of algebras
 -- 1 = (1,1)
 -- (a1,b1) * (a2,b2) = (a1*a2, b1*b2)
diff --git a/Math/Algebras/TensorAlgebra.hs b/Math/Algebras/TensorAlgebra.hs
--- a/Math/Algebras/TensorAlgebra.hs
+++ b/Math/Algebras/TensorAlgebra.hs
@@ -5,6 +5,8 @@
 -- |A module defining the tensor algebra, symmetric algebra, exterior (or alternating) algebra, and tensor coalgebra
 module Math.Algebras.TensorAlgebra where
 
+import Prelude hiding ( (*>) )
+
 import qualified Data.List as L
 
 import Math.Algebras.VectorSpace
@@ -195,7 +197,7 @@
     mult xy = nf $ xy >>= (\(Ext i xs, Ext j ys) -> signedMerge 1 (0,[]) (i,xs) (j,ys))
         where signedMerge s (k,zs) (i,x:xs) (j,y:ys) =
                   case compare x y of
-                  EQ -> zero
+                  EQ -> zerov
                   LT -> signedMerge s (k+1,x:zs) (i-1,xs) (j,y:ys)
                   GT -> let s' = if even i then s else -s -- we had to commute y past x:xs, with i sign changes
                         in signedMerge s' (k+1,y:zs) (i,x:xs) (j-1,ys)
diff --git a/Math/Algebras/TensorProduct.hs b/Math/Algebras/TensorProduct.hs
--- a/Math/Algebras/TensorProduct.hs
+++ b/Math/Algebras/TensorProduct.hs
@@ -1,10 +1,12 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE NoMonomorphismRestriction #-}
 
 -- |A module defining direct sum and tensor product of vector spaces
 module Math.Algebras.TensorProduct where
 
+import Prelude hiding ( (*>) )
+
 import Math.Algebras.VectorSpace
 
 infix 7 `te`, `tf`
@@ -38,12 +40,12 @@
 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
+    p1' (Right b) = zerov
 
 -- |Projection onto right summand from direct sum
 p2 :: (Eq k, Num k, Ord b) => Vect k (DSum a b) -> Vect k b
 p2 = linear p2' where
-    p2' (Left a) = zero
+    p2' (Left a) = zerov
     p2' (Right b) = return b
 
 -- |The product of two linear functions (with the same source).
@@ -83,7 +85,7 @@
 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 [])
+    where sum = foldl add zerov
 
 
 -- tensor isomorphisms
diff --git a/Math/Algebras/VectorSpace.hs b/Math/Algebras/VectorSpace.hs
--- a/Math/Algebras/VectorSpace.hs
+++ b/Math/Algebras/VectorSpace.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE NoMonomorphismRestriction #-}
 {-# OPTIONS_HADDOCK prune #-}
@@ -6,12 +6,16 @@
 -- |A module defining the type and operations of free k-vector spaces over a basis b (for a field k)
 module Math.Algebras.VectorSpace where
 
+import Prelude hiding ( (*>) )
+
+import Control.Applicative hiding ( (<*), (*>) )
+import Control.Monad (ap)
 import qualified Data.List as L
 import qualified Data.Set as S -- only needed for toSet
 
 infixr 7 *>
 infixl 7 <*
-infixl 6 <+>, <->
+infixl 6 <+>, <->, <<+>>, <<->>
 
 
 -- |Given a field type k and a basis type b, Vect k b is the type of the free k-vector space over b.
@@ -48,24 +52,24 @@
 terms (V ts) = ts
 
 -- |Return the coefficient of the specified basis element in a vector
+coeff :: (Num k, Eq b) => b -> Vect k b -> k
 coeff b v = sum [k | (b',k) <- terms v, b' == b]
 
 -- |Remove the term for a specified basis element from a vector
-removeTerm b v = v <-> coeff b v *> return b
-
--- Deprecated
-zero = V []
+removeTerm :: (Eq k, Num k, Ord b) => b -> Vect k b -> Vect k b
+removeTerm b (V ts) = V $ filter ((/=b) . fst) ts
+-- v <-> coeff b v *> return b
 
 -- |The zero vector
 zerov :: Vect k b
 zerov = V []
 
 -- |Addition of vectors
-add :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b
+add :: (Eq k, Num k, Ord b) => 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, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b
+(<+>) :: (Eq k, Num k, Ord b) => Vect k b -> Vect k b -> Vect k b
 (<+>) = add
 
 addmerge ((a,x):ts) ((b,y):us) =
@@ -77,23 +81,20 @@
 addmerge [] us = us
 
 -- |Sum of a list of vectors
-sumv :: (Ord b, Eq k, Num k) => [Vect k b] -> Vect k b
+sumv :: (Eq k, Num k, Ord b) => [Vect k b] -> Vect k b
 sumv = foldl (<+>) zerov
 
--- Deprecated
-neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts
-
 -- |Negation of a vector
 negatev :: (Eq k, Num k) => Vect k b -> Vect k b
 negatev (V ts) = V $ map (\(b,x) -> (b,-x)) ts
 
 -- |Subtraction of vectors
-(<->) :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b
+(<->) :: (Eq k, Num k, Ord b) => Vect k b -> Vect k b -> Vect k b
 (<->) u v = u <+> negatev v
 
 -- |Scalar multiplication (on the left)
 smultL :: (Eq k, Num k) => k -> Vect k b -> Vect k b
-smultL 0 _ = zero -- V []
+smultL 0 _ = zerov -- V []
 smultL k (V ts) = V [(ei,k*xi) | (ei,xi) <- ts]
 
 -- |Same as smultL. Mnemonic is \"multiply through (from the left)\"
@@ -102,7 +103,7 @@
 
 -- |Scalar multiplication on the right
 smultR :: (Eq k, Num k) => Vect k b -> k -> Vect k b
-smultR _ 0 = zero -- V []
+smultR _ 0 = zerov -- V []
 smultR (V ts) k = V [(ei,xi*k) | (ei,xi) <- ts]
 
 -- |Same as smultR. Mnemonic is \"multiply through (from the right)\"
@@ -119,7 +120,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, Eq k, Num k) => Vect k b -> Vect k b
+nf :: (Eq k, Num k, Ord b) => 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
@@ -136,7 +137,7 @@
 --
 -- In the mathematical sense, this can be regarded as a functor from the category Set (of sets) to the category k-Vect
 -- (of k-vector spaces). In Haskell, instead of Set we have Hask, the category of Haskell types. However, for our purposes
--- it is helpful to identify Hask with Set, but identifying a Haskell type with its set of inhabitants.
+-- it is helpful to identify Hask with Set, by identifying a Haskell type with its set of inhabitants.
 --
 -- The type constructor (Vect k) gives the action of the functor on objects in the category,
 -- taking a set (type) to a free k-vector space. fmap gives the action of the functor on arrows in the category,
@@ -149,6 +150,16 @@
     fmap f (V ts) = V [(f b, x) | (b,x) <- ts]
 -- Note that if f is not order-preserving, then we need to call "nf" afterwards
 
+-- From GHC 7.10, Monad has Applicative as a superclass, so we must define an instance.
+-- It doesn't particularly make sense for Vect k.
+-- (Although given Vect k b, we could represent the dual space as Vect k (b -> ()),
+-- and then have a use for <*>.)
+instance Num k => Applicative (Vect k) where
+    pure = return
+    -- pure b = V [(b,1)]
+    (<*>) = ap
+    -- V fs <*> V xs = V [(f x, a*b) | (f,a) <- fs, (x,b) <- xs]
+
 -- |Given a field k, the type constructor (Vect k) is a monad, the \"free k-vector space monad\".
 --
 -- In order to understand this, it is probably easiest to think of a free k-vector space as a kind of container,
@@ -172,7 +183,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, Eq k, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b
+linear :: (Eq k, Num k, Ord b) => (a -> Vect k b) -> Vect k a -> Vect k b
 linear f v = nf $ v >>= f
 
 newtype EBasis = E Int deriving (Eq,Ord)
@@ -193,7 +204,7 @@
 
 -- |Wrap an element of the field k to an element of the trivial k-vector space
 wrap :: (Eq k, Num k) => k -> Vect k ()
-wrap 0 = zero
+wrap 0 = zerov
 wrap x = V [( (),x)]
 
 -- |Unwrap an element of the trivial k-vector space to an element of the field k
@@ -223,6 +234,20 @@
 
 (f <<+>> g) v = f v <+> g v
 
+(f <<->> g) v = f v <-> g v
+
 zerof v = zerov
 
 sumf fs = foldl (<<+>>) zerof fs
+
+
+-- Lens
+coeffLens :: (Ord b, Eq k, Num k, Functor f) => b -> (k -> f k) -> (Vect k b -> f (Vect k b))
+coeffLens b = lens (coeff b) (setter b)
+    where setter b = \(V ts) k -> (k *> return b) <+> (V $ filter ((/=b) . fst) ts)
+          lens getter setter f a = fmap (setter a) (f (getter a))
+-- Can be used with lens-family, for example
+-- e1 ^. coeffLens (E 2) --> 0
+-- e1 & coeffLens (E 2) .~ 2 --> e1+2e2
+-- e1 & coeffLens (E 1) %~ (+2) --> 3e1
+
diff --git a/Math/Combinatorics/CombinatorialHopfAlgebra.hs b/Math/Combinatorics/CombinatorialHopfAlgebra.hs
--- a/Math/Combinatorics/CombinatorialHopfAlgebra.hs
+++ b/Math/Combinatorics/CombinatorialHopfAlgebra.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2012, David Amos. All rights reserved.
+-- Copyright (c) 2012-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction, ScopedTypeVariables, DeriveFunctor #-}
 
@@ -34,6 +34,8 @@
 -- Lie Algebras and Hopf Algebras
 -- Michiel Hazewinkel, Nadiya Gubareni, V.V.Kirichenko
 
+import Prelude hiding ( (*>) )
+
 import Data.List as L
 import Data.Maybe (fromJust)
 import qualified Data.Set as S
@@ -198,17 +200,17 @@
 
 -- |Convert an element of SSym represented in the monomial basis to the fundamental basis
 ssymMtoF :: (Eq k, Num k) => Vect k SSymM -> Vect k SSymF
-ssymMtoF = linear ssymMtoF'
-    where ssymMtoF' (SSymM u) = sumv [mu (set,po) u v *> return (SSymF v) | v <- set, po u v]
-              where set = L.permutations u
-                    po = weakOrder
+ssymMtoF = linear ssymMtoF' where
+    ssymMtoF' (SSymM u) = sumv [mu (set,po) u v *> return (SSymF v) | v <- set, po u v]
+        where set = L.permutations u
+              po = weakOrder
 
 -- |Convert an element of SSym represented in the fundamental basis to the monomial basis
 ssymFtoM :: (Eq k, Num k) => Vect k SSymF -> Vect k SSymM
-ssymFtoM = linear ssymFtoM'
-    where ssymFtoM' (SSymF u) = sumv [return (SSymM v) | v <- set, po u v]
-              where set = L.permutations u
-                    po = weakOrder
+ssymFtoM = linear ssymFtoM' where
+    ssymFtoM' (SSymF u) = sumv [return (SSymM v) | v <- set, po u v]
+        where set = L.permutations u
+              po = weakOrder
 
 -- (p,q)-shuffles: permutations of [1..p+q] having at most one descent, at position p
 -- denoted S^{(p,q)} in Aguiar&Sottile
@@ -457,16 +459,16 @@
 
 -- |Convert an element of YSym represented in the monomial basis to the fundamental basis
 ysymMtoF :: (Eq k, Num k) => Vect k YSymM -> Vect k (YSymF ())
-ysymMtoF = linear ysymMtoF'
-    where ysymMtoF' (YSymM t) = sumv [mu (set,po) t s *> return (YSymF s) | s <- set]
-              where po = tamariOrder
-                    set = tamariUpSet t -- [s | s <- trees (nodecount t), t `tamariOrder` s]
+ysymMtoF = linear ysymMtoF' where
+    ysymMtoF' (YSymM t) = sumv [mu (set,po) t s *> return (YSymF s) | s <- set]
+        where po = tamariOrder
+              set = tamariUpSet t -- [s | s <- trees (nodecount t), t `tamariOrder` s]
 
 -- |Convert an element of YSym represented in the fundamental basis to the monomial basis
 ysymFtoM :: (Eq k, Num k) => Vect k (YSymF ()) -> Vect k YSymM
-ysymFtoM = linear ysymFtoM'
-    where ysymFtoM' (YSymF t) = sumv [return (YSymM s) | s <- tamariUpSet t]
-                            -- sumv [return (YSymM s) | s <- trees (nodecount t), t `tamariOrder` s]
+ysymFtoM = linear ysymFtoM' where
+    ysymFtoM' (YSymF t) = sumv [return (YSymM s) | s <- tamariUpSet t]
+                       -- sumv [return (YSymM s) | s <- trees (nodecount t), t `tamariOrder` s]
 
 
 instance (Eq k, Num k) => Algebra k YSymM where
@@ -565,13 +567,13 @@
 
 -- |Convert an element of QSym represented in the monomial basis to the fundamental basis
 qsymMtoF :: (Eq k, Num k) => Vect k QSymM -> Vect k QSymF
-qsymMtoF = linear qsymMtoF'
-    where qsymMtoF' (QSymM alpha) = sumv [(-1) ^ (length beta - length alpha) *> return (QSymF beta) | beta <- refinements alpha]
+qsymMtoF = linear qsymMtoF' where
+    qsymMtoF' (QSymM alpha) = sumv [(-1) ^ (length beta - length alpha) *> return (QSymF beta) | beta <- refinements alpha]
 
 -- |Convert an element of QSym represented in the fundamental basis to the monomial basis
 qsymFtoM :: (Eq k, Num k) => Vect k QSymF -> Vect k QSymM
-qsymFtoM = linear qsymFtoM'
-    where qsymFtoM' (QSymF alpha) = sumv [return (QSymM beta) | beta <- refinements alpha] -- ie beta <- up-set of alpha
+qsymFtoM = linear qsymFtoM' where
+    qsymFtoM' (QSymF alpha) = sumv [return (QSymM beta) | beta <- refinements alpha] -- ie beta <- up-set of alpha
 
 instance (Eq k, Num k) => Algebra k QSymF where
     unit x = x *> return (QSymF [])
@@ -592,13 +594,13 @@
 -- the above induces Hopf algebra structure on quasi-symmetric functions via
 -- m_alpha -> sum [product (zipWith (^) (map x_ is) alpha | is <- combinationsOf k [] ] where k = length alpha
 
-xvars n = [glexvar ("x" ++ show i) | i <- [1..n] ]
+-- xvars n = [glexvar ("x" ++ show i) | i <- [1..n] ]
 
 -- |@qsymPoly n is@ is the quasi-symmetric polynomial in n variables for the indices is. (This corresponds to the
 -- monomial basis for QSym.) For example, qsymPoly 3 [2,1] == x1^2*x2+x1^2*x3+x2^2*x3.
 qsymPoly :: Int -> [Int] -> GlexPoly Q String
 qsymPoly n is = sum [product (zipWith (^) xs' is) | xs' <- combinationsOf r xs]
-    where xs = xvars n
+    where xs = [glexvar ("x" ++ show i) | i <- [1..n] ]
           r = length is
 
 
@@ -776,6 +778,7 @@
 
 -- "inverse" for descendingTree
 -- These are the maps called gamma in Loday.pdf
+-- or are they? - these give the min and max inverse images in the lexicographic order, rather than the weak order?
 minPerm t = minPerm' (lrCountTree t)
     where minPerm' E = []
           minPerm' (T l (lc,rc) r) = minPerm' l ++ [lc+rc+1] ++ map (+lc) (minPerm' r)
@@ -810,10 +813,9 @@
 
 -- The composition of [1..n] obtained by treating each descent as a cut
 descentComposition [] = []
-descentComposition xs = dc $ zipWith (>) xs (tail xs) ++ [False]
-    where dc bs = case break id bs of
-                  (ls,r:rs) -> (length ls + 1) : dc rs
-                  (ls,[]) -> [length ls]
+descentComposition xs = descComp 0 xs where
+    descComp c (x1:x2:xs) = if x1 < x2 then descComp (c+1) (x2:xs) else (c+1) : descComp 0 (x2:xs)
+    descComp c [x] = [c+1]
 
 -- |Given a permutation of [1..n], its descents are those positions where the next number is less than the previous number.
 -- For example, the permutation [2,3,5,1,6,4] has descents from 5 to 1 and from 6 to 4. The descents can be regarded as cutting
diff --git a/Math/Combinatorics/Digraph.hs b/Math/Combinatorics/Digraph.hs
--- a/Math/Combinatorics/Digraph.hs
+++ b/Math/Combinatorics/Digraph.hs
@@ -11,7 +11,7 @@
 import qualified Data.Map as M
 import qualified Data.Set as S
 
-toSet = S.toList . S.fromList
+import Math.Core.Utils (picks, toSet)
 
 -- |A digraph is represented as DG vs es, where vs is the list of vertices, and es is the list of edges.
 -- Edges are directed: an edge (u,v) means an edge from u to v.
@@ -50,8 +50,8 @@
     | otherwise = digraphIsos' [] vsa vsb
     where digraphIsos' xys [] [] = [xys]
           digraphIsos' xys (x:xs) ys =
-              concat [ digraphIsos' ((x,y):xys) xs (L.delete y ys)
-                     | y <- ys, isCompatible (x,y) xys]
+              concat [ digraphIsos' ((x,y):xys) xs ys'
+                     | (y,ys') <- picks ys, isCompatible (x,y) xys]
           isCompatible (x,y) xys = and [ ((x,x') `elem` esa) == ((y,y') `elem` esb)
                                       && ((x',x) `elem` esa) == ((y',y) `elem` esb)
                                        | (x',y') <- xys ]
@@ -76,8 +76,8 @@
                                          (x',y') <- xys]
           dfs xys [] [] = [xys]
           dfs xys (x:xs) ys =
-              concat [ dfs ((x,y):xys) xs (L.delete y ys)
-                     | y <- ys, isCompatible (x,y) xys]
+              concat [ dfs ((x,y):xys) xs ys'
+                     | (y,ys') <- picks ys, isCompatible (x,y) xys]
 
 -- For DAGs, can almost certainly do better than the above by using the height partition
 -- However see remarks in Poset on orderIsos:
@@ -111,8 +111,8 @@
           dfs xys [] [] = [xys]
           dfs xys ([]:las) ([]:lbs) = dfs xys las lbs
           dfs xys ((x:xs):las) (ys:lbs) =
-              concat [ dfs ((x,y):xys) (xs:las) (L.delete y ys : lbs)
-                     | y <- ys, isCompatible (x,y) xys]
+              concat [ dfs ((x,y):xys) (xs:las) (ys' : lbs)
+                     | (y,ys') <- picks ys, isCompatible (x,y) xys]
           isCompatible (x,y) xys =
               let preds_x = M.findWithDefault [] x predsA
                   preds_y = M.findWithDefault [] y predsB
@@ -169,7 +169,7 @@
           dfs xys [] [] = [xys]
           dfs xys ([]:sls) ([]:tls) = dfs xys sls tls
           dfs xys ((x:xs):sls) (ys:tls) =
-              concat [ dfs ((x,y):xys) (xs:sls) (L.delete y ys : tls) | y <- ys]
+              concat [ dfs ((x,y):xys) (xs:sls) (ys' : tls) | (y,ys') <- picks ys]
               -- not applying any compatibility condition yet
 
 
diff --git a/Math/Combinatorics/FiniteGeometry.hs b/Math/Combinatorics/FiniteGeometry.hs
--- a/Math/Combinatorics/FiniteGeometry.hs
+++ b/Math/Combinatorics/FiniteGeometry.hs
@@ -1,9 +1,11 @@
--- Copyright (c) David Amos, 2008-2011. All rights reserved.
+-- Copyright (c) David Amos, 2008-2015. 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.
 module Math.Combinatorics.FiniteGeometry where
 
+import Prelude hiding ( (*>) )
+
 import Data.List as L
 import qualified Data.Set as S
 
@@ -189,7 +191,7 @@
     es = L.sort [ [Left x, Right b] | b <- lines, x <- closurePG b]
 -- Could also consider incidence structure between points and planes, etc
 
--- incidenceAuts (incidenceGraphPG n fq) == PGL(n,fq) * auts fq
+-- incidenceAuts (incidenceGraphPG n fq) == PGL(n+1,fq) * auts fq
 -- For example, incidenceAuts (incidenceGraphPG 2 f4) =
 -- PGL(3,f4) * auts f4
 -- where PGL(3,f4)/PSL(3,f4) == f4* (multiplicative group of f4),
diff --git a/Math/Combinatorics/Graph.hs b/Math/Combinatorics/Graph.hs
--- a/Math/Combinatorics/Graph.hs
+++ b/Math/Combinatorics/Graph.hs
@@ -14,7 +14,7 @@
 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
+import qualified Math.Algebra.Group.SchreierSims as SS
 
 -- Main source: Godsil & Royle, Algebraic Graph Theory
 
@@ -105,7 +105,7 @@
 
 -- |c n is the cyclic graph on n vertices
 c :: (Integral t) => t -> Graph t
-c n = graph (vs,es) where
+c n | n >= 3 = graph (vs,es) where
     vs = [1..n]
     es = L.insert [1,n] [[i,i+1] | i <- [1..n-1]]
 -- automorphism group is D2n
@@ -165,7 +165,12 @@
            [2,7],[7,3],[3,8],[8,4],[4,9],[9,5],[5,10],[10,6],[6,11],[11,2] ]
 
 
+-- Prisms are regular, vertex-transitive, but not edge-transitive unless n == 1, 2, 4.
+-- (prism 2 ~= q 2, prism 4 ~= q 3)
+prism :: Int -> Graph (Int,Int)
+prism n = k 2 `cartProd` c n
 
+
 -- convert a graph to have [1..n] as vertices
 to1n (G vs es) = graph (vs',es') where
     mapping = M.fromList $ zip vs [1..] -- the mapping from vs to [1..n]
@@ -223,11 +228,22 @@
 lineGraph' (G vs es) = graph (es, [ [ei,ej] | ei <- es, ej <- dropWhile (<= ei) es, ei `intersect` ej /= [] ])
 
 
+-- For example cartProd (c m) (c n) is a wireframe for a torus
+-- cartProd (q m) (q n) `isGraphIso` q (m+n)
+-- Godsil and Royle p154
+cartProd (G vs es) (G vs' es') = G us [e | e@[u,u'] <- combinationsOf 2 us, u `adj` u' ]
+    where us = [(v,v') | v <- vs, v' <- vs']
+          eset = S.fromList es
+          eset' = S.fromList es'
+          adj (x1,y1) (x2,y2) = x1 == x2 && L.sort [y1,y2] `S.member` eset'
+                             || y1 == y2 && L.sort [x1,x2] `S.member` eset
+
+
 -- SIMPLE PROPERTIES OF GRAPHS
 
-order g = length (vertices g)
+order = length . vertices
 
-size g = length (edges g)
+size = length . edges
 
 -- also called degree
 valency (G vs es) v = length $ filter (v `elem`) es
@@ -297,16 +313,24 @@
 
 -- circumference = max cycle - Bollobas p104
 
+-- Vertices that are not in the same component as the start vertex all go into a final cell
+distancePartition g@(G vs es) v = distancePartitionS vs (S.fromList es) v
 
-distancePartition g v = distancePartition' S.empty (S.singleton v) where
-    distancePartition' interior boundary
-        | S.null boundary = []
-        | otherwise = let interior' = S.union interior boundary
-                          boundary' = foldl S.union S.empty [S.fromList (nbrs g x) | x <- S.toList boundary] S.\\ interior'
-                      in S.toList boundary : distancePartition' interior' boundary'
+distancePartitionS vs eset v = distancePartition' (S.singleton v) (S.delete v (S.fromList vs)) where
+    distancePartition' boundary exterior
+        | S.null boundary = if S.null exterior then [] else [S.toList exterior] -- graph may not have been connected
+        | otherwise = let (boundary', exterior') = S.partition (\v -> any (`S.member` eset) [L.sort [u,v] | u <- S.toList boundary]) exterior
+                      in S.toList boundary : distancePartition' boundary' exterior'
 
 -- the connected component to which v belongs
-component g v = L.sort $ concat $ distancePartition g v
+component g v = component' S.empty (S.singleton v) where
+    component' interior boundary
+        | S.null boundary = S.toList interior
+        | otherwise = let interior' = S.union interior boundary
+                          boundary' = foldl S.union S.empty [S.fromList (nbrs g x) | x <- S.toList boundary] S.\\ interior'
+                      in component' interior' boundary'
+-- TODO: This can almost certainly be made more efficient.
+-- nbrs is O(n), and this calls it for each vertex in the component, so it is O(n^2)
 
 -- |Is the graph connected?
 isConnected :: (Ord t) => Graph t -> Bool
@@ -317,8 +341,8 @@
     where components' [] = []
           components' (v:vs) = let c = component g v in c : components' (vs LS.\\ c)
 
--- MORE GRAPHS
 
+-- MORE GRAPHS
 
 -- Generalized Johnson graph, Godsil & Royle p9
 -- Also called generalised Kneser graph, http://en.wikipedia.org/wiki/Kneser_graph
@@ -356,7 +380,7 @@
       ++ (map . map) Right [ [i, (i+k) `mod` n] | i <- [0..n-1] ]
 
 petersen2 = gp 5 2
-prism n = gp n 1
+prism' n = gp n 1
 durer = gp 6 2
 mobiusKantor = gp 8 3
 dodecahedron2 = gp 10 2
diff --git a/Math/Combinatorics/GraphAuts.hs b/Math/Combinatorics/GraphAuts.hs
--- a/Math/Combinatorics/GraphAuts.hs
+++ b/Math/Combinatorics/GraphAuts.hs
@@ -1,20 +1,26 @@
--- Copyright (c) David Amos, 2009. All rights reserved.
+-- Copyright (c) David Amos, 2009-2014. All rights reserved.
 
+{-# LANGUAGE NoMonomorphismRestriction, TupleSections, DeriveFunctor #-}
+
 module Math.Combinatorics.GraphAuts (isVertexTransitive, isEdgeTransitive,
-                                     isArcTransitive, is2ArcTransitive, is3ArcTransitive, isnArcTransitive,
+                                     isArcTransitive, is2ArcTransitive, is3ArcTransitive, is4ArcTransitive, isnArcTransitive,
                                      isDistanceTransitive,
-                                     graphAuts, incidenceAuts,
+                                     graphAuts, incidenceAuts, graphAuts7, graphAuts8, incidenceAuts2,
+                                     isGraphAut, isIncidenceAut,
                                      graphIsos, incidenceIsos,
                                      isGraphIso, isIncidenceIso) where
 
-import Data.Either (lefts)
+import Data.Either (lefts, rights, partitionEithers)
 import qualified Data.List as L
 import qualified Data.Map as M
 import qualified Data.Set as S
 import Data.Maybe
+import Data.Ord (comparing)
+import qualified Data.Foldable as Foldable
+import qualified Data.Sequence as Seq
 
 import Math.Common.ListSet
-import Math.Core.Utils (combinationsOf, pairs)
+import Math.Core.Utils (combinationsOf, intersectAsc, pairs, picks, (^-))
 import Math.Combinatorics.Graph
 -- import Math.Combinatorics.StronglyRegularGraph
 -- import Math.Combinatorics.Hypergraph -- can't import this, creates circular dependency
@@ -55,7 +61,8 @@
     orbitP auts v == v:vs && -- isVertexTransitive g
     orbitP stab n == n:ns
     where auts = graphAuts g
-          stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
+          stab = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts
+          -- stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
           n:ns = nbrs g v
 
 -- execution time of both of the above is dominated by the time to calculate the graph auts, so their performance is similar
@@ -77,7 +84,7 @@
 
 -- note that a graph with triangles can't be 3-arc transitive, etc, because an aut can't map a self-crossing arc to a non-self-crossing arc
 
--- |A graph is n-arc-transitive is its automorphism group is transitive on n-arcs. (An n-arc is an ordered sequence (v0,...,vn) of adjacent vertices, with crossings allowed but not doubling back.)
+-- |A graph is n-arc-transitive if its automorphism group is transitive on n-arcs. (An n-arc is an ordered sequence (v0,...,vn) of adjacent vertices, with crossings allowed but not doubling back.)
 isnArcTransitive :: (Ord t) => Int -> Graph t -> Bool
 isnArcTransitive _ (G [] []) = True
 isnArcTransitive n g@(G (v:vs) es) =
@@ -85,7 +92,8 @@
     orbit (->^) a stab == a:as
     -- closure [a] [ ->^ h | h <- stab] == a:as
     where auts = graphAuts g
-          stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
+          stab = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts
+          -- stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
           a:as = findArcs g v n
 
 is2ArcTransitive :: (Ord t) => Graph t -> Bool
@@ -94,6 +102,10 @@
 is3ArcTransitive :: (Ord t) => Graph t -> Bool
 is3ArcTransitive g = isnArcTransitive 3 g
 
+-- The incidence graphs of the projective planes PG(2,Fq) are 4-arc-transitive
+is4ArcTransitive :: (Ord t) => Graph t -> Bool
+is4ArcTransitive g = isnArcTransitive 4 g
+
 -- Godsil & Royle 66-7
 -- |A graph is distance transitive if given any two ordered pairs of vertices (u,u') and (v,v') with d(u,u') == d(v,v'),
 -- there is an automorphism of the graph that takes (u,u') to (v,v')
@@ -103,34 +115,30 @@
     | isConnected g =
         orbitP auts v == v:vs && -- isVertexTransitive g
         length stabOrbits == diameter g + 1 -- the orbits under the stabiliser of v coincide with the distance partition from v
-    | otherwise = error "isDistanceTransitive: only defined for connected graphs"
+    | otherwise = error "isDistanceTransitive: only implemented for connected graphs"
     where auts = graphAuts g
-          stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
+          stab = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts
+          -- stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
           stabOrbits = let os = orbits stab in os ++ map (:[]) ((v:vs) L.\\ concat os) -- include fixed point orbits
 
 
 -- GRAPH AUTOMORPHISMS
 
--- !! Note, in the literature the following is just called the intersection of two partitions
--- !! Refinement actually refers to the process of refining to an equitable partition
-
--- refine one partition by another
-refine p1 p2 = filter (not . null) $ refine' p1 p2
--- Refinement preserves ordering within cells but not between cells
--- eg the cell [1,2,3,4] could be refined to [2,4],[1,3]
-
--- refine, but leaving null cells in
--- we use this in the graphAuts functions when comparing two refinements to check that they split in the same way
-refine' p1 p2 = concat [ [c1 `intersect` c2 | c2 <- p2] | c1 <- p1]
-
-
+-- |Is the permutation an automorphism of the graph?
+isGraphAut :: Ord t => Graph t -> Permutation t -> Bool
 isGraphAut (G vs es) h = all (`S.member` es') [e -^ h | e <- es]
     where es' = S.fromList es
 -- this works best on sparse graphs, where p(edge) < 1/2
 -- if p(edge) > 1/2, it would be better to test on the complement of the graph
 
-
-
+-- |Is the permutation an automorphism of the incidence structure represented by the graph?
+-- (Note that an incidence graph colours points as Left, blocks as Right, and a permutation
+-- that swaps points and blocks, even if it is an automorphism of the graph, does not represent
+-- an automorphism of the incidence structure. Instead, a point-block crossover is called a duality.)
+isIncidenceAut :: (Ord p, Ord b) => Graph (Either p b) -> Permutation (Either p b) -> Bool
+isIncidenceAut (G vs es) h = all (`S.member` es') [e ->^ h | e <- es]
+    -- using ->^ instead of -^ excludes dualities, since each edge is of the form [Left p, Right b]
+    where es' = S.fromList es
 
 -- Calculate a map consisting of neighbour lists for each vertex in the graph
 -- If a vertex has no neighbours then it is left out of the map
@@ -143,260 +151,482 @@
 -- ALTERNATIVE VERSIONS OF GRAPH AUTS
 -- (showing how we got to the final version)
 
--- return all graph automorphisms, using naive depth first search
-graphAuts1 (G vs es) = dfs [] vs vs
-    where dfs xys (x:xs) ys =
-              concat [dfs ((x,y):xys) xs (L.delete y ys) | y <- ys, isCompatible (x,y) xys]
-          dfs xys [] [] = [fromPairs xys]
-          isCompatible (x,y) xys = and [([x',x] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x',y') <- xys]
-          es' = S.fromList es
+data SearchTree a = T Bool a [SearchTree a] deriving (Eq, Ord, Show, Functor)
+-- The boolean indicates whether or not this is a terminal / solution node
 
--- return generators for graph automorphisms
--- (using Lemma 9.1.1 from Seress p203 to prune the search tree)
-graphAuts2 (G vs es) = graphAuts' [] vs
-    where graphAuts' us (v:vs) =
-              let uus = zip us us
-              in concat [take 1 $ dfs ((v,w):uus) vs (v : L.delete w vs) | w <- vs, isCompatible (v,w) uus]
-              ++ graphAuts' (v:us) vs
-              -- stab us == transversal for stab (v:us) ++ stab (v:us)  (generators thereof)
-          graphAuts' _ [] = [] -- we're not interested in finding the identity element
-          dfs xys (x:xs) ys =
-              concat [dfs ((x,y):xys) xs (L.delete y ys) | y <- ys, isCompatible (x,y) xys]
-          dfs xys [] [] = [fromPairs xys]
-          isCompatible (x,y) xys = and [([x',x] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x',y') <- xys]
-          es' = S.fromList es
+leftDepth (T _ _ []) = 1
+leftDepth (T _ _ (t:ts)) = 1 + leftDepth t
 
--- Now using distance partitions
--- Note that because of the use of distance partitions, this is only valid for connected graphs
-graphAuts3 g@(G vs es) = graphAuts' [] [vs] where
-    graphAuts' us ((x:ys):pt) =
-        let px = refine' (ys : pt) (dps M.! x)
-            p y = refine' ((x : L.delete y ys) : pt) (dps M.! y)
-            uus = zip us us
-            p' = L.sort $ filter (not . null) $ px
-        in concat [take 1 $ dfs ((x,y):uus) px (p y) | y <- ys]
-        ++ graphAuts' (x:us) p'
-    graphAuts' us ([]:pt) = graphAuts' us pt
-    graphAuts' _ [] = []
-    dfs xys p1 p2
-        | map length p1 /= map length p2 = []
-        | otherwise =
-            let p1' = filter (not . null) p1
-                p2' = filter (not . null) p2
-            in if all isSingleton p1'
-               then let xys' = xys ++ zip (concat p1') (concat p2')
-                    in if isCompatible xys' then [fromPairs' xys'] else []
-                    -- we shortcut the search when we have all singletons, so must check isCompatible to avoid false positives
-               else let (x:xs):p1'' = p1'
-                        ys:p2'' = p2'
-                    in concat [dfs ((x,y):xys)
-                                   (refine' (xs : p1'') (dps M.! x))
-                                   (refine' ((L.delete y ys):p2'') (dps M.! y))
-                                   | y <- ys]
-    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
-    dps = M.fromList [(v, distancePartition g v) | v <- vs]
+leftWidths (T _ _ []) = []
+leftWidths (T _ _ ts@(t:_)) = length ts : leftWidths t
+
+graphAutsEdgeSearchTree (G vs es) = dfs [] vs vs where
+    dfs xys (x:xs) yys = T False xys [dfs ((x,y):xys) xs ys | (y,ys) <- picks yys, isCompatible xys (x,y)]
+    dfs xys [] [] = T True xys []
+    isCompatible xys (x',y') = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys]
     es' = S.fromList es
 
+graphAuts1 = map fromPairs . terminals . graphAutsEdgeSearchTree
+
+terminals (T False _ ts) = concatMap terminals ts
+terminals (T True xys _) = [xys]
+
+-- Using Lemma 9.1.1 from Seress p203 to prune the search tree
+-- Because auts form a group, it is sufficient to expand only each leftmost branch of the tree in full.
+-- For every other branch, it is sufficient to find a single representative, since the other elements
+-- can then be obtained by multiplication in the group (using the leftmost elements).
+-- In effect, we are finding a transversal generating set.
+-- Note however, that this transversal generating set is relative to whatever base order the tree uses,
+-- so for clarity, the tree should use natural vertex order.
+transversalTerminals (T False _ (t:ts)) = concatMap (take 1 . transversalTerminals) ts ++ transversalTerminals t
+-- transversalTerminals (T False _ (t:ts)) = transversalTerminals t ++ concatMap (take 1 . transversalTerminals) ts
+transversalTerminals (T True xys _) = [xys]
+transversalTerminals _ = []
+
+graphAuts2 = filter (/=1) . map fromPairs . transversalTerminals . graphAutsEdgeSearchTree
+-- init because last is identity
+
 isSingleton [_] = True
 isSingleton _ = False
 
+intersectCells p1 p2 = concat [ [c1 `intersectAsc` c2 | c2 <- p2] | c1 <- p1]
+-- Intersection preserves ordering within cells but not between cells
+-- eg the cell [1,2,3,4] could be refined to [2,4],[1,3]
 
--- Now we try to use generators we've already found at a given level to save us having to look for others
--- For example, if we have found (1 2)(3 4) and (1 3 2), then we don't need to look for something taking 1 -> 4
-graphAuts4 g@(G vs es) = graphAuts' [] [vs] where
-    graphAuts' us p@((x:ys):pt) =
-        -- let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
-        let p' = L.sort $ refine (ys:pt) (dps M.! x)
-        in level us p x ys []
-        ++ graphAuts' (x:us) p'
-    graphAuts' us ([]:pt) = graphAuts' us pt
-    graphAuts' _ [] = []
-    level us p@(ph:pt) x (y:ys) hs =
-        let px = refine' (L.delete x ph : pt) (dps M.! x)
-            py = refine' (L.delete y ph : pt) (dps M.! y)
-            uus = zip us us
-        in case dfs ((x,y):uus) px py of
-           []  -> level us p x ys hs
-           h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
-    level _ _ _ [] _ = []
-    dfs xys p1 p2
-        | map length p1 /= map length p2 = []
-        | otherwise =
-            let p1' = filter (not . null) p1
-                p2' = filter (not . null) p2
-            in if all isSingleton p1'
-               then let xys' = xys ++ zip (concat p1') (concat p2')
-                    in if isCompatible xys' then [fromPairs' xys'] else []
-               else let (x:xs):p1'' = p1'
-                        ys:p2'' = p2'
-                    in concat [dfs ((x,y):xys)
-                                   (refine' (xs : p1'') (dps M.! x))
-                                   (refine' ((L.delete y ys):p2'') (dps M.! y))
-                                   | y <- ys]
+
+graphAutsDistancePartitionSearchTree g@(G vs es) = dfs [] ([vs],[vs]) where
+    dfs xys (srcPart,trgPart)
+        | all isSingleton srcPart =
+             let xys' = zip (concat srcPart) (concat trgPart)
+             in T (isCompatible xys') (xys++xys') []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x:xs):srcCells = srcPart
+                          yys   :trgCells = trgPart
+                          srcPart' = intersectCells (xs : srcCells) (dps M.! x)
+                      in T False xys -- the L.sort in the following line is so that we traverse vertices in natural order
+                         [dfs ((x,y):xys) ((unzip . L.sort) (zip (filter (not . null) srcPart') (filter (not . null) trgPart')))
+                         | (y,ys) <- picks yys,
+                           let trgPart' = intersectCells (ys : trgCells) (dps M.! y),
+                           map length srcPart' == map length trgPart']
     isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
-    dps = M.fromList [(v, distancePartition g v) | v <- vs]
     es' = S.fromList es
+    dps = M.fromAscList [(v, distancePartitionS vs es' v) | v <- vs]
 
--- contrary to first thought, you can't stop when a level is null - eg kb 2 3, the third level is null, but the fourth isn't
+graphAuts3 = filter (/=1) . map fromPairs . transversalTerminals . graphAutsDistancePartitionSearchTree
 
+-- Whereas transversalTerminals produced a transversal generating set, here we produce a strong generating set.
+-- In particular, if we have already found (3 4), and then we find (1 2 3),
+-- then there is no need to look for (1 3 ...) or (1 4 ...), since it is clear that such elements exist
+-- as products of those we have already found.
+strongTerminals = strongTerminals' [] where
+    strongTerminals' gs (T False xys ts) =
+        case listToMaybe $ reverse $ filter (\(x,y) -> x /= y) xys of -- the first vertex that isn't fixed
+        Nothing -> L.foldl' (\hs t -> strongTerminals' hs t) gs ts
+        Just (x,y) -> if y `elem` (x .^^ gs)
+                      then gs
+                      -- Since we're not on the leftmost spine, we can stop as soon as we find one new element
+                      else find1New gs ts
+                      -- else L.foldl' (\hs t -> if hs /= gs then hs else strongTerminals' hs t) gs ts
+    strongTerminals' gs (T True xys []) = fromPairs xys : gs
+    find1New gs (t:ts) = let hs = strongTerminals' gs t
+                         in if take 1 gs /= take 1 hs -- we know a new element would be placed at the front
+                            then hs
+                            else find1New gs ts
+    find1New gs [] = gs
 
+-- |Given a graph g, @graphAuts g@ returns a strong generating set for the automorphism group of g.
+graphAuts :: (Ord a) => Graph a -> [Permutation a]
+graphAuts = filter (/=1) . strongTerminals . graphAutsDistancePartitionSearchTree
 
--- an example for equitable partitions
--- this is a graph whose distance partition (from any vertex) can be refined to an equitable partition
-eqgraph = G vs es where
-    vs = [1..14]
-    es = L.sort $ [[1,14],[2,13]] ++ [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, v1+1 == v2 || v1+3 == v2 && even v2]
 
--- refine a partition to give an equitable partition
-toEquitable g cells = L.sort $ toEquitable' [] cells where
-    toEquitable' ls (r:rs) =
-        let (lls,lrs) = L.partition isSingleton $ map (splitNumNbrs r) ls
-            -- so the lrs split, and the lls didn't
-            rs' = concatMap (splitNumNbrs r) rs
-        in if isSingleton r -- then we know it won't split further, so can remove it from further processing
-           then r : toEquitable' (concat lls) (concat lrs ++ rs')
-           else toEquitable' (r : concat lls) (concat lrs ++ rs')
-    toEquitable' ls [] = ls
-    splitNumNbrs t c = map (map snd) $ L.groupBy (\x y -> fst x == fst y) $ L.sort
-                    [ (length ((nbrs_g M.! v) `intersect` t), v) | v <- c]
-    nbrs_g = M.fromList [(v, nbrs g v) | v <- vertices g]
+-- Using colourings (M.Map vertex colour, M.Map colour [vertex]), in place of partitions ([[vertex]])
+-- This turns out to be slower than using partitions.
+-- Updating the colour partition incrementally seems to be much less efficient than just recalculating it each time
+-- (Recalculating each time is O(n), incrementally updating is O(n^2)?)
+graphAutsDistanceColouringSearchTree g@(G vs es) = dfs [] unitCol unitCol where
+    unitCol = (M.fromList $ map (,[]) vs, M.singleton [] vs) -- "unit colouring"
+    dfs xys srcColouring@(srcVmap,srcCmap) trgColouring@(trgVmap,trgCmap)
+        -- ( | M.map length srcCmap /= M.map length trgCmap = T False xys [] )
+        | all isSingleton (M.elems srcCmap) = -- discrete colouring
+             let xys' = zip (concat $ M.elems srcCmap) (concat $ M.elems trgCmap)
+             in T (isCompatible xys') (reverse xys'++xys) []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x,c) = M.findMin srcVmap
+                          (xVmap,xCmap) = dcs M.! x
+                          ys = trgCmap M.! c
+                          srcVmap' = M.delete x (intersectColouring srcVmap xVmap)
+                          srcCmap' = colourPartition srcVmap'
+                          -- srcCmap' = M.fromAscList [(c1++c2, cell) | (c1,srcCell) <- M.assocs srcCmap, (c2,xCell) <- M.assocs xCmap,
+                          --                                            let cell = L.delete x (intersectAsc srcCell xCell),
+                          --                                            (not . null) cell]
+                      in T False xys
+                         [dfs ((x,y):xys) (srcVmap',srcCmap') (trgVmap',trgCmap')
+                         | y <- ys,
+                           let (yVmap,yCmap) = dcs M.! y,
+                           let trgVmap' = M.delete y (intersectColouring trgVmap yVmap),
+                           let trgCmap' = colourPartition trgVmap',
+                           -- let trgCmap' = M.fromAscList [(c1++c2, cell) | (c1,trgCell) <- M.assocs trgCmap, (c2,yCell) <- M.assocs yCmap,
+                           --                                                let cell = L.delete y (intersectAsc trgCell yCell),
+                           --                                                (not . null) cell],
+                           M.map length srcCmap' == M.map length trgCmap' ]
+    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+    es' = S.fromList es
+    dcs = M.fromAscList [(v, distanceColouring v) | v <- vs]
+    distanceColouring u = let dp = distancePartitionS vs es' u
+                              vmap = M.fromList [(v,[c]) | (cell,c) <- zip dp [0..], v <- cell]
+                              cmap = M.fromList $ zip (map (:[]) [0..]) dp
+                          in (vmap, cmap)
 
+{-
+-- If we are going to recalculate the colour partition each time anyway,
+-- then we don't need to carry it around, and can simplify the code
+graphAutsDistanceColouringSearchTree g@(G vs es) = dfs [] initCol initCol where
+    initCol = M.fromList $ map (,[]) vs
+    dfs xys srcCol trgCol
+        | M.map length srcPart /= M.map length trgPart = T False xys []
+        | all isSingleton (M.elems srcPart) =
+             let xys' = zip (concat $ M.elems srcPart) (concat $ M.elems trgPart)
+             in T (isCompatible xys') (reverse xys'++xys) []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x,c) = M.findMin srcCol
+                          ys = trgPart M.! c
+                          srcCol' = M.delete x $ intersectColouring srcCol (dcs M.! x)
+                      in T False xys
+                         [dfs ((x,y):xys) srcCol' trgCol'
+                         | y <- ys,
+                           let trgCol' = M.delete y (intersectColouring trgCol (dcs M.! y))]
+        where srcPart = colourPartition srcCol
+              trgPart = colourPartition trgCol
+    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+    es' = S.fromList es
+    dcs = M.fromAscList [(v, distanceColouring v) | v <- vs]
+    distanceColouring u = M.fromList [(v,[c]) | (cell,c) <- zip (distancePartitionS vs es' u) [0..], v <- cell]
+-}
+distanceColouring (G vs es) u = M.fromList [(v,[c]) | (cell,c) <- zip (distancePartitionS vs es' u) [0..], v <- cell]
+    where es' = S.fromList es
 
--- try to refine two partitions in parallel, failing if they become mismatched
-toEquitable2 nbrs_g psrc ptrg = unzip $ L.sort $ toEquitable' [] (zip psrc ptrg) where
-    toEquitable' ls (r:rs) =
-        let ls' = map (splitNumNbrs nbrs_g r) ls
-            (lls,lrs) = L.partition isSingleton $ map fromJust ls'
-            rs' = map (splitNumNbrs nbrs_g r) rs
-        in if any isNothing ls' || any isNothing rs'
-           then []
-           else
-               {- if (isSingleton . fst) r
-               then r : toEquitable' (concat lls) (concat lrs ++ concatMap fromJust rs')
-               else -} toEquitable' (r : concat lls) (concat lrs ++ concatMap fromJust rs')
-    toEquitable' ls [] = ls
+intersectColouring c1 c2 = M.intersectionWith (++) c1 c2
 
-splitNumNbrs nbrs_g (t_src,t_trg) (c_src,c_trg) =
-    let src_split = L.groupBy (\x y -> fst x == fst y) $ L.sort
-                    [ (length ((nbrs_g M.! v) `intersect` t_src), v) | v <- c_src]
-        trg_split = L.groupBy (\x y -> fst x == fst y) $ L.sort
-                    [ (length ((nbrs_g M.! v) `intersect` t_trg), v) | v <- c_trg]
-    in if map length src_split == map length trg_split
-       && map (fst . head) src_split == map (fst . head) trg_split
-       then Just $ zip (map (map snd) src_split) (map (map snd) trg_split)
-       else Nothing
-       -- else error (show (src_split, trg_split)) -- for debugging
+colourPartition c = L.foldr (\(k,v) m -> M.insertWith (++) v [k] m) M.empty (M.assocs c)
 
--- Now, every time we intersect two partitions, refine to an equitable partition
 
--- |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 = 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 "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 []
-              ++ graphAuts' (x:us) p'
-          graphAuts' us ([]:pt) = graphAuts' us pt
-          graphAuts' _ [] = []
-          level us p@(ph:pt) x (y:ys) hs =
-              let px = refine' (L.delete x ph : pt) (dps M.! x)
-                  py = refine' (L.delete y ph : pt) (dps M.! y)
-                  uus = zip us us
-              in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
-                 []  -> level us p x ys hs
-                 h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
-          level _ _ _ [] _ = []
-          dps = M.fromList [(v, distancePartition g v) | v <- vs]
-          es' = S.fromList es
-          nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
+-- Based on McKay’s Canonical Graph Labeling Algorithm, by Stephen G. Hartke and A. J. Radcliffe
+-- (http://www.math.unl.edu/~aradcliffe1/Papers/Canonical.pdf)
 
-dfsEquitable (dps,es',nbrs_g) xys p1 p2 = dfs xys p1 p2 where
-    dfs xys p1 p2
-        | map length p1 /= map length p2 = []
-        | otherwise =
-            let p1' = filter (not . null) p1
-                p2' = filter (not . null) p2
-                (p1e,p2e) = toEquitable2 nbrs_g p1' p2'
-            in if null p1e
-               then []
-               else
-                   if all isSingleton p1e
-                   then let xys' = xys ++ zip (concat p1e) (concat p2e)
-                        in if isCompatible xys' then [fromPairs' xys'] else []
-                   else let (x:xs):p1'' = p1e
-                            ys:p2'' = p2e
-                        in concat [dfs ((x,y):xys)
-                                       (refine' (xs : p1'') (dps M.! x))
-                                       (refine' ((L.delete y ys):p2'') (dps M.! y))
-                                       | y <- ys]
-    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+equitableRefinement g@(G vs es) p = equitableRefinement' (S.fromList es) p
 
+equitableRefinement' edgeset partition = go partition where
+    go cells = let splits = L.zip (L.inits cells) (L.tails cells)
+                   shatterPairs = [(L.zip ci counts,ls,rs) | (ls,ci:rs) <- splits, cj <- cells,
+                                                             let counts = map (nbrCount cj) ci, isShatter counts]
+               in case shatterPairs of -- by construction, the lexicographic least (i,j) comes first
+                  [] -> cells
+                  (vcs,ls,rs):_ -> let fragments = shatter vcs 
+                                   in go (ls ++ fragments ++ rs)
+    isShatter (c:cs) = any (/= c) cs
+    shatter vcs = map (map fst) $ L.groupBy (\x y -> snd x == snd y) $ L.sortBy (comparing snd) $ vcs
+    -- Memoizing here results in about 10% speed improvement. Not worth it for loss of generality (ie requiring HasTrie instances)
+    -- nbrCount = memo2 nbrCount'
+    -- How many neighbours in cell does vertex have
+    nbrCount cell vertex = length (filter (isEdge vertex) cell)
+    isEdge u v = L.sort [u,v] `S.member` edgeset
 
+equitablePartitionSearchTree g@(G vs es) p = dfs [] p where
+    dfs bs p = let p' = equitableRefinement' es' p in
+               if all isSingleton p'
+               then T True (p',bs) []
+               else T False (p',bs) [dfs (b:bs) p'' | (b,p'') <- splits [] p']
+    -- For now, we just split the first non-singleton cell we find
+    splits ls (r:rs) | isSingleton r = splits (r:ls) rs
+                     | otherwise = let ls' = reverse ls in [(x, ls' ++ [x]:xs:rs) | (x,xs) <- picks r]
+    es' = S.fromList es
+
+
+{-
+-- Using Data.Sequence instead of list for the partitions
+-- Makes no difference to speed (in fact slightly slower)
+equitableRefinementSeq' edgeset partition = go partition where
+    go cells = let splits = Seq.zip (Seq.inits cells) (Seq.tails cells)
+                   shatterPairs = [(L.zip ci counts,ls,rs') | (ls,rs) <- Foldable.toList splits, (not . Seq.null) rs, let ci Seq.:< rs' = Seq.viewl rs,
+                                                              cj <- Foldable.toList cells,
+                                                              let counts = map (nbrCount cj) ci, isShatter counts]
+               in case shatterPairs of -- by construction, the lexicographic least (i,j) comes first
+                  [] -> cells
+                  (vcs,ls,rs):_ -> let fragments = Seq.fromList (shatter vcs) 
+                                   in go (ls Seq.>< fragments Seq.>< rs)
+    isShatter (c:cs) = any (/= c) cs
+    shatter vcs = map (map fst) $ L.groupBy (\x y -> snd x == snd y) $ L.sortBy (comparing snd) $ vcs
+    -- How many neighbours in cell does vertex have
+    nbrCount cell vertex = length (filter (isEdge vertex) cell)
+    isEdge u v = L.sort [u,v] `S.member` edgeset
+
+equitablePartitionSeqSearchTree g@(G vs es) p = dfs [] (Seq.fromList p) where
+    dfs bs p = let p' = equitableRefinementSeq' es' p in
+               if Foldable.all isSingleton p'
+               then T True (Foldable.toList p',bs) []
+               else T False (Foldable.toList p',bs) [dfs (b:bs) p'' | (b,p'') <- splits p']
+    -- For now, we just split the first non-singleton cell we find
+    splits cells = case Seq.findIndexL (not . isSingleton) cells of
+                   Just i -> let (ls,rs) = Seq.splitAt i cells
+                                 r Seq.:< rs' = Seq.viewl rs
+                             in [(x, ls Seq.>< ([x] Seq.<| xs Seq.<| rs')) | (x,xs) <- picks r]
+                   Nothing -> error "Not possible, as we know there are non-singleton cells"
+    es' = S.fromList es
+-}
+
+-- In this version, whenever we have an equitable partition, we separate out all the singleton cells and put them to one side.
+-- (Since the partition is equitable, singleton cells have already done any work they are going to do in shattering other cells,
+-- so they will no longer play any part.)
+-- This seems to result in about 20% speedup.
+equitablePartitionSearchTree2 g@(G vs es) p = dfs [] ([],p) where
+    dfs bs (ss,cs) = let (ss',cs') = L.partition isSingleton $ equitableRefinement' es' cs
+                         ss'' = ss++ss'
+                     in case cs' of
+                        [] -> T True (ss'',bs) []
+                        -- We just split the first non-singleton cell
+                        -- c:cs'' -> T False (ss''++cs',bs) [dfs (x:bs) (ss'',[x]:xs:cs'') | (x,xs) <- picks c]
+                        c:cs'' -> T False (cs'++ss'',bs) [dfs (x:bs) (ss'',[x]:xs:cs'') | (x,xs) <- picks c]
+    es' = S.fromList es
+-- TODO: On the first level, we can use a stronger partitioning function (eg distance partitions, + see nauty manual, vertex invariants)
+
+equitableDistancePartitionSearchTree g@(G vs es) p = dfs [] p where
+    dfs bs p = let p' = equitableRefinement' es' p in
+               if all isSingleton p'
+               then T True (p',bs) []
+               else T False (p',bs) [dfs (b:bs) p'' | (b,p'') <- splits [] p']
+    -- For now, we just split the first non-singleton cell we find
+    splits ls (r:rs) | isSingleton r = splits (r:ls) rs
+                     | otherwise = [(x, p'') | let ls' = reverse ls,
+                                               (x,xs) <- picks r,
+                                               let p' = ls' ++ [x]:xs:rs,
+                                               let p'' = filter (not . null) (intersectCells p' (dps M.! x))]
+    es' = S.fromList es
+    dps = M.fromAscList [(v, distancePartitionS vs es' v) | v <- vs]
+
+
+{-
+-- This is just fmap (\(p,bs) -> (p,bs,trace p)) t
+equitablePartitionTracedSearchTree g@(G vs es) trace p = dfs [] p where
+    dfs bs p = let p' = equitableRefinement' es' p
+               in if all isSingleton p'
+                  then T True (p',bs,trace p') []
+                  else T False (p',bs,trace p') [dfs (b:bs) p'' | (b,p'') <- splits [] p']
+    -- For now, we just split the first non-singleton cell we find
+    splits ls (r:rs) | isSingleton r = splits (r:ls) rs
+                     | otherwise = let ls' = reverse ls in [(x, ls' ++ [x]:xs:rs) | (x,xs) <- picks r]
+    es' = S.fromList es
+-}
+
+-- Intended as a node invariant
+trace1 p = map (\xs@(x:_) -> (x, length xs)) $ L.group $ L.sort $ map length p
+
+equitablePartitionGraphSearchTree g@(G vs es) = equitablePartitionSearchTree g unitPartition
+    where unitPartition = [vs]
+
+-- The incidence graph has vertices that are coloured left (points) or right (blocks).
+-- We are not interested in dualities (automorphisms that swap points and blocks), so we look for colour-preserving automorphisms
+equitablePartitionIncidenceSearchTree g@(G vs es) = equitablePartitionSearchTree g lrPartition
+    where (lefts, rights) = partitionEithers vs
+          lrPartition = [map Left lefts, map Right rights]
+
+leftLeaf (T False _ (t:ts)) = leftLeaf t
+leftLeaf (T True (p,bs) []) = (concat p, reverse bs)
+{-
+leftSpine (T False x (t:ts)) = x : leftSpine t
+leftSpine (T True x []) = [x]
+-}
+allLeaves (T False _ ts) = concatMap allLeaves ts
+allLeaves (T True (p,bs) []) = [(concat p, reverse bs)]
+
+{-
+partitionTransversals tree = [fromPairs (zip canonical partition) | partition <- findTransversals tree] where
+    (_,canonical) = leftLeaf tree
+    findTransversals (T False _ (t:ts)) = concatMap (take 1 . findTransversals) ts ++ findTransversals t
+    findTransversals (T True (_,partition) []) = [concat partition]
+
+graphAuts5 = partitionTransversals . equitablePartitionGraphSearchTree
+-}
+-- NOT WORKING
+partitionBSGS0 g@(G vs es) t = (bs, findLevels t) where
+    (p1,bs) = leftLeaf t
+    g1 = fromPairs $ zip p1 vs
+    g1' = g1^-1
+    es1 = S.fromList $ edges $ fmap (.^ g1) g -- the edges of the isomorph corresponding to p1. (S.fromList makes it unnecessary to call nf.)
+    findLevels (T True (partition,_) []) = []
+    findLevels (T False (partition,_) (t:ts)) =
+        let hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+        in findLevel v hs (zip vs ts)
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut
+    -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.)
+    find1New (T True (partition,_) []) = let h = fromPairs $ zip (concat partition) vs
+                                             g' = fmap (.^ h) g
+                                         in if all (`S.member` es1) (edges g') then [h*g1'] else []
+    -- isAut h = all (`S.member` es') [e -^ h | e <- es]
+    -- es' = S.fromList es
+
+-- Given a partition search tree, return a base and strong generating set for graph automorphism group.
+partitionBSGS g@(G vs es) t = (bs, findLevels t) where
+    (canonical,bs) = leftLeaf t
+    findLevels (T True (partition,_) []) = []
+    findLevels (T False (partition,_) (t:ts)) =
+        let hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+        in findLevel v hs (zip vs ts)
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs -- TODO: Memoize this orbit
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- Some leaf nodes correspond to different isomorphs of the graph, and hence don't yield automorphisms
+    find1New (T True (partition,_) []) = let h = fromPairs $ zip canonical (concat partition)
+                                         in filter isAut [h]
+    isAut h = all (`S.member` es') [e -^ h | e <- es]
+    es' = S.fromList es
+-- The tree for g1 has leaf nodes of two different isomorphs, as does the tree for incidenceGraphPG 2 f8
+
+-- Returns auts as Right, different isomorphs as Left
+-- (Must be used with the tree which doesn't put singletons to end)
+partitionBSGS3 g@(G vs es) t = (bs, findLevels t) where
+    (p1,bs) = leftLeaf t
+    findLevels (T True (partition,_) []) = []
+    findLevels (T False (partition,_) (t:ts)) =
+        let hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+        in findLevel v hs (zip vs ts)
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ rights hs
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut
+    -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.)
+    find1New (T True (partition,_) []) = let h = fromPairs $ zip p1 (concat partition)
+                                         in if isAut h then [Right h] else [Left h]
+    isAut h = all (`S.member` es') [e -^ h | e <- es]
+    es' = S.fromList es
+-- TODO: I think we are only justified in doing find1New (ie only finding 1) if we *do* find an aut.
+-- If we don't, we should potentially keep looking in that subtree
+-- (See section 6 of paper. If we find isomorphic leaves, then the two subtrees of their common parent are isomorphic,
+-- so no need to continue searching the second.)
+
+
+-- This is using a node invariant to do more pruning.
+-- However, seems to be much slower on very regular graphs (where perhaps there is no pruning to be done)
+-- (This suggests that perhaps using fmap is not good - perhaps a space leak?)
+-- (Or perhaps it's just that calculating and comparing the node invariants is expensive)
+-- TODO: Perhaps use something simpler, like just the number of cells in the partition
+partitionBSGS2 g@(G vs es) t = (bs, findLevels t') where
+    t' = fmap (\(p,bs) -> (p,bs,trace1 p)) t
+    trace1 = length -- the number of cells in the partition
+    (canonical,bs) = leftLeaf t
+    findLevels (T True (partition,_,_) []) = []
+    findLevels (T False (partition,_,_) (t:ts)) =
+        let (T _ (_,_,trace) _) = t
+            hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+            vts = filter (\(_,T _ (_,_,trace') _) -> trace == trace') $ zip vs ts
+        in findLevel v hs vts
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut
+    -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.)
+    -- (The graph g1, below, shows a simple example where this will happen.)
+    find1New (T True (partition,_,_) []) = let h = fromPairs $ zip canonical (concat partition)
+                                           in filter isAut [h]
+    isAut h = all (`S.member` es') [e -^ h | e <- es]
+    es' = S.fromList es
+
+
+graphAuts7 g = (partitionBSGS g) (equitablePartitionGraphSearchTree g)
+
+-- This is faster on kneser graphs, but slower on incidenceGraphPG
+graphAuts8 g = (partitionBSGS g) (equitableDistancePartitionSearchTree g [vertices g])
+
+-- This is a graph where the node invariant should cause pruning.
+-- The initial equitable partition will be [[1..8],[9,10]], because it can do no better than distinguish by degree
+-- However, vertices 1..4 and vertices 5..8 are in fact different (there is no aut that takes one set to the other),
+-- so the subtrees starting 1..4 have a different invariant to those starting 5..8
+g1 = G [1..10] [[1,2],[1,3],[1,9],[2,4],[2,10],[3,4],[3,9],[4,10],[5,6],[5,8],[5,9],[6,7],[6,10],[7,8],[7,9],[8,10]]
+
+g1' = nf $ fmap (\x -> if x <= 4 then x+4 else if x <= 8 then x-4 else x) g1
+-- G [1..10] [[1,2],[1,4],[1,9],[2,3],[2,10],[3,4],[3,9],[4,10],[5,6],[5,7],[5,9],[6,8],[6,10],[7,8],[7,9],[8,10]]
+
+g2 = G [1..12] [[1,2],[1,4],[1,11],[2,3],[2,12],[3,4],[3,11],[4,12],[5,6],[5,8],[5,11],[6,9],[6,12],[7,8],[7,10],[7,11],[8,12],[9,10],[9,11],[10,12]]
+
+-- NOT WORKING: This fails to find the isomorphism between g1 and g1' above.
+-- Instead of using left leaf, we need to find the canonical isomorph, as described in the paper.
+-- (In a graph where not all leaves lead to automorphisms, we might happen to end up with non-isomorphic left leaves)
+maybeGraphIso g1 g2 = let (vs1,_) = (leftLeaf . equitablePartitionGraphSearchTree) g1
+                          (vs2,_) = (leftLeaf . equitablePartitionGraphSearchTree) g2
+                          f = M.fromList (zip vs1 vs2)
+                      in if length vs1 == length vs2 && (nf . fmap (f M.!)) g1 == g2 then Just f else Nothing
+
+
 -- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH
 
--- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
+-- This code is nearly identical to the corresponding graphAuts code, with two exceptions:
+-- 1. We start by partitioning into lefts and rights.
+-- This avoids left-right crossover auts, which while they are auts of the graph,
+-- are not auts of the incidence structure
+-- 2. When labelling the nodes, we filter out Right blocks, and unLeft the Left points
+incidenceAutsDistancePartitionSearchTree g@(G vs es) = dfs [] (lrPart, lrPart) where
+    dfs xys (srcPart,trgPart)
+        | all isSingleton srcPart =
+             let xys' = zip (concat srcPart) (concat trgPart)
+             in T (isCompatible xys') (unLeft $ xys++xys') []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x:xs):srcCells = srcPart
+                          yys   :trgCells = trgPart
+                          srcPart' = intersectCells (xs : srcCells) (dps M.! x)
+                      in T False (unLeft xys) -- the L.sort in the following line is so that we traverse vertices in natural order
+                         [dfs ((x,y):xys) ((unzip . L.sort) (zip (filter (not . null) srcPart') (filter (not . null) trgPart')))
+                         | (y,ys) <- picks yys,
+                           let trgPart' = intersectCells (ys : trgCells) (dps M.! y),
+                           map length srcPart' == map length trgPart']
+    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+    (lefts, rights) = partitionEithers vs
+    lrPart = [map Left lefts, map Right rights] -- Partition the vertices into left and right, to exclude crossover auts
+    unLeft xys = [(x,y) | (Left x, Left y) <- xys] -- also filters out Rights
+    es' = S.fromList es
+    dps = M.fromList [(v, distancePartitionS vs es' v) | v <- vs]
 
 -- |Given the incidence graph of an incidence structure between points and blocks
 -- (for example, a set system),
--- @incidenceAuts g@ returns generators for the automorphism group of the incidence structure.
+-- @incidenceAuts g@ returns a strong generating set 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.
--- 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 = 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)
+incidenceAuts = filter (/= p []) . strongTerminals . incidenceAutsDistancePartitionSearchTree
 
-incidenceAutsCon g@(G vs es) 
-    | isConnected g = map points (incidenceAuts' [] [vs])
-    | 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)
-              let p' = L.sort $ refine (ys:pt) (dps M.! x)
-              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' _ [] = []
-          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)
-                  uus = zip us us
-              in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
-                 []  -> level us p x ys hs
-                 h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
-          level _ _ _ _ _ = [] -- includes the case where y matches Right _, which can only occur on first level, before we've distance partitioned
-          dps = M.fromList [(v, distancePartition g v) | v <- vs]
-          es' = S.fromList es
-          nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
 
+-- TODO: Filter out rights, map unLeft - to bs and gs
+incidenceAuts2 g = (partitionBSGS g) (equitablePartitionIncidenceSearchTree g)
+    where unLeft (Left x) = x
+          -- map (\g -> fromPairs . map (\(Left x, Left y) -> (x,y)) . filter (\(x,y) -> isLeft x) . toPairs) gs
 
+
 -- GRAPH ISOMORPHISMS
 
 -- !! 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
     | length cs1 /= length cs2 = []
@@ -404,9 +634,8 @@
     where cs1 = map (inducedSubgraph g1) (components g1)
           cs2 = map (inducedSubgraph g2) (components g2)
           graphIsos' (ci:cis) cjs =
-              [iso ++ iso' | cj <- cjs,
+              [iso ++ iso' | (cj,cjs') <- picks cjs,
                              iso <- graphIsosCon ci cj,
-                             let cjs' = L.delete cj cjs,
                              iso' <- graphIsos' cis cjs']
           graphIsos' [] [] = [[]]
 
@@ -428,12 +657,16 @@
                      else let (x:xs):p1'' = p1'
                               ys:p2'' = p2'
                           in concat [dfs ((x,y):xys)
-                                         (refine' (xs : p1'') (dps1 M.! x))
-                                         (refine' ((L.delete y ys):p2'') (dps2 M.! y))
-                                         | y <- ys]
+                                         (intersectCells (xs : p1'') (dps1 M.! x))
+                                         (intersectCells (ys': p2'') (dps2 M.! y))
+                                         | (y,ys') <- picks ys]
           isCompatible xys = and [([x,x'] `S.member` es1) == (L.sort [y,y'] `S.member` es2) | (x,y) <- xys, (x',y') <- xys, x < x']
-          dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
-          dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
+          dps1 = M.fromAscList [(v, distancePartitionS vs1 es1 v) | v <- vs1]
+          dps2 = M.fromAscList [(v, distancePartitionS vs2 es2 v) | v <- vs2]
+          -- dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
+          -- dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
+          vs1 = vertices g1
+          vs2 = vertices g2
           es1 = S.fromList $ edges g1
           es2 = S.fromList $ edges g2
 
@@ -445,10 +678,7 @@
 -- !! then the cost of calculating distancePartitions may not be warranted
 -- !! (see Math.Combinatorics.Poset: orderIsos01 versus orderIsos)
 
--- !! deprecate
-isIso g1 g2 = (not . null) (graphIsos g1 g2)
 
-
 -- the following differs from graphIsos in only two ways
 -- we avoid Left, Right crossover isos, by insisting that a Left is taken to a Left (first two lines)
 -- we return only the action on the Lefts, and unLeft it
@@ -461,9 +691,8 @@
     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 ++ iso' | (cj,cjs') <- picks cjs,
                              iso <- incidenceIsosCon ci cj,
-                             let cjs' = L.delete cj cjs,
                              iso' <- incidenceIsos' cis cjs']
           incidenceIsos' [] [] = [[]]
 
@@ -484,9 +713,9 @@
                      else let (x:xs):p1'' = p1'
                               ys:p2'' = p2'
                           in concat [dfs ((x,y):xys)
-                                         (refine' (xs : p1'') (dps1 M.! x))
-                                         (refine' ((L.delete y ys):p2'') (dps2 M.! y))
-                                         | y <- ys]
+                                         (intersectCells (xs : p1'') (dps1 M.! x))
+                                         (intersectCells (ys': p2'') (dps2 M.! y))
+                                         | (y,ys') <- picks ys]
           isCompatible xys = and [([x,x'] `S.member` es1) == (L.sort [y,y'] `S.member` es2) | (x,y) <- xys, (x',y') <- xys, x < x']
           dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
           dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
@@ -498,63 +727,5 @@
      Graph (Either p1 b1) -> Graph (Either p2 b2) -> Bool
 isIncidenceIso g1 g2 = (not . null) (incidenceIsos g1 g2)
 
-{-
-removeGens x gs = removeGens' [] gs where
-    baseOrbit = x .^^ gs
-    removeGens' ls (r:rs) =
-        if x .^^ (ls++rs) == baseOrbit
-        then removeGens' ls rs
-        else removeGens' (r:ls) rs
-    removeGens' ls [] = reverse ls
--- !! reverse is probably pointless
-
-
--- !! DON'T THINK THIS IS WORKING PROPERLY
--- eg graphAutsSGSNew $ toGraph ([1..7],[[1,3],[2,3],[3,4],[4,5],[4,6],[4,7]])
--- returns [[[1,2]],[[5,6]],[[5,7,6]],[[6,7]]]
--- whereas [[6,7]] was a Schreier generator, so shouldn't have been listed
-
--- Using Schreier generators to seed the next level
--- At the moment this is slower than the above
--- (This could be modified to allow us to start the search with a known subgroup)
-graphAutsNew g@(G vs es) = graphAuts' [] [] [vs] where
-    graphAuts' us hs p@((x:ys):pt) =
-        let ys' = ys L.\\ (x .^^ hs) -- don't need to consider points which can already be reached from Schreier generators
-            hs' = level us p x ys' []
-            p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
-            reps = cosetRepsGx (hs'++hs) x
-            schreierGens = removeGens x $ schreierGeneratorsGx (x,reps) (hs'++hs)
-        in hs' ++ graphAuts' (x:us) schreierGens p'
-    graphAuts' us hs ([]:pt) = graphAuts' us hs pt
-    graphAuts' _ _ [] = []
-    level us p@(ph:pt) x (y:ys) hs =
-        let px = refine' (L.delete x ph : pt) (dps M.! x)
-            py = refine' (L.delete y ph : pt) (dps M.! y)
-            uus = zip us us
-        in if map length px /= map length py
-           then level us p x ys hs
-           else case dfs ((x,y):uus) (filter (not . null) px) (filter (not . null) py) of
-                []  -> level us p x ys hs
-                h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
-                -- if h1 = (1 2)(3 4), and h2 = (1 3 2), then we can remove 4 too
-    level _ _ _ [] _ = []
-    dfs xys p1 p2
-        | map length p1 /= map length p2 = []
-        | otherwise =
-            let p1' = filter (not . null) p1
-                p2' = filter (not . null) p2
-            in if all isSingleton p1'
-               then let xys' = xys ++ zip (concat p1') (concat p2')
-                    in if isCompatible xys' then [fromPairs' xys'] else []
-               else let (x:xs):p1'' = p1'
-                        ys:p2'' = p2'
-                    in concat [dfs ((x,y):xys)
-                                   (refine' (xs : p1'') (dps M.! x))
-                                   (refine' ((L.delete y ys):p2'') (dps M.! y))
-                                   | y <- ys]
-    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
-    dps = M.fromList [(v, distancePartition g v) | v <- vs]
-    es' = S.fromList es
--}
 
 
diff --git a/Math/Combinatorics/IncidenceAlgebra.hs b/Math/Combinatorics/IncidenceAlgebra.hs
--- a/Math/Combinatorics/IncidenceAlgebra.hs
+++ b/Math/Combinatorics/IncidenceAlgebra.hs
@@ -1,10 +1,12 @@
--- Copyright (c) 2011, David Amos. All rights reserved.
+-- Copyright (c) 2011-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction #-}
 
 
 module Math.Combinatorics.IncidenceAlgebra where
 
+import Prelude hiding ( (*>) )
+
 import Math.Core.Utils
 
 import Math.Combinatorics.Digraph
@@ -103,11 +105,11 @@
     -- |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.
-    unit 0 = zero -- so that sum works
+    unit 0 = zerov -- so that sum works
     -- unit x = x *> sumv [return (Iv (a,a)) | a <- poset] -- the delta function
     -- but we can't know from the types alone which poset we are working in
     mult = linear mult'
-        where mult' (Iv poset (a,b), Iv _ (c,d)) = if b == c then return (Iv poset (a,d)) else zero
+        where mult' (Iv poset (a,b), Iv _ (c,d)) = if b == c then return (Iv poset (a,d)) else zerov
 
 -- So multiplication in the incidence algebra is about composition of intervals
 
diff --git a/Math/Combinatorics/Matroid.hs b/Math/Combinatorics/Matroid.hs
--- a/Math/Combinatorics/Matroid.hs
+++ b/Math/Combinatorics/Matroid.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2011, David Amos. All rights reserved.
+-- Copyright (c) 2011-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE NoMonomorphismRestriction, DeriveFunctor #-}
 
@@ -6,6 +6,8 @@
 module Math.Combinatorics.Matroid where
 
 -- Source: Oxley, Matroid Theory (second edition)
+
+import Prelude hiding ( (*>) )
 
 import Math.Core.Utils
 import Math.Core.Field hiding (f7)
diff --git a/Math/Combinatorics/Poset.hs b/Math/Combinatorics/Poset.hs
--- a/Math/Combinatorics/Poset.hs
+++ b/Math/Combinatorics/Poset.hs
@@ -7,7 +7,7 @@
 
 import Math.Common.ListSet as LS -- set operations on strictly ascending lists
 import Math.Core.Utils -- for set/multiset operations on ordered lists
-import Math.Algebra.Field.Base
+import Math.Core.Field
 import Math.Combinatorics.FiniteGeometry
 import Math.Algebra.LinearAlgebra
 
@@ -72,9 +72,13 @@
 
 -- LATTICE OF (POSITIVE) DIVISORS OF N
 
-divides a b = b `mod` a == 0
+divides a b = b `rem` a == 0
 
-divisors n | n >= 1 = [a | a <- [1..n], a `divides` n]
+divisors n = toSet [ d' | d <- takeWhile (\d -> d*d <= n) [1..],
+                          let (q,r) = n `quotRem` d, r == 0,
+                          d' <- [d,q] ]
+-- The toSet call sorts, and deduplicates if n is a square
+-- divisors n | n >= 1 = [a | a <- [1..n], a `divides` n]
 
 -- |posetD n is the lattice of (positive) divisors of n
 posetD :: Int -> Poset Int
@@ -108,7 +112,7 @@
 partitions [] = [[]]
 partitions [x] = [[[x]]]
 partitions (x:xs) = let ps = partitions xs in
-    map ([x]:) ps ++ [ (x:cell):(L.delete cell p) | p <- ps, cell <- p]
+    map ([x]:) ps ++ [ (x:cell):p' | p <- ps, (cell,p') <- picks p]
 -- if the input is sorted, then so is the output
 
 isRefinement a b = and [or [acell `isSubset` bcell | bcell <- b] | acell <- a]
@@ -123,6 +127,7 @@
 
 
 -- LATTICE OF INTERVAL PARTITIONS OF [1..N] ORDERED BY REFINEMENT
+-- Interval partitions of [1..n] correspond to compositions of n
 
 intervalPartitions xs = filter (all isInterval) (partitions xs)
 
@@ -138,16 +143,12 @@
 
 -- LATTICE OF INTEGER PARTITIONS OF N ORDERED BY REFINEMENT
 
-integerPartitions1 n = ips (reverse [1..n]) n
-    where ips [] 0 = [[]]
-          ips [] _ = []
-          ips (x:xs) n | x > n     = ips xs n
-                       | otherwise = map (x:) (ips (x:xs) (n-x)) ++ ips xs n
-
 -- For example, integerPartitions 5 -> [ [5], [4,1], [3,2], [3,1,1], [2,2,1], [2,1,1,1], [1,1,1,1,1] ]
-integerPartitions n = dfs ([],n,n)
-    where dfs (xs, 0, _) = [reverse xs]
-          dfs (xs, r, i) = concatMap dfs [ (i':xs, r-i', i') | i' <- reverse [1..min r i] ]
+integerPartitions n | n >= 0 = ips n n where
+    ips 0 _ = [[]]
+    ips _ 0 = []
+    ips n m | m <= n = map (m:) (ips (n-m) m) ++ ips n (m-1)
+            | otherwise = ips n n
 
 isIPRefinement ys xs = dfs xs ys
     where dfs (x:xs) (y:ys) | x < y = False
@@ -187,7 +188,8 @@
 -- 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 :: (Eq fq, FiniteField fq) => Int -> [fq] -> Poset [[fq]]
+-- Example usage: posetL 2 f3
+posetL :: (Eq fq, Num fq) => Int -> [fq] -> Poset [[fq]]
 posetL n fq = Poset ( subspaces fq n, isSubspace ) 
 
 
@@ -259,22 +261,26 @@
     and [ x `poa` y == f x `pob` f y | x <- seta, y <- seta ]
 
 -- Find all order isomorphisms between two posets
--- This algorithm is faster to find out whether or not there are any
+-- This is the most naive algorithm, and should not be used on larger posets
+-- For example, already on the following, this takes forever compared to almost instant for orderIsos:
+-- > head $ orderIsos01 (posetD $ 8*9*25*49) (posetD $ 4*27*25*121)
 orderIsos01 (Poset (seta,poa)) (Poset (setb,pob))
     | length seta /= length setb = []
     | otherwise = orderIsos' [] seta setb
     where orderIsos' xys [] [] = [xys]
           orderIsos' xys (x:xs) ys =
-              concat [ orderIsos' ((x,y):xys) xs (L.delete y ys)
-                     | y <- ys, and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ]
+              concat [ orderIsos' ((x,y):xys) xs ys'
+                     | (y,ys') <- picks ys,
+                       and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ]
 
 -- |Are the two posets order-isomorphic?
-isOrderIso :: (Eq a, Eq b) => Poset a -> Poset b -> Bool
-isOrderIso poseta posetb = (not . null) (orderIsos01 poseta posetb)
+isOrderIso :: (Ord a, Ord b) => Poset a -> Poset b -> Bool
+isOrderIso poseta posetb = (not . null) (orderIsos poseta posetb)
 
--- Find all order isomorphisms between two posets
--- This algorithm is faster to find all isomorphisms, if there are many
--- (It may be that it is faster to find any, for large posets, but the break-even point seems to be quite big)
+-- For small posets, it may be that the up-front cost of calculating the hasseDigraph and heightPartitionDAG
+-- are not justified.
+-- |Find all order isomorphisms between two posets
+orderIsos :: (Ord a, Ord b) => Poset a -> Poset b -> [[(a,b)]]
 orderIsos posetA@(Poset (_,poa)) posetB@(Poset (_,pob))
     | map length heightPartA /= map length heightPartB = []
     | otherwise = dfs [] heightPartA heightPartB
@@ -283,8 +289,9 @@
           dfs xys [] [] = [xys]
           dfs xys ([]:las) ([]:lbs) = dfs xys las lbs
           dfs xys ((x:xs):las) (ys:lbs) =
-              concat [ dfs ((x,y):xys) (xs:las) (L.delete y ys : lbs)
-                     | y <- ys, and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ]
+              concat [ dfs ((x,y):xys) (xs:las) (ys' : lbs)
+                     | (y,ys') <- picks ys,
+                       and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ]
 -- A variant on this algorithm would use the Hasse digraph rather than the partial order in the test on the last line
 -- This might be faster, depending how expensive the partial order comparison function is
 -- In effect though, it would then be a DAG isomorphism function
diff --git a/Math/Combinatorics/StronglyRegularGraph.hs b/Math/Combinatorics/StronglyRegularGraph.hs
--- a/Math/Combinatorics/StronglyRegularGraph.hs
+++ b/Math/Combinatorics/StronglyRegularGraph.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2008, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, David Amos. All rights reserved.
 
 -- |A module defining various strongly regular graphs, including the Clebsch, Hoffman-Singleton, Higman-Sims, and McLaughlin graphs.
 --
@@ -8,6 +8,8 @@
 --
 -- Strongly regular graphs are highly symmetric, and have large automorphism groups.
 module Math.Combinatorics.StronglyRegularGraph where
+
+import Prelude hiding ( (*>) )
 
 import qualified Data.List as L
 import Data.Maybe (isJust)
diff --git a/Math/CommutativeAlgebra/Polynomial.hs b/Math/CommutativeAlgebra/Polynomial.hs
--- a/Math/CommutativeAlgebra/Polynomial.hs
+++ b/Math/CommutativeAlgebra/Polynomial.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2011, David Amos. All rights reserved.
+-- Copyright (c) 2011-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances, DeriveFunctor #-}
 
@@ -12,6 +12,8 @@
 --
 -- > [t,u,v,x,y,z] = map glexvar ["t","u","v","x","y","z"]
 module Math.CommutativeAlgebra.Polynomial where
+
+import Prelude hiding ( (*>) )
 
 import Math.Core.Field
 import Math.Core.Utils (toSet)
diff --git a/Math/Core/Field.hs b/Math/Core/Field.hs
--- a/Math/Core/Field.hs
+++ b/Math/Core/Field.hs
@@ -2,7 +2,8 @@
 
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 
--- |A module defining the field Q of rationals and the small finite fields F2, F3, F4, F5, F7, F8, F9, F11, F13, F16, F17, F19, F23, F25.
+-- |A module defining the field Q of rationals and the small finite fields (Galois fields)
+-- F2, F3, F4, F5, F7, F8, F9, F11, F13, F16, F17, F19, F23, F25.
 --
 -- Given a prime power q, Fq is the type representing elements of the field (eg @F4@),
 -- fq is a list of the elements of the field, beginning 0,1,... (eg @f4@),
@@ -380,10 +381,10 @@
 
 instance Num F9 where
     F9 x + F9 y = F9 $ z1 + z0
-        where z = x+y; z1 = (z .&. 0xff00) `mod` 0x300; z0 = (z .&. 0xff) `mod` 3
+        where z = x+y; z1 = (z .&. 0xff00) `rem` 0x300; z0 = (z .&. 0xff) `rem` 3
     negate (F9 x) = F9 $ z1 + z0
-        where z = 0x303 - x; z1 = (z .&. 0xff00) `mod` 0x300; z0 = (z .&. 0xff) `mod` 3
-    F9 x * F9 y = F9 $ ((z2 + z1) `mod` 0x300) + ((z2 + z0) `mod` 3) 
+        where z = 0x303 - x; z1 = (z .&. 0xff00) `rem` 0x300; z0 = (z .&. 0xff) `rem` 3
+    F9 x * F9 y = F9 $ ((z2 + z1) `rem` 0x300) + ((z2 + z0) `rem` 3) 
         where z = x*y; z2 = z .&. 0xff0000; z1 = z .&. 0xff00; z0 = z .&. 0xff
         -- Explanation: We are substituting x^2 = x+1.
         -- We could do z2 `shiftR` 8 and z2 `shiftR` 16
@@ -470,10 +471,10 @@
 
 instance Num F25 where
     F25 x + F25 y = F25 $ z1 + z0
-        where z = x+y; z1 = (z .&. 0xff00) `mod` 0x500; z0 = (z .&. 0xff) `mod` 5
+        where z = x+y; z1 = (z .&. 0xff00) `rem` 0x500; z0 = (z .&. 0xff) `rem` 5
     negate (F25 x) = F25 $ z1 + z0
-        where z = 0x505 - x; z1 = (z .&. 0xff00) `mod` 0x500; z0 = (z .&. 0xff) `mod` 5
-    F25 x * F25 y = F25 $ ((z2 + z1) `mod` 0x500) + ((3*z2 + z0) `mod` 5) 
+        where z = 0x505 - x; z1 = (z .&. 0xff00) `rem` 0x500; z0 = (z .&. 0xff) `rem` 5
+    F25 x * F25 y = F25 $ ((z2 + z1) `rem` 0x500) + ((3*z2 + z0) `rem` 5) 
         where z = x*y; z2 = z .&. 0xff0000; z1 = z .&. 0xff00; z0 = z .&. 0xff
         -- Explanation: We are substituting x^2 = x+3.
         -- We could do z2 `shiftR` 8 and z2 `shiftR` 16
diff --git a/Math/Core/Utils.hs b/Math/Core/Utils.hs
--- a/Math/Core/Utils.hs
+++ b/Math/Core/Utils.hs
@@ -26,6 +26,25 @@
     GT -> y : setUnionAsc (x:xs) ys
 setUnionAsc xs ys = xs ++ ys
 
+setUnionDesc :: Ord a => [a] -> [a] -> [a]
+setUnionDesc (x:xs) (y:ys) =
+    case compare x y of
+    GT -> x : setUnionDesc xs (y:ys)
+    EQ -> x : setUnionDesc xs ys
+    LT -> y : setUnionDesc (x:xs) ys
+setUnionDesc xs ys = xs ++ ys
+
+-- |The (multi-)set intersection of two ascending lists. If both inputs are strictly increasing,
+-- then the output is the set intersection and is strictly increasing. If both inputs are weakly increasing,
+-- then the output is the multiset intersection (with multiplicity), and is weakly increasing.
+intersectAsc :: Ord a => [a] -> [a] -> [a]
+intersectAsc (x:xs) (y:ys) =
+    case compare x y of
+    LT -> intersectAsc xs (y:ys)
+    EQ -> x : intersectAsc xs ys
+    GT -> intersectAsc (x:xs) ys
+intersectAsc _ _ = []
+
 -- |The multiset sum of two ascending lists. If xs and ys are ascending, then multisetSumAsc xs ys == sort (xs++ys).
 -- The code does not check that the lists are ascending.
 multisetSumAsc :: Ord a => [a] -> [a] -> [a]
@@ -36,7 +55,7 @@
     GT -> y : multisetSumAsc (x:xs) ys
 multisetSumAsc xs ys = xs ++ ys
 
--- |The multiset sum of two descending lists. If xs and ys are descending, then multisetSumDesc xs ys == sort (xs++ys).
+-- |The multiset sum of two descending lists. If xs and ys are descending, then multisetSumDesc xs ys == sortDesc (xs++ys).
 -- The code does not check that the lists are descending.
 multisetSumDesc :: Ord a => [a] -> [a] -> [a]
 multisetSumDesc (x:xs) (y:ys) =
@@ -79,6 +98,33 @@
     GT -> isSubMultisetAsc (x:xs) ys
 isSubMultisetAsc [] ys = True
 isSubMultisetAsc xs [] = False
+
+-- |Is the element in the ascending list?
+--
+-- With infinite lists, this can fail to terminate.
+-- For example, elemAsc 1 [1/2,3/4,7/8..] would fail to terminate.
+-- However, with a list of Integer, this will always terminate.
+elemAsc :: Ord a => a -> [a] -> Bool
+elemAsc x (y:ys) = case compare x y of
+                   LT -> False
+                   EQ -> True
+                   GT -> elemAsc x ys
+-- or x `elemAsc` ys = x `elem` takeWhile (<= x) ys
+
+-- |Is the element not in the ascending list? (With infinite lists, this can fail to terminate.)
+notElemAsc :: Ord a => a -> [a] -> Bool
+notElemAsc x (y:ys) = case compare x y of
+                      LT -> True
+                      EQ -> False
+                      GT -> notElemAsc x ys
+
+
+-- From Conor McBride
+-- http://stackoverflow.com/questions/12869097/splitting-list-into-a-list-of-possible-tuples/12872133#12872133
+-- |Return all the ways to \"pick one and leave the others\" from a list
+picks :: [a] -> [(a,[a])]
+picks [] = []
+picks (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- picks xs]
 
 
 pairs (x:xs) = map (x,) xs ++ pairs xs
diff --git a/Math/NumberTheory/Factor.hs b/Math/NumberTheory/Factor.hs
--- a/Math/NumberTheory/Factor.hs
+++ b/Math/NumberTheory/Factor.hs
@@ -1,16 +1,72 @@
 -- Copyright (c) 2006-2011, David Amos. All rights reserved.
 
+{-# LANGUAGE BangPatterns #-}
+
 -- |A module for finding prime factors.
 module Math.NumberTheory.Factor (module Math.NumberTheory.Prime,
-                                 pfactors) where
+                                 pfactors, ppfactors, pfactorsTo, ppfactorsTo) where
 
-import Math.NumberTheory.Prime
+import Control.Arrow (second, (&&&))
 import Data.Either (lefts)
-import Data.List (zip4)
+import Data.List as L
+import Math.Core.Utils (multisetSumAsc)
+import Math.NumberTheory.Prime
 
--- Cohen, A Course in Computational Algebraic Number Theory, p488
+-- |List the prime factors of n (with multiplicity). For example:
+-- >>> pfactors 60
+-- [2,2,3,5]
+--
+-- This says that 60 = 2 * 2 * 3 * 5
+-- 
+-- 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]
+pfactors n | n > 0 = pfactors' n $ takeWhile (< 10000) primes
+           | n < 0 = -1 : pfactors' (-n) (takeWhile (< 10000) primes)
+    where pfactors' n (d:ds) | n == 1 = []
+                             | n < d*d = [n]
+                             | r == 0 = d : pfactors' q (d:ds)
+                             | otherwise = pfactors' n ds
+                             where (q,r) = quotRem n d
+          pfactors' n [] = pfactors'' n
+          pfactors'' n = if isMillerRabinPrime n then [n]
+                         else let d = findFactorParallelECM n -- findFactorECM n
+                              in multisetSumAsc (pfactors'' d) (pfactors'' (n `div` d))
 
+-- |List the prime power factors of n. For example:
+-- >>> ppfactors 60
+-- [(2,2),(3,1),(5,1)]
+--
+-- This says that 60 = 2^2 * 3^1 * 5^1
+ppfactors :: Integer -> [(Integer,Int)]
+ppfactors = map (head &&& length) . L.group . pfactors
+-- ppfactors = map (\xs -> (head xs, length xs)) . L.group . pfactors
 
+-- |Find the prime factors of all numbers up to n. Thus @pfactorsTo n@ is equivalent to @[(m, pfactors m) | m <- [1..n]]@,
+-- except that the results are not returned in order. For example:
+-- >>> pfactorsTo 10
+-- [(8,[2,2,2]),(4,[2,2]),(6,[3,2]),(10,[5,2]),(2,[2]),(9,[3,3]),(3,[3]),(5,[5]),(7,[7]),(1,[])]
+--
+-- @pfactorsTo n@ is significantly faster than @map pfactors [1..n]@ for larger n.
+pfactorsTo n = pfactorsTo' (1,[]) primes where
+    pfactorsTo' (!m,!qs) ps@(ph:pt) | m' > n = [(m,qs)]
+                                    | otherwise = pfactorsTo' (m',ph:qs) ps ++ pfactorsTo' (m,qs) pt
+        where m' = m*ph
+-- We avoid a reverse call, because it does make a noticeable difference to the speed.
+
+-- |Find the prime power factors of all numbers up to n. Thus @ppfactorsTo n@ is equivalent to @[(m, ppfactors m) | m <- [1..n]]@,
+-- except that the results are not returned in order. For example:
+-- >>> ppfactorsTo 10
+-- [(8,[(2,3)]),(4,[(2,2)]),(6,[(3,1),(2,1)]),(10,[(5,1),(2,1)]),(2,[(2,1)]),(9,[(3,2)]),(3,[(3,1)]),(5,[(5,1)]),(7,[(7,1)]),(1,[])]
+--
+-- @ppfactorsTo n@ is significantly faster than @map ppfactors [1..n]@ for larger n.
+ppfactorsTo = map (second (map (head &&& length) . L.group)) . pfactorsTo
+
+
+-- Cohen, A Course in Computational Algebraic Number Theory, p488
+
 -- return (u,v,d) s.t ua+vb = d, with d = gcd a b
 extendedEuclid a b
     | b == 0 = (1,0,a)
@@ -18,7 +74,6 @@
                       (s,t,d) = extendedEuclid b r -- s*b+t*r == d
                   in (t,s-q*t,d)                   -- s*b+t*(a-q*b) == d
 
-
 -- ELLIPTIC CURVE ARITHMETIC
 
 data EllipticCurve = EC Integer Integer Integer deriving (Eq, Show)
@@ -66,7 +121,6 @@
                               (Left _, _) -> p'
                               (_, Left _) -> q'
 
-
 -- ELLIPTIC CURVE FACTORISATION
 
 -- We choose an elliptic curve E over Zn, and a point P on the curve
@@ -74,7 +128,6 @@
 -- What we are hoping is that at some stage we will fail because we can't invert an element in Zn
 -- This will lead to finding a non-trivial factor of n
 
-
 discriminantEC a b = 4 * a * a * a + 27 * b * b
 
 -- perform a sequence of scalar multiplications in the elliptic curve, hoping for a bailout
@@ -105,32 +158,7 @@
     -- the filter is because d might be a multiple of n,
     -- for example if the problem was that the discriminant was divisible by n
 
-
--- |List the prime factors of n (with multiplicity).
--- 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]
-pfactors n | n > 0 = pfactors' n $ takeWhile (< 10000) primes
-           | n < 0 = -1 : pfactors' (-n) (takeWhile (< 10000) primes)
-    where pfactors' n (d:ds) | n == 1 = []
-                             | n < d*d = [n]
-                             | r == 0 = d : pfactors' q (d:ds)
-                             | otherwise = pfactors' n ds
-                             where (q,r) = quotRem n d
-          pfactors' n [] = pfactors'' n
-          pfactors'' n = if isMillerRabinPrime n then [n]
-                         else let d = findFactorParallelECM n -- findFactorECM n
-                              in merge (pfactors'' d) (pfactors'' (n `div` d))
-
-merge (x:xs) (y:ys) =
-    case compare x y of
-    LT -> x : merge xs (y:ys)
-    EQ -> x : y : merge xs ys
-    GT -> y : merge (x:xs) ys
-merge xs ys = xs ++ ys
-
+-- TESTING MULTIPLE CURVES IN PARALLEL
 
 -- Cohen p489
 -- find inverse of as mod n in parallel, or a non-trivial factor of n
@@ -144,7 +172,7 @@
 
 parallelEcAdd n ecs ps1 ps2 =
     case parallelInverse n (zipWith f ps1 ps2) of
-    Right invs -> Right [g ec p1 p2 inv | (ec,p1,p2,inv) <- zip4 ecs ps1 ps2 invs]
+    Right invs -> Right [g ec p1 p2 inv | (ec,p1,p2,inv) <- L.zip4 ecs ps1 ps2 invs]
     Left d -> Left d
     where f Inf pt = 1
           f pt Inf = 1
diff --git a/Math/NumberTheory/Prime.hs b/Math/NumberTheory/Prime.hs
--- a/Math/NumberTheory/Prime.hs
+++ b/Math/NumberTheory/Prime.hs
@@ -11,16 +11,24 @@
 
 
 isTrialDivisionPrime n
-    | n > 1 = isNotDivisibleBy primes
+    | n > 1 = not $ any (\p -> n `rem` p == 0) (takeWhile (\p -> p*p <= n) primes)
     | otherwise = False
-    where isNotDivisibleBy (d:ds) | d*d > n         = True
-                                  | n `rem` d == 0  = False
-                                  | otherwise       = isNotDivisibleBy ds
 
 -- |A (lazy) list of the primes
 primes :: [Integer]
-primes = 2 : 3 : filter isTrialDivisionPrime (concat [ [m6-1,m6+1] | m6 <- [6,12..] ])
+primes = 2 : filter isPrime [3,5..] where
+    isPrime n = not $ any (\p -> n `rem` p == 0) (takeWhile (\p -> p*p <= n) primes)
 
+{-
+-- This is just marginally faster, but less elegant
+primes2 :: [Integer]
+primes2 = 2 : 3 : 5 : 7 : filter isPrime
+    (concat [ [m30+11,m30+13,m30+17,m30+19,m30+23,m30+29,m30+31,m30+37] | m30 <- [0,30..] ])
+    where isPrime n = not $ any (\p -> n `rem` p == 0) (takeWhile (\p -> p*p <= n) primes2')
+          primes2' = drop 3 primes2
+-}
+
+{-
 -- initial version. This isn't going to be very good if n has any "large" prime factors (eg > 10000)
 pfactors1 n | n > 0 = pfactors' n primes
             | n < 0 = -1 : pfactors' (-n) primes
@@ -29,7 +37,7 @@
                              | r == 0 = d : pfactors' q (d:ds)
                              | otherwise = pfactors' n ds
                              where (q,r) = quotRem n d
-
+-}
 
 -- MILLER-RABIN TEST
 -- Cohen, A Course in Computational Algebraic Number Theory, p422
@@ -59,7 +67,8 @@
 -- 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 `rem` n) (if even t then y else x*y `rem` n) (t `div` 2)
+          powerMod' x y t = let (q,r) = t `quotRem` 2
+                            in powerMod' (x*x `rem` n) (if r == 0 then y else x*y `rem` n) q
 
 isMillerRabinPrime' n
     | n >= 4 =
@@ -106,6 +115,7 @@
             where n6 = (n `div` 6) * 6
                   candidates = dropWhile (<= n) $ concat [ [m6+1,m6+5] | m6 <- [n6, n6+6..] ]
 
+{-
 -- slightly better version. This is okay so long as n has at most one "large" prime factor (> 10000)
 -- if it has more, it does at least tell you, via an error message, that it has run into difficulties
 pfactors2 n | n > 0 = pfactors' n $ takeWhile (< 10000) primes
@@ -116,4 +126,4 @@
                              | otherwise = pfactors' n ds
                              where (q,r) = quotRem n d
           pfactors' n [] = if isMillerRabinPrime n then [n] else error ("pfactors2: can't factor " ++ show n)
-
+-}
diff --git a/Math/NumberTheory/QuadraticField.hs b/Math/NumberTheory/QuadraticField.hs
--- a/Math/NumberTheory/QuadraticField.hs
+++ b/Math/NumberTheory/QuadraticField.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2011, David Amos. All rights reserved.
+-- Copyright (c) 2011-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances, OverlappingInstances #-}
 
@@ -16,7 +16,7 @@
 -- > i * sqrt(-3)
 module Math.NumberTheory.QuadraticField where
 
-import Prelude hiding (sqrt)
+import Prelude hiding (sqrt, (*>) )
 
 import Data.List as L
 import Math.Core.Field
diff --git a/Math/Projects/ChevalleyGroup/Classical.hs b/Math/Projects/ChevalleyGroup/Classical.hs
--- a/Math/Projects/ChevalleyGroup/Classical.hs
+++ b/Math/Projects/ChevalleyGroup/Classical.hs
@@ -1,6 +1,8 @@
--- Copyright (c) 2008, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Classical where
+
+import Prelude hiding ( (*>) )
 
 import Math.Algebra.Field.Base
 import Math.Algebra.Field.Extension hiding ( (<+>), (<*>) )
diff --git a/Math/Projects/ChevalleyGroup/Exceptional.hs b/Math/Projects/ChevalleyGroup/Exceptional.hs
--- a/Math/Projects/ChevalleyGroup/Exceptional.hs
+++ b/Math/Projects/ChevalleyGroup/Exceptional.hs
@@ -1,7 +1,9 @@
--- Copyright (c) 2008, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Exceptional where
 
+import Prelude hiding ( (*>) )
+
 import Data.List as L
 
 -- import Math.Algebra.Field.Base
@@ -165,7 +167,7 @@
 
 -- Unit imaginary octonions form one orbit under the action of G2
 
--- [alpha', beta', gamma' generate G2(3) as a permutation group on 702 points (the number of unit imaginary octonions over F3)
+-- [alpha', beta', gamma'] generate G2(3) as a permutation group on 702 points (the number of unit imaginary octonions over F3)
 -- Interestingly, http://brauer.maths.qmul.ac.uk/Atlas/v3/exc/G23/ doesn't seem to have this permutation representation
 
 
diff --git a/Math/Projects/MiniquaternionGeometry.hs b/Math/Projects/MiniquaternionGeometry.hs
--- a/Math/Projects/MiniquaternionGeometry.hs
+++ b/Math/Projects/MiniquaternionGeometry.hs
@@ -1,6 +1,8 @@
--- Copyright (c) David Amos, 2009. All rights reserved.
+-- Copyright (c) David Amos, 2009-2015. All rights reserved.
 
 module Math.Projects.MiniquaternionGeometry where
+
+import Prelude hiding ( (*>) )
 
 import qualified Data.List as L
 
diff --git a/Math/Projects/RootSystem.hs b/Math/Projects/RootSystem.hs
--- a/Math/Projects/RootSystem.hs
+++ b/Math/Projects/RootSystem.hs
@@ -1,6 +1,8 @@
--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) David Amos, 2008-2015. All rights reserved.
 
 module Math.Projects.RootSystem where
+
+import Prelude hiding ( (*>) )
 
 import Data.Ratio
 import qualified Data.List as L
diff --git a/Math/QuantumAlgebra/OrientedTangle.hs b/Math/QuantumAlgebra/OrientedTangle.hs
--- a/Math/QuantumAlgebra/OrientedTangle.hs
+++ b/Math/QuantumAlgebra/OrientedTangle.hs
@@ -1,10 +1,12 @@
--- Copyright (c) David Amos, 2010. All rights reserved.
+-- Copyright (c) David Amos, 2010-2015. All rights reserved.
 
 {-# LANGUAGE TypeFamilies, EmptyDataDecls #-}
 
 
 module Math.QuantumAlgebra.OrientedTangle where
 
+import Prelude hiding ( (*>), (<*>) )
+
 import Math.Algebra.Field.Base
 import Math.Algebras.LaurentPoly -- hiding (lvar, q, q')
 
@@ -26,7 +28,7 @@
 data OrientedTangle
 
 -- In GHCi 6.12.1, we appear to be limited to 8 value constructors within an associated data family
-instance Category OrientedTangle where
+instance MCategory OrientedTangle where
     data Ob OrientedTangle = OT [Oriented] deriving (Eq,Ord,Show)
     data Ar OrientedTangle = IdT [Oriented]
                            | CapT HorizDir
@@ -54,7 +56,7 @@
     target (SeqT as) = target (last as)
     a >>> b | target a == source b = SeqT [a,b]
 
-instance TensorCategory OrientedTangle where
+instance Monoidal OrientedTangle where
     tunit = OT []
     tob (OT as) (OT bs) = OT (as++bs)
     tar a b = ParT [a,b]
@@ -64,11 +66,11 @@
 idV = id
 idV' = id
 
-evalV  = \(E i, E j) -> if i + j == 0 then return () else zero
-evalV' = \(E i, E j) -> if i + j == 0 then return () else zero
+evalV  = \(E i, E j) -> if i + j == 0 then return () else zerov
+evalV' = \(E i, E j) -> if i + j == 0 then return () else zerov
 
-coevalV  m = foldl (<+>) zero [e i `te` e (-i) | i <- [1..m] ]
-coevalV' m = foldl (<+>) zero [e (-i) `te` e i | i <- [1..m] ]
+coevalV  m = foldl (<+>) zerov [e i `te` e (-i) | i <- [1..m] ]
+coevalV' m = foldl (<+>) zerov [e (-i) `te` e i | i <- [1..m] ]
 
 lambda m = q' ^ m -- q^-m
 
diff --git a/Math/QuantumAlgebra/QuantumPlane.hs b/Math/QuantumAlgebra/QuantumPlane.hs
--- a/Math/QuantumAlgebra/QuantumPlane.hs
+++ b/Math/QuantumAlgebra/QuantumPlane.hs
@@ -53,7 +53,7 @@
     powers (Aq20 m) = powers m
 
 instance Algebra (LaurentPoly Q) (Aq20 String) where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- V []
     unit x = V [(munit,x)] where munit = Aq20 (NCM 0 [])
     mult x = x''' where
         x' = mult $ fmap ( \(Aq20 a, Aq20 b) -> (a,b) ) x -- unwrap and multiply
@@ -75,7 +75,7 @@
     powers (Aq02 m) = powers m
 
 instance Algebra (LaurentPoly Q) (Aq02 String) where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- V []
     unit x = V [(munit,x)] where munit = Aq02 (NCM 0 [])
     mult x = x''' where
         x' = mult $ fmap ( \(Aq02 a, Aq02 b) -> (a,b) ) x -- unwrap and multiply
@@ -99,7 +99,7 @@
     powers (M2q m) = powers m
 
 instance Algebra (LaurentPoly Q) (M2q String) where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- V []
     unit x = V [(munit,x)] where munit = M2q (NCM 0 [])
     mult x = x''' where
         x' = mult $ fmap ( \(M2q a, M2q b) -> (a,b) ) x -- unwrap and multiply
@@ -188,7 +188,7 @@
     powers (SL2q m) = powers m
 
 instance Algebra (LaurentPoly Q) (SL2q String) where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- V []
     unit x = V [(munit,x)] where munit = SL2q (NCM 0 [])
     mult x = x''' where
         x' = mult $ fmap ( \(SL2q a, SL2q b) -> (a,b) ) x -- unwrap and multiply
diff --git a/Math/QuantumAlgebra/Tangle.hs b/Math/QuantumAlgebra/Tangle.hs
--- a/Math/QuantumAlgebra/Tangle.hs
+++ b/Math/QuantumAlgebra/Tangle.hs
@@ -1,4 +1,4 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, TypeFamilies, FlexibleInstances, EmptyDataDecls #-}
 
@@ -15,7 +15,7 @@
 import Math.Algebra.Field.Base
 import Math.Algebras.LaurentPoly
 
-import Math.QuantumAlgebra.TensorCategory
+import Math.QuantumAlgebra.TensorCategory hiding (Vect)
 
 
 instance Mon [a] where
@@ -25,7 +25,7 @@
 -- type TensorAlgebra k a = Vect k [a]
 
 instance (Eq k, Num k, Ord a) => Algebra k [a] where
-    unit 0 = zero -- V []
+    unit 0 = zerov -- V []
     unit x = V [(munit,x)]
     mult = nf . fmap (\(a,b) -> a `mmult` b)
 
@@ -37,7 +37,7 @@
 
 data Tangle
 
-instance Category Tangle where
+instance MCategory Tangle where
     data Ob Tangle = OT Int deriving (Eq,Ord,Show)
     data Ar Tangle = IdT Int
                    | CapT
@@ -71,7 +71,7 @@
 --    a >>> b | target a == source b = SeqT a b
     a >>> b | target a == source b = SeqT [a,b]
 
-instance TensorCategory Tangle where
+instance Monoidal Tangle where
     tunit = OT 0
     tob (OT a) (OT b) = OT (a+b)
 --    tar a b = ParT a b
@@ -93,7 +93,7 @@
 cup :: [Oriented] -> TangleRep [Oriented]
 cup [Plus, Minus] = (-q'^2) *> return []
 cup [Minus, Plus] = return []
-cup _ = zero
+cup _ = zerov
 
 -- also called xminus
 over :: [Oriented] -> TangleRep [Oriented]
diff --git a/Math/QuantumAlgebra/TensorCategory.hs b/Math/QuantumAlgebra/TensorCategory.hs
--- a/Math/QuantumAlgebra/TensorCategory.hs
+++ b/Math/QuantumAlgebra/TensorCategory.hs
@@ -1,14 +1,13 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2014, David Amos. All rights reserved.
 
-{-# LANGUAGE TypeFamilies, EmptyDataDecls #-}
+{-# LANGUAGE TypeFamilies, EmptyDataDecls, MultiParamTypeClasses #-}
 
--- |A module defining classes and example instances of categories and tensor categories
+-- |A module defining classes and example instances of categories, monoidal categories and braided categories
 module Math.QuantumAlgebra.TensorCategory where
 
-import Math.Algebra.Group.PermutationGroup
-
+import Data.List as L
 
-class Category c where
+class MCategory c where
     data Ob c :: *
     data Ar c :: *
     id_ :: Ob c -> Ar c
@@ -18,73 +17,169 @@
 -- whereas we want the objects to be values of a single type.
 
 
+class (MCategory a, MCategory b) => MFunctor a b where
+    fob :: Ob a -> Ob b -- functor on objects
+    far :: Ar a -> Ar b -- functor on arrows
+
+-- We could also define tensor functors and braided functors, which are just functors which commute appropriately
+-- with tensor and braiding operations
+
+
 -- Kassel p282
--- The following is actually definition of a strict tensor category
-class Category c => TensorCategory c where
+-- The following is actually definition of a _strict_ monoidal category
+-- Also called tensor category
+class MCategory c => Monoidal c where
     tunit :: Ob c
     tob :: Ob c -> Ob c -> Ob c -- tensor product of objects
     tar :: Ar c -> Ar c -> Ar c -- tensor product of arrows
 
-class TensorCategory c => StrictTensorCategory c where {}
+class Monoidal c => StrictMonoidal c where {}
 -- we want to be able to declare some tensor categories as strict
 
-class TensorCategory c => WeakTensorCategory c where
-    assoc :: Ob c -> Ob c -> Ob c -- (u `tob` v) `tob` w -> u `tob` (v `tob` w)
-    lunit :: Ob c -> Ob c         -- unit `tob` v -> v
-    runit :: Ob c -> Ob c         -- v `tob` unit -> v
+class Monoidal c => WeakMonoidal c where
+    assoc :: Ob c -> Ob c -> Ob c -> Ar c -- assoc u v w is an arrow (natural transformation?): (u `tob` v) `tob` w -> u `tob` (v `tob` w)
+    lunit :: Ob c -> Ar c                 -- lunit v is an arrow (isomorphism): tunit `tob` v -> v
+    runit :: Ob c -> Ar c                 -- runit v is an arrow (isomorphism): v `tob` tunit -> v
 
 {-
-instance (TensorCategory c, Eq (Ar c), Show (Ar c)) => Num (Ar c) where
+instance (Monoidal c, Eq (Ar c), Show (Ar c)) => Num (Ar c) where
     (*) = tar
 -}
 
+class Monoidal c => Braided c where
+    twist :: Ob c -> Ob c -> Ar c
+    -- twist v w is a map from v tensor w to w tensor v
+    -- twist must be natural, and satisfy certain commutative diagrams - Kock 161, 169
 
--- SYMMETRIC GROUPOID
+class Braided c => Symmetric c where {}
+-- if twist satisfies twist v w >>> twist w v == id_ (v tensor w), then the category is symmetric
 
-data SymmetricGroupoid
 
-instance Category SymmetricGroupoid where
-    data Ob SymmetricGroupoid = OS Int deriving (Eq,Ord,Show)
-    data Ar SymmetricGroupoid = AS Int (Permutation Int) deriving (Eq,Ord,Show)
-    id_ (OS n) = AS n 1
-    source (AS n _) = OS n
-    target (AS n _) = OS n
-    AS m g >>> AS n h | m == n = AS m (g*h)
+-- SIMPLEX CATEGORIES
 
-instance TensorCategory SymmetricGroupoid where
-    tunit = OS 0
-    tob (OS m) (OS n) = OS (m+n)
-    tar (AS m g) (AS n h) = AS (m+n) (g * h~^k)
-        where k = p [[1..m+n]] ^ m
---    tar (AS m g) (AS n h) = AS (m+n) (fromPairs $ toPairs g ++ map (\(x,y)->(x+m,y+m)) (toPairs h))
+-- Kock, Frobenius Algebras ..., p178-9
+-- The skeleton of FinOrd (finite ordered sets)
+-- The objects are the finite ordinals n == [0..n-1]
+-- The arrows are the order-preserving maps
+data FinOrd
 
+instance MCategory FinOrd where
+    data Ob FinOrd = FinOrdOb Int deriving (Eq,Ord,Show)
+    -- FinOrdOb n represents the oriented simplex n == [0..n-1]
+    data Ar FinOrd = FinOrdAr Int Int [Int] deriving (Eq,Ord,Show)
+    -- FinOrdAr s t fs represents the order-preserving map, zip [0..s-1] fs.
+    -- For example FinOrdAr 3 2 [0,0,1] represents the map 0 -> 0, 1 -> 0, 2 -> 1
+    id_ (FinOrdOb n) = FinOrdAr n n [0..n-1]
+    source (FinOrdAr s _ _) = FinOrdOb s
+    target (FinOrdAr _ t _) = FinOrdOb t
+    FinOrdAr sf tf fs >>> FinOrdAr sg tg gs | tf == sg = FinOrdAr sf tg [let j = fs !! i in gs !! j | i <- [0..sf-1] ]
 
+instance Monoidal FinOrd where
+    tunit = FinOrdOb 0
+    tob (FinOrdOb m) (FinOrdOb n) = FinOrdOb (m+n)
+    tar (FinOrdAr sf tf fs) (FinOrdAr sg tg gs) = FinOrdAr (sf+sg) (tf+tg) (fs ++ map (+tf) gs)
+
+finOrdAr s t fs | s == length fs && minimum fs >= 0 && maximum fs < t && isOrderPreserving fs
+    = FinOrdAr s t fs
+    where isOrderPreserving (f1:f2:fs) = f1 <= f2 && isOrderPreserving (f2:fs)
+          isOrderPreserving _ = True
+
+
+-- The skeleton of FinSet
+-- The objects are the finite cardinals n == {0..n-1} (with no order)
+-- The arrows are the maps
+data FinCard
+
+instance MCategory FinCard where
+    data Ob FinCard = FinCardOb Int deriving (Eq,Ord,Show)
+    -- FinCardOb n represents the unoriented simplex n == {0..n-1}
+    data Ar FinCard = FinCardAr Int Int [Int] deriving (Eq,Ord,Show)
+    -- FinCardAr s t fs represents the map, zip [0..s-1] fs.
+    -- For example FinCardAr 3 2 [0,1,0] represents the map 0 -> 0, 1 -> 1, 2 -> 0
+    id_ (FinCardOb n) = FinCardAr n n [0..n-1]
+    source (FinCardAr s _ _) = FinCardOb s
+    target (FinCardAr _ t _) = FinCardOb t
+    FinCardAr sf tf fs >>> FinCardAr sg tg gs | tf == sg = FinCardAr sf tg [let j = fs !! i in gs !! j | i <- [0..sf-1] ]
+
+instance Monoidal FinCard where
+    tunit = FinCardOb 0
+    tob (FinCardOb m) (FinCardOb n) = FinCardOb (m+n)
+    tar (FinCardAr sf tf fs) (FinCardAr sg tg gs) = FinCardAr (sf+sg) (tf+tg) (fs ++ map (+tf) gs)
+
+finCardAr s t fs | s == length fs && minimum fs >= 0 && maximum fs < t -- for finite cardinals, the map doesn't have to be order-preserving
+    = FinCardAr s t fs
+
+-- Finite permutations form a subcategory of FinCard
+-- having as objects the finite cardinals n == {0..n-1}
+-- and as arrows the bijections (== permutations)
+finPerm fs | L.sort fs == [0..n-1] = FinCardAr n n fs
+    where n = length fs
+-- (Note that these are permutations of [0..n-1], rather than [1..n])
+
+
+-- This is the forgetful functor FinOrd -> FinCard (FinSet)
+instance MFunctor FinOrd FinCard where
+    fob (FinOrdOb n) = FinCardOb n
+    far (FinOrdAr s t fs) = FinCardAr s t fs
+
+
 -- BRAID CATEGORY
 
 data Braid
 
-instance Category Braid where
-    data Ob Braid = OB Int deriving (Eq,Ord,Show)
-    data Ar Braid = AB Int [Int] deriving (Eq,Ord,Show)
-    id_ (OB n) = AB n []
-    source (AB n _) = OB n
-    target (AB n _) = OB n
-    AB m is >>> AB n js | m == n = AB m (is ++ js)
+instance MCategory Braid where
+    data Ob Braid = BraidOb Int deriving (Eq,Ord,Show)
+    data Ar Braid = BraidAr Int [Int] deriving (Eq,Ord,Show)
+    id_ (BraidOb n) = BraidAr n []
+    source (BraidAr n _) = BraidOb n
+    target (BraidAr n _) = BraidOb n
+    BraidAr m is >>> BraidAr n js | m == n = BraidAr m (cancel (reverse is) js)
+        where cancel (x:xs) (y:ys) = if x+y == 0 then cancel xs ys else reverse xs ++ x:y:ys
+              cancel xs ys = reverse xs ++ ys
 
-s n i | 0 < i && i < n = AB n [i]
+t n 0 = BraidAr n [] -- the identity braid
+t n i |  0 < i && i < n = BraidAr n [i]
+      | -n < i && i < 0 = BraidAr n [i]
+-- The generators of B_n are [t n i | i <- [1..n-1]]
 
-instance TensorCategory Braid where
-    tunit = OB 0
-    tob (OB a) (OB b) = OB (a+b)
-    tar (AB m is) (AB n js) = AB (m+n) (is ++ map (+m) js)
+-- The inverses of the braid generators
+t' n i | 0 < i && i < n = BraidAr n [-i]
 
+instance Monoidal Braid where
+    tunit = BraidOb 0
+    tob (BraidOb m) (BraidOb n) = BraidOb (m+n)
+    tar (BraidAr m is) (BraidAr n js) = BraidAr (m+n) (is ++ map (+m) js)
 
+instance Braided Braid where
+    twist (BraidOb m) (BraidOb n) = BraidAr (m+n) $ concat [[i..i+n-1] | i <- [m,m-1..1]]
 
+-- Note that in FinCard we consider the objects as [0..n-1], whereas in Braid we consider them as [1..n], so that s_i twists [i,i+1]
+instance MFunctor Braid FinCard where
+    fob (BraidOb n) = FinCardOb n
+    far (BraidAr n ss) = foldr (>>>) (id_ (FinCardOb n)) [finPerm ([0..ti-1] ++ [ti+1,ti] ++ [ti+2..n-1]) | si <- ss, let ti = abs si - 1]
 
+
+-- VECT
+
+data Vect k
+
+instance Num k => MCategory (Vect k) where
+    data Ob (Vect k) = VectOb Int deriving (Eq,Ord,Show)
+    data Ar (Vect k) = VectAr Int Int [[Int]] deriving (Eq,Ord,Show)
+    id_ (VectOb n) = VectAr n n idMx where idMx = [[if i == j then 1 else 0 | j <- [1..n]] | i <- [1..n]]
+    source (VectAr m _ _) = VectOb m
+    target (VectAr _ n _) = VectOb n
+    VectAr r c xss >>> VectAr r' c' yss | c == r' = undefined -- matrix multiplication
+
+-- functor from FinPerm to Vect k
+
+
+-- 2-COBORDISMS
+
 data Cob2
 -- works very similar to Tangle category
 
-instance Category Cob2 where
+instance MCategory Cob2 where
     data Ob Cob2 = O Int deriving (Eq,Ord,Show)
     data Ar Cob2 = Id Int
                  | Unit
@@ -111,7 +206,7 @@
     target (Seq a b) = target b
     a >>> b | target a == source b = Seq a b
 
-instance TensorCategory Cob2 where
+instance Monoidal Cob2 where
     tunit = O 0
     tob (O a) (O b) = O (a+b)
     tar a b = Par a b
diff --git a/Math/Test/TAlgebras/TOctonions.hs b/Math/Test/TAlgebras/TOctonions.hs
--- a/Math/Test/TAlgebras/TOctonions.hs
+++ b/Math/Test/TAlgebras/TOctonions.hs
@@ -1,6 +1,8 @@
--- Copyright (c) 2011, David Amos. All rights reserved.
+-- Copyright (c) 2011-2015, David Amos. All rights reserved.
 
 module Math.Test.TAlgebras.TOctonions where
+
+import Prelude hiding ( (*>) )
 
 import Test.QuickCheck
 
diff --git a/Math/Test/TAlgebras/TStructures.hs b/Math/Test/TAlgebras/TStructures.hs
--- a/Math/Test/TAlgebras/TStructures.hs
+++ b/Math/Test/TAlgebras/TStructures.hs
@@ -6,6 +6,8 @@
 
 module Math.Test.TAlgebras.TStructures where
 
+import Prelude hiding ( (*>) )
+
 -- import Test.QuickCheck
 -- don't actually need, as we don't define any Arbitrary instances here
 
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
@@ -1,9 +1,11 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE EmptyDataDecls, ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies, RankNTypes #-}
 
 module Math.Test.TAlgebras.TTensorProduct where
+
+import Prelude hiding ( (*>) )
 
 import Test.QuickCheck
 import Math.Algebras.VectorSpace
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
@@ -1,10 +1,12 @@
--- Copyright (c) 2010, David Amos. All rights reserved.
+-- Copyright (c) 2010-2015, David Amos. All rights reserved.
 
 {-# LANGUAGE FlexibleInstances, NoMonomorphismRestriction, ScopedTypeVariables, GeneralizedNewtypeDeriving #-}
 
 
 module Math.Test.TAlgebras.TVectorSpace where
 
+import Prelude hiding ( (*>) )
+
 import Test.QuickCheck
 import Math.Algebras.VectorSpace
 import Math.Algebras.TensorProduct
@@ -17,8 +19,8 @@
 prop_AddGrp (x,y,z) =
     x <+> (y <+> z) == (x <+> y) <+> z && -- associativity
     x <+> y == y <+> x                 && -- commutativity
-    x <+> zero == x                    && -- identity
-    x <+> neg x == zero                   -- inverse
+    x <+> zerov == x                   && -- identity
+    x <+> negatev x == zerov              -- inverse
 
 prop_VecSp (a,b,x,y,z) =
     prop_AddGrp (x,y,z) &&
@@ -51,9 +53,9 @@
 
 
 prop_Linear f (a,x,y) =
-    f (x <+> y) == f x <+> f y &&
-    f zero == zero             &&
-    f (neg x) == neg (f x)     &&
+    f (x <+> y) == f x <+> f y       &&
+    f zerov == zerov                 &&
+    f (negatev x) == negatev (f x)   &&
     f (a *> x) == a *> f x
 
 prop_LinearQn f (a,x,y) = prop_Linear f (a,x,y)
@@ -73,8 +75,12 @@
 -- DIRECT SUM
 
 {-
+instance Num k => Alternative (Vect k) where
+    (<|>) = mplus
+    empty = mzero
+
 instance Num k => MonadPlus (Vect k) where
-    mzero = zero
+    mzero = zerov
     mplus (V xs) (V ys) = V (xs++ys) -- need to call nf afterwards
 -}
 
@@ -96,7 +102,7 @@
 linfun avbs = linear f where
     f a = case lookup a avbs of
           Just vb -> vb
-          Nothing -> zero
+          Nothing -> zerov
 
 
 prop_Product (f',g',x) =
@@ -174,7 +180,7 @@
     ((a, b) -> Vect k c) -> Vect k (Either a b) -> Vect k c
 bilinear f = linear f . tensor
 
-dot = bilinear (\(a,b) -> if a == b then return () else zero)
+dot = bilinear (\(a,b) -> if a == b then return () else zerov)
 
 polymult = bilinear (\(E i, E j) -> return (E (i+j)))
 
diff --git a/Math/Test/TCombinatorics/TGraphAuts.hs b/Math/Test/TCombinatorics/TGraphAuts.hs
--- a/Math/Test/TCombinatorics/TGraphAuts.hs
+++ b/Math/Test/TCombinatorics/TGraphAuts.hs
@@ -6,10 +6,14 @@
 import Math.Core.Field hiding (f7)
 import Math.Core.Utils (combinationsOf)
 import Math.Algebra.Group.PermutationGroup as P
+import Math.Algebra.Group.RandomSchreierSims as SS
 import Math.Combinatorics.Graph as G
 import Math.Combinatorics.GraphAuts
 import Math.Combinatorics.Matroid as M
+import Math.Combinatorics.FiniteGeometry
 
+import qualified Math.Algebra.Field.Extension as F
+
 import Test.HUnit
 
 
@@ -20,7 +24,18 @@
     testlistIncidenceAutsOrder,
     testlistGraphIsos,
     testlistIsGraphIso,
-    testlistIncidenceIsos
+    testlistIncidenceIsos,
+    testlistVertexTransitive,
+    testlistNotVertexTransitive,
+    testlistEdgeTransitive,
+    testlistNotEdgeTransitive,
+    testlistArcTransitive,
+    testlistNotArcTransitive,
+    testlist2ArcTransitive,
+    testlistNot2ArcTransitive,
+    testlist4ArcTransitive,
+    testlistDistanceTransitive,
+    testlistNotDistanceTransitive
     ]
 
 
@@ -125,3 +140,121 @@
                           (G [Left 3, Left 4, Right 4] [[Left 4, Right 4]])
                           [[(1,4),(2,3)]]
     ]
+
+testcaseVertexTransitive (desc, graph) = TestCase $
+    assertBool ("isVertexTransitive " ++ desc) $ isVertexTransitive graph
+
+testlistVertexTransitive = TestList $ map testcaseVertexTransitive [
+    -- because we're mapping, these all have to be of same type, hence Graph Int
+    ("(q 3)", q 3),
+    ("(q 4)", q 4),
+    ("petersen", G.to1n petersen)
+    ]
+
+testcaseNotVertexTransitive (desc, graph) = TestCase $
+    assertBool ("not isVertexTransitive " ++ desc) $ (not . isVertexTransitive) graph
+
+testlistNotVertexTransitive = TestList $ map testcaseNotVertexTransitive [
+    ("(kb 2 3)", G.to1n $ kb 2 3),
+    ("(kb 3 4)", G.to1n $ kb 3 4),
+    ("regular not vertex transitive" , G [(1::Int)..8] [[1,2],[1,3],[1,8],[2,3],[2,4],[3,5],[4,5],[4,6],[5,7],[6,7],[6,8],[7,8]])
+    ]
+
+
+testcaseEdgeTransitive (desc, graph) = TestCase $
+    assertBool ("isEdgeTransitive " ++ desc) $ isEdgeTransitive graph
+
+testlistEdgeTransitive = TestList $ map testcaseEdgeTransitive [
+    ("(kb 2 3)", kb 2 3),
+    ("(kb 3 4)", kb 3 4)
+    ]
+
+testcaseNotEdgeTransitive (desc, graph) = TestCase $
+    assertBool ("not isEdgeTransitive " ++ desc) $ (not . isEdgeTransitive) graph
+
+testlistNotEdgeTransitive = TestList $ map testcaseNotEdgeTransitive [
+    ("pyramid 4", pyramid 4),
+    ("pyramid 5", pyramid 5),
+    ("prism 3", G.to1n $ prism 3),
+    ("prism 5", G.to1n $ prism 5)
+    ]
+    where pyramid n = let G vs es = c n in graph (0:vs, [[0,v] | v <- vs] ++ es)
+
+
+testcaseArcTransitive (desc, graph) = TestCase $
+    assertBool ("isArcTransitive " ++ desc) $ isArcTransitive graph
+
+testlistArcTransitive = TestList $ map testcaseArcTransitive [
+    -- Godsil and Royle, p60 - j v k i is arc-transitive
+    ("(j 4 2 0)", j 4 2 0),
+    ("(j 5 2 0)", j 5 2 0),
+    ("(j 5 2 1)", j 5 2 1)
+    ]
+
+testcaseNotArcTransitive (desc, graph) = TestCase $
+    assertBool ("not isArcTransitive " ++ desc) $ (not . isArcTransitive) graph
+
+testlistNotArcTransitive = TestList $ map testcaseNotArcTransitive [
+    ("kb 3 2", kb 3 2),
+    ("kb 4 3", kb 4 3)
+    ]
+
+
+testcase2ArcTransitive (desc, graph) = TestCase $
+    assertBool ("is2ArcTransitive " ++ desc) $ is2ArcTransitive graph
+
+testlist2ArcTransitive = TestList $ map testcase2ArcTransitive [
+    -- Godsil and Royle, p60 - j (2k+1) k 0 is 2-arc-transitive
+    ("(j 3 1 0)", j 3 1 0),
+    ("(j 5 2 0)", j 5 2 0),
+    ("(j 7 3 0)", j 7 3 0)
+    ]
+
+testcaseNot2ArcTransitive (desc, graph) = TestCase $
+    assertBool ("not is2ArcTransitive " ++ desc) $ (not . is2ArcTransitive) graph
+
+testlistNot2ArcTransitive = TestList $ map testcaseNot2ArcTransitive [
+    -- because a 2-arc can be two sides of a triangle, or not, so they are not all alike
+    ("octahedron", octahedron),
+    ("icosahedron", icosahedron)
+    ]
+
+
+testcase4ArcTransitive (desc, graph) = TestCase $
+    assertBool ("is4ArcTransitive " ++ desc) $ is4ArcTransitive graph
+
+testlist4ArcTransitive = TestList $ map testcase4ArcTransitive [
+    -- Godsil and Royle, p80-1
+    ("PG(2,F2)", G.to1n $ incidenceGraphPG 2 f2),
+    ("PG(2,F3)", G.to1n $ incidenceGraphPG 2 f3),
+    ("PG(2,F4)", G.to1n $ incidenceGraphPG 2 f4)
+    ]
+
+
+testcaseDistanceTransitive (desc, graph) = TestCase $
+    assertBool ("isDistanceTransitive " ++ desc) $ isDistanceTransitive graph
+
+testlistDistanceTransitive = TestList $ map testcaseDistanceTransitive [
+    ("(kb 3 3)", G.to1n $ kb 3 3),
+    ("(kb 4 4)", G.to1n $ kb 4 4),
+    ("(q 3)", G.to1n $ q 3),
+    ("(q 4)", G.to1n $ q 4),
+    ("petersen", G.to1n $ petersen),
+    -- Godsil and Royle, p67 - j v k (k-1) and j (2k+1) (k+1) 0 are distance-transitive
+    ("(j 3 2 1)", G.to1n $ j 3 2 1),
+    ("(j 4 2 1)", G.to1n $ j 4 2 1),
+    ("(j 5 3 2)", G.to1n $ j 5 3 2)
+    -- ("(j 3 2 0)", G.to1n $ j 3 2 0),
+    -- ("(j 5 3 0)", G.to1n $ j 5 3 0),
+    -- ("(j 7 4 0)", G.to1n $ j 7 4 0)
+    ]
+
+testcaseNotDistanceTransitive (desc, graph) = TestCase $
+    assertBool ("not isDistanceTransitive " ++ desc) $ (not . isDistanceTransitive) graph
+
+testlistNotDistanceTransitive = TestList $ map testcaseNotDistanceTransitive [
+    ("(prism 3)", prism 3),
+    -- not prism 4, which is the cube
+    ("(prism 5)", prism 5)
+    ]
+
diff --git a/Math/Test/TNumberTheory/TPrimeFactor.hs b/Math/Test/TNumberTheory/TPrimeFactor.hs
--- a/Math/Test/TNumberTheory/TPrimeFactor.hs
+++ b/Math/Test/TNumberTheory/TPrimeFactor.hs
@@ -11,6 +11,7 @@
 
 
 testlistPrimeFactor = TestList [
+    testcasePrimesList,
     testlistSmallPrimes,
     testlistMillerRabin,
     testlistMersennePrimes,
@@ -30,11 +31,14 @@
     testlistFactorOrder
     ]
 
+primesTo100 = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
 
+testcasePrimesList = TestCase $ assertEqual "primes list"
+    primesTo100 (takeWhile (<100) primes)
+
 testlistSmallPrimes = TestList [
     TestCase $ assertEqual "small primes"
-               [False,True,True,False,True,False,True,False,False,False]
-               (map isPrime [1..10]),
+               primesTo100 (filter isPrime [1..100]),
     TestCase $ assertBool "negative primes" (all notPrime [-10..0])
     ]