#############################################################################
##
#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, their
##  root  systems,  and  some special functions to  compute Kazhdan - Lusztig
##  polynomials and left cells.  
##

#############################################################################
##
#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 );'
##
CartanMat := function ( t, n )
    local  An, r, a, i, j, sr, cos, sin;
    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 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 = 6  then
            return List( [ 1 .. 6 ], function ( i )
                    return List( [ 1 .. 6 ], function ( j )
                            return a[i][j];
                        end );
                end );
        elif n = 7  then
            return List( [ 1 .. 7 ], function ( i )
                    return List( [ 1 .. 7 ], function ( j )
                            return a[i][j];
                        end );
                end );
        elif n = 8  then
            return a;
        else
            return "funny rank";
        fi;
    elif t = "H"  then
        r := (1 + ER( 5 )) / 2;
        if n = 3  then
            return [ [ 2, -1 * r, 0 ], [ -1 * r, 2, -1 ], [ 0, -1, 2 ] ];
        elif n = 4  then
            return 
             [ [ 2, -1 * r, 0, 0 ], [ -1 * r, 2, -1, 0 ], [ 0, -1, 2, -1 ], 
                [ 0, 0, -1, 2 ] ];
        else
            return "funny rank";
        fi;
    elif t = "I2" then
        cos := ( E( 2 * n ) + E( 2 * n ) ^ -1 ) / 2;
        sin := - E( 4 ) * ( E( 2 * n ) - E( 2 * n ) ^ -1 ) / 2; 
        sr := [ [ 1, 0 ], [ - cos, sin ] ]; 
        a := [ ];
        for i  in [ 1 .. 2 ]  do 
            a[ i ] := [ ]; 
            for j  in [ 1 .. 2 ]  do 
                a[ i ][ j ] := 2 * (sr[ i ] * sr[ j ]) / (sr[ i ] * sr[ i ]);
            od;
        od;
        return  a;
    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, i, j, n, m;
    n := Length( a );
    m := Length( b );
    c := [  ];
    for i  in [ 1 .. n + m ]  do
        c[i] := [  ];
        for j  in [ 1 .. n + m ]  do
            if i <= n and j <= n  then
                c[i][j] := a[i][j];
            elif i > n and j > n  then
                c[i][j] := b[i - n][j - n];
            else
                c[i][j] := 0;
            fi;
        od;
    od;
    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, i, k, a, e, g;
    l := Length( A );
    e := A ^ 0;
    g := [  ];
    for k  in [ 1 .. l ]  do
        a := [  ];
        for i  in [ 1 .. l ]  do
            a[i] := Copy( e[i] );
            a[i][k] := e[i][k] - A[k] * e[i];
        od;
        g[k] := a;
    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, neua;
    R := A ^ 0;
    if ForAll( Iterated( A, Concatenation ), IsInt )  then
        for a  in R  do
            v := A * a;
            for i  in [ 1 .. Length( A ) ]  do
                r := 1;
                while a - r * R[i] in R  do
                    r := r + 1;
                od;
                if r - v[i] > 1  then
                    neua := a + R[i];
                    if not neua in R  then
                        Add( R, neua );
                    fi;
                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
                        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  i, j, p, l, R, perm, n;
    R := Concatenation( r, -1 * r );
    n := Length( R );
    p := [  ];
    for i  in [ 1 .. Length( mats ) ]  do
        perm := r * mats[i];
        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;

WeylGroupOps := rec( Print := function( W )  Print( "Weyl( ",
             W.cartan," )" );  end, Size := function( W )  return 
             Size( Group( W.permgens, () ) ); end );

#############################################################################
##
#F  Weyl( <mat> )  . . . . . . . . . . . . . . . . create a record containing
#F  data on  roots and  the Weyl group determined by the Cartan matrix <mat>.
##  Typically, <mat> can be taken as the result of the function 'CartanMat'.
##
##  'Weyl' returns a record with the following entries: 
##  cartan    : the Cartan matrix <mat>,
##  dim       : the size of <mat>
##  degree    : the number of positive roots
##  N         : the total number of roots
##  roots     : the root vectors
##  matgens   : the matrices of the simple reflections
##  permgens  : the permutations on the root vectors
##  parameter : [1 .. 1]
##
Weyl := function( A )
    local  r, m;
    if not IsMat( A ) then
       return "funny input";
    fi;
    r := Rootsystem( A );
    m := SimpleReflectionMatrices( A );
    return rec( isDomain := true, isWeylGroup := true,
         cartan := A, dim := Length( A ), degree := 2 * Length( r ),
         N := Length( r ), roots := r, matgens := m,
         permgens := PermRepresentationRoots( m, r ), 
         parameter := List( [ 1 .. Length( A ) ], i -> 1 ),
         operations := WeylGroupOps );
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  r, l;
    l := 0;
    for r  in [ 1 .. W.N ]  do
        if r ^ w > W.N  then
            l := l + 1;
        fi;
    od;
    return l;
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  l, i, a;
    l := Length( w );
    if l = 0  then
        return ();
    elif l = 1  then
        return W.permgens[w[1]];
    else
        a := W.permgens[w[1]];
        for i  in [ 2 .. l ]  do
            a := a * W.permgens[w[i]];
        od;
        return a;
    fi;
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 i ^ ww <= W.N  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 )
    local  i, l, ww;
    if Length( w ) <= 1  then
        return w;
    else
        ww := ();
        for i  in w  do
            ww := ww * W.permgens[i];
        od;
        l := [  ];
        while ww <> ()  do
            i := 1;
            while i ^ ww <= W.N  do
                i := i + 1;
            od;
            Add( l, i );
            ww := W.permgens[i] * ww;
        od;
        return l;
    fi;
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 i ^ ( w ^ -1 ) > W.N 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  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 )
    local  p, l, r, s, i, j, x;
    r := Copy( W.permgens );
    for i  in r  do
        for s  in W.permgens  do
            x := s * i * s;
            if not x in r  then
                Add( r, x );
            fi;
        od;
    od;
    p := [  ];
    for i  in r  do
        j := 1;
        while j ^ i <> j + W.N  do
            j := j + 1;
        od;
        Add( p, j );
    od;
    l := [  ];
    for i  in [ 1 .. W.N ]  do
        l[ p[ i ] ] := WeylWordPerm( W, r[ i ] );
    od;
    return l;
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 r ^ ( a ^ -1 ) <= W.N   then
                    x := a * W.permgens[r];
                    if not x in en  then
                        j := 1;
                        while j <= Length( J ) and J[ j ] ^ x <= W.N 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 r ^ w > W.N  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( Group( W.permgens, () ) );
    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. 
##
WeylConjugacyClassesA := function ( n )
    local  part, p, i, z, l, reps, w;

    part := Partitions( n + 1 );
    reps := [  ];
    for p  in part  do
        w := [  ];
        i := 1;
        for l  in p  do
            z := 1;
            while z < l  do
                Add( w, i );
                z := z + 1;
                i := i + 1;
            od;
            i := i + 1;
        od;
        Add( reps, w );
    od;
    return reps;
end;

WeylConjugacyClassesB := function ( n )
   local  z, l, part, i, j, p, reps, w;

   part:= PartitionTuples(n, 2);
   reps:= [];
   for p  in part  do
      p[2]:= Reversed(p[2]);
      w:= [];
      i:= 1;
      for l  in p[2]  do
          for j  in [1..i]  do
              Add(w, i + 1 - j);
          od;
          for j  in [2..i]  do
              Add(w, j);
          od;
          i:= i + 1;
          z:= 1;
          while z < l  do
              Add(w, i);
              z:= z + 1;
              i:= i + 1;
          od;
      od;
      i:= i + 1;
      for l  in p[1]  do
          z:= 1;
          while z < l  do
              Add(w, i);
              z:= z + 1;
              i:= i + 1;
          od;
          i:= i + 1;
      od;
      Add(reps, w);
   od;
   return reps;
end;

WeylConjugacyClasses := function ( W )
    local   y,  bc,  f1,  wdJ,  G,  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 J[r] ^ ww <= W.N  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 s ^ j > W.N  then
                      if s ^ (y ^ (-1 * 1)) > W.N  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 s ^ (y ^ (-1 * 1)) > W.N 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( List( J, function ( i )
                    return W.permgens[i];
                end ), () ) ), 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" );
    G := Group( W.permgens, () );
    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( G, 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 s ^ w <= W.N  do
           s := s + 1;
       od;
       if s ^ y > W.N  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  KazhdanLusztigPolynomial( <W>, <y>, <w>, <u> ) Kazhdan-Lusztig polynomial
##
##  'KazhdanLustzigPolynomial'   returns   the  Kazhdan - Lusztig  polynomial 
##  in the indeterminate <u> corresponding to the elments <y> and <w>  (given 
##  as reduced expressions for <y> and <w> of the Weyl group <W>. 
##
KazhdanLusztigPolynomial := function ( W, y, w, u )
    local  Mue, kl, l;
      Mue := function ( W, x, y, lx, ly )
          local  i;
          if ly = lx + 1  then
              return 1;
          fi; 
          for i  in W.permgens  do
              if i ^ x <= W.N  and  i ^ y > W.N  then
                  return 0;
              fi;
          od;
          return -1;
      end;
    kl := function ( W, y, w )
          local  ly, lv, lw, i, j, x, s, v, c, z, m, p, w0;
          if not Bruhat( W, y, w )  then
              return  0 * u;
          fi;
          ly := WeylLengthPerm( W, y );
          lw := WeylLengthPerm( W, w );
          if lw - ly <= 2  then
              return  u ^ 0;
          fi; 
          s := 1;
          while s <= W.dim and (s ^ y > W.N or s ^ w <= W.N)  do
              s := s + 1;
          od;
          if s <= W.dim  then
              return kl( W, W.permgens[s] * y, w );
          fi;
          if not IsBound( W.all )  then
              W.all := List( WeylElements( W ), i -> 
                                 List( i, j -> PermWeylWord( W, j ) ) );
          fi;
          s := 1;
          while s ^ w <= W.N  do
              s := s + 1;
          od;
          v := W.permgens[s] * w;
          lv := lw - 1;
          p := kl( W, W.permgens[s] * y, v ) + u * kl( W, y, v );
          i := lv - 1;
          while lw - i <= 2 * Degree( p )  do
              j := 1;
              while
               j <= Length( W.all[i + 1] ) and lw - i <= 2 * Degree ( p )  and 
                    p.coefficients[(lw - i) / 2 + 1] > 0  do
                  z := W.all[i + 1][j];
                  if s ^ z > W.N  and Bruhat( W, y, z ) and 
                                              Bruhat( W, z, v )  then
                      m := Mue( W, z, v, i, lv );
                      if m = -1  then
                          x := kl( W, z, v );
                          if Degree( x ) * 2 = lv - i - 1  then
                              m := x.coefficients[(lv - i + 1) / 2];
                          else
                              m := 0;
                          fi;
                      fi;
                      if m <> 0  then
                          p := p - m * u ^ ( ( lv - i + 1)/2 ) * kl( W, y, z );
                      fi;
                  fi;
                  j := j + 1;
              od;
              i := i - 2;
          od;
          return p;
      end;
    return kl( W, PermWeylWord( W, y ), PermWeylWord( W, w ) );
end;

#############################################################################
##
#F  KLCoefficient( <W>, <y>, <w>, k )  . . . . . . . . . . . . . . . . . . .
#F  . . . . . . . . . .  the k-th coefficient of a Kazhdan-Lusztig polynomial
##
KLCoefficient := function ( W, y, w, k )
    local  Mue, m, lw, ly, s, sy, v, lv, z, i, j, m1, lg;
     Mue := function ( W, x, y, lx, ly )
          local  i;
          if ly = lx + 1  then
              return 1;
          fi; 
          for i  in W.permgens  do
              if i ^ x <= W.N  and  i ^ y > W.N  then
                  return 0;
              fi;
          od;
          return -1;
      end;
    lw := WeylLengthPerm( W, w );
    ly := WeylLengthPerm( W, y );
    if k < 0 or not IsInt( k ) or ly >= lw or 2 * k > lw - ly - 1 
               or not Bruhat( W, y, w )  then
        return 0;
    elif k = 0  then
        return 1;
    fi;
    s := 1;
    while s <= W.dim and (s ^ y > W.N or s ^ w <= W.N)  do
        s := s + 1;
    od;
    if s <= W.dim  then
        return KLCoefficient( W, W.permgens[s] * y, w, k );
    fi;
    s := 1;
    while s ^ w <= W.N  do
        s := s + 1;
    od;
    v := W.permgens[s] * w;
    lv := lw - 1;
    sy := W.permgens[s] * y;
    m := KLCoefficient( W, sy, v, k ) + KLCoefficient( W, y, v, (k - 1) );
    if m = 0  then
        return 0;
    fi;
    if not IsBound( W.all )  then
        W.all := List( WeylElements( W ), i -> 
                                 List( i, j -> PermWeylWord( W, j ) ) );
    fi;
    i := lv - 1;
    while i >= ly  do
        lg := Length( W.all[i + 1] );
        j := 1;
        while j <= lg  do
            z := W.all[i + 1][j];
            if s ^ z > W.N and Bruhat( W, y, z ) and Bruhat( W, z, v )  then
                m1 := Mue( W, z, v, i, lv );
                if m1 = -1  then
                    m1 := KLCoefficient( W, z, v, (lv - i - 1) / 2 );
                fi;
                if m1 > 0  then
                    m := m - m1 * KLCoefficient( W, y, z, 
                                                 k - (lv - i + 1) / 2 );
                    if m = 0  then
                        return 0;
                    fi;
                fi;
            fi;
            j := j + 1;
        od;
        i := i - 2;
    od;
    return m;
end;

#############################################################################
##
#F  WeylMueMat( <W>, <list> )  . . . .  the matrix of leading coefficients of
#F  Kazhdan-Lusztig polynomials  of elements in a given list of reduced words
##
##  The elements in <list> must be ordered by increasing length. E.g., <list>
##  is the result of 'Iterated(WeylElements(<W>),Concatenation)'.
##
WeylMueMat := function ( W, c )
    local  m, i, j, k, ll, x, li, lj, n, p, d, s, w0;
    w0 := PermWeylWord( W, LongestWeylWord( W ) );
    if not IsBound( W.all )  then
        W.all := List( WeylElements( W ), i -> 
                                 List( i, j -> PermWeylWord( W, j ) ) );
    fi;
    n := Length( c );
    ll := List( c, function ( i ) return PermWeylWord( W, i ); end );
    m := List( [ 1 .. n ], function ( i ) return [  ]; end );
    for k  in [ 1 .. Length( c[n] ) - Length( c[1] ) + 1 ]  do
        for i  in [ 1 .. n ]  do
            li := ll[i];
            for j  in [ 1 .. i ]  do
                lj := ll[j];
                d := Length( c[i] ) - Length( c[j] );
                if d = k - 1  then
                     if Length( c[i] ) + Length( c[j] ) > W.N  then
                          m[i][j] := KLCoefficient( W, w0 * li, w0 * lj, 
                                                         ( d - 1 ) / 2  );
                     else
                          m[i][j] := KLCoefficient( W, lj, li, ( d - 1) / 2 );
                     fi;
                fi;
            od;
        od;
    od;
    return  m;
end;

#############################################################################
##
#F  DecomposedLeftCells( <W>, <c>, <mue> )  . . . . decompose into left cells
##
##  'DecomposedLeftCells' returns a list of pairs.The first component of each
##  pair  consists of those elements in the list <c>  of reduced words in the 
##  Weyl group  <W>  which  lie  in one  left cell  C, the  second  component 
##  consists  of the  corresponding  matrix of highest coefficients mue(y,w),
##  where  y, w  are in C. We require that the  elements in  <c>  are reduced
##  words in the standard generators ordered by incresing lengths.
##
DecomposedLeftCells := function ( W, c, mue )
    local  Lleq, s1, s2, d, ce, i, j, m, rest;
    Lleq := function ( W, mue, x, y )
          local  i, ex, ey;
          if x = y  then
              return false;
          elif x < y and mue[y][x] = 0  then
              return false;
          elif x > y and mue[x][y] = 0  then
              return false;
          else
              ex := PermWeylWord( W, c[x] );
              ey := PermWeylWord( W, c[y] );
              for i  in [ 1 .. W.dim ]  do
                  if i ^ ex > W.N  and  i ^ ey <= W.N  then
                      return true;
                  fi;
              od;
              return false;
          fi;
      end;
    ce := [  ];
    rest := [ 1 .. Length( c ) ];
    while rest <> [  ]  do
        s1 := [ rest[1] ];
        for j  in s1  do
            for i  in [ 1 .. Length( c ) ]  do
                if not i in s1 and Lleq( W, mue, j, i )  then
                    Add( s1, i );
                fi;
            od;
        od;
        Sort( s1 );
        s2 := [ rest[1] ];
        for j  in s2  do
            for i  in [ 1 .. Length( c ) ]  do
                if not i in s2 and Lleq( W, mue, i, j )  then
                    Add( s2, i );
                fi;
            od;
        od;
        Sort( s2 );
        d := Intersection( s1, s2 );
        Sort( d );
        m := [  ];
        for i  in [ 1 .. Length( d ) ]  do
            m[i] := [  ];
            for j  in [ 1 .. i ]  do
                m[i][j] := mue[d[i]][d[j]];
            od;
        od;
        Add( ce, [ List( d, function ( i ) return c[i]; end ), m ] );
        SubtractSet( rest, d );
    od;
    return ce;
end;

#############################################################################
##
#F  LeftCells( <W> )  . . . . . . . . . . . . . . . . . . the left cells of W
##
##  'LeftCells' returns a list of pairs. The first  component  of  each  pair
##  consists of the reduced words in the Weyl group <W> which lie in one left
##  cell C,  the  second  component  consists  of the corresponding matrix of 
##  highest coefficients mue(y,w), where y,w are in C.
##
LeftCells := function ( arg )
    local  W, Rw, t, i, j, max, ze, rw, x;
    Rw := function ( W, w )
          local  s, r, ww, l;
          l := Length( w );
          ww := PermWeylWord( W, w );
          r := [  ];
          for s  in [ 1 .. W.dim ]  do
              if WeylLengthPerm( W, ww * W.permgens[s] ) < l  then
                  Add( r, s );
              fi;
          od;
          return r;
      end;
    W := arg[1];
    t := Combinations( [ 1 .. W.dim ] );
    rw := [  ];
    for j  in [ 1 .. Length( t ) ]  do
        rw[j] := [  ];
    od;
    if not IsBound( W.all )  then
        W.all := WeylElements( W );
    fi;
    if  Length( arg ) = 2  then
        max := arg[2];
    else
        max := Sum( W.all, Length );
    fi;
    for j  in W.all  do
        for i  in j  do
            Add( rw[Position( t, Rw( W, i ) )], i );
        od;
    od;
    W.all := List( W.all, i -> List( i, j -> PermWeylWord( W, j ) ) );
    ze := [  ];
    for i  in rw  do
        if Length(i) <= max then
            Print( "#I  R(w) = ", t[Position( rw, i )],
                                     " :  Size = ", Length( i ) , "\c" );
            x := DecomposedLeftCells( W, i, WeylMueMat( W, i ) );
            if Length( x ) = 1  then
                Print( ",  ", Length( x ), " new cell \n" );
            else
                Print( ",  ", Length( x ), " new cells \n" );
            fi;
            Append( ze, x );
        fi;
    od;
    Unbind( W.all );
    return ze;
end;

#############################################################################
##
#F  LeftCellRepresentation( <W> , <cell>, <v> )  . . . . . . . representation
#F  of  the  Hecke  algebra  associated  to  a  left  cell  of  a  Weyl group
##
##  'LeftCellRepresentation' returns a list of matrices giving the left  cell 
##  representation of the  Hecke algebra with parameter <v>^2 associated with
##  the  Weyl group  <W>. The argument <cell> is a pair with  first component
##  a list of reduced  words which  form  a  left cell, and  second component
##  the corresponding matrix of mue's. Typically, <cell> is the result of the
##  function 'WeylMueMat'.
##
LeftCellRepresentation := function ( W, cell, v )
    local  s, i, j, rep, t, f1, f2, pc, c, mue, eins, null;
    eins := v ^ 0;
    null := 0 * v;
    c := cell[1];
    pc := List( c, function ( i )
            return PermWeylWord( W, i );
        end );
    mue := cell[2];
    rep := [  ];
    for s  in W.permgens  do
        f1 := [  ];
        for i  in pc  do
            if
             WeylLengthPerm( W, s * i ) > WeylLengthPerm( W, i )
                    then
                Add( f1, true );
            else
                Add( f1, false );
            fi;
        od;
        f2 := List( pc, function ( i )
                return Position( pc, s * i );
            end );
        t := [  ];
        for i  in [ 1 .. Length( c ) ]  do
            t[i] := [  ];
            for j  in [ 1 .. Length( c ) ]  do
                if i = j  then
                    if f1[i]  then
                        t[i][i] := v ^ 2;
                    else
                        t[i][i] := -1 * eins;
                    fi;
                elif i > j  then
                    if f2[j] = i  then
                        t[i][j] := v;
                    else
                        t[i][j] := null;
                    fi;
                else
                    if f1[j] and mue[j][i] <> 0 and not f1[i]  then
                        t[i][j] 
                         := Sum( List( [ 1 .. mue[j][i] ], function ( x )
                                    return eins;
                                end ) ) * v;
                    else
                        t[i][j] := null;
                    fi;
                fi;
            od;
        od;
        Add( rep, t );
    od;
    return rep;
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, W1, i, j, k, cp;
    closure := function ( W, l )
          local  i, j, n, r, nn, x, p;
          if l = [  ]  then
              return [  ];
          fi;
          r := Concatenation( W.roots, -1 * W.roots );
          n := List( l, function ( i )
                  return r[i];
              end );
          n := Set( Concatenation( n, -1 * 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 );
    W1 := Group( W.permgens, () );
    t := [  ];
    n := 0;
    for i  in te  do
        if ForAll( t, function ( j )
                  return Size( i ) <> Size( j ) 
                    or PermGroupOps.RepresentativeSet( W1, 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( W1, 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( W1, 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 Iterated( t, Concatenation );
end;

#############################################################################
##
#F  Hecke( <W>, <arg> )  . .  the Hecke algebra corresponding to a Weyl group 
##
##  'Hecke'  adds  to the record of the Weyl group  <W>  a function  T  which 
##  produces,  for each element  w in W, the corresponding basis element T(w)
##  in the Hecke algebra H associated to W.  These basis elements satisfy the
##  rules T(s)^2 = <q[s]> T(1) + ( <q[s]> - 1 )T(s),  for standard generators 
##  s, and  T(vw) = T(v) T(w), if l(vw) = l(v) + l(w). 
##  Typically, <W> is the result of the function 'Weyl'.
##  <arg> specifies the index paramters q[s]: if <arg> is a single number or
##  polynomial, then all parameters are assumed to be equal to this element.
##  Alternatively, <arg> may be a list of such numbers or polynomials.
##  Examples:  W:=Weyl(CartanMat("B",3)); q:=Indeterminate("q");
##             Hecke(W,1);         or          Hecke(W,[q,q^2,q^2]);
##  
Hecke:=function(W,arg)
  local i,j,ListOrder,HeckeAdd,HeckeProd1,HeckeProd,HeckeOps,q;
  if not IsList(arg) then
    q:=List([1..W.dim],i->arg);
  else
    q:=arg;
    for i in [1..W.dim] do
      for j in [i+1..W.dim] do
        if q[i]<>q[j] and 
                  not IsInt(OrderPerm(W.permgens[i]*W.permgens[j])/2) then
          Error("funny parameters");
        fi;
      od;
    od;
  fi;

  ListOrder:=function(l1,l2)
    if Length(l1)<>Length(l2) then
      return Length(l1)<Length(l2);
    else
      return l1<=l2;
    fi;
  end;

  HeckeAdd:=function(x,y)
    local a,z,i,j;
    z:=[]; i:=1; j:=1;
    if x[1]=[0*q[1],[]] then
      if y[1]<>[0*q[1],[]] then 
        i:=i+1;
      else 
        return [[0*q[1],[]]];
      fi;
    elif y[1]=[0*q[1],[]] then 
      j:=j+1;
    fi;
    while i<=Length(x) and j<=Length(y) do
      if ListOrder(x[i][2],y[j][2]) then
        if x[i][2]=y[j][2] then
          a:=x[i][1]+y[j][1];
          if a<>0*q[1] then
            Add(z,[a,x[i][2]]);
          fi;
          i:=i+1; j:=j+1;
        else
          Add(z,x[i]); i:=i+1;
        fi;
      else
        Add(z,y[j]); j:=j+1;
      fi;
    od;
    if i>Length(x) then
      Append(z,y{[j..Length(y)]});
    elif j>Length(y) then
      Append(z,x{[i..Length(x)]});
    fi;
    if z=[] then
      return  [[0*q[1],[]]];
    else
      return z;
    fi;
  end;

  HeckeProd1:=function(s,y)
    local  pjy,jy,z,s;
    z:=[[0*q[1],[]]];
    for jy in y do
      pjy:=PermWeylWord(W,jy[2]);
      if s^pjy<=W.N then
        z:=HeckeAdd(z,[[jy[1],WeylWordPerm(W,W.permgens[s]*pjy)]]);
      else
        if q[s]-1=0*q[s] then
          z:=HeckeAdd(z,[[jy[1],WeylWordPerm(W,W.permgens[s]*pjy)]]);
        else
          z:=HeckeAdd(z,[[q[s]*jy[1],WeylWordPerm(W,W.permgens[s]*pjy)],
                         [(q[s]-1)*jy[1],jy[2]]]);
        fi;
      fi;
    od;
    return z;
  end;

  HeckeProd:=function(x,y)
    local  ix,iz,z,s,i,n,neu;
    neu:=[[0*q[1],[]]];
    for ix in x do
      n:=Length(ix[2]);
      if n=0 then
        z:=[];
        for i in y do
          Add(z,[ix[1]*i[1],i[2]]);
        od;
      else
        z:=HeckeProd1(ix[2][n],y);
        for s in Reversed([1..n-1]) do
          z:=HeckeProd1(ix[2][s],z);
        od;
        for iz in z do
          iz[1]:=ix[1]*iz[1];
        od;
      fi;
      neu:=HeckeAdd(neu,z);
    od;
    return neu;
  end;

  HeckeOps:=rec(
      \+:=function(h1,h2)
         return rec(basrep:=HeckeAdd(h1.basrep,h2.basrep),
                    operations:=h1.operations);
      end,
     \-:=function(h1,h2)
        return h1+rec(basrep:=[[-q[1]^0,[]]],operations:=h1.operations)*h2;
      end,
     \*:=function(h1,h2)
        if not IsRec(h1) then
          h1:=rec(basrep:=[[h1*q[1]^0,[]]],operations:=h2.operations);
        else
          if IsPolynomial(h1) then
            h1:=rec(basrep:=[[h1,[]]],operations:=h2.operations);
          fi;
        fi;
        if h1.basrep[1][1]=0*q[1] then
          return rec(basrep:=[[0*q[1],[]]],operations:=h2.operations);
        elif h2.basrep[1][1]=0*q[1] then
          return rec(basrep:=[[0*q[1],[]]],operations:=h1.operations);
        else
          return rec(basrep:=HeckeProd(h1.basrep,h2.basrep),
                     operations:=h1.operations);
        fi;
     end,
    \^:=function(h,n)
       local i,p;
       if n<0 then
         if Length(h.basrep)=1 then
           if n=-1 then
             p:=[[h.basrep[1][1]^(-1),[]]];
             for i in Reversed(h.basrep[1][2]) do
               p:=HeckeProd(p,[[q[i]^(-1)-1,[]],[q[i]^(-1),[i]]]);
             od;
             return rec(basrep:=p,operations:=h.operations);
           else
             return (h^(-1))^(-n);
           fi;
         else
           Error("negative exponent");
         fi;
       elif n=0 then
         return rec(basrep:=[[q[1]^0,[]]],operations:=h.operations);
       else
         p:=h.basrep;
         for i in [2..n] do
           p:=HeckeProd(p,h.basrep);
         od;
         return rec(basrep:=p,operations:=h.operations);
       fi;
    end,

    Print:=function(h)
      local i,b;
      for i in [1..Length(h.basrep)] do
        b:=h.basrep[i];
        if i=1 then
          if IsRat(b[1]) then
            Print(b[1],"*T(",b[2],")");
          else
            Print("(",b[1],")*T(",b[2],")");
          fi;
        else
          if b[1]<>0*q[1] then
            if IsRat(b[1]) then
              if b[1]>0 then
                Print("+",b[1],"*T(",b[2],")");
              else
                Print(b[1],"*T(",b[2],")");
              fi;
            else
              Print("+(",b[1],")*T(",b[2],")");
            fi;
          fi;
        fi;
      od;
    end);
  W.T:=function(w)
         return rec(basrep:=[[q[1]^0,ReducedWeylWord(W,w)]],
                    operations:=HeckeOps);
       end;
  W.parameter:=q;
end;

#############################################################################
##
#F  HeckeReflectionRepresentation( <W> )  . . . . . . . . . . . . . . . . . . 
#F  . . . . . . . . . . . .  the reflection representation of a Hecke algebra
##
##  <W> must be a result of the function 'Weyl' 
##
HeckeReflectionRepresentation:=function(W)
  local t,i,j,n,m,a,C,ord,v;
  if Length(Set(W.parameter))>1 then
    Error("not yet implemted");
  fi;

  ord:=function(W,w)
  local pw,n,x;
    pw:=PermWeylWord(W,w);
    n:=1;
    x:=pw;
    while x<>()do
      x:=pw*x;
      n:=n+1;
    od;
    return n;
  end;

  v:=W.parameter[1];
  C:=[];
  for i in [1..W.dim] do
    C[i]:=[];
    for j in [1..i-1] do
      m:=ord(W,[i,j]);
      if m=2 then
        C[i][j]:=0*v;
        C[j][i]:=0*v;
      elif m=3 then
        C[i][j]:=v^0;
        C[j][i]:=v;
      elif m=4 then
        C[i][j]:=2*v^0;
        C[j][i]:=v;
      elif m=6 then
        C[i][j]:=3*v^0;
        C[j][i]:=v;
      fi;
    od;
    C[i][i]:=v+v^0;
  od;
  t:=[];
  for i in [1..W.dim] do
    a:=[];
    for j in [1..W.dim] do
      a[j]:=List([1..W.dim],function(k) return 0*v; end);
      a[j][j]:=v;
      a[j][i]:=a[j][i]-C[i][j];
    od;
    Add(t,a);
  od;
  return t;
end;
  
#############################################################################
##
#F  CheckHeckeDefiningRelations( <W> , <t> )  . . . . . . . . . . . . . check
#F  the  defining  relations  of  a  Hecke algebra for a given representation
##  
##  'CheckHeckeDefiningRelations' returns true or false, according to whether
##  a given  set  <t>  of  matrices  corresponding to the standard generators
##  of the  Weyl group  <W>  defines a representation of the associated Hecke 
##  algebra   or not. 
## 
CheckHeckeDefiningRelations:=function(W,t)
  local i,j,k,e,m,x,y,q;
  q:=W.parameter;
  for i in [1..W.dim] do
    e:=q[i]*t[1]^0;
    if not t[i]^2=e+(q[i]-1)*t[i] then
      Print("#I  Error in ",Ordinal(i)," quadradic relation\n");
      return false;
    fi;
    for j in [i+1..W.dim] do
      x:=W.permgens[i]*W.permgens[j];
      y:=x;
      m:=1;
      while y<>()do
        m:=m+1;
        y:=y*x;
      od;
      x:=t[i];
      for k in [2..m] do
        if k mod 2=0 then
          x:=x*t[j];
        else
          x:=x*t[i];
        fi;
      od;
      y:=t[j];
      for k in [2..m] do
        if k mod 2=0 then
          y:=y*t[i];
        else
          y:=y*t[j];
        fi;
      od;
      if not x=y then
        Print("#I  Error in (",i,",",j,") homogenous relation\n");
        return false;
      fi;
    od;
  od;
  return true;
end;
  
#############################################################################
##
#F  WeylClassPolynomials( <W>, <reps>, <w> )  . . . . . the class polynomials
##
##  'WeylClassPolynomials' returns  the  class polynomials of the element <w>
##  with  respect  to  representatives  <reps>  of  minimal  lengths  in  the
##  conjugacy  classes of the Weyl group <W>. Typically, <reps> is the result
##  of the  'WeylConjugacyClasses' and  <w>  is  a  reduced expression in the 
##  standard generators of <W>.

CycleType:=function(p,l)
  local i,x,c;
  c:=[];
  while l<>[] do
    x:=CyclePermInt(p,l[1]);
    Add(c,Length(x));
    SubtractSet(l,x);
  od;
  Sort(c);
  return c;
end;
  
RefChar:=function(W,p)
  local s,c,x;
  c:=0;
  for s in [1..W.dim] do
    x:=s^p;
    if x>W.N then
      c:=c-W.roots[x-W.N][s];
    else
      c:=c+W.roots[x][s];
    fi;
  od;
  return c;
end;

WeylClassPolynomials:=function(W,reps,elm)
  local orb,orb1,sim,G,alt,b,fertig,w,min,ll,l,class,pa,C,pr,pw,q,o,i,j,p,cl;

  if not IsBound(W.T) then
    Hecke(W,1);
  fi; 
  q:=W.parameter;

  orb:=function(w)
    local bahn,s,j,y,yy;
    bahn:=[PermWeylWord(W,w)];
    for j in bahn do
      for s in [1..W.dim] do
        y:=W.permgens[s]*j;
          if s^j>W.N then
            if s^(y^-1)>W.N then 
              return [WeylWordPerm(W,y*W.permgens[s]),WeylWordPerm(W,y),s];
            else
              yy:=y*W.permgens[s];
              if not yy in bahn then
                Add(bahn,yy);
              fi;
            fi;
          else
            yy:=y*W.permgens[s];
            if s^(y^-1)>W.N and not yy in bahn then
              Add(bahn,yy);
            fi;
          fi; 
       od;
     od;
    return false;
  end;

  orb1:=function(w)
    local bahn,s,j,y,yy;
    bahn:=[PermWeylWord(W,w)];
    for j in bahn do
      for s in [1..W.dim] do
        y:=W.permgens[s]*j;
          if s^j>W.N then
            if s^(y^-1)<=W.N then 
              yy:=y*W.permgens[s];
              if not yy in bahn then
                Add(bahn,yy);
              fi;
              if Length(bahn)>1000 then
                return bahn;
              fi;
            fi;
          else
            yy:=y*W.permgens[s];
            if s^(y^-1)>W.N and not yy in bahn then
              Add(bahn,yy);
            fi;
            if Length(bahn)>1000 then
              return bahn;
            fi;
          fi; 
       od;
     od;
    return bahn;
  end;
    
  sim:=function(a,pa,C,b)
    local i,c,pb;
    pb:=PermWeylWord(W,b);
    c:=PermGroupOps.RepresentativeConjugationElements(G,pb,pa);
    return c<>false and PermGroupOps.ElementProperty(C,i->
           WeylLengthPerm(W,c*i*pa)=WeylLengthPerm(W,c*i)+Length(a))<>false;
  end;

  G:=Group(W.permgens,());
  pr:=List(reps,i->PermWeylWord(W,i));
  class:=List([1..Length(reps)],i->[Length(reps[i]),RefChar(W,pr[i])]);

  alt:=W.T(elm);
  min:=0*W.T([]);
  #Print("#I  ");
  while alt<>0*W.T([]) do
    l:=Maximum(List(alt.basrep,i->Length(i[2])));
    ll:=Filtered(alt.basrep,i->Length(i[2])=l); 
  #  Print(Length(alt.basrep)," \c");
    for w in ll do
      o:=orb(w[2]);
      if o=false then
        pw:=PermWeylWord(W,w[2]);
        cl:=[Length(w[2]),RefChar(W,pw)];
        i:=1;
        while cl<>class[i] or not IsConjugate(G,pw,pr[i]) do
          i:=i+1;
        od;
        min:=min+w[1]*W.T(reps[i]);
      else
        i:=1;
        while i<=Length(alt.basrep) and 
                              Length(o[1])<>Length(alt.basrep[i][2]) do
          i:=i+1;
        od;
        if i>Length(alt.basrep) then
          alt:=alt+w[1]*q[o[3]]*W.T(o[1]);
        else
          fertig:=false;
          j:=i;
          while not fertig and j<=Length(alt.basrep) and 
                               Length(o[1])=Length(alt.basrep[j][2]) do
            if alt.basrep[j][2]=o[1] then
              alt:=alt+w[1]*q[o[3]]*W.T(alt.basrep[j][2]);
              fertig:=true;
            fi;
            j:=j+1;
          od;
          if not fertig then 
            b:=orb1(o[1]);
            j:=i;
            while not fertig and j<=Length(alt.basrep) and 
                                 Length(o[1])=Length(alt.basrep[j][2]) do
              if PermWeylWord(W,alt.basrep[j][2]) in b then
                alt:=alt+w[1]*q[o[3]]*W.T(alt.basrep[j][2]);
                fertig:=true;
              fi;
              j:=j+1;
            od;
          fi;
          if not fertig then
            j:=i;
            pa:=PermWeylWord(W,o[1]);
            C:=Centralizer(G,pa);
            while not fertig and j<=Length(alt.basrep) and 
                                 Length(o[1])=Length(alt.basrep[j][2]) do
              if sim(o[1],pa,C,alt.basrep[j][2]) then
                alt:=alt+w[1]*q[o[3]]*W.T(alt.basrep[j][2]);
                fertig:=true;
              fi;
              j:=j+1;
            od;
          fi;
          if not fertig then
            alt:=alt+w[1]*q[o[3]]*W.T(o[1]);
          fi;
        fi;
        i:=1;
        while i<=Length(alt.basrep) and 
                              Length(o[2])<>Length(alt.basrep[i][2]) do
          i:=i+1;
        od;
        if i>Length(alt.basrep) then
          alt:=alt+w[1]*(q[o[3]]-1)*W.T(o[2]);
        else
          fertig:=false;
          j:=i;
          while not fertig and j<=Length(alt.basrep) and 
                               Length(o[2])=Length(alt.basrep[j][2]) do
            if alt.basrep[j][2]=o[2] then
              alt:=alt+w[1]*(q[o[3]]-1)*W.T(alt.basrep[j][2]);
              fertig:=true;
            fi;
            j:=j+1;
          od;
          if not fertig then 
            b:=orb1(o[2]);
            j:=i;
            while not fertig and j<=Length(alt.basrep) and 
                                 Length(o[2])=Length(alt.basrep[j][2]) do
              if PermWeylWord(W,alt.basrep[j][2]) in b then
                alt:=alt+w[1]*(q[o[3]]-1)*W.T(alt.basrep[j][2]);
                fertig:=true;
              fi;
              j:=j+1;
            od;
          fi;
          if not fertig then
            j:=i;
            pa:=PermWeylWord(W,o[2]);
            C:=Centralizer(G,pa);
            while not fertig and j<=Length(alt.basrep) and 
                                 Length(o[2])=Length(alt.basrep[j][2]) do
              if sim(o[2],pa,C,alt.basrep[j][2]) then
                alt:=alt+w[1]*(q[o[3]]-1)*W.T(alt.basrep[j][2]);
                fertig:=true;
              fi;
              j:=j+1;
            od;
          fi;
          if not fertig then
            alt:=alt+w[1]*(q[o[3]]-1)*W.T(o[2]);
          fi;
        fi;
      fi;
      alt:=alt-w[1]*W.T(w[2]);
    od;
  od;

  p:=0*[1..Length(reps)];
  for i in min.basrep do
    p[Position(reps,i[2])]:=i[1];
  od;
  #Print("\n");
  return p;
end;

#############################################################################
##
#F  CharHeckeRepresentation( <rep> , <elts> ) . . . . . . . . . . . . . . . .
#F  . . . . . . . . . . . the character of a representation on given elements
##  
##  <rep> is a list of matrices, <elts> a list of reduced expressions. 
## 
CharHeckeRepresentation:=function(rep,elts)
  local char,e,i,c,n,x;
  char:=[];
  n:=Length(rep[1]);
  x:=First(Iterated(rep[1],Concatenation),i->i<>0*i);
  for e in elts do
    if e=[] then
      Add(char,n*x^0);
    else
      c:=rep[e[1]];
      for i in [2..Length(e)] do
        c:=c*rep[e[i]];
      od;
      Add(char,Sum([1..n],i->c[i][i]));
    fi;
  od;
  return char;
end;

#############################################################################
##
#F  HeckeCharTable( <W>, <rep>, <elts>) . .character table of a Hecke algebra
##  
HeckeCharTable:=function(W,c,irr)
  local C,G,pc;
  C:=rec(name:=ConcatenationString("H()"));
  C.identifier:=C.name;
  C.cartan:=W.cartan;
  C.parameter:=W.parameter;
  G:=Group(W.permgens,());
  C.order:=Size(G);
  C.size:=Size(G);
  pc:=List(c,i->PermWeylWord(W,i));
  C.classtext:=c;
  C.orders:=List(pc,i->OrderPerm(i));
  C.centralizers:=List(pc,i->Size(Centralizer(G,i)));
  C.classes:=List(C.centralizers,i->C.size/i);
  C.powermap:=[];
  C.irreducibles:=irr; 
  C.operations:=CharTableOps;
  return C;
end;