subroutine bisection_min ( f, a_init, b_init, x_tol, a, m, b, it, nf )

!*****************************************************************************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:
!
!    02 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.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) a
  real ( kind = rk8 ) a_init
  real ( kind = rk8 ) b
  real ( kind = rk8 ) b_init
  real ( kind = rk8 ) d
  real ( kind = rk8 ) e
  real ( kind = rk8 ), external :: f
  real ( kind = rk8 ) fa
  real ( kind = rk8 ) fb
  real ( kind = rk8 ) fd
  real ( kind = rk8 ) fe
  real ( kind = rk8 ) fm
  integer it
  real ( kind = rk8 ) m
  integer nf
  real ( kind = rk8 ) x_tol

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

  do 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 ) then

      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.5D+00 * ( a + m )
      fd = f ( d )

      e = 0.5D+00 * ( 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 ) then
        b = m
        fb = fm
        m = d
        fm = fd
!
!  Right midpoint is less than or equal to fm.
!  Shrink interval to [m,right midpoint,b].
!
      else if ( fe <= fm ) then
        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 if

    end if

    it = it + 1

    if ( abs ( b - a ) < x_tol ) then
      exit
    end if

  end do

  return
end