```-- Copyright (c) David Amos, 2008. All rights reserved.

module Math.Algebra.Group.SchreierSims where

import qualified Data.List as L
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.Common.ListSet (toListSet)

-- COSET REPRESENTATIVES FOR STABILISER OF A POINT

-- Given a group G = <gs>, and a point x, find coset representatives for Gx (stabiliser of x) in G
-- In other words, for each x' in the orbit of x under G, we find a g <- G taking x to x'
-- The code is similar to the code for calculating orbits, but modified to keep track of the group elements that we used to get there
cosetRepsGx gs x = cosetRepsGx' gs M.empty (M.singleton x 1) where
cosetRepsGx' gs interior boundary
| M.null boundary = interior
| otherwise =
let interior' = M.union interior boundary
boundary' = M.fromList [(p .^ g, h*g) | g <- gs, (p,h) <- M.toList boundary] M.\\ interior'
in cosetRepsGx' gs interior' boundary'

-- SCHREIER GENERATORS

-- toSet xs = (map head . group . sort) xs

-- Generators for Gx, the stabiliser of x, given that G is generated by gs, and rs is a set of coset representatives for Gx in G.
-- Schreier's Lemma states that if H < G = <S>, and R is a set of coset reps for H in G
-- then H is generated by { rs(rs)*^-1 | r <- R, s <- S } (where * means "the coset representative of")
-- In particular, with H = Gx, this gives us a way of finding a set of generators for Gx
schreierGeneratorsGx (x,rs) gs = L.nub \$ filter (/= 1) [schreierGenerator r g | r <- M.elems rs,  g <- gs]
where schreierGenerator r g = let h = r * g
h' = rs M.! (x .^ h)
in h * inverse h'

-- SCHREIER-SIMS ALGORITHM

sift bts g = sift' bts g where
sift' [] g = if g == 1 then Nothing else Just g
sift' ((b,t):bts) g = case M.lookup (b .^ g) t of
Nothing -> Just g
Just h -> sift' bts (g * inverse h)

findBase gs = minimum \$ concatMap supp gs
{-
-- Find base and strong generating set using Schreier-Sims algorithm
bsgs gs | all (/= 1) gs = map fst \$ ss [newLevel gs] []

newLevel s = let b = findBase s
t = cosetRepsGx s b
in ((b,t),s)

let bts = map fst goods
sgs = schreierGeneratorsGx (b,t) s
siftees = filter isJust \$ map (sift bts) sgs
in  if null siftees
else let Just h = head siftees
in if null goods
then ss (newLevel [h] : bad : bads) []
else let ((b_,t_),s_) = head goods
s' = h:s_
t' = cosetRepsGx s' b_
in ss (((b_,t'),s') : bad : bads) (tail goods)
ss [] goods = goods
-}

-- strong generating set, with implied base from the Ord instance
sgs gs = toListSet \$ concatMap snd \$ ss bs gs
where bs = toListSet \$ concatMap supp gs

-- Find base and strong generating set using Schreier-Sims algorithm
-- !! This function is poorly named - it actually finds you a base and sets of transversals
-- This version guarantees to use bases in order
bsgs gs = bsgs' bs gs
where bs = toListSet \$ concatMap supp gs

-- This version lets you pass in bases in the order you want them (or [], and it will find its own)
bsgs' bs gs = map fst \$ ss bs gs

-- For example, bsgs (_A 5) uses [1,2,3] as the bases, but bsgs' [] (_A 5) uses [1,3,2]

newLevel (b:bs) s = (bs, newLevel' b s)
newLevel [] s = ([], newLevel' b s) where b = findBase s
newLevel' b s = ((b,t),s) where t = cosetRepsGx s b

ss bs gs = ss' bs' [level] []
where (bs',level) = newLevel bs \$ filter (/=1) gs

let bts = map fst goods
sgs = schreierGeneratorsGx (b,t) s
siftees = filter isJust \$ map (sift bts) sgs
in  if null siftees
else let Just h = head siftees
in if null goods
then let (bs', level) = newLevel bs [h]
in ss' bs' (level : bad : bads) []
else let ((b_,t_),s_) = head goods
s' = h:s_
t' = cosetRepsGx s' b_
in ss' bs (((b_,t'),s') : bad : bads) (tail goods)
ss' _ [] goods = goods
{-
extendbsgs [] g = bsgs [g]
extendbsgs (((b,t),s):bts) g = ss (((b,t),g:s):bts) []

bsgs' gs = map fst \$ foldl extendbsgs [] gs
-}

-- The above is written for simplicity.
-- Its efficiency could be improved by incrementally updating the transversals,
-- and keeping track of Schreier generators we have already tried.
-- (Remember to add new Schreier generators every time the generating set or transversal is augmented.)

-- USING THE SCHREIER-SIMS TRANSVERSALS

isMemberBSGS bts g = isNothing \$ sift bts g

-- By Lagrange's thm, every g <- G can be written uniquely as g = r_m ... r_1 (Seress p56)
-- Note that we have to reverse the list of coset representatives
eltsBSGS bts = map (product . reverse) (cartProd ts)
where ts = map (M.elems . snd) bts

cartProd (set:sets) = [x:xs | x <- set, xs <- cartProd sets]
cartProd [] = [[]]

orderBSGS bts = product (map (toInteger . M.size . snd) bts)

-- |Given generators for a group, determine whether a permutation is a member of the group, using Schreier-Sims algorithm
isMember :: (Ord t, Show t) => [Permutation t] -> Permutation t -> Bool
isMember gs h = isMemberBSGS (bsgs gs) h

-- |Given generators for a group, return a (sorted) list of all elements of the group, using Schreier-Sims algorithm
elts :: (Ord t, Show t) => [Permutation t] -> [Permutation t]
elts [] = 
elts gs = eltsBSGS \$ bsgs gs

-- |Given generators for a group, return the order of the group (the number of elements), using Schreier-Sims algorithm
order :: (Ord t, Show t) => [Permutation t] -> Integer
order [] = 1
order gs = orderBSGS \$ bsgs gs

isSubgp hs gs = all (isMemberBSGS gs') hs
where gs' = bsgs gs

isNormal hs gs = hs `isSubgp` gs
&& all (isMemberBSGS hs') [h~^g | h <- hs, g <- gs]
where hs' = bsgs hs

index gs hs = order gs `div` order hs

-- given list of generators, try to find a shorter list
reduceGens gs = fst \$ reduceGensBSGS (filter (/=1) gs)

reduceGensBSGS (g:gs) = reduceGens' ([g],bsgs [g]) gs where
reduceGens' (gs,bsgsgs) (h:hs) =
if isMemberBSGS bsgsgs h
then reduceGens' (gs,bsgsgs) hs
else reduceGens' (h:gs, bsgs \$ h:gs) hs
reduceGens' (gs,bsgsgs) [] = (reverse gs,bsgsgs)
reduceGensBSGS [] = ([],[])

-- normal closure of H in G
-- for efficiency, should be called with gs and hs already reduced sets of generators
normalClosure gs hs = reduceGens \$ hs ++ [h ~^ g | h <- hs, g <- gs']
where gs' = gs ++ map inverse gs

-- commutator gp of H and K
commutatorGp hs ks = normalClosure (hsks) [h^-1 * k^-1 * h * k | h <- hs', k <- ks']
where hs' = reduceGens hs
ks' = reduceGens ks
hsks = reduceGens (hs' ++ ks')
-- no point processing more potential generators than we have to

-- derived subgroup (or commutator subgroup)
derivedSubgp gs = normalClosure gs' [g^-1 * h^-1 * g * h | g <- gs', h <- gs']
where gs' = reduceGens gs
-- == commutatorGp gs gs

{-
isPerfect gs = order gs == order (derivedSubgp gs)
-- We compare orders rather than the generators themselves, because derivedSubgp will usually find different generators

derivedSeries gs = derivedSeries' (gs, order gs) where
derivedSeries' ([],1) = [[]]
derivedSeries' (hs, orderhs) =
let hs' = derivedSubgp hs
orderhs' = order hs'
in if orderhs' == orderhs
then [hs]
else hs : derivedSeries' (hs',orderhs')

lowerCentralSeries gs = lowerCentralSeries' (gs, order gs) where
lowerCentralSeries' ([],1) = [[]]
lowerCentralSeries' (hs, orderhs) =
let hs' = commutatorGp gs hs
orderhs' = order hs'
in if orderhs' == orderhs
then [hs]
else hs : lowerCentralSeries' (hs',orderhs')
-}

```