HaskellForMaths 0.1.8 → 0.1.9

raw patch · 10 files changed

+324/−317 lines, 10 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• HaskellForMaths.cabal +2/−2
• Math/Algebra/Group/PermutationGroup.hs +74/−87
• Math/Algebra/LinearAlgebra.hs +12/−12
• Math/Combinatorics/Design.hs +49/−113
• Math/Combinatorics/FiniteGeometry.hs +29/−10
• Math/Combinatorics/Graph.hs +2/−0
• Math/Combinatorics/GraphAuts.hs +144/−86
• Math/Test/TDesign.hs +5/−4
• Math/Test/TFiniteGeometry.hs +6/−3
• Math/Test/TGraph.hs +1/−0

@@ -1,7 +1,7 @@    Name:                HaskellForMaths
-   Version:             0.1.8
+   Version:             0.1.9
    Category:            Math
-   Description:         Math library - combinatorics, group theory, commutative algebra, non-commutative algebra
+   Description:         A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended for educational purposes, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra
    License:             BSD3
    License-file:        license.txt

@@ -1,4 +1,4 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) David Amos, 2008-2009. All rights reserved.
 
 module Math.Algebra.Group.PermutationGroup where
 
@@ -126,18 +126,7 @@ orbitB gs b = orbit (-^) b gs
 orbitE gs b = orbit (-^) b gs
 
-{-
--- orbit of a vertex / point
-x .^^ gs = closure [x] [ .^g | g <- gs]
-orbitV gs x = closure [x] [ .^g | g <- gs]
-orbitP gs x = closure [x] [ .^g | g <- gs]
 
--- orbit of an edge / block
-b -^^ gs = closure [b] [ -^g | g <- gs]
-orbitE gs b = closure [b] [ -^g | g <- gs]
-orbitB gs b = closure [b] [ -^g | g <- gs]
--}
-
 action xs f = fromPairs [(x, f x) | x <- xs]
 
 
@@ -285,6 +274,8 @@ -- 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
+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
@@ -323,6 +314,45 @@           conjClasses' (h:hs) = let c = conjClass gs h in c : conjClasses' (hs L.\\ c)
 -}
 
+
+-- given list of generators, try to find a shorter list
+reduceGens (1:gs) = reduceGens gs
+reduceGens (g:gs) = reduceGens' ([g], eltsS [g]) gs where
+    reduceGens' (gs,eltsgs) (h:hs) =
+        if h `S.member` eltsgs
+        then reduceGens' (gs,eltsgs) hs
+        else reduceGens' (h:gs, eltsS $ h:gs) hs
+    reduceGens' (gs,_) [] = reverse gs
+
+
+-- SUBGROUPS
+
+isSubgp hs gs = all (`S.member` gs') hs
+    where gs' = eltsS gs
+
+-- The following is similar to the "cyclic extension" method - Holt p385
+-- However, Holt only looks at normal cyclic extensions (ie, by an elt of prime order), and so only finds solvable subgps
+
+-- |Return the subgroups of a group. Only suitable for use on small groups (eg < 100 elts)
+subgps :: (Ord a, Show a) => [Permutation a] -> [[Permutation a]]
+subgps gs = [] : subgps' S.empty [] (map (:[]) hs) where
+    hs = filter isMinimal $ elts gs
+    subgps' found ls (r:rs) =
+        let ks = elts r in
+        if ks `S.member` found
+        then subgps' found ls rs
+        else r : subgps' (S.insert ks found) (r:ls) rs
+    subgps' found [] [] = []
+    subgps' found ls [] = subgps' found [] [l ++ [h] | l <- reverse ls, h <- hs, last l < h]
+
+-- g is the minimal elt in the cyclic subgp it generates
+isMinimal 1 = False
+isMinimal g = all (g <=) primitives -- g == minimum primitives
+    where powers = takeWhile (/=1) $ tail $ iterate (*g) 1
+          n = orderElt g -- == length powers + 1
+          primitives = filter (\h -> orderElt h == n) powers
+
+
 -- centralizer of a subgroup or a set of elts
 -- the centralizer of H in G is the set of elts of G which commute with all elts of H
 centralizer gs hs = [k | k <- elts gs, all (\h -> h*k == k*h) hs]
@@ -344,16 +374,6 @@ -- setwise stabiliser of a set
 setStab gs xs = [g | g <- elts gs, xs -^ g == xs]
 
-
--- given list of generators, try to find a shorter list
-reduceGens (1:gs) = reduceGens gs
-reduceGens (g:gs) = reduceGens' ([g], eltsS [g]) gs where
-    reduceGens' (gs,eltsgs) (h:hs) =
-        if h `S.member` eltsgs
-        then reduceGens' (gs,eltsgs) hs
-        else reduceGens' (h:gs, eltsS $ h:gs) hs
-    reduceGens' (gs,_) [] = reverse gs
-
 -- normal closure of H in G
 normalClosure gs hs = reduceGens $ hs ++ [h ~^ g | h <- hs, g <- gs ++ map inverse gs]
 
@@ -368,87 +388,54 @@ derivedSubgp gs = commutatorGp gs gs
 
 
--- ACTIONS ON COSETS AND SUBGROUPS (QUOTIENT GROUPS)
-
-isSubgp hs gs = all (isMember gs) hs
+-- ACTION ON COSETS, QUOTIENT GROUPS
 
-isNormal hs gs = isSubgp hs gs && all (isMember hs) [h~^g | h <- hs, g <- gs]
+xs -*- ys = toListSet [x*y | x <- xs, y <- ys]
 
+xs -*  y  = L.sort [x*y | x <- xs] -- == xs -*- [y]
+x   *- ys = L.sort [x*y | y <- ys] -- == [x] -*- ys
 
--- action of a group on cosets by right multiplication
--- (hs should be all elts, not just generators)
-hs **^ g = L.sort [h*g | h <- hs]
+-- |isNormal gs ks returns True if \<ks\> is normal in \<gs\>.
+-- Note, it is caller's responsibility to ensure that \<ks\> is a subgroup of \<gs\> (ie that each k is in \<gs\>).
+isNormal :: (Ord a, Show a) => [Permutation a] -> [Permutation a] -> Bool
+isNormal gs ks = all (== ks') [ (g^-1) *- ks' -* g | g <- gs]
+    where ks' = elts ks
 
--- Cosets are disjoint, which leads to Lagrange's theorem
+-- |Return the normal subgroups of a group. Only suitable for use on small groups (eg < 100 elts)
+normalSubgps :: (Ord a, Show a) => [Permutation a] -> [[Permutation a]]
+normalSubgps gs = filter (isNormal gs) (subgps gs)
 
--- cosets gs hs = closure [hs] [ **^ g | g <- gs]
-cosets gs hs = orbit (**^) hs gs
--- the group acts transitively on cosets of a subgp, so this gives all cosets
--- hs #^^ gs = orbit (#^) gs hs
+isSimple gs = length (normalSubgps gs) == 2
 
-cosetAction gs hs =
-    let _H = elts hs
-        cosets_H = cosets gs _H
-    in toSn [action cosets_H (**^ g) | g <- gs]
-    -- in toSn $ map (induced (**^) cosets_H) gs
+-- Note: caller must ensure that hs is a subgp of gs
+cosets gs hs = orbit (-*) hs' gs
+    where hs' = elts hs
 
--- if H normal in G, then each element within a given coset gives rise to the same action on other cosets,
--- and we get a well defined multiplication Hx * Hy = Hxy (where it doesn't depend on which coset rep we chose)
-quotientGp gs hs
-    | hs `isNormal` gs = gens $ cosetAction gs hs
-    | otherwise = error "quotientGp: not well defined unless H normal in G"
--- the call to gens removes identity and duplicates
+-- |quotientGp gs ks returns \<gs\> / \<ks\>
+quotientGp :: (Ord a, Show a) => [Permutation a] -> [Permutation a] -> [Permutation Int]
+quotientGp gs ks
+    | ks `isNormal` gs = gens $ toSn [action cosetsK (-* g) | g <- gs]
+    | otherwise = error "quotientGp: not well defined unless ks normal in gs"
+    where cosetsK = cosets gs ks
 
-gs // hs = quotientGp gs hs
+-- |Synonym for quotientGp
+(//) :: (Ord a, Show a) => [Permutation a] -> [Permutation a] -> [Permutation Int]
+gs // ks = quotientGp gs ks
 
 
--- action of group on a subgroup by conjugation
--- (hs should be all elts, not just generators)
-hs ~~^ g = L.sort [h ~^ g | h <- hs]
+-- action of group element on a subset by conjugation
+xs ~~^ g = L.sort [x ~^ g | x <- xs]
 
--- don't think that this is necessarily transitive on isomorphic subgps
--- conjugateSubgps gs hs = closure [hs] [ ~~^ g | g <- gs]
-conjugateSubgps gs hs = orbit (~~^) hs gs
--- hs ~~^^ gs = orbit (~~^) gs hs
+conjugateSubgps gs hs = orbit (~~^) hs' gs
+    where hs' = elts hs
+-- not necessarily transitive on isomorphic subgps - eg a gp with an outer aut
 
 subgpAction gs hs =
-    let _H = elts hs
-        conjugates_H = conjugateSubgps gs _H
-    in toSn [action conjugates_H (~~^ g) | g <- gs]
-    -- in toSn $ map (induced (~~^) conjugates_H) gs
+    let conjugatesH = conjugateSubgps gs hs
+    in toSn [action conjugatesH (~~^ g) | g <- gs]
 
 
 -- in cube gp, the subgps all appear to correspond to stabilisers of subsets, or of blocks
 
 
-
-{-
-OLDER VERSIONS
-
-
--- the orbit of a point or block under the action of a set of permutations
-orbit action x gs = S.toList $ orbitS action x gs
-
-orbitS action x gs = orbit' S.empty (S.singleton x) where
-    orbit' interior boundary
-        | S.null boundary = interior
-        | otherwise =
-            let interior' = S.union interior boundary
-                boundary' = S.fromList [p `action` g | g <- gs, p <- S.toList boundary] S.\\ interior'
-            in orbit' interior' boundary'
-
--- orbit of a point
--}
-{-
--- the induced action of g on a set of blocks
--- Note: the set of blocks must be closed under the action of g, otherwise we will get an error in fromPairs
--- To ensure that it is closed, generate the blocks as the orbit of a starting block
-inducedAction bs g = fromPairs [(b, b -^ g) | b <- bs]
-
-induced action bs g = fromPairs [(b, b `action` g) | b <- bs]
-
-inducedB bs g = induced (-^) bs g
--}
--- elts gs = orbit (*) 1 gs
--- eltsS gs = orbitS (*) 1 gs
 

@@ -9,7 +9,7 @@ --
 -- The mnemonic for many of the arithmetic operations is that the number of angle brackets
 -- 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.
+-- v \<*\>\> m is multiplication of a vector on the left by a matrix on the right.
 module Math.Algebra.LinearAlgebra where
 
 import qualified Data.List as L
@@ -26,50 +26,50 @@ 
 -- vector operations
 
--- |u <+> v returns the sum u+v of vectors
+-- |u \<+\> v returns the sum u+v of vectors
 (<+>) :: (Num a) => [a] -> [a] -> [a]
 u <+> v = zipWith (+) u v
 
--- |u <-> v returns the difference u-v of vectors
+-- |u \<-\> v returns the difference u-v of vectors
 (<->) :: (Num a) => [a] -> [a] -> [a]
 u <-> v = zipWith (-) u v
 
--- |k *> v returns the product k*v of the scalar k and the vector v
+-- |k *\> v returns the product k*v of the scalar k and the vector v
 (*>) :: (Num a) => a -> [a] -> [a]
 k *> v = map (k*) v
 
--- |u <.> v returns the dot product of vectors (also called inner or scalar product)
+-- |u \<.\> v returns the dot product of vectors (also called inner or scalar product)
 (<.>) :: (Num a) => [a] -> [a] -> a
 u <.> v = sum (zipWith (*) u v)
 
--- |u <*> v returns the tensor product of vectors (also called outer or matrix product)
+-- |u \<*\> v returns the tensor product of vectors (also called outer or matrix product)
 (<*>) :: (Num a) => [a] -> [a] -> [[a]]
 u <*> v = [ [a*b | b <- v] | a <- u]
 
 
 -- matrix operations
 
--- |a <<+>> b returns the sum a+b of matrices
+-- |a \<\<+\>\> b returns the sum a+b of matrices
 (<<+>>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
 a <<+>> b = (zipWith . zipWith) (+) a b
 
--- |a <<->> b returns the difference a-b of matrices
+-- |a \<\<-\>\> b returns the difference a-b of matrices
 (<<->>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
 a <<->> b = (zipWith . zipWith) (-) a b
 
--- |a <<*>> b returns the product a*b of matrices
+-- |a \<\<*\>\> b returns the product a*b of matrices
 (<<*>>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
 a <<*>> b = [ [u <.> v | v <- L.transpose b] | u <- a]
  
--- |k *> m returns the product k*m of the scalar k and the matrix m
+-- |k *\> m returns the product k*m of the scalar k and the matrix m
 (*>>) :: (Num a) => a -> [[a]] -> [[a]]
 k *>> m = (map . map) (k*) m
 
--- |m <<*> v is multiplication of a vector by a matrix on the left
+-- |m \<\<*\> v is multiplication of a vector by a matrix on the left
 (<<*>) :: (Num a) => [[a]] -> [a] -> [a]
 m <<*> v = map (<.> v) m
 
--- |v <*>> m is multiplication of a vector by a matrix on the right
+-- |v \<*\>\> m is multiplication of a vector by a matrix on the right
 (<*>>) :: (Num a) => [a] -> [[a]] -> [a]
 v <*>> m = map (v <.>) (L.transpose m)
 

@@ -7,13 +7,13 @@ import qualified Data.Map as M
 import qualified Data.Set as S
 
-import Math.Common.ListSet (symDiff)
+import Math.Common.ListSet (intersect, symDiff)
 import Math.Algebra.Field.Base
 import Math.Algebra.Field.Extension
 import Math.Algebra.Group.PermutationGroup hiding (elts, order, isMember)
 import Math.Algebra.Group.SchreierSims as SS
 import Math.Combinatorics.Graph as G hiding (to1n)
-import Math.Combinatorics.GraphAuts (refine, isSingleton, graphAuts, incidenceAuts, removeGens)
+import Math.Combinatorics.GraphAuts (refine', isSingleton, graphAuts, incidenceAuts) -- , removeGens)
 import Math.Combinatorics.FiniteGeometry
 
 -- Cameron & van Lint, Designs, Graphs, Codes and their Links
@@ -77,7 +77,7 @@         case reverse (takeWhile (isJust . snd) [(t, findlambda t d) | t <- [0..k] ]) of
         [] -> Nothing
         (t,Just lambda):_ -> Just (t,(v,k,lambda))
-
+-- Note that a 0-(v,k,lambda) design just means that there are lambda blocks, all of size k, with no other regularity
 
 isStructure t d = isJust $ tDesignParams t d
 
@@ -213,106 +213,6 @@ -- The incidence graph is a bipartite graph, so the distance function naturally partitions points from blocks
 
 
--- !! Not strictly correct to stop when a level is empty
--- For example, Fano plane, then if you fix 1,2, you can't move 3, but you can move 4
--- However, the reason it doesn't seem to cause any problem is that in fact it's almost certainly okay to just find the first level
--- (since we know the group is transitive on the points)
--- (In the graph case, non-transitive graphs like kb 2 3 require us to continue after null levels)
-
--- !! No, not okay to stop at first level
--- The strong generators of one level don't necessarily generate the whole group
--- For example, suppose the group is Sn. We might by chance find a transversal consisting entirely of elts of An
--- We might only finally get Sn when looking for transversal for n-1, and finding (n-1 n)
-
--- The distance partition of the incidence graph of a design will always look like this:
--- [x], [b | x `elem` b], xs \\ [x], [b | x `notElem` b],
--- When we refine it, we're crossing it with a similar partition but with a different x
--- This won't have much effect on the point cell, but will hopefully split the block cells in half
-
-
--- !! superceded by Math.Combinatorics.GraphAuts.incidenceAuts, leading to designAuts above
--- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
-designAuts2 d@(D xs bs) = map points (designAuts' [] [vs]) where
-    n = length xs
-    g@(G vs es) = incidenceGraph d
-    points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
-    designAuts' us p@((x@(Left _):ys):pt) =
-        let p' = L.sort $ filter (not . null) $ refine (ys:pt) (dps M.! x)
-        in level us p x ys []
-        ++ designAuts' (x:us) p'
-    designAuts' us ([]:pt) = designAuts' us pt
-    designAuts' _ (((Right _):_):_) = [] -- if we fix all the points, then the blocks must be fixed too 
-    -- designAuts' _ [] = []
-    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 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 _ _ _ [] _ = []
-    level _ _ _ _ _ = [] -- includes the case where y matches Right _, which can only occur on first level, before we've distance partitioned
-    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, G.distancePartition g v) | v <- vs]
-    es' = S.fromList es
-
-
-designAutsNew d@(D xs bs) = map points (designAuts' [] [] [vs]) where
-    n = length xs
-    g@(G vs es) = incidenceGraph d
-    -- xs' = take n vs -- the points of the design as vertices of the incidence graph
-    points (P m) = P $ M.fromList [(k,v) | (Left k, Left v) <- M.toList m]
-    designAuts' us hs p@((x@(Left _):ys):pt) =
-        let p' = L.sort $ filter (not . null) $ refine (ys:pt) (dps M.! x)
-            hs' = level us p x (ys L.\\ (x .^^ hs)) []
-            reps = cosetRepsGx (hs'++hs) x
-            schreierGens = removeGens x $ schreierGeneratorsGx (x,reps) (hs'++hs)
-        in hs' ++ designAuts' (x:us) schreierGens p'
-    designAuts' us hs ([]:pt) = designAuts' us hs pt
-    designAuts' us hs p@((x@(Right _):ys):pt) = []
-    designAuts' _ _ [] = []
-    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 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 _ _ _ [] _ = []
-    level _ _ _ _ _ = [] -- includes the case where y matches Right _, which can only occur on first level, before we've distance partitioned
-    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, G.distancePartition g v) | v <- vs]
-    es' = S.fromList es
-
-
 -- MATHIEU GROUPS AND WITT DESIGNS
 
 alphaL2_23 = p [[-1],[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]]                      -- t -> t+1
@@ -340,7 +240,7 @@ m22sgs = filter (\g -> 0.^g == 0) m23sgs
 -- order 443520
 
--- !! The above assume that the base is [-1,0,..], which isn't guaranteed
+-- sgs uses the base implied by the Ord instance, which will be [-1,0,..]
 
 
 -- Steiner system S(5,8,24)
@@ -352,27 +252,63 @@ s_5_8_24 = design ([-1..22], octad -^^ l2_23)
 -- S(5,8,24) constructed as the image of a single octad under the action of PSL(2,23)
 -- 759 blocks ( (24 `choose` 5) `div` (8 `choose` 5) )
+-- Automorphism group is M24
 
 s_4_7_23 = derivedDesign s_5_8_24 (-1)
 -- 253 blocks ( (23 `choose` 4) `div` (7 `choose` 4) )
+-- Automorphism group is M23
 
 s_3_6_22 = derivedDesign s_4_7_23 0
 -- 77 blocks
+-- Automorphism group is M22
 
+-- Derived design of s_3_6_22 is PG(2,F4)
+
+
+-- An alternative construction
+s_5_8_24' = D xs bs where
+    xs = [1..24]
+    bs = sift [] (combinationsOf 8 xs)
+    sift ls (r:rs) = if all ((<=4) . length) [r `intersect` l | l <- ls]
+                     then r : sift (r:ls) rs 
+                     else sift ls rs
+    sift ls [] = []
+
+
 -- Could test that m22sgs are all designAuts of s_3_6_22
 
-{-
-s_5_8_24' = octad -^^ m24 -- [alphaL2_23, betaL2_23, gammaL2_23]
--- S(5,8,24) constructed as the image of a single octad under the action of PSL(2,23)
--- 759 blocks ( (24 `choose` 5) `div` (8 `choose` 5) )
 
+-- S(5,6,12) and M12
 
-s_4_7_23' = [xs | (-1):xs <- s_5_8_24]
--- 253 blocks ( (23 `choose` 4) `div` (7 `choose` 4) )
+alphaL2_11 = p [[-1],[0,1,2,3,4,5,6,7,8,9,10]]          -- t -> t+1
+betaL2_11  = p [[-1],[0],[1,3,9,5,4],[2,6,7,10,8]]      -- t -> 3*t
+gammaL2_11 = p [[-1,0],[1,10],[2,5],[3,7],[4,8],[6,9]]  -- t -> -1/t
 
-s_3_6_22' = [xs | 0:xs <- s_4_7_23]
--- 77 blocks
--}
+l2_11 = [alphaL2_11, betaL2_11, gammaL2_11]
+
+deltaM12 = p [[-1],[0],[1],[2,10],[3,4],[5,9],[6,7],[8]]
+-- Conway&Sloane p271, 327
+
+hexad = [0,1,3,4,5,9]
+-- the squares (quadratic residues) in F11
+-- http://en.wikipedia.org/wiki/Steiner_system
+
+s_5_6_12 = design ([-1..10], hexad -^^ l2_11)
+-- S(5,6,12) constructed as the image of a single hexad under the action of PSL(2,11)
+-- 132 blocks ( (12 `choose` 5) `div` (6 `choose` 5) )
+-- Automorphism group is M12
+
+s_4_5_11 = derivedDesign s_5_6_12 (-1)
+-- 66 blocks
+-- Automorphism group is M11
+
+m12 = [alphaL2_11, betaL2_11, gammaL2_11, deltaM12]
+
+m12sgs = sgs m12
+-- order 95040
+
+m11sgs = filter (\g -> (-1).^g == -1) m12sgs
+-- order 7920
 
 
 {-

@@ -7,6 +7,8 @@ import Data.List as L
 import qualified Data.Set as S
 
+import Math.Common.ListSet (toListSet)
+
 import Math.Algebra.Field.Base
 import Math.Algebra.Field.Extension hiding ( (<+>) ) -- , (*>) )
 import Math.Algebra.LinearAlgebra -- hiding ( det )
@@ -14,6 +16,7 @@ import Math.Combinatorics.Graph
 import Math.Combinatorics.GraphAuts -- for use in GHCi
 import Math.Algebra.Group.PermutationGroup -- for use in GHCi
+import Math.Algebra.Group.SchreierSims as SS -- for use in GHCi
 
 -- |ptsAG n fq returns the points of the affine geometry AG(n,Fq), where fq are the elements of Fq
 ptsAG :: (FiniteField a) => Int -> [a] -> [[a]]
@@ -36,14 +39,25 @@ ispnf _ = False
 
 -- closure of points in AG(n,Fq)
--- result is sorted
+-- given p1, .., pk, we're looking for all a1 p1 + ... + ak pk, s.t. a1 + ... + ak = 1
+-- if m is the matrix with p1, .., pk as rows, and vs are the vectors [a1, .., ak]
+-- then this is the same as [v <*>> m | v <- vs] == [m' <<*> v | v <- vs]
 closureAG ps =
-    let multipliers = [ (1 - sum xs) : xs | xs <- ptsAG (k-1) fq ] -- k-vectors over fq whose sum is 1
-    in S.toList $ S.fromList [foldl1 (<+>) $ zipWith (*>) m ps | m <- multipliers]
-    where n = length $ head ps -- the dimension of the space we're working in
-          k = length ps        -- the dimension of the flat
+    let vs = [ (1 - sum xs) : xs | xs <- ptsAG (k-1) fq ] -- k-vectors over fq whose sum is 1
+    in toListSet [m' <<*> v | v <- vs]
+    where k = length ps -- the dimension of the flat (assuming ps are independent)
+          m' = L.transpose ps
           fq = eltsFq undefined
+-- toListSet call sorts the result, and also removes duplicates in case the points weren't independent
 
+{-
+closureAG ps =
+    let vs = [ (1 - sum xs) : xs | xs <- ptsAG (k-1) fq ] -- k-vectors over fq whose sum is 1
+    in toListSet [foldl1 (<+>) $ zipWith (*>) v ps | v <- vs]
+    where k = length ps        -- the dimension of the flat
+          fq = eltsFq undefined
+-}
+
 lineAG [p1,p2] = L.sort [ p1 <+> (c *> dp) | c <- fq ] where
     dp = p2 <-> p1
     fq = eltsFq undefined
@@ -51,11 +65,13 @@ -- closure of points in PG(n,Fq)
 -- take all linear combinations of the points (ie the subspace generated by the points, considered as points in Fq ^(n+1) )
 -- then discard all which aren't in PNF (thus dropping back into PG(n,Fq))
-closurePG ps = L.sort $ filter ispnf $ map (<*>> ps) $ ptsAG k fq where
+closurePG ps = toListSet $ filter ispnf $ map (<*>> ps) $ ptsAG k fq where
     k = length ps
     fq = eltsFq undefined
-
+-- toListSet call sorts the result, and also removes duplicates in case the points weren't independent
 
+linePG [p1,p2] = toListSet $ filter ispnf [(a *> p1) <+> (b *> p2) | a <- fq, b <- fq]
+    where fq = eltsFq undefined
 
 -- van Lint & Wilson, p325, 332
 qtorial n q | n >= 0 = product [(q^i - 1) `div` (q-1) | i <- [1..n]]
@@ -136,13 +152,12 @@ -- among all pairs of distinct points, select those which are the first two in the line they generate
 linesAG1 n fq = [ [x,y] | [x,y] <- combinationsOf 2 (ptsAG n fq),
                           [x,y] == take 2 (lineAG [x,y]) ]
-
+-- the point of the condition at the end is to avoid listing the same line more than once
 
 -- almost certainly not as efficient as linesAG, because requires lineAG/closureAG call
 -- a line in AG(n,fq) is a translation (x) of a line through the origin (y)
 linesAG2 n fq = [ [x,z] | x <- ptsAG n fq, y <- ptsPG (n-1) fq,
                           z <- [x <+> y], [x,z] == take 2 (closureAG [x,z]) ]
--- the point of the condition at the end is to avoid listing the same line more than once
 
 
 -- INCIDENCE GRAPH
@@ -164,7 +179,6 @@ -- The full group is called PGammaL(3,f4)
 
 
-
 -- |Incidence graph of AG(n,fq), considered as an incidence structure between points and lines
 incidenceGraphAG :: (Ord a, FiniteField a) => Int -> [a] -> Graph (Either [a] [[a]])
 incidenceGraphAG n fq = G vs es where
@@ -185,8 +199,13 @@ 
 orderAff n q = q^n * orderGL n q
 
+orderPGL n q = orderGL n q `div` (q-1)
 
 -- NOTE:
 -- AG(n,F2) is degenerate:
 -- Every pair of points is a line, so it is the complete graph on 2^n points
 -- And as such has aut group S(2^n)
+
+
+-- Heawood graph = incidence graph of Fano plane
+heawood = to1n $ incidenceGraphPG 2 f2

@@ -211,6 +211,8 @@ 
 valencies g@(G vs es) = map (head &&& length) $ L.group $ L.sort $ map (valency g) vs
 
+valencyPartition g@(G vs es) = map (map snd) $ L.groupBy (\x y -> fst x == fst y) [(valency g v, v) | v <- vs]
+
 regularParam g =
     case valencies g of
     [(v,_)] -> Just v

@@ -5,6 +5,7 @@ import qualified Data.List as L
 import qualified Data.Map as M
 import qualified Data.Set as S
+import Data.Maybe
 
 import Math.Common.ListSet
 import Math.Combinatorics.Graph
@@ -87,14 +88,22 @@           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 = concat [ [c1 `intersect` c2 | c2 <- p2] | c1 <- p1]
+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]
 
+
 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
@@ -130,8 +139,8 @@ -- Now using distance partitions
 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)
+        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]
@@ -150,8 +159,8 @@                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))
+                                   (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]
@@ -163,18 +172,17 @@ 
 -- 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
--- |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 g@(G vs es) = graphAuts' [] [vs] where
+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 $ 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)
+        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
@@ -191,8 +199,8 @@                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))
+                                   (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]
@@ -201,27 +209,129 @@ -- 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
 
 
+
+-- 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]
+
+
+-- 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
+
+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
+
+-- Now, every time we intersect two partitions, refine to an equitable partition
+-- |Given a graph g, graphAuts g returns a strong generating set for the automorphism group of g.
+graphAuts :: (Ord a) => Graph a -> [Permutation a]
+graphAuts g@(G vs es) = graphAuts' [] (toEquitable g $ valencyPartition g) 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]
+
+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']
+
+
 -- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH
 
 -- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
 incidenceAuts g@(G vs es) = map points (incidenceAuts' [] [vs]) where
     points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
-    -- points (P m) = P $ M.fromList [(k,v) | (Left k, Left v) <- M.toList m]
     incidenceAuts' us p@((x@(Left _):ys):pt) =
-        let p' = L.sort $ filter (not . null) $ refine (ys:pt) (dps M.! x)
+        -- 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)
+        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
+        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]
+
+
+-- GRAPH ISOMORPHISMS
+
+-- !! not yet using equitable partitions, so could probably be more efficient
+
+graphIsos g1 g2 = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2) | v2 <- vertices g2] where
+    v1 = head $ vertices g1
     dfs xys p1 p2
         | map length p1 /= map length p2 = []
         | otherwise =
@@ -229,21 +339,23 @@                 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 []
+                    in if isCompatible xys' then [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))
+                                   (refine' (xs : p1'') (dps1 M.! x))
+                                   (refine' ((L.delete y ys):p2'') (dps2 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
-
-
+    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]
+    es1 = S.fromList $ edges g1
+    es2 = S.fromList $ edges g2
 
+isIso g1 g2 = (not . null) (graphIsos g1 g2)
 
 
+{-
 removeGens x gs = removeGens' [] gs where
     baseOrbit = x .^^ gs
     removeGens' ls (r:rs) =
@@ -266,15 +378,15 @@     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)
+            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)
+        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
@@ -294,66 +406,12 @@                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))
+                                   (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
-
-
--- GRAPH ISOMORPHISMS
-
-graphIsos g1 g2 = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2) | v2 <- vertices g2] where
-    v1 = head $ vertices g1
-    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 [xys'] else []
-               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]
-    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]
-    es1 = S.fromList $ edges g1
-    es2 = S.fromList $ edges g2
-
-isIso g1 g2 = (not . null) (graphIsos g1 g2)
-
--- graphAuts3 g = map fromPairs $ graphIsos g g
-
-
-
-
-{-
-graphAuts2 (G vs es) = graphAuts' [] 1 (bsgsSym vs) where
-    graphAuts' bs g ((b,t):bts) = concat [graphAuts' (b:bs) (h*g) bts | h <- M.elems t, isCompatible (b:bs) (h*g)]
-    -- has to be h*g not g*h - not quite sure why
-    graphAuts' _ g [] = [g]
-    isCompatible (b:bs) g = and [(e `S.member` es') == ((e -^ g) `S.member` es') | e <- [ [b',b] | b' <- bs] ] -- if bs ordered then b' < b
-    es' = S.fromList es
-
-graphAutsSGS2 (G vs es) = transversals [] (bsgsSym vs) where
-    transversals bs ((b,t):bts) = let t' = concat [take 1 $ dfs (b:bs) h bts | h <- tail (M.elems t), isCompatible (b:bs) h]
-                                  in t' ++ transversals (b:bs) bts
-    transversals _ [] = []
-    dfs bs g ((b,t):bts) = concat [dfs (b:bs) (h*g) bts | h <- M.elems t, isCompatible (b:bs) (h*g)]
-    dfs _ g [] = [g]
-    isCompatible (b:bs) g = and [(e `S.member` es') == ((e -^ g) `S.member` es') | e <- [ [b',b] | b' <- bs] ] -- if bs ordered then b' < b
-    es' = S.fromList es
-
--- base and strong generating set for Sym(xs)
-bsgsSym xs = [(x, t x) | x <- init xs]
-    where t x = M.fromList $ (x,p []) : [(y, p [[x,y]]) | y <- dropWhile (<= x) xs]
-
-bsgs_S n = bsgsSym [1..n]
 -}
+
 

@@ -6,14 +6,10 @@ 
 import qualified Data.List as L
 
--- import PermutationGroup as P
 import Math.Algebra.Field.Base
 import Math.Algebra.Field.Extension
--- import Graph as G
--- import StronglyRegularGraph as SRG
 import Math.Combinatorics.Design as D
 import Math.Algebra.Group.SchreierSims as SS
--- import LinearAlgebra hiding ( (^-) )
 
 
 factorial n = product [1..n]
@@ -30,7 +26,12 @@     ,designParams (pg2 f2) == Just (2,(7,3,1))
     ,designParams (pg2 f3) == Just (2,(13,4,1))
     ,designParams (pg2 f4) == Just (2,(21,5,1))
+    ,designParams s_3_6_22 == Just (3,(22,6,1))
+    ,designParams (derivedDesign s_3_6_22 1) == Just (2,(21,5,1)) -- it is PG(2,F4)
+    ,designParams s_4_5_11 == Just (4,(11,5,1))
+    ,designParams (derivedDesign (derivedDesign s_4_5_11 0) 1) == Just (2,(9,3,1)) -- it is AG(2,F3)
     ]
+
 
 designAutTest = all (uncurry (==)) designAutTests
 

@@ -23,7 +23,10 @@     ,numFlatsPG 3 4 1 == length (flatsPG 3 f4 1)
     ,numFlatsPG 3 4 2 == length (flatsPG 3 f4 2)
     ,numFlatsPG 3 4 3 == length (flatsPG 3 f4 3)
-    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f2) == orderAff 2 2 * degree f2 
-    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f3) == orderAff 2 3 * degree f3 
-    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f4) == orderAff 2 4 * degree f4 
+    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f2) == orderAff 2 2 * toInteger (degree f2) 
+    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f3) == orderAff 2 3 * toInteger (degree f3)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphAG 2 f4) == orderAff 2 4 * toInteger (degree f4)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f2) == orderPGL 3 2 * toInteger (degree f2)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f3) == orderPGL 3 3 * toInteger (degree f3)
+    ,(orderSGS $ incidenceAuts $ incidenceGraphPG 2 f4) == orderPGL 3 4 * toInteger (degree f4)
     ]

@@ -60,6 +60,7 @@ -- by virtue of not even being vertex- or edge-transitive
 -- It is actually rather hard to find graphs which are vertex- and edge-transitive but not arc-transitive, but here is one
 -- Doyle, "A 27-vertex graph that is vertex-transitive and edge-transitive but not 1-transitive"
+-- http://en.wikipedia.org/wiki/Holt_graph
 doyleGraph = G gs es where
     relations = knuthBendix [("aaaaaaaaa",""), ("ccccccccc",""), ("aaaaaa","ccc"), ("cccccc","aaa"), ("ccccccccac","aaaa"), ("aaaaaaaaca","cccc")]
     gs = L.sort $ nfs ("ac",relations) -- all elements of the group, reduced to normal form