# include <cmath>
# include <iomanip>
# include <iostream>

using namespace std;

# include "bisection.hpp"

//****************************************************************************80

void bisection ( double &a, double &b, double tol, double f ( double x ),
  int &it )

//****************************************************************************80
//
//  Purpose:
//
//    bisection() carries out the bisection method to seek a root of F(X) = 0.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    06 December 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double &a, &b: two points at which the function differs in sign.  
//
//    double tol: the interval size tolerance.  Once |B-A| < TOL,
//    the iteration will stop.
//
//    double f(double x): the name of the routine that evaluates the function.
//
//  Output:
//
//    double &a, &b,: the new endpoints that constitute an
//    change of sign interval no larger than TOL.
//
//    int &it: the number of bisections.
//
{
  double c;
  double fa;
  double fb;
  double fc;

  it = 0;

  fa = f ( a );
  if ( fa == 0.0 )
  {
    b = a;
    return;
  }

  fb = f ( b );
  if ( fb == 0.0 )
  {
    a = b;
    return;
  }

  if ( 0.0 < fa * fb )
  {
    printf ( "bisection(): [A,B] not a change of sign interval!\n" );
    exit ( 1 );
  }

  while ( tol < fabs ( b - a ) )
  {
    c = ( a + b ) / 2.0;
    fc = f ( c );
    it = it + 1;

    if ( 100 < it )
    {
      printf ( "bisection(): Too many iterations!\n" );
      exit ( 1 );
    }

    if ( fc == 0.0 )
    {
      a = c;
      b = c;
      return;
    }
    else if ( 0.0 < fc * fa )
    {
      a = c;
      fa = fc;
    }
    else
    {
      b = c;
      fb = fc;
    }

  }

  return;
}