function bisect ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% bisect() carries out the bisection method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  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.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real t
  real x
  real x1
  real x2
%
%  Initialization.
%
  ierror = 0
  k = 0
%
%  Evaluate the function at the starting points.
%
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )
%
%  Set XPOS and XNEG to the X values for which F(X) is positive
%  and negative, respectively.
%
  if ( 0.0 <= fx .and. fx1 <= 0.0 )

  else if ( fx <= 0.0 .and. 0.0 <= fx1 )
    t = x
    x = x1
    x1 = t
    t = fx
    fx = fx1
    fx1 = t
  else
    ierror = 1
    return
  end if
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

    if ( kmax < k )
      ierror = 2
      exit
    end if
%
%  Set the increment.
%
    dx = 0.5 * ( x1 - x )
%
%  Update the iterate and function values.
%
    x2 = x + dx
    fx2 = f ( x2, 0 )

    if ( 0.0 <= fx2 )
      x = x2
      fx = fx2
    else
      x1 = x2
      fx1 = fx2
    end if

  end do

  return
end
function brent ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************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.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d
  real dx
  real e
  real em
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real p
  real q
  real r
  real s
  real t
  real tol
  real x
  real x1
  real x2
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )
  x2 = x1
  fx2 = fx1
%
%  Iteration loop:
%
  do

    k = k + 1

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

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

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

    tol = 2.0 * epsilon ( 1.0 ) * abs ( x1 ) + abserr
    em = 0.5 * ( x2 - x1 )

    if ( abs ( em ) <= tol .or. fx1 == 0.0 )
      x = x1
      return
    end if

    if ( abs ( e ) < tol .or. 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

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

      s = e
      e = d

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

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

  end do

  return
end
function c8_muller ( func, fatol, itmax, x1, x2, x3, xatol, xrtol, ...
  xnew, fxnew )

%*****************************************************************************80
%
%% c8_muller() carries out Muller's method, using C8 arithmetic.
%
%  Discussion:
%
%    "C8" arithmetic is complex double precision arithmetic.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gisela Engeln-Muellges, Frank Uhlig,
%    Numerical Algorithms with C,
%    Springer, 1996,
%    ISBN: 3-540-60530-4,
%    LC: QA297.E56213.
%
%  Input:
%
%    external FUNC, the name of the routine that evaluates the function.
%    FUNC should have the form:
%      function func ( x, fx )
%      complex fx
%      complex x
%
%    real FATOL, the absolute error tolerance for F(X).
%
%    integer ITMAX, the maximum number of steps allowed.
%
%    complex ( kind = ck ) X1, X2, X3, three distinct points to start the
%    iteration.
%
%    real XATOL, XRTOL, absolute and relative
%    error tolerances for the root.
%
%  Output:
%
%    complex ( kind = ck ) XNEW, the estimated root.
%
%    complex ( kind = ck ) FXNEW, the value of the function at XNEW.
%
  implicit none

  integer, parameter :: ck = kind ( ( 1.0, 1.0 ) )
  integer, parameter :: rk = kind ( 1.0 )

  complex ( kind = ck ) a
  complex ( kind = ck ) b
  complex ( kind = ck ) c
  complex ( kind = ck ) c8_temp
  complex ( kind = ck ) discrm
  real fatol
  complex ( kind = ck ) fminus
  complex ( kind = ck ) fplus
  external func
  complex ( kind = ck ) fxmid
  complex ( kind = ck ) fxnew
  complex ( kind = ck ) fxold
  integer iterate
  integer itmax
  real x_ave
  complex ( kind = ck ) x_inc
  complex ( kind = ck ) x1
  complex ( kind = ck ) x2
  complex ( kind = ck ) x3
  real xatol
  complex ( kind = ck ) xlast
  complex ( kind = ck ) xmid
  complex ( kind = ck ) xminus
  complex ( kind = ck ) xnew
  complex ( kind = ck ) xold
  complex ( kind = ck ) xplus
  real xrtol

  xnew = x1
  xmid = x2
  xold = x3

  call func ( xnew, fxnew )
  call func ( xmid, fxmid )
  call func ( xold, fxold )
'
  fprintf ( 1, ' '
  fprintf ( 1, 'c8_muller():'
  fprintf ( 1, '  Muller''s method (complex root version)'
  fprintf ( 1, ' '
  write ( *, '(a)' ) ...
    '  Iteration     x_real              x_imag' // ...
    '             ||fx||           ||disc||'
  fprintf ( 1, ' '

  iterate = -2
  write ( *, '(i6,f20.10,f20.10,f20.10)' ) iterate, xold, abs ( fxold )
  iterate = -1
  write ( *, '(i6,f20.10,f20.10,f20.10)' ) iterate, xmid, abs ( fxmid )
  iterate = 0
  write ( *, '(i6,f20.10,f20.10,f20.10)' ) iterate, xnew, abs ( fxnew )

  if ( abs ( fxnew ) < fatol )
    fprintf ( 1, ' '
    fprintf ( 1, 'c8_muller():'
    fprintf ( 1, '  |F(X)| is below the tolerance.'
    return
  end if

  do
%
%  You may need to swap (XMID,FXMID) and (XNEW,FXNEW).
%
    if ( abs ( fxmid ) <= abs ( fxnew ) )

      c8_temp = xnew
      xnew = xmid
      xmid = c8_temp

      c8_temp = fxnew
      fxnew = fxmid
      fxmid = c8_temp

    end if

    xlast = xnew
    iterate = iterate + 1

    if ( itmax < iterate )
      fprintf ( 1, ' '
      fprintf ( 1, 'c8_muller():'
      fprintf ( 1, '  Maximum number of steps taken.'
      exit
    end if

    a =  ( ( xmid - xnew ) * ( fxold - fxnew ) ...
         - ( xold - xnew ) * ( fxmid - fxnew ) )

    b = ( ( xold - xnew )**2 * ( fxmid - fxnew ) ...
        - ( xmid - xnew )**2 * ( fxold - fxnew ) )

    c = ( ( xold - xnew ) * ( xmid - xnew ) * ( xold - xmid ) * fxnew )

    xold = xmid
    xmid = xnew
%
%  Apply the quadratic formula to get roots XPLUS and XMINUS.
%
    discrm = b ** 2 - 4.0 * a * c

    if ( a == 0.0 )
      fprintf ( 1, ' '
      fprintf ( 1, 'c8_muller():'
      fprintf ( 1, '  The algorithm has broken down.'
      fprintf ( 1, '  The quadratic coefficient A is zero.'
      exit
    end if

    xplus = xnew + ( ( - b + sqrt ( discrm ) ) / ( 2.0 * a ) )

    call func ( xplus, fplus )

    xminus = xnew + ( ( - b - sqrt ( discrm ) ) / ( 2.0 * a ) )

    call func ( xminus, fminus )
%
%  Choose the root with smallest function value.
%
    if ( abs ( fminus ) < abs ( fplus ) )
      xnew = xminus
    else
      xnew = xplus
    end if

    fxold = fxmid
    fxmid = fxnew
    call func ( xnew, fxnew )
    write ( *, '(i6,f20.10,f20.10,f20.10,f20.10)' ) ...
      iterate, xnew, abs ( fxnew ), abs ( discrm )
%
%  Check for convergence.
%
    x_ave = abs ( xnew + xmid + xold ) / 3.0
    x_inc = xnew - xmid

    if ( abs ( x_inc ) <= xatol )
      fprintf ( 1, ' '
      fprintf ( 1, 'c8_muller():'
      fprintf ( 1, '  Absolute convergence of the X increment.'
      exit
    end if

    if ( abs ( x_inc ) <= xrtol * x_ave )
      fprintf ( 1, ' '
      fprintf ( 1, 'c8_muller():'
      fprintf ( 1, '  Relative convergence of the X increment.'
      exit
    end if

    if ( abs ( fxnew ) <= fatol )
      fprintf ( 1, ' '
      fprintf ( 1, 'c8_muller():'
      fprintf ( 1, '  Absolute convergence of |F(X)|.'
      exit
    end if

  end do

  return
end
function cap_phi_03 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% cap_phi_03() implements the Traub capital Phi(0,3) method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 232.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: if IERROR = 0, X is an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    z = 1.0 - 2.0 * fx * d2fx / dfx**2
%
%  Set the increment.
%
    if ( 0.0 < z )
      dx = - 2.0 * ( fx / dfx ) / ( 1.0 + sqrt ( z ) )
    else
      dx = - fx / dfx
    end if
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function cap_phi_21 ( x, x1, x2, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% cap_phi_21() implements the Traub capital PHI(2,1) function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964.
%
%  Input:
%
%    real X, X1, X2: three distinct points that start the method.
%
%    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, X2: X is an approximation which satisfies 
%    abs ( F(X) ) < ABSERR, and X1 and X2 are the previous estimates.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    Output, integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d1
  real d2
  real det
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real x
  real x1
  real x2
  real z
%
%  Initialization.
%
  ierror = 0
  k = -2
  fx2 = f ( x2, 0 )
  k = -1
  fx1 = f ( x1, 0 )
  k = 0
  fx = f ( x, 0 )

  if ( x1 == x2 )
    ierror = 3
    return
  end if

  d1 = ( fx1 - fx2 ) / ( x1 - x2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    d2 = d1

    if ( x == x1 )
      ierror = 3
      return
    end if

    d1 = ( fx - fx1 ) / ( x - x1 )

    if ( x == x2 )
      ierror = 3
      return
    end if

    d2 = ( d1 - d2 ) / ( x - x2 )

    z = d1 + ( x - x1 ) * d2
    det = z**2 - 4.0 * fx * d2
    det = max ( det, 0.0 )
%
%  Set the increment.
%
    dx = - 2.0 * fx / ( z + sqrt ( det ) )
%
%  Remember current data for next step.
%
    x2 = x1
    fx2 = fx1

    x1 = x
    fx1 = fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function chebyshev ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% chebyshev() implements Chebyshev's method.
%
%  Discussion:
%
%    This is also known as the Euler-Chebyshev method.
%
%    x(n+1) = x(n) - ( f(x(n)) / f'(x(n)) ) * ( 1 + 0.5 * L(x(n)) )
%
%    where
%
%      L(x) = f(x) * f''(x) / f'(x)**2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, an estimate for the root of the equation.
%
%    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: if IERROR = 0, X is an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = d2fx * fx / dfx**2

    if ( 1.0 + 0.5 * u == 0.0 )
      ierror = 4
      return
    end if
%
%  Set the increment.
%
    dx = - ( fx / dfx ) * ( 1.0 + 0.5 * u )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function dagger_e12 ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% dagger_e12() implements the dagger E 1,2 algorithm.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 234.
%
%  Input:
%
%    real X, X1, distinct points that start the method.
%
%    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 an approximation to a root of the equation
%    which satisfies abs ( F(X) ) < ABSERR, and X1 is the
%    previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d
  real dfx
  real dfx1
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real x1
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )

  dfx1 = f ( x1, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    if ( x - x1 == 0.0 )
      ierror = 3
      return
    end if

    d = ( dfx - dfx1 ) / ( x - x1 )
%
%  Set the increment.
%
    dx = - ( fx / dfx ) - 0.5 * ( fx / dfx )**2 * d / dfx
%
%  Remember current data for next step.
%
    x1 = x
    dfx1 = dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function e3 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% e3() implements the Traub E3 method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 232.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real v
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    v = d2fx / ( 2.0 * dfx )
%
%  Set the increment.
%
    dx = - ( fx / dfx ) * ( 1.0 + v * u )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function e4 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% e4() implements the Traub E4 method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 232.
%
%  Input:
%
%    real X, an estimate for the root of the equation.
%
%    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:
%
%    Ireal X: an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real d3fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real v
  real w
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
  d3fx = f ( x, 3 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    v = d2fx / ( 2.0 * dfx )
    w = d3fx / ( 6.0 * dfx )
%
%  Set the increment.
%
    dx = - u * ( 1.0 + u * ( v + u * ( 2.0 * v**2 - w ) ) )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )
    d3fx = f ( x, 3 )

  end do

  return
end
function euler ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% euler() implements the Euler method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: f IERROR = 0, X is an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    z = dfx**2 - 2.0 * fx * d2fx
%
%  Set the increment.
%
    if ( 0.0 < z )

      if ( dfx + sqrt ( z ) == 0.0 )
        ierror = 3
        return
      end if

      dx = - 2.0 * fx / ( dfx + sqrt ( z ) )

    else

      if ( dfx == 0.0 )
        ierror = 3
        return
      end if

      dx = - fx / dfx

    end if
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function halley1 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% halley1() implements Halley's method.
%
%  Discussion:
%
%    x(n+1) = x(n) - ( f(x(n)) / f'(x(n)) ) / ( 1 - 0.5 * L(x(n)) )
%
%    where
%
%      L(x) = f(x) * f''(x) / f'(x)**2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = d2fx * fx / dfx**2

    if ( 2.0 - u == 0.0 )
      ierror = 4
      return
    end if
%
%  Set the increment.
%
    dx = - ( fx / dfx ) / ( 1.0 - 0.5 * u )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function halley2 ( x, x1, x2, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% halley2() implements Halley's method, with finite differences.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, X1, X2: three distinct points that start the method.
%
%    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, X2: X is an approximate root of the equation
%    which satisfies abs ( F(X) ) < ABSERR, and X1 and X2 are the
%    previous estimates.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d1
  real d2
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real x
  real x1
  real x2
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )
  fx2 = f ( x2, 0 )

  d1 = ( fx1 - fx2 ) / ( x1 - x2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    d2 = d1

    if ( x == x1 )
      ierror = 3
      return
    end if

    d1 = ( fx - fx1 ) / ( x - x1 )

    if ( x == x2 )
      ierror = 3
      return
    end if

    d2 = ( d1 - d2 ) / ( x - x2 )

    if ( d1 == 0.0 )
      ierror = 3
      return
    end if

    z = d1 - fx1 * d2 / d1

    if ( z == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx / z
%
%  Remember current data for next step.
%
    x2 = x1
    fx2 = fx1
    x1 = x
    fx1 = fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function halley_super ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% halley_super() implements the super Halley method.
%
%  Discussion:
%
%    x(n+1) = x(n) - 0.5 * ( f(x(n)) / f'(x(n)) )
%      * ( 2 - L(x(n)) ) / ( 1 - L(x(n)) )
%
%    where
%
%      L(x) = f(x) * f''(x) / f'(x)**2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = d2fx * fx / dfx**2

    if ( 1.0 - u == 0.0 )
      ierror = 4
      return
    end if
%
%  Set the increment.
%
    dx = - 0.5 * ( fx / dfx ) * ( 2.0 - u ) / ( 1.0 - u )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function hansen ( x, beta, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% hansen() implements the Hansen and Patrick method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eldon Hansen, Merrell Patrick,
%    A Family of Root Finding Methods,
%    Numerische Mathematik,
%    Volume 27, 1977, pages 257 - 269.
%
%  Input:
%
%    real X, an estimate for the root of the equation.
%
%    real BETA, ???
%
%    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: if IERROR = 0, X is an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real beta
  real bot
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )

  write ( *, '(a,i4,2g14.6)' ) 'KXFXDFX', k, x, fx
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    z = dfx**2 - ( beta + 1.0 ) * fx * d2fx
    z = max ( z, 0.0 )

    bot = beta * dfx + sqrt ( z )

    if ( bot == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - ( beta + 1.0 ) * fx / bot
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

    write ( *, '(a,i4,3g14.6)' ) 'HANSEN - KXFXDFX', k, x, fx, dfx

  end do

  return
end
function jarratt ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% jarratt() implements the Jarratt fourth order method.
%
%  Discussion:
%
%    Jarratt's method is of fourth order.
%
%    x(n+1) = x(n) - 1/2 * ( f(x(n) / f'(x(n)) )
%      + f(x(n)) / ( f'(x(n)) - 3*f'( x(n)-(2/3)*f(x(n))/f'(x(n)) ) )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    P Jarratt,
%    Some fourth-order multipoint iterative methods for solving equations,
%    Mathematics of Computation,
%    Volume 20, Number 95, 1966, pages 434 - 437.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfz
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    z = x - 2.0 * fx / ( 3.0 * dfx )
    dfz = f ( z, 1 )

    if ( dfx - 3.0 * dfz == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = ( - 0.5 / dfx + 1.0 / ( dfx - 3.0 * dfz ) ) * fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function jarratt2 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% jarratt2() implements the inverse-free Jarratt fourth order method.
%
%  Discussion:
%
%    Jarratt's inverse-free method is of fourth order.
%
%    x(n+1) = x(n) - u(x(n)) + (3/4) * u(x(n)) * h(x(n)) * ( 1-(3/2)*h(x(n)) )
%
%    where
%
%      u(x) = f(x) / f'(x)
%      h(x) = ( f'(x-(2/3)*u(x)) - f'(x) ) / f'(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    P Jarratt,
%    Some fourth-order multipoint iterative methods for solving equations,
%    Mathematics of Computation,
%    Volume 20, Number 95, 1966, pages 434 - 437.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  real hx
  integer ierror
  integer k
  integer kmax
  real x
  real ux
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    ux = fx / dfx

    hx = ( f ( x - 2.0 / 3.0 * ux, 1 ) - dfx ) / dfx
%
%  Set the increment.
%
    dx = - ux + 0.75 * ux * hx * ( 1.0 - 1.5 * hx )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function king ( x, beta, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% king() implements a family of fourth order methods.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard King,
%    A family of fourth order methods,
%    SIAM Journal on Numerical Analysis,
%    Volume 10, 1973, pages 876 - 879.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    real BETA, a parameter in the algorithm.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real beta
  real dfx
  real dx
  real, external :: f
  real fw
  real fx
  integer ierror
  integer k
  integer kmax
  real w
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    w = x - fx / dfx
    fw = f ( w, 0 )

    if ( fx + ( beta - 2.0 ) * fw == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - ( fx / dfx ) - ( fw / dfx ) ...
      * ( fx + beta * fw ) / ( fx + ( beta - 2.0 ) * fw )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function laguerre ( x0, ipoly, abserr, kmax, f, x, ierror, k )

%*****************************************************************************80
%
%% laguerre() implements the Laguerre rootfinding method for polynomials.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eldon Hansen, Merrell Patrick,
%    A Family of Root Finding Methods,
%    Numerische Mathematik,
%    Volume 27, 1977, pages 257 - 269.
%
%  Input:
%
%    real X0.
%    An initial estimate for the root of the equation.
%
%    integer IPOLY: the polynomial degree of the function.
%    IPOLY must be at least 2.
%
%    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: the estimated solution, if IERROR=0.
%
%    integer IERROR: error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K: the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real beta
  real bot
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer ipoly
  integer k
  integer kmax
  real x
  real x0
  real z
%
%  Check.
%
  if ( ipoly < 2 )
    ierror = 1
    fprintf ( 1, ' '
    fprintf ( 1, 'laguerre(): Fatal error%'
    fprintf ( 1, '  IPOLY < 2.'
    return
  end if
%
%  Initialization.
%
  ierror = 0
  x = x0

  beta = 1.0 / real ( ipoly - 1, kind = rk )

  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    z = dfx**2 - ( beta + 1.0 ) * fx * d2fx
    z = max ( z, 0.0 )

    bot = beta * dfx + sqrt ( z )

    if ( bot == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - ( beta + 1.0 ) * fx / bot
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function midpoint ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% midpoint() implements the midpoint method.
%
%  Discussion:
%
%    The midpoint method is of order 3.
%
%    x(n+1) = x(n) - f(x(n)) / g(x(n))
%
%    where
%
%    g(x) = f'( x - f(x) / ( 2 * f'(x) ) )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 164.
%
%  Input:
%
%    real X, the point that starts the method.
%
%    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: an approximate root satisfying abs ( F(X) ) < ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real arg
  real dfx
  real dx
  real, external :: f
  real fx
  real gx
  integer ierror
  integer k
  integer kmax
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0

  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    arg = x - fx / ( 2.0 * dfx )
    gx = f ( arg, 1 )

    if ( gx == 0.0 )
      ierror = 4
      return
    end if
%
%  Set the increment.
%
    dx = - fx / gx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function newton ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% newton() implements Newton's method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx / dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function newton_mod ( x, nsub, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% newton_mod() implements the modified Newton method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 236.
%
%  Input:
%
%    real X: n estimate for the root of the equation.
%
%    integer NSUB, the number of steps in the sub-iteration.
%    The derivative is only evaluated before the first step,
%    and after every subsequent NSUB steps.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  integer nsub
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( mod ( k - 1, nsub ) == 0 )

      dfx = f ( x, 1 )

      if ( dfx == 0.0 )
        ierror = 3
        return
      end if

    end if
%
%  Set the increment.
%
    dx = - fx / dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function newton_sec ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% newton_sec() implements the Newton - secant method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 236.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  real fxu
  integer ierror
  integer k
  integer kmax
  real x
  real x1
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    x1 = x - fx / dfx
    fxu = f ( x1, 0 )

    if ( fxu - fx == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - ( fx / dfx ) * ( 1.0 - fxu / ( fxu - fx ) )
%
%  Remember current data for next step.
%
    x1 = x
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function ostrowski_sqrt ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% ostrowski_sqrt() implements the Ostrowski square root method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eldon Hansen, Merrell Patrick,
%    A Family of Root Finding Methods,
%    Numerische Mathematik,
%    Volume 27, 1977, pages 257 - 269.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    real BETA, ???
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real bot
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    z = dfx**2 - fx * d2fx
%
%  Set the increment.
%
    if ( 0.0 < z )
      bot = sqrt ( z )
    else
      bot = dfx
    end if

    dx = - fx / bot
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function perp_e_12 ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% perp_e_12() implements the Traub E 1,2 algorithm.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 234.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    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 an approximate root
%    satisfiying abs ( F(X) ) < ABSERR, and X1 is the previous
%    estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfx1
  real dx
  real, external :: f
  real fx
  real fx1
  integer ierror
  integer k
  integer kmax
  real x
  real x1
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  fx1 = f ( x1, 0 )
  dfx1 = f ( x1, 1 )

  if ( dfx1 == 0.0 )
    ierror = 3
    return
  end if
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    if ( dfx1 == 0.0 )
      ierror = 3
      return
    end if

    if ( fx == fx1 )
      ierror = 3
      return
    end if

    z = 2.0 / dfx + 1.0 / dfx1 - 3.0 * ( x - x1 ) / ( fx - fx1 )
%
%  Set the increment.
%
    dx = - fx / dfx + fx**2 * z / ( fx - fx1 )
%
%  Remember current data for next step.
%
    x1 = x
    fx1 = fx
    dfx1 = dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function perp_e_21 ( x, x1, x2, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% perp_e_21() implements the Traub perp E 21 method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 233.
%
%  Input:
%
%    real X, X1, X2: three distinct points that start the method.
%
%    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, X2: X is an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 and X2 are the previous estimates.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d
  real d1
  real d2
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real x
  real x1
  real x2
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )
  fx2 = f ( x2, 0 )

  if ( x1 == x2 )
    ierror = 3
    return
  end if

  d1 = ( fx1 - fx2 ) / ( x1 - x2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    d2 = d1

    if ( x == x1 )
      ierror = 3
      return
    end if

    if ( x == x2 )
      ierror = 3
      return
    end if

    d1 = ( fx - fx1 ) / ( x - x1 )
    d = ( fx - fx2 ) / ( x - x2 )
%
%  Set the increment.
%
    dx = - fx * ( 1.0 / d1 + 1.0 / d - 1.0 / d2 )
%
%  Remember current data for next step.
%
    x2 = x1
    fx2 = fx1
    x1 = x
    fx1 = fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function phi_12 ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% phi_12() implements the Traub capital PHI(1,2) method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 233.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    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 an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real c
  real d
  real dfx
  real dfx1
  real dx
  real, external :: f
  real fx
  real fx1
  real h
  integer ierror
  integer k
  integer kmax
  real x
  real x1
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )
  dfx1 = f ( x1, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    c = fx - fx1

    if ( x == x1 )
      ierror = 3
      return
    end if

    d = ( fx - fx1 ) / ( x - x1 )

    h = 1.0 / c * ( 1.0 / dfx - 1.0 / d ) - ( fx1 / c**2 ) * ...
      ( 1.0 / dfx + 1.0 / dfx1 - 2.0 / d )
%
%  Set the increment.
%
    dx = - ( fx / dfx ) + fx**2 * h
%
%  Remember current data for next step.
%
    x1 = x
    fx1 = fx
    dfx1 = dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function psi_21 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% psi_21() implements the Traub PSI 2,1 method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 232.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real d3fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real v
  real w
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )
    d3fx = f ( x, 3 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    v = d2fx / ( 2.0 * dfx )
    w = d3fx / ( 6.0 * dfx )
%
%  Set the increment.
%
    dx = - u * ( v - u * ( v**2 - w ) ) / ( v - u * ( 2.0 * v**2 - w ) )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function psi_i2 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% psi_12() implements the Traub PSI 1,2 method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 232.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real d3fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real v
  real w
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )
    d3fx = f ( x, 3 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    v = d2fx / ( 2.0 * dfx )
    w = d3fx / ( 6.0 * dfx )
%
%  Set the increment.
%
    dx = - u / ( 1.0 - u * v + ( v**2 - w ) * u**2 )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function r8_muller ( x, x1, x2, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% r8_muller() implements Muller's method for real arithmetic.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, X1, X2: three distinct points that start the method.
%
%    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, X2: X is an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 and X2 are the previous estimates.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real a
  real abserr
  real b
  real c
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real q
  real term
  real x
  real x1
  real x2
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx2 = f ( x2, 0 )
  fx1 = f ( x1, 0 )
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    q = ( x - x1 ) / ( x1 - x2 )

    a = q * fx - q * ( 1.0 + q ) * fx1 + q**2 * fx2
    b = ( 2.0 * q + 1.0 ) * fx - ( 1.0 + q )**2 * fx1 + q**2 * fx2
    c = ( 1.0 + q ) * fx

    term = b**2 - 4.0 * a * c
    term = max ( term, 0.0 )
    term = sqrt ( term )
    if ( b < 0.0 )
      term = - term
    end if
%
%  Set the increment.
%
    dx = - ( x - x1 ) * 2.0 * c / ( b + term )
%
%  Remember current data for next step.
%
    x2 = x1
    fx2 = fx1
    x1 = x
    fx1 = fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function red_cap_phi_04 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% red_cap_phi_04() implements the Traub reduced capital PHI(0,4) method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 232.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real d3fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real v
  real w
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )
    d3fx = f ( x, 3 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    v = d2fx / ( 2.0 * dfx )
    w = d3fx / ( 6.0 * dfx )
    z = 1.0 - 4.0 * u * ( v - u * w )
%
%  Set the increment.
%
    if ( 0.0 < z )
      dx = - 2.0 * u / ( 1.0 + sqrt ( z ) )
    else
      dx = - fx / dfx
    end if
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function regula ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% regula() implements the Regula Falsi method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    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 an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real t
  real x
  real x1
  real x2
%
%  Initialization.
%
  ierror = 0
  k = 0

  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )

  if ( fx < 0.0  )
    t = x
    x = x1
    x1 = t
    t = fx
    fx = fx1
    fx1 = t
  end if
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

    if ( kmax < k )
      ierror = 2
      return
    end if
%
%  Set the increment.
%
    dx = - fx1 * ( x - x1 ) / ( fx - fx1 )
%
%  Update the iterate and function values.
%
    x2 = x1 + dx
    fx2 = f ( x2, 0 )

    if ( 0.0 <= fx2 )
      x = x2
      fx = fx2
    else
      x1 = x2
      fx1 = fx2
    end if

  end do

  return
end
function rhein1 ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% rhein1() implements the Rheinboldt bisection - secant method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Werner Rheinboldt,
%    Algorithms for finding zeros of a function,
%    UMAP Journal,
%    Volume 2, Number 1, 1981, pages 43 - 72.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    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 an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dx
  real em
  real, external :: f
  logical forced
  real fx
  real fx1
  real fx2
  integer i
  integer ierror
  integer k
  integer kmax
  real p
  real q
  real t
  real x
  real x1
  real x2
%
%  Initialization.
%
  i = 0
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )

  if ( abs ( fx1 ) < abs ( fx ) )
    t = x
    x = x1
    x1 = t
    t = fx
    fx = fx1
    fx1 = t
  end if

  x2 = x1
  fx2 = fx1
  t = 0.5 * abs ( x - x1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

    if ( kmax < k )
      ierror = 2
      return
    end if
%
%  Force ABS ( FX ) <= ABS ( FX1 ).
%
    if ( abs ( fx1 ) < abs ( fx ) )

      x2 = x
      x = x1
      x1 = x2

      fx2 = fx
      fx = fx1
      fx1 = fx2

    end if

    em = 0.5 * ( x1 - x )
%
%  Compute the numerator and denominator for secant step.
%
    p = ( x - x2 ) * fx
    q = fx2 - fx

    if ( p < 0.0 )
      p = - p
      q = - q
    end if
%
%  Save the old minimum residual point.
%
    x2 = x
    fx2 = fx
%
%  Test for forced bisection.
%
    i = i + 1
    forced = .false.

    if ( 3 < i )
      if ( t < 8.0 * abs ( em ) )
        forced = .true.
      else
        i = 0
        t = em
      end if
    end if
%
%  Set the increment.
%
    if ( forced )
      dx = em
    else if ( p <= abs ( q ) * abserr )
      dx = sign ( 1.0, em ) * abserr
    else if ( p < q * em )
      dx = p / q
    else
      dx = em
    end if
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
%
%  Preserve the change of sign interval.
%
    if ( sign ( 1.0, fx ) == sign ( 1.0, fx1 ) )
      x1 = x2
      fx1 = fx2
    end if

  end do

  return
end
function rhein2 ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% rhein2() implements the Rheinboldt bisection/secant/inverse quadratic method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Werner Rheinboldt,
%    Algorithms for Finding Zeros of a Function,
%    UMAP Journal,
%    Volume 2, Number 1, 1981, pages 43 - 72.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    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 an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dx
  real em
  real, external :: f
  logical forced
  real fx
  real fx1
  real fx2
  integer i
  integer ierror
  integer k
  integer kmax
  real piq
  real ps
  real qiq
  real qs
  real s
  real stpmin
  real t
  real u
  real v
  real w
  real x
  real x1
  real x2
%
%  Initialization.
%
  i = 0
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )

  if ( abs ( fx1 ) < abs ( fx ) )
    s = x
    x = x1
    x1 = s
    s = fx
    fx = fx1
    fx1 = s
  end if

  x2 = x1
  fx2 = fx1
  k = 0
  t = 0.5 * abs ( x - x1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    i = i + 1

    em = 0.5 * ( x1 - x )
%
%  Compute the numerator and denominator for secant step.
%
    if ( 2.0 * abs ( x2 - x ) < abs ( x1 - x ) )
      ps = ( x - x2 ) * fx
      qs = fx2 - fx
    else
      ps = ( x - x1 ) * fx
      qs = fx1 - fx
    end if

    if ( ps < 0.0 )
      ps = - ps
      qs = - qs
    end if
%
%  Compute the numerator and denominator for inverse quadratic.
%
    piq = 0.0
    qiq = 0.0

    if ( x1 /= x2 )
      u = fx / fx2
      v = fx2 / fx1
      w = fx / fx1
      piq = u * ( 2.0 * em * v * ( v - w ) - ( x - x2 ) * ( w - 1.0 ) )
      qiq = ( u - 1.0 ) * ( v - 1.0 ) * ( w - 1.0 )

      if ( 0.0 < piq )
        qiq =  - qiq
      end if

      piq = abs ( piq )
    end if
%
%  Save the old minimum residual point.
%
    x2 = x
    fx2 = fx

    stpmin = ( abs ( x ) + abs ( em ) + 1.0 ) * abserr
%
%  Choose bisection, secant or inverse quadratic step.
%
    forced = .false.

    if ( 3 < i )
      if ( t < 8.0 * abs ( em ) )
        forced = .true.
      else
        i = 0
        t = em
      end if
    end if
%
%  Set the increment.
%
    if ( forced )
      dx = em
    else if ( piq < 1.5 * em * qiq .and. ...
      abs ( qiq ) * stpmin < abs ( piq ) )
      dx = piq / qiq
    else if ( ps < qs * em .and. abs ( qs ) * stpmin < abs ( ps ) )
      dx = ps / qs
    else
      dx = em
    end if
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
%
%  Set the new X1 as either X1 or X2, depending on whether
%  F(X1) or F(X2) has the opposite sign from F(X).
%
    if ( sign ( 1.0, fx ) == sign ( 1.0, fx1 ) )
      x1 = x2
      fx1 = fx2
    end if
%
%  Force ABS ( FX ) <= ABS ( FX1 ).
%
    if ( abs ( fx1 ) < abs ( fx ) )
      s = x
      x = x1
      x1 = s
      s = fx
      fx = fx1
      fx1 = s
    end if

  end do

  return
end
function script_e2 ( x, m, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% script_e2() implements the Traub script E - 2 function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 234.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    integer M, the multiplicity of the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  integer m
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - real ( m, kind = rk ) * fx / dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function script_e3 ( x, m, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% script_e3() implements the Traub script E - 3 function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 235.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    integer M, the multiplicity of the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real a2
  real abserr
  real d2fx
  real dfx
  real dx
  real em
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  integer m
  real u
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  em = real ( m, kind = rk )
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    a2 = d2fx / ( 2.0 * dfx )
%
%  Set the increment.
%
    dx = - em * u * ( ( 3.0 - em ) / 2.0 + em * a2 * u )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function script_e4 ( x, m, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% script_e4() implements the Traub script E - 4 function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 235.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    integer M, the multiplicity of the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real a2
  real a3
  real abserr
  real d2fx
  real d3fx
  real dfx
  real dx
  real em
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  integer m
  real u
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  em = real ( m, kind = rk )
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )
    d3fx = f ( x, 3 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    a2 = d2fx / ( 2.0 * dfx )
    a3 = d3fx / ( 6.0 * dfx )
%
%  Set the increment.
%
    dx = - em * u * ( ( em**2 - 6.0 * em + 11.0 ) / 6.0 ...
       + em * ( 2.0 - em ) * a2 * u ...
       + em**2 * ( 2.0 * a2**2 - a3 ) * u**2 )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function secant ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% secant() implements the secant method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    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 an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dx
  real, external :: f
  real fx
  real fx1
  integer ierror
  integer k
  integer kmax
  real x
  real x1
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( fx - fx1 == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx * ( x - x1 ) / ( fx - fx1 )
%
%  Remember current data for next step.
%
    x1 = x
    fx1 = fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function secant_extended ( x, x1, x2, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% secant_extended() carries out the extended secant algorithm.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, X1, X3: three distinct points that start the method.
%
%    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, X3: X is an approximate root satisfying
%    wabs ( F(X) ) < ABSERR, and X1 and X2 are two previous estimates.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d1
  real d2
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real x
  real x1
  real x2
%
%  Initialization.
%
  ierror = 0
  k = 0

  fx1 = f ( x1, 0 )
  fx2 = f ( x2, 0 )
  d1 = ( fx1 - fx2 ) / ( x1 - x2 )

  if ( d1 == 0.0 )
    ierror = 3
    return
  end if

  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    d2 = d1

    if ( x == x1 )
      ierror = 3
      return
    end if

    d1 = ( fx - fx1 ) / ( x - x1 )

    if ( d1 == 0.0 )
      ierror = 3
      return
    end if

    if ( fx2 == fx1 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx / d1 + ( fx * fx1 / ( fx - fx2 ) ) ...
      * ( 1.0 / d1 - 1.0 / d2 )
%
%  Remember current data for next step.
%
    x2 = x1
    fx2 = fx1
    x1 = x
    fx1 = fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function star_e12 ( x, x1, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% star_e12() implements the Traub *E12 method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 234.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    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 an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d
  real dfx
  real dfx1
  real dx
  real, external :: f
  real fx
  real fx1
  integer ierror
  integer k
  integer kmax
  real u
  real x
  real x1
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )

  fx1 = f ( x1, 0 )
  dfx1 = f ( x1, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx

    if ( x == x1 )
      ierror = 3
      return
    end if

    d = ( fx - fx1 ) / ( x - x1 )
    z = 2.0 * dfx + dfx1 - 3.0 * d
%
%  Set the increment.
%
    dx = - u - u**2 * z / ( dfx * ( x - x1 ) )
%
%  Remember current data for next step.
%
    x1 = x
    fx1 = fx
    dfx1 = dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function star_e21 ( x, x1, x2, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% star_e21() implements the Traub *E21 method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 234.
%
%  Input:
%
%    real X, X1, X2: three distinct points that start the method.
%
%    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, X2: X is an approximate root satisfying
%    abs ( F(X) ) < ABSERR, X1 and X2 are the previous estimates.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d
  real d1
  real d2
  real dx
  real, external :: f
  real fx
  real fx1
  real fx2
  integer ierror
  integer k
  integer kmax
  real x
  real x1
  real x2
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  fx1 = f ( x1, 0 )
  fx2 = f ( x2, 0 )
  d1 = ( fx1 - fx2 ) / ( x1 - x2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    d2 = d1

    if ( x == x1 )
      ierror = 3
      return
    end if

    d1 = ( fx - fx1 ) / ( x - x1 )

    if ( x == x2 )
      ierror = 3
      return
    end if

    d = ( fx - fx2 ) / ( x - x2 )

    if ( d1 + d - d2 == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx / ( d1 + d - d2 )
%
%  Remember current data for next step.
%
    x2 = x1
    fx2 = fx1

    x1 = x
    fx1 = fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function steffenson ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% steffenson() implements Steffenson's method.
%
%  Discussion:
%
%    Steffenson's method is of order 2.
%
%    x(n+1) = x(n) - f(x(n)) / g(x(n))
%
%    where
%
%    g(x) = ( f(x+f(x)) - f(x) ) / f(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 178.
%
%  Input:
%
%    real X: the point that starts the method.
%
%    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: an approximate root satisfying abs ( F(X) ) < ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dx
  real, external :: f
  real fx
  real gx
  integer ierror
  integer k
  integer kmax
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0

  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    gx = ( f ( x + fx, 0 ) - fx ) / fx

    if ( gx == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx / gx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function stirling ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% stirling() implements Stirling's method.
%
%  Discussion:
%
%    Stirling's method is of order 2.
%
%    x(n+1) = x(n) - f(x(n)) / f'( x(n) - f(x(n) )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: the point that starts the method.
%
%    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: an approximate root satisfying abs ( F(X) ) < ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dx
  real, external :: f
  real fx
  real gx
  integer ierror
  integer k
  integer kmax
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0

  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    gx = f ( x - fx, 1 )

    if ( gx == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx / gx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function t14 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% t14() implements the Traub fourteenth function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 237.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfxu
  real dfz
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    dfxu = f ( x - ( fx / dfx ), 1 )

    if ( dfxu == 0.0 )
      ierror = 3
      return
    end if

    z = x - 0.25 * ( fx / dfx + fx / dfxu  )
    dfz = f ( z, 1 )
%
%  Set the increment.
%
    dx = - ( fx / dfx + ( fx / dfxu ) + 4.0 * ( fx / dfz ) ) / 6.0
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function t15 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% t15() implements the Traub fifteenth function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 238.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfxu
  real dfz
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    dfxu = f ( x - ( fx / dfx ), 1 )

    if ( dfxu == 0.0 )
      ierror = 3
      return
    end if

    z = x - 2.0 * ( 2.0 * u + fx / dfxu ) / 9.0
    dfz = f ( z, 1 )

    if ( dfz == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - 0.25 * ( fx / dfx + 3.0 * fx / dfz )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function t16 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% t16() implements the Traub sixteenth function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 238.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfxu
  real dfz
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    dfxu = f ( x - ( fx / ( 3.0 * dfx ) ), 1 )

    if ( dfxu == 0.0 )
      ierror = 3
      return
    end if

    z = x - 2.0 * fx / ( 3.0 * dfxu )
    dfz = f ( z, 1 )

    if ( dfz == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - 0.25 * ( fx / dfx + 3.0 * fx / dfz )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function t_family1 ( x, c, d, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% t_family1() implements the Traub first family of iterations.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, pages 236 - 237.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    real C, D, ???
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real c
  real d
  real dfx
  real dfz
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    z = x - d * fx / dfx
    dfz = f ( z, 1 )

    if ( dfz == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - ( c / dfx + ( 1.0 - c ) / dfz ) * fx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function t_family2 ( x, a, b, c, d, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% t_family2() implements the Traub second family of iterations.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 236 - 237.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    real A, B, C, D, ???
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real a
  real abserr
  real b
  real c
  real d
  real dfx
  real dfz
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    z = x - d * fx / dfx
    dfz = f ( z, 1 )
%
%  Set the increment.
%
    dx = - fx * ( b * dfx - c * dfz ) / ( a * dfx**2 )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function te11f ( x, x1, m, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% te11f() implements the Traub *E - 1,1(f) function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 235.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    integer M, the multiplicity of the root.
%
%    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 an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d
  real dx
  real expon
  real, external :: f
  real fnorm
  real fx
  real fx1
  real gx
  real gx1
  integer ierror
  integer k
  integer kmax
  integer m
  real x
  real x1
%
%  Initialization.
%
  ierror = 0
  k = 0
  expon = 1.0 / real ( m, kind = rk )
  fx1 = f ( x1, 0 )

  fnorm = abs ( fx1 )
  gx1 = sign ( 1.0, fx1 ) * fnorm**expon
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    fnorm = abs ( fx )
    gx = sign ( 1.0, fx ) * fnorm**expon

    if ( x == x1 )
      ierror = 3
      return
    end if

    d = ( gx - gx1 ) / ( x - x1 )
%
%  Set the increment.
%
    dx = - gx / d
%
%  Remember current data for next step.
%
    gx1 = gx
    x1 = x
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function te2u ( x, em, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% te2u() implements the Traub E - 2(u) function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 235.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    real EM, an estimate of the multiplicity of the root.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real em
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    if ( dfx**2 - fx * d2fx == 0.0 )
      ierror = 3
      return
    end if

    em = dfx**2 / ( dfx**2 - fx * d2fx )
    em = max ( em, 1.0 )
%
%  Set the increment.
%
    dx = - em * fx / dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function tphi1u ( x, em, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% tphi1u() implements the Traub phi - 1,1(u) function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 235.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    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:  an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    real EM, an estimate of the multiplicity of the root.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real em
  real, external :: f
  real fu
  real fu1
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real u1
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  u1 = f ( x, 0 ) / f ( x, 1 )
  x = x - u1
  fx = f ( x, 0 )
  dfx = f ( x, 1 )

  fu1 = f ( u1, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    fu = f ( u, 0 )

    em = ( u - u1 ) / ( fu - fu1 )
    em = max ( em, 1.0 )
%
%  Set the increment
%
    dx = - em * fx / dfx
%
%  Remember current data for next step.
%
    u1 = u
    fu1 = fu
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function traub1 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% traub1() implements the Traub first method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 236.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfz
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    z = x - fx / dfx
    dfz = f ( z, 1 )

    if ( dfz == 0.0 )
      ierror = 3
      return
    end if
%
%  Set the increment.
%
    dx = - fx / dfz
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function traub4 ( x, nsub, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% traub4() implements the Traub fourth method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 236.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    integer NSUB, the number of steps in the sub-iteration.
%    The derivatives are only evaluated before the first step,
%    and after every subsequent NSUB steps.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  integer nsub
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( mod ( k - 1, nsub ) == 0 )

      dfx = f ( x, 1 )
      d2fx = f ( x, 2 )

      if ( dfx == 0.0 )
        ierror = 3
        return
      end if

      if ( dfx - d2fx * fx / dfx == 0.0 )
        ierror = 3
        return
      end if

    end if
%
%  Set the increment.
%
    dx = - fx / ( dfx - d2fx * fx / dfx )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function traub8 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% traub8() implements the Traub eighth function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfxu
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    dfxu = f ( x - 2.0 * fx / ( 3.0 * dfx ), 1 )
%
%  Set the increment.
%
    dx = - 4.0 * fx / ( dfx + 3.0 * dfxu )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function traub9 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% traub9() implements the Traub ninth method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 237.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  real fxu
  integer ierror
  integer k
  integer kmax
  real u
  real x
  real xu
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    xu = x - fx / dfx
    fxu = f ( xu, 0 )

    if ( 2.0 * fxu - fx == 0.0 )
      ierror = 4
      return
    end if
%
%  Set the increment.
%
    dx = - fx / dfx + u * fxu / ( 2.0 * fxu - fx )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function traub_ostrowski ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% traub_ostrowski() implements the Traub-Ostrowksi method.
%
%  Discussion:
%
%    The Traub-Ostrowski method is of order 4.
%
%    x(n+1) = x(n) - u(x(n)) * ( f ( x(n) - u(x(n)) - f(x(n)) ) /
%      ( 2 * f ( x(n) - u(x(n)) - f(x(n)) )
%
%    where
%
%    u(x) = f(x) / f'(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 184.
%
%  Input:
%
%    real X: the point that starts the method.
%
%    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: an approximate root satisfying abs ( F(X) ) < ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  real fx1
  integer ierror
  integer k
  integer kmax
  real ux
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0

  fx = f ( x, 0 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    dfx = f ( x, 1 )

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    ux = f ( x, 0 ) / dfx

    fx1 = f ( x - ux, 0 )

    if ( 2.0 * fx1 - fx == 0.0 )
      ierror = 4
      return
    end if
%
%  Set the increment.
%
    dx = - ux * ( fx1 - fx ) / ( 2.0 * fx1 - fx )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )

  end do

  return
end
function tt1f ( x, a, rho, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% tt1f() implements the Traub type 1 functions 10 and 11.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 237.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    real A, ???
%
%    real RHO, ???
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real a
  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  real fxu
  real fz
  integer ierror
  integer k
  integer kmax
  real rho
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0

  if ( rho == 0.0 )
    ierror = 3
    return
  end if

  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    z = x + rho * fx / dfx
    fxu = f ( z, 0 )

    z = x - fxu / ( rho**2 * dfx )
    fz = f ( z, 0 )
%
%  Set the increment.
%
    dx = - ( a * fz + fxu / rho**2 ) / dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function tthip ( x, em, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% tthip() implements the Traub third function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joseph Traub,
%    Iterative Methods for the Solution of Equations,
%    ISBN: 0828403120,
%    LC: QA297.T7.
%    Prentice Hall, 1964, page 235.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    real EM, an estimate of the multiplicity of the root.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real em
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    if ( fx == 0.0 )
      em = 1.0
    else if ( log ( abs ( fx / dfx ) ) == 0.0 )
      em = 1.0
    else
      em = log ( abs ( fx ) ) / log ( abs ( fx / dfx ) )
    end if

    em = max ( em, 1.0 )
%
%  Set the increment.
%
    dx = - em * fx / dfx
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function vandev ( x, em, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% vandev() implements the Van de Vel iteration.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Hugo vandeVel,
%    A method for computing a root of a single nonlinear equation,
%    including its multiplicity,
%    Computing,
%    Volume 14, 1975, pages 167 - 171.
%
%  Input:
%
%    real X: an estimate for the root.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    real EM, an estimate of the multiplicity of the root.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dfz
  real dx
  real em
  real em1
  real, external :: f
  real fx
  real fz
  integer ierror
  integer k
  integer kmax
  real u
  real u1
  real x
  real z
%
%  Initialization.
%
  ierror = 0
  k = 0
  em1 = 1.0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx
    z = x - em1 * fx / dfx
    fz = f ( z, 0 )
    dfz = f ( z, 1 )

    if ( dfz == 0.0 )
      ierror = 3
      return
    end if

    u1 = fz / dfz

    if ( u - u1 /= 0.0 )
      em = em1 * u / ( u - u1 )
    else
      em = em1
    end if

    em = max ( em, 1.0 )
%
%  Set the increment.
%
    dx = - em1 * fx / dfx - em * fz / dfz
%
%  Remember current data for next step.
%
    em1 = em
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function vandev2 ( x, em, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% vandev2() implements the improved Van de Vel iteration.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard King,
%    Improving the van de Vel root-finding method,
%    Computing,
%    Volume 30, 1983, pages 373 - 378.
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which
%    abs ( F(X) ) <= ABSERR.
%
%    real EM, an estimate of the multiplicity of the root.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real dfx
  real dx
  real, external :: f
  real fx
  real em
  real em1
  integer ierror
  integer k
  integer kmax
  real u
  real u1
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  em1 = 1.0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )

  if ( dfx == 0.0 )
    ierror = 3
    return
  end if

  u1 = fx / dfx

  x = x - fx / dfx
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = fx / dfx

    if ( u1 - u /= 0.0 )
      em = em1 * u1 / ( u1 - u )
    else
      em = em1
    end if

    em = max ( em, 1.0 )
%
%  Set the increment.
%
    dx = - em * fx / dfx
%
%  Remember current data for next step.
%
    em1 = em
    u1 = u
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )

  end do

  return
end
function whittaker ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% whittaker() implements the convex acceleeration of Whittaker's method.
%
%  Discussion:
%
%    x(n+1) = x(n) - ( f(x(n)) / f'(x(n)) ) * ( 1 - 0.5 * L(x(n)) )
%
%    where
%
%      L(x) = f(x) * f''(x) / f'(x)**2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = d2fx * fx / dfx**2
%
%  Set the increment.
%
    dx = - ( fx / dfx ) * ( 1.0 - 0.5 * u )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end
function whittaker2 ( x, abserr, kmax, f, ierror, k )

%*****************************************************************************80
%
%% whittaker2() implements the double convex acceleeration of Whittaker's method.
%
%  Discussion:
%
%    x(n+1) = x(n) - (1/4) * ( f(x(n)) / f'(x(n)) ) *
%      ( 2 - L(x(n)) + ( 2 + L(x(n)) ) / ( 1 - L(x(n)) + 0.5 * L(x(n))**2 ) )
%
%    where
%
%      L(x) = f(x) * f''(x) / f'(x)**2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%
  implicit none

  integer, parameter :: rk = kind ( 1.0 )

  real abserr
  real d2fx
  real dfx
  real dx
  real, external :: f
  real fx
  integer ierror
  integer k
  integer kmax
  real u
  real x
%
%  Initialization.
%
  ierror = 0
  k = 0
  fx = f ( x, 0 )
  dfx = f ( x, 1 )
  d2fx = f ( x, 2 )
%
%  Iteration loop:
%
  do
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      exit
    end if

    k = k + 1

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

    if ( dfx == 0.0 )
      ierror = 3
      return
    end if

    u = d2fx * fx / dfx**2

    if ( 1.0 - u + 0.5 * u * u  == 0.0 )
      ierror = 4
      return
    end if
%
%  Set the increment.
%
    dx = - 0.25 * ( fx / dfx ) ...
      * ( 2.0 - u + ( 2.0 + u ) / ( 1.0 - u + 0.5 * u * u ) )
%
%  Update the iterate and function values.
%
    x = x + dx
    fx = f ( x, 0 )
    dfx = f ( x, 1 )
    d2fx = f ( x, 2 )

  end do

  return
end