# compute HL polynomial in n variables: gensym(n);comphl1([partition]);

#read "N:/My Documents/soft/alcoves/hall-litt-a.txt":
with(combinat);



Perm2Length:=proc(p) local i,j,l;
l:=0;
for i from 1 to nops(p)-1 do
 for j from i+1 to nops(p) do
  if p[j]<p[i] then l:=l+1 fi
 od
od;
l
end;


gensym:=proc(n) local i,b0;global bc,sbc;
bc:=permute(n);
sbc:=[seq(Perm2Length(bc[i]),i=1..nops(bc))];
print(`number of elements`,nops(bc));
end;

lchf:=proc(n,nn) local l,i,j,m1;
l:=[];
for i from n to 1 by -1 do
 for j from n+1 to nn do l:=[[i,j],op(l)] od od;
m1:=nops(l);
l,m1;
end;

#lt listed as [3,2,1]
lch:=proc(lt,n) local i,lc,ll,r;
lc:=[];ll:=[[0,0]];
for i from nops(lt) to 1 by -1 do
 r:=lchf(lt[i],n);
 lc:=[op(lc),op(r[1])];
 ll:=[op(ll),[ll[-1][1]+r[2],lt[i]]];
od;
lc,ll;
end;

decl:=proc(w,p) if w[p[1]]>w[p[2]] then true else false fi;
end;


rev:=proc(s) local i; [seq(s[nops(s)+1-i],i=1..nops(s))] end;

tr:=proc(p,t) local aux,p1;
p1:=p;aux:=p1[t[1]];p1[t[1]]:=p1[t[2]];p1[t[2]]:=aux;p1
end;


#ll[1] limit in chain of perms;ll[2] height of column; typical [[0,0],[4,2],[7,1]] for [[1,4],[1,3][2,4],[2,3],[1,4],[1,3],[1,2]]
#returns list of tableaux for each folding plus related parameters
fold:=proc(lc,w,ll) local i,j,ww,lfld,lw,llw,fld,ready,lt,t,w1,k,kk;
#print(lc,w,ll);
lfld:=[[]];ww:=w;llw:=[[w]];
i:=1;lw:=[w];fld:=[];ready:=false;
while not ready do
 if i<=nops(lc) then
  if decl(ww,lc[i]) then
   ww:=tr(ww,lc[i]);
   lw:=[op(lw),ww];
   fld:=[op(fld),i];lfld:=[op(lfld),fld];llw:=[op(llw),lw];
  fi;
  i:=i+1;
 fi;
 if i>nops(lc) then
  if fld=[] then ready:=true
  else
   i:=fld[-1]+1;fld:=subsop(-1=NULL,fld);lw:=subsop(-1=NULL,lw);
   ww:=lw[-1];
  fi
 fi;
od;
#print(lfld);
lt:=[];
for i from 1 to nops(lfld) do
 t:=[];fld:=lfld[i];kk:=0;w1:=w;lw:=llw[i];
 for j from 1 to nops(ll)-1 do
  t:=[op(t),[seq(w1[k],k=1..ll[j+1][2])]];
  while (kk+1<=nops(fld)) and (fld[kk+1]<=ll[j+1][1]) do kk:=kk+1 od;
  if kk<=nops(fld) then w1:=lw[kk+1] fi;
 od;
 lt:=[op(lt),[t,Perm2Length(lw[-1]),nops(fld)]];
od;
#print(`----`,lt);
lt;
end;

#USE tab();
#lam as [2,1,1,0]; 
comphl1:=proc(lam) local lam1,lc,ll,n,i,hl,j,lt,k,nhl,r,lpt;global bc,sbc;
n:=nops(lam);lam1:=cnjpart(lam);
r:=lch(lam1,n);lc:=r[1];ll:=r[2];
#print(lc,ll);
hl:=array(1..80000);nhl:=0;lpt:=table();
for i from 1 to nops(bc) do
 j:=1;while (j<=n-1) and ((lam[j]<>lam[j+1]) or (bc[i][j]<bc[i][j+1])) do j:=j+1 od;
 if j=n then
   lt:=fold(lc,bc[i],ll);
   for j from 1 to nops(lt) do
    k:=lpt[op(lt[j][1])];
#print(w1,pt);
    if op(0,k)<>`Integer` then 
     nhl:=nhl+1;
     if nhl mod 200=0 then print(nhl,i) fi;
     lpt[op(lt[j][1])]:=nhl;
     hl[nhl]:=[t^((sbc[i]+lt[j][2]-lt[j][3])/2)*(1-t)^lt[j][3],lt[j][1]]
    else hl[k][1]:=expand(hl[k][1]+t^((sbc[i]+lt[j][2]-lt[j][3])/2)*(1-t)^lt[j][3])
    fi;
   od;
  fi;
od;
for i from 1 to nhl do hl[i][1]:=factor(hl[i][1]) od;
[seq(hl[i],i=1..nhl)];
end;

red:=proc(e) local e1;
e1:=expand(e);
while rem(e1,t,t)=0 do e1:=simplify(e1/t) od;
while rem(e1,1-t,t)=0 do e1:=simplify(e1/(1-t)) od;
e1;
end;

red1:=proc(e) local e1;
e1:=expand(e);
while rem(e1,t,t)=0 do e1:=simplify(e1/t) od;
while rem(e1,1-t,t)=0 do e1:=simplify(e1/(1-t)) od;
if member(e1,{2-t,t^2-t+1,t^2-2*t+2,2,3-2*t,t+1,t-t^2+1}) then e1:=1 fi;
e1;
end;

#group by tableaux with equivalence classes of rows; only record stuff if t-polynomial if very bad (reduce some 2-term stuff too)
comphlred:=proc(lam) local lam1,f,nl,ct,lc,ll,n,i,j,lt,k,r,lpt,dif,inv,iv,rt,nfl;global bc,sbc,hl,nhl;
n:=nops(lam);lam1:=cnjpart(lam);
r:=lch(lam1,n);lc:=r[1];ll:=r[2];
#print(lc,ll);
hl:=array(1..800000);nhl:=0;nfl:=0;lpt:=table();
for i from 1 to nops(bc) do
 j:=1;while (j<=n-1) and ((lam[j]<>lam[j+1]) or (bc[i][j]<bc[i][j+1])) do j:=j+1 od;
 if j=n then
   lt:=fold(lc,bc[i],ll);
   nfl:=nfl+nops(lt);
   for j from 1 to nops(lt) do
    rt:=conjtab(lt[j][1]);
    ct:=ordtab(rt);
    dif:=differ(lt[j][1]);
#    inv:=dinv(rt);
    k:=lpt[op(ct)];
#print(w1,pt);
    if op(0,k)<>`Integer` then 
     nhl:=nhl+1;
     if nhl mod 2000=0 then print(nhl,i) fi;
     lpt[op(ct)]:=nhl;
     hl[nhl]:=[t^((sbc[i]+lt[j][2]-lt[j][3])/2)*(1-t)^lt[j][3],ct,dif]
    else hl[k][1]:=expand(hl[k][1]+t^((sbc[i]+lt[j][2]-lt[j][3])/2)*(1-t)^lt[j][3]);
         hl[k][3]:=min(dif,hl[k][3]);
#hl[k][4]:=min(inv,hl[k][4]);
    fi;
   od;
  fi;
od;
nl:=[];
for i from 1 to nhl do 
# if expand(hl[i][1]-t^hl[i][4]*(1-t)^hl[i][3])<>0 then nl:=[op(nl),[factor(hl[i][1]),hl[i][4],hl[i][3],hl[i][2]]] fi 
inv:=dinv(hl[i][2]);if expand(hl[i][1]-t^inv*(1-t)^hl[i][3])<>0 then nl:=[op(nl),[factor(hl[i][1]),inv,hl[i][3],hl[i][2]]] fi 
od;
print(`terms`,nhl,`compression factor`,nfl/nhl*1.0,`errors`,nl);
end;

#conjugate a tableau in Japanese form
conjtab:=proc(t0) local t,i,t1,r;
t1:=[];t:=t0;
while t<>[] do
 r:=[];
 for i from nops(t) to 1 by -1 do
  r:=[t[i][1],op(r)];if nops(t[i])=1 then t:=subsop(i=NULL,t) else t[i]:=subsop(1=NULL,t[i]) fi;
 od;
# if pl=nops(r) then t1[-1]:={op(t1[-1]),r} else t1:=[op(t1),{r}];pl:=nops(r) fi;
 t1:=[op(t1),r]
od;
t1;
end;

cnjpart:=proc(p0) local p,p1,i;
p:=p0;while p[-1]=0 do p:=subsop(-1=NULL,p) od;
p1:=[];
while p<>[] do
 p1:=[op(p1),nops(p)];
 while p<>[] and p[-1]=1 do p:=subsop(-1=NULL,p) od;
 for i from 1 to nops(p) do p[i]:=p[i]-1 od;
od;
p1;
end;












switch2:=proc(mat) local i,j,c,m1,rev;
c:=links2(mat[1],mat[2]);m1:=mat;rev:=false;
for j from 1 to c do
 if mat[1,j]>mat[2,j] then m1[1][j]:=mat[2,j];m1[2][j]:=mat[1,j];rev:=true fi;
od;
m1,rev;
end;

switchall:=proc(mat0) local mat,i,j,rev,r;
i:=nops(mat0)-1;rev:=true;mat:=mat0;
while rev and i>=1 do
 rev:=false;
 for j from 1 to i do
  r:=switch2([mat[j],mat[j+1]]);
  if r[2] then rev:=true;mat[j]:=r[1][1];mat[j+1]:=r[1][2] fi
 od;
 i:=i-1;
od;
mat;
end;

links2:=proc(r1,r2) local i,t;
t:=true;i:=1;
while i<nops(r1) and t do if crossint([r1[i],r2[i]],[r1[i+1],r2[i+1]]) then i:=i+1 else t:=false fi od;
i;
end;

crossint:=proc(a0,b0) local a,b;
a:=[min(a0[1],a0[2]),max(a0[1],a0[2])];
b:=[min(b0[1],b0[2]),max(b0[1],b0[2])];
if a[2]>=b[2] and b[2]>=a[1] and a[1]>=b[1] then true else false fi;
end;

ordtab:=proc(t) local pl,i,j,t1,mat;
t1:=[];pl:=0;mat:=[];
for i from 1 to nops(t) do
 if nops(t[i])=pl then mat:=[op(mat),t[i]]
 else 
  if nops(mat)=1 then t1:=[op(t1),op(mat)] else if nops(mat)>1 then t1:=[op(t1),op(switchall(mat))] fi fi;
  mat:=[t[i]];pl:=nops(t[i]);
 fi;
od;
if nops(mat)=1 then t1:=[op(t1),op(mat)] else if nops(mat)>1 then t1:=[op(t1),op(switchall(mat))] fi fi;
t1;
end;

#number of inversions for tableaux by rows, NE justified
dinv:=proc(t) local i,i0,k,n;
n:=0;
for i from 2 to nops(t) do
 for k from 0 to nops(t[i])-1 do
  for i0 from 1 to i-1 do
   if t[i0][nops(t[i0])-k]>t[i][nops(t[i])-k] then
     if (k=nops(t[i])-1) or (t[i][nops(t[i])-k-1]>t[i0][nops(t[i0])-k]) then n:=n+1 fi
   fi
  od
 od
od;
n;
end;

#number of descents in rows for tableaux by columns, NE justified
differ:=proc(tt) local n,i,j;
n:=0;
for i from 1 to nops(tt)-1 do
 for j from 1 to nops(tt[i]) do
  if tt[i][j]<>tt[i+1][j] then n:=n+1 fi
 od
od;
n;
end;






nxt:=proc(p,pos) local i,j;
i:=pos[1];j:=pos[2];
if i<p[j] then i:=i+1 else
 if j<nops(p) then j:=j+1;i:=1 else i:=0;j:=0 fi
fi;
[i,j];
end;

prev:=proc(p,pos) local i,j;
i:=pos[1];j:=pos[2];
if i>1 then i:=i-1 else
 if j>1 then j:=j-1;i:=p[j] else i:=0;j:=0 fi
fi;
[i,j];
end;

stat:=proc(tt) local s,n,i,j,k;
n:=0;
for i from 1 to nops(tt)-1 do
 for j from 1 to nops(tt[i+1]) do
  if tt[i][j]<>tt[i+1][j] then n:=n+1 fi
 od
od;
s:=(1-t)^n;
n:=0;
for i from 1 to nops(tt) do
 for j from 1 to nops(tt[i])-1 do
  for k from j+1 to nops(tt[i]) do
   if (tt[i][j]>tt[i][k]) and ((i=nops(tt)) or (k>nops(tt[i+1])) or (tt[i][j]<tt[i+1][k])) then n:=n+1 fi;
  od
 od
od;
s*t^n
end;

tn:=proc(n) local i,s;
s:=1;
for i from 1 to n-1 do s:=s+t^i od;
s;
end;

fn:=proc(n)  local i,s;
s:=1;
for i from 1 to n do s:=s*tn(i) od;
s;
end;


tab2mon:=proc(t) local m,i,j;
m:=1;
for i from 1 to nops(t) do
 for j from 1 to nops(t[i]) do
  m:=m*cat(x,t[i][j]);
 od;
od;
end;

stat:=proc(tt) local s,n,i,j,k;
n:=0;
for i from 1 to nops(tt)-1 do
 for j from 1 to nops(tt[i+1]) do
  if tt[i][j]<>tt[i+1][j] then n:=n+1 fi
 od
od;
s:=(1-t)^n;
n:=0;
for i from 1 to nops(tt) do
 for j from 1 to nops(tt[i])-1 do
  for k from j+1 to nops(tt[i]) do
   if (tt[i][j]>tt[i][k]) and ((i=nops(tt)) or (k>nops(tt[i+1])) or (tt[i][j]<tt[i+1][k])) then n:=n+1 fi;
  od
 od
od;
s*t^n
end;

genhhl:=proc(p,n) local lt,la,i,j,adm,ct,cp,cp1,pp,hhl,nhhl,lpt,k,s,den,p1,rep,poss,la1,err;global nhl,hl;
cp:=[1,1];p1:=cnjpart(p);den:=1;rep:=0;
for i from 1 to nops(p) do
 if i=1 or p[i]=p[i-1] then rep:=rep+1
 else den:=den*fn(rep);rep:=1
 fi;
od;
den:=den*fn(rep);
lt:=[];la:=[{seq(i,i=1..n)}];ct:=[seq([0$p1[j]],j=1..nops(p1))];
while cp<>[0,0] do
 ct[cp[2]][cp[1]]:=la[-1][1];la[-1]:=subsop(1=NULL,la[-1]);
 cp1:=nxt(p1,cp);poss:=true;la1:=[];
 while cp1<>[0,0] and poss do
  adm:={seq(ct[cp1[2]][i],i=1..cp1[1]-1)};
  if cp1[2]>1 then pp:=ct[cp1[2]-1][cp1[1]];adm:=adm union {seq(ct[cp1[2]-1][i],i=1..cp1[1]-1)} else pp:=1 fi;
  adm:={seq(i,i=pp..n)} minus adm;
  if nops(adm)>0 then
   ct[cp1[2]][cp1[1]]:=adm[1];adm:=subsop(1=NULL,adm);la1:=[op(la1),adm];
   cp1:=nxt(p1,cp1);
  else poss:=false;
  fi;
 od;
 if poss then
  la:=[op(la),op(la1)];
  lt:=[op(lt),ct];
  cp:=[p1[-1],nops(p1)];
 else
  cp:=prev(p1,cp);la:=subsop(-1=NULL,la);
 fi;
 while cp<>[0,0] and la[-1]={} do la:=subsop(-1=NULL,la);cp:=prev(p1,cp) od;
od;

hhl:=array(1..800000);nhhl:=0;lpt:=table();
for i from 1 to nops(lt) do
 ct:=ordtab(conjtab(rev(lt[i])));
 k:=lpt[op(ct)];s:=stat(lt[i]);
 if op(0,k)<>`Integer` then 
  nhhl:=nhhl+1;
  if nhhl mod 2000=0 then print(nhhl) fi;
  lpt[op(ct)]:=nhhl;
  hhl[nhhl]:=[s,ct]
 else hhl[k][1]:=expand(hhl[k][1]+s)
 fi;
od;

err:=[];
if nhl=nhhl then
 for i from 1 to nhl do
  k:=lpt[op(hl[i][2])];
  if op(0,k)<>`Integer` then ERROR(`--- element not found`,hl[i][2])
  else if expand(hhl[k][1]-hl[i][1]*den)<>0 then err:=[op(err),[factor(hl[i][1]),factor(hhl[k][1]),hl[i][2]]] fi
  fi
 od
else ERROR(`different number of terms`);
fi;
print(`***`,err);
end;

#example of running: gensym(5);comphlred([3,2,2,0,0]);genhhl([3,2,2],5); - need to get [] in both cases as "error terms"
print(`for Maple 9,10; run gensym first`);