diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal --- a/HaskellForMaths.cabal +++ b/HaskellForMaths.cabal @@ -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, diff --git a/Math/Algebra/Commutative/GBasis.hs b/Math/Algebra/Commutative/GBasis.hs --- a/Math/Algebra/Commutative/GBasis.hs +++ b/Math/Algebra/Commutative/GBasis.hs @@ -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] + + diff --git a/Math/Algebra/LinearAlgebra.hs b/Math/Algebra/LinearAlgebra.hs --- a/Math/Algebra/LinearAlgebra.hs +++ b/Math/Algebra/LinearAlgebra.hs @@ -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 diff --git a/Math/Algebras/Commutative.hs b/Math/Algebras/Commutative.hs --- a/Math/Algebras/Commutative.hs +++ b/Math/Algebras/Commutative.hs @@ -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 diff --git a/Math/Algebras/Structures.hs b/Math/Algebras/Structures.hs --- a/Math/Algebras/Structures.hs +++ b/Math/Algebras/Structures.hs @@ -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 diff --git a/Math/Algebras/TensorAlgebra.hs b/Math/Algebras/TensorAlgebra.hs --- a/Math/Algebras/TensorAlgebra.hs +++ b/Math/Algebras/TensorAlgebra.hs @@ -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 diff --git a/Math/Algebras/VectorSpace.hs b/Math/Algebras/VectorSpace.hs --- a/Math/Algebras/VectorSpace.hs +++ b/Math/Algebras/VectorSpace.hs @@ -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 diff --git a/Math/Combinatorics/FiniteGeometry.hs b/Math/Combinatorics/FiniteGeometry.hs --- a/Math/Combinatorics/FiniteGeometry.hs +++ b/Math/Combinatorics/FiniteGeometry.hs @@ -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 diff --git a/Math/Combinatorics/Graph.hs b/Math/Combinatorics/Graph.hs --- a/Math/Combinatorics/Graph.hs +++ b/Math/Combinatorics/Graph.hs @@ -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 diff --git a/Math/Combinatorics/GraphAuts.hs b/Math/Combinatorics/GraphAuts.hs --- a/Math/Combinatorics/GraphAuts.hs +++ b/Math/Combinatorics/GraphAuts.hs @@ -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 diff --git a/Math/Combinatorics/Matroid.hs b/Math/Combinatorics/Matroid.hs new file mode 100644 --- /dev/null +++ b/Math/Combinatorics/Matroid.hs @@ -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 + diff --git a/Math/Combinatorics/Poset.hs b/Math/Combinatorics/Poset.hs --- a/Math/Combinatorics/Poset.hs +++ b/Math/Combinatorics/Poset.hs @@ -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) diff --git a/Math/Combinatorics/StronglyRegularGraph.hs b/Math/Combinatorics/StronglyRegularGraph.hs --- a/Math/Combinatorics/StronglyRegularGraph.hs +++ b/Math/Combinatorics/StronglyRegularGraph.hs @@ -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] diff --git a/Math/CommutativeAlgebra/GroebnerBasis.hs b/Math/CommutativeAlgebra/GroebnerBasis.hs new file mode 100644 --- /dev/null +++ b/Math/CommutativeAlgebra/GroebnerBasis.hs @@ -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 + diff --git a/Math/CommutativeAlgebra/Polynomial.hs b/Math/CommutativeAlgebra/Polynomial.hs new file mode 100644 --- /dev/null +++ b/Math/CommutativeAlgebra/Polynomial.hs @@ -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)] diff --git a/Math/Core/Field.hs b/Math/Core/Field.hs new file mode 100644 --- /dev/null +++ b/Math/Core/Field.hs @@ -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 + + diff --git a/Math/Core/Utils.hs b/Math/Core/Utils.hs new file mode 100644 --- /dev/null +++ b/Math/Core/Utils.hs @@ -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] diff --git a/Math/Test/TCombinatorics/TMatroid.hs b/Math/Test/TCombinatorics/TMatroid.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TCombinatorics/TMatroid.hs @@ -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 + ] diff --git a/Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs b/Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs @@ -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,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-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+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-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] [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,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) + ] diff --git a/Math/Test/TCore/TField.hs b/Math/Test/TCore/TField.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TCore/TField.hs @@ -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)