(* :Title: Weyl.m *)

(* :Context: Projects`Weyl.m *)

(* :Author: Donald E. Taylor *)

(* :Summary:
   This package defines functions to compute with finite Weyl groups.
   It uses special properties of Weyl groups whenever it can and thus
   avoids doing general group theory.
*)

(* :Copyright: © 1999 by D. E. Taylor *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 3.0 *)

(* :History:
   Version 1.0 September 1999
*)

(* :Keywords: Weyl group, Cartan matrix, Lie algebra, 
   highest weight module, reflection *)

(* :Sources:
   N. Bourbaki, Groupes et algebres de Lie, Chapitres 4-6, Hermann, 1968.
   H. Freudenthal and H. de Vries, Linear Lie groups, Academic Press, 1969.
   J.F. Humphreys, Introduction to Lie algebras and representation theory,
   Springer Verlag, 1972.
   R.V. Moody and J. Patera, Fast recursion formula for weight multiplicities,
   Bull. Amer. Math. Soc. 7 (1982), 237-242.
   J. Stembridge, A Maple package for root systems and
   Coxeter groups. Maple Share Library, 1998.
*)

(* :Warnings:
  The package is restricted to finite Weyl groups and the functions are not
  guaranteed to work with Cartan matrices of types H3, H4 and I2.
*)

(* :Limitations:
   The function CoxeterForm is designed to work with the Cartan matrices
   generated internally by this package and will not necessarily return
   the correct result for user-defined Cartan matrices.  The functions
   which depend on CoxeterForm are OrbitLength, WeylDimension and
   WeightMultiplicity.
*)

(* :Discussion:
   A (finite) Weyl group is defined by its Cartan matrix.  In this package 
   we do not generate elements of the group itself but refer back to the 
   Cartan matrix and the basic reflections which it defines.  The roots, 
   weights and reflections are written with respect to the basis of simple 
   roots determined by the Cartan matrix.
     The primary goal of the package is to do simple weight space 
   calculations.
*)

(* :Requirements:
   none
*)

(* :Examples:
   c = CartanMatrix["G",2]
   fw = FundamentalWeights[c]
   WeylDimension[c,fw[[1]]+fw[[2]]
   WeightMultiplicity[c,fw[[1]]+fw[[2]],fw[[1]]] 
*)

BeginPackage["Projects`Weyl`"]

(* usage messages for the exported functions and the context itself *)

Weyl::usage = 
    "This package deals with finite Weyl groups only.  Coxeter groups
of types H and I are not fully supported"

CartanMatrix::usage = 
    "CartanMatrix[type,rank] returns the Cartan matrix of the given
(simple) type and rank"

MatrixBlocks::usage =
    "MatrixBlocks[M] returns a partition of {1..Length[n]} representing
the connected components of the graph defined by joining i to j whenever
i > j and M[[i,j]] != 0"

CartanType::usage =
    "CartanType[C] returns a list of pairs consisting of the string type
of each component of C together with the indices giving the location of
the component in C"

WeylOrder::usage =
    "WeylOrder[C] is the order of the Weyl group with Cartan matrix C"

Degrees::usage =
    "Degrees[C] returns the sorted list of degrees of the basic invariants
of the Weyl group with Cartan matrix C"

Exponents::usage =
    "Exponents[C] returns the sorted list of exponents of the Weyl group
with Cartan matrix C"

NumberOfReflections::usage =
    "NumberOfReflections[C] is the number of reflections in the Weyl group
with Cartan matrix C"

PositiveRoots::usage =
    "PositiveRoots[C] returns the list of positive roots in the root system
determined by the Cartan matrix C"

Reflection::usage =
    "Reflection[C,k] returns the k-th simple reflection matrix corresponding
to the Cartan matrix C"

SimpleReflections::usage =
    "SimpleReflections[C] is the list of simple reflections of the
Weyl group with Cartan matrix C"

FundamentalWeights::usage =
    "FundamentalWeights[C] is the list of fundamental weights of the Weyl 
group with Cartan matrix C.  It coincides with the inverse of C"

DominantWeights::usage =
    "DominantWeights[C,lambda] is the list of dominant
weights <= lambda for the root system with Cartan matrix C,
where lambda is itself a dominant weight"

Orbit::usage =
    "Orbit[C,v] returns the orbit of the vector v under the
action of the Weyl group with Cartan matrix C.
    Orbit[C,v,subs] uses only the reflections indexed by subs"

OrbitLength::usage = 
    "OrbitLength[C,v,options] is the length of the orbit of the vector
v under the action of the Weyl group with Cartan matrix C.
    OrbitLength[C,v,subs,options] uses only the roots and reflections
indexed by subs"

Extended::usage = "Extended is an option for OrbitLength that specifies
whether or not to include -1 in the group defining the action"

Saturation::usage =
    "Saturation[C,lambda] returns all weights of the irreducible module
V(lambda) with highest weight lambda for the Lie algebra over the complex
numbers with Cartan matrix C"

CoxeterForm::usage =
    "CoxeterForm[C] returns a matrix describing the bilinear form
associated with the Cartan matrix C"

WeylDimension::usage =
    "WeylDimension[C,lambda] returns the dimension of the representation of
the irreducible module with highest weight lambda for the Lie algebra with
Cartan matrix C"

DominantRepresentative::usage =
    "DominantRepresentative[C,v] returns the unique element in the W-orbit
of v which lies in the fundamental Weyl chamber, where W is the Weyl
group with Cartan matrix C.
    DominantRepresentative[C,v,subs] uses only the roots and reflections
indexed by subs"

DominantSign::usage =
    "DominantSign[C,v] is the 0 if no element of the Weyl group W of the
Cartan matrix C fixes v, otherwise it is the sign of the unique element
of W that takes v to a dominant weight"

DominantQ::usage =
    "DominantQ[C,v] yields True if the vector v is a dominant weight with
respect to the Cartan matrix C"

WeightMultiplicity::usage =
    "WeightMultiplicity[C,lambda,mu] returns the dimension of the weight
space of mu in the irreducible module with heighest weight lambda for the
Lie algebra with Cartan matrix C"

TensorProduct::usage =
    "TensorProduct[C,mu,xi] is a list of pairs {m,lambda} where m is the
multiplicity of the irreducible module V(lambda) in the tensor product
of V(mu) and V(xi).  It is best if the dimension of V(mu) is small."

CoxeterNumber::usage =
    "CoxeterNumber[C] is the Coxeter number of the Weyl group with Cartan
matrix C"

(* error messages for the exported objects *)

CartanMatrix::badtype = "The type `` is not supported"

CartanType::inf = "Not a finite Weyl group"
CartanType::cycle = "Dynkin diagram of `` contains a cycle"

CoxeterForm::unrec =
    "`` is not an irreducible Cartan matrix"

Degrees::unknown = "Unknown Cartan type ``"

DominantWeights::notdom = "The weight `` is not dominant"
DominantWeights::incompat = "Matrix and weight have incompatible dimensions"

Begin["`Private`"]

CartanMatrix[ type_String, rank_Integer ] :=
  Module[{ m = Table[Switch[i-j,-1,-1,0,2,1,-1,_,0],{i,rank},{j,rank}] },
  (Switch[ type, 
  "A"|"a", Null,
  "B"|"b", m[[2,1]] = -2,
  "C"|"c", m[[1,2]] = -2,
  "D"|"d", m[[1,2]] = m[[2,1]] = 0;
           m[[1,3]] = m[[3,1]] = -1,
  "E"|"e", m[[1,2]] = m[[2,1]] = m[[2,3]] = m[[3,2]] = 0;
           m[[1,3]] = m[[3,1]] = m[[2,4]] = m[[4,2]] = -1,
  "F"|"f", m[[3,2]] = -2,
  "G"|"g", m[[2,1]] = -3,
  "H"|"h", m[[1,2]] = m[[2,1]] = -GoldenRatio,
  "I2"|"i2",
    m = With[ {z = Exp[2Pi I/rank],
    	       w = Exp[2Pi I(rank+1)/rank]+Exp[-2Pi I(rank+1)/rank]},
      Switch[ rank,
        3, {{2,-1},{-1,2}},
        4, {{2,-1},{-2,2}},
        6, {{2,-1},{-3,2}},
        _, If[EvenQ[rank],{{2,-1},{-2-z-z^-1,2}},{{2,w},{w,2}}]
      ]
    ],
  _,       Message[CartanMatrix::badtype,type];Return[$Failed]
  ];  m) /;
  (Switch[ type,
  "A"|"a", 1 <= rank,
  "B"|"b", 2 <= rank,
  "C"|"c", 2 <= rank,
  "D"|"d", 4 <= rank,
  "E"|"e", 6 <= rank && rank <= 8,
  "F"|"f", rank == 4,
  "G"|"g", rank == 2,
  _, True])
  ]

(* Non-crystallographic systems -- not implemented
Hn[ k_ ] :=
  Module[ {m = An[k]},
    m[[1,2]] = m[[2,1]] = -GoldenRatio;
    m
  ]

I2[ k_ ] :=
  With[ {z = Exp[2Pi I/k], w = Exp[2Pi I(k+1)/k]+Exp[-2Pi I(k+1)/k]},
    Switch[ k,
      3, {{2,-1},{-1,2}},
      4, {{2,-1},{-2,2}},
      6, {{2,-1},{-3,2}},
      _, If[EvenQ[k],{{2,-1},{-2-z-z^-1,2}},{{2,w},{w,2}}]
    ]
  ]
*)

(* Version 1
MatrixBlocks[M_?MatrixQ] := 
  Module[{i,j,k,l,p,n = Length[M],X},
    X = Range[n];
    For[ i = 2, i <= n, i++,
      For[ j = 1, j < i, j++,
        If[ M[[i,j]] == 0, Continue[] ];
        p = X[[{i,j}]];
        k = Min[p];
        l = Max[p];
        If[ k == l, Message[MatrixBlocks::cycle,M]; Return[$Failed] ];
        X = X /. { l->k };
        ]
      ];
    Map[Flatten[Position[X,#]]&,Union[X]]
    ]
*)

(* Version 2
MatrixBlocks[M_?MatrixQ] := 
  Module[{f,res},
    With[{n=Length[M],c=M},
      f[X_,p_] := With[{q = X[[p]]},With[{k=Min[q],l=Max[q]},
        If[k==l, Message[MatrixBlocks::cycle,c];Throw[$Failed], X/.{l->k}]]];
      res = Catch[Fold[f,Range[n],
        Select[Join @@ Table[Table[{i,j},{i,j-1}],{j,2,n}],
        c[[Sequence @@ #]] != 0&]]]
      ];
    If[res === $Failed, res, Map[Flatten[Position[res,#]]&,Union[res]]]
    ]
*)

MatrixBlocks[M_?MatrixQ] := 
  Module[{f,res},
    With[{n=Length[M]},
      f[X_,p_] := With[{q = X[[p]]},With[{k=Min[q],l=Max[q]}, X/.{l->k}]];
      res = Fold[f,Range[n],
        Select[Join @@ Table[{i,j},{j,2,n},{i,j-1}],
        M[[Sequence @@ #]] != 0&]]
      ];
    Flatten[Position[res,#]]& /@ Union[res]
    ]

CartanType[c_?MatrixQ] := 
    Map[simpleCartanType,Map[c[[#,#]]&,MatrixBlocks[c]]]

simpleCartanType[c_?MatrixQ] :=
  Module[{n=Length[c],edges,strong,gamma,deg,ends,fork,s,m},
    If[ n == 1, Return[{"A",1}] ];
    edges = Select[Join @@ Table[{i,j}, {i,n}, {j,n}],
        c[[Sequence @@ #]] < 0 &];
    If[ 2n != Length[edges]+2, Message[CartanType::cycle,c];Return[$Failed] ];
    strong = Select[edges,c[[Sequence @@ #]] < -1&];
    If[ Length[strong] > 1, Message[CartanType::inf];Return[$Failed] ];
    strong = Flatten[strong];
    m = If[ Length[strong] == 0, -1, c[[Sequence @@ strong]] ];
    With[{et = Transpose[edges]}, (* List of neighbours *)
      gamma = Last[et][[Flatten[Position[First[et],#]]]]& /@ Range[n]];
    deg = Length /@ gamma;
    ends = Flatten[Position[deg,1]];
    Switch[ Length[ends],
      2, If[ m == -1, Return[{"A",n}] ];
         s = Intersection[ends,strong];
         Switch[ Length[s],
           0, If[ n == 4 && m == -2, Return[{"F",4}] ],
           1, If[ m == -2,
               If[ strong[[1]] == s[[1]], Return[{"C",n}], Return[{"B",n}]],
               ],
           2, If[ m == -2, Return[{"C",2}] ];
              If[ m == -3, Return[{"G",2}] ]
           ],
      3, s = Length[Intersection[ ends,
               gamma[[First[Flatten[Position[deg,3]]]]] ]];
         If[ s == 1 && n < 9, Return[{"E",n}] ];
         If[ s > 1 && n > 3,  Return[{"D",n}] ]
      ];
    Message[CartanType::inf];
    $Failed
    ]

PositiveRoots[ c_?MatrixQ ] :=
  With[{n=Length[c]},
    With[{b=IdentityMatrix[n],
      t = Select[Join @@ Table[{i,j}, {i,1,n-1}, {j,i+1,n}],
          c[[Sequence @@ #]] < 0 &]},
    Join @@ allLayers[c,{b},Apply[Function[{x,y},b[[x]]+b[[y]]],t,{1}]]
  ]]

allLayers[ c_, layer_, {} ] := layer
allLayers[ c_, layer_, r_ ] := 
  Module[{f},
    With[{b=First[layer],h=Length[layer],ndx=Range[Length[c]]},
      f[R_,beta_] := Module[{bC = beta.c},
        Join[R,beta + b[[#]]& /@ Select[ndx,
          Function[i,With[{ip=bC[[i]],alpha=b[[i]]},
            (ip < 0 || (h > ip && beta[[i]] > ip &&
             MemberQ[layer[[h-ip]],beta-(ip+1)alpha])) &&
             FreeQ[R,beta+alpha] ]]]]]
      ];
    allLayers[c,Append[layer,r],Fold[f,{},r]]
  ]

Reflection[c_?MatrixQ, k_Integer] :=
  Module[{M = IdentityMatrix[Length[c]], t = Transpose[c]},
    M[[k]] = M[[k]] - t[[k]]; Transpose[M]
    ]
Reflection[c_?MatrixQ, lst_List] := Reflection[c,#]& /@ lst

SimpleReflections[c_?MatrixQ] := Reflection[c,Range[Length[c]]]

FundamentalWeights[c_?MatrixQ] := Inverse[c]

DominantWeights[ c_?MatrixQ, lambda_?VectorQ ] := 
  Module[{ R, dom, ndx, u, v },

  If[ Length[c] != Length[lambda],
    Message[DominantWeights::incompat]; Return[$Failed]];
  If[ Min[lambda . c] < 0,
    Message[DominantWeights::notdom,lambda]; Return[$Failed]];
  
  R = PositiveRoots[c];
  dom = { lambda };
  ndx = 0;
  While[ ndx < Length[dom], ndx++;
    u = dom[[ndx]];
    dom = Join[dom,Select[ (u - #&) /@ R,
      Function[v, Min[v . c ] >= 0 && FreeQ[dom,v]]]]
    ];
  dom
  ]

Orbit[ c_?MatrixQ, v_?VectorQ ] := orb[SimpleReflections[c],v]
Orbit[ c_?MatrixQ, v_?VectorQ, subs_List ] := 
                               orb[SimpleReflections[c][[subs]],v]

orb[ gen_List, v_?VectorQ ] := 
  Module[{o={v},ndx=0},
    While[ ndx < Length[o],
      ndx++;
      u = o[[ndx]];
      o = Join[o, Select[ (u . #&) /@ gen, FreeQ[o,#]& ]]
      ];
    o
    ]

Saturation[c_?MatrixQ,lambda_?VectorQ] :=
    Join @@ (Orbit[c,#]& /@ DominantWeights[c, lambda])

Options[OrbitLength] = {Extended->False}

OrbitLength[c_?MatrixQ, v_?VectorQ, subs_List, opts___?OptionQ ] :=
  Module[{u = DominantRepresentative[c,v,subs],
  	  b = v . CoxeterForm[c],d,extend,res},
    extend = Extended /. {opts} /. Options[OrbitLength];
    d = Select[subs, b[[#]] == 0& ];
    res = WeylOrder[c[[subs,subs]]]/WeylOrder[c[[d,d]]];
    If[extend && u != DominantRepresentative[c,-v,subs], 2 res, res]
    ]
OrbitLength[c_?MatrixQ, v_?VectorQ, opts___?OptionQ ] :=
  OrbitLength[c,v,Range[Length[c]],opts]

(* This relies completely on the form of the Cartan
matrix produced by this package *)
irrRootNorms[c_?MatrixQ] :=
  With[ { n = Length[c] },
    If[ c == Transpose[c], Table[1,{n}],
      If[ c[[1,2]] == -1,
        Switch[ c[[2,1]],
          -1, {1/2,1/2,1,1},
          -2, {1/2,Table[1,{n-1}]},
          -3, {1/2,3/2},
           _, Message[CoxeterForm::unrec, c ]; Return[$Failed]
          ],
        {1,Table[1/2,{n-1}]}
        ]
      ]
    ]

CoxeterForm[c_?MatrixQ] :=
  c . DiagonalMatrix[Flatten[Map[irrRootNorms,
        Map[c[[#,#]]&,MatrixBlocks[c]]]]]

WeylDimension[c_?MatrixQ, lambda_?VectorQ] := 
  With[{ R = PositiveRoots[c], b = CoxeterForm[c] },
    With[{ delta = (Plus @@ R) / 2 },
      With[{ s = (lambda+delta) . b, t = delta . b },
        Times @@ ((s . #)/(t . #)&) /@ R
        ]
      ]
    ]

DominantRepresentative[c_?MatrixQ, u_?VectorQ, subs_List] :=
  Module[{ i, v = u, vc = u . c },
    With[{ S = SimpleReflections[c] },
      While[ True,
        i = Select[ subs, vc[[#]] < 0&, 1 ];
        If[Length[i] == 0, Break[]];
        v = v . S[[First[i]]];
        vc = v . c
        ]
      ];
    v
    ]
DominantRepresentative[c_?MatrixQ, u_?VectorQ] :=
  DominantRepresentative[c,u,Range[Length[c]]]

DominantQ[c_?MatrixQ,v_?VectorQ] := Min[ v . c ] >= 0

DominantSign[c_?MatrixQ, u_?VectorQ] :=
  Module[{ i, v = u, vc = u . c, s },
    With[{ S = SimpleReflections[c], ndx = Range[Length[c]] },
      s = 1;
      While[ True,
        i = Select[ ndx, vc[[#]] < 0&, 1 ];
        If[Length[i] == 0, Break[]];
        v = v . S[[First[i]]];
        vc = v . c;
        s = -s
        ];
      s = If[ Length[Select[ndx, vc[[#]] == 0&, 1]] == 0, s, 0]
      ];
    s
    ]

(* From Stembridge *)
WeightMultiplicity[ c_?MatrixQ, lambda_?VectorQ, mu_?VectorQ ] :=
  Module[{u,v,wts,R,delta,n,ip,wl,mults,i,J,m,k,uu,j,orblen,ll,beta,ipb,ldip},
    u   = DominantRepresentative[c,mu];
    wts = DominantWeights[c,lambda];
    If[ FreeQ[wts,u], Return[0] ];

    wts = First /@ Sort[Thread[List[wts,Apply[Plus, wts, 1]]],
          Last[#1] > Last[#2]&];
    n  = First[First[Position[wts,u]]];
    ll = Range[Length[c]];
    ip = CoxeterForm[c];
    R  = PositiveRoots[c];
    delta = (Plus @@ R) / 2;
    ldip = (lambda+delta) . ip . (lambda+delta);
    wl = wts . c;
    mults[1] = 1;
    For[ i = 2, i <= n, i++,
      With[{alpha=wl[[i]]},
        J = Select[ll, alpha[[#]] == 0& ]];
      v = wts[[i]];
      m = 0;
      For[ k = 1, k <= Length[R], k++,
        beta = R[[k]];
        With[{alpha=beta . c},
          If[ Min[Map[alpha[[#]]&, J]] < 0, Continue[] ]];
        u = v + beta;
        uu = DominantRepresentative[c,u];
        If[ FreeQ[wts,uu], Continue[] ];
        j = First[First[Position[wts,uu]]];
        orblen = OrbitLength[c,beta,J,Extended->True];
        ipb = ip . beta;
        While[True,
          m += mults[j] orblen (u . ipb);
          u += beta;
          uu = DominantRepresentative[c,u];
          If[ FreeQ[wts,uu], Break[]];
          j = First[First[Position[wts,uu]]]
          ]
        ];
      mults[i] = m/(ldip - (v+delta).ip.(v+delta))
      ];
    mults[n]
    ]

(* Unfinished - from Stembridge
TensorProductWeights[c_?MatrixQ,u_?VectorQ,v_?VectorQ] :=
  Module[{wts,wl,cow,ip,coRoots,lb,n,i,k},
    ip = CoxeterForm[c];
    R = PositiveRoots[c];
    wts = DominantWeights[c,u+v];
    coRoots = Table[2 r[[i]]/(r[[i]] . ip . r[[i]]),{i,Length[c]}];
    cow = coRoots . c;
    wl = wts . c;
    
    ]
*)

TensorProduct[c_?MatrixQ,u_?VectorQ,v_?VectorQ] :=
  Module[{sat,delta,ndx,dr,mul,wts},
    delta = (Plus @@ PositiveRoots[c]) / 2;
    sat = Saturation[c,u];
    sgn = DominantSign[c,v+delta+#]& /@ sat;
    ndx = Select[Range[Length[sat]],sgn[[#]] != 0&];
    sat = sat[[ndx]];
    sgn = sgn[[ndx]];
    dr  = DominantRepresentative[c,v+delta+#]-delta& /@ sat;
    mul = MapThread[Times,{sgn,WeightMultiplicity[c,u,#]& /@ sat}];
    wts = Union[dr];
    mul = Map[Apply[Plus,mul[[#]]]&,Map[Flatten[Position[dr,#]]&,wts]];
    Transpose[{mul,wts}]
    ]

Degrees[c_?MatrixQ] := Flatten[ irrDegrees /@ CartanType[c] ]

irrDegrees[{t_,n_}] := 
  Switch[ t,
    "A", Table[i,{i,2,n+1}],
    "B"|"C", Table[2i, {i,n}],
    "D", Append[Table[2i,{i,n-1}],n],
    "E",
      Switch[ n,
        6, {2,5,6,8,9,12},
        7, {2,6,8,10,12,14,18},
        8, {2,8,12,14,18,20,24,30}
        ],
    "F", {2,6,8,12},
    "G", {2,6},
     _, Message[Degrees::unknown,t]; Return[$Failed]
     ]

Exponents[c_?MatrixQ] := Degrees[c]-1

WeylOrder[{}] := 1
WeylOrder[c_?MatrixQ] := Times @@ Degrees[c]

CoxeterNumber[c_?MatrixQ] := LCM @@ Max /@ irrDegrees /@ CartanType[c]

NumberOfReflections[c_?MatrixQ] := Plus @@ Exponents[c]

End[]

Protect[CartanMatrix,CartanType,CoxeterForm,CoxeterNumber,Degrees,DominantQ,
      DominantSign,DominantRepresentative,DominantWeights,Exponents,
      FundamentalWeights,MatrixBlocks,NumberOfReflections,Orbit,OrbitLength,
      PositiveRoots,Reflection,Saturation,SimpleReflections,TensorProduct
      WeightMultiplicity,WeylDimension,WeylOrder]

EndPackage[ ]