# include <cmath>
# include <cstdlib>

using namespace std;

# include "bisection_min.hpp"

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

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 )

//****************************************************************************80
//
//  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:
//
//    05 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;
}