#run this first to initialize data about root system
# w=weights, r=roots, dec=decompositions of s_\alpha, b=simple roots, 
# rh=rho, rs=name of root system
prepare:=proc(rs0) global w,r,dec,b,rh,rs;
rs:=rs0;w:=weights(rs);
rh:=rho(rs);
b:=base(rs);
root_dec();
print(`ready`);
end;

#outputs lambda-chain and the list of levels
lchain:=proc(lam) local i,k,v,kk,n,av,a,comp,ll,ord,aux,lc,lev;global w,r,rs;
comp:=proc(a,b) local i;
i:=1;
while a[i]=b[i] do i:=i+1 od;
if a[i]<b[i] then 0 else 1 fi;
end;
n:=nops(w);
lc:=[];v:=[];lev:=[];
for i from 1 to nops(r) do
  av:=2*r[i]/iprod(r[i],r[i]);k:=iprod(lam,av);
  if k>0 then
   lc:=[op(lc),seq(r[i],kk=0..k-1)];
   lev:=[op(lev),seq(kk,kk=0..k-1)];
   a:=[seq(iprod(w[kk],av)/k,kk=1..n)]; 
   v:=[op(v),seq([kk/k,op(a)],kk=0..k-1)];
  fi;
od;
ll:=nops(lc);ord:=false;
while (not ord) and (ll>1) do
 ord:=true;
 for i from 1 to ll-1 do
  if comp(v[i],v[i+1])=1 then
   ord:=false;
   aux:=v[i];v:=subsop(i=v[i+1],v);v:=subsop(i+1=aux,v);
   aux:=lc[i];lc:=subsop(i=lc[i+1],lc);lc:=subsop(i+1=aux,lc);
   aux:=lev[i];lev:=subsop(i=lev[i+1],lev);lev:=subsop(i+1=aux,lev);
  fi;
 od;
 ll:=ll-1;
od;
lc,lev;
end;

#computes decomposition of s_\alpha for every positive root; only needed in initialization
root_dec:=proc() local i,v,kk,n,av,a,comp,ll,ord,aux,sr,k,sr0;global r,dec,w,b,rh,rs;
comp:=proc(a,b) local i;
i:=1;
while a[i]=b[i] do i:=i+1 od;
if a[i]<b[i] then 0 else 1 fi;
end;
r:=pos_roots(rs);n:=nops(w);
v:=[];
for i from 1 to nops(r) do
   av:=2*r[i]/iprod(r[i],r[i]);k:=iprod(rh,av);
   a:=[seq(iprod(w[kk],av)/k,kk=1..n)]; 
   v:=[op(v),a];
od;
ll:=nops(r);ord:=false;
while (not ord) and (ll>1) do
 ord:=true;
 for i from 1 to ll-1 do
  if comp(v[i],v[i+1])=1 then
   ord:=false;
   aux:=v[i];v:=subsop(i=v[i+1],v);v:=subsop(i+1=aux,v);
   aux:=r[i];r:=subsop(i=r[i+1],r);r:=subsop(i+1=aux,r);
  fi;
 od;
 ll:=ll-1;
od;
sr0:=[r[1],seq(reflect(seq(r[i],i=1..kk-1),r[kk]),kk=2..nops(r))];sr:=[];
for i from 1 to nops(sr0) do
 k:=1;while b[k]<>sr0[i] do k:=k+1 od;sr:=[op(sr),k];
od;
dec:=[seq(reduce([seq(sr[i],i=1..kk-1),sr[kk],seq(sr[kk-i],i=1..kk-1)],rs),kk=1..nops(sr))];
end;

upbruhat1:=proc(r0,p) local i,d;global r,dec,rs;
 i:=1;while r[i]<>r0 do i:=i+1 od;
 d:=reduce([op(p),op(dec[i])],rs);
 if nops(d)=nops(p)+1 then true,d else false,[] fi;
end;

#outputs all foldings of a lambda-chain and the list of their keys
fold:=proc(lc) local i,j,u,r0,f,p,pos,fld,nn,sp,aux,key;
f:=[[]];nn:=nops(lc);key:=[[]];
p:=[];sp:=[p];
pos:=1;fld:=[];
while (pos<=nn) and (nops(f)<20000) do
 r0:=lc[pos];
 u:=upbruhat1(r0,p);
 if u[1] then 
  fld:=[op(fld),pos];
  p:=u[2];
  sp:=[op(sp),p];f:=[op(f),fld];key:=[op(key),p];
%print(f);
  if irem(nops(f),2000)=0 then print(nops(f)) fi;
 fi;
 pos:=pos+1;
 while (pos>nn) and (not (fld=[])) do
  pos:=fld[-1]+1;p:=sp[-2];
  fld:=subsop(-1=NULL,fld);sp:=subsop(-1=NULL,sp);
 od;
od;
f,key;
end;

#gives weight in coordinate form
weight:=proc(lc,lev,fld,lam) local i,nf,a,w,n;global rs;
nf:=nops(fld);w:=lam;
for i from nf to 1 by -1 do
 a:=lc[fld[i]];
 w:=reflect(a,w)+lev[fld[i]]*a;
od;
weight_coords(w,rs);
end;

#converts adm. subsequence to list of roots (folded chain)
admfold:=proc(lc,fld) local i,nn,a,flc,j,indf,ii;
nn:=nops(lc);flc:=[seq([lc[i],lc[i]],i=1..nn)];
#fold chain
for i from nops(fld) to 1 by -1 do
 ii:=fld[i];a:=lc[ii];
 flc[ii][2]:=-a;
#better algorithm by computing r_j_1...r_j_i; later do (r_j_1...r_j_s(\rho),\alpha_p^\vee)
 for j from ii+1 to nn do
  if flc[j][1]=flc[j][2] then indf:=false else indf:=true fi;
  flc[j][1]:=reflect(a,flc[j][1]);
  if indf then flc[j][2]:=-flc[j][1] else flc[j][2]:=flc[j][1] fi;
 od;
od;
flc;
end;

tstlrfold:=proc(nu,fld,key,lc) 
local n,i,j,flc,nn,indf,a1,a2,aa,mup,mp,lrch,tstneg,kr,found;global b,r;
tstneg:=proc(kr,a) local aa;global r;
aa:=reflect(op(kr),a);
if member(aa,r) then false else true fi
end;
n:=nops(b);nn:=nops(lc);
kr:=[seq(b[key[nops(key)+1-i]],i=1..nops(key))];
flc:=admfold(lc,fld);
i:=1;lrch:=true;
while lrch and (i<=n) do
 aa:=0;a1:=0;a2:=0;found:=false;
 for j from 1 to nn do
  if (flc[j][1]=b[i]) or (flc[j][1]=-b[i]) then 
   found:=true;
   if flc[j][2]=-b[i] then indf:=false else indf:=true fi;
   if flc[j][1]=b[i] then 
    a1:=aa;if flc[j][2]=b[i] then aa:=aa+1 fi;
   else a1:=aa-1;aa:=a1 
   fi;
   if a1>a2 then a2:=a1 fi;
  fi;
 od;
 if found then if indf then a1:=a1+1 else if tstneg(kr,b[i]) then a1:=a1-1 fi fi fi;
 mup:=a1;if a1>a2 then mp:=a1 else mp:=a2 fi;
 if iprod(nu,2*b[i]/iprod(b[i],b[i]))+mup<mp then lrch:=false else i:=i+1 fi;
od;
lrch;
end;

#inserts a number into a sequence preserving order.
insfold:=proc(l,a) local i,ll;
ll:=[];
i:=1;while (i<=nops(l)) and (l[i]<a) do ll:=[op(ll),l[i]];i:=i+1 od;
ll:=[op(ll),a];
while i<=nops(l) do ll:=[op(ll),l[i]];i:=i+1 od;
ll;
end;

fp:=proc(q,lc,fld) 
local i,j,flc,ii,nn,indf,a1,a2,aa,pos,posmax,inf,mp,fld1,found;
global b,r;
nn:=nops(lc);inf:=999999;fld1:=fld;
pos:=[];
flc:=admfold(lc,fld);
i:=q;
 aa:=0;a1:=0;a2:=0;found:=false;
 for j from 1 to nn do
  if (flc[j][1]=b[i]) or (flc[j][1]=-b[i]) then 
   found:=true;
   if flc[j][2]=-b[i] then indf:=false else indf:=true fi;
   if flc[j][1]=b[i] then 
    a1:=aa;if flc[j][2]=b[i] then aa:=aa+1 fi;
   else a1:=aa-1;aa:=a1 
   fi;
   if a1>a2 then a2:=a1;posmax:=nops(pos)+1 fi;
   pos:=[op(pos),j];
  fi;
 od;
 if found and indf then a1:=a1+1 fi;
 if a1>a2 then posmax:=inf;mp:=a1 else mp:=a2 fi;
if mp>0 then 
 if posmax<inf then
  i:=1;while fld1[i]<>pos[posmax] do i:=i+1 od; fld1:=subsop(i=NULL,fld1);
  insfold(fld1,pos[posmax-1]);
 else insfold(fld1,pos[-1])
 fi;
else 0
fi;
end;

#tstneg:=proc(kr,a) local aa;global r;
#aa:=reflect(op(kr),a);
#if member(aa,r) then false else true fi
#end;

#needs weight in vector form; returns both folding and weight
ep:=proc(q,lc,fld,we) 
local i,j,flc,ii,nn,indf,a1,a2,aa,pos,posmax,mp,fld1;
global b;
nn:=nops(lc);fld1:=fld;
pos:=[];
#kr:=[seq(b[key[nops(key)+1-i]],i=1..nops(key))];
flc:=admfold(lc,fld);
i:=q;
 aa:=0;a1:=0;a2:=0;
# found:=false;
 for j from 1 to nn do
  if (flc[j][1]=b[i]) or (flc[j][1]=-b[i]) then 
#   found:=true;
#   if flc[j][2]=-b[i] then indf:=false else indf:=true fi;
   if flc[j][1]=b[i] then 
    a1:=aa;if flc[j][2]=b[i] then aa:=aa+1 fi;
   else a1:=aa-1;aa:=a1 
   fi;
   if a1>=a2 then a2:=a1;posmax:=nops(pos)+1 fi;
   pos:=[op(pos),j];
  fi;
 od;
 a1:=iprod(we,2*b[q]/iprod(b[q],b[q]));
# if found then if indf then a1:=a1+1 else if tstneg(kr,b[i]) then a1:=a1-1 fi fi fi;
#print(pos,posmax,a1,a2);
if a1>=a2 then 0,0 else
 i:=1;while fld1[i]<>pos[posmax] do i:=i+1 od; fld1:=subsop(i=NULL,fld1);
 if posmax=nops(pos) then fld1,we+b[q] else insfold(fld1,pos[posmax+1]),we+b[q] fi;
fi;
end;

#needs weight of fld2 in vector form; returns both pair of foldings and weight of new fld1,fld2; if empty result, then return 0,0
fp2:=proc(q,lc1,fld1,lc2,fld2,w1,w2) 
local i,j,flc,ii,nn,a1,a2,aa,lrch,pos1,posmax1,inf,
mp1,mup1,mup2,pos2,posmax2,mp2,fld11,fld21;global b,r;
nn:=nops(lc1);inf:=999999;fld11:=fld1;
pos1:=[];
flc:=admfold(lc1,fld1);
i:=q;
 aa:=0;a1:=0;a2:=0;
#found:=false;
 for j from 1 to nn do
  if (flc[j][1]=b[i]) or (flc[j][1]=-b[i]) then 
#   found:=true;
#   if flc[j][2]=-b[i] then indf:=false else indf:=true fi;
   if flc[j][1]=b[i] then 
    a1:=aa;if flc[j][2]=b[i] then aa:=aa+1 fi;
   else a1:=aa-1;aa:=a1 
   fi;
   if a1>a2 then a2:=a1;posmax1:=nops(pos1)+1 fi;
   pos1:=[op(pos1),j];
  fi;
 od;
# if found and indf then a1:=a1+1 fi;
 mup1:=iprod(w1,2*b[q]/iprod(b[q],b[q]));

 if mup1>a2 then posmax1:=inf;mp1:=mup1 else mp1:=a2 fi;

nn:=nops(lc2);fld21:=fld2;
pos2:=[];
flc:=admfold(lc2,fld2);
i:=q;
 aa:=0;a1:=0;a2:=0;
#found:=false;
 for j from 1 to nn do
  if (flc[j][1]=b[i]) or (flc[j][1]=-b[i]) then 
#   found:=true;
#   if flc[j][2]=-b[i] then indf:=false else indf:=true fi;
   if flc[j][1]=b[i] then 
    a1:=aa;if flc[j][2]=b[i] then aa:=aa+1 fi;
   else a1:=aa-1;aa:=a1 
   fi;
   if a1>a2 then a2:=a1;posmax2:=nops(pos2)+1 fi;
   pos2:=[op(pos2),j];
  fi;
 od;
# if found then if indf then a1:=a1+1 else ...;
 mup2:=iprod(w2,2*b[q]/iprod(b[q],b[q]));

 if mup2>a2 then posmax2:=inf;mp2:=mup2 else mp2:=a2 fi;

if mp1+mup2-mp2>0 then

if mp1>0 then 
 if posmax1<inf then
  i:=1;while fld11[i]<>pos1[posmax1] do i:=i+1 od; fld11:=subsop(i=NULL,fld11);
  [insfold(fld11,pos1[posmax1-1]),fld21],w1-b[q],w2;
 else [insfold(fld11,pos1[-1]),fld21],w1-b[q],w2
 fi;
else 0,0,0
fi;

else 

if mp2>0 then 
 if posmax2<inf then
  i:=1;while fld21[i]<>pos2[posmax2] do i:=i+1 od; fld21:=subsop(i=NULL,fld21);
  [fld11,insfold(fld21,pos2[posmax2-1])],w1,w2-b[q];
 else [fld11,insfold(fld21,pos2[-1])],w1,w2-b[q]
 fi;
else 0,0,0
fi;

fi;
end;

repf:=proc(lc1,fld1,sf) local i,r;
r:=fld1;i:=1;
while (i<=nops(sf)) and (r<>0) do
 r:=fp(sf[i],lc1,r);i:=i+1;
od;
r;
end;


#needs weight of fld1,fld2 in vector form; returns foldings and weights of new fld1,fld2; if empty result, then return 0,0,0
ep2:=proc(q,lc1,fld1,lc2,fld2,w1,w2) 
local i,j,flc,ii,nn,a1,a2,aa,lrch,pos1,posmax1,inf,
mp1,mup1,mup2,pos2,posmax2,mp2,fld11,fld21;global b,r;
nn:=nops(lc1);inf:=999999;fld11:=fld1;
pos1:=[];
flc:=admfold(lc1,fld1);
i:=q;
 aa:=0;a1:=0;a2:=0;
#found:=false;
 for j from 1 to nn do
  if (flc[j][1]=b[i]) or (flc[j][1]=-b[i]) then 
#   found:=true;
#   if flc[j][2]=-b[i] then indf:=false else indf:=true fi;
   if flc[j][1]=b[i] then 
    a1:=aa;if flc[j][2]=b[i] then aa:=aa+1 fi;
   else a1:=aa-1;aa:=a1 
   fi;
   if a1>=a2 then a2:=a1;posmax1:=nops(pos1)+1 fi;
   pos1:=[op(pos1),j];
  fi;
 od;
# if found and indf then a1:=a1+1 fi;
 mup1:=iprod(w1,2*b[q]/iprod(b[q],b[q]));

 if mup1>a2 then posmax1:=inf;mp1:=mup1 else mp1:=a2 fi;

nn:=nops(lc2);fld21:=fld2;
pos2:=[];
flc:=admfold(lc2,fld2);
i:=q;
 aa:=0;a1:=0;a2:=0;
#found:=false;
 for j from 1 to nn do
  if (flc[j][1]=b[i]) or (flc[j][1]=-b[i]) then 
#   found:=true;
#   if flc[j][2]=-b[i] then indf:=false else indf:=true fi;
   if flc[j][1]=b[i] then 
    a1:=aa;if flc[j][2]=b[i] then aa:=aa+1 fi;
   else a1:=aa-1;aa:=a1 
   fi;
   if a1>=a2 then a2:=a1;posmax2:=nops(pos2)+1 fi;
   pos2:=[op(pos2),j];
  fi;
 od;
# if found then if indf then a1:=a1+1 else ...;
 mup2:=iprod(w2,2*b[q]/iprod(b[q],b[q]));

 if mup2>a2 then posmax2:=inf;mp2:=mup2 else mp2:=a2 fi;

if mp1+mup2-mp2>=0 then

if mup1=mp1 then 0,0,0 else

#print(fld11,pos1,posmax1);

 i:=1;while fld11[i]<>pos1[posmax1] do i:=i+1 od; fld11:=subsop(i=NULL,fld11);
 if posmax1=nops(pos1) then [fld11,fld21],w1+b[q],w2 
                       else [insfold(fld11,pos1[posmax1+1]),fld21],w1+b[q],w2 fi;
fi;

else

if mup2=mp2 then 0,0,0 else
 i:=1;while fld21[i]<>pos2[posmax2] do i:=i+1 od; fld21:=subsop(i=NULL,fld21);
 if posmax2=nops(pos2) then [fld11,fld21],w1,w2+b[q]
                       else [fld11,insfold(fld21,pos2[posmax2+1])],w1,w2+b[q] fi;
fi;

fi;

end;

repe:=proc(lc1,fld1,sf,we) local i,r,res;
r:=fld1;i:=1;w:=we;res:=r,w;
while (i<=nops(sf)) and (res[1]<>0) do
 res:=ep(sf[i],lc1,res[1],res[2]);if res[1]<>0 then r:=res[1] fi;i:=i+1;
od;
r;
end;


#fld={f1,f2}, apply root ops on pair until get to the bottom of crystal
gotobot:=proc(lc1,lc2,lev1,lev2,fld,lam1,lam2) local ww,mup2,mup1,i,q,n, ready,r,rr,s;global b,w;
s:=[];q:=1;ready:=false;n:=nops(b);
ww:=weight(lc1,lev1,fld[1],lam1);
mup1:=sum('ww[i]*w[i]','i'=1..n);
ww:=weight(lc2,lev2,fld[2],lam2);
mup2:=sum('ww[i]*w[i]','i'=1..n);
r:=fld,mup1,mup2;
while not ready do
 rr:=fp2(q,lc1,r[1][1],lc2,r[1][2],r[2],r[3]);
 if rr[1]<>0 then 
  s:=[op(s),q];r:=rr;q:=1;
 else q:=q+1; if q=n+1 then ready:=true fi;
 fi;
od;
s,r[1];
end;

gotobot1:=proc(lc,fld) local i,q,n, ready,r,rr,s;global b;
s:=[];q:=1;r:=fld;ready:=false;n:=nops(b);
while not ready do
 rr:=fp(q,lc,r);
 if rr<>0 then 
  s:=[op(s),q];r:=rr;q:=1;
 else q:=q+1; if q=n+1 then ready:=true fi;
 fi;
od;
s;
end;

#fld={f1,f2}, apply root ops on pair until get to the bottom of crystal
gototop:=proc(lc1,lc2,lev1,lev2,fld,lam1,lam2) local ww,mup2,mup1,i,q,n, ready,r,rr,s;global b,w;
s:=[];q:=1;ready:=false;n:=nops(b);
ww:=weight(lc1,lev1,fld[1],lam1);
mup1:=sum('ww[i]*w[i]','i'=1..n);
ww:=weight(lc2,lev2,fld[2],lam2);
mup2:=sum('ww[i]*w[i]','i'=1..n);
r:=fld,mup1,mup2;
while not ready do
 rr:=ep2(q,lc1,r[1][1],lc2,r[1][2],r[2],r[3]);
 if rr[1]<>0 then 
  s:=[op(s),q];r:=rr;q:=1;
 else q:=q+1; if q=n+1 then ready:=true fi;
 fi;
od;
s,r[1];
end;

gototop1:=proc(lc,fld,we) local i,q,n, ready,r,rr,s,w;global b;
s:=[];q:=1;r:=fld;ready:=false;n:=nops(b);w:=we;
while not ready do
 rr:=ep(q,lc,r,w);
 if rr[1]<>0 then 
  s:=[op(s),q];r:=rr[1];w:=rr[2];q:=1;
 else q:=q+1; if q=n+1 then ready:=true fi;
 fi;
od;
s;
end;

#computes key as a reduced word
compkey:=proc(lc,fld) local i,j,p;global rs,r,dec;
p:=[];
for i from 1 to nops(fld) do 
 j:=1;while r[j]<>lc[fld[i]] do j:=j+1 od;
 p:=reduce([op(p),op(dec[j])],rs);
od;
p;
end;

#involution on simple roots given by w0
w0s:=proc(k) local r,i,le;global b,rs;
le:=longest_elt(rs);
le:=[seq(b[le[i]],i=1..nops(le))];
r:=-reflect(op(le),b[k]);
i:=1;while b[i]<>r do i:=i+1 od;
i;
end;


#action of a simple reflection on a folding
sp:=proc(q,lc,fld,we) local i,f1,w,we1;global b;
w:=iprod(we,2*b[q]/iprod(b[q],b[q]));f1:=fld;we1:=we;
if w<0 then for i from 1 to -w do f1:=ep(q,lc,f1,we1);we1:=f1[2];f1:=f1[1] od
else if w>0 then for i from 1 to w do f1:=fp(q,lc,f1) od fi fi;
f1;
end;

#action of a permutation on a folding
actw:=proc(r,lc,fld,we) local i,f1,we1;global b;
f1:=fld;we1:=we;
for i from nops(r) to 1 by -1 do f1:=sp(r[i],lc,f1,we1);we1:=reflect(b[r[i]],we1) od;
f1;
end;

#action of a simple reflection on a pair of foldings
sp2:=proc(q,lc1,fld1,lc2,fld2,w1,w2) local i,f1,f2,res,we2,we1,w;global b;
w:=iprod(w1+w2,2*b[q]/iprod(b[q],b[q]));
f1:=fld1;f2:=fld2;we1:=w1;we2:=w2;
if w<0 then for i from 1 to -w do 
     res:=ep2(q,lc1,f1,lc2,f2,we1,we2);
     we1:=res[2];we2:=res[3];f1:=res[1][1];f2:=res[1][2] od
else if w>0 then for i from 1 to w do res:=fp2(q,lc1,f1,lc2,f2,we1,we2);we1:=res[2];we2:=res[3];
     f1:=res[1][1];f2:=res[1][2] od fi fi;
f1,f2,we1,we2;
end;

#for the moment returns just weights
#action of a permutation on a pair of foldings; 
actw2:=proc(r,lc1,fld1,lc2,fld2,w1,w2) local i,f1,f2,we1,we2,res,tr;global b;
f1:=fld1;f2:=fld2;we1:=w1;we2:=w2;
for i from nops(r) to 1 by -1 do 
 res:=sp2(r[i],lc1,f1,lc2,f2,we1,we2);
 we1:=res[3];we2:=res[4];
 f1:=res[1];f2:=res[2];
od;
[f1,f2];
end;


#action of involution S on a folding
invol:=proc(lc,fld) local lr,i;
lr:=gotobot1(lc,fld);
lr:=[seq(w0s(lr[nops(lr)+1-i]),i=1..nops(lr))];
repf(lc,[],lr)
end;

commutator:=proc(lc1,lc2,lev2,fld2,p2,lam2) local i,lr1,sfld1,fld1,sf2,fld11,sf1,fld21,p1,a,t,lr2;global rs;
sfld1:=gotobot(lc1,lc2,lev2,[[],fld2],lam2)[2][1];
lr1:=gotobot1(lc1,sfld1);
lr1:=[seq(w0s(lr1[nops(lr1)+1-i]),i=1..nops(lr1))];
#test if fld1, fld2 obtained by inverse sequences of root operators
#lr2:=[seq(lr1[nops(lr1)+1-i],i=1..nops(lr1))];
#fld21:=repf(lc2,[],lr2);
#if fld21=fld2 then print(`**** Checked ****`,fld2,lr1) else ERROR(`NOT INVERSE`,fld1,fld2) fi;
fld1:=repf(lc1,[],lr1);
[[],fld1];
end;