// DK.C --- Durand Kerner 法で複素係数代数方程式を解く

#include <iostream.h>
#include <iomanip.h>    // setprecision() のため
#include <math.h>
#include <complex.h>    // GNU ANSI C++ Library の comples

#define MAXN 100

complex<double> polynomial(int, complex<double> *, complex<double>);
complex<double> bunbo(int, complex<double> *, int);

main()
{
  int i, k, n;
  complex<double> a[MAXN+1], x[MAXN], newx[MAXN], dx[MAXN];
  complex<double> g, I(0,1);
  double r0, max, pi;

  pi = 4 * atan(1.0);

  // 表示は 16 桁とする
  cout << setprecision(16);
  // 方程式の入力
  cout << "次数 n=";
  cin >> n;
  if (n > MAXN || n <= 0) {
    cerr << "次数は" << MAXN << "以下の自然数として下さい。" << endl;
    exit(0);
  }
  for (i = 0; i <= n; i++) {
    cout << (n-i) << "次の係数を入力してください: ";
    cin >> a[i];
  }
  cout << "monic にします。" << endl;
  for (i = 1; i <= n; i++)
    a[i] /= a[0];
  a[0] = 1;
  cout << "修正した係数" << endl;
  for (i = 0; i <= n; i++)
    cout << "a[" << i << "]=" << a[i] << endl;

  // Aberth の初期値
  g = - a[1] / n;
  cout << "根の重心" << g << endl;
  max = 0;
  for (i = 1; i <= n; i++)
    if (abs(a[i]) > max)
      max = abs(a[i]);
  cout << "max|a_i|=" << max << endl;
  r0 = abs(g) + 1 + max;
  cout << "根は重心中心で、半径 r0=" << r0 << "の円盤内にある" << endl;

  for (i = 0; i < n; i++) {
    double theta;
    theta = 2 * i * pi / n + pi / (2 * n);
    x[i] = g + r0 * exp(I * theta);
    cout << x[i] << endl;
  }

  /* 反復 */
  for (k = 1; k <= 1000; k++) {
    double error;
    cout << "第" << k << "反復" << endl;
    for (i = 0; i < n; i++) {
      dx[i] = polynomial(n, a, x[i]) / bunbo(n, x, i);
      newx[i] = x[i] - dx[i];
    }
    for (i = 0; i < n; i++) {
      x[i] = newx[i];
      cout << x[i] << endl;
    }
    error = 0.0;
    for (i = 0; i < n; i++)
      error += abs(dx[i]);
    if (error < 1e-12)
      break;
  }
  return 0;
}

// Horner 法 (いわゆる組立除法で多項式の値を計算する)
complex<double> polynomial(int n,
			   complex<double> *a,
			   complex<double> z)
{
  int i;
  complex<double> w;
  w = a[0];
  for (i = 1; i <= n; i++)
    w = (w * z + a[i]);
  return w;
}

complex<double> bunbo(int n,
		      complex<double> *z,
		      int i)
{
  int j;
  complex<double> w(1,0);
  for (j = 0; j < n; j++)
    if (j != i)
      w *= (z[i] - z[j]);
  return w;
}