HaskellForMaths 0.1.9 → 0.2.0

raw patch · 10 files changed

+668/−22 lines, 10 filesdep +QuickCheckdep +arraydep +randomPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, array, random

API changes (from Hackage documentation)

Files

• HaskellForMaths.cabal +6/−4
• Math/Algebra/Group/RandomSchreierSims.hs +135/−0
• Math/Algebra/Group/SchreierSims.hs +6/−2
• Math/Combinatorics/Design.hs +1/−1
• Math/Combinatorics/LatinSquares.hs +122/−0
• Math/Combinatorics/StronglyRegularGraph.hs +16/−12
• Math/Projects/ChevalleyGroup/Classical.hs +1/−0
• Math/Projects/MiniquaternionGeometry.hs +372/−0
• Math/Projects/Rubik.hs +4/−1
• Math/Test/TPermutationGroup.hs +5/−2

@@ -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

@@ -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+-}

@@ -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
 

@@ -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).

@@ -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

@@ -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
 

@@ -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

@@ -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]+-}

@@ -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.

@@ -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