/* Cubic equation solution. Real coefficients case.

   int Cubic(double *x,double a,double b,double c);
   Parameters:
   x - solution array (size 3). On output:
       3 real roots -> then x is filled with them;
       1 real + 2 complex -> x[0] is real, x[1] is real part of 
                             complex roots, x[2] - non-negative 
                             imaginary part.
   a, b, c - coefficients: x^3 + ax^2 + bx + c = 0.
   Returns: 3 - 3 real roots;
            1 - 1 real root + 2 complex;
            2 - 1 real root + complex roots imaginary part is zero 
                (i.e. 2 real roots). 
*/

#include <math.h>   /* for sqrt(), fabs(), pow(), cos(), acos(). */
#define M_PI (3.141592653589793)
#define M_2PI (2.*M_PI)

int Cubic(double *x,double a,double b,double c) {
 double q,r,r2,q3;
 q=(a*a-3.*b)/9.; r=(a*(2.*a*a-9.*b)+27.*c)/54.;
 r2=r*r; q3=q*q*q;
 if(r2<q3) {
  double t=acos(r/sqrt(q3));
  a/=3.; q=-2.*sqrt(q);
  x[0]=q*cos(t/3.)-a;
  x[1]=q*cos((t+M_2PI)/3.)-a;
  x[2]=q*cos((t-M_2PI)/3.)-a;
  return(3);
 }
 else {
  double aa,bb;
  if(r<0.) r=-r;
  aa=-pow(r+sqrt(r2-q3),1./3.);
  if(aa!=0.) bb=q/aa;
  else bb=0.;
  a/=3.; q=aa+bb; r=aa-bb;
  x[0]=q-a;
  x[1]=-0.5*q-a;
  x[2]=(sqrt(3.)*0.5)*fabs(r);
  if(x[2]==0.) return(2);
  return(1);
 }
}