Permalink
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
389 lines (302 sloc) 14.2 KB
-- Copyright (c) Yuchen Pei, 2017. (Added positive roots and highest elements etc.)
-- Copyright (c) David Amos, 2008-2015. All rights reserved.
--module Math.Combinatorics.RootSystem where
module RootSystem where
import Prelude hiding ( (*>) )
import Data.Ratio
import Data.List
import Data.Maybe
import qualified Data.Set as S
import Math.Algebra.LinearAlgebra
type Q = Rational
data Type = A | B | C | D | E | F | G deriving Show
type SimpleSystem = [[Q]]
-- 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
--basisElt' :: Int -> Int -> [Int] -- 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 -> SimpleSystem
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 B 1 = simpleSystem A 1
simpleSystem C n | n >= 2 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [2 *> e n]
where e = basisElt n
simpleSystem C 1 = simpleSystem A 1
simpleSystem D n | n >= 4 = [e i <-> e (i+1) | i <- [1..n-1]] ++ [e (n-1) <+> e n]
where e = basisElt n
simpleSystem D 3 = simpleSystem A 3
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
simpleSystem t n = error $ "Invalid root system of type " ++ (show t) ++ " and rank " ++ (show n) ++ "."
dimensionOfHostSpace :: Type -> Int -> Int
dimensionOfHostSpace t n = length $ head $ simpleSystem t n
-- Humphreys p3
-- Reflection corresponding to a root
-- s alpha beta = s_\alpha \beta
s :: [Q] -> [Q] -> [Q]
s alpha beta = beta <-> (dynkinIndex alpha beta) *> alpha
-- |The indices of the simple roots in the reduced decomposition of the longest elements
longestElementIndex :: SimpleSystem -> [Int]
longestElementIndex ss = (+1) <$> fromJust <$> flip elemIndex ss <$> longestElement ss
-- |The reduced decomposition of the longest elements
longestElement :: SimpleSystem -> [[Q]]
longestElement ss = longestElement' [] posRoots
where
posRoots = positiveRoots ss
longestElement' :: [[Q]] -> [[Q]] -> [[Q]]
longestElement' xs rs =
let ys = (S.fromList ss) `S.intersection` (S.fromList rs) in
if S.null ys
then xs
else let alpha = (head (S.toList ys)) in
longestElement' (alpha:xs) (s alpha <$> rs)
-- |all positive roots
positiveRoots :: SimpleSystem -> [[Q]]
positiveRoots ss = (positiveRoots' $ cartanMatrix' ss) <<*>> ss
-- |return root indices of all positive roots
positiveRoots' :: [[Q]] -> [[Q]]
positiveRoots' cm = positiveRoots'' [] (iMx $ length cm)
where positiveRoots'' :: [[Q]] -> [[Q]] -> [[Q]]
positiveRoots'' pr npr
| null npr = pr
| otherwise = positiveRoots'' (pr ++ npr) (S.toList . S.fromList $ mconcat $ zipWith newRootIndices npr (pIndex pr cm <$> npr))
-- |calculate the successors of ri in the simple string
-- |given root index ri and its p-vector pi
newRootIndices :: [Q] -> [Q] -> [[Q]]
newRootIndices ri pi = (ri <+>) <$> (go pi [] 1)
where go :: [Q] -> [Int] -> Int -> [[Q]]
go [] ys n = basisElt n <$> ys
go (x:xs) ys k = go xs (if x == 0 then ys else ys ++ [k]) (k + 1)
-- |given the already-found positive roots allRootIndices,
-- |the transposed Cartan Matrix cm and the root index, calculate its p-vector
pIndex :: [[Q]] -> [[Q]] -> [Q] -> [Q]
pIndex allRootIndices cm rootIndex = (mIndex rootIndex allRootIndices) <-> (dynkinIndex' rootIndex cm)
-- |given the root index and all already-found positive roots, calculate its m-vector
mIndex :: [Q] -> [[Q]] -> [Q]
mIndex rootIndex allRootIndices =
if isMultBasis rootIndex
then 2 *> rootIndex
else maxV $ filter isMultBasis [rootIndex <-> r | r <- allRootIndices]
maxV :: Ord a => [[a]] -> [a]
maxV xs = maximum <$> (transpose xs)
sumV :: Num a => [[a]] -> [a]
sumV = foldl1 (<+>)
cartanMatrix :: SimpleSystem -> [[Q]]
cartanMatrix ss = [dynkinIndex r <$> ss | r <- ss]
cartanMatrix' :: SimpleSystem -> [[Q]]
cartanMatrix' ss = transpose $ cartanMatrix ss
isMultBasis :: [Q] -> Bool
isMultBasis v = (length $ filter (/=0) v) == 1
dynkinRow :: [Q] -> SimpleSystem -> [Q]
dynkinRow r ss = dynkinIndex r <$> ss
dynkinIndex :: [Q] -> [Q] -> Q
dynkinIndex r s = 2 * (r <.> s) / (r <.> r)
dynkinIndex' :: [Q] -> [[Q]] -> [Q]
dynkinIndex' ri cm = ri <*>> cm
isInWeylChamber :: SimpleSystem -> [Q] -> Bool
isInWeylChamber ss r = all (\t -> t <.> r >= 0) ss
-- 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
-- |Auxilieary function to transform an Int to a Root system for quick check
int2TypeInt :: Int -> (Type, Int)
int2TypeInt n
| n > 108 = int2TypeInt (n `mod` 108 + 1)
| n > 105 = (E, n - 100)
| n > 96 = (G, 2)
| n > 88 = (F, 4)
| n > 85 = (E, n - 80)
| n > 62 = (D, n - 60)
| n > 40 = (C, n - 40)
| n > 20 = (B, n - 20)
| n > 0 = (A, n)
| otherwise = int2TypeInt (n `mod` 108 + 1)
-- | test the number of positive roots is half that of all roots
prop_positiveRoots :: Int -> Bool
prop_positiveRoots n = prop_positiveRoots' t m
where (t, m) = int2TypeInt n
prop_positiveRoots' :: Type -> Int -> Bool
prop_positiveRoots' t n = (length $ positiveRoots (simpleSystem t n)) * 2 == numRoots t n
-- | test the positive roots and negative roots form all the roots
prop_positiveRoots1 :: Int -> Bool
prop_positiveRoots1 n = prop_positiveRoots1' t m
where (t, m) = int2TypeInt n
prop_positiveRoots1' :: Type -> Int -> Bool
prop_positiveRoots1' t n =
let pr = positiveRoots (simpleSystem t n) in
(S.fromList $ pr ++ ((-1) *>> pr)) == (S.fromList $ allRoots (simpleSystem t n))
-- Given a simple system, return the full root system
-- The closure of a set of roots under reflection
allRoots :: SimpleSystem -> [[Q]]
allRoots ss = S.toList $ closure S.empty (S.fromList ss) where
closure interior boundary
| S.null boundary = interior
| otherwise =
let interior' = S.union interior boundary
boundary' = S.fromList [s alpha beta | alpha <- ss, beta <- S.toList boundary] S.\\ interior'
in closure interior' boundary'
-- | Test that the length of the longest element is half of the number of roots
prop_longestElement :: Int -> Bool
prop_longestElement n = let (t, m) = int2TypeInt n in
numRoots t m == 2 * (length $ longestElement $ simpleSystem t m)
-- Old code in HaskellForMaths
{--
-- 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]] ]
-- 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)]
-}
--}