double bisect1 ( double a, double b )
//
//  BISECT1 is a program to carry out bisection to solve F(X) = 0.
//  The function must be declared as:
//
//    double f ( double x );
//
//  Because we call FABS() the program needs the "#include <cmath>" statement.
//
{
  double c;

  while ( true )
  {
//
//  C = midpoint.
//
    c = ( a + b ) / 2.0;
//
//  Is F(C) extremely close to 0?
//
    if ( fabs ( f ( c ) ) <= 0.000001 )
    {
      break;
    }
//
//  If F(C) is opposite in sign to F(A), C replaces B.
//
    if ( f ( c ) * f ( a ) < 0.0 )
    {
      b = c;
    }
    else
    {
      a = c;
    }
  }
  return c;
}