I have an application that requires finding the roots of a cubic polynomial, so I programmed a closed-form cubic root finder based on the algorithm on Wikipedia. The formula is mathematically impressive, and it works fairly well, but there is a potential problem with it.
The three roots of a cubic polynomial consist of a real root and two other roots that can be either real or complex. The problem is that numerical roundoff error can cause the real roots to appear to have a negligible but nonzero imaginary part. My application requires the smallest positive real root, so I need to identify the real roots.
Since we know that one root must be real, we can take the root with the smallest imaginary part and zero that component, making the root real. The problem is that for the other two roots, some threshold must be selected for determining whether the root is real or imaginary.
Using standard 64-bit “Doubles”, I calculate the absolute value of the ratio of the imaginary part to the real part, and set the threshold to 1e-9. That may be a bit conservative, but it seems to work for my application.
For higher reliability, I would like to compute the roots using higher numerical precision than standard “Doubles”. I looked into BigDecimal, but it has no square root or cube root functions, which are needed for the algorithm. I guess I could roll my own square root and cube root functions, but I am surprised that one is not available in the standard library.
Any suggestions for an alternative to BigDecimal or where to find the code for square root and cube root functions?
Beat me to it – yeah, Spire provides several high-precision numeric types, and type classes providing various functionality. I don’t do a lot with numerics myself, but that’s where I would look first.
OK, thanks. Spire looks impressive, and maybe I should learn how to use it. For now though, does anyone know if it provides sqrt and cbrt functions for BigDecimal? Also, can the “Real” type be used more or less as a plug-in replacement for “Double”?
If your cubic equation has real-valued coefficients, then your genuinely complex roots will be conjugate. So you could partition the roots by the sign of their imaginary parts and pair the complex roots off by seeing which of the two imaginary parts are closest to the negative of the other imaginary part.
The real parts of the two conjugate roots across the partition will also be closest together too.
Yes, my coefficients are real valued. So are you saying that the complex roots, if there are any, must be a conjugate pair, with identical real parts and opposite imaginary parts? Or could there be tiny differences? And if there are still differences, don’t I still need a threshold?
Yes. Yes (numerical error). No - just choose the ones whose real parts are closest and whose negated imaginary parts are closest. As a final tiebreaker, substitute just the real parts back into the cubic and choose the one with the least magnitude residual error.
Yes, I saw some of that stackoverflow history, but I guess I was hoping that one of the wizards here could suggest some simpler trick.
After reading the replies here and thinking about it a bit, I realize that the problem not as “complex” as I thought it was (pun intended).
The plot of a cubic polynomial is like two parabolas in opposite directions spliced together into a smooth curve that can have a local minimum and a local maximum. The roots are where it crosses the x axis, which can occur at one, two, or three places.
If it crosses at two places, that is a rare special case where it crosses the x axis exactly at either the local minimum or the local maximum. In that case, two real roots will be identical. For roots that have a negligible but nonzero imaginary component, I should be able to figure out what is going on by looking at the difference between the real components. If the real components are close together, the roots are at or near the local min or max.
I don’t think I need the higher precision numerics after all.