var('x,y,z,w')
R=x+y-z-w
s=x-y+z-w
t=x-y-z+w
h(x,y,z,w)=R^2+s^2+t^2+R^2*s^2+R^2*t^2+s^2*t^2+R*s*t
r(a,b,c,d)=0
for i in range(7):
    g(x,y,z)=h(x,y,z,0)
    q(a,b,c)=0
    for n in range(8):
        f(x,y)=g(x,y,0)
        p(a,b)=0
        for m in range(11):
            p(a,b)=p(a,b)+f(a,0)*b^m
            f(x,y)=expand((f(x,y)-f(x+y,0))/x/y)
        q(a,b,c)=q(a,b,c)+p(a,b)*c^n
        g(x,y,z)=expand((g(x,y,z)-p(x+y+z,x*y+x*z+y*z))/x/y/z)
    r(a,b,c,d)=r(a,b,c,d)+q(a,b,c)*d^i
    h(x,y,z,w)=expand((h(x,y,z,w)-q(x+y+z+w,x*y+x*z+x*w+y*z+y*w+z*w,x*y*z+x*y*w+x*z*w+y*z*w))/x/y/z/w)
print (r(a,b,c,d))
#The linear combinations of the roots of a quartic to be solved are all those with two coefficients 1 and two coefficients -1, normalized to first coefficient 1 by scalar multiplication not affecting the square so the squares are all the squares and are permuted by permutation of the variables. Thus h is symmetric and a symmetric polynomial program writes h in terms of the elementary symmetric polynomials whence the cubic with the squares for roots is known. Solving it yields our linear combinations which along with the first elementary symmetric polynomial add to 4x. Solving the cubic is similar with two linear combinations, the two with distinct third roots of unity first one normalized to 1 by scalar multiplication. Here the scalar multiplications do not affect the cube so the two cubes are all the cubes and are permuted. Thus the quadratic with the two cubes for roots has symmetric coefficients and a program evaluates its coefficients. Our program embeds the three variable program in lines 9 through 17 and attacks our symmetric polynomial with the last variable set equal to 0 with solution q which is used to redefine h in line 19. Setting w=0 vanishes the numerator, hence the even division by w. The other divisions are even because the numerator is symmetric. Line 18 adjusts oue solution for the subtraction with the d^i factor adjusting for the i earlier divisions by the last elementary symmetric polynomial. Algorithm programs for 2 through 6 variables are in cocalc.com/share. The output yields resolvent x^3-(3a^2-8b)x^2+(3a^4-16ba^2+16b^2+16ac-64d)x-(a^3-4ab+8c)^2 for x^4-ax^3+bx^2-cx+d and the cubic output yields resolvent x^2-(2a^3-9ab +27c)x+(a^2-3b)^3 for x^3-ax^2+bx-c.