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

module Math.Algebra.Group.StringRewriting where

import Data.List as L
import Data.Maybe (catMaybes)

-- REWRITING

-- |Given a list of rewrite rules of the form (left,right), and a word,
-- rewrite it by repeatedly replacing any left substring in the word by the corresponding right
rewrite :: (Eq a) => [([a], [a])] -> [a] -> [a]
rewrite rules word = rewrite' rules word where
rewrite' (r:rs) xs =
case rewrite1 r xs of
Nothing -> rewrite' rs xs
Just ys -> rewrite' rules ys
rewrite' [] xs = xs

rewrite1 (l,r) xs =
case xs `splitSubstring` l of
Nothing -> Nothing
Just (a,b) -> Just (a++r++b)

-- given a string x and a substring b, find if possible (a,c) such that xs = abc
splitSubstring xs ys = splitSubstring' [] xs where
splitSubstring' ls [] = Nothing
splitSubstring' ls (r:rs) =
if ys `L.isPrefixOf` (r:rs)
then Just (reverse ls, drop (length ys) (r:rs))
else splitSubstring' (r:ls) rs
-- there might be a more efficient way to do this

-- KNUTH-BENDIX

-- given two strings x,y, find if possible a,b,c with x=ab y=bc
findOverlap xs ys = findOverlap' [] xs ys where
findOverlap' as [] cs = Nothing -- (reverse as, [], cs)
findOverlap' as (b:bs) cs =
if (b:bs) `L.isPrefixOf` cs
then Just (reverse as, b:bs, drop (length (b:bs)) cs)
else findOverlap' (b:as) bs cs
-- there might be a more efficient way to do this

-- note that findOverlap "abab" "abab" won't find the partial overlap ("ab","ab","ab")

-- Knuth-Bendix algorithm
-- http://en.wikipedia.org/wiki/Knuth-Bendix_algorithm
-- Given a set of rules (assumed already reduced with respect to each other)
-- return a confluent rewrite system
knuthBendix1 rules = knuthBendix' rules pairs where
pairs = [(lri,lrj) | lri <- rules, lrj <- rules, lri /= lrj]
knuthBendix' rules [] = rules -- should reduce in some way
knuthBendix' rules ( ((li,ri),(lj,rj)) : ps) =
case findOverlap li lj of
Nothing -> knuthBendix' rules ps
Just (a,b,c) -> case ordpair (rewrite rules (ri++c)) (rewrite rules (a++rj)) of
Nothing -> knuthBendix' rules ps -- they both reduce to the same thing
Just rule' -> let rules' = reduce rule' rules
ps' = ps ++ [(rule',rule) | rule <- rules'] ++ [(rule,rule') | rule <- rules']
in knuthBendix' (rule':rules') ps'
-- the new rule comes from seeing that
-- a ++ b ++ c == l1 ++ c -> r1 ++ c (by rule 1)
-- a ++ b ++ c == a ++ l2 -> a ++ r2 (by rule 2)
reduce rule@(l,r) rules = filter (\(l',r') -> not (L.isInfixOf l l')) rules
-- [rule' | rule'@(l',r') <- rules, not (l `L.isInfixOf` l')]

ordpair x y =
case shortlex x y of
LT -> Just (y,x)
EQ -> Nothing
GT -> Just (x,y)

shortlex x y = compare (length x, x) (length y, y)

-- for groups, where "letters" will take the form Either a a, we will want a different order, because we will want x^-1 -> x^3 to be the right way round

-- An optimisation - keep the rules ordered smallest first, and process the pairs smallest first
-- Appears to be significantly faster on average
knuthBendix2 rules = map snd \$ knuthBendix' rules' pairs where
rules' = L.sort \$ map sizedRule rules
pairs = L.sort [sizedPair sri srj | sri <- rules', srj <- rules', sri /= srj]
knuthBendix' rules [] = rules
knuthBendix' rules ( (s,(li,ri),(lj,rj)) : ps) =
case findOverlap li lj of
Nothing -> knuthBendix' rules ps
Just (a,b,c) -> case ordpair (rewrite (map snd rules) (ri++c)) (rewrite (map snd rules) (a++rj)) of
Nothing -> knuthBendix' rules ps -- they both reduce to the same thing
Just rule' -> let rules' = reduce (snd rule') rules
-- ps' = L.sort \$ ps ++ [sizedPair rule' rule | rule <- rules'] ++ [sizedPair rule rule' | rule <- rules']
ps' = merge ps \$ merge [sizedPair rule' rule | rule <- rules'] [sizedPair rule rule' | rule <- rules']
in knuthBendix' (L.insert rule' rules') ps'
reduce rule@(l,r) rules = filter (\(s',(l',r')) -> not (L.isInfixOf l l')) rules
-- reduce rule@(l,r) rules = [rule' | rule'@(s',(l',r')) <- rules, not (l `L.isInfixOf` l')]
ordpair x y =
let lx = length x; ly = length y in
case compare (lx,x) (ly,y) of
LT -> Just (ly,(y,x)); EQ -> Nothing; GT -> Just (lx,(x,y))
sizedRule (rule@(l,r)) = (length l, rule)
sizedPair (s1,r1) (s2,r2) = (s1+s2,r1,r2)

-- merge two ordered lists
merge (x:xs) (y:ys) =
case compare x y of
LT -> x : merge xs (y:ys)
GT -> y : merge (x:xs) ys
EQ -> error "" -- shouldn't happen in our case
merge xs ys = xs++ys

-- Another optimisation - at the stage where we remove some rules, we remove corresponding pairs too
-- Seems to perform about 25% faster on large problems (eg Coxeter groups A4-12, B4-12)
knuthBendix3 rules = knuthBendix' rules' pairs (length rules' + 1) where
rules' = L.sort \$ zipWith (\i (l,r) -> (length l,i,(l,r)) ) [1..] rules
pairs = L.sort [sizedPair ri rj | ri <- rules', rj <- rules', ri /= rj]
knuthBendix' rules [] k = map (\(s,i,r) -> r) rules
knuthBendix' rules ( (s,(i,j),((li,ri),(lj,rj))) : ps) k =
case findOverlap li lj of
Nothing -> knuthBendix' rules ps k
Just (a,b,c) -> case ordpair k (rewrite (map third rules) (ri++c)) (rewrite (map third rules) (a++rj)) of
Nothing -> knuthBendix' rules ps k -- they both reduce to the same thing
Just rule'@(_,_,(l,r)) ->
let (outrules,inrules) = L.partition (\(s',i',(l',r')) -> L.isInfixOf l l') rules
removedIndices = map second outrules
ps' = [p | p@(s,(i,j),(ri,rj)) <- ps, i `notElem` removedIndices, j `notElem` removedIndices]
ps'' = merge ps' \$ merge [sizedPair rule' rule | rule <- inrules] [sizedPair rule rule' | rule <- inrules]
in knuthBendix' (L.insert rule' inrules) ps'' (k+1)
ordpair k x y =
let lx = length x; ly = length y in
case compare (lx,x) (ly,y) of
LT -> Just (ly,k,(y,x)); EQ -> Nothing; GT -> Just (lx,k,(x,y))
second (s,i,r) = i
third (s,i,r) = r
sizedPair (si,i,ri) (sj,j,rj) = (si+sj,(i,j),(ri,rj))

-- |Implementation of the Knuth-Bendix algorithm. Given a list of relations, return a confluent rewrite system.
-- The algorithm is not guaranteed to terminate.
knuthBendix :: (Ord a) => [([a], [a])] -> [([a], [a])]
knuthBendix relations = knuthBendix3 (reduce [] rules) where
rules = catMaybes [ordpair x y | (x,y) <- relations]
reduce ls (r:rs) = reduce (r: reduce' r ls) (reduce' r rs)
reduce ls [] = ls
reduce' r rules = catMaybes [ordpair (rewrite [r] lhs) (rewrite [r] rhs) | (lhs,rhs) <- rules]

-- |Given generators and a confluent rewrite system, return (normal forms of) all elements
nfs :: (Ord a) => ([a], [([a], [a])]) -> [[a]]
nfs (gs,rs) = nfs' [[]] where
nfs' [] = [] -- we have run out of words - this monoid is finite
nfs' ws = let ws' = [g:w | g <- gs, w <- ws, not (any (`L.isPrefixOf` (g:w)) (map fst rs))]
in ws ++ nfs' ws'

-- |Given generators and relations, return (normal forms of) all elements
elts :: (Ord a) => ([a], [([a], [a])]) -> [[a]]
elts (gs,rs) = nfs (gs, knuthBendix rs)

-- PRESENTATIONS FOR SOME STANDARD GROUPS
-- Would like to add a few more to this list

newtype SGen = S Int deriving (Eq,Ord)

instance Show SGen where
show (S i) = "s" ++ show i

s_ i = S i
s1 = s_ 1
s2 = s_ 2
s3 = s_ 3

-- D L Johnson, Presentations of Groups, p62

-- symmetric group, generated by adjacent transpositions
_S n = (gs, r ++ s ++ t) where
gs = map s_ [1..n-1]
r = [([s_ i, s_ i],[]) | i <- [1..n-1]]
s = [(concat \$ replicate 3 [s_ i, s_ (i+1)],[]) | i <- [1..n-2]]
t = [([s_ i, s_ j, s_ i, s_ j],[]) | i <- [1..n-1], j <- [i+2..n-1]]

-- braid presentation for Sn
_S' n = (gs, r ++ s ++ t) where
gs = map s_ [1..n-1]
r = [([s_ i, s_ i], []) | i <- [1..n-1]]
s = [([s_ (i+1), s_ i, s_ (i+1)], [s_ i, s_ (i+1), s_ i] ) | i <- [1..n-2]]
t = [([s_ i, s_ j, s_ i, s_ j], []) | i <- [1..n-1], j <- [i+2..n-1]]

-- http://en.wikipedia.org/wiki/Triangle_group
-- triangle groups - Johnson p127ff
tri l m n = ("abc", [("aa",""),("bb",""),("cc",""),("ab" ^ l,""),("bc" ^ n,""),("ca" ^ m,"" )])
where xs ^ i = concat \$ replicate i xs

-- von Dyck groups - Johnson p121ff
-- The subgroup of index 2 in the triangle group consisting of elts that preserve the orientation of the triangle
_D l m n = ("xy", [("x" ^ l,""), ("y" ^ m,""), ("xy" ^ n,"")])
where xs ^ i = concat \$ replicate i xs

-- Degenerate cases: n == 1 => cyclic group
-- l,2,2, l>=2 -> n-gon bipyramid - dihedral group
-- Spherical case: 1/l+1/m+1/n > 1
-- 3,3,2 -> tetrahedron; 4,3,2 -> octahedron; 5,3,2 -> icosahedron
-- Euclidean case: 1/l+1/m+1/n == 1
-- 3,3,3  4,4,2  6,3,2
-- Hyperbolic case: 1/l+1/m+1/n < 1

```