A simple program to solve quadratic equations with

#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;
}