What are the best methods for solving cubic and quartic equations by computer programs?

We know that there are closed form formulas for real roots of degree 4 and 3 polynomials, but people sometimes advise to use numerical (e.g. Newton) methods anyway. They claim that closed formulas (Cardano or Viéte...) lack robustness.

What is the best approach to this problem, for a programmer?

• Which library (possibly open source) do you recommend?
• Iterative or algebraic approach?
• What are the problems with algebraic methods?
-
"They claim that closed formulas lacks robustness. "- because people are often sloppy that they pay no attention to things like subtractive cancellation. Even the quadratic formula is not trivial to implement as a stable formula. –  Guess who it is. Jul 26 '12 at 11:41
"How about "speed" concerns?" - you certainly have to do testing yourself; it depends. FWIW, iterative methods are not very useful if you don't have good starting points or bounds, so that's one point that complicates their use. –  Guess who it is. Jul 26 '12 at 11:42
Note that "algebraic" approaches require drawing roots, which isn't an elementary operation on a computer; the number of elementary operations required for a well-chosen numerical method for a cubic or quartic equation may not be all that much more than what's required for drawing roots, and the closed formulas typically contain several roots. –  joriki Jul 26 '12 at 11:58
Considering that even such a simple operation as division is often implemented in software - as a newton iteration - the lines between all these techniques are blurry. –  Peter Sheldrick Jul 26 '12 at 12:12
Although you can probably stabilize the closed-form representations of the roots, I suspect that in doing so you would likely end up with no fewer operations, anyway. The general representations for cubics and quartics are pretty heinous anyways. –  Arkamis Jul 26 '12 at 12:45

I might as well expand on the comments I gave.

Firstly, note that square roots and cube roots, and even division in itself are computed through Newton-Raphson or some other iterative method, so the dichotomy of "iterative" versus "algebraic" methods in implementing these on a computer are more or less moot; they're all iterative.

Secondly, by virtue of any of these things being iterative, there is no surefire way to say whether directly applying the quadratic formula (to use the simplest case) or using Newton-Raphson on $ax^2+bx+c=0$ would be quicker. (Also, compilers do the darndest things, so always do some testing.) The little snag that these people who keep pushing "iterative" over "algebraic" approaches apparently forget to talk about is that Newton-Raphson only works out nicely if you have a good starting point. If you're very much concerned with efficiency, you can't afford to be a slouch in choosing the starting point, and that's not always easy to do automatically.

It must be noted, however, that there are certain modifications of Newton-Raphson that do not require accurate starting values. These usually fall under the rubric of "simultaneous iteration methods", the prototypical example of which is the (Weierstrass-)Durand-Kerner method. This is effectively the application of Newton-Raphson to the Vieta formulae relating the roots and coefficients of a polynomial. (See this for more details.)

But, getting back to the topic at hand, even the numerically sound implementation of the classical algebraic formulae is not completely trivial, and most people are careless, writing lousy code that does not take into account things like subtractive cancellation. To again use the quadratic formula as an example:

$$x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$$

consider the case where $b^2\gg4ac$, with which the square root of the discriminant is nearly the same in magnitude as $b$. Then, the computation corresponding to the plus sign involves what is called subtractive cancellation, where one is subtracting nearly equal quantities to obtain a tiny result. This is a good way to lose significant figures in your results.

The right way starts by first considering the "Citardauq" formula:

$$x=\frac{2c}{-b\pm\sqrt{b^2-4ac}}$$

In that case, one of the roots is accurately computed, and the other one is the one prone to subtractive cancellation. The good news is that the two formulae are complementary: one is able to accurately compute the root that the other method is not too good at. Thus, if one wants to stably implement the "algebraic approach", one should first compute

$$q=-\frac12\left(b+\mathrm{sign}(b)\sqrt{b^2-4ac}\right)$$

(assuming that the coefficients are all real; if the coefficients are complex, things are slightly more intricate), and after which, the two roots of $ax^2+bx+c=0$ are computed as $q/a$ and $c/q$. Similar stabilized rearrangements are possible for the cubic formula (both the trigonometric form that is needed for casus irreducibilis and the Cardano form that you can use when only one of the roots is real), and then one has to note that the solution of a quartic hinges on the solution of a (resolvent) cubic and a quadratic.

In short: you still have to do the testing for determining which is "best" for your application.

-
Thank you for answer. I am aware of problems that comes with numerical competitions, but I'm not good at solving them. That's why I would prefer tested library that computes those roots, rather then implementing it myself. Do you know any? –  Piotr Jul 26 '12 at 17:29
I don't know of any computing environments that implement the specialized methods for quadratics, cubic, and quartics. Most of them just have an iterative routine that works for polynomials of all degrees, so they don't bother with "special cases". –  Guess who it is. Jul 27 '12 at 0:22
If you reformualte these problems as eigenvalue computations of small matrices than the Fortran library LAPACK is a solid choice. –  Peter Sheldrick Jul 27 '12 at 1:33
@Peter mentions one possible iterative route; they are pretty good with respect to accuracy, but are slightly expensive with respect to effort. –  Guess who it is. Jul 27 '12 at 11:11