# include <math.h>
# include <stdbool.h>
# include <stdlib.h>

# include "bisection_min.h"

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

void bisection_min ( double f ( double x ), double a_init, double b_init, 
  double x_tol, double *a, double *m, double *b, int *it, int *nf )

/******************************************************************************/
/*
  Purpose:

    bisection_min() uses bisection to find the minimum of a unimodal function.

  Discussion:

    The function f(x) is assumed to be strictly unimodal:
    to the left of the single minimizer x* in [a,b], f(x) is strictly 
    decreasing, and to the right it is strictly increasing.

    If f(x) is only weakly unimodal, allowing intervals of constant value,
    this procedure is liable to failure.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    04 March 2026

  Author:

    John Burkardt

  Input:

    real f(x): evaluates the unimodal function.

    real a_init, b_init: the initial interval.

    real x_tol: bisection stops when the interval is smaller than x_tol.

  Output:

    real a, m, b: the left, middle, and right points of the final interval.

    integer it: the number of iterations.

    integer nf: the number of function evaluations.
*/
{
  double d;
  double e;
  double fd;
  double fe;
  double fm;

  *a = a_init;
  *m = 0.5 * ( a_init + b_init );
  *b = b_init;
  *it = 0;
  *nf = 0;

  while ( true )
  {
/*
  Start the iteration by evaluating f at the midpoint.
  By strict unimodality, we must have fm < fa and fm < fb.
*/
    if ( it == 0 )
    {
      fm = f ( *m );
      *nf = *nf + 3;
    }
/*
  Evaluate f at the left and right midpoints
  and consider how to proceed.
*/
    else
    {
      d = 0.5 * ( *a + *m );
      fd = f ( d );

      e = 0.5 * ( *m + *b );
      fe = f ( e );

      *nf = *nf + 2;
/*
  Left midpoint is less than or equal to fm.
  Shrink interval to [a, left midpoint, m].
*/
      if ( fd <= fm )
      {
        *b = *m;
        *m = d;
        fm = fd;
      }
/*
  Right midpoint is less than or equal to fm.
  Shrink interval to [m,right midpoint,b].
*/
      else if ( fe <= fm )
      {
        *a = *m;
        *m = e;
        fm = fe;
      }
/*
  fm is strictly smaller than right or left midpoints.
  Shrink interval to [left midpoint,m,right midpoint].
*/
      else
      {
        *a = d;
        *b = e;
      }
    }

    *it = *it + 1;

    if ( fabs ( *b - *a ) < x_tol )
    {
      break;
    }
  }

  return;
}