#  Initializes the symmetric group on 4 letters
S4:=[
[1,2,3,4],[1,2,4,3],[1,3,2,4],[1,3,4,2],[1,4,2,3],[1,4,3,2],
[2,1,3,4],[2,1,4,3],[2,3,1,4],[2,3,4,1],[2,4,1,3],[2,4,3,1],
[3,1,2,4],[3,1,4,2],[3,2,1,4],[3,2,4,1],[3,4,1,2],[3,4,2,1],
[4,1,2,3],[4,1,3,2],[4,2,1,3],[4,2,3,1],[4,3,1,2],[4,3,2,1]]:

SYMMETRIZE:=proc(PPOLY)
# This procedure takes a polynomial and symmetrizes it, 
# returning a normalized version.
#
local nterms,i,Term,EXP,COEFF,jj,T,SYMPOLY,Poly,pi,POLY;

POLY:=2*expand(PPOLY):

nterms:=nops(POLY):
SYMPOLY:=0:
for i from 1 to nterms do
EXP:=[0,0,0,0]:
Term:=op(i,POLY):			#
COEFF:=op(1,Term):			#   This takes 
  for jj from 2 to nops(Term) do	#   the ith term 
   T:= op(jj,Term):			#   of POLY
   if (nops(T)=1) then EXP[jj-1]:=1:	#   and computes its
      else EXP[jj-1]:= op(2,T):		#   exponent vector
   fi:					#
  od:					#
  if (nops(Term)=3) then		#
   EXP:=[EXP[],0]:  fi:			#
 Poly:=0:
  for pi in S4 do
   Poly:= Poly + s^EXP[pi[1]] * t^EXP[pi[2]] * u^EXP[pi[3]] * v^EXP[pi[4]];
  od:
SYMPOLY:=SYMPOLY+COEFF*Poly:
od:
					# Normalize the symmetrized
SYMPOLY:=SYMPOLY/igcd(coeffs(SYMPOLY))	# polynomial
end: