/* Merhly method */
/* maehly_ehrich_aberth.c */
#include <stdio.h>
#include <math.h>
#include <complex.h>
#define NUMAX 20
#define EPS 1.0e-12 
#define N 5 /* f(z) の次数 */
#define r 3.875  /* 零点を含む重心中心の円の半径の上限 */
double complex c = 1.0; /* f(z)の最高次の係数 */
double complex f(double complex z);
double complex df(double complex z);

int main(void)
{
  int i,j,nu; 
  double res; /* 残差 */
  double complex z[N],z0[N],w;
  

  /* 初期値 */
  for(j = 0; j < N;j++){
    z[j] =0.6 + r * cexp(I * M_PI * (2.0 * j / N + 1.0 / (2 * N)));
  }


  res = 0.0;
  for(j = 0;j < N;j++)
    if(res < cabs(f(z[j]))) res = cabs(f(z[j]));

  nu = 0;
  printf("%2d %10.3e\n",nu,res);
  for(j = 0;j < N;j++)
    printf("  %d %20.15f %20.15f\n",j,creal(z[j]),cimag(z[j])); 

  while(res > EPS && nu < NUMAX){
    nu++;
    for(j = 0;j < N;j++){
      z0[j] = z[j];
    }
    for(j = 0;j < N;j++){
      w = 0.0;
      for(i = 0;i < N;i++){
	if (j == i) continue;
	w += 1.0 / (z0[j] - z0[i]);
      }
      z[j] = z0[j] - 1.0 / (df(z0[j])/f(z0[j]) - w);
    }
  
    res = 0.0;
    for(j = 0;j < N;j++)
      if(cabs(f(z[j])) > res) res = cabs(f(z[j]));

    printf("%2d %10.3e\n",nu,res);
 
    for(j = 0;j < N;j++)
      printf("  %d %20.15f %20.15f\n",j,creal(z[j]),cimag(z[j])); 
  } /* end while */

  return 0;
}


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

  w = -50 + z * (80 + z * (-37 + z * (9 + z * (-3 + z))));
  
  return w;
}

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

  w = 80 + z * (-74 + z * (27 + z * (-12 + 5 * z)));
  
  return w;
}