packages feed

HaskellForMaths 0.2.0 → 0.2.1

raw patch · 9 files changed

+365/−27 lines, 9 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

HaskellForMaths.cabal view
@@ -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,
Math/Algebra/Group/RandomSchreierSims.hs view
@@ -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
Math/Algebra/Group/SchreierSims.hs view
@@ -12,7 +12,7 @@ 
 -- COSET REPRESENTATIVES FOR STABILISER OF A POINT
 
--- Given a group G = <gs>, and a point x, find coset representatives for Gx (stabiliser of x) in G
+-- Given a group G = <gs>, 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] []
+ Math/Algebra/Group/Subquotients.hs view
@@ -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 <H^G> (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 <H^G> `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' _ [] = []+
Math/Combinatorics/LatinSquares.hs view
@@ -1,4 +1,4 @@-+-- Copyright (c) David Amos, 2009. All rights reserved.  module Math.Combinatorics.LatinSquares where
Math/Projects/MiniquaternionGeometry.hs view
@@ -1,4 +1,4 @@-+-- Copyright (c) David Amos, 2009. All rights reserved.  module Math.Projects.MiniquaternionGeometry where
Math/Projects/Rubik.hs view
@@ -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)
+ Math/Test/TSubquotients.hs view
@@ -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
Math/Test/TestAll.hs view
@@ -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