R.<z4,z3,z2,z1,z0,y4,y3,y2,y1,y0,x,zeta5,sqrtm10m2sqrt5,sqrtm10p2sqrt5,sqrt5,rt5of2> = PolynomialRing(QQ,['z4','z3','z2','z1','z0','y4','y3','y2','y1','y0','x','zeta5','sqrtm10m2sqrt5','sqrtm10p2sqrt5','sqrt5','rt5of2'],order='degrevlex(5),degrevlex(5),lex(6)')
Izeta = R.ideal(sum([zeta5^i for i in range(5)]))
Ifullnum = R.ideal([zeta5 - 1/8*sqrtm10p2sqrt5*sqrt5 - 1/8*sqrtm10p2sqrt5 - 1/4*sqrt5 + 1/4, sqrtm10m2sqrt5 - 1/2*sqrtm10p2sqrt5*sqrt5 - 1/2*sqrtm10p2sqrt5, sqrtm10p2sqrt5^2 - 2*sqrt5 + 10, sqrt5^2 - 5])

f = x^5 + 100*x^2 + 1000
elt = SymmetricFunctionAlgebra(QQ, basis='elementary')
sigma1 = elt([1]).expand(5).subs(x0=z0,x1=z1,x2=z2,x3=z3,x4=z4)
sigma2 = elt([2]).expand(5).subs(x0=z0,x1=z1,x2=z2,x3=z3,x4=z4)
sigma3 = elt([3]).expand(5).subs(x0=z0,x1=z1,x2=z2,x3=z3,x4=z4)
sigma4 = elt([4]).expand(5).subs(x0=z0,x1=z1,x2=z2,x3=z3,x4=z4)
sigma5 = elt([5]).expand(5).subs(x0=z0,x1=z1,x2=z2,x3=z3,x4=z4)
sigmas = [1,sigma1,sigma2,sigma3,sigma4,sigma5]
rootrels = [sigmas[i] - (-1)^i * f.coefficient({x:f.degree()-i}) for i in range(1,len(sigmas))]
Idecomp = R.ideal(rootrels)
zs = [z0,z1,z2,z3,z4]
ys = [y0,y1,y2,y3,y4]
symp = sum([zs[i]^2 * (zs[(i+1)%5]*zs[(i+4)%5] + zs[(i+2)%5]*zs[(i+3)%5]) for i in range(5)])
symq = sum([zs[i] * zs[(i+1)%5] for i in range(5)])
Idecp = Idecomp + R.ideal(symp-2000)
Idecq = Izeta + Idecp + R.ideal(symq-20*(2*(1 + zeta5 + zeta5^4) - 1))
lagrange = [sum([zs[i]*zeta5^((i*j)%5) for i in range(5)]) for j in range(5)]
Iresolv0 = Izeta + Idecq + R.ideal([ys[i]-lagrange[i] for i in range(5)]) + R.ideal([y0, y1^2*y3 - 1000, y1*y3^3 + 5000, y2^5 + 400000, y2*y4^2 + 1000, y1*y4 - 100, y3^2*y4 + 500, y2*y3 + 100, y3*y4^3 + 5000, y4^5 + 50000, y1^5 + 200000, y1^3*y2 - 20000, y3^5 - 25000, y2^3*y4 - 20000, y1*y2^2 + 2000])
B0 = Iresolv0.groebner_basis()
Iresolv = R.ideal(B0) + Ifullnum
B = Iresolv.groebner_basis()
Iresolv2 = R.ideal(B) + R.ideal(rt5of2 + y1/10)
B2 = Iresolv2.groebner_basis()