title: C++ solving for quartic roots (fourth order polynomial)
body: <p>As part of my project I need to solve a quartic polynomial in a <strong>closed form</strong> in C++.</p>
  
  <blockquote>
    <p>A*x<sup>4</sup> + B*x<sup>3</sup> + C*x<sup>2</sup> + D*x + E = 0</p>
  </blockquote>
  
  <p>I found several links to this end. One of them is <a href="http://www.gamedev.net/topic/451048-best-way-of-solving-a-polynomial-of-the-fourth-degree/" rel="nofollow">here</a>. But it computes all roots, while I want just real roots. The algorithms mainly use Ferrari's method to reduce the order.</p>
  
  <pre><code>bool solveQuartic(double a, double b, double c, double d, double e, double &amp;root)
  {
  // I switched to this method, and it seems to be more numerically stable.
  // http://www.gamedev.n...topic_id=451048 
  
  // When a or (a and b) are magnitudes of order smaller than C,D,E
  // just ignore them entirely. This seems to happen because of numerical
  // inaccuracies of the line-circle algorithm. I wanted a robust solver,
  // so I put the fix here instead of there.
  if(a == 0.0 || abs(a/b) &lt; 1.0e-5 || abs(a/c) &lt; 1.0e-5 || abs(a/d) &lt; 1.0e-5)
      return solveCubic(b, c, d, e, root);
  
  double B = b/a, C = c/a, D = d/a, E = e/a;
  double BB = B*B;
  double I = -3.0*BB*0.125 + C;
  double J = BB*B*0.125 - B*C*0.5 + D;
  double K = -3*BB*BB/256.0 + C*BB/16.0 - B*D*0.25 + E;
  
  double z;
  bool foundRoot2 = false, foundRoot3 = false, foundRoot4 = false, foundRoot5 = false;
  if(solveCubic(1.0, I+I, I*I - 4*K, -(J*J), z))
  {
      double value = z*z*z + z*z*(I+I) + z*(I*I - 4*K) - J*J;
  
      double p = sqrt(z);
      double r = -p;
      double q = (I + z - J/p)*0.5;
      double s = (I + z + J/p)*0.5;
  
      bool foundRoot = false, foundARoot;
      double aRoot;
      foundRoot = solveQuadratic(1.0, p, q, root);
      root -= B/4.0;
  
      foundARoot = solveQuadratic(1.0, r, s, aRoot);
      aRoot -= B/4.0;
      if((foundRoot &amp;&amp; foundARoot &amp;&amp; ((aRoot &lt; root &amp;&amp; aRoot &gt;= 0.0) 
          || root &lt; 0.0)) || (!foundRoot &amp;&amp; foundARoot)) 
      {
          root = aRoot;
          foundRoot = true;
      }
  
      foundARoot = solveQuadraticOther(1.0, p, q, aRoot);
      aRoot -= B/4.0;
      if((foundRoot &amp;&amp; foundARoot &amp;&amp; ((aRoot &lt; root &amp;&amp; aRoot &gt;= 0.0) 
          || root &lt; 0.0)) || (!foundRoot &amp;&amp; foundARoot)) 
      {
          root = aRoot;
          foundRoot = true;
      }
  
      foundARoot = solveQuadraticOther(1.0, r, s, aRoot);
      aRoot -= B/4.0;
      if((foundRoot &amp;&amp; foundARoot &amp;&amp; ((aRoot &lt; root &amp;&amp; aRoot &gt;= 0.0) 
          || root &lt; 0.0)) || (!foundRoot &amp;&amp; foundARoot)) 
      {
          root = aRoot;
          foundRoot = true;
      }
      return foundRoot;
  }
  return false;
  }
  </code></pre>
  
  <p>This uses solveCubic() which gives both real and imaginary solution:</p>
  
  <pre><code>bool solveCubic(double &amp;a, double &amp;b, double &amp;c, double &amp;d, double &amp;root)
  {
  if(a == 0.0 || abs(a/b) &lt; 1.0e-6)
      return solveQuadratic(b, c, d, root);
  
  double B = b/a, C = c/a, D = d/a;
  
  double Q = (B*B - C*3.0)/9.0, QQQ = Q*Q*Q;
  double R = (2.0*B*B*B - 9.0*B*C + 27.0*D)/54.0, RR = R*R;
  
  // 3 real roots
  if(RR&lt;QQQ)
  {
      /* This sqrt and division is safe, since RR &gt;= 0, so QQQ &gt; RR,    */
      /* so QQQ &gt; 0.  The acos is also safe, since RR/QQQ &lt; 1, and      */
      /* thus R/sqrt(QQQ) &lt; 1.                                     */
      double theta = acos(R/sqrt(QQQ));
      /* This sqrt is safe, since QQQ &gt;= 0, and thus Q &gt;= 0             */
      double r1, r2, r3;
      r1 = r2 = r3 = -2.0*sqrt(Q);
      r1 *= cos(theta/3.0);
      r2 *= cos((theta+2*PI)/3.0);
      r3 *= cos((theta-2*PI)/3.0);
  
      r1 -= B/3.0;
      r2 -= B/3.0;
      r3 -= B/3.0; 
  
      root = 1000000.0;
  
      if(r1 &gt;= 0.0) root = r1;
      if(r2 &gt;= 0.0 &amp;&amp; r2 &lt; root) root = r2;
      if(r3 &gt;= 0.0 &amp;&amp; r3 &lt; root) root = r3;
  
      return true;
  }
  // 1 real root
  else
  {
      double A2 = -pow(fabs®+sqrt(RR-QQQ),1.0/3.0);
      if (A2!=0.0) {
          if (R&lt;0.0) A2 = -A2; 
          root = A2 + Q/A2; 
      }
      root -= B/3.0;
      return true;
  }
  }
  </code></pre>
  
  <p>Here are some links that explain the code. <a href="http://www.1728.org/cubic2.htm" rel="nofollow">solveCubic</a> and <a href="http://www.1728.org/quartic2.htm" rel="nofollow">solveQuartic</a></p>
  
  <p>Is there anyone who can modify the code to solve a quartic polynomial for real roots?</p>
  
  <p>I want to implement it as efficiently as possible. BTW I'd appreciate if someone introduced a useful library for this purpose like LAPACK (it seems to be unable to calculate the roots of a quartic polynomial directly).</p>