The correct way to solve quadratic equations is none of the above. For the correct algorithm as applied in the following function, see
numerical recipes in C and the Wolfram site at
http://mathworld.wolfram.com/QuadraticEquation.html
#include <limits>
bool quadratic(double a, double b, double c, double& x1, double& x2, double y, double z)
{
double delta=0.0e+00, q=0.0e+00, bsign, y1, y2;
const double eps = std::numeric_limits()<double>.epsilon();
if ( fabs(a) <= eps ) { return false; }
delta=b*b-4.0*a*c;
if ( delta<0.0e+00 ) { return false; }
if ( delta>=0.0e+00 && delta <= eps )
{ x1=(-1.0*b/(2.0*a)); x2=x1; return true; }
bsign=1.0e+00;
if ( fabs(b)>1.0e-16 ) { bsign=b/fabs(b); }
else { if ( b<0.0e+00 ) { bsign=-1.0e+00; } }
q=(-0.5)*(b+bsign*sqrt(delta));
y1=q/a; y2=c/q;
if (( y<=y1 )&&( y1<=z ))
{
x1=y1; x2=y2;
}
else if (( y<=y2 )&&( y2<=z ))
{
x1=y2; x2=y1;
}
else
{
x1=y1; x2=y2;
}
return true;
}