# include <complex.h>
# include <math.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>
#
# include "zero_muller.h"
#
/******************************************************************************/

void zero_muller ( double complex func ( double complex x ), double fatol,
  int itmax, double complex x1, double complex x2, double complex x3, 
  double xatol, double xrtol, double complex *xnew, double complex *fxnew )

/******************************************************************************/
/*
  Purpose:

    zero_muller() carries out Muller's method, using complex arithmetic.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 March 2024

  Author:

    John Burkardt

  Reference:

    Gisela Engeln-Muellges, Frank Uhlig,
    Numerical Algorithms with C,
    Springer, 1996,
    ISBN: 3-540-60530-4,
    LC: QA297.E56213.

  Input:

    double complex func ( double complex x ): evaluates the function.
 
    double FATOL, the absolute error tolerance for F(X).

    integer ITMAX, the maximum number of steps allowed.

    double complex X1, X2, X3, three distinct points to start the
    iteration.

    double XATOL, XRTOL, absolute and relative
    error tolerances for the root.

  Output:

    double complex *xnew, the estimated root.

    double complex *fxnew, the value of the function at xnew.
*/
{
  double complex a;
  double complex b;
  double complex c;
  double complex c8_temp;
  double complex discrm;
  double complex fminus;
  double complex fplus;
  double complex fxmid;
  double complex fxold;
  int iterate;
  double x_ave;
  double complex x_inc;
  double complex xmid;
  double complex xminus;
  double complex xold;
  double complex xplus;
 
  *xnew = x1;
  xmid = x2;
  xold = x3;

  *fxnew = func ( *xnew );
  fxmid = func ( xmid );
  fxold = func ( xold );

  printf ( "\n" );
  printf ( "zero_muller():\n" );
  printf ( "  Muller's root-finding method (complex root version)\n" );
  printf ( "\n" );
  printf ( "  Iteration     x_real              x_imag             ||fx||           ||disc||\n" );
  printf ( "\n" );

  iterate = -2;
  printf ( "%6d%20.10f%20.10f%20.10f\n", 
    iterate, creal ( xold ), cimag ( xold ), cabs ( fxold ) );
  iterate = -1;
  printf ( "%6d%20.10f%20.10f%20.10f\n", 
    iterate, creal ( xmid ), cimag ( xmid ), cabs ( fxmid ) );
  iterate = 0;
  printf ( "%6d%20.10f%20.10f%20.10f\n", 
    iterate, creal ( *xnew ), cimag ( *xnew ), cabs ( *fxnew ) );

  if ( cabs ( *fxnew ) < fatol )
  {
    printf ( "\n" );
    printf ( "zero_muller():\n" );
    printf ( "  |F(X)| is below the tolerance.\n" );
    return;
  }

  while ( true )
  {
/*
  You may need to swap (XMID,FXMID) and (XNEW,FXNEW).
*/
    if ( cabs ( fxmid ) <= cabs ( *fxnew ) )
    {
      c8_temp = *xnew;
      *xnew = xmid;
      xmid = c8_temp;

      c8_temp = *fxnew;
      *fxnew = fxmid;
      fxmid = c8_temp;
    }

    iterate = iterate + 1;

    if ( itmax < iterate )
    {
      printf ( "\n" );
      printf ( "zero_muller(): Warning!\n" );
      printf ( "  Maximum number of steps taken.\n" );
      break;
    }

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

    b = ( ( xold - *xnew ) * ( xold - *xnew ) * ( fxmid - *fxnew ) 
        - ( xmid - *xnew ) * ( xmid - *xnew ) * ( 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 * b - 4.0 * a * c;

    if ( a == 0.0 )
    {
      printf ( "\n" );
      printf ( "zero_muller(): Warning!\n" );
      printf ( "  The algorithm has broken down.\n" );
      printf ( "  The quadratic coefficient A is zero.\n" );
      break;
    }

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

    fplus = func ( xplus );

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

    fminus = func ( xminus );
/*
  Choose the root with smallest function value.
*/
    if ( cabs ( fminus ) < cabs ( fplus ) )
    {
      *xnew = xminus;
    }
    else
    {
      *xnew = xplus;
    }

    fxold = fxmid;
    fxmid = *fxnew;
    *fxnew = func ( *xnew );

    printf ( "%6d%20.10f%20.10f%20.10f%20.10f\n", 
      iterate, creal ( *xnew ), cimag ( *xnew ), cabs ( *fxnew ), cabs ( discrm ) );
/*
  Check for convergence.
*/
    x_ave = cabs ( *xnew + xmid + xold ) / 3.0;
    x_inc = *xnew - xmid;

    if ( cabs ( x_inc ) <= xatol )
    {
      printf ( "\n" );
      printf ( "zero_muller():\n" );
      printf ( "  Absolute convergence of the X increment.\n" );
      break;
    }

    if ( cabs ( x_inc ) <= xrtol * x_ave )
    {
      printf ( "\n" );
      printf ( "zero_muller():\n" );
      printf ( "  Relative convergence of the X increment.\n" );
      break;
    }

    if ( cabs ( *fxnew ) <= fatol )
    {
      printf ( "\n" );
      printf ( "zero_muller():\n" );
      printf ( "  Absolute convergence of |F(X)|.\n" );
      break;
    }

  }

  return;
}