function x = golden_search ( f, a, b, xtol )

%*****************************************************************************80
%
%% golden_search seeks a minimizer of a function in an interval [a,b].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function handle f: the function for which a minimizer bracket
%    is sought.
%
%    real a, b: two distinct points which define the initial search interval.
%
%    real xtol: a positive tolerance.  When the search interval is smaller
%    than xtol, the search will halt.
%
%  Output:
%
%    real x: a point such that a < x < b, so that f(x) is estimated to
%    be the minimum value of f() in the interval.
%
  print_stuff = true;

  phi = ( sqrt ( 5.0 ) + 1.0 ) / 2.0;

  it = 0;
  dx = abs ( b - a );

  while ( true )

    it = it + 1;
%
%  Compute C and D so that:
%  A < C < D < B
%
    c = b - ( b - a ) / phi;
    d = a + ( b - a ) / phi;

    if ( xtol < abs ( c - d ) <= xtol )
      break
    end
%
%  If F(C) < F(D), new interval is [A,C,D]
%  Otherwise, new interval is [C,D,B].
%
    dxold = dx;
    if ( f ( c ) < f ( d ) )
      b = d;
      dx = abs ( b - a );
      alpha = dx / dxold;
      if ( print_stuff )
        fprintf ( 1, 'it=%d, f(%20.16g)=%g, dx = %g  alpha = %g\n', it, c, f(c), dx, alpha );
      end
    else
      a = c;
      dx = abs ( b - a );
      alpha = dx / dxold;
      if ( print_stuff )
        fprintf ( 1, 'it=%d, f(%20.16g)=%g, dx = %g  alpha = %g\n', it, d, f(d), dx, alpha );
      end
    end

  end

  x = ( b + a ) / 2.0;

  return
end

