# include <math.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>

# include "zero_laguerre.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  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.

    integer degree: the polynomial degree, at least 2.

    double ABSERR: an error tolerance.

    integer 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;
}