interface(quiet=true):

nops(Disc);

% Initializes Sn = S4
%
%
%

S3:=[
[1,2,3],[1,3,2],[2,1,3],[2,3,1],[3,1,2],[3,2,1]]:

vars:=[s,t,u]:

Poly := proc(Coef,EXP)  
 local pi,poly:
 poly:=0:
 if EXP=[10,10,10] then 
  Coef*s^10*t^10*u^10
 else 
  for pi in S3 do
   poly:= poly + s^EXP[pi[1]] * t^EXP[pi[2]] * u^EXP[pi[3]];
  od:
  Coef*poly*nops(poly)/6
 fi:
end:

PART:=proc(exs)
#
#    This takes an exponent vector and sorts it into a partition
#
local II,sorted,try,j,i,big,jj;

 II:={1,2,3}:	sorted:=[0,0,0]:	try:=[]:
 for j in II do  try:= [try[],exs[j]]: od:
 big:=max(try[]):	 sorted[1]:=big:
 for jj from 1 to 3 do 
   if (member(jj,II) and big=exs[jj]) then i:=jj: fi:   od:
 II:= II minus {i}:	try:=[]:
 for j in II do  try:= [try[],exs[j]]: od:
 big:=max(try[]):	sorted[2]:=big:
 for jj from 1 to 3 do 
    if (member(jj,II) and big=exs[jj]) then i:=jj: fi:	od:
 II:= II minus {i}:	try:=[]:
 for j in II do  try:= [try[],exs[j]]: od:
 big:=max(try[]):	sorted[3]:=big:
 for jj from 1 to 3 do 
   if (member(jj,II) and big=exs[jj]) then i:=jj: fi:	od:
 II:= II minus {i}:
 sorted[4]:=exs[op(1,II)]:
 sorted
end:

deconstruct := proc(Term)
  local coef,jj,T,exlist:
 exlist:=[0,0,0]:
 coef:=op(1,Term):
  for jj from 2 to nops(Term) do
   T:= op(jj,Term):
   if (nops(T)=1) then exlist[jj-1]:=1:
      else exlist[jj-1]:= op(2,T):
   fi;
  od:
  if (nops(Term)=3) then exlist:=[exlist[],0]:  fi:

[coef,exlist]

end:

D1:= Disc:


TERMS:={}:

while not(D1=0) do

T:=deconstruct(op(1,D1));
TERMS:=TERMS union {T}:

D1:= D1-Poly(T[1],T[2]):
od:

terms:={}:
for j from 1 to nops(TERMS) do
terms:= terms union {[op(j,TERMS)[1],PART(op(j,TERMS)[2])]}:
od:

lprint(terms);

quit;