FUNCTION INSERT

/* Insert a new point adjacent to the current minimum.
 * the point may lie to the left or right of x[index].
 * Return TRUE if the new Y is a new minimum.
 */
int insert(double X, double Y){
   assert(npoints > 2);
   assert((x[index-1]< X) && (X < x[index+1]));
   assert((0 < index) && (index < npoints-1));
   assert(X != x[index]);
   if(X > x[index])
		return insert_right(X, Y);
   else
 		return insert_left(X, Y);
}


BISECTION

// bisect the interval between elements k-1 and k
STATE bisect_left(double (*f)(double), int k){
	assert((k > 0) && (k < npoints-1) && (npoints > 2));
   double X = (x[k] + x[k-1])/2;
   double Y = f(X);
   insert_left(X, Y);
   return PARA_FIT;
}

// bisect the interval between elements k and k+1
STATE bisect_right(double (*f)(double), int k){
   assert((k > 0) && (k < npoints-1) && (npoints > 2));
   double X = (x[k] + x[k+1])/2;
   double Y = f(X);
   insert_right(X, Y);
   return PARA_FIT;
}


// bisect the largest interval adjacent to the current minimum
STATE bisect(double (*f)(double)){
   assert((npoints > 2) && (index > 0) && (index < npoints-1));
   double dx0
    = x[index] - x[index-1];
   double dx1 = x[index+1] - x[index];
   if(dx0 > dx1)
		return bisect_left(f, index);
   else
      return bisect_right(f, index);
}


 THE INTERPOLATION ROUTINES

/* do a parabolic fit of the three points centered at index
 * If the fit produces a new minimum, repeat.
 * Else do a bisection.
 */
STATE para_fit(double (*f)(double), double &dx){
   #define MAX_PARA 1
   assert(npoints >= 3);
   assert((index > 0) && (index < npoints-1));

   static int para_count = MAX_PARA;
   double X, Y, Yguess;
   double dx0 = x[index] - x[index-1];
   double dx1 = x[index+1] - x[index];
   double span = dx0 + dx1;
	double m0 = (y[index] - y[index-1])/dx0;
	double m1 = (y[index+1] - y[index])/dx1;
	double m  = (y[index+1] - y[index-1])/span;
	double c  = (m1 - m0)/span;
   double oldy = y[index];
   // compute the estimated minimum
	X = (x[index-1] + x[index+1] - m/c)/2.0;
   dx = X - x[index];
   Yguess = y[index-1] - c*X*X;
   Y = f(X);


   // make sure results are in range
   // don't insert if not
   if((X < x[index-1]) || (X > x[index+1]))
      return BISECT;
   if(!insert(X, Y))
      return BISECT;
   if(!(--para_count)){
      para_count = MAX_PARA;
      return BISECT;
   }
   if((Yguess < oldy) && (Y > oldy))
      return FLIP;
   else
      return PARA_FIT;
}

// Flip the step to the opposite side of index
STATE flip(double (*f)(double), double dx){
   double X = x[index] - 2*dx;
   double Y = f(X);
   insert(X, Y);
   return PARA_FIT;
}


// do a "wide" parabolic fit of the three points x0, x2, and x4
void para_wide(double (*f)(double)){
   assert(npoints ==5);
   double X, Yguess, Y;

	double m0 = (y[2] - y[0])/(x[2] - x[0]);
	double m1 = (y[4] - y[2])/(x[4] - x[2]);
	double m  = (y[4] - y[0])/(x[4] - x[0]);
	double c  = (m1 - m0)/(x[4] - x[0]);

	X = (x[0] + x[4] - m/c)/2.0;
   Yguess = y[0] - c*X*X;
   Y = f(X);
   insert(X, Y);
}


// Fit straight lines through the outer two pairs of points,
// and find the intersection.
// return the results in X and Y
void line_fit(double (*f)(double)){
   double X, Yguess, Y;
   assert(npoints ==5);

	double m0 = (y[1] - y[0])/(x[1] - x[0]);
	double m3 = (y[4] - y[3])/(x[4] - x[3]);
   cout << m0 << ' ' << m3 << endl;
	X = (m3*x[3] + m0*x[0] - (y[3]-y[0]))/(m3-m0);
   Yguess = y[0] + m0*(X-x[0]);
   Y = f(X);
   insert(X, Y);
}