include/crbn/solver/cubic.hpp

#ifndef __cubic_hpp__
#define __cubic_hpp__

#include <crbn/basic/scalar.hpp>

// Solve cubic equation (only real roots are considered)
inline int cubic_r(const float& a, const float& b, const float& c, const float& d, float& x0, float& x1, float& x2) {
  float a0, a1, a2, a3, q, r, d, s, t;

  a3 = a; a2 = b; a1 = c; a0 = d; 

  a2 /= a3;
  a1 /= a3;
  a0 /= a3;

  q = ((3.0 * a1) - (a2 * a2)) / 9.0;
  r = ((9.0 * a1 * a2) - (27.0 * a0) - (2.0 * a2 * a2 * a2)) / 54.0;
  
  d = (q * q * q) + (r * r); 
  
  if(d > 0.0f) {
      d = sqrt(d);

      s = cbrt(r + d);
      t = cbrt(r - d);
      
      x0 = x1 = x2 = -(a2 / 3.0) + s + t;   
      
      return 1;
  }

  if (d < 0.0f) {
    float theta = acos(r / sqrt(-(q * q * q)));
    s = 2.0 * sqrt(-q);
    x0 = (s * cos(theta / 3.0)) - (a2 / 3.0);
    x1 = (s * cos((theta + (2.0 * M_PI)) / 3.0)) - (a2 / 3.0);
    x2 = (s * cos((theta - (4.0 * M_PI)) / 3.0)) - (a2 / 3.0);
    
    if(x0 > x1) {
      t = x0;
      x0 = x1;
      x1 = t;
    }
    if(x0 > x2) {
      t = x0;
      x0 = x2;
      x2 = t;
    }
    if(x1 > x2) {
      t = x1;
      x1 = x2;
      x2 = t;
    }
    
    return 3;
  }

  s = cbrt(r);
  
  x0 = -(a2 / 3.0) + (2.0 * s);   
  x1 = x2 = -(a2 / 3.0) - s;
  
  if(x0 > x1) {
    x2 = x0;
    x0 = x1;
    x1 = x2;
  }
  
  return 2;
}

#endif

---