function [ x, x1, ierror, k ] = brent ( x, x1, abserr, kmax, f )

%*****************************************************************************80
%
%% brent() implements the Brent bisection-based zero finder.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization without Derivatives,
%    Dover, 2002,
%    ISBN: 0-486-41998-3,
%    LC: QA402.5.B74.
%
%  Input:
%
%    real X, X1.  Two points defining
%    the interval in which the search will take place.  F(X) and F(X1)
%    should have opposite signs.
%
%    real ABSERR, an error tolerance.
%
%    integer KMAX, the maximum number of iterations allowed.
%
%    real external F, the name of the routine that
%    evaluates the function or its derivatives, of the form
%      function f ( x, ider )
%
%  Output:
%
%    real X, X1.  X is the best estimate for
%    the root, and X1 is a recently computed point such that F changes
%    sign in the interval between X and X1.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%

%
%  Initialization.
%
  ierror = 0;
  k = 0;
  fx = f ( x, 0 );
  fx1 = f ( x1, 0 );
  x2 = x1;
  fx2 = fx1;
%
%  Iteration loop:
%
  while ( true )

    k = k + 1;

    if ( kmax < k );
      ierror = 2;
      return;
    end

    if ( ( 0.0 < fx1 ) == ( 0.0 < fx2 ) )
      x2 = x;
      fx2 = fx;
      d = x1 - x;
      e = x1 - x;
    end

    if ( abs ( fx2 ) < abs ( fx1 ) )
      t = x1;
      x1 = x2;
      x2 = t;
      t = fx1;
      fx1 = fx2;
      fx2 = t;
    end

    tol = 2.0 * eps * abs ( x1 ) + abserr;
    em = 0.5 * ( x2 - x1 );

    if ( abs ( em ) <= tol || fx1 == 0.0 )
      x = x1;
      return
    end

    if ( abs ( e ) < tol || abs ( fx ) <= abs ( fx1 ) )

      d = em;
      e = em;

    else

      s = fx1 / fx;

      if ( x == x2 )
        p = 2.0 * em * s;
        q = 1.0 - s;
      else
        q = fx / fx2;
        r = fx1 / fx2;
        p = s * ( 2.0 * em * q * ( q - r ) - ( x1 - x ) * ( r - 1.0 ) );
        q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
      end

      if ( 0 < p )
        q = - q;
      else
        p = - p;
      end

      s = e;
      e = d;

      if ( 2.0 * p < 3.0 * em * q - abs ( tol * q ) || p < abs ( 0.5 * s * q ) )
        d = p / q;
      else
        d = em;
        e = em;
      end

    end
%
%  Set the increment.
%
    if ( tol < abs ( d ) )
      dx = + d;
    elseif ( 0.0 < em )
      dx = + tol;
    else
      dx = - tol;
    end
%
%  Remember current data for next step.
%
    x = x1;
    fx = fx1;
%
%  Update the iterate and function values.
%
    x1 = x1 + dx;
    fx1 = f ( x1, 0 );

  end

  return
end
