// DK3.C -- DKA 法で代数方程式の根を求める
//
//   数学科 6701 号室のワークステーション環境でのコンパイル
//     (グラフィックス・ライブラリィ GLSC を使っているので
//      情報科学センターではコンパイルできない。悪しからず。)
//     g++ -o DK3 DK3.C -I/usr/local/include -lglscd -lX11 -lsocket
//   あるいは
//     ccmg+ DK3.C
//
//  Win32 環境下の glscwin でも (若干の修正の下に) 動作させられることが
//  分かったので、画面に表示するメッセージを変更しました (2000/10/18)。

#include <iostream.h>
#include <iomanip.h>   // setprecision() のため
#include <math.h>
#include <complex.h>   // complex<double> のため

// GLSC のヘッダーファイルをインクルード
extern "C" {
#define G_DOUBLE
#include <glsc.h>
}

// プロトタイプ宣言
complex<double> polynomial(int, complex<double> *, complex<double>);
complex<double> bunbo(int, complex<double> *, int);

// 解きたい代数方程式の次数の最高値
#define MAXN 100

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, R, max, pi;

  // 数学定数 (円周率) の準備
  pi = 4 * atan(1.0);

  // 表\示の桁数を 16 桁にする
  cout << setprecision(16);

  // 方程式の入力
  cout << "次数 n を入力してください (1≦n≦" << MAXN << "): ";
  cin >> n;
  if (n > MAXN || n <= 0) {
    cerr << "次数は" << MAXN << "以下の自然数として下さい。" << endl;
    exit(0);
  }
  for (i = 0; i <= n; i++) {
    cout << (n-i) << "次の係数を入力してください: ";
    cin >> a[i];
  }

  // 多項式を最高次の係数で割って monic にする
  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;

  cout << "円の半径 (分からなければ上の値を指定してください): ";
  R = r0;
  cin >> r0;
  if (r0 > R)
    R = r0;
  cout << "図は根の重心を中心として半径 " << R
#ifdef __CYGWIN__
       << "の円が表\示できるようにします。" << endl;
#else
       << "の円が表示できるようにします。" << endl;
#endif

  // Aberth の初期値
  cout << "初期値" << 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;
  }

  // グラフィックス・ライブラリィ GLSC の初期化
  g_init("DKGRAPH", 120.0, 120.0);
  g_device(G_BOTH);
  // 座標系の指定
  g_def_scale(0,
	      real(g) - 1.1 * R, real(g) + 1.1 * R,
	      imag(g) - 1.1 * R, imag(g) + 1.1 * R,
	      10.0, 10.0, 100.0, 100.0);
  g_sel_scale(0);
  // 線種、マーカーの指定
  g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
  g_sel_line(0);
  g_def_marker(0, G_RED, 2, G_MARKER_CIRC);
  g_sel_marker(0);

  // 初期値と初期値を置く円を描く
  g_circle(real(g), imag(g), r0, G_YES, G_NO);
  for (i = 0; i < n; i++)
    g_marker(real(x[i]), imag(x[i]));

  // 反復
  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];
      // 図示する
      g_move(real(x[i]), imag(x[i]));
      g_plot(real(newx[i]), imag(newx[i]));
      g_marker(real(newx[i]), imag(newx[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]);
    cout << "変化量=" << error << endl;
    // 変化量が小さければ収束したと判断して反復を終了する。
    if (error < 1e-12)
      break;
  }

  // GLSC 終了の処理
  cout << "終了には" << endl;
  cout << " グラフィックスのウィンドウをクリック (X の場合)。" << endl;
  cout << " グラフィックスのウィンドウをクローズ (Windows の場合)。" << endl;
  g_sleep(-1.0);
  g_term();

  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;
}

//  Π (ｚi-ｚj) の計算
// j≠i
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;
}