HaskellForMaths 0.3.3 → 0.3.4

raw patch · 20 files changed

+2900/−114 lines, 20 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• HaskellForMaths.cabal +7/−2
• Math/Algebra/Commutative/GBasis.hs +12/−0
• Math/Algebra/LinearAlgebra.hs +10/−1
• Math/Algebras/Commutative.hs +6/−1
• Math/Algebras/Structures.hs +6/−2
• Math/Algebras/TensorAlgebra.hs +0/−1
• Math/Algebras/VectorSpace.hs +7/−2
• Math/Combinatorics/FiniteGeometry.hs +28/−8
• Math/Combinatorics/Graph.hs +32/−10
• Math/Combinatorics/GraphAuts.hs +131/−61
• Math/Combinatorics/Matroid.hs +1146/−0
• Math/Combinatorics/Poset.hs +5/−8
• Math/Combinatorics/StronglyRegularGraph.hs +18/−18
• Math/CommutativeAlgebra/GroebnerBasis.hs +282/−0
• Math/CommutativeAlgebra/Polynomial.hs +368/−0
• Math/Core/Field.hs +466/−0
• Math/Core/Utils.hs +64/−0
• Math/Test/TCombinatorics/TMatroid.hs +148/−0
• Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs +75/−0
• Math/Test/TCore/TField.hs +89/−0

@@ -1,5 +1,5 @@    Name:                HaskellForMaths
-   Version:             0.3.3
+   Version:             0.3.4
    Category:            Math
    Description:         A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended for educational purposes, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra
@@ -31,6 +31,9 @@         Math/Test/TCombinatorics/TPoset.hs
         Math/Test/TCombinatorics/TDigraph.hs
         Math/Test/TCombinatorics/TIncidenceAlgebra.hs
+        Math/Test/TCombinatorics/TMatroid.hs
+        Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs
+        Math/Test/TCore/TField.hs
 
    Library
      Build-Depends:     base >= 2 && < 5, containers, array, random, QuickCheck
@@ -50,8 +53,10 @@         Math.Combinatorics.Graph, Math.Combinatorics.GraphAuts, Math.Combinatorics.StronglyRegularGraph,
         Math.Combinatorics.Design, Math.Combinatorics.FiniteGeometry, Math.Combinatorics.Hypergraph,
         Math.Combinatorics.LatinSquares, Math.Combinatorics.Poset, Math.Combinatorics.IncidenceAlgebra,
-        Math.Combinatorics.Digraph,
+        Math.Combinatorics.Digraph, Math.Combinatorics.Matroid,
         Math.Common.IntegerAsType, Math.Common.ListSet,
+        Math.CommutativeAlgebra.Polynomial, Math.CommutativeAlgebra.GroebnerBasis,
+        Math.Core.Utils, Math.Core.Field,
         Math.Projects.RootSystem,
         Math.Projects.Rubik, Math.Projects.MiniquaternionGeometry,
         Math.Projects.ChevalleyGroup.Classical, Math.Projects.ChevalleyGroup.Exceptional,

@@ -269,3 +269,15 @@     else u : merge (t:ts) us
 merge ts us = ts ++ us -- one of them is null
 -}
+
+-- OPERATIONS ON IDEALS
+
+-- Cox et al, p181
+-- Geometric interpretation: V(I+J) = V(I) `intersect` V(J)
+sumI fs gs = gb $ fs ++ gs
+
+-- Cox et al, p183
+-- Geometric interpretation: V(I.J) = V(I) `union` V(J)
+productI fs gs = gb [f * g | f <- fs, g <- gs]
+
+

@@ -61,7 +61,7 @@ (<<*>>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
 a <<*>> b = [ [u <.> v | v <- L.transpose b] | u <- a]
  
--- |k *\> m returns the product k*m of the scalar k and the matrix m
+-- |k *\>\> m returns the product k*m of the scalar k and the matrix m
 (*>>) :: (Num a) => a -> [[a]] -> [[a]]
 k *>> m = (map . map) (k*) m
 
@@ -186,6 +186,15 @@     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
+
+isZero v = all (==0) v
+
+-- inSpanRE m v returns whether the vector v is in the span of the matrix m, where m is required to be in row echelon form
+inSpanRE ((1:xs):bs) (y:ys) = inSpanRE (map tail bs) (if y == 0 then ys else ys <-> y *> xs)
+inSpanRE ((0:xs):bs) (y:ys) = if y == 0 then inSpanRE (xs : map tail bs) ys else False
+inSpanRE _ ys = isZero ys
+
+rank m = length $ filter (not . isZero) $ rowEchelonForm m
 
 -- kernel of a matrix
 -- returns basis for vectors v s.t m <<*> v == 0

@@ -63,7 +63,11 @@ type GlexPoly k v = Vect k (GlexMonomial v)  -+-- |glexVar creates a variable in the algebra of commutative polynomials with Glex term ordering.+-- For example, the following code creates variables called x, y and z:+--+-- > [x,y,z] = map glexVar ["x","y","z"] :: GlexPoly Q String+glexVar :: (Num k) => v -> GlexPoly k v glexVar v = V [(Glex 1 [(v,1)], 1)]  @@ -115,6 +119,7 @@  infixl 7 %% +-- |(%%) reduces a polynomial with respect to a list of polynomials. (%%) :: (Fractional k, Ord b, Show b, Algebra k b, DivisionBasis b)      => Vect k b -> [Vect k b] -> Vect k b f %% gs = r where (_,r) = quotRemMP f gs

@@ -70,13 +70,17 @@     unit = wrap     -- unit 0 = zero -- V []     -- unit x = V [( (),x)]-    mult (V [( ((),()), x)]) = V [( (),x)]+    mult = linear mult' where mult' ((),()) = return ()+    -- mult (V [( ((),()), x)]) = V [( (),x)]+    -- mult (V []) = zerov  instance Num k => Coalgebra k () where     counit = unwrap     -- counit (V []) = 0     -- counit (V [( (),x)]) = x-    comult (V [( (),x)]) = V [( ((),()), x)]+    comult = linear comult' where comult' () = return ((),())+    -- comult (V [( (),x)]) = V [( ((),()), x)]+    -- comult (V []) = zerov  unit' :: (Num k, Algebra k b) => Trivial k -> Vect k b unit' = unit . unwrap -- where unwrap = counit :: Num k => Trivial k -> k

@@ -22,7 +22,6 @@ -- to the infinite direct sum: -- -- T(V) = k ⊕ V ⊕ V⊗V ⊕ V⊗V⊗V ⊕ ...--- where 1 is the unit vector space k =~ Vect k () data TensorAlgebra a = TA Int [a] deriving (Eq,Ord)  instance Show a => Show (TensorAlgebra a) where

@@ -84,7 +84,7 @@ smultL 0 _ = zero -- V [] smultL k (V ts) = V [(ei,k*xi) | (ei,xi) <- ts] --- |Same as smultL. Mnemonic is "multiply through (from the left)"+-- |Same as smultL. Mnemonic is \"multiply through (from the left)\" (*>) :: (Num k) => k -> Vect k b -> Vect k b (*>) = smultL @@ -93,7 +93,7 @@ smultR _ 0 = zero -- V [] smultR (V ts) k = V [(ei,xi*k) | (ei,xi) <- ts] --- |Same as smultR. Mnemonic is "multiply through (from the right)"+-- |Same as smultR. Mnemonic is \"multiply through (from the right)\" (<*) :: (Num k) => Vect k b -> k -> Vect k b (<*) = smultR @@ -130,6 +130,11 @@     V ts >>= f = V $ concat [ [(b,y*x) | let V us = f a, (b,y) <- us] | (a,x) <- ts]     -- Note that as we can't assume Ord a in the Monad instance, we need to call "nf" afterwards +-- |A linear map between vector spaces A and B can be defined by giving its action on the basis elements of A.+-- The action on all elements of A then follows by linearity.+--+-- If we have A = Vect k a, B = Vect k b, and f :: a -> Vect k b is a function from the basis elements of A into B,+-- then @linear f@ is the linear map that this defines by linearity. linear :: (Ord b, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b linear f v = nf $ v >>= f 

@@ -18,13 +18,16 @@ import Math.Algebra.Group.PermutationGroup -- for use in GHCi
 import Math.Algebra.Group.SchreierSims as SS -- for use in GHCi
 
+-- !! The following two functions previously required (FiniteField a) as context
+-- but this has been temporarily removed to enable them to work with Math.Core.Field
+
 -- |ptsAG n fq returns the points of the affine geometry AG(n,Fq), where fq are the elements of Fq
-ptsAG :: (FiniteField a) => Int -> [a] -> [[a]]
+ptsAG :: Int -> [a] -> [[a]]
 ptsAG 0 fq = [[]]
 ptsAG n fq = [x:xs | x <- fq, xs <- ptsAG (n-1) fq]
 
 -- |ptsPG n fq returns the points of the projective geometry PG(n,Fq), where fq are the elements of Fq
-ptsPG :: (FiniteField a) => Int -> [a] -> [[a]]
+ptsPG :: Num a => Int -> [a] -> [[a]]
 ptsPG 0 _ = [[1]]
 ptsPG n fq = map (0:) (ptsPG (n-1) fq) ++ map (1:) (ptsAG n fq)
 
@@ -42,6 +45,9 @@ -- given p1, .., pk, we're looking for all a1 p1 + ... + ak pk, s.t. a1 + ... + ak = 1
 -- if m is the matrix with p1, .., pk as rows, and vs are the vectors [a1, .., ak]
 -- then this is the same as [v <*>> m | v <- vs] == [m' <<*> v | v <- vs]
+
+-- |Given a list of points in AG(n,Fq), return their closure, the smallest flat containing them
+closureAG :: (Ord a, FiniteField a) => [[a]] -> [[a]]
 closureAG ps =
     let vs = [ (1 - sum xs) : xs | xs <- ptsAG (k-1) fq ] -- k-vectors over fq whose sum is 1
     in toListSet [m' <<*> v | v <- vs]
@@ -65,6 +71,9 @@ -- closure of points in PG(n,Fq)
 -- take all linear combinations of the points (ie the subspace generated by the points, considered as points in Fq ^(n+1) )
 -- then discard all which aren't in PNF (thus dropping back into PG(n,Fq))
+
+-- |Given a set of points in PG(n,Fq), return their closure, the smallest flat containing them
+closurePG :: (Ord a, FiniteField a) => [[a]] -> [[a]]
 closurePG ps = toListSet $ filter ispnf $ map (<*>> ps) $ ptsAG k fq where
     k = length ps
     fq = eltsFq undefined
@@ -119,8 +128,13 @@     starColumn r = replicate (r-1) Star ++ replicate (k+1-r) Zero
 
 
--- |flatsPG n fq k returns the k-flats in PG(n,Fq), where fq are the elements of Fq
-flatsPG :: (FiniteField a) => Int -> [a] -> Int -> [[[a]]]
+-- flatsPG :: (FiniteField a) => Int -> [a] -> Int -> [[[a]]]
+
+-- |@flatsPG n fq k@ returns the k-flats in PG(n,Fq), where fq are the elements of Fq.
+-- The returned flats are represented as matrices in reduced row echelon form,
+-- the rows of which are the points that generate the flat.
+-- The full set of points in the flat can be recovered by calling 'closurePG'
+flatsPG :: (Num a) => Int -> [a] -> Int -> [[[a]]]
 flatsPG n fq k = concatMap substStars $ rrefs (n+1) (k+1) where
     substStars (r:rs) = [r':rs' | r' <- substStars' r, rs' <- substStars rs]
     substStars [] = [[]]
@@ -131,20 +145,26 @@ 
 
 -- Flats in AG(n,Fq) are just the flats in PG(n,Fq) which are not "at infinity"
--- |flatsAG n fq k returns the k-flats in AG(n,Fq), where fq are the elements of Fq
-flatsAG :: (FiniteField a) => Int -> [a] -> Int -> [[[a]]]
+-- flatsAG :: (FiniteField a) => Int -> [a] -> Int -> [[[a]]]
+
+-- |flatsAG n fq k returns the k-flats in AG(n,Fq), where fq are the elements of Fq.
+flatsAG :: (Num a) => Int -> [a] -> Int -> [[[a]]]
 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)
 
 
+-- linesPG :: (FiniteField a) => Int -> [a] -> [[[a]]]
+
 -- |The lines (1-flats) in PG(n,fq)
-linesPG :: (FiniteField a) => Int -> [a] -> [[[a]]]
+linesPG :: (Num a) => Int -> [a] -> [[[a]]]
 linesPG n fq = flatsPG n fq 1
 
+-- linesAG :: (FiniteField a) => Int -> [a] -> [[[a]]]
+
 -- |The lines (1-flats) in AG(n,fq)
-linesAG :: (FiniteField a) => Int -> [a] -> [[[a]]]
+linesAG :: (Num a) => Int -> [a] -> [[[a]]]
 linesAG n fq = flatsAG n fq 1
 
 

@@ -11,7 +11,7 @@ import qualified Data.Set as S
 import Control.Arrow ( (&&&) )
 
-import Math.Common.ListSet
+import Math.Common.ListSet as LS
 import Math.Algebra.Group.PermutationGroup hiding (fromDigits, fromBinary)
 import Math.Algebra.Group.SchreierSims as SS
 
@@ -63,6 +63,7 @@ -- graph (vs,es) checks that vs and es are valid before returning the graph.
 graph :: (Ord t) => ([t], [[t]]) -> Graph t
 graph (vs,es) | isGraph vs es = G vs es
+--              | otherwise = error ( "graph " ++ show (vs,es) )
 -- isValid g = g where g = G vs es
 
 toGraph (vs,es) | isGraph vs' es' = G vs' es' where
@@ -123,26 +124,29 @@     es = [[i,j] | i <- [1..n-1], j <- [i+1..n]] -- == combinationsOf 2 [1..n]
 -- automorphism group is Sn
 
--- The complete bipartite graph on m and n vertices
--- kb :: (Integral t) => t -> t -> Graph t
+-- |kb m n is the complete bipartite graph on m and n vertices
+kb :: (Integral t) => t -> t -> Graph t
 kb m n = to1n $ kb' m n
 
--- The complete bipartite graph on m left and n right vertices
--- kb :: (Integral t) => t -> t -> Graph (Either t t)
+-- |kb' m n is the complete bipartite graph on m left and n right vertices
+kb' :: (Integral t) => t -> t -> Graph (Either t t)
 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)
 
+-- |q k is the graph of the k-cube
+q :: (Integral t) => Int -> Graph t
+q k = fromBinary $ q' k
+
+q' :: (Integral t) => Int -> Graph [t]
 q' k = graph (vs,es) where
     vs = sequence $ replicate k [0,1] -- ptsAn k f2
     es = [ [u,v] | [u,v] <- combinationsOf 2 vs, hammingDistance u v == 1 ]
     hammingDistance as bs = length $ filter id $ zipWith (/=) as bs
 -- can probably type-coerce this to be Graph [F2] if required
 
-q k = fromBinary $ q' k
 
-
 tetrahedron = k 4
 
 cube = q 3
@@ -197,6 +201,7 @@     es' = (map . map) fromBinary' es
 -}
 
+petersen :: Graph [Integer]
 petersen = graph (vs,es) where
     vs = combinationsOf 2 [1..5]
     es = [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, disjoint v1 v2]
@@ -207,9 +212,16 @@ 
 -- NEW GRAPHS FROM OLD
 
-complement (G vs es) = graph (vs,es') where es' = combinationsOf 2 vs \\ es
+complement :: (Ord t) => Graph t -> Graph t
+complement (G vs es) = graph (vs,es') where es' = combinationsOf 2 vs LS.\\ es
 -- es' = [e | e <- combinationsOf 2 vs, e `notElem` es]
 
+-- |The restriction of a graph to a subset of the vertices
+restriction :: (Eq a) => Graph a -> [a] -> Graph a
+restriction g@(G vs es) us = G us (es `restrict` us)
+    where es `restrict` us = [e | e@[i,j] <- es, i `elem` us, j `elem` us]
+
+inducedSubgraph :: (Eq a) => Graph a -> [a] -> Graph a
 inducedSubgraph g@(G vs es) us = G us (es `restrict` us)
     where es `restrict` us = [e | e@[i,j] <- es, i `elem` us, j `elem` us]
 
@@ -236,8 +248,12 @@     [(v,_)] -> Just v
     _       -> Nothing
 
+-- |A graph is regular if all vertices have the same valency (degree)
+isRegular :: (Eq t) => Graph t -> Bool
 isRegular g = isJust $ regularParam g
 
+-- |A 3-regular graph is called a cubic graph
+isCubic :: (Eq t) => Graph t -> Bool
 isCubic g = regularParam g == Just 3
 
 
@@ -255,7 +271,8 @@         | otherwise = bfs (nodes ++ [(w:z:zs) | w <- nbrs g z, w `notElem` zs])
     bfs [] = []
 
--- length of the shortest path from x to y
+-- |Within a graph G, the distance d(u,v) between vertices u, v is length of the shortest path from u to v
+distance :: (Eq a) => Graph a -> a -> a -> Int
 distance g x y =
     case findPaths g x y of
     [] -> -1 -- infinite
@@ -296,11 +313,16 @@                       in S.toList boundary : distancePartition' interior' boundary'
 
 -- the connected component to which v belongs
-component g v = concat $ distancePartition g v
+component g v = L.sort $ concat $ distancePartition g v
 
+-- |Is the graph connected?
+isConnected :: (Ord t) => Graph t -> Bool
 isConnected g@(G (v:vs) es) = length (component g v) == length (v:vs)
 isConnected (G [] []) = True
 
+components g = components' (vertices g)
+    where components' [] = []
+          components' (v:vs) = let c = component g v in c : components' (vs LS.\\ c)
 
 -- MORE GRAPHS
 

@@ -20,10 +20,14 @@ 
 -- TRANSITIVITY PROPERTIES OF GRAPHS
 
+-- |A graph is vertex-transitive if its automorphism group acts transitively on the vertices. Thus, given any two distinct vertices, there is an automorphism mapping one to the other.
+isVertexTransitive :: (Ord t) => Graph t -> Bool
 isVertexTransitive (G [] []) = True -- null graph is trivially vertex transitive
 isVertexTransitive g@(G (v:vs) es) = orbitV auts v == v:vs where
     auts = graphAuts g
 
+-- |A graph is edge-transitive if its automorphism group acts transitively on the edges. Thus, given any two distinct edges, there is an automorphism mapping one to the other.
+isEdgeTransitive :: (Ord t) => Graph t -> Bool
 isEdgeTransitive (G _ []) = True
 isEdgeTransitive g@(G vs (e:es)) = orbitE auts e == e:es where
     auts = graphAuts g
@@ -32,6 +36,8 @@ -- unlike edges/blocks, arcs are directed, so the action on them does not sort
 
 -- Godsil & Royle 59-60
+-- |A graph is arc-transitive (or flag-transitive) if its automorphism group acts transitively on arcs. (An arc is an ordered pair of adjacent vertices.)
+isArcTransitive :: (Ord t) => Graph t -> Bool
 isArcTransitive (G _ []) = True -- empty graphs are trivially arc transitive
 isArcTransitive g@(G vs es) = orbit (->^) a auts == a:as where
 -- isArcTransitive g@(G vs es) = closure [a] [ ->^ h | h <- auts] == a:as where
@@ -64,6 +70,8 @@ 
 -- note that a graph with triangles can't be 3-arc transitive, etc, because an aut can't map a self-crossing arc to a non-self-crossing arc
 
+-- |A graph is n-arc-transitive is its automorphism group is transitive on n-arcs. (An n-arc is an ordered sequence (v0,...,vn) of adjacent vertices, with crossings allowed but not doubling back.)
+isnArcTransitive :: (Ord t) => Int -> Graph t -> Bool
 isnArcTransitive _ (G [] []) = True
 isnArcTransitive n g@(G (v:vs) es) =
     orbitP auts v == v:vs && -- isVertexTransitive g
@@ -73,11 +81,16 @@           stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
           a:as = findArcs g v n
 
+is2ArcTransitive :: (Ord t) => Graph t -> Bool
 is2ArcTransitive g = isnArcTransitive 2 g
 
+is3ArcTransitive :: (Ord t) => Graph t -> Bool
 is3ArcTransitive g = isnArcTransitive 3 g
 
 -- Godsil & Royle 66-7
+-- |A graph is distance transitive if given any two ordered pairs of vertices (u,u') and (v,v') with d(u,u') == d(v,v'),
+-- there is an automorphism of the graph that takes (u,u') to (v,v')
+isDistanceTransitive :: (Ord t) => Graph t -> Bool
 isDistanceTransitive (G [] []) = True
 isDistanceTransitive g@(G (v:vs) es)
     | isConnected g =
@@ -147,6 +160,7 @@           es' = S.fromList es
 
 -- Now using distance partitions
+-- Note that because of the use of distance partitions, this is only valid for connected graphs
 graphAuts3 g@(G vs es) = graphAuts' [] [vs] where
     graphAuts' us ((x:ys):pt) =
         let px = refine' (ys : pt) (dps M.! x)
@@ -267,26 +281,32 @@        -- else error (show (src_split, trg_split)) -- for debugging
 
 -- Now, every time we intersect two partitions, refine to an equitable partition
--- |Given a graph g, graphAuts g returns a strong generating set for the automorphism group of g.
+-- |Given a graph g, @graphAuts g@ returns a strong generating set for the automorphism group of g.
+--
+-- Note that the implementation is currently only valid for connected graphs
 graphAuts :: (Ord a) => Graph a -> [Permutation a]
-graphAuts g@(G vs es) = graphAuts' [] (toEquitable g $ valencyPartition g) where
-    graphAuts' us p@((x:ys):pt) =
-        let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
-        in level us p x ys []
-        ++ graphAuts' (x:us) p'
-    graphAuts' us ([]:pt) = graphAuts' us pt
-    graphAuts' _ [] = []
-    level us p@(ph:pt) x (y:ys) hs =
-        let px = refine' (L.delete x ph : pt) (dps M.! x)
-            py = refine' (L.delete y ph : pt) (dps M.! y)
-            uus = zip us us
-        in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
-           []  -> level us p x ys hs
-           h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
-    level _ _ _ [] _ = []
-    dps = M.fromList [(v, distancePartition g v) | v <- vs]
-    es' = S.fromList es
-    nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
+graphAuts g@(G vs es)
+    | isConnected g = graphAuts' [] (toEquitable g $ valencyPartition g)
+    | otherwise = error "graphAuts: only implemented for connected graphs"
+    where graphAuts' us p@((x:ys):pt) =
+              let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
+              in level us p x ys []
+              ++ graphAuts' (x:us) p'
+          graphAuts' us ([]:pt) = graphAuts' us pt
+          graphAuts' _ [] = []
+          level us p@(ph:pt) x (y:ys) hs =
+              let px = refine' (L.delete x ph : pt) (dps M.! x)
+                  py = refine' (L.delete y ph : pt) (dps M.! y)
+                  uus = zip us us
+              in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
+                 []  -> level us p x ys hs
+                 h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
+          level _ _ _ [] _ = []
+          dps = M.fromList [(v, distancePartition g v) | v <- vs]
+          es' = S.fromList es
+          nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
+-- To handle disconnected graphs, you not only need to find auts of each component,
+-- you also need to find auts that swap components
 
 dfsEquitable (dps,es',nbrs_g) xys p1 p2 = dfs xys p1 p2 where
     dfs xys p1 p2
@@ -313,60 +333,110 @@ -- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH
 
 -- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
-incidenceAuts g@(G vs es) = map points (incidenceAuts' [] [vs]) where
-    points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
-    incidenceAuts' us p@((x@(Left _):ys):pt) =
-        -- let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
-        let p' = L.sort $ refine (ys:pt) (dps M.! x)
-        in level us p x ys []
-        ++ incidenceAuts' (x:us) p'
-    incidenceAuts' us ([]:pt) = incidenceAuts' us pt
-    incidenceAuts' _ (((Right _):_):_) = [] -- if we fix all the points, then the blocks must be fixed too 
-    -- incidenceAuts' _ [] = []
-    level us p@(ph:pt) x (y@(Left _):ys) hs =
-        let px = refine' (L.delete x ph : pt) (dps M.! x)
-            py = refine' (L.delete y ph : pt) (dps M.! y)
-            uus = zip us us
-        in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
-           []  -> level us p x ys hs
-           h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
-    level _ _ _ _ _ = [] -- includes the case where y matches Right _, which can only occur on first level, before we've distance partitioned
-    dps = M.fromList [(v, distancePartition g v) | v <- vs]
-    es' = S.fromList es
-    nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
+-- |Given the incidence graph of an incidence structure between points and blocks
+-- (for example, a set system),
+-- @incidenceAuts g@ returns a strong generating set for the automorphism group of the incidence structure.
+-- The generators are represented as permutations of the points.
+-- The incidence graph should be represented with the points on the left and the blocks on the right.
+--
+-- Note that the implementation is currently only valid for connected incidence graphs
+incidenceAuts :: (Ord p, Ord b) => Graph (Either p b) -> [Permutation p]
+incidenceAuts g@(G vs es) 
+    | isConnected g = map points (incidenceAuts' [] [vs])
+    | otherwise = error "incidenceAuts: only implemented for connected incidence graphs"
+    where points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
+          incidenceAuts' us p@((x@(Left _):ys):pt) =
+              -- let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
+              let p' = L.sort $ refine (ys:pt) (dps M.! x)
+              in level us p x ys []
+              ++ incidenceAuts' (x:us) p'
+          incidenceAuts' us ([]:pt) = incidenceAuts' us pt
+          incidenceAuts' _ (((Right _):_):_) = [] -- if we fix all the points, then the blocks must be fixed too 
+          -- incidenceAuts' _ [] = []
+          level us p@(ph:pt) x (y@(Left _):ys) hs =
+              let px = refine' (L.delete x ph : pt) (dps M.! x)
+                  py = refine' (L.delete y ph : pt) (dps M.! y)
+                  uus = zip us us
+              in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
+                 []  -> level us p x ys hs
+                 h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
+          level _ _ _ _ _ = [] -- includes the case where y matches Right _, which can only occur on first level, before we've distance partitioned
+          dps = M.fromList [(v, distancePartition g v) | v <- vs]
+          es' = S.fromList es
+          nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
 
 
 -- GRAPH ISOMORPHISMS
 
 -- !! not yet using equitable partitions, so could probably be more efficient
 
-graphIsos g1 g2 = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2) | v2 <- vertices g2] where
-    v1 = head $ vertices g1
-    dfs xys p1 p2
-        | map length p1 /= map length p2 = []
-        | otherwise =
-            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
+-- graphIsos :: (Ord a, Ord b) => Graph a -> Graph b -> [[(a,b)]]
+graphIsos g1 g2 
+    | isConnected g1 && isConnected g2
+        = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2) | v2 <- vertices g2]
+    | otherwise = error "graphIsos: only implemented for connected graphs"
+    where v1 = head $ vertices g1
+          dfs xys p1 p2
+              | 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
 
+
 -- !! If we're only interested in seeing whether or not two graphs are iso,
 -- !! then the cost of calculating distancePartitions may not be warranted
 -- !! (see Math.Combinatorics.Poset: orderIsos01 versus orderIsos)
+isGraphIso g1 g2 = (not . null) (graphIsos g1 g2)
+
+-- !! deprecate
 isIso g1 g2 = (not . null) (graphIsos g1 g2)
 
+
+-- the following differs from graphIsos in only two ways
+-- we avoid Left, Right crossover isos, by insisting that a Left is taken to a Left (first two lines)
+-- we return only the action on the Lefts, and unLeft it
+-- incidenceIsos :: (Ord p1, Ord b1, Ord p2, Ord b2) =>
+--     Graph (Either p1 b1) -> Graph (Either p2 b2) -> [[(p1,p2)]]
+incidenceIsos g1 g2
+    | isConnected g1 && isConnected g2
+        = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2) | v2@(Left _) <- vertices g2]
+    | otherwise = error "incidenceIsos: only implemented for connected graphs"
+    where v1@(Left _) = head $ vertices g1
+          dfs xys p1 p2
+              | 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 [[(x,y) | (Left x, Left y) <- 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
+
+isIncidenceIso g1 g2 = (not . null) (incidenceIsos g1 g2)
 
 {-
 removeGens x gs = removeGens' [] gs where

@@ -0,0 +1,1146 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++{-# LANGUAGE NoMonomorphismRestriction, DeriveFunctor #-}++-- |A module providing functions to construct and investigate (small, finite) matroids.+module Math.Combinatorics.Matroid where++-- Source: Oxley, Matroid Theory (second edition)++import Math.Core.Utils+import Math.Core.Field hiding (f7)++import Math.Common.ListSet as LS -- set operations on strictly ascending lists+-- import Math.Algebra.Field.Base hiding (Q, F2, F3, F5, F7, F11, f2, f3, f5, f7, f11)+import Math.Algebra.LinearAlgebra hiding (rank, (*>))++import qualified Math.Combinatorics.Graph as G -- hiding (combinationsOf, restriction, component, isConnected)+import Math.Combinatorics.FiniteGeometry+import Math.Combinatorics.GraphAuts+import Math.Algebra.Group.PermutationGroup hiding (closure)++import Math.Algebras.VectorSpace hiding (dual)+import Math.Algebras.Structures+import Math.Algebras.Commutative++import Data.List as L+import qualified Data.Map as M+import qualified Data.Set as S++implies p q = q || not p+exists = not . null+unique [x] = x++shortlex xs ys = case compare (length xs) (length ys) of+    LT -> LT+    EQ -> compare xs ys+    GT -> GT++isShortlex xs = foldcmpl (\x1 x2 -> shortlex x1 x2 /= GT) xs++toShortlex xs = map snd $ L.sort [(length x, x) | x <- xs]++isClutter ss = and [ (s1 `LS.isSubset` s2) `implies` (s1 == s2) | s1 <- ss, s2 <- ss ]+-- note that this definition would allow duplicates++-- single-element deletions of xs+-- if xs is sorted, then so is reverse (deletions xs)+deletions xs = zipWith (++) (inits xs) (tail $ tails xs)++closedUnderSubsets xss = and [xs' `S.member` xss' | xs <- xss, xs' <- deletions xs]+    where xss' = S.fromList xss+++-- |The data structure that we use to store the bases of the matroid+data TrieSet a = TS [(a, TrieSet a)] deriving (Eq,Ord,Functor)+-- Note that in a trie we would normally have a Bool at each node saying whether or not the node is terminal+-- (ie corresponds to a member and not just a prefix of a member).+-- However, since we intend to use the trie to store the bases, and they are all the same length,+-- we can use lack of children to detect that a node is terminal.++-- We could try storing the children in a map rather than a list++-- for debugging+tsshow (TS xts) = "TS [" ++ concatMap (\(x,t) -> "(" ++ show x ++ "," ++ tsshow t ++ ")") xts ++ "]"++instance Show a => Show (TrieSet a) where+    show = show . tstolist++tsempty = TS []++tsinsert (x:xs) (TS ts) =+    case L.lookup x ts of+    Nothing -> let t = tsinsert xs (TS [])+               in TS $ L.insert (x,t) ts+    Just t -> let t' = tsinsert xs t+              in TS $ L.insert (x,t') $ L.delete (x,t) ts+tsinsert [] t = t++tsmember (x:xs) (TS ts) =+    case lookup x ts of+    Nothing -> False+    Just t -> tsmember xs t+tsmember [] (TS []) = True -- the node has no children, and hence is terminal+tsmember [] _ = False++-- xs is a subset of a member of t+tssubmember (x:xs) (TS ts) = or [ case compare x y of+                                  LT -> False+                                  EQ -> tssubmember xs t+                                  GT -> tssubmember (x:xs) t+                                | (y,t) <- ts ]+tssubmember [] _ = True++tstolist (TS []) = [[]]+tstolist (TS xts) = concatMap (\(x,t) -> map (x:) (tstolist t)) xts ++tsfromlist = foldl' (flip tsinsert) tsempty+++-- |A datatype to represent a matroid. @M es bs@ is the matroid whose elements are @es@, and whose bases are @bs@.+-- The normal form is for the @es@ to be in order, for each of the @bs@ individually to be in order.+-- (So the TrieSet should have the property that any path from the root to a leaf is strictly increasing).+data Matroid a = M [a] (TrieSet a) deriving (Eq,Show,Functor)++-- |Return the elements over which the matroid is defined.+elements :: Matroid t -> [t]+elements (M es bs) = es++-- |Return all the independent sets of a matroid, in shortlex order.+indeps :: (Ord a) => Matroid a -> [[a]]+indeps m = bfs [ ([],es) ]+    where es = elements m+          bfs ( (ls,rs) : nodes ) =+              let ls' = reverse ls in+              if isIndependent m ls'+              then ls' : bfs ( nodes ++ successors (ls,rs) )+              else bfs nodes+          bfs [] = []+          successors (ls,rs) = [ (r:ls, rs') | r:rs' <- L.tails rs ]++isIndependent :: (Ord a) => Matroid a -> [a] -> Bool+isIndependent (M es bs) xs = xs `tssubmember` bs++isDependent :: (Ord a) => Matroid a -> [a] -> Bool+isDependent m = not . isIndependent m++-- |Are these the independent sets of a matroid? (The sets must individually be ordered.)+isMatroidIndeps :: (Ord a) => [[a]] -> Bool+isMatroidIndeps is =+    [] `elem` is &&+    closedUnderSubsets is &&+    and [ (l1 < l2) `implies` exists [e | e <- i2 LS.\\ i1, L.insert e i1 `elem` is]+        | i1 <- is, let l1 = length i1, i2 <- is, let l2 = length i2 ]++-- |Construct a matroid from its elements and its independent sets.+fromIndeps :: (Ord a) => [a] -> [[a]] -> Matroid a+fromIndeps es is = fromBases es bs+    where bs = dfs [] [([],es)]+          dfs bs ( node@(ls,rs) : nodes ) =+              let succs = successors node+              in if null succs then dfs (ls:bs) nodes else dfs bs (succs ++ nodes)+          dfs ls [] = let r = length $ last ls -- we know that the first one we found is a true base+                      in map reverse $ filter (\b -> length b == r) ls+                      -- we might have had null succs simply because we ran out of vectors to add+          successors (ls,rs) = [ (r:ls, rs') | r:rs' <- L.tails rs, (r:ls) `S.member` is' ]+          is' = S.fromList $ map reverse is++-- this seems to be just slightly slower+fromIndeps1 es is = fromBases es bs+    where b = greedy [] es -- first find any basis+          greedy ls (r:rs) = if (r:ls) `S.member` ris'+                             then greedy (r:ls) rs+                             else greedy ls rs+          greedy ls [] = reverse ls+          ris' = S.fromList $ map reverse is+          bs = closure S.empty (S.singleton b) -- now find all other bases by passing to "neighbouring" bases+          closure interior boundary =+              if S.null boundary+              then S.toList interior+              else let interior' = interior `S.union` boundary+                       boundary' = S.fromList [ b' | b <- S.toList boundary,+                                                     x <- b, y <- es LS.\\ b,+                                                     let b' = L.insert y (L.delete x b),+                                                     b' `S.notMember` interior',+                                                     b' `S.member` is' ]+                   in closure interior' boundary'+          is' = S.fromList is+-- The basis properties imply that the set of all bases is connected under the "neighbour" relation+++-- |Given a matrix, represented as a list of rows, number the columns [1..],+-- and construct the matroid whose independent sets correspond to those sets of columns which are linearly independent+-- (or in case there are repetitions, those multisets of columns which are sets, and which are linearly independent).+vectorMatroid :: (Fractional k) => [[k]] -> Matroid Int+vectorMatroid = vectorMatroid' . L.transpose++-- |Given a list of vectors (or rows of a matrix), number the vectors (rows) [1..], and construct the matroid whose independent sets+-- correspond to those sets of vectors (rows) which are linearly independent+-- (or in case there are repetitions, those multisets which are sets, and which are linearly independent).+vectorMatroid' :: (Fractional k) => [[k]] -> Matroid Int+vectorMatroid' vs = fromBases (map fst vs') bs+    where vs' = zip [1..] vs+          bs = dfs [] [([],[],vs')]+          dfs ls (r@(i,ref,es) : rs) =+              let succs = successors r in+              if null succs then dfs (i:ls) (succs ++ rs) else dfs ls (succs ++ rs)+          dfs ls [] = let r = length $ last ls -- we know that the first one we found is a true base+                      in map reverse $ filter (\b -> length b == r) ls+                      -- we might have had null succs simply because we ran out of vectors to add+          successors (i,ref,es) = [(i',ref',es') | (j,e):es' <- L.tails es,+                                                   not (inSpanRE ref e),+                                                   let ref' = rowEchelonForm (ref ++ [e]), -- is this really better than e:ref?+                                                   let i' = j : i ]++-- |Given the edges of an undirected graph, number the edges [1..], and construct the matroid whose independent sets+-- correspond to those sets of edges which contain no cycle. The bases therefore correspond to maximal forests within the graph.+-- The edge set is allowed to contain loops or parallel edges.+cycleMatroid :: (Ord a) => [[a]] -> Matroid Int+cycleMatroid es = fromBases (map fst es') bs+    where es' = zip [1..] es+          bs = dfs [] [([], M.empty, es')]+          dfs ls (r@(i,ref,es) : rs) =+              let succs = successors r+              in if null succs then dfs (i:ls) (succs ++ rs) else dfs ls (succs ++ rs)+          dfs ls [] = let r = length $ last ls -- we know that the first one we found is a true base+                      in map reverse $ filter (\b -> length b == r) ls+                      -- we might have had null succs simply because we ran out of edges to add+          successors (i, reps, (j,[u,v]):es' ) =+              if u == v+              then successors (i, reps, es')+              else case (M.lookup u reps, M.lookup v reps) of+                   (Nothing, Nothing) -> (j:i, M.insert u u $ M.insert v u reps, es') : successors (i, reps, es')+                                         -- neither of these vertices has been seen before, so add this edge as a new tree in the forest+                   (Just u', Nothing) -> (j:i, M.insert v u' reps, es') : successors (i, reps, es')+                                         -- we have seen u before but not v, so add v to the u-tree+                   (Nothing, Just v') -> (j:i, M.insert u v' reps, es') : successors (i, reps, es')+                                         -- we have seen v before but not u, so add u to the v-tree+                   (Just u', Just v') -> if u' == v'+                                         then successors (i,reps,es')+                                              -- u and v are already in the same tree, so adding this edge would make a cycle+                                         else (j:i, M.map (\w -> if w == v' then u' else w) reps, es') : successors (i, reps, es')+                                              -- u and v are in different trees, so join the trees together+          successors (_, _, []) = []++-- A version of cycle matroid that retains the original edges rather than relabeling with integers.+-- Not really valid if there are parallel edges (because they can't be distinguished in the output).+-- Required for "markedfcim" below.+cycleMatroid' es = fmap lookupEdge $ cycleMatroid es+    where table = M.fromList $ zip [1..] es+          lookupEdge = (M.!) table+++-- |Given a matroid over an arbitrary type, relabel to obtain a matroid over the integers.+to1n :: (Ord a) => Matroid a -> Matroid Int+to1n m = fmap to1n' m+    where es = elements m+          table = M.fromList $ zip es [1..]+          to1n' = (M.!) table+++-- ISOMORPHISMS AND AUTOMORPHISMS++incidenceGraphB m = G.G vs' es'+    where es = elements m+          bs = bases m+          vs' = map Left es ++ map Right bs+          es' = L.sort [ [Left e, Right b] | b <- bs, e <- b ]+-- incidence graph for the matroid considered as an incidence structure between elements and bases++incidenceGraphH m = G.G vs' es'+    where es = elements m+          hs = L.sort $ hyperplanes m+          vs' = map Left es ++ map Right hs+          es' = L.sort [ [Left e, Right h] | h <- hs, e <- h ]+-- incidence graph for the matroid considered as an incidence structure between elements and hyperplanes+-- for "sparse" matroids, there are likely to be fewer hyperplanes than bases (why?)+-- incidenceGraphH may not be connected - eg u 2 4++-- (So a rank 3 or higher matroid, provided it has more than one hyperplane, certainly has a connected incidenceGraphH)++matroidIsos m1 m2 = incidenceIsos (incidenceGraphH m1) (incidenceGraphH m2)++-- |Are the two matroids isomorphic?+isMatroidIso :: (Ord a, Ord b) => Matroid a -> Matroid b -> Bool+isMatroidIso m1 m2 = isIncidenceIso (incidenceGraphH m1) (incidenceGraphH m2)++-- |Return the automorphisms of the matroid.+matroidAuts :: (Ord a) => Matroid a -> [Permutation a]+matroidAuts m+    | G.isConnected hgraph = incidenceAuts hgraph+    | G.isConnected bgraph = incidenceAuts bgraph+    | otherwise = error "matroidAuts: incidence graph is not connected"+    where hgraph = incidenceGraphH m+          bgraph = incidenceGraphB m+-- Note that the results aren't always what one intuitively expects from the geometric representation.+-- This is because geometric representations suggest additional structure beyond matroid structure.+-- For example, for the Vamos matroid v8,+-- it returns auts which are not "geometric" auts of the geometric representation.+-- Matroids are really combinatorial objects, not geometric.+-- In the case of v8, the key is that the only thing the auts have to preserve are the planes, not the lines+++-- CIRCUITS++-- |A circuit in a matroid is a minimal dependent set.+isCircuit :: (Ord a) => Matroid a -> [a] -> Bool+isCircuit m c =+    isDependent m c &&+    all (isIndependent m) (deletions c)++-- |Return all circuits for the given matroid, in shortlex order.+circuits :: (Ord a) => Matroid a -> [[a]]+circuits m = toShortlex $ dfs S.empty [L.insert e b | b <- bs, e <- es LS.\\ b]+    where es = elements m+          bs = bases m+          dfs vs (c:cs) | c `S.member` vs = dfs vs cs+                        | otherwise = let cs' = successors c+                                          vs' = S.insert c vs+                                      in if null cs' then c : dfs vs' cs else dfs vs' (cs' ++ cs)+          dfs _ [] = []+          successors c = [c' | c' <- deletions c, isDependent m c' ]++-- Oxley p10+-- |Are the given sets the circuits of some matroid?+isMatroidCircuits :: (Ord a) => [[a]] -> Bool+isMatroidCircuits cs =+    [] `notElem` cs &&+    and [(c1 `LS.isSubset` c2) `implies` (c1 == c2) | c1 <- cs, c2 <- cs] &&+    and [ exists [c3 | c3 <- cs, c3 `LS.isSubset` c12']+        | c1 <- cs, c2 <- cs, c1 /= c2,+          e <- c1 `LS.intersect` c2, let c12' = L.delete e (c1 `LS.union` c2)] ++-- |Reconstruct a matroid from its elements and circuits.+fromCircuits :: (Ord a) => [a] -> [[a]] -> Matroid a+fromCircuits es cs = fromBases es bs+    where b = greedy [] es -- first find any basis+          greedy ls (r:rs) = let ls' = ls ++ [r] in+                             if isIndep ls'+                             then greedy ls' rs+                             else greedy ls rs+          greedy ls [] = ls+          bs = closure S.empty (S.singleton b) -- now find all other bases by passing to "neighbouring" bases+          closure interior boundary =+              if S.null boundary+              then S.toList interior+              else let interior' = interior `S.union` boundary+                       boundary' = S.fromList [ b' | b <- S.toList boundary,+                                                     x <- b, y <- es LS.\\ b,+                                                     let b' = L.insert y (L.delete x b),+                                                     b' `S.notMember` interior',+                                                     isIndep b' ]+                   in closure interior' boundary'+          isIndep xs = not (any (`LS.isSubset` xs) cs)+++-- |An element e in a matroid M is a loop if {e} is a circuit of M.+isLoop :: (Ord a) => Matroid a -> a -> Bool+isLoop m e = isCircuit m [e]+-- isLoop (M es is) e = [e] `notElem` is++-- |Elements f and g in a matroid M are parallel if {f,g} is a circuit of M.+isParallel :: (Ord a) => Matroid a -> a -> a -> Bool+isParallel m f g = isCircuit m [f,g]++-- |A matroid is simple if it has no loops or parallel elements+isSimple :: (Ord a) => Matroid a -> Bool+isSimple m = all ( (>2) . length ) (circuits m)+++-- BASES++-- |A base or basis in a matroid is a maximal independent set.+isBase :: (Ord a) => Matroid a -> [a] -> Bool+isBase (M es bs) b = b `tsmember` bs++-- |Return all bases for the given matroid+bases :: (Ord a) => Matroid a -> [[a]]+bases (M es bs) = tstolist bs++-- |Are the given sets the bases of some matroid?+isMatroidBases :: (Ord a) => [[a]] -> Bool+isMatroidBases bs =+    (not . null) bs &&+    and [ exists [y | y <- b2 LS.\\ b1, L.insert y (L.delete x b1) `elem` bs]+        | b1 <- bs, b2 <- bs, x <- b1 LS.\\ b2 ]++-- |Reconstruct a matroid from its elements and bases.+fromBases :: (Ord a) => [a] -> [[a]] -> Matroid a+fromBases es bs = M es (tsfromlist bs)+-- The elements are required because a loop does not appear in any basis.++-- Oxley p17+-- |Given a matroid m, a basis b, and an element e, @fundamentalCircuit m b e@ returns the unique circuit contained in b union {e},+-- which is called the fundamental circuit of e with respect to b.+fundamentalCircuit :: (Ord a) => Matroid a -> [a] -> a -> [a]+fundamentalCircuit m b e = unique [c | c <- circuits m, c `LS.isSubset` be]+    where be = L.insert e b++uniformMatroid m n | m <= n = fromBases es bs+    where es = [1..n]+          bs = combinationsOf m es++-- |The uniform matroid U m n is the matroid whose independent sets are all subsets of [1..n] with m or fewer elements.+u :: Int -> Int -> Matroid Int+u = uniformMatroid+++-- RANK FUNCTION++-- Oxley p103, 3.1.14+restriction1 m xs = fromBases xs bs'+    where bs = bases m+          is' = toShortlex [b `LS.intersect` xs | b <- bs]+          r = length $ last is'+          bs' = dropWhile ( (< r) . length ) is'++-- |The restriction of a matroid to a subset of its elements+restriction :: (Ord a) => Matroid a -> [a] -> Matroid a+restriction m@(M es bs) xs = M xs bs'+    where (_,bs') = balance $ prune bs+          prune (TS yts) = let (ins, outs) = L.partition (\(y,t) -> y `elem` xs) yts+                               ins' = [(y, prune t) | (y,t) <- ins]+                               outs' = concat [ zts | (y,t) <- outs, let TS zts = prune t ]+                           in TS $ ins' ++ outs'+          balance (TS yts) = let dyt's = [(d',(y,t')) | (y,t) <- yts, let (d',t') = balance t]+                                 d = maximum $ 0 : map fst dyt's+                             in (d+1, TS $ toSet [(y,t') | (d',(y,t')) <- dyt's, d' == d])+-- we have to make the toSet call *after* we balance, otherwise two trees may appear unequal because of undergrowth that is going to be removed++-- !! Need thorough testing to prove that restriction == restriction1+++-- |Given a matroid m, @rankfun m@ is the rank function on subsets of its element set+rankfun :: (Ord a) => Matroid a -> [a] -> Int+rankfun m xs = (length . head . bases) (restriction m xs)+-- no danger of head [], because bases must be non-null++-- |The rank of a matroid is the cardinality of a basis+rank :: (Ord a) => Matroid a -> Int+rank m = length $ head $ bases m+-- rank m@(M es bs) = rankfun m es++-- Oxley p23+-- |Reconstruct a matroid from its elements and rank function+fromRankfun :: (Ord a) => [a] -> ([a] -> Int) -> Matroid a+fromRankfun es rkf = fromBases es bs+    where b = greedy 0 [] es -- first find any basis+          greedy rk ls (r:rs) = let ls' = ls ++ [r] in+                                if rkf ls' == rk+1+                                then greedy (rk+1) ls' rs+                                else greedy rk ls rs+          greedy _ ls [] = ls+          rk = rkf b+          isBasis b' = rkf b' == rk+          bs = closure S.empty (S.singleton b) S.empty -- now find all other bases by passing to "neighbouring" bases+          closure interior boundary exterior =+              if S.null boundary+              then S.toList interior+              else let interior' = interior `S.union` boundary+                       candidates = S.fromList [ b' | b <- S.toList boundary,+                                                      x <- b, y <- es LS.\\ b,+                                                      let b' = L.insert y (L.delete x b),+                                                      b' `S.notMember` interior',+                                                      b' `S.notMember` exterior ]+                       (boundary', exterior') = S.partition isBasis candidates+                   in closure interior' boundary' (S.union exterior exterior')+-- The purpose of keeping track of exterior points we have encountered+-- is to avoid making repeat calls to the rank function rkf+++-- CLOSURE OPERATOR AND FLATS++-- |Given a matroid m, @closure m@ is the closure operator on subsets of its element set+closure :: (Ord a) => Matroid a -> [a] -> [a]+closure m xs = [x | x <- es, x `elem` xs || rankfun m (L.insert x xs) == rankxs]+    where es = elements m+          rankxs = rankfun m xs+-- The intuition is that closure xs is all elements within the span of xs++-- |Reconstruct a matroid from its elements and closure operator+fromClosure :: (Ord a) => [a] -> ([a] -> [a]) -> Matroid a+fromClosure es cl = fromBases es bs+    where b = greedy (cl []) [] es -- first find any basis+          greedy span ls (r:rs) = let ls' = ls ++ [r] in+                                  if r `notElem` span -- r is independent relative to ls+                                  then greedy (cl ls') ls' rs+                                  else greedy span ls rs+          greedy _ ls [] = ls+          rk = length b+          isBasis b' = cl b' == es+          bs = closure S.empty (S.singleton b) S.empty -- now find all other bases by passing to "neighbouring" bases+          closure interior boundary exterior =+              if S.null boundary+              then S.toList interior+              else let interior' = interior `S.union` boundary+                       candidates = S.fromList [ b' | b <- S.toList boundary,+                                                      x <- b, y <- es LS.\\ b,+                                                      let b' = L.insert y (L.delete x b),+                                                      b' `S.notMember` interior',+                                                      b' `S.notMember` exterior ]+                       (boundary', exterior') = S.partition isBasis candidates+                   in closure interior' boundary' (S.union exterior exterior')+-- The purpose of keeping track of exterior points we have encountered+-- is to avoid making repeat calls to the closure operator cl++{-+-- Not quite sure why the following is so much slower than the above+-- May just be because Data.Set is compiled, and this would be comparable if also compiled+fromClosure2 es cl = M es bs+    where b = greedy (cl []) [] es -- first find any basis+          greedy span ls (r:rs) = let ls' = ls ++ [r] in+                                  if r `notElem` span -- r is independent relative to ls+                                  then greedy (cl ls') ls' rs+                                  else greedy span ls rs+          greedy _ ls [] = ls+          rk = length b+          isBasis b' = cl b' == es+          bs = closure tsempty [b] tsempty -- now find all other bases by passing to "neighbouring" bases+          closure interior boundary exterior =+              if null boundary+              then interior+              else let interior' = foldl' (flip tsinsert) interior boundary+                       candidates = [ b' | b <- boundary,+                                           x <- b, y <- es LS.\\ b,+                                           let b' = L.insert y (L.delete x b),+                                           not (b' `tsmember` interior'),+                                           not (b' `tsmember` exterior) ]+                       (boundary', exterior') = L.partition isBasis candidates+                   in closure interior' boundary' (foldl' (flip tsinsert) exterior exterior')+-}++-- |A flat in a matroid is a closed set, that is a set which is equal to its own closure+isFlat :: (Ord a) => Matroid a -> [a] -> Bool+isFlat m xs = closure m xs == xs++flats1 m = [xs | xs <- powersetbfs es, isFlat m xs]+    where es = elements m+-- first, inefficient, implementation++-- given xs, a flat in m, return the flats which cover xs+-- these have the property that they partition es \\ xs+coveringFlats m xs = coveringFlats' (es LS.\\ xs)+    where es = elements m+          coveringFlats' (y:ys) = let zs = closure m (L.insert y xs)+                                  in zs : coveringFlats' (ys LS.\\ zs)+          coveringFlats' [] = []++-- since we are dealing with finite matroids, the lattice of flats is finite, so it has a minimal element+minimalFlat m = head $ filter (isFlat m) $ powersetbfs $ elements m++-- |The flats of a matroid are its closed sets. They form a lattice under inclusion.+flats :: (Ord a) => Matroid a -> [[a]]+flats m = flats' S.empty [minimalFlat m]+    where flats' ls (r:rs) = if r `S.member` ls+                             then flats' ls rs+                             else flats' (S.insert r ls) (rs ++ coveringFlats m r)+          flats' ls [] = toShortlex $ S.toList ls+          -- this is just breadth-first search++-- isMatroidFlats+-- Oxley p31-32+{-+isMatroidFlats es fs =+    es `elem` fs &&+    [ (f1 `LS.intersect` f2) `elem` fs | f1 <- fs, f2 <- fs ] &&+    -- for all flats f, the minimal flats containing f partition es-f+-}++-- |Reconstruct a matroid from its flats. (The flats must be given in shortlex order.)+fromFlats :: (Ord a) => [[a]] -> Matroid a+fromFlats fs | isShortlex fs = fromFlats' fs+             | otherwise = error "fromFlats: flats must be in shortlex order"++fromFlats' fs = fromClosure es cl+    where es = last fs -- es is a flat, and last in shortlex order+          cl xs = head [f | f <- fs, xs `LS.isSubset` f] -- the first flat is minimal, because of shortlex order+-- !! we can probably do better (efficiency-wise)+-- eg by constructing the lattice as a DAG, and climbing up it++-- |A subset of the elements in a matroid is spanning if its closure is all the elements+isSpanning :: (Ord a) => Matroid a -> [a] -> Bool+isSpanning m xs = closure m xs == es+    where es = elements m++-- |A hyperplane is a flat whose rank is one less than that of the matroid+isHyperplane :: (Ord a) => Matroid a -> [a] -> Bool+isHyperplane m xs = isFlat m xs && rankfun m xs == rank m - 1++hyperplanes1 m = [h | h <- flats m, rankfun m h == rk - 1]+    where rk = rank m++-- Oxley p65: h is a hyperplane iff its complement is a cocircuit+hyperplanes :: (Ord a) => Matroid a -> [[a]]+hyperplanes m = toShortlex $ map complement $ cocircuits m+    where es = elements m+          complement cc = es LS.\\ cc+-- This appears to be faster than hyperplanes1++-- Oxley p70+isMatroidHyperplanes :: (Ord a) => [a] -> [[a]] -> Bool+isMatroidHyperplanes es hs =+    es `notElem` hs &&+    isClutter hs &&+    and [ exists [h3 | h3 <- hs, h12e `LS.isSubset` h3] | (h1,h2) <- pairs hs,+                                                          e <- es LS.\\ (LS.union h1 h2),+                                                          let h12e = L.insert e (LS.intersect h1 h2) ]+-- Note that contrary to what one might initially think+-- it does *not* follow from the third condition that every element is in some hyperplane+-- since there might be only one hyperplane, eg fromBases [1,2] [[2]] has [1] as it's only hyperplane++-- Haven't actually proven that the following is always correct+-- but it seems to work+fromHyperplanes1 es hs = fromFlats $ closure S.empty (S.fromList hs)+    where closure interior boundary =+              if S.null boundary+              then (toShortlex $ S.toList interior) ++ [es]+              else let interior' = S.union interior boundary+                       boundary' = S.fromList [ f1 `LS.intersect` f2 | (f1,f2) <- pairs (S.toList boundary) ]+                                   S.\\ interior'+                   in closure interior' boundary'++-- |Reconstruct a matroid from its elements and hyperplanes+fromHyperplanes :: (Ord a) => [a] -> [[a]] -> Matroid a+fromHyperplanes es hs = fromCocircuits es $ map complement hs+    where fromCocircuits es = dual . fromCircuits es+          complement xs = es LS.\\ xs+++-- GEOMETRIC REPRESENTATION++-- |Given a list of points in k^n, number the points [1..], and construct the matroid whose independent sets+-- correspond to those sets of points which are affinely independent.+--+-- A multiset of points in k^n is said to be affinely dependent if it contains two identical points,+-- or three collinear points, or four coplanar points, or ... - and affinely independent otherwise.+affineMatroid :: (Fractional k) => [[k]] -> Matroid Int+affineMatroid vs = vectorMatroid' $ map (1:) vs++-- |fromGeoRep returns a matroid from a geometric representation consisting of dependent flats of various ranks.+-- Given lists of dependent rank 0 flats (loops), rank 1 flats (points), rank 2 flats (lines) and rank 3 flats (planes),+-- @fromGeoRep loops points lines planes@ returns the matroid having these as dependent flats.+-- Note that if all the elements lie in the same plane, then this should still be listed as an argument.+fromGeoRep :: (Ord a) => [[a]] -> [[a]] -> [[a]] -> [[a]] -> Matroid a+fromGeoRep loops points lines planes =+    fromCircuits es $ minimal $+        loops +++        concatMap (combinationsOf 2) points +++        concatMap (combinationsOf 3) lines +++        concatMap (combinationsOf 4) planes +++        combinationsOf 5 es+    where es = toSet $ concat loops ++ concat points ++ concat lines ++ concat planes++-- Note that we don't check that the inputs are valid++-- xss assumed to be in shortlex order+minimal xss = minimal' [] xss+    where minimal' ls (r:rs) = if any (`LS.isSubset` r) ls+                               then minimal' ls rs+                               else minimal' (r:ls) rs+          minimal' ls [] = reverse ls++-- |A simple matroid has no loops or parallel elements, hence its geometric representation has no loops or dependent points.+-- @simpleFromGeoRep lines planes@ returns the simple matroid having these dependent flats.+simpleFromGeoRep :: (Ord a) => [[a]] -> [[a]] -> Matroid a+simpleFromGeoRep lines planes = fromGeoRep [] []  lines planes++-- Oxley p37-8+isSimpleGeoRep lines planes =+    all ( (<= 1) . length ) [ l1 `LS.intersect` l2 | (l1,l2) <- pairs lines ] &&+    all ( \i -> length i <= 2 || i `elem` lines ) [ p1 `LS.intersect` p2 | (p1,p2) <- pairs planes ] &&+    and [ any (u `LS.isSubset`) planes | (l1,l2) <- pairs lines, length (l1 `LS.intersect` l2) == 1, let u = l1 `LS.union` l2 ] &&+    and [ length i == 1 || i == l | l <- lines, p <- planes, let i = l `LS.intersect` p ]+++isCircuitHyperplane m xs = isCircuit m xs && isHyperplane m xs++-- |List the circuit-hyperplanes of a matroid.+circuitHyperplanes :: (Ord a) => Matroid a -> [[a]]+circuitHyperplanes m = [ h | h <- hyperplanes m, isCircuit m h ] ++-- Oxley p39+-- |Given a matroid m, and a set of elements b which is both a circuit and a hyperplane in m,+-- then @relaxation m b@ is the matroid which is obtained by adding b as a new basis.+-- This corresponds to removing b from the geometric representation of m.+relaxation :: (Ord a) => Matroid a -> [a] -> Matroid a+relaxation m b+    | isCircuitHyperplane m b = fromBases es bs+    | otherwise = error "relaxation: not a circuit-hyperplane"+    where es = elements m+          bs = b : bases m +++-- TRANSVERSAL MATROIDS++ex161 = [ [1,2,6], [3,4,5,6], [2,3], [2,4,6] ]++-- the edges of the bipartite graph+transversalGraph as = [(Left x, Right i) | (a,i) <- zip as [1..], x <- a]++partialMatchings es = dfs [(S.empty, [], es)]+    where dfs (node@(vs,ls,rs): nodes) = ls : dfs (successors node ++ nodes)+          dfs [] = []+          successors (vs,ls,rs) =  [ (S.insert u $ S.insert v vs, L.insert r ls, rs')+                                   | r@(u,v):rs' <- L.tails rs, u `S.notMember` vs, v `S.notMember` vs ]++-- |Given a set of elements es, and a sequence as = [a1,...,am] of subsets of es,+-- return the matroid whose independent sets are the partial transversals of the as.+transversalMatroid :: (Ord a) => [a] -> [[a]] -> Matroid a+transversalMatroid es as = fromBases es bs+    where is@(i:_) = reverse $ toShortlex $ toSet $ (map . map) (unLeft . fst) $ partialMatchings (transversalGraph as)+          unLeft (Left x) = x+          l = length i+          bs = reverse $ takeWhile ( (== l) . length ) is+-- In this case, as is called a presentation of the matroid+-- Note that there may be partial transversals even if there are no transversals+++-- Not obvious how to efficiently find the bases without finding the independent sets,+-- since neighbouring partial transversals might arise in different ways+-- (Although Oxley p93 seems to be saying that they don't)+{-+transversalMatroid2 es as = fromBases es bs+    where es' = transversalGraph as+          b = greedy [] es' -- first find any basis+          greedy ls (r@(u,v):rs) = if u `notElem` map fst ls && v `notElem` map snd ls+                                   then greedy (r:ls) rs+                                   else greedy ls rs+          greedy ls [] = ls+-- now choose the subset of the as indicated by the Right part of b+-- and seek all full transversals+          bs = closure S.empty (S.singleton b) -- now find all other bases by passing to "neighbouring" bases+          closure interior boundary =+              if S.null boundary+              then S.toList interior+              else let interior' = interior `S.union` boundary+                       boundary' = S.fromList [ b' | b <- S.toList boundary,+                                                     x <- b, y <- es LS.\\ b,+                                                     let b' = L.insert y (L.delete x b),+                                                     b' `S.notMember` interior',+                                                     cl b' == es ] -- b' is spanning, and same size as b, hence must be basis+                                   -- S.\\ interior'+                   in closure interior' boundary'+-}+++-- DUALITY++-- |The dual matroid+dual :: (Ord a) => Matroid a -> Matroid a+dual m = fromBases es bs'+    where es = elements m+          bs = bases m+          bs' = L.sort $ map (es LS.\\) bs++isCoindependent m xs = isIndependent (dual m) xs++isCobase m xs = isBase (dual m) xs+-- quicker but less clear to calculate this directly++isCocircuit m xs = isCircuit (dual m) xs++cocircuits :: (Ord a) => Matroid a -> [[a]]+cocircuits m = circuits (dual m)++isColoop m e = isLoop (dual m) e++isCoparallel m f g = isParallel (dual m) f g+++-- MINORS++deletion :: (Ord a) => Matroid a -> [a] -> Matroid a+deletion m xs = restriction m (es LS.\\ xs)+    where es = elements m++(\\\) = deletion++contraction :: (Ord a) => Matroid a -> [a] -> Matroid a+contraction m xs = dual (deletion (dual m) xs)++(///) = contraction+++-- CONNECTIVITY++-- Oxley p120+-- |A matroid is (2-)connected if, for every pair of distinct elements, there is a circuit containing both+isConnected :: (Ord a) => Matroid a -> Bool+isConnected m = and [any (pair `LS.isSubset`) cs | pair <- combinationsOf 2 es]+    where es = elements m+          cs = circuits m++component m x = closure S.empty (S.singleton x)+    where cs = circuits m+          closure interior boundary =+              if S.null boundary+              then S.toList interior+              else let interior' = S.union interior boundary+                       boundary' = S.fromList (concat [c | c <- cs, (not . null) (LS.intersect c (S.toList boundary)) ])+                                   S.\\ interior'+                   in closure interior' boundary'++-- |The direct sum of two matroids+dsum :: (Ord a, Ord b) => Matroid a -> Matroid b -> Matroid (Either a b)+dsum m1 m2 = fromBases es bs+    where es = map Left (elements m1) ++ map Right (elements m2)+          bs = [map Left b1 ++ map Right b2 | b1 <- bases m1, b2 <- bases m2]+++-- REPRESENTABILITY++-- |@matroidPG n fq@ returns the projective geometry PG(n,Fq), where fq is a list of the elements of Fq+matroidPG :: (Fractional a) => Int -> [a] -> Matroid Int+matroidPG n fq = vectorMatroid' $ ptsPG n fq++-- |@matroidAG n fq@ returns the affine geometry AG(n,Fq), where fq is a list of the elements of Fq+matroidAG :: (Fractional a) => Int -> [a] -> Matroid Int+matroidAG n fq = vectorMatroid' $ ptsAG n fq+++-- Oxley p182+-- |Given a matroid m, the fundamental-circuit incidence matrix relative to a base b+-- has rows indexed by the elements of b, and columns indexed by the elements not in b.+-- The bi, ej entry is 1 if bi is in the fundamental circuit of ej relative to b, and 0 otherwise.+fundamentalCircuitIncidenceMatrix :: (Ord a, Num k) => Matroid a -> [a] -> [[k]]+fundamentalCircuitIncidenceMatrix m b = L.transpose $ fundamentalCircuitIncidenceMatrix' m b++fundamentalCircuitIncidenceMatrix' m b =+    [ [if e `elem` fundamentalCircuit m b e' then 1 else 0 | e <- b]+    | e' <- elements m LS.\\ b ]++fcim = fundamentalCircuitIncidenceMatrix+fcim' = fundamentalCircuitIncidenceMatrix'++-- Then fcim w4 [1..4] == L.transpose (fcim (dual w4) [5..8])+++-- Given a matrix of 0s and 1s, return a matrix of 0s, 1s and *s+-- where 0 -> 0, and 1 -> 1 if it is the first 1 in either a row or column, 1 -> * otherwise+markNonInitialRCs mx = mark (replicate w False) mx+    where w = length $ head mx -- the width+          mark cms (r:rs) = let (cms', r') = mark' False [] cms [] r in r' : mark cms' rs+          mark _ [] = []+          mark' rm cms' (cm:cms) ys (x:xs)+              | x == 0 = mark' rm (cm:cms') cms (Zero:ys) xs+              | x == 1 = if rm && cm+                         then mark' True (True:cms') cms (Star:ys) xs+                         else mark' True (True:cms') cms (One:ys) xs+          mark' _ cms' [] ys [] = (reverse cms', reverse ys)++-- Given a matrix of 0s, 1s and *s, return all distinct matrices that can be obtained+-- by substituting non-zero elements of Fq for the *s+substStars mx fq = substStars' mx+    where fq' = tail fq -- non-zero elts of fq+          substStars' (r:rs) = [r':rs' | r' <- substStars'' r, rs' <- substStars' rs]+          substStars' [] = [[]]+          substStars'' (Zero:xs) = map (0:) $ substStars'' xs+          substStars'' (One:xs) = map (1:) $ substStars'' xs+          substStars'' (Star:xs) = [x':xs' | x' <- fq', xs' <- substStars'' xs]+          substStars'' [] = [[]]++starSubstitutionsV fq' (Zero:xs) = map (0:) $ starSubstitutionsV fq' xs+starSubstitutionsV fq' (One:xs) = map (1:) $ starSubstitutionsV fq' xs+starSubstitutionsV fq' (Star:xs) = [x':xs' | x' <- fq', xs' <- starSubstitutionsV fq' xs]+starSubstitutionsV _ [] = [[]]+++-- Oxley p184-5+-- Note that the particular representations you get depend on which basis is used+-- (Perhaps we should let you pass in the basis to use)+representations1 fq m = [ L.transpose d | d <- substStars dhash fq, let mx = ir ++ d,+                                          to1n m == (vectorMatroid' $ map snd $ L.sort $ zip (b ++ b') mx) ]+    where b = head $ bases m+          b' = elements m LS.\\ b+          r = length b -- rank of the matroid+          ir = idMx r -- identity matrix+          dhash = markNonInitialRCs $ fcim' m b++-- edges of the fundamental circuit incidence graph+fcig m b = [ [e,e'] | e <- b, e' <- elements m LS.\\ b, e `elem` fundamentalCircuit m b e' ]++-- the fcim of m relative to b, with 1s and *s marked+markedfcim m b = mark b b' (fcim m b)+    where b' = elements m LS.\\ b -- the elements of b are the row labels, those of b' are the column labels+          entries = fcig m b -- the [row,column] coordinates of the non-zero entries in the fcim+          ones = head $ bases $ cycleMatroid' entries -- a set of entries which we can take to be 1+          stars = entries LS.\\ ones -- the set of entries which we are then still free to assign+          mark (i:is) js (r:rs) = (mark' i js r) : mark is js rs+          mark [] _ [] = []+          mark' i (j:js) (x:xs)+              | x == 0 = Zero : mark' i js xs+              | x == 1 = (if [i,j] `elem` stars then Star else One) : mark' i js xs+          mark' _ [] [] = []+-- The markedfcim will sometimes do better than markNonInitialRCs, and is best possible++-- Oxley p184-5+representations2 fq m = [ L.transpose mx | d <- substStars dhash' fq, let mx = ir ++ d,+                                           m' == (vectorMatroid' $ map snd $ L.sort $ zip (b ++ b') mx) ]+    where m' = to1n m+          es = elements m+          b = head $ bases m+          b' = es LS.\\ b+          r = length b -- rank of the matroid+          ir = idMx r -- identity matrix+          dhash' = L.transpose $ markedfcim m b++-- In the following, we check the representations of the restrictions as we add a column at a time+-- This enables us to early out if a potential representation doesn't work on a restriction++-- |Find representations of the matroid m over fq. Specifically, this function will find one representative+-- of each projective equivalence class of representation.+representations :: (Fractional fq, Ord a) => [fq] -> Matroid a -> [[[fq]]]+representations fq m = map L.transpose $ representations' (reverse $ zip b ir) (zip b' dhash')+    where fq' = tail fq -- fq \ {0}+          b = head $ bases m+          b' = elements m LS.\\ b+          r = length b -- rank of the matroid+          ir = idMx r -- identity matrix+          dhash' = L.transpose $ markedfcim m b+          representations' ls ((i,r):rs) = concat+              [ representations' ((i,r'):ls) rs+              | r' <- starSubstitutionsV fq' r,+                let (is,vs) = unzip $ L.sortBy cmpfst ((i,r'):ls),+                to1n (restriction m is) == (vectorMatroid' vs) ]+          representations' ls [] = [map snd $ reverse ls]++-- |Is the matroid representable over Fq? For example, to find out whether a matroid m is binary, evaluate @isRepresentable f2 m@.+isRepresentable :: (Fractional fq, Ord a) => [fq] -> Matroid a -> Bool+isRepresentable fq m = (not . null) (representations fq m)++-- |A binary matroid is a matroid which is representable over F2+isBinary :: (Ord a) => Matroid a -> Bool+isBinary = isRepresentable f2++-- |A ternary matroid is a matroid which is representable over F3+isTernary :: (Ord a) => Matroid a -> Bool+isTernary = isRepresentable f3+++-- CONSTRUCTIONS+-- The next three functions not very thoroughly tested++data LMR a b = L a | Mid | R b deriving (Eq, Ord, Show)++-- Oxley p252+seriesConnection (m1,p1) (m2,p2)+    | not (isLoop m1 p1) && not (isColoop m1 p1) && not (isLoop m2 p2) && not (isColoop m2 p2) =+        fromCircuits es cs+    | otherwise = error "not yet implemented"+        where es = map L (elements m1 LS.\\ [p1]) ++ [Mid] ++ map R (elements m2 LS.\\ [p2])+              cs = (map . map) L (circuits $ m1 \\\ [p1]) +++                   (map . map) R (circuits $ m2 \\\ [p2]) +++                   [ map L (L.delete p1 c1) ++ [Mid] ++ map R (L.delete p2 c2)+                   | c1 <- circuits m1, p1 `elem` c1, c2 <- circuits m2, p2 `elem` c2]++parallelConnection (m1,p1) (m2,p2)+    | not (isLoop m1 p1) && not (isColoop m1 p1) && not (isLoop m2 p2) && not (isColoop m2 p2) =+        fromCircuits es cs+    | otherwise = error "not yet implemented"+        where es = map L (elements m1 LS.\\ [p1]) ++ [Mid] ++ map R (elements m2 LS.\\ [p2])+              cs = (map . map) L (circuits $ m1 \\\ [p1]) +++                   [ map L (L.delete p1 c1) ++ [Mid] | c1 <- circuits m1, p1 `elem` c1 ] +++                   (map . map) R (circuits $ m2 \\\ [p2]) +++                   [ [Mid] ++ map R (L.delete p2 c2) | c2 <- circuits m2, p2 `elem` c2 ] +++                   [ map L (L.delete p1 c1) ++ map R (L.delete p2 c2)+                   | c1 <- circuits m1, p1 `elem` c1, c2 <- circuits m2, p2 `elem` c2 ]++twoSum (m1,p1) (m2,p2)+    | not (isLoop m1 p1) && not (isColoop m1 p1) && not (isLoop m2 p2) && not (isColoop m2 p2) =+        fromCircuits es cs+    | otherwise = error "not yet implemented"+        where es = map L (elements m1 LS.\\ [p1]) ++ [Mid] ++ map R (elements m2 LS.\\ [p2])+              cs = (map . map) L (circuits $ m1 \\\ [p1]) +++                   (map . map) R (circuits $ m2 \\\ [p2]) +++                   [ map L (L.delete p1 c1) ++ map R (L.delete p2 c2)+                   | c1 <- circuits m1, p1 `elem` c1, c2 <- circuits m2, p2 `elem` c2]+-- Note: The only different from seriesConnection is that we don't have Mid in the last set++-- Oxley p427+matroidUnion m1 m2 = fromBases es bs+    where es = LS.union (elements m1) (elements m2)+          is = toShortlex $ toSet [ LS.union b1 b2 | b1 <- bases m1, b2 <- bases m2 ]+          r = length $ last is+          bs = dropWhile ( (< r) . length ) is++++-- SOME INTERESTING MATROIDS++-- |The Fano plane F7 = PG(2,F2)+f7 :: Matroid Int+f7 = fromGeoRep [] [] [[1,2,3],[1,4,7],[1,5,6],[2,4,6],[2,5,7],[3,4,5],[3,6,7]] [[1..7]]++-- Oxley p36, fig 1.12:+-- cycleMatroid [[1,2],[1,3],[2,3],[3,4],[2,4],[1,4]] == restriction fanoPlane [1..6]++-- |F7-, the relaxation of the Fano plane by removal of a line+f7m :: Matroid Int+f7m = relaxation f7 [2,4,6]++-- Oxley p39+-- |The Pappus configuration from projective geometry+pappus :: Matroid Int+pappus = fromGeoRep [] [] [[1,2,3],[1,5,7],[1,6,8],[2,4,7],[2,6,9],[3,4,8],[3,5,9],[4,5,6],[7,8,9]] [[1..9]]++-- more logical relabeling+-- pappus' = fromGeoRep [] [] [[1,2,3],[1,5,9],[1,6,8],[2,4,9],[2,6,7],[3,4,8],[3,5,7],[4,5,6],[7,8,9]] [[1..9]]++-- |Relaxation of the Pappus configuration by removal of a line+nonPappus :: Matroid Int+nonPappus = relaxation pappus [7,8,9]+-- fromGeoRep [] [] [[1,2,3],[1,5,7],[1,6,8],[2,4,7],[2,6,9],[3,4,8],[3,5,9],[4,5,6]] [[1..9]]++-- |The Desargues configuration+desargues :: Matroid Int+desargues = fromGeoRep [] [] [[1,2,5],[1,3,6],[1,4,7],[2,3,8],[2,4,9],[3,4,10],[5,6,8],[5,7,9],[6,7,10],[8,9,10]]+                              [[1,2,3,5,6,8],[1,2,4,5,7,9],[1,3,4,6,7,10],[2,3,4,8,9,10],[5,6,7,8,9,10]]+-- desargues == cycleMatroid (combinationsOf 2 [1..5])+-- (ie Desargues = M(K5) )+-- interestingly, although these are all the dependent flats that are evident from the diagram,+-- there are also some rank 3 flats consisting of a line and an "antipodal" point+-- eg, in terms of K5, [1,8,9,10] corresponds to [[1,2],[3,4],[3,5],[4,5]]++-- Oxley p71+vamosMatroid1 = fromHyperplanes [1..8] (hs4 ++ hs3)+    where hs4 = [ [1,2,3,4], [1,4,5,6], [2,3,5,6], [1,4,7,8], [2,3,7,8] ]+          hs3 = [ h3 | h3 <- combinationsOf 3 [1..8], null [h4 | h4 <- hs4, h3 `LS.isSubset` h4] ]+vamosMatroid = fromGeoRep [] [] [] [[1,2,3,4],[1,4,5,6],[2,3,5,6],[1,4,7,8],[2,3,7,8]]++-- |The Vamos matroid V8. It is not representable over any field.+v8 :: Matroid Int+v8 = vamosMatroid+-- v8 is self-dual (isomorphic to its own dual)++-- Oxley p188+-- |P8 is a minor-minimal matroid that is not representable over F4, F8, F16, ... .+-- It is Fq-representable if and only if q is not a power of 2.+p8 :: Matroid Int+p8 = vectorMatroid $+    ( [ [1,0,0,0,  0, 1, 1,-1],+        [0,1,0,0,  1, 0, 1, 1],+        [0,0,1,0,  1, 1, 0, 1],+        [0,0,0,1, -1, 1, 1, 0] ] :: [[F3]] )++p8' = fromGeoRep [] [] [] [ [1,2,3,8], [1,2,4,7], [1,3,4,6], [1,4,5,8], [1,5,6,7], [2,3,4,5], [2,3,6,7], [2,5,6,8], [3,5,7,8], [4,6,7,8] ]++-- |P8- is a relaxation of P8. It is Fq-representable if and only if q >= 4.+p8m :: Matroid Int+p8m = relaxation p8 [2,3,6,7]++-- |P8-- is a relaxation of P8-. It is a minor-minimal matroid that is not representable over F4.+-- It is Fq-representable if and only if q >= 5.+p8mm :: Matroid Int+p8mm = relaxation p8m [1,4,5,8]+++-- Oxley p317+-- r-spoked wheel graph+wheelGraph r = G.G vs es+    where vs = [0..r]+          es = [ [0,i] | i <- [1..r] ] ++ [ [i,i+1] | i <- [1..r-1] ] ++ [ [1,r] ]+-- for normal form, should L.sort es++mw4 = cycleMatroid $ G.edges $ wheelGraph 4+-- has [5,6,7,8] as unique circuit-hyperplane++w4' = relaxation mw4 $ unique $ circuitHyperplanes mw4 -- [5,6,7,8]++-- Oxley p183+w4 = fromGeoRep [] [] [[1,2,5],[1,4,8],[2,3,6],[3,4,7]] [[1,2,3,5,6],[1,2,4,5,8],[1,3,4,7,8],[2,3,4,6,7]]+++-- BINARY MATROIDS++-- Oxley p344+isBinary2 m = all (even . length) [ c `LS.intersect` cc | c <- circuits m, cc <- cocircuits m ]+++-- RANK POLYNOMIAL+-- Godsil & Royle p356++[x,y] = map glexVar ["x","y"] :: [GlexPoly Integer String]++-- first naive version+rankPoly1 m = sum [ x^(rm - r a) * y^(rm' - r' a') | a <- powersetdfs es, let a' = es LS.\\ a ]+    where es = elements m+          rm = rank m+          r = rankfun m+          m' = dual m+          rm' = rank m'+          r' = rankfun m'++-- |Given a matroid m over elements es, the rank polynomial is a polynomial r(x,y),+-- which is essentially a generating function for the subsets of es, enumerated by size and rank.+-- It is efficiently calculated using deletion and contraction.+--+-- It has the property that r(0,0) is the number of bases in m, r(1,0) is the number of independent sets,+-- r(0,1) is the number of spanning sets. It can also be used to derive the chromatic polynomial of a graph,+-- the weight enumerator of a linear code, and more.+rankPoly :: (Ord a) => Matroid a -> GlexPoly Integer String+rankPoly m+    | null es = 1+    | isLoop m e = (1+y) * rankPoly (m \\\ [e]) -- deletion+    | isColoop m e = (1+x) * rankPoly (m /// [e]) -- contraction+    | otherwise = rankPoly (m \\\ [e]) + rankPoly (m /// [e])+    where es = elements m+          e = head es++numBases m = unwrap $ rankPoly m `bind` (\v -> case v of "x" -> 0; "y" -> 0)++numIndeps m = unwrap $ rankPoly m `bind` (\v -> case v of "x" -> 1; "y" -> 0)++numSpanning m = unwrap $ rankPoly m `bind` (\v -> case v of "x" -> 0; "y" -> 1)++-- It is then possible to derive the chromatic poly of a graph from the rankPoly of its cycle matroid, etc+++-- Oxley p586++-- How many independent sets are there of each size+indepCounts m = map length $ L.groupBy eqfst $ [(length i, i) | i <- indeps m]+-- relying on shortlex order++{-+-- The following tries to calculate the indep counts+-- as the Hilbert series of the Stanley-Reisner ring.+-- However, it's not giving the right answer.+indepCounts2 m = take r $ indepCounts' (circuits m)+    where n = length (elements m) -- the number of "variables"+          r = rank m+          n' = toInteger n+          indepCounts' [] = map (\k -> (k+n'-1) `choose` (n'-1)) [0..]+              -- the number of ways of choosing n-1 separators, so as to leave a product of k powers+          indepCounts' (c:cs) =+              let d = length c -- the "degree" of c+                  cs' = reduce [] $ toShortlex $ toSet $ map (LS.\\ c) cs -- the quotient of cs by c+              in indepCounts' cs <-> (replicate d 0 ++ indepCounts' cs')+          reduce ls (r:rs) = if any (`LS.isSubset` r) ls then reduce ls rs else reduce (r:ls) rs+          reduce ls [] = reverse ls+          (x:xs) <+> (y:ys) = (x+y) : (xs <+> ys)+          xs <+> ys = xs ++ ys -- one of them is null+          xs <-> ys = xs <+> map negate ys +-}+++-- Whitney numbers of the second kind+-- The number of flats of each rank+whitney2nd m = map length $ L.groupBy eqfst $ L.sort [(rankfun m f, f) | f <- flats m]++-- Whitney numbers of the first kind+-- The number of subsets (of the element set) of each rank+whitney1st m = alternatingSign $ map length $ L.groupBy eqfst $ L.sort [(rankfun m x, x) | x <- powersetdfs (elements m)]+    where alternatingSign (x:xs) = x : alternatingSign (map negate xs)+          alternatingSign [] = []++++-- TODO++-- 1. Sort out the isomorphism / automorphism code in the case where the incidence graph isn't connected++-- 2. We could generate the geometric representation from a matroid (provided its rank <= 4)+-- geoRep m = filter (isDependent m) (flats m)++-- 3. isMatroidFlats+

@@ -5,6 +5,7 @@  module Math.Combinatorics.Poset where +import Math.Common.ListSet as LS -- set operations on strictly ascending lists import Math.Algebra.Field.Base import Math.Combinatorics.FiniteGeometry import Math.Algebra.LinearAlgebra@@ -85,6 +86,7 @@ powerset [] = [[]] powerset (x:xs) = let p = powerset xs in p ++ map (x:) p +{- -- subset test for sorted lists isSubset (x:xs) (y:ys) =     case compare x y of@@ -93,10 +95,11 @@     GT -> isSubset (x:xs) ys isSubset [] _ = True isSubset _ [] = False+-}  -- |posetB n is the lattice of subsets of [1..n] ordered by inclusion posetB :: Int -> Poset [Int]-posetB n = Poset ( powerset [1..n], isSubset )+posetB n = Poset ( powerset [1..n], LS.isSubset )   -- LATTICE OF PARTITIONS OF [1..N] ORDERED BY REFINEMENT@@ -135,15 +138,9 @@ -- LATTICE OF SUBSPACES OF Fq^n  subspaces fq n = [] : concatMap (flatsPG (n-1) fq) [0..n-1]--- the PG(n-1,Fq) is the set of subspaces of fq^n--isZero v = all (==0) v+-- note that flatsPG returns the subspaces as a matrix of row vectors in reduced row echelon form  -- inSpanRE m v returns whether the vector v is in the span of the matrix m, where m is required to be in row echelon form-inSpanRE ((1:xs):bs) (y:ys) = inSpanRE (map tail bs) (if y == 0 then ys else ys <-> y *> xs)-inSpanRE ((0:xs):bs) (y:ys) = if y == 0 then inSpanRE (xs : map tail bs) ys else False-inSpanRE _ ys = isZero ys- isSubspace s1 s2 = all (inSpanRE s2) s1  -- This is the projective geometry PG(n,q)

@@ -51,9 +51,9 @@ 
 -- Triangular graph - van Lint & Wilson p262
 -- http://mathworld.wolfram.com/TriangularGraph.html
-t m = G.to1n $ t' m
+t' m = G.to1n $ t m
 
-t' m | m >= 4 = graph (vs,es) where
+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
@@ -61,9 +61,9 @@ 
 -- Lattice graph - van Lint & Wilson p262
 -- http://mathworld.wolfram.com/LatticeGraph.html
-l2 m = G.to1n $ l2' m
+l2' m = G.to1n $ l2 m
 
-l2' m | m >= 2 = graph (vs,es) where
+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)
@@ -79,9 +79,9 @@ -- CLEBSCH GRAPH
 -- van Lint & Wilson, p263
 
-clebsch = G.to1n clebsch'
+clebsch' = G.to1n clebsch
 
-clebsch' = graph (vs,es) where
+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]
 
@@ -118,9 +118,9 @@ plane +^^ gs = orbit (+^) plane gs
 -- plane +^^ gs = closure [plane] [ +^ g | g <- gs ]
 
-hoffmanSingleton = G.to1n hoffmanSingleton'
+hoffmanSingleton' = G.to1n hoffmanSingleton
 
-hoffmanSingleton' = graph (vs,es) where
+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
@@ -129,7 +129,7 @@ 
 -- induced action of A7 on Hoffman-Singleton graph
 inducedA7 g = fromPairs [(v, v ~^ g) | v <- vs] where
-    vs = vertices hoffmanSingleton'
+    vs = vertices hoffmanSingleton
     (Left h) ~^ g = Left (h +^ g)
     (Right t) ~^ g = Right (t -^ g)
 
@@ -140,9 +140,9 @@ -- van Lint & Wilson p266-7
 -- (also called Sims-Gewirtz graph)
 
-gewirtz = G.to1n gewirtz'
+gewirtz' = G.to1n gewirtz
 
-gewirtz' = graph (vs,es) where
+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.
@@ -154,10 +154,10 @@ 
 data DesignVertex = C | P Integer | B [Integer] deriving (Eq,Ord,Show)
 
-higmanSimsGraph = G.to1n higmanSimsGraph'
+higmanSimsGraph' = G.to1n higmanSimsGraph
 
 -- Cameron & van Lint, p107
-higmanSimsGraph' = graph (vs,es) where
+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]
@@ -170,7 +170,7 @@ -- induced action of M22 on Higman-Sims graph
 inducedM22 g = fromPairs [(v, v ~^ g) | v <- vs] where
     -- G vs _ = higmanSimsGraph'
-    vs = vertices higmanSimsGraph'
+    vs = vertices higmanSimsGraph
     (B b) ~^ g = B (b -^ g)
     (P p) ~^ g = P (p .^ g)
     C ~^ _ = C
@@ -219,9 +219,9 @@                  -- edges between us and its complement are switched
 
 -- Godsil & Royle p259
-schlafli = G.to1n schlafli'
+schlafli' = G.to1n schlafli
 
-schlafli' = graph (vs,es') where
+schlafli = graph (vs,es') where
     g = lineGraph $ k 8
     v:vs = vertices g
     es = edges g
@@ -234,9 +234,9 @@ -- http://people.csse.uwa.edu.au/gordon/constructions/mclaughlin/
 -- http://mathworld.wolfram.com/McLaughlinGraph.html
 
-mcLaughlin = G.to1n mcLaughlin'
+mcLaughlin' = G.to1n mcLaughlin
 
-mcLaughlin' = graph (vs',es') where
+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]

@@ -0,0 +1,282 @@+-- Copyright (c) David Amos, 2011. All rights reserved.++{-# LANGUAGE TupleSections, NoMonomorphismRestriction #-}++-- |A module providing an efficient implementation of the Buchberger algorithm for calculating the (reduced) Groebner basis for an ideal,+-- together with some straightforward applications.+module Math.CommutativeAlgebra.GroebnerBasis where++import Data.List as L+import qualified Data.IntMap as IM+import qualified Data.Set as S++import Math.Core.Utils+import Math.Core.Field+import Math.Algebras.VectorSpace+import Math.Algebras.Structures+import Math.CommutativeAlgebra.Polynomial++-- 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 d = tgcd (lt f) (lt g)+            in (lt g `tdiv` d) *-> f - (lt f `tdiv` d) *-> 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)+++-- Initial, naive version+-- 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' (l:ls) rs = let l' = l %% (rs ++ ls) in+                        if l' == 0 then reduce' ls rs else reduce' ls (toMonic l':rs)+    reduce' [] rs = L.sort rs+-- 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 mcoprime (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 = mlcm (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) `mdivides` l+    gb' gs [] _ = gs++xs ! i = xs !! (i-1) -- in other words, index the list from 1 not 0+++-- Doesn't result in any speedup+gb2a fs = reduce $ gb' fs' (pairs [1..s]) s+    where fs' = IM.fromList $ zip [1..] $ filter (/= 0) fs+          s = IM.size fs'+          gb' gs ((i,j):ps) t =+              let fi = gs IM.! i; fj = gs IM.! j in+              if mcoprime (lm fi) (lm fj) || criterion fi fj+              then gb' gs ps t+              else let h = sPoly fi fj %% IM.elems gs in+                   if h == 0+                   then gb' gs ps t+                   else let t' = t+1+                            gs' = IM.insert t' h gs+                            ps' = ps ++ map (,t') [1..t]+                        in gb' gs' ps' t'+              where criterion fi fj = let l = mlcm (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 IM.! k) `mdivides` l+          gb' gs [] _ = IM.elems gs+++-- version of gb2 where we eliminate pairs as they're created, rather than as they're processed+gb3 fs = reduce $ gb1' [] fs [] 0+    where+    gb1' gs (f:fs) ps t = let ps' = updatePairs gs ps f t+                          in gb1' (gs ++ [f]) fs ps' (t+1)+    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 `mdivides` mlcm (lm (gs!i)) (lm (gs!j)))] +     ++ [ (i,t+1) | (gi,i) <- zip gs [1..t],+                    not (mcoprime (lm gi) (lm f)),+                    not (criterion (mlcm (lm gi) (lm f)) i) ]+        where criterion l i = any (`mdivides` 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")+gb4 fs = reduce $ gb1' [] fs' [] 0+    where fs' = reverse $ L.sort $ filter (/=0) fs+          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 gs ps f t =+              let oldps = [p | p@(l,(i,j)) <- ps, not (lm f `mdivides` l)]+                  newps = sortBy (flip cmpfst)+                          [ (l,(i,t+1)) | (gi,i) <- zip gs [1..t], let l = mlcm (lm gi) (lm f),+                                          not (mcoprime (lm gi) (lm f)),+                                          not (criterion l i) ]+              in mergeBy (flip cmpfst) oldps newps+              where criterion l i = any (`mdivides` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` map snd ps]+++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+++-- 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++-- It is only for Lex ordering that sugar seems to give a substantial improvement+-- gb4 is fine for Grevlex++-- !! can probably get rid of the Ord k requirement in the following++-- |Given a list of polynomials over a field, return a Groebner basis for the ideal generated by the polynomials.+gb :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>+    [Vect k m] -> [Vect k m]+gb fs =+    let fs' = reverse $ 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 gs ps gk k =+        let newps = [let l = mlcm (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,+                       let ((sik,tik),_) = newps ! i, let ((sjk,tjk),_) = newps ! j,+                       not ( (tik `mproperlydivides` tij) && (tjk `mproperlydivides` tij) ) ] -- sloppy variant+            newps' = discard1 [] newps+            newps'' = sortBy (flip cmpSug) $ discard2 [] $ sortBy (flip cmpNormal) newps'+        in mergeBy (flip cmpSug) ps' newps''+        where+            discard1 ls (r@((_sik,tik),(i,_k)):rs) =+                if lm (gs!i) `mcoprime` 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 `mdivides` tjk)) rs+            discard2 ls [] = ls+-- 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 = mdeg h + max (deg f - mdeg (lm f)) (deg g - mdeg (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)+++{-+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+-}++-- OPERATIONS ON IDEALS++memberGB f gs = f %% gs == 0++-- |@memberI f gs@ returns whether f is in the ideal generated by gs+memberI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>+    Vect k m -> [Vect k m] -> Bool+memberI f gs = memberGB f (gb gs)++-- Cox et al, p181+-- |Given ideals I and J, their sum is defined as I+J = {f+g | f \<- I, g \<- J}.+--+-- If fs and gs are generators for I and J, then @sumI fs gs@ returns generators for I+J.+--+-- The geometric interpretation is that the variety of the sum is the intersection of the varieties,+-- ie V(I+J) = V(I) intersect V(J)+sumI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>+     [Vect k m] -> [Vect k m] -> [Vect k m]+sumI fs gs = gb $ fs ++ gs++-- Cox et al, p183+-- |Given ideals I and J, their product I.J is the ideal generated by all products {f.g | f \<- I, g \<- J}.+--+-- If fs and gs are generators for I and J, then @productI fs gs@ returns generators for I.J.+--+-- The geometric interpretation is that the variety of the product is the union of the varieties,+-- ie V(I.J) = V(I) union V(J)+productI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>+     [Vect k m] -> [Vect k m] -> [Vect k m]+productI fs gs = gb [f * g | f <- fs, g <- gs]++-- Cox et al, p185-6+-- |The intersection of ideals I and J is the set of all polynomials which belong to both I and J.+--+-- If fs and gs are generators for I and J, then @intersectI fs gs@ returns generators for the intersection of I and J+--+-- The geometric interpretation is that the variety of the intersection is the union of the varieties,+-- ie V(I intersect J) = V(I) union V(J).+--+-- The reason for prefering the intersection over the product is that the intersection of radical ideals is radical,+-- whereas the product need not be.+intersectI :: (Fractional k, Ord k, Monomial m, Ord m) =>+     [Vect k m] -> [Vect k m] -> [Vect k m]+intersectI fs gs =+    let t = toElimFst $ return $ (mvar "t" :: Glex String)+        hs = map ((t *) . toElimSnd) fs ++ map (((1-t) *) . toElimSnd) gs+    in eliminateFst hs++toElimFst = fmap (\m -> Elim2 m munit)+toElimSnd = fmap (\m -> Elim2 munit m)+isElimFst = (/= munit) . (\(Elim2 m _) -> m) . lm+fromElimSnd = fmap (\(Elim2 _ m) -> m)+eliminateFst = map fromElimSnd . dropWhile isElimFst . gb+++-- Cox et al, p193-4+-- |Given ideals I and J, their quotient is defined as I:J = {f | f \<- R, f.g is in I for all g in J}.+--+-- If fs and gs are generators for I and J, then @quotientI fs gs@ returns generators for I:J.+--+-- The ideal quotient is the algebraic analogue of the Zariski closure of a difference of varieties.+-- V(I:J) contains the Zariski closure of V(I)-V(J), with equality if k is algebraically closed and I is a radical ideal.+quotientI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>+    [Vect k m] -> [Vect k m] -> [Vect k m]+quotientI _ [] = [1]+quotientI fs gs = foldl1 intersectI $ map (quotientP fs) gs+-- quotientI fs gs = foldl intersectI [1] $ map (quotientP fs) gs++quotientP fs g = map ( // g ) $ intersectI fs [g]+    where h // g = let ([u],_) = quotRemMP h [g] in u+

@@ -0,0 +1,368 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances #-}++-- |A module defining the algebra of commutative polynomials over a field k.+-- Polynomials are represented as the free k-vector space with the monomials as basis.+-- +-- A monomial ordering is required to specify how monomials are to be ordered.+-- The Lex, Glex, and Grevlex monomial orders are defined, with the possibility to add others.+module Math.CommutativeAlgebra.Polynomial where++import Math.Core.Field+import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures++-- |In order to work with monomials, we need to be able to multiply them and divide them.+-- Multiplication is defined by the Mon (monoid) class. Division is defined in this class.+-- The functions here are primarily intended for internal use only.+class (Eq m, Show m, Mon m) => Monomial m where+    -- mdivmaybe :: m -> m -> Maybe m+    mdivides :: m -> m -> Bool+    mdiv :: m -> m -> m+    mgcd :: m -> m -> m+    mlcm :: m -> m -> m+    mcoprime :: m -> m -> Bool+    mdeg :: m -> Int++-- mlcm m1 m2 = let m = mgcd m1 m2 in mmult m1 (mdiv m2 m)++-- mdividesne m1 m2 = m1 /= m2 && mdivides m1 m2++mproperlydivides m1 m2 = m1 /= m2 && mdivides m1 m2+++-- |We want to be able to construct monomials over any set of variables that we choose.+-- Although we will often use String as the type of our variables,+-- it is useful to define polymorphic types for monomials.+class MonomialConstructor m where+    mvar :: v -> m v+    mindices :: m v -> [(v,Int)]++-- |@var v@ creates a variable in the vector space of polynomials.+-- For example, if we want to work in Q[x,y,z], we might define:+--+-- > [x,y,z] = map var ["x","y","z"] :: [GlexPoly Q String]+--+-- Notice that, in general, it is necessary to provide a type annotation so that+-- the compiler knows which field and which term order to use.+var :: (Num k, MonomialConstructor m) => v -> Vect k (m v)+var = return . mvar++-- class MonomialOrder m where+--     isGraded :: m -> Bool+++-- MONOMIALS++-- |The underlying implementation of monomials in variables of type v. Most often, we will be interested in MonImpl String,+-- with the variable \"x\" represented by M 1 [(\"x\",1)]. However, other types can be used instead.+--+-- No Ord instance is defined for MonImpl v, so it cannot be used as the basis for a free vector space of polynomials.+-- Instead, several different newtype wrappers are provided, corresponding to different monomial orderings.+data MonImpl v = M Int [(v,Int)] deriving (Eq)+-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons++instance Show v => Show (MonImpl v) where+    show (M _ []) = "1"+    show (M _ xis) = concatMap (\(x,i) -> if i==1 then showVar x else showVar x ++ "^" ++ show i) xis+        where showVar x = filter ( /= '"' ) (show x) -- in case v == String++instance (Ord v) => Mon (MonImpl v) where+    munit = M 0 []+    mmult (M si xis) (M sj yjs) = M (si+sj) $ addmerge xis yjs++instance (Ord v, Show v) => Monomial (MonImpl v) where+{-+    mdivmaybe m1 m2 =+        let m@(M s xis) = mdiv m1 m2+        in if s < 0 || any (< 0) (map snd xis)+           then Nothing+           else Just m+-}+    mdivides (M si xis) (M sj yjs) = si <= sj && mdivides' xis yjs where+        mdivides' ((x,i):xis) ((y,j):yjs) =+            case compare x y of+            LT -> False+            GT -> mdivides' ((x,i):xis) yjs+            EQ -> if i<=j then mdivides' xis yjs else False+        mdivides' [] _ = True+        mdivides' _ [] = False+    mdiv (M si xis) (M sj yjs) = M (si-sj) $ addmerge xis $ map (\(y,j) -> (y,-j)) yjs+    -- we don't check that the result has no negative indices+    mgcd (M _ xis) (M _ yjs) = mgcd' 0 [] xis yjs+        where mgcd' s zks ((x,i):xis) ((y,j):yjs) =+                  case compare x y of+                  LT -> mgcd' s zks xis ((y,j):yjs)+                  GT -> mgcd' s zks ((x,i):xis) yjs+                  EQ -> let k = min i j in mgcd' (s+k) ((x,k):zks) xis yjs+              mgcd' s zks _ _ = M s (reverse zks)+    mlcm (M si xis) (M sj yjs) = mlcm' 0 [] xis yjs+        where mlcm' s zks ((x,i):xis) ((y,j):yjs) =+                  case compare x y of+                  LT -> mlcm' (s+i) ((x,i):zks) xis ((y,j):yjs)+                  GT -> mlcm' (s+j) ((y,j):zks) ((x,i):xis) yjs+                  EQ -> let k = max i j in mlcm' (s+k) ((x,k):zks) xis yjs+              mlcm' s zks xis yjs = let zks' = xis ++ yjs; s' = sum (map snd zks') -- either xis or yjs is null+                                    in M (s+s') (reverse zks ++ zks')+    mcoprime (M _ xis) (M _ yjs) = mcoprime' xis yjs+        where mcoprime' ((x,i):xis) ((y,j):yjs) =+                  case compare x y of+                  LT -> mcoprime' xis ((y,j):yjs)+                  GT -> mcoprime' ((x,i):xis) yjs+                  EQ -> False+              mcoprime' _ _ = True+    -- mcoprime m1 m2 = mgcd m1 m2 == munit+    mdeg (M s _) = s++instance MonomialConstructor MonImpl where+    mvar v = M 1 [(v,1)]+    mindices (M si xis) = xis+++-- LEX ORDER++-- |A type representing monomials with Lex ordering.+--+-- Lex stands for lexicographic ordering.+-- For example, in Lex ordering, monomials up to degree two would be ordered as follows: x^2+xy+xz+x+y^2+yz+y+z^2+z+1.+newtype Lex v = Lex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving++instance Show v => Show (Lex v) where+    show (Lex m) = show m++instance Ord v => Ord (Lex v) where+    compare (Lex (M si xis)) (Lex (M sj yjs)) = compare' xis yjs+        where compare' ((x,i):xis) ((y,j):yjs) =+                  case compare x y of+                  LT -> LT+                  GT -> GT+                  EQ -> case compare i j of+                        LT -> GT+                        GT -> LT+                        EQ -> compare' xis yjs+              compare' [] [] = EQ+              compare' _ [] = LT+              compare' [] _ = GT+        -- unfortunately we can't use the following, because we want [] sorted after everything, not before+        -- compare [(x,-i) | (x,i) <- xis] [(y,-j) | (y,j) <- yjs]++-- instance MonomialOrder Lex where isGraded _ = False++-- |A type representing polynomials with Lex term ordering.+type LexPoly k v = Vect k (Lex v)++-- |@lexvar v@ creates a variable in the algebra of commutative polynomials over Q with Lex term ordering.+-- It is provided as a shortcut, to avoid having to provide a type annotation, as with @var@.+-- For example, the following code creates variables called x, y and z:+--+-- > [x,y,z] = map lexvar ["x","y","z"]+lexvar :: v -> LexPoly Q v+lexvar v = return $ Lex $ M 1 [(v,1)]++instance (Num k, Ord v, Show v) => Algebra k (Lex v) where+    unit x = x *> return munit+    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy+++-- GLEX ORDER++-- |A type representing monomials with Glex ordering.+--+-- Glex stands for graded lexicographic. Thus monomials are ordered first by degree, then by lexicographic order.+-- For example, in Glex ordering, monomials up to degree two would be ordered as follows: x^2+xy+xz+y^2+yz+z^2+x+y+z+1.+newtype Glex v = Glex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving++instance Show v => Show (Glex v) where+    show (Glex m) = show m++instance Ord v => Ord (Glex v) where+    compare (Glex (M si xis)) (Glex (M sj yjs)) =+        compare (-si, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- yjs])++-- instance MonomialOrder Glex where isGraded _ = True++-- |A type representing polynomials with Glex term ordering.+type GlexPoly k v = Vect k (Glex v)++-- |@glexvar v@ creates a variable in the algebra of commutative polynomials over Q with Glex term ordering.+-- It is provided as a shortcut, to avoid having to provide a type annotation, as with @var@.+-- For example, the following code creates variables called x, y and z:+--+-- > [x,y,z] = map glexvar ["x","y","z"]+glexvar :: v -> GlexPoly Q v+glexvar v = return $ Glex $ M 1 [(v,1)]++instance (Num k, Ord v, Show v) => Algebra k (Glex v) where+    unit x = x *> return munit+    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy+++-- GREVLEX ORDER++-- |A type representing monomials with Grevlex ordering.+--+-- Grevlex stands for graded reverse lexicographic. Thus monomials are ordered first by degree, then by reverse lexicographic order.+-- For example, in Grevlex ordering, monomials up to degree two would be ordered as follows: x^2+xy+y^2+xz+yz+z^2+x+y+z+1.+--+-- In general, Grevlex leads to the smallest Groebner bases.+newtype Grevlex v = Grevlex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving++instance Show v => Show (Grevlex v) where+    show (Grevlex m) = show m++instance Ord v => Ord (Grevlex v) where+    compare (Grevlex (M si xis)) (Grevlex (M sj yjs)) =+        compare (-si, reverse xis) (-sj, reverse yjs)++-- instance MonomialOrder Grevlex where isGraded _ = True++-- |A type representing polynomials with Grevlex term ordering.+type GrevlexPoly k v = Vect k (Grevlex v)++-- |@grevlexvar v@ creates a variable in the algebra of commutative polynomials over Q with Grevlex term ordering.+-- It is provided as a shortcut, to avoid having to provide a type annotation, as with @var@.+-- For example, the following code creates variables called x, y and z:+--+-- > [x,y,z] = map grevlexvar ["x","y","z"]+grevlexvar :: v -> GrevlexPoly Q v+grevlexvar v = return $ Grevlex $ M 1 [(v,1)]++instance (Num k, Ord v, Show v) => Algebra k (Grevlex v) where+    unit x = x *> return munit+    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy+++-- ELIMINATION ORDER++data Elim2 a b = Elim2 !a !b deriving (Eq)++instance (Ord a, Ord b) => Ord (Elim2 a b) where+    compare (Elim2 a1 b1) (Elim2 a2 b2) = compare (a1,b1) (a2,b2)++instance (Show a, Show b) => Show (Elim2 a b) where+    show (Elim2 ma mb) = case (show ma, show mb) of+                       ("1","1") -> "1"+                       (ma',"1") -> ma'+                       ("1",mb') -> mb'+                       (ma',mb') -> ma' ++ mb'++instance (Mon a, Mon b) => Mon (Elim2 a b) where+    munit = Elim2 munit munit+    mmult (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mmult a1 a2) (mmult b1 b2)++instance (Monomial a, Monomial b) => Monomial (Elim2 a b) where+    -- mdivmaybe :: Ord v => m v -> m v -> Maybe (m v)+    mdivides (Elim2 a1 b1) (Elim2 a2 b2) = mdivides a1 a2 && mdivides b1 b2+    mdiv (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mdiv a1 a2) (mdiv b1 b2)+    mgcd (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mgcd a1 a2) (mgcd b1 b2)+    mlcm (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mlcm a1 a2) (mlcm b1 b2)+    mcoprime (Elim2 a1 b1) (Elim2 a2 b2) = mcoprime a1 a2 && mcoprime b1 b2+    mdeg (Elim2 a b) = mdeg a + mdeg b++instance (Num k, Ord a, Mon a, Ord b, Mon b) => Algebra k (Elim2 a b) where+    unit x = x *> return munit+    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy+++-- VARIABLE SUBSTITUTION++-- In effect, we have (Num k, Monomial m) => Monad (\v -> Vect k (m v)), with return = var, and (>>=) = bind.+-- However, we can't express this directly in Haskell, firstly because of the Ord v constraint,+-- secondly because Haskell doesn't support type functions.+bind :: (MonomialConstructor m, Num k, Ord a, Show a, Algebra k a) =>+    Vect k (m v) -> (v -> Vect k a) -> Vect k a+V ts `bind` f = sum [c *> product [f x ^ i | (x,i) <- mindices m] | (m, c) <- ts] ++flipbind f = linear (\m -> product [f x ^ i | (x,i) <- mindices m])+++-- DIVISION ALGORITHM FOR POLYNOMIALS++lt (V (t:ts)) = t -- leading term+lm = fst . lt     -- leading monomial+lc = snd . lt     -- leading coefficient++-- deg :: (Num k, Monomial m, MonomialOrder m) => Vect k m -> Int+deg (V []) = -1+deg f = maximum $ [mdeg m | (m,c) <- terms f]+{-+deg f | isGraded (lm f) = mdeg (lm f)+      | otherwise = maximum $ [mdeg m | (m,c) <- terms f]+-}+-- the true degree of the polynomial, not the degree of the leading term+-- required for sugar strategy when computing Groebner basis++toMonic 0 = 0+toMonic f = (1 / lc f) *> f++-- tdivmaybe (m1,x1) (m2,x2) = fmap (\m -> (m,x1/x2)) $ mdivmaybe m1 m2++tdivides (m1,x1) (m2,x2) = mdivides m1 m2++tdiv (m1,x1) (m2,x2) = (mdiv m1 m2, x1/x2)++tgcd (m1,_) (m2,_) = (mgcd m1 m2, 1)+-- tlcm (m1,_) (m2,_) = (mlcm m1 m2, 1)++tmult (m,c) (m',c') = (mmult m m',c*c')++infixl 7 *->+t *-> V ts = V $ map (tmult 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 `tdivides` lt h+        then let t = V [lt h `tdiv` lt g]+                 h' = h - t*g+                 u' = u+t+             in quotRemMP' h' (reverse us' ++ u':us, r)+        else divisionStep h (gs,u:us',us,r)+    divisionStep h ([],us',[],r) =+        let (lth,h') = splitlt h+        in quotRemMP' h' (reverse us', r+lth)+    splitlt (V (t:ts)) = (V [t], V ts)++{-+-- version using mdivmaybe - surprisingly, this seems to be slower+quotRemMP2 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) =+        case tdivmaybe (lt h) (lt g) of+        Nothing -> divisionStep h (gs,u:us',us,r)+        Just t -> let h' = h - V [t] * g+                      u' = u + V [t]+                  in quotRemMP' h' (reverse us' ++ u':us, r)+    divisionStep h ([],us',[],r) =+        let (lth,h') = splitlt h+        in quotRemMP' h' (reverse us', r+lth)+    splitlt (V (t:ts)) = (V [t], V ts)+-}++infixl 7 %%++-- |@f %% gs@ is the reduction of a polynomial f with respect to a list of polynomials gs.+-- In the case where the gs are a Groebner basis for an ideal I,+-- then @f %% gs@ is the equivalence class representative of f in R/I,+-- and is zero if and only if f is in I.+(%%) :: (Fractional k, Monomial m, Ord m, Algebra k m) =>+     Vect k m -> [Vect k m] -> Vect k m+f %% gs = r where (_,r) = quotRemMP f gs+++-- |As a convenience, a partial instance of Fractional is defined for polynomials.+-- The instance is well-defined only for scalars, and gives an error if used on other values.+-- The purpose of this is to allow entry of fractional scalars, in expressions such as @x/2@.+-- On the other hand, an expression such as @2/x@ will return an error.+instance (Fractional k, Monomial m, Ord m, Algebra k m) => Fractional (Vect k m) where+    recip (V [(m,c)]) | m == munit = V [(m,1/c)]+                      | otherwise = error "Polynomial recip: only defined for scalars"+    fromRational x = V [(munit, fromRational x)]

@@ -0,0 +1,466 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++{-# LANGUAGE GeneralizedNewtypeDeriving #-}++-- |A module defining the field Q of rationals and the small finite fields F2, F3, F4, F5, F7, F8, F9, F11, F13, F16, F17, F19, F23, F25.+--+-- Given a prime power q, Fq is the type representing elements of the field (eg @F4@),+-- fq is a list of the elements of the field, beginning 0,1,... (eg @f4@),+-- and for prime power fields, aq is a primitive element, which generates the multiplicative group (eg @a4@).+--+-- The design philosophy is that fq, the list of elements, represents the field.+-- Thus, many functions elsewhere in the library expect to take fq as an argument,+-- telling them which field to work over.+module Math.Core.Field where++import Data.Ratio+import Data.Bits+import Data.List as L++-- |Q is just the rationals, but with a better show function than the Prelude version+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+++-- The following implementations of the prime fields are only slightly faster than the versions in Math.Algebra.Field.Base++-- |F2 is a type for the finite field with 2 elements+newtype F2 = F2 Int deriving (Eq,Ord)++instance Show F2 where+    show (F2 x) = show x++instance Num F2 where+    F2 x + F2 y = F2 $ (x+y) .&. 1 -- `mod` 2+    negate x = x+    F2 x * F2 y = F2 $ x*y+    fromInteger n = F2 $ fromInteger n `mod` 2+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F2 where+    recip (F2 0) = error "F2.recip 0"+    recip (F2 1) = F2 1+    fromRational _ = error "F2.fromRational: not well defined"++-- |f2 is a list of the elements of F2+f2 :: [F2]+f2 = map fromInteger [0..1] -- :: [F2]+++-- |F3 is a type for the finite field with 3 elements+newtype F3 = F3 Int deriving (Eq,Ord)++instance Show F3 where+    show (F3 x) = show x++instance Num F3 where+    F3 x + F3 y = F3 $ (x+y) `mod` 3+    negate (F3 0) = F3 0+    negate (F3 x) = F3 $ 3 - x+    F3 x * F3 y = F3 $ (x*y) `mod` 3+    fromInteger n = F3 $ fromInteger n `mod` 3+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F3 where+    recip (F3 0) = error "F3.recip 0"+    recip (F3 x) = F3 x+    fromRational _ = error "F3.fromRational: not well defined"++-- |f3 is a list of the elements of F3+f3 :: [F3]+f3 = map fromInteger [0..2] -- :: [F3]+++-- |F5 is a type for the finite field with 5 elements+newtype F5 = F5 Int deriving (Eq,Ord)++instance Show F5 where+    show (F5 x) = show x++instance Num F5 where+    F5 x + F5 y = F5 $ (x+y) `mod` 5+    negate (F5 0) = F5 0+    negate (F5 x) = F5 $ 5 - x+    F5 x * F5 y = F5 $ (x*y) `mod` 5+    fromInteger n = F5 $ fromInteger n `mod` 5+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F5 where+    recip (F5 0) = error "F5.recip 0"+    recip (F5 x) = F5 $ (x^3) `mod` 5+    fromRational _ = error "F5.fromRational: not well defined"++-- |f5 is a list of the elements of F5+f5 :: [F5]+f5 = map fromInteger [0..4]+++-- |F7 is a type for the finite field with 7 elements+newtype F7 = F7 Int deriving (Eq,Ord)++instance Show F7 where+    show (F7 x) = show x++instance Num F7 where+    F7 x + F7 y = F7 $ (x+y) `mod` 7+    negate (F7 0) = F7 0+    negate (F7 x) = F7 $ 7 - x+    F7 x * F7 y = F7 $ (x*y) `mod` 7+    fromInteger n = F7 $ fromInteger n `mod` 7+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F7 where+    recip (F7 0) = error "F7.recip 0"+    recip (F7 x) = F7 $ (x^5) `mod` 7+    fromRational _ = error "F7.fromRational: not well defined"++-- |f7 is a list of the elements of F7+f7 :: [F7]+f7 = map fromInteger [0..6]+++-- |F11 is a type for the finite field with 11 elements+newtype F11 = F11 Int deriving (Eq,Ord)++instance Show F11 where+    show (F11 x) = show x++instance Num F11 where+    F11 x + F11 y = F11 $ (x+y) `mod` 11+    negate (F11 0) = F11 0+    negate (F11 x) = F11 $ 11 - x+    F11 x * F11 y = F11 $ (x*y) `mod` 11+    fromInteger n = F11 $ fromInteger n `mod` 11+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F11 where+    recip (F11 0) = error "F11.recip 0"+    recip (F11 x) = F11 $ (x^9) `mod` 11+    fromRational _ = error "F11.fromRational: not well defined"++-- |f11 is a list of the elements of F11+f11 :: [F11]+f11 = map fromInteger [0..10]+++-- |F13 is a type for the finite field with 13 elements+newtype F13 = F13 Int deriving (Eq,Ord)++instance Show F13 where+    show (F13 x) = show x++instance Num F13 where+    F13 x + F13 y = F13 $ (x+y) `mod` 13+    negate (F13 0) = F13 0+    negate (F13 x) = F13 $ 13 - x+    F13 x * F13 y = F13 $ (x*y) `mod` 13+    fromInteger n = F13 $ fromInteger n `mod` 13+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F13 where+    recip (F13 0) = error "F13.recip 0"+    recip (F13 x) = F13 $ (x5*x5*x) `mod` 13 where x5 = x^5 `mod` 13 -- 12^11 would overflow Int+    fromRational _ = error "F13.fromRational: not well defined"++-- |f13 is a list of the elements of F13+f13 :: [F13]+f13 = map fromInteger [0..12]+++-- |F17 is a type for the finite field with 17 elements+newtype F17 = F17 Int deriving (Eq,Ord)++instance Show F17 where+    show (F17 x) = show x++instance Num F17 where+    F17 x + F17 y = F17 $ (x+y) `mod` 17+    negate (F17 0) = F17 0+    negate (F17 x) = F17 $ 17 - x+    F17 x * F17 y = F17 $ (x*y) `mod` 17+    fromInteger n = F17 $ fromInteger n `mod` 17+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F17 where+    recip (F17 0) = error "F17.recip 0"+    recip (F17 x) = F17 $ (x5^3) `mod` 17 where x5 = x^5 `mod` 17 -- 16^15 would overflow Int+    fromRational _ = error "F17.fromRational: not well defined"++-- |f17 is a list of the elements of F17+f17 :: [F17]+f17 = map fromInteger [0..16]+++-- |F19 is a type for the finite field with 19 elements+newtype F19 = F19 Int deriving (Eq,Ord)++instance Show F19 where+    show (F19 x) = show x++instance Num F19 where+    F19 x + F19 y = F19 $ (x+y) `mod` 19+    negate (F19 0) = F19 0+    negate (F19 x) = F19 $ 19 - x+    F19 x * F19 y = F19 $ (x*y) `mod` 19+    fromInteger n = F19 $ fromInteger n `mod` 19+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F19 where+    recip (F19 0) = error "F17.recip 0"+    recip (F19 x) = F19 $ (x4^4*x) `mod` 19 where x4 = x^4 `mod` 19 -- 18^17 would overflow Int+    fromRational _ = error "F19.fromRational: not well defined"++-- |f19 is a list of the elements of F19+f19 :: [F19]+f19 = map fromInteger [0..18]+++-- |F23 is a type for the finite field with 23 elements+newtype F23 = F23 Int deriving (Eq,Ord)++instance Show F23 where+    show (F23 x) = show x++instance Num F23 where+    F23 x + F23 y = F23 $ (x+y) `mod` 23+    negate (F23 0) = F23 0+    negate (F23 x) = F23 $ 23 - x+    F23 x * F23 y = F23 $ (x*y) `mod` 23+    fromInteger n = F23 $ fromInteger n `mod` 23+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F23 where+    recip (F23 0) = error "F23.recip 0"+    recip (F23 x) = F23 $ (x5^4*x) `mod` 23 where x5 = x^5 `mod` 23 -- 22^21 would overflow Int+    fromRational _ = error "F23.fromRational: not well defined"++-- |f23 is a list of the elements of F23+f23 :: [F23]+f23 = map fromInteger [0..22]+++-- The following implementations of the prime power fields are significantly faster than the versions in Math.Algebra.Field.Extension++-- |F4 is a type for the finite field with 4 elements.+-- F4 is represented as the extension of F2 by an element a4 satisfying x^2+x+1 = 0+newtype F4 = F4 Int deriving (Eq,Ord)++instance Show F4 where+    show (F4 0x00) = "0"+    show (F4 0x01) = "1"+    show (F4 0x10) = "a4"+    show (F4 0x11) = "a4+1" -- == a4^2++-- |a4 is a primitive element for F4 as an extension over F2. a4 satisfies x^2+x+1 = 0.+a4 :: F4+a4 = F4 0x10++instance Num F4 where+    F4 x + F4 y = F4 $ (x+y) .&. 0x11+    negate x = x+    F4 x * F4 y = let z = x*y in+                  if z `testBit` 8+                  then F4 ((z + 0x11) .&. 0x11) -- this is replacing x^2 by x+1+                  else F4 z+    fromInteger n = F4 $ fromInteger n .&. 1+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F4 where+    recip (F4 0) = error "F4.recip 0"+    recip (F4 1) = F4 1+    recip (F4 x) = F4 (x `xor` 1)+    fromRational _ = error "F4.fromRational: not well defined"++-- |f4 is a list of the elements of F4+f4 :: [F4]+f4 = L.sort $ 0 : powers a4++powers x | x /= 0 = 1 : takeWhile (/=1) (iterate (*x) x)+++-- |F8 is a type for the finite field with 8 elements.+-- F8 is represented as the extension of F2 by an element a8 satisfying x^3+x+1 = 0+newtype F8 = F8 Int deriving (Eq,Ord)++instance Show F8 where+    show (F8 0x0) = "0"+    show (F8 0x1) = "1"+    show (F8 0x10) = "a8"+    show (F8 0x11) = "a8+1"+    show (F8 0x100) = "a8^2"+    show (F8 0x101) = "a8^2+1"+    show (F8 0x110) = "a8^2+a8"+    show (F8 0x111) = "a8^2+a8+1"++-- |a8 is a primitive element for F8 as an extension over F2. a8 satisfies x^3+x+1 = 0.+a8 :: F8+a8 = F8 0x10++instance Num F8 where+    F8 x + F8 y = F8 $ (x+y) .&. 0x111+    negate x = x+    F8 x * F8 y = F8 $ ((z43 `shiftR` 8) + (z43 `shiftR` 12) + z) .&. 0x111+        where z = x*y; z43 = z .&. 0xff000; -- z210 = z .&. 0xfff+        -- Explanation: We are making the substitution x^3 = x+1, x^4 = x^2+x+    fromInteger n = F8 $ fromInteger n .&. 0x1+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F8 where+    recip (F8 0) = error "F8.recip 0"+    recip x = x^6+    fromRational _ = error "F8.fromRational: not well defined"++-- |f8 is a list of the elements of F8+f8 :: [F8]+f8 = L.sort $ 0 : powers a8+++-- |F9 is a type for the finite field with 9 elements.+-- F9 is represented as the extension of F3 by an element a9 satisfying x^2+2x+2 = 0+newtype F9 = F9 Int deriving (Eq,Ord)++instance Show F9 where+    show (F9 0x00) = "0"+    show (F9 0x01) = "1"+    show (F9 0x02) = "2"+    show (F9 0x100) = "a9"+    show (F9 0x101) = "a9+1"+    show (F9 0x102) = "a9+2"+    show (F9 0x200) = "2a9"+    show (F9 0x201) = "2a9+1"+    show (F9 0x202) = "2a9+2"++-- |a9 is a primitive element for F9 as an extension over F3. a9 satisfies x^2+2x+2 = 0.+a9 :: F9+a9 = F9 0x100++instance Num F9 where+    F9 x + F9 y = F9 $ z1 + z0+        where z = x+y; z1 = (z .&. 0xff00) `mod` 0x300; z0 = (z .&. 0xff) `mod` 3+    negate (F9 x) = F9 $ z1 + z0+        where z = 0x303 - x; z1 = (z .&. 0xff00) `mod` 0x300; z0 = (z .&. 0xff) `mod` 3+    F9 x * F9 y = F9 $ ((z2 + z1) `mod` 0x300) + ((z2 + z0) `mod` 3) +        where z = x*y; z2 = z .&. 0xff0000; z1 = z .&. 0xff00; z0 = z .&. 0xff+        -- Explanation: We are substituting x^2 = x+1.+        -- We could do z2 `shiftR` 8 and z2 `shiftR` 16+        -- However, because 0x100 `mod` 3 == 1, we don't need to +    fromInteger n = F9 $ fromInteger n `mod` 3+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F9 where+    recip (F9 0) = error "F9.recip 0"+    recip x = x^7+    fromRational _ = error "F9.fromRational: not well defined"++-- |f9 is a list of the elements of F9+f9 :: [F9]+f9 = L.sort $ 0 : powers a9+++-- |F16 is a type for the finite field with 16 elements.+-- F16 is represented as the extension of F2 by an element a16 satisfying x^4+x+1 = 0+newtype F16 = F16 Int deriving (Eq,Ord)++instance Show F16 where+    show (F16 0x0) = "0"+    show (F16 0x1) = "1"+    show (F16 0x10) = "a16"+    show (F16 0x11) = "a16+1"+    show (F16 0x100) = "a16^2"+    show (F16 0x101) = "a16^2+1"+    show (F16 0x110) = "a16^2+a16"+    show (F16 0x111) = "a16^2+a16+1"+    show (F16 0x1000) = "a16^3"+    show (F16 0x1001) = "a16^3+1"+    show (F16 0x1010) = "a16^3+a16"+    show (F16 0x1011) = "a16^3+a16+1"+    show (F16 0x1100) = "a16^3+a16^2"+    show (F16 0x1101) = "a16^3+a16^2+1"+    show (F16 0x1110) = "a16^3+a16^2+a16"+    show (F16 0x1111) = "a16^3+a16^2+a16+1"++-- |a16 is a primitive element for F16 as an extension over F2. a16 satisfies x^4+x+1 = 0.+a16 :: F16+a16 = F16 0x10++instance Num F16 where+    F16 x + F16 y = F16 $ (x+y) .&. 0x1111+    negate x = x+    F16 x * F16 y = F16 $ ((z654 `shiftR` 12) + (z654 `shiftR` 16) + z) .&. 0x1111+        where z = x*y; z654 = z .&. 0xfff0000; -- z3210 = z .&. 0xffff+        -- Explanation: We are making the substitution x^4 = x+1 (and also for x^5, x^6)+    fromInteger n = F16 $ fromInteger n .&. 0x1+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F16 where+    recip (F16 0) = error "F16.recip 0"+    recip x = x^14+    fromRational _ = error "F16.fromRational: not well defined"++-- |f16 is a list of the elements of F16+f16 :: [F16]+f16 = L.sort $ 0 : powers a16+++-- |F25 is a type for the finite field with 25 elements.+-- F25 is represented as the extension of F5 by an element a25 satisfying x^2+4x+2 = 0+newtype F25 = F25 Int deriving (Eq,Ord)++instance Show F25 where+    show (F25 x) = case ( (x .&. 0xff00) `shiftR` 8, x .&. 0xff ) of+                   (0,x0) -> show x0+                   (1,0) -> "a25"+                   (1,x0) -> "a25+" ++ show x0+                   (x1,0) -> show x1 ++ "a25"+                   (x1,x0) -> show x1 ++ "a25+" ++ show x0++-- |a25 is a primitive element for F25 as an extension over F5. a25 satisfies x^2+4x+2 = 0.+a25 :: F25+a25 = F25 0x100++instance Num F25 where+    F25 x + F25 y = F25 $ z1 + z0+        where z = x+y; z1 = (z .&. 0xff00) `mod` 0x500; z0 = (z .&. 0xff) `mod` 5+    negate (F25 x) = F25 $ z1 + z0+        where z = 0x505 - x; z1 = (z .&. 0xff00) `mod` 0x500; z0 = (z .&. 0xff) `mod` 5+    F25 x * F25 y = F25 $ ((z2 + z1) `mod` 0x500) + ((3*z2 + z0) `mod` 5) +        where z = x*y; z2 = z .&. 0xff0000; z1 = z .&. 0xff00; z0 = z .&. 0xff+        -- Explanation: We are substituting x^2 = x+3.+        -- We could do z2 `shiftR` 8 and z2 `shiftR` 16+        -- However, because 0x100 `mod` 5 == 1, we don't need to +    fromInteger n = F25 $ fromInteger n `mod` 5+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"++instance Fractional F25 where+    recip (F25 0) = error "F25.recip 0"+    recip x = x^23+    fromRational _ = error "F25.fromRational: not well defined"++-- |f25 is a list of the elements of F25+f25 :: [F25]+f25 = L.sort $ 0 : powers a25++

@@ -0,0 +1,64 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++{-# LANGUAGE NoMonomorphismRestriction, TupleSections #-}++-- |A module of simple utility functions which are used throughout the rest of the library+module Math.Core.Utils where++import Data.List as L+import qualified Data.Set as S+++toSet = S.toList . S.fromList++pairs (x:xs) = map (x,) xs ++ pairs xs+pairs [] = []++ordpair x y | x < y     = (x,y)+            | otherwise = (y,x)+++-- fold a comparison operator through a list+foldcmpl p (x1:x2:xs) = p x1 x2 && foldcmpl p (x2:xs)+foldcmpl _ _ = True+-- This can be expressed as a pure fold:+-- foldcmpl cmp (x:xs) = snd $ foldl (\(bool,x') x -> (bool && cmp x' x, x)) (True,x)+-- foldcmpl _ [] = True+-- However, that is less efficient, as we can't abort as soon as we fail+-- (What about using the Maybe monad?)++-- for use with L.sortBy+cmpfst x y = compare (fst x) (fst y)++-- for use with L.groupBy+eqfst x y = (==) (fst x) (fst y)+++fromBase b xs = foldl' (\n x -> n * b + x) 0 xs++-- |Given a set @xs@, represented as an ordered list, @powersetdfs xs@ returns the list of all subsets of xs, in lex order+powersetdfs :: [a] -> [[a]]+powersetdfs xs = map reverse $ dfs [ ([],xs) ]+    where dfs ( (ls,rs) : nodes ) = ls : dfs (successors (ls,rs) ++ nodes)+          dfs [] = []+          successors (ls,rs) = [ (r:ls, rs') | r:rs' <- L.tails rs ]++-- |Given a set @xs@, represented as an ordered list, @powersetbfs xs@ returns the list of all subsets of xs, in shortlex order+powersetbfs :: [a] -> [[a]]+powersetbfs xs = map reverse $ bfs [ ([],xs) ]+    where bfs ( (ls,rs) : nodes ) = ls : bfs ( nodes ++ successors (ls,rs) )+          bfs [] = []+          successors (ls,rs) = [ (r:ls, rs') | r:rs' <- L.tails rs ]+++-- |Given a positive integer @k@, and a set @xs@, represented as a list,+-- @combinationsOf k xs@ returns all k-element subsets of xs.+-- The result will be in lex order, relative to the order of the xs.+combinationsOf :: Int -> [a] -> [[a]]+combinationsOf 0 _ = [[]]+combinationsOf _ [] = []+combinationsOf k (x:xs) | k > 0 = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs++-- |@choose n k@ is the number of ways of choosing k distinct elements from an n-set+choose :: (Integral a) => a -> a -> a+choose n k = product [n-k+1..n] `div` product [1..k]

@@ -0,0 +1,148 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++module Math.Test.TCombinatorics.TMatroid where++import Test.HUnit++import Math.Combinatorics.Matroid++import Math.Core.Field hiding (f7)+-- import Math.Algebra.Field.Base hiding (f7)+import Math.Core.Utils (combinationsOf)++testlistMatroid = TestList [+    testlistKnownEquivalent,+    testlistFromCircuits,+    testlistFromBases,+    testlistFromRankfun,+    testlistFromClosure,+    testlistFromFlats,+    testlistFromHyperplanes,+    testlistRepresentable+    ]+++elts = elements+++-- Oxley p8+ex112 = vectorMatroid' [[1,0],[0,1],[0,0],[1,0],[1,1::Q]]++-- Oxley p11 - leads to same matroid as ex112+ex118 = cycleMatroid [ [1,2],[1,3],[3,3],[1,2],[2,3]]++-- Oxley p19 - two non-isomorphic graphs giving rise to the same matroid+fig13a = cycleMatroid [ [1,1],[2,3],[2,4],[2,5],[6,7]]++fig13b = cycleMatroid [ [1,1],[1,2],[2,3],[3,4],[4,5]]++-- Oxley p34, fig 1.8+ex154 = fromGeoRep [[3]] [[1,4]] [[1,2,5],[2,4,5]] [[1..5]]++-- Oxley p46+ex163a = cycleMatroid [ [1,2],[1,2],[2,3],[2,3],[3,4],[3,4],[1,4] ]+ex163b = transversalMatroid [1..7] [ [1,2,7],[3,4,7],[5,6,7] ]+++{-+-- not currently used+ex16 = affineMatroid [ [0,0],[0,1],[0,2],[1,0],[1,1],[2,0::Q]]++ex17 = affineMatroid [ [0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,1,0::Q]]++-}++++testcaseEquivalent desc m1 m2 = TestCase (assertEqual desc m1 m2)++testlistKnownEquivalent = TestList [+    testcaseEquivalent "ex112=118" ex112 ex118,+    testcaseEquivalent "fig13a=fig13b" fig13a fig13b,+    testcaseEquivalent "ex112=ex154" ex112 ex154,+    testcaseEquivalent "desargues = M(K5)" desargues (cycleMatroid (combinationsOf 2 [1..5])), -- Oxley p39+    testcaseEquivalent "ex163a=ex163b" ex163a ex163b+    ]+++exampleList = [+    ("ex112", ex112),+    ("fig13a", fig13a),+    ("f7", f7),+    ("pappus",pappus),+    ("nonPappus",nonPappus),+    ("desargues",desargues),+    ("v8",v8)+    ]+-- if we're really just going to use the same examples for each test, should use map++testcaseFromCircuits (desc, m) = TestCase (assertEqual ("fromCircuits " ++ desc) m (fromCircuits (elts m) (circuits m)))++testlistFromCircuits = TestList $ map testcaseFromCircuits exampleList+++testcaseFromBases (desc, m) = TestCase (assertEqual ("fromBases " ++ desc) m (fromBases (elts m) (bases m)))++testlistFromBases = TestList $ map testcaseFromBases exampleList+++testcaseFromRankfun (desc, m) = TestCase (assertEqual ("fromRankfun " ++ desc) m (fromRankfun (elts m) (rankfun m)))++testlistFromRankfun = TestList $ map testcaseFromRankfun exampleList++ +testcaseFromClosure (desc, m) = TestCase (assertEqual ("fromClosure " ++ desc) m (fromClosure (elts m) (closure m)))++testlistFromClosure = TestList $ map testcaseFromClosure exampleList+++testcaseFromFlats (desc, m) = TestCase (assertEqual ("fromFlats " ++ desc) m (fromFlats (flats m)))++testlistFromFlats = TestList $ map testcaseFromFlats exampleList+++testcaseFromHyperplanes (desc, m) = TestCase (assertEqual ("fromHyperplanes " ++ desc) m (fromHyperplanes (elts m) (hyperplanes m)))++testlistFromHyperplanes = TestList $ map testcaseFromHyperplanes exampleList+++testcaseIsRepresentable desc fq m = TestCase (assertBool ("isRepresentable " ++ desc) (isRepresentable fq m))+testcaseNotRepresentable desc fq m = TestCase (assertBool ("notRepresentable " ++ desc) (not (isRepresentable fq m)))++testlistRepresentable = TestList [+    testcaseNotRepresentable "f2 v8" f2 v8, -- Oxley p84+    testcaseNotRepresentable "f3 v8" f3 v8,+    testcaseNotRepresentable "f4 v8" f4 v8,+    testcaseNotRepresentable "f5 v8" f5 v8,++    testcaseIsRepresentable "f2 f7" f2 f7, -- Oxley p187+    testcaseNotRepresentable "f3 f7" f3 f7,+    testcaseIsRepresentable "f4 f7" f4 f7,+    testcaseNotRepresentable "f5 f7" f5 f7,+    testcaseNotRepresentable "f2 f7m" f2 f7m,+    testcaseIsRepresentable "f3 f7m" f3 f7m,+    testcaseNotRepresentable "f4 f7m" f4 f7m,+    testcaseIsRepresentable "f5 f7m" f5 f7m,++    testcaseNotRepresentable "f2 p8" f2 p8, -- Oxley p189-90+    testcaseIsRepresentable "f3 p8" f3 p8,+    testcaseNotRepresentable "f4 p8" f4 p8,+    testcaseIsRepresentable "f5 p8" f5 p8,+    testcaseNotRepresentable "f2 p8m" f2 p8m,+    testcaseNotRepresentable "f3 p8m" f3 p8m,+    testcaseIsRepresentable "f4 p8m" f4 p8m,+    testcaseIsRepresentable "f5 p8m" f5 p8m,+    testcaseNotRepresentable "f2 p8mm" f2 p8mm,+    testcaseNotRepresentable "f3 p8mm" f3 p8mm,+    testcaseNotRepresentable "f4 p8mm" f4 p8mm,+    testcaseIsRepresentable "f5 p8mm" f5 p8mm,++    testcaseNotRepresentable "f2 u24" f2 (u 2 4), -- Oxley p193+    testcaseNotRepresentable "f3 u25" f3 (u 2 5),+    testcaseNotRepresentable "f3 u35" f3 (u 3 5),+    testcaseNotRepresentable "f4 u26" f4 (u 2 6),+    testcaseNotRepresentable "f4 u46" f4 (u 4 6),+    testcaseNotRepresentable "f5 u27" f5 (u 2 7),+    testcaseNotRepresentable "f5 u57" f5 (u 5 7)+    -- Oxley does mention other excluded minors for other fields, but we have to stop somewhere+    ]

@@ -0,0 +1,75 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++module Math.Test.TCommutativeAlgebra.TGroebnerBasis where++import Test.HUnit++import Math.Core.Field+import Math.Algebras.VectorSpace+import Math.CommutativeAlgebra.Polynomial+import Math.CommutativeAlgebra.GroebnerBasis++-- Sources:+-- Eisenbud, Commutative Algebra with a View Toward Algebraic Geometry+-- [IVA] - Cox, Little, O'Shea: Ideals, Varieties and Algorithms, 2nd ed+-- (Note that I have the 5th printing, which is known to have some misprints)+-- [UAG] - Cox, Little, O'Shea: Using Algebraic Geometry+-- Schenck, Computational Algebraic Geometry++testlistGroebnerBasis = TestList [+    testlistLexGb,+    testlistGlexGb,+    testlistIntersectI,+    testlistQuotientI+    ]+++data Var = X | Y | Z deriving (Eq,Ord)++instance Show Var where+    show X = "x"+    show Y = "y"+    show Z = "z"++[x,y,z] = map grevlexvar [X,Y,Z]+++testcaseGb desc input output = TestCase (assertEqual desc output (gb input))++testlistLexGb =+    let [x,y,z] = map lexvar ["x","y","z"] in+    TestList [+        testcaseGb "Lex <x^2, xy+y^2>" [x^2,x*y+y^2] [x^2,x*y+y^2,y^3], -- Eisenbud p339+        testcaseGb "Lex <x^2+y^2+z^2-1,x^2+z^2-y,x-z>" [x^2+y^2+z^2-1,x^2+z^2-y,x-z] [x-z,y-2*z^2,z^4+1/2*z^2-1/4], -- IVA p93-4+        testcaseGb "Lex <x^2+y^2+z^2-1,x*y*z-1>" [x^2+y^2+z^2-1,x*y*z-1] [x+y^3*z+y*z^3-y*z,y^4*z^2+y^2*z^4-y^2*z^2+1], -- IVA p116+        testcaseGb "Lex <x*y-4,y^2-(x^3-1)>" [x*y-4,y^2-(x^3-1)] [x-1/16*y^4-1/16*y^2,y^5+y^3-64] -- IVA p117, misprint corrected+    ]++testlistGlexGb =+    let [x,y,z] = map glexvar ["x","y","z"] in+    TestList [+        testcaseGb "Glex <x*z-y^2,x^3-z^2>" [x*z-y^2,x^3-z^2] [y^6-z^5,x*y^4-z^4,x^2*y^2-z^3,x^3-z^2,x*z-y^2], -- IVA p93+        testcaseGb "Glex <x^2+y,x*y+x>" [x^2+y,x*y+x] [x^2+y,x*y+x,y^2+y] -- Schenck p54+    ]+++testcaseIntersectI desc i1 i2 iout =+    TestCase (assertEqual ("intersectI " ++ desc) iout (intersectI i1 i2))++testlistIntersectI =+    let [x,y] = map grevlexvar ["x","y"] in+    TestList [+        testcaseIntersectI "[x^2*y] [x*y^2]" [x^2*y] [x*y^2] [x^2*y^2], -- IVA p186+        testcaseIntersectI "IVA 186/2" [(x+y)^4*(x^2+y)^2*(x-5*y)] [(x+y)*(x^2+y)^3*(x+3*y)]+                                          [(x+y)^4*(x^2+y)^3*(x-5*y)*(x+3*y)] -- IVA p186+    ]++testcaseQuotientI desc i j q =+    TestCase (assertEqual ("quotientI " ++ desc) q (i `quotientI` j))++testlistQuotientI =+    let [x,y,z] = map grevlexvar ["x","y","z"] in+    TestList [+        testcaseQuotientI "[x*z, y*z] [z]" [x*z, y*z] [z] [x,y], -- IVA p192+        testcaseQuotientI "[y^2, z^2] [y*z]" [y^2, z^2] [y*z] [y,z] -- Schenck p56 (in passing)+    ]

@@ -0,0 +1,89 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++module Math.Test.TCore.TField where++import Test.QuickCheck+import Math.Core.Field++prop_Field (x,y,z) =+    (x+y)+z == x+(y+z)                          &&  -- associativity of addition+    (x+0) == x && (0+x) == x                    &&  -- identity for addition+    x+(-x) == 0 && (-x)+x == 0                  &&  -- additive inverse+    x+y == y+x                                  &&  -- commutativity of addition+    (x*y)*z == x*(y*z)                          &&  -- associativity of multiplication+    (x*1) == x && (1*x) == x                    &&  -- identity for multiplication+    (x == 0 || (x*(1/x) == 1 && (1/x)*x == 1))  &&  -- multiplicative inverse+    x*y == y*x                                  &&  -- commutativity of multiplication+    x*(y+z) == x*y + x*z && (x+y)*z == x*z + y*z    -- distributivity of multiplication over addition++instance Arbitrary F2 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F3 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F5 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F7 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F11 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F13 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F17 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F19 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F23 where+    arbitrary = do {n <- arbitrary; return (fromInteger n)}++instance Arbitrary F4 where+    arbitrary = do {n <- arbitrary; return (f4 !! mod (fromInteger n) 4)}++instance Arbitrary F8 where+    arbitrary = do {n <- arbitrary; return (f8 !! mod (fromInteger n) 8)}++instance Arbitrary F9 where+    arbitrary = do {n <- arbitrary; return (f9 !! mod (fromInteger n) 9)}++instance Arbitrary F16 where+    arbitrary = do {n <- arbitrary; return (f16 !! mod (fromInteger n) 16)}++instance Arbitrary F25 where+    arbitrary = do {n <- arbitrary; return (f25 !! mod (fromInteger n) 25)}+++test = do putStrLn "Testing F2..."+          quickCheck (prop_Field :: (F2,F2,F2) -> Bool)+          putStrLn "Testing F3..."+          quickCheck (prop_Field :: (F3,F3,F3) -> Bool)+          putStrLn "Testing F5..."+          quickCheck (prop_Field :: (F5,F5,F5) -> Bool)+          putStrLn "Testing F7..."+          quickCheck (prop_Field :: (F7,F7,F7) -> Bool)+          putStrLn "Testing F11..."+          quickCheck (prop_Field :: (F11,F11,F11) -> Bool)+          putStrLn "Testing F13..."+          quickCheck (prop_Field :: (F13,F13,F13) -> Bool)+          putStrLn "Testing F17..."+          quickCheck (prop_Field :: (F17,F17,F17) -> Bool)+          putStrLn "Testing F19..."+          quickCheck (prop_Field :: (F19,F19,F19) -> Bool)+          putStrLn "Testing F23..."+          quickCheck (prop_Field :: (F23,F23,F23) -> Bool)+          putStrLn "Testing F4..."+          quickCheck (prop_Field :: (F4,F4,F4) -> Bool)+          putStrLn "Testing F8..."+          quickCheck (prop_Field :: (F8,F8,F8) -> Bool)+          putStrLn "Testing F9..."+          quickCheck (prop_Field :: (F9,F9,F9) -> Bool)+          putStrLn "Testing F16..."+          quickCheck (prop_Field :: (F16,F16,F16) -> Bool)+          putStrLn "Testing F25..."+          quickCheck (prop_Field :: (F25,F25,F25) -> Bool)