A simple program to solve quadratic equations with

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();
   // equation is: a*x^2 + b*x + c = 0; a cannot be zero.
   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));
 
   // the roots of the equation are:
   y1=q/a; y2=c/q;
 
   // find if any of the roots is in the given [y,z] interval
   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;
}