#include <iostream.h>
#include <iomanip.h>
#include <math.h>


//============= Begin bisect() ===================
//
// Solves the single nonlinear equation
//               f(x) = 0
// on the user-specified interval [a,b] by the method of bisection.
// Receives pointer to function f.

double bisect(double (*f)(double), double a, double b, 
              double tol, int &error)
{
   //
   // Solve an equation f(x) = 0 for x using the method of bisection. 
   //
   double f_a, f_b, f_c, c;
   
   f_a = f(a);  f_b = f(b);  c = (a+b)/2.0;       // Call f twice.
   
   if (f_a*f_b > 0.0)
       error = 1;
   else if (b < a)           // Check input data for errors.
   { 
       error = 2;
   }
   else if (b - a <= tol)
       return c;
   else
   {
       while ( b-a > tol )     // Note: only need to call f() ONCE
       {                       //       per loop iteration.
            f_c = f(c);
            if ( f_c*f_a <= 0 )       // [a, c] must contain a root.
            {
                  b   = c;
                  f_b = f_c;
            }
            else                     // [c, b] must contain a root.
            {
                  a   = c;
                  f_a = f_c;
            } 
            c = (a+b)/2;
       }
       return c;
   }
}

//============= End of bisect() ==================