#############################################################################
##
#A  weylgrp.g                  GAP Library                       Meinolf Geck
##
#Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##
##  This file  contains general functions that  deal with  Weyl groups and their
##  root  systems
##
# Meinolf Geck 16-12-93:
# added function 'HeighestShortRoot'

# Jean Michel 1-10-93:
# The introduction of Weyl subgroups means that an expression like
#       while i ^ ww <= W.N  do
# has to be rewritten either as
#       while W.rootIndices[i] ^ ww <= W.parentN  do
# or as
#       while i ^ (ww ^ W.weylFusion) <= W.N  do
#  (for programs to work with Weyl Subgroups)

#############################################################################
##
#F  CartanMat( <type>, <n> )  . . . .  a Cartan matrix of given type and rank
##
##  'CartanMat' returns  the  Cartan matrix  of type  <type> and rank <n>.
##  E.g., 'CartanMat( "F", 4 );'
##  For type I2(n), a third argument is needed: 'CartanMat( "I", 2, n )'.
##
CartanMat := function ( arg )
    local tr, t, n, An, r, a, i, j;
    if Length(arg)=1 then
      tr:=arg[1];
    else
      tr:=arg;
    fi;
    t:=tr[1];n:=tr[2];
    An := function ( n )
          local  a, i, j;
          a := [  ];
          for i  in [ 1 .. n ]  do
              a[i] := [  ];
              for j  in [ 1 .. n ]  do
                  if i = j  then
                      a[i][j] := 2;
                  elif j = i + 1 or j = i - 1  then
                      a[i][j] := -1;
                  else
                      a[i][j] := 0;
                  fi;
              od;
          od;
          return a;
      end;
    if Length(tr)=3 then
        if t <> "I" or n<>2 then return "funny type";fi;
	n:=tr[3];
	if n=2 then 
	   return [[2,0],[0,2]];
	elif n mod 2=0 then
           return  [[2,-1],[-2-E(n)-E(n)^-1,2]];
	else 
	   return  [[2,-E(2*n)-E(2*n)^-1],[-E(2*n)-E(2*n)^-1,2]];
        fi;
    fi;
    if n < 1 or not IsInt( n )  then
        return "funny rank";
    fi;
    if t = "A" or t = "a"  then
        return An( n );
    elif t = "B" or t = "b"  then
        if n < 2  then
            return "funny rank";
        fi;
        a := An( n );
        a[1][2] := -2;
        return a;
    elif t = "C" or t = "c"  then
        if n < 2  then
            return "funny rank";
        fi;
        a := An( n );
        a[2][1] := -2;
        return a;
    elif t = "D" or t = "d"  then
        if n < 3  then
            return "funny rank";
        fi;
        a := An( n );
        a[1][3] := -1;
        a[3][1] := -1;
        a[1][2] := 0;
        a[2][1] := 0;
        return a;
    elif t = "G" or t = "g"  then
        if n <> 2  then
            return "funny rank";
        fi;
        a := An( 2 );
        a[2][1] := -3;
        return a;
    elif t = "F" or t = "f"  then
        if n <> 4  then
            return "funny rank";
        fi;
        a := An( 4 );
        a[3][2] := -2;
        return a;
    elif t = "E" or t = "e"  then
        a := [ [ 2, 0, -1, 0, 0, 0, 0, 0 ], [ 0, 2, 0, -1, 0, 0, 0, 0 ], 
            [ -1, 0, 2, -1, 0, 0, 0, 0 ], [ 0, -1, -1, 2, -1, 0, 0, 0 ], 
            [ 0, 0, 0, -1, 2, -1, 0, 0 ], [ 0, 0, 0, 0, -1, 2, -1, 0 ], 
            [ 0, 0, 0, 0, 0, -1, 2, -1 ], [ 0, 0, 0, 0, 0, 0, -1, 2 ] ];
        if n in [6..8] then
            return a {[1..n]}{[1..n]};
        else
            return "funny rank";
        fi;
    elif t = "H"  then
        r := (1 + ER( 5 )) / 2;
        a := [ [ 2, -1 * r, 0, 0 ], [ -1 * r, 2, -1, 0 ], [ 0, -1, 2, -1 ], 
                [ 0, 0, -1, 2 ] ];
        if n in [2..4] then
            return a {[1..n]}{[1..n]};
        else
            return "funny rank";
        fi;
    else
        return "funny type";
    fi;
end;

#############################################################################
##
#F  DirectSumCartanMat( <mat1>, <mat2> )  . direct sum of two Cartan Matrices
##
##  'DirectSumCartanMat' returns the block diagonal direct sum of the  Cartan
##  matrices <mat1>  and <mat2>.
##
DirectSumCartanMat := function ( a, b )
    local  c, n, m;
    n := Length( a );
    m := Length( b );
    c := List([1..n+m],x->[1..n+m]*0);
    c{[1..n]}{[1..n]}:=a;
    c{[n+1..n+m]}{[n+1..n+m]}:=b;
    return c;
end;

#############################################################################
##
#F  SimpleReflectionMatrices( <mat> )  . . matrices of the simple reflections
##
##  'SimpleReflectionMatrices'  returns the   matrices  (in row   convention)
##  of the simple  reflections of  the Weyl group  determined  by  the Cartan
##  matrix  <mat> with  respect to  the  basis consisting of the  fundamental
##  root vectors corresponding to the rows of <mat>.
##  (This function is only used in the function 'Weyl'.)
##
SimpleReflectionMatrices := function( A )
   local  l, k, e, g;
    l := Length( A );
    e := A ^ 0;
    g := List([ 1 .. l ], x->Copy(e));
    for k  in [ 1 .. l ]  do
            g[k]{[1..l]}[k] := e[k]-A[k];
    od;
    return g;
end;

#############################################################################
##
#F  Rootsystem( <mat> )  . . . .the root system determined by a Cartan matrix
##
##  'Rootsystem' returns the positive roots of the root system determined  by
##  the Cartan matrix <mat>,  given  by  their coefficients when expressed as
##  linear combinations of fundamental roots. Thus the fundamental roots  are
##  the standard basis vectors e_i, corresponding to the rows i=1..n of <mat>.
##  The roots  are  ordered  by  increasing  height.  Note  that  we  use the 
##  convention that <mat>[i,j]=2(e_i,e_j)/(e_i,e_i).
##  (This function is only used in the function 'Weyl'.)
##
Rootsystem := function ( A )
    local   R, a, b, t, i, v, r;
    R := A ^ 0;
    if ForAll( Concatenation( A ), IsInt )  then
        for a  in R  do
            v := A * a;
            for i  in [ 1 .. Length( A ) ]  do
                if (v[i]<0 or a-(v[i]+1)*R[i] in R) and not a+R[i] in R then
                        Add( R, a+R[i] );
                fi;
            od;
        od;
        return R;
    else
        t := SimpleReflectionMatrices( A );
        for b  in R  do
            for a  in t  do
                if Position( R, b ) <> Position( t, a )  then
                    v := b * a;
                    if not v in R  then
	#    Print("m",Position(t,a),"*(",pqv(b),")->(",pqv(v),")\n");
                        Add( R, v );
                    fi;
                fi;
            od;
        od;
        return  R;
    fi;
end;
#############################################################################
##
#F  PermRepresentationRoots( <mats> , <roots> )  . . . . . . . . . . . . . .
#F  . . . . . . . . . . . .the permutation representation on the root vectors
##
##  'PermRepresentationRoots'   returns  the  permutations  induced  by   the
##  fundamental   reflections   <mats>   of   a  given   Weyl  group  on  the 
##  corresponding root sytem <roots>.  If  there  are 2N roots in total, then 
##  the first N roots are positive and the last  N roots negative.
##  (This function is only used in the function 'Weyl'.)
##
PermRepresentationRoots := function( mats, r)
    local  p, l, R, perm, perms;
    R := Concatenation( r, -r );
    l:=Sortex(R)^-1;
    p := [  ];
    perms:=r*mats;
    for perm in perms do
        Add(p, Sortex(Concatenation( perm, -perm ))*l);
        perm := Concatenation( perm, -1 * perm );
        l := [  ];
        for i  in [ 1 .. n ]  do
            Add( l, Position( R, perm[i] ) );
        od;
        Add( p, PermList( l ) );
    od;
    return p;
end;

#############################################################################
##
#F  DecomposedMat( <mat> ) . . . . . . . . . Find if the square matrix mat
#F  admits a block decomposition.
##  
##  Define a graph G with vertices [1..Length(mat)] and with an edge
##  between i and j if either mat[i][j] or mat[j][i] is non-zero.
##  DecomposedMat return a list of lists l such that l[1],l[2], etc..
##  are the vertices in each connected component of G. In other words,
##  mat{l[1]}{l[1]},mat{l[2]}{l[2]},etc... are blocks of the matrix mat
## 
##  This routine is general in scope and has nothing to do in this
##  particular file.
##
DecomposedMat:=function(mat)local r,sets,todo,set,x,done;
 r:=Length(mat);
 sets:=[];
 todo:=[1..r];
 while Length(todo)>0 do
   set:=[todo[1]];done:=[];
   while Length(done)<Length(set) do
     x:=Difference(set,done);
     x:=x[1];
     AddSet(done,x);
     UniteSet(set,Filtered([1..r],y->mat[x][y]<>0));
   od;
   SubtractSet(todo,set);
   Add(sets,set);
 od;
 return sets;
 end;

##

#############################################################################
##
#F  CartanType( <C>)  . . . . . . . type of a Cartan Matrix
##
##  Returns the type of the Cartan Matrix C, i.e. a list  l which for
##  each irreducible component of C contains a pair such that
##  if l[i]=[t,v] and we put C_i:=CartanMat(t,Length(v)) then
##  C{Concatenation(v_1,...,v_r)}=DirectSumCartanMat(C_1,...,C_r)
##  (excepted that for an irreducible component of type I2(n) the
##   pair is replaced by a triple ["I",v,n] )
##
CartanType:=function(C)local type,s,t, m,r,x,list,branch,l,p;
   type:=[];
   for s in DecomposedMat(C) do
     m:=C{s}{s};
     r:=Length(m);
     list:=[1];
     repeat
       x:=Filtered([1..r],i->m[list[Length(list)]][i]<>0 and not i in list);
       if Length(x)>0 then Add(list,x[1]);fi;
     until Length(x)=0;
     repeat
       x:=Filtered([1..r],i->m[list[1]][i]<>0 and not i in list);
       if Length(x)>0 then list:=Concatenation([x[1]],list);fi;
     until Length(x)=0;
     l:=Length(list);
     if l=r then # types A,B,C,F,G,H,I
       if not ForAll(Concatenation(m),IsInt) then
	  if r>2 then 
	    if m[list[1]][list[2]]=-1 then list:=Reversed(list);fi;
	    t:=["H",list];
	  else
	    x:=Maximum(List(Concatenation(m),NofCyc));
	    if x mod 2<>0 and m[list[1]][list[2]]<>m[2][1] then x:=2*x; fi;
	    if m[list[1]][list[2]]<>m[list[2]][list[1]] and
	       m[list[1]][list[2]]<>-1 then list:=Reversed(list);fi;
	    if x=5 then t:=["H",list];
	    else t:= ["I",list,x];
	    fi;
	  fi;
       elif -3 in Concatenation(m) then
	    if m[list[1]][list[2]]<>-1 then list:=Reversed(list);fi;
	    t:= ["G",list];
       elif -2 in Concatenation(m) then
	    if m[list[1]][list[2]]<>-2 and m[list[2]][list[1]]<>-2 then
	      list:=Reversed(list);
	    fi;
	    if m[list[1]][list[2]]<>-2 and m[list[2]][list[1]]<>-2 then
	      if m[list[3]][list[2]]<>-2 then list:=Reversed(list);fi;
	      t:=["F",list];
	    elif m[list[1]][list[2]]=-2 then t:= ["B",list];
	    else t:= ["C",list];
	    fi;
       else t:= ["A",list];
       fi;
     else # types D,E
       branch:=Filtered(list,
                        x->Number([1..r],i->m[x][i]<>0 and not i in list)<>0);
       if Length(branch)<>1 then Error("not a Cartan matrix");fi;
       repeat
	x:=Filtered([1..r],i->m[branch[Length(branch)]][i]<>0
	                      and not i in branch and not i in list);
	if Length(x)>0 then Add(branch,x[1]);fi;
       until Length(x)=0;
       if Length(branch)+l<>r+1 then Error("not a Cartan matrix");fi;
       p:=Position(list,branch[1]);
       if Length(branch)>=l+1-p then 
	 x:=list{[p..l]};
	 list:=Concatenation(list{[1..p-1]},branch);
	 branch:=x;
	 l:=Length(list);
       fi;
       if Length(branch)>=p then
	 x:=list{[1..p]};
	 list:=Concatenation(Reversed(branch),list{[p+1..l]});
	 p:=Length(branch);
	 branch:=Reversed(x);
	 l:=Length(list);
       fi;
       if p-1 > l-p then list:=Reversed(list);p:=l+1-p;fi;
       if p=2 then t:=["D",Concatenation([branch[2]],list)];
       else  
         t:=["E",Concatenation([list[1]],[branch[2]],list{[2..Length(list)]})];
       fi;
     fi;
     t[2]:=s{t[2]};
     Add(type,t);
    od;
    return type;
    end;

#############################################################################
##
#F  WeylName( <type> ) . . . . . . . . . Find the type of a Weyl group
##  
##  Gives a string which describes the type of W (like "A2B3" for the
##  direct product of a group of type A2 by one of type B3, or like "I2(7)").
##  Call as Weylname(W.type);
##
WeylName:=function(type)local set,res,t;
  res:=[];
  for t in type do 
    Append(res,t[1]);
    Append(res,String(Length(t[2])));
    if Length(t)=3 then
       Append(res,"(");
       Append(res,String(t[3]));
       Append(res,")");
    fi;
  od;
  return res;
  end;

#############################################################################
##
#F  ReflectionDegreesWeyl( <W>)  . . . . . . . degrees as a reflection group
##
##  Returns the degrees of the Crystallographic reflection group W,
##  as an increasing list of integers.
##  Returns currently false if the Weyl Group is not Crystallographic
##  
ReflectionDegreesWeyl:= function(W)local p;
    if not IsBound(W.degrees) then
      if ForAll( Concatenation( W.cartan ), IsInt )  then
# Crystallographic group: general algorithm
        p:=List(Collected(List(W.roots{[1..W.N]},Sum)),x->x[2]);
        W.degrees:=Reversed(1+List([1..Maximum(p)],x->Number(p,y->y>=x)));
      else 
# else use the classification
        if Length(W.type)=1 then
          if W.type[1][1]="H" then 
            if Length(W.type[1][2])=4 then W.degrees:=[2,12,20,30];
              elif Length(W.type[1][2])=3 then W.degrees:=[2,6,10];
              else W.degrees:=[2,5];
            fi;
          else # type I2
             W.degrees:=[2,W.type[1][3]];
          fi;
        else
          W.degrees:=[];
          for p in W.type do
            Append(W.degrees,ReflectionDegreesWeyl(WeylSubgroup(W,p[2])));
          od;
          Sort(W.degrees);
        fi;
      fi;
    fi;
    return W.degrees;
    end;

#############################################################################
##
#V  WeylGroupOps  . . . . . . operation record for Weyl group category
##
##  'WeylGroupOps'  is  the  operation  record for   Weyl  groups.  It
##  contains   the  functions  for   domain  operation,    e.g.,    'Size'
##  as well  as  the functions  for group    operations, e.g., 'Centralizer'.
##
##  'WeylGroupOps' is initially a copy of 'PermGroupOps', thus Weyl groups
##  inherit the default group functions, e.g., 'DerivedSubgroup' and 'Index'.
##  However  'WeylGroupOps'   overlays  some of  those    functions with more
##  efficient ones, e.g., 'Size'.
##
WeylGroupOps := Copy( PermGroupOps );

WeylGroupOps.Size:= W->Product(ReflectionDegreesWeyl(W));

WeylGroupOps.Print:= function(W)
    Print( "Weyl( ", W.cartan);
    if not ForAll(W.parameter,x->x=1) then
     Print(",",W.parameter);
    fi;
    Print( ")");
    end;

WeylGroupOps.CharTable:=function(W)
  local t,p,res,rank,Res,tbl;
  if IsBound(W.charTable) and IsCharTable(W.charTable) then
    return W.charTable;
  fi;
  W.name:=WeylName(W.type);
  Res:=[];
  for t in W.type do
  rank:=Length(t[2]);
  p:=MappingPermListList([1..rank],t[2])^-1;
  if t[1]="A" then 
    tbl:=HeckeCharTableA(rank,W.parameter[1^p]);
  elif t[1]="B" or t[1]="C" then
    tbl:=HeckeCharTableB(rank,W.parameter{[1^p,2^p]});
  elif t[1]="D" then
    tbl:=HeckeCharTableD(rank,W.parameter[1^p]);
  elif t[1]="E" then
    tbl:=HeckeCharTableE(rank,W.parameter[1^p]);
  elif t[1]="F" then
    tbl:=HeckeCharTableF4(W.parameter{[2^p,3^p]});
  elif t[1]="G" then
    tbl:=HeckeCharTableG2(W.parameter{[1^p,2^p]});
  elif t[1]="H" then
    tbl:=HeckeCharTableH(rank,W.parameter[1^p]);
  elif t[1]="I" then
    if t[3] mod 2 =0 then
      tbl:=HeckeCharTableDihedral(t[3],W.parameter{[1^p,2^p]});
    else
      tbl:=HeckeCharTableDihedral(t[3],W.parameter[1^p]);
    fi;
  fi;
  tbl.classtext:=List(tbl.classtext,x->OnTuples(x,p^-1));
  IsCharTable(tbl);
  Add(Res,tbl);
  od;
  if Length(Res)=1 then 
    W.charTable:=Res[1];
  else
    res:=Res[1];
    for tbl in Res{[2..Length(Res)]} do
      res:=HeckeCharTableDirectProduct(res,tbl,0);
    od;
    W.charTable:=res;
  fi;
  if not IsBound(W.charTable.identifier) then
    W.charTable.identifier:=W.charTable.name;
  fi;
  if not IsBound(W.charTable.size) then
    W.charTable.size:=W.charTable.order;
  fi;
  return W.charTable;
end;

# the next function is absolutely necessary so that fusion maps are OK
  WeylGroupOps.ConjugacyClasses:=W->List(WeylConjugacyClasses(W),
     x->ConjugacyClass(W,PermWeylWord(W,x)));

#############################################################################
##
#F  Weyl( <mat> , [ <parameter>] )  . . . . . . . . . . create a group record
#F  for the Weyl group determined by the Cartan matrix <mat>.
##  Typically, <mat> can be taken as the result of the function 'CartanMat'.
##  The group is given as a permutation group on the set of roots.
##
##  'Weyl' returns a record with the following entries: 
##    permutation group entries:
##  isDomain, isGroup, identity, generators, operations, isPermGroup,
##  isFinite, "1", ... "n" (n=number of generators)
##    in addition
##  isWeylGroup: true
##  cartan    : the Cartan matrix <mat>,
##  dim       : the size of <mat>
##  degree    : the total number of roots
##  N         : the number of positive roots
##  roots     : the root vectors
##  matgens   : the matrices of the simple reflections
##  permgens  : the permutations on the root vectors
##  parameter : the parameters of the Hecke algebra ([1,..,1] by default)
##  type      : the type, as return by CartanType
##  W.rootIndices, W.weylFusion, W.parentN are meaningful for Weyl subgroups;
##
Weyl := function( arg)
    local  A,r, m, W, p;

    A:=arg[1];
    if not IsMat( A ) then
       return "funny input";
    fi;
    r := Rootsystem( A );
    m := SimpleReflectionMatrices( A );
    W:=Group(PermRepresentationRoots(m,r),());
    W.isWeylGroup := true;
    W.cartan := A;
    W.dim := Length( A );
    W.N := Length( r );
    W.degree := 2 * W.N;
    W.rootIndices := [1..W.degree];
    W.weylFusion:=();
    W.parentN:=W.N;
    W.roots := Concatenation(r,-r);
    W.matgens := m;
    W.permgens := W.generators;
    W.permref:=WeylReflectionsPerm(W);
    if Length(arg)=2 then Hecke(W,arg[2]);
    elif Length(arg)=1 then W.parameter := List( [ 1 .. W.dim ], i -> 1 );
    else Error("syntax Weyl(cartanmat) or Weyl(cartanmat, Heckeparameters)");
    fi;
    W.operations:=WeylGroupOps;
    W.type:=CartanType(W.cartan);
    return W;
end;

#############################################################################
##
#F  WeylLengthPerm( <W> , <w> )  . . . . . . . length of a permutation w in W
##
##  'WeylLengthPerm'  returns  the  length  of the  permutation  w in  W as a 
##  reduced expression in the standard generators.
##  E.g., <w> can be taken as the result of the function 'PermWeylWord'.
##  
WeylLengthPerm := function ( W, w )local res,r;
# it's a pity the line below is 5 times less efficient than the next 5 lines
#    return Number(W.rootIndices{[1..W.N]},r->r ^ w > W.parentN);
#
    res:=0;
    for r in W.rootIndices{[1..W.N]} do
      if r ^ w > W.parentN then res:=res+1; fi;
    od;
    return res;
end;

#############################################################################
##
#F  PermWeylWord( <W> , <w> )  . . . . . . .  convert a word to a permutation
##
##  'PermWeylWord'  returns the  permutation  on  the root vectors determined
##  by  the  element  w  of  W,  given  as  a  list  of integers representing
##  the standard generators.
##  
PermWeylWord := function ( W, w )
    local  r, a;
    a := ();
    for r  in W.permgens{w} do
	a := a * r;
    od;
    return a;
end;

#############################################################################
##
#F  WeylWordPerm( <W> , <w> )  . . .  convert a permutation to a reduced word
##
##  'WeylWordPerm'  returns  a  reduced word  in the standard generators of W 
##  determined by the permutation w on the root vectors.
##  
WeylWordPerm := function( W, w )
    local  i, l, ww;
    l := [  ];
    ww := w;
    while ww <> ()  do
        i := 1;
        while W.rootIndices[i] ^ ww <= W.parentN  do
            i := i + 1;
        od;
        Add( l, i );
        ww := W.permgens[i] * ww;
    od;
    return l;
end;

#############################################################################
##
#F  ReducedWeylWord( <W> , w )  . . . . . . . . . . a reduced word for w in W
##
##  'ReducedWeylWord'   returns a reduced expression for the element w, given 
##  as  an  arbitrary  list of  integers  where  each  entry  i  in this list
##  represents the  i-th standard  generator of W.
##
ReducedWeylWord := function( W, w )
    return WeylWordPerm(W,PermWeylWord(W,w));
end;

#############################################################################
##
#F  LongestWeylWord( <W> )  . . . . . . . . . . . .  the longest element in W
##
## 'LongestWeylWord' returns a reduced expression in the  standard generators
##  of the unique longest element of the Weyl group <W>.
##
LongestWeylWord := function( W )
    local  i, w, fertig;
    fertig := false;
    w := ();
    while not fertig  do
        i := 1;
        while i <= W.dim and W.rootIndices[i] ^ ( w ^ -1 ) > W.parentN do
            i := i + 1;
        od;
        if i > W.dim  then
            fertig := true;
        else
            w := w * W.permgens[i];
        fi;
    od;
    return WeylWordPerm( W, w );
end;

#############################################################################
##
#F  WeylReflectionsPerm( <W> )  . . . . . . . . . . . .  the reflections in W
##
## 'WeylReflections'  returns  the set of reflections in the the Weyl
##  group <W>, as a list of permutations. The i-th entry in this list
##  is the reflection along the i-th root in <W>.roots.
##
WeylReflectionsPerm := function(W)
    local  i, j, r, x, w, res;
    w := ();
    W.permref := Copy(W.permgens);
    while true do
      i := 1;
      while i <= W.dim and W.rootIndices[i] ^ ( w ^ -1 ) > W.parentN do
	  i := i + 1;
      od;
      if i > W.dim  then
	  return W.permref;
      else
	  x:=w*W.permgens[i]*w^-1;
	  j := 1;
	  while W.rootIndices[j] ^ x <> W.rootIndices[j] + W.parentN  do
	      j := j + 1;
	  od;
	  W.permref[j]:=x;
	  w := w * W.permgens[i];
      fi;
    od;
end;

#############################################################################
##
#F  WeylReflections( <W> )  . . . . . . . . . . . . . .  the reflections in W
##
## 'WeylReflections'  returns  a  list of reduced expressions in the standard
##  generators for the set of reflections in the the Weyl group <W>. The i-th
##  entry in this list is the reflection along the i-th root in <W>.roots.
##
WeylReflections := function( W )
    return List( W.permref, x->WeylWordPerm(W,x));
end;
  
#############################################################################
##
#F  WeylRightCosetRepresentatives( <W>, <I>, <J> )  . . . . . . . . . . . . . 
#F  . . . . . . . . . . . . . . .  distinguished right coset representatives
##
##  'WeylRightCosetRepresentatives'  returns a list of reduced words  in  the
##  Weyl group <W>  that are  distinguished right   coset representatives for
##  the right cosets W(I)/W(J) where  W(K) is the  subgroup  generated by the
##  simple reflections corresponding to the subset K of [1..n] and  n  is the
##  rank of <W>.
##
WeylRightCosetRepresentatives := function( W, I, J ) 
    local  s, a, r, i, j, e, ea, en, x, n;
    e := [ () ];
    n := 1;
    ea := [ () ];
    for i  in [ 1 .. W.N ]  do
      en := [  ];
      for a  in ea  do
        for r  in I  do
          if W.rootIndices[r] ^ ( a ^ -1 ) <= W.parentN   then
            x := a * W.permgens[r];
            if not x in en  then
              j := 1;
              while j<=Length(J) and W.rootIndices[J[j]]^x<=W.parentN do
                j := j + 1;
              od;
              if j > Length( J ) then
                Add( en, x );
                n := n + 1;
              fi;
            fi;
         fi;
       od;
     od;
     ea := en;
     Append( e, ea );
   od;
   Print( "#I  Number of cosets = ", n, "\n" );
   s := [  ];
   for i  in e  do
     Add( s, WeylWordPerm( W, i ) );
   od;
   return s;
end;

#############################################################################
##
#F  WeylCosetPermRepresentation( <W>, <J> ) . . . . . . . . . . . . . . . . .
#F  . . . .  permutation representation on the cosets of a parabolic subgroup
##
##  'WeylCosetPermRepresentation'  returns  the  list of permutations induced
##  by the standard generators of the Weyl group  <W>  on the cosets  of  the
##  parabolic subgroup generated by the elements in the set <J>. The   cosets 
##  are   in   the   order   determined   by   the  result  of  the  function
##  'WeylRightCosetRepresentatives( <W>, <I>, <J> )' where I={1,...,n}.
##
WeylCosetPermRepresentation := function( W, J )
    local  wjdj, n, D, i, per, per1, s, d, j;
    wjdj := function ( W, D, J, w )
          local  w, r, c;
          if w = ()  then
              return [ (), () ];
          else
              if w in D  then
                  return [ (), w ];
              fi;
              for r  in J  do
                  if W.rootIndices[r] ^ w > W.parentN  then
                      c := wjdj( W, D, J, W.permgens[r] * w );
                      return [ W.permgens[r] * c[1], c[2] ];
                  fi;
              od;
          fi;
      end;
    D := List( WeylRightCosetRepresentatives( W, [ 1 .. W.dim ], J ), 
                           function ( i ) return PermWeylWord( W, i ); end );
    n := Length( D );
    s := [  ];
    for i  in W.permgens  do
        per := [  ];
        for d  in D  do
            Add( per, wjdj( W, D, J, d * i ) );
        od;
        per1 := [  ];
        for d  in [ 1 .. n ]  do
            j := 1;
            while per[d][2] <> D[j]  do
                j := j + 1;
            od;
            per1[d] := j;
        od;
        Add( s, PermList( per1 ) );
    od;
    return s;
end;

#############################################################################
##
#F  WeylElements( <W> )  . . . . . . . . . . . . . . . . .  all elements of W
##
##  'WeylElements' returns a list of which the  i - th entry is the list   of
##  reduced words of length  i - 1 in the Weyl group <W>.
##
WeylElements := function( W )
    local  all, i, l, l1, w;
    l1 := Elements( W );
    l := List( [ 1 .. W.N + 1 ], function ( i )
            return [  ];
        end );
    Print( "#I  Order = ", Length( l1 ), "\n" );
    for i  in l1  do
        w := WeylWordPerm( W, i );
        Add( l[Length( w ) + 1], w );
    od;
    for i  in l  do
        Sort( i );
    od;
    return l;
end;

#############################################################################
##
#F  WeylConjugacyClasses( <W> )  . . . . . . . . . . . . .  conjugacy classes
##
##  'WeylConjugacyClasses' returns a list of representatives of the conjugacy
##  classes of the Weyl group  <W>.  Each  element in this list is given as a
##  word in the standard generators,  and  it has  the property that it is of
##  minimal length in its conjugacy class. 
##

WeylConjugacyClasses:=function(W)local t,p,res,rank,res,cl;
  W.name:=WeylName(W.type);
  res:=[];
  for t in W.type do
    rank:=Length(t[2]);
    p:=MappingPermListList([1..rank],t[2]);
    if t[1]="A" then  cl:=WeylConjugacyClassesA(rank);
    elif t[1]="B" then  cl:=WeylConjugacyClassesB(rank);
    elif t[1]="D" then  cl:=WeylConjugacyClassesD(rank);
    elif t[1]="E" then  cl:=WeylConjugacyClassesE(rank);
    elif t[1]="F" then  cl:=WeylConjugacyClassesF4();
    elif t[1]="G" then  cl:=WeylConjugacyClassesG2();
    elif t[1]="H" then  cl:=WeylConjugacyClassesH(rank);
    elif t[1]="I" then  cl:=WeylConjugacyClassesDihedral(t[3]);
    fi;
    cl:=List(cl,x->OnTuples(x,p));
    Add(res,cl);
  od;
  return List(Cartesian(res),x->Concatenation(x));
  end;

#############################################################################
##
#F  WeylConjugacyClasses( <W> )  . . . . . . . . . . . . .  conjugacy classes
##
##  'WeylConjugacyClasses' returns a list of representatives of the conjugacy
##  classes of the Weyl group  <W>.  Each  element in this list is given as a
##  word in the standard generators,  and  it has  the property that it is of
##  minimal length in its conjugacy class. 
##

OldWeylConjugacyClasses := function ( W )
    local   y,  bc,  f1,  wdJ,  z,  y,  rep,  nb,  a,  b,  r,  x,  m,
            c,  J,  dd,  i,  j,  s,  n;

    if W.cartan = CartanMat( "A", W.dim )  then
        return WeylConjugacyClassesA( W.dim );
    elif W.cartan = CartanMat( "B", W.dim )  then
        return WeylConjugacyClassesB( W.dim );
    fi;  
    wdJ := function ( w, J )
       local  ww, r, c;
       ww := Copy( w );
       if ww = ()  then
         return ();
       else
         while true  do
           r := 1;
           while r <= Length( J ) and W.rootIndices[J[r]]^ ww <= W.parentN  do
             r := r + 1;
           od;
           if r > Length( J )  then
             return ww;
           else
             ww := W.permgens[J[r]] * ww;
           fi;
         od;
       fi;
    end;
    f1 := function ( w )
      local  bahn, s, j, y, yy;
      bahn := [ w ];
      for j  in bahn  do
        for s  in [ 1 .. W.dim ]  do
          y := W.permgens[s] * j;
          if W.rootIndices[s]^ j > W.parentN  then
            if W.rootIndices[s] ^ (y ^ (-1 * 1)) > W.parentN  then
              return true;
            else
              yy := y * W.permgens[s];
              if not yy in bahn  then
                Add( bahn, yy );
              fi;
            fi;
          else
            yy := y * W.permgens[s];
                      if W.rootIndices[s]^ (y ^ (-1 * 1)) > W.parentN and not yy in bahn  then
                          Add( bahn, yy );
                      fi;
                  fi;
              od;
          od;
          return false;
      end;
    dd := [ () ];
    J := [ 1 .. W.dim ];
    n := W.dim;
    while n > Int( W.dim / 2 ) -1  do
        Unbind( J[n] );
        c := WeylRightCosetRepresentatives( W, Concatenation( J, [n] ), J );
        m := Length( dd );
        for x  in c  do
            i := PermWeylWord( W, x );
            if i <> ()  then
                for j  in [ 1 .. m ]  do
                    Add( dd, i * dd[j] );
                od;
            fi;
        od;
        for s  in J  do
       dd := Filtered(dd,i->WeylLengthPerm( W, wdJ( i * W.permgens[s], J ) ) 
                      >= WeylLengthPerm( W, i ) );
     od;
     n := n - 1;
   od;
   a := [ Elements( Group( W.permgens{J}, () ) ), dd ];
   b := [  ];
   Print( "#I  ", Length(a[1])*Length(a[2]), " elements to consider\n" );
   for i  in a[1]  do
     for j  in a[2]  do
       x := i * j;
       if not f1( x )  then
         Add( b, x );
       fi;
     od;
   od;
   Print( "#I  Still ", Length( b ), " elements to consider\n" );
   bc := List( b, i -> OrderPerm( i ) );
   rep := [ 1 .. Length( b ) ];
   nb := [  ];
   while rep <> [  ]  do
     x := [  ];
     for i  in rep  do
       if  bc[rep[1]] = bc[i]  
           and WeylLengthPerm(W,b[rep[1]]) = WeylLengthPerm(W,b[i])
           and IsConjugate( W, b[rep[1]], b[i] ) then
                Add( x, i );
        fi;
      od;
      SubtractSet( rep, x );
      y := List( x, i -> WeylLengthPerm( W, b[i] ) );
      m := Minimum( y );
      y := Filtered( [ 1 .. Length( x ) ], i -> y[i] = m );
      y := List( y, i -> WeylWordPerm( W, b[x[i]] ) );
      Sort( y, function ( x, y )
              if Length( x ) <> Length( y )  then
                return Length( x ) < Length( y );
              else
                return x <= y;
              fi;
          end );
      Add( nb, y[1] );
    od;
    Print( "#I  ", Length( nb ), " conjugacy classes found \n" );
    Sort( nb, function ( x, y )
          if Length( x ) <> Length( y )  then
            return Length( x ) < Length( y );
          else
            return x <= y;
          fi;
      end );
    return nb;
end;


#############################################################################
##
#F  Bruhat( <W>, <y>, <w> )  . . . . . . . . . . . . . . Bruhat partial order
##
##  'Bruhat'  returns true, if the element  <y>  is less than or equal to the 
##  element <w> of the Weyl group <W>, and false otherwise. Both <y> and  <w>
##  must be given as permutations on the root vectors of <W>.
##
Bruhat := function ( W, y, w )
    local  s, lw, ly;
    if y = ()  then 
       return  true;
    elif w = ()  then 
       return  y = w;
    else
       s := 1;
       while W.rootIndices[s] ^ w <= W.parentN  do
           s := s + 1;
       od;
       if W.rootIndices[s] ^ y > W.parentN  then
           return Bruhat( W, W.permgens[s] * y, W.permgens[s] * w );
       else
           return Bruhat( W, y, W.permgens[s] * w );
       fi;
    fi;
end;
   
#############################################################################
##
#F  PermWeylMat( <W> , <mat> ) . . . . . . . . . Find the permutation of
#F  the roots effected by an element of the Weyl group W given in
#F  matrix form
##  
## 
PermWeylMat:=function(W,mat)local r,res;
  res:=[];
  for r in W.roots do Add(res,Position(W.roots,r*mat));od;
  return PermList(res);
  end;

#############################################################################
##
#F  MatWeylPerm( <W> , <w> ) . . . . . . . . . Find the matrix in the
#F  reflection representation for an element of the Weyl group W given as
#F  a permutation.
##  
## 
MatWeylPerm:=function(W,p)local x,res;
# return W.roots{OnTuples([1..W.dim],p^W.weylFusion)};
  return W.roots{List([1..W.dim],x->x^(p^W.weylFusion))};
  end;

#############################################################################
##
#F  PrintDynkinDiagram( <type>) . . . . . . . . .
#F  Prints the Dynkin diagram of some Coxeter group
##  
##  This is a purely descriptive routine which shows the labeling of
##  the Dynkin diagram corresponding to the given data. The data should
##  conform to the format of W.type for a Weyl group; that is a list
##  of triples: a string (type), a permutation of [1..rank],
##  and for type I2(n) an additional argument specifying n.
##  Call as:
##                PrintDynkinDiagram(W.type);
##
PrintDynkinDiagram:=function(t)local i,type,v,r,n,t,a;
  for a in t do
  type:=a[1];v:=a[2];r:=Length(v);
  if Length(a)=3 then Print("I2(",a[3],")    ",v[1]," - ",v[2],"\n");
  elif type="B" then
    Print("B",r,"    ",v[1]," < ",v[2]);
    for i in [3..r] do Print(" - ",v[i]);od;Print("\n");
  elif type="C" then
    Print("C",r,"    ",v[1]," > ",v[2]);
    for i in [3..r] do Print(" - ",v[i]);od;Print("\n");
  elif type="A" then
    Print("A",r,"    ",v[1]);
    for i in [2..r] do Print(" - ",v[i]);od;Print("\n");
  elif type="D" then
    Print("D",String(r,2),"      ",v[2],"\n");
    Print("         |\n");
    Print("     ",v[1]);for i in [3..r] do Print(" - ",v[i]);od;Print("\n");
  elif type="E" then
    Print("E",r,"          ",v[2],"             \n");
    Print("            |\n");
    Print("    ",v[1]);for i in [3..r] do Print(" - ",v[i]);od;Print("\n");
  elif type="G" then
    Print("G2  ",v[1]," < ",v[2]," \n");
  elif type="F" then
    Print("F4    ",v[1]," - ",v[2]," > ",v[3]," - ",v[4],"\n");
  elif type="H" then
    Print("        5 \n");
    Print("H",r,"    ",v[1]," - ",v[2]," - ",v[3]);
    if r=4 then Print(" - ",v[4]);fi;
    Print("\n");
  fi;
  od;
  end;

#############################################################################
##
#F  CartanCoeff( <a>, <b> )  . . . given alpha=W.roots[a] and beta=W.roots[b],
#F                compute the coefficient 2<alpha, beta>/<alpha, alpha>
##
##  This allow to find the geometry (the Cartan matrix) of any subset of roots
##  of W.
##
CartanCoeff:=function(W,a,b)local p;
  p:=PositionProperty(W.roots[a],x->x<>0);
  return (W.roots[b][p]-W.roots[b^(W.permref[a]^W.weylFusion)][p])
	    /W.roots[a][p];
end;

#############################################################################
##
#F  WeylSubgroup( <W>, <r> )  . . . given a Weyl group W and a set of roots
#F   given as their index in the roots of W, return the Weyl subgroup
#F   generated by the reflections w.r.t. these roots.
##
##
##  A Weyl Subgroup is a permutation subgroup, and otherwise has the same
##  fields as a Weyl group. The difference is that the following fields
##  contain useful values:
##
##  rootIndices: the index of the roots in the parent's roots
##  weylFusion: a permutation such that 
##              OnTuples(W.rootIndices,W.weylFusion)=[1..W.degree]
##  parentN:  the number of positive roots of the parent
##
WeylSubgroup:=function(W,genref)
  local H,simpleroots,permrefs,i,j,s,t,r,gens,p,q;
  gens:=1+List(genref-1,i->i mod W.N);
  if ForAll(genref, x->x<=W.dim) then simpleroots:=genref;
    gens:=W.permgens{simpleroots};
    H:=Subgroup(W,gens);
  else 
# 
# For the algorithm used, see Dyer's thesis.
# 
    permrefs:=W.permref;
    H:=Subgroup(W,permrefs{gens});
    simpleroots:=Filtered(Concatenation(Orbits(H,W.rootIndices{genref})),
                 y->y<=W.N);
    simpleroots:=Filtered(simpleroots,j->
            Number(simpleroots,i->i^permrefs[j]>W.N)=1);
    gens:=permrefs{simpleroots};
  fi;
  if IsIdentical(H,W) then 
    H:=Copy(H); H.parent:=W;
  fi;
  H.cartan:=[];
  for i in [1..Length(simpleroots)] do
    H.cartan[i]:=[];
    for j in [1..Length(simpleroots)] do
      H.cartan[i][j]:=CartanCoeff(W,simpleroots[i],simpleroots[j]);
    od;
  od;
  H.isWeylGroup := true;
  H.dim:=Length(simpleroots);
  H.matgens := SimpleReflectionMatrices(H.cartan);
  H.permgens := gens;
  H.parentN:=W.N;
  r:=Rootsystem(H.cartan);
  H.N:=Length(r);
  H.degree:=2*H.N;
  H.roots := Concatenation(r,-r);
  r:=H.roots*W.roots{simpleroots};
  H.rootIndices:=ListBlist([1..W.degree],BlistList(W.roots,r));
  H.rootIndices:=
       Permuted(H.rootIndices,Sortex(W.roots{H.rootIndices})*Sortex(r)^-1);
  H.weylFusion:=MappingPermListList(H.rootIndices,[1..H.degree]);
  if not ForAll(W.parameter,x->x=1) then
     if not ForAll(H.rootIndices{[1..H.dim]},x->x in [1..W.dim]) then
       Error(Concatenation(
             "cannot handle Hecke subalgebras generated by non-simple roots\n",
             "try setting the parameters of the parent to 1 (e.g. Hecke(W,1))")
            );
     else
       H.parameter:=W.parameter{H.rootIndices{[1..H.dim]}};
     fi;
  else
     H.parameter := List( [ 1 .. H.dim ], i -> 1 );
  fi;
  H.type:=CartanType(H.cartan);
  H.permref:=WeylReflectionsPerm(H);
  H.operations:=WeylGroupOps;
  return H;
end;

#############################################################################
##
#F  ReducedInWeylCoset( <W>, <w> ) . . . given w, an automorphism of W given as
#F     a permutation of its roots,  or an element of the parent group of W,
#F    return the element in the right coset wW which sends all simple roots
#F    of W to positive roots.
##
ReducedInWeylCoset:=function(W,w)local res,v,s;
   v:=w;
   while WeylLengthPerm(W,v)>0 do
     for s in W.permgens do
      if WeylLengthPerm(W,v*s)<WeylLengthPerm(W,v) then
       v:=v*s;
      fi;
     od;
   od;
   return v;
end;

#############################################################################
##
#F  ReflectionChar( <W>, <w> ) ... The reflection character on w
##
##  'ReflectionChar' returns  the value of the reflection character of W at
##  the element w.

ReflectionChar:=function(W,p)
  return Sum([1..W.dim],s->W.roots[s^(p^W.weylFusion)][s]);
end;

########################################################################
##
#F  HighestShortRoot( <W> ) . . . . . . . . . . . . . highest short root 
##
##  'HighestShortRoot' computes the unique short root of maximal height.
##  Note that if all roots have the same lenght then this is the  unique
##  root of maximal height. The latter can be extracted by  W.roots[W.N]
##
HighestShortRoot:=function(W)
  local l,s,r,i,k,q;
  if Length(W.type)>1 then
    Print("#W root system is not indecomposable \n");
    return false;
  fi;
  l:=Combinations([1..W.dim],2);
  for k in l do
    if W.cartan[k[1]][k[2]]<>0 then
      q:=W.cartan[k[1]][k[2]]/W.cartan[k[2]][k[1]];
      if q<>1 then 
        if q>1 then
          r:=Orbit(W,k[1]);
        else 
          r:=Orbit(W,k[2]);
        fi;
        s:=List(r,i->Sum(W.roots[i]));
        return W.roots[r[Position(s,Maximum(s))]];
      fi;
    fi;
  od;
  return W.roots[W.N];
end;

#############################################################################
##
#F  ParametersCentralizers( <W> ) . . parameters for the conjugacy classes of
##  centralizers  of  semisimple elements  in a  Chevalley group with a given
##  Weyl group
##
##  'ParametersCentralizers' returns a  list   of   pairs  which parametrizes 
##  the  conjugacy classes  of  centralizers  of  semisimple  elements  in  a 
##  Chevalley group with Weyl group <W>. Each pair consists of a subset which
##  contains fundamental roots or the highest root in the root system of <W>,
##  and an element in the normalizer in <W> of the subgroup generated by  the
##  reflections along the roots in the given subset. 
## 
ParametersCentralizers := function ( W )
    local  closure, x, n, t, te, i, j, k, cp;
    closure := function ( W, l )
          local  i, j, n, r, nn, x, p;
          if l = [  ]  then
              return [  ];
          fi;
          r := W.roots;
          n := List( l, function ( i )
                  return r[i];
              end );
          n := Set( Concatenation( n, -n ) );
          nn := Copy( l );
          for i  in l  do
              if i > W.N  then
                  Add( nn, i - W.N );
              else
                  Add( nn, i + W.N );
              fi;
          od;
          for i  in n  do
              for j  in n  do
                  x := i + j;
                  p := Position( r, x );
                  if p <> false and not x in n  then
                      Add( n, x );
                      Add( nn, p );
                      if not -1 * x in n  then
                          Add( n, -1 * x );
                          if p > W.N  then
                              Add( nn, p - W.N );
                          else
                              Add( nn, p + W.N );
                          fi;
                      fi;
                  fi;
              od;
          od;
          Sort( nn );
          return nn;
      end;
    te := Combinations( [ 1 .. W.dim ] );
    Sort( te, function ( i, j )
          if Length( i ) = Length( j )  then
              return i <= j;
          else
              return Length( i ) < Length( j );
          fi;
      end );
    t := [  ];
    n := 0;
    for i  in te  do
        if ForAll( t, function ( j )
                  return Size( i ) <> Size( j ) 
                    or PermGroupOps.RepresentativeSet( W, j, i ) = false;
              end )  then
            Add( t, i );
            n := n + 1;
        fi;
    od;
    for i  in te  do
        if i <> [ 1 .. W.dim ]  then
            j := Concatenation( i, [ W.degree ] );
            if ForAll( t, function ( k )
                      return Size( k ) <> Size( j ) 
                        or PermGroupOps.RepresentativeSet( W, j, k ) = false;
                  end )  then
                Add( t, j );
                n := n + 1;
            fi;
        fi;
    od;
    Sort( t, function ( i, j )
          return IsSubset( closure( W, j ), closure( W, i ) );
      end );
    Print( "#I  Number of Psi's = ", n, "\c" );
    n := 0;
    cp := [  ];
    for i  in t  do
        x := ConjugacyClasses( PermGroupOps.StabilizerSet( W, i ) );
        n := n + Length( x );
        Print( "\n#I  ", Length( x ), " conjugacy classes " );
        Add( cp, List( x, function ( j )
                return WeylWordPerm( W, j.representative );
            end ) );
    od;
    Print( "\n", "#I  Number of Pairs = ", n, "" );
    t := List( [ 1 .. Length( t ) ], function ( i )
            return List( cp[i], function ( j )
                    return [ t[i], j ];
                end );
        end );
    Print( "\n" );
    return Concatenation( t );
end;
