As suggested in my comment: you can eliminate x1,…,xt from the ideal I=⟨yi−ei(p1(x1,…,xt),…,pm(x1,…,xt)):i=1,…,m⟩ in the ring Q[x1,…,xt,y1,…,ym], where ei are the elementary symmetric polynomials in m variables.
For efficiency reasons we compute the elimination ideal I∩Q[y1,…,ym] manually, using a Groebner basis G of I with respect to a block ordering where the x's are the greatest and the y's have a weighted degrevlex ordering where yk has degree k (simulating the degree of ek), so that substituting yi=ei(x1,…,xt) into the polynomials in G∩Q[y1,…,ym] yields symmetric dependencies of minimal degree.
def symmetric_dependencies(p):
R = p[0].parent()
m = len(p)
S = PolynomialRing(R.base_ring(),
names=list(R.gens()) + ['y{}'.format(k+1) for k in range(m)],
order=TermOrder('degrevlex', R.ngens()) + \
TermOrder('wdegrevlex',list(range(1,m+1))))
p = [S(f) for f in p]
X = S.gens()[0:R.ngens()]
Y = S.gens()[R.ngens():]
E = SymmetricFunctions(QQ).elementary()
Es = [E[k+1].expand(m) for k in range(m)]
I = S.ideal([Y[k] - Es[k](p) for k in range(m)])
G = I.groebner_basis()
GY = [g for g in G if all(v in Y for v in g.lt().variables())]
E_subs = {Y[k] : Es[k] for k in range(m)}
dependencies = [f.subs(E_subs) for f in GY]
#assert all(f(p) == 0 for f in dependencies)
return dependencies
Example:
sage: R.<x> = QQ[]
sage: [min(f.degree() for f in symmetric_dependencies([x, x^2, x^(2*k)])) for k in range(1,8)]
[4, 5, 8, 11, 12, 12, 12]
It seems to be slightly faster than your function on this example.