subroutine laguerre ( x0, degree, abserr, kmax, f, x, ierror, k )

c*********************************************************************72
c
cc laguerre() implements the Laguerre rootfinding method for polynomials.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 March 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Eldon Hansen, Merrell Patrick,
c    A Family of Root Finding Methods,
c    Numerische Mathematik,
c    Volume 27, 1977, pages 257 - 269.
c
c  Input:
c
c    double precision X0.
c    An initial estimate for the root of the equation.
c
c    integer degree: the polynomial degree of the function, at least 2.
c
c    double precision ABSERR: an error tolerance.
c
c    integer KMAX: the maximum number of iterations allowed.
c
c    double precision external F: the name of the routine that
c    evaluates the function or its derivatives, of the form
c      function f ( x, ider )
c
c  Output:
c
c    double precision X: the estimated solution, if IERROR=0.
c
c    integer IERROR: error indicator.
c    0, no error occurred.
c    nonzero, an error occurred, and the iteration was halted.
c
c    integer K: the number of steps taken.
c
      implicit none

      double precision abserr
      double precision beta
      double precision bot
      double precision d2fx
      integer degree
      double precision dfx
      double precision dx
      double precision, external :: f
      double precision fx
      integer ierror
      integer k
      integer kmax
      double precision x
      double precision x0
      double precision z
c
c  Initialization.
c
      ierror = 0
      x = x0

      beta = 1.0D+00 / dble ( degree - 1 )

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

        k = k + 1

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

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

        bot = beta * dfx + sqrt ( z )

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

      end do

      return
      end