##########################################################################
#    Contribution to the Chevie package
#    Meinolf Geck -Jean Michel   May 1995
#    
#    This file contains functions implementing the algorithm described
#    in our paper ``On ``good'' elements in the conjugacy classes of
#     finite Coxeter groups and their eigenvalues on the irreducible
#     representations of Iwahori-Hecke algebras.''
#       Proc.\ London Math.\ Soc. (to appear).
#
#    At the end of the file, an example shows the instructions used
#    to get the character table for type E8
##########################################################################
# The following function transforms a polynomial to double list
# [degrees of non 0 coeffs, coeffs]
poltocouple:=function(p)local c,ind;
  c:= [p.valuation..p.valuation+Length(p.coefficients)-1];
  ind:=Filtered([1..Length(p.coefficients)],i->p.coefficients[i]<>0);
  return [c{ind},p.coefficients{ind}];
  end;

# The following function keeps only those eigenvalues which can contribute
# to $\chi(T_w)$, i.e. those $v^a$ such that $a/2$ is integral if $\chi$
# is rational or such that $a$ is integral otherwise. The second argument
# is a list of pairs |[i,j]| such that $\chi_i$ is Galois-conjugate to
# $\chi_j$. |eig| contains for each $\chi$ just
# the list $n_{\chi,1},\ldots,n_{\chi,r}$ such that the eigenvalues are
# $v^{2 n_{\chi,1}},\ldots,v^{2n_{\chi,r}}$

reduce:=function(eig,gal)
  return List([1..Length(eig)],i->Filtered(eig[i],
	      function(x) if i in Concatenation(gal)
	                  then return IsInt(2*x);
			  else return IsInt(x);fi;end));
  end;

# Now is a series of routines giving relations 1--7. The result of all these
# routines is always a list of vectors $[v_1,\ldots,v_r]$ such that each $v_i$
# is of same length $1+A$ where |A=Sum(eig,Length)|
# (the number of variables $a_{\chi,i}$) and represents a relation
#  $v_{A+1} + \sum_{j=1}^{j=A} a_j v_i=0$ (where $a_j$ is the $j$th variable
#  $a_{\chi,i}$).

# The next function returns the relations given by the value
# $\chi_(T_w)_{v\mapsto 1}$. The first argument is the index of the
# class of $w$ in the character table of $W$.
#
relsvalue:=function(index,W,eig)local i,v,res,model,ti;
  InfoChevie("using relations coming from character values\n");
  ti:=CharTable(W);
  res:=[];model:=List(eig,x->List(x,y->0));
  for i in [1..Length(eig)] do
    v:=Copy(model);v[i]:=v[i]+1;
    if Length(eig[i])>0 then
       v:=Concatenation(v);
       Add(v,-ti.irreducibles[i][index]);
       Add(res,v);
    fi;
  od;
  return res;
end;

# The next function returns the relations coming from Curtis-Alvis duality.
# The first argument is the length of $w$.

relsdual:=function(len,W,eig)local dual,model,res,i,c,v,ti,j;
  InfoChevie("using relations coming from Curtis-Alvis duality\n");
  ti:=CharTable(W);
  res:=[]; model:=List(eig,x->List(x,y->0));
  dual:=List(ti.irreducibles,x->Position(ti.irreducibles,
    List([1..Length(x)],i->x[i]*(-1)^Length(ti.classtext[i]))));
  InfoChevie("dual",dual,"\n");
  for i in [1..Length(ti.classtext)] do
   for c in [1..Length(eig[i])] do
     v:=Copy(model);
     v[i][c]:=-1;
     j:=Position(eig[dual[i]],len-eig[i][c]);
     if j<>false then v[dual[i]][j]:=v[dual[i]][j]+(-1)^len;fi;
# the addition is to handle correctly a self-dual coefficient
# the 'if' is because integrality condition may have been stronger for dual
     v:=Concatenation(v);Add(v,0);Add(res,v);
   od;
  od;
  return res;
end;

# The next function takes an argument |gal| as |reduce| does and returns
# relations expressing that $\chi_i$ and $\chi_j$ are Galois-conjugate

relsgal:=function(eig,gal)local model,c,v,res,x;
   InfoChevie("using relations coming from Galois-conjugation\n");
   res:=[]; model:=List(eig,x->List(x,y->0));
   for x in gal do
   for c in [1..Length(eig[x[1]])] do
     v:=Copy(model);
     v[x[1]][c]:=-1; v[x[2]][c]:=(-1)^(2*eig[x[1]][c]);
     v:=Concatenation(v);
     Add(v,0);
     Add(res,v);
   od;
  od;
  return res;
end;

# The next function (used in the one just after) returns the values on $T_w$
# of the exterior powers of the reflection character of the Hecke
# algebra of $W$ as polynomials in $v$.
# The first argument is the |CoxeterWord| representing $w$, the second
# argument is the Hecke algebra of $W$.
#
refchars:=function(w,H)local l,p,T,q,r;
 q:=H.parameter[1];T:=X(q.domain);r:=CoxeterGroup(H).semisimpleRank;
 if Length(w)=0 then l:=IdentityMat(r);
 else l:=Product(HeckeReflectionRepresentation(H){w});fi;
 p:=DeterminantMat(IdentityMat(r)*T-l);
 return List([1..r+1],i->p.coefficients[i]*q^(Length(w)*(i-r))*(-1)^(r+1-i));
end;

# The next function returns the indices of the exterior powers of the reflexion
# character in the list of characters of W
refinds:=function(W)local ti,refls;
  ti:=CharTable(W);
  refls:=Position(ti.irreducibles,List(CoxeterConjugacyClasses(W),
		  w->ReflectionCharValue(W,PermCoxeterWord (W,w))));
  refls:=List([W.semisimpleRank,W.semisimpleRank-1..1],
              i->AntiSymmetricParts(ti,ti.irreducibles{[refls]},i)[1]);
  refls:=List(refls,x->Position(ti.irreducibles,x));
  Add(refls,PositionProperty(ti.irreducibles,x->ForAll(x,y->y=1)));
  return refls;
end;

# The next function returns the relations expressing the value on $T_w$ of
# the exterior powers of the reflection character.
# The first argument is the |CoxeterWord| of $w$, and the second is the
# Hecke algebra.
#
relsrefl:=function(w,H,eig)
  local i,j,m,c,v,res,refs,model,pos,refind,ti;
  InfoChevie("using the values of the exterior power of reflection character\n");
  ti:=CharTable(H);
  refs:=refchars(w,H);
  refind:=refinds(CoxeterGroup(H));
  res:=[]; model:=List(eig,x->List(x,y->0));
  for i in [1..Length(refind)] do
    m:=poltocouple(refs[i]);
    if IsBound(H.sqrtParameter[1]) then m[1]:=m[1]/2;fi;
    v:=Copy(model);
    c:=IdentityMat(Length(eig[refind[i]]));
    for j in [1..Length(eig[refind[i]])] do
     v[refind[i]]:=c[j];
     pos:=Position(m[1],eig[refind[i]][j]);
     if pos=false then
       Add(res,Concatenation(Concatenation(v),[0]));
     else
       Add(res,Concatenation(Concatenation(v),[-m[2][pos]]));
       m[2][pos]:=0;
     fi;
    od;
    if Number(m[2],x->x<>0)>0 then 
     InfoChevie(" refs[",i,"]:",refs[i]," eig[",refind[i],"]",eig[refind[i]],"\n");
     Error("value of refl. chars");fi;
  od;
  return res;
end;
    
# The next function returns relations expressing
# $\varphi(T_w)_{v\mapsto \exp(\pi/d)}=0$ for virtual characters orthogonal
# to $\Phi_d$-projectives.
#  The third argument is a list of vectors $v_i$ such that $v_i$ is the
# list of coefficients of some $\varphi$ on the irreducibles.

relsd:=function(d,eig,null)local v,i,res,n,model;
  InfoChevie("using relations from the Phi_",d,"-decomposition matrix\n");
  res:=[];
  model:=List(eig,k->List(2*k,x->E(2*d)^x));
  for i in null do
     v:=Concatenation(List([1..Length(model)],k->model[k]*i[k]));
     Add(v,0);
     n:=NofCyc(v);
     v:=TransposedMat(List(v,x->CoeffsCyc(x,n)));
     v:=Filtered(v,x->Number(x,y->y<>0)>0);
     Append(res,v);
  od;
  return res;
end;

# The next function returns the relations coming from 
# $\sum_\chi \chi(T_w)d_\chi=0$. The second argument is the Hecke algebra.
#
relsdegg:=function(eig,H)local offs,v,gendeg,res,i,j,k,degs,factor;
  InfoChevie("using relations coming from generic degrees\n");
  offs:=List(eig,Length);
  offs:=Concatenation([0],List([1..Length(offs)],i->Sum(offs{[1..i]})));
  v:=SchurElements(H);
  gendeg:=List(v,x->v[1]/x);
  res:=List([1..1+4*CoxeterGroup(H).N],x->[1..offs[Length(offs)]]*0);
  for i in [1..Length(eig)] do
   degs:=poltocouple(gendeg[i]);
   if IsBound(H.sqrtParameter[1]) then degs[1]:=degs[1]/2;fi;
   for j in [1..Length(eig[i])] do
     res{2*(eig[i][j]+degs[1])}[j+offs[i]]:=degs[2];
   od;
  od;
  res:=Filtered(res,x->Number(x,y->y<>0)>0);
  for v in res do Add(v,0);od;
  return res;
end;

# To save space we apply each new set of relations as soon as possible.
# We represent our current knowledge as a record |M| with 3 fields:
# \begin{itemize}
#  \item |M.known| is the indices of the variables $a_j$ that we know already.
#  \item |M.values| is the values of the variables $a_j$ that we know already.
#  \item |M.relations| is a set of vectors representing the relations on
#     the remaining variables.
# \end{itemize}
# The next function takes a bunch of new relations and modifies |M| accordingly.
# It uses the routine  |TriangulizeMat| to find all completely known basis
# vectors resulting from |M.relations|, and then supresses the corresponding
# columns from |M.relations|, adding entries instead to |M.known| and
# |M.values|.
# It returns |true| iff at the end of the process all values are known.

apply:=function(newrels,M)local compl,i,newval,newind;
  if Length(newrels)=0 then return false;fi;
  compl:=Filtered([1..Length(newrels[1])-1],i->not i in M.known);
  if Length(M.known)>0 then
    newrels:=List(newrels,
                x->Concatenation(x{compl},[x[Length(x)]+x{M.known}*M.values]));
  fi;
  Append(M.relations,newrels);
  InfoChevie("starting triangulization\n");
  TriangulizeMat(M.relations);
  M.relations:=M.relations{[1..RankMat(M.relations)]};
  if Length(M.relations)>0 and 
     PositionProperty(M.relations[Length(M.relations)],y->y<>0)=
              Length(M.relations[1])
  then Error("contradictory relations");fi;
  newval:=Filtered(M.relations,x->Number(x{[1..Length(x)-1]},y->y<>0)=1);
  M.relations:=Filtered(M.relations,x->Number(x{[1..Length(x)-1]},y->y<>0)<>1);
  newind:=List(newval,x->PositionProperty(x,y->y<>0));
  M.relations:=List(M.relations,x->Concatenation(
	x{Filtered([1..Length(x)-1],y->not y in newind)},[x[Length(x)]]));
  Append(M.known,compl{newind});
  Append(M.values,List(newval,x->-x[Length(x)]));
  SortParallel(M.known,M.values);
  InfoChevie("known:",Length(M.known)," unknown:");
  if Length(M.relations)>0 then
     InfoChevie(Length(M.relations[1])-1," rank:",RankMat(M.relations),"\n");
     return false;
  else InfoChevie("0\n");return true;fi;
end;

# Now we put everything together to get a routine which returns all
# character values on $T_w$. The arguments are:
# \begin{itemize}
#  \item The index of the class of $w$ in the character table of $W$.
#  \item The Hecke algebra of $W$.
#  \item A list |null| which should be such that |null[e]| is bound
#        precisely when $\Phi_e$ divides the Poincar\'e polynomial, and
#        then |null[e]| contains the kernel of the $\Phi_e$-decomposition
#        matrix s explained in |relsd|.
#  \item A list |gal| as explained in |reduce|.
# \end{itemize}
#
# The result is the list of character values on $T_w$ as polynomials in $v$.
#
getchar:=function(ind,H,null,gal)local eig,W,M,values,e,ti,w,d;
  InfoChevie("trying element ",ind," ...\n");

# The next function returns the polynomial described by |eig| and a record
# |M| such that |M.known| is everybody and |M.relations| is empty.
# Actually we accept that |M.know| is not everybody and then give the
# arbitrary value 999 to unknown entries, which is quite distinctive
# and allows us at a glance to see what is known or not in
# a partially known character value.
#
  values:=function(M,eig)local vals,offs;
    vals:=999+[1..Sum(eig,Length)]*0;
    vals{M.known}:=M.values;
    offs:=List(eig,Length);
    offs:=Concatenation([0],List([1..Length(offs)],i->Sum(offs{[1..i]})));
    vals:=List([1..Length(eig)],i->vals{[offs[i]+1..offs[i+1]]});
    return 
      List([1..Length(vals)],function(i)
	if Length(vals[i])=0 then return 0;
	else return vals[i]*List(eig[i],function(n)
              if IsInt(n) then return H.parameter[1]^n;
                          else return H.sqrtParameter[1]^(2*n);
              fi;end);
	fi;end);
    end;

  W:=CoxeterGroup(H);
  w:=CoxeterConjugacyClasses(W)[ind];
  d:=OrderPerm(PermCoxeterWord (W,w));
  eig:=List(HeckeCharValuesGood(H,w),x->poltocouple(x)[1]/(2*d));
  eig:=reduce(eig,gal);
  M:=rec(known:=[],values:=[],relations:=[]);
  if apply(relsvalue(ind,W,eig),M) then return values(M,eig);fi;
  if apply(relsgal(eig,gal),M) then return values(M,eig);fi;
  if apply(relsrefl(w,H,eig),M) then return values(M,eig);fi;
  for e in Filtered([1..Length(null)],j->IsBound(null[j])) do
    if apply(relsd(e,eig,null[e]),M) then return values(M,eig);fi;
  od;
  if apply(relsdegg(eig,H),M) then return values(M,eig);fi;
  if apply(relsdual(Length(w),W,eig),M) then return values(M,eig);fi;
  InfoChevie("not enough relations\n");return values(M,eig);
  end;

##########################################################################
##
#F RelationsDecMats( <arg> ) . . . . . . .  compute relations derived from
## . . . . . . . . . . . . . . . . . . Phi_e-modular decomposition numbers
##
## 'RelationsDecMats' returns a list (with holes) which contains at the
## e-th position a set of vectors which yield a linear  relation among
## the Phi_e-modular reductions of the irreducible characters.
## <arg> is either just a Cartan matrix (in which case e runs over all
## possible values for which Phi_e has multiplicity at least two in the
## Poincare polynomials) or a Cartan matrix and a list of specified 
## values for e.
##
## Examples: gap> RelationsDecMats(CartanMat("E",6));   # or
##           gap> RelationsDecMats(CartanMat("E",6), [2,3]); 
##
## The relations are obtained by induction of all possible linear
## independent relations from the maximal (proper) parabolic subalgebras.
##
RelationsDecMats:=function(arg)
  local NullSpaceMat,Koerper,HV,V,W,HW,HW1,ls,ind,e,cartan,
                                           null,rels,v,b,l,base,f,i,j;
  NullSpaceMat:=function(mat)
    local a,i,j,k,piv,null,opiv,n,m,bas,x;
    null:=0*mat[1][1];
    a:=TransposedMat(mat);
    TriangulizeMat(a);
    m:=Length(Filtered(a,i->i<>0*i));
    n:=Length(a[1]);
    piv:=[ ];
    for i in [1..m] do
      j:=i;
      while a[i][j]=null do
        j:=j+1;
      od;
      piv[i]:=j;
    od;
    opiv:=[ ];
    for i in [1..n] do
      if not i in piv then
        Add(opiv,i);
      fi;
    od;
    bas:=[ ];
    for i in [1..n-m] do
      x:=List([1..n],j->null);
      x[opiv[i]]:=null^0;
      for k in [1..m] do
        x[piv[k]]:=-a[k][opiv[i]];
      od;
      IsVector(x);
      Add(bas,x);
    od;
    return bas;
  end;
  
  Koerper:=function(e,min)
    local i,q;
    i:=1;
    while true do
      q:=1+2*e*i;
      if q>min and IsPrime(q) then
        return [i,q];
      fi;
      i:=i+1;
    od;
  end; 

  if Length(arg)=1 then
    cartan:=arg[1];
    W:=CoxeterGroup(cartan);
    ls:=[];
    for e in ReflectionDegrees(W) do
      Append(ls,Filtered([2..e],i->IsInt(e/i)));
    od;
    ls:=Set(ls);
    #ls:=Filtered(Set(ls),i->Number(ls,j->j=i)>1);
  else
    cartan:=arg[1];
    W:=CoxeterGroup(cartan);
    ls:=arg[2];
  fi;
  InfoChevie("\nModuli: ",ls,"\nComputing induce/restrict matrices ");
  ind:=[];
  for i in [1..W.semisimpleRank] do
    InfoChevie("+\c");
    V:=ReflectionSubgroup(W,Difference([1..W.semisimpleRank],[i]));
    ind[i]:=TransposedMat(InductionTable(V,W).scalar);
  od;
  null:=[];
  for e in ls do
    v:=Koerper(e,2*Maximum(Factors(Size(W))));
    f:=[[],[]];
    for j in [1..v[2]] do
      f[1][j]:=j-(v[2]+1)/2;
      f[2][j]:=f[1][j]*Z(v[2])^0;
    od;
    HW:=Hecke(CoxeterGroup(cartan),Z(v[2])^(2*v[1]),Z(v[2])^v[1]);
    HW1:=Hecke(CoxeterGroup(cartan),E(2*e)^2,E(2*e));
    InfoChevie("\ne = ",e,"; \c");
    InfoChevie(" using finite field of order ",v[2],"\n");
    rels:=[];
    for i in [1..W.semisimpleRank] do
      HV:=HeckeSubAlgebra(HW,Difference([1..W.semisimpleRank],[i]));
      base:=NullSpaceMat(CharTable(HV).irreducibles*Z(v[2])^0);
      HV:=HeckeSubAlgebra(HW1,Difference([1..W.semisimpleRank],[i]));
      for b in base do 
        for j in [1..Length(b)] do
          b[j]:=f[1][Position(f[2],b[j])];
        od;
        if b*CharTable(HV).irreducibles=0*b then
          Add(rels,Sum(List([1..Length(b)],j->b[j]*ind[i][j])));
        else
          InfoChevie("wrong relation \n");
        fi;
      od;
      TriangulizeMat(rels);
      rels:=Filtered(rels,j->j<>0*j);
      InfoChevie("No. of relations after using subalgebra of type ",
                                 CartanName(CoxeterGroup(HV))," : ",Length(rels),"\n");
    od;
    if rels<>[] then
      null[e]:=Copy(rels);
    fi;
  od;
  return null;
end;

# return indices of the cuspidal classes in the list of classes of W 
CuspidalClasses:=function(W) local cl;
  cl:=CoxeterConjugacyClasses(W);
  return Filtered([1..Length(cl)],i->Set(cl[i])=[1..W.semisimpleRank]);
end;

# The following lines use the above routines to write all the character
# values on cuspidal classes of the Hecke algebra of $E_8$.

# define $v$
v:=X(Rationals); v.name:="v";

# define the hecke algebra.
W:=CoxeterGroup("E",7);
H:=Hecke(W,v^2,v);

# Define the galois-relations
#gal:=[[62,105],[63,106]];
gal:=[];
gal:=[[59,60]];

# And go!
r:=RelationsDecMats(W.cartan);
l:=List(CuspidalClasses(W),i->getchar(i,H,r,gal));