aboutsummaryrefslogtreecommitdiff
path: root/Math/Combinatorics/RootSystem.hs
blob: 716819cd640834bcb2d8a4a944bfdff5f450b35e (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
-- 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)]
-}
--}