/* 
bisection method for root finding.
of function f(x) in closed interval [a,b].
*/

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <gsl/gsl_math.h>

using std::endl;
using std::cout;

// define function:

double func_f(double x){
  return pow(x,3.0) + 1.0;
}


int main() 
{    
   int N= 0;
   int N_limit = 1000;           // max number of trys

   double a = -100.0;   // lower bound for [a,b]
   double b = 3.0;      // upper bound for [a,b]

   double tol = 1.0e-5;          // tolerance for finding zero

   if (func_f(a) == 0.0) {
     cout << a << " is a root." << endl;
     return 0;
   }

   if (func_f(b) == 0.0) {
     cout << b << " is a root." << endl;
     return 0;
   }

   if (GSL_SIGN(func_f(a)) * GSL_SIGN(func_f(b)) > 0.0){
     cout << "[" << a << "," << b << "] do not bracket a root." << endl;
     return 0;
   }
   
   double x_low = a;
   double x_high = b;
   double x_mid;

   while ((fabs(x_high - x_low) >= tol) && (N < N_limit)){
     x_mid = 0.5*(x_low + x_high);
     if ( GSL_SIGN(func_f(x_mid)) * GSL_SIGN(func_f(x_high)) < 0){
       x_low = x_mid;
     }
     else {
       x_high = x_mid;
     }
     ++N;
   }

   if(N > N_limit)
     cout << "iteration limit exceeded" << endl;
   else
     cout << "estimated root = " << x_mid << " with " << N << " iterations." << endl;
   
   return 0;
}