calc_l_weight:= proc(term)
   local i1,i2,i3,i4,i5,i6,i0;
   i0 := degree(term,lambda[0]);
   i1 := degree(term,lambda[1]);
   i2 := degree(term,lambda[2]);
   i3 := degree(term,lambda[3]);
   i4 := degree(term,lambda[4]);
   return(i0*15 + i1*12 + i2*9 + i3*6 + i4*3);
end:
chop_order2 := (Y,Z)-> select(X->calc_l_weight(X)<Z, expand(Y)):
#
pl_ult:=proc(lev, size)
local All,Good,Bad, i, P,Blocks,len, Perms:
All:=seq(i, i=1..lev-1); Good:=seq(i*n, i=1..(lev-1)/n); Bad:={All} minus {Good};
P:=remove(X -> nops(X)>size, remove(has, combinat[partition](lev-1), Bad)): 
Blocks:=[]:
for i from 1 to nops(P) do len:=nops(op(i,P)): Blocks:=[op(Blocks), [seq(0,j=1..size-len),op(op(i,P))] ]: od:
Perms:=[]:
for i from 1 to nops(Blocks) do Perms:=[op(Perms), op(combinat[permute](op(i,Blocks),size)) ]: od:
RETURN(Perms): 
end proc:
#
Expand_Product:=proc(Term, lev, size, subs_set, havezeta)
local i,j,k, div,term, T,count, S,wt, Result,perms,W,test, deg;
if(type(Term, `+`)) then print(error_not_product); RETURN(error_not_product); fi:
div:=remove(has, Term, si): # print(div); 
term:=Term/div: count:=1: # print(term);   

if   type(term, `^`) then  
          if(op(2,term)<>size) then print(error_weird_power); RETURN(error_weird_power):  
          else for i from 1 to size do T[i]:=expand(eval(subs(si=sigma,op(1,term)))): count:=count+1: od:
          fi:
elif type(term, `*`) then for i from 1 to nops(term) do
  if type(op(i,term),`^`) then  
    for j from count to count+op(2,op(i,term))-1 do ; 
      T[j]:=expand(eval(subs(si=sigma,op(1,op(i,term))))): count:=count+1: 
    od:
  elif op(0,op(i,term))=si or op(0,op(i,term))=diff then
    T[count]:=expand(eval(subs(si=sigma,op(i,term)))): count:=count+1:    
  else print(error_weird_term_part): RETURN(error_weird_term_part):
  fi:
od: fi:
if count<>size+1 then print(error_count): RETURN(error_count): fi: 
T[1]:=expand( chop_order2( div*T[1], lev) ): 
for i from 1 to size do if( nops(T[i])=1 ) then if(T[i]=0) then RETURN(0):  fi: fi: od:

if havezeta=yes then for i from 1 to size do 
  deg:=degree(T[i],zeta): # print(degree=deg); 
  T[i]:= expand(subs(zeta^2=-zeta-1, subs(
         [seq( zeta^i=1, i=3..3*ceil(deg/3), 3), seq( zeta^i=zeta, i=4..4+3*ceil(deg/3), 3), seq( zeta^i=zeta^2, i=5..5+3*ceil(deg/3), 3) ]
         , T[i] ))): # print(degree=degree(T[i],zeta)); 
od: fi:

for i from 1 to size do # print(nops(T[i])); 
  T[i]:=select(X->subs(subs_set,X)<>0, T[i]): # print(nops(T[i])); 
od:

for i from 1 to size do for j from 0 to lev-1 do S[i,j]:=0: od: od:
for i from 1 to size do # print(S||i); 
  if type(T[i],`+`) then for j from 1 to nops(T[i]) do # if modp(j,1000)=0 then print(j); fi;
                           wt:=calc_l_weight( op(j,T[i]) ):
                           if(wt<lev) then S[i,wt]:=S[i,wt] + op(j,T[i]): fi:
                         od:
  else wt:=calc_l_weight( T[i] ): if(wt<lev) then S[i,wt]:=S[i,wt] + T[i]: fi: 
  fi: 
od: 
for i from 1 to size do test[i]:=add( S[i,j], j=0..lev-1); if T[i]-test[i]<>0 then print(error5); print([i,T[i]]); fi; od:    

Result:=0: wt:='wt': 
for wt from (lev) by -3 to 4 do # print(wt);   
  perms:=pl_ult(wt, size): # print(perms);  
  for i from 1 to nops(perms) do # print(op(i,perms));
    for j from 1 to size do W[j]:=op(j,op(i,perms)): od:  
    Result:=Result + subs(subs_set, mul( S[i,W[i]], i=1..size ) ): 
  od:
od:    
Result:=Result + subs(subs_set,mul(S[i,0], i=1..size)):

RETURN(expand(Result)): 
end proc: