/* 複素Laguerre法 */
/* complex_laguerre.c */
#include <stdio.h>
#include <complex.h>
#define N 3
#define NMAX 10
#define EPS 1.0e-14 
double complex f(double complex z);  /* f(z) */
double complex df(double complex z);  /* f'(z) */
double complex ddf(double complex z);  /* f''(z) */

int main(void)
{
  int n;
  double absfz,absfw; 
  double complex z,fz,dfz,ddfz;
  double complex w,fw,dfw,ddfw;  

  z = -1.0 + 1.0 * I; /* 初期値 */
  w = -1.0 + 1.0 * I; /* 初期値 */  
  n = 0;
  absfz = cabs(f(z));
  absfw = cabs(f(w));  
  printf("%2s %20s %19s %10s  %20s %19s %10s\n","n","Re(z_n)  ","Im(z_n)  ","|f(z_n)|","Re(w_n)  ","Im(w_n)  ","|f(w_n)|");
printf("%2d %20.15f %19.15f %10.3e  %20.15f %19.15f %10.3e\n",n,creal(z),cimag(z),absfz,creal(w),cimag(w),absfw); 
 
    while((absfz > EPS || absfw > EPS) && n < NMAX)
    {
      n++;
      fz = f(z);
      dfz = df(z);
      ddfz = ddf(z);      
      z = z - N * fz /(dfz+csqrt((N-1)*(N-1)*dfz*dfz-N*(N-1)*fz*ddfz));
      fw = f(w);
      dfw = df(w);
      ddfw = ddf(w);      
      w = w - N * fw /(dfw-csqrt((N-1)*(N-1)*dfw*dfw-N*(N-1)*fw*ddfw));
      absfz=cabs(f(z));
      absfw=cabs(f(w));
      printf("%2d %20.15f %19.15f %10.3e  %20.15f %19.15f %10.3e\n",n,creal(z),cimag(z),absfz,creal(w),cimag(w),absfw); 
    }
  
  return 0;
}

double complex f(double complex z)
{
  double complex w;

  w = (z * z - 2.0) * z - 5.0;

  return w;
}


double complex df(double complex z)
{
  double complex w;

  w = 3.0 * z * z - 2.0;

  return w;
}

double complex ddf(double complex z)
{
  double complex w;

  w = 6.0 * z;

  return w;
}