#include <iostream>
#include <cmath>
using namespace std;


void bisection(double, double, double, int);
double f(double);


int main()
{
  int imax;          // maximum number of iterations
  double a, b;       // left and right ends of the original interval
  double epsilon;    // convergence criterion
	
  // obtain the input data
  cout << "Enter the limits of the original search interval, a and b: ";
  cin  >> a >> b;
  cout << "Enter the convergence criteria: ";
  cin  >> epsilon;
  cout << "Enter the maximum number of iterations allowed: ";
  cin  >> imax;

  bisection(a, b, epsilon, imax);

  return 0;
}

// A bisection function that finds roots of a function
// The interval a < x < b is known to contain a root of f(x). The estimate
// of the root is successively improved by finding in which half of the interval
// the root lies and then replacing the original interval by that half-interval.
//
void bisection(double a, double b, double epsilon, int imax)
{
  int i;             // current iteration counter
  double x1, x2, x3; // left, right, and midpoint of current interval
  double f1, f2, f3; // function evaluated at these points
  double width;      // width of original interval = (b - a)
  double curwidth;   // width of current interval = (x3 - x1)
	
  // echo back the passed input data
  cout << "\nThe original search interval is from " << a << " to " << b << endl;
  cout << "The convergence criterion is: interval < "  << epsilon << endl;
  cout << "The maximum number of iterations allowed is " << imax << endl;

  // calculate the root
  x1 = a;
  x3 = b;
  f1 = f(x1);
  f3 = f(x3);
  width = (b - a);

  // verify there is a root in the interval
  if (f1 * f3 > 0.0)
    cout << "\nNo root in the original interval exists" << endl;
  else
  {
    for (i = 1; i <= imax; i++)
    {
      // find which half of the interval contains the root
      x2 = (x1 + x3) / 2.0;
      f2 = f(x2);
      if (f1 * f2 <= 0.0)  // root is in left half interval
      {
        curwidth = (x2 - x1) / 2.0;
        f3 = f2;
        x3 = x2;
      }
      else  // root is in right half interval
      {
        curwidth = (x3 - x2) / 2.0;
        f1 = f2;
        x1 = x2;
      }
      if (curwidth < epsilon)
      {
   	cout << "\nA root at x = " << x2 << " was found " 
    	     << "in " << i << " iterations" << endl;
	cout << "The value of the function is " << f2 << endl;
	return;
      }
    }
  }
  cout << "\nAfter " << imax << " iterations, no root was found "
       << "within the convergence criterion" << endl;

  return;
}

// function to evaluate f(x)
double f(double x)
{
	const double PI = 3.14;

	return (exp(-x) - sin(0.5 * PI * x));
}