#by Viktor Petrov 2006.

with(combinat,multinomial):

initchow:=proc(gid) 
global group_id;
#highestweight; 
#local i; 
group_id:=gid; 
#highestweight:=0:
#for i from 1 to nops(weights(group_id)) do 
#highestweight:=highestweight+weights(group_id)[i];od: 
diagram(gid); 
end proc:

gennext:=proc(p,r) 
local newp, i, m; 
m:=nops(p); newp:=p; i:=m; 
while p[i]=0 do i:=i-1; od; 
newp[i]:=0; 
newp[i-1]:=p[i-1]+1; 
newp[m]:=p[i]-1; 
newp; 
end proc:

gennextcomb:=proc(p,n,k) 
local newp, i, j, m; 
m:=nops(p); 
newp:=p;
i:=k; 
while p[i]=n-m+i do i:=i-1; od;
for j from i to m do newp[j]:=p[i]+j-i+1; od; 
newp; 
end proc:

triang:=proc(p) option remember;
 local m, k, ind, i, j, q, s, res, newp; 
 k:=nops(p);
 if p[k]=0 then return 0; fi;
 if k=1 then return 1; fi; res:=0; ind:=[]; 
 for i to k-1 do if triangA[i,k]<>0 then ind:=[op(ind),i]; fi; od; 
 m:=nops(ind);
 if m=0 then if p[k]=1 then return triang(p[1..k-1]); fi; return 0; fi;
 q:=[seq(0,i=1..m)]; q[m]:=p[k]-1;
 while true do s:=multinomial(p[k]-1,op(q)); 
 for i to m do s:=s*triangA[ind[i],k]^(q[i]); od;
 newp:=p[1..k-1]; for i to m do newp[ind[i]]:=p[ind[i]]+q[i]; od;
 res:=res+s*triang(newp);
 if q[1]=p[k]-1 then break; fi; q:=gennext(q); od; res; 
end proc:

generatelist:=proc(w,u) local p,n,m,i,S,k,D1,D2,v; 
 n:=nops(w); m:=nops(u); S:=[];
 D2:=pos_roots(group_id); 
 D1:={seq(reflect(seq(base(group_id)[u[j]],j=1..nops(u)),-D2[i]),i=1..nops(D2))} intersect {op(D2)}; 
 p:=[seq(i,i=1..m)]; while true do k:=m+1; D2:={};
 for i to m while k=m+1 do 
 v:=reflect(seq(base(group_id)[w[p[j]]],j=1..i-1),base(group_id)[w[p[i]]]); 
 if (v in D2) or not (v in D1) then k:=i; fi; D2:=D2 union {v}; od; 
 if k=m+1 then S:=[op(S),p]; k:=m; fi; if p[1]=n-m+1 then break; fi; 
 p:=gennextcomb(p,n,k); od; S;
end proc:

inittriangA:=proc(w) global triangA; local c,i,j; 
 c:=cartan_matrix(group_id); triangA:=matrix(nops(w),nops(w));
 forget(triang);
 for i to nops(w) do for j to nops(w) do
 if i<j then triangA[i,j]:=-c[w[j],w[i]]; else triangA[i,j]:=0; fi; od; od;
 triangA;
end proc:

chowcoef:=proc(w,u,v) global triangA; local l1,l2,i,j,k,p,res;
 inittriangA(w);
 l1:=generatelist(w,u); l2:=generatelist(w,v); res:=0;
 for i to nops(l1) do for j to nops(l2) do p:=[seq(0,i=1..nops(w))];
 for k to nops(u) do p[l1[i][k]]:=p[l1[i][k]]+1; od;
 for k to nops(v) do p[l2[j][k]]:=p[l2[j][k]]+1; od;
 res:=res+triang(p); od; od; res;
end proc:

steencoef:=proc(w,u,p,k) local l, q, x, res, i, j; 
 if k>nops(u) then return 0; fi;
 inittriangA(w); l:=generatelist(w,u); res:=0; for i to nops(l) do
 q:=[seq(i,i=1..k)];
 while true do x:=[seq(0,i=1..nops(w))];
 for j to nops(u) do x[l[i][j]]:=1; od; for j to k do x[l[i][q[j]]]:=p; od;
 res:=(res+triang(x)) mod p; if q[1]=nops(u)-k+1 then break; fi;
 q:=gennextcomb(q,nops(u),k); od; od; res;
end proc: