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

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;

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;