#define TNNearlyZero  1E-20

double UserFunction(double X,double (*FuncAddr)(double));

double _abs(double x);

void TestInput(
			   double   LeftEndpoint,
			   double	RightEndpoint,
			   double	Tol,
			   int		MaxIter,
			   double	*Answer,
			   double	*fAnswer,
			   int      *Error,
			   int 	    *Found,
			   double (*FuncPtr)(double))

/*--------------------------------------------------------------*/
/*- Input:  LeftEndpoint, RightEndpoint, Tol, MaxIter          -*/
/*- Output: Answer, fAnswer, Error, Found                      -*/
/*-                                                            -*/
/*- This function tests the input data for errors.  If         -*/
/*- LeftEndpoint > RightEndpoint, Tol <= 0, or MaxIter < 0,    -*/
/*- ) an error is returned.  If one the of the endpoints       -*/
/*- (LeftEndpoint, RightEndpoint) is a root, ) Found = TRUE    -*/
/*- and Answer and fAnswer are returned.                       -*/
/*--------------------------------------------------------------*/


{

  double  yLeft, yRight; /* The values of the function at the endpoints. */

  yLeft = UserFunction(LeftEndpoint, FuncPtr);
  yRight = UserFunction(RightEndpoint, FuncPtr);

  if ( _abs(yLeft) <= TNNearlyZero )
  {
	*Answer = LeftEndpoint;
	*fAnswer = UserFunction(*Answer, FuncPtr);
	*Found = 1;
  }

  if ( _abs(yRight) <= TNNearlyZero )
  {
	*Answer = RightEndpoint;
	*fAnswer = UserFunction(*Answer, FuncPtr);
	*Found = 1;
  }

  if (!(*Found))   /* Test for errors */
  {
	if ( yLeft * yRight > 0 )
	  *Error = 2;
	if ( Tol <= 0 )
	  *Error = 3;
	if ( MaxIter < 0 )
	  *Error = 4;
  }
} /* void Tests */

void Initial_Converge_Bisect(double LeftEndpoint,
			 int *Iter,
			 double *yLeft,
			 double (*FuncPtr)(double))
/* Initialize variables. */
{
  *Iter = 0;
  *yLeft = UserFunction(LeftEndpoint, FuncPtr);
} /* function Initial_Converge_Bisect */

int TestForRoot_Converge_Bisect(double X, double OldX, double Y, double Tol)
/*--------------------------------------------------------------------*/
/*-  These are the stopping criteria.  Four different ones are       -*/
/*-  provided.  If you wish to change the active criteria, simply    -*/
/*-  comment off the current criteria (including the appropriate OR) -*/
/*-  and remove the comment brackets from the criteria (including    -*/
/*-  the appropriate OR) you wish to be active.                      -*/
/*--------------------------------------------------------------------*/
{
  return                       /*---------------------------*/
	(_abs(Y) <= TNNearlyZero)          /*- Y = 0                   -*/
									  /*-                         -*/
		   ||                         /*-                         -*/
									  /*-                         -*/
	(_abs(X - OldX) < _abs(OldX * Tol)) /*- Relative change in X    -*/
									  /*-                         -*/
									  /*-                         -*/
/*         ||                      */ /*-                         -*/
/*                                 */ /*-                         -*/
/*  (_abs(X - OldX) < Tol)          */ /*- Absolute change in X    -*/
/*                                 */ /*-                         -*/
/*         ||                      */ /*-                         -*/
/*                                 */ /*-                         -*/
/*  (_abs(Y) <= Tol)                */ /*- Absolute change in Y    -*/
		   ;	 				      /*---------------------------*/

/*-----------------------------------------------------------------------*/
/*- The first criteria simply checks to see if ( the value of the         -*/
/*- function is zero.  You should probably always keep this criteria    -*/
/*- active.                                                             -*/
/*-                                                                     -*/
/*- The second criteria checks the relative error in x. This criteria   -*/
/*- evaluates the fractional change in x between interations. Note      -*/
/*- that x has been multiplied through the inequality to avoid divide   -*/
/*- by zero errors.                                                     -*/
/*-                                                                     -*/
/*- The third criteria checks the absolute difference in x between      -*/
/*- iterations.                                                         -*/
/*-                                                                     -*/
/*- The fourth criteria checks the absolute difference between          -*/
/*- the value of the function and zero.                                 -*/
/*-----------------------------------------------------------------------*/

} /* function TestForRoot_Converge_Bisect */



void Converge_Bisect(double   *LeftEndpoint,
			  double   *RightEndpoint,
			  double   Tol,
			  int	   *Found,
			  int	   MaxIter,
			  double   *Answer,
			  double   *fAnswer,
			  int	   *Iter,
			  int	   *Error,
			  double (*FuncPtr)(double))

/*------------------------------------------------------------------*/
/*- Input:  LeftEndpoint, RightEndpoint, Tol, MaxIter              -*/
/*- Output: Found, Answer, fAnswer, Iter, Error                    -*/
/*-                                                                -*/
/*- This void applies the bisection method to find a          -*/
/*- root to TNTargetF(x) on the interval [LeftEndpoint,            -*/
/*- RightEndpoint]. The root must be found within MaxIter          -*/
/*- iterations to a tolerance of Tol.  If a root is found, ) it -*/
/*- is returned in Answer, the value of the function at the        -*/
/*- approximated root is returned in fAnswer (should be close to   -*/
/*- zero), and the number of iterations is returned in Iter.       -*/
/*------------------------------------------------------------------*/

{
  double yLeft, MidPoint, yMidPoint;


  /* void Converge_Bisect */
  Initial_Converge_Bisect(*LeftEndpoint, Iter, &yLeft, FuncPtr);
  while ( !(*Found) && (*Iter < MaxIter) )
  {
	(*Iter)++;
	MidPoint = (*LeftEndpoint + *RightEndpoint) / 2;
	yMidPoint = UserFunction(MidPoint, FuncPtr);
	*Found = TestForRoot_Converge_Bisect(MidPoint, *LeftEndpoint, yMidPoint, Tol);

	if ( (yLeft * yMidPoint) < 0 )
	   *RightEndpoint = MidPoint;
	else
	  {
		*LeftEndpoint = MidPoint;
		yLeft = yMidPoint;
	  }
  }
  *Answer = MidPoint;
  *fAnswer = yMidPoint;
  if ( *Iter >= MaxIter )
	*Error = 1;
} /* function Converge_Bisect */

double UserFunction(double X, double (*ProcAddr)(double))
{
 return ProcAddr(X);
}

double _abs(double x)
{
 return ( x < 0 ? -x : x);
}

void Bisect(  double LeftEnd,
			  double RightEnd,
			  double Tol,
			  int	 MaxIter,
			  double *Answer,
			  double *fAnswer,
			  int	 *Iter,
			  int    *Error,
			  double (*FuncPtr)(double))

{
  int  Found;

  /* function Bisect */
  Found = 0;
  TestInput(LeftEnd, RightEnd, Tol, MaxIter,
			Answer, fAnswer, Error, &Found, FuncPtr);
  if ( !(*Error) && !Found )  /* i.e. no error */
	Converge_Bisect(&LeftEnd, &RightEnd, Tol, &Found, MaxIter,
			 Answer, fAnswer, Iter, Error, FuncPtr);
} /* function Bisect */
