function [ x, ierror, k ] = zero_laguerre ( x0, degree, abserr, kmax, f )

%*****************************************************************************80
%
%% zero_laguerre() implements the Laguerre rootfinding method for polynomials.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 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.
%
%    integer degree: the polynomial degree of the function, at least 2.
%
%    real abserr: an error tolerance.
%
%    integer kmax: the maximum number of iterations allowed.
%
%    function fx = f ( x, ider ): evaluates the function or
%    first or second derivative (ider = 0, 1, 2)
%
%  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.
%

%
%  Initialization.

  x = x0;
  ierror = 0;
  k = 0;

  beta = 1.0 / ( degree - 1 );
%
%  Iteration loop:
%
  while ( true )

    fx = f ( x, 0 );
    dfx = f ( x, 1 );
    d2fx = f ( x, 2 );
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      break
    end

    k = k + 1;

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

    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
%
%  Set the increment.
%
    dx = - ( beta + 1.0 ) * fx / bot;
%
%  Update the iterate and function values.
%
    x = x + dx;

  end

  return
end
