


elem := (i,vrs) -> sum(subsets(vrs,i), product);
expandElem := (P,vrs,els) -> (
    if P == 0 then return 0;
    c := coefficients(P,Variables=>vrs);
    M := c#0_(0,0); C := c#1_(0,0);
    e := append((first exponents M)_(apply(vrs,index)),0);
    ee := apply(#vrs, i -> e#i - e#(i+1));
    if any(ee, i->i<0) then error "nonsymmetric polynomial";
    Q := P - C * product(#vrs, i -> (elem(i+1,vrs))^(ee#i));
    sub(C,ring first els) * product(#vrs, i -> (els#i)^(ee#i)) + expandElem(Q,vrs,els)
    )


F3 = product apply(subsets(gens R, 3), set -> sum set)
rep = expandElem(F3,gens R,gens S)


f0 = (x^6 + t*x^5 + t)

f1 = (x^4 + t) * (x^2 + t*x + t)

f2 = (x + 1) * (x^5 + t*x + t)

f3 = (x^2 + t) * (x^2 + t*x + t)* (x^2 + t^2*x + t)

f5 = (x + 1) * (x + t) * (x^4 + t)

f7 = (x + 1) * (x + t) * (x^2 + t) * (x^2 + tx + t)

f9 = (x + 1) * (x + t) * (x + t^2) * (x^3 + t*x^2 + t) 

f15 = (x + 1) * (x + t) * (x + t^2) * (x + t^3) * (x^2 + t) 

f27 = (x + 1) * (x + t) * (x + t^2) * (x + t^3) * (x + t^4)

pols = [f1, f2, f3, f5, f7, f15, f27]
for f in pols do
    c0 = coefficient(x^0, f)
    c1 = coefficient(x^1, f)
    c2 = coefficient(x^2, f)
    c3 = coefficient(x^3, f)
    c4 = coefficient(x^4, f)
    result = substitute(rep, {e1 => -c4, e2 => c3, e3 => -c2, e4 => c1, e5 => -c0})
end    
