HaskellForMaths 0.1.6 → 0.1.7

raw patch · 13 files changed

+283/−105 lines, 13 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• HaskellForMaths.cabal +1/−1
• Math/Algebra/Commutative/MPoly.hs +4/−4
• Math/Algebra/Commutative/Monomial.hs +3/−3
• Math/Algebra/Field/Base.hs +22/−6
• Math/Algebra/Field/Extension.hs +52/−25
• Math/Algebra/Group/PermutationGroup.hs +43/−13
• Math/Algebra/Group/SchreierSims.hs +1/−0
• Math/Combinatorics/Design.hs +10/−11
• Math/Combinatorics/FiniteGeometry.hs +75/−13
• Math/Combinatorics/Graph.hs +24/−27
• Math/Combinatorics/GraphAuts.hs +46/−0
• Math/Combinatorics/StronglyRegularGraph.hs +1/−1
• Math/Common/ListSet.hs +1/−1

@@ -1,5 +1,5 @@    Name:                HaskellForMaths
-   Version:             0.1.6
+   Version:             0.1.7
    Category:            Math
    Description:         Math library - combinatorics, group theory, commutative algebra, non-commutative algebra
    License:             BSD3

@@ -15,9 +15,9 @@ 
 -- MULTIVARIATE POLYNOMIALS
 
--- |Type for multivariate polynomials
--- |ord is a phantom type defining how terms are ordered, r is the type of the ring we are working over
--- |For example, a common choice will be MPoly Grevlex Q, meaning polynomials over Q with the grevlex term ordering
+-- |Type for multivariate polynomials.
+-- ord is a phantom type defining how terms are ordered, r is the type of the ring we are working over.
+-- For example, a common choice will be MPoly Grevlex Q, meaning polynomials over Q with the grevlex term ordering
 newtype MPoly ord r = MP [(Monomial ord,r)] deriving (Eq)
 -- deriving instance (Ord (Monomial ord), Ord r) => Ord (MPoly ord r)
 -- standalone deriving supported from GHC 6.8
@@ -76,7 +76,7 @@     recip _ = error "MPoly.recip: only supported for (non-zero) constants or monomials"
 
 -- |Create a variable with the supplied name.
--- |By convention, variable names should usually be a single letter followed by none, one or two digits
+-- By convention, variable names should usually be a single letter followed by none, one or two digits.
 var :: String -> MPoly Grevlex Q
 var v = MP [(Monomial $ M.singleton v 1, 1)] :: MPoly Grevlex Q
 

@@ -39,9 +39,9 @@ -- |Phantom type representing grevlex term ordering
 data Grevlex
 
--- |Phantom type for an elimination term ordering
--- |In the ordering, xis come before yjs come before zks, but within the xis, or yjs, or zks, grevlex ordering is used 
-data Elim -- a term order for elimination
+-- |Phantom type for an elimination term ordering.
+-- In the ordering, xis come before yjs come before zks, but within the xis, or yjs, or zks, grevlex ordering is used 
+data Elim
 
 
 diffs a b = M.elems m where Monomial m = a/b

@@ -10,7 +10,7 @@ 
 -- RATIONALS
 
--- Rationals with a better show function
+-- |Q is just the rationals, but with a better show function than the Prelude version
 newtype Q = Q Rational deriving (Eq,Ord,Num,Fractional)
 
 instance Show Q where
@@ -53,7 +53,7 @@                    in Fp $ u `mod` p
                    where p = value (undefined :: n)
 
-class FiniteField fq where
+class Fractional fq => FiniteField fq where
     eltsFq :: fq -> [fq]  -- return all elts of the field
     basisFq :: fq -> [fq] -- return an additive basis for the field (as Z-module)
 
@@ -71,17 +71,33 @@ char fq = head [p | p <- [2..], length fq `mod` p == 0]
 
 
+-- |F2 is a type for the finite field with 2 elements
 type F2 = Fp T2
-f2 = map fromInteger [0..1] :: [F2]
 
+-- |f2 lists the elements of F2
+f2 :: [F2]
+f2 = map fromInteger [0..1] -- :: [F2]
+
+-- |F3 is a type for the finite field with 3 elements
 type F3 = Fp T3
-f3 = map fromInteger [0..2] :: [F3]
 
+-- |f3 lists the elements of F3
+f3 :: [F3]
+f3 = map fromInteger [0..2] -- :: [F3]
+
+-- |F5 is a type for the finite field with 5 elements
 type F5 = Fp T5
-f5 = map fromInteger [0..4] :: [F5]
 
+-- |f5 lists the elements of F5
+f5 :: [F5]
+f5 = map fromInteger [0..4] -- :: [F5]
+
+-- |F7 is a type for the finite field with 7 elements
 type F7 = Fp T7
-f7 = map fromInteger [0..6] :: [F7]
+
+-- |f7 lists the elements of F7
+f7 :: [F7]
+f7 = map fromInteger [0..6] -- :: [F7]
 
 type F11 = Fp T11
 f11 = map fromInteger [0..10] :: [F11]

@@ -4,6 +4,8 @@ 
 module Math.Algebra.Field.Extension where
 
+import Data.List as L (elemIndex)
+
 import Math.Common.IntegerAsType
 import Math.Algebra.Field.Base
 
@@ -16,23 +18,26 @@ x = UP [0,1] :: UPoly Integer
 
 instance (Show a, Num a) => Show (UPoly a) where
-    show (UP []) = "0"
-    show (UP as) = let powers = filter ( (/=0) . fst ) $ zip as [0..]
-                       c:cs = concatMap showTerm powers
-                   in if c == '+' then cs else c:cs
-        where showTerm (a,i) = showCoeff a ++ showPower a i
-              showCoeff a = case show a of
-                            "1" -> "+"
-                            "-1" -> "-"
-                            '-':cs -> '-':cs
-                            cs -> '+':cs
-              showPower a i | i == 0 = case show a of
-                                       "1" -> "1"
-                                       "-1" -> "1"
-                                       otherwise -> ""
-                            | i == 1 = "x"
-                            | i > 1  = "x^" ++ show i
+    -- show (UP []) = "0"
+    show (UP as) = showUP "x" as
 
+showUP _ [] = "0"
+showUP v as = let powers = filter ( (/=0) . fst ) $ zip as [0..]
+                  c:cs = concatMap showTerm powers
+              in if c == '+' then cs else c:cs
+    where showTerm (a,i) = showCoeff a ++ showPower a i
+          showCoeff a = case show a of
+                        "1" -> "+"
+                        "-1" -> "-"
+                        '-':cs -> '-':cs
+                        cs -> '+':cs
+          showPower a i | i == 0 = case show a of
+                                   "1" -> "1"
+                                   "-1" -> "1"
+                                   otherwise -> ""
+                        | i == 1 = v -- "x"
+                        | i > 1  = v ++ "^" ++ show i -- "x^" ++ show i
+
 instance Num a => Num (UPoly a) where
 	UP as + UP bs = toUPoly $ as <+> bs
 	negate (UP as) = UP $ map negate as
@@ -98,7 +103,9 @@ data ExtensionField k poly = Ext (UPoly k) deriving (Eq,Ord)
 
 instance Num k => Show (ExtensionField k poly) where
-    show (Ext f) = show f
+    -- show (Ext f) = show f
+    -- show (Ext (UP [])) = "0"
+    show (Ext (UP as)) = showUP "a" as
 
 instance (Num k, Fractional k, PolynomialAsType k poly) => Num (ExtensionField k poly) where
     Ext x + Ext y = Ext $ (x+y) `modUP` pvalue (undefined :: (k,poly))
@@ -112,7 +119,7 @@                         (u,v,1) = extendedEuclidUP f g -- so u*f + v*g == 1. (We know the gcd is 1 as g is irreducible)
                     in Ext $ u `modUP` g
 
-instance (Num k, FiniteField k, PolynomialAsType k poly) => FiniteField (ExtensionField k poly) where
+instance (FiniteField k, PolynomialAsType k poly) => FiniteField (ExtensionField k poly) where
     eltsFq _ = map Ext (polys (d-1) fp) where
         fp = eltsFq (undefined :: k)
         d = deg $ pvalue (undefined :: (k,poly))
@@ -134,45 +141,58 @@ instance PolynomialAsType F2 ConwayF4 where pvalue _ = convert $ x^2+x+1
 type F4 = ExtensionField F2 ConwayF4
 f4 = map Ext (polys 2 f2) :: [F4]
-x4 = embed x :: F4
+a4 = embed x :: F4
 
 data ConwayF8
 instance PolynomialAsType F2 ConwayF8 where pvalue _ = convert $ x^3+x+1
 type F8 = ExtensionField F2 ConwayF8
 f8 = map Ext (polys 3 f2) :: [F8]
-x8 = embed x :: F8
+a8 = embed x :: F8
 
 data ConwayF9
 instance PolynomialAsType F3 ConwayF9 where pvalue _ = convert $ x^2+2*x+2
 type F9 = ExtensionField F3 ConwayF9
 f9 = map Ext (polys 2 f3) :: [F9]
-x9 = embed x :: F9
+a9 = embed x :: F9
 
 data ConwayF16
 instance PolynomialAsType F2 ConwayF16 where pvalue _ = convert $ x^4+x+1
 type F16 = ExtensionField F2 ConwayF16
 f16 = map Ext (polys 4 f2) :: [F16]
-x16 = embed x :: F16
+a16 = embed x :: F16
 
 data ConwayF25
 instance PolynomialAsType F5 ConwayF25 where pvalue _ = convert $ x^2+4*x+2
 type F25 = ExtensionField F5 ConwayF25
 f25 = map Ext (polys 2 f5) :: [F25]
-x25 = embed x :: F25
+a25 = embed x :: F25
 
 data ConwayF27
 instance PolynomialAsType F3 ConwayF27 where pvalue _ = convert $ x^3+2*x+1
 type F27 = ExtensionField F3 ConwayF27
 f27 = map Ext (polys 3 f3) :: [F27]
-x27 = embed x :: F27
+a27 = embed x :: F27
 
 data ConwayF32
 instance PolynomialAsType F2 ConwayF32 where pvalue _ = convert $ x^5+x^2+1
 type F32 = ExtensionField F2 ConwayF32
 f32 = map Ext (polys 5 f2) :: [F32]
-x32 = embed x :: F32
+a32 = embed x :: F32
 
 
+-- generator for the automorphism group of fq, as applied to an element of fq
+frobeniusAut x = x ^ p where
+    p = char $ eltsFq x
+
+-- the degree of fq as an extension over fp
+-- (hence also, the order of the automorphism group of fq)
+degree fq = n where
+    q = length fq
+    p = char fq
+    Just n = L.elemIndex q $ iterate (*p) 1
+
+
+
 -- QUADRATIC EXTENSIONS OF Q
 
 data Sqrt a = Sqrt a
@@ -205,3 +225,10 @@ 
 type QSqrtMinus5 = ExtensionField Q (Sqrt (M TMinus1 T5))
 sqrtminus5 = embed x :: QSqrtMinus5
+
+
+-- conjugation automorphism of quadratic field
+-- conjugate of a + b sqrt d is a - b sqrt d
+conjugate :: ExtensionField Q (Sqrt d) -> ExtensionField Q (Sqrt d)
+conjugate (Ext (UP [a,b])) = Ext (UP [a,-b])
+conjugate x = x -- the zero or constant cases

@@ -34,18 +34,24 @@ supp (P g) = M.keys g
 -- (This is guaranteed not to contain fixed points provided the permutations have been constructed using the supplied constructors)
 
--- |x .^ g returns the image of a vertex or point x under the action of the permutation g
+-- |x .^ g returns the image of a vertex or point x under the action of the permutation g.
+-- The dot is meant to be a mnemonic for point or vertex.
 (.^) :: (Ord k) => k -> Permutation k -> k
 x .^ P g = case M.lookup x g of
            Just y  -> y
            Nothing -> x -- if x `notElem` supp (P g), then x is not moved
 
+-- |b -^ g returns the image of an edge or block b under the action of the permutation g
+-- The dash is meant to be a mnemonic for edge or line or block.
+(-^) :: (Ord t) => [t] -> Permutation t -> [t]
+xs -^ g = L.sort [x .^ g | x <- xs]
+
 -- construct a permutation from cycles
 fromCycles cs = fromPairs $ concatMap fromCycle cs
     where fromCycle xs = zip xs (rotateL xs)
 
--- |Construct a permutation from a list of cycles
--- |For example, p [[1,2,3],[4,5]] returns the permutation that sends 1 to 2, 2 to 3, 3 to 1, 4 to 5, 5 to 4
+-- |Construct a permutation from a list of cycles.
+-- For example, p [[1,2,3],[4,5]] returns the permutation that sends 1 to 2, 2 to 3, 3 to 1, 4 to 5, 5 to 4
 p :: (Ord a) => [[a]] -> Permutation a
 p cs = fromCycles cs
 -- can't specify in pointfree style because of monomorphism restriction
@@ -83,7 +89,8 @@ instance (Ord a, Show a) => Fractional (Permutation a) where
     recip = inverse
 
--- |g ~^ h returns the conjugate of g by h
+-- |g ~^ h returns the conjugate of g by h.
+-- The tilde is meant to a mnemonic, because conjugacy is an equivalence relation.
 (~^) :: (Ord t, Show t) => Permutation t -> Permutation t -> Permutation t
 g ~^ h = h^-1 * g * h
 
@@ -93,11 +100,6 @@ 
 -- ORBITS
 
--- |b -^ g returns the image of an edge or block b under the action of g
-(-^) :: (Ord t) => [t] -> Permutation t -> [t]
-xs -^ g = L.sort [x .^ g | x <- xs]
-
-
 closureS xs fs = closure' S.empty (S.fromList xs) where
     closure' interior boundary
         | S.null boundary = interior
@@ -110,14 +112,20 @@ 
 orbit action x gs = closure [x] [ (`action` g) | g <- gs]
 
+-- |x .^^ gs returns the orbit of the point or vertex x under the action of the gs
+(.^^) :: (Ord a) => a -> [Permutation a] -> [a]
 x .^^ gs = orbit (.^) x gs
+
 orbitP gs x = orbit (.^) x gs
 orbitV gs x = orbit (.^) x gs
 
--- orbit of a block
+-- |b -^^ gs returns the orbit of the block or edge b under the action of the gs
+(-^^) :: (Ord t) => [t] -> [Permutation t] -> [[t]]
 b -^^ gs = orbit (-^) b gs
+
 orbitB gs b = orbit (-^) b gs
 orbitE gs b = orbit (-^) b gs
+
 {-
 -- orbit of a vertex / point
 x .^^ gs = closure [x] [ .^g | g <- gs]
@@ -131,7 +139,6 @@ -}
 
 action xs f = fromPairs [(x, f x) | x <- xs]
--- probably supercedes the three following functions
 
 
 -- find all the orbits of a group
@@ -205,21 +212,41 @@     mapping = M.fromList $ zip _X [1..] -- the mapping from _X to [1..n]
     toSn' g = fromPairs' $ map (\(x,x') -> (mapping M.! x, mapping M.! x')) $ toPairs g
 
+-- Given a permutation over lists of small positive integers, such as [1,2,3],
+-- return a permutation over the integers obtained by interpreting the lists as digits.
+-- For example, [1,2,3] -> 123.
+fromDigits g = fromPairs [(fromDigits' x, fromDigits' y) | (x,y) <- toPairs g]
 
+fromDigits' xs = f (reverse xs) where
+    f (x:xs) = x + 10 * f xs
+    f [] = 0
+
+-- Given a permutation over lists of 0s and 1s,
+-- return the permutation obtained by interpreting these as binary digits.
+-- For example, [1,1,0] -> 6.
+fromBinary g = fromPairs [(fromBinary' x, fromBinary' y) | (x,y) <- toPairs g]
+
+fromBinary' xs = f (reverse xs) where
+    f (x:xs) = x + 2 * f xs
+    f [] = 0
+
+
+
+
 -- INVESTIGATING GROUPS
 -- Functions to investigate groups in various ways
 -- Most of these functions will only be efficient for small groups (say |G| < 10000)
 -- For larger groups we will need to use Schreier-Sims and associated algorithms
 
 -- |Given generators for a group, return a (sorted) list of all elements of the group.
--- |Implemented using a naive closure algorithm, so only suitable for small groups (|G| < 10000)
+-- Implemented using a naive closure algorithm, so only suitable for small groups (|G| < 10000)
 elts :: (Num a, Ord a) => [a] -> [a]
 elts gs = closure [1] [ (*g) | g <- gs]
 
 eltsS gs = closureS [1] [ (*g) | g <- gs]
 
 -- |Given generators for a group, return the order of the group (the number of elements).
--- |Implemented using a naive closure algorithm, so only suitable for small groups (|G| < 10000)
+-- Implemented using a naive closure algorithm, so only suitable for small groups (|G| < 10000)
 order :: (Num a, Ord a) => [a] -> Int
 order gs = S.size $ eltsS gs -- length $ elts gs
 
@@ -280,6 +307,9 @@ conjClass gs h = closure [h] [ (~^ g) | g <- gs]
 -- conjClass gs h = h ~^^ gs
 
+-- |conjClassReps gs returns a conjugacy class representatives and sizes for the group generated by gs.
+-- This implementation is only suitable for use with small groups (|G| < 10000).
+conjClassReps :: (Ord t, Show t) => [Permutation t] -> [(Permutation t, Int)]
 conjClassReps gs = conjClassReps' (elts gs) where
     conjClassReps' (h:hs) =
         let cc = conjClass gs h in (h, length cc) : conjClassReps' (hs \\ cc)

@@ -142,6 +142,7 @@ 
 -- |Given generators for a group, return a (sorted) list of all elements of the group, using Schreier-Sims algorithm
 elts :: (Ord t, Show t) => [Permutation t] -> [Permutation t]
+elts [] = [1]
 elts gs = eltsBSGS $ bsgs gs
 
 -- |Given generators for a group, return the order of the group (the number of elements), using Schreier-Sims algorithm

@@ -12,8 +12,8 @@ import Math.Algebra.Field.Extension
 import Math.Algebra.Group.PermutationGroup hiding (elts, order, isMember)
 import Math.Algebra.Group.SchreierSims as SS
-import Math.Combinatorics.Graph as G hiding (to1n, combinationsOf)
-import Math.Combinatorics.GraphAuts (refine, isSingleton, graphAuts, removeGens)
+import Math.Combinatorics.Graph as G hiding (to1n)
+import Math.Combinatorics.GraphAuts (refine, isSingleton, graphAuts, incidenceAuts, removeGens)
 import Math.Combinatorics.FiniteGeometry
 
 -- Cameron & van Lint, Designs, Graphs, Codes and their Links
@@ -21,11 +21,6 @@ 
 {-
 set xs = map head $ group $ sort xs
-
--- subsets of size k (returned in ascending order)
-combinationsOf 0 _ = [[]]
-combinationsOf _ [] = []
-combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs
 -}
 isSubset xs ys = all (`elem` ys) xs
 
@@ -200,14 +195,17 @@ isDesignAut (D xs bs) g | supp g `isSubset` xs = all (`S.member` bs') [b -^ g | b <- bs]
     where bs' = S.fromList bs
 
-incidenceGraph (D xs bs) = graph (vs,es) where
+incidenceGraph (D xs bs) = G vs es where -- graph (vs,es) where
     vs = L.sort $ map Left xs ++ map Right bs
     es = L.sort [ [Left x, Right b] | x <- xs, b <- bs, x `elem` b ]
 
+
+designAuts d = incidenceAuts $ incidenceGraph d
+
 -- We find design auts by finding graph auts of the incidence graph of the design
 -- In a square design, we need to watch out for graph auts which are mapping points <-> blocks
 designAuts1 d = filter (/=1) $ map points $ graphAuts $ incidenceGraph d where
-    points (P m) = P $ M.fromList [(x,y) | (Left x, Left y) <- M.toList m]
+    points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h]
      -- This implicitly filters out (Right x, Right y) action on blocks,
      -- and also (Left x, Right y) auts taking points to blocks.
      -- The filter (/=1) is to remove points <-> blocks auts
@@ -232,11 +230,12 @@ -- This won't have much effect on the point cell, but will hopefully split the block cells in half
 
 
+-- !! superceded by Math.Combinatorics.GraphAuts.incidenceAuts, leading to designAuts above
 -- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
-designAuts d@(D xs bs) = map points (designAuts' [] [vs]) where
+designAuts2 d@(D xs bs) = map points (designAuts' [] [vs]) where
     n = length xs
     g@(G vs es) = incidenceGraph d
-    points (P m) = P $ M.fromList [(k,v) | (Left k, Left v) <- M.toList m]
+    points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
     designAuts' us p@((x@(Left _):ys):pt) =
         let p' = L.sort $ filter (not . null) $ refine (ys:pt) (dps M.! x)
         in level us p x ys []

@@ -1,29 +1,23 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) David Amos, 2008-2009. All rights reserved.
 
 module Math.Combinatorics.FiniteGeometry where
 
 import Data.List as L
 import qualified Data.Set as S
--- import qualified Data.Map as M -- not really required
+
 import Math.Algebra.Field.Base
 import Math.Algebra.Field.Extension hiding ( (<+>) )
 import Math.Algebra.LinearAlgebra -- hiding ( det )
 
--- import PermutationGroup
--- import SchreierSims as SS
-
-
--- !! This should really live somewhere else
--- subsets of size k
-combinationsOf 0 _ = [[]]
-combinationsOf _ [] = []
-combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs
-
-
+import Math.Combinatorics.Graph
 
+-- |ptsAG n fq returns the points of the affine geometry AG(n,Fq), where fq are the elements of Fq
+ptsAG :: (FiniteField a) => Int -> [a] -> [[a]]
 ptsAG 0 fq = [[]]
 ptsAG n fq = [x:xs | x <- fq, xs <- ptsAG (n-1) fq]
 
+-- |ptsPG n fq returns the points of the projective geometry PG(n,Fq), where fq are the elements of Fq
+ptsPG :: (FiniteField a) => Int -> [a] -> [[a]]
 ptsPG 0 _ = [[1]]
 ptsPG n fq = map (0:) (ptsPG (n-1) fq) ++ map (1:) (ptsAG n fq)
 
@@ -37,6 +31,8 @@ ispnf (1:xs) = True
 ispnf _ = False
 
+-- closure of points in AG(n,Fq)
+-- result is sorted
 closureAG ps =
     let multipliers = [ (1 - sum xs) : xs | xs <- ptsAG (k-1) fq ] -- k-vectors over fq whose sum is 1
     in S.toList $ S.fromList [foldl1 (<+>) $ zipWith (*>) m ps | m <- multipliers]
@@ -99,6 +95,9 @@     oneColumn r = replicate (r-1) Zero ++ One : replicate (k-r) Zero
     starColumn r = replicate (r-1) Star ++ replicate (k+1-r) Zero
 
+
+-- |flatsPG n fq k returns the k-flats in PG(n,Fq), where fq are the elements of Fq
+flatsPG :: (FiniteField a) => Int -> [a] -> Int -> [[[a]]]
 flatsPG n fq k = concatMap substStars $ rrefs (n+1) (k+1) where
     substStars (r:rs) = [r':rs' | r' <- substStars' r, rs' <- substStars rs]
     substStars [] = [[]]
@@ -109,9 +108,72 @@ 
 
 -- Flats in AG(n,Fq) are just the flats in PG(n,Fq) which are not "at infinity"
+-- |flatsAG n fq k returns the k-flats in AG(n,Fq), where fq are the elements of Fq
+flatsAG :: (FiniteField a) => Int -> [a] -> Int -> [[[a]]]
 flatsAG n fq k = [map tail (r : map (r <+>) rs) | r:rs <- flatsPG n fq k, head r == 1]
 -- The head r == 1 condition is saying that we want points which are in the "finite" part of PG(n,Fq), not points at infinity
 -- The reason we add r to each of the rs is to bring them into the "finite" part
 -- (If you don't do this, it can lead to incorrect results, for example some of the flats having the same closure)
 
 
+-- |The lines (1-flats) in PG(n,fq)
+linesPG :: (FiniteField a) => Int -> [a] -> [[[a]]]
+linesPG n fq = flatsPG n fq 1
+
+-- |The lines (1-flats) in AG(n,fq)
+linesAG :: (FiniteField a) => Int -> [a] -> [[[a]]]
+linesAG n fq = flatsAG n fq 1
+
+
+-- less efficient but perhaps more intuitive
+-- a line in AG(n,fq) is a translation (x) of a line through the origin (y)
+linesAG1 n fq = [ [x,z] | x <- ptsAG n fq, y <- ptsPG (n-1) fq,
+                          z <- [x <+> y], [x,z] == take 2 (closureAG [x,z]) ]
+-- the point of the condition at the end is to avoid listing the same line more than once
+
+
+-- INCIDENCE GRAPH
+
+-- |Incidence graph of PG(n,fq), considered as an incidence structure between points and lines
+incidenceGraphPG :: (Ord a, FiniteField a) => Int -> [a] -> Graph (Either [a] [[a]])
+incidenceGraphPG n fq = G vs es where
+    points = ptsPG n fq
+    lines = linesPG n fq
+    vs = L.sort $ map Left points ++ map Right lines
+    es = L.sort [ [Left x, Right b] | b <- lines, x <- closurePG b]
+-- Could also consider incidence structure between points and planes, etc
+
+-- incidenceAuts (incidenceGraphPG n fq) == PGL(n,fq) * auts fq
+-- For example, incidenceAuts (incidenceGraphPG 2 f4) =
+-- PGL(3,f4) * auts f4
+-- where PGL(3,f4)/PSL(3,f4) == f4* (multiplicative group of f4),
+-- and auts f4 == { 1, x -> x^2 } (the Frobenius aut)
+-- The full group is called PGammaL(3,f4)
+
+
+
+-- |Incidence graph of AG(n,fq), considered as an incidence structure between points and lines
+incidenceGraphAG :: (Ord a, FiniteField a) => Int -> [a] -> Graph (Either [a] [[a]])
+incidenceGraphAG n fq = G vs es where
+    points = ptsAG n fq
+    lines = linesAG n fq
+    vs = L.sort $ map Left points ++ map Right lines
+    es = L.sort [ [Left x, Right b] | b <- lines, x <- closureAG b]
+
+-- incidenceAuts (incidenceGraphAG n fq) == Aff(n,fq) * auts fq
+-- where Aff(n,fq), the affine group, is the semi-direct product GL(n,fq) * (fq^n)+
+-- where (fq^n)+ is the additive group of translations
+-- Each elt of Aff(n,fq) is of the form x -> ax + b, where a <- GL(n,fq), b <- (fq^n)+
+
+orderGL n q = product [q^n - q^i | i <- [0..n-1] ]
+-- for the first row, we can choose any vector except zero, hence q^n-1
+-- for the second row, we can choose any vector except a multiple of the first, hence q^n-q
+-- etc
+
+orderAff n q = q^n * orderGL n q
+
+
+-- NOTE:
+-- AG(n,F2) is degenerate:
+-- Every pair of points is a line, so it is the complete graph on 2^n points
+-- And as such as aut group S(2^n)

@@ -9,7 +9,7 @@ import Control.Arrow ( (&&&) )
 
 import Math.Common.ListSet
-import Math.Algebra.Group.PermutationGroup
+import Math.Algebra.Group.PermutationGroup hiding (fromDigits, fromBinary)
 import Math.Algebra.Group.SchreierSims as SS
 
 -- Main source: Godsil & Royle, Algebraic Graph Theory
@@ -24,8 +24,8 @@ powerset [] = [[]]
 powerset (x:xs) = let p = powerset xs in p ++ map (x:) p
 
--- |combinationsOf k xs returns the subsets of xs of size k
--- |If xs is in ascending order, then the returned list is in ascending order
+-- |combinationsOf k xs returns the subsets of xs of size k.
+-- If xs is in ascending order, then the returned list is in ascending order
 combinationsOf :: (Integral t) => t -> [a] -> [[a]]
 combinationsOf 0 _ = [[]]
 combinationsOf _ [] = []
@@ -34,9 +34,9 @@ 
 -- GRAPH
 
--- |Datatype for graphs, represented as a list of vertices and a list of edges
--- |Both the list of vertices and the list of edges, and also the 2-element lists representing the edges,
--- |are required to be in ascending order, without duplicates
+-- |Datatype for graphs, represented as a list of vertices and a list of edges.
+-- Both the list of vertices and the list of edges, and also the 2-element lists representing the edges,
+-- are required to be in ascending order, without duplicates.
 data Graph a = G [a] [[a]] deriving (Eq,Ord,Show)
 
 -- we require that vs, es, and each individual e are sorted
@@ -45,7 +45,8 @@ isGraph vs es = isSetSystem vs es && all ( (==2) . length) es
 
 -- |Safe constructor for graph from lists of vertices and edges.
--- |graph (vs,es) checks that vs and es are valid before returning the graph.
+-- graph (vs,es) checks that vs and es are valid before returning the graph.
+graph :: (Ord t) => ([t], [[t]]) -> Graph t
 graph (vs,es) | isGraph vs es = G vs es
 -- isValid g = g where g = G vs es
 
@@ -88,25 +89,25 @@ nullGraph :: Graph Int -- type signature needed
 nullGraph = G [] []
 
--- |The cyclic graph on n vertices
+-- |c n is the cyclic graph on n vertices
 c :: (Integral t) => t -> Graph t
 c n = graph (vs,es) where
     vs = [1..n]
     es = L.insert [1,n] [[i,i+1] | i <- [1..n-1]]
 -- automorphism group is D2n
 
--- |The complete graph on n vertices
+-- |k n is the complete graph on n vertices
 k :: (Integral t) => t -> Graph t
 k n = graph (vs,es) where
     vs = [1..n]
     es = [[i,j] | i <- [1..n-1], j <- [i+1..n]] -- == combinationsOf 2 [1..n]
 -- automorphism group is Sn
 
--- |The complete bipartite graph on m and n vertices
+-- The complete bipartite graph on m and n vertices
 -- kb :: (Integral t) => t -> t -> Graph t
 kb m n = to1n $ kb' m n
 
--- |The complete bipartite graph on m left and n right vertices
+-- The complete bipartite graph on m left and n right vertices
 -- kb :: (Integral t) => t -> t -> Graph (Either t t)
 kb' m n = graph (vs,es) where
     vs = map Left [1..m] ++ map Right [1..n]
@@ -119,22 +120,13 @@     hammingDistance as bs = length $ filter id $ zipWith (/=) as bs
 -- can probably type-coerce this to be Graph [F2] if required
 
+q k = fromBinary $ q' k
+{-
 -- note, this definition only in versions >0.1.3
-q k = gmap (\v -> v <.> pows2) (q' k) where
+q'' k = gmap (\v -> v <.> pows2) (q' k) where
     pows2 = reverse $ take k $ iterate (*2) 1
     u <.> v = sum $ zipWith (*) u v
     gmap f (G vs es) = G (map f vs) ((map . map) f es)
-
-{-
--- definitions in versions <= 0.1.3
-q k = let vs = zip [0..] (powerset [1..k])
-          es = [ [i,j] | (i,iset) <- vs, (j,jset) <- vs, i < j, length (iset `symDiff` jset) == 1 ]
-      in graph (map fst vs,es)
-
-q' k = let us = powerset $ map (2^) [0..k-1]
-           vs = [0..2^k-1] -- == L.sort $ map sum us
-           es = L.sort [ L.sort [sum u, sum v] | [u,v] <- combinationsOf 2 us, length (u `symDiff` v) == 1 ]
-       in graph (vs, es)
 -}
 
 tetrahedron = k 4
@@ -169,7 +161,7 @@     vs' = M.elems mapping
     es' = [map (mapping M.!) e | e <- es] -- the edges will already be sorted correctly by construction
 
--- Given a graph with vertices which are lists of small integers, eg [1,2,3]
+-- |Given a graph with vertices which are lists of small integers, eg [1,2,3],
 -- return a graph with vertices which are the numbers obtained by interpreting these as digits, eg 123.
 -- The caller is responsible for ensuring that this makes sense (eg that the small integers are all < 10)
 fromDigits :: Integral a => Graph [a] -> Graph a
@@ -177,9 +169,14 @@     vs' = map fromDigits' vs
     es' = (map . map) fromDigits' es
 
-fromDigits' xs = f (reverse xs) where
-    f (x:xs) = x + 10 * f xs
-    f [] = 0
+-- |Given a graph with vertices which are lists of 0s and 1s,
+-- return a graph with vertices which are the numbers obtained by interpreting these as binary digits.
+-- For example, [1,1,0] -> 6.
+fromBinary :: Integral a => Graph [a] -> Graph a
+fromBinary (G vs es) = graph (vs',es') where
+    vs' = map fromBinary' vs
+    es' = (map . map) fromBinary' es
+
 
 -- this definition only in versions >0.1.3
 petersen = graph (vs,es) where

@@ -146,6 +146,7 @@             in if all isSingleton p1'
                then let xys' = xys ++ zip (concat p1') (concat p2')
                     in if isCompatible xys' then [fromPairs' xys'] else []
+                    -- we shortcut the search when we have all singletons, so must check isCompatible to avoid false positives
                else let (x:xs):p1'' = p1'
                         ys:p2'' = p2'
                     in concat [dfs ((x,y):xys)
@@ -162,6 +163,8 @@ 
 -- Now we try to use generators we've already found at a given level to save us having to look for others
 -- For example, if we have found (1 2)(3 4) and (1 3 2), then we don't need to look for something taking 1 -> 4
+-- |Given a graph g, graphAuts g returns a strong generating set for the automorphism group of g.
+graphAuts :: (Ord a) => Graph a -> [Permutation a]
 graphAuts g@(G vs es) = graphAuts' [] [vs] where
     graphAuts' us p@((x:ys):pt) =
         let p' = L.sort $ filter (not . null) $ refine (ys:pt) (dps M.! x)
@@ -196,6 +199,49 @@     es' = S.fromList es
 
 -- contrary to first thought, you can't stop when a level is null - eg kb 2 3, the third level is null, but the fourth isn't
+
+
+-- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH
+
+-- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
+incidenceAuts g@(G vs es) = map points (incidenceAuts' [] [vs]) where
+    points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
+    -- points (P m) = P $ M.fromList [(k,v) | (Left k, Left v) <- M.toList m]
+    incidenceAuts' us p@((x@(Left _):ys):pt) =
+        let p' = L.sort $ filter (not . null) $ refine (ys:pt) (dps M.! x)
+        in level us p x ys []
+        ++ incidenceAuts' (x:us) p'
+    incidenceAuts' us ([]:pt) = incidenceAuts' us pt
+    incidenceAuts' _ (((Right _):_):_) = [] -- if we fix all the points, then the blocks must be fixed too 
+    -- incidenceAuts' _ [] = []
+    level us p@(ph:pt) x (y@(Left _):ys) hs =
+        let px = refine (L.delete x ph : pt) (dps M.! x)
+            py = refine (L.delete y ph : pt) (dps M.! y)
+            uus = zip us us
+        in case dfs ((x,y):uus) px py of
+           []  -> level us p x ys hs
+           h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
+    level _ _ _ _ _ = [] -- includes the case where y matches Right _, which can only occur on first level, before we've distance partitioned
+    dfs xys p1 p2
+        | map length p1 /= map length p2 = []
+        | otherwise =
+            let p1' = filter (not . null) p1
+                p2' = filter (not . null) p2
+            in if all isSingleton p1'
+               then let xys' = xys ++ zip (concat p1') (concat p2')
+                    in if isCompatible xys' then [fromPairs' xys'] else []
+               else let (x:xs):p1'' = p1'
+                        ys:p2'' = p2'
+                    in concat [dfs ((x,y):xys)
+                                   (refine (xs : p1'') (dps M.! x))
+                                   (refine ((L.delete y ys):p2'') (dps M.! y))
+                                   | y <- ys]
+    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+    dps = M.fromList [(v, distancePartition g v) | v <- vs]
+    es' = S.fromList es
+
+
+
 
 
 removeGens x gs = removeGens' [] gs where

@@ -15,7 +15,7 @@ import Math.Combinatorics.Design as D
 import Math.Algebra.LinearAlgebra -- hiding (t)
 import Math.Algebra.Field.Base -- for F2
-import Math.Combinatorics.FiniteGeometry hiding (combinationsOf)
+import Math.Combinatorics.FiniteGeometry
 
 -- Sources
 -- Godsil & Royle, Algebraic Graph Theory

@@ -1,4 +1,4 @@-
+-- Copyright (c) David Amos, 2008. All rights reserved.
 
 module Math.Common.ListSet where