R. = 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()