/* Quadratic equation solution. Real coefficients case.

   int Quadratic(double *x,double a,double b,double c);
   Parameters:
   x - solution array (size 2). On output:
       2 real roots -> then x is filled with them;
       2 complex-conjugate roots -> x[0] is real part, 
         x[1] is non-negative imaginary part. 
       other cases -> x[0] unique valid root if (-1) returned,
                     no valid roots otherwise.
   a, b, c - coefficients: ax^2 + bx + c = 0. 
   Returns: 2 - 2 real and distince roots;
            1 - 1 real distinct root (x[0]=x[1]);
            0 - 2 complex roots;
            -1 - one real root in case a==0;
            -2 - no roots in case a=0, b=0;
            -3 - infinite number of roots (a=b=c=0). 
*/

#include <math.h>   /* for sqrt() */

int Quadratic(double *x,double a,double b,double c) {
 double d;
 /* degenerated cases */
 if(a==0.) {
  if(b==0.) {
   if(c==0.) return(-3);
   return(-2);
  }
  x[0]=-c/b; return(-1);
 }
 /* the main case */
 d=b*b-4.*a*c;       /* the discriminant */
 /* one distinct root */
 if(d==0.) {
  x[0]=x[1]=-b/(2.*a); return(1);
 }
 /* conjugate complex roots */
 if(d<0.) {
  double t=0.5/a;
  x[0]=-b*t; x[1]=sqrt(-d)*t;
  return(0);
 }
 /* 2 real roots: avoid subtraction of 2 close numbers */
 if(b>=0.) d=(-0.5)*(b+sqrt(d));
 else d=(-0.5)*(b-sqrt(d));
 x[0]=d/a; x[1]=c/d;
 return(2);
}