Listing 1   The worldÕs best root finder

double iterate(double x0, double x2, double (*f)(double), 
double eps, int imax, nt & error)
{
  double x1,y0,y1,y2,b,c,temp,y10,y20,y21,xm,ym;
  double xmlast = x0;
  error = NONE;
  y0 = f(x0);
  if (y0 ==0.0) return x0;
  
  y2 = f(x2);
  if (y2 ==0.0) return x2;

  if(y2 * y0 > 0.0)
  {
    error = BAD_DATA;
    return x0;
  }

  for (int i = 0; i < imax; i++)
  {
    x1 = 0.5 * (x2 + x0);
    y1 = f(x1);
    if (y1 == 0.0) return x1;
    if (abs (x1 - x0) < eps) return x1;
    if (y1 * y0 > 0.0)
    {
      temp = x0;
      x0 = x2;
      x2 = temp;
      temp = y0;
      y0 = y2;
      y2 = temp;
    }
    y10 = y1 - y0;
    y21 = y2 - y1;
    y20 = y2 - y0;
    if (y2 * y20 <2.0 * y1 * y10)
    {
      x2 = x1;
      y2 = y1;
    }
    else
    {
      b = (x1 - x0) / y10;   
      c = (y10 -y21) / (y21 * y20); 
      xm = x0 - b * y0 * (1.0-c * y1);
      ym = f(xm);
      if (ym == 0.0) return xm;
      if (abs (xm - xmlast) < eps)
      {
        return xm;
      }
      xmlast = xm;
      if (ym * y0 < 0.0)
      {
        x2 =x m;
        y2 = ym;
      }
      else
      {
        x0 = xm;
        y0 = ym;
        x2 = x1;
        y2 = y1;
      }
    }
  }
  error = NO_CONVERGENCE;
  return x1;
}