Additionally to what malex said, here is proc which also works
as a certificate whether the inout is symmetric.
The problem of expressing a symmetric polynomial in terms of
the elementary polynomials is just rewriting process.
See for instance
Cox, Little, O'Shea: Ideals, varieties and algorithms,
Ch. 7, especially p. 314
there you will also find the notice that
Gauss already has answered this question.
Now with the help of Groebner basis calculations, this rewriting
process is automatized by computing the normalform.
First we need a procedure to generate the elementary symmetric polynomials.
Code:
proc elemsymp(list #)
"USAGE: elemsymp(); or elemsymp(k); k int, 0<=k<=nvars(basering)
RETURN: ideal I; containing the elementary symmetric polynomials
I[k] = s_k, 1<=k<=nvars, deg(s_k) = k
poly, the elementary symmetric polynomial s_k, 0<=k<=nvars(basering)
EXAMPLE: example elemsymp; shows an example
"
{
int n = nvars(basering);
int i,k;
if (size(#))
{
k=#[1];
if (k<0 or k>n) {ERROR("// -- proc elemsymp: argument outside range");}
if (k==0) {return(1);}
}
def d = basering;
list rl = ringlist(d);
rl[2]= insert(rl[2],"@Z@",n);
def rnew = ring(rl);
setring rnew;
poly f=1;
for (i=1;i<=n;i++)
{
f = f*(var(n+1)-var(i));
}
ideal I = coeffs(f,var(n+1));
for (i=1;i<=n;i++)
{
I[i] = (-1)^(n-i+1)*I[i];
}
I = I[n..1];
setring d;
ideal I = fetch(rnew,I);
if (k) {return(I[k]);}
else { return(I);}
}
When called without any argument, this proc returns an ideal such
that slot i (1<=i<=nvars(basering)) contains the elementary polynomial of degree i. (So the constant polynomial is ommited.)
If this proc is called with an integer i , (0<=i<=nvars), then it returns the i-th symmetric polynomial.
Code:
> ring r=0,(x,y,z),dp;
> elemsymp();
_[1]=x+y+z
_[2]=xy+xz+yz
_[3]=xyz
> elemsymp(0);
1
> elemsymp(1);
x+y+z
> elemsymp(2);
xy+xz+yz
> elemsymp(3);
xyz
Now the rewriting process:
Proposition 4 p. 314 says:
[list=] Let e_1,...,e_n the elementary symmetric polynomials in the variables x_1,...,x_n.
Introduce new variables y_1,...,y_n, such the monomial in the old
variables x_1,...,x_n are greater h´than those in y_i
Form the ideal I = <e_1 - y_1,..., e_n - y_n>,
Compute a Grobener basis G of I,
and determine for the given polynomial f the normalform w.r.t. G.
Then f is symmetric if and only if g only depends on the variables y_1,...,y_n.
Moreover, in this case, g is the desired polynomial.
This is, substituting the variable y_i by the poylnomial e_i
gives f = g(e_1,...,g_n).
[/list=]
Here is the proc which does this, it uses the proc
substitute from
poly.lib.
Code:
proc reprsympoly(poly f)
"USAGE: reprsympoly(f); f poly
RETURN: poly 0, if f is not symmetric, otherwise
the polynomial which expresses f in the elementary polynomials,
ASSUME: global polynomial ordering of the basering
NOTE: the polynomial is considered symmetric w.r.t. to the variables
of the basering
EXAMPLE: example reprsympoly(f); shows an example
"
{
int n = nvars(basering);
def d = basering;
ideal Isym = elemsymp();
ideal maxid1 = maxideal(1);
list rl = ringlist(d);
for (int i=1;i<=n;i++)
{
rl[2]= insert(rl[2],"@"+string(var(i)),n+i-1);
}
def rnew = ring(rl);
setring rnew;
ideal I = fetch(d,Isym);
for (i=1;i<=n;i++)
{
I[i] = -var(n+i) +I[i];
}
ideal Istd = std(I);
poly f = fetch(d,f);
poly G = NF(f,Istd);
ideal maxid1 = maxideal(1);
ideal I1 = maxid1[1..n];
ideal I2 = maxid1[n+1..2*n];
ideal varG = variables(G);
if (substitute(G,I1,ideal(0:n))!=G) // poly.lib
{
return(0);
}
else
{
G = substitute(G,I2,I1);
}
setring d;
poly G = fetch(rnew,G);
return(G);
}
If the input polynomial is not symmetric, then 0 is return,
For the convenience of the user, the returned polynomial is in the original variables.
Code:
> LIB "poly.lib"; // for subsitute
> ideal I = elemsymp();
> I;
I[1]=x+y+z
I[2]=xy+xz+yz
I[3]=xyz
> poly f = -7*I[1]^3 +4*I[3] - 11 * I[2] + 1; // malex example
> reprsympoly(f); // leading coeff is the right one
-7x3-11y+4z+1
> reprsympoly(xy+z); // not symmetric
0
> reprsympoly(x3y+x3z+xy3+xz3+y3z+yz3);
x2y-2y2-xz
> substitute(_,maxideal(1),elemsymp());
x3y+xy3+x3z+y3z+xz3+yz3
> reprsympoly((x-y)^2*(x-z)^2*(y-z)^2);
x2y2-4x3z-4y3+18xyz-27z2
> factorize(substitute(_,maxideal(1),elemsymp()));
[1]:
_[1]=1
_[2]=x-z
_[3]=-x+y
_[4]=y-z
[2]:
1,2,2,2
Further comments:
1.) As usuals, the certificat "symmetric" is w.r.t. the ambiant ring.
This means, in a ring with three variables the polynomial
x+y is not (considered as) symmetric.
Cleary, it is not invariant under the permutations involving z.
Code:
> string(basering);
(0),(x,y,z),(dp(3),C)
> reprsympoly(x+y);
0
2.) The proc reprsympoly does not need a substitution
y_1 <- e_1,...,y_n <- e_n
to verify the correctness of the computation
You may do this afterward, as shown in the examples above,
but in general the substitution may take longer than computing
the representation itself.
The following timings are taken on a slow machine, to demonstrate the difference.
Code:
> poly F = (((x-y)^2*(x-z)^2*(y-z)^2)^21);
> timer=0;
> timer = 1;
> poly H = reprsympoly(F); // find the representation
//used time: 37.99 sec
> timer;
38
> timer = 1;
> poly S = substitute(H,maxideal(1),elemsymp()); // go to check the result
//used time: 64.68 sec
> timer;
65
> F == S;
1
//-------------------
Ch. Gorzel