#
# bruhat_order - Bruhat ordering of a finite Coxeter group
#
# Calling sequence:
#  bruhat_order(R,<options>);
#
# Parameters:
#  R = a root system data structure
# The optional arguments are any of the following, in any order:
#  m = a nonnegative integer
#  w = a list of integers representing an element of W(R)
#  Q = an unassigned name
#
# IMPORTANT: This requires the coxeter package (version 2.0 or later)
#
# Let R be a finite root system and W(R) the associated reflection group
# (a Coxeter group). The Bruhat ordering of W(R) is the transitive closure
# of the relations w < t*w for all w in W(R) and all reflections t in W(R)
# such that l(w) < l(t*w), where l() denotes the length function.
#
# If R is the name of a finite root system, bruhat_order(R) returns the
# covering relation of a poset isomorphic to the Bruhat ordering of W(R).
#
# If an optional group element w is specified, the subinterval of the Bruhat
# order from w to w0 (the longest element) is returned.
#
# If an optional integer m is specified, the poset is truncated at rank m.
#
# If an unassigned name is specified, then this name will be assigned the
# set of covering edges of the poset that belong to the (left) weak order;
# i.e., the relations between pairs that differ by a simple reflection.
#
# Examples:
#  with(coxeter): with(posets):
#
#  P:=bruhat_order(A3,'Q');
#  wk:=table([green=Q]):
#  plot_poset(P,ecolor=wk);
#  eulerian(P),lattice(P);    # the Bruhat order is Eulerian/not a lattice
#
#  P:=bruhat_order(H3,5,'Q');
#  wk:=table([green=Q]):
#  plot_poset(P,ecolor=wk,stretch=2.0);
#
# Use the antiautomorphism w -> w_0*w to extract the subinterval of the
# Bruhat order of E6 from [] to [1,2,3,4,5,6]:
#  w:=[op(longest_elt(E6)),1,2,3,4,5,6];
#  P:=bruhat_order(E6,w);
#  plot_poset(P,proportional);
#  isom(P,J({},6));
#
# Determine all isomorphism classes of subintervals of length 3 in the
# Bruhat ordering of B3:
#  P:=bruhat_order(B3): F:=filter(P);
#  L:=[seq(seq(seq(subinterval(P,[i,j]),i=F[r-3]),j=F[r]),r=4..nops(F))]:
#  rm_isom(L);
#
bruhat_order:=proc(R) local S,coS,m,n,v0,P,Q,new,old,p,q,i,j,c,v,EPS;
  S:=coxeter['base'](R);
  coS:=[seq(2*v/coxeter['iprod'](v,v),v=S)];
  if type(S,'list'('polynom'('rational'))) then
    EPS:=0 else EPS:=`coxeter/default`['epsilon'] fi;
  v0:=coxeter['interior_pt'](S);
  v:=-v0; m:=coxeter['num_refl'](R);
  for i from 2 to nargs do
    if type(args[i],'list') then
      v:=coxeter['reflect'](seq(S[j],j=args[i]),-v0)
      elif type(args[i],'integer') then m:=args[i]
      elif type(args[i],'name') then Q:=NULL; n:=i
    fi
  od;
  P:=NULL; new:=[v]; q:=0;
  while m>0 do
    p:=q; q:=q+nops(new);
    old:=new; new:=[];
    for i to nops(old) do
      c:=`bruhat_order/covers`(old[i],v0,S,coS);
      for v in c do
        if `bruhat_order/new`(v,new,EPS,'j') then
          new:=[op(new),v]; j:=nops(new) fi;
        P:=P,[p+i,q+j]
      od;
      if type(n,'integer') then for c to nops(S) do
        v:=coxeter['iprod'](coS[c],old[i]);
        if v<0 then v:=old[i]-v*S[c] else next fi;
        `bruhat_order/new`(v,new,EPS,'j'); Q:=Q,[p+i,q+j]
      od fi
    od;
    m:=m-1; if new=[] then break fi;
  od;
  if type(n,'integer') then assign(args[n],{Q}) fi; {P}
end;

`bruhat_order/covers`:=proc(u,v0,S,coS) local w,alive,v,i;
  coxeter['vec2fc'](u,S,'w'); v:=v0;
  i:=nops(w); alive:=[];
  while i>0 do
    alive:=map((x,y)->[x,coxeter['iprod'](y,x)],alive,coS[w[i]]);
    alive:=[v,op(map(proc(x,y) if x[2]>0 then x[1]-x[2]*y fi end,
      alive,S[w[i]]))];
    v:=coxeter['reflect'](S[w[i]],v); i:=i-1
  od; alive
end;

`bruhat_order/new`:=proc(v,L,EPS) local i;
  if EPS=0 then RETURN(not member(v,L,args[4])) fi;
  for i to nops(L) do if
    convert(map(abs,[coeffs(L[i]-v)]),`+`)<EPS then
    assign(args[4],i); RETURN(false)
  fi od; true
end;