//	Bisect4.cc

//	Recursive function
//			Example: Find the point of zero by bisection

//	Version 4: func(x) and epsilon as parameters
//		   more robust wrt. intervall bounds


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

double f(const double x)		// declaration and
  { return sin(x) - 0.5*x ; }		//   definition of function f(x)

double g(const double x)		// declaration and
  { return -(x-1.234567)*(x+0.987654) ; }//   definition of function g(x)

					// declaration of Bisect3
double Bisect3(double (*func)(double), const double a, const double b, 
		const double eps);

main()
{
 const double EPS = 1e-6;		// accuracy constant
 double a,b,x0,fa,fb;
 
 
 cout << endl;
 cout << " Determine point of zero in [a,b] by bisection " << endl;

 cout << "f(a) > 0,  a : "; cin >> a;
 cout << "f(b) < 0,  b : "; cin >> b;
 cout << endl;
 
 fa = f(a);  fb = f(b);
 
 if ( fabs(fa)<EPS || fabs(fb)<EPS)	// solution is a or/and b
  {
   if ( fabs(fa)<EPS )			// solution is a
    {
      x0 = a; 
      cout << endl << " point of zero = " << x0 << endl;
    }
   if ( fabs(fb)<EPS )			// solution is b
    {
      x0 = b; 
      cout << endl << " point of zero = " << x0 << endl;
    }
  }
 else if ( fa*fb > 0.0)			// no solution in [a,b]
  {
   cout << endl;
   cout << "There is potentially no solution for " 
        << "f(x) = 0  in [" << a << "," << b << "]" << endl;
  }
 else
  {
   if ( fa < 0.0 ) 			// f(a) < 0 && f(b) > 0
    {
     cout << "I have to swap a and b" << endl;
     x0 = Bisect3(f,b,a,EPS);		// call recursive function
    }
   else					// f(a) > 0 && f(b) < 0
    {
     x0 = Bisect3(f,a,b,EPS);		// call recursive function
    }
   cout << endl << " point of zero = " << x0 << endl;
  }
 cout << endl;
}

//	---------------------------------------------------------------
//	Recursive  function Bisect3
//	---------------------------------------------------------------
//	definition of Bisect3

double Bisect3(double (*func)(double), const double a, const double b, 
		const double eps)
{
 double x0, fc, c = (a+b)/2;
 
 fc = func(c);
 if ( fabs(fc) < eps )
  {
   x0 = c;			// end of recursion
  }
 else if (  fc > 0.0 )
  {
   x0 = Bisect3(func,c,b,eps);	// search in right intervall
  }
 else 			// i.e., fc < 0.0
  {
   x0 = Bisect3(func,a,c,eps);	// search in left intervall
  } 

 return x0;			// return with solution
}