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

    module Math.Projects.RootSystem where

    import Data.Ratio
    import qualified Data.List as L
    import qualified Data.Set as S

    import Math.Algebra.LinearAlgebra
    import Math.Algebra.Group.PermutationGroup hiding (elts, order, closure)
    import Math.Algebra.Group.SchreierSims as SS
    import Math.Algebra.Group.StringRewriting as SG

    import Math.Algebra.Field.Base -- for Q

    data Type = A | B | C | D | E | F | G

    -- Humphreys, Reflection Groups and Coxeter Groups

    -- SIMPLE SYSTEMS
    -- sometimes called fundamental systems

    -- The ith basis vector in K^n
    basisElt :: Int -> Int -> [Q] -- this type signature determines all the rest
    basisElt n i = replicate (i-1) 0 ++ 1 : replicate (n-i) 0
    -- We need to work over the rationals to ensure that arithmetic is exact
    -- So long as our simple systems are rational, then reflection matrices are rational

    -- A simple system is like a basis for the root system (see Humphreys p8 for full definition)
    -- simpleSystem :: Type -> Int -> [[Q]]
    simpleSystem A n | n >= 1 = [e i <-> e (i+1) | i <- [1..n]]
        where e = basisElt (n+1)
    simpleSystem B n | n >= 2 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [e n]
        where e = basisElt n
    simpleSystem C n | n >= 2 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [2 *> e n]
        where e = basisElt n
    simpleSystem D n | n >= 4 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [e (n-1) <+> e n]
        where e = basisElt n
    simpleSystem E n | n `elem` [6,7,8] = take n simpleroots
        where e = basisElt 8
              simpleroots = ((1/2) *> (e 1 <-> e 2 <-> e 3 <-> e 4 <-> e 5 <-> e 6 <-> e 7 <+> e 8))
                          : (e 1 <+> e 2)
                          : [e (i-1) <-> e (i-2) | i <- [3..8]]
    simpleSystem F 4 = [e 2 <-> e 3, e 3 <-> e 4, e 4, (1/2) *> (e 1 <-> e 2 <-> e 3 <-> e 4)]
        where e = basisElt 4
    simpleSystem G 2 = [e 1 <-> e 2, ((-2) *> e 1) <+> e 2 <+> e 3]
        where e = basisElt 3

    -- ROOT SYSTEMS
    -- Calculating the full root system from the fundamental roots

    -- Humphreys p3
    -- Weyl group element corresponding to a root
    -- w r is the reflection in the hyperplane orthogonal to r
    w r s = s <-> (2 * (s <.> r) / (r <.> r)) *> r

    -- Given a simple system, return the full root system
    -- The closure of a set of roots under reflection
    closure rs = S.toList $ closure' S.empty (S.fromList rs) where
        closure' interior boundary
            | S.null boundary = interior
            | otherwise =
                let interior' = S.union interior boundary
                    boundary' = S.fromList [w r s | r <- rs, s <- S.toList boundary] S.\\ interior'
                in closure' interior' boundary'

    -- WEYL GROUP
    -- The finite reflection group generated by the root system

    -- Generators of the Weyl group as permutation group on the roots
    weylPerms t n =
        let rs = simpleSystem t n
            xs = closure rs
            toPerm r = fromPairs [(x, w r x) | x <- xs]
        in map toPerm rs

    -- Generators of the Weyl group as a matrix group
    weylMatrices t n = map wMx (simpleSystem t n)

    -- The Weyl group element corresponding to a root, represented as a matrix
    wMx r = map (w r) [e i | i <- [1..d]] -- matrix for reflection in hyperplane orthogonal to r
        where d = length r -- dimension of the space
              e = basisElt d
    -- the images of the basis elts form the columns of the matrix
    -- however, reflection matrices are symmetric, so they also form the rows

    -- CARTAN MATRIX, DYNKIN DIAGRAM, COXETER SYSTEM

    cartanMatrix t n = [[2 * (ai <.> aj) / (ai <.> ai) | aj <- roots] | ai <- roots]
        where roots = simpleSystem t n
    -- Note: The Cartan matrices for A, D, E systems are symmetric.
    -- Those of B, C, F, G are not
    -- Carter, Simple Groups of Lie Type, p44-5 gives the expected answers
    -- They agree with our answers except for G2, which is the transpose
    -- (So probably Carter defines the roots of G2 the other way round to Humphreys)

    -- set the diagonal entries of (square) matrix mx to constant c
    setDiag c mx@((x:xs):rs) = (c:xs) : zipWith (:) (map head rs) (setDiag c $ map tail rs)
    setDiag _ [[]] = [[]]

    -- Carter, Segal, Macdonald p17-18
    -- given a Cartan matrix, derive the corresponding matrix describing the Dynkin diagram
    -- nij = Aij * Aji, nii = 0
    dynkinFromCartan aij = setDiag 0 $ (zipWith . zipWith) (*) aij (L.transpose aij)

    dynkinDiagram t n = dynkinFromCartan $ cartanMatrix t n

    -- given the Dynkin diagram nij, derive the coefficients mij of the Coxeter group <si | si^2, (sisj)^mij> (so mii == 1)
    -- using nij = 4 cos^2 theta_ij
    -- nij == 0 <=> theta = pi/2
    -- nij == 1 <=> theta = pi/3
    -- nij == 2 <=> theta = pi/4
    -- nij == 3 <=> theta = pi/6
    coxeterFromDynkin nij = setDiag 1 $ (map . map) f nij
        where f 0 = 2; f 1 = 3; f 2 = 4; f 3 = 6

    -- The mij coefficients of the Coxeter group <si | si^2, (sisj)^mij>, as a matrix
    coxeterMatrix t n = coxeterFromDynkin $ dynkinDiagram t n

    -- Given the matrix of coefficients mij, return the Coxeter group <si | si^2, (sisj)^mij>
    -- We assume but don't check that mii == 1 and mij == mji
    fromCoxeterMatrix mx = (gs,rs) where
        n = length mx
        gs = map s_ [1..n]
        rs = rules mx 1
        rules [] _ = []
        rules ((1:xs):rs) i = ([s_ i, s_ i],[]) : [powerRelation i j m | (j,m) <- zip [i+1..] xs] ++ rules (map tail rs) (i+1)
        powerRelation i j m = (concat $ replicate m [s_ i, s_ j],[])

    -- Another presentation for the Coxeter group, using braid relations
    fromCoxeterMatrix2 mx = (gs,rs) where
        n = length mx
        gs = map s_ [1..n]
        rs = rules mx 1
        rules [] _ = []
        rules ((1:xs):rs) i = ([s_ i, s_ i],[]) : [braidRelation i j m | (j,m) <- zip [i+1..] xs] ++ rules (map tail rs) (i+1)
        braidRelation i j m = (take m $ cycle [s_ j, s_ i], take m $ cycle [s_ i, s_ j])

    coxeterPresentation t n = fromCoxeterMatrix $ coxeterMatrix t n

    eltsCoxeter t n = SG.elts $ fromCoxeterMatrix2 $ coxeterMatrix t n
    -- it's just slightly faster to use the braid presentation

    poincarePoly t n = map length $ L.group $ map length $ eltsCoxeter t n

    -- LIE ALGEBRAS

    elemMx n i j = replicate (i-1) z ++ e j : replicate (n-i) z
        where z = replicate n 0
              e = basisElt n

    lieMult x y = x*y - y*x

    -- for gluing matrices together
    (+|+) = zipWith (++) -- glue two matrices together side by side
    (+-+) = (++)         -- glue two matrices together above and below

    form D n = (zMx n +|+ idMx n)
                      +-+
              (idMx n +|+ zMx n)
    form C n = (2 : replicate (2*n) 0) :
               (map (0:) (form D n))
    form B n = let id' = (-1) *>> idMx n
               in (zMx n +|+ idMx n)
                         +-+
                   (id'  +|+ zMx n)

    -- TESTING
    -- The expected values of the root system, number of roots, order of Weyl group
    -- for comparison against the calculated values

    -- !! Not yet got root systems for E6,7,8, F4

    -- Humphreys p41ff

    -- The full root system
    -- L.sort (rootSystem t n) == closure (simpleSystem t n)
    -- rootSystem :: Type -> Int -> [[QQ]]
    rootSystem A n | n >= 1 = [e i <-> e j | i <- [1..n+1], j <- [1..n+1], i /= j]
        where e = basisElt (n+1)
    rootSystem B n | n >= 2 = shortRoots ++ longRoots
        where e = basisElt n
              shortRoots = [e i | i <- [1..n]]
                        ++ [[] <-> e i | i <- [1..n]]
              longRoots  = [e i <+> e j | i <- [1..n], j <- [i+1..n]]
                        ++ [e i <-> e j | i <- [1..n], j <- [i+1..n]]
                        ++ [[] <-> e i <+> e j | i <- [1..n], j <- [i+1..n]]
                        ++ [[] <-> e i <-> e j | i <- [1..n], j <- [i+1..n]]
    rootSystem C n | n >= 2 = longRoots ++ shortRoots
        where e = basisElt n
              longRoots  = [2 *> e i | i <- [1..n]]
                        ++ [[] <-> (2 *> e i) | i <- [1..n]]
              shortRoots = [e i <+> e j | i <- [1..n], j <- [i+1..n]]
                        ++ [e i <-> e j | i <- [1..n], j <- [i+1..n]]
                        ++ [[] <-> e i <+> e j | i <- [1..n], j <- [i+1..n]]
                        ++ [[] <-> e i <-> e j | i <- [1..n], j <- [i+1..n]]
    rootSystem D n | n >= 4 =
        [e i <+> e j | i <- [1..n], j <- [i+1..n]]
        ++ [e i <-> e j | i <- [1..n], j <- [i+1..n]]
        ++ [[] <-> e i <+> e j | i <- [1..n], j <- [i+1..n]]
        ++ [[] <-> e i <-> e j | i <- [1..n], j <- [i+1..n]]
        where e = basisElt n
    rootSystem G 2 = shortRoots ++ longRoots
        where e = basisElt 3
              shortRoots = [e i <-> e j | i <- [1..3], j <- [1..3], i /= j]
              longRoots = concatMap (\r-> [r,[] <-> r]) [2 *> e i <-> e j <-> e k | i <- [1..3], [j,k] <- [[1..3] L.\\ [i]] ]

    -- numRoots t n == length (closure $ simpleSystem t n)
    numRoots A n = n*(n+1)
    numRoots B n = 2*n*n
    numRoots C n = 2*n*n
    numRoots D n = 2*n*(n-1)
    numRoots E 6 = 72
    numRoots E 7 = 126
    numRoots E 8 = 240
    numRoots F 4 = 48
    numRoots G 2 = 12

    -- The order of the Weyl group
    -- orderWeyl t n == S.order (weylPerms t n)
    orderWeyl A n = factorial (n+1)
    orderWeyl B n = 2^n * factorial n
    orderWeyl C n = 2^n * factorial n
    orderWeyl D n = 2^(n-1) * factorial n
    orderWeyl E 6 = 2^7 * 3^4 * 5
    orderWeyl E 7 = 2^10 * 3^4 * 5 * 7
    orderWeyl E 8 = 2^14 * 3^5 * 5^2 * 7
    orderWeyl F 4 = 2^7 * 3^2
    orderWeyl G 2 = 12

    factorial n = product [1..toInteger n]

    {-
    -- now moved to TRootSystem
    test1 = all (\(t,n) -> orderWeyl t n == L.genericLength (eltsCoxeter t n))
        [(A,3),(A,4),(A,5),(B,3),(B,4),(B,5),(C,3),(C,4),(C,5),(D,4),(D,5),(F,4),(G,2)]

    test2 = all (\(t,n) -> orderWeyl t n == SS.order (weylPerms t n))
        [(A,3),(A,4),(A,5),(B,3),(B,4),(B,5),(C,3),(C,4),(C,5),(D,4),(D,5),(E,6),(F,4),(G,2)]
    -}