packages feed

HaskellForMaths (empty) → 0.1

raw patch · 28 files changed

+4894/−0 lines, 28 filesdep +basedep +containerssetup-changed

Dependencies added: base, containers

Files

+ HaskellForMaths.cabal view
@@ -0,0 +1,24 @@+   Name:                HaskellForMaths
+   Version:             0.1
+   Description:         Math library - combinatorics, group theory, commutative algebra, non-commutative algebra
+   License:             BSD3
+   License-file:        license.txt
+   Author:              David Amos
+   Maintainer:          haskellformaths-at-googlemail-dot-com
+   Homepage:            http://www.polyomino.f2s.com/haskellformathsv2/HaskellForMathsv2.html
+   Build-Type:          Simple
+   Cabal-Version:       >=1.2
+   
+   Library
+     Build-Depends:     base >=3 && < 4, containers
+     Exposed-modules:
+        Math.Algebra.LinearAlgebra, Math.Algebra.Commutative.MPoly, Math.Algebra.Commutative.GBasis,
+        Math.Algebra.Field.Base, Math.Algebra.Field.Extension,
+        Math.Algebra.Group.PermutationGroup, Math.Algebra.Group.SchreierSims, Math.Algebra.Group.StringRewriting,
+        Math.Algebra.NonCommutative.NCPoly, Math.Algebra.NonCommutative.GSBasis, Math.Algebra.NonCommutative.TensorAlgebra,
+        Math.Combinatorics.Graph, Math.Combinatorics.GraphAuts, Math.Combinatorics.StronglyRegularGraph,
+        Math.Combinatorics.Design, Math.Combinatorics.FiniteGeometry, Math.Combinatorics.Hypergraph,
+        Math.Common.IntegerAsType, Math.Common.ListSet,
+        Math.Projects.RootSystem, Math.Projects.ChevalleyGroup.Classical, Math.Projects.ChevalleyGroup.Exceptional,
+        Math.Projects.KnotTheory.LaurentMPoly, Math.Projects.KnotTheory.TemperleyLieb, Math.Projects.KnotTheory.IwahoriHecke
+     ghc-options:       -w
+ Math/Algebra/Commutative/GBasis.hs view
@@ -0,0 +1,267 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+module Math.Algebra.Commutative.GBasis where
+
+import Data.List
+import qualified Data.Map as M
+
+import Math.Algebra.Commutative.Monomial
+import Math.Algebra.Commutative.MPoly
+
+-- Sources:
+-- Cox, Little, O'Shea: Ideals, Varieties and Algorithms
+-- Giovini, Mora, Niesi, Robbiano, Traverso, "One sugar cube please, or Selection strategies in the Buchberger algorithm"
+
+
+
+sPoly f g = let h = lcmT (lt f) (lt g)
+            in h `divT` lt f .* f - h `divT` lt g .* g
+-- The point about the s-poly is that it cancels out the leading terms of the two polys, exposing their second terms
+
+isGB fs = all (\h -> h %% fs == 0) (pairWith sPoly fs)
+
+
+-- Cox p87
+gb1 fs = gb' fs (pairWith sPoly fs) where
+    gb' gs (h:hs) = let h' = h %% gs in
+                        if h' == 0 then gb' gs hs else gb' (h':gs) (hs ++ map (sPoly h') gs)
+    gb' gs [] = gs
+
+-- [f xi xj | xi <- xs, xj <- xs, i < j]
+pairWith f (x:xs) = map (f x) xs ++ pairWith f xs
+pairWith _ [] = []
+
+
+-- Cox p89-90
+reduce gs = reduce' [] gs where
+    reduce' gs' (g:gs) | g' == 0   = reduce' gs' gs
+                       | otherwise = reduce' (g':gs') gs
+                       where g' = g %% (gs'++gs)
+    reduce' gs' [] = reverse $ sort $ map toMonic gs'
+-- the reverse means that when using an elimination order, the elimination ideal will be at the end
+
+
+-- Cox et al p106-7
+-- No need to calculate an spoly fi fj if
+-- 1. the lm fi and lm fj are coprime, or
+-- 2. there exists some fk, with (i,k) (j,k) already considered, and lm fk divides lcm (lm fi) (lm fj) 
+-- some slight inefficiencies from looking up fi, fj repeatedly
+gb2 fs = reduce $ gb' fs (pairs [1..s]) s where
+    s = length fs
+    gb' gs ((i,j):ps) t =
+        let fi = gs!i; fj = gs!j in
+        if coprimeM (lm fi) (lm fj) || criterion fi fj
+        -- if lcmM (lm fi) (lm fj) == lm fi * lm fj || criterion fi fj
+        then gb' gs ps t
+        else let h = sPoly fi fj %% gs in
+             if h == 0 then gb' gs ps t else gb' (gs++[h]) (ps ++ [ (i,t+1) | i <- [1..t] ]) (t+1)
+        where
+            criterion fi fj = let l = lcmM (lm fi) (lm fj) in any (test l) [1..t]
+            test l k = k `notElem` [i,j]
+                    && ordpair i k `notElem` ps
+                    && ordpair j k `notElem` ps
+                    && lm (gs!k) `dividesM` l
+    gb' gs [] _ = gs
+
+
+pairs (x:xs) = map (\y->(x,y)) xs ++ pairs xs
+pairs [] = []
+
+xs ! i = xs !! (i-1) -- in other words, index the list from 1 not 0
+
+ordpair x y | x < y     = (x,y)
+            | otherwise = (y,x)
+
+-- version of gb2 where we eliminate pairs as they're created, rather than as they're processed
+gb2b fs = reduce $ gb1' [] fs [] 0 where
+    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where
+        ps' = updatePairs gs ps f t
+    gb1' ls [] ps t = gb2' ls ps t
+    gb2' gs ((i,j):ps) t =
+        let h = sPoly (gs!i) (gs!j) %% gs in
+        if h == 0
+        then gb2' gs ps t
+        else let ps' = updatePairs gs ((i,j):ps) h t in gb2' (gs++[h]) ps' (t+1)
+    gb2' gs [] _ = gs
+    updatePairs gs ps f t =
+        [p | p@(i,j) <- ps,
+             not (lm f `dividesM` lcmM (lm (gs!i)) (lm (gs!j)))] 
+     ++ [ (i,t+1) | (gi,i) <- zip gs [1..t],
+                    not (coprimeM (lm gi) (lm f)),
+                    not (criterion (lcmM (lm gi) (lm f)) i) ]
+        where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` ps]
+
+
+-- Cox et al 108
+-- 1. list smallest fs first, as more likely to reduce
+-- 2. order the pairs with smallest lcm fi fj first ("normal selection strategy")
+gb3b fs =
+    let fs' = sort $ filter (/=0) fs
+    in reduce $ gb1' [] fs' [] 0 where
+    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where
+        ps' = updatePairs gs ps f t
+    gb1' ls [] ps t = gb2' ls ps t
+    gb2' gs ((l,(i,j)):ps) t =
+        let h = sPoly (gs!i) (gs!j) %% gs in
+        if h == 0
+        then gb2' gs ps t
+        else let ps' = updatePairs gs ((l,(i,j)):ps) h t in gb2' (gs++[h]) ps' (t+1)
+    gb2' gs [] _ = gs
+    updatePairs :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [(Monomial ord, (Int,Int))] -> (MPoly ord r) -> Int -> [(Monomial ord, (Int,Int))]
+    updatePairs gs ps f t =
+        mergeBy cmpFst
+        [p | p@(l,(i,j)) <- ps,
+             not (lm f `dividesM` l)]
+        $ sortBy cmpFst
+        [ (l,(i,t+1)) | (gi,i) <- zip gs [1..t], l <- [lcmM (lm gi) (lm f)],
+                    not (coprimeM (lm gi) (lm f)),
+                    not (criterion l i) ]
+        where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` map snd ps]
+
+cmpFst (a,_) (b,_) = compare a b
+
+mergeBy cmp (t:ts) (u:us) =
+    case cmp t u of
+    LT -> t : mergeBy cmp ts (u:us)
+    EQ -> t : mergeBy cmp ts (u:us)
+    GT -> u : mergeBy cmp (t:ts) us
+mergeBy _ ts us = ts ++ us -- one of them is null
+
+
+-- naive implementation of "sugar strategy"
+gb4b fs =
+    let fs' = sort $ filter (/=0) fs
+    in reduce $ gb1' [] fs' [] 0 where
+    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where
+        ps' = updatePairs gs ps f t
+    gb1' ls [] ps t = gb2' ls ps t
+    gb2' gs ((sl,(i,j)):ps) t =
+        let h = sPoly (gs!i) (gs!j) %% gs in
+        if h == 0
+        then gb2' gs ps t
+        else let ps' = updatePairs gs ((sl,(i,j)):ps) h t in gb2' (gs++[h]) ps' (t+1)
+    gb2' gs [] _ = gs
+    updatePairs :: (Ord (Monomial ord), Fractional r) =>
+                   [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> (MPoly ord r) -> Int -> [((Int,Monomial ord), (Int,Int))]
+    updatePairs gs ps f t =
+        mergeBy cmpFst
+        [p | p@((s,l),(i,j)) <- ps,
+             not (lm f `dividesM` l)]
+        $ sortBy cmpFst
+        [ ((s,l),(i,t+1)) | (gi,i) <- zip gs [1..t], l <- [lcmM (lm gi) (lm f)], s <- [sugar gi f l],
+                    not (coprimeM (lm gi) (lm f)),
+                    not (criterion l i) ]
+        where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` map snd ps]
+
+
+-- Giovini et al
+-- The point of sugar is, given fi, fj, to give an upper bound on the degree of sPoly fi fj without having to calculate it
+-- We can then select by preference pairs with lower sugar, expecting therefore that the s-polys will have lower degree
+gb fs =
+    -- let fs' = sort $ filter (/=0) fs
+    let fs' = sort $ map toMonic $ filter (/=0) fs
+    in reduce $ gb1' [] fs' [] 0 where
+    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1)
+        where ps' = updatePairs gs ps f (t+1)
+    gb1' ls [] ps t = gb2' ls ps t
+    gb2' gs (p@(_,(i,j)):ps) t =
+        if h == 0
+        then gb2' gs ps t
+        else gb2' (gs++[h]) ps' (t+1)
+        where h = toMonic $ sPoly (gs!i) (gs!j) %% gs
+              ps' = updatePairs gs (p:ps) h (t+1)
+    gb2' gs [] _ = gs
+    updatePairs :: (Ord (Monomial ord), Fractional r) =>
+                   [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> (MPoly ord r) -> Int -> [((Int,Monomial ord), (Int,Int))]
+    updatePairs gs ps gk k =
+        let newps = [let l = lcmM (lm gi) (lm gk) in ((sugar gi gk l, l), (i,k)) | (gi,i) <- zip gs [1..k-1] ]
+            ps' = [p | p@((sij,tij),(i,j)) <- ps, ((sik,tik),_) <- [newps ! i], ((sjk,tjk),_) <- [newps ! j],
+                       not ( (tik `properlyDividesM` tij) && (tjk `properlyDividesM` tij) ) ] -- sloppy variant
+            newps' = discard1 [] newps
+            newps'' = sortBy cmpSug $ discard2 [] $ sortBy cmpNormal newps'
+        in mergeBy cmpSug ps' newps''
+        where
+            discard1 ls (r@((_sik,tik),(i,_k)):rs) =
+                if lm (gs!i) `coprimeM` lm gk
+                -- then discard [l | l@((_,tjk),_) <- ls, tjk /= tik] [r | r@((_,tjk),_) <- ls, tjk /= tik]
+                then discard1 (filter (\((_,tjk),_) -> tjk /= tik) ls) (filter (\((_,tjk),_) -> tjk /= tik) rs)
+                else discard1 (r:ls) rs
+            discard1 ls [] = ls
+            discard2 ls (r@((_sik,tik),(i,k)):rs) = discard2 (r:ls) $ filter (\((_sjk,tjk),(j,k')) -> not (k == k' && tik `dividesM` tjk)) rs
+            discard2 ls [] = ls
+-- The type annotation on updatePairs appears to be required
+-- The two calls to toMonic are designed to prevent coefficient explosion, but it is unproven that they are effective
+
+-- sugar of sPoly f g, where h = lcm (lt f) (lt g)
+-- this is an upper bound on deg (sPoly f g)
+sugar f g h = degM h + max (deg f - degM (lm f)) (deg g - degM (lm g))
+
+cmpNormal ((s1,t1),(i1,j1)) ((s2,t2),(i2,j2)) = compare (t1,j1) (t2,j2)
+
+cmpSug ((s1,t1),(i1,j1)) ((s2,t2),(i2,j2)) = compare (s1,t1,j1) (s2,t2,j2)
+
+
+-- earlier version of gb3b
+gb3 fs =
+    let gs = sort $ filter (/=0) fs
+        ps = sortBy cmpFst $ pairWith (\(i,fi) (j,fj) -> (lcmM (lm fi) (lm fj), (i,j)) ) $ zip [1..] gs
+    in reduce $ gb' gs ps s
+    where
+        s = length fs
+        gb' :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [(Monomial ord, (Int,Int))] -> Int -> [MPoly ord r]
+        gb' gs ((l,(i,j)):ps) t =
+            let fi = gs!i; fj = gs!j in
+            if coprimeM (lm fi) (lm fj) || any (test l) [1..t]
+            then gb' gs ps t
+            else let h = sPoly fi fj %% gs in
+                 if h == 0
+                 then gb' gs ps t
+                 else let ps' = mergeBy cmpFst ps $ sortBy cmpFst $ zip [lcmM (lm h) (lm fi) | fi <- gs] [(i,t+1) | i <- [1..t]]
+                 -- else let ps' = mergeBy cmpFst ps $ zip [lcmM (lm h) (lm fi) | fi <- gs] [(i,t+1) | i <- [1..t]]
+                      in gb' (gs++[h]) ps' (t+1)
+            where
+                test l k = k `notElem` [i,j]
+                        && ordpair i k `notElem` map snd ps
+                        && ordpair j k `notElem` map snd ps
+                        && lm (gs!k) `dividesM` l
+        gb' gs [] _ = gs
+-- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm
+
+
+
+-- earlier version of gb4b
+gb4 fs =
+    let gs = sort $ filter (/=0) fs
+        ps = sortBy cmpFst $ pairWith (\(i,fi) (j,fj) -> let l = lcmM (lm fi) (lm fj) in ((sugar fi fj l, l), (i,j)) ) $ zip [1..] gs
+    in reduce $ gb' gs ps s
+    where
+        s = length fs
+        gb' :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> Int -> [MPoly ord r]
+        gb' gs (((s,l),(i,j)):ps) t =
+            let fi = gs!i; fj = gs!j in
+            if coprimeM (lm fi) (lm fj) || any (test l) [1..t]
+            then gb' gs ps t
+            else let h = sPoly fi fj %% gs in
+                 if h == 0
+                 then gb' gs ps t
+                 else let ps' = mergeBy cmpFst ps $ sortBy cmpFst $ zip [let l = lcmM (lm fi) (lm h) in (sugar fi h l, l) | fi <- gs] [(i,t+1) | i <- [1..t]]
+                      in gb' (gs++[h]) ps' (t+1)
+            where
+                test l k = k `notElem` [i,j]
+                        && ordpair i k `notElem` map snd ps
+                        && ordpair j k `notElem` map snd ps
+                        && lm (gs!k) `dividesM` l
+        gb' gs [] _ = gs
+-- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm
+
+
+
+{-
+merge (t:ts) (u:us) =
+    if t <= u
+    then t : merge ts (u:us)
+    else u : merge (t:ts) us
+merge ts us = ts ++ us -- one of them is null
+-}
+ Math/Algebra/Commutative/MPoly.hs view
@@ -0,0 +1,197 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+module Math.Algebra.Commutative.MPoly where
+
+import qualified Data.Map as M
+import Data.List as L
+import Control.Arrow (first, second)
+import Data.Ratio (denominator)
+
+import Math.Algebra.Field.Base
+import Math.Algebra.Commutative.Monomial
+
+
+-- MULTIVARIATE POLYNOMIALS
+
+newtype MPoly ord r = MP [(Monomial ord,r)] deriving (Eq)
+-- deriving instance (Ord (Monomial ord), Ord r) => Ord (MPoly ord r)
+-- standalone deriving supported from GHC 6.8
+
+instance (Ord (Monomial ord), Ord r) => Ord (MPoly ord r) where
+    compare (MP ts) (MP us) = compare ts us
+
+instance (Show r, Num r) => Show (MPoly ord r) where
+    show (MP []) = "0"
+    show (MP ts) =
+        let (c:cs) = concatMap showTerm ts
+        in if c == '+' then cs else c:cs
+        where showTerm (m,c) =
+                  case show c of
+                  "1" -> "+" ++ show m
+                  "-1" -> "-" ++ show m
+                  cs@(x:_) -> (if x == '-' then cs else '+':cs) ++ (if m == 1 then "" else show m)
+
+
+instance (Ord (Monomial ord), Num r) => Num (MPoly ord r) where
+    MP ts + MP us = MP (mergeTerms ts us)
+    negate (MP ts) = MP $ map (second negate) ts
+    MP ts * MP us = MP $ collect $ sortBy cmpTerm $ [(g*h,c*d) | (g,c) <- ts, (h,d) <- us]
+    {-
+    -- The following appears to be slightly slower, perhaps because sortBy is compiled
+    MP (t@(g,c):ts) * MP (u@(h,d):us) =
+        let MP vs = MP ts * MP us
+        in MP $ mergeTerms ((g*h,c*d):vs) $ mergeTerms [(g*h,c*d) | (h,d) <- us] [(g*h,c*d) | (g,c) <- ts]
+    _ * _ = MP []
+    -}
+    fromInteger 0 = MP []
+    fromInteger n = MP [(fromInteger 1, fromInteger n)]
+
+cmpTerm (a,c) (b,d) = case compare a b of EQ -> EQ; GT -> LT; LT -> GT -- in mpolys we put "larger" terms first
+
+-- inputs in descending order
+mergeTerms (t@(g,c):ts) (u@(h,d):us) =
+    case compare g h of
+    GT -> t : mergeTerms ts (u:us)
+    LT -> u : mergeTerms (t:ts) us
+    EQ -> if e == 0 then mergeTerms ts us else (g,e) : mergeTerms ts us
+    where e = c + d
+mergeTerms ts us = ts ++ us -- one of them is null
+
+collect (t1@(g,c):t2@(h,d):ts)
+    | g == h = collect $ (g,c+d):ts
+    | c == 0  = collect $ t2:ts
+    | otherwise = t1 : collect (t2:ts)
+collect ts = ts
+
+-- Fractional instance so that we can enter fractional coefficients
+-- Only lets us divide by field elements (with unit monomial), not any other polynomials
+instance (Ord (Monomial ord), Fractional r) => Fractional (MPoly ord r) where
+    recip (MP [(m,c)]) = MP [(recip m, recip c)]
+    -- recip (MP [(m,c)]) | m == fromInteger 1 = MP [(m, recip c)]
+    recip _ = error "MPoly.recip: only supported for (non-zero) constants or monomials"
+
+
+var v = MP [(Monomial $ M.singleton v 1, 1)] :: MPoly Grevlex Q
+a = var "a"
+b = var "b"
+c = var "c"
+d = var "d"
+s = var "s"
+t = var "t"
+u = var "u"
+v = var "v"
+w = var "w"
+x = var "x"
+y = var "y"
+z = var "z"
+
+x_ i = var ("x" ++ show i)
+x0 = x_ 0
+x1 = x_ 1
+x2 = x_ 2
+x3 = x_ 3
+
+
+-- convertMP :: Ord (Monomial ord') => MPoly ord k -> MPoly ord' k
+convertMP (MP ts) = MP $ sortBy cmpTerm $ map (first convertM) ts
+
+toLex :: MPoly ord k -> MPoly Lex k
+toLex = convertMP
+
+toGlex :: MPoly ord k -> MPoly Glex k
+toGlex = convertMP
+
+toGrevlex :: MPoly ord k -> MPoly Grevlex k
+toGrevlex = convertMP
+
+toElim :: MPoly ord k -> MPoly Elim k
+toElim = convertMP
+
+
+varLex v = toLex $ var v
+
+varElim v = toElim $ var v
+
+
+-- DIVISION ALGORITHM
+
+lt (MP (t:ts)) = t
+
+lm = fst . lt
+
+deg 0 = -1
+deg (MP ts) = maximum [degM m | (m,c) <- ts]
+-- the true degree of the polynomial, not the degree of the leading term
+-- required for sugar strategy when computing Groebner basis
+
+mulT (m,c) (m',c') = (m*m',c*c')
+
+divT (m,c) (m',c') = (m/m',c/c')
+
+dividesT (m,_) (m',_) = dividesM m m'
+
+properlyDividesT (m,_) (m',_) = dividesM m m' && m /= m'
+
+lcmT (m,c) (m',c') = (lcmM m m',1)
+
+
+infixl 7 .*
+t .* MP ts = MP $ map (mulT t) ts -- preserves term order
+
+
+-- given f, gs, find as, r such that f = sum (zipWith (*) as gs) + r, with r not divisible by any g
+quotRemMP f gs = quotRemMP' f (replicate n 0, 0) where
+    n = length gs
+    quotRemMP' 0 (us,r) = (us,r)
+    quotRemMP' h (us,r) = divisionStep h (gs,[],us,r)
+    divisionStep h (g:gs,us',u:us,r) =
+        if lt g `dividesT` lt h
+        then let t = MP [lt h `divT` lt g]
+                 h' = h - t*g
+                 u' = u+t
+             in quotRemMP' h' (reverse us' ++ u':us, r)
+        else divisionStep h (gs,u:us',us,r)
+    divisionStep h ([],us',[],r) =
+        let (lth,h') = splitlt h
+        in quotRemMP' h' (reverse us', r+lth)
+    splitlt (MP (t:ts)) = (MP [t], MP ts)
+
+infixl 7 %%
+f %% gs = r where (_,r) = quotRemMP f gs
+
+-- div and mod by single mpoly
+divModMP f g = (q,r) where ([q],r) = quotRemMP f [g]
+
+divMP f g = q where ([q],_) = quotRemMP f [g]
+
+modMP f g = r where (_,r) = quotRemMP f [g]
+
+
+-- OTHER STUFF
+
+-- injection of field elements into polynomial ring
+inject 0 = MP []
+inject c = MP [(fromInteger 1, c)]
+
+toMonic 0 = 0
+toMonic (MP ts@((_,c):_))
+    | c == 1 = MP ts
+    | otherwise = MP $ map (second (/c)) ts
+
+-- multiply out all denominators
+toZ (MP ts) = MP $ map (second (*c)) ts where c = fromInteger $ foldl lcm 1 $ [denominator c | (m,Q c) <- ts]
+
+-- substitute terms for variables in an MPoly
+-- eg subst [(x,a),(y,a+b),(z,c^2)] (x*y+z) -> a*(a+b)+c^2
+subst vts (MP us) = sum [inject c * substM m | (m,c) <- us] where
+    substM (Monomial m) = product [substV v ^ i | (v,i) <- M.toList m]
+    substV v =
+        let v' = MP [(Monomial $ M.singleton v 1, 1)] in
+        case L.lookup v' vts of
+        Just t -> t
+        Nothing -> v' -- no substitute, so keep as is
+
+support (MP ts) = [MP [(m,1)] | m <- reverse $ L.sort $ support' ts]
+    where support' ts = foldl L.union [] [supportM m | (m,c) <- ts]
+ Math/Algebra/Field/Base.hs view
@@ -0,0 +1,147 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+module Math.Algebra.Field.Base where
+
+import Data.Ratio
+import Math.Common.IntegerAsType
+
+
+-- RATIONALS
+
+-- Rationals with a better show function
+newtype Q = Q Rational deriving (Eq,Ord,Num,Fractional)
+
+instance Show Q where
+    show (Q x) | b == 1    = show a
+               | otherwise = show a ++ "/" ++ show b
+               where a = numerator x
+                     b = denominator x
+
+
+numeratorQ (Q x) = Data.Ratio.numerator x
+denominatorQ (Q x) = Data.Ratio.denominator x
+
+
+-- PRIME FIELDS
+
+-- extendedEuclid a b returns (u,v,d) such that u*a + v*b = d
+extendedEuclid a b | a >= 0 && b >= 0 = extendedEuclid' a b [] where
+    extendedEuclid' d 0 qs = let (u,v) = unwind 1 0 qs in (u,v,d)
+    extendedEuclid' a b qs = let (q,r) = quotRem a b in extendedEuclid' b r (q:qs)
+    unwind u v [] = (u,v)
+    unwind u v (q:qs) = unwind v (u-v*q) qs
+
+
+newtype Fp n = Fp Integer deriving (Eq,Ord)
+
+instance Show (Fp n) where
+    show (Fp x) = show x
+
+instance IntegerAsType n => Num (Fp n) where
+    Fp x + Fp y = Fp $ (x+y) `mod` p where p = value (undefined :: n)
+    negate (Fp 0) = Fp 0
+    negate (Fp x) = Fp $ p - x where p = value (undefined :: n)
+    Fp x * Fp y = Fp $ (x*y) `mod` p where p = value (undefined :: n)
+    fromInteger m = Fp $ m `mod` p where p = value (undefined :: n)
+
+-- n must be prime - could perhaps use a type to guarantee this
+instance IntegerAsType n => Fractional (Fp n) where
+    recip 0 = error "Fp.recip 0"
+    recip (Fp x) = let (u,v,1) = extendedEuclid x p -- so ux+vp = 1. (We know the gcd is 1 as p prime)
+                   in Fp $ u `mod` p
+                   where p = value (undefined :: n)
+
+class FiniteField fq where
+    eltsFq :: fq -> [fq]  -- return all elts of the field
+    basisFq :: fq -> [fq] -- return an additive basis for the field (as Z-module)
+
+instance IntegerAsType p => FiniteField (Fp p) where
+    eltsFq _ = map fromInteger [0..p'-1] where p' = value (undefined :: p)
+    basisFq _ = [fromInteger 1]
+
+primitiveElt fq = head [x | x <- tail fq, length (powers x) == q-1] where
+    q = length fq
+
+powers x | x /= 0 = 1 : takeWhile (/=1) (iterate (*x) x)
+
+
+-- characteristic of a finite field
+char fq = head [p | p <- [2..], length fq `mod` p == 0]
+
+
+type F2 = Fp T2
+f2 = map fromInteger [0..1] :: [F2]
+
+type F3 = Fp T3
+f3 = map fromInteger [0..2] :: [F3]
+
+type F5 = Fp T5
+f5 = map fromInteger [0..4] :: [F5]
+
+type F7 = Fp T7
+f7 = map fromInteger [0..6] :: [F7]
+
+type F11 = Fp T11
+f11 = map fromInteger [0..10] :: [F11]
+
+type F13 = Fp T13
+f13 = map fromInteger [0..12] :: [F13]
+
+type F17 = Fp T17
+f17 = map fromInteger [0..16] :: [F17]
+
+type F19 = Fp T19
+f19 = map fromInteger [0..18] :: [F19]
+
+type F23 = Fp T23
+f23 = map fromInteger [0..22] :: [F23]
+
+type F29 = Fp T29
+f29 = map fromInteger [0..28] :: [F29]
+
+type F31 = Fp T31
+f31 = map fromInteger [0..30] :: [F31]
+
+type F37 = Fp T37
+f37 = map fromInteger [0..36] :: [F37]
+
+type F41 = Fp T41
+f41 = map fromInteger [0..40] :: [F41]
+
+type F43 = Fp T43
+f43 = map fromInteger [0..42] :: [F43]
+
+type F47 = Fp T47
+f47 = map fromInteger [0..46] :: [F47]
+
+type F53 = Fp T53
+f53 = map fromInteger [0..52] :: [F53]
+
+type F59 = Fp T59
+f59 = map fromInteger [0..58] :: [F59]
+
+type F61 = Fp T61
+f61 = map fromInteger [0..60] :: [F61]
+
+type F67 = Fp T67
+f67 = map fromInteger [0..66] :: [F67]
+
+type F71 = Fp T71
+f71 = map fromInteger [0..70] :: [F71]
+
+type F73 = Fp T73
+f73 = map fromInteger [0..72] :: [F73]
+
+type F79 = Fp T79
+f79 = map fromInteger [0..78] :: [F79]
+
+type F83 = Fp T83
+f83 = map fromInteger [0..82] :: [F83]
+
+type F89 = Fp T89
+f89 = map fromInteger [0..88] :: [F89]
+
+type F97 = Fp T97
+f97 = map fromInteger [0..96] :: [F97]
+ Math/Algebra/Field/Extension.hs view
@@ -0,0 +1,207 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+module Math.Algebra.Field.Extension where
+
+import Math.Common.IntegerAsType
+import Math.Algebra.Field.Base
+
+
+-- UNIVARIATE POLYNOMIALS
+
+newtype UPoly a = UP [a] deriving (Eq,Ord)
+-- the list [a_0, a_1, ..., a_n] represents the polynomial a_0 + a_1 x + ... + a_n x^n
+
+x = UP [0,1] :: UPoly Integer
+
+instance (Show a, Num a) => Show (UPoly a) where
+    show (UP []) = "0"
+    show (UP as) = let powers = filter ( (/=0) . fst ) $ zip as [0..]
+                       c:cs = concatMap showTerm powers
+                   in if c == '+' then cs else c:cs
+        where showTerm (a,i) = showCoeff a ++ showPower a i
+              showCoeff a = case show a of
+                            "1" -> "+"
+                            "-1" -> "-"
+                            '-':cs -> '-':cs
+                            cs -> '+':cs
+              showPower a i | i == 0 = case show a of
+                                       "1" -> "1"
+                                       "-1" -> "1"
+                                       otherwise -> ""
+                            | i == 1 = "x"
+                            | i > 1  = "x^" ++ show i
+
+instance Num a => Num (UPoly a) where
+	UP as + UP bs = toUPoly $ as <+> bs
+	negate (UP as) = UP $ map negate as
+	UP as * UP bs = toUPoly $ as <*> bs
+	fromInteger 0 = UP []
+	fromInteger a = UP [fromInteger a]
+
+toUPoly as = UP (reverse (dropWhile (== 0) (reverse as)))
+
+(a:as) <+> (b:bs) = (a+b) : (as <+> bs)
+as <+> [] = as
+[] <+> bs = bs
+
+[] <*> _ = []
+_ <*> [] = []
+(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
+-- > convert (x^2+3*x+2) :: UPoly F2
+-- x^2+x
+-- > convert (x^2+3*x+2) :: UPoly F3
+-- x^2+2
+
+
+-- DIVISION ALGORITHM
+
+-- degree
+deg (UP as) = length as
+
+-- leading term
+lt (UP as) = last as
+
+monomial a i = UP $ replicate i 0 ++ [a]
+
+-- quotRem for UPolys over a field
+quotRemUP :: (Num k, Fractional k) => UPoly k -> UPoly k -> (UPoly k, UPoly k)
+quotRemUP f g = qr 0 f where
+    qr q r = if deg r < deg_g
+             then (q,r)
+             else let s = monomial (lt r / lt_g) (deg r - deg_g)
+                  in qr (q+s) (r-s*g)
+    deg_g = deg g
+    lt_g = lt g
+
+
+modUP f g = snd $ quotRemUP f g
+
+-- extendedEuclidUP f g returns (u,v,d) such that u*f + v*g = d
+extendedEuclidUP f g = extendedEuclidUP' f g [] where
+    extendedEuclidUP' d 0 qs = let (u,v) = unwind 1 0 qs in (u,v,d)
+    extendedEuclidUP' f g qs = let (q,r) = quotRemUP f g in extendedEuclidUP' g r (q:qs)
+    unwind u v [] = (u,v)
+    unwind u v (q:qs) = unwind v (u-v*q) qs
+
+
+-- EXTENSION FIELDS
+
+class PolynomialAsType k poly where
+    pvalue :: (k,poly) -> UPoly k
+
+data ExtensionField k poly = Ext (UPoly k) deriving (Eq,Ord)
+
+instance Num k => Show (ExtensionField k poly) where
+    show (Ext f) = show f
+
+instance (Num k, Fractional k, PolynomialAsType k poly) => Num (ExtensionField k poly) where
+    Ext x + Ext y = Ext $ (x+y) `modUP` pvalue (undefined :: (k,poly))
+    Ext x * Ext y = Ext $ (x*y) `modUP` pvalue (undefined :: (k,poly))
+    negate (Ext x) = Ext $ negate x
+    fromInteger x = Ext $ fromInteger x
+
+instance (Num k, Fractional k, PolynomialAsType k poly) => Fractional (ExtensionField k poly) where
+    recip 0 = error "ExtensionField.recip 0"
+    recip (Ext f) = let g = pvalue (undefined :: (k,poly))
+                        (u,v,1) = extendedEuclidUP f g -- so u*f + v*g == 1. (We know the gcd is 1 as g is irreducible)
+                    in Ext $ u `modUP` g
+
+instance (Num k, FiniteField k, PolynomialAsType k poly) => FiniteField (ExtensionField k poly) where
+    eltsFq _ = map Ext (polys (d-1) fp) where
+        fp = eltsFq (undefined :: k)
+        d = deg $ pvalue (undefined :: (k,poly))
+    basisFq _ = map embed $ take (d-1) $ iterate (*x) 1 where
+        d = deg $ pvalue (undefined :: (k,poly))
+
+embed f = Ext (convert f)
+
+
+-- PRIME POWER FINITE FIELDS
+
+polys d fp = map toUPoly $ polys' d where
+    polys' 0 = [[]]
+    polys' d = [x:xs | x <- fp, xs <- polys' (d-1)]
+
+-- Conway polynomials from Holt, Handbook of Computational Group Theory, p60
+
+data ConwayF4
+instance PolynomialAsType F2 ConwayF4 where pvalue _ = convert $ x^2+x+1
+type F4 = ExtensionField F2 ConwayF4
+f4 = map Ext (polys 2 f2) :: [F4]
+x4 = embed x :: F4
+
+data ConwayF8
+instance PolynomialAsType F2 ConwayF8 where pvalue _ = convert $ x^3+x+1
+type F8 = ExtensionField F2 ConwayF8
+f8 = map Ext (polys 3 f2) :: [F8]
+x8 = embed x :: F8
+
+data ConwayF9
+instance PolynomialAsType F3 ConwayF9 where pvalue _ = convert $ x^2+2*x+2
+type F9 = ExtensionField F3 ConwayF9
+f9 = map Ext (polys 2 f3) :: [F9]
+x9 = embed x :: F9
+
+data ConwayF16
+instance PolynomialAsType F2 ConwayF16 where pvalue _ = convert $ x^4+x+1
+type F16 = ExtensionField F2 ConwayF16
+f16 = map Ext (polys 4 f2) :: [F16]
+x16 = embed x :: F16
+
+data ConwayF25
+instance PolynomialAsType F5 ConwayF25 where pvalue _ = convert $ x^2+4*x+2
+type F25 = ExtensionField F5 ConwayF25
+f25 = map Ext (polys 2 f5) :: [F25]
+x25 = embed x :: F25
+
+data ConwayF27
+instance PolynomialAsType F3 ConwayF27 where pvalue _ = convert $ x^3+2*x+1
+type F27 = ExtensionField F3 ConwayF27
+f27 = map Ext (polys 3 f3) :: [F27]
+x27 = embed x :: F27
+
+data ConwayF32
+instance PolynomialAsType F2 ConwayF32 where pvalue _ = convert $ x^5+x^2+1
+type F32 = ExtensionField F2 ConwayF32
+f32 = map Ext (polys 5 f2) :: [F32]
+x32 = embed x :: F32
+
+
+-- QUADRATIC EXTENSIONS OF Q
+
+data Sqrt a = Sqrt a
+
+-- n should be square-free
+instance IntegerAsType n => PolynomialAsType Q (Sqrt n) where
+    pvalue _ = convert $ x^2 - fromInteger (value (undefined :: n))
+
+type QSqrt2 = ExtensionField Q (Sqrt T2)
+sqrt2 = embed x :: QSqrt2
+
+type QSqrt3 = ExtensionField Q (Sqrt T3)
+sqrt3 = embed x :: QSqrt3
+
+type QSqrt5 = ExtensionField Q (Sqrt T5)
+sqrt5 = embed x :: QSqrt5
+
+type QSqrt7 = ExtensionField Q (Sqrt T7)
+sqrt7 = embed x :: QSqrt7
+
+type QSqrtMinus1 = ExtensionField Q (Sqrt TMinus1)
+i = embed x :: QSqrtMinus1
+
+
+type QSqrtMinus2 = ExtensionField Q (Sqrt (M TMinus1 T2))
+sqrtminus2 = embed x :: QSqrtMinus2
+
+type QSqrtMinus3 = ExtensionField Q (Sqrt (M TMinus1 T3))
+sqrtminus3 = embed x :: QSqrtMinus3
+
+type QSqrtMinus5 = ExtensionField Q (Sqrt (M TMinus1 T5))
+sqrtminus5 = embed x :: QSqrtMinus5
+ Math/Algebra/Group/PermutationGroup.hs view
@@ -0,0 +1,313 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Algebra.Group.PermutationGroup where
+
+import qualified Data.List as L
+import qualified Data.Map as M
+import qualified Data.Set as S
+
+import Math.Common.ListSet (toListSet, union) -- a version of union which assumes the arguments are ascending sets (no repeated elements)
+
+rotateL (x:xs) = xs ++ [x]
+
+
+-- PERMUTATIONS
+
+newtype Permutation a = P (M.Map a a) deriving (Eq,Ord)
+
+fromPairs xys | isValid   = fromPairs' xys
+              | otherwise = error "Not a permutation"
+    where (xs,ys) = unzip xys
+          (xs',ys') = (L.sort xs, L.sort ys)
+          isValid = xs' == ys' && all ((==1) . length) (L.group xs') -- ie the domain and range are the same, and are *sets*
+
+fromPairs' xys = P $ M.fromList $ filter (uncurry (/=)) xys
+-- we remove fixed points, so that the derived Eq instance works as expected
+
+toPairs (P g) = M.toList g
+
+fromList xs = fromPairs $ zip xs (L.sort xs)
+-- for example, fromList [2,3,1] is [[1,3,2]] - because the 1 moved to the 3 position
+
+-- the support of a permutation is the points it moves (returned in ascending order)
+supp (P g) = M.keys g
+-- (This is guaranteed not to contain fixed points provided the permutations have been constructed using the supplied constructors)
+
+-- image of x under action of g
+x .^ P g = case M.lookup x g of
+           Just y  -> y
+           Nothing -> x -- if x `notElem` supp (P g), then x is not moved
+
+-- construct a permutation from cycles
+fromCycles cs = fromPairs $ concatMap fromCycle cs
+    where fromCycle xs = zip xs (rotateL xs)
+
+-- as we will use fromCycles a lot, we provide a shorthand for it
+p cs = fromCycles cs
+-- can't specify in pointfree style because of monomorphism restriction
+
+-- convert a permutation to cycles
+toCycles g = toCycles' $ supp g
+    where toCycles' ys@(y:_) = let c = cycleOf g y in c : toCycles' (ys L.\\ c)
+          toCycles' [] = []
+
+cycleOf g x = cycleOf' x [] where
+    cycleOf' y ys = let y' = y .^ g in if y' == x then reverse (y:ys) else cycleOf' y' (y:ys)
+
+instance (Ord a, Show a) => Show (Permutation a) where
+    show g = show (toCycles g)
+
+parity g = let cs = toCycles g in (length (concat cs) - length cs) `mod` 2
+-- parity' g = length (filter (even . length) $ toCycles g) `mod` 2
+
+sign g = (-1)^(parity g)
+
+orderElt g = foldl lcm 1 $ map length $ toCycles g
+-- == order [g]
+
+instance (Ord a, Show a) => Num (Permutation a) where
+    g * h = fromPairs' [(x, x .^ g .^ h) | x <- supp g `union` supp h]
+    -- signum = sign -- doesn't work, complains about no (+) instance
+    fromInteger 1 = P $ M.empty
+
+inverse (P g) = P $ M.fromList $ map (\(x,y)->(y,x)) $ M.toList g
+
+g ^- n = inverse g ^ n
+
+instance (Ord a, Show a) => Fractional (Permutation a) where
+    recip = inverse
+
+-- conjugation
+h ~^ g = g^-1 * h * g
+
+-- commutator
+comm g h = g^-1 * h^-1 * g * h
+
+
+-- ORBITS
+
+-- action on blocks
+xs -^ g = L.sort [x .^ g | x <- xs]
+
+-- 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
+x .^^ gs = orbit (.^) x gs
+orbitP gs x = orbit (.^) x gs
+
+-- orbit of a block
+b -^^ gs = orbit (-^) b gs
+orbitB gs b = orbit (-^) b gs
+
+-- 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
+
+-- find all the orbits of a group
+-- (as we typically work with transitive groups, this is more useful for studying induced actions)
+-- (Note that of course this won't find orbits of points which are fixed by all elts of G)
+orbits gs = let xs = foldl union [] $ map supp gs in orbits' xs
+    where orbits' [] = []
+          orbits' (x:xs) = let o = x .^^ gs in o : orbits' (xs L.\\ o)
+
+
+-- GROUPS
+-- Some standard sequences of groups, and constructions of new groups from old
+
+-- Cn, cyclic group of order n
+_C n | n >= 2 = [p [[1..n]]]
+
+-- D2n, dihedral group of order 2n, symmetry group of n-gon
+-- For example, _D 8 == _D2 4 == symmetry group of square
+_D n | r == 0 = _D2 q where (q,r) = n `quotRem` 2
+
+_D2 n | n >= 3 = [a,b] where
+    a = p [[1..n]]                            -- rotation
+    b = fromPairs $ [(i,n+1-i) | i <- [1..n]] -- reflection
+
+-- Sn, symmetric group on [1..n]
+_S n | n >= 3 = [s,t] where
+    s = p [[1..n]]
+    t = p [[1,2]]
+
+-- An, alternating group on [1..n]
+_A n | n == 3 = [t]
+     | n > 3 = [s,t] where
+    s | odd n  = p [[3..n]]
+      | even n = p [[1,2], [3..n]]
+    t = p [[1,2,3]]
+
+
+-- Cartesian product of groups
+-- Given generators for H and K, acting on sets X and Y respectively,
+-- return generators for H*K, acting on the disjoint union X+Y (== Either X Y)
+cp hs ks =
+    [P $ M.fromList $ map (\(x,x') -> (Left x,Left x')) $ M.toList h' | P h' <- hs] ++
+    [P $ M.fromList $ map (\(y,y') -> (Right y,Right y')) $ M.toList k' | P k' <- ks]
+
+-- Wreath product of groups
+-- Given generators for H and K, acting on sets X and Y respectively,
+-- return generators for H wr K, acting on X*Y (== (X,Y))
+-- (Cameron, Combinatorics, p229-230; Cameron, Permutation Groups, p11-12)
+wr hs ks =
+    let _X = S.toList $ foldl S.union S.empty [M.keysSet h' | P h' <- hs] -- set on which H acts
+        _Y = S.toList $ foldl S.union S.empty [M.keysSet k' | P k' <- ks] -- set on which K acts
+        -- Then the wreath product acts on cartesian product X * Y,
+        -- regarded as a fibre bundle over Y of isomorphic copies of X
+        _B = [P $ M.fromList $ map (\(x,x') -> ((x,y),(x',y))) $ M.toList h' | P h' <- hs, y <- _Y]
+        -- bottom group B applies the action of H within each fibre
+        _T = [P $ M.fromList [((x,y),(x,y')) | x <- _X, (y,y') <- M.toList k'] | P k' <- ks]
+        -- top group T uses the action of K to permute the fibres
+    in _B ++ _T -- semi-direct product of B and T
+
+
+-- embed group elts into Sn - ie, convert so that the set acted on is [1..n]
+toSn gs = [toSn' g | g <- gs] where
+    _X = foldl union [] $ map supp gs   -- the set on which G acts
+    mapping = M.fromList $ zip _X [1..] -- the mapping from _X to [1..n]
+    toSn' g = fromPairs' $ map (\(x,x') -> (mapping M.! x, mapping M.! x')) $ toPairs g
+
+
+-- INVESTIGATING GROUPS
+-- Functions to investigate groups in various ways
+-- Most of these functions will only be efficient for small groups (say |G| < 10000)
+-- For larger groups we will need to use Schreier-Sims and associated algorithms
+
+-- Given generators, list all elements of a group
+-- Note, result is guaranteed to be in order, which we use on occasion
+elts gs = orbit (*) 1 gs
+
+eltsS gs = orbitS (*) 1 gs
+
+order gs = S.size $ eltsS gs -- length $ elts gs
+
+isMember gs h = h `S.member` eltsS gs -- h `elem` elts gs
+
+
+-- given the elts of a group, find generators
+gens hs = gens' [] (S.singleton 1) hs where
+    gens' gs eltsG (h:hs) = if h `S.member` eltsG then gens' gs eltsG hs else gens' (h:gs) (eltsS $ h:gs) hs
+    gens' gs _ [] = reverse gs
+
+
+
+-- conjClass gs h = orbit (~^) gs h
+
+-- Conjugacy class - should only be used for small groups
+h ~^^ gs = orbit (~^) h gs
+
+conjClass gs h = h ~^^ gs
+
+-- This is just the orbits under conjugation. Can we generalise "orbits" to help us here?
+conjClasses gs = conjClasses' (elts gs)
+    where conjClasses' [] = []
+          conjClasses' (h:hs) = let c = conjClass gs h in c : conjClasses' (hs L.\\ c)
+
+
+-- 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]
+
+-- the centre of G is the set of elts of G which commute with all other elts
+centre gs = centralizer gs gs
+
+-- normaliser of a subgroup
+-- the normaliser of H in G is {g <- G | g^-1Hg == H}
+-- it is a subgroup of G, and H is a normal subgroup of it: H <|= N_G(H) <= G
+normalizer gs hs = [g | g <- elts gs, all (\h -> h~^g `elem` elts hs) hs]
+
+-- stabilizer of a point
+stabilizer gs x = [g | g <- elts gs, x .^ g == x]
+
+-- pointwise stabiliser of a set
+ptStab gs xs = [g | g <- elts gs, and [x .^ g == x | x <- xs] ]
+
+-- 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]
+
+-- commutator gp of H and K
+commutatorGp hs ks = normalClosure (hsks) [h^-1 * k^-1 * h * k | h <- hs', k <- ks']
+    where hs' = reduceGens hs
+          ks' = reduceGens ks
+          hsks = reduceGens (hs' ++ ks')
+          -- no point processing more potential generators than we have to
+
+-- derived subgroup
+derivedSubgp gs = commutatorGp gs gs
+
+
+-- ACTIONS ON COSETS AND SUBGROUPS (QUOTIENT GROUPS)
+
+isSubgp hs gs = all (isMember gs) hs
+
+isNormal hs gs = isSubgp hs gs && all (isMember hs) [h~^g | h <- hs, g <- gs]
+
+
+-- 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]
+
+-- Cosets are disjoint, which leads to Lagrange's theorem
+
+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
+
+cosetAction gs hs =
+    let _H = elts hs
+        cosets_H = cosets gs _H
+    in toSn $ map (induced (**^) cosets_H) gs
+
+-- 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
+
+gs // hs = quotientGp gs hs
+
+
+-- action of group on a subgroup by conjugation
+-- (hs should be all elts, not just generators)
+hs ~~^ g = L.sort [h ~^ g | h <- hs]
+
+-- don't think that this is necessarily transitive on isomorphic subgps
+conjugateSubgps gs hs = orbit (~~^) hs gs
+-- hs ~~^^ gs = orbit (~~^) gs hs
+
+subgpAction gs hs =
+    let _H = elts hs
+        conjugates_H = conjugateSubgps gs _H
+    in toSn $ map (induced (~~^) conjugates_H) gs
+
+
+-- in cube gp, the subgps all appear to correspond to stabilisers of subsets, or of blocks
+
+ Math/Algebra/Group/SchreierSims.hs view
@@ -0,0 +1,209 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Algebra.Group.SchreierSims where
+
+import qualified Data.List as L
+import Data.Maybe (isNothing, isJust)
+import qualified Data.Map as M
+import Math.Algebra.Group.PermutationGroup hiding (elts, order, gens, isMember, isSubgp, isNormal, reduceGens, normalClosure, commutatorGp, derivedSubgp)
+
+
+-- COSET REPRESENTATIVES FOR STABILISER OF A POINT
+
+-- Given a group G = <gs>, and a point x, find coset representatives for Gx (stabiliser of x) in G
+-- In other words, for each x' in the orbit of x under G, we find a g <- G taking x to x'
+-- The code is similar to the code for calculating orbits, but modified to keep track of the group elements that we used to get there
+cosetRepsGx gs x = cosetRepsGx' gs M.empty (M.singleton x 1) where
+    cosetRepsGx' gs interior boundary
+        | M.null boundary = interior
+        | otherwise =
+            let interior' = M.union interior boundary
+                boundary' = M.fromList [(p .^ g, h*g) | g <- gs, (p,h) <- M.toList boundary] M.\\ interior'
+            in cosetRepsGx' gs interior' boundary'
+
+
+-- SCHREIER GENERATORS
+
+-- toSet xs = (map head . group . sort) xs
+
+-- Generators for Gx, the stabiliser of x, given that G is generated by gs, and rs is a set of coset representatives for Gx in G. 
+-- Schreier's Lemma states that if H < G = <S>, and R is a set of coset reps for H in G
+-- then H is generated by { rs(rs)*^-1 | r <- R, s <- S } (where * means "the coset representative of")
+-- In particular, with H = Gx, this gives us a way of finding a set of generators for Gx
+schreierGeneratorsGx (x,rs) gs = L.nub $ filter (/= 1) [schreierGenerator r g | r <- M.elems rs,  g <- gs]
+    where schreierGenerator r g = let h = r * g
+                                      h' = rs M.! (x .^ h)
+                                  in h * inverse h'
+
+
+-- SCHREIER-SIMS ALGORITHM
+
+sift bts g = sift' bts g where
+    sift' [] g = if g == 1 then Nothing else Just g
+    sift' ((b,t):bts) g = case M.lookup (b .^ g) t of
+                          Nothing -> Just g
+                          Just h -> sift' bts (g * inverse h)
+
+findBase gs = minimum $ concatMap supp gs
+{-
+-- Find base and strong generating set using Schreier-Sims algorithm
+bsgs gs | all (/= 1) gs = map fst $ ss [newLevel gs] []
+
+newLevel s = let b = findBase s
+                 t = cosetRepsGx s b
+             in ((b,t),s)
+
+ss (bad@((b,t),s):bads) goods =
+    let bts = map fst goods
+        sgs = schreierGeneratorsGx (b,t) s
+        siftees = filter isJust $ map (sift bts) sgs
+    in  if null siftees
+        then ss bads (bad:goods)
+        else let Just h = head siftees
+             in if null goods
+                then ss (newLevel [h] : bad : bads) []
+                else let ((b_,t_),s_) = head goods
+                         s' = h:s_
+                         t' = cosetRepsGx s' b_
+                     in ss (((b_,t'),s') : bad : bads) (tail goods)
+ss [] goods = goods
+-}
+-- Find base and strong generating set using Schreier-Sims algorithm
+-- This version guarantees to use bases in order
+bsgs gs = bsgs' bs gs
+    where bs = (map head . L.group . L.sort) $ concatMap supp gs
+
+-- This version lets you pass in bases in the order you want them (or [], and it will find its own)
+bsgs' bs gs = map fst $ ss bs gs
+
+-- For example, bsgs (_A 5) uses [1,2,3] as the bases, but bsgs' [] (_A 5) uses [1,3,2]
+
+newLevel (b:bs) s = (bs, newLevel' b s)
+newLevel [] s = ([], newLevel' b s) where b = findBase s
+newLevel' b s = ((b,t),s) where t = cosetRepsGx s b
+
+ss bs gs = ss' bs' [level] []
+    where (bs',level) = newLevel bs gs
+
+ss' bs (bad@((b,t),s):bads) goods =
+    let bts = map fst goods
+        sgs = schreierGeneratorsGx (b,t) s
+        siftees = filter isJust $ map (sift bts) sgs
+    in  if null siftees
+        then ss' bs bads (bad:goods)
+        else let Just h = head siftees
+             in if null goods
+                then let (bs', level) = newLevel bs [h]
+                     in ss' bs' (level : bad : bads) []
+                else let ((b_,t_),s_) = head goods
+                         s' = h:s_
+                         t' = cosetRepsGx s' b_
+                     in ss' bs (((b_,t'),s') : bad : bads) (tail goods)
+ss' _ [] goods = goods
+{-
+extendbsgs [] g = bsgs [g]
+extendbsgs (((b,t),s):bts) g = ss (((b,t),g:s):bts) []
+
+bsgs' gs = map fst $ foldl extendbsgs [] gs
+-}
+
+-- The above is written for simplicity.
+-- Its efficiency could be improved by incrementally updating the transversals,
+-- and keeping track of Schreier generators we have already tried.
+-- (Remember to add new Schreier generators every time the generating set or transversal is augmented.)
+
+
+-- USING THE SCHREIER-SIMS TRANSVERSALS
+
+isMemberBSGS bts g = isNothing $ sift bts g
+
+-- By Lagrange's thm, every g <- G can be written uniquely as g = r_m ... r_1 (Seress p56)
+-- Note that we have to reverse the list of coset representatives
+eltsBSGS bts = map (product . reverse) (cartProd ts)
+	where ts = map (M.elems . snd) bts
+
+cartProd (set:sets) = [x:xs | x <- set, xs <- cartProd sets]
+cartProd [] = [[]]
+
+orderBSGS bts = product (map (toInteger . M.size . snd) bts)
+
+
+isMember gs h = isMemberBSGS (bsgs gs) h
+
+elts gs = eltsBSGS $ bsgs gs
+
+order [] = 1
+order gs = orderBSGS $ bsgs gs
+
+isSubgp hs gs = all (isMemberBSGS gs') hs
+    where gs' = bsgs gs
+
+isNormal hs gs = hs `isSubgp` gs
+              && all (isMemberBSGS hs') [h~^g | h <- hs, g <- gs]
+    where hs' = bsgs hs
+
+index gs hs = order gs `div` order hs
+
+
+-- strong generating set
+-- sgs gs = filter (/=1) $ concatMap (M.elems . snd) $ bsgs gs
+-- sgs gs = concatMap snd $ ss [newLevel gs] []
+sgs gs = L.nub $ concatMap snd $ ss bs gs -- bs' [level] []
+    where bs = (map head . L.group . L.sort) $ concatMap supp gs
+          -- (bs',level) = newLevel bs gs
+-- !! Note, not properly tested - results not in expected order
+-- the sgs calculated during bsgs may be a smaller set (for example, the whole transversal could be powers of a single generator).
+
+
+-- given list of generators, try to find a shorter list
+reduceGens gs = fst $ reduceGensBSGS (filter (/=1) gs)
+
+reduceGensBSGS (g:gs) = reduceGens' ([g],bsgs [g]) gs where
+    reduceGens' (gs,bsgsgs) (h:hs) =
+        if isMemberBSGS bsgsgs h
+        then reduceGens' (gs,bsgsgs) hs
+        else reduceGens' (h:gs, bsgs $ h:gs) hs
+    reduceGens' (gs,bsgsgs) [] = (reverse gs,bsgsgs)
+reduceGensBSGS [] = ([],[])
+
+-- normal closure of H in G
+-- for efficiency, should be called with gs and hs already reduced sets of generators 
+normalClosure gs hs = reduceGens $ hs ++ [h ~^ g | h <- hs, g <- gs']
+    where gs' = gs ++ map inverse gs
+
+-- commutator gp of H and K
+commutatorGp hs ks = normalClosure (hsks) [h^-1 * k^-1 * h * k | h <- hs', k <- ks']
+    where hs' = reduceGens hs
+          ks' = reduceGens ks
+          hsks = reduceGens (hs' ++ ks')
+          -- no point processing more potential generators than we have to
+
+-- derived subgroup (or commutator subgroup)
+derivedSubgp gs = normalClosure gs' [g^-1 * h^-1 * g * h | g <- gs', h <- gs']
+    where gs' = reduceGens gs
+    -- == commutatorGp gs gs
+
+{-
+isPerfect gs = order gs == order (derivedSubgp gs)
+-- We compare orders rather than the generators themselves, because derivedSubgp will usually find different generators
+
+
+derivedSeries gs = derivedSeries' (gs, order gs) where
+    derivedSeries' ([],1) = [[]]
+    derivedSeries' (hs, orderhs) =
+        let hs' = derivedSubgp hs
+            orderhs' = order hs'
+        in if orderhs' == orderhs
+           then [hs]
+           else hs : derivedSeries' (hs',orderhs')
+
+lowerCentralSeries gs = lowerCentralSeries' (gs, order gs) where
+    lowerCentralSeries' ([],1) = [[]]
+    lowerCentralSeries' (hs, orderhs) =
+        let hs' = commutatorGp gs hs
+            orderhs' = order hs'
+        in if orderhs' == orderhs
+           then [hs]
+           else hs : lowerCentralSeries' (hs',orderhs')
+-}
+
+ Math/Algebra/Group/StringRewriting.hs view
@@ -0,0 +1,188 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Algebra.Group.StringRewriting where
+
+import Data.List as L
+import Data.Maybe (catMaybes)
+
+-- REWRITING
+
+rewrite rules word = rewrite' rules word where
+    rewrite' (r:rs) xs =
+        case rewrite'' r xs of
+        Nothing -> rewrite' rs xs
+        Just ys -> rewrite' rules ys
+    rewrite' [] xs = xs
+rewrite'' (l,r) xs =
+    case xs `splitSubstring` l of
+    Nothing -> Nothing
+    Just (a,b) -> Just (a++r++b)
+
+-- given a string x and a substring b, find if possible (a,c) such that xs = abc
+splitSubstring xs b = splitSubstring' [] xs where
+    splitSubstring' ls [] = Nothing
+    splitSubstring' ls (r:rs) =
+        if b `L.isPrefixOf` (r:rs)
+        then Just (reverse ls, drop (length b) (r:rs))
+        else splitSubstring' (r:ls) rs
+-- there might be a more efficient way to do this
+
+
+-- KNUTH-BENDIX
+
+-- given two strings x,y, find if possible a,b,c with x=ab y=bc
+findOverlap xs ys = findOverlap' [] xs ys where
+    findOverlap' as [] cs = Nothing -- (reverse as, [], cs)
+    findOverlap' as (b:bs) cs =
+        if (b:bs) `L.isPrefixOf` cs
+        then Just (reverse as, b:bs, drop (length (b:bs)) cs)
+        else findOverlap' (b:as) bs cs
+-- there might be a more efficient way to do this
+
+-- note that findOverlap "abab" "abab" won't find the partial overlap ("ab","ab","ab")
+
+-- Knuth-Bendix algorithm
+-- http://en.wikipedia.org/wiki/Knuth-Bendix_algorithm
+-- Given a set of rules (assumed already reduced with respect to each other)
+-- return a confluent rewrite system
+knuthBendix1 rules = knuthBendix' rules pairs where
+    pairs = [(lri,lrj) | lri <- rules, lrj <- rules, lri /= lrj]
+    knuthBendix' rules [] = rules -- should reduce in some way
+    knuthBendix' rules ( ((li,ri),(lj,rj)) : ps) =
+        case findOverlap li lj of
+        Nothing -> knuthBendix' rules ps
+        Just (a,b,c) -> case ordpair (rewrite rules (ri++c)) (rewrite rules (a++rj)) of
+                        Nothing -> knuthBendix' rules ps -- they both reduce to the same thing
+                        Just rule' -> let rules' = reduce rule' rules
+                                          ps' = ps ++ [(rule',rule) | rule <- rules'] ++ [(rule,rule') | rule <- rules']
+                                      in knuthBendix' (rule':rules') ps'
+                    	-- the new rule comes from seeing that
+                    	-- a ++ b ++ c == l1 ++ c -> r1 ++ c (by rule 1)
+                    	-- a ++ b ++ c == a ++ l2 -> a ++ r2 (by rule 2)
+    reduce rule@(l,r) rules = filter (\(l',r') -> not (L.isInfixOf l l')) rules
+        -- [rule' | rule'@(l',r') <- rules, not (l `L.isInfixOf` l')]
+-- !! Possible efficiency improvement is to somehow prune the pairs queue, as or after we prune the rules
+
+ordpair x y =
+    case shortlex x y of
+    LT -> Just (y,x)
+    EQ -> Nothing
+    GT -> Just (x,y)
+
+shortlex x y = compare (length x, x) (length y, y)
+
+-- for groups, where "letters" will take the form Either a a, we will want a different order, because we will want x^-1 -> x^3 to be the right way round
+
+
+-- An optimisation - keep the rules ordered smallest first, and process the pairs smallest first
+-- Appears to be significantly faster on average
+knuthBendix2 rules = map snd $ knuthBendix' rules' pairs where
+    rules' = L.sort $ map sizedRule rules
+    pairs = L.sort [sizedPair sri srj | sri <- rules', srj <- rules', sri /= srj]
+    knuthBendix' rules [] = rules
+    knuthBendix' rules ( (s,(li,ri),(lj,rj)) : ps) =
+        case findOverlap li lj of
+        Nothing -> knuthBendix' rules ps
+        Just (a,b,c) -> case ordpair (rewrite (map snd rules) (ri++c)) (rewrite (map snd rules) (a++rj)) of
+                        Nothing -> knuthBendix' rules ps -- they both reduce to the same thing
+                        Just rule' -> let rules' = reduce (snd rule') rules
+                                          -- ps' = L.sort $ ps ++ [sizedPair rule' rule | rule <- rules'] ++ [sizedPair rule rule' | rule <- rules']
+                                          ps' = merge ps $ merge [sizedPair rule' rule | rule <- rules'] [sizedPair rule rule' | rule <- rules']
+                                     in knuthBendix' (L.insert rule' rules') ps'
+    reduce rule@(l,r) rules = filter (\(s',(l',r')) -> not (L.isInfixOf l l')) rules
+    -- reduce rule@(l,r) rules = [rule' | rule'@(s',(l',r')) <- rules, not (l `L.isInfixOf` l')]
+    ordpair x y =
+        let lx = length x; ly = length y in
+            case compare (lx,x) (ly,y) of
+            LT -> Just (ly,(y,x)); EQ -> Nothing; GT -> Just (lx,(x,y))
+    sizedRule (rule@(l,r)) = (length l, rule)
+    sizedPair (s1,r1) (s2,r2) = (s1+s2,r1,r2)
+
+-- merge two ordered lists
+merge (x:xs) (y:ys) =
+    case compare x y of
+    LT -> x : merge xs (y:ys)
+    GT -> y : merge (x:xs) ys
+    EQ -> error "" -- shouldn't happen in our case
+merge xs ys = xs++ys
+
+-- Another optimisation - at the stage where we remove some rules, we remove corresponding pairs too
+-- Seems to perform about 25% faster on large problems (eg Coxeter groups A4-12, B4-12)
+knuthBendix3 rules = knuthBendix' rules' pairs (length rules' + 1) where
+    rules' = L.sort $ zipWith (\i (l,r) -> (length l,i,(l,r)) ) [1..] rules
+    pairs = L.sort [sizedPair ri rj | ri <- rules', rj <- rules', ri /= rj]
+    knuthBendix' rules [] k = map (\(s,i,r) -> r) rules
+    knuthBendix' rules ( (s,(i,j),((li,ri),(lj,rj))) : ps) k =
+        case findOverlap li lj of
+        Nothing -> knuthBendix' rules ps k
+        Just (a,b,c) -> case ordpair k (rewrite (map third rules) (ri++c)) (rewrite (map third rules) (a++rj)) of
+                        Nothing -> knuthBendix' rules ps k -- they both reduce to the same thing
+                        Just rule'@(_,_,(l,r)) ->
+                            let (outrules,inrules) = L.partition (\(s',i',(l',r')) -> L.isInfixOf l l') rules
+                                removedIndices = map second outrules
+                                ps' = [p | p@(s,(i,j),(ri,rj)) <- ps, i `notElem` removedIndices, j `notElem` removedIndices]
+                                ps'' = merge ps' $ merge [sizedPair rule' rule | rule <- inrules] [sizedPair rule rule' | rule <- inrules]
+                            in knuthBendix' (L.insert rule' inrules) ps'' (k+1)
+    ordpair k x y =
+        let lx = length x; ly = length y in
+            case compare (lx,x) (ly,y) of
+            LT -> Just (ly,k,(y,x)); EQ -> Nothing; GT -> Just (lx,k,(x,y))
+    second (s,i,r) = i
+    third (s,i,r) = r
+    sizedPair (si,i,ri) (sj,j,rj) = (si+sj,(i,j),(ri,rj))
+
+
+-- Version of knuthBendix that makes sure the initial rules are reduced with respect to each other
+knuthBendix rules = knuthBendix3 (reduce [] rules) where
+    reduce ls (r:rs) = reduce (r: reduce' r ls) (reduce' r rs)
+    reduce ls [] = ls
+    reduce' r rules = catMaybes [ordpair (rewrite [r] lhs) (rewrite [r] rhs) | (lhs,rhs) <- rules]
+
+
+-- given generators and rules/relations, list all normal forms
+-- the rules are assumed to be a confluent rewrite system
+nfs (gs,rs) = nfs' [[]] where
+    nfs' [] = [] -- we have run out of words - this is a finite semigroup
+    nfs' ws = let ws' = [g:w | g <- gs, w <- ws, not (any (`L.isPrefixOf` (g:w)) (map fst rs))]
+              in ws ++ nfs' ws'
+
+elts (gs,rs) = nfs (gs, knuthBendix rs)
+
+
+-- PRESENTATIONS FOR SOME STANDARD GROUPS
+-- Would like to add a few more to this list
+
+newtype S = S Int deriving (Eq,Ord)
+
+instance Show S where
+    show (S i) = "s" ++ show i
+
+s_ i = S i
+s1 = s_ 1
+s2 = s_ 2
+s3 = s_ 3
+
+-- D L Johnson, Presentations of Groups, p62
+
+-- symmetric group, generated by adjacent transpositions
+_S n = (gs, r ++ s ++ t) where
+    gs = map s_ [1..n-1]
+    r = [([s_ i, s_ i],[]) | i <- [1..n-1]]
+    s = [(concat $ replicate 3 [s_ i, s_ (i+1)],[]) | i <- [1..n-2]]
+    t = [([s_ i, s_ j, s_ i, s_ j],[]) | i <- [1..n-1], j <- [i+2..n-1]]
+
+-- http://en.wikipedia.org/wiki/Triangle_group
+-- triangle group
+tri l m n = ("abc", [("aa",""),("bb",""),("cc",""),("ab" ^ l,""),("bc" ^ n,""),("ca" ^ m,"" )])
+    where xs ^ i = concat $ replicate i xs
+
+-- von Dyck groups
+-- The subgroup of index 2 of the triangle group of elts that preserve the orientation of the triangle
+_D l m n = ("xy", [("x" ^ l,""), ("y" ^ m,""), ("xy" ^ n,"")])
+    where xs ^ i = concat $ replicate i xs
+
+-- So 2,3,3 -> tetrahedron; 2,3,4 -> cube/octahedron; 2,3,5 -> dodecahedron/icosahedron
+-- 2,2,n, n>=2 -> n-gon bipyramid
+-- Other values correspond to Euclidean or Hyperbolic groups, which are infinite
+
+
+ Math/Algebra/LinearAlgebra.hs view
@@ -0,0 +1,192 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+module Math.Algebra.LinearAlgebra where
+
+import qualified Data.List as L
+import Math.Algebra.Field.Base -- not actually used in this module
+
+
+infixr 8 *>, *>>
+infixr 7 <<*>
+infixl 7 <.>, <*>, <<*>>, <*>>
+infixl 6 <+>, <->, <<+>>, <<->>
+
+-- The mnemonic for these operations is that the number of angle brackets on each side indicates the dimension of the argument on that side
+
+u <+> v = zipWith (+) u v
+
+u <-> v = zipWith (-) u v
+
+-- scalar multiplication
+k *> v = map (k*) v
+
+k *>> m = (map . map) (k*) m
+
+-- dot product of vectors (also called inner or scalar product)
+u <.> v = sum (zipWith (*) u v)
+
+-- tensor product of vectors (also called outer or matrix product)
+u <*> v = [ [a*b | b <- v] | a <- u]
+
+
+-- matrix operations
+
+a <<+>> b = (zipWith . zipWith) (+) a b
+
+a <<->> b = (zipWith . zipWith) (-) a b
+
+a <<*>> b = [ [u <.> v | v <- L.transpose b] | u <- a]
+ 
+-- action on the left
+m <<*> v = map (<.> v) m
+
+-- action on the right
+v <*>> m = map (v <.>) (L.transpose m)
+
+
+fMatrix n f = [[f i j | j <- [1..n]] | i <- [1..n]] 
+
+-- version with indices from zero
+fMatrix' n f = [[f i j | j <- [0..n-1]] | i <- [0..n-1]] 
+
+
+-- idMx n = fMatrix n (\i j -> if i == j then 1 else 0)
+
+idMx n = idMxs !! n where
+    idMxs = map snd $ iterate next (0,[])
+    next (j,m) = (j+1, (1 : replicate j 0) : map (0:) m)
+
+jMx n = replicate n (replicate n 1)
+
+zMx n = replicate n (replicate n 0)
+
+{-
+-- VECTORS
+
+data Vector d k = V [k] deriving (Eq,Ord,Show) 
+
+instance (IntegerAsType d, Num k) => Num (Vector d k) where
+    V a + V b = V $ a <+> b
+    V a - V b = V $ a <-> b
+    negate (V a) = V $ map negate a
+    fromInteger 0 = V $ replicate d' 0 where d' = fromInteger $ value (undefined :: d)
+
+V v <>> M m = V $ v <*>> m
+
+M m <<> V v = V $ m <<*> v
+
+k |> V v = V $ k *> v
+-}
+
+-- MATRICES
+
+{-
+-- Square matrices of dimension d over field k
+data Matrix d k = M [[k]] deriving (Eq,Ord,Show)
+
+instance (IntegerAsType d, Num k) => Num (Matrix d k) where
+    M a + M b = M $ a <<+>> b
+    M a - M b = M $ a <<->> b
+    negate (M a) = M $ (map . map) negate a
+    M a * M b = M $ a <<*>> b
+    fromInteger 0 = M $ zMx d' where d' = fromInteger $ value (undefined :: d)
+    fromInteger 1 = M $ idMx d' where d' = fromInteger $ value (undefined :: d)
+
+instance (IntegerAsType d, Fractional a) => Fractional (Matrix d a) where
+	recip (M a) = case inverse a of
+		Nothing -> error "Matrix.recip: matrix is singular"
+		Just a' -> M a'
+-}
+
+inverse m =
+    let d = length m -- the dimension
+        i = idMx d
+        m' = zipWith (++) m i
+        i1 = inverse1 m'
+        i2 = inverse2 i1
+    in if length i1 == d
+       then Just i2
+       else Nothing
+
+-- given (M|I), use row operations to get to (U|A), where U is upper triangular with 1s on diagonal
+inverse1 [] = []
+inverse1 ((x:xs):rs) =
+    if x /= 0
+    then let r' = (1/x) *> xs
+         in (1:r') : inverse1 [ys <-> y *> r' | (y:ys) <- rs]
+    else case filter (\r' -> head r' /= 0) rs of
+         [] -> [] -- early termination, which will be detected in calling function
+         r:_ -> inverse1 (((x:xs) <+> r) : rs)
+-- This is basically row echelon form
+
+-- given (U|A), use row operations to get to M^-1
+inverse2 [] = []
+inverse2 ((1:r):rs) = inverse2' r rs : inverse2 rs where
+    inverse2' xs [] = xs
+    inverse2' (x:xs) ((1:r):rs) = inverse2' (xs <-> x *> r) rs
+-- This is basically reduced row echelon form
+
+xs ! i = xs !! (i-1) -- ie, a 1-based list lookup instead of 0-based
+
+rowEchelonForm [] = []
+rowEchelonForm ((x:xs):rs) =
+    if x /= 0
+    then let r' = (1/x) *> xs
+         in (1:r') : map (0:) (rowEchelonForm [ys <-> y *> r' | (y:ys) <- rs])
+    else case filter (\r' -> head r' /= 0) rs of
+         [] -> map (0:) (rowEchelonForm $ xs : map tail rs)
+         r:_ -> rowEchelonForm (((x:xs) <+> r) : rs)
+rowEchelonForm zs@([]:_) = zs
+
+reducedRowEchelonForm m = reverse $ reduce $ reverse $ rowEchelonForm m where
+    reduce (r:rs) = let r':rs' = reduceStep (r:rs) in r' : reduce rs' -- is this scanl or similar?
+    reduce [] = []
+    reduceStep ((1:xs):rs) = (1:xs) : [ 0: (ys <-> y *> xs) | y:ys <- rs]
+    reduceStep rs@((0:_):_) = zipWith (:) (map head rs) (reduceStep $ map tail rs)
+    reduceStep rs = rs
+
+-- kernel of a matrix
+-- returns basis for vectors v s.t m <<*> v == 0
+kernel m = kernelRRE $ reducedRowEchelonForm m
+
+kernelRRE m =
+    let nc = length $ head m -- the number of columns
+        is = findLeadingCols 1 (L.transpose m) -- these are the indices of the columns which have a leading 1
+        js = [1..nc] L.\\ is
+        freeCols = let m' = take (length is) m -- discard zero rows
+                   in zip is $ L.transpose [map (negate . (!j)) m' | j <- js]
+        boundCols = zip js (idMx $ length js)
+    in L.transpose $ map snd $ L.sort $ freeCols ++ boundCols
+    where
+    findLeadingCols i (c@(1:_):cs) = i : findLeadingCols (i+1) (map tail cs)
+    findLeadingCols i (c@(0:_):cs) = findLeadingCols (i+1) cs
+    findLeadingCols _ _ = []
+
+m ^- n = recip m ^ n
+
+-- t (M m) = M (L.transpose m)
+
+det [[x]] = x
+det ((x:xs):rs) =
+    if x /= 0
+    then let r' = (1/x) *> xs
+         in x * det [ys <-> y *> r' | (y:ys) <- rs]
+    else case filter (\r' -> head r' /= 0) rs of
+         [] -> 0
+         r:_ -> det (((x:xs) <+> r) : rs)
+
+
+{-
+class IntegerAsType a where
+    value :: a -> Integer
+
+data Z
+instance IntegerAsType Z where
+    value _ = 0
+
+data S a
+instance IntegerAsType a => IntegerAsType (S a) where
+    value _ = value (undefined :: a) + 1
+-}
+ Math/Algebra/NonCommutative/GSBasis.hs view
@@ -0,0 +1,66 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Algebra.NonCommutative.GSBasis where
+
+import Data.List as L
+
+import Math.Algebra.NonCommutative.NCPoly
+
+
+-- given two monomials f g, find if possible a,b,c with f=ab g=bc
+findOverlap (M xs) (M ys) = findOverlap' [] xs ys where
+    findOverlap' as [] cs = Nothing -- (reverse as, [], cs)
+    findOverlap' as (b:bs) cs =
+        if (b:bs) `L.isPrefixOf` cs
+        then Just (M $ reverse as, M $ b:bs, M $ drop (length (b:bs)) cs)
+        else findOverlap' (b:as) bs cs
+
+-- given two monomials f g, find if possible l,r with g = lfr
+-- findInclusion (M xs) (M ys) = findInclusion' 
+
+sPoly f@(NP ((xs,c):_)) g@(NP ((ys,d):_)) =
+    case findOverlap xs ys of
+    Just (l,m,r) -> f * NP [(r,d)] - NP [(l,c)] * g
+    Nothing -> 0
+sPoly _ _ = 0 -- !! shouldn't reach this
+-- The point about the s-poly is that it cancels out the leading terms of the two polys, exposing their second terms
+
+
+gb1 fs = gb' fs [sPoly fi fj | fi <- fs, fj <- fs, fi /= fj] where -- unlike the commutative case, we take sPolys both ways round
+    gb' gs (h:hs) = let h' = h %% gs in
+                    if h' == 0 then gb' gs hs else gb' (h':gs) (hs ++ [sPoly h' g | g <- gs] ++ [sPoly g h' | g <- gs])
+    gb' gs [] = gs
+
+reduce gs = reduce' [] gs where
+    reduce' gs' (g:gs) | g' == 0   = reduce' gs' gs
+                       | otherwise = reduce' (g':gs') gs
+                       where g' = g %% (gs'++gs)
+    reduce' gs' [] = reverse $ sort $ gs'
+
+gb fs = map toMonic $ reduce $ gb1 fs
+
+gb' fs = reduce $ gb1 fs
+
+
+gb2 fs = gb' fs [(fi,fj) | fi <- fs, fj <- fs, fi /= fj] where -- unlike the commutative case, we take sPolys both ways round
+    gb' gs ((fi,fj):pairs) =
+        let h = sPoly fi fj %% gs in
+        if h == 0 then gb' gs pairs else gb' (h:gs) (pairs ++ [(h,g) | g <- gs] ++ [(g,h) | g <- gs])
+    gb' gs [] = gs
+
+gb2' fs = gb' fs [(fi,fj) | fi <- fs, fj <- fs, fi /= fj] where -- unlike the commutative case, we take sPolys both ways round
+    gb' gs ((fi,fj):pairs) =
+        let h = sPoly fi fj %% gs in
+        if h == 0 then gb' gs pairs else (fi,fj,sPoly fi fj,h) : gb' (h:gs) (pairs ++ [(h,g) | g <- gs] ++ [(g,h) | g <- gs])
+    gb' gs [] = [] -- gs
+
+
+-- Monomial basis for the quotient algebra, where gs are the generators, rs the relations
+mbasisQA gs rs = mbasisQA' [1] where
+    mbasisQA' [] = [] -- the quotient ring has a finite monomial basis
+    mbasisQA' ms = let ms' = [g*m | g <- gs, m <- ms, g*m %% rs == g*m] -- ie, not reducible
+                   in ms ++ mbasisQA' ms'
+{-
+isGB fs = all (\h -> h %% fs == 0) (pairWith sPoly fs)
+-}
+
+ Math/Algebra/NonCommutative/NCPoly.hs view
@@ -0,0 +1,207 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Algebra.NonCommutative.NCPoly where
+
+import Data.List as L
+import Math.Algebra.Field.Base
+
+
+-- (NON-COMMUTATIVE) MONOMIALS
+
+newtype Monomial v = M [v] deriving (Eq)
+
+instance Ord v => Ord (Monomial v) where
+    compare (M xs) (M ys) = compare (length xs,xs) (length ys,ys)
+-- Glex ordering
+
+instance Show v => Show (Monomial v) where
+    show (M xs) | null xs = "1"
+                | otherwise = concatMap show xs
+-- !! we can do better - we should show "xx" as "x^2"
+
+instance (Eq v, Show v) => Num (Monomial v) where
+    M xs * M ys = M (xs ++ ys)
+    fromInteger 1 = M []
+
+-- try to find l, r such that a = lbr
+divM (M a) (M b) = divM' [] a where
+    divM' ls (r:rs) =
+        if b `L.isPrefixOf` (r:rs)
+        then Just (M $ reverse ls, M $ drop (length b) (r:rs))
+        else divM' (r:ls) rs
+    divM' _ [] = Nothing
+
+
+-- (NON-COMMUTATIVE) POLYNOMIALS
+
+newtype NPoly r v = NP [(Monomial v,r)] deriving (Eq)
+
+instance (Ord r, Ord v) => Ord (NPoly r v) where
+    compare (NP ts) (NP us) = compare ts us
+
+instance (Show r, Eq v, Show v) => Show (NPoly r v) where
+    show (NP []) = "0"
+    show (NP ts) =
+        let (c:cs) = concatMap showTerm ts
+        in if c == '+' then cs else c:cs
+        where showTerm (m,a) =
+                  case show a of
+                  "1" -> "+" ++ show m
+                  "-1" -> "-" ++ show m
+                  -- cs@(x:_) -> (if x == '-' then cs else '+':cs) ++ (if m == 1 then "" else show m)
+                  cs -> showCoeff cs ++ (if m == 1 then "" else show m)
+              showCoeff (c:cs) = if any (`elem` ['+','-']) cs
+                                 then "+(" ++ c:cs ++ ")"
+                                 else if c == '-' then c:cs else '+':c:cs
+
+
+instance (Ord v, Show v, Num r) => Num (NPoly r v) where
+    NP ts + NP us = NP (mergeTerms ts us)
+    negate (NP ts) = NP $ map (\(m,c) -> (m,-c)) ts
+    NP ts * NP us = NP $ collect $ L.sortBy cmpTerm $ [(g*h,c*d) | (g,c) <- ts, (h,d) <- us]
+    fromInteger 0 = NP []
+    fromInteger n = NP [(fromInteger 1, fromInteger n)]
+
+cmpTerm (a,c) (b,d) = case compare a b of EQ -> EQ; GT -> LT; LT -> GT -- in mpolys we put "larger" terms first
+
+-- inputs in descending order
+mergeTerms (t@(g,c):ts) (u@(h,d):us) =
+    case cmpTerm t u of
+    LT -> t : mergeTerms ts (u:us)
+    GT -> u : mergeTerms (t:ts) us
+    EQ -> if e == 0 then mergeTerms ts us else (g,e) : mergeTerms ts us
+    where e = c + d
+mergeTerms ts us = ts ++ us -- one of them is null
+
+collect (t1@(g,c):t2@(h,d):ts)
+    | g == h = collect $ (g,c+d):ts
+    | c == 0  = collect $ t2:ts
+    | otherwise = t1 : collect (t2:ts)
+collect ts = ts
+
+
+-- Fractional instance so that we can enter fractional coefficients
+-- Only lets us divide by field elements (with unit monomial), not any other polynomials
+instance (Ord v, Show v, Fractional r) => Fractional (NPoly r v) where
+    recip (NP [(1,c)]) = NP [(1, recip c)]
+    recip _ = error "NPoly.recip: only supported for (non-zero) constants"
+
+
+-- SOME VARIABLES (INDETERMINATES)
+-- The idea is that you define your own type of indeterminates as required, along the same lines as this
+
+data Var = X | Y | Z deriving (Eq,Ord)
+
+instance Show Var where
+    show X = "x"
+    show Y = "y"
+    show Z = "z"
+
+var v = NP [(M [v], 1)]
+
+x = var X :: NPoly Q Var
+y = var Y :: NPoly Q Var
+z = var Z :: NPoly Q Var
+
+
+-- DIVISION ALGORITHM
+
+lm (NP ((m,c):ts)) = m
+lc (NP ((m,c):ts)) = c
+lt (NP (t:ts)) = NP [t]
+
+-- given f, gs, find ls, rs, f' such that f = sum (zipWith3 (*) ls gs rs) + f', with f' not divisible by any g
+quotRemNP f gs | all (/=0) gs = quotRemNP' f (replicate n (0,0), 0)
+               | otherwise = error "quotRemNP: division by zero"
+    where
+    n = length gs
+    quotRemNP' 0 (lrs,f') = (lrs,f')
+    quotRemNP' h (lrs,f') = divisionStep h (gs,[],lrs,f')
+    divisionStep h (g:gs, lrs', (l,r):lrs, f') =
+        case lm h `divM` lm g of
+        Just (l',r') -> let l'' = NP [(l',lc h / lc g)]
+                            r'' = NP [(r',1)]
+                            h' = h - l'' * g * r''
+                        in quotRemNP' h' (reverse lrs' ++ (l+l'',r+r''):lrs, f')
+        Nothing -> divisionStep h (gs,(l,r):lrs',lrs,f')
+    divisionStep h ([],lrs',[],f') =
+        let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
+        in quotRemNP' (h-lth) (reverse lrs', f'+lth)
+
+-- It is only marginally (5-10%) more space/time efficient not to track the (lazily unevaluated) factors
+remNP f gs | all (/=0) gs = remNP' f 0
+-- let result = remNP' f 0 in if result == remNP2 f gs then result else error ("remNP2 " ++ show f ++ " " ++ show gs)
+           | otherwise = error "remNP: division by zero"
+    where
+    n = length gs
+    remNP' 0 f' = f'
+    remNP' h f' = divisionStep h gs f'
+    divisionStep h (g:gs) f' =
+        case lm h `divM` lm g of
+        Just (l',r') -> let l'' = NP [(l',lc h / lc g)]
+                            r'' = NP [(r',1)]
+                            h' = h - l'' * g * r''
+                        in remNP' h' f'
+        Nothing -> divisionStep h gs f'
+    divisionStep h [] f' =
+        let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
+        in remNP' (h-lth) (f'+lth)
+
+infixl 7 %%
+-- f %% gs = r where (_,r) = quotRemNP f gs
+f %% gs = remNP f gs
+
+-- !! Not sure if the following is valid
+-- The idea is to avoid dividing by lc g, because sometimes our coefficient ring is not a field
+-- Passes all the knot theory tests
+-- However, it may be that if we ever get a non-invertible element at the front, we are in trouble anyway
+remNP2 f gs | all (/=0) gs = remNP' f 0
+           | otherwise = error "remNP: division by zero"
+    where
+    n = length gs
+    remNP' 0 f' = f'
+    remNP' h f' = divisionStep h gs f'
+    divisionStep h (g:gs) f' =
+        case lm h `divM` lm g of
+        Just (l',r') -> let l'' = NP [(l',1)] -- NP [(l',lc h / lc g)]
+                            r'' = NP [(r',1)]
+                            lcg = inject (lc g)
+                            lch = inject (lc h)
+                            -- h' = h - l'' * g * r''
+                            h' = lcg * h - lch * l'' * g * r''
+                        in remNP' h' (lcg * f') -- must multiply f' by lcg too (otherwise get incorrect results, eg tlBasis 4)
+        Nothing -> divisionStep h gs f'
+    divisionStep h [] f' =
+        let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
+        in remNP' (h-lth) (f'+lth)
+
+
+-- OTHER STUFF
+
+toMonic 0 = 0
+toMonic (NP ts@((_,c):_))
+    | c == 1 = NP ts
+    | otherwise = NP $ map (\(m,d)->(m,d/c)) ts
+
+-- injection of field elements into polynomial ring
+inject 0 = NP []
+inject c = NP [(fromInteger 1, c)]
+
+-- substitute terms for variables in an NPoly
+-- eg subst [(x,a),(y,a+b),(z,c^2)] (x*y+z) -> a*(a+b)+c^2
+subst vts (NP us) = sum [inject c * substM m | (m,c) <- us] where
+    substM (M xs) = product [substV x | x <- xs]
+    substV v =
+        let v' = NP [(M [v], 1)] in
+        case L.lookup v' vts of
+        Just t -> t
+        Nothing -> error ("subst: no substitute supplied for " ++ show v')
+
+
+-- INVERTIBLE
+-- To support algebras which have invertible elements
+
+class Invertible a where
+    inv :: a -> a
+
+x ^- k = inv x ^ k
+ Math/Algebra/NonCommutative/TensorAlgebra.hs view
@@ -0,0 +1,132 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Algebra.NonCommutative.TensorAlgebra where
+
+import Math.Algebra.Field.Base
+import Math.Algebra.NonCommutative.NCPoly hiding (X)
+import Math.Algebra.NonCommutative.GSBasis
+
+
+-- TENSOR ALGEBRA
+-- Tensor product satisfies the universal property that any multilinear map from the cartesian product can be factored through the tensor product
+
+-- The tensor algebra is the free algebra on the basis elts of the vector space
+
+data Basis = E Int deriving (Eq,Ord)
+
+instance Show Basis where
+    show (E i) = 'e': show i
+
+e_ i = NP [(M [E i], 1)] :: NPoly Q Basis
+
+e1 = e_ 1
+e2 = e_ 2
+e3 = e_ 3
+e4 = e_ 4
+
+-- given an elt of the tensor algebra, return the dimension of the vector space it's defined over
+dim (NP ts) = maximum $ 0 : [i | (M bs,c) <- ts, E i <- bs]
+
+-- Monomial basis for tensor algebra over k^n - infinite
+tensorBasis n = mbasisQA [e_ i | i <- [1..n]] []
+
+
+-- EXTERIOR ALGEBRA
+-- Exterior product satisfies the universal property that any alternating multilinear map from the cartesian product can be factored through the exterior product
+
+-- Exterior algebra over k^n is tensor algebra over k^n quotiented by these relations
+extRelations n = [e_ i * e_ i | i <- [1..n] ] ++
+                 [e_ i * e_ j + e_ j * e_ i | i <- [1..n], j <- [i+1..n] ]
+
+extnf t = t %% (extRelations $ dim t)
+
+-- Monomial basis for exterior algebra over k^n - finite
+exteriorBasis n = mbasisQA [e_ i | i <- [1..n]] $ extRelations n
+
+
+-- SYMMETRIC ALGEBRA
+-- Symmetric product satisfies the universal property that any symmetric multilinear map from the cartesian product can be factored through the symmetric product
+
+-- Symmetric algebra over k^n is tensor algebra over k^n quotiented by these relations
+symRelations n = [e_ i * e_ j - e_ j * e_ i | i <- [1..n], j <- [i+1..n] ]
+
+symnf t = t %% (symRelations $ dim t)
+
+-- Monomial basis for symmetric algebra over k^n - infinite
+symmetricBasis n = mbasisQA [e_ i | i <- [1..n]] $ symRelations n
+
+
+-- WEYL ALGEBRAS
+-- http://en.wikipedia.org/wiki/Weyl_algebra
+-- Coutinho, A Primer of Algebraic D-modules, ch1
+
+-- Given a symplectic form w, represented by
+-- [0  I]
+-- [-I 0]
+-- on R^2n
+-- The Weyl algebra is the tensor algebra quotiented by < u*v-v*u-w(u,v) >
+-- It has a natural interpretation as an operator algebra in which
+-- e_1 .. e_i .. e_n correspond to x_i (the "multiply by x_i" operator),
+-- e_n+1 .. e_n+i .. e_2*n correspond to d_x_i (the "differentiate wrt x_i" operator)
+
+-- Weyl algebra W(V) is a "quantization" of the Symmetric algebra Sym(V)
+
+weylRelations n = [e_ j * e_ i - e_ i * e_ j | i <- [1..2*n], j <- [i+1..2*n], j /= i+n ] ++
+                  [e_ (i+n) * e_ i - e_ i * e_ (i+n) - 1 | i <- [1..n] ]
+
+weylnf n t = t %% (weylRelations n)
+
+weylBasis n = mbasisQA [e_ i | i <- [1..2*n]] $ weylRelations n
+
+
+-- Explicit construction of Weyl algebra in terms of d_x_i and x_i operators
+
+data WeylGens = X Int | D Int deriving (Eq,Ord)
+
+instance Show WeylGens where
+    show (D i) = 'd': show i
+    show (X i) = 'x': show i
+
+d_ i = NP [(M [D i], 1)] :: NPoly Q WeylGens
+x_ i = NP [(M [X i], 1)] :: NPoly Q WeylGens
+
+d1 = d_ 1
+d2 = d_ 2
+d3 = d_ 3
+
+x1 = x_ 1
+x2 = x_ 2
+x3 = x_ 3
+
+comm p q = p*q - q*p
+
+delta i j = if i == j then 1 else 0
+
+weylRelations' n = [comm (x_ i) (x_ j) | i <- [1..n], j <- [i+1..n] ] ++
+                   [comm (d_ i) (d_ j) | i <- [1..n], j <- [i+1..n] ] ++
+                   [comm (d_ i) (x_ j) - delta i j | i <- [1..n], j <- [1..n] ]
+
+weylnf' f@(NP ts) = f %% weylRelations' n where
+    n = maximum $ 0 : [i | (M bs,c) <- ts, X i <- bs] ++ [i | (M bs,c) <- ts, D i <- bs]
+
+weylBasis' n = mbasisQA (map x_ [1..n] ++ map d_ [1..n]) (weylRelations' n)
+
+{-
+-- HEISENBERG ALGEBRA
+
+data Heisenberg = D | U deriving (Eq,Ord)
+
+instance Show Heisenberg where
+    show D = "d"
+    show U = "u"
+
+d = NP [(M [D], 1)] :: NPoly Q Heisenberg
+u = NP [(M [U], 1)] :: NPoly Q Heisenberg
+
+
+heisenberg = [u*d-d*u-1]
+
+
+-- Monomial basis for Heisenberg algebra - infinite
+hBasis = mbasisQA [d,u] (gb heisenberg)
+-}
+ Math/Combinatorics/Design.hs view
@@ -0,0 +1,437 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Combinatorics.Design where
+
+import Data.Maybe (fromJust, isJust)
+import qualified Data.List as L
+import qualified Data.Map as M
+import qualified Data.Set as S
+
+import Math.Common.ListSet (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, combinationsOf)
+import Math.Combinatorics.GraphAuts (refine, isSingleton, graphAuts, removeGens)
+import Math.Combinatorics.FiniteGeometry
+
+-- Cameron & van Lint, Designs, Graphs, Codes and their Links
+
+
+{-
+set xs = map head $ group $ sort xs
+
+-- subsets of size k (returned in ascending order)
+combinationsOf 0 _ = [[]]
+combinationsOf _ [] = []
+combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs
+-}
+isSubset xs ys = all (`elem` ys) xs
+
+
+-- DESIGNS
+
+data Design a = D [a] [[a]] deriving (Eq,Ord,Show)
+
+design (xs,bs) | isValid d = d where d = D xs bs
+
+toDesign (xs,bs) = D xs' bs' where
+    xs' = L.sort xs
+    bs' = L.sort $ map L.sort bs -- in fact don't require that the blocks are in order
+
+isValid (D xs bs) = (xs == L.sort xs || error "design: points are not in order")
+                 && (all (\b -> b == L.sort b) bs || error "design: blocks do not have points in order")
+-- could also check that each block is a subset of xs, etc
+
+points (D xs bs) = xs
+
+blocks (D xs bs) = bs
+
+
+-- FINDING DESIGN PARAMETERS
+
+noRepeatedBlocks (D xs bs) = all ( (==1) . length ) $ L.group $ L.sort bs
+
+
+-- Note that the design parameters functions don't check no repeated blocks, so they're also valid for t-structures
+
+-- given t and a t-(v,k,lambda) design, return (v,k,lambda)
+tDesignParams t d =
+    case findvk d of
+    Nothing -> Nothing
+    Just (v,k) ->
+        case findlambda t d of
+        Nothing -> Nothing
+        Just lambda -> Just (v,k,lambda)
+
+findvk (D xs bs) =
+    let k:ls = map length bs
+    in if all (==k) ls then Just (v,k) else Nothing
+    where v = length xs
+
+findlambda t (D xs bs) =
+    let lambda:ls = [length [b | b <- bs, ts `isSubset` b] | ts <- combinationsOf t xs]
+    in if all (==lambda) ls then Just lambda else Nothing
+
+-- given (xs,bs), return design parameters t-(v,k,lambda) with t maximal
+designParams d =
+    case findvk d of
+    Nothing -> Nothing
+    Just (v,k) ->
+        case reverse (takeWhile (isJust . snd) [(t, findlambda t d) | t <- [0..k] ]) of
+        [] -> Nothing
+        (t,Just lambda):_ -> Just (t,(v,k,lambda))
+
+
+isStructure t d = isJust $ tDesignParams t d
+
+isDesign t d = noRepeatedBlocks d && isStructure t d
+
+is2Design d = isDesign 2 d
+
+-- square 2-design (more often called "symmetric" in the literature)
+isSquare d@(D xs bs) = is2Design d && length xs == length bs
+
+
+-- incidence matrix of a design
+-- (rows and columns indexed by blocks and points respectively)
+-- (this follows Cameron & van Lint, though elsewhere in the literature it is sometimes the other way round)
+incidenceMatrix (D xs bs) = [ [if x `elem` b then 1 else 0 | x <- xs] | b <- bs]
+
+
+-- SOME FAMILIES OF DESIGNS
+
+-- the following is trivially a k-(v,k,lambda) design
+subsetDesign v k = design (xs,bs) where
+    xs = [1..v]
+    bs = combinationsOf k xs
+
+-- Cameron & van Lint, p30
+-- the pair design on n points is the complete graph on n points considered as a 2-(n,2,1) design
+pairDesign n = D vs es where
+    graph = G.k n
+    vs = vertices graph
+    es = edges graph
+
+-- the affine plane AG(2,Fq) - a 2-(q^2,q,1) design
+ag2 fq = design (points, lines) where
+    points = ptsAG 2 fq
+    lines = map line $ tail $ ptsPG 2 fq
+    line [a,b,c] = [ [x,y] | [x,y] <- points, a*x+b*y+c==0 ]
+
+-- the projective plane PG(2,Fq) - a square 2-(q^2+q+1,q+1,1) design
+pg2 fq = design (points, lines) where
+    points = ptsPG 2 fq
+    lines = map line points
+    line u = [v | v <- points, u <.> v == 0]
+    u <.> v = sum (zipWith (*) u v)
+-- Remember that the points and lines of PG(2,Fp) are really the lines and planes of AG(3,Fp).
+-- A line in AG(3,Fp) defines a plane orthogonal to it.
+
+
+-- The points and i-flats of PG(n,fq), 1<=i<=n-1, form a 2-design
+-- For i==1, this is a 2-((q^(n+1)-1)/(q-1),q+1,1) design
+-- For i==n-1, this is a 2-((q^(n+1)-1)/(q-1),(q^n-1)/(q-1),(q^(n-1)-1)/(q-1)) design
+-- Cameron & van Lint, p8
+flatsDesignPG n fq k = design (points, blocks) where
+    points = ptsPG n fq
+    blocks = map closurePG $ flatsPG n fq k -- the closurePG replaces the generators of the flat by the list of points of the flat
+
+-- The projective point-hyperplane design is also denoted PG(n,q)
+pg n fq = flatsDesignPG n fq (n-1)
+
+-- (Cameron & van Lint don't actually state that this is a design except when k == n-1)
+flatsDesignAG n fq k = design (points, blocks) where
+    points = ptsAG n fq
+    blocks = map closureAG $ flatsAG n fq k -- the closureAG replaces the generators of the flat by the list of points of the flat
+
+-- The affine point-hyperplane design is also denoted AG(n,q)
+-- It a 2-(q^n,q^(n-1),(q^(n-1)-1)/(q-1)) design
+-- Cameron & van Lint, p17
+ag n fq = flatsDesignAG n fq (n-1)
+
+
+
+-- convert a design to be defined over the set [1..n]
+to1n (D xs bs) = (D xs' bs') where
+    mapping = M.fromList $ zip xs [1..] -- the mapping from vs to [1..n]
+    xs' = M.elems mapping
+    bs' = [map (mapping M.!) b | b <- bs] -- the blocks will already be sorted correctly by construction
+
+
+-- Cameron & van Lint p10
+paleyDesign fq | length fq `mod` 4 == 3 = design (xs,bs) where
+    xs = fq
+    qs = set [x^2 | x <- xs] L.\\ [0] -- the non-zero squares in Fq
+    bs = [L.sort (map (x+) qs) | x <- xs]
+
+fanoPlane = paleyDesign f7
+-- isomorphic to PG(2,F2)
+
+
+-- NEW DESIGNS FROM OLD
+
+-- Dual of a design. Cameron & van Lint p11
+dual (D xs bs) = design (bs, map beta xs) where
+    beta x = filter (x `elem`) bs
+
+-- Derived design relative to a point. Cameron & van Lint p11
+-- Derived design of a t-(v,k,lambda) is a t-1-(v-1,k-1,lambda) design.
+derivedDesign (D xs bs) p = design (xs L.\\ [p], [b L.\\ [p] | b <- bs, p `elem` b])
+
+-- Residual design relative to a point. Cameron & van Lint p13
+-- Point-residual of a t-(v,k,lambda) is a t-1-(v-1,k,mu).
+pointResidual (D xs bs) p = design (xs L.\\ [p], [b | b <- bs, p `notElem` b])
+
+-- Complementary design relative to a point. Cameron & van Lint p13
+-- Complement of a t-(v,k,lambda) is a t-(v,v-k,mu).
+complementaryDesign (D xs bs) = design (xs, [xs L.\\ b | b <- bs])
+
+-- Residual design relative to a block. Cameron & van Lint p13
+-- This is only a design if (xs,bs) is a square design
+-- It may have repeated blocks - but if so, residuals of the complement will not
+-- Block-residual of a 2-(v,k,lambda) is a 2-(v-k,k-lambda,lambda).
+blockResidual d@(D xs bs) b | isSquare d = design (xs L.\\ b, [b' L.\\ b | b' <- bs, b' /= b])
+
+
+-- DESIGN AUTOMORPHISMS
+
+isDesignAut (D xs bs) g | supp g `isSubset` xs = all (`S.member` bs') [b -^ g | b <- bs]
+    where bs' = S.fromList bs
+
+incidenceGraph (D xs bs) = graph (vs,es) where
+    vs = L.sort $ map Left xs ++ map Right bs
+    es = L.sort [ [Left x, Right b] | x <- xs, b <- bs, x `elem` b ]
+
+-- We find design auts by finding graph auts of the incidence graph of the design
+-- In a square design, we need to watch out for graph auts which are mapping points <-> blocks
+designAuts1 d = filter (/=1) $ map points $ graphAuts $ incidenceGraph d where
+    points (P m) = P $ M.fromList [(x,y) | (Left x, Left y) <- M.toList m]
+     -- This implicitly filters out (Right x, Right y) action on blocks,
+     -- and also (Left x, Right y) auts taking points to blocks.
+     -- The filter (/=1) is to remove points <-> blocks auts
+
+-- 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
+
+
+-- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
+designAuts d@(D xs bs) = map points (designAuts' [] [vs]) where
+    n = length xs
+    g@(G vs es) = incidenceGraph d
+    points (P m) = P $ M.fromList [(k,v) | (Left k, Left v) <- M.toList m]
+    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
+betaL2_23  = p [[-1],[0],[1,2,4,8,16,9,18,13,3,6,12],[5,10,20,17,11,22,21,19,15,7,14]]                  -- t -> 2*t
+gammaL2_23 = p [[-1,0],[1,22],[2,11],[3,15],[4,17],[5,9],[6,19],[7,13],[8,20],[10,16],[12,21],[14,18]]  -- t -> -1/t
+
+l2_23 = [alphaL2_23, betaL2_23, gammaL2_23]
+
+-- Mathieu group M24
+-- Conway and Sloane p274ff
+-- This is the automorphism group of the extended binary Golay code G24
+-- or alternatively of the unique Steiner system S(5,8,24) (which consists of the weight 8 codewords of the above)
+
+deltaM24 = p [[-1],[0],[1,18,4,2,6],[3],[5,21,20,10,7],[8,16,13,9,12],[11,19,22,14,17],[15]]
+-- this is t -> t^3 / 9 (for t a quadratic residue), t -> 9 t^3 (t a non-residue)
+
+m24 = [alphaL2_23, betaL2_23, gammaL2_23, deltaM24]
+
+m24sgs = sgs m24
+-- order 244823040
+
+m23sgs = filter (\g -> (-1).^g == -1) m24sgs
+-- order 10200960
+
+m22sgs = filter (\g -> 0.^g == 0) m23sgs
+-- order 443520
+
+-- !! The above assume that the base is [-1,0,..], which isn't guaranteed
+
+
+-- Steiner system S(5,8,24)
+
+octad = [0,1,2,3,4,7,10,12]
+-- Conway&Sloane p276 - this is a weight 8 codeword from Golay code G24
+
+
+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) )
+
+s_4_7_23 = derivedDesign s_5_8_24 (-1)
+-- 253 blocks ( (23 `choose` 4) `div` (7 `choose` 4) )
+
+s_3_6_22 = derivedDesign s_4_7_23 0
+-- 77 blocks
+
+-- 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_4_7_23' = [xs | (-1):xs <- s_5_8_24]
+-- 253 blocks ( (23 `choose` 4) `div` (7 `choose` 4) )
+
+s_3_6_22' = [xs | 0:xs <- s_4_7_23]
+-- 77 blocks
+-}
+
+
+{-
+-- WITT DESIGNS
+-- S(5,8,24) AND S(5,6,12)
+
+-- Let D be a square 2-design.
+-- An n-arc is a set of n points of D, no three of which are contained in a block
+arcs n (D xs bs) = map reverse $ dfs n [] xs where
+    dfs 0 ys _ = [ys]
+    dfs i ys xs = concat [dfs (i-1) (x:ys) (dropWhile (<=x) xs) | x <- xs, isCompatible (x:ys)]
+    isCompatible ys = all ((<=2) . length) [ys `L.intersect` b | b <- bs]
+
+tangents (D xs bs) arc = [b | b <- bs, length (arc `L.intersect` b) == 1]
+
+-- !! NOT QUITE AS EXPECTED
+-- Cameron van Lint implies that ovals should have n = 1+(k-1)/lambda, whereas I'm finding that they're one bigger than that
+-- eg length $ ovals $ ag2 f3 should be 54
+-- But ag2 f3 isn't a *square* design
+ovals d =
+    let Just (_,k,lambda) = tDesignParams 2 d
+        (q,r) = (k-1) `quotRem` lambda
+        n = 2+q -- == 1+(k-1)/lambda
+    in if r == 0
+       then [arc | arc <- arcs n d, arc == L.sort (concat $ map (L.intersect arc) $ tangents d arc)] -- each point has a unique tangent
+       else []
+
+hyperovals d =
+    let Just (_,k,lambda) = tDesignParams 2 d
+        (q,r) = k `quotRem` lambda
+        n = 1+q -- == 1+k/lambda
+    in if r == 0
+       then filter (null . tangents d) $ arcs n d
+       else []
+
+-- Cameron & van Lint, p22
+-- s_5_8_24 = [length (intersect (head h) (head s)) | h <- [h1,h2,h3], s <- [s1,s2,s3]] where
+s_5_8_24 = design (points,lines) where
+    points = map Left xs ++ map Right [1,2,3]
+    lines = [map Left b ++ map Right [1,2,3] | b <- bs] ++ -- line plus three points at infinity
+            [map Left h ++ map Right [2,3] | h <- h1] ++ -- hyperoval plus two points at infinity
+            [map Left h ++ map Right [1,3] | h <- h2] ++
+            [map Left h ++ map Right [1,2] | h <- h3] ++
+            [map Left s ++ map Right [1] | s <- s1] ++ -- Baer subplanes plus one point at infinity
+            [map Left s ++ map Right [2] | s <- s2] ++
+            [map Left s ++ map Right [3] | s <- s3] ++
+            [map Left (l1 `symDiff` l2) | l1 <- bs, l2 <- bs, l1 < l2]
+    d@(D xs bs) = pg2 f4
+    hs = hyperovals d
+    [h1,h2,h3] = evenClasses hs
+    [s2,s1,s3] = oddClasses baerSubplanes -- we have to number the ss so that if h <- hi, s <- sj, then |h intersect s| is even <=> i == j
+    evenClasses (h:hs) = let (ys,ns) = partition (even . length . L.intersect h) hs in (h:ys) : evenClasses ns
+    evenClasses [] = []
+    oddClasses (h:hs) = let (ys,ns) = partition (odd . length . L.intersect h) hs in (h:ys) : oddClasses ns
+    oddClasses [] = []
+    baerSubplanes = [s | s <- baerSubplanes', and [length (L.intersect s b) `elem` [1,3] | b <- bs] ] 
+    baerSubplanes' = map reverse $ dfs 7 [] xs where
+        dfs 0 ys _ = [ys]
+        dfs i ys xs = concat [dfs (i-1) (x:ys) (dropWhile (<=x) xs) | x <- xs, isCompatible (x:ys)]
+        isCompatible ys = all ((<=3) . length) [ys `L.intersect` b | b <- bs]
+-}
+ Math/Combinatorics/FiniteGeometry.hs view
@@ -0,0 +1,117 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Combinatorics.FiniteGeometry where
+
+import Data.List as L
+import qualified Data.Set as S
+-- import qualified Data.Map as M -- not really required
+import Math.Algebra.Field.Base
+import Math.Algebra.Field.Extension hiding ( (<+>) )
+import Math.Algebra.LinearAlgebra -- hiding ( det )
+
+-- import PermutationGroup
+-- import SchreierSims as SS
+
+
+-- !! This should really live somewhere else
+-- subsets of size k
+combinationsOf 0 _ = [[]]
+combinationsOf _ [] = []
+combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs
+
+
+
+ptsAG 0 fq = [[]]
+ptsAG n fq = [x:xs | x <- fq, xs <- ptsAG (n-1) fq]
+
+ptsPG 0 _ = [[1]]
+ptsPG n fq = map (0:) (ptsPG (n-1) fq) ++ map (1:) (ptsAG n fq)
+
+
+-- "projective normal form"
+pnf (0:xs) = 0 : pnf xs
+pnf (1:xs) = 1 : xs
+pnf (x:xs) = 1 : map (* x') xs where x' = recip x
+
+ispnf (0:xs) = ispnf xs
+ispnf (1:xs) = True
+ispnf _ = False
+
+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
+          fq = eltsFq undefined
+
+
+-- 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
+    k = length ps
+    fq = eltsFq undefined
+
+
+
+-- van Lint & Wilson, p325, 332
+qtorial n q | n >= 0 = product [(q^i - 1) `div` (q-1) | i <- [1..n]]
+
+-- van Lint & Wilson, p326
+qnomial n k q = (n `qtorial` q) `div` ( (k `qtorial` q) * ((n-k) `qtorial` q) )
+
+-- Cameron, p129
+numFlatsPG n q k = qnomial (n+1) (k+1) q -- because it's the number of subspaces in AG n+1 
+
+-- Cameron, p137
+numFlatsAG n q k = q^(n-k) * qnomial n k q
+
+
+qtorials q = scanl (*) 1 [(q^i - 1) `div` (q-1) | i <- [1..]]
+
+qnomials q = iterate succ [1] where
+    succ xs = L.zipWith3 (\l r c -> l+c*r) (0:xs) (xs++[0]) (iterate (*q) 1)
+    -- succ xs = zipWith (+) (0:xs) $ zipWith (*) (xs++[0]) $ iterate (*q) 1 
+-- This amounts to saying
+-- [n+1,k]_q = [n,k-1]_q + q^k [n,k]_q
+-- Cameron, Combinatorics, p126
+
+
+
+-- FLATS VIA REDUCED ROW ECHELON FORMS
+-- Suggested by Cameron p125
+
+
+data ZeroOneStar = Zero | One | Star deriving (Eq)
+
+instance Show ZeroOneStar where
+    show Zero = "0"
+    show One  = "1"
+    show Star = "*"
+
+-- reduced row echelon forms
+rrefs n k = map (rref 1 1) (combinationsOf k [1..n]) where
+    rref r c (x:xs) =
+        if c == x
+        then zipWith (:) (oneColumn r) (rref (r+1) (c+1) xs)
+        else zipWith (:) (starColumn r) (rref r (c+1) (x:xs))
+    rref _ c [] = replicate k (replicate (n+1-c) Star)
+    oneColumn r = replicate (r-1) Zero ++ One : replicate (k-r) Zero
+    starColumn r = replicate (r-1) Star ++ replicate (k+1-r) Zero
+
+flatsPG n fq k = concatMap substStars $ rrefs (n+1) (k+1) where
+    substStars (r:rs) = [r':rs' | r' <- substStars' r, rs' <- substStars rs]
+    substStars [] = [[]]
+    substStars' (Star:xs) = [x':xs' | x' <- fq, xs' <- substStars' xs]
+    substStars' (Zero:xs) = map (0:) $ substStars' xs
+    substStars' (One:xs) = map (1:) $ substStars' xs
+    substStars' [] = [[]]
+
+
+-- Flats in AG(n,Fq) are just the flats in PG(n,Fq) which are not "at infinity"
+flatsAG n fq k = [map tail (r : map (r <+>) rs) | r:rs <- flatsPG n fq k, head r == 1]
+-- The head r == 1 condition is saying that we want points which are in the "finite" part of PG(n,Fq), not points at infinity
+-- The reason we add r to each of the rs is to bring them into the "finite" part
+-- (If you don't do this, it can lead to incorrect results, for example some of the flats having the same closure)
+
+
+ Math/Combinatorics/Graph.hs view
@@ -0,0 +1,256 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Combinatorics.Graph where
+
+import qualified Data.List as L
+import Data.Maybe (isJust)
+import qualified Data.Map as M
+import qualified Data.Set as S
+import Control.Arrow ( (&&&) )
+
+import Math.Common.ListSet
+import Math.Algebra.Group.PermutationGroup
+import Math.Algebra.Group.SchreierSims as SS
+
+-- Main source: Godsil & Royle, Algebraic Graph Theory
+
+
+-- COMBINATORICS
+-- Some functions we'll use
+
+set xs = map head $ L.group $ L.sort xs
+
+-- subsets of a set (returned in "binary" order)
+powerset [] = [[]]
+powerset (x:xs) = let p = powerset xs in p ++ map (x:) p
+
+-- subsets of size k (returned in ascending order)
+combinationsOf 0 _ = [[]]
+combinationsOf _ [] = []
+combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs
+
+
+-- GRAPH
+
+data Graph a = G [a] [[a]] deriving (Eq,Ord,Show)
+
+-- we require that vs, es, and each individual e are sorted
+isSetSystem xs bs = isListSet xs && isListSet bs && all isListSet bs && all (`isSubset` xs) bs
+
+isGraph vs es = isSetSystem vs es && all ( (==2) . length) es
+
+graph (vs,es) | isGraph vs es = G vs es
+-- isValid g = g where g = G vs es
+
+toGraph (vs,es) | isGraph vs' es' = G vs' es' where
+    vs' = L.sort vs
+    es' = L.sort $ map L.sort es
+-- note that calling isListSet on a sorted list still checks that there are no duplicates
+
+vertices (G vs _) = vs
+
+edges (G _ es) = es
+
+
+-- OTHER REPRESENTATIONS
+
+-- incidence matrix of a graph
+-- (rows and columns indexed by edges and vertices respectively)
+-- (warning: in the literature it is often the other way round)
+incidenceMatrix (G vs es) = [ [if v `elem` e then 1 else 0 | v <- vs] | e <- es]
+
+fromIncidenceMatrix m = graph (vs,es) where
+    n = L.genericLength $ head m
+    vs = [1..n]
+    es = L.sort $ map edge m
+    edge row = [v | (1,v) <- zip row vs]
+
+adjacencyMatrix (G vs es) =
+    [ [if L.sort [i,j] `S.member` es' then 1 else 0 | j <- vs] | i <- vs]
+    where es' = S.fromList es
+
+fromAdjacencyMatrix m = graph (vs,es) where
+    n = L.genericLength m
+    vs = [1..n]
+    es = es' 1 m
+    es' i (r:rs) = [ [i,j] | (j,1) <- drop i (zip vs r)] ++ es' (i+1) rs
+    es' _ [] = []
+
+-- SOME SIMPLE FAMILIES OF GRAPHS
+
+nullGraph :: Graph Int -- type signature needed
+nullGraph = G [] []
+
+-- cyclic graph
+c n = graph (vs,es) where
+    vs = [1..n]
+    es = L.insert [1,n] [[i,i+1] | i <- [1..n-1]]
+-- automorphism group is D2n
+
+-- complete graph
+k n = graph (vs,es) where
+    vs = [1..n]
+    es = [[i,j] | i <- [1..n-1], j <- [i+1..n]] -- == combinationsOf 2 [1..n]
+-- automorphism group is Sn
+
+-- complete bipartite graph
+kb m n = to1n $ kb' m n
+
+kb' m n = graph (vs,es) where
+    vs = map Left [1..m] ++ map Right [1..n]
+    es = [ [Left i, Right j] | i <- [1..m], j <- [1..n] ]
+-- automorphism group is Sm*Sn (plus a flip if m==n)
+
+-- k-cube
+q k = let vs = zip [0..] (powerset [1..k])
+          es = [ [i,j] | (i,iset) <- vs, (j,jset) <- vs, i < j, length (iset `symDiff` jset) == 1 ]
+      in graph (map fst vs,es)
+
+q' k = let us = powerset $ map (2^) [0..k-1]
+           vs = [0..2^k-1] -- == L.sort $ map sum us
+           es = L.sort [ L.sort [sum u, sum v] | [u,v] <- combinationsOf 2 us, length (u `symDiff` v) == 1 ]
+       in graph (vs, es)
+
+tetrahedron = k 4
+
+cube = q 3
+
+octahedron = graph (vs,es) where
+    vs = [1..6]
+    es = combinationsOf 2 vs L.\\ [[1,6],[2,5],[3,4]]
+
+dodecahedron = toGraph (vs,es) where
+    vs = [1..20]
+    es = [ [1,2],[2,3],[3,4],[4,5],[5,1],
+           [6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15],[15,6],
+           [16,17],[17,18],[18,19],[19,20],[20,16],
+           [1,6],[2,8],[3,10],[4,12],[5,14],
+           [7,16],[9,17],[11,18],[13,19],[15,20] ]
+
+icosahedron = toGraph (vs,es) where
+    vs = [1..12]
+    es = [ [1,2],[1,3],[1,4],[1,5],[1,6],
+           [2,3],[3,4],[4,5],[5,6],[6,2],
+           [7,12],[8,12],[9,12],[10,12],[11,12],
+           [7,8],[8,9],[9,10],[10,11],[11,7],
+           [2,7],[7,3],[3,8],[8,4],[4,9],[9,5],[5,10],[10,6],[6,11],[11,2] ]
+
+
+
+-- 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]
+    vs' = M.elems mapping
+    es' = [map (mapping M.!) e | e <- es] -- the edges will already be sorted correctly by construction
+
+
+-- NEW GRAPHS FROM OLD
+
+complement (G vs es) = graph (vs,es') where es' = combinationsOf 2 vs \\ es
+-- es' = [e | e <- combinationsOf 2 vs, e `notElem` es]
+
+lineGraph g = to1n $ lineGraph' g
+
+lineGraph' (G vs es) = graph (es, [ [ei,ej] | ei <- es, ej <- dropWhile (<= ei) es, ei `intersect` ej /= [] ])
+
+petersen = complement $ lineGraph $ k 5
+
+
+-- SIMPLE PROPERTIES OF GRAPHS
+
+order g = length (vertices g)
+
+size g = length (edges g)
+
+-- also called degree
+valency (G vs es) v = length $ filter (v `elem`) es
+
+valencies g@(G vs es) = map (head &&& length) $ L.group $ L.sort $ map (valency g) vs
+
+regularParam g =
+    case valencies g of
+    [(v,_)] -> Just v
+    _       -> Nothing
+
+isRegular g = isJust $ regularParam g
+
+isCubic g = regularParam g == Just 3
+
+
+nbrs (G vs es) v = [u | [u,v'] <- es, v == v']
+                ++ [w | [v',w] <- es, v == v']
+-- if the graph is valid, then the neighbours will be returned in ascending order
+
+
+-- find paths from x to y using bfs
+-- by definition, a path is a subgraph isomorphic to a "line" - it can't have self-crossings
+-- (a walk allows self-crossings, a trail allows self-crossings but no edge reuse)
+findPaths g@(G vs es) x y = map reverse $ bfs [ [x] ] where
+    bfs ((z:zs) : nodes)
+        | z == y    = (z:zs) : bfs nodes
+        | otherwise = bfs (nodes ++ [(w:z:zs) | w <- nbrs g z, w `notElem` zs])
+    bfs [] = []
+
+-- length of the shortest path from x to y
+distance g x y =
+    case findPaths g x y of
+    [] -> -1 -- infinite
+    p:ps -> length p - 1
+
+-- diameter of a graph is maximum distance between two distinct vertices
+diameter g@(G vs es)
+    | isConnected g = maximum $ map maxDistance vs
+    | otherwise = -1
+    where maxDistance v = length (distancePartition g v) - 1
+
+-- find cycles starting at x
+-- by definition, a cycle is a subgraph isomorphic to a cyclic graph - it can't have self-crossings
+-- (a circuit allows self-crossings but not edge reuse)
+findCycles g@(G vs es) x = [reverse (x:z:zs) | z:zs <- bfs [ [x] ], z `elem` nbrsx, length zs > 1] where
+    nbrsx = nbrs g x
+    bfs ((z:zs) : nodes) = (z:zs) : bfs (nodes ++ [ w:z:zs | w <- nbrs g z, w `notElem` zs])
+    bfs [] = []
+
+-- girth of a graph is the size of the smallest cycle it contains
+-- Note: If graph contains no cycles, we return -1, representing infinity
+girth g@(G vs es) = minimum' $ map minCycle vs where
+    minimum' xs = let (zs,nzs) = L.partition (==0) xs in if null nzs then -1 else minimum nzs
+    minCycle v = case findCycles g v of
+                 [] -> 0
+                 c:cs -> length c - 1 -- because v occurs twice in c, as startpoint and endpoint
+
+-- circumference = max cycle - Bollobas p104
+
+
+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'
+
+-- the connected component to which v belongs
+component g v = concat $ distancePartition g v
+
+isConnected g@(G (v:vs) es) = length (component g v) == length (v:vs)
+isConnected (G [] []) = True
+
+
+-- MORE GRAPHS
+
+
+-- Generalized Johnson graph, Godsil & Royle p9
+j v k i | v >= k && k >= i
+    = graph (vs,es) where
+        vs = combinationsOf k [1..v]
+        es = [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, length (v1 `intersect` v2) == i ]
+-- j v k i is isomorphic to j v (v-k) (v-2k+i), so may as well have v >= 2k
+
+kneser v k | v >= 2*k = j v k 0
+
+johnson v k | v >= 2*k = j v k (k-1)
+
+petersen1 = to1n $ j 5 2 0
+
+
+
+ Math/Combinatorics/GraphAuts.hs view
@@ -0,0 +1,311 @@+-- Copyright (c) David Amos, 2009. All rights reserved.
+
+module Math.Combinatorics.GraphAuts where
+
+import qualified Data.List as L
+import qualified Data.Map as M
+import qualified Data.Set as S
+
+import Math.Common.ListSet
+import Math.Combinatorics.Graph
+-- import Math.Combinatorics.StronglyRegularGraph
+-- import Math.Combinatorics.Hypergraph -- can't import this, creates circular dependency
+import Math.Algebra.Group.PermutationGroup
+import Math.Algebra.Group.SchreierSims as SS
+
+
+-- The code for finding automorphisms - "graphAuts" - follows later on in file
+
+
+-- TRANSITIVITY PROPERTIES OF GRAPHS
+
+isVertexTransitive (G [] []) = True -- null graph is trivially vertex transitive
+isVertexTransitive g@(G (v:vs) es) = orbitP auts v == v:vs where
+    auts = graphAuts g
+
+isEdgeTransitive (G _ []) = True
+isEdgeTransitive g@(G vs (e:es)) = orbitB auts e == e:es where
+    auts = graphAuts g
+
+arc ->^ g = map (.^ g) arc
+-- unlike blocks, arcs are directed, so the action on them does not sort
+
+-- Godsil & Royle 59-60
+isArcTransitive (G _ []) = True -- empty graphs are trivially arc transitive
+isArcTransitive g@(G vs es) = orbit (->^) a auts == a:as where
+    a:as = L.sort $ es ++ map reverse es
+    auts = graphAuts g
+
+isArcTransitive' g@(G (v:vs) es) =
+    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
+          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
+
+
+-- then k n, kb n n, q n, other platonic solids, petersen graph, heawood graph, pappus graph, desargues graph are all arc-transitive
+
+
+-- find arcs of length l from x using dfs - results returned in order
+-- an arc is a sequence of vertices connected by edges, no doubling back, but self-crossings allowed
+findArcs g@(G vs es) x l = map reverse $ dfs [ ([x],0) ] where
+    dfs ( (z1:z2:zs,l') : nodes)
+        | l == l'   = (z1:z2:zs) : dfs nodes
+        | otherwise = dfs $ [(w:z1:z2:zs,l'+1) | w <- nbrs g z1, w /= z2] ++ nodes
+    dfs ( ([z],l') : nodes)
+        | l == l'   = [z] : dfs nodes
+        | otherwise = dfs $ [([w,z],l'+1) | w <- nbrs g z] ++ nodes
+    dfs [] = []
+
+-- 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
+
+isnArcTransitive _ (G [] []) = True
+isnArcTransitive n g@(G (v:vs) es) =
+    orbitP auts v == v:vs && -- isVertexTransitive g
+    orbit (->^) a stab == a:as
+    where auts = graphAuts g
+          stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
+          a:as = findArcs g v n
+
+is2ArcTransitive g = isnArcTransitive 2 g
+
+is3ArcTransitive g = isnArcTransitive 3 g
+
+-- Godsil & Royle 66-7
+isDistanceTransitive (G [] []) = True
+isDistanceTransitive g@(G (v:vs) es)
+    | 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"
+    where auts = graphAuts g
+          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
+
+-- refine one partition by another
+refine p1 p2 = concat [ [c1 `intersect` c2 | c2 <- p2] | c1 <- p1]
+-- Refinement preserves ordering within cells but not between cells
+-- eg the cell [1,2,3,4] could be refined to [2,4],[1,3]
+
+
+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
+
+
+-- 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
+
+-- 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
+
+-- 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)
+            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 []
+               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
+
+isSingleton [_] = True
+isSingleton _ = False
+
+
+-- 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
+graphAuts 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)
+        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]
+    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
+
+-- 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
+
+
+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
+
+
+-- 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]
+-}
+
+ Math/Combinatorics/Hypergraph.hs view
@@ -0,0 +1,214 @@+-- Copyright (c) David Amos, 2009. All rights reserved.
+
+module Math.Combinatorics.Hypergraph where
+
+import qualified Data.List as L
+import Math.Common.ListSet
+import Math.Combinatorics.Graph hiding (incidenceMatrix)
+import Math.Algebra.Group.PermutationGroup (orbitB, p) -- needed for construction of Coxeter group
+
+-- not used in this module, only in GHCi
+import Math.Algebra.Field.Base
+import Math.Algebra.Field.Extension
+import Math.Combinatorics.Design hiding (incidenceMatrix, incidenceGraph, dual, isSubset, fanoPlane)
+
+
+-- set system or hypergraph
+data Hypergraph a = H [a] [[a]] deriving (Eq,Ord,Show)
+
+hypergraph xs bs | isSetSystem xs bs = H xs bs
+
+toHypergraph xs bs = H xs' bs' where
+    xs' = L.sort xs
+    bs' = L.sort $ map L.sort bs
+-- this still doesn't guarantee that all bs are subset of xs
+
+
+-- uniform hypergraph - all blocks are same size
+isUniform h@(H xs bs) = isSetSystem xs bs && same (map length bs)
+
+same (x:xs) = all (==x) xs
+same [] = True
+
+
+
+fromGraph (G vs es) = H vs es
+fromDesign (D xs bs) = H xs (L.sort bs)
+-- !! should insist that designs have blocks in order
+-- !! dual probably doesn't guarantee this at present
+
+{-
+dual (H xs bs) = toHypergraph (bs, map beta xs) where
+    beta x = filter (x `elem`) bs
+-}
+
+
+
+-- INCIDENCE GRAPH
+
+data Incidence a = P a | B [a] deriving (Eq, Ord, Show)
+
+-- compare Design, where we just use Left, Right
+
+-- Also called the Levi graph
+incidenceGraph (H xs bs) = G vs es where
+    vs = map P xs ++ map B bs
+    es = L.sort [ [P x, B b] | b <- bs, x <- b]
+
+
+-- INCIDENCE MATRIX
+
+-- !! why are we doing this the other way round to the literature ??
+
+
+-- incidence matrix of a hypergraph
+-- (rows and columns indexed by edges and vertices respectively)
+-- (warning: in the literature it is often the other way round)
+incidenceMatrix (H vs es) = [ [if v `elem` e then 1 else 0 | v <- vs] | e <- es]
+
+fromIncidenceMatrix m = H vs es where
+    n = L.genericLength $ head m
+    vs = [1..n]
+    es = L.sort $ map edge m
+    edge row = [v | (1,v) <- zip row vs]
+
+
+
+
+-- isTwoGraph
+
+
+-- We can represent various incidence structures as hypergraphs,
+-- by identifying the lines with the sets of points that they contain
+
+isPartialLinearSpace h@(H ps ls) =
+    isSetSystem ps ls &&
+    all ( (<=1) . length ) [filter (pair `isSubset`) ls | pair <- combinationsOf 2 ps]
+    -- any two points are incident with at most one line
+
+-- Godsil & Royle, p79
+isProjectivePlane h@(H ps ls) =
+    isSetSystem ps ls &&
+    all ( (==1) . length) [intersect l1 l2 | [l1,l2] <- combinationsOf 2 ls] && -- any two lines meet in a unique point
+    all ( (==1) . length) [ filter ([p1,p2] `isSubset`) ls | [p1,p2] <- combinationsOf 2 ps] -- any two points lie in a unique line
+
+-- a projective plane with a triangle
+-- this is a weak non-degeneracy condition, which eliminates all points on the same line, or all lines through the same point
+isProjectivePlaneTri h@(H ps ls) =
+    isProjectivePlane h && any triangle (combinationsOf 3 ps)
+    where triangle t@[p1,p2,p3] =
+                   (not . null) [l | l <- ls, [p1,p2] `isSubset` l, p3 `notElem` l] && -- there is a line containing p1,p2 but not p3
+                   (not . null) [l | l <- ls, [p1,p3] `isSubset` l, p2 `notElem` l] &&
+                   (not . null) [l | l <- ls, [p2,p3] `isSubset` l, p1 `notElem` l] 
+
+-- a projective plane with a quadrangle
+-- this is a stronger non-degeneracy condition
+isProjectivePlaneQuad h@(H ps ls) =
+    isProjectivePlane h && any quadrangle (combinationsOf 4 ps)
+    where quadrangle q = all (not . collinear) (combinationsOf 3 q) -- no three points collinear
+          collinear ps = any (ps `isSubset`) ls
+
+
+-- > isProjectivePlaneQuad $ fromDesign $ pg2 f2
+-- True
+
+
+-- GENERALIZED QUADRANGLES
+
+-- Godsil & Royle p81
+isGeneralizedQuadrangle h@(H ps ls) =
+    isPartialLinearSpace h &&
+    all (\(l,p) -> unique [p' | p' <- l, collinear (pair p p')]) [(l,p) | l <- ls, p <- ps, p `notElem` l] &&
+    -- given any line l and point p not on l, there is a unique point p' on l with p and p' collinear
+    any (not . collinear) (powerset ps) && -- there are non collinear points
+    any (not . concurrent) (powerset ls) -- there are non concurrent lines
+    where unique xs = length xs == 1
+          pair x y = if x < y then [x,y] else [y,x]
+          collinear ps = any (ps `isSubset`) ls
+          concurrent ls = any (\p -> all (p `elem`) ls) ps
+
+
+grid m n = H ps ls where
+    ps = [(i,j) | i <- [1..m], j <- [1..n] ]
+    ls = L.sort $ [ [(i,j) | i <- [1..m] ] | j <- [1..n] ] -- horizontal lines
+               ++ [ [(i,j) | j <- [1..n] ] | i <- [1..m] ] -- vertical lines
+
+dualGrid m n = fromGraph $ kb m n
+-- the lines of the grid are the points of the dual, and the points of the grid are the lines of the dual
+
+isGenQuadrangle' h = diameter g == 4 && girth g == 8 -- !! plus non-degeneracy conditions
+    where g = incidenceGraph h
+
+
+-- CONFIGURATIONS
+
+-- http://en.wikipedia.org/wiki/Projective_configuration
+isConfiguration h@(H ps ls) =
+    isUniform h && -- a set system, with each line incident with the same number of points
+    same [length (filter (p `elem`) ls) | p <- ps] -- each point is incident with the same number of lines
+
+
+
+
+fanoPlane = toHypergraph [1..7] [[1,2,4],[2,3,5],[3,4,6],[4,5,7],[5,6,1],[6,7,2],[7,1,3]]
+
+heawoodGraph = incidenceGraph fanoPlane
+
+
+desarguesConfiguration = H xs bs where
+    xs = combinationsOf 2 [1..5]
+    bs = [ [x | x <- xs, x `isSubset` b] | b <- combinationsOf 3 [1..5] ]
+
+desarguesGraph = incidenceGraph desarguesConfiguration
+
+
+pappusConfiguration = H xs bs where
+    xs = [1..9]
+    bs = L.sort [ [1,2,3], [4,5,6], [7,8,9], [1,5,9], [1,6,8], [2,4,9], [3,4,8], [2,6,7], [3,5,7] ]
+
+pappusGraph = incidenceGraph pappusConfiguration
+
+
+
+-- !! no particular reason why the following is here rather than elsewhere
+{-
+triples = combinationsOf 3 [1..7]
+
+heptads = [ [a,b,c,d,e,f,g] | a <- triples,
+                              b <- triples, a < b, meetOne b a,
+                              c <- triples, b < c, all (meetOne c) [a,b],
+                              d <- triples, c < d, all (meetOne d) [a,b,c],
+                              e <- triples, d < e, all (meetOne e) [a,b,c,d],
+                              f <- triples, e < f, all (meetOne f) [a,b,c,d,e],
+                              g <- triples, f < g, all (meetOne g) [a,b,c,d,e,f],
+                              foldl intersect [1..7] [a,b,c,d,e,f,g] == [] ]
+    where meetOne x y = length (intersect x y) == 1
+    -- each pair of triples meet in exactly one point, and there is no point in all of them - Godsil & Royle p69
+    -- (so these are the projective planes over 7 points)
+-}
+-- Godsil & Royle p69
+coxeterGraph = G vs es where
+    g = p [[1..7]]
+    vs = L.sort $ concatMap (orbitB [g]) [[1,2,4],[3,5,7],[3,6,7],[5,6,7]]
+    es = [ e | e@[v1,v2] <- combinationsOf 2 vs, disjoint v1 v2]
+
+-- is this the incidence graph of a hypergraph involving heptads over triples?
+
+
+-- edges of K6
+duads = combinationsOf 2 [1..6]
+
+-- 1-factors of K6
+-- 15 different ways to pick three disjoint duads from [1..6]
+synthemes = [ [d1,d2,d3] | d1 <- duads,
+                           d2 <- duads, d2 > d1, disjoint d1 d2,
+                           d3 <- duads, d3 > d2, disjoint d1 d3, disjoint d2 d3 ]
+
+-- Tutte 8-cage
+tutteCoxeterGraph = incidenceGraph $ H duads synthemes
+
+
+-- Also known as line graph
+intersectionGraph (H xs bs) = G vs es where
+    vs = bs
+    es = [pair | pair@[b1,b2] <- combinationsOf 2 bs, not (disjoint b1 b2)]
+ Math/Combinatorics/StronglyRegularGraph.hs view
@@ -0,0 +1,280 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Combinatorics.StronglyRegularGraph where
+
+import qualified Data.List as L
+import Data.Maybe (isJust)
+import qualified Data.Map as M
+import qualified Data.Set as S
+
+import Math.Common.ListSet
+import Math.Algebra.Group.PermutationGroup hiding (P)
+import Math.Algebra.Group.SchreierSims as SS
+import Math.Combinatorics.Graph as G hiding (G)
+import Math.Combinatorics.GraphAuts
+import Math.Combinatorics.Design as D
+import Math.Algebra.LinearAlgebra -- hiding (t)
+import Math.Algebra.Field.Base -- for F2
+import Math.Combinatorics.FiniteGeometry hiding (combinationsOf)
+
+-- Sources
+-- Godsil & Royle, Algebraic Graph Theory
+-- Cameron & van Lint, Designs, Graphs, Codes and their Links
+-- van Lint & Wilson, A Course in Combinatorics, 2nd ed
+
+
+-- STRONGLY REGULAR GRAPHS
+
+-- strongly regular graphs
+srgParams g =
+    let vs = vertices g
+        n = length vs
+        es = edges g
+        es' = combinationsOf 2 vs \\ es -- the non-edges
+        k:ks = map (valency g) vs
+        lambda:ls = map (length . commonNbrs) es  -- common neighbours of adjacent vertices
+        mu:ms = map (length . commonNbrs) es' -- common neighbours of non-adjacent vertices
+        commonNbrs [v1,v2] = nbrs g v1 `intersect` nbrs g v2
+    in if all (==k) ks && all (==lambda) ls && all (==mu) ms
+       then Just (n,k,lambda,mu)
+       else Nothing
+
+isSRG g = isJust $ srgParams g
+
+
+-- SIMPLE EXAMPLES
+
+-- Triangular graph - van Lint & Wilson p262
+-- http://mathworld.wolfram.com/TriangularGraph.html
+t m = G.to1n $ t' m
+
+t' m | m >= 4 = graph (vs,es) where
+    vs = combinationsOf 2 [1..m]
+    es = [ [v,v'] | v <- vs, v' <- dropWhile (<= v) vs, not (disjoint v v')]
+-- This is just lineGraph (k m), by another name
+
+petersen = complement $ t 5
+
+-- Lattice graph - van Lint & Wilson p262
+-- http://mathworld.wolfram.com/LatticeGraph.html
+l2 m = G.to1n $ l2' m
+
+l2' m | m >= 2 = graph (vs,es) where
+    vs = [ (i,j) | i <- [1..m], j <- [1..m] ]
+    es = [ [v,v'] | v@(i,j) <- vs, v'@(i',j') <- dropWhile (<= v) vs, i == i' || j == j']
+-- This is lineGraph (kb m m)
+-- Automorphism group is Sm * Sm * C2
+-- via i -> ig, j -> jg, i <-> j
+
+paleyGraph fq | length fq `mod` 4 == 1 = graph (vs,es) where
+    vs = fq
+    qs = set [x^2 | x <- vs] \\ [0] -- the non-zero squares in Fq
+    es = [ [x,y] | x <- vs, y <- vs, x < y, (x-y) `elem` qs]
+
+
+-- CLEBSCH GRAPH
+-- van Lint & Wilson, p263
+
+clebsch = G.to1n clebsch'
+
+clebsch' = graph (vs,es) where
+    vs = L.sort $ filter (even . length) $ powerset [1..5]
+    es = [ [v,v'] | v <- vs, v' <- dropWhile (<= v) vs, length (symDiff v v') == 4]
+
+-- Alternative construction from Cameron & van Lint p106
+clebsch2 = graph (vs,es) where
+    D xs bs = pairDesign 5
+    vs = [C] ++ [P x | x <- xs] ++ [B b | b <- bs]
+    es = L.sort $ [ [B a, B b]  | a <- bs, b <- dropWhile (<=a) bs, disjoint a b]
+               ++ [ [P p, B b] | b <- bs, p <- b]
+               ++ [ [C, P p]   | p <- xs ]
+
+
+-- HOFFMAN-SINGLETON GRAPH
+-- Cameron, Permutation Groups, p79ff
+-- Godsil & Royle, p92ff
+-- Aut group is U3(5).2 (Atlas p34)
+
+triples = combinationsOf 3 [1..7]
+
+heptads = [ [a,b,c,d,e,f,g] | a <- triples,
+                              b <- triples, a < b, meetOne b a,
+                              c <- triples, b < c, all (meetOne c) [a,b],
+                              d <- triples, c < d, all (meetOne d) [a,b,c],
+                              e <- triples, d < e, all (meetOne e) [a,b,c,d],
+                              f <- triples, e < f, all (meetOne f) [a,b,c,d,e],
+                              g <- triples, f < g, all (meetOne g) [a,b,c,d,e,f],
+                              foldl intersect [1..7] [a,b,c,d,e,f,g] == [] ]
+    where meetOne x y = length (intersect x y) == 1
+    -- each pair of triples meet in exactly one point, and there is no point in all of them - Godsil & Royle p69
+    -- (so these are the projective planes over 7 points)
+
+plane +^ g = L.sort [line -^ g | line <- plane]
+
+plane +^^ gs = orbit (+^) plane gs
+
+hoffmanSingleton = G.to1n hoffmanSingleton'
+
+hoffmanSingleton' = graph (vs,es) where
+    h = head heptads
+    hs = h +^^ _A 7 -- an A7 orbit of a heptad
+    vs = map Left hs ++ map Right triples
+    es = [ [Left h, Right t] | h <- hs, t <- triples, t `elem` h]
+      ++ [ [Right t, Right t'] | t <- triples, t' <- dropWhile (<= t) triples, t `disjoint` t']
+
+-- induced action of A7 on Hoffman-Singleton graph
+inducedA7 g = fromPairs [(v, v ~^ g) | v <- vs] where
+    vs = vertices hoffmanSingleton'
+    (Left h) ~^ g = Left (h +^ g)
+    (Right t) ~^ g = Right (t -^ g)
+
+hsA7 = toSn $ map inducedA7 $ _A 7
+
+
+-- GEWIRTZ GRAPH
+-- van Lint & Wilson p266-7
+-- (also called Sims-Gewirtz graph)
+
+gewirtz = G.to1n gewirtz'
+
+gewirtz' = graph (vs,es) where
+    vs = [xs | xs <- blocks s_3_6_22, 22 `notElem` xs]
+    -- The 21 blocks of S(3,6,22) which contain 22 are the lines of PG(2,4) (projective plane over F4)
+    -- The 56 blocks which don't are hyperovals in this plane. They form a 2-(21,6,4) design.
+    es = [ [v,v'] | v <- vs, v' <- dropWhile (<= v) vs, length (v `intersect` v') == 0]
+
+
+-- HIGMAN-SIMS GRAPH
+-- Aut group is HS.2, where HS is the Higman-Sims sporadic simple group
+
+data DesignVertex = C | P Integer | B [Integer] deriving (Eq,Ord,Show)
+
+higmanSimsGraph = G.to1n higmanSimsGraph'
+
+-- Cameron & van Lint, p107
+higmanSimsGraph' = graph (vs,es) where
+    D xs bs = s_3_6_22
+    vs = [C] ++ [P x | x <- xs] ++ [B b | b <- bs]
+    es = L.sort $ [ [B a, B b]  | a <- bs, b <- dropWhile (<=a) bs, disjoint a b]
+               ++ [ [P p, B b] | b <- bs, p <- b]
+               ++ [ [C, P p]   | p <- xs ]
+    -- s_3_6_22' = blocks s_3_6_22
+
+-- There is an induced action of M22 on Higman Sims graph
+
+-- induced action of M22 on Higman-Sims graph
+inducedM22 g = fromPairs [(v, v ~^ g) | v <- vs] where
+    -- G vs _ = higmanSimsGraph'
+    vs = vertices higmanSimsGraph'
+    (B b) ~^ g = B (b -^ g)
+    (P p) ~^ g = P (p .^ g)
+    C ~^ _ = C
+
+higmanSimsM22 = toSn $ map inducedM22 $ m22sgs
+-- all (isGraphAut higmanSimsGraph) higmanSimsM22
+
+-- M22 is one point stabilizer (of C)
+
+-- HS.2, where HS is Higman-Sims sporadic group
+_HS2 = SS.reduceGens $ graphAuts higmanSimsGraph
+-- (It will actually find 11 strong generators, but the first 4 are sufficient to generate the group)
+
+_HS = SS.derivedSubgp _HS2
+
+
+
+-- SYMPLECTIC GRAPHS
+
+-- Godsil & Royle p242
+sp2 r = graph (vs,es) where
+    vs = tail $ ptsAG (2*r) f2 -- all non-zero pts in F2^2r
+    es = [ [u,v] | [u,v] <- combinationsOf 2 vs, u <*>> n <.> v == 1] -- uT N v == 1, ie vectors adjacent if non-orthogonal
+    n = fMatrix (2*r) (\i j -> if abs (i-j) == 1 && even (max i j) then 1 else 0) -- matrix defining a symplectic form
+
+sp n | even n = sp2 (n `div` 2)
+
+
+-- TWO GRAPHS AND SWITCHING
+
+-- SCHLAFLI GRAPH
+-- An srg(27,16,10,8)
+-- Has geometric interpretation in terms of 27 lines on general cubic surface in projective 3-space
+-- Aut group is G.2 where G = U4(2) = S4(3) (Atlas p26)
+-- (G.2 is also the Weyl group of E6 - don't know if there's any connection)
+
+-- Godsil & Royle p254ff
+switch g us | us `D.isSubset` vs = graph (vs, L.sort switchedes) where
+    vs = vertices g
+    us' = vs L.\\ us -- complement of us in vs
+    es = edges g
+    es' = S.fromList es
+    switchedes = [e | e@[v1,v2] <- es, (v1 `elem` us) == (v2 `elem` us)]
+                 -- edges within us or its complement are unchanged
+              ++ [ L.sort [v1,v2] | v1 <- us, v2 <- us', L.sort [v1,v2] `S.notMember` es']
+                 -- edges between us and its complement are switched
+
+-- Godsil & Royle p259
+schlafli = G.to1n schlafli'
+
+schlafli' = graph (vs,es') where
+    g = lineGraph $ k 8
+    v:vs = vertices g
+    es = edges g
+    gswitched = switch g (nbrs g v) -- switch off the vertex v
+    es' = edges gswitched
+
+
+-- MCLAUGHLIN GRAPH
+-- Aut group is McL.2, where McL is the McLaughlin sporadic simple group
+-- http://people.csse.uwa.edu.au/gordon/constructions/mclaughlin/
+-- http://mathworld.wolfram.com/McLaughlinGraph.html
+
+mcLaughlin = G.to1n mcLaughlin'
+
+mcLaughlin' = graph (vs',es') where
+    D xs bs = s_4_7_23
+    vs = map P xs ++ map B bs
+    es = [ [P x, B b] | x <- xs, b <- bs, x `notElem` b]
+      ++ [ [B b1, B b2] | b1 <- bs, b2 <- bs, b1 < b2, length (b1 `intersect` b2) == 1]
+    g276 = graph (vs,es)
+    g276switched = switch g276 (nbrs g276 (P 0))
+    P 0 : vs' = vs -- drop P 0 as it's now not connected
+    es' = edges g276switched
+
+_McL2 = SS.reduceGens $ graphAuts mcLaughlin
+-- finds 14 auts - but takes half an hour (interpreted) to do so
+-- in fact just the first 2 are sufficient to generate the group
+
+_McL = SS.derivedSubgp $ _McL2
+
+{-
+-- TWO GRAPH ON 276 VERTICES
+-- Has Conway's .3 as automorphism group
+
+-- Godsil & Royle p260ff
+
+twoGraph276 =
+    let nt = D.incidenceMatrix s_4_7_23
+        n = L.transpose nt -- Godsil & Royle do incidence matrix the other way round to us
+        s = L.transpose $
+            (j 23 23 <<->> i 23)      +|+ (j 23 253 <<->> 2 *>> n)
+                                      ++
+            (j 253 23 <<->> 2 *>> nt) +|+ (nt <<*>> n <<->> 5 *>> i 253 <<->> 2 *>> j 253 253)
+        a = (map . map) (`div` 2) (j 276 276 <<->> i 276 <<->> s)
+    in fromAdjacencyMatrix a
+    where j r c = replicate r (replicate c 1)
+          i = idMx
+          (+|+) = zipWith (++)
+
+-- Its automorphism group *as a two-graph* is .3 (Co3)
+-- But its aut group as a graph is only M23
+
+twoGraph276' = graph (vs,es) where
+    D xs bs = s_4_7_23
+    vs = map P xs ++ map B bs
+    es = [ [P x, B b] | x <- xs, b <- bs, x `notElem` b]
+      ++ [ [B b1, B b2] | b1 <- bs, b2 <- bs, b1 < b2, length (b1 `intersect` b2) == 1]
+-- !! This isn't isomorphic to twoGraph276
+-- (Perhaps it is in the same switching class though)
+-- We can obtain McLaughlin graph from this by switching in neighbourhood of P 0
+-}
+ Math/Common/IntegerAsType.hs view
@@ -0,0 +1,98 @@+-- Copyright (c) David Amos, 2009. All rights reserved.
+
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+module Math.Common.IntegerAsType where
+
+class IntegerAsType a where
+    value :: a -> Integer
+
+-- multiplication of IntegerAsType
+data M a b = M a b
+instance (IntegerAsType a, IntegerAsType b) => IntegerAsType (M a b) where
+    value _ = value (undefined :: a) * value (undefined :: b)
+
+data TMinus1
+instance IntegerAsType TMinus1 where value _ = -1
+
+data TZero
+instance IntegerAsType TZero where value _ = 0
+
+data TOne
+instance IntegerAsType TOne where value _ = 1
+
+data T2
+instance IntegerAsType T2 where value _ = 2
+
+data T3
+instance IntegerAsType T3 where value _ = 3
+
+data T5
+instance IntegerAsType T5 where value _ = 5
+
+data T7
+instance IntegerAsType T7 where value _ = 7
+
+data T11
+instance IntegerAsType T11 where value _ = 11
+
+data T13
+instance IntegerAsType T13 where value _ = 13
+
+data T17
+instance IntegerAsType T17 where value _ = 17
+
+data T19
+instance IntegerAsType T19 where value _ = 19
+
+data T23
+instance IntegerAsType T23 where value _ = 23
+
+data T29
+instance IntegerAsType T29 where value _ = 29
+
+data T31
+instance IntegerAsType T31 where value _ = 31
+
+data T37
+instance IntegerAsType T37 where value _ = 37
+
+data T41
+instance IntegerAsType T41 where value _ = 41
+
+data T43
+instance IntegerAsType T43 where value _ = 43
+
+data T47
+instance IntegerAsType T47 where value _ = 47
+
+data T53
+instance IntegerAsType T53 where value _ = 53
+
+data T59
+instance IntegerAsType T59 where value _ = 59
+
+data T61
+instance IntegerAsType T61 where value _ = 61
+
+data T67
+instance IntegerAsType T67 where value _ = 67
+
+data T71
+instance IntegerAsType T71 where value _ = 71
+
+data T73
+instance IntegerAsType T73 where value _ = 73
+
+data T79
+instance IntegerAsType T79 where value _ = 79
+
+data T83
+instance IntegerAsType T83 where value _ = 83
+
+data T89
+instance IntegerAsType T89 where value _ = 89
+
+data T97
+instance IntegerAsType T97 where value _ = 97
+
+ Math/Common/ListSet.hs view
@@ -0,0 +1,54 @@+
+
+module Math.Common.ListSet where
+
+import Data.List (group,sort)
+-- versions of Data.List functions which assume that the lists are ascending sets (no repeated elements)
+
+toListSet xs = map head $ group $ sort xs
+
+isListSet (x1:x2:xs) = x1 < x2 && isListSet (x2:xs)
+isListSet _ = True
+
+
+union (x:xs) (y:ys) =
+    case compare x y of
+    LT -> x : union xs (y:ys)
+    EQ -> x : union xs ys
+    GT -> y : union (x:xs) ys
+union xs ys = xs ++ ys -- one of them is null
+
+intersect (x:xs) (y:ys) =
+    case compare x y of
+    LT -> intersect xs (y:ys)
+    EQ -> x : intersect xs ys
+    GT -> intersect (x:xs) ys
+intersect _ _ = []
+
+(x:xs) \\ (y:ys) =
+    case compare x y of
+    LT -> x : (xs \\ (y:ys))
+    EQ -> xs \\ ys
+    GT -> (x:xs) \\ ys
+[] \\ _ = []
+xs \\ [] = xs
+
+symDiff (x:xs) (y:ys) =
+    case compare x y of
+    LT -> x : symDiff xs (y:ys)
+    EQ -> symDiff xs ys
+    GT -> y : symDiff (x:xs) ys
+symDiff xs ys = xs ++ ys -- one of them is null
+
+disjoint xs ys = null (intersect xs ys)
+
+isSubset (x:xs) (y:ys) =
+    case compare x y of
+    LT -> False
+    EQ -> isSubset xs ys
+    GT -> isSubset (x:xs) ys
+isSubset [] _ = True
+isSubset _ [] = False
+
+-- Note that an ListSet.elem turned out to be slower than Data.List.elem
+-- (Perhaps because it's slower when x `notElem` xs)
+ Math/Projects/ChevalleyGroup/Classical.hs view
@@ -0,0 +1,127 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Projects.ChevalleyGroup.Classical where
+
+import Math.Algebra.Field.Base
+import Math.Algebra.Field.Extension hiding ( (<+>), (<*>) )
+import Math.Algebra.LinearAlgebra
+
+import Math.Algebra.Group.PermutationGroup
+import Math.Algebra.Group.SchreierSims as SS
+
+import Math.Combinatorics.FiniteGeometry
+
+
+numPtsAG n q = q^n
+
+numPtsPG n q = (q^(n+1)-1) `div` (q-1)
+
+
+
+-- LINEAR GROUPS
+
+-- SL(n,Fq) is generated by elementary transvections
+sl n fq = [elemTransvection n (r,c) l | r <- [1..n], c <- [1..n], r /= c, l <- fq']
+    where fq' = basisFq undefined -- tail fq
+    -- Carter p68 - x_r(t1) x_r(t2) == x_r(t1+t2) - this is true in general, not just in this case
+
+elemTransvection n (r,c) l = fMatrix n (\i j -> if i == j then 1 else if (i,j) == (r,c) then l else 0)
+
+-- PSL(n,Fq) == A(n,Fq) == SL(n,Fq)/Z
+l n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- sl n fq]
+    where ps = ptsPG (n-1) fq
+
+orderL n q = ( q^(n*(n-1) `div` 2) * product [ q^i-1 | i <- [n,n-1..2] ] )
+              `div` gcd (q-1) n
+
+
+-- SYMPLECTIC GROUPS
+-- Carter p186 and 181-3
+
+sp2 n fq =
+    [_I <<+>> t *>> (e i j <<->> e (-j) (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+    [_I <<->> t *>> (e (-i) (-j) <<->> e j i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+    [_I <<+>> t *>> (e i (-j) <<+>> e j (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+    [_I <<+>> t *>> (e (-i) j <<+>> e (-j) i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++ -- Carter expresses this slightly differently
+    [_I <<+>> t *>> e i (-i) | i <- [1..n], t <- fq' ] ++
+    [_I <<+>> t *>> e (-i) i | i <- [1..n], t <- fq' ]
+    where
+        fq' = basisFq undefined -- tail fq -- multiplicative group
+        _I = idMx (2*n)
+        e i j = e' (if i > 0 then i else n-i) (if j > 0 then j else n-j)
+        e' i j = fMatrix (2*n) (\k l -> if (k,l) == (i,j) then 1 else 0)
+
+-- PSp2n(Fq) == Cn(Fq) == Sp2n(Fq)/Z
+s2 n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- sp2 n fq]
+    where ps = ptsPG (2*n-1) fq
+
+s n fq | even n = s2 (n `div` 2) fq
+
+
+orderS2 n q = (q^n^2 * product [ q^i-1 | i <- [2*n,2*n-2..2] ]) `div` gcd (q-1) 2
+
+orderS n q | even n = orderS2 (n `div` 2) q
+
+
+-- ORTHOGONAL GROUPS
+
+-- Carter p185 and 178-9
+-- Omega2n(q) - commutator subgroup of O2n(q)
+omegaeven n fq =
+    [_I <<+>> t *>> (e i j <<->> e (-j) (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+    [_I <<->> t *>> (e (-i) (-j) <<->> e j i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+    [_I <<+>> t *>> (e i (-j) <<->> e j (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+    [_I <<->> t *>> (e (-i) j <<->> e (-j) i) | i <- [1..n], j <- [i+1..n], t <- fq' ]
+    where
+        fq' = basisFq undefined -- tail fq -- multiplicative group
+        _I = idMx (2*n)
+        e i j = e' (if i > 0 then i else n-i) (if j > 0 then j else n-j)
+        e' i j = fMatrix (2*n) (\k l -> if (k,l) == (i,j) then 1 else 0)
+
+
+-- O+2n(Fq)  Artin/Conway notation (Atlas, pxii)
+-- Dn(Fq)    Chevalley group
+d n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- omegaeven n fq]
+    where ps = ptsPG (2*n-1) fq
+
+
+-- Carter p186-8
+-- Omega2n+1(q)
+omegaodd n fq
+    | char fq /= 2 =
+        [_I <<+>> t *>> (e i j <<->> e (-j) (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+        [_I <<->> t *>> (e (-i) (-j) <<->> e j i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+        [_I <<+>> t *>> (e i (-j) <<->> e j (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+        [_I <<->> t *>> (e (-i) j <<->> e (-j) i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+        [_I <<+>> t *>> (2 *>> e i 0 <<->> e 0 (-i)) <<->> (t^2) *>> e i (-i) | i <- [1..n], t <- fq' ] ++
+        [_I <<->> t *>> (2 *>> e (-i) 0 <<->> e 0 i) <<->> (t^2) *>> e (-i) i | i <- [1..n], t <- fq' ]
+    | char fq == 2 =
+        [_I <<+>> t *>> (e i j <<->> e (-j) (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+        [_I <<->> t *>> (e (-i) (-j) <<->> e j i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+        [_I <<+>> t *>> (e i (-j) <<->> e j (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++ -- !! Carter has a + in place of a - here
+        [_I <<->> t *>> (e (-i) j <<->> e (-j) i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
+        [_I <<+>> (t^2) *>> e i (-i) | i <- [1..n], t <- fq' ] ++
+        [_I <<+>> (t^2) *>> e (-i) i | i <- [1..n], t <- fq' ]
+    where
+        fq' = basisFq undefined -- tail fq -- multiplicative group
+        _I = idMx (2*n+1)
+        e i j = e' (if i >= 0 then i else n-i) (if j >= 0 then j else n-j)
+        e' i j = fMatrix' (2*n+1) (\k l -> if (k,l) == (i,j) then 1 else 0)
+
+-- O2n+1(Fq)  Artin/Conway notation
+-- Bn(Fq)     Chevalley group
+b n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- omegaodd n fq]
+    where ps = ptsPG (2*n) fq
+
+
+o n fq | even n = d (n `div` 2) fq
+       | odd  n = b (n `div` 2) fq
+
+-- The orthogonal groups aren't transitive on PG(n-1,Fq),
+-- so the above permutation representation actually splits into smaller representations on the orbits
+-- eg map length $ orbits $ o 7 f3 -> [364,378,351]
+-- which is the first three permutation representations listed at http://brauer.maths.qmul.ac.uk/Atlas/v3/clas/O73/
+
+
+-- UNITARY GROUPS
+-- The unitary group U(n+1,q) is the twisted Chevalley group 2An(q)
+ Math/Projects/ChevalleyGroup/Exceptional.hs view
@@ -0,0 +1,180 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Projects.ChevalleyGroup.Exceptional where
+
+import Data.List as L
+
+import Math.Algebra.Field.Base
+import Math.Algebra.Field.Extension hiding ( (<+>), (<*>) )
+import Math.Algebra.LinearAlgebra
+
+import Math.Algebra.Group.PermutationGroup
+import Math.Algebra.Group.SchreierSims as SS
+
+import Math.Combinatorics.FiniteGeometry (ptsAG)
+-- import ClassicalChevalleyGroup (ptsAG)
+
+
+-- Follows Conway's notation
+
+-- The octonion xinf + x0 i0 + x1 i1 + ... + x6 i6
+-- is represented as O [(-1,xinf),(0,x0),(1,x1),...,(6,x6)]
+-- where a list element may be omitted if one if the coefficient is zero
+
+
+newtype Octonion k = O [(Int,k)] deriving (Eq, Ord)
+
+i0 = O [(0,1)] :: Octonion Q
+i1 = O [(1,1)] :: Octonion Q
+i2 = O [(2,1)] :: Octonion Q
+i3 = O [(3,1)] :: Octonion Q
+i4 = O [(4,1)] :: Octonion Q
+i5 = O [(5,1)] :: Octonion Q
+i6 = O [(6,1)] :: Octonion Q
+
+fromList as = O $ filter ((/=0) . snd) $ zip [-1..6] as
+
+toList (O xs) = toList' xs [-1..6] where
+    toList' ((i,a):xs) (j:js) =
+        if i == j then a : toList' xs js else 0 : toList' ((i,a):xs) js
+    toList' [] (j:js) = 0 : toList' [] js
+    toList' _ [] = []
+
+expose (O ts) = ts
+
+instance Show k => Show (Octonion k) where
+    show (O []) = "0"
+    show (O ts) = let c:cs = concatMap showTerm ts
+                  in if c == '+' then cs else c:cs
+        where showTerm (i,a) = showCoeff a ++ showImag a i
+              showCoeff a = case show a of
+                            "1" -> "+"
+                            "-1" -> "-"
+                            '-':cs -> '-':cs
+                            cs -> '+':cs
+              showImag a i | i == -1 = case show a of
+                                       "1" -> "1"
+                                       "-1" -> "1"
+                                       otherwise -> ""
+                           | otherwise = "i" ++ show i
+
+instance (Ord k, Num k) => Num (Octonion k) where -- Ord k not strictly required, but keeps nf simpler
+    O ts + O us = O $ nf $ ts ++ us
+    negate (O ts) = O $ map (\(i,a) -> (i,-a)) ts
+    O ts * O us = O $ nf [m t u | t <- ts, u <- us]
+    fromInteger 0 = O []
+    fromInteger n = O [(-1, fromInteger n)]
+
+nf ts = nf' $ L.sort ts where
+    nf' ((i1,a1):(i2,a2):ts) =
+        if i1 == i2
+        then if a1+a2 == 0 then nf' ts else nf' ((i1,a1+a2):ts)
+        else (i1,a1) : nf' ((i2,a2):ts)
+    nf' ts = ts
+
+m (-1,a) (i,b) = (i,a*b)
+m (i,a) (-1,b) = (i,a*b)
+m (i,a) (j,b) =
+    case (j-i) `mod` 7 of
+    0 -> (-1,-a*b)
+    1 -> ( (i+3) `mod` 7, a*b)	-- i_n+1 * i_n+2 == i_n+4
+    2 -> ( (i+6) `mod` 7, a*b)	-- i_n+2 * i_n+4 == i_n+1
+    3 -> ( (i+1) `mod` 7, -a*b)	-- i_n+1 * i_n+4 == -i_n+2
+    4 -> ( (i+5) `mod` 7, a*b)	-- i_n+4 * i_n+1 == i_n+2
+    5 -> ( (i+4) `mod` 7, -a*b)	-- i_n+4 * i_n+2 == -i_n+1
+    6 -> ( (i+2) `mod` 7, -a*b)	-- i_n+2 * i_n+1 == -i_n+4
+
+
+conj (O ts) = O $ map (\(i,a) -> if i == -1 then (i,a) else (i,-a)) ts
+
+sqnorm (O ts) = sum [a^2 | (i,a) <- ts]
+
+instance (Ord k, Num k, Fractional k) => Fractional (Octonion k) where
+    recip x = let O x' = conj x
+                  xx' = sqnorm x
+              in O $ map (\(i,a) -> (i,a/xx')) x'
+
+
+isOrthogonal (O ts) (O us) = dot ts us == 0 where
+    dot ((i,a):ts) ((j,b):us) =
+        case compare i j of
+        EQ -> a*b + dot ts us
+        LT -> dot ts ((j,b):us)
+        GT -> dot ((i,a):ts) us
+    dot _ _ = 0
+
+antiCommutes x y = x*y + y*x == 0
+
+-- anti-commuting and being orthogonal appear to be equivalent for unit imaginary octonions,
+-- provided we're not in characteristic 2
+
+
+-- OCTONIONS OVER FINITE FIELDS
+
+{-
+octonions fq = map O $ octonions' [-1..6] where
+    octonions' (i:is) = [if a == 0 then ts else (i,a):ts | a <- fq, ts <- octonions' is]
+    octonions' [] = [[]]
+-}
+octonions fq = map fromList $ ptsAG 8 fq
+
+isUnit x = sqnorm x == 1
+
+unitImagOctonions fq = filter isUnit $ map (fromList . (0:)) $ ptsAG 7 fq
+
+
+-- given the images of i0, i1, i2, return the automorphism
+-- the inputs must be pure imaginary unit octonions
+-- and we must have isOrthogonal i0 i1, isOrthogonal i0 i2, isOrthogonal i1 i2, and isOrthogonal (i0*i1) i2
+autFrom i0' i1' i2' =
+    let 0:r0 = toList i0'
+        0:r1 = toList i1'
+        0:r2 = toList i2'
+        0:r3 = toList $ i0'*i1'
+        0:r4 = toList $ i1'*i2'
+        0:r5 = toList $ i0'*(i1'*i2')
+        0:r6 = toList $ i0'*i2'
+    in [r0,r1,r2,r3,r4,r5,r6]
+
+
+x %^ g =
+    let a:as = toList x
+    in fromList $ a : (as <*>> g)
+
+
+-- G2(3)
+
+alpha3 = autFrom (O [(1,1::F3)]) (O [(2,1)]) (O [(3,1)])
+beta3 = autFrom (O [(0,1::F3)]) (O [(2,1)]) (O [(4,1)])
+
+gamma3s = [x | x <- unitImagOctonions f3, isOrthogonal (O [(0,1)]) x, isOrthogonal (O [(1,1)]) x, isOrthogonal (O [(3,1)]) x]
+
+gamma3 = autFrom (O [(0,1::F3)]) (O [(1,1)]) (O [(2,1),(4,1),(5,1),(6,1)])
+
+alpha3' = fromPairs [(x, x %^ alpha3) | x <- unitImagOctonions f3]
+beta3' = fromPairs [(x, x %^ beta3) | x <- unitImagOctonions f3]
+gamma3' = fromPairs [(x, x %^ gamma3) | x <- unitImagOctonions f3]
+
+-- These three together generate a group of order 4245696, which is therefore the whole of G2(3)
+-- (But takes nearly 10 minutes to construct the BSGS in interpreter)
+
+-- Unit imaginary octonions form one orbit under the action of G2
+
+-- [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
+
+
+-- G2(4)
+
+alpha4 = autFrom (O [(1,1::F4)]) (O [(2,1)]) (O [(3,1)])
+beta4 = autFrom (O [(0,1::F4)]) (O [(2,1)]) (O [(4,1)])
+
+gamma4s = [x | x <- unitImagOctonions f4, isOrthogonal (O [(0,1)]) x, isOrthogonal (O [(1,1)]) x, isOrthogonal (O [(3,1)]) x]
+
+gamma4 = autFrom (O [(0,1::F4)]) (O [(1,1)]) (O [(5,embed x),(6,embed $ 1+x)])
+
+alpha4' = fromPairs [(x, x %^ alpha4) | x <- unitImagOctonions f4]
+beta4' = fromPairs [(x, x %^ beta4) | x <- unitImagOctonions f4]
+gamma4' = fromPairs [(x, x %^ gamma4) | x <- unitImagOctonions f4]
+
+-- Haven't checked whether these generate whole group - can be expected to run a long time
+ Math/Projects/KnotTheory/IwahoriHecke.hs view
@@ -0,0 +1,127 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -XFlexibleInstances #-}
+
+module Math.Projects.KnotTheory.IwahoriHecke where
+
+-- import qualified Data.Map as M
+
+import Math.Algebra.Field.Base
+import Math.Algebra.NonCommutative.NCPoly as NP hiding (x,y,z)
+import Math.Algebra.NonCommutative.GSBasis
+
+import Math.Projects.KnotTheory.LaurentMPoly as LP hiding (z)
+import Math.Projects.KnotTheory.Braid
+
+
+-- IWAHORI-HECKE ALGEBRAS
+
+data IwahoriHeckeGens = T Int deriving (Eq,Ord)
+
+instance Show IwahoriHeckeGens where
+    show (T i) = 't': show i
+
+t_ i = NP [(M [T i], 1)] :: NPoly LPQ IwahoriHeckeGens
+
+t1 = t_ 1
+t2 = t_ 2
+t3 = t_ 3
+t4 = t_ 4
+
+q = LP.var "q" -- :: LPQ
+z = LP.var "z" -- :: LPQ
+
+q' = NP.inject q :: NPoly LPQ IwahoriHeckeGens
+z' = NP.inject z :: NPoly LPQ IwahoriHeckeGens
+
+-- inverses
+instance Invertible (NPoly LPQ IwahoriHeckeGens) where
+    inv (NP [(M [T i], 1)]) = (t_ i - z') / q'
+
+-- x ^- n = inv x ^ n
+
+-- Iwahori-Hecke algebra Hn(q,z), generated by n-1 elts t1..t_n-1, together with relations
+ihRelations n =
+    [t_ i * t_ j - t_ j * t_ i | i <- [1..n-1], j <- [i+2..n-1] ] ++
+    [t_ i * t_ j * t_ i - t_ j * t_ i * t_ j | i <- [1..n-1], j <- [1..n-1], abs (i-j) == 1 ] ++
+    [(t_ i)^2 - z' * t_ i - q' | i <- [1..n-1] ]
+
+-- given an elt of the Temperley-Lieb algebra, return the dimension it's defined over (ie the number of points)
+dimIH (NP ts) = 1 + maximum (0 : [i | (M bs,c) <- ts, T i <- bs])
+
+-- Reduce to normal form
+ihnf f = f %% (gb $ ihRelations $ dimIH f)
+
+-- Monomial basis for Iwahori-Hecke algebra (as quotient of free algebra by Iwahori-Hecke relations)
+ihBasis n = mbasisQA [t_ i | i <- [1..n-1]] (gb $ ihRelations n)
+
+
+-- OCNEANU TRACE
+
+-- Trace function on single terms
+tau' 1 (1,c) = c
+tau' 2 (1,c) = c * (1-q)/z
+tau' 2 (m,c) = c
+tau' n (m,c) = case m `divM` M [T (n-1)] of
+               Just (l,r) -> tau (n-1) $ NP [(l*r,c)] -- have to call tau not tau', because l*r might be reducible
+               Nothing -> tau' (n-1) (m,c*(1-q)/z)
+
+-- Trace function on polynomials, by linearity
+-- Given an elt of Iwahori-Hecke algebra, returns an elt of Q[q,z]
+tau n f | dimIH f <= n = let NP ts = ihnf f in sum [tau' n t | t <- ts]
+
+
+fromBraid f = ihnf (NP.subst skeinRelations f) where
+    skeinRelations = concat [ [(s_ i, t_ i), (s_ (-i), (t_ i - z') / q')] | i <- [1..] ]
+
+
+-- HOMFLY polynomial, also called Jones-Conway polynomial
+-- Kassel, Turaev version, as poly in x,y
+-- Satisfies skein relation x P_L_+(x,y) - x^-1 P_L_-(x,y) == y P_L_0(x,y)
+-- n is the index of the Iwahori-Hecke algebra we're working in (ie the number of strings in the braid)
+-- f is the braid expressed in the Iwahori-Hecke generators ti
+homfly n f = LP.subst [(q,1/x^2),(z,y/x)] $ tau n $ fromBraid f
+
+i = LP.var "i" :: LPQ
+l = LP.var "l" :: LPQ
+m = LP.var "m" :: LPQ
+
+-- HOMFLY polynomial (for braid f over n strings)
+-- Lickorish version, as poly in l,m
+-- Satisfies skein relation l P_L_+(l,m) + l^-1 P_L_-(l,m) + m P_L_0(l,m) == 0
+homfly' n f =
+    let f' = LP.subst [(x,i^3*l),(y,i*m)] (homfly n f)
+    in reduceLP f' (i^2+1)
+
+-- Closer to Thistlethwaite notation
+homfly'' n f = sum $ zipWith (*) (map LP.inject $ coeffs (m^2) (homfly' n f)) (iterate (*(m'^2)) 1)
+    where m' = LP.var "m" :: LaurentMPoly LPQ
+    -- where m' = LP [(LM $ M.singleton "m" 1, 1)] :: LaurentMPoly LPQ
+
+-- express an lpoly as a upoly over one of its variables
+coeffs v 0 = []
+coeffs v f = let (f',c) = quotRemLP f v in c : coeffs v f'
+
+
+-- Jones polynomial (for braid f over m strings)
+-- from the HOMFLY polynomial via the substitution x -> 1/t, y -> t^1/2 - t^-1/2
+-- The reason the code is so complicated is that y == t^1/2 - t^-1/2 can appear in the denominator,
+-- and we have to cancel it out with the numerator by hand because our code can't do it for us
+jones' m f = let f' = homfly m f
+                 n = d*f'
+                 d = denominatorLP f'
+                 subs = [(x,1/t),(y,t^^^(1/2)-1/t^^^(1/2))]
+                 -- subs = [(x,1/t^2),(y,t-1/t)]
+                 n' = LP.subst subs n
+                 d' = LP.subst subs d
+                 nn = nd*n'
+                 nd = denominatorLP n'
+                 dn = dd*d'
+                 dd = denominatorLP d'
+                 (q,r) = quotRemLP nn dn
+             -- in if r == 0 then halfExponents' (dd/nd * q) else error ""
+             in if r == 0 then (dd/nd * q) else error ""
+
+-- Alexander polynomial is the substitution x -> 1, y -> t^1/2 - t^-1/2
+
+
+ Math/Projects/KnotTheory/LaurentMPoly.hs view
@@ -0,0 +1,194 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+module Math.Projects.KnotTheory.LaurentMPoly where
+
+import qualified Data.Map as M
+import Data.List as L
+
+import Math.Algebra.Field.Base
+
+
+-- It would be possible to refactor this and other code into common code for semigroup rings
+-- But there are enough small fiddly differences that it's easier not to
+
+
+-- LAURENT MONOMIALS
+
+newtype LaurentMonomial = LM (M.Map String Q) deriving (Eq)
+-- We allow exponents to be rationals, because we want to support rings such as Z[q^1/2,q^-1/2] in connection with Hecke algebras
+
+-- the multidegree - not sure how meaningful this is if we have negative indices too
+degLM (LM m) = sum $ M.elems m
+
+-- Glex ordering
+instance Ord LaurentMonomial where
+    compare a b = let ds = M.elems m where LM m = a/b in
+                  case compare (sum ds) 0 of
+                  GT -> GT
+                  LT -> LT
+                  EQ -> if null ds then EQ else
+                        if head ds > 0 then GT else LT
+
+instance Show LaurentMonomial where
+    show (LM a) | M.null a = "1"
+                | otherwise = concatMap showVar $ M.toList a
+                where showVar (v,1) = v
+                      showVar (v,i) = v ++ "^" ++ show i
+
+instance Num LaurentMonomial where
+    LM a * LM b = LM $ M.filter (/=0) $ M.unionWith (+) a b
+    fromInteger 1 = LM M.empty
+
+instance Fractional (LaurentMonomial) where
+    recip (LM m) = LM $ M.map negate m
+
+
+-- untested
+-- numeratorLM (LM a) = LM $ M.filter (>0) a
+
+denominatorLM (LM a) = recip $ LM $ M.filter (<0) a
+
+-- not valid for arguments with negative exponents (because 0 won't trump -i)
+lcmLM (LM a) (LM b) = LM $ M.unionWith max a b
+
+-- not tested -- for arguments with non-zero denominators
+-- gcdLM (LM a) (LM b) = LM $ M.intersectionWith min a b
+
+divLM a b = let LM c = a/b in if all (>=0) (M.elems c) then Just (LM c) else Nothing
+
+
+-- LAURENT POLYNOMIALS
+
+newtype LaurentMPoly r = LP [(LaurentMonomial,r)] deriving (Eq,Ord)
+
+instance Show r => Show (LaurentMPoly r) where
+    show (LP []) = "0"
+    show (LP ts) =
+        let (c:cs) = concatMap showTerm (reverse ts) -- we show Laurent polys with smallest terms first
+        in if c == '+' then cs else c:cs
+        where showTerm (m,c) =
+                  case show c of
+                  "1" -> "+" ++ show m
+                  "-1" -> "-" ++ show m
+                  -- cs@(x:_) -> (if x == '-' then cs else '+':cs) ++ (if m == 1 then "" else show m)
+                  cs -> showCoeff cs ++ (if m == 1 then "" else show m)
+              showCoeff (c:cs) = if any (`elem` ['+','-']) cs
+                                 then "+(" ++ c:cs ++ ")"
+                                 else if c == '-' then c:cs else '+':c:cs
+              -- we don't attempt sign reversal within brackets in case we have expressions like t^-1 inside the brackets
+
+instance Num r => Num (LaurentMPoly r) where
+    LP ts + LP us = LP (mergeTerms ts us)
+    negate (LP ts) = LP $ map (\(m,c)->(m,-c)) ts
+    LP ts * LP us = LP $ collect $ sortBy cmpTerm $ [(g*h,c*d) | (g,c) <- ts, (h,d) <- us]
+    fromInteger 0 = LP []
+    fromInteger n = LP [(fromInteger 1, fromInteger n)]
+
+cmpTerm (a,c) (b,d) = case compare a b of EQ -> EQ; GT -> LT; LT -> GT
+-- we have to put largest terms first so that quotRem works
+
+-- inputs in descending order
+mergeTerms (t@(g,c):ts) (u@(h,d):us) =
+    case cmpTerm t u of
+    LT -> t : mergeTerms ts (u:us)
+    GT -> u : mergeTerms (t:ts) us
+    EQ -> if e == 0 then mergeTerms ts us else (g,e) : mergeTerms ts us
+    where e = c + d
+mergeTerms ts us = ts ++ us -- one of them is null
+
+collect (t1@(g,c):t2@(h,d):ts)
+    | g == h = collect $ (g,c+d):ts
+    | c == 0  = collect $ t2:ts
+    | otherwise = t1 : collect (t2:ts)
+collect ts = ts
+
+-- Fractional instance so that we can enter fractional coefficients
+-- Only lets us divide by single terms, not any other polynomials
+instance Fractional r => Fractional (LaurentMPoly r) where
+    recip (LP [(m,c)]) = LP [(recip m, recip c)]
+    recip _ = error "LaurentMPoly.recip: only supported for (non-zero) constants or monomials"
+
+
+lm (LP ((m,c):ts)) = m
+lc (LP ((m,c):ts)) = c
+lt (LP ((m,c):ts)) = LP [(m,c)]
+
+quotRemLP f g
+    | g == 0 = error "quotRemLP: division by zero"
+    | denominatorLP f /= 1 || denominatorLP g /= 1 = error "quotRemLP: negative exponents"
+    | otherwise = quotRemLP' f (0,0)
+    where
+    quotRemLP' 0 (q,r) = (q,r)
+    quotRemLP' h (q,r) =
+        case lm h `divLM` lm g of
+        Just m -> let t = LP [(m, lc h / lc g)]
+                  in quotRemLP' (h-t*g) (q+t,r)
+        Nothing -> let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
+                   in quotRemLP' (h-lth) (q, r+lth)
+
+
+-- g must be a binomial without negative exponents - eg i^2+1
+reduceLP f g@(LP [_,_]) =
+    let fn = f * fd
+        fd = denominatorLP f
+        (_,rn) = quotRemLP fn g
+        (_,rd) = quotRemLP fd g
+    in rn / rd
+
+
+var v = LP [(LM $ M.singleton v 1, 1)]
+
+t = var "t" :: LaurentMPoly Q
+x = var "x" :: LaurentMPoly Q
+y = var "y" :: LaurentMPoly Q
+z = var "z" :: LaurentMPoly Q
+
+
+denominatorLP (LP ts) = LP [(m',1)] where
+    m' = foldl lcmLM 1 [denominatorLM m | (m,c) <- ts]
+
+{-
+-- not tested for terms with non-zero denominator
+gcdTermsLP (LP ts) = LP [(m',1)] where
+    m' = foldl gcdLM 1 [m | (m,c) <- ts]
+-}
+
+-- injection of field elements into polynomial ring
+inject 0 = LP []
+inject c = LP [(fromInteger 1, c)]
+
+sqrtvar v = LP [(LM $ M.singleton v (1/2), 1)]
+
+-- substitute terms for variables in an MPoly
+-- eg subst [(x,a),(y,a+b),(z,c^2)] (x*y+z) -> a*(a+b)+c^2
+subst vts (LP us) = sum [inject c * substM m | (m,c) <- us] where
+    substM (LM m) = product [substV v ^^^ i | (v,i) <- M.toList m]
+    substV v =
+        let v' = var v in
+        case L.lookup v' vts of
+        Just t -> t
+        Nothing -> v' -- no substitute, so keep as is
+
+f ^^^ i | denominatorQ i == 1 = f ^^ numeratorQ i -- exponent is an integer
+        | otherwise = case f of
+                      LP [(LM m,1)] -> LP [(LM $ M.map (*i) m ,1)]-- base is a monomial
+                      otherwise     -> error ("(^^^): Cannot calculate " ++ show f ++ " ^^^ " ++ show i)
+{-
+-- halve all indices - useful when we really want to be working over k[t^1/2,t^-1/2]
+halfExponents (LP ts) =
+    if any odd (concatMap (\(LM m,c) -> M.elems m) ts)
+    then error ("halfExponents: " ++ show (LP ts))
+    else LP $ map (\(LM m, c) -> (LM $ M.filter (/=0) $ M.map (`div` 2) m, c)) ts
+
+halfExponents' f@(LP ts) =
+    let f'@(LP us) = LP $ map (\(LM m, c) -> (LM $ M.map (`div` 2) m, c)) ts
+    in if any (==0) (concatMap (\(LM m,c) -> M.elems m) us)
+       then Left f
+       else Right f'
+
+quarterExponents' f@(LP ts) =
+    let f'@(LP us) = LP $ map (\(LM m, c) -> (LM $ M.map (`div` 4) m, c)) ts
+    in if any (==0) (concatMap (\(LM m,c) -> M.elems m) us)
+       then Left f
+       else Right f'
+-}
+ Math/Projects/KnotTheory/TemperleyLieb.hs view
@@ -0,0 +1,85 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -XFlexibleInstances #-}
+
+module Math.Projects.KnotTheory.TemperleyLieb where
+
+import Data.List ( (\\) )
+
+import Math.Algebra.Field.Base
+import Math.Algebra.NonCommutative.NCPoly as NP
+import Math.Algebra.NonCommutative.GSBasis
+
+import Math.Projects.KnotTheory.LaurentMPoly as LP
+import Math.Projects.KnotTheory.Braid
+
+
+-- TEMPERLEY-LIEB ALGEBRAS
+
+data TemperleyLiebGens = E Int deriving (Eq,Ord)
+
+instance Show TemperleyLiebGens where
+    show (E i) = 'e': show i
+
+e_ i = NP [(M [E i], 1)] :: NPoly LPQ TemperleyLiebGens
+
+-- d is the value of a closed loop
+d = LP.var "d"
+d' = NP.inject d :: NPoly LPQ TemperleyLiebGens
+
+e1 = e_ 1
+e2 = e_ 2
+e3 = e_ 3
+e4 = e_ 4
+
+-- Temperley-Lieb algebra An(d), generated by n-1 elts e1..e_n-1, together with relations
+tlRelations n =
+    [e_ i * e_ j - e_ j * e_ i | i <- [1..n-1], j <- [i+2..n-1] ] ++
+    [e_ i * e_ j * e_ i - e_ i | i <- [1..n-1], j <- [1..n-1], abs (i-j) == 1 ] ++
+    [(e_ i)^2 - d' * e_ i | i <- [1..n-1] ]
+
+-- given an elt of the Temperley-Lieb algebra, return the dimension it's defined over (ie the number of points)
+dimTL (NP ts) = 1 + maximum (0 : [i | (M bs,c) <- ts, E i <- bs])
+
+-- Reduce to normal form
+tlnf f = f %% (gb $ tlRelations $ dimTL f)
+
+-- Monomial basis for Temperley-Lieb algebra (as quotient of free algebra by Temperley-Lieb relations)
+tlBasis n = mbasisQA [e_ i | i <- [1..n-1]] (gb $ tlRelations n)
+
+
+-- trace function
+-- the trace of an elt is d^k, where k is the number of loops in its closure (ie join the top and bottom of the diagram to make an annulus)
+-- this is clearly the same as the number of cycles of the elt when thought of as an elt of Sn, with ei mapped to the transposition (i i+1)
+tr' n (M g) = d ^ ( -1 + length (orbits g [1..n]) ) where
+    image i [] = i
+    image i (E j : es) | i == j    = image (i+1) es
+                       | i == j+1  = image (i-1) es
+                       | otherwise = image i es
+    orbits g [] = []
+    orbits g (i:is) = let i' = orbit i [] in i' : orbits g ((i:is) \\ i')
+    orbit j js = let j' = image j g in if j' `elem` (j:js) then reverse (j:js) else orbit j' (j:js)
+-- Note, some authors define the trace so that tr 1 == 1.
+-- That is the same as this trace except for a factor of d^(n-1)
+
+tr n f@(NP ts) = sum [c * tr' n m | (m,c) <- ts]
+
+
+-- JONES POLYNOMIAL
+
+a = LP.var "a"
+a' = NP.inject a :: NPoly LPQ TemperleyLiebGens
+
+-- Convert a braid to Temperley-Lieb algebra using Skein relation
+fromBraid f = tlnf (NP.subst skeinRelations f) where
+    skeinRelations = concat [ [(s_ i, 1/a' * e_ i + a'), (s_ (-i), a' * e_ i + 1/a')] | i <- [1..] ]
+
+
+-- Jones polynomial
+-- n the number of strings, f the braid
+jones n f = let kauffman = LP.subst [(d, - a^2 - 1/a^2)] $ tr n (fromBraid f)
+                j = (-a)^^(-3 * writhe f) * kauffman
+            -- in halfExponents $ halfExponents $ LP.subst [(a,1/t)] j
+            -- in quarterExponents' $ LP.subst [(a,1/t)] j
+            in LP.subst [(a,1/t^^^(1/4))] j
+
+ Math/Projects/RootSystem.hs view
@@ -0,0 +1,253 @@+-- Copyright (c) David Amos, 2008. All rights reserved.
+
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+module Math.Projects.RootSystem where
+
+import Data.Ratio
+import qualified Data.List as L
+import qualified Data.Set as S
+
+import Math.Algebra.LinearAlgebra
+import Math.Algebra.Group.PermutationGroup hiding (elts, order)
+import Math.Algebra.Group.SchreierSims as SS
+import Math.Algebra.Group.StringRewriting as SG
+
+import Math.Algebra.Field.Base -- for Q
+
+data Type = A | B | C | D | E | F | G
+
+-- Humphreys, Reflection Groups and Coxeter Groups
+
+
+-- SIMPLE SYSTEMS
+-- sometimes called fundamental systems
+
+-- The ith basis vector in K^n
+basisElt :: Int -> Int -> [Q] -- this type signature determines all the rest
+basisElt n i = replicate (i-1) 0 ++ 1 : replicate (n-i) 0
+-- We need to work over the rationals to ensure that arithmetic is exact
+-- So long as our simple systems are rational, then reflection matrices are rational
+
+-- A simple system is like a basis for the root system (see Humphreys p8 for full definition)
+-- simpleSystem :: Type -> Int -> [[Q]]
+simpleSystem A n | n >= 1 = [e i <-> e (i+1) | i <- [1..n]]
+    where e = basisElt (n+1)
+simpleSystem B n | n >= 2 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [e n]
+    where e = basisElt n
+simpleSystem C n | n >= 2 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [2 *> e n]
+    where e = basisElt n
+simpleSystem D n | n >= 4 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [e (n-1) <+> e n]
+    where e = basisElt n
+simpleSystem E n | n `elem` [6,7,8] = take n simpleroots
+    where e = basisElt 8
+          simpleroots = ((1/2) *> (e 1 <-> e 2 <-> e 3 <-> e 4 <-> e 5 <-> e 6 <-> e 7 <+> e 8))
+                      : (e 1 <+> e 2)
+                      : [e (i-1) <-> e (i-2) | i <- [3..8]]
+simpleSystem F 4 = [e 2 <-> e 3, e 3 <-> e 4, e 4, (1/2) *> (e 1 <-> e 2 <-> e 3 <-> e 4)]
+    where e = basisElt 4
+simpleSystem G 2 = [e 1 <-> e 2, ((-2) *> e 1) <+> e 2 <+> e 3]
+    where e = basisElt 3
+
+
+-- ROOT SYSTEMS
+-- Calculating the full root system from the fundamental roots
+
+-- Humphreys p3
+-- Weyl group element corresponding to a root
+-- w r is the reflection in the hyperplane orthogonal to r
+w r s = s <-> (2 * (s <.> r) / (r <.> r)) *> r
+
+-- Given a simple system, return the full root system
+-- The closure of a set of roots under reflection
+closure rs = S.toList $ closure' S.empty (S.fromList rs) where
+    closure' interior boundary
+        | S.null boundary = interior
+        | otherwise =
+            let interior' = S.union interior boundary
+                boundary' = S.fromList [w r s | r <- rs, s <- S.toList boundary] S.\\ interior'
+            in closure' interior' boundary'
+
+
+-- WEYL GROUP
+-- The finite reflection group generated by the root system
+
+-- Generators of the Weyl group as permutation group on the roots
+weylPerms t n =
+    let rs = simpleSystem t n
+        xs = closure rs
+        toPerm r = fromPairs [(x, w r x) | x <- xs]
+    in map toPerm rs
+    
+-- Generators of the Weyl group as a matrix group
+weylMatrices t n = map wMx (simpleSystem t n)
+
+-- The Weyl group element corresponding to a root, represented as a matrix
+wMx r = map (w r) [e i | i <- [1..d]] -- matrix for reflection in hyperplane orthogonal to r
+    where d = length r -- dimension of the space
+          e = basisElt d
+-- the images of the basis elts form the columns of the matrix
+-- however, reflection matrices are symmetric, so they also form the rows
+
+
+-- CARTAN MATRIX, DYNKIN DIAGRAM, COXETER SYSTEM
+
+cartanMatrix t n = [[2 * (ai <.> aj) / (ai <.> ai) | aj <- roots] | ai <- roots]
+    where roots = simpleSystem t n
+-- Note: The Cartan matrices for A, D, E systems are symmetric.
+-- Those of B, C, F, G are not
+-- Carter, Simple Groups of Lie Type, p44-5 gives the expected answers
+-- They agree with our answers except for G2, which is the transpose
+-- (So probably Carter defines the roots of G2 the other way round to Humphreys)
+
+-- set the diagonal entries of (square) matrix mx to constant c
+setDiag c mx@((x:xs):rs) = (c:xs) : zipWith (:) (map head rs) (setDiag c $ map tail rs)
+setDiag _ [[]] = [[]]
+
+-- Carter, Segal, Macdonald p17-18
+-- given a Cartan matrix, derive the corresponding matrix describing the Dynkin diagram
+-- nij = Aij * Aji, nii = 0
+dynkinFromCartan aij = setDiag 0 $ (zipWith . zipWith) (*) aij (L.transpose aij)
+
+dynkinDiagram t n = dynkinFromCartan $ cartanMatrix t n
+
+-- given the Dynkin diagram nij, derive the coefficients mij of the Coxeter group <si | si^2, (sisj)^mij> (so mii == 1)
+-- using nij = 4 cos^2 theta_ij
+-- nij == 0 <=> theta = pi/2
+-- nij == 1 <=> theta = pi/3
+-- nij == 2 <=> theta = pi/4
+-- nij == 3 <=> theta = pi/6
+coxeterFromDynkin nij = setDiag 1 $ (map . map) f nij
+    where f 0 = 2; f 1 = 3; f 2 = 4; f 3 = 6
+
+-- The mij coefficients of the Coxeter group <si | si^2, (sisj)^mij>, as a matrix
+coxeterMatrix t n = coxeterFromDynkin $ dynkinDiagram t n
+
+
+-- Given the matrix of coefficients mij, return the Coxeter group <si | si^2, (sisj)^mij>
+-- We assume but don't check that mii == 1 and mij == mji
+fromCoxeterMatrix mx = (gs,rs) where
+    n = length mx
+    gs = map s_ [1..n]
+    rs = rules mx 1
+    rules [] _ = []
+    rules ((1:xs):rs) i = ([s_ i, s_ i],[]) : [powerRelation i j m | (j,m) <- zip [i+1..] xs] ++ rules (map tail rs) (i+1)
+    powerRelation i j m = (concat $ replicate m [s_ i, s_ j],[])
+
+-- Another presentation for the Coxeter group, using braid relations
+fromCoxeterMatrix2 mx = (gs,rs) where
+    n = length mx
+    gs = map s_ [1..n]
+    rs = rules mx 1
+    rules [] _ = []
+    rules ((1:xs):rs) i = ([s_ i, s_ i],[]) : [braidRelation i j m | (j,m) <- zip [i+1..] xs] ++ rules (map tail rs) (i+1)
+    braidRelation i j m = (take m $ cycle [s_ j, s_ i], take m $ cycle [s_ i, s_ j])
+
+
+
+coxeterPresentation t n = fromCoxeterMatrix $ coxeterMatrix t n
+
+eltsCoxeter t n = SG.elts $ fromCoxeterMatrix2 $ coxeterMatrix t n
+-- it's just slightly faster to use the braid presentation
+
+poincarePoly t n = map length $ L.group $ map length $ eltsCoxeter t n
+
+
+-- LIE ALGEBRAS
+
+elemMx n i j = replicate (i-1) z ++ e j : replicate (n-i) z
+    where z = replicate n 0
+          e = basisElt n
+
+
+lieMult x y = x*y - y*x
+
+-- for gluing matrices together
+(+|+) = zipWith (++) -- glue two matrices together side by side
+(+-+) = (++)         -- glue two matrices together above and below
+
+form D n = (zMx n +|+ idMx n)
+                  +-+
+          (idMx n +|+ zMx n)
+form C n = (2 : replicate (2*n) 0) :
+           (map (0:) (form D n))
+form B n = let id' = (-1) *>> idMx n
+           in (zMx n +|+ idMx n)
+                     +-+
+               (id'  +|+ zMx n)
+
+
+-- TESTING
+-- The expected values of the root system, number of roots, order of Weyl group
+-- for comparison against the calculated values
+
+-- !! Not yet got root systems for E6,7,8, F4
+
+-- Humphreys p41ff
+
+-- The full root system
+-- L.sort (rootSystem t n) == closure (simpleSystem t n)
+-- rootSystem :: Type -> Int -> [[QQ]]
+rootSystem A n | n >= 1 = [e i <-> e j | i <- [1..n+1], j <- [1..n+1], i /= j]
+    where e = basisElt (n+1)
+rootSystem B n | n >= 2 = shortRoots ++ longRoots
+    where e = basisElt n
+          shortRoots = [e i | i <- [1..n]]
+                    ++ [[] <-> e i | i <- [1..n]]
+          longRoots  = [e i <+> e j | i <- [1..n], j <- [i+1..n]]
+                    ++ [e i <-> e j | i <- [1..n], j <- [i+1..n]]
+                    ++ [[] <-> e i <+> e j | i <- [1..n], j <- [i+1..n]]
+                    ++ [[] <-> e i <-> e j | i <- [1..n], j <- [i+1..n]]
+rootSystem C n | n >= 2 = longRoots ++ shortRoots
+    where e = basisElt n
+          longRoots  = [2 *> e i | i <- [1..n]]
+                    ++ [[] <-> (2 *> e i) | i <- [1..n]]
+          shortRoots = [e i <+> e j | i <- [1..n], j <- [i+1..n]]
+                    ++ [e i <-> e j | i <- [1..n], j <- [i+1..n]]
+                    ++ [[] <-> e i <+> e j | i <- [1..n], j <- [i+1..n]]
+                    ++ [[] <-> e i <-> e j | i <- [1..n], j <- [i+1..n]]
+rootSystem D n | n >= 4 =
+    [e i <+> e j | i <- [1..n], j <- [i+1..n]]
+    ++ [e i <-> e j | i <- [1..n], j <- [i+1..n]]
+    ++ [[] <-> e i <+> e j | i <- [1..n], j <- [i+1..n]]
+    ++ [[] <-> e i <-> e j | i <- [1..n], j <- [i+1..n]]
+    where e = basisElt n
+rootSystem G 2 = shortRoots ++ longRoots
+    where e = basisElt 3
+          shortRoots = [e i <-> e j | i <- [1..3], j <- [1..3], i /= j]
+          longRoots = concatMap (\r-> [r,[] <-> r]) [2 *> e i <-> e j <-> e k | i <- [1..3], [j,k] <- [[1..3] L.\\ [i]] ]
+
+
+-- numRoots t n == length (closure $ simpleSystem t n)
+numRoots A n = n*(n+1)
+numRoots B n = 2*n*n
+numRoots C n = 2*n*n
+numRoots D n = 2*n*(n-1)
+numRoots E 6 = 72
+numRoots E 7 = 126
+numRoots E 8 = 240    
+numRoots F 4 = 48
+numRoots G 2 = 12
+
+-- The order of the Weyl group
+-- orderWeyl t n == S.order (weylPerms t n)
+orderWeyl A n = factorial (n+1)
+orderWeyl B n = 2^n * factorial n
+orderWeyl C n = 2^n * factorial n
+orderWeyl D n = 2^(n-1) * factorial n
+orderWeyl E 6 = 2^7 * 3^4 * 5
+orderWeyl E 7 = 2^10 * 3^4 * 5 * 7
+orderWeyl E 8 = 2^14 * 3^5 * 5^2 * 7
+orderWeyl F 4 = 2^7 * 3^2
+orderWeyl G 2 = 12
+
+
+factorial n = product [1..toInteger n]
+
+
+
+test1 = all (\(t,n) -> orderWeyl t n == L.genericLength (eltsCoxeter t n))
+    [(A,3),(A,4),(A,5),(B,3),(B,4),(B,5),(C,3),(C,4),(C,5),(D,4),(D,5),(F,4),(G,2)]
+
+test2 = all (\(t,n) -> orderWeyl t n == SS.order (weylPerms t n))
+    [(A,3),(A,4),(A,5),(B,3),(B,4),(B,5),(C,3),(C,4),(C,5),(D,4),(D,5),(E,6),(F,4),(G,2)]
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple
+main = defaultMain
+ license.txt view
@@ -0,0 +1,10 @@+Copyright (c) 2009, David Amos
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 
+Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 
+Neither the name of Haskell for Maths nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. 
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+