diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal --- a/HaskellForMaths.cabal +++ b/HaskellForMaths.cabal @@ -1,5 +1,5 @@ Name: HaskellForMaths - Version: 0.1.9 + Version: 0.2.0 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 @@ -22,19 +22,21 @@ Math/Test/TestAll.hs Library - Build-Depends: base >=3 && < 4, containers + Build-Depends: base >=3 && < 4, containers, array, random, QuickCheck Exposed-modules: Math.Algebra.LinearAlgebra, Math.Algebra.Commutative.Monomial, Math.Algebra.Commutative.MPoly, Math.Algebra.Commutative.GBasis, Math.Algebra.Field.Base, Math.Algebra.Field.Extension, - Math.Algebra.Group.PermutationGroup, Math.Algebra.Group.SchreierSims, Math.Algebra.Group.StringRewriting, + Math.Algebra.Group.PermutationGroup, Math.Algebra.Group.SchreierSims, + Math.Algebra.Group.RandomSchreierSims, 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.Combinatorics.LatinSquares, Math.Common.IntegerAsType, Math.Common.ListSet, Math.Projects.RootSystem, - Math.Projects.Rubik, + Math.Projects.Rubik, Math.Projects.MiniquaternionGeometry, Math.Projects.ChevalleyGroup.Classical, Math.Projects.ChevalleyGroup.Exceptional, Math.Projects.KnotTheory.Braid, Math.Projects.KnotTheory.LaurentMPoly, Math.Projects.KnotTheory.TemperleyLieb, Math.Projects.KnotTheory.IwahoriHecke diff --git a/Math/Algebra/Group/RandomSchreierSims.hs b/Math/Algebra/Group/RandomSchreierSims.hs new file mode 100644 --- /dev/null +++ b/Math/Algebra/Group/RandomSchreierSims.hs @@ -0,0 +1,135 @@ + + +module Math.Algebra.Group.RandomSchreierSims where + +import System.Random +import Data.List as L +import qualified Data.Map as M +import Data.Maybe + +import Control.Monad +import Data.Array.MArray +import Data.Array.IO +import System.IO.Unsafe + +import Math.Common.ListSet (toListSet) +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 + +-- Holt p69-71 +-- Product replacement algorithm for generating uniformly distributed random elts of a black box group + +initProdRepl :: (Ord a, Show a) => [Permutation a] -> IO (Int, IOArray Int (Permutation a)) +initProdRepl gs = + let n = length gs + r = max 10 n + xs = (1:) $ take r $ concat $ repeat gs + in do xs' <- newListArray (0,r) xs + replicateM_ 60 $ nextProdRepl (r,xs') -- perform initial mixing + return (r,xs') + +nextProdRepl :: (Ord a, Show a) => (Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a)) +nextProdRepl (r,xs) = + do s <- randomRIO (1,r) + t <- randomRIO (1,r) + u <- randomRIO (0,3 :: Int) + out <- updateArray xs s t u + return out + +updateArray xs s t u = + let (swap,invert) = quotRem u 2 in + if s == t + then return Nothing + else do + x_0 <- readArray xs 0 + x_s <- readArray xs s + x_t <- readArray xs t + let x_s' = mult (swap,invert) x_s x_t + x_0' = mult (swap,0) x_0 x_s' + writeArray xs 0 x_0' + writeArray xs s x_s' + return (Just x_0') + where mult (swap,invert) a b = case (swap,invert) of + (0,0) -> a * b + (0,1) -> a * b^-1 + (1,0) -> b * a + (1,1) -> b^-1 * a + + +-- Holt p97-8 +-- Random Schreier-Sims algorithm, for finding strong generating set of permutation group + +-- It's possible that the following code can be improved by introducing levels only as we need them? + +-- |Given generators for a permutation group, return a strong generating set. +-- The result is calculated using random Schreier-Sims algorithm, so has a small (\<10^-6) chance of being incomplete. +-- The sgs is relative to the base implied by the Ord instance. +sgs :: (Ord a, Show a) => [Permutation a] -> [Permutation a] +sgs gs = toListSet $ concatMap snd $ rss gs + +rss gs = unsafePerformIO $ + do (r,xs) <- initProdRepl gs + rss' (r,xs) (initLevels gs) 0 + +rss' (r,xs) levels i + | i == 25 = return levels -- stop if we've had 25 successful sifts in a row + | otherwise = do g <- nextProdRepl (r,xs) + let (changed,levels') = updateLevels levels g + rss' (r,xs) levels' (if changed then 0 else i+1) +-- if we currently have an sgs for a subgroup of the group, then it must have index >= 2 +-- so the chance of a random elt sifting to identity is <= 1/2 + +initLevels gs = [((b,M.singleton b 1),[]) | b <- bs] + where bs = toListSet $ concatMap supp gs + +updateLevels levels Nothing = (False,levels) -- not strictly correct to increment count on a Nothing +updateLevels levels (Just g) = + case sift (map fst levels) g of + Nothing -> (False, levels) + -- Just 1 -> error "Just 1" + Just g' -> (True, updateLevels' [] levels g' (minsupp g')) + +updateLevels' ls (r@((b,t),s):rs) h b' = + if b == b' + then reverse ls ++ ((b, cosetRepsGx (h:s) b), h:s) : rs + else updateLevels' (r:ls) rs h b' +-- updateLevels' ls [] h b' = error $ "updateLevels: " ++ show (ls,[],h,b') + +-- used the following in debugging +-- orderLevels levels = product $ [if M.null t then 1 else toInteger (M.size t) | ((b,t),s) <- levels] + + +-- 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 + +-- |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 +isMemberSGS gs h = let bts = baseTransversalsSGS gs in isNothing $ sift bts h + + +{- +-- Alternative where we carry on with Schreier-Sims when we finish Random Schreier-Sims, just to make sure +-- !! Unfortunately, doesn't appear to work - perhaps ss' doesn't like finding empty levels +sgs2 gs = toListSet $ concatMap snd $ rss2 gs + +rss2 gs = unsafePerformIO $ + do (r,xs) <- initProdRepl gs + levels <- rss' (r,xs) (initLevels gs) 0 + return $ ss' bs (reverse levels) [] + where bs = toListSet $ concatMap supp gs +-} 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 @@ -41,10 +41,12 @@ -- SCHREIER-SIMS ALGORITHM sift bts g = sift' bts g where - sift' [] g = if g == 1 then Nothing else Just g + 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 findBase gs = minimum $ concatMap supp gs {- @@ -72,7 +74,9 @@ -} --- strong generating set, with implied base from the Ord instance +-- |Given generators for a permutation group, return a strong generating set. +-- The result is calculated using Schreier-Sims algorithm, and is relative to the base implied by the Ord instance +sgs :: (Ord a, Show a) => [Permutation a] -> [Permutation a] sgs gs = toListSet $ concatMap snd $ ss bs gs where bs = toListSet $ concatMap supp gs diff --git a/Math/Combinatorics/Design.hs b/Math/Combinatorics/Design.hs --- a/Math/Combinatorics/Design.hs +++ b/Math/Combinatorics/Design.hs @@ -118,7 +118,7 @@ -- the projective plane PG(2,Fq) - a square 2-(q^2+q+1,q+1,1) design pg2 fq = design (points, lines) where points = ptsPG 2 fq - lines = map line points + lines = L.sort $ map line points line u = [v | v <- points, u <.> v == 0] u <.> v = sum (zipWith (*) u v) -- Remember that the points and lines of PG(2,Fp) are really the lines and planes of AG(3,Fp). diff --git a/Math/Combinatorics/LatinSquares.hs b/Math/Combinatorics/LatinSquares.hs new file mode 100644 --- /dev/null +++ b/Math/Combinatorics/LatinSquares.hs @@ -0,0 +1,122 @@ + + +module Math.Combinatorics.LatinSquares where + +import qualified Data.List as L +import qualified Data.Set as S +import qualified Data.Map as M + +-- import Math.Combinatorics.FiniteGeometry +import Math.Combinatorics.Design +import Math.Algebra.Field.Base +import Math.Algebra.Field.Extension +import Math.Algebra.LinearAlgebra (fMatrix') +import Math.Combinatorics.Graph +import Math.Combinatorics.GraphAuts +import Math.Combinatorics.StronglyRegularGraph + + +-- LATIN SQUARES + +findLatinSqs xs = findLatinSqs' 1 [xs] where + n = length xs + findLatinSqs' i rows + | i == n = [reverse rows] + | otherwise = concat [findLatinSqs' (i+1) (row:rows) + | row <- findRows (L.transpose rows) [] xs] + findRows (col:cols) ls rs = concat [findRows cols (r:ls) (L.delete r rs) + | r <- rs, r `notElem` col] + findRows [] ls _ = [reverse ls] + +isLatinSq rows = all isOneOfEach rows && all isOneOfEach cols where + cols = L.transpose rows + +isOneOfEach xs = length xs == S.size (S.fromList xs) + + +-- The incidence graph of a latin square +-- It is distance-regular +-- Godsil & Royle p69 +incidenceGraphLS l = graph (vs,es) where + n = length l -- the order of the latin square + vs = [ (i, j, l ! (i,j)) | i <- [1..n], j <- [1..n] ] + es = [ [v1,v2] | [v1@(i,j,lij), v2@(i',j',lij')] <- combinationsOf 2 vs, i == i' || j == j' || lij == lij' ] + m ! (i,j) = m !! (i-1) !! (j-1) + +incidenceGraphLS' l = graph (vs,es) where + n = length l -- the order of the latin square + vs = [ (i, j) | i <- [1..n], j <- [1..n] ] + es = [ [v1,v2] | [v1@(i,j), v2@(i',j')] <- combinationsOf 2 vs, i == i' || j == j' || l' M.! (i,j) == l' M.! (i',j') ] + l' = M.fromList [ ( (i,j), l !! (i-1) !! (j-1) ) | i <- [1..n], j <- [1..n] ] +-- vertices are grid positions +-- adjacent if they're in the same row, same column, or have the same entry + + +-- ORTHOGONAL AND MUTUALLY ORTHOGONAL LATINS SQUARES + +isOrthogonal greeks latins = isOneOfEach pairs + where pairs = zip (concat greeks) (concat latins) + +findMOLS k lsqs = findMOLS' k [] lsqs where + findMOLS' 0 ls _ = [reverse ls] + findMOLS' i ls (r:rs) = + if all (isOrthogonal r) ls + then findMOLS' (i-1) (r:ls) rs ++ findMOLS' i ls rs + else findMOLS' i ls rs + findMOLS' _ _ [] = [] + +isMOLS (greek:latins) = all (isOrthogonal greek) latins && isMOLS latins +isMOLS [] = True + +-- MOLS from a projective plane +fromProjectivePlane (D xs bs) = map toLS parallelClasses where + k = [x | [0,1,x] <- xs] -- the field we're working over + n = length k -- the order of the projective plane + parallelClasses = drop 2 $ L.groupBy (\l1 l2 -> head l1 == head l2) bs + -- The classes of parallel lines + -- Each line has its ideal point at its head + -- The first two classes have [0,0,1] and [0,1,0] as ideal points, and hence consist of horizontal and vertical lines + toLS ls = let grid = M.fromList [ ((x,y),i) | (i, [0,1,mu]:ps) <- zip [1..] ls, [1,x,y] <- ps] + in fMatrix' n (\i j -> grid M.! (k !! i, k !! j)) + + +-- ORTHOGONAL ARRAYS +-- Godsil & Royle p224 + +isOA (k,n) rows = + length rows == k && + all ( (== n^2) . length ) rows && + all isOneOfEach [zip ri rj | [ri,rj] <- combinationsOf 2 rows ] + +-- An OA(3,n) from a latin square +fromLS l = + [ concat [replicate n i | i <- [1..n] ] -- row numbers + , concat (replicate n [1..n]) -- column numbers + , concat l -- entries + ] + where n = length l -- the order of the latin square + +fromMOLS mols = + (concat [replicate n i | i <- [1..n] ]) : -- row numbers + (concat (replicate n [1..n]) ) : -- column numbers + map concat mols -- entries for each lsq + where n = length $ head mols -- the order of the latin squares + +-- The graph defined by an OA(k,n) +-- It is strongly regular with parameters ( n^2, (n-1)k, n-2+(k-1)(k-2), k(k-1) ) +-- Godsil & Royle p225 +graphOA rows = graph (vs,es) where + vs = L.transpose rows -- the vertices are the columns of the OA + es = [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, or (zipWith (==) v1 v2) ] + -- two vertices are adjacent if they agree in any position + +-- Expected SRG parameters +srgParamsOA (k,n) = Just ( n^2, (n-1)*k, n-2+(k-1)*(k-2), k*(k-1) ) + +-- eg srgParams (4,4) == srgParams $ graphOA $ init $ fromMOLS $ fromProjectivePlane $ pg2 f4 + + +-- Todo: +-- Would like a way to find out to what extent two sets of MOLS are really the same, +-- eg can one be obtained from the other by row or column reordering (with renumbering) +-- This might provide a proof of the distinctness of phi, omega, omegaD, psi diff --git a/Math/Combinatorics/StronglyRegularGraph.hs b/Math/Combinatorics/StronglyRegularGraph.hs --- a/Math/Combinatorics/StronglyRegularGraph.hs +++ b/Math/Combinatorics/StronglyRegularGraph.hs @@ -26,18 +26,22 @@ -- STRONGLY REGULAR GRAPHS -- strongly regular graphs -srgParams g = - let vs = vertices g - n = length vs - es = edges g - es' = combinationsOf 2 vs \\ es -- the non-edges - k:ks = map (valency g) vs - lambda:ls = map (length . commonNbrs) es -- common neighbours of adjacent vertices - mu:ms = map (length . commonNbrs) es' -- common neighbours of non-adjacent vertices - commonNbrs [v1,v2] = nbrs g v1 `intersect` nbrs g v2 - in if all (==k) ks && all (==lambda) ls && all (==mu) ms - then Just (n,k,lambda,mu) - else Nothing +srgParams g + | null es = error "srgParams: not defined for null graph" + | null es' = error "srgParams: not defined for complete graph" + | otherwise = + if all (==k) ks && all (==lambda) ls && all (==mu) ms + then Just (n,k,lambda,mu) + else Nothing + where vs = vertices g + n = length vs + es = edges g + es' = combinationsOf 2 vs \\ es -- the non-edges + k:ks = map (valency g) vs + lambda:ls = map (length . commonNbrs) es -- common neighbours of adjacent vertices + mu:ms = map (length . commonNbrs) es' -- common neighbours of non-adjacent vertices + commonNbrs [v1,v2] = (nbrs_g M.! v1) `intersect` (nbrs_g M.! v2) + nbrs_g = M.fromList [ (v, nbrs g v) | v <- vs ] isSRG g = isJust $ srgParams g diff --git a/Math/Projects/ChevalleyGroup/Classical.hs b/Math/Projects/ChevalleyGroup/Classical.hs --- a/Math/Projects/ChevalleyGroup/Classical.hs +++ b/Math/Projects/ChevalleyGroup/Classical.hs @@ -21,6 +21,7 @@ -- LINEAR GROUPS -- SL(n,Fq) is generated by elementary transvections +-- sl :: FiniteField k => Int -> k -> [[[k]]] sl n fq = [elemTransvection n (r,c) l | r <- [1..n], c <- [1..n], r /= c, l <- fq'] where fq' = basisFq undefined -- tail fq -- Carter p68 - x_r(t1) x_r(t2) == x_r(t1+t2) - this is true in general, not just in this case diff --git a/Math/Projects/MiniquaternionGeometry.hs b/Math/Projects/MiniquaternionGeometry.hs new file mode 100644 --- /dev/null +++ b/Math/Projects/MiniquaternionGeometry.hs @@ -0,0 +1,372 @@ + + +module Math.Projects.MiniquaternionGeometry where + +import qualified Data.List as L + +import Math.Common.ListSet as LS + +import Math.Algebra.Field.Base +import Math.Combinatorics.FiniteGeometry (pnf, ispnf, orderPGL) +import Math.Combinatorics.Graph (combinationsOf) +import Math.Combinatorics.GraphAuts +import Math.Algebra.Group.PermutationGroup hiding (order) +import qualified Math.Algebra.Group.SchreierSims as SS +import Math.Algebra.Group.RandomSchreierSims +import Math.Combinatorics.Design as D +import Math.Algebra.LinearAlgebra -- ( (<.>), (<+>) ) + +import Math.Projects.ChevalleyGroup.Classical + +import Test.QuickCheck + + + +-- Sources: +-- Miniquaternion Geometry, Room & Kirkpatrick +-- Survey of Non-Desarguesian Planes, Charles Weibel + + +-- F9, defined by adding sqrt of -1 to F3. (The Conway poly for F9 is not so convenient for us here) +data F9 = F9 F3 F3 deriving (Eq,Ord) + +instance Show F9 where + show (F9 0 0) = "0" + show (F9 0 1) = "e" + show (F9 0 2) = "-e" + show (F9 1 0) = "1" + show (F9 1 1) = "1+e" + show (F9 1 2) = "1-e" + show (F9 2 0) = "-1" + show (F9 2 1) = "-1+e" + show (F9 2 2) = "-1-e" + +e = F9 0 1 -- sqrt of -1 + +instance Num F9 where + F9 a1 b1 + F9 a2 b2 = F9 (a1+a2) (b1+b2) + F9 a1 b1 * F9 a2 b2 = F9 (a1*a2-b1*b2) (a1*b2+a2*b1) + negate (F9 a b) = F9 (negate a) (negate b) + fromInteger n = F9 (fromInteger n) 0 + +f9 = [F9 a b | a <- f3, b <- f3] + +w = 1-e -- a primitive element - generates the multiplicative group + +conj (F9 a b) = F9 a (-b) +-- This is just the Frobenius aut x -> x^3 + +norm (F9 a b) = a^2 + b^2 +-- == x * conj x + +instance Fractional F9 where + recip x@(F9 a b) = F9 (a/n) (-b/n) where n = norm x + +instance FiniteField F9 where + basisFq _ = [1,e] + + +-- J9, or Q, defined by modifying the multiplication in F9 +data J9 = J9 F9 deriving (Eq,Ord) + +instance Show J9 where + show (J9 (F9 0 0)) = "0" + show (J9 (F9 0 1)) = "-j" + show (J9 (F9 0 2)) = "j" + show (J9 (F9 1 0)) = "1" + show (J9 (F9 1 1)) = "-k" + show (J9 (F9 1 2)) = "i" + show (J9 (F9 2 0)) = "-1" + show (J9 (F9 2 1)) = "-i" + show (J9 (F9 2 2)) = "k" + +squaresF9 = [1,w^2,w^4,w^6] -- and 0, but not needed here + +instance Num J9 where + J9 x + J9 y = J9 (x+y) + 0 * _ = 0 + _ * 0 = 0 + J9 x * J9 y = + if y `elem` squaresF9 + then J9 (x*y) + else J9 (conj x * y) + negate (J9 x) = J9 (negate x) + fromInteger n = J9 (fromInteger n) + +i = J9 w +j = J9 (w^6) -- == i-1 +k = J9 (w^7) -- == i+1 + +j9 = [J9 x | x <- f9] + + +-- the aut of J9 that sends i to x +autJ9 x = fromPairs [ (a+b*i, a+b*x) | a <- [0,1,-1], b <- [1,-1] ] + +autA = autJ9 (-i) -- sends i -> -i +autB = autJ9 (-k) -- sends j -> -j +autC = autJ9 (-j) -- sends k -> -k + +autsJ9 = [autA, autC] +-- these two auts generate the group, which is isomorphic to S3 +-- indeed, the auts permute the pairs {i,-i}, {j,-j}, {k,-k} + + +conj' (J9 x) = J9 (conj x) +-- Note that conj' x == x .^ autB + + +isAut k sigma = and [sigma x + sigma y == sigma (x+y) | x <- k, y <- k] + && and [sigma x * sigma y == sigma (x*y) | x <- k, y <- k] + + +isReal x = x `elem` [0,1,-1] +isComplex = not . isReal + +instance Fractional J9 where + recip 0 = error "J9.recip: 0" + recip x | isReal x = x + | otherwise = -x + +instance FiniteField J9 where + basisFq _ = [1,i,j,k] + eltsFq _ = j9 + + +-- Near fields + +prop_NearField (a,b,c) = + a+(b+c) == (a+b)+c && -- addition is associative + a+b == b+a && -- addition is commutative + a+0 == a && -- additive identity + a+(-a) == 0 && -- additive inverse + a*(b*c) == (a*b)*c && -- multiplication is associative + a*1 == a && 1*a == a && -- multiplicative identity + (a+b)*c == a*c + b*c && -- right-distributivity + a*0 == 0 + +instance Arbitrary F9 where + arbitrary = do x <- arbitrary :: Gen Int + return (f9 !! (x `mod` 9)) + coarbitrary = undefined -- !! only required if we want to test functions over the type + +instance Arbitrary J9 where + arbitrary = do x <- arbitrary :: Gen Int + return (j9 !! (x `mod` 9)) + coarbitrary = undefined -- !! only required if we want to test functions over the type + +prop_NearFieldF9 (a,b,c) = prop_NearField (a,b,c) where + types = (a,b,c) :: (F9,F9,F9) + +prop_NearFieldJ9 (a,b,c) = prop_NearField (a,b,c) where + types = (a,b,c) :: (J9,J9,J9) + + + +-- PROJECTIVE PLANES + +ptsPG2 r = [ [0,0,1] ] ++ [ [0,1,x] | x <- r ] ++ [ [1,x,y] | x <- r, y <- r ] +-- if r is sorted, then so is the result + +orthogonalLinesPG2 xs = L.sort [ [x | x <- xs, u <.> x == 0] | u <- xs ] + +rightLinesPG2 r = + [ [0,0,1] : [ [0,1,x] | x <- r] ] ++ -- line at infinity + [ [0,0,1] : [ [1,x,y] | y <- r] | x <- r ] ++ -- vertical lines + [ [0,1,a] : [ [1,x,y] | x <- r, y <- [x*a+b] ] | a <- r, b <- r ] -- slope multiplies on the right +-- if r is sorted, then so is the result, and each line in the result + +leftLinesPG2 r = + [ [0,0,1] : [ [0,1,x] | x <- r] ] ++ -- line at infinity + [ [0,0,1] : [ [1,x,y] | y <- r] | x <- r ] ++ -- vertical lines + [ [0,1,a] : [ [1,x,y] | x <- r, y <- [a*x+b] ] | a <- r, b <- r ] -- slope multiplies on the left + + +-- Projective plane PG2(F9) +phi = design (xs,bs) where + xs = ptsPG2 f9 + bs = orthogonalLinesPG2 xs -- L.sort [ [x | x <- xs, u <.> x == 0] | u <- xs ] + +-- Then the collineations of phi consist of projective transformations, +-- together with a conjugacy collineation induced by the Frobenius aut + +-- alternative construction of PG2(F9) - gives same result +phi' = design (xs,bs) where + xs = ptsPG2 f9 + bs = rightLinesPG2 f9 + + +collineationsPhi = l 3 f9 ++ [fieldAut] where + D xs bs = phi + fieldAut = fromPairs [ (x , map conj x) | x <- xs ] +-- in general, this would be PSigmaL(n,Fq), whereas we want PGammaL(n,Fq). However, for F9 they coincide. +-- order 84913920 + + +liftToGraph (D xs bs) g = fromPairs $ [(Left x, Left (x .^ g)) | x <- xs] ++ [(Right b, Right (b -^ g)) | b <- bs] + + + +-- This construction appears to produce a projective plane +-- (However, Room & Kirkpatrick point out that it's not really well-defined +-- - if we had chosen different quasi-homogeneous coords, we would have got different results) +-- However, it's not the same as either omega or omegaD below +omega0 = design (xs,bs) where + xs = ptsPG2 j9 + bs = orthogonalLinesPG2 xs -- L.sort [ [x | x <- xs, u <.> x == 0] | u <- xs ] + + +-- Room & Kirkpatrick, p103 +omega = design (xs,bs) where + xs = ptsPG2 j9 + bs = rightLinesPG2 j9 + +-- another construction that produces same result (but slower) +omega2 = design (xs,bs) where + xs = ptsPG2 j9 + bs = [ l | [p,q] <- combinationsOf 2 xs, l <- [line p q], [p,q] == take 2 l] + line p q = toListSet $ filter ispnf [(a *> p) <+> (b *> q) | a <- j9, b <- j9] + + +-- Room & Kirkpatrick, p107, p114 +collineationsOmega = + [r] + ++ [s rho sigma | rho <- j9 \\ [0], sigma <- j9 \\ [0], rho == 1 || sigma == 1] + ++ [t delta epsilon | delta <- j9, epsilon <- j9, delta * epsilon == 0] -- for generators sufficient to have only one non-zero + ++ [u] + ++ [a lambda | lambda <- autsJ9] where + D xs bs = omega + fromMatrix m = fromPairs [ (x, pnf (x <*>> m)) | x <- xs] + r = fromMatrix [[1,0,0],[0,0,1],[0,1,0]] -- reflect in the line x = y in the affine subplane + s rho sigma = fromPairs $ [([1,x,y], [1,x*rho,y*sigma]) | x <- j9, y <- j9] + ++ [([0,1,mu],[0,1,(recip rho)*mu*sigma]) | mu <- j9] + ++ [([0,0,1],[0,0,1])] -- leaves "Y" fixed + -- fromMatrix [[1,0,0],[0,rho,0],[0,0,sigma]] -- scale x,y -> rho x, sigma y + t delta epsilon = fromMatrix [[1,delta,epsilon],[0,1,0],[0,0,1]] -- translation x,y -> x+delta, y+epsilon + u = fromPairs $ [([1,x,y], [1,x+y,x-y]) | x <- j9, y <- j9] + ++ [([0,1,mu],[0,1,-mu]) | mu <- filter isComplex j9] + ++ [([0,1,0],[0,1,1]), ([0,1,1],[0,1,0]), ([0,1,-1],[0,0,1]), ([0,0,1],[0,1,-1])] + -- fromMatrix [[1,0,0],[0,1,-1],[0,1,1]] + a lambda = fromPairs [ (x, map (.^ lambda) x) | x <- xs] +-- order 311040 +-- (which means this is also the plane constructed in Weibel?) + + +-- dual plane of omega +omegaD = design (xs,bs) where + xs = ptsPG2 j9 + bs = leftLinesPG2 j9 + +omegaD1 = D.to1n $ dual omega +-- need proof omega /~= omegaD + +omegaD2 = design (xs,bs) where + xs = ptsPG2 j9 + bs = [ l | [p,q] <- combinationsOf 2 xs, l <- [line p q], [p,q] == take 2 l] + line p q = toListSet $ filter ispnf [(p <* a) <+> (q <* b) | a <- j9, b <- j9] + +us <* x = map (*x) us + + +-- Room and Kirkpatrick p130 +psi = design (xs,bs) where + xs = ptsPG2 j9 + isReal x = all (`elem` [0,1,-1]) x + xrs = ptsPG2 [0,1,-1] -- the thirteen real points, a copy of PG2(F3) within psi + bs = toListSet [line p q | p <- xrs, q <- xs, q /= p] + line p q = L.sort $ p : [pnf ( (p <* a) <+> q) | a <- j9] + + +-- Room & Kirkpatrick p137 +psi2 = design (xs,bs) where + xs = ptsPG2 j9 + bs = L.sort $ + [ [0,0,1] : [ [0,1,x] | x <- j9] ] ++ -- line at infinity, z=0 + [ [0,0,1] : [ [1,kappa,y] | y <- j9] | kappa <- j9 ] ++ -- vertical lines x = kappa + [ [0,1,m] : [ [1,x,m*x+kappa] | x <- j9 ] | m <- [0,1,-1], kappa <- j9 ] ++ -- lines with real slope + [ [0,1,kappa] : [ [1,x,kappa*(x-r)+s] | x <- j9 ] | r <- [0,1,-1], s <- [0,1,-1], kappa <- j9 \\ [0,1,-1] ] + -- lines with complex slope + +-- Room & Kirkpatrick p134-6 +collineationsPsi = realProjectivities -- real transvections, generating real projectivities + ++ [a lambda | lambda <- autsJ9] where + D xs bs = psi + n = 3 + realTransvections = [elemTransvection n (r,c) l | r <- [1..n], c <- [1..n], r /= c, l <- [1]] + realProjectivities = [fromPairs $ [(x, pnf (x <*>> m)) | x <- xs] | m <- realTransvections] + a lambda = fromPairs [ (x, map (.^ lambda) x) | x <- xs] +-- order 33696 + + +-- The order of a projective plane +order (D xs bs) = length (head bs) - 1 + +isProjectivePlane pi = designParams pi == Just (2,(q^2+q+1,q,1)) + where q = order pi + + +collinear (D xs bs) ys = (not . null) [b | b <- bs, all (`elem` b) ys] + +-- assume p1..4 are distinct +isQuadrangle plane ps@[p1,p2,p3,p4] = + all (not . collinear plane) (combinationsOf 3 ps) + + +concurrent (D xs bs) ls = (not . null) [x | x <- xs, all (x `elem`) ls] + +isQuadrilateral plane ls@[l1,l2,l3,l4] = + all (not . concurrent plane) (combinationsOf 3 ls) + + +isOval pi ps = length ps == order pi+1 + && all (not . collinear pi) (combinationsOf 3 ps) + +findOvals1 pi = findOvals' 0 ([], points pi) where + n = order pi + findOvals' i (ls,rs) + | i == n+1 = [reverse ls] + | otherwise = concatMap (findOvals' (i+1)) + [ (r:ls, rs') | r:rs' <- L.tails rs, all (not . collinear pi) (map (r:) (combinationsOf 2 ls)) ] +-- if we have a function to quickly generate the line through two points, +-- then we just need to see whether the third point is on it, which is much faster than testing collinearity + +findQuadrangles pi = findQuadrangles' 0 ([], points pi) where + findQuadrangles' i (ls,rs) + | i == 4 = [reverse ls] + | otherwise = concatMap (findQuadrangles' (i+1)) + [ (r:ls, rs') | r:rs' <- L.tails rs, all (not . collinear pi) (map (r:) (combinationsOf 2 ls)) ] + + +findOvals pi@(D xs bs) = findOvals' 0 ([],xs) bs where + n = order pi + findOvals' i (ls,rs) bs + | i == n+1 = [reverse ls] + | otherwise = concat + [let rls = reverse (r:ls) + (notchords, chords) = L.partition (\b -> length (rls `LS.intersect` b) < 2) bs + rs'' = foldl (\\) rs' chords + -- if any line is already a chord, remove remaining points on it from further consideration + in findOvals' (i+1) (r:ls, rs'') notchords + | r:rs' <- L.tails rs] + +-- Todo: +-- Code that shows that phi is Desarguesian, and omega, omegaD and psi are not +{- +-- !! NOT WORKING +-- finds apparent counterexamples in phi too +findNonDesarguesian pi@(D xs bs) = + [ [p,x,y,z,x',y',z',k,l,m] | p <- xs, + x <- xs \\ [p], + y <- xs \\ [p,x], + z <- xs \\ [p,x,y], + (not . collinear pi) [x,y,z], + x' <- line p x \\ L.sort [p,x], + y' <- line p y \\ L.sort [p,y], + z' <- line p z \\ L.sort [p,z], + (not . collinear pi) [x',y',z'], + k <- line x y `intersect` line x' y', -- will only have one element + l <- line x z `intersect` line x' z', + m <- line y z `intersect` line y' z', + (not . collinear pi) [k,l,m] ] + where line p q = head [b | b <- bs, p `elem` b, q `elem` b] +-} diff --git a/Math/Projects/Rubik.hs b/Math/Projects/Rubik.hs --- a/Math/Projects/Rubik.hs +++ b/Math/Projects/Rubik.hs @@ -3,7 +3,8 @@ module Math.Projects.Rubik where import Math.Algebra.Group.PermutationGroup -import Math.Algebra.Group.SchreierSims +import Math.Algebra.Group.SchreierSims as SS +import Math.Algebra.Group.RandomSchreierSims as RSS -- 11 12 13 @@ -22,5 +23,7 @@ r = p [[41,43,49,47],[42,46,48,44],[ 3,13,57,33],[ 6,16,54,36],[ 9,19,51,39]] u = p [[11,13,19,17],[12,16,18,14],[ 1,21,51,41],[ 2,22,52,42],[ 3,23,53,43]] 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. diff --git a/Math/Test/TPermutationGroup.hs b/Math/Test/TPermutationGroup.hs --- a/Math/Test/TPermutationGroup.hs +++ b/Math/Test/TPermutationGroup.hs @@ -8,8 +8,10 @@ import Math.Algebra.Group.PermutationGroup as P import Math.Algebra.Group.SchreierSims as SS +import Math.Algebra.Group.RandomSchreierSims as RSS import Math.Combinatorics.Graph import Math.Combinatorics.GraphAuts +import Math.Projects.Rubik import Test.QuickCheck hiding (choose) @@ -74,7 +76,7 @@ choose n m | m <= n = product [m+1..n] `div` product [1..n-m] -test = and [sgsTest, ssTest, ccTest] +test = and [sgsTest, ssTest, ccTest, rubikTest] sgsTest = all (uncurry (==)) sgsTests @@ -86,8 +88,9 @@ [let _G = toSn (_S 3 `dp` _S 3) in (sgsOrder _G, SS.order _G) ] ++ [let _G = toSn (_C 3 `wr` _S 3) in (sgsOrder _G, SS.order _G) ] ++ [let _G = toSn (_S 3 `wr` _C 3) in (sgsOrder _G, SS.order _G) ] - where sgsOrder = orderTGS . tgsFromSgs . sgs + where sgsOrder = orderTGS . tgsFromSgs . SS.sgs +rubikTest = orderSGS (RSS.sgs rubikCube) == 43252003274489856000 ssTest = all (uncurry (==)) ssTests