# include <cmath>
# include <complex>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "zero_muller.hpp"

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

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

//****************************************************************************80
//
//  Purpose:
//
//    zero_muller() carries out Muller's method, using complex arithmetic.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 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:
//
//    complex <double> func ( complex <double> x ): evaluates the function.
//
//    double FATOL, the absolute error tolerance for F(X).
//
//    integer ITMAX, the maximum number of steps allowed.
//
//    complex <double> X1, X2, X3, three distinct points to start the
//    iteration.
//
//    double XATOL, XRTOL, absolute and relative
//    error tolerances for the root.
//
//  Output:
//
//    complex <double> &xnew, the estimated root.
//
//    complex <double> &fxnew, the value of the function at xnew.
//
{
  complex <double> a;
  complex <double> b;
  complex <double> c;
  complex <double> c8_temp;
  complex <double> discrm;
  complex <double> fminus;
  complex <double> fplus;
  complex <double> fxmid;
  complex <double> fxold;
  int iterate;
  double x_ave;
  complex <double> x_inc;
  complex <double> xmid;
  complex <double> xminus;
  complex <double> xold;
  complex <double> xplus;
 
  xnew = x1;
  xmid = x2;
  xold = x3;

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

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

  iterate = -2;
  cout << "  " << iterate
       << "  " << real ( xold )
       << "  " << imag ( xold )
       << "  " << abs ( fxold ) << "\n";
  iterate = -1;
  cout << "  " << iterate
       << "  " << real ( xmid )
       << "  " << imag ( xmid )
       << "  " << abs ( fxmid ) << "\n";
  iterate = 0;
  cout << "  " << iterate
       << "  " << real ( xnew )
       << "  " << imag ( xnew )
       << "  " << abs ( fxnew ) << "\n";

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

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

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

    iterate = iterate + 1;

    if ( itmax < iterate )
    {
      cout << "\n";
      cout << "zero_muller(): Warning!\n";
      cout << "  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 )
    {
      cout << "\n";
      cout << "zero_muller(): Warning!\n";
      cout << "  The algorithm has broken down.\n";
      cout << "  The quadratic coefficient A is zero.\n";
      break;
    }

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

    fplus = func ( xplus );

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

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

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

    cout << "  " << iterate
         << "  " << real ( xnew )
         << "  " << imag ( xnew )
         << "  " << abs ( fxnew )
         << "  " << abs ( discrm ) << "\n";
//
//  Check for convergence.
//
    x_ave = abs ( xnew + xmid + xold ) / 3.0;
    x_inc = xnew - xmid;

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

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

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

  }

  return;
}