,sxrjd / 2 / sqeq

#include <math.h>

// this function returns:
// 0 - not equation (a=0, b=0)
// 1 - linear           x1=x2
// 2 - real equal roots x1=x2
// 3 - real two roots   x1,x2
// 4 - complex real x1, imaginary x2
//
int sqeq(double a, double b, double c, double *x1, double *x2){
    if( a == 0 )
      if( b == 0 ) return 0;
      else{
          *x2 = *x1 = -c/b;
          return 1;
      }

double disc = b*b - 4*a*c;
    if( disc == 0 ){
        *x2 = *x1 = -b/(2*a);
        return 2;
    }
    if( disc > 0 ){
        *x1 = (-b - sqrt(disc))/(2*a);
        *x2 = (-b + sqrt(disc))/(2*a);
        return 3;
    }
    *x1 = -b/(2*a);
    *x2 = sqrt(-disc)/(2*a);
    return 4;
}