###########################################################################
# Contribution to the Chevie Package
#
#  This file contaoins some supplementary programs for working with
#  braids
#    (C) Jean MICHEL 1995-1996
###########################################################################
##
#F  orbitcircperm(w,F) computes the equivalence class of w in the braid monoid
##    for the equivalence relation where x * F(y) and y * x are equivalent
##    F can be omitted if trivial
##
##   (this is sometimes the whole F-conjugacy class of w)
##
orbitcircperm:=function(arg)local res,current,v,w,s,c,W,B,F,b,gens;
  b:=arg[1];if Length(arg)=2 then F:=arg[2];else F:=x->x;fi;
  W:=CoxeterGroup(b);B:=Braid(W); res:=[b]; current:=1;
  while current<=Length(res) do
    c:=res[current];
    if c.pw0>0 then gens:=List(W.generators,B);
    else gens:=List(LeftDescentSet(W,c.elm[1]),B);
    fi;
    for s in gens do
      v:=s^-1*c*F(s);
      if not v in res then Add(res,v);fi;
    od;
    current:=current+1;
  od;
  return res;
end;

#############################################################################
##
#F  leftconjugations(w) computes all elements w' of the braid monoid such
##    that there exists x reduced, x * w' = w * x 
##
##   (according to Deligne, the whole conjugacy class is reached by a
##    succession of such operations)
##
leftconjugations:=function(w)local v,x,B,W;
  W:=CoxeterGroup(w);B:=Braid(W);
  for x in Elements(W) do
    v:=B(x)^-1*w*B(x);
    if v.pw0>=0 then Print(w,"->(",B(x),")",v,"\n"); fi;
  od;
end;

#############################################################################
##
#F  decompositions( W, x, y)   x,y are elements of the CoxeterGroup group W
##                       given as permutations of the roots
##
##  returns as a list of pairs all decompositions w=xy where l(x)+l(y)=l(w)
##
decompositions:=function(W,w)local res,i,v,rest,x,new;
  rest:=[[(),w]];
  res:=[];
  while Length(rest)>0 do
    Append(res,rest);
    new:=[];
    for x in rest do
      Append(new,List(LeftDescentSet(W,x[2]),
		      i->[x[1]*W.generators[i],W.generators[i]*x[2]]));
    od;
    rest:=Set(new);
  od;
  return res;
end;

########################################################################
##
#F  alpha(b) returns the longest reduced braid dividing the element b of
##           the braid monoid
##
alpha:=function(b)local W;
  W:=CoxeterGroup(b);
  if b.pw0>0 then return Braid(W)(LongestCoxeterElement(W));
  elif Length(b.elm)>0 then return  Braid(W)(b.elm[1]);
  else return Braid(W)();
  fi;
end;

########################################################################
##
#F  gcd(a,b) returns the (left) gcd of elements a,b of the braid monoid
##
gcd:=function(a,b)local x,y,res,W;
 W:=CoxeterGroup(b);
 res:=Braid(W)();
 repeat
  x:=LeftDescentSet(W,PermBraid(alpha(a)));
  y:=LeftDescentSet(W,PermBraid(alpha(b)));
  x:=Intersection(x,y);
  y:=Braid(W)(LongestCoxeterElement(ReflectionSubgroup(W,x)));
  res:=res*y;a:=y^-1*a;b:=y^-1*b;
 until Length(x)=0;
 return res;
end;

########################################################################
##
#F  rev(a) returns the 'reversed' of element a of the braid monoid
##
rev:=a->Braid(CoxeterGroup(a))(Reversed(WordBraid(a)));

########################################################################
##
#F  lcm(a,b) returns the (left) lcm of elements a,b of the braid monoid
##
lcm:=function(a,b)local W,w0;
  W:=CoxeterGroup(a);
  w0:=Braid(W)(LongestCoxeterElement(W));
  return rev(w0*gcd(rev(a^-1*w0),rev(b^-1*w0))^-1);
end;

########################################################################
##
#F  allpaths(b,lim) returns all elements of the centralizer of b in
##                  the braid monoid obtained as paths of length <=lim
##                  in the graph of cyclic shift conjugacy for element b
##
allpaths:=function(b,lim)
  local i,p,newpaths,res,oldpaths,res,edges,current,v,w,s,c,W,B,F,b,gens;
  W:=CoxeterGroup(b);B:=Braid(W);
  res:=[b];edges:=[];
  current:=1;
  while current<=Length(res) do
    c:=res[current];
    if c.pw0>0 then gens:=List(W.generators,B);
    else gens:=List(LeftDescentSet(W,c.elm[1]),B);
    fi;
    for s in gens do
      v:=s^-1*c*s;
      Add(edges,[c,s,v]);
      if not v in res then Add(res,v);fi;
    od;
    current:=current+1;
  od;
  oldpaths:=[[b,B()]];res:=[];
  for i in [1..lim] do
    newpaths:=[];
    for p in oldpaths do
       Append(newpaths,List(Filtered(edges,x->x[1]=p[1]),e->[e[3],p[2]*e[2]]));
    od;
    newpaths:=Set(newpaths);
    for p in newpaths do
      if p[1]=b then 
	   Print(p[2],"\n");
	   if not p[2] in res then Add(res,p[2]);fi;
      fi;
    od;
    oldpaths:=newpaths;
  od;
  return res;
end;

########################################################################
##
#F  puiss(n,x,F) let x be an object (monoid element) on which F acts.
##               this function computes  (xF)^n 
##
puiss:=function(n,x,F)  
  if n=0 then return x^0;
  else return x* F(puiss(n-1,x,F));
  fi;
end;

########################################################################
##
#F  roots(wF,d) this function returns all d-th roots of pi in the
##              Coxeter coset WF using a very bestial algorithm
##
roots:=function(WF,d)local W,e,pi,F,p,phiorder;
  phiorder:=function(w,p)
    return First(DivisorsInt(OrderPerm(w*p)),i->(w*p)^i=p^i);
  end;
  if IsBound(WF.phi) then W:=CoxeterGroup(WF);p:=WF.phi;F:=Frobenius(WF);
		     else W:=WF;p:=();F:=x->x;
  fi;
  e:=CoxeterElementsLength(W,2*W.N/d);
  e:=List(Filtered(e,x->phiorder(x,p)=d),Braid(W));
  pi:=Braid(W)(LongestCoxeterElement(W))^2;
  e:=Filtered(e,x->puiss(d,x,F)=pi);
  return e;
end;

#test:=function(WF)local F,r,d,rr,n,pi,b,r,W,red,bad,rbad;
#  if IsBound(WF.phi) then W:=CoxeterGroup(WF);F:=Frobenius(WF);
#  else W:=WF;F:=x->x;fi;
#  r:=RegularNumbers(WF);
#  r:=r{[2..Length(r)]};
#  pi:=Braid(W)(LongestCoxeterElement(W))^2;
#  for d in r do
#    rr:=roots(WF,d);
#    Print("d=",d," reduced roots:",Length(rr));
#    rr:=orbitcircperm(rr[1],F);
#    Print(" orbitcirc:", Length(rr));
#    bad:=0;
#    rbad:=0;
#    red:=0;
#    for b in rr do
#      if Length(b.elm)<2 then red:=red+1;fi;
#      if puiss(d,b,F)<>pi then 
#          if Length(b.elm)<2 then rbad:=rbad+1;fi;
#	  bad:=bad+1;
#      fi;
#    od;
#    Print(" reduced:",red," bad:",bad," rbad:",rbad,"\n");
#  od;
#end;