zero_brent <- function ( a, b, t, f )

#*****************************************************************************80
#
## zero_brent() seeks the 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 * EPS * abs ( C ) + 2 * T, where EPS is the machine epsilon.
#
#  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 ), 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 function calls.
#
{
  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
    }

    tol <- 2.0 * .Machine$double.eps * abs ( sb ) + t
    m <- 0.5 * ( c - sb )

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

    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 )
      }

      if ( 0.0 < p )
      {
        q <- - q
      }
      else
      {
        p <- - p
      }

      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
      }

    }

    sa <- sb
    fa <- fb

    if ( tol < abs ( d ) )
    {
      sb <- sb + d
    }
    else if ( 0.0 < m )
    {
      sb <- sb + tol
    }
    else
    {
      sb <- sb - tol
    }

    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
    }

  }

  value <- sb

  output <- c ( value, calls )

  return ( output )
}