function [ a, m, b, it, nf ] = bisection_min ( f, a_init, b_init, x_tol )

%*****************************************************************************80
%
%% 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:
%
%    01 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.
%
  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 )

      fa = f ( a );
      fm = f ( m );
      fb = f ( b );
      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;
        fb = fm;
        m = d;
        fm = fd;
%
%  Right midpoint is less than or equal to fm.
%  Shrink interval to [m,right midpoint,b].
%
      elseif ( fe <= fm )
        a = m;
        fa = fm;
        m = e;
        fm = fe;
%
%  fm is strictly smaller than right or left midpoints.
%  Shrink interval to [left midpoint,m,right midpoint].
%
      else
        a = d;
        fa = fd;
        b = e;
        fb = fe;
      end

    end

    it = it + 1;

    if ( abs ( b - a ) < x_tol )
      break
    end

  end

  return
end

