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

Question

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?

Answer 1

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.

Answer 2

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.)

Answer 3

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).

Answer 4

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.