Algorithm for closed-form polynomial root finding

Question

I'm looking for a robust algorithm (or a paper describing an algorithm) that can find roots of polynomials (ideally up to the 4th debree, but anything will do) using a closed-form solution. I'm only interested in the real roots.

My first take on solving quadratic equations involved this (I also have code in similar style for cubics / quartics, but let's focus on quadratics right now):

/**
 *  @brief a simple quadratic equation solver
 *
 *  With double-precision floating-point, this reaches 1e-12 worst-case and 1e-15 average
 *  precision of the roots (the value of the function in the roots). The roots can be however
 *  quite far from the true roots, up to 1e-10 worst-case and 1e-18 average absolute difference
 *  for cases when two roots exist. If only a single root exists, the worst-case precision is
 *  1e-13 and average-case precision is 1e-18.
 *
 *  With single-precision floating-point, this reaches 1e-3 worst-case and 1e-7 average
 *  precision of the roots (the value of the function in the roots). The roots can be however
 *  quite far from the true roots, up to 1e-1 worst-case and 1e-10 average absolute difference
 *  for cases when two roots exist. If only a single root exists, the worst-case precision is
 *  1e+2 (!) and average-case precision is 1e-2. Do not use single-precision floating point,
 *  except if pressed by time.
 *
 *  All the precision measurements are scaled by the maximum absolute coefficient value.
 *
 *  @tparam T is data type of the arguments (default double)
 *  @tparam b_sort_roots is root sorting flag (if set, the roots are
 *      given in ascending (not absolute) value; default true)
 *  @tparam n_2nd_order_coeff_log10_thresh is base 10 logarithm of threshold
 *      on the first coefficient (if below threshold, the equation is a linear one; default -6)
 *  @tparam n_zero_discriminant_log10_thresh is base 10 logarithm of threshold
 *      on the discriminant (if below negative threshold, the equation does not
 *      have a real root, if below threshold, the equation has just a single solution; default -6)
 */
template <class T = double, const bool b_sort_roots = true,
    const int n_2nd_order_coeff_log10_thresh = -6,
    const int n_zero_discriminant_log10_thresh = -6>
class CQuadraticEq {
protected:
    T a; /**< @brief the 2nd order coefficient */
    T b; /**< @brief the 1st order coefficient */
    T c; /**< @brief 0th order coefficient */
    T p_real_root[2]; /**< @brief list of the roots (real parts) */
    //T p_im_root[2]; // imaginary part of the roots
    size_t n_real_root_num; /**< @brief number of real roots */

public:
    /**
     *  @brief default constructor; solves for roots of \f$ax^2 + bx + c = 0\f$
     *
     *  This finds roots of the given equation. It tends to find two identical roots instead of one, rather
     *  than missing one of two different roots - the number of roots found is therefore orientational,
     *  as the roots might have the same value.
     *
     *  @param[in] _a is the 2nd order coefficient
     *  @param[in] _b is the 1st order coefficient
     *  @param[in] _c is 0th order coefficient
     */
    CQuadraticEq(T _a, T _b, T _c) // ax2 + bx + c = 0
        :a(_a), b(_b), c(_c)
    {
        T _aa = fabs(_a);
        if(_aa < f_Power_Static(10, n_2nd_order_coeff_log10_thresh)) { // otherwise division by a yields large numbers, this is then more precise
            p_real_root[0] = -_c / _b;
            //p_im_root[0] = 0;
            n_real_root_num = 1;
            return;
        }
        // a simple linear equation

        if(_aa < 1) { // do not divide always, that makes it worse
            _b /= _a;
            _c /= _a;
            _a = 1;

            // could copy the code here and optimize away division by _a (optimizing compiler might do it for us)
        }
        // improve numerical stability if the coeffs are very small

        const double f_thresh = f_Power_Static(10, n_zero_discriminant_log10_thresh);
        double f_disc = _b * _b - 4 * _a * _c;
        if(f_disc < -f_thresh) // only really negative
            n_real_root_num = 0; // only two complex roots
        else if(/*fabs(f_disc) < f_thresh*/f_disc <= f_thresh) { // otherwise gives problems for double root situations
            p_real_root[0] = T(-_b / (2 * _a));
            n_real_root_num = 1;
        } else {
            f_disc = sqrt(f_disc);
            int i = (b_sort_roots)? ((_a > 0)? 0 : 1) : 0; // produce sorted roots, if required
            p_real_root[i] = T((-_b - f_disc) / (2 * _a));
            p_real_root[1 - i] = T((-_b + f_disc) / (2 * _a));
            //p_im_root[0] = 0;
            //p_im_root[1] = 0;
            n_real_root_num = 2;
        }
    }

    /**
     *  @brief gets number of real roots
     *  @return Returns number of real roots (0 to 2).
     */
    size_t n_RealRoot_Num() const
    {
        _ASSERTE(n_real_root_num >= 0);
        return n_real_root_num;
    }

    /**
     *  @brief gets value of a real root
     *  @param[in] n_index is zero-based index of the root
     *  @return Returns value of the specified root.
     */
    T f_RealRoot(size_t n_index) const
    {
        _ASSERTE(n_index < 2 && n_index < n_real_root_num);
        return p_real_root[n_index];
    }

    /**
     *  @brief evaluates the equation for a given argument
     *  @param[in] f_x is value of the argument \f$x\f$
     *  @return Returns value of \f$ax^2 + bx + c\f$.
     */
    T operator ()(T f_x) const
    {
        T f_x2 = f_x * f_x;
        return f_x2 * a + f_x * b + c;
    }
};

The code is horrible, and I hate all the thresholds. But for random equations with roots in the [-100, 100] interval, this is not so bad:

root response precision 1e-100: 6315 cases
root response precision 1e-19: 2 cases
root response precision 1e-17: 2 cases
root response precision 1e-16: 6 cases
root response precision 1e-15: 6333 cases
root response precision 1e-14: 3765 cases
root response precision 1e-13: 241 cases
root response precision 1e-12: 3 cases
2-root solution precision 1e-100: 5353 cases
2-root solution precision 1e-19: 656 cases
2-root solution precision 1e-18: 4481 cases
2-root solution precision 1e-17: 2312 cases
2-root solution precision 1e-16: 455 cases
2-root solution precision 1e-15: 68 cases
2-root solution precision 1e-14: 7 cases
2-root solution precision 1e-13: 2 cases
1-root solution precision 1e-100: 3022 cases
1-root solution precision 1e-19: 38 cases
1-root solution precision 1e-18: 197 cases
1-root solution precision 1e-17: 68 cases
1-root solution precision 1e-16: 7 cases
1-root solution precision 1e-15: 1 cases

Note that this precision is relative to the magnitude of the coefficients, which is typically in the 10^6 range (so finally the precision is far from perfect, but probably mostly usable). Without the thresholds, however, it is near to useless.

I have tried using multiple precision arithmetics, which generally works well, but tends to reject many of the roots simply because the coefficients of the polynomial are not multiple precision and some polynomials cannot be exactly represented (if there is a double root in a 2nd degree polynomial, it mostly either splits it to two roots (which I wouldn't mind) or says that there is no root whatsoever). If I want to recover perhaps even slightly imprecise roots, my code gets complicated and full of thresholds.

So far, I've tried using CCmath, but either I can't use it correctly, or the precision is really bad. Also, it uses iterative (not closed-form) solver in plrt().

I have tried using GNU scientific library gsl_poly_solve_quadratic() but that seems to be a naive approach, and not very numerically stable.

Using std::complex numbers naively also turned out to be a really bad idea, as both the precision and speed can be bad (especially with cubic / quartic equations where the code is heavy with transcendental functions).

Is recovering the roots as complex numbers the only way to go? Then no roots are missed and the user can select how precise the roots need to be (and thus ignore small imaginary components in less precise roots).

at school, we used Bairstow's method to extract roots link Maybe it's what you need... — iksess
@iksess Bairstow's method you propose is fairly interesting (thanks), but it uses Newton's method, which is not closed-form. — the swine
Sorry ! My english is not good enough to understand what is "closed-form" — iksess
@iksess No problemo. Closed-form is an expression that can be evaluated in a finite number of steps, which generally cannot be said about the Newton's method, as it might not converge at all. On the other hand, the code I posted is closed-form, as it does not even loop, and therefore always finishes in a fixed number of steps. — the swine

Ian Cook Ian Cook · Accepted Answer · 2015-01-15T22:43:25

This isn't really answering your question but I think you can improve on what you've got since you currently have a 'loss of significance' problem when b^2 >> ac. In such cases, you end up with a formula along the lines of (-b + (b + eps))/(2 * a) where the cancellation of the b's can lose many significant figures from eps.

The correct way of handling this is to use the 'normal' equation for roots of a quadratic for one root and the lesser known 'alternative' or 'upside down' equation for the other root. Which way round you take them depends on the sign of _b.

A change to your code along this lines of the following should reduce the errors resulting from this.

if( _b > 0 ) {
    p_real_root[i] = T((-_b - f_disc) / (2 * _a));
    p_real_root[1 - i] = T((2 * _c) / (-_b - f_disc));
}
else{
    p_real_root[i] = T((2 * _c) / (-_b + f_disc));
    p_real_root[1 - i] = T((-_b + f_disc) / (2 * _a));
}

Algorithm for closed-form polynomial root finding

1 Answers