The Euler method for quartic equations

Consider the (depressed) quartic equation

    z4 + pz2 + qz + r = 0

for complex coefficients. The following solution method is due to Euler [1729]. The basic idea is
to seek the roots as four plain linear combinations of three square roots. Employing the ansatz
z = u + v + w yields

z2 = (u2 + v2 + w2) + 2(uv + uw + vw)

z4 = (u2 + v2 + w2)2 + 4(u2 + v2 + w2)(uv + uw + vw) +
       4(u2v2 + u2w2 + v2w2) + 8uvw(u + v + w)

which after substitution into the original equation results in

    (u2 + v2 + w2)2 + p(u2 + v2 + w2) + (4(u2 + v2 + w2) + 2p)(uv + uw + vw) +
    4(u2v2 + u2w2 + v2w2) + (8uvw + q)(u + v + w) + r = 0

By first requiring

(1)    4(u2 + v2 + w2) + 2p = 0

(3)    8uvw + q = 0

and subsequently exploiting (1), the remaining equation becomes

(2)    0 = (u2 + v2 + w2)2 + p(u2 + v2 + w2) + 4(u2v2 + u2w2 + v2w2) + r
            = (−p/2)2 + p(−p/2) + 4(u2v2 + u2w2 + v2w2) + r
            = 4(u2v2 + u2w2 + v2w2) − p2/4 + r

Introducing t1 = u2, t2 = v2, t3 = w2, equations (1–3) are reformulated as

(1.2)    t1 + t2 + t3 = −p/2
(2.2)    t1t2 + t1t3 + t2t3 = p2/16 − r/4
(3.2)    t1t2t3 = q2/64

According to the Viète formulas, t1, t2, t3 are the roots to the cubic resolvent

(4)    t3 + (p/2)t2 + (p2/16 − r/4)tq2/64 = 0

Thus, the following possibilities arise u = ∓√t1, v = ∓√t2, w = ∓√t3 from which eight linear
combinations for the ansatz (u + v + w) could be formed. Regardless of the choices of sign,
the product uvw becomes either √t1t2t3 or -√t1t2t3. Consequently, it is convenient to choose

    u = √t1,   v = √t2,   w = √t3

and use (3) as the discriminator between case 1 (8uvw + q = 0) or case 2 (-8uvw + q = 0).
For case 1, the appropriate linear combinations are u+v+w, u-v-w, -u+v-w, -u-v+w
(each product of the three terms is uvw). For case 2, the appropriate linear combinations are
-u-v-w, -u+v+w, u-v+w, u+v-w (each product of the three terms is -uvw).

For the case q ≠ 0, the cubic resolvent (4) is depressed by the change of variable t = sp/6 into

    s3 + 3bs − 2d = 0
where
    b = −r/12 − p2/144
    d = q2/128 − pr/48 + p3/1728

the roots of which yield

    u = √(s1p/6),   v = √(s2p/6),   w = √(s3p/6)

The roots to the original quartic equation become
if  |8uvw + q| < |8uvwq|
    z1 =   u + v + w
    z2 =   uvw
    z3 = −u + vw
    z4 = −uv + w
else
    z1 = −uvw
    z2 = −u + v + w
    z3 =   uv + w
    z4 =   u + vw

For the biquadratic case q = 0, the cubic resolvent (4) is factorized into

    (t2 + (p/2)t + p2/16 − r/4)t = 0

the roots of which are
    t1 = −p/4 + √(r/4)
    t2 = −p/4 − √(r/4)
    t3 = 0
yielding w = 0 and

    u = √(−p + 2√r)/2,   v = √(−p − 2√r)/2

The roots to the original quartic equation simply become
    z1 =   u + v
    z2 =   uv
    z3 = −u + v
    z4 = −uv