diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal --- a/HaskellForMaths.cabal +++ b/HaskellForMaths.cabal @@ -1,5 +1,5 @@ Name: HaskellForMaths - Version: 0.2.0 + Version: 0.2.1 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 @@ -19,6 +19,7 @@ Math/Test/TGraph.hs, Math/Test/TNonCommutativeAlgebra.hs, Math/Test/TPermutationGroup.hs, + Math/Test/TSubquotients.hs, Math/Test/TestAll.hs Library @@ -29,7 +30,8 @@ Math.Algebra.Commutative.MPoly, Math.Algebra.Commutative.GBasis, Math.Algebra.Field.Base, Math.Algebra.Field.Extension, Math.Algebra.Group.PermutationGroup, Math.Algebra.Group.SchreierSims, - Math.Algebra.Group.RandomSchreierSims, Math.Algebra.Group.StringRewriting, + Math.Algebra.Group.RandomSchreierSims, Math.Algebra.Group.Subquotients, + Math.Algebra.Group.StringRewriting, Math.Algebra.NonCommutative.NCPoly, Math.Algebra.NonCommutative.GSBasis, Math.Algebra.NonCommutative.TensorAlgebra, Math.Combinatorics.Graph, Math.Combinatorics.GraphAuts, Math.Combinatorics.StronglyRegularGraph, Math.Combinatorics.Design, Math.Combinatorics.FiniteGeometry, Math.Combinatorics.Hypergraph, diff --git a/Math/Algebra/Group/RandomSchreierSims.hs b/Math/Algebra/Group/RandomSchreierSims.hs --- a/Math/Algebra/Group/RandomSchreierSims.hs +++ b/Math/Algebra/Group/RandomSchreierSims.hs @@ -1,4 +1,4 @@ - +-- Copyright (c) David Amos, 2009. All rights reserved. module Math.Algebra.Group.RandomSchreierSims where @@ -16,16 +16,7 @@ import Math.Algebra.Group.PermutationGroup import Math.Algebra.Group.SchreierSims (sift, cosetRepsGx, ss') -{- --- all the imports below used only for testing -import Math.Algebra.Field.Base -import Math.Algebra.Field.Extension -import Math.Projects.ChevalleyGroup.Classical -import Math.Projects.ChevalleyGroup.Exceptional --} - - testProdRepl = do (r,xs) <- initProdRepl $ _D 10 hs <- replicateM 20 $ nextProdRepl (r,xs) mapM_ print hs @@ -114,8 +105,10 @@ -- recover the base tranversals from the sgs. gs must be an sgs -baseTransversalsSGS gs = [let hs = [h | h <- gs, b <= minsupp h] in (b, cosetRepsGx hs b) | b <- bs] - where bs = toListSet $ concatMap supp gs +-- baseTransversalsSGS gs = [let hs = [h | h <- gs, b <= minsupp h] in (b, cosetRepsGx hs b) | b <- bs] +baseTransversalsSGS gs = [let hs = filter ( (b <=) . minsupp ) gs in (b, cosetRepsGx hs b) | b <- bs] + where bs = toListSet $ map minsupp gs + -- where bs = toListSet $ concatMap supp gs -- |Given a strong generating set gs, isMemberSGS gs is a membership test for the group isMemberSGS :: (Ord a, Show a) => [Permutation a] -> Permutation a -> Bool diff --git a/Math/Algebra/Group/SchreierSims.hs b/Math/Algebra/Group/SchreierSims.hs --- a/Math/Algebra/Group/SchreierSims.hs +++ b/Math/Algebra/Group/SchreierSims.hs @@ -12,7 +12,7 @@ -- COSET REPRESENTATIVES FOR STABILISER OF A POINT --- Given a group G = , and a point x, find coset representatives for Gx (stabiliser of x) in G +-- Given a group G = , and a point x, find (right) coset representatives for Gx (stabiliser of x) in G -- In other words, for each x' in the orbit of x under G, we find a g <- G taking x to x' -- The code is similar to the code for calculating orbits, but modified to keep track of the group elements that we used to get there cosetRepsGx gs x = cosetRepsGx' gs M.empty (M.singleton x 1) where @@ -40,15 +40,18 @@ -- SCHREIER-SIMS ALGORITHM -sift bts g = sift' bts g where - sift' _ 1 = Nothing - sift' ((b,t):bts) g = case M.lookup (b .^ g) t of - Nothing -> Just g - -- Nothing -> sift' bts g -- if we allow empty levels - Just h -> sift' bts (g * inverse h) - sift' [] g = if g == 1 then Nothing else Just g +-- Given a list of right transversals for a stabiliser chain, sift a group element through it +-- Note, this version assumes the base is non-redundant +sift _ 1 = Nothing +sift ((b,t):bts) g = case M.lookup (b .^ g) t of + Nothing -> Just g + -- Nothing -> sift bts g -- if we allow redundant levels + Just h -> sift bts (g * inverse h) +sift [] g = Just g -- g == 1 case already caught above -findBase gs = minimum $ concatMap supp gs + +-- findBase gs = minimum $ concatMap supp gs +findBase gs = minimum $ map minsupp gs {- -- Find base and strong generating set using Schreier-Sims algorithm bsgs gs | all (/= 1) gs = map fst $ ss [newLevel gs] [] diff --git a/Math/Algebra/Group/Subquotients.hs b/Math/Algebra/Group/Subquotients.hs new file mode 100644 --- /dev/null +++ b/Math/Algebra/Group/Subquotients.hs @@ -0,0 +1,196 @@ +-- Copyright (c) David Amos, 2009. All rights reserved. + +{-# LANGUAGE NoMonomorphismRestriction #-} +-- Because unRight defined point-free + +module Math.Algebra.Group.Subquotients where + +import qualified Data.List as L +import qualified Data.Map as M + +import Math.Common.ListSet +import Math.Algebra.Group.PermutationGroup hiding (ptStab, normalClosure) +import Math.Algebra.Group.SchreierSims (cosetRepsGx) +import Math.Algebra.Group.RandomSchreierSims + + +-- Source: Seress, Permutation Group Algorithms + + +isLeft (Left _) = True +isLeft (Right _) = False + +isRight (Right _) = True +isRight (Left _) = False + + +unRight = fromPairs . map (\(Right a, Right b) -> (a,b)) . toPairs + +restrictLeft g = fromPairs [(a,b) | (Left a, Left b) <- toPairs g] +-- note that this is doing a filter - taking only the left part of the action - and a map, unLefting + + +-- pointwise stabiliser of xs +ptStab gs delta = map unRight $ dropWhile (isLeft . minsupp) $ sgs gs' where + gs' = [ (fromPairs . map (\(a,b) -> (lr a, lr b)) . toPairs) g | g <- gs] + lr x = if x `elem` delta then Left x else Right x + + +{- +-- !! NEXT TWO FUNCTIONS NOT TESTED +-- Need some meaningful examples of homomorphisms +-- eg Sn -> Sym(k-subsets of n) +-- restrict to a transitive constituent +-- blocks + +-- Given generators gs for a group G, and f : G -> H a homomorphism, +-- return the "semi-diagonal" subgroup [(f g, g) | g <- gs] of f(G) * G +homomorphismConstruction :: (Ord a, Ord b) => [Permutation a] -> (Permutation a -> Permutation b) -> [Permutation (Either b a)] +homomorphismConstruction gs f = [lift g | g <- gs] where + lift g = fromPairs $ [(Right x, Right y) | (x,y) <- toPairs g] ++ [(Left x', Left y') | (x',y') <- toPairs (f g)] + +ker gs f = ks where + gbar = homomorphismConstruction gs f + gs' = sgs gbar + ks' = dropWhile (\h -> isLeft $ minsupp h) gs' -- !! should filter isRight - sgs might not be in order + ks = map unRight ks' + unRight = fromPairs . map (\(Right a, Right b) -> (a,b)) . toPairs +-} + + +isTransitive gs = length (orbits gs) == 1 + + +-- TRANSITIVE CONSTITUENTS + +{- +-- find largest composition factor of a group which is not transitive +-- we do this by taking the smallest orbit delta, +-- then constructing the homomorphism G -> Sym(delta) +-- and returning the kernel and the image +factorNotTransitive gs = transitiveConstituentHomomorphism' gs delta where + delta = smallest $ orbits gs + sizeSorted lists = map snd $ L.sort $ [(length l, l) | l <- lists] + smallest = head . sizeSorted +-} + +-- Seress p81 +-- A transitive constituent homomorphism is the restriction of G <= Sym(omega) to an orbit delta <= omega +-- This function returns the kernel and the image +transitiveConstituentHomomorphism gs delta + | delta == closure delta [(.^ g) | g <- gs] -- delta is closed under action of gs, hence a union of orbits + = transitiveConstituentHomomorphism' gs delta + +transitiveConstituentHomomorphism' gs delta = (ker, im) where + gs' = sgs $ map (fromPairs . map (\(a,b) -> (lr a, lr b)) . toPairs) gs + -- as delta is a transitive constituent, we will always have a and b either both Left or both Right + lr x = if x `elem` delta then Left x else Right x + ker = map unRight $ dropWhile (isLeft . minsupp) gs' -- pointwise stabiliser of delta + im = map restrictLeft $ takeWhile (isLeft . minsupp) gs' -- restriction of the action to delta + + +-- BLOCKS OF IMPRIMITIVITY + +-- Holt p83ff (and also Seress p107ff) +-- Find a minimal block containing ys. ys are assumed to be sorted. +minimalBlock gs ys@(y1:yt) = minimalBlock' p yt gs where + xs = foldl union [] $ map supp gs + p = M.fromList $ [(yi,y1) | yi <- ys] ++ [(x,x) | x <- xs \\ ys] + minimalBlock' p (q:qs) (h:hs) = + let r = p M.! q -- representative of class containing q + k = p M.! (q .^ h) -- rep of class (q^h) + l = p M.! (r .^ h) -- rep of class (r^h) + in if k /= l -- then we need to merge the classes + then let p' = M.map (\x -> if x == l then k else x) p + qs' = qs ++ [l] + in minimalBlock' p' (q:qs') hs + else minimalBlock' p (q:qs) hs + minimalBlock' p (q:qs) [] = minimalBlock' p qs gs + minimalBlock' p [] _ = + let reps = toListSet $ M.elems p + in L.sort [ filter (\x -> p M.! x == r) xs | r <- reps ] +-- Because the support of the permutations is not constrained to be [1..n], we have to use a map instead of an array +-- This probably affects the complexity, but isn't a problem in practice + +-- Find all block systems +blockSystems gs + | isTransitive gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- xs ] + | otherwise = error "blockSystems: not transitive" + where x:xs = foldl union [] $ map supp gs + + +-- More efficient version if we have an sgs +blockSystemsSGS gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- rs ] + where x:xs = foldl union [] $ map supp gs + hs = filter (\g -> x < minsupp g) gs -- sgs for stabiliser Gx + os = orbits hs + rs = map head os ++ (xs \\ L.sort (concat os)) -- orbit representatives, including singleton cycles +-- Perhaps we could have a function which just returns orbit reps for stabiliser + +-- eg for D 10, the stabiliser of 1 is [[2,6],[3,5]] - we need to make sure we don't forget 4 + +-- If we didn't have an SGS, we could try to randomly generate a few elts of stabiliser Gx, as that would still be better than nothing +-- see Holt RandomStab function + + +isPrimitive gs = null (blockSystems gs) + +isPrimitiveSGS gs = null (blockSystemsSGS gs) + +-- There are other optimisations we haven't done +-- see Holt p86 + + +blockHomomorphism gs bs + | bs == closure bs [(-^ g) | g <- gs] -- bs is closed under action of gs + = blockHomomorphism' gs bs + +blockHomomorphism' gs bs = (ker,im) where + gs' = sgs $ map lr gs + lr g = fromPairs $ [(Left b, Left $ b -^ g) | b <- bs] ++ [(Right x, Right y) | (x,y) <- toPairs g] + ker = map unRight $ dropWhile (isLeft . minsupp) gs' -- stabiliser of the blocks + im = map restrictLeft $ takeWhile (isLeft . minsupp) gs' -- restriction to the action on blocks + +-- Note that there is a slightly more efficient way to calculate block homomorphism, +-- but requires change of base algorithm which we haven't implemented yet + + +-- NORMAL CLOSURE + +-- Seress 115 +-- Given G, H < Sym(Omega) return (the normal closure) +normalClosure gs hs = map unRight $ dropWhile (isLeft . minsupp) $ sgs ks where + xs = foldl union [] $ map supp $ gs ++ hs + ds = map diag gs -- {(g,g) | g <- G} + diag g = fromPairs $ concat [ [(Left x, Left y) , (Right x, Right y)] | (x,y) <- toPairs g] + hsR = map inR hs -- {(1,h) | h <- H} + inR h = fromPairs [(Right x, Right y) | (x,y) <- toPairs h] + ks = ds ++ hsR + +-- Seress 116 +-- Given G, H < Sym(Omega) return `intersection` G +intersectionNormalClosure gs hs = map unRight $ dropWhile (isLeft . minsupp) $ sgs ks where + xs = foldl union [] $ map supp $ gs ++ hs + ds = map diag gs -- {(g,g) | g <- G} + diag g = fromPairs $ concat [ [(Left x, Left y) , (Right x, Right y)] | (x,y) <- toPairs g] + hsL = map inL hs -- {(h,1) | h <- H} + inL h = fromPairs [(Left x, Left y) | (x,y) <- toPairs h] + ks = ds ++ hsL + + +-- CENTRALISER IN THE SYMMETRIC GROUP + +-- Centralizer of G in Sym(X) - transitive case +centralizerSymTrans gs = filter (/= 1) $ centralizerSymTrans' [] fix_g_a where + xs@(a:_) = foldl union [] $ map supp gs + ss = sgs gs + g_a = dropWhile ( (==a) . minsupp ) ss -- pt stabiliser of a + fix_g_a = xs \\ (foldl union [] $ map supp g_a) -- the pts fixed by stabiliser of a + reps_a = cosetRepsGx gs a + -- xs = M.keys reps_a + centralizingElt b = fromPairs [ let g = reps_a M.! x in (x, b .^ g) | x <- xs ] + centralizerSymTrans' ls (r:rs) = + let c = centralizingElt r + in c : centralizerSymTrans' (c:ls) (rs \\ orbitP (c:ls) a) + centralizerSymTrans' _ [] = [] + diff --git a/Math/Combinatorics/LatinSquares.hs b/Math/Combinatorics/LatinSquares.hs --- a/Math/Combinatorics/LatinSquares.hs +++ b/Math/Combinatorics/LatinSquares.hs @@ -1,4 +1,4 @@ - +-- Copyright (c) David Amos, 2009. All rights reserved. module Math.Combinatorics.LatinSquares where diff --git a/Math/Projects/MiniquaternionGeometry.hs b/Math/Projects/MiniquaternionGeometry.hs --- a/Math/Projects/MiniquaternionGeometry.hs +++ b/Math/Projects/MiniquaternionGeometry.hs @@ -1,4 +1,4 @@ - +-- Copyright (c) David Amos, 2009. All rights reserved. module Math.Projects.MiniquaternionGeometry where diff --git a/Math/Projects/Rubik.hs b/Math/Projects/Rubik.hs --- a/Math/Projects/Rubik.hs +++ b/Math/Projects/Rubik.hs @@ -2,11 +2,14 @@ module Math.Projects.Rubik where -import Math.Algebra.Group.PermutationGroup +import Math.Algebra.Group.PermutationGroup hiding (_D) import Math.Algebra.Group.SchreierSims as SS import Math.Algebra.Group.RandomSchreierSims as RSS +import Math.Algebra.Group.Subquotients +-- Rubik's cube + -- 11 12 13 -- 14 U 16 -- 17 18 19 @@ -25,5 +28,88 @@ d = p [[31,33,39,37],[32,36,38,34],[ 7,47,57,27],[ 8,48,58,28],[ 9,49,59,29]] rubikCube = [f,b,l,r,u,d] +-- In Singmaster notation these would be capital letters. --- In Singmaster notation these would be capital letters. +[cornerFaces,edgeFaces] = orbits rubikCube + +(kerCornerFaces,imCornerFaces) = transitiveConstituentHomomorphism rubikCube cornerFaces +-- kernel is the elts which fix all corner faces +-- image is the action restricted to the corner faces + +(kerEdgeFaces,imEdgeFaces) = transitiveConstituentHomomorphism rubikCube edgeFaces +-- kernel is the elts which fix all edge faces +-- image is the action restricted to the edge faces + +[cornerBlocks] = blockSystems imCornerFaces +[edgeBlocks] = blockSystems imEdgeFaces + +(kerCornerBlocks,imCornerBlocks) = blockHomomorphism imCornerFaces cornerBlocks +-- kernel is elts which fix all the corners as blocks, with order 3^7 +-- (Whenever you twist one corner you must untwist another +-- - so the action on 7 corners determines the 8th) +-- image is the action on the corners as blocks, which is S8 of order 20160 + +(kerEdgeBlocks,imEdgeBlocks) = blockHomomorphism imEdgeFaces edgeBlocks +-- kernel is elts which fix all the edges as blocks, with order 2^11 +-- (Whenever you flip one edge, you must flip another edge +-- - so the action on 11 edges determines the 12th) +-- image is the action on the edges as blocks, which is S12 of order 479001600 + +-- Note that orderSGS imCornerFaces * orderSGS imEdgeFaces == 2 * orderSGS (sgs rubikCube) +-- This is because you can't operate on corners and edges totally independently +-- If you swap two corners, you must also swap two edges + +-- See also +-- http://www.gap-system.org/Doc/Examples/rubik.html + +-- (Note that the kernel of the corner constituent homomorphism /= image of edge constituent homomorphism +-- For example, [[36,38],[48,58]] is in the latter, but not the former because it's not in the Rubik group +-- ie there is an elt in the Rubik group which does just that to the edges, but may do some things to the corners) + + +-- Rubik's revenge (4*4*4 cube) + +-- 1 2 3 4 +-- 5 6 7 8 +-- 9 10 11 12 +-- 13 14 15 16 +-- 101 102 103 104 201 202 203 204 301 302 303 304 401 402 403 404 +-- 105 106 107 108 205 206 207 208 305 306 307 308 405 406 407 408 +-- 109 110 111 112 209 210 211 212 309 310 311 312 409 410 411 412 +-- 113 114 115 116 213 214 215 216 313 314 315 316 413 414 415 416 +-- 501 502 503 504 +-- 505 506 507 508 +-- 509 510 511 512 +-- 513 514 515 516 + +_U = p [[1,13,16,4],[2,9,15,8],[3,5,14,12],[6,10,11,7], + [101,201,301,401],[102,202,302,402],[103,203,303,403],[104,204,304,404]] +_u = p [[105,205,305,405],[106,206,306,406],[107,207,307,407],[108,208,308,408]] +_d = p [[109,209,309,409],[110,210,310,410],[111,211,311,411],[112,212,312,412]] +_D = p [[113,213,313,413],[114,214,314,414],[115,215,315,415],[116,216,316,416], + [501,504,516,513],[502,508,515,509],[503,512,514,505],[506,507,511,510]] + +bf = p [[1,304,516,113],[2,308,515,109],[3,312,514,105],[4,316,513,101], + [5,303,512,114],[6,307,511,110],[7,311,510,106],[8,315,509,102], + [9,302,508,115],[10,306,507,111],[11,310,506,107],[12,314,505,103], + [13,301,504,116],[14,305,503,112],[15,309,502,108],[16,313,501,104], + [201,204,216,213],[202,208,215,209],[203,212,214,205],[206,207,211,210], + [401,413,416,404],[402,409,415,408],[403,405,414,412],[406,410,411,407]] + +_R = _U ~^ bf +_r = _u ~^ bf +_l = _d ~^ bf +_L = _D ~^ bf + +ud = _U * _u * _d * _D + +_B = _R ~^ ud +_b = _r ~^ ud +_f = _l ~^ ud +_F = _L ~^ ud + +-- Note that orderSGS $ sgs [_U,_u,_d,_D,bf] comes out much too large, +-- because it includes rotations of the whole cube (24) +-- and exchanges of indistinguishable centre faces (24 for each of 6 colours) +-- So we have to divide by 24^7 / 2. +-- (The /2 is because we can only have even permutations when exchanging indistinguishable centres) diff --git a/Math/Test/TSubquotients.hs b/Math/Test/TSubquotients.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TSubquotients.hs @@ -0,0 +1,56 @@ +-- Copyright (c) David Amos, 2009. All rights reserved. + +module Math.Test.TSubquotients where + +import Math.Algebra.Group.PermutationGroup hiding (ptStab, normalClosure) +-- import Math.Algebra.Group.SchreierSims (cosetRepsGx) +-- import Math.Algebra.Group.RandomSchreierSims + +import Math.Algebra.Group.Subquotients + +import qualified Math.Algebra.Group.PermutationGroup as P -- for testing + + +test = and [testPtStab, testTransitiveConstituentHomomorphism, testBlockSystems, testBlockHomomorphism, + testNormalClosure, testCentralizerSymTrans] + +-- TESTS + +testPtStab = + let gs = [p [[1..5],[6,7]], p [[1,2],[6..10]] ] -- S 5 * S 5 + gs16 = ptStab gs [1,6] + in orderSGS gs16 == 24*24 && orbitP gs16 1 == [1] && orbitP gs16 6 == [6] + +testTransitiveConstituentHomomorphism = + let gs1 = [p [[1..5],[6,7]], p [[1,2],[6..10]] ] -- S 5 * S 5 + (ker1,im1) = transitiveConstituentHomomorphism gs1 [1..5] + gs2 = [p [[1,2,3],[4,5,6]], p [[1,2],[4,5]]] -- note that the two halves don't move independently + (ker2,im2) = transitiveConstituentHomomorphism gs2 [1,2,3] + in orderSGS ker1 == 120 && orderSGS im1 == 120 && + null ker2 && orderSGS im2 == 6 + + +testBlockSystems = + blockSystems [p [[1..3]], p [[4..6]], p [[1,4],[2,5],[3,6]]] == [ [[1,2,3],[4,5,6]] ] && + blockSystems (_D 10) == [] && + blockSystems (_D 12) == [ [[1,3,5],[2,4,6]], [[1,4],[2,5],[3,6]] ] + +testBlockHomomorphism = + let (ker14,im14) = blockHomomorphism (_D 12) [[1,4],[2,5],[3,6]] + (ker135,im135) = blockHomomorphism (_D 12) [[1,3,5],[2,4,6]] + in (orderSGS ker14, orderSGS im14) == (2,6) && + (orderSGS ker135, orderSGS im135) == (6,2) + +-- !! Need to improve this test +testNormalClosure = + let gs = [p [[1,4]]] + hs = [p [[1,2,3]]] + in elts (P.normalClosure gs hs) == elts (normalClosure gs hs) +-- == _A 4 + + +testCentralizerSymTrans = + let gs = [ p [[1..5]] ] + hs = [ p [[1,2],[3,4]], p [[1,3],[2,4]] ] + in centralizerSymTrans gs == gs + && centralizerSymTrans hs == hs diff --git a/Math/Test/TestAll.hs b/Math/Test/TestAll.hs --- a/Math/Test/TestAll.hs +++ b/Math/Test/TestAll.hs @@ -3,6 +3,7 @@ import Math.Test.TGraph import Math.Test.TDesign import Math.Test.TPermutationGroup +import Math.Test.TSubquotients import Math.Test.TFiniteGeometry import Math.Test.TCommutativeAlgebra import Math.Test.TNonCommutativeAlgebra @@ -14,6 +15,7 @@ [Math.Test.TGraph.test ,Math.Test.TDesign.test ,Math.Test.TPermutationGroup.test + ,Math.Test.TSubquotients.test ,Math.Test.TFiniteGeometry.test ,Math.Test.TCommutativeAlgebra.test ,Math.Test.TField.test