The tedious way is to expand $(x_1 - x_2)^2 (x_2 - x_3)^2 (x_1 - x_3)^2$ out completely and then write it in terms of $x_1 + x_2 + x_3$ and so forth. You are guaranteed that this is possible by the fundamental theorem of symmetric polynomials, the proof of which even gives an algorithm for doing this, but it's a pain to do by hand (although it's not a bad exercise in algebraic manipulation).
A less tedious way is to argue as follows. We will first work under the assumption that $p = 0$. Now, $q, r$ are polynomials of degrees $2, 3$, and the discriminant is a polynomial of degree $6$, so the discriminant must be a linear combination of the two monomials $q^3, r^2$. Thus we can write
$$\Delta = a q^3 + b r^2$$
for two constants $a, b$, where $\Delta$ is the discriminant. Setting $q = -1, r = 0$ we get the polynomial $x^3 - x = 0$ with roots $0, 1, -1$. We compute that the discriminant is equal to $4$, from which it follows that $a = -4$.
Setting $q = 0, r = -1$ we get the polynomial $x^3 - 1 = 0$ with roots $1, \omega, \omega^2$ where $\omega$ is a primitive third root of unity. Using the identity
$$(\omega - 1)^2 = \omega^2 - 2 \omega + 1 = - 3 \omega$$
we compute that the discriminant is equal to $-27$, from which it follows that $b = -27$. Thus
$$\Delta = -4 q^3 - 27 r^2.$$
To get from here to an arbitrary choice of $x_1, x_2, x_3$, apply the above formula to the polynomial with roots $x_1 - \frac{p}{3}, x_2 - \frac{p}{3}, x_3 - \frac{p}{3}$ and note that subtracting the same constant from each of the three roots doesn't change the discriminant.
The problem occurs when you are taking a cube root: in the Wolfram version, when you compute $$T = \sqrt[3]{R - \sqrt{D}} = \sqrt[3]{-\frac{35}{27}-\sqrt{\frac{50}{27}}}$$ and in corresponding places in the other versions of the cubic formula.
The expression inside the cube root is negative (approximately $-2.66$), and the cubic formula works if you take the real cube root (approximately $-1.39$). But writing $(R - \sqrt D)^{1/3}$ in many computer algebra systems instead computes the principal cube root: the complex root with largest real part. This will be the real cube root for a positive number, but here, it gives us approximately $0.69 - 1.2 i$, and that is the source of your error.
I can't identify the language you're using by looking at it, so I don't know what the best way to avoid this problem is. In Mathematica, there's a built-in CubeRoot command, which will always take the real root of a real number.
You can always make it go with some conditionals: define $T$ to be
- $(R - \sqrt{D})^{1/3}$, if $R - \sqrt D \ge 0$, and
- $-(\sqrt D - R)^{1/3}$, if $R - \sqrt D < 0$.
(And do a similar thing everywhere else you take a cube root.)
The above will work for real coefficients; for complex coefficients (or even when $D$ is negative), we want to be more careful, because then a real root might not exist.
The idea is that we ran into trouble here, not because we weren't taking the real root necessarily, but because we took two cube roots that "didn't match up with each other". We can rewrite the cubic formula in a few different ways so that we only take one cube root, and express the other cube root we'd have to take in terms of the first. Then the issue doesn't occur.
This gives us a better way to solve the problem in code: after computing $S = (R + \sqrt D)^{1/3}$, set $T = -Q/S$. The justification is this: since $S^3 = R + \sqrt D$ and $T^3 = R - \sqrt D$, we have $S^3 T^3 = R^2 - D = -Q^3$. Naively, we can cancel the cube roots and assume $ST = -Q$, and this works out.
See also this question and its answer.
Best Answer
Using the formula $\cos(3\theta) = 4\cos^3(\theta) - 3\cos(\theta)$ is actually quite useful, but first we have to do $x=2y$, then we get $$8y^3 - 6y - 1 = 0$$ Now we have the 4:-3 ratio and we can substitute $y=\cos(\theta)$ to get, from the above formula, $$2\cos(3\theta) = 1$$ and from there we can get cosine form of the roots.