# include <cmath>
# include <cstdlib>

# include "zero_laguerre.hpp"

//****************************************************************************80

void zero_laguerre ( double x0, int degree, double abserr, int kmax, 
  double f ( double x, int ider ), double &x, int &ierror, int &k )

//****************************************************************************80
//
//  Purpose:
//
//    zero_laguerre() implements the Laguerre rootfinding method for polynomials.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    30 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:
//
//    double x0: initial estimate for a root.
//
//    int degree: the polynomial degree, at least 2.
//
//    double ABSERR: an error tolerance.
//
//    int KMAX: the maximum number of iterations allowed.
//
//    double f ( double x, int ider ): evaluates the function or its first 
//    two derivatives.
//
//  Output:
//
//    double X: the estimated solution, if IERROR=0.
//
//    int IERROR: error indicator.
//    0, no error occurred.
//    nonzero, an error occurred, and the iteration was halted.
//
//    int K: the number of steps taken.
//
{
  double beta;
  double bot;
  double d2fx;
  double dfx;
  double dx;
  double fx;
  double z;
//
//  Initialization.
//
  x = x0;
  ierror = 0;
  k = 0;

  beta = 1.0 / ( double ) ( 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 ( fabs ( fx ) <= abserr )
    {
      break;
    }

    k = k + 1;

    if ( kmax < k )
    {
      ierror = 2;
      break;
    }

    z = dfx * dfx - ( beta + 1.0 ) * fx * d2fx;
    z = fmax ( z, 0.0 );

    bot = beta * dfx + sqrt ( z );

    if ( bot == 0.0 )
    {
      ierror = 3;
      break;
    }
//
//  Set the increment.
//
    dx = - ( beta + 1.0 ) * fx / bot;
//
//  Update the iterate and function values.
//
    x = x + dx;
  }

  return;
}