# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>
//
//  We declare all our functions here, so anyone can use them.
//
double bisect2 ( double a, double b, double f ( double x ) );
double f1 ( double x );
double f2 ( double x );
double f3 ( double x );

using namespace std;

int main ( )
{
  double a;
  double b;
  double c;
//
//  Function #1.
//
  a = - 10.0;
  b = + 10.0;
  c = bisect2 ( a, b, f1 );

  cout << "\n";
  cout << "BISECT2 returned solution estimate C = " << c << "\n";
  cout << "F1(C) = " << f1 ( c ) << "\n";
//
//  Function #2.
//
  a = 0.0;
  b = 3.0;
  c = bisect2 ( a, b, f2 );

  cout << "\n";
  cout << "BISECT2 returned solution estimate C = " << c << "\n";
  cout << "F2(C) = " << f2 ( c ) << "\n";
//
//  Function #3.
//
  a = 2.0;
  b = 3.0;
  c = bisect2 ( a, b, f3 );

  cout << "\n";
  cout << "BISECT2 returned solution estimate C = " << c << "\n";
  cout << "F3(C) = " << f3 ( c ) << "\n";
  return 0;
}
double bisect2 ( double a, double b, double f ( double x ) )
//
//  BISECT2 is a bisection solver.  The name of the function is an input argument.
//  BISECT2 thinks the function is called "F", but you can call it anything,
//  as long as it takes a double input and returns a double output.
//
{
  double c;
//
//  Make sure [A,B] is a change of sign interval.
//
  if ( 0 < f ( a ) * f ( b ) )
  {
    cerr << "\n";
    cerr << "BISECT2 - Fatal error!\n";
    cerr << "  F(A) = " << f(a) << "\n";
    cerr << "  F(B) = " << f(b) << "\n";
    cerr << "  [" << a << ", " << b << "] is not a change of sign interval.\n";
  }

  while ( true )
  {
//
//  C = midpoint.
//
    c = ( a + b ) / 2.0;
//
//  Is F(C) extremely close to 0?
//
    if ( fabs ( f ( c ) ) <= 0.000001 )
    {
      break;
    }
//
//  If F(C) is opposite in sign to F(A), C replaces B.
//
    if ( f ( c ) * f ( a ) < 0.0 )
    {
      b = c;
    }
    else
    {
      a = c;
    }
  }
  return c;
}
double f1 ( double x )
//
//  F1 is our "cosine" function.
{
  double value;

  value = cos ( x ) - x;

  return value;
}
double f2 ( double x )
//
//  F2 is our "satellite" function.
//
{
  double value;

  value = 1.082 * cos ( x ) + sin ( x );

  return value;
}
double f3 ( double x ) 
//
//  F3 is our "pole" function.
//
{
  double value;

  value = pow ( x, 4 ) + 2 * pow ( x, 3 ) - 14 * pow ( x, 2 ) + 2 * x + 1;

  return value;
}