What is a simple way to find real roots of a (cubic) polynomial?

Question

this seems like an obvious question to me, but I couldn't find it anywhere on SO. I have a cubic polynomial and I need to find real roots of the function. What is THE way of doing this?

I have found several closed form formulas for roots of a cubic function, but all of them use either complex numbers or lots of goniometric functions and I don't like them (and also don't know which one to choose).

I need something simple; faster is better; and I know that I will eventually need to solve polynomials of higher order, so having a numerical solver would maybe help too. I know I could use some library to do the hard work for me, but lets say I want to do this as an exercise.

I'm coding in C, so no import magic_poly_solver, please.

Bonus question: How do I find only roots inside a given interval?

Answer 1

For a cubic polynomial there are closed form solutions, but they are not particularly well suited for numerical calculus.

I'd do the following for the cubic case: any cubic polynomial has at least one real root, you can find it easily with Newton's method. Then, you use deflation to get the remaining quadratic polynomial to solve, see my answer there for how to do this latter step correctly.

One word of caution: if the discriminant is close to zero, there will be a numerically multiple real root, and Newton's method will miserably fail. Moreover, since to the vicinity of the root, the polynomial is like (x - x0)^2, you'll lose half your significant digits (since P(x) will be < epsilon as soon as x - x0 < sqrt(epsilon)). So you may want to rule this out and use the closed form solution in this particular case, or solve the derivative polynomial.

If you want to find roots in a given interval, check Sturm's theorem.

A more general (complex) algorithm for generic polynomial solving is Jenkins-Traub algorithm. This is clearly overkill here, but it works well on cubics. Usually, you use a third-party implementation.

Since you do C, using the GSL is surely your best bet.

Another generic method is to find the eigenvalues of the companion matrix with eg. balanced QR decomposition, or reduction to Householder form. This is the approach taken by GSL.