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? 

*How about "speed" concerns?  

 A: 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.
