title: Algorithm for closed-form polynomial root finding
body: <p>I'm looking for a <em>robust</em> 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.</p>
  
  <p>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):</p>
  
  <pre><code>/**
   *  @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 &lt;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&gt;
  class CQuadraticEq {
  protected:
      T a; /**&lt; @brief the 2nd order coefficient */
      T b; /**&lt; @brief the 1st order coefficient */
      T c; /**&lt; @brief 0th order coefficient */
      T p_real_root[2]; /**&lt; @brief list of the roots (real parts) */
      //T p_im_root[2]; // imaginary part of the roots
      size_t n_real_root_num; /**&lt; @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 &lt; 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 &lt; 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 &lt; -f_thresh) // only really negative
              n_real_root_num = 0; // only two complex roots
          else if(/*fabs(f_disc) &lt; f_thresh*/f_disc &lt;= 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 &gt; 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 &gt;= 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 &lt; 2 &amp;&amp; n_index &lt; 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;
      }
  };
  </code></pre>
  
  <p>The code is horrible, and I hate all the thresholds. But for random equations with roots in the <code>[-100, 100]</code> interval, this is not so bad:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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.</p>
  
  <p>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.</p>
  
  <p>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 <code>plrt()</code>.</p>
  
  <p>I have tried using GNU scientific library <code>gsl_poly_solve_quadratic()</code> but that seems to be a naive approach, and not very numerically stable.</p>
  
  <p>Using <code>std::complex</code> 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).</p>
  
  <p>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).</p>