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

void modregfalsi(double, double, double, int); // function prototype
double f(double); // function prototype

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;

  modregfalsi(a, b, epsilon, imax);

  return 0;
}

// A modified regula falsi function that finds roots of a function
// The maximum number of iterations permitted is imax. The convergence
// criterion is the fractional size of the search interval (x3 - x1) / (b - a)
// is less than epsilon. A relaxation factor RELAX is used
void modregfalsi(double a, double b, double epsilon, int imax)
{
  const double RELAX = 0.9;  // the relaxation factor

  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 
       << "\nThe convergence criterion is: interval < "  << epsilon 
       << "\nThe maximum number of iterations allowed is " << imax << endl;

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

  // iterations
  for (i = 1; i <= imax; i++)
  {
    curwidth = (x3 - x1) / width;
    x2 = x1 - width * curwidth * f1 / (f3 - f1);
    f2 = f(x2);
    if (fabs(curwidth) < epsilon) // root is found
    {
      cout << "\nA root at x = " << x2 << " was found " 
           << "in " << i << " iterations" << endl;
      cout << "The value of the function is " << f2 << endl;
      return;
    }
    else  // check for left and right crossing
    {
      if(f1 * f2 < 0.0) // check for crossing on the left
      {
	 x3 = x2;
	 f3 = f2;
	 f1 = RELAX * f1;
      }
      else if (f2 * f3 < 0.0)  // check for crossing on the right
      {
	x1 = x2;
	f1 = f2;
	f3 = RELAX * f3;
      }
      else  // no crossing in the interval
      {
	cout << "The search for a root has failed due to no root in the interval\n"
	     << "In step " << i << " of the iteration the function does not change sign"              << endl;
      }	
    }
  }
  cout << "\nAfter " << imax << " iterations, no root was found "
       << "within the convergence criterion\n" 
       << "The search for a root has failed due to excessive iterations\n"
       << "after the maximum number of " << imax << " iterations" << endl;
		
  return;
}

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

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