Consider the arbitrary quartic equation \[ ax^4 + bx^3 + cx^2 + dx + e = 0 \] for real numbers $a$, $b$, $c$, $d$, $e$ with $a\neq0$. By the fundamental theorem of algebra this equation has four roots $x_1$, $x_2$, $x_3$, $x_4$ over the complex numbers. Using the factor theorem gives the factorization \[ ax^4 + bx^3 + cx^2 + dx + e = a(x-x_1)(x-x_2)(x-x_3)(x-x_4) . \] Expanding out the right-hand side gives \[ ax^4 - a(x_1+x_2+x_3+x_4)x^3 + a(x_1x_2+x_1x_3+x_1x_4+x_2x_3+x_2x_4+x_3x_4)x^2 - a(x_1x_2x_3+x_1x_2x_4+x_1x_3x_4+x_2x_3x_4)x + ax_1x_2x_3x_4 , \] and equating coefficients with the original expression gives the following system of equations: \begin{align*} b &= - a(x_1+x_2+x_3+x_4) \\ c &= a(x_1x_2+x_1x_3+x_1x_4+x_2x_3+x_2x_4+x_3x_4) \\ d &= - a(x_1x_2x_3+x_1x_2x_4+x_1x_3x_4+x_2x_3x_4) \\ e &= ax_1x_2x_3x_4 \end{align*} At this point, we're basically stuck; this is a complicated non-linear system that we want to solve for $x_1$, $x_2$, $x_3$, $x_4$.

## Some Sleight-of-hand

It turns out that it is now helpful to introduce the variables $y_1$, $y_2$, $y_3$ using the following definitions: \begin{align*} y_1 &= a(x_1 + x_2 - x_3 - x_4) \\ y_2 &= a(x_1 - x_2 + x_3 - x_4) \\ y_3 &= a(x_1 - x_2 - x_3 + x_4) \end{align*} And for the real bit of magic, I claim that the following identities hold: \begin{align*} y_1^2 + y_2^2 + y_3^2 &= 3b^2 - 8ac \\ y_1^2y_2^2 + y_1^2y_3^2 + y_2^2y_3^2 &= 3b^4 + 16a^2c^2 + 16a^2bd - 16ab^2c - 64a^3e \\ y_1y_2y_3 &= -b^3+4abc-8a^2d \end{align*} The real trick, of course, is where these identities came from—but once someone gives you them, it is not necessary to know how they were derived to check that they are true. It is straightforward, although tedious, to expand out the left sides using the definitions of $y_1$, $y_2$, and $y_3$, and to expand out the right sides using the expressions for $b$, $c$, $d$, and $e$. The result will be equations in terms of $x_1$, $x_2$, $x_3$, $x_4$, and $a$, and one merely needs to check that both sides are identical. [A real magician never reveals their secrets, but if you're curious I've also written a bit about how the identities were computed.]

## Finding the $y_i$

What make the $y_i$ so useful is that one can solve for them using the above identities. It's not immediately obvious how one would go about doing this, since again we have a complicated non-linear system of equations to solve. The trick is to consider a cubic equation which has as its roots $y_1^2$, $y_2^2$, and $y_3^2$; for example, $(y-y_1^2)(y-y_2^2)(y-y_3^2)$. Expanding this out gives \[ y^3 - (y_1^2+y_2^2+y_3^2)y^2 + (y_1^2y_2^2+y_1^2y_3^2+y_2^2y_3^2)y - (y_1y_2y_3)^2 , \] and the left-hand side of the identities make an appearance! Rewriting this using the right-hand sides gives \[ y^3 - (3b^2 - 8ac)y^2 + (3b^4 + 16a^2c^2 + 16a^2bd - 16ab^2c - 64a^3e)y - (-b^3+4abc-8a^2d)^2 , \] and now one can solve this using the cubic formula, and therefore find the roots $y_1^2$, $y_2^2$, and $y_3^2$. Taking square roots gives $y_1$, $y_2$, and $y_3$, except now a problem presents itself: since every number has two square roots, how do you know which one to take?

Actually, the problem isn't as bad as it seems. As long as one ensures that the values for $y_i$ chosen satisfy the third identity $y_1y_2y_3 = -b^3+4abc-8a^2d$, choosing different values for the square roots will just end up causing the $x_i$ to be labelled differently. So you don't have a completely free choice of square roots (in particular, one can't just use principal square roots), but 4 of the 8 possible square root selections will give correct answers.

## Finding the $x_i$

Now that the $y_i$ have been found we are in a great position to use them to find the $x_i$ and thus solve the quartic. Consider the definitions of the $y_i$ augmented with the expression for $b$: \begin{align*} b &= - a(x_1+x_2+x_3+x_4) \\ y_1 &= a(x_1 + x_2 - x_3 - x_4) \\ y_2 &= a(x_1 - x_2 + x_3 - x_4) \\ y_3 &= a(x_1 - x_2 - x_3 + x_4) \end{align*} We know $a$, $b$, $y_1$, $y_2$, and $y_3$, and wish to solve this system for $x_1$, $x_2$, $x_3$, and $x_4$. Since this is a nonsingular linear system, linear algebra allows us to do this! The solution (calculated using Gaussian elimation, for example) is as follows: \begin{align*} x_1 &= (-b+y_1+y_2+y_3)/(4a) \\ x_2 &= (-b+y_1-y_2-y_3)/(4a) \\ x_3 &= (-b-y_1+y_2-y_3)/(4a) \\ x_4 &= (-b-y_1-y_2+y_3)/(4a) \end{align*} Or more concisely, the solution can be expressed as \[ x = \frac{-b\pm(y_1\pm y_2)\pm y_3}{4a} \] where all choices of the $\pm$ signs are chosen with the last two equivalent.

Although this is a perfectly legitimate solution of the quartic, it relies on one "manually" choosing values for square roots so that $\sqrt{y_1^2}\sqrt{y_2^2}\sqrt{y_3^2} = -b^3+4abc-8a^2d$ is satisfied. For simplicity, I've also derived a slight modification which only uses principal square roots.

## References

Galois Theory by David A. Cox, in particular see Chapter 12.1C

Galois Theory by Harold M. Edwards, in particular see Section 17

Elements of Abstract Algebra by Allan Clark, in particular see Article 148

Abstract Algebra by Davis S. Dummit and Richard M. Foote, in particular see Chapter 14.7