HaskellForMaths 0.4.5 → 0.4.6

raw patch · 46 files changed

+1222/−579 lines, 46 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/Field/Extension.hs +23/−5
• Math/Algebra/Group/PermutationGroup.hs +24/−5
• Math/Algebra/Group/SchreierSims.hs +1/−1
• Math/Algebra/LinearAlgebra.hs +3/−1
• Math/Algebras/AffinePlane.hs +2/−2
• Math/Algebras/Commutative.hs +12/−2
• Math/Algebras/GroupAlgebra.hs +3/−1
• Math/Algebras/LaurentPoly.hs +1/−1
• Math/Algebras/Matrix.hs +7/−5
• Math/Algebras/NonCommutative.hs +4/−2
• Math/Algebras/Octonions.hs +3/−1
• Math/Algebras/Quaternions.hs +5/−3
• Math/Algebras/Structures.hs +4/−2
• Math/Algebras/TensorAlgebra.hs +3/−1
• Math/Algebras/TensorProduct.hs +6/−4
• Math/Algebras/VectorSpace.hs +44/−19
• Math/Combinatorics/CombinatorialHopfAlgebra.hs +28/−26
• Math/Combinatorics/Digraph.hs +8/−8
• Math/Combinatorics/FiniteGeometry.hs +4/−2
• Math/Combinatorics/Graph.hs +37/−13
• Math/Combinatorics/GraphAuts.hs +484/−313
• Math/Combinatorics/IncidenceAlgebra.hs +5/−3
• Math/Combinatorics/Matroid.hs +3/−1
• Math/Combinatorics/Poset.hs +31/−24
• Math/Combinatorics/StronglyRegularGraph.hs +3/−1
• Math/CommutativeAlgebra/Polynomial.hs +3/−1
• Math/Core/Field.hs +8/−7
• Math/Core/Utils.hs +47/−1
• Math/NumberTheory/Factor.hs +62/−34
• Math/NumberTheory/Prime.hs +18/−8
• Math/NumberTheory/QuadraticField.hs +2/−2
• Math/Projects/ChevalleyGroup/Classical.hs +3/−1
• Math/Projects/ChevalleyGroup/Exceptional.hs +4/−2
• Math/Projects/MiniquaternionGeometry.hs +3/−1
• Math/Projects/RootSystem.hs +3/−1
• Math/QuantumAlgebra/OrientedTangle.hs +9/−7
• Math/QuantumAlgebra/QuantumPlane.hs +4/−4
• Math/QuantumAlgebra/Tangle.hs +6/−6
• Math/QuantumAlgebra/TensorCategory.hs +138/−43
• Math/Test/TAlgebras/TOctonions.hs +3/−1
• Math/Test/TAlgebras/TStructures.hs +2/−0
• Math/Test/TAlgebras/TTensorProduct.hs +3/−1
• Math/Test/TAlgebras/TVectorSpace.hs +15/−9
• Math/Test/TCombinatorics/TGraphAuts.hs +134/−1
• Math/Test/TNumberTheory/TPrimeFactor.hs +6/−2

@@ -1,5 +1,5 @@    Name:                HaskellForMaths
-   Version:             0.4.5
+   Version:             0.4.6
    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 as an educational resource, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra

@@ -1,9 +1,11 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) David Amos, 2008-2015. All rights reserved.
 
 {-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, ScopedTypeVariables, EmptyDataDecls, FlexibleInstances #-}
 
 module Math.Algebra.Field.Extension where
 
+import Prelude hiding ( (<*>) )
+
 import Data.Ratio
 import Data.List as L (elemIndex)
 
@@ -41,9 +43,9 @@                         | i > 1  = v ++ "^" ++ show i -- "x^" ++ show i
 
 instance (Eq a, Num a) => Num (UPoly a) where
-    UP as + UP bs = toUPoly $ as <+> bs
+    UP as + UP bs = UP $ as <+> bs
     negate (UP as) = UP $ map negate as
-    UP as * UP bs = toUPoly $ as <*> bs
+    UP as * UP bs = UP $ as <*> bs
     fromInteger 0 = UP []
     fromInteger a = UP [fromInteger a]
     abs _ = error "Prelude.Num.abs: inappropriate abstraction"
@@ -51,14 +53,30 @@ 
 toUPoly as = UP (reverse (dropWhile (== 0) (reverse as)))
 
-(a:as) <+> (b:bs) = (a+b) : (as <+> bs)
+-- The fussiness of the code is to avoid adding trailing zeroes, eg [3] <+> [-3]
+-- Otherwise we would have to normalise after every addition
 as <+> [] = as
 [] <+> bs = bs
+-- (a:as) <+> (b:bs) = (a+b) : (as <+> bs)
+(a:as) <+> (b:bs) = let c = a+b
+                        cs = as <+> bs
+                    in if c == 0 && null cs then [] else c:cs
 
+-- The fussiness of the code is to avoid adding trailing zeroes.
+-- Note that since we call <+>, we rely on it having similar properties.
 [] <*> _ = []
+_ <*> [] = [] -- to avoid [0,1] <*> [] -> [0]
+(a:as) <*> bs = if null as then map (a*) bs else map (a*) bs <+> (0 : as <*> bs)
+
+-- > let valid xs = null xs || last xs /= 0
+-- > quickCheck (\as bs -> not (valid as) || not (valid bs) || valid (as <*> bs))
+
+{-
+-- The following definition introduces unnecessary trailing zeroes, eg [3] <*> [2] -> [6,0]
+[] <*> _ = []
 _ <*> [] = []
 (a:as) <*> (b:bs) = [a*b] <+> (0 : map (a*) bs) <+> (0 : map (*b) as) <+> (0 : 0 : as <*> bs)
-
+-}
 
 convert (UP as) = toUPoly $ map fromInteger as
 -- Can be used with type annotations to construct polynomials over other types, eg

@@ -31,6 +31,9 @@ -- |A type for permutations, considered as functions or actions which can be performed on an underlying set.
 newtype Permutation a = P (M.Map a a) deriving (Eq,Ord)
 
+-- Can't make a Functor instance because we need an Ord instance
+fmapP f = fromPairs . map (\(x,y) -> (f x, f y)) . toPairs
+
 -- |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
@@ -69,7 +72,8 @@ xs -^ g = L.sort [x .^ g | x <- xs]
 
 -- construct a permutation from cycles
-fromCycles cs = fromPairs $ concatMap fromCycle cs
+-- fromCycles cs = fromPairs $ concatMap fromCycle cs
+fromCycles cs = product $ map (fromPairs . fromCycle) cs
     where fromCycle xs = zip xs (rotateL xs)
 
 -- convert a permutation to cycles
@@ -94,7 +98,7 @@ 
 -- |The Num instance is what enables us to write @g*h@ for the product of group elements and @1@ for the group identity.
 -- Unfortunately we can't of course give sensible definitions for the other functions declared in the Num typeclass.
-instance (Ord a, Show a) => Num (Permutation a) where
+instance Ord a => Num (Permutation a) where
     g * h = fromPairs' [(x, x .^ g .^ h) | x <- supp g `union` supp h]
     -- signum = sign -- doesn't work, complains about no (+) instance
     fromInteger 1 = P $ M.empty
@@ -104,13 +108,13 @@     signum _ = error "(Permutation a).signum: not applicable"
 
 -- |The HasInverses instance is what enables us to write @g^-1@ for the inverse of a group element.
-instance (Ord a, Show a) => HasInverses (Permutation a) where
+instance Ord a => HasInverses (Permutation a) where
     inverse (P g) = P $ M.fromList $ map (\(x,y)->(y,x)) $ M.toList g
 
 
 -- |g ~^ h returns the conjugate of g by h, that is, h^-1*g*h.
 -- The tilde is meant to a mnemonic, because conjugacy is an equivalence relation.
-(~^) :: (Ord a, Show a) => Permutation a -> Permutation a -> Permutation a
+(~^) :: Ord a => Permutation a -> Permutation a -> Permutation a
 g ~^ h = h^-1 * g * h
 
 -- commutator
@@ -119,6 +123,7 @@ 
 -- ORBITS
 
+{-
 closureS xs fs = closure' S.empty (S.fromList xs) where
     closure' interior boundary
         | S.null boundary = interior
@@ -126,6 +131,12 @@             let interior' = S.union interior boundary
                 boundary' = S.fromList [f x | x <- S.toList boundary, f <- fs] S.\\ interior'
             in closure' interior' boundary'
+-}
+closureS xs fs = closure' S.empty xs where
+    closure' interior (x:xs)
+        | S.member x interior = closure' interior xs
+        | otherwise = closure' (S.insert x interior) ([f x | f <- fs] ++ xs)
+    closure' interior [] = interior
 
 closure xs fs = S.toList $ closureS xs fs
 
@@ -293,11 +304,19 @@ -- For example, sgs (_A 5) == [[[1,2,3]],[[2,4,5]],[[3,4,5]]]
 -- So we need all three to generate the first transversal, then the last two to generate the second transversal, etc
 
--- |Given a strong generating set, return the order of the group it generates
+-- |Given a strong generating set, return the order of the group it generates.
+-- Note that the SGS is assumed to be relative to the natural order of the points on which the group acts.
 orderSGS :: (Ord a) => [Permutation a] -> Integer
 orderSGS sgs = product $ map (L.genericLength . fundamentalOrbit) bs where
     bs = toListSet $ map minsupp sgs
     fundamentalOrbit b = b .^^ filter ( (b <=) . minsupp ) sgs
+
+-- !! Needs more testing
+-- |Given a base and strong generating set, return the order of the group it generates.
+orderBSGS :: (Ord a) => ([a],[Permutation a]) -> Integer
+orderBSGS (bs,sgs) = go 1 bs sgs where
+    go n [] _ = n
+    go n (b:bs) gs = go (n * L.genericLength (b .^^ gs)) bs (filter (\g -> b .^ g == b) gs)
 
 
 -- MORE INVESTIGATIONS

@@ -6,7 +6,7 @@ import Data.Maybe (isNothing, isJust)
 import qualified Data.Set as S
 import qualified Data.Map as M
-import Math.Algebra.Group.PermutationGroup hiding (elts, order, gens, isMember, isSubgp, isNormal, reduceGens, normalClosure, commutatorGp, derivedSubgp)
+import Math.Algebra.Group.PermutationGroup hiding (elts, order, orderBSGS, gens, isMember, isSubgp, isNormal, reduceGens, normalClosure, commutatorGp, derivedSubgp)
 import Math.Common.ListSet (toListSet)
 import Math.Core.Utils hiding (elts)
 

@@ -1,4 +1,4 @@--- Copyright (c) 2008-2012, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, David Amos. All rights reserved.
 
 -- |A module providing elementary operations involving scalars, vectors, and matrices
 -- over a ring or field. Vectors are represented as [a], matrices as [[a]].
@@ -9,6 +9,8 @@ -- on each side indicates the dimension of the argument on that side. For example,
 -- v \<*\>\> m is multiplication of a vector on the left by a matrix on the right.
 module Math.Algebra.LinearAlgebra where
+
+import Prelude hiding ( (*>), (<*>) )
 
 import qualified Data.List as L
 import Math.Core.Field -- not actually used in this module

@@ -1,4 +1,4 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} @@ -38,7 +38,7 @@ instance Show v => Show (SL2 v) where show (SL2 m) = show m  instance Algebra Q (SL2 ABCD) where -- to do this for Num k instead of Q we would need a,b,c,d defined for Num k-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)] where munit = SL2 (Glex 0 [])     mult x = x''' where         x' = mult $ fmap ( \(SL2 a, SL2 b) -> (a,b) ) x -- perform the multiplication in GlexPoly

@@ -1,10 +1,16 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} --- |A module defining the algebra of commutative polynomials over a field k+-- |A module defining the algebra of commutative polynomials over a field k.+--+-- Most users should probably use Math.CommutativeAlgebra.Polynomial instead, which is basically the same thing+-- but more fully-featured. This module will probably be deprecated at some point, but remains for now because+-- it has a simpler implementation which may be more helpful for people wanting to understand the code. module Math.Algebras.Commutative where +import Prelude hiding ( (*>) )+ import Math.Algebra.Field.Base hiding (powers) import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct@@ -35,6 +41,10 @@ instance Functor GlexMonomial where     fmap f (Glex si xis) = Glex si [(f x, i) | (x,i) <- xis] -- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards++instance Applicative GlexMonomial where+    pure = return+    (<*>) = ap  -- GlexMonomial is the free commutative monoid, and hence a monad instance Monad GlexMonomial where

@@ -1,4 +1,4 @@--- Copyright (c) 2010-2012, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances #-} -- ScopedTypeVariables@@ -13,6 +13,8 @@ -- Elements of the group algebra consist of arbitrary K-linear combinations of elements of G. -- For example, @p [[1,2,3]] + 2 * p [[1,2],[3,4]]@ module Math.Algebras.GroupAlgebra (GroupAlgebra, p) where++import Prelude hiding ( (*>) )  import Math.Core.Field import Math.Core.Utils hiding (elts)

@@ -32,7 +32,7 @@     mmult (LM si xis) (LM sj yjs) = LM (si+sj) $ addmerge xis yjs  instance (Eq k, Num k) => Algebra k LaurentMonomial where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)]      mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts)     -- mult (V ts) = nf $ V [(a `mmult` b, x) | (T a b, x) <- ts]

@@ -1,10 +1,12 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances #-}   module Math.Algebras.Matrix where +import Prelude hiding ( (*>) )+ import Math.Algebra.Field.Base import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct@@ -45,9 +47,9 @@ toMat2 [[a,b],[c,d]] = sum $ zipWith (\x e -> unit x * return e) [a,b,c,d] [E2 1 1, E2 1 2, E2 2 1, E2 2 2] -- fromMat2 -toEB2 [x,y] = foldl add zero $ zipWith (\x e -> x `smultL` return e) [x,y] [E 1, E 2]+toEB2 [x,y] = foldl add zerov $ zipWith (\x e -> x `smultL` return e) [x,y] [E 1, E 2] -toEB xs = foldl add zero $ zipWith (\x e -> x `smultL` return e) xs (map E [1..])+toEB xs = foldl add zerov $ zipWith (\x e -> x `smultL` return e) xs (map E [1..])   @@ -58,7 +60,7 @@ instance (Eq k, Num k) => Coalgebra k Mat2' where     counit (V ts) = sum [xij * delta i j | (E2' i j, xij) <- ts]     -- comult (V ts) = V $ concatMap (\(E2' i j,xij) -> [(T (E2' i k) (E2' k j), xij) | k <- [1..2]]) ts-    comult = linear (\(E2' i j) -> foldl (<+>) zero [return (E2' i k, E2' k j) | k <- [1..2]])+    comult = linear (\(E2' i j) -> foldl (<+>) zerov [return (E2' i k, E2' k j) | k <- [1..2]]) -- In other words -- counit (a b) = (1 0) --        (c d)   (0 1)@@ -74,7 +76,7 @@ -- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere  instance (Eq k, Num k) => Algebra k M3 where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(E3 i i, x) | i <- [1..3] ]     -- mult (V ts) = nf $ V $ map (\((E3 i j, E3 k l), x) -> (E3 i l, delta j k * x)) ts     mult = linear mult' where

@@ -1,10 +1,12 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}  -- |A module defining the algebra of non-commutative polynomials over a field k module Math.Algebras.NonCommutative where +import Prelude hiding ( (*>) )+ import Math.Algebra.Field.Base hiding (powers) import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct@@ -30,7 +32,7 @@     mmult (NCM i xs) (NCM j ys) = NCM (i+j) (xs++ys)  instance (Eq k, Num k, Ord v) => Algebra k (NonComMonomial v) where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)]     mult = nf . fmap (\(a,b) -> a `mmult` b) 

@@ -1,4 +1,4 @@--- Copyright (c) 2011, David Amos. All rights reserved.+-- Copyright (c) 2011-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances, NoMonomorphismRestriction #-} @@ -7,6 +7,8 @@ -- The octonions are the algebra defined by the basis {1,i0,i1,i2,i3,i4,i5,i6}, -- where each i_n * i_n = -1, and i_n+1 * i_n+2 = i_n+4 (where the indices are modulo 7). module Math.Algebras.Octonions where++import Prelude hiding ( (*>) )  import Math.Core.Field import Math.Algebras.VectorSpace

@@ -1,4 +1,4 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances, NoMonomorphismRestriction #-} @@ -7,6 +7,8 @@ -- The quaternions are the algebra defined by the basis {1,i,j,k}, where i^2 = j^2 = k^2 = ijk = -1 module Math.Algebras.Quaternions where +import Prelude hiding ( (*>) )+ import Math.Core.Field import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct@@ -68,7 +70,7 @@     sqnorm :: Vect k a -> k  -- |If an algebra has a conjugation operation, then it has multiplicative inverses,--- via 1/x = conj x / sqnorm x+-- via 1\/x = conj x \/ sqnorm x instance (Eq k, Fractional k, Ord a, Show a, HasConjugation k a) => Fractional (Vect k a) where     recip 0 = error "recip 0"     recip x = (1 / sqnorm x) *> conj x@@ -167,7 +169,7 @@ instance (Eq k, Num k) => Coalgebra k (Dual HBasis) where     counit = unwrap . linear counit'         where counit' (Dual One) = return ()-              counit' _ = zero+              counit' _ = zerov     comult = linear comult'         where comult' (Dual One) = return (Dual One, Dual One) <+>                   (-1) *> ( return (Dual I, Dual I) <+> return (Dual J, Dual J) <+> return (Dual K, Dual K) )

@@ -1,4 +1,4 @@--- Copyright (c) David Amos, 2010-2012. All rights reserved.+-- Copyright (c) David Amos, 2010-2015. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction #-} {-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}@@ -8,6 +8,8 @@ -- - specifically algebra, coalgebra, bialgebra, Hopf algebra, module, comodule module Math.Algebras.Structures where +import Prelude hiding ( (*>) )+ import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct @@ -103,7 +105,7 @@     mult = linear mult'         where mult' (Left a1, Left a2) = i1 $ mult $ return (a1,a2)               mult' (Right b1, Right b2) = i2 $ mult $ return (b1,b2)-              mult' _ = zero+              mult' _ = zerov -- This is the product algebra, which is the product in the category of algebras -- 1 = (1,1) -- (a1,b1) * (a2,b2) = (a1*a2, b1*b2)

@@ -5,6 +5,8 @@ -- |A module defining the tensor algebra, symmetric algebra, exterior (or alternating) algebra, and tensor coalgebra module Math.Algebras.TensorAlgebra where +import Prelude hiding ( (*>) )+ import qualified Data.List as L  import Math.Algebras.VectorSpace@@ -195,7 +197,7 @@     mult xy = nf $ xy >>= (\(Ext i xs, Ext j ys) -> signedMerge 1 (0,[]) (i,xs) (j,ys))         where signedMerge s (k,zs) (i,x:xs) (j,y:ys) =                   case compare x y of-                  EQ -> zero+                  EQ -> zerov                   LT -> signedMerge s (k+1,x:zs) (i-1,xs) (j,y:ys)                   GT -> let s' = if even i then s else -s -- we had to commute y past x:xs, with i sign changes                         in signedMerge s' (k+1,y:zs) (i,x:xs) (j-1,ys)

@@ -1,10 +1,12 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE NoMonomorphismRestriction #-}  -- |A module defining direct sum and tensor product of vector spaces module Math.Algebras.TensorProduct where +import Prelude hiding ( (*>) )+ import Math.Algebras.VectorSpace  infix 7 `te`, `tf`@@ -38,12 +40,12 @@ p1 :: (Eq k, Num k, Ord a) => Vect k (DSum a b) -> Vect k a p1 = linear p1' where     p1' (Left a) = return a-    p1' (Right b) = zero+    p1' (Right b) = zerov  -- |Projection onto right summand from direct sum p2 :: (Eq k, Num k, Ord b) => Vect k (DSum a b) -> Vect k b p2 = linear p2' where-    p2' (Left a) = zero+    p2' (Left a) = zerov     p2' (Right b) = return b  -- |The product of two linear functions (with the same source).@@ -83,7 +85,7 @@ tf :: (Eq k, Num k, Ord a', Ord b') => (Vect k a -> Vect k a') -> (Vect k b -> Vect k b')    -> Vect k (Tensor a b) -> Vect k (Tensor a' b') tf f g (V ts) = sum [x *> te (f $ return a) (g $ return b) | ((a,b), x) <- ts]-    where sum = foldl add zero -- (V [])+    where sum = foldl add zerov   -- tensor isomorphisms

@@ -1,4 +1,4 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE NoMonomorphismRestriction #-} {-# OPTIONS_HADDOCK prune #-}@@ -6,12 +6,16 @@ -- |A module defining the type and operations of free k-vector spaces over a basis b (for a field k) module Math.Algebras.VectorSpace where +import Prelude hiding ( (*>) )++import Control.Applicative hiding ( (<*), (*>) )+import Control.Monad (ap) import qualified Data.List as L import qualified Data.Set as S -- only needed for toSet  infixr 7 *> infixl 7 <*-infixl 6 <+>, <->+infixl 6 <+>, <->, <<+>>, <<->>   -- |Given a field type k and a basis type b, Vect k b is the type of the free k-vector space over b.@@ -48,24 +52,24 @@ terms (V ts) = ts  -- |Return the coefficient of the specified basis element in a vector+coeff :: (Num k, Eq b) => b -> Vect k b -> k coeff b v = sum [k | (b',k) <- terms v, b' == b]  -- |Remove the term for a specified basis element from a vector-removeTerm b v = v <-> coeff b v *> return b---- Deprecated-zero = V []+removeTerm :: (Eq k, Num k, Ord b) => b -> Vect k b -> Vect k b+removeTerm b (V ts) = V $ filter ((/=b) . fst) ts+-- v <-> coeff b v *> return b  -- |The zero vector zerov :: Vect k b zerov = V []  -- |Addition of vectors-add :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b+add :: (Eq k, Num k, Ord b) => Vect k b -> Vect k b -> Vect k b add (V ts) (V us) = V $ addmerge ts us  -- |Addition of vectors (same as add)-(<+>) :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b+(<+>) :: (Eq k, Num k, Ord b) => Vect k b -> Vect k b -> Vect k b (<+>) = add  addmerge ((a,x):ts) ((b,y):us) =@@ -77,23 +81,20 @@ addmerge [] us = us  -- |Sum of a list of vectors-sumv :: (Ord b, Eq k, Num k) => [Vect k b] -> Vect k b+sumv :: (Eq k, Num k, Ord b) => [Vect k b] -> Vect k b sumv = foldl (<+>) zerov --- Deprecated-neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts- -- |Negation of a vector negatev :: (Eq k, Num k) => Vect k b -> Vect k b negatev (V ts) = V $ map (\(b,x) -> (b,-x)) ts  -- |Subtraction of vectors-(<->) :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b -> Vect k b+(<->) :: (Eq k, Num k, Ord b) => Vect k b -> Vect k b -> Vect k b (<->) u v = u <+> negatev v  -- |Scalar multiplication (on the left) smultL :: (Eq k, Num k) => k -> Vect k b -> Vect k b-smultL 0 _ = zero -- V []+smultL 0 _ = zerov -- V [] smultL k (V ts) = V [(ei,k*xi) | (ei,xi) <- ts]  -- |Same as smultL. Mnemonic is \"multiply through (from the left)\"@@ -102,7 +103,7 @@  -- |Scalar multiplication on the right smultR :: (Eq k, Num k) => Vect k b -> k -> Vect k b-smultR _ 0 = zero -- V []+smultR _ 0 = zerov -- V [] smultR (V ts) k = V [(ei,xi*k) | (ei,xi) <- ts]  -- |Same as smultR. Mnemonic is \"multiply through (from the right)\"@@ -119,7 +120,7 @@  -- |Convert an element of Vect k b into normal form. Normal form consists in having the basis elements in ascending order, -- with no duplicates, and all coefficients non-zero-nf :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b+nf :: (Eq k, Num k, Ord b) => Vect k b -> Vect k b nf (V ts) = V $ nf' $ L.sortBy compareFst ts where     nf' ((b1,x1):(b2,x2):ts) =         case compare b1 b2 of@@ -136,7 +137,7 @@ -- -- In the mathematical sense, this can be regarded as a functor from the category Set (of sets) to the category k-Vect -- (of k-vector spaces). In Haskell, instead of Set we have Hask, the category of Haskell types. However, for our purposes--- it is helpful to identify Hask with Set, but identifying a Haskell type with its set of inhabitants.+-- it is helpful to identify Hask with Set, by identifying a Haskell type with its set of inhabitants. -- -- The type constructor (Vect k) gives the action of the functor on objects in the category, -- taking a set (type) to a free k-vector space. fmap gives the action of the functor on arrows in the category,@@ -149,6 +150,16 @@     fmap f (V ts) = V [(f b, x) | (b,x) <- ts] -- Note that if f is not order-preserving, then we need to call "nf" afterwards +-- From GHC 7.10, Monad has Applicative as a superclass, so we must define an instance.+-- It doesn't particularly make sense for Vect k.+-- (Although given Vect k b, we could represent the dual space as Vect k (b -> ()),+-- and then have a use for <*>.)+instance Num k => Applicative (Vect k) where+    pure = return+    -- pure b = V [(b,1)]+    (<*>) = ap+    -- V fs <*> V xs = V [(f x, a*b) | (f,a) <- fs, (x,b) <- xs]+ -- |Given a field k, the type constructor (Vect k) is a monad, the \"free k-vector space monad\". -- -- In order to understand this, it is probably easiest to think of a free k-vector space as a kind of container,@@ -172,7 +183,7 @@ -- -- If we have A = Vect k a, B = Vect k b, and f :: a -> Vect k b is a function from the basis elements of A into B, -- then @linear f@ is the linear map that this defines by linearity.-linear :: (Ord b, Eq k, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b+linear :: (Eq k, Num k, Ord b) => (a -> Vect k b) -> Vect k a -> Vect k b linear f v = nf $ v >>= f  newtype EBasis = E Int deriving (Eq,Ord)@@ -193,7 +204,7 @@  -- |Wrap an element of the field k to an element of the trivial k-vector space wrap :: (Eq k, Num k) => k -> Vect k ()-wrap 0 = zero+wrap 0 = zerov wrap x = V [( (),x)]  -- |Unwrap an element of the trivial k-vector space to an element of the field k@@ -223,6 +234,20 @@  (f <<+>> g) v = f v <+> g v +(f <<->> g) v = f v <-> g v+ zerof v = zerov  sumf fs = foldl (<<+>>) zerof fs+++-- Lens+coeffLens :: (Ord b, Eq k, Num k, Functor f) => b -> (k -> f k) -> (Vect k b -> f (Vect k b))+coeffLens b = lens (coeff b) (setter b)+    where setter b = \(V ts) k -> (k *> return b) <+> (V $ filter ((/=b) . fst) ts)+          lens getter setter f a = fmap (setter a) (f (getter a))+-- Can be used with lens-family, for example+-- e1 ^. coeffLens (E 2) --> 0+-- e1 & coeffLens (E 2) .~ 2 --> e1+2e2+-- e1 & coeffLens (E 1) %~ (+2) --> 3e1+

@@ -1,4 +1,4 @@--- Copyright (c) 2012, David Amos. All rights reserved.+-- Copyright (c) 2012-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction, ScopedTypeVariables, DeriveFunctor #-} @@ -34,6 +34,8 @@ -- Lie Algebras and Hopf Algebras -- Michiel Hazewinkel, Nadiya Gubareni, V.V.Kirichenko +import Prelude hiding ( (*>) )+ import Data.List as L import Data.Maybe (fromJust) import qualified Data.Set as S@@ -198,17 +200,17 @@  -- |Convert an element of SSym represented in the monomial basis to the fundamental basis ssymMtoF :: (Eq k, Num k) => Vect k SSymM -> Vect k SSymF-ssymMtoF = linear ssymMtoF'-    where ssymMtoF' (SSymM u) = sumv [mu (set,po) u v *> return (SSymF v) | v <- set, po u v]-              where set = L.permutations u-                    po = weakOrder+ssymMtoF = linear ssymMtoF' where+    ssymMtoF' (SSymM u) = sumv [mu (set,po) u v *> return (SSymF v) | v <- set, po u v]+        where set = L.permutations u+              po = weakOrder  -- |Convert an element of SSym represented in the fundamental basis to the monomial basis ssymFtoM :: (Eq k, Num k) => Vect k SSymF -> Vect k SSymM-ssymFtoM = linear ssymFtoM'-    where ssymFtoM' (SSymF u) = sumv [return (SSymM v) | v <- set, po u v]-              where set = L.permutations u-                    po = weakOrder+ssymFtoM = linear ssymFtoM' where+    ssymFtoM' (SSymF u) = sumv [return (SSymM v) | v <- set, po u v]+        where set = L.permutations u+              po = weakOrder  -- (p,q)-shuffles: permutations of [1..p+q] having at most one descent, at position p -- denoted S^{(p,q)} in Aguiar&Sottile@@ -457,16 +459,16 @@  -- |Convert an element of YSym represented in the monomial basis to the fundamental basis ysymMtoF :: (Eq k, Num k) => Vect k YSymM -> Vect k (YSymF ())-ysymMtoF = linear ysymMtoF'-    where ysymMtoF' (YSymM t) = sumv [mu (set,po) t s *> return (YSymF s) | s <- set]-              where po = tamariOrder-                    set = tamariUpSet t -- [s | s <- trees (nodecount t), t `tamariOrder` s]+ysymMtoF = linear ysymMtoF' where+    ysymMtoF' (YSymM t) = sumv [mu (set,po) t s *> return (YSymF s) | s <- set]+        where po = tamariOrder+              set = tamariUpSet t -- [s | s <- trees (nodecount t), t `tamariOrder` s]  -- |Convert an element of YSym represented in the fundamental basis to the monomial basis ysymFtoM :: (Eq k, Num k) => Vect k (YSymF ()) -> Vect k YSymM-ysymFtoM = linear ysymFtoM'-    where ysymFtoM' (YSymF t) = sumv [return (YSymM s) | s <- tamariUpSet t]-                            -- sumv [return (YSymM s) | s <- trees (nodecount t), t `tamariOrder` s]+ysymFtoM = linear ysymFtoM' where+    ysymFtoM' (YSymF t) = sumv [return (YSymM s) | s <- tamariUpSet t]+                       -- sumv [return (YSymM s) | s <- trees (nodecount t), t `tamariOrder` s]   instance (Eq k, Num k) => Algebra k YSymM where@@ -565,13 +567,13 @@  -- |Convert an element of QSym represented in the monomial basis to the fundamental basis qsymMtoF :: (Eq k, Num k) => Vect k QSymM -> Vect k QSymF-qsymMtoF = linear qsymMtoF'-    where qsymMtoF' (QSymM alpha) = sumv [(-1) ^ (length beta - length alpha) *> return (QSymF beta) | beta <- refinements alpha]+qsymMtoF = linear qsymMtoF' where+    qsymMtoF' (QSymM alpha) = sumv [(-1) ^ (length beta - length alpha) *> return (QSymF beta) | beta <- refinements alpha]  -- |Convert an element of QSym represented in the fundamental basis to the monomial basis qsymFtoM :: (Eq k, Num k) => Vect k QSymF -> Vect k QSymM-qsymFtoM = linear qsymFtoM'-    where qsymFtoM' (QSymF alpha) = sumv [return (QSymM beta) | beta <- refinements alpha] -- ie beta <- up-set of alpha+qsymFtoM = linear qsymFtoM' where+    qsymFtoM' (QSymF alpha) = sumv [return (QSymM beta) | beta <- refinements alpha] -- ie beta <- up-set of alpha  instance (Eq k, Num k) => Algebra k QSymF where     unit x = x *> return (QSymF [])@@ -592,13 +594,13 @@ -- the above induces Hopf algebra structure on quasi-symmetric functions via -- m_alpha -> sum [product (zipWith (^) (map x_ is) alpha | is <- combinationsOf k [] ] where k = length alpha -xvars n = [glexvar ("x" ++ show i) | i <- [1..n] ]+-- xvars n = [glexvar ("x" ++ show i) | i <- [1..n] ]  -- |@qsymPoly n is@ is the quasi-symmetric polynomial in n variables for the indices is. (This corresponds to the -- monomial basis for QSym.) For example, qsymPoly 3 [2,1] == x1^2*x2+x1^2*x3+x2^2*x3. qsymPoly :: Int -> [Int] -> GlexPoly Q String qsymPoly n is = sum [product (zipWith (^) xs' is) | xs' <- combinationsOf r xs]-    where xs = xvars n+    where xs = [glexvar ("x" ++ show i) | i <- [1..n] ]           r = length is  @@ -776,6 +778,7 @@  -- "inverse" for descendingTree -- These are the maps called gamma in Loday.pdf+-- or are they? - these give the min and max inverse images in the lexicographic order, rather than the weak order? minPerm t = minPerm' (lrCountTree t)     where minPerm' E = []           minPerm' (T l (lc,rc) r) = minPerm' l ++ [lc+rc+1] ++ map (+lc) (minPerm' r)@@ -810,10 +813,9 @@  -- The composition of [1..n] obtained by treating each descent as a cut descentComposition [] = []-descentComposition xs = dc $ zipWith (>) xs (tail xs) ++ [False]-    where dc bs = case break id bs of-                  (ls,r:rs) -> (length ls + 1) : dc rs-                  (ls,[]) -> [length ls]+descentComposition xs = descComp 0 xs where+    descComp c (x1:x2:xs) = if x1 < x2 then descComp (c+1) (x2:xs) else (c+1) : descComp 0 (x2:xs)+    descComp c [x] = [c+1]  -- |Given a permutation of [1..n], its descents are those positions where the next number is less than the previous number. -- For example, the permutation [2,3,5,1,6,4] has descents from 5 to 1 and from 6 to 4. The descents can be regarded as cutting

@@ -11,7 +11,7 @@ import qualified Data.Map as M import qualified Data.Set as S -toSet = S.toList . S.fromList+import Math.Core.Utils (picks, toSet)  -- |A digraph is represented as DG vs es, where vs is the list of vertices, and es is the list of edges. -- Edges are directed: an edge (u,v) means an edge from u to v.@@ -50,8 +50,8 @@     | otherwise = digraphIsos' [] vsa vsb     where digraphIsos' xys [] [] = [xys]           digraphIsos' xys (x:xs) ys =-              concat [ digraphIsos' ((x,y):xys) xs (L.delete y ys)-                     | y <- ys, isCompatible (x,y) xys]+              concat [ digraphIsos' ((x,y):xys) xs ys'+                     | (y,ys') <- picks ys, isCompatible (x,y) xys]           isCompatible (x,y) xys = and [ ((x,x') `elem` esa) == ((y,y') `elem` esb)                                       && ((x',x) `elem` esa) == ((y',y) `elem` esb)                                        | (x',y') <- xys ]@@ -76,8 +76,8 @@                                          (x',y') <- xys]           dfs xys [] [] = [xys]           dfs xys (x:xs) ys =-              concat [ dfs ((x,y):xys) xs (L.delete y ys)-                     | y <- ys, isCompatible (x,y) xys]+              concat [ dfs ((x,y):xys) xs ys'+                     | (y,ys') <- picks ys, isCompatible (x,y) xys]  -- For DAGs, can almost certainly do better than the above by using the height partition -- However see remarks in Poset on orderIsos:@@ -111,8 +111,8 @@           dfs xys [] [] = [xys]           dfs xys ([]:las) ([]:lbs) = dfs xys las lbs           dfs xys ((x:xs):las) (ys:lbs) =-              concat [ dfs ((x,y):xys) (xs:las) (L.delete y ys : lbs)-                     | y <- ys, isCompatible (x,y) xys]+              concat [ dfs ((x,y):xys) (xs:las) (ys' : lbs)+                     | (y,ys') <- picks ys, isCompatible (x,y) xys]           isCompatible (x,y) xys =               let preds_x = M.findWithDefault [] x predsA                   preds_y = M.findWithDefault [] y predsB@@ -169,7 +169,7 @@           dfs xys [] [] = [xys]           dfs xys ([]:sls) ([]:tls) = dfs xys sls tls           dfs xys ((x:xs):sls) (ys:tls) =-              concat [ dfs ((x,y):xys) (xs:sls) (L.delete y ys : tls) | y <- ys]+              concat [ dfs ((x,y):xys) (xs:sls) (ys' : tls) | (y,ys') <- picks ys]               -- not applying any compatibility condition yet  

@@ -1,9 +1,11 @@--- Copyright (c) David Amos, 2008-2011. All rights reserved.
+-- Copyright (c) David Amos, 2008-2015. All rights reserved.
 
 -- |Constructions of the finite geometries AG(n,Fq) and PG(n,Fq), their points, lines and flats,
 -- together with the incidence graphs between points and lines.
 module Math.Combinatorics.FiniteGeometry where
 
+import Prelude hiding ( (*>) )
+
 import Data.List as L
 import qualified Data.Set as S
 
@@ -189,7 +191,7 @@     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
+-- incidenceAuts (incidenceGraphPG n fq) == PGL(n+1,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),

@@ -14,7 +14,7 @@ import Math.Common.ListSet as LS
 import Math.Core.Utils
 import Math.Algebra.Group.PermutationGroup hiding (fromDigits, fromBinary)
-import Math.Algebra.Group.SchreierSims as SS
+import qualified Math.Algebra.Group.SchreierSims as SS
 
 -- Main source: Godsil & Royle, Algebraic Graph Theory
 
@@ -105,7 +105,7 @@ 
 -- |c n is the cyclic graph on n vertices
 c :: (Integral t) => t -> Graph t
-c n = graph (vs,es) where
+c n | n >= 3 = graph (vs,es) where
     vs = [1..n]
     es = L.insert [1,n] [[i,i+1] | i <- [1..n-1]]
 -- automorphism group is D2n
@@ -165,7 +165,12 @@            [2,7],[7,3],[3,8],[8,4],[4,9],[9,5],[5,10],[10,6],[6,11],[11,2] ]
 
 
+-- Prisms are regular, vertex-transitive, but not edge-transitive unless n == 1, 2, 4.
+-- (prism 2 ~= q 2, prism 4 ~= q 3)
+prism :: Int -> Graph (Int,Int)
+prism n = k 2 `cartProd` c n
 
+
 -- convert a graph to have [1..n] as vertices
 to1n (G vs es) = graph (vs',es') where
     mapping = M.fromList $ zip vs [1..] -- the mapping from vs to [1..n]
@@ -223,11 +228,22 @@ lineGraph' (G vs es) = graph (es, [ [ei,ej] | ei <- es, ej <- dropWhile (<= ei) es, ei `intersect` ej /= [] ])
 
 
+-- For example cartProd (c m) (c n) is a wireframe for a torus
+-- cartProd (q m) (q n) `isGraphIso` q (m+n)
+-- Godsil and Royle p154
+cartProd (G vs es) (G vs' es') = G us [e | e@[u,u'] <- combinationsOf 2 us, u `adj` u' ]
+    where us = [(v,v') | v <- vs, v' <- vs']
+          eset = S.fromList es
+          eset' = S.fromList es'
+          adj (x1,y1) (x2,y2) = x1 == x2 && L.sort [y1,y2] `S.member` eset'
+                             || y1 == y2 && L.sort [x1,x2] `S.member` eset
+
+
 -- SIMPLE PROPERTIES OF GRAPHS
 
-order g = length (vertices g)
+order = length . vertices
 
-size g = length (edges g)
+size = length . edges
 
 -- also called degree
 valency (G vs es) v = length $ filter (v `elem`) es
@@ -297,16 +313,24 @@ 
 -- circumference = max cycle - Bollobas p104
 
+-- Vertices that are not in the same component as the start vertex all go into a final cell
+distancePartition g@(G vs es) v = distancePartitionS vs (S.fromList es) v
 
-distancePartition g v = distancePartition' S.empty (S.singleton v) where
-    distancePartition' interior boundary
-        | S.null boundary = []
-        | otherwise = let interior' = S.union interior boundary
-                          boundary' = foldl S.union S.empty [S.fromList (nbrs g x) | x <- S.toList boundary] S.\\ interior'
-                      in S.toList boundary : distancePartition' interior' boundary'
+distancePartitionS vs eset v = distancePartition' (S.singleton v) (S.delete v (S.fromList vs)) where
+    distancePartition' boundary exterior
+        | S.null boundary = if S.null exterior then [] else [S.toList exterior] -- graph may not have been connected
+        | otherwise = let (boundary', exterior') = S.partition (\v -> any (`S.member` eset) [L.sort [u,v] | u <- S.toList boundary]) exterior
+                      in S.toList boundary : distancePartition' boundary' exterior'
 
 -- the connected component to which v belongs
-component g v = L.sort $ concat $ distancePartition g v
+component g v = component' S.empty (S.singleton v) where
+    component' interior boundary
+        | S.null boundary = S.toList interior
+        | otherwise = let interior' = S.union interior boundary
+                          boundary' = foldl S.union S.empty [S.fromList (nbrs g x) | x <- S.toList boundary] S.\\ interior'
+                      in component' interior' boundary'
+-- TODO: This can almost certainly be made more efficient.
+-- nbrs is O(n), and this calls it for each vertex in the component, so it is O(n^2)
 
 -- |Is the graph connected?
 isConnected :: (Ord t) => Graph t -> Bool
@@ -317,8 +341,8 @@     where components' [] = []
           components' (v:vs) = let c = component g v in c : components' (vs LS.\\ c)
 
--- MORE GRAPHS
 
+-- MORE GRAPHS
 
 -- Generalized Johnson graph, Godsil & Royle p9
 -- Also called generalised Kneser graph, http://en.wikipedia.org/wiki/Kneser_graph
@@ -356,7 +380,7 @@       ++ (map . map) Right [ [i, (i+k) `mod` n] | i <- [0..n-1] ]
 
 petersen2 = gp 5 2
-prism n = gp n 1
+prism' n = gp n 1
 durer = gp 6 2
 mobiusKantor = gp 8 3
 dodecahedron2 = gp 10 2

@@ -1,20 +1,26 @@--- Copyright (c) David Amos, 2009. All rights reserved.
+-- Copyright (c) David Amos, 2009-2014. All rights reserved.
 
+{-# LANGUAGE NoMonomorphismRestriction, TupleSections, DeriveFunctor #-}
+
 module Math.Combinatorics.GraphAuts (isVertexTransitive, isEdgeTransitive,
-                                     isArcTransitive, is2ArcTransitive, is3ArcTransitive, isnArcTransitive,
+                                     isArcTransitive, is2ArcTransitive, is3ArcTransitive, is4ArcTransitive, isnArcTransitive,
                                      isDistanceTransitive,
-                                     graphAuts, incidenceAuts,
+                                     graphAuts, incidenceAuts, graphAuts7, graphAuts8, incidenceAuts2,
+                                     isGraphAut, isIncidenceAut,
                                      graphIsos, incidenceIsos,
                                      isGraphIso, isIncidenceIso) where
 
-import Data.Either (lefts)
+import Data.Either (lefts, rights, partitionEithers)
 import qualified Data.List as L
 import qualified Data.Map as M
 import qualified Data.Set as S
 import Data.Maybe
+import Data.Ord (comparing)
+import qualified Data.Foldable as Foldable
+import qualified Data.Sequence as Seq
 
 import Math.Common.ListSet
-import Math.Core.Utils (combinationsOf, pairs)
+import Math.Core.Utils (combinationsOf, intersectAsc, pairs, picks, (^-))
 import Math.Combinatorics.Graph
 -- import Math.Combinatorics.StronglyRegularGraph
 -- import Math.Combinatorics.Hypergraph -- can't import this, creates circular dependency
@@ -55,7 +61,8 @@     orbitP auts v == v:vs && -- isVertexTransitive g
     orbitP stab n == n:ns
     where auts = graphAuts g
-          stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
+          stab = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts
+          -- stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
           n:ns = nbrs g v
 
 -- execution time of both of the above is dominated by the time to calculate the graph auts, so their performance is similar
@@ -77,7 +84,7 @@ 
 -- note that a graph with triangles can't be 3-arc transitive, etc, because an aut can't map a self-crossing arc to a non-self-crossing arc
 
--- |A graph is n-arc-transitive is its automorphism group is transitive on n-arcs. (An n-arc is an ordered sequence (v0,...,vn) of adjacent vertices, with crossings allowed but not doubling back.)
+-- |A graph is n-arc-transitive if its automorphism group is transitive on n-arcs. (An n-arc is an ordered sequence (v0,...,vn) of adjacent vertices, with crossings allowed but not doubling back.)
 isnArcTransitive :: (Ord t) => Int -> Graph t -> Bool
 isnArcTransitive _ (G [] []) = True
 isnArcTransitive n g@(G (v:vs) es) =
@@ -85,7 +92,8 @@     orbit (->^) a stab == a:as
     -- closure [a] [ ->^ h | h <- stab] == a:as
     where auts = graphAuts g
-          stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
+          stab = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts
+          -- stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
           a:as = findArcs g v n
 
 is2ArcTransitive :: (Ord t) => Graph t -> Bool
@@ -94,6 +102,10 @@ is3ArcTransitive :: (Ord t) => Graph t -> Bool
 is3ArcTransitive g = isnArcTransitive 3 g
 
+-- The incidence graphs of the projective planes PG(2,Fq) are 4-arc-transitive
+is4ArcTransitive :: (Ord t) => Graph t -> Bool
+is4ArcTransitive g = isnArcTransitive 4 g
+
 -- Godsil & Royle 66-7
 -- |A graph is distance transitive if given any two ordered pairs of vertices (u,u') and (v,v') with d(u,u') == d(v,v'),
 -- there is an automorphism of the graph that takes (u,u') to (v,v')
@@ -103,34 +115,30 @@     | isConnected g =
         orbitP auts v == v:vs && -- isVertexTransitive g
         length stabOrbits == diameter g + 1 -- the orbits under the stabiliser of v coincide with the distance partition from v
-    | otherwise = error "isDistanceTransitive: only defined for connected graphs"
+    | otherwise = error "isDistanceTransitive: only implemented for connected graphs"
     where auts = graphAuts g
-          stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
+          stab = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts
+          -- stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
           stabOrbits = let os = orbits stab in os ++ map (:[]) ((v:vs) L.\\ concat os) -- include fixed point orbits
 
 
 -- GRAPH AUTOMORPHISMS
 
--- !! Note, in the literature the following is just called the intersection of two partitions
--- !! Refinement actually refers to the process of refining to an equitable partition
-
--- refine one partition by another
-refine p1 p2 = filter (not . null) $ refine' p1 p2
--- Refinement preserves ordering within cells but not between cells
--- eg the cell [1,2,3,4] could be refined to [2,4],[1,3]
-
--- refine, but leaving null cells in
--- we use this in the graphAuts functions when comparing two refinements to check that they split in the same way
-refine' p1 p2 = concat [ [c1 `intersect` c2 | c2 <- p2] | c1 <- p1]
-
-
+-- |Is the permutation an automorphism of the graph?
+isGraphAut :: Ord t => Graph t -> Permutation t -> Bool
 isGraphAut (G vs es) h = all (`S.member` es') [e -^ h | e <- es]
     where es' = S.fromList es
 -- this works best on sparse graphs, where p(edge) < 1/2
 -- if p(edge) > 1/2, it would be better to test on the complement of the graph
 
-
-
+-- |Is the permutation an automorphism of the incidence structure represented by the graph?
+-- (Note that an incidence graph colours points as Left, blocks as Right, and a permutation
+-- that swaps points and blocks, even if it is an automorphism of the graph, does not represent
+-- an automorphism of the incidence structure. Instead, a point-block crossover is called a duality.)
+isIncidenceAut :: (Ord p, Ord b) => Graph (Either p b) -> Permutation (Either p b) -> Bool
+isIncidenceAut (G vs es) h = all (`S.member` es') [e ->^ h | e <- es]
+    -- using ->^ instead of -^ excludes dualities, since each edge is of the form [Left p, Right b]
+    where es' = S.fromList es
 
 -- Calculate a map consisting of neighbour lists for each vertex in the graph
 -- If a vertex has no neighbours then it is left out of the map
@@ -143,260 +151,482 @@ -- ALTERNATIVE VERSIONS OF GRAPH AUTS
 -- (showing how we got to the final version)
 
--- return all graph automorphisms, using naive depth first search
-graphAuts1 (G vs es) = dfs [] vs vs
-    where dfs xys (x:xs) ys =
-              concat [dfs ((x,y):xys) xs (L.delete y ys) | y <- ys, isCompatible (x,y) xys]
-          dfs xys [] [] = [fromPairs xys]
-          isCompatible (x,y) xys = and [([x',x] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x',y') <- xys]
-          es' = S.fromList es
+data SearchTree a = T Bool a [SearchTree a] deriving (Eq, Ord, Show, Functor)
+-- The boolean indicates whether or not this is a terminal / solution node
 
--- return generators for graph automorphisms
--- (using Lemma 9.1.1 from Seress p203 to prune the search tree)
-graphAuts2 (G vs es) = graphAuts' [] vs
-    where graphAuts' us (v:vs) =
-              let uus = zip us us
-              in concat [take 1 $ dfs ((v,w):uus) vs (v : L.delete w vs) | w <- vs, isCompatible (v,w) uus]
-              ++ graphAuts' (v:us) vs
-              -- stab us == transversal for stab (v:us) ++ stab (v:us)  (generators thereof)
-          graphAuts' _ [] = [] -- we're not interested in finding the identity element
-          dfs xys (x:xs) ys =
-              concat [dfs ((x,y):xys) xs (L.delete y ys) | y <- ys, isCompatible (x,y) xys]
-          dfs xys [] [] = [fromPairs xys]
-          isCompatible (x,y) xys = and [([x',x] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x',y') <- xys]
-          es' = S.fromList es
+leftDepth (T _ _ []) = 1
+leftDepth (T _ _ (t:ts)) = 1 + leftDepth t
 
--- Now using distance partitions
--- Note that because of the use of distance partitions, this is only valid for connected graphs
-graphAuts3 g@(G vs es) = graphAuts' [] [vs] where
-    graphAuts' us ((x:ys):pt) =
-        let px = refine' (ys : pt) (dps M.! x)
-            p y = refine' ((x : L.delete y ys) : pt) (dps M.! y)
-            uus = zip us us
-            p' = L.sort $ filter (not . null) $ px
-        in concat [take 1 $ dfs ((x,y):uus) px (p y) | y <- ys]
-        ++ graphAuts' (x:us) p'
-    graphAuts' us ([]:pt) = graphAuts' us pt
-    graphAuts' _ [] = []
-    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 []
-                    -- 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)
-                                   (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]
+leftWidths (T _ _ []) = []
+leftWidths (T _ _ ts@(t:_)) = length ts : leftWidths t
+
+graphAutsEdgeSearchTree (G vs es) = dfs [] vs vs where
+    dfs xys (x:xs) yys = T False xys [dfs ((x,y):xys) xs ys | (y,ys) <- picks yys, isCompatible xys (x,y)]
+    dfs xys [] [] = T True xys []
+    isCompatible xys (x',y') = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys]
     es' = S.fromList es
 
+graphAuts1 = map fromPairs . terminals . graphAutsEdgeSearchTree
+
+terminals (T False _ ts) = concatMap terminals ts
+terminals (T True xys _) = [xys]
+
+-- Using Lemma 9.1.1 from Seress p203 to prune the search tree
+-- Because auts form a group, it is sufficient to expand only each leftmost branch of the tree in full.
+-- For every other branch, it is sufficient to find a single representative, since the other elements
+-- can then be obtained by multiplication in the group (using the leftmost elements).
+-- In effect, we are finding a transversal generating set.
+-- Note however, that this transversal generating set is relative to whatever base order the tree uses,
+-- so for clarity, the tree should use natural vertex order.
+transversalTerminals (T False _ (t:ts)) = concatMap (take 1 . transversalTerminals) ts ++ transversalTerminals t
+-- transversalTerminals (T False _ (t:ts)) = transversalTerminals t ++ concatMap (take 1 . transversalTerminals) ts
+transversalTerminals (T True xys _) = [xys]
+transversalTerminals _ = []
+
+graphAuts2 = filter (/=1) . map fromPairs . transversalTerminals . graphAutsEdgeSearchTree
+-- init because last is identity
+
 isSingleton [_] = True
 isSingleton _ = False
 
+intersectCells p1 p2 = concat [ [c1 `intersectAsc` c2 | c2 <- p2] | c1 <- p1]
+-- Intersection preserves ordering within cells but not between cells
+-- eg the cell [1,2,3,4] could be refined to [2,4],[1,3]
 
--- 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
-graphAuts4 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)
-        let p' = L.sort $ refine (ys:pt) (dps M.! x)
-        in level us p x ys []
-        ++ graphAuts' (x:us) p'
-    graphAuts' us ([]:pt) = graphAuts' us pt
-    graphAuts' _ [] = []
-    level us p@(ph:pt) x (y: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 _ _ _ [] _ = []
-    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]
+
+graphAutsDistancePartitionSearchTree g@(G vs es) = dfs [] ([vs],[vs]) where
+    dfs xys (srcPart,trgPart)
+        | all isSingleton srcPart =
+             let xys' = zip (concat srcPart) (concat trgPart)
+             in T (isCompatible xys') (xys++xys') []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x:xs):srcCells = srcPart
+                          yys   :trgCells = trgPart
+                          srcPart' = intersectCells (xs : srcCells) (dps M.! x)
+                      in T False xys -- the L.sort in the following line is so that we traverse vertices in natural order
+                         [dfs ((x,y):xys) ((unzip . L.sort) (zip (filter (not . null) srcPart') (filter (not . null) trgPart')))
+                         | (y,ys) <- picks yys,
+                           let trgPart' = intersectCells (ys : trgCells) (dps M.! y),
+                           map length srcPart' == map length trgPart']
     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
+    dps = M.fromAscList [(v, distancePartitionS vs es' v) | v <- vs]
 
--- 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
+graphAuts3 = filter (/=1) . map fromPairs . transversalTerminals . graphAutsDistancePartitionSearchTree
 
+-- Whereas transversalTerminals produced a transversal generating set, here we produce a strong generating set.
+-- In particular, if we have already found (3 4), and then we find (1 2 3),
+-- then there is no need to look for (1 3 ...) or (1 4 ...), since it is clear that such elements exist
+-- as products of those we have already found.
+strongTerminals = strongTerminals' [] where
+    strongTerminals' gs (T False xys ts) =
+        case listToMaybe $ reverse $ filter (\(x,y) -> x /= y) xys of -- the first vertex that isn't fixed
+        Nothing -> L.foldl' (\hs t -> strongTerminals' hs t) gs ts
+        Just (x,y) -> if y `elem` (x .^^ gs)
+                      then gs
+                      -- Since we're not on the leftmost spine, we can stop as soon as we find one new element
+                      else find1New gs ts
+                      -- else L.foldl' (\hs t -> if hs /= gs then hs else strongTerminals' hs t) gs ts
+    strongTerminals' gs (T True xys []) = fromPairs xys : gs
+    find1New gs (t:ts) = let hs = strongTerminals' gs t
+                         in if take 1 gs /= take 1 hs -- we know a new element would be placed at the front
+                            then hs
+                            else find1New gs ts
+    find1New gs [] = gs
 
+-- |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 = filter (/=1) . strongTerminals . graphAutsDistancePartitionSearchTree
 
--- an example for equitable partitions
--- this is a graph whose distance partition (from any vertex) can be refined to an equitable partition
-eqgraph = G vs es where
-    vs = [1..14]
-    es = L.sort $ [[1,14],[2,13]] ++ [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, v1+1 == v2 || v1+3 == v2 && even v2]
 
--- refine a partition to give an equitable partition
-toEquitable g cells = L.sort $ toEquitable' [] cells where
-    toEquitable' ls (r:rs) =
-        let (lls,lrs) = L.partition isSingleton $ map (splitNumNbrs r) ls
-            -- so the lrs split, and the lls didn't
-            rs' = concatMap (splitNumNbrs r) rs
-        in if isSingleton r -- then we know it won't split further, so can remove it from further processing
-           then r : toEquitable' (concat lls) (concat lrs ++ rs')
-           else toEquitable' (r : concat lls) (concat lrs ++ rs')
-    toEquitable' ls [] = ls
-    splitNumNbrs t c = map (map snd) $ L.groupBy (\x y -> fst x == fst y) $ L.sort
-                    [ (length ((nbrs_g M.! v) `intersect` t), v) | v <- c]
-    nbrs_g = M.fromList [(v, nbrs g v) | v <- vertices g]
+-- Using colourings (M.Map vertex colour, M.Map colour [vertex]), in place of partitions ([[vertex]])
+-- This turns out to be slower than using partitions.
+-- Updating the colour partition incrementally seems to be much less efficient than just recalculating it each time
+-- (Recalculating each time is O(n), incrementally updating is O(n^2)?)
+graphAutsDistanceColouringSearchTree g@(G vs es) = dfs [] unitCol unitCol where
+    unitCol = (M.fromList $ map (,[]) vs, M.singleton [] vs) -- "unit colouring"
+    dfs xys srcColouring@(srcVmap,srcCmap) trgColouring@(trgVmap,trgCmap)
+        -- ( | M.map length srcCmap /= M.map length trgCmap = T False xys [] )
+        | all isSingleton (M.elems srcCmap) = -- discrete colouring
+             let xys' = zip (concat $ M.elems srcCmap) (concat $ M.elems trgCmap)
+             in T (isCompatible xys') (reverse xys'++xys) []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x,c) = M.findMin srcVmap
+                          (xVmap,xCmap) = dcs M.! x
+                          ys = trgCmap M.! c
+                          srcVmap' = M.delete x (intersectColouring srcVmap xVmap)
+                          srcCmap' = colourPartition srcVmap'
+                          -- srcCmap' = M.fromAscList [(c1++c2, cell) | (c1,srcCell) <- M.assocs srcCmap, (c2,xCell) <- M.assocs xCmap,
+                          --                                            let cell = L.delete x (intersectAsc srcCell xCell),
+                          --                                            (not . null) cell]
+                      in T False xys
+                         [dfs ((x,y):xys) (srcVmap',srcCmap') (trgVmap',trgCmap')
+                         | y <- ys,
+                           let (yVmap,yCmap) = dcs M.! y,
+                           let trgVmap' = M.delete y (intersectColouring trgVmap yVmap),
+                           let trgCmap' = colourPartition trgVmap',
+                           -- let trgCmap' = M.fromAscList [(c1++c2, cell) | (c1,trgCell) <- M.assocs trgCmap, (c2,yCell) <- M.assocs yCmap,
+                           --                                                let cell = L.delete y (intersectAsc trgCell yCell),
+                           --                                                (not . null) cell],
+                           M.map length srcCmap' == M.map length trgCmap' ]
+    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+    es' = S.fromList es
+    dcs = M.fromAscList [(v, distanceColouring v) | v <- vs]
+    distanceColouring u = let dp = distancePartitionS vs es' u
+                              vmap = M.fromList [(v,[c]) | (cell,c) <- zip dp [0..], v <- cell]
+                              cmap = M.fromList $ zip (map (:[]) [0..]) dp
+                          in (vmap, cmap)
 
+{-
+-- If we are going to recalculate the colour partition each time anyway,
+-- then we don't need to carry it around, and can simplify the code
+graphAutsDistanceColouringSearchTree g@(G vs es) = dfs [] initCol initCol where
+    initCol = M.fromList $ map (,[]) vs
+    dfs xys srcCol trgCol
+        | M.map length srcPart /= M.map length trgPart = T False xys []
+        | all isSingleton (M.elems srcPart) =
+             let xys' = zip (concat $ M.elems srcPart) (concat $ M.elems trgPart)
+             in T (isCompatible xys') (reverse xys'++xys) []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x,c) = M.findMin srcCol
+                          ys = trgPart M.! c
+                          srcCol' = M.delete x $ intersectColouring srcCol (dcs M.! x)
+                      in T False xys
+                         [dfs ((x,y):xys) srcCol' trgCol'
+                         | y <- ys,
+                           let trgCol' = M.delete y (intersectColouring trgCol (dcs M.! y))]
+        where srcPart = colourPartition srcCol
+              trgPart = colourPartition trgCol
+    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+    es' = S.fromList es
+    dcs = M.fromAscList [(v, distanceColouring v) | v <- vs]
+    distanceColouring u = M.fromList [(v,[c]) | (cell,c) <- zip (distancePartitionS vs es' u) [0..], v <- cell]
+-}
+distanceColouring (G vs es) u = M.fromList [(v,[c]) | (cell,c) <- zip (distancePartitionS vs es' u) [0..], v <- cell]
+    where es' = S.fromList es
 
--- try to refine two partitions in parallel, failing if they become mismatched
-toEquitable2 nbrs_g psrc ptrg = unzip $ L.sort $ toEquitable' [] (zip psrc ptrg) where
-    toEquitable' ls (r:rs) =
-        let ls' = map (splitNumNbrs nbrs_g r) ls
-            (lls,lrs) = L.partition isSingleton $ map fromJust ls'
-            rs' = map (splitNumNbrs nbrs_g r) rs
-        in if any isNothing ls' || any isNothing rs'
-           then []
-           else
-               {- if (isSingleton . fst) r
-               then r : toEquitable' (concat lls) (concat lrs ++ concatMap fromJust rs')
-               else -} toEquitable' (r : concat lls) (concat lrs ++ concatMap fromJust rs')
-    toEquitable' ls [] = ls
+intersectColouring c1 c2 = M.intersectionWith (++) c1 c2
 
-splitNumNbrs nbrs_g (t_src,t_trg) (c_src,c_trg) =
-    let src_split = L.groupBy (\x y -> fst x == fst y) $ L.sort
-                    [ (length ((nbrs_g M.! v) `intersect` t_src), v) | v <- c_src]
-        trg_split = L.groupBy (\x y -> fst x == fst y) $ L.sort
-                    [ (length ((nbrs_g M.! v) `intersect` t_trg), v) | v <- c_trg]
-    in if map length src_split == map length trg_split
-       && map (fst . head) src_split == map (fst . head) trg_split
-       then Just $ zip (map (map snd) src_split) (map (map snd) trg_split)
-       else Nothing
-       -- else error (show (src_split, trg_split)) -- for debugging
+colourPartition c = L.foldr (\(k,v) m -> M.insertWith (++) v [k] m) M.empty (M.assocs c)
 
--- Now, every time we intersect two partitions, refine to an equitable partition
 
--- |Given a graph g, @graphAuts g@ returns generators for the automorphism group of g.
--- If g is connected, then the generators will be a strong generating set.
-graphAuts :: (Ord a) => Graph a -> [Permutation a]
-graphAuts g = autsWithinComponents ++ isosBetweenComponents
-    where cs = map (inducedSubgraph g) (components g)
-          -- autsWithinComponents = concatMap graphAutsCon cs
-          autsWithinComponents = concatMap graphAuts4 cs
-          isosBetweenComponents = map swapFromIso $ concat [take 1 (graphIsos ci cj) | (ci,cj) <- pairs cs]
-          swapFromIso xys = fromPairs (xys ++ map swap xys)
-          swap (x,y) = (y,x)
--- Using graphAuts4 instead of graphAutsCon as latter appears to have a bug, eg
--- > graphAuts4 $ G [1..3] [[1,2],[2,3]]
--- [[[1,3]]]
--- > graphAutsCon $ G [1..3] [[1,2],[2,3]]
--- []
 
--- Automorphisms of a connected graph
-graphAutsCon g@(G vs es)
-    | isConnected g = graphAuts' [] (toEquitable g $ valencyPartition g)
-    | otherwise = error "graphAutsCon: graph is not connected"
-    where graphAuts' us p@((x:ys):pt) =
-              let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
-              in level us p x ys []
-              ++ graphAuts' (x:us) p'
-          graphAuts' us ([]:pt) = graphAuts' us pt
-          graphAuts' _ [] = []
-          level us p@(ph:pt) x (y: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 dfsEquitable (dps,es',nbrs_g) ((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 _ _ _ [] _ = []
-          dps = M.fromList [(v, distancePartition g v) | v <- vs]
-          es' = S.fromList es
-          nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
+-- Based on McKay’s Canonical Graph Labeling Algorithm, by Stephen G. Hartke and A. J. Radcliffe
+-- (http://www.math.unl.edu/~aradcliffe1/Papers/Canonical.pdf)
 
-dfsEquitable (dps,es',nbrs_g) xys p1 p2 = dfs xys p1 p2 where
-    dfs xys p1 p2
-        | map length p1 /= map length p2 = []
-        | otherwise =
-            let p1' = filter (not . null) p1
-                p2' = filter (not . null) p2
-                (p1e,p2e) = toEquitable2 nbrs_g p1' p2'
-            in if null p1e
-               then []
-               else
-                   if all isSingleton p1e
-                   then let xys' = xys ++ zip (concat p1e) (concat p2e)
-                        in if isCompatible xys' then [fromPairs' xys'] else []
-                   else let (x:xs):p1'' = p1e
-                            ys:p2'' = p2e
-                        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']
+equitableRefinement g@(G vs es) p = equitableRefinement' (S.fromList es) p
 
+equitableRefinement' edgeset partition = go partition where
+    go cells = let splits = L.zip (L.inits cells) (L.tails cells)
+                   shatterPairs = [(L.zip ci counts,ls,rs) | (ls,ci:rs) <- splits, cj <- cells,
+                                                             let counts = map (nbrCount cj) ci, isShatter counts]
+               in case shatterPairs of -- by construction, the lexicographic least (i,j) comes first
+                  [] -> cells
+                  (vcs,ls,rs):_ -> let fragments = shatter vcs 
+                                   in go (ls ++ fragments ++ rs)
+    isShatter (c:cs) = any (/= c) cs
+    shatter vcs = map (map fst) $ L.groupBy (\x y -> snd x == snd y) $ L.sortBy (comparing snd) $ vcs
+    -- Memoizing here results in about 10% speed improvement. Not worth it for loss of generality (ie requiring HasTrie instances)
+    -- nbrCount = memo2 nbrCount'
+    -- How many neighbours in cell does vertex have
+    nbrCount cell vertex = length (filter (isEdge vertex) cell)
+    isEdge u v = L.sort [u,v] `S.member` edgeset
 
+equitablePartitionSearchTree g@(G vs es) p = dfs [] p where
+    dfs bs p = let p' = equitableRefinement' es' p in
+               if all isSingleton p'
+               then T True (p',bs) []
+               else T False (p',bs) [dfs (b:bs) p'' | (b,p'') <- splits [] p']
+    -- For now, we just split the first non-singleton cell we find
+    splits ls (r:rs) | isSingleton r = splits (r:ls) rs
+                     | otherwise = let ls' = reverse ls in [(x, ls' ++ [x]:xs:rs) | (x,xs) <- picks r]
+    es' = S.fromList es
+
+
+{-
+-- Using Data.Sequence instead of list for the partitions
+-- Makes no difference to speed (in fact slightly slower)
+equitableRefinementSeq' edgeset partition = go partition where
+    go cells = let splits = Seq.zip (Seq.inits cells) (Seq.tails cells)
+                   shatterPairs = [(L.zip ci counts,ls,rs') | (ls,rs) <- Foldable.toList splits, (not . Seq.null) rs, let ci Seq.:< rs' = Seq.viewl rs,
+                                                              cj <- Foldable.toList cells,
+                                                              let counts = map (nbrCount cj) ci, isShatter counts]
+               in case shatterPairs of -- by construction, the lexicographic least (i,j) comes first
+                  [] -> cells
+                  (vcs,ls,rs):_ -> let fragments = Seq.fromList (shatter vcs) 
+                                   in go (ls Seq.>< fragments Seq.>< rs)
+    isShatter (c:cs) = any (/= c) cs
+    shatter vcs = map (map fst) $ L.groupBy (\x y -> snd x == snd y) $ L.sortBy (comparing snd) $ vcs
+    -- How many neighbours in cell does vertex have
+    nbrCount cell vertex = length (filter (isEdge vertex) cell)
+    isEdge u v = L.sort [u,v] `S.member` edgeset
+
+equitablePartitionSeqSearchTree g@(G vs es) p = dfs [] (Seq.fromList p) where
+    dfs bs p = let p' = equitableRefinementSeq' es' p in
+               if Foldable.all isSingleton p'
+               then T True (Foldable.toList p',bs) []
+               else T False (Foldable.toList p',bs) [dfs (b:bs) p'' | (b,p'') <- splits p']
+    -- For now, we just split the first non-singleton cell we find
+    splits cells = case Seq.findIndexL (not . isSingleton) cells of
+                   Just i -> let (ls,rs) = Seq.splitAt i cells
+                                 r Seq.:< rs' = Seq.viewl rs
+                             in [(x, ls Seq.>< ([x] Seq.<| xs Seq.<| rs')) | (x,xs) <- picks r]
+                   Nothing -> error "Not possible, as we know there are non-singleton cells"
+    es' = S.fromList es
+-}
+
+-- In this version, whenever we have an equitable partition, we separate out all the singleton cells and put them to one side.
+-- (Since the partition is equitable, singleton cells have already done any work they are going to do in shattering other cells,
+-- so they will no longer play any part.)
+-- This seems to result in about 20% speedup.
+equitablePartitionSearchTree2 g@(G vs es) p = dfs [] ([],p) where
+    dfs bs (ss,cs) = let (ss',cs') = L.partition isSingleton $ equitableRefinement' es' cs
+                         ss'' = ss++ss'
+                     in case cs' of
+                        [] -> T True (ss'',bs) []
+                        -- We just split the first non-singleton cell
+                        -- c:cs'' -> T False (ss''++cs',bs) [dfs (x:bs) (ss'',[x]:xs:cs'') | (x,xs) <- picks c]
+                        c:cs'' -> T False (cs'++ss'',bs) [dfs (x:bs) (ss'',[x]:xs:cs'') | (x,xs) <- picks c]
+    es' = S.fromList es
+-- TODO: On the first level, we can use a stronger partitioning function (eg distance partitions, + see nauty manual, vertex invariants)
+
+equitableDistancePartitionSearchTree g@(G vs es) p = dfs [] p where
+    dfs bs p = let p' = equitableRefinement' es' p in
+               if all isSingleton p'
+               then T True (p',bs) []
+               else T False (p',bs) [dfs (b:bs) p'' | (b,p'') <- splits [] p']
+    -- For now, we just split the first non-singleton cell we find
+    splits ls (r:rs) | isSingleton r = splits (r:ls) rs
+                     | otherwise = [(x, p'') | let ls' = reverse ls,
+                                               (x,xs) <- picks r,
+                                               let p' = ls' ++ [x]:xs:rs,
+                                               let p'' = filter (not . null) (intersectCells p' (dps M.! x))]
+    es' = S.fromList es
+    dps = M.fromAscList [(v, distancePartitionS vs es' v) | v <- vs]
+
+
+{-
+-- This is just fmap (\(p,bs) -> (p,bs,trace p)) t
+equitablePartitionTracedSearchTree g@(G vs es) trace p = dfs [] p where
+    dfs bs p = let p' = equitableRefinement' es' p
+               in if all isSingleton p'
+                  then T True (p',bs,trace p') []
+                  else T False (p',bs,trace p') [dfs (b:bs) p'' | (b,p'') <- splits [] p']
+    -- For now, we just split the first non-singleton cell we find
+    splits ls (r:rs) | isSingleton r = splits (r:ls) rs
+                     | otherwise = let ls' = reverse ls in [(x, ls' ++ [x]:xs:rs) | (x,xs) <- picks r]
+    es' = S.fromList es
+-}
+
+-- Intended as a node invariant
+trace1 p = map (\xs@(x:_) -> (x, length xs)) $ L.group $ L.sort $ map length p
+
+equitablePartitionGraphSearchTree g@(G vs es) = equitablePartitionSearchTree g unitPartition
+    where unitPartition = [vs]
+
+-- The incidence graph has vertices that are coloured left (points) or right (blocks).
+-- We are not interested in dualities (automorphisms that swap points and blocks), so we look for colour-preserving automorphisms
+equitablePartitionIncidenceSearchTree g@(G vs es) = equitablePartitionSearchTree g lrPartition
+    where (lefts, rights) = partitionEithers vs
+          lrPartition = [map Left lefts, map Right rights]
+
+leftLeaf (T False _ (t:ts)) = leftLeaf t
+leftLeaf (T True (p,bs) []) = (concat p, reverse bs)
+{-
+leftSpine (T False x (t:ts)) = x : leftSpine t
+leftSpine (T True x []) = [x]
+-}
+allLeaves (T False _ ts) = concatMap allLeaves ts
+allLeaves (T True (p,bs) []) = [(concat p, reverse bs)]
+
+{-
+partitionTransversals tree = [fromPairs (zip canonical partition) | partition <- findTransversals tree] where
+    (_,canonical) = leftLeaf tree
+    findTransversals (T False _ (t:ts)) = concatMap (take 1 . findTransversals) ts ++ findTransversals t
+    findTransversals (T True (_,partition) []) = [concat partition]
+
+graphAuts5 = partitionTransversals . equitablePartitionGraphSearchTree
+-}
+-- NOT WORKING
+partitionBSGS0 g@(G vs es) t = (bs, findLevels t) where
+    (p1,bs) = leftLeaf t
+    g1 = fromPairs $ zip p1 vs
+    g1' = g1^-1
+    es1 = S.fromList $ edges $ fmap (.^ g1) g -- the edges of the isomorph corresponding to p1. (S.fromList makes it unnecessary to call nf.)
+    findLevels (T True (partition,_) []) = []
+    findLevels (T False (partition,_) (t:ts)) =
+        let hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+        in findLevel v hs (zip vs ts)
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut
+    -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.)
+    find1New (T True (partition,_) []) = let h = fromPairs $ zip (concat partition) vs
+                                             g' = fmap (.^ h) g
+                                         in if all (`S.member` es1) (edges g') then [h*g1'] else []
+    -- isAut h = all (`S.member` es') [e -^ h | e <- es]
+    -- es' = S.fromList es
+
+-- Given a partition search tree, return a base and strong generating set for graph automorphism group.
+partitionBSGS g@(G vs es) t = (bs, findLevels t) where
+    (canonical,bs) = leftLeaf t
+    findLevels (T True (partition,_) []) = []
+    findLevels (T False (partition,_) (t:ts)) =
+        let hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+        in findLevel v hs (zip vs ts)
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs -- TODO: Memoize this orbit
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- Some leaf nodes correspond to different isomorphs of the graph, and hence don't yield automorphisms
+    find1New (T True (partition,_) []) = let h = fromPairs $ zip canonical (concat partition)
+                                         in filter isAut [h]
+    isAut h = all (`S.member` es') [e -^ h | e <- es]
+    es' = S.fromList es
+-- The tree for g1 has leaf nodes of two different isomorphs, as does the tree for incidenceGraphPG 2 f8
+
+-- Returns auts as Right, different isomorphs as Left
+-- (Must be used with the tree which doesn't put singletons to end)
+partitionBSGS3 g@(G vs es) t = (bs, findLevels t) where
+    (p1,bs) = leftLeaf t
+    findLevels (T True (partition,_) []) = []
+    findLevels (T False (partition,_) (t:ts)) =
+        let hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+        in findLevel v hs (zip vs ts)
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ rights hs
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut
+    -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.)
+    find1New (T True (partition,_) []) = let h = fromPairs $ zip p1 (concat partition)
+                                         in if isAut h then [Right h] else [Left h]
+    isAut h = all (`S.member` es') [e -^ h | e <- es]
+    es' = S.fromList es
+-- TODO: I think we are only justified in doing find1New (ie only finding 1) if we *do* find an aut.
+-- If we don't, we should potentially keep looking in that subtree
+-- (See section 6 of paper. If we find isomorphic leaves, then the two subtrees of their common parent are isomorphic,
+-- so no need to continue searching the second.)
+
+
+-- This is using a node invariant to do more pruning.
+-- However, seems to be much slower on very regular graphs (where perhaps there is no pruning to be done)
+-- (This suggests that perhaps using fmap is not good - perhaps a space leak?)
+-- (Or perhaps it's just that calculating and comparing the node invariants is expensive)
+-- TODO: Perhaps use something simpler, like just the number of cells in the partition
+partitionBSGS2 g@(G vs es) t = (bs, findLevels t') where
+    t' = fmap (\(p,bs) -> (p,bs,trace1 p)) t
+    trace1 = length -- the number of cells in the partition
+    (canonical,bs) = leftLeaf t
+    findLevels (T True (partition,_,_) []) = []
+    findLevels (T False (partition,_,_) (t:ts)) =
+        let (T _ (_,_,trace) _) = t
+            hs = findLevels t
+            -- TODO: It might be better to use the b that is added in t to find the cell that splits
+            cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split
+            vts = filter (\(_,T _ (_,_,trace') _) -> trace == trace') $ zip vs ts
+        in findLevel v hs vts
+    findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs
+                                   then findLevel v hs vts
+                                   else let h = find1New t' in findLevel v (h++hs) vts
+    findLevel _ hs [] = hs
+    find1New (T False _ ts) = take 1 $ concatMap find1New ts
+    -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut
+    -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.)
+    -- (The graph g1, below, shows a simple example where this will happen.)
+    find1New (T True (partition,_,_) []) = let h = fromPairs $ zip canonical (concat partition)
+                                           in filter isAut [h]
+    isAut h = all (`S.member` es') [e -^ h | e <- es]
+    es' = S.fromList es
+
+
+graphAuts7 g = (partitionBSGS g) (equitablePartitionGraphSearchTree g)
+
+-- This is faster on kneser graphs, but slower on incidenceGraphPG
+graphAuts8 g = (partitionBSGS g) (equitableDistancePartitionSearchTree g [vertices g])
+
+-- This is a graph where the node invariant should cause pruning.
+-- The initial equitable partition will be [[1..8],[9,10]], because it can do no better than distinguish by degree
+-- However, vertices 1..4 and vertices 5..8 are in fact different (there is no aut that takes one set to the other),
+-- so the subtrees starting 1..4 have a different invariant to those starting 5..8
+g1 = G [1..10] [[1,2],[1,3],[1,9],[2,4],[2,10],[3,4],[3,9],[4,10],[5,6],[5,8],[5,9],[6,7],[6,10],[7,8],[7,9],[8,10]]
+
+g1' = nf $ fmap (\x -> if x <= 4 then x+4 else if x <= 8 then x-4 else x) g1
+-- G [1..10] [[1,2],[1,4],[1,9],[2,3],[2,10],[3,4],[3,9],[4,10],[5,6],[5,7],[5,9],[6,8],[6,10],[7,8],[7,9],[8,10]]
+
+g2 = G [1..12] [[1,2],[1,4],[1,11],[2,3],[2,12],[3,4],[3,11],[4,12],[5,6],[5,8],[5,11],[6,9],[6,12],[7,8],[7,10],[7,11],[8,12],[9,10],[9,11],[10,12]]
+
+-- NOT WORKING: This fails to find the isomorphism between g1 and g1' above.
+-- Instead of using left leaf, we need to find the canonical isomorph, as described in the paper.
+-- (In a graph where not all leaves lead to automorphisms, we might happen to end up with non-isomorphic left leaves)
+maybeGraphIso g1 g2 = let (vs1,_) = (leftLeaf . equitablePartitionGraphSearchTree) g1
+                          (vs2,_) = (leftLeaf . equitablePartitionGraphSearchTree) g2
+                          f = M.fromList (zip vs1 vs2)
+                      in if length vs1 == length vs2 && (nf . fmap (f M.!)) g1 == g2 then Just f else Nothing
+
+
 -- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH
 
--- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts
+-- This code is nearly identical to the corresponding graphAuts code, with two exceptions:
+-- 1. We start by partitioning into lefts and rights.
+-- This avoids left-right crossover auts, which while they are auts of the graph,
+-- are not auts of the incidence structure
+-- 2. When labelling the nodes, we filter out Right blocks, and unLeft the Left points
+incidenceAutsDistancePartitionSearchTree g@(G vs es) = dfs [] (lrPart, lrPart) where
+    dfs xys (srcPart,trgPart)
+        | all isSingleton srcPart =
+             let xys' = zip (concat srcPart) (concat trgPart)
+             in T (isCompatible xys') (unLeft $ xys++xys') []
+             -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible.
+             -- However, we do need to check that the xys' are edge-compatible with each other.
+        | otherwise = let (x:xs):srcCells = srcPart
+                          yys   :trgCells = trgPart
+                          srcPart' = intersectCells (xs : srcCells) (dps M.! x)
+                      in T False (unLeft xys) -- the L.sort in the following line is so that we traverse vertices in natural order
+                         [dfs ((x,y):xys) ((unzip . L.sort) (zip (filter (not . null) srcPart') (filter (not . null) trgPart')))
+                         | (y,ys) <- picks yys,
+                           let trgPart' = intersectCells (ys : trgCells) (dps M.! y),
+                           map length srcPart' == map length trgPart']
+    isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
+    (lefts, rights) = partitionEithers vs
+    lrPart = [map Left lefts, map Right rights] -- Partition the vertices into left and right, to exclude crossover auts
+    unLeft xys = [(x,y) | (Left x, Left y) <- xys] -- also filters out Rights
+    es' = S.fromList es
+    dps = M.fromList [(v, distancePartitionS vs es' v) | v <- vs]
 
 -- |Given the incidence graph of an incidence structure between points and blocks
 -- (for example, a set system),
--- @incidenceAuts g@ returns generators for the automorphism group of the incidence structure.
+-- @incidenceAuts g@ returns a strong generating set for the automorphism group of the incidence structure.
 -- The generators are represented as permutations of the points.
 -- The incidence graph should be represented with the points on the left and the blocks on the right.
--- If the incidence graph is connected, then the generators will be a strong generating set.
 incidenceAuts :: (Ord p, Ord b) => Graph (Either p b) -> [Permutation p]
-incidenceAuts g = autsWithinComponents ++ isosBetweenComponents
-    where cs = map (inducedSubgraph g) (components g)
-          autsWithinComponents = concatMap incidenceAutsCon cs
-          isosBetweenComponents = map swapFromIso $ concat [take 1 (incidenceIsos ci cj) | (ci,cj) <- pairs cs]
-          swapFromIso xys = fromPairs (xys ++ map swap xys)
-          swap (x,y) = (y,x)
+incidenceAuts = filter (/= p []) . strongTerminals . incidenceAutsDistancePartitionSearchTree
 
-incidenceAutsCon g@(G vs es) 
-    | isConnected g = map points (incidenceAuts' [] [vs])
-    | otherwise = error "incidenceAutsCon: graph is not connected"
-    where points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
-          incidenceAuts' us p@((x@(Left _):ys):pt) =
-              -- let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
-              let p' = L.sort $ 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 dfsEquitable (dps,es',nbrs_g) ((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
-          dps = M.fromList [(v, distancePartition g v) | v <- vs]
-          es' = S.fromList es
-          nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]
 
+-- TODO: Filter out rights, map unLeft - to bs and gs
+incidenceAuts2 g = (partitionBSGS g) (equitablePartitionIncidenceSearchTree g)
+    where unLeft (Left x) = x
+          -- map (\g -> fromPairs . map (\(Left x, Left y) -> (x,y)) . filter (\(x,y) -> isLeft x) . toPairs) gs
 
+
 -- GRAPH ISOMORPHISMS
 
 -- !! not yet using equitable partitions, so could probably be more efficient
 
-
 -- graphIsos :: (Ord a, Ord b) => Graph a -> Graph b -> [[(a,b)]]
 graphIsos g1 g2
     | length cs1 /= length cs2 = []
@@ -404,9 +634,8 @@     where cs1 = map (inducedSubgraph g1) (components g1)
           cs2 = map (inducedSubgraph g2) (components g2)
           graphIsos' (ci:cis) cjs =
-              [iso ++ iso' | cj <- cjs,
+              [iso ++ iso' | (cj,cjs') <- picks cjs,
                              iso <- graphIsosCon ci cj,
-                             let cjs' = L.delete cj cjs,
                              iso' <- graphIsos' cis cjs']
           graphIsos' [] [] = [[]]
 
@@ -428,12 +657,16 @@                      else let (x:xs):p1'' = p1'
                               ys:p2'' = p2'
                           in concat [dfs ((x,y):xys)
-                                         (refine' (xs : p1'') (dps1 M.! x))
-                                         (refine' ((L.delete y ys):p2'') (dps2 M.! y))
-                                         | y <- ys]
+                                         (intersectCells (xs : p1'') (dps1 M.! x))
+                                         (intersectCells (ys': p2'') (dps2 M.! y))
+                                         | (y,ys') <- picks ys]
           isCompatible xys = and [([x,x'] `S.member` es1) == (L.sort [y,y'] `S.member` es2) | (x,y) <- xys, (x',y') <- xys, x < x']
-          dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
-          dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
+          dps1 = M.fromAscList [(v, distancePartitionS vs1 es1 v) | v <- vs1]
+          dps2 = M.fromAscList [(v, distancePartitionS vs2 es2 v) | v <- vs2]
+          -- dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
+          -- dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
+          vs1 = vertices g1
+          vs2 = vertices g2
           es1 = S.fromList $ edges g1
           es2 = S.fromList $ edges g2
 
@@ -445,10 +678,7 @@ -- !! then the cost of calculating distancePartitions may not be warranted
 -- !! (see Math.Combinatorics.Poset: orderIsos01 versus orderIsos)
 
--- !! deprecate
-isIso g1 g2 = (not . null) (graphIsos g1 g2)
 
-
 -- the following differs from graphIsos in only two ways
 -- we avoid Left, Right crossover isos, by insisting that a Left is taken to a Left (first two lines)
 -- we return only the action on the Lefts, and unLeft it
@@ -461,9 +691,8 @@     where cs1 = map (inducedSubgraph g1) (filter (not . null . lefts) $ components g1)
           cs2 = map (inducedSubgraph g2) (filter (not . null . lefts) $ components g2)
           incidenceIsos' (ci:cis) cjs =
-              [iso ++ iso' | cj <- cjs,
+              [iso ++ iso' | (cj,cjs') <- picks cjs,
                              iso <- incidenceIsosCon ci cj,
-                             let cjs' = L.delete cj cjs,
                              iso' <- incidenceIsos' cis cjs']
           incidenceIsos' [] [] = [[]]
 
@@ -484,9 +713,9 @@                      else let (x:xs):p1'' = p1'
                               ys:p2'' = p2'
                           in concat [dfs ((x,y):xys)
-                                         (refine' (xs : p1'') (dps1 M.! x))
-                                         (refine' ((L.delete y ys):p2'') (dps2 M.! y))
-                                         | y <- ys]
+                                         (intersectCells (xs : p1'') (dps1 M.! x))
+                                         (intersectCells (ys': p2'') (dps2 M.! y))
+                                         | (y,ys') <- picks ys]
           isCompatible xys = and [([x,x'] `S.member` es1) == (L.sort [y,y'] `S.member` es2) | (x,y) <- xys, (x',y') <- xys, x < x']
           dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
           dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
@@ -498,63 +727,5 @@      Graph (Either p1 b1) -> Graph (Either p2 b2) -> Bool
 isIncidenceIso g1 g2 = (not . null) (incidenceIsos g1 g2)
 
-{-
-removeGens x gs = removeGens' [] gs where
-    baseOrbit = x .^^ gs
-    removeGens' ls (r:rs) =
-        if x .^^ (ls++rs) == baseOrbit
-        then removeGens' ls rs
-        else removeGens' (r:ls) rs
-    removeGens' ls [] = reverse ls
--- !! reverse is probably pointless
-
-
--- !! DON'T THINK THIS IS WORKING PROPERLY
--- eg graphAutsSGSNew $ toGraph ([1..7],[[1,3],[2,3],[3,4],[4,5],[4,6],[4,7]])
--- returns [[[1,2]],[[5,6]],[[5,7,6]],[[6,7]]]
--- whereas [[6,7]] was a Schreier generator, so shouldn't have been listed
-
--- Using Schreier generators to seed the next level
--- At the moment this is slower than the above
--- (This could be modified to allow us to start the search with a known subgroup)
-graphAutsNew g@(G vs es) = graphAuts' [] [] [vs] where
-    graphAuts' us hs p@((x:ys):pt) =
-        let ys' = ys L.\\ (x .^^ hs) -- don't need to consider points which can already be reached from Schreier generators
-            hs' = level us p x ys' []
-            p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
-            reps = cosetRepsGx (hs'++hs) x
-            schreierGens = removeGens x $ schreierGeneratorsGx (x,reps) (hs'++hs)
-        in hs' ++ graphAuts' (x:us) schreierGens p'
-    graphAuts' us hs ([]:pt) = graphAuts' us hs pt
-    graphAuts' _ _ [] = []
-    level us p@(ph:pt) x (y: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 if map length px /= map length py
-           then level us p x ys hs
-           else case dfs ((x,y):uus) (filter (not . null) px) (filter (not . null) py) of
-                []  -> level us p x ys hs
-                h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
-                -- if h1 = (1 2)(3 4), and h2 = (1 3 2), then we can remove 4 too
-    level _ _ _ [] _ = []
-    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
--}
 
 

@@ -1,10 +1,12 @@--- Copyright (c) 2011, David Amos. All rights reserved.+-- Copyright (c) 2011-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction #-}   module Math.Combinatorics.IncidenceAlgebra where +import Prelude hiding ( (*>) )+ import Math.Core.Utils  import Math.Combinatorics.Digraph@@ -103,11 +105,11 @@     -- |Note that we are not able to give a generic definition of unit for the incidence algebra,     -- because it depends on which poset we are working in,     -- and that information is encoded at the value level rather than the type level. See unitIA.-    unit 0 = zero -- so that sum works+    unit 0 = zerov -- so that sum works     -- unit x = x *> sumv [return (Iv (a,a)) | a <- poset] -- the delta function     -- but we can't know from the types alone which poset we are working in     mult = linear mult'-        where mult' (Iv poset (a,b), Iv _ (c,d)) = if b == c then return (Iv poset (a,d)) else zero+        where mult' (Iv poset (a,b), Iv _ (c,d)) = if b == c then return (Iv poset (a,d)) else zerov  -- So multiplication in the incidence algebra is about composition of intervals 

@@ -1,4 +1,4 @@--- Copyright (c) 2011, David Amos. All rights reserved.+-- Copyright (c) 2011-2015, David Amos. All rights reserved.  {-# LANGUAGE NoMonomorphismRestriction, DeriveFunctor #-} @@ -6,6 +6,8 @@ module Math.Combinatorics.Matroid where  -- Source: Oxley, Matroid Theory (second edition)++import Prelude hiding ( (*>) )  import Math.Core.Utils import Math.Core.Field hiding (f7)

@@ -7,7 +7,7 @@  import Math.Common.ListSet as LS -- set operations on strictly ascending lists import Math.Core.Utils -- for set/multiset operations on ordered lists-import Math.Algebra.Field.Base+import Math.Core.Field import Math.Combinatorics.FiniteGeometry import Math.Algebra.LinearAlgebra @@ -72,9 +72,13 @@  -- LATTICE OF (POSITIVE) DIVISORS OF N -divides a b = b `mod` a == 0+divides a b = b `rem` a == 0 -divisors n | n >= 1 = [a | a <- [1..n], a `divides` n]+divisors n = toSet [ d' | d <- takeWhile (\d -> d*d <= n) [1..],+                          let (q,r) = n `quotRem` d, r == 0,+                          d' <- [d,q] ]+-- The toSet call sorts, and deduplicates if n is a square+-- divisors n | n >= 1 = [a | a <- [1..n], a `divides` n]  -- |posetD n is the lattice of (positive) divisors of n posetD :: Int -> Poset Int@@ -108,7 +112,7 @@ partitions [] = [[]] partitions [x] = [[[x]]] partitions (x:xs) = let ps = partitions xs in-    map ([x]:) ps ++ [ (x:cell):(L.delete cell p) | p <- ps, cell <- p]+    map ([x]:) ps ++ [ (x:cell):p' | p <- ps, (cell,p') <- picks p] -- if the input is sorted, then so is the output  isRefinement a b = and [or [acell `isSubset` bcell | bcell <- b] | acell <- a]@@ -123,6 +127,7 @@   -- LATTICE OF INTERVAL PARTITIONS OF [1..N] ORDERED BY REFINEMENT+-- Interval partitions of [1..n] correspond to compositions of n  intervalPartitions xs = filter (all isInterval) (partitions xs) @@ -138,16 +143,12 @@  -- LATTICE OF INTEGER PARTITIONS OF N ORDERED BY REFINEMENT -integerPartitions1 n = ips (reverse [1..n]) n-    where ips [] 0 = [[]]-          ips [] _ = []-          ips (x:xs) n | x > n     = ips xs n-                       | otherwise = map (x:) (ips (x:xs) (n-x)) ++ ips xs n- -- For example, integerPartitions 5 -> [ [5], [4,1], [3,2], [3,1,1], [2,2,1], [2,1,1,1], [1,1,1,1,1] ]-integerPartitions n = dfs ([],n,n)-    where dfs (xs, 0, _) = [reverse xs]-          dfs (xs, r, i) = concatMap dfs [ (i':xs, r-i', i') | i' <- reverse [1..min r i] ]+integerPartitions n | n >= 0 = ips n n where+    ips 0 _ = [[]]+    ips _ 0 = []+    ips n m | m <= n = map (m:) (ips (n-m) m) ++ ips n (m-1)+            | otherwise = ips n n  isIPRefinement ys xs = dfs xs ys     where dfs (x:xs) (y:ys) | x < y = False@@ -187,7 +188,8 @@ -- This is the projective geometry PG(n,q) -- |posetL n fq is the lattice of subspaces of the vector space Fq^n, ordered by inclusion. -- Subspaces are represented by their reduced row echelon form.-posetL :: (Eq fq, FiniteField fq) => Int -> [fq] -> Poset [[fq]]+-- Example usage: posetL 2 f3+posetL :: (Eq fq, Num fq) => Int -> [fq] -> Poset [[fq]] posetL n fq = Poset ( subspaces fq n, isSubspace )   @@ -259,22 +261,26 @@     and [ x `poa` y == f x `pob` f y | x <- seta, y <- seta ]  -- Find all order isomorphisms between two posets--- This algorithm is faster to find out whether or not there are any+-- This is the most naive algorithm, and should not be used on larger posets+-- For example, already on the following, this takes forever compared to almost instant for orderIsos:+-- > head $ orderIsos01 (posetD $ 8*9*25*49) (posetD $ 4*27*25*121) orderIsos01 (Poset (seta,poa)) (Poset (setb,pob))     | length seta /= length setb = []     | otherwise = orderIsos' [] seta setb     where orderIsos' xys [] [] = [xys]           orderIsos' xys (x:xs) ys =-              concat [ orderIsos' ((x,y):xys) xs (L.delete y ys)-                     | y <- ys, and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ]+              concat [ orderIsos' ((x,y):xys) xs ys'+                     | (y,ys') <- picks ys,+                       and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ]  -- |Are the two posets order-isomorphic?-isOrderIso :: (Eq a, Eq b) => Poset a -> Poset b -> Bool-isOrderIso poseta posetb = (not . null) (orderIsos01 poseta posetb)+isOrderIso :: (Ord a, Ord b) => Poset a -> Poset b -> Bool+isOrderIso poseta posetb = (not . null) (orderIsos poseta posetb) --- Find all order isomorphisms between two posets--- This algorithm is faster to find all isomorphisms, if there are many--- (It may be that it is faster to find any, for large posets, but the break-even point seems to be quite big)+-- For small posets, it may be that the up-front cost of calculating the hasseDigraph and heightPartitionDAG+-- are not justified.+-- |Find all order isomorphisms between two posets+orderIsos :: (Ord a, Ord b) => Poset a -> Poset b -> [[(a,b)]] orderIsos posetA@(Poset (_,poa)) posetB@(Poset (_,pob))     | map length heightPartA /= map length heightPartB = []     | otherwise = dfs [] heightPartA heightPartB@@ -283,8 +289,9 @@           dfs xys [] [] = [xys]           dfs xys ([]:las) ([]:lbs) = dfs xys las lbs           dfs xys ((x:xs):las) (ys:lbs) =-              concat [ dfs ((x,y):xys) (xs:las) (L.delete y ys : lbs)-                     | y <- ys, and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ]+              concat [ dfs ((x,y):xys) (xs:las) (ys' : lbs)+                     | (y,ys') <- picks ys,+                       and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ] -- A variant on this algorithm would use the Hasse digraph rather than the partial order in the test on the last line -- This might be faster, depending how expensive the partial order comparison function is -- In effect though, it would then be a DAG isomorphism function

@@ -1,4 +1,4 @@--- Copyright (c) 2008, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, David Amos. All rights reserved.
 
 -- |A module defining various strongly regular graphs, including the Clebsch, Hoffman-Singleton, Higman-Sims, and McLaughlin graphs.
 --
@@ -8,6 +8,8 @@ --
 -- Strongly regular graphs are highly symmetric, and have large automorphism groups.
 module Math.Combinatorics.StronglyRegularGraph where
+
+import Prelude hiding ( (*>) )
 
 import qualified Data.List as L
 import Data.Maybe (isJust)

@@ -1,4 +1,4 @@--- Copyright (c) 2011, David Amos. All rights reserved.+-- Copyright (c) 2011-2015, David Amos. All rights reserved.  {-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances, DeriveFunctor #-} @@ -12,6 +12,8 @@ -- -- > [t,u,v,x,y,z] = map glexvar ["t","u","v","x","y","z"] module Math.CommutativeAlgebra.Polynomial where++import Prelude hiding ( (*>) )  import Math.Core.Field import Math.Core.Utils (toSet)

@@ -2,7 +2,8 @@  {-# LANGUAGE GeneralizedNewtypeDeriving #-} --- |A module defining the field Q of rationals and the small finite fields F2, F3, F4, F5, F7, F8, F9, F11, F13, F16, F17, F19, F23, F25.+-- |A module defining the field Q of rationals and the small finite fields (Galois fields)+-- F2, F3, F4, F5, F7, F8, F9, F11, F13, F16, F17, F19, F23, F25. -- -- Given a prime power q, Fq is the type representing elements of the field (eg @F4@), -- fq is a list of the elements of the field, beginning 0,1,... (eg @f4@),@@ -380,10 +381,10 @@  instance Num F9 where     F9 x + F9 y = F9 $ z1 + z0-        where z = x+y; z1 = (z .&. 0xff00) `mod` 0x300; z0 = (z .&. 0xff) `mod` 3+        where z = x+y; z1 = (z .&. 0xff00) `rem` 0x300; z0 = (z .&. 0xff) `rem` 3     negate (F9 x) = F9 $ z1 + z0-        where z = 0x303 - x; z1 = (z .&. 0xff00) `mod` 0x300; z0 = (z .&. 0xff) `mod` 3-    F9 x * F9 y = F9 $ ((z2 + z1) `mod` 0x300) + ((z2 + z0) `mod` 3) +        where z = 0x303 - x; z1 = (z .&. 0xff00) `rem` 0x300; z0 = (z .&. 0xff) `rem` 3+    F9 x * F9 y = F9 $ ((z2 + z1) `rem` 0x300) + ((z2 + z0) `rem` 3)          where z = x*y; z2 = z .&. 0xff0000; z1 = z .&. 0xff00; z0 = z .&. 0xff         -- Explanation: We are substituting x^2 = x+1.         -- We could do z2 `shiftR` 8 and z2 `shiftR` 16@@ -470,10 +471,10 @@  instance Num F25 where     F25 x + F25 y = F25 $ z1 + z0-        where z = x+y; z1 = (z .&. 0xff00) `mod` 0x500; z0 = (z .&. 0xff) `mod` 5+        where z = x+y; z1 = (z .&. 0xff00) `rem` 0x500; z0 = (z .&. 0xff) `rem` 5     negate (F25 x) = F25 $ z1 + z0-        where z = 0x505 - x; z1 = (z .&. 0xff00) `mod` 0x500; z0 = (z .&. 0xff) `mod` 5-    F25 x * F25 y = F25 $ ((z2 + z1) `mod` 0x500) + ((3*z2 + z0) `mod` 5) +        where z = 0x505 - x; z1 = (z .&. 0xff00) `rem` 0x500; z0 = (z .&. 0xff) `rem` 5+    F25 x * F25 y = F25 $ ((z2 + z1) `rem` 0x500) + ((3*z2 + z0) `rem` 5)          where z = x*y; z2 = z .&. 0xff0000; z1 = z .&. 0xff00; z0 = z .&. 0xff         -- Explanation: We are substituting x^2 = x+3.         -- We could do z2 `shiftR` 8 and z2 `shiftR` 16

@@ -26,6 +26,25 @@     GT -> y : setUnionAsc (x:xs) ys setUnionAsc xs ys = xs ++ ys +setUnionDesc :: Ord a => [a] -> [a] -> [a]+setUnionDesc (x:xs) (y:ys) =+    case compare x y of+    GT -> x : setUnionDesc xs (y:ys)+    EQ -> x : setUnionDesc xs ys+    LT -> y : setUnionDesc (x:xs) ys+setUnionDesc xs ys = xs ++ ys++-- |The (multi-)set intersection of two ascending lists. If both inputs are strictly increasing,+-- then the output is the set intersection and is strictly increasing. If both inputs are weakly increasing,+-- then the output is the multiset intersection (with multiplicity), and is weakly increasing.+intersectAsc :: Ord a => [a] -> [a] -> [a]+intersectAsc (x:xs) (y:ys) =+    case compare x y of+    LT -> intersectAsc xs (y:ys)+    EQ -> x : intersectAsc xs ys+    GT -> intersectAsc (x:xs) ys+intersectAsc _ _ = []+ -- |The multiset sum of two ascending lists. If xs and ys are ascending, then multisetSumAsc xs ys == sort (xs++ys). -- The code does not check that the lists are ascending. multisetSumAsc :: Ord a => [a] -> [a] -> [a]@@ -36,7 +55,7 @@     GT -> y : multisetSumAsc (x:xs) ys multisetSumAsc xs ys = xs ++ ys --- |The multiset sum of two descending lists. If xs and ys are descending, then multisetSumDesc xs ys == sort (xs++ys).+-- |The multiset sum of two descending lists. If xs and ys are descending, then multisetSumDesc xs ys == sortDesc (xs++ys). -- The code does not check that the lists are descending. multisetSumDesc :: Ord a => [a] -> [a] -> [a] multisetSumDesc (x:xs) (y:ys) =@@ -79,6 +98,33 @@     GT -> isSubMultisetAsc (x:xs) ys isSubMultisetAsc [] ys = True isSubMultisetAsc xs [] = False++-- |Is the element in the ascending list?+--+-- With infinite lists, this can fail to terminate.+-- For example, elemAsc 1 [1/2,3/4,7/8..] would fail to terminate.+-- However, with a list of Integer, this will always terminate.+elemAsc :: Ord a => a -> [a] -> Bool+elemAsc x (y:ys) = case compare x y of+                   LT -> False+                   EQ -> True+                   GT -> elemAsc x ys+-- or x `elemAsc` ys = x `elem` takeWhile (<= x) ys++-- |Is the element not in the ascending list? (With infinite lists, this can fail to terminate.)+notElemAsc :: Ord a => a -> [a] -> Bool+notElemAsc x (y:ys) = case compare x y of+                      LT -> True+                      EQ -> False+                      GT -> notElemAsc x ys+++-- From Conor McBride+-- http://stackoverflow.com/questions/12869097/splitting-list-into-a-list-of-possible-tuples/12872133#12872133+-- |Return all the ways to \"pick one and leave the others\" from a list+picks :: [a] -> [(a,[a])]+picks [] = []+picks (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- picks xs]   pairs (x:xs) = map (x,) xs ++ pairs xs

@@ -1,16 +1,72 @@ -- Copyright (c) 2006-2011, David Amos. All rights reserved. +{-# LANGUAGE BangPatterns #-}+ -- |A module for finding prime factors. module Math.NumberTheory.Factor (module Math.NumberTheory.Prime,-                                 pfactors) where+                                 pfactors, ppfactors, pfactorsTo, ppfactorsTo) where -import Math.NumberTheory.Prime+import Control.Arrow (second, (&&&)) import Data.Either (lefts)-import Data.List (zip4)+import Data.List as L+import Math.Core.Utils (multisetSumAsc)+import Math.NumberTheory.Prime --- Cohen, A Course in Computational Algebraic Number Theory, p488+-- |List the prime factors of n (with multiplicity). For example:+-- >>> pfactors 60+-- [2,2,3,5]+--+-- This says that 60 = 2 * 2 * 3 * 5+-- +-- The algorithm uses trial division to find small factors,+-- followed if necessary by the elliptic curve method to find larger factors.+-- The running time increases with the size of the second largest prime factor of n.+-- It can find 10-digit prime factors in seconds, but can struggle with 20-digit prime factors.+pfactors :: Integer -> [Integer]+pfactors n | n > 0 = pfactors' n $ takeWhile (< 10000) primes+           | n < 0 = -1 : pfactors' (-n) (takeWhile (< 10000) primes)+    where pfactors' n (d:ds) | n == 1 = []+                             | n < d*d = [n]+                             | r == 0 = d : pfactors' q (d:ds)+                             | otherwise = pfactors' n ds+                             where (q,r) = quotRem n d+          pfactors' n [] = pfactors'' n+          pfactors'' n = if isMillerRabinPrime n then [n]+                         else let d = findFactorParallelECM n -- findFactorECM n+                              in multisetSumAsc (pfactors'' d) (pfactors'' (n `div` d)) +-- |List the prime power factors of n. For example:+-- >>> ppfactors 60+-- [(2,2),(3,1),(5,1)]+--+-- This says that 60 = 2^2 * 3^1 * 5^1+ppfactors :: Integer -> [(Integer,Int)]+ppfactors = map (head &&& length) . L.group . pfactors+-- ppfactors = map (\xs -> (head xs, length xs)) . L.group . pfactors +-- |Find the prime factors of all numbers up to n. Thus @pfactorsTo n@ is equivalent to @[(m, pfactors m) | m <- [1..n]]@,+-- except that the results are not returned in order. For example:+-- >>> pfactorsTo 10+-- [(8,[2,2,2]),(4,[2,2]),(6,[3,2]),(10,[5,2]),(2,[2]),(9,[3,3]),(3,[3]),(5,[5]),(7,[7]),(1,[])]+--+-- @pfactorsTo n@ is significantly faster than @map pfactors [1..n]@ for larger n.+pfactorsTo n = pfactorsTo' (1,[]) primes where+    pfactorsTo' (!m,!qs) ps@(ph:pt) | m' > n = [(m,qs)]+                                    | otherwise = pfactorsTo' (m',ph:qs) ps ++ pfactorsTo' (m,qs) pt+        where m' = m*ph+-- We avoid a reverse call, because it does make a noticeable difference to the speed.++-- |Find the prime power factors of all numbers up to n. Thus @ppfactorsTo n@ is equivalent to @[(m, ppfactors m) | m <- [1..n]]@,+-- except that the results are not returned in order. For example:+-- >>> ppfactorsTo 10+-- [(8,[(2,3)]),(4,[(2,2)]),(6,[(3,1),(2,1)]),(10,[(5,1),(2,1)]),(2,[(2,1)]),(9,[(3,2)]),(3,[(3,1)]),(5,[(5,1)]),(7,[(7,1)]),(1,[])]+--+-- @ppfactorsTo n@ is significantly faster than @map ppfactors [1..n]@ for larger n.+ppfactorsTo = map (second (map (head &&& length) . L.group)) . pfactorsTo+++-- Cohen, A Course in Computational Algebraic Number Theory, p488+ -- return (u,v,d) s.t ua+vb = d, with d = gcd a b extendedEuclid a b     | b == 0 = (1,0,a)@@ -18,7 +74,6 @@                       (s,t,d) = extendedEuclid b r -- s*b+t*r == d                   in (t,s-q*t,d)                   -- s*b+t*(a-q*b) == d - -- ELLIPTIC CURVE ARITHMETIC  data EllipticCurve = EC Integer Integer Integer deriving (Eq, Show)@@ -66,7 +121,6 @@                               (Left _, _) -> p'                               (_, Left _) -> q' - -- ELLIPTIC CURVE FACTORISATION  -- We choose an elliptic curve E over Zn, and a point P on the curve@@ -74,7 +128,6 @@ -- What we are hoping is that at some stage we will fail because we can't invert an element in Zn -- This will lead to finding a non-trivial factor of n - discriminantEC a b = 4 * a * a * a + 27 * b * b  -- perform a sequence of scalar multiplications in the elliptic curve, hoping for a bailout@@ -105,32 +158,7 @@     -- the filter is because d might be a multiple of n,     -- for example if the problem was that the discriminant was divisible by n ---- |List the prime factors of n (with multiplicity).--- The algorithm uses trial division to find small factors,--- followed if necessary by the elliptic curve method to find larger factors.--- The running time increases with the size of the second largest prime factor of n.--- It can find 10-digit prime factors in seconds, but can struggle with 20-digit prime factors.-pfactors :: Integer -> [Integer]-pfactors n | n > 0 = pfactors' n $ takeWhile (< 10000) primes-           | n < 0 = -1 : pfactors' (-n) (takeWhile (< 10000) primes)-    where pfactors' n (d:ds) | n == 1 = []-                             | n < d*d = [n]-                             | r == 0 = d : pfactors' q (d:ds)-                             | otherwise = pfactors' n ds-                             where (q,r) = quotRem n d-          pfactors' n [] = pfactors'' n-          pfactors'' n = if isMillerRabinPrime n then [n]-                         else let d = findFactorParallelECM n -- findFactorECM n-                              in merge (pfactors'' d) (pfactors'' (n `div` d))--merge (x:xs) (y:ys) =-    case compare x y of-    LT -> x : merge xs (y:ys)-    EQ -> x : y : merge xs ys-    GT -> y : merge (x:xs) ys-merge xs ys = xs ++ ys-+-- TESTING MULTIPLE CURVES IN PARALLEL  -- Cohen p489 -- find inverse of as mod n in parallel, or a non-trivial factor of n@@ -144,7 +172,7 @@  parallelEcAdd n ecs ps1 ps2 =     case parallelInverse n (zipWith f ps1 ps2) of-    Right invs -> Right [g ec p1 p2 inv | (ec,p1,p2,inv) <- zip4 ecs ps1 ps2 invs]+    Right invs -> Right [g ec p1 p2 inv | (ec,p1,p2,inv) <- L.zip4 ecs ps1 ps2 invs]     Left d -> Left d     where f Inf pt = 1           f pt Inf = 1

@@ -11,16 +11,24 @@   isTrialDivisionPrime n-    | n > 1 = isNotDivisibleBy primes+    | n > 1 = not $ any (\p -> n `rem` p == 0) (takeWhile (\p -> p*p <= n) primes)     | otherwise = False-    where isNotDivisibleBy (d:ds) | d*d > n         = True-                                  | n `rem` d == 0  = False-                                  | otherwise       = isNotDivisibleBy ds  -- |A (lazy) list of the primes primes :: [Integer]-primes = 2 : 3 : filter isTrialDivisionPrime (concat [ [m6-1,m6+1] | m6 <- [6,12..] ])+primes = 2 : filter isPrime [3,5..] where+    isPrime n = not $ any (\p -> n `rem` p == 0) (takeWhile (\p -> p*p <= n) primes) +{-+-- This is just marginally faster, but less elegant+primes2 :: [Integer]+primes2 = 2 : 3 : 5 : 7 : filter isPrime+    (concat [ [m30+11,m30+13,m30+17,m30+19,m30+23,m30+29,m30+31,m30+37] | m30 <- [0,30..] ])+    where isPrime n = not $ any (\p -> n `rem` p == 0) (takeWhile (\p -> p*p <= n) primes2')+          primes2' = drop 3 primes2+-}++{- -- initial version. This isn't going to be very good if n has any "large" prime factors (eg > 10000) pfactors1 n | n > 0 = pfactors' n primes             | n < 0 = -1 : pfactors' (-n) primes@@ -29,7 +37,7 @@                              | r == 0 = d : pfactors' q (d:ds)                              | otherwise = pfactors' n ds                              where (q,r) = quotRem n d-+-}  -- MILLER-RABIN TEST -- Cohen, A Course in Computational Algebraic Number Theory, p422@@ -59,7 +67,8 @@ -- power_mod b t n == b^t mod n power_mod b t n = powerMod' b 1 t     where powerMod' x y 0 = y-          powerMod' x y t = powerMod' (x*x `rem` n) (if even t then y else x*y `rem` n) (t `div` 2)+          powerMod' x y t = let (q,r) = t `quotRem` 2+                            in powerMod' (x*x `rem` n) (if r == 0 then y else x*y `rem` n) q  isMillerRabinPrime' n     | n >= 4 =@@ -106,6 +115,7 @@             where n6 = (n `div` 6) * 6                   candidates = dropWhile (<= n) $ concat [ [m6+1,m6+5] | m6 <- [n6, n6+6..] ] +{- -- slightly better version. This is okay so long as n has at most one "large" prime factor (> 10000) -- if it has more, it does at least tell you, via an error message, that it has run into difficulties pfactors2 n | n > 0 = pfactors' n $ takeWhile (< 10000) primes@@ -116,4 +126,4 @@                              | otherwise = pfactors' n ds                              where (q,r) = quotRem n d           pfactors' n [] = if isMillerRabinPrime n then [n] else error ("pfactors2: can't factor " ++ show n)-+-}

@@ -1,4 +1,4 @@--- Copyright (c) 2011, David Amos. All rights reserved.+-- Copyright (c) 2011-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances, OverlappingInstances #-} @@ -16,7 +16,7 @@ -- > i * sqrt(-3) module Math.NumberTheory.QuadraticField where -import Prelude hiding (sqrt)+import Prelude hiding (sqrt, (*>) )  import Data.List as L import Math.Core.Field

@@ -1,6 +1,8 @@--- Copyright (c) 2008, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Classical where
+
+import Prelude hiding ( (*>) )
 
 import Math.Algebra.Field.Base
 import Math.Algebra.Field.Extension hiding ( (<+>), (<*>) )

@@ -1,7 +1,9 @@--- Copyright (c) 2008, David Amos. All rights reserved.
+-- Copyright (c) 2008-2015, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Exceptional where
 
+import Prelude hiding ( (*>) )
+
 import Data.List as L
 
 -- import Math.Algebra.Field.Base
@@ -165,7 +167,7 @@ 
 -- Unit imaginary octonions form one orbit under the action of G2
 
--- [alpha', beta', gamma' generate G2(3) as a permutation group on 702 points (the number of unit imaginary octonions over F3)
+-- [alpha', beta', gamma'] generate G2(3) as a permutation group on 702 points (the number of unit imaginary octonions over F3)
 -- Interestingly, http://brauer.maths.qmul.ac.uk/Atlas/v3/exc/G23/ doesn't seem to have this permutation representation
 
 

@@ -1,6 +1,8 @@--- Copyright (c) David Amos, 2009. All rights reserved.+-- Copyright (c) David Amos, 2009-2015. All rights reserved.  module Math.Projects.MiniquaternionGeometry where++import Prelude hiding ( (*>) )  import qualified Data.List as L 

@@ -1,6 +1,8 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) David Amos, 2008-2015. All rights reserved.
 
 module Math.Projects.RootSystem where
+
+import Prelude hiding ( (*>) )
 
 import Data.Ratio
 import qualified Data.List as L

@@ -1,10 +1,12 @@--- Copyright (c) David Amos, 2010. All rights reserved.+-- Copyright (c) David Amos, 2010-2015. All rights reserved.  {-# LANGUAGE TypeFamilies, EmptyDataDecls #-}   module Math.QuantumAlgebra.OrientedTangle where +import Prelude hiding ( (*>), (<*>) )+ import Math.Algebra.Field.Base import Math.Algebras.LaurentPoly -- hiding (lvar, q, q') @@ -26,7 +28,7 @@ data OrientedTangle  -- In GHCi 6.12.1, we appear to be limited to 8 value constructors within an associated data family-instance Category OrientedTangle where+instance MCategory OrientedTangle where     data Ob OrientedTangle = OT [Oriented] deriving (Eq,Ord,Show)     data Ar OrientedTangle = IdT [Oriented]                            | CapT HorizDir@@ -54,7 +56,7 @@     target (SeqT as) = target (last as)     a >>> b | target a == source b = SeqT [a,b] -instance TensorCategory OrientedTangle where+instance Monoidal OrientedTangle where     tunit = OT []     tob (OT as) (OT bs) = OT (as++bs)     tar a b = ParT [a,b]@@ -64,11 +66,11 @@ idV = id idV' = id -evalV  = \(E i, E j) -> if i + j == 0 then return () else zero-evalV' = \(E i, E j) -> if i + j == 0 then return () else zero+evalV  = \(E i, E j) -> if i + j == 0 then return () else zerov+evalV' = \(E i, E j) -> if i + j == 0 then return () else zerov -coevalV  m = foldl (<+>) zero [e i `te` e (-i) | i <- [1..m] ]-coevalV' m = foldl (<+>) zero [e (-i) `te` e i | i <- [1..m] ]+coevalV  m = foldl (<+>) zerov [e i `te` e (-i) | i <- [1..m] ]+coevalV' m = foldl (<+>) zerov [e (-i) `te` e i | i <- [1..m] ]  lambda m = q' ^ m -- q^-m 

@@ -53,7 +53,7 @@     powers (Aq20 m) = powers m  instance Algebra (LaurentPoly Q) (Aq20 String) where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)] where munit = Aq20 (NCM 0 [])     mult x = x''' where         x' = mult $ fmap ( \(Aq20 a, Aq20 b) -> (a,b) ) x -- unwrap and multiply@@ -75,7 +75,7 @@     powers (Aq02 m) = powers m  instance Algebra (LaurentPoly Q) (Aq02 String) where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)] where munit = Aq02 (NCM 0 [])     mult x = x''' where         x' = mult $ fmap ( \(Aq02 a, Aq02 b) -> (a,b) ) x -- unwrap and multiply@@ -99,7 +99,7 @@     powers (M2q m) = powers m  instance Algebra (LaurentPoly Q) (M2q String) where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)] where munit = M2q (NCM 0 [])     mult x = x''' where         x' = mult $ fmap ( \(M2q a, M2q b) -> (a,b) ) x -- unwrap and multiply@@ -188,7 +188,7 @@     powers (SL2q m) = powers m  instance Algebra (LaurentPoly Q) (SL2q String) where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)] where munit = SL2q (NCM 0 [])     mult x = x''' where         x' = mult $ fmap ( \(SL2q a, SL2q b) -> (a,b) ) x -- unwrap and multiply

@@ -1,4 +1,4 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE MultiParamTypeClasses, TypeFamilies, FlexibleInstances, EmptyDataDecls #-} @@ -15,7 +15,7 @@ import Math.Algebra.Field.Base import Math.Algebras.LaurentPoly -import Math.QuantumAlgebra.TensorCategory+import Math.QuantumAlgebra.TensorCategory hiding (Vect)   instance Mon [a] where@@ -25,7 +25,7 @@ -- type TensorAlgebra k a = Vect k [a]  instance (Eq k, Num k, Ord a) => Algebra k [a] where-    unit 0 = zero -- V []+    unit 0 = zerov -- V []     unit x = V [(munit,x)]     mult = nf . fmap (\(a,b) -> a `mmult` b) @@ -37,7 +37,7 @@  data Tangle -instance Category Tangle where+instance MCategory Tangle where     data Ob Tangle = OT Int deriving (Eq,Ord,Show)     data Ar Tangle = IdT Int                    | CapT@@ -71,7 +71,7 @@ --    a >>> b | target a == source b = SeqT a b     a >>> b | target a == source b = SeqT [a,b] -instance TensorCategory Tangle where+instance Monoidal Tangle where     tunit = OT 0     tob (OT a) (OT b) = OT (a+b) --    tar a b = ParT a b@@ -93,7 +93,7 @@ cup :: [Oriented] -> TangleRep [Oriented] cup [Plus, Minus] = (-q'^2) *> return [] cup [Minus, Plus] = return []-cup _ = zero+cup _ = zerov  -- also called xminus over :: [Oriented] -> TangleRep [Oriented]

@@ -1,14 +1,13 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2014, David Amos. All rights reserved. -{-# LANGUAGE TypeFamilies, EmptyDataDecls #-}+{-# LANGUAGE TypeFamilies, EmptyDataDecls, MultiParamTypeClasses #-} --- |A module defining classes and example instances of categories and tensor categories+-- |A module defining classes and example instances of categories, monoidal categories and braided categories module Math.QuantumAlgebra.TensorCategory where -import Math.Algebra.Group.PermutationGroup-+import Data.List as L -class Category c where+class MCategory c where     data Ob c :: *     data Ar c :: *     id_ :: Ob c -> Ar c@@ -18,73 +17,169 @@ -- whereas we want the objects to be values of a single type.  +class (MCategory a, MCategory b) => MFunctor a b where+    fob :: Ob a -> Ob b -- functor on objects+    far :: Ar a -> Ar b -- functor on arrows++-- We could also define tensor functors and braided functors, which are just functors which commute appropriately+-- with tensor and braiding operations++ -- Kassel p282--- The following is actually definition of a strict tensor category-class Category c => TensorCategory c where+-- The following is actually definition of a _strict_ monoidal category+-- Also called tensor category+class MCategory c => Monoidal c where     tunit :: Ob c     tob :: Ob c -> Ob c -> Ob c -- tensor product of objects     tar :: Ar c -> Ar c -> Ar c -- tensor product of arrows -class TensorCategory c => StrictTensorCategory c where {}+class Monoidal c => StrictMonoidal c where {} -- we want to be able to declare some tensor categories as strict -class TensorCategory c => WeakTensorCategory c where-    assoc :: Ob c -> Ob c -> Ob c -- (u `tob` v) `tob` w -> u `tob` (v `tob` w)-    lunit :: Ob c -> Ob c         -- unit `tob` v -> v-    runit :: Ob c -> Ob c         -- v `tob` unit -> v+class Monoidal c => WeakMonoidal c where+    assoc :: Ob c -> Ob c -> Ob c -> Ar c -- assoc u v w is an arrow (natural transformation?): (u `tob` v) `tob` w -> u `tob` (v `tob` w)+    lunit :: Ob c -> Ar c                 -- lunit v is an arrow (isomorphism): tunit `tob` v -> v+    runit :: Ob c -> Ar c                 -- runit v is an arrow (isomorphism): v `tob` tunit -> v  {--instance (TensorCategory c, Eq (Ar c), Show (Ar c)) => Num (Ar c) where+instance (Monoidal c, Eq (Ar c), Show (Ar c)) => Num (Ar c) where     (*) = tar -} +class Monoidal c => Braided c where+    twist :: Ob c -> Ob c -> Ar c+    -- twist v w is a map from v tensor w to w tensor v+    -- twist must be natural, and satisfy certain commutative diagrams - Kock 161, 169 --- SYMMETRIC GROUPOID+class Braided c => Symmetric c where {}+-- if twist satisfies twist v w >>> twist w v == id_ (v tensor w), then the category is symmetric -data SymmetricGroupoid -instance Category SymmetricGroupoid where-    data Ob SymmetricGroupoid = OS Int deriving (Eq,Ord,Show)-    data Ar SymmetricGroupoid = AS Int (Permutation Int) deriving (Eq,Ord,Show)-    id_ (OS n) = AS n 1-    source (AS n _) = OS n-    target (AS n _) = OS n-    AS m g >>> AS n h | m == n = AS m (g*h)+-- SIMPLEX CATEGORIES -instance TensorCategory SymmetricGroupoid where-    tunit = OS 0-    tob (OS m) (OS n) = OS (m+n)-    tar (AS m g) (AS n h) = AS (m+n) (g * h~^k)-        where k = p [[1..m+n]] ^ m---    tar (AS m g) (AS n h) = AS (m+n) (fromPairs $ toPairs g ++ map (\(x,y)->(x+m,y+m)) (toPairs h))+-- Kock, Frobenius Algebras ..., p178-9+-- The skeleton of FinOrd (finite ordered sets)+-- The objects are the finite ordinals n == [0..n-1]+-- The arrows are the order-preserving maps+data FinOrd +instance MCategory FinOrd where+    data Ob FinOrd = FinOrdOb Int deriving (Eq,Ord,Show)+    -- FinOrdOb n represents the oriented simplex n == [0..n-1]+    data Ar FinOrd = FinOrdAr Int Int [Int] deriving (Eq,Ord,Show)+    -- FinOrdAr s t fs represents the order-preserving map, zip [0..s-1] fs.+    -- For example FinOrdAr 3 2 [0,0,1] represents the map 0 -> 0, 1 -> 0, 2 -> 1+    id_ (FinOrdOb n) = FinOrdAr n n [0..n-1]+    source (FinOrdAr s _ _) = FinOrdOb s+    target (FinOrdAr _ t _) = FinOrdOb t+    FinOrdAr sf tf fs >>> FinOrdAr sg tg gs | tf == sg = FinOrdAr sf tg [let j = fs !! i in gs !! j | i <- [0..sf-1] ] +instance Monoidal FinOrd where+    tunit = FinOrdOb 0+    tob (FinOrdOb m) (FinOrdOb n) = FinOrdOb (m+n)+    tar (FinOrdAr sf tf fs) (FinOrdAr sg tg gs) = FinOrdAr (sf+sg) (tf+tg) (fs ++ map (+tf) gs)++finOrdAr s t fs | s == length fs && minimum fs >= 0 && maximum fs < t && isOrderPreserving fs+    = FinOrdAr s t fs+    where isOrderPreserving (f1:f2:fs) = f1 <= f2 && isOrderPreserving (f2:fs)+          isOrderPreserving _ = True+++-- The skeleton of FinSet+-- The objects are the finite cardinals n == {0..n-1} (with no order)+-- The arrows are the maps+data FinCard++instance MCategory FinCard where+    data Ob FinCard = FinCardOb Int deriving (Eq,Ord,Show)+    -- FinCardOb n represents the unoriented simplex n == {0..n-1}+    data Ar FinCard = FinCardAr Int Int [Int] deriving (Eq,Ord,Show)+    -- FinCardAr s t fs represents the map, zip [0..s-1] fs.+    -- For example FinCardAr 3 2 [0,1,0] represents the map 0 -> 0, 1 -> 1, 2 -> 0+    id_ (FinCardOb n) = FinCardAr n n [0..n-1]+    source (FinCardAr s _ _) = FinCardOb s+    target (FinCardAr _ t _) = FinCardOb t+    FinCardAr sf tf fs >>> FinCardAr sg tg gs | tf == sg = FinCardAr sf tg [let j = fs !! i in gs !! j | i <- [0..sf-1] ]++instance Monoidal FinCard where+    tunit = FinCardOb 0+    tob (FinCardOb m) (FinCardOb n) = FinCardOb (m+n)+    tar (FinCardAr sf tf fs) (FinCardAr sg tg gs) = FinCardAr (sf+sg) (tf+tg) (fs ++ map (+tf) gs)++finCardAr s t fs | s == length fs && minimum fs >= 0 && maximum fs < t -- for finite cardinals, the map doesn't have to be order-preserving+    = FinCardAr s t fs++-- Finite permutations form a subcategory of FinCard+-- having as objects the finite cardinals n == {0..n-1}+-- and as arrows the bijections (== permutations)+finPerm fs | L.sort fs == [0..n-1] = FinCardAr n n fs+    where n = length fs+-- (Note that these are permutations of [0..n-1], rather than [1..n])+++-- This is the forgetful functor FinOrd -> FinCard (FinSet)+instance MFunctor FinOrd FinCard where+    fob (FinOrdOb n) = FinCardOb n+    far (FinOrdAr s t fs) = FinCardAr s t fs++ -- BRAID CATEGORY  data Braid -instance Category Braid where-    data Ob Braid = OB Int deriving (Eq,Ord,Show)-    data Ar Braid = AB Int [Int] deriving (Eq,Ord,Show)-    id_ (OB n) = AB n []-    source (AB n _) = OB n-    target (AB n _) = OB n-    AB m is >>> AB n js | m == n = AB m (is ++ js)+instance MCategory Braid where+    data Ob Braid = BraidOb Int deriving (Eq,Ord,Show)+    data Ar Braid = BraidAr Int [Int] deriving (Eq,Ord,Show)+    id_ (BraidOb n) = BraidAr n []+    source (BraidAr n _) = BraidOb n+    target (BraidAr n _) = BraidOb n+    BraidAr m is >>> BraidAr n js | m == n = BraidAr m (cancel (reverse is) js)+        where cancel (x:xs) (y:ys) = if x+y == 0 then cancel xs ys else reverse xs ++ x:y:ys+              cancel xs ys = reverse xs ++ ys -s n i | 0 < i && i < n = AB n [i]+t n 0 = BraidAr n [] -- the identity braid+t n i |  0 < i && i < n = BraidAr n [i]+      | -n < i && i < 0 = BraidAr n [i]+-- The generators of B_n are [t n i | i <- [1..n-1]] -instance TensorCategory Braid where-    tunit = OB 0-    tob (OB a) (OB b) = OB (a+b)-    tar (AB m is) (AB n js) = AB (m+n) (is ++ map (+m) js)+-- The inverses of the braid generators+t' n i | 0 < i && i < n = BraidAr n [-i] +instance Monoidal Braid where+    tunit = BraidOb 0+    tob (BraidOb m) (BraidOb n) = BraidOb (m+n)+    tar (BraidAr m is) (BraidAr n js) = BraidAr (m+n) (is ++ map (+m) js) +instance Braided Braid where+    twist (BraidOb m) (BraidOb n) = BraidAr (m+n) $ concat [[i..i+n-1] | i <- [m,m-1..1]] +-- Note that in FinCard we consider the objects as [0..n-1], whereas in Braid we consider them as [1..n], so that s_i twists [i,i+1]+instance MFunctor Braid FinCard where+    fob (BraidOb n) = FinCardOb n+    far (BraidAr n ss) = foldr (>>>) (id_ (FinCardOb n)) [finPerm ([0..ti-1] ++ [ti+1,ti] ++ [ti+2..n-1]) | si <- ss, let ti = abs si - 1] ++-- VECT++data Vect k++instance Num k => MCategory (Vect k) where+    data Ob (Vect k) = VectOb Int deriving (Eq,Ord,Show)+    data Ar (Vect k) = VectAr Int Int [[Int]] deriving (Eq,Ord,Show)+    id_ (VectOb n) = VectAr n n idMx where idMx = [[if i == j then 1 else 0 | j <- [1..n]] | i <- [1..n]]+    source (VectAr m _ _) = VectOb m+    target (VectAr _ n _) = VectOb n+    VectAr r c xss >>> VectAr r' c' yss | c == r' = undefined -- matrix multiplication++-- functor from FinPerm to Vect k+++-- 2-COBORDISMS+ data Cob2 -- works very similar to Tangle category -instance Category Cob2 where+instance MCategory Cob2 where     data Ob Cob2 = O Int deriving (Eq,Ord,Show)     data Ar Cob2 = Id Int                  | Unit@@ -111,7 +206,7 @@     target (Seq a b) = target b     a >>> b | target a == source b = Seq a b -instance TensorCategory Cob2 where+instance Monoidal Cob2 where     tunit = O 0     tob (O a) (O b) = O (a+b)     tar a b = Par a b

@@ -1,6 +1,8 @@--- Copyright (c) 2011, David Amos. All rights reserved.+-- Copyright (c) 2011-2015, David Amos. All rights reserved.  module Math.Test.TAlgebras.TOctonions where++import Prelude hiding ( (*>) )  import Test.QuickCheck 

@@ -6,6 +6,8 @@  module Math.Test.TAlgebras.TStructures where +import Prelude hiding ( (*>) )+ -- import Test.QuickCheck -- don't actually need, as we don't define any Arbitrary instances here 

@@ -1,9 +1,11 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE EmptyDataDecls, ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies, RankNTypes #-}  module Math.Test.TAlgebras.TTensorProduct where++import Prelude hiding ( (*>) )  import Test.QuickCheck import Math.Algebras.VectorSpace

@@ -1,10 +1,12 @@--- Copyright (c) 2010, David Amos. All rights reserved.+-- Copyright (c) 2010-2015, David Amos. All rights reserved.  {-# LANGUAGE FlexibleInstances, NoMonomorphismRestriction, ScopedTypeVariables, GeneralizedNewtypeDeriving #-}   module Math.Test.TAlgebras.TVectorSpace where +import Prelude hiding ( (*>) )+ import Test.QuickCheck import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct@@ -17,8 +19,8 @@ prop_AddGrp (x,y,z) =     x <+> (y <+> z) == (x <+> y) <+> z && -- associativity     x <+> y == y <+> x                 && -- commutativity-    x <+> zero == x                    && -- identity-    x <+> neg x == zero                   -- inverse+    x <+> zerov == x                   && -- identity+    x <+> negatev x == zerov              -- inverse  prop_VecSp (a,b,x,y,z) =     prop_AddGrp (x,y,z) &&@@ -51,9 +53,9 @@   prop_Linear f (a,x,y) =-    f (x <+> y) == f x <+> f y &&-    f zero == zero             &&-    f (neg x) == neg (f x)     &&+    f (x <+> y) == f x <+> f y       &&+    f zerov == zerov                 &&+    f (negatev x) == negatev (f x)   &&     f (a *> x) == a *> f x  prop_LinearQn f (a,x,y) = prop_Linear f (a,x,y)@@ -73,8 +75,12 @@ -- DIRECT SUM  {-+instance Num k => Alternative (Vect k) where+    (<|>) = mplus+    empty = mzero+ instance Num k => MonadPlus (Vect k) where-    mzero = zero+    mzero = zerov     mplus (V xs) (V ys) = V (xs++ys) -- need to call nf afterwards -} @@ -96,7 +102,7 @@ linfun avbs = linear f where     f a = case lookup a avbs of           Just vb -> vb-          Nothing -> zero+          Nothing -> zerov   prop_Product (f',g',x) =@@ -174,7 +180,7 @@     ((a, b) -> Vect k c) -> Vect k (Either a b) -> Vect k c bilinear f = linear f . tensor -dot = bilinear (\(a,b) -> if a == b then return () else zero)+dot = bilinear (\(a,b) -> if a == b then return () else zerov)  polymult = bilinear (\(E i, E j) -> return (E (i+j))) 

@@ -6,10 +6,14 @@ import Math.Core.Field hiding (f7) import Math.Core.Utils (combinationsOf) import Math.Algebra.Group.PermutationGroup as P+import Math.Algebra.Group.RandomSchreierSims as SS import Math.Combinatorics.Graph as G import Math.Combinatorics.GraphAuts import Math.Combinatorics.Matroid as M+import Math.Combinatorics.FiniteGeometry +import qualified Math.Algebra.Field.Extension as F+ import Test.HUnit  @@ -20,7 +24,18 @@     testlistIncidenceAutsOrder,     testlistGraphIsos,     testlistIsGraphIso,-    testlistIncidenceIsos+    testlistIncidenceIsos,+    testlistVertexTransitive,+    testlistNotVertexTransitive,+    testlistEdgeTransitive,+    testlistNotEdgeTransitive,+    testlistArcTransitive,+    testlistNotArcTransitive,+    testlist2ArcTransitive,+    testlistNot2ArcTransitive,+    testlist4ArcTransitive,+    testlistDistanceTransitive,+    testlistNotDistanceTransitive     ]  @@ -125,3 +140,121 @@                           (G [Left 3, Left 4, Right 4] [[Left 4, Right 4]])                           [[(1,4),(2,3)]]     ]++testcaseVertexTransitive (desc, graph) = TestCase $+    assertBool ("isVertexTransitive " ++ desc) $ isVertexTransitive graph++testlistVertexTransitive = TestList $ map testcaseVertexTransitive [+    -- because we're mapping, these all have to be of same type, hence Graph Int+    ("(q 3)", q 3),+    ("(q 4)", q 4),+    ("petersen", G.to1n petersen)+    ]++testcaseNotVertexTransitive (desc, graph) = TestCase $+    assertBool ("not isVertexTransitive " ++ desc) $ (not . isVertexTransitive) graph++testlistNotVertexTransitive = TestList $ map testcaseNotVertexTransitive [+    ("(kb 2 3)", G.to1n $ kb 2 3),+    ("(kb 3 4)", G.to1n $ kb 3 4),+    ("regular not vertex transitive" , G [(1::Int)..8] [[1,2],[1,3],[1,8],[2,3],[2,4],[3,5],[4,5],[4,6],[5,7],[6,7],[6,8],[7,8]])+    ]+++testcaseEdgeTransitive (desc, graph) = TestCase $+    assertBool ("isEdgeTransitive " ++ desc) $ isEdgeTransitive graph++testlistEdgeTransitive = TestList $ map testcaseEdgeTransitive [+    ("(kb 2 3)", kb 2 3),+    ("(kb 3 4)", kb 3 4)+    ]++testcaseNotEdgeTransitive (desc, graph) = TestCase $+    assertBool ("not isEdgeTransitive " ++ desc) $ (not . isEdgeTransitive) graph++testlistNotEdgeTransitive = TestList $ map testcaseNotEdgeTransitive [+    ("pyramid 4", pyramid 4),+    ("pyramid 5", pyramid 5),+    ("prism 3", G.to1n $ prism 3),+    ("prism 5", G.to1n $ prism 5)+    ]+    where pyramid n = let G vs es = c n in graph (0:vs, [[0,v] | v <- vs] ++ es)+++testcaseArcTransitive (desc, graph) = TestCase $+    assertBool ("isArcTransitive " ++ desc) $ isArcTransitive graph++testlistArcTransitive = TestList $ map testcaseArcTransitive [+    -- Godsil and Royle, p60 - j v k i is arc-transitive+    ("(j 4 2 0)", j 4 2 0),+    ("(j 5 2 0)", j 5 2 0),+    ("(j 5 2 1)", j 5 2 1)+    ]++testcaseNotArcTransitive (desc, graph) = TestCase $+    assertBool ("not isArcTransitive " ++ desc) $ (not . isArcTransitive) graph++testlistNotArcTransitive = TestList $ map testcaseNotArcTransitive [+    ("kb 3 2", kb 3 2),+    ("kb 4 3", kb 4 3)+    ]+++testcase2ArcTransitive (desc, graph) = TestCase $+    assertBool ("is2ArcTransitive " ++ desc) $ is2ArcTransitive graph++testlist2ArcTransitive = TestList $ map testcase2ArcTransitive [+    -- Godsil and Royle, p60 - j (2k+1) k 0 is 2-arc-transitive+    ("(j 3 1 0)", j 3 1 0),+    ("(j 5 2 0)", j 5 2 0),+    ("(j 7 3 0)", j 7 3 0)+    ]++testcaseNot2ArcTransitive (desc, graph) = TestCase $+    assertBool ("not is2ArcTransitive " ++ desc) $ (not . is2ArcTransitive) graph++testlistNot2ArcTransitive = TestList $ map testcaseNot2ArcTransitive [+    -- because a 2-arc can be two sides of a triangle, or not, so they are not all alike+    ("octahedron", octahedron),+    ("icosahedron", icosahedron)+    ]+++testcase4ArcTransitive (desc, graph) = TestCase $+    assertBool ("is4ArcTransitive " ++ desc) $ is4ArcTransitive graph++testlist4ArcTransitive = TestList $ map testcase4ArcTransitive [+    -- Godsil and Royle, p80-1+    ("PG(2,F2)", G.to1n $ incidenceGraphPG 2 f2),+    ("PG(2,F3)", G.to1n $ incidenceGraphPG 2 f3),+    ("PG(2,F4)", G.to1n $ incidenceGraphPG 2 f4)+    ]+++testcaseDistanceTransitive (desc, graph) = TestCase $+    assertBool ("isDistanceTransitive " ++ desc) $ isDistanceTransitive graph++testlistDistanceTransitive = TestList $ map testcaseDistanceTransitive [+    ("(kb 3 3)", G.to1n $ kb 3 3),+    ("(kb 4 4)", G.to1n $ kb 4 4),+    ("(q 3)", G.to1n $ q 3),+    ("(q 4)", G.to1n $ q 4),+    ("petersen", G.to1n $ petersen),+    -- Godsil and Royle, p67 - j v k (k-1) and j (2k+1) (k+1) 0 are distance-transitive+    ("(j 3 2 1)", G.to1n $ j 3 2 1),+    ("(j 4 2 1)", G.to1n $ j 4 2 1),+    ("(j 5 3 2)", G.to1n $ j 5 3 2)+    -- ("(j 3 2 0)", G.to1n $ j 3 2 0),+    -- ("(j 5 3 0)", G.to1n $ j 5 3 0),+    -- ("(j 7 4 0)", G.to1n $ j 7 4 0)+    ]++testcaseNotDistanceTransitive (desc, graph) = TestCase $+    assertBool ("not isDistanceTransitive " ++ desc) $ (not . isDistanceTransitive) graph++testlistNotDistanceTransitive = TestList $ map testcaseNotDistanceTransitive [+    ("(prism 3)", prism 3),+    -- not prism 4, which is the cube+    ("(prism 5)", prism 5)+    ]+

@@ -11,6 +11,7 @@   testlistPrimeFactor = TestList [+    testcasePrimesList,     testlistSmallPrimes,     testlistMillerRabin,     testlistMersennePrimes,@@ -30,11 +31,14 @@     testlistFactorOrder     ] +primesTo100 = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97] +testcasePrimesList = TestCase $ assertEqual "primes list"+    primesTo100 (takeWhile (<100) primes)+ testlistSmallPrimes = TestList [     TestCase $ assertEqual "small primes"-               [False,True,True,False,True,False,True,False,False,False]-               (map isPrime [1..10]),+               primesTo100 (filter isPrime [1..100]),     TestCase $ assertBool "negative primes" (all notPrime [-10..0])     ]