# include <float.h>
# include <math.h>
# include <stdio.h>
# include <stdlib.h>

# include "bisection.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    05 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;

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

  *it = 0;

  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 ( 0.0 < fc * fa )
    {
      *a = c;
      fa = fc;
    }
    else
    {
      *b = c;
      fb = fc;
    }

  }

  return;
}