Flocke developed an algorithm for obtaining all the zeros of cubic and quartic polynomials. The key to accuracy is scaling the polynomials so that all coefficients in absolute value are bounded by unity. A recent book by Boyd [1] contains a chapter on the zeros of cubic and a chapter on the zeros of quartic polynomials.
For the cubic solver, the algorithm considers six classes based on the sign of each coefficient. For each class, one can find a good starting value for Newton’s method to converge fast to the desired root.
For the quartic solver, the algorithm uses the cubic solver to find the three stationary points. In the cases of real roots, Newton’s method is used to get the extreme real root. Deflation is used to get the other roots. If the quartic has only complex roots, a sixth degree polynomial is constructed to obtain the real parts of those roots and a quadratic polynomial is used to get the imaginary parts.
Examples showing the efficiency and accuracy of the solver for both the cubic and quartic polynomials are given. It would be interesting to compare this algorithm to the one suggested by Boyd.