//	Bisect6.cc

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

//	Version 6: All three Bisection functions witzn the same name
//

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

const double EPS = 1e-6;		// accuracy constant

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)

double h(const double x)		// declaration and
  { return 3.0 - exp(x) ; }		//   definition of function h(x)

double t(const double x)		// declaration and
  { return 1-x*x ; }			//   definition of function t(x)
  
//	The three versions of Bisect differ by its signature
					//	declaration of Bisect
double Bisect(const double a, const double b, const double eps);
					// declaration of Bisect2
double Bisect(const double a, const double b);
					// declaration of Bisect3
double Bisect(double (* const func)(double), const double a,  
		const double b, const double eps);

main()
{
 double a,b,x0,fa,fb;
 char   choice;
 double (*ff)(double);
 
 
 cout << endl;
 cout << " Determine point of zero in [a,b] by bisection " << endl;

 cout << "  f(x) :=  sin(x) - x/2" << endl;
 cout << "  g(x) := (1.234567-x)*(x+0.987654)" << endl;
 cout << "  h(x) :=  3.0 - exp(x)" << endl;
 cout << "  t(x) :=  1-x*x" << endl;

 cout << endl << "Which function do you prefer ? "; 
 cin  >> choice;
 
 switch (choice)
 {
   case 'f':  ff = f; break;
   case 'g':  ff = g; break;
   case 't':  ff = t; break;
   default: 
   	choice = 'h';
	cout <<  " no correct choice. h(x) is used." << endl;
   case 'h':  ff = h; break;
 }

 cout << " " << choice << "(a) > 0,  a : "; cin >> a;
 cout << " " << choice << "(b) < 0,  b : "; cin >> b;
 cout << endl;
 
 fa = ff(a);  fb = ff(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 " << choice 
        << "(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 = Bisect(ff,b,a,EPS);		// Bisect(3) is called
    }
   else					// f(a) > 0 && f(b) < 0
    {
     x0 = Bisect(ff,a,b,EPS);		// Bisect(3) is called
    }
   cout << endl << " point of zero = " << x0 << endl;
  }
 cout << endl;
}

//	---------------------------------------------------------------
//	Recursive  function Bisect
//	---------------------------------------------------------------
//	definition of Bisect

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

 return x0;			// return the solution
}

//	---------------------------------------------------------------
//	Recursive  function Bisect2
//	---------------------------------------------------------------

double Bisect(const double a, const double b)	// definition of Bisect2
{
 double x0, fc, c = (a+b)/2;
 
 fc = f(c);
 if ( fabs(fc) < EPS )
  {
   x0 = c;			// end of recursion
  }
 else if (  fc > 0.0 )
  {
   x0 = Bisect(c,b);		// search in right intervall
  }
 else 				// i.e., fc < 0.0
  {
   x0 = Bisect(a,c);		// search in left intervall
  } 

 return x0;			// return with solution
}

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

double Bisect(double (* const 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 = Bisect(func,c,b,eps);	// search in right intervall
  }
 else 			// i.e., fc < 0.0
  {
   x0 = Bisect(func,a,c,eps);	// search in left intervall
  } 

 return x0;			// return with solution
}