function [ value, calls ] = zero_brent ( a, b, t, f )

%*****************************************************************************80
%
%% zero_brent() seeks a root of a function F(X) in an interval [A,B].
%
%  Discussion:
%
%    The interval [A,B] must be a change of sign interval for F.
%    That is, F(A) and F(B) must be of opposite signs.  Then
%    assuming that F is continuous implies the existence of at least
%    one value C between A and B for which F(C) = 0.
%
%    The location of the zero is determined to within an accuracy
%    of 4 * EPSILON * abs ( C ) + 2 * T, where EPSILON is the machine epsilon.
%
%    Thanks to Thomas Secretin for pointing out a transcription error in the
%    setting of the value of P, 11 February 2013.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 July 2021
%
%  Author:
%
%    Original FORTRAN77 version by Richard Brent.
%    This 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 change of sign interval.
%
%    real T, a positive error tolerance.
%
%    real value = F ( x ), the name of a user-supplied
%    function which evaluates the function whose zero is being sought.
%
%  Output:
%
%    real VALUE, the estimated value of a zero of the function F.
%
%    integer CALLS: the number of calls to F.
%
  calls = 0;
%
%  Make local copies of A and B.
%
  sa = a;
  fa = f ( sa );
  calls = calls + 1;

  sb = b;
  fb = f ( sb );
  calls = calls + 1;

  c = sa;
  fc = fa;
  e = sb - sa;
  d = e;

  while ( true )

    if ( abs ( fc ) < abs ( fb ) )

      sa = sb;
      sb = c;
      c = sa;
      fa = fb;
      fb = fc;
      fc = fa;

    end

    tol = 2.0 * eps * abs ( sb ) + t;
    m = 0.5 * ( c - sb );

    if ( abs ( m ) <= tol || fb == 0.0 )
      break
    end

    if ( abs ( e ) < tol || abs ( fa ) <= abs ( fb ) )

      e = m;
      d = e;

    else

      s = fb / fa;

      if ( sa == c )

        p = 2.0 * m * s;
        q = 1.0 - s;

      else

        q = fa / fc;
        r = fb / fc;
        p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
        q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );

      end

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

      s = e;
      e = d;

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

    end

    sa = sb;
    fa = fb;

    if ( tol < abs ( d ) )
      sb = sb + d;
    elseif ( 0.0 < m )
      sb = sb + tol;
    else
      sb = sb - tol;
    end

    fb = f ( sb );
    calls = calls + 1;

    if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
      c = sa;
      fc = fa;
      e = sb - sa;
      d = e;
    end

  end

  value = sb;

  return
end
