Software | University at Albany

Formal root polynomials in hyperbolic Schubert calculus

This Maple package checks the type G_2 Yang-Baxter relation for the elements h_i^Y(lambda) in Definition 3.1 of the paper "Towards generalized cohomology Schubert calculus via formal root polynomials" with my coauthor, K. Zainoulline.

The Maple package and a brief documentation

#Maple program checking the YB relation in type G_2 for the elements h_i^Y(lambda) in Definition 3.1 of our paper "Towards generalized cohomology Schubert calculus via formal root polynomials"

with(combinat);

#the hyperbolic formal group law (a=mu_1, b=mu_2)
f:=proc(x,y) (x+y-a*x*y)/(1+b*x*y) end;

#generates all subsets of {1,...,6} in order to expand the product of h_i(lambda) over the 6 positive roots of G_2
p6:=powerset(6);

#implements Y_i^2=a Y_i
sqr:=proc(s)
local s1, i, coef;
coef:=1;
s1:=[s[1]];
for i from 2 to nops(s) do
if s[i] <> s[i-1] then s1:=[op(s1), s[i]]
else coef:=coef*a
fi
od;
[coef, s1]
end;

#calculates a term in the product of h_i(lambda) corresponding to a subset "s" in "p6"
#a product "Y_1.Y_2.Y_1.Y_2" is represented as Y[1,2,1,2], with the appropriate coefficient (product of y_i=y[i], and possibly a, b)
term:=proc(s)
local a, b, i, ys, r;
if s=[] then 1
else
ys:=[];
for i to nops(s) do
if s[i] mod 2=0 then ys:=[op(ys), 2]
else ys:=[op(ys), 1]
fi
od;
r:=sqr(ys);
(-1)^nops(s)*product('y[s[i]]', 'i'=1 .. nops(s))*r[1]*Y[op(r[2])]
fi
end;

#calculates a term in the reversed product of h_i(lambda) corresponding to a subset "s" in "p6"
revterm:=proc(s)
local a, b, i, ys, r;
if s=[] then 1
else
ys:=[];
for i to nops(s) do
if s[i] mod 2=1 then ys:=[op(ys), 2]
else ys:=[op(ys), 1]
fi
od;
r:=sqr(ys);
(-1)^nops(s)*
product('y[7-s[i]]', 'i'=1 .. nops(s))*r[1]*Y[op(r[2])]
fi
end;

#expresses the variables y_i in terms of y_1=x and y_6=z corresponding to the simple roots alpha, beta
#we use the reflection order
#alpha=(1,-1,0),
#3 alpha+beta=(2,-1,-1),
#2 alpha+beta=(1,0,-1),
#3 alpha+2 beta=(1,1,-2),
#alpha+beta=(0,1,-1),
#beta=(-1,2,-1)
sbs:=proc(e)
simplify(subs({y[1]=x, y[6]=z,
y[2]=simplify(f(x, f(x, f(x,z)))),
y[3]=simplify(f(x, f(x, z))),
y[4]=simplify(f(x, f(x, f(x, f(z,z))))),
y[5]=simplify(f(x,z))}, e))
end;

#calculates the left-hand side in the Yang-Baxter relation
ls:=sum('term([op(p6[i])])','i'=1..nops(p6));

#calculates the right-hand side in the Yang-Baxter relation
rs:=sum('revterm([op(p6[i])])','i'=1..nops(p6));

#calculates the difference of the two sides, after applying the corresponding twisted braid relation in the Demazure algebra (cf. Leclerc, arXiv:1505.05097, Example 4.12 (ii))
dif:=expand(subs(Y[2,1,2,1,2,1]=Y[1,2,1,2,1,2]+4*b*Y[1,2,1,2]-4*b*Y[2,1,2,1]
+3*b^2*Y[1,2]-3*b^2*Y[2,1],ls-rs));

y[6] y[2] y[1] a Y[2, 1] - y[6] y[5] y[4] y[1] Y[2, 1, 2, 1]
+ y[6] y[3] y[1] a Y[2, 1]
- y[6] y[5] y[2] y[1] Y[2, 1, 2, 1]
+ y[5] y[3] y[2] a Y[1, 2] - y[6] y[4] y[3] y[1] a^2 Y[2, 1]
+ y[6] y[5] y[1] a Y[2, 1]
- y[5] y[4] y[3] y[2] Y[1, 2, 1, 2]
- 4 y[1] y[2] y[3] y[4] y[5] y[6] b Y[1, 2, 1, 2]
+ 4 y[1] y[2] y[3] y[4] y[5] y[6] b Y[2, 1, 2, 1]
- 3 y[1] y[2] y[3] y[4] y[5] y[6] b^2 Y[1, 2]
+ 3 y[1] y[2] y[3] y[4] y[5] y[6] b^2 Y[2, 1]
- y[6] y[5] y[3] y[1] a^2 Y[2, 1]
- y[6] y[4] y[2] y[1] a^2 Y[2, 1] + y[6] y[5] y[3] a Y[2, 1]
+ y[3] y[4] y[5] y[6] Y[1, 2, 1, 2] - y[6] y[3] Y[2, 1]
- y[1] y[3] y[4] y[5] y[6] a Y[1, 2, 1, 2]
+ y[5] y[4] y[2] a Y[1, 2]
+ y[1] y[4] y[5] y[6] Y[1, 2, 1, 2]
- y[1] y[2] y[4] y[5] y[6] a Y[1, 2, 1, 2]
+ y[5] y[6] Y[1, 2] - y[1] y[5] y[6] a Y[1, 2]
+ y[1] y[2] y[5] y[6] Y[1, 2, 1, 2] - y[5] y[4] Y[1, 2]
- y[3] y[5] y[6] a Y[1, 2] + y[1] y[3] y[5] y[6] a^2 Y[1, 2]
- y[1] y[2] y[3] y[5] y[6] a Y[1, 2, 1, 2]
+ y[1] y[6] Y[1, 2] - y[1] y[2] y[6] a Y[1, 2]
+ y[3] y[6] Y[1, 2] - y[1] y[3] y[6] a Y[1, 2]
+ y[1] y[2] y[3] y[6] Y[1, 2, 1, 2]
- y[1] y[4] y[6] a Y[1, 2] + y[1] y[2] y[4] y[6] a^2 Y[1, 2]
- y[3] y[4] y[6] a Y[1, 2] + y[1] y[3] y[4] y[6] a^2 Y[1, 2]
- y[1] y[2] y[3] y[4] y[6] a Y[1, 2, 1, 2]
- y[6] y[5] Y[2, 1] + y[1] y[2] Y[1, 2] + y[2] y[3] Y[2, 1]
+ y[6] y[5] y[4] y[3] y[1] a Y[2, 1, 2, 1]
+ y[1] y[4] Y[1, 2] - y[1] y[2] y[4] a Y[1, 2]
+ y[3] y[4] Y[1, 2] - y[1] y[3] y[4] a Y[1, 2]
+ y[1] y[2] y[3] y[4] Y[1, 2, 1, 2] + y[2] y[5] Y[2, 1]
- y[6] y[5] y[4] y[3] Y[2, 1, 2, 1]
- y[2] y[3] y[5] a Y[2, 1] + y[4] y[5] Y[2, 1]
- y[2] y[4] y[5] a Y[2, 1] + y[6] y[4] y[3] a Y[2, 1]
+ y[2] y[3] y[4] y[5] Y[2, 1, 2, 1] - y[3] y[4] Y[2, 1]
- y[2] y[3] Y[1, 2] - y[1] y[4] Y[2, 1]
+ y[1] y[3] y[4] a Y[2, 1] - y[1] y[2] Y[2, 1]
- y[5] y[2] Y[1, 2] + y[1] y[2] y[4] a Y[2, 1]
- y[1] y[2] y[3] y[4] Y[2, 1, 2, 1]
+ y[6] y[5] y[4] y[2] y[1] a Y[2, 1, 2, 1]
+ y[6] y[4] y[1] a Y[2, 1] - y[6] y[1] Y[2, 1]
+ y[1] y[2] y[3] y[4] y[6] a Y[2, 1, 2, 1]
- y[1] y[2] y[3] y[6] Y[2, 1, 2, 1]
+ y[1] y[2] y[3] y[5] y[6] a Y[2, 1, 2, 1]

#calculating the difference of the two sides in the Yang-Baxter relation after expressing all the y_i in terms of y_1 and y_6 corresponding to the simple roots
#as we can see from the above expression of "dif", only the coefficients of "Y_1.Y_2.Y_1.Y_2", "Y_2.Y_1.Y_2.Y_1", "Y_1.Y_2", and "Y_2.Y_1" need to be considered; so we only calculate these, and we get "0" each time, as expected.

sbs(coeff(dif,Y[1,2,1,2]));

sbs(coeff(dif,Y[2,1,2,1]));

0
sbs(coeff(dif,Y[1,2]));

sbs(coeff(dif,Y[2,1]));

Hall-Littlewood polynomials

This Maple package computes Hall-Littlewood polynomials of type A and C based on Schwer's formula. It also allows to experiment with compressing this formula to one in terms of tableaux, cf. the Haglund-Haiman-Loehr formula. C. L. was partially supported by National Science Foundation grant DMS-0701044.

The package for type A and a brief documentation

# 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;

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;

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`);

The package for type C and a brief documentation

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

#read `N:/My Documents/soft/alcoves/hall-litt-c.txt`:
#read `C:/Documents and Settings/lenart/My Documents/soft/alcoves/hall-litt-c.txt`:

print(`for Maple 9,10; run genhyp first`);

Perm2LengthB:=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 cmp(p[j],p[i]) then l:=l+1 fi
od
od;
l
end;

ListCode:=proc(n)
if (n=1) then
RETURN([[0], [1]])
fi;
RETURN(map(proc(code, n) local i;
op(code); seq([%, i], i=0..2*n-1)
end, ListCode(n-1), n))
end;

Code2Bar:=proc(code)
local
bar, # the result...
rd, # a reduced decomposition...
i, # variable for loop...
j; # variable for sequence...
rd:=Code2Rd(code);
bar:=[seq(j, j=1..nops(code))];
for i in rd do
if (i=0) then
bar:=subsop(1=-bar[1], bar)
else
bar:=subsop(i=bar[i+1], i+1=bar[i], bar)
fi
od;
RETURN(bar)
end;

Bar2Length:=proc(p) local s,i;s:=Perm2LengthB(p);for i from 1 to nops(p) do if p[i]<0 then s:=s+nops(p)+p[i]+1 fi od;s; end;

Code2Rd:=proc(code)
local
i, j; # variables for sequence...
[seq( [seq(-j, j=-i+1..-1), seq(j, j=0..i-1)] , i=1..nops(code))];
RETURN([seq(op(1..code[i], %[i]), i=1..nops(code))])
end;

genhyp:=proc(n) local i,b0;global bc,sbc;
b0:=ListCode(n);bc:=[seq(Code2Bar(b0[i]),i=1..nops(b0))];
sbc:=[seq(Bar2Length(bc[i]),i=1..nops(bc))];
print(`number of elements`,nops(bc));
end;

lchf:=proc(n,nn) local l,i,j,m1,m2,m3,m4;
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);
for i from n to 1 by -1 do
for j from i to 1 by -1 do l:=[[i,j],op(l)] od od;
m2:=nops(l);
for i from n to 1 by -1 do
for j from nn to n+1 by -1 do l:=[[j,i],op(l)] od od;
m3:=nops(l);
for i from n to 1 by -1 do
for j from i-1 to 1 by -1 do l:=[[i,j],op(l)] od od;
m4:=nops(l);
l,m4-m3,m4;
#l,m4-m3,m4-m2,m4-m1,m4;
#l,m4;
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]],[ll[-1][1]+r[3],lt[i]]];
# ll:=[op(ll),[ll[-1][1]+r[2],lt[i]],[ll[-1][1]+r[3],lt[i]],[ll[-1][1]+r[4],lt[i]],[ll[-1][1]+r[5],lt[i]]];
# ll:=[op(ll),[ll[-1][1]+r[2],lt[i]]];
od;
lc,ll;
end;

cmp:=proc(i,j) if i*j>0 then if i<j then true else false fi else if i>0 then true else false fi fi end;

decllong:=proc(w,p) local l1,l2;l1:=Bar2Length(w);l2:=Bar2Length(tr(w,p)); if l1>l2 then true else false fi end;

#### - take care that for ei+ej, entries [i,j] appear as [j,i]
decl:=proc(w,p) if p[1]<p[2] then cmp(w[p[2]],w[p[1]])
else if p[1]=p[2] then if w[p[1]]<0 then true else false fi
else cmp(-w[p[1]],w[p[2]]) fi fi
end;

tr:=proc(w,p) local w1;if p[1]<p[2] then tr1(w,p) else if p[1]=p[2] then w1:=w;w1[p[1]]:=-w1[p[1]];w1
else w1:=tr1(w,p);w1[p[1]]:=-w1[p[1]];w1[p[2]]:=-w1[p[2]];w1 fi fi end;

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

tr1:=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];
#if nops(lfld)>20000 then ERROR() fi;
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,Bar2Length(lw[-1]),nops(fld),fld]];
#*****lt:=[op(lt),[t,Bar2Length(lw[-1]),nops(fld),fld,lw[1],lw[-1]]];
od;
#print(`----`,lt);
lt;
end;

#returns only fld
fold1:=proc(lc,w) 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];
if nops(lfld)>22000 then ERROR() fi;
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;
[seq([lfld[i],Bar2Length(llw[i][-1])],i=1..nops(lfld))];
end;

#returns only fld; uses increasing condition on chains
fold2:=proc(lc,w) 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 not 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;
[seq([lfld[i],Bar2Length(llw[i][-1])],i=1..nops(lfld))];
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;

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

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

die1:=rand(1..2);with(combinat);

compress:=proc(p0) local p,i,j,aux,p1,pa;
p:=[seq(abs(p0[i]),i=1..nops(p0))];pa:=p;
for i from 1 to nops(p)-1 do
for j from i+1 to nops(p) do
if p[i]>p[j] then aux:=p[i];p[i]:=p[j];p[j]:=aux fi;
od
od;
p1:=[0$nops(p)];
for i from 1 to nops(p) do
j:=1;while pa[i]<>p[j] do j:=j+1 od;
p1[i]:=j*sign(p0[i]);
od;
p1
end;

#tests sums for last column only, length k; w is a permutation, of length n, corresponding to the last column; it's randomly generated;
tstlast:=proc(k,n) local w,lw,lc,ll,s,i,j,r,lt,ii,p1,p2;
r:=lch([k],n);lc:=r[1];ll:=r[2];
for i from 1 to ll[2][1] do lc:=subsop(1=NULL,lc) od;#ll[3][1]:=ll[3][1]-ll[2][1];ll:=subsop(2=NULL,ll);
print(lc,ll);
s:=0;
w:=randperm(n);
for i from 1 to n do j:=die1();if j=2 then w[i]:=-w[i] fi od;
print(w,nops(lc));
p1:=[seq(w[i],i=1..k)];p2:=compress([seq(w[i],i=k+1..n)]);
lw:=Bar2Length(w);
lt:=fold1(lc,w);
#print(lt);
for j from 1 to nops(lt) do
s:=s+expand(t^((lw+lt[j][2]-nops(lt[j][1]))/2)*(1-t)^nops(lt[j][1]))
od;
if s<>t^(Perm2LengthB(p1)+Bar2Length(p2)) then ERROR(`-----`,factor(s),nops(lt)) fi;
print(`***`,factor(s),nops(lt));
end;

weight:=proc(lc,lev,fld,p,lam) local i,nf,a,w,al,j,n,ps,w1;
nf:=nops(fld);w:=lam;n:=nops(lam);
for i from nf to 1 by -1 do
a:=lc[fld[i]];
if a[1]<a[2] then al:=[seq(0,j=1..a[1]-1),1,seq(0,j=a[1]+1..a[2]-1),-1,seq(0,j=a[2]+1..n)];ps:=w[a[1]]-w[a[2]]
else if a[1]>a[2] then al:=[seq(0,j=1..a[2]-1),1,seq(0,j=a[2]+1..a[1]-1),1,seq(0,j=a[1]+1..n)];ps:=w[a[2]]+w[a[1]]
else al:=[seq(0,j=1..a[1]-1),2,seq(0,j=a[1]+1..n)];ps:=w[a[1]] fi fi;
w:=w+(lev[fld[i]]-ps)*al;
od;
w1:=[0$n];for i from 1 to n do w1[abs(p[i])]:=sign(p[i])*w[i] od;
w1;
end;

refl:=proc(a,b,a1);
if b<0 then if a1=a then -b else if a1=-b then a else a1 fi fi
else if a1=a then -b else if a1=b then -a else a1 fi fi
fi;
end;

posdir:=proc(lc,fld,w)
local n,i,j,flc,ii,nn,a,b,a1,b1,np,x;
nn:=nops(lc);flc:=[[0,0]$nn];
for i from 1 to nn do
if lc[i][1]<lc[i][2] then flc[i]:=[lc[i][1],-lc[i][2]] else flc[i]:=[lc[i][2],lc[i][1]] fi
od;
#fold chain
for i from nops(fld) to 1 by -1 do
ii:=fld[i];a:=flc[ii][1];b:=flc[ii][2];
flc[ii]:=[0,0];
#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]<>0 then
a1:=sign(flc[j][1])*refl(a,b,abs(flc[j][1]));b1:=sign(flc[j][2])*refl(a,b,abs(flc[j][2]));
if abs(a1)<abs(b1) then flc[j]:=[a1,b1] else flc[j]:=[b1,a1] fi;
fi;
od;
od;
np:=0;
for i from 1 to nn do
a:=flc[i][1];
if a<>0 then
b:=flc[i][2];
x:=[sign(a)*w[abs(a)],sign(b)*w[abs(b)]]; if abs(x[1])>abs(x[2]) then x:=[x[2],x[1]] fi;
if x[1]>0 then np:=np+1 fi
fi;
od;
np;
end;

levels:=proc(lc) local i,j,lev;
lev:=[1$nops(lc)];
for i from 2 to nops(lc) do
j:=i-1;
while (j>=1) and (lc[j]<>lc[i]) do j:=j-1 od;
if j>=1 then lev[i]:=lev[j]+1 fi;
od;
lev;
end;

isdom:=proc(w) local i;
i:=1;
while (i<=nops(w)-1) and (w[i]>=w[i+1]) do i:=i+1 od;
if i=nops(w) and w[-1]>=0 then true else false fi;
end;

#checks if schwer's formula agrees with ram's formula
compschwer:=proc(lam) local lev,lam1,lc,n,i,hl,j,lt,k,nhl,r,lpt,w,nz,hl1,lpt1,nhl1,lc1,lev1,ps,kk,nfl;global bc,sbc;
n:=nops(lam);lam1:=cnjpart(lam);nz:=0;i:=n;while lam[i]=0 do i:=i-1;nz:=nz+1 od;nfl:=0;
r:=lch(lam1,n);lc:=r[1];lev:=levels(lc);
lc1:=[[2,1],[3,1],[4,1],[1,1],[1,4],[1,3],[2,1],[3,2],[3,1],[4,2],[4,1],[2,2],[2,1],[1,1],[2,4],[1,4],[3,2],[3,1],[2,1]];lev1:=levels(lc1);
hl:=array(1..1100000);nhl:=0;lpt:=table();hl1:=array(1..1100000);nhl1:=0;lpt1:=table();
for i from 1 to nops(bc) do
j:=1;while (j<=n-1) and ((lam[j]<>lam[j+1]) or (cmp(bc[i][j],bc[i][j+1]))) do j:=j+1 od;
if j=n then
while lam[j]=0 and bc[i][j]>0 do j:=j-1 od;
if lam[j]<>0 then
lt:=fold1(lc,bc[i]);
for j from 1 to nops(lt) do
nfl:=nfl+1;
if nfl mod 10000=0 then print(nfl/1000,i) fi;
#if sbc[i]+lt[j][2]-lt[j][3]<0 then ERROR(lc,bc[i],ll,nops(lt),j) fi;
w:=weight(lc,lev,lt[j][1],bc[i],lam);
if isdom(w) then
k:=lpt[op(w)];
#print(w1,pt);
if op(0,k)<>`Integer` then
nhl:=nhl+1;
lpt[op(w)]:=nhl;
hl[nhl]:=[t^((sbc[i]+lt[j][2]-nops(lt[j][1]))/2)*(1-t)^nops(lt[j][1]),w]
else hl[k][1]:=expand(hl[k][1]+t^((sbc[i]+lt[j][2]-nops(lt[j][1]))/2)*(1-t)^nops(lt[j][1]));
fi;
fi;
od;
lt:=fold1(lc1,bc[i]);
for j from 1 to nops(lt) do
w:=weight(lc1,lev1,lt[j][1],bc[i],lam);
if isdom(w) then
ps:=sum((lam[kk]+w[kk])*(n-kk+1/2),kk=1..n);k:=lpt1[op(w)];
#print(w1,pt);
if op(0,k)<>`Integer` then
nhl1:=nhl1+1;
lpt1[op(w)]:=nhl1;
hl1[nhl1]:=[t^(ps+nz^2-n^2+sbc[i]-nops(lt[j][1])-posdir(lc1,lt[j][1],bc[i]))*(1-t)^nops(lt[j][1]),w]
else hl1[k][1]:=expand(hl1[k][1]+t^(ps+nz^2-n^2+sbc[i]-nops(lt[j][1])-posdir(lc1,lt[j][1],bc[i]))*(1-t)^nops(lt[j][1]));
fi;
fi;
od;
fi;
fi;
od;
if nhl<>nhl1 then ERROR(`different number of weights`,nhl,nhl1) fi;
for i from 1 to nhl do
k:=lpt1[op(hl[i][2])];if op(0,k)<>`Integer` then ERROR(`missing weight in Schwer`,hl[i]) fi;
if hl[i][1]<>hl1[k][1] then ERROR(`different coefficients`,hl[i],hl1[k]) fi;
od;
print(`*** Checked ***`);
end;

#compression in Schwer's formula
compschwer1:=proc(lam) local hl1,lc,ll,n,i,hl,j,lt,k,nhl,r,lpt,nfl,hl2,hl3,lev,nz,w,ps;global bc,sbc;
n:=nops(lam);nz:=0;i:=n;while lam[i]=0 do i:=i-1;nz:=nz+1 od;
lc:=[[2,1],[3,1],[4,1],[1,1],[1,4],[1,3],[2,1],[3,2],[3,1],[4,2],[4,1],[2,2],[2,1],[1,1],[2,4],[1,4],[3,2],[3,1],[2,1]];
ll:=[[0,0],[0,1],[6,1],[7,2],[16,2],[19,3],[19,3]];
lev:=levels(lc);
hl:=array(1..80000000);nhl:=0;lpt:=table();nfl:=0;
for i from 1 to nops(bc) do
j:=1;while (j<=n-1) and ((lam[j]<>lam[j+1]) or (cmp(bc[i][j],bc[i][j+1]))) do j:=j+1 od;
if j=n then
while lam[j]=0 and bc[i][j]>0 do j:=j-1 od;
if lam[j]<>0 then
lt:=fold(lc,bc[i],ll);
nfl:=nfl+nops(lt);
for j from 1 to nops(lt) do
w:=weight(lc,lev,lt[j][-1],bc[i],lam);
if isdom(w) then
ps:=sum((lam[kk]+w[kk])*(n-kk+1/2),kk=1..n);
k:=lpt[op(lt[j][1])];
if op(0,k)<>`Integer` then
nhl:=nhl+1;
if nhl mod 2000=0 then print(nhl,i) fi;
lpt[op(lt[j][1])]:=nhl;
hl[nhl]:=[t^(ps+nz^2-n^2+sbc[i]-nops(lt[j][-1])-posdir(lc,lt[j][-1],bc[i]))*(1-t)^nops(lt[j][-1]),lt[j][1],1]
else hl[k][1]:=expand(hl[k][1]+t^(ps+nz^2-n^2+sbc[i]-nops(lt[j][-1])-posdir(lc,lt[j][-1],bc[i]))*(1-t)^nops(lt[j][-1]));
hl[k][3]:=hl[k][3]+1;
fi;
fi;
od;
fi;
fi;
od;
hl3:=[];
hl1:=[];hl2:=[];k:=0;
for i from 1 to nhl do
# hl3:=[op(hl3),[factor(hl[i][1]),hl[i][3],hl[i][2]]];
r:=red(hl[i][1]); if r<>1 then k:=k+1;
hl2:=[op(hl2),[factor(hl[i][1]),hl[i][3],hl[i][2]]]
#****hl[i][1]:=factor(hl[i][1]);
#hl1:=[op(hl1),[factor(hl[i][1]),hl[i][3],hl[i][2]]]
#hl1:=[op(hl1),[r+t-1,hl[i][3],hl[i][2]]]
fi od;
print(`compression factor`,nfl/nhl*1.0,`terms`,nhl,`of which`,k,`that is`,k/nhl*100.0,`per cent in nonstandard form`);
hl2;
end;

#USE tab();
#lam as [3,2,1,0,0,0], lam1 is the transpose - only nonzero column lengths
comphl1:=proc(lam) local hl1,lam1,lc,ll,n,i,hl,j,lt,k,nhl,r,lpt,nfl,hl2,hl3;global bc,sbc;
n:=nops(lam);lam1:=cnjpart(lam);
r:=lch(lam1,n);lc:=r[1];ll:=r[2];
hl:=array(1..800000000);nhl:=0;lpt:=table();nfl:=0;
for i from 1 to nops(bc) do
j:=1;while (j<=n-1) and ((lam[j]<>lam[j+1]) or (cmp(bc[i][j],bc[i][j+1]))) do j:=j+1 od;
if j=n then
while lam[j]=0 and bc[i][j]>0 do j:=j-1 od;
if lam[j]<>0 then
lt:=fold(lc,bc[i],ll);
nfl:=nfl+nops(lt);
for j from 1 to nops(lt) do
#if sbc[i]+lt[j][2]-lt[j][3]<0 then ERROR(lc,bc[i],ll,nops(lt),j) fi;
k:=lpt[op(lt[j][1])];
#print(w1,pt);
if op(0,k)<>`Integer` then
nhl:=nhl+1;
if nhl mod 2000=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],1]
#*****hl[nhl]:=[t^((sbc[i]+lt[j][2]-lt[j][3])/2)*(1-t)^lt[j][3],lt[j][1],1,[t^((sbc[i]+lt[j][2]-lt[j][3])/2)*(1-t)^lt[j][3],lt[j][-3],lt[j][-2],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]);hl[k][3]:=hl[k][3]+1;
#*****if hl[k][3]<=13 then hl[k]:=[op(hl[k]),[t^((sbc[i]+lt[j][2]-lt[j][3])/2)*(1-t)^lt[j][3],lt[j][-3],lt[j][-2],lt[j][-1]]] fi
fi;
od;
fi;
fi;
od;
hl3:=[];
hl1:=[];hl2:=[];k:=0;
for i from 1 to nhl do
## hl3:=[op(hl3),[factor(hl[i][1]),hl[i][3],hl[i][2]]];
r:=red(hl[i][1]); if r<>1 then k:=k+1;
hl2:=[op(hl2),[factor(hl[i][1]),hl[i][3],hl[i][2]]]
# if redt(r+t-1)<>1 then hl2:=[op(hl2),[factor(hl[i][1]),hl[i][3],hl[i][2]]]
##****hl[i][1]:=factor(hl[i][1]);
##hl1:=[op(hl1),[factor(hl[i][1]),hl[i][3],hl[i][2]]]
##hl1:=[op(hl1),[r+t-1,hl[i][3],hl[i][2]]]
#fi
fi od;
print(`compression factor`,nfl/nhl*1.0,`terms`,nhl,`of which`,k,`that is`,k/nhl*100.0,`per cent in nonstandard form`);
nhl,nfl/nhl*1.0,hl2,seq(hl[i],i=1..nhl);
end;

rev:=proc(s) local i; [seq(s[nops(s)+1-i],i=1..nops(s))] 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;s:=1;
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 cmp(tt[i][k],tt[i][j]) and ((i=nops(tt)) or (k>nops(tt[i+1])) or cmp(tt[i][j],tt[i+1][k])) then n:=n+1 fi;
od
od
od;
s*t^n
end;

tst:=proc(r) local tab,s,s1,i;
for i from 1 to nops(r) do
tab:=r[i][3];
if tab[1]=tab[2] and tab[3]=tab[4] and tab[5]=tab[6] then
s:=stat(rev(tab));s1:=powt(r[i][1]);
if expand(s-s1)<>0 then print(s,factor(r[i][1]),r[i][2],tab) fi
fi;
od;
end;

q-weight multiplicities for the superalgebras spo(2n,M)

This Maple package computes Lusztig's q-weight multiplicities for the semisimple Lie algebras of type B-D, as well as the q-weight multiplicities for the Lie superalgebras spo(2n,M) introduced in a joint paper with C. Lecouvey. C. L. was partially supported by National Science Foundation grant DMS-0701044.

The package and a brief documentation

# NOTE: for Maple 9 and 10.
# works well for osp(2m+1,2n) with m<=4 and n<=3
# NOTE: in a pair (lam1,lam2), the first partition is the orthogonal one, and the second the symplectic one
# istyp(lam1,lam2) - verifies if the dominant pair (lam1,lam2) is typical for lam1 of type B
# istyp2(lam1,lam2) - verifies if the dominant pair (lam1,lam2) is typical for lam1 of type D
# isfte(lam1,lam2) - verifies if the dominant pair (lam1,lam2) corresponds to a finite dimensional irrep.
# rhobc(m,n), rhodc(m,n) - the corresponding rho=rho^+-rho^-
# klmb, klmc, klmd(lam,mu) - Lusztig's q-weight multiplicities
# num(m,n) - computes the Kac numerator for type B/C; has to be executed before any subsequent klm with partitions of length m,n
# num2(m,n) - computes the Kac numerator for type D/C; has to be executed before any subsequent klm2 with partitions of length m,n
# klm(lam1,lam2,mu1,mu2) - computes the corresponding q-weight multiplicity for type B/C
# klm2(lam1,lam2,mu1,mu2) - computes the corresponding q-weight multiplicity for type D/C
# num0, klm0, istyp0 - similar procedures for gl(m,n)
# alldom(lam1,lam2) - computes the weights of the irreps of g^0 in the irrep of g of type B/C with weight (lam1,lam2); has to be
# executed before klmp
# klmp(lam1,lam2,mu1,mu2) - computes the different q-analog sum m(lambda,gamma) K_q^{g_0}(gamma,mu); good for verifications
# (for q=1 need to get the same result as klm); NOTE: klmp always produces a polynomial in N[q]
# verifb,verifc,verifd(w) - for verification purposes only; compares all the weight multiplicities in an irrep of classical algebra
# obtained with Stembridge's package (needs to be loaded) with the ones obtained via the procedure for generating partitions
# in this package; w are the coordinates of a dominant weight in the basis of fundamental weights (omega_1=e_1,omega_2=e1+e2,...)

# EXAMPLE in the paper:
#num(3,3); (should return the number 70592, i.e. the number of terms in the expansion of the Kac numerator).
#istyp([3,2,0],[5,4,4]);
#isfte([3,2,0],[5,4,4]);
#klm([3,2,0],[5,4,4],[1,1,0],[3,2,1]); (runs around 21 seconds on a 2.8GHz processor)
#klm([3,2,0],[10,9,9],[1,1,0],[8,7,6]); (this is the stabilized K)

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[i]>p[j] then l:=l+1 fi
od
od;
l
end;

ListCode:=proc(n)
if (n=1) then
RETURN([[0], [1]])
fi;
RETURN(map(proc(code, n) local i;
op(code); seq([%, i], i=0..2*n-1)
end, ListCode(n-1), n))
end;

Code2Length:=proc(code)
RETURN(convert(code, `+`))
end;

Code2Rd:=proc(code)
local
i, j; # variables for sequence...
[seq( [seq(-j, j=-i+1..-1), seq(j, j=0..i-1)] , i=1..nops(code))];
RETURN([seq(op(1..code[i], %[i]), i=1..nops(code))])
end;

cmbc:=proc(n,a) local a1,ch,j,i,c;
a1:=0;c:=[];
while a1<=a do
ch:=choose(a+2*n-3-a1,2*n-3);
c:=[op(c),seq([a1/2,a1,seq(ch[i][j]-j+a1,j=1..2*n-3),a],i=1..nops(ch))];
a1:=a1+2;
od;
c;
end;

cmba:=proc(n,a) local a1,ch,j,i,c;
ch:=choose(a+n-2,n-2);
c:=[seq([0,seq(ch[i][j]-j,j=1..n-2),a],i=1..nops(ch))];
c;
end;

cmbb:=proc(n,a) local a1,ch,j,i,c;
a1:=0;c:=[];
while a1<=a do
ch:=choose(a+2*n-3-a1,2*n-3);
c:=[op(c),seq([0,a1,seq(ch[i][j]-j+a1,j=1..2*n-3),a],i=1..nops(ch))];
a1:=a1+1;
od;
c;
end;

cmbd:=proc(n,a) local a1,ch,j,i,c;
ch:=choose(a+2*n-3,2*n-3);
c:=[seq([0,seq(ch[i][j]-j,j=1..2*n-3),a],i=1..nops(ch))];
c;
end;

genpd:=proc(lam) local i,n,c,s,j,ex,lam1,cc;
n:=nops(lam);
if n=2 then if (lam[1]>=abs(lam[2])) and ((lam[1]+lam[2]) mod 2 =0) then ex:=q^(lam[1]) else ex:=0 fi
else
i:=1;s:=0;while (i<=n) and (s>=0) do s:=s+lam[i];i:=i+1 od;
if (s<0) or (s mod 2=1) then ex:=0
else
ex:=0;
c:=cmbd(n,lam[1]);
for i from 1 to nops(c) do
cc:=c[i];
lam1:=subsop(1=NULL,lam);
for j from 1 to n-1 do
lam1[j]:=lam1[j]+2*cc[2*j]-cc[2*j-1]-cc[2*j+1];
od;
ex:=ex+expand(genpd(lam1)*q^(cc[2*n-1]))
od;
fi fi;
ex;
end;

genpc:=proc(lam) local i,n,c,s,j,ex,lam1,cc;
n:=nops(lam);
if n=1 then if (lam[1]>=0) and (lam[1] mod 2 =0) then ex:=q^(lam[1]/2) else ex:=0 fi
else
i:=1;s:=0;while (i<=n) and (s>=0) do s:=s+lam[i];i:=i+1 od;
if (s<0) or (s mod 2=1) then ex:=0
else
ex:=0;
c:=cmbc(n,lam[1]);
for i from 1 to nops(c) do
cc:=c[i];
lam1:=subsop(1=NULL,lam);
for j from 1 to n-1 do
lam1[j]:=lam1[j]+2*cc[2*j+1]-cc[2*j]-cc[2*j+2];
od;
ex:=ex+expand(genpc(lam1)*q^(cc[2*n]-cc[1]))
od;
fi fi;
ex;
end;

genpa:=proc(lam) local i,n,c,s,j,ex,lam1,cc;
n:=nops(lam);
if n=2 then if (lam[1]>=0) then ex:=q^(lam[1]) else ex:=0 fi
else
i:=1;s:=0;while (i<=n) and (s>=0) do s:=s+lam[i];i:=i+1 od;
if (s<0) or (s>0) then ex:=0
else
ex:=0;
c:=cmba(n,lam[1]);
for i from 1 to nops(c) do
cc:=c[i];
lam1:=subsop(1=NULL,lam);
for j from 1 to n-1 do
lam1[j]:=lam1[j]+cc[j+1]-cc[j];
od;
ex:=ex+expand(genpa(lam1)*q^(cc[n]))
od;
fi fi;
ex;
end;

genpb:=proc(lam) local i,n,c,s,j,ex,lam1,cc;
n:=nops(lam);
if n=1 then if lam[1]>=0 then ex:=q^lam[1] else ex:=0 fi
else
i:=1;s:=0;while (i<=n) and (s>=0) do s:=s+lam[i];i:=i+1 od;
if s<0 then ex:=0
else
ex:=0;
c:=cmbb(n,lam[1]);
for i from 1 to nops(c) do
cc:=c[i];
lam1:=subsop(1=NULL,lam);
for j from 1 to n-1 do
lam1[j]:=lam1[j]+2*cc[2*j+1]-cc[2*j]-cc[2*j+2];
od;
ex:=ex+expand(genpb(lam1)*q^cc[2*n])
od;
fi fi;
ex;
end;

acta:=proc(p,w) local i,w1,n;
n:=nops(w);w1:=[0$n];
for i from 1 to n do w1[p[i]]:=w[i] od;
w1;
end;

actb:=proc(p,w) local i,w1,n;
n:=nops(w);w1:=[0$n];
for i from 1 to n do w1[abs(p[i])]:=sign(p[i])*w[i] od;
w1;
end;

rhoaa:=proc(m,n) local i,s;s:=[];for i from 1 to m do s:=[op(s),n] od;
for i from 1 to n do s:=[op(s),-m] od;[op(rhoa(m)),op(rhoa(n))]+s end;
rhoa:=proc(n) local i;[seq((n-1)/2-i+1,i=1..n)] end;
rhoc:=proc(n) local i;[seq(n+1-i,i=1..n)] end;
rhob:=proc(n) local i;[seq(n-i+1/2,i=1..n)] end;
rhod:=proc(n) local i;[seq(n-i,i=1..n)] end;
rhobc:=proc(m,n) local i,s;s:=[];for i from 1 to m do s:=[op(s),0] od;
for i from 1 to n do s:=[op(s),-m-1/2] od;[op(rhob(m)),op(rhoc(n))]+s end;
rhodc:=proc(m,n) local i,s;s:=[];for i from 1 to m do s:=[op(s),0] od;
for i from 1 to n do s:=[op(s),-m] od;[op(rhod(m)),op(rhoc(n))]+s end;

istyp:=proc(lam1,lam2) local m,n,i,j,w,f;
m:=nops(lam1);n:=nops(lam2);
w:=[op(lam1),op(lam2)]+rhobc(m,n);f:=false;
j:=1;
while (j<=n) and (not f) do
i:=1;
while (i<=m) and (not f) do
if (w[i]=w[j+m]) or (w[i]=-w[j+m]) then f:=true;#print(i,j,w)
fi;
i:=i+1
od;
j:=j+1;
od;
not f;
end;

istyp0:=proc(lam1,lam2) local m,n,i,j,w,f;
m:=nops(lam1);n:=nops(lam2);
w:=[op(lam1),op(lam2)]+rhoaa(m,n);f:=false;
j:=1;
while (j<=n) and (not f) do
i:=1;
while (i<=m) and (not f) do
if (w[i]=-w[j+m]) then f:=true;#print(i,j,w)
fi;
i:=i+1
od;
j:=j+1;
od;
not f;
end;

isfte:=proc(lam1,lam2) local i,j,m;
i:=1;m:=nops(lam1);
while (i<=m) and (lam1[i]>0) do i:=i+1 od;
if i-1<=lam2[-1] then true else false fi;
end;

istyp2:=proc(lam1,lam2) local m,n,i,j,w,f;
m:=nops(lam1);n:=nops(lam2);
w:=[op(lam1),op(lam2)]+rhodc(m,n);f:=false;
i:=1;
while (i<=m) and (not f) do
j:=1;
while (j<=n) and (not f) do
if (w[i]=w[j+m]) or (w[i]=-w[j+m]) then f:=true;#print(i,j,w)
fi;
j:=j+1
od;
i:=i+1
od;
not f;
end;

klmam:=proc(lam,mu) local s,i,n,mur,lamr,no,ss;global aam,saam;
n:=nops(lam);
ss:=sum(lam[i]-mu[i],i=1..n);
if ss <>0 then 0 else
s:=0;
mur:=rhoa(n)+mu;#+[ss/n$n];
lamr:=rhoa(n)+lam;
#print(mur);
for i from 1 to nops(aam) do
no:=genpa(acta(aam[i],lamr)-mur);
#if no>0 then print(bc[i],actb(bc[i],lamr)-mur,sbc[i]*no) fi;
s:=s+saam[i]*no od;
s;
fi;
end;

klman:=proc(lam,mu) local s,i,n,mur,lamr,no,ss;global aan,saan;
n:=nops(lam);
ss:=sum(lam[i]-mu[i],i=1..n);
if ss <>0 then 0 else
s:=0;
mur:=rhoa(n)+mu;#+[ss/n$n];
lamr:=rhoa(n)+lam;
#print(mur);
for i from 1 to nops(aan) do
no:=genpa(acta(aan[i],lamr)-mur);
#if no>0 then print(bc[i],actb(bc[i],lamr)-mur,sbc[i]*no) fi;
s:=s+saan[i]*no od;
s;
fi;
end;

klmc:=proc(lam,mu) local s,i,n,mur,lamr,no;global bc,sbc;
n:=nops(lam);
s:=0;mur:=rhoc(n)+mu;lamr:=rhoc(n)+lam;
#print(mur);
for i from 1 to nops(bc) do
no:=genpc(actb(bc[i],lamr)-mur);
#if no>0 then print(bc[i],actb(bc[i],lamr)-mur,sbc[i]*no) fi;
s:=s+sbc[i]*no od;
s;
end;

klmb:=proc(lam,mu) local s,i,n,mur,lamr,no;global bb,sbb;
n:=nops(lam);
s:=0;mur:=rhob(n)+mu;lamr:=rhob(n)+lam;
#print(mur);
for i from 1 to nops(bb) do
no:=genpb(actb(bb[i],lamr)-mur);
#if subs(q=1,no)>0 then
#print(b[i],actb(bb[i],lamr)-mur,sbb[i]*no);
#fi;
s:=s+sbb[i]*no od;
s;
end;

klmd:=proc(lam,mu) local s,i,n,mur,lamr,no;global bd,sbd;
n:=nops(lam);
s:=0;mur:=rhod(n)+mu;lamr:=rhod(n)+lam;
#print(mur);
for i from 1 to nops(bd) do
no:=genpd(actb(bd[i],lamr)-mur);
#if no>0 then print(bd[i],actb(bd[i],lamr)-mur,sbd[i]*no) fi;
s:=s+sbd[i]*no od;
s;
end;

klmc1:=proc(lam,mu,k) local s,i,n,mur,lamr,no;global bc,sbc;
n:=nops(lam);
s:=0;mur:=rhoc(n)+mu+[(-k-1/2)$n];lamr:=rhoc(n)+lam+[(-k-1/2)$n];
#print(mur);
for i from 1 to nops(bc) do
no:=genpc(actb(bc[i],lamr)-mur);
#if no>0 then print(bc[i],actb(bc[i],lamr)-mur,sbc[i]*no) fi;
s:=s+sbc[i]*no od;
s;
end;

gensymm:=proc(n) local i,b0;global aam,saam;
aam:=[op(permute(n))];
saam:=[seq((-1)^Perm2Length(aam[i]),i=1..nops(aam))];
print(`number of elements`,nops(aam));
end;

gensymn:=proc(n) local i,b0;global aan,saan;
aan:=[op(permute(n))];
saan:=[seq((-1)^Perm2Length(aan[i]),i=1..nops(aan))];
print(`number of elements`,nops(aan));
end;

genhypb:=proc(n) local i,b0;global bb,sbb;
b0:=ListCode(n);bb:=[seq(Code2Bar(b0[i]),i=1..nops(b0))];
sbb:=[seq((-1)^Code2Length(b0[i]),i=1..nops(bb))];
print(`number of elements`,nops(bb));
end;

genhypc:=proc(n) local i,b0;global bc,sbc;
b0:=ListCode(n);bc:=[seq(Code2Bar(b0[i]),i=1..nops(b0))];
sbc:=[seq((-1)^Code2Length(b0[i]),i=1..nops(bc))];
print(`number of elements`,nops(bc));
end;

genhypd:=proc(n) local i,b0,bc1,c,j;global bd,sbd;
b0:=ListCode(n);bc1:=[seq(Code2Bar(b0[i]),i=1..nops(b0))];sbd:=[];bd:=[];
for i from 1 to nops(bc1) do
c:=0;
for j from 1 to n do if bc1[i][j]<0 then c:=c+1 fi od;
if c mod 2 =0 then bd:=[op(bd),bc1[i]];sbd:=[op(sbd),(-1)^Code2Length(b0[i])] fi
od;
print(`number of elements`,nops(bd));
end;

#reinitialize each time because of array of pointers lp
#lw=list of weights (kappa1,2); lc=list of coefficients;lp=list of pointers
num:=proc(m,n) global lp,lc,lw,nn,lim;local i,lc1,j,k,w,w1,pt,nn0;
lim:=1000000;
lc:=array(1..lim);lw:=array(1..lim);lc1:=array(1..lim);lp:=table();
nn:=1;
genhypb(m);genhypc(n);
lw[1]:=[0$(m+n)];
lp[op(lw[1])]:=1;
lc[1]:=1;
for i from 1 to m do
for j from 1 to n do
for k from 1 to nn do lc1[k]:=lc[k] od;
w:=[0$(i-1),-1,0$(m-i),0$(j-1),1,0$(n-j)];
nn0:=nn;
for k from 1 to nn0 do
w1:=w+lw[k];
#print(w1);
pt:=lp[op(w1)];
#print(w1,pt);
if op(0,pt)=`Integer` then
lc[pt]:=lc[pt]+lc1[k] else nn:=nn+1;if nn>lim then ERROR(`too long array`) fi;lc[nn]:=lc1[k];lw[nn]:=w1;lp[op(w1)]:=nn; fi;
od;
for k from 1 to nn do lc1[k]:=lc[k] od;
w:=[0$(i-1),1,0$(m-i),0$(j-1),1,0$(n-j)];
nn0:=nn;
for k from 1 to nn0 do
w1:=w+lw[k];
pt:=lp[op(w1)];
if op(0,pt)=`Integer` then
lc[pt]:=lc[pt]+lc1[k] else nn:=nn+1;if nn>lim then ERROR(`too long array`) fi;lc[nn]:=lc1[k];lw[nn]:=w1;lp[op(w1)]:=nn; fi;
od;
od od;
for j from 1 to n do
for k from 1 to nn do lc1[k]:=lc[k] od;
w:=[0$m,0$(j-1),1,0$(n-j)];
nn0:=nn;
for k from 1 to nn0 do
w1:=w+lw[k];
pt:=lp[op(w1)];
if op(0,pt)=`Integer` then
lc[pt]:=lc[pt]+lc1[k] else nn:=nn+1;if nn>lim then ERROR(`too long array`) fi;lc[nn]:=lc1[k];lw[nn]:=w1;lp[op(w1)]:=nn; fi;
od;
od;
end;

num0:=proc(m,n) global lp,lc,lw,nn,lim;local i,lc1,j,k,w,w1,pt,nn0;
lim:=1000000;
lc:=array(1..lim);lw:=array(1..lim);lc1:=array(1..lim);lp:=table();
nn:=1;
gensymm(m);gensymn(n);
lw[1]:=[0$(m+n)];
lp[op(lw[1])]:=1;
lc[1]:=1;
for i from 1 to m do
for j from 1 to n do
for k from 1 to nn do lc1[k]:=lc[k] od;
w:=[0$(i-1),-1,0$(m-i),0$(j-1),1,0$(n-j)];
nn0:=nn;
for k from 1 to nn0 do
w1:=w+lw[k];
#print(w1);
pt:=lp[op(w1)];
#print(w1,pt);
if op(0,pt)=`Integer` then
lc[pt]:=lc[pt]+lc1[k] else nn:=nn+1;if nn>lim then ERROR(`too long array`) fi;lc[nn]:=lc1[k];lw[nn]:=w1;lp[op(w1)]:=nn; fi;
od;
od od;
end;

#reinitialize each time because of array of pointers lp
#lw=list of weights (kappa1,2); lc=list of coefficients;lp=list of pointers
num2:=proc(m,n) global lp,lc,lw,nn,lim;local i,lc1,j,k,w,w1,pt,nn0;
lim:=1000000;
lc:=array(1..lim);lw:=array(1..lim);lc1:=array(1..lim);lp:=table();
nn:=1;
genhypd(m);genhypc(n);
lw[1]:=[0$(m+n)];
lp[op(lw[1])]:=1;
lc[1]:=1;
for i from 1 to m do
for j from 1 to n do
for k from 1 to nn do lc1[k]:=lc[k] od;
w:=[0$(i-1),-1,0$(m-i),0$(j-1),1,0$(n-j)];
nn0:=nn;
for k from 1 to nn0 do
w1:=w+lw[k];
#print(w1);
pt:=lp[op(w1)];
#print(w1,pt);
if op(0,pt)=`Integer` then
lc[pt]:=lc[pt]+lc1[k] else nn:=nn+1;if nn>lim then ERROR(`too long array`) fi;lc[nn]:=lc1[k];lw[nn]:=w1;lp[op(w1)]:=nn; fi;
od;
for k from 1 to nn do lc1[k]:=lc[k] od;
w:=[0$(i-1),1,0$(m-i),0$(j-1),1,0$(n-j)];
nn0:=nn;
for k from 1 to nn0 do
w1:=w+lw[k];
pt:=lp[op(w1)];
if op(0,pt)=`Integer` then
lc[pt]:=lc[pt]+lc1[k] else nn:=nn+1;if nn>lim then ERROR(`too long array`) fi;lc[nn]:=lc1[k];lw[nn]:=w1;lp[op(w1)]:=nn; fi;
od;
od od;
end;

klm:=proc(lam1,lam2,mu1,mu2) local s,i,j,m,n,k1,k2,sm1,sm2,pm1,pm2,nn1,nn2,km1,km2,pt,lim,x;global nn,lc,lw;
#x:=[];
s:=0;m:=nops(lam1);n:=nops(lam2);
lim:=300000;
sm1:=array(1..lim);sm2:=array(1..lim);pm1:=table();pm2:=table();
nn1:=0;nn2:=0;
for i from 1 to nn do
if i mod 1000 =0 then print(i) fi;
k2:=[seq(lw[i][j],j=m+1..m+n)];
pt:=pm2[op(k2)];
if op(0,pt)=`Integer` then
km2:=sm2[pt];
else
nn2:=nn2+1;
km2:=klmc1(lam2,mu2+k2,m);
sm2[nn2]:=km2;
pm2[op(k2)]:=nn2;
fi;
if km2<>0 then
k1:=[seq(lw[i][j],j=1..m)];pt:=pm1[op(k1)];
if op(0,pt)=`Integer` then
km1:=sm1[pt];
else
nn1:=nn1+1;
km1:=klmb(lam1,mu1+k1);
sm1[nn1]:=km1;
pm1[op(k1)]:=nn1;
fi;
#if km1*km2<>0 then x:=[op(x),[lc[i],k1,k2,km1,km2]] fi;
#if subs(q=1,km1)<0 or subs(q=1,km2)<0 then print(k1,k2,km1,km2) fi;
s:=s+expand(lc[i]*q^sum(k2[j],j=1..n)*km1*km2)
fi;
od;
sort(s);
#s,x;
end;

klm0:=proc(lam1,lam2,mu1,mu2) local s,i,j,m,n,k1,k2,sm1,sm2,pm1,pm2,nn1,nn2,km1,km2,pt,lim,x;global nn,lc,lw;
#x:=[];
s:=0;m:=nops(lam1);n:=nops(lam2);
lim:=300000;
sm1:=array(1..lim);sm2:=array(1..lim);pm1:=table();pm2:=table();
nn1:=0;nn2:=0;
for i from 1 to nn do
if i mod 1000 =0 then print(i) fi;
k2:=[seq(lw[i][j],j=m+1..m+n)];
pt:=pm2[op(k2)];
if op(0,pt)=`Integer` then
km2:=sm2[pt];
else
nn2:=nn2+1;
km2:=klman(lam2,mu2+k2);
sm2[nn2]:=km2;
pm2[op(k2)]:=nn2;
fi;
if km2<>0 then
k1:=[seq(lw[i][j],j=1..m)];pt:=pm1[op(k1)];
if op(0,pt)=`Integer` then
km1:=sm1[pt];
else
nn1:=nn1+1;
km1:=klmam(lam1,mu1+k1);
sm1[nn1]:=km1;
pm1[op(k1)]:=nn1;
fi;
#if km1*km2<>0 then x:=[op(x),[lc[i],k1,k2,km1,km2]] fi;
#if subs(q=1,km1)<0 or subs(q=1,km2)<0 then print(k1,k2,km1,km2) fi;
s:=s+expand(lc[i]*q^sum(k2[j],j=1..n)*km1*km2)
fi;
od;
sort(s);
#s,x;
end;

klm2:=proc(lam1,lam2,mu1,mu2) local s,i,j,m,n,k1,k2,sm1,sm2,pm1,pm2,nn1,nn2,km1,km2,pt,lim;global nn,lc,lw;
s:=0;m:=nops(lam1);n:=nops(lam2);
lim:=300000;
sm1:=array(1..lim);sm2:=array(1..lim);pm1:=table();pm2:=table();
nn1:=0;nn2:=0;
for i from 1 to nn do
if i mod 1000 =0 then print(i) fi;
k2:=[seq(lw[i][j],j=m+1..m+n)];
pt:=pm2[op(k2)];
if op(0,pt)=`Integer` then
km2:=sm2[pt];
else
nn2:=nn2+1;
km2:=klmc1(lam2,mu2+k2,m);
sm2[nn2]:=km2;
pm2[op(k2)]:=nn2;
fi;
if km2<>0 then
k1:=[seq(lw[i][j],j=1..m)];
pt:=pm1[op(k1)];
if op(0,pt)=`Integer` then
km1:=sm1[pt];
else
nn1:=nn1+1;
km1:=klmd(lam1,mu1+k1);
sm1[nn1]:=km1;
pm1[op(k1)]:=nn1;
fi;
s:=s+expand(lc[i]*q^sum(k2[j],j=1..n)*km1*km2)
fi;
od;
sort(s);
end;

dist:=proc(m) local a,i,j,x,na;global nn,lw;
a:=array(1..1000000);na:=0;
for i from 1 to nn do
x:=[seq(lw[i][j],j=1..m)];
j:=1;while (j<=na) and (a[j]<>x) do j:=j+1 od;
if j>na then na:=na+1;a[na]:=x fi
od;
na;
end;

dist1:=proc(m,n) local a,i,j,x,na;global nn,lw;
a:=array(1..1000000);na:=0;
for i from 1 to nn do
x:=[seq(lw[i][j],j=m+1..m+n)];
j:=1;while (j<=na) and (a[j]<>x) do j:=j+1 od;
if j>na then na:=na+1;a[na]:=x fi
od;
na;
end;

klmold:=proc(lam1,lam2,mu1,mu2) local s,i,j,m,n,k1,k2;global nn,lc,lw;
s:=0;m:=nops(lam1);n:=nops(lam2);
for i from 1 to nn do
k1:=[seq(lw[i][j],j=1..m)];k2:=[seq(lw[i][j],j=m+1..m+n)];
s:=s+expand(lc[i]*q^sum(k1[j],j=1..m)*klmb(lam1,mu1+k1)*klmc(lam2,mu2+k2))
od;
s;
end;

verifc:=proc(sw) local i,w,n,wm,lam,t,c,j,m,k,mu,bc;
n:=nops(sw);w:=weights(C.n);
lam:=sum(sw[i]*w[n+1-i],i=1..n);
wm:=weight_mults(lam,C.n);
lam:=[0$n];
for i from 1 to n do
for j from 1 to i do lam[j]:=lam[j]+sw[i] od
od;
for i from 1 to nops(wm) do
t:=op(i,wm);
if op(0,t)=`M` then c:=1 else c:=op(1,t);t:=op(2,t) fi;
mu:=[0$n];
for j from 1 to n do
m:=op(n+1-j,t);
for k from 1 to j do mu[k]:=mu[k]+m od;
od;
m:=subs(q=1,klmc(lam,mu));if m=c then print(`***`,c*t) else ERROR(`---`,c*t,m) fi;
od;
end;

verifa:=proc(sw) local i,w,n,wm,lam,t,c,j,m,k,mu,bc;
n:=nops(sw)+1;w:=weights(cat(A,n-1));
lam:=sum(sw[i]*w[n-i],i=1..n-1);
wm:=weight_mults(lam,cat(A,n-1));
lam:=[0$n];
for i from 1 to n-1 do
for j from 1 to i do lam[j]:=lam[j]+sw[i] od
od;
for i from 1 to nops(wm) do
t:=op(i,wm);
if op(0,t)=`M` then c:=1 else c:=op(1,t);t:=op(2,t) fi;
mu:=[0$n];
for j from 1 to n-1 do
m:=op(n-j,t);
for k from 1 to j do mu[k]:=mu[k]+m od;
od;
m:=subs(q=1,klmam(lam,mu));if m=c then print(`***`,c*t) else ERROR(`---`,c*t,m) fi;
od;
end;

verifb:=proc(sw) local i,w,n,wm,lam,t,c,j,m,k,mu,bc;
n:=nops(sw);w:=weights(B.n);
lam:=sum(sw[i]*w[n+1-i],i=1..n);
wm:=weight_mults(lam,B.n);
if sw[n] mod 2 <>0 then ERROR(`spin column`) fi;
lam:=[0$n];
for i from 1 to n do
for j from 1 to i do if i=n then lam[j]:=lam[j]+sw[i]/2 else lam[j]:=lam[j]+sw[i] fi od
od;
for i from 1 to nops(wm) do
t:=op(i,wm);
if op(0,t)=`M` then c:=1 else c:=op(1,t);t:=op(2,t) fi;
mu:=[0$n];
for j from 1 to n do
m:=op(n+1-j,t);
for k from 1 to j do if j=n then mu[k]:=mu[k]+m/2 else mu[k]:=mu[k]+m fi od;
od;
m:=subs(q=1,klmb(lam,mu));if m=c then print(`***`,c*t) else ERROR(`---`,c*t,m) fi;
od;
end;

verifd:=proc(sw) local i,w,n,wm,lam,t,c,j,m,k,mu,bc;
n:=nops(sw);w:=weights(D.n);
lam:=sum(sw[i]*w[n+1-i],i=1..n);
wm:=weight_mults(lam,D.n);
lam:=[0$n];
for i from 1 to n-2 do
for j from 1 to i do lam[j]:=lam[j]+sw[i] od
od;
for j from 1 to n-1 do lam[j]:=lam[j]+sw[n-1]/2 od;lam[n]:=lam[n]-sw[n-1]/2;
for j from 1 to n do lam[j]:=lam[j]+sw[n]/2 od;
for i from 1 to nops(wm) do
t:=op(i,wm);
if op(0,t)=`M` then c:=1 else c:=op(1,t);t:=op(2,t) fi;
mu:=[0$n];
for j from 1 to n-2 do
m:=op(n+1-j,t);
for k from 1 to j do mu[k]:=mu[k]+m od;
od;
for k from 1 to n-1 do mu[k]:=mu[k]+op(2,t)/2 od;mu[n]:=mu[n]-op(2,t)/2;
for k from 1 to n do mu[k]:=mu[k]+op(1,t)/2 od;
m:=subs(q=1,klmd(lam,mu));if m=c then print(`***`,c*t) else ERROR(`---`,c*t,m) fi;
od;
end;

isdom:=proc(d,m,n) local i,j;
i:=1;
while (i<=m-1) and (d[i]>=d[i+1]) do i:=i+1 od;
if (i=m) and (d[m]>=0) then
j:=1;
while (j<=n-1) and (d[m+j]>=d[m+j+1]) do j:=j+1 od;
if (j=n) and (d[m+n]>=0) then true else false fi
else false
fi;
end;

alldom:=proc(lam1,lam2) global lw,bb,bc,nn,nd,ld;local i,j,k,l1,l2,lr1,lr2,r1,r2,pd,m,n,d,rr;
pd:=table();
m:=nops(lam1);n:=nops(lam2);r1:=rhob(m);r2:=rhoc(n);lr1:=lam1+r1+[(0)$m];lr2:=lam2+r2+[(-m-1/2)$n];ld:=[];nd:=0;rr:=rhobc(m,n);
for j from 1 to nops(bb) do
l1:=actb(bb[j],lr1);
for k from 1 to nops(bc) do
l2:=actb(bc[k],lr2);
for i from 1 to nn do
d:=[op(l1),op(l2)]-rr-lw[i];
if isdom(d,m,n) then
if op(0,pd[op(d)])<>`Integer` then nd:=nd+1;ld:=[op(ld),d];pd[op(d)]:=nd;
fi;
fi;
od;
od;
od;
print(nd);
end;

mult:=proc(lam1,lam2,g) global bb,bc,lp,lc,sbb,sbc;local i,j,mm,l1,l2,m,n,r1,r2,w,lr1,lr2,pd,rr;
mm:=0;m:=nops(lam1);n:=nops(lam2);r1:=rhob(m);r2:=rhoc(n);lr1:=lam1+r1+[(0)$m];lr2:=lam2+r2+[(-m-1/2)$n];rr:=rhobc(m,n);
for i from 1 to nops(bb) do
l1:=actb(bb[i],lr1);
for j from 1 to nops(bc) do
l2:=actb(bc[j],lr2);
w:=[op(l1),op(l2)]-rr-g;
if op(0,lp[op(w)])=`Integer` then pd:=lp[op(w)];mm:=mm+sbb[i]*sbc[j]*lc[pd];
#print(i,j,g,lw[pd],sbb[i]*sbc[j]*lc[pd]);
fi;
od;
od;
mm;
end;

klmp:=proc(lam1,lam2,mu1,mu2) global ld,nd;local i,j,m,n,s;
m:=nops(lam1);n:=nops(lam2);
s:=0;
for i from 1 to nd do
if i mod 100=0 then print(i) fi;
#print(ld[i],mult(lam1,lam2,ld[i]),klmb([seq(ld[i][j],j=1..m)],mu1),klmc([seq(ld[i][j],j=m+1..m+n)],mu2));
s:=s+expand(mult(lam1,lam2,ld[i])*klmb([seq(ld[i][j],j=1..m)],mu1)*klmc([seq(ld[i][j],j=m+1..m+n)],mu2));
od;
sort(s);
end;

The package Alcove_Path for crystals

This Maple package, based on the alcove path model and using Stembridge's Coxeter and Weyl packages, is useful for investigating the combinatorics of irreducible crystals corresponding to the complex simple Lie algebras. C. L. was partially supported by National Science Foundation grants DMS-0403029.

Documentation

"Alcove_Path" : A package for the alcove path model based on Stembridge's Maple packages Coxeter/Weyl

This package, based on the alcove path model, is useful for investigating the combinatorics of irreducible crystals corresponding to the complex simple Lie algebras.

The alcove path model was introduced in my joint papers with A. Postnikov entitled "Affine Weyl groups in K-theory and representation theory" and "A combinatorial model for crystals of Kac-Moody algebras"; both are available on my webpage and on arXiv. I refer the user to the short survey paper "A new combinatorial model in representation theory" on my webpage for the concepts that are mentioned below. In the alcove path model, vertices of crystals are labeled by certain finite sets of positive integers, called admissible subsets (in type A, there is an easy weight-preserving bijection between admissible subsets and semistandard Young tableaux, described in a forthcoming publication); in the discussion below, we will identify a vertex of a crystal with the corresponding admissible subset.

The alcove path model can be thought of as a discrete counterpart to the Littelmann path model. The connection to the latter is discussed in the second joint paper mentioned above. The advantage of having this discrete model is that more combinatorics can be developed. For instance, I gave a combinatorial realization of Lusztig's involution on crystals, in finite types, in the paper "On the combinatorics of crystal graphs, I. Lusztig's involution" (also available on my webpage). This paper also contains a type-independent analog of jeu de taquin for tableaux; this is not yet implemented, but has several applications. A forthcoming paper entitled "On the combinatorics of
crystal graphs II. Crystal multiplication" will contain a combinatorial realization of the commutator in the category of crystals (for finite types) that has been recently studied by Henriques and Kamnitzer. This package was first used for experiments related to this paper. Projects include the connection between the alcove path model and the tableaux models in classical types, plus the Kyoto path model in the affine types.

Stembridge's packages, on which this package is based, have to be loaded into Maple first (available from www.math.lsa.umich.edu/~jrs), then the "alcove_path" file is loaded. The data format in Stembridge's packages is used. The procedure "prepare" (see below) has to be run each time before switching to a new root system.

Procedures

prepare - initializes the basic structures related to a root system

lchain - constructs a particular lambda-chain of roots

fold - constructs all admissible subsets for a given lambda-chain

weight - computes the weight of an admissible subset

admfold - constructs the admissible folding corresponding to an admissible subset

tstlrfold - tests whether an admissible subset indexes an irreducible crystal in the tensor product of two irreducible crystals

fp - constructs the result of applying a root operator f_p to an admissible subset

repf - applies a sequence of root operators f_p to an admissible subset

ep - constructs the result of applying a root operator e_p to an admissible subset

repe - applies a sequence of root operators e_p to an admissible subset

fp2 - constructs the result of applying a root operator f_p to a pair of admissible subsets (in a tensor product of irreducible crystals)

ep2 - constructs the result of applying a root operator e_p to a pair of admissible subsets (in a tensor product of irreducible crystals)

gototop1 - constructs a sequence of root operators e_p corresponding to a path from a vertex to the highest weight vertex in an irreducible crystal

gototop - constructs a sequence of root operators e_p corresponding to a path from a vertex to the highest weight vertex in a tensor product of two irreducible crystals

gotobot1 - constructs a sequence of root operators f_p corresponding to a path from a vertex to the lowest weight vertex in an irreducible crystal

gotobot - constructs a sequence of root operators f_p corresponding to a path from a vertex to the lowest weight vertex in a tensor product of two irreducible crystals

compkey - computes the Weyl group element w(J) for an admissible subset J (this is an analog of the Lascoux-Schuetzenberger key, and is relevant to the corresponding Demazure modules)

sp - constructs the result of applying a simple reflection to a vertex of an irreducible crystal

actw - constructs the result of applying a sequence of simple reflections to a vertex of an irreducible crystal

sp2 - constructs the result of applying a simple reflection to a pair of admissible subsets (in a tensor product of two irreducible crystals)

actw2 - constructs the result of applying a sequence of simple reflections to a pair of admissible subsets (in a tensor product of two irreducible crystals)

invol - constructs the result of applying Lusztig's involution to a vertex of an irreducible crystal

commutator - realizes the commutator in the category of crystals studied by Henriques and Kamnitzer

Description

PROCEDURE: prepare - initializes the basic structures related to a root system

CALLING SEQUENCE: prepare(rs)

PARAMETERS: rs = root system type (cf. Stembridge's packages)

OUTPUT: no output, but some invisible global parameters are initialized (b = simple roots, w = fundamental weights etc).

PROCEDURE: lchain - constructs a particular lambda-chain (of roots)

CALLING SEQUENCE: lchain(lam)

PARAMETERS: lam = dominant weight indexing an irreducible representation (cf. Stembridge's format)

OUTPUT: lambda-chain of roots; sequence of "levels" corresponding to it

PROCEDURE: fold - constructs all admissible subsets for a given lambda-chain

CALLING SEQUENCE: fold(lc)

PARAMETERS: lc = lambda-chain constructed by "lchain"

OUTPUT: sequence of admissible subsets (each one is a sequence J_i=[p, q, ...]); corresponding sequence of Weyl group elements [w(J_1), w(J_2), ...], each one being in the form of a reduced word - see "compkey"

PROCEDURE: weight - computes the weight of an admissible subset

CALLING SEQUENCE: weight(lc,lev,fld,lam)

PARAMETERS: lc = lambda-chain of roots

lev = sequence of "levels" associated to a lambda-chain (output of lchain)

fld = admissible subset

lam = corresponding dominant weight

OUTPUT: weight of fld, denoted \mu(fld), in the form [a_1,...,a_n], where a_i is the coefficient of the ith fundamental weight in \mu(fld)

PROCEDURE: admfold - constructs the admissible folding corresponding to an admissible subset

CALLING SEQUENCE: admfold(lc,fld)

PARAMETERS: lc, fld = same as the ones with the same name in "weight"

OUTPUT: sequence of pairs of roots (see definition of "admissible folding" in the paper cited at the beginning of the file)

PROCEDURE: tstlrfold - tests whether an admissible subset fld indexes an irreducible crystal in the tensor product of two irreducible crystals

CALLING SEQUENCE: tstlrfold(nu,fld,key,lc)

PARAMETERS: lc, fld = same as the ones with the same name in "weight" (correspond to the second irreducible crystal in the tensor product)

key = Weyl group element w(fld) in the form of a reduced word, see "compkey"; this is computed by "compkey" and is also an output of "fold"

nu = dominant weight corresponding to the first irreducible crystal in the tensor product

NOTE: a highest weight vertex of an irreducible crystal in a tensor product of two irreducible crystals corresponds to a pair of admissible subsets of the form [[], fld]

OUTPUT: true/false

PROCEDURE: fp - constructs the result of applying a root operator f_p to an admissible subset

CALLING SEQUENCE: fp(p,lc,fld)

PARAMETERS: lc, fld = same as the ones with the same name in "weight"

p = index of root operator "f"

OUTPUT: f_p(fld) or 0, if the former is not defined

PROCEDURE: repf - applies a sequence of root operators f_p to an admissible subset

CALLING SEQUENCE: repf(lc,fld,sf)

PARAMETERS: lc, fld = same as the ones with the same name in "weight"

sf= sequence of indices of root operators "f" (applied from left to right)

OUTPUT: admissible subset obtained by applying the given sequence of root operators, or 0 if some root operator is undefined

PROCEDURE: ep - constructs the result of applying a root operator e_p to an admissible subset

CALLING SEQUENCE: ep(p,lc,fld,w)

PARAMETERS: p, lc, fld = same as the ones with the same name in "fp"

w = weight of fld

OUTPUT: e_p(fld) or 0, if the former is undefined; weight of e_p(fld) or 0

PROCEDURE: repe - applies a sequence of root operators e_p to an admissible subset

CALLING SEQUENCE: repe(lc,fld,sf,w)

PARAMETERS: lc, fld, sf = same as in "repf"

w = weight of fld

OUTPUT: same as in "repf"

PROCEDURE: fp2 - constructs the result of applying a root operator f_p to a pair of admissible subsets (in a tensor product of irreducible crystals)

CALLING SEQUENCE: fp2(p,lc1,fld1,lc2,fld2,w1,w2)

PARAMETERS: p = index of root operator "f"

lc1, lc2 = the two lambda-chains corresponding to the two irreducible crystals in the tensor product

fld1,fld2 = admissible subsets corresponding to lc1 and lc2

w1,w2 = weights of fld1 and fld2

OUTPUT: f_p([fld1,fld2]) or 0, if the former is undefined; weights of the two admissible subsets in f_p([fld1,fld2]) or 0, 0

PROCEDURE: ep2 - constructs the result of applying a root operator e_p to a pair of admissible subsets (in a tensor product of irreducible crystals)

CALLING SEQUENCE: ep2(p,lc1,fld1,lc2,fld2,w1,w2)

PARAMETERS: same as in "fp2"

OUTPUT: e_p([fld1,fld2]) or 0, if the former is undefined; weights of the two admissible subsets in e_p([fld1,fld2]) or 0, 0

PROCEDURE: gototop1 - constructs a sequence of root operators e_p corresponding to a path from a vertex to the highest weight vertex in an irreducible crystal

CALLING SEQUENCE: gototop1(lc,fld,w)

PARAMETERS: lc, fld = same as the ones with the same name in "weight"

w = weight of fld

OUTPUT: sequence of indices of root operators "e_p" applied from left to right

PROCEDURE: gototop - constructs a sequence of root operators e_p corresponding to a path from a vertex to the highest weight vertex in a tensor product of two irreducible crystals (see NOTE in "tstlrfold")

CALLING SEQUENCE: gototop(lc1,lc2,lev1,lev2,fld,lam1,lam2)

PARAMETERS: lc1, lc2 = the two lambda-chains corresponding to the two irreducible crystals in the tensor product

lev1, lev2 = the two sequences of "levels" corresponding to lc1 and lc2 (output of "lchain")

fld = pair of admissible subsets

lam1, lam2 = dominant weights corresponding to the two irreducible crystals

OUTPUT: sequence of indices of root operators "e_p" applied from left to right; highest weight pair in the tensor product

PROCEDURE: gotobot1 - constructs a sequence of root operators f_p corresponding to a path from a vertex to the lowest weight vertex in an irreducible crystal

CALLING SEQUENCE: gotobot1(lc,fld)

PARAMETERS: same as in "gototop1"

OUTPUT: sequence of indices of root operators "f_p" applied from left to right

PROCEDURE: gotobot - constructs a sequence of root operators f_p corresponding to a path from a vertex to the lowest weight vertex in a tensor product of two irreducible crystals

CALLING SEQUENCE: gotobot(lc1,lc2,lev1,lev2,fld,lam1,lam2)

PARAMETERS: same as in gototop

OUTPUT: sequence of indices of root operators "f_p" applied from left to right; lowest weight pair in the tensor product

PROCEDURE: compkey - computes the Weyl group element w(J) for an admissible subset J (this is an analog of the Lascoux-Schuetzenberger key, and is relevant to the corresponding Demazure module)

CALLING SEQUENCE: compkey(lc,fld)

PARAMETERS: same as the ones with the same name in "weight"

OUTPUT: the Weyl group element w(fld) in the form of a reduced word

PROCEDURE: sp - constructs the result of applying a simple reflection to a vertex of an irreducible crystal

CALLING SEQUENCE: sp(p,lc,fld,w)

PARAMETERS: lc, fld = same as the ones with the same name in "weight"

p = index of the corresponding simple reflection

w = weight of fld

OUTPUT: s_p(fld)

PROCEDURE: actw - constructs the result of applying a sequence of simple reflections to a vertex of an irreducible crystal

CALLING SEQUENCE: actw(r,lc,fld,w)

PARAMETERS: lc, fld, w = same as the ones with the same name in "sp"

r = sequence of indices of simple reflections (applied from right to
left)

OUTPUT: admissible subset obtained by applying the simple reflections
in r to fld

PROCEDURE: sp2 - constructs the result of applying a simple reflection to a pair of admissible subsets (in a tensor product of two irreducible crystals)

CALLING SEQUENCE: sp2(p,lc1,fld1,lc2,fld2,w1,w2)

PARAMETERS: p = index of simple reflection applied

lc1, lc2, fld1, fld2, w1, w2 = same as the ones with the same name in "fp2"

OUTPUT: s_p([fld1,fld2]); weights of the two admissible subsets in s_p([fld1,fld2])

PROCEDURE: actw2 - constructs the result of applying a sequence of simple reflections to a pair of admissible subsets (in a tensor product of two irreducible crystals)

CALLING SEQUENCE: actw2(r,lc1,fld1,lc2,fld2,w1,w2)

PARAMETERS: lc1, lc2, fld1, fld2, w1, w2 = same as the ones with the same name in "fp2"

r = same as in "actw"

OUTPUT: pair of admissible subsets obtained by applying the simple reflections in r to [fld1,fld2]

PROCEDURE: invol - constructs the result of applying Lusztig's involution to a vertex of an irreducible crystal

CALLING SEQUENCE: invol(lc,fld)

PARAMETERS: same as the ones with the same name in "weight"

OUTPUT: admissible subset obtained by applying Lusztig's involution to fld

PROCEDURE: commutator - realizes the commutator in the category of crystals studied by Henriques and Kamnitzer

CALLING SEQUENCE: commutator(lc1,lc2,lev2,fld2,p2,lam2)

PARAMETERS: lc1, lc2 = lambda-chains for two irreducible crystals A and B

lev2 = the level sequence corresponding to the crystal B (output of "lchain")

lam2 = the dominant weight corresponding to the crystal B

fld2 = an admissible subset for lc2, such that [[],fld2] is a highest weight vertex for an irreducible crystal in the tensor product A x B

p2 = w(fld2), see "compkey"; this Weyl group element is computed by compkey and is also an output of "fold"

OUTPUT: the pair [[],fld1] in the tensor product B x A that corresponds to [[],fld2] in A x B via the commutator studied by Henriques-Kamnitzer. Here fld1 is an admissible subset for the lambda-chain lc1.

2006 C. Lenart

The package

#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;

References related to clustering software

C. Lenart, Software for classification, Studia Univ. "Babes-Bolyai", Mathematica 36 no. 3 (1991), 41-49.
C. Lenart, Topographical data management system, Studia Univ. "Babes-Bolyai", Mathematica 36 no. 3 (1991), 51-59.

Abstract related to clustering software

Between 1989 and 1992, I developed two computer packages, consisting of several programs written in Pascal and C languages. The first one is a package for clustering, while the second one is a package for processing geological data (including map processing).