function root_rc ( x, fx, ferr, xerr, q )

!*****************************************************************************80
!
!! root_rc() solves a single nonlinear equation using reverse communication.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 January 2013
!
!  Author:
!
!    Original FORTRAN77 version by Gaston Gonnet.
!    This version by John Burkardt.
!
!  Reference:
!     
!    Gaston Gonnet,
!    On the Structure of Zero Finders,
!    BIT Numerical Mathematics,
!    Volume 17, Number 2, June 1977, pages 170-183.
!
!  Parameters:
!
!    Input, real ( kind = rk ) X, an estimate for the root.  On the first
!    call, this must be a value chosen by the user.  Thereafter, it may
!    be a value chosen by the user, or the value of ROOT returned on the
!    previous call to the function.
!
!    Input, real ( kind = rk ) FX, the value of the function at X.
!
!    Output, real ( kind = rk ) FERR, the smallest value of F encountered.
!
!    Output, real ( kind = rk ) XERR, the width of the change-in-sign interval,
!    if one was encountered.
!
!    Input/output, real ( kind = rk ) Q(9), storage needed by the function.
!    Before the first call, the user must set Q(1) to 0.
!
!    Output, real ( kind = rk ) ROOT_RC, an improved estimate for the root.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) d
  real ( kind = rk ) decr
  real ( kind = rk ) ferr
  real ( kind = rk ) fx
  integer i
  real ( kind = rk ) p
  real ( kind = rk ) q(9)
  real ( kind = rk ) r
  real ( kind = rk ) r8_sign
  real ( kind = rk ) root_rc
  real ( kind = rk ) u
  real ( kind = rk ) v
  real ( kind = rk ) w
  real ( kind = rk ) x
  real ( kind = rk ) xerr
  real ( kind = rk ) z
!
!  If we found an exact zero, there is nothing more to do.
!
  if ( fx == 0.0D+00 ) then
    ferr = 0.0D+00
    xerr = 0.0D+00
    root_rc = x
    return
  end if

  ferr = abs ( fx )
!
!  If this is the first time, initialize, estimate the first root, and exit.
!
  if ( q(1) == 0.0D+00 ) then
    q(1) = fx
    q(2) = x
    q(3:9) = 0.0D+00
    root_rc = x + fx
    xerr = huge ( xerr )
    return
  end if
!
!  This is not the first call.
!
  q(9) = q(9) + 1.0D+00
!
!  Check for too many iterations.
!
  if ( 80.0D+00 < q(9) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ROOT_RC - Fatal error!'
    write ( *, '(a,i6)' ) '  Number of iterations = ', int ( q(9) )
    stop
  end if
!
!  Check for a repeated X value.
!
  if ( ( 2.0 <= q(9) .and. x == q(4) ) .or. x == q(2) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ROOT_RC - Fatal error!'
    write ( *, '(a,i6)' ) '  Value of X has been input before.'
    stop
  end if
!
!  Push X -> A -> B -> C
!
  do i = 6, 3, -1
    q(i) = q(i-2)
  end do
  q(1) = fx
  q(2) = x
!
!  If we have a change-in-sign interval, store the opposite value.
!
  if ( r8_sign ( q(1) ) /= r8_sign ( q(3) ) ) then
    q(7) = q(3)
    q(8) = q(4)
  end if
!
!  Calculate XERR.
!
  if ( q(7) /= 0.0D+00 ) then
    xerr = abs ( q(8) - q(2) )
  else
    xerr = huge ( xerr )
  end if
!
!  If more than 30 iterations, and we have change-in-sign interval, bisect.
!
  if ( 30.0D+00 < q(9) .and. q(7) /= 0.0D+00 ) then
    root_rc = q(2) + ( q(8) - q(2) ) / 2.0D+00
    return
  end if

  v = ( q(3) - q(1) ) / ( q(4) - q(2) )
!
!  If 3 or more points, try Muller.
!
  if ( q(5) /= 0.0D+00 ) then
    u = ( q(5) - q(3) ) / ( q(6) - q(4) )
    w = q(4) - q(2)
    z = ( q(6) - q(2) ) / w
    r = ( z + 1.0D+00 ) * v - u

    if ( r /= 0.0D+00 ) then
      p = 2.0D+00 * z * q(1) / r
      d = 2.0D+00 * p / ( w * r ) * ( v - u )
      if ( -1.0D+00 <= d ) then
        root_rc = q(2) - p / ( 1.0D+00 + sqrt ( 1.0D+00 + d ) )
        if ( q(7) == 0.0D+00 .or. &
             ( q(2) < root_rc .and. root_rc < q(8) ) .or. &
             ( q(8) < root_rc .and. root_rc < q(2) ) ) then
          return
        end if
      end if
    end if
  end if
!
!  Try the secant step.
!
  if ( q(1) /= q(3) .or. q(7) == 0.0D+00 ) then
    if ( q(1) == q(3) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'ROOT_RC - Fatal error!'
      write ( *, '(a)' ) '  Cannot apply any method.'
      stop
    end if
    decr = q(1) / v
    if ( abs ( decr ) * 4.6D+18 < abs ( q(2) ) ) then
      decr = 1.74D-18 * abs ( q(2) ) * r8_sign ( decr )
    end if
    root_rc = q(2) - decr
    if ( q(7) == 0.0D+00 .or. &
        ( q(2) < root_rc .and. root_rc < q(8) ) .or. &
        ( q(8) < root_rc .and. root_rc < q(2) ) ) then
      return
    end if
  end if
!
!  Apply bisection.
!
  root_rc = q(2) + ( q(8) - q(2) ) / 2.0D+00

  return
end
function r8_sign ( x )

!*****************************************************************************80
!
!! r8_sign() returns the sign of an R8.
!
!  Discussion:
!
!    value = -1 if X < 0;
!    value = +1 if X => 0.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 March 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) X, the number whose sign is desired.
!
!    Output, real ( kind = rk ) R8_SIGN, the sign of X:
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) r8_sign
  real ( kind = rk ) x

  if ( x < 0.0D+00 ) then
    r8_sign = -1.0D+00
  else
    r8_sign = +1.0D+00
  end if

  return
end