# Quartic Solution on Wikipedia special cases problem $S=0$ how to "change the choice of cubic root"?

So, I've posted a question regarding Wikipedia's quartic page. This was from the first question.

I'm trying to implement the general quartic solution for use in a ray tracer, but I'm having some trouble. The solvers I've found do cause some strange false negatives leaving holes in the tori I'm testing with.

Most implementations use the depressed quartic solutions, I don't understand the math involved and can't figure out why I'm having false non-intersections (link to layman explanation would be great). So I'm trying to implement the general solution at this wikipedia page. I got the stuff up until the special cases implemented, but at that point I have an issue.

With lots of rays being traced most of the special cases become common. I've found a set of coefficients that Wolfram Alpha tells me has two real roots, but my code was just returning NaN, further searching I found my S was coming up as $$\sqrt{-4.9 \times 10^{-11}}$$ Floating point precision error, means this should equate to 0, so I need the special case for S=0, it says we need to "change choice of cubic root in Q" but it does not explain how to do this. I did try changing the sign of Q when S=0, but that doesn't work. Does anyone know what this means and how I can do it?

• For most practical purposes, the "exact" solution of the quartic in radicals is not very useful: besides being complicated, it can have severe roundoff errors. A more useful idea might be to isolate the roots in intervals using Sturm's theorem and then approximate them using Newton-Raphson. Jun 17, 2015 at 22:21
• @RobertIsrael interesting point, I've already seen some round-off errors as differences between my discriminant and WolframAlpha's, fixed by going from calculating the 16 terms outright, to using the $\Delta_1^2-4\Delta_0^3=-27\Delta$ identity. Do you have a good link for the Strum's/Newton-Raphson? Jun 18, 2015 at 2:04
• @MatrixPeckham en.wikipedia.org/wiki/Root-finding_algorithm is a good starting point. I second the recommendation that you may as well use a rootfinding algorithm - especially since whatever square root/cube root algorithm you're finding is doing that behind the scenes anyway. Jun 18, 2015 at 4:54

Try this version. Given,

$$x^4+ax^3+bx^2+cx+d=0$$

then,

$$x_{1,2} = -\tfrac{1}{4}a+\tfrac{1}{2}\sqrt{u}\pm\tfrac{1}{4}\sqrt{3a^2-8b-4u+\frac{-a^3+4ab-8c}{\sqrt{u}}}\tag1$$

$$x_{3,4} = -\tfrac{1}{4}a-\tfrac{1}{2}\sqrt{u}\pm\tfrac{1}{4}\sqrt{3a^2-8b-4u-\frac{-a^3+4ab-8c}{\sqrt{u}}}\tag2$$

where,

$$u = \frac{3a^2-8b}{12} +\frac{1}{3}\left(v_1^{1/3}+\frac{b^2 - 3 a c + 12 d}{v_1^{1/3}}\right)$$

and $v_1$ is any non-zero root of the quadratic,

$$v^2 + (-2 b^3 + 9 a b c - 27 c^2 - 27 a^2 d + 72 b d)v + (b^2 - 3 a c + 12 d)^3 = 0$$

P.S. This is essentially the method used by Mathematica, though much simplified for aesthetics.

• @MatrixPeckham: How does this method fare on your ray tracer implementation? Jun 18, 2015 at 14:22
• It might work, but after adding complex numbers to my existing code I think I have it working, but need to test more. Jun 19, 2015 at 3:17

If $\alpha$ is a root of $x^3 = A$, then so too is $\gamma \alpha$, where $\gamma$ is a complex cube root of unity. This follows since $(\gamma\alpha)^3 = \gamma^3\alpha^3 = \alpha^3 = A$.

The nontrivial values of these are given by $\gamma = -\frac{1}{2} + \frac{i\sqrt{3}}{2}$ or $\gamma = -\frac{1}{2} - \frac{i\sqrt{3}}{2}$.

In your case, you need to replace $Q$ with $\gamma Q = \left(-\frac{1}{2} + \frac{i\sqrt{3}}{2}\right)Q$ in this case, if I'm reading Wikipedia's algorithm correctly.

(Implementation might be easier if you used the note that this case is always accompanied by the depressed quartic being biquadratic, in which case the solutions follow from applying the quadratic formula.)

• This would require Imaginary numbers, which my implementation language does not natively support, The biquadratic method requires a re-write of the rest of my code... I may implement the imaginary numbers just to make everything consistent with the overall formula. Jun 18, 2015 at 2:11
• Even for cubic equations, you can't get by without complex numbers. That is, if an irreducible (over the rationals) cubic has three real roots, those roots can't be expressed using only real radicals. This would extend to a quartic that is the product of a linear factor with rational root and one of those cubics. See en.wikipedia.org/wiki/Casus_irreducibilis Jun 18, 2015 at 22:36
• That is true, however my need is for real roots only (ray-torus intersection) so I was hoping to scrape by without, but I have since implemented complex numbers and I think my code is working, however I don't know for sure, I need some more test cases with known solutions. Jun 19, 2015 at 3:13

As noted in a comment by Robert Israel, if you are using floating point numbers, it's better to use a root-finding algorithm than the exact formula for quartic roots, as root finding algorithms tend to have much better numerical stability, and avoid issues like this one. They also tend to work much more generally (i.e., they won't be limited to just a single degree of polynomial, or just polynomials of degree < 5).

If you are interested in the various ways of representing the quartic solution symbolically, the SymPy code goes through a few methods and tries to pick the best one for the given polynomial. There are also some useful references there for further reading.