Fastest numerical solution of a real cubic polynomial?

Question

R question: Looking for the fastest way to NUMERICALLY solve a bunch of arbitrary cubics known to have real coeffs and three real roots. The polyroot function in R is reported to use Jenkins-Traub's algorithm 419 for complex polynomials, but for real polynomials the authors refer to their earlier work. What are the faster options for a real cubic, or more generally for a real polynomial?

Answer 1

The numerical solution for doing this many times in a reliable, stable manner, involve: (1) Form the companion matrix, (2) find the eigenvalues of the companion matrix.

You may think this is a harder problem to solve than the original one, but this is how the solution is implemented in most production code (say, Matlab).

For the polynomial:

p(t) = c0 + c1 * t + c2 * t^2 + t^3

the companion matrix is:

[[0 0 -c0],[1 0 -c1],[0 1 -c2]]

Find the eigenvalues of such matrix; they correspond to the roots of the original polynomial.

For doing this very fast, download the singular value subroutines from LAPACK, compile them, and link them to your code. Do this in parallel if you have too many (say, about a million) sets of coefficients.

Notice that the coefficient of t^3 is one, if this is not the case in your polynomials, you will have to divide the whole thing by the coefficient and then proceed.

Good luck.

Edit: Numpy and octave also depend on this methodology for computing the roots of polynomials. See, for instance, this link.