function [ x, fx, calls ] = glomin ( a, b, c, m, e, t, f )

%*****************************************************************************80
%
%% glomin() seeks a global minimum of a function F(X) in an interval [A,B].
%
%  Discussion:
%
%    This function assumes: 
%    * F(X) is twice continuously differentiable over [A,B];
%    * F''(X) <= M for all X in [A,B];
%    * the user can supply the value of this upper bound M.
%
%    Thanks to Hans Bieshaar for supplying several corrections to the code,
%    28 May 2021
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    28 May 2021
%
%  Author:
%
%    Original FORTRAN77 version by Richard Brent.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization Without Derivatives,
%    Dover, 2002,
%    ISBN: 0-486-41998-3,
%    LC: QA402.5.B74.
%
%  Input:
%
%    real A, B, the endpoints of the interval.
%    It must be the case that A < B.
%
%    real C, an initial guess for the global
%    minimizer.  If no good guess is known, C = A or B is acceptable.
%
%    real M, the bound on the second derivative.
%
%    real E, a positive tolerance, a bound for the
%    absolute error in the evaluation of F(X) for any X in [A,B].
%
%    real T, a positive error tolerance.
%
%    function value = F ( x ), the name of a user-supplied
%    function whose global minimum is being sought.
%
%  Output:
%
%    real X, the estimated value of the abscissa
%    for which F attains its global minimum value in [A,B].
%
%    real FX, the value F(X).
%
%    integer CALLS: the number of function calls.
%
  calls = 0;
  a0 = b;
  x = a0;
  a2 = a;
  y0 = f ( b );
  calls = calls + 1;
  yb = y0;
  y2 = f ( a );
  calls = calls + 1;
  y = y2;

  if ( y0 < y )
    y = y0;
  else
    x = a;
  end

  if ( m <= 0.0 || b <= a )
    fx = y;
    return
  end

  m2 = 0.5 * ( 1.0 + 16.0 * eps ) * m;

  if ( c <= a || b <= c )
    c = 0.5 * ( a + b );
  end

  y1 = f ( c );
  calls = calls + 1;
  k = 3;
  d0 = a2 - c;
  h = 9.0 / 11.0;

  if ( y1 < y )
    x = c;
    y = y1;
  end

  while ( true )

    d1 = a2 - a0;
    d2 = c - a0;
    z2 = b - a2;
    z0 = y2 - y1;
    z1 = y2 - y0;
    r = d1 * d1 * z0 - d0 * d0 * z1;
    p = r;
    qs = 2.0 * ( d0 * z1 - d1 * z0 );
    q = qs;
%
%  Loop control corrected by Hans Bieshaar, 28 May 2021.
%
    force_first = true;

    if ( 100000 < k && y < y2 )
      k = mod ( 1611 * k, 1048576 );
      q = 1.0;
      r = ( b - a ) * 0.00001 * k;
      force_first = false;
    end

    while ( r < z2 || force_first )

      force_first = false;

      if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) < ...
        z2 * m2 * r * ( z2 * q - r ) )
        a3 = a2 + r / q;
        y3 = f ( a3 );
        calls = calls + 1;

        if ( y3 < y )
          x = a3;
          y = y3;
        end
      end

      k = mod ( 1611 * k, 1048576 );
      q = 1.0;
      r = ( b - a ) * 0.00001 * k;

    end

    r = m2 * d0 * d1 * d2;
    s = sqrt ( ( ( y2 - y ) + t ) / m2 );
    h = 0.5 * ( 1.0 + h );
    p = h * ( p + 2.0 * r * s );
%
%  Correction supplied by Hans Bieshaar, 27 May 2021.
%
    q = r + 0.5 * qs;
    r = - 0.5 * ( d0 + ( z0 + 2.01 * e ) / ( d0 * m2 ) );

    if ( r < s || d0 < 0.0 )
      r = a2 + s;
    else
      r = a2 + r;
    end

    if ( 0.0 < p * q )
      a3 = a2 + p / q;
    else
      a3 = r;
    end

    while ( true )

      a3 = max ( a3, r );

      if ( b <= a3 )
        a3 = b;
        y3 = yb;
      else
        y3 = f ( a3 );
        calls = calls + 1;
      end

      if ( y3 < y )
        x = a3;
        y = y3;
      end

      d0 = a3 - a2;

      if ( a3 <= r )
        break
      end

      p = 2.0 * ( y2 - y3 ) / ( m * d0 );

      if ( ( 1.0 + 9.0 * eps ) * d0 <= abs ( p ) )
        break
      end

      if ( 0.5 * m2 * ( d0 * d0 + p * p ) <= ...
        ( y2 - y ) + ( y3 - y ) + 2.0 * t )
        break
      end

      a3 = 0.5 * ( a2 + a3 );
      h = 0.9 * h;

    end

    if ( b <= a3 )
      break
    end

    a0 = c;
    c = a2;
    a2 = a3;
    y0 = y1;
    y1 = y2;
    y2 = y3;

  end

  fx = y;
 
  return
end
 
