Magma includes functions for working with symmetric polynomials.

ElementarySymmetricPolynomial(P, k) : RngMPol, RngIntElt -> RngMPolElt

Given a polynomial ring P of rank n, and an integer k with 1 <= k <= n, return the k-th elementary symmetric polynomial of P.

IsSymmetric(f, S) : RngMPolElt, RngMPol -> BoolElt, RngMPolElt

Given a polynomial f from a polynomial ring P of rank n, return whether f is a symmetric polynomial of P (i.e., is symmetric in all the n variables of P). If the answer is true, a polynomial g from a new polynomial ring of rank n is returned such that f = g(e_1, ..., e_n), where e_i is the i-th elementary symmetric polynomial of P. If g is desired to be a member of a particular polynomial ring S of rank n (to obtain predetermined names of variables, for example), then S may also be passed.

> P<a, b, c, d> := PolynomialRing(RationalField(), 4, "grevlex");
> f :=
> a^2*b^2*c*d + a^2*b*c^2*d + a*b^2*c^2*d + a^2*b*c*d^2 + a*b^2*c*d^2 + 
>     a*b*c^2*d^2 - a^2*b^2*c - a^2*b*c^2 - a*b^2*c^2 - a^2*b^2*d -
>     3*a^2*b*c*d - 3*a*b^2*c*d - a^2*c^2*d - 3*a*b*c^2*d - b^2*c^2*d -
>     a^2*b*d^2 - a*b^2*d^2 - a^2*c*d^2 - 3*a*b*c*d^2 - b^2*c*d^2 -
>     a*c^2*d^2 - b*c^2*d^2 + a + b + c + d;
> // Check orbit under Sym(4) has size one:
> #(f^Sym(4));
1
> Q<e1, e2, e3, e4> := PolynomialRing(RationalField(), 4);
> l, E := IsSymmetric(f, Q); 
> l;
true
> E;
e1 - e2*e3 + e2*e4