/*
  Author:  Pate Williams (c) 1997

  Complex roots of a squarefree polynomial by
  the modified Newton's method. Uses the
  standard C++ library class complex.
*/

#include <iostream.h>
#include <math.h>
#include <complex>
using namespace std;

const int SIZE = 8192;

typedef complex<double> dcomplex;

dcomplex Evaluate(int n, dcomplex x, dcomplex *P)
//evaluation of a polynomial using Horner's method
{
  dcomplex value = P[n];

  for (int i = n - 1; i >= 0; i--)
    value = value * x + P[i];
  return value;
}

void Roots(bool real, double epsilon, int degreeP, dcomplex *P, dcomplex *z)
{
  int c, count = 0, f = real ? 1 : 0, i, n = degreeP;
  double m, m1;
  dcomplex Px, PPx, alpha, beta, dx, v, v1, x, x0(1.3, 0.314159), y;
  dcomplex A[SIZE], B[SIZE], PP[SIZE], Q[SIZE], QP[SIZE];

  epsilon = fabs(epsilon);
  for (i = 0; i <= n; i++) Q[i] = P[i];
  //calculate the derivative of the polynomial
  for (i = 1; i <= n; i++) {
    PP[i - 1] = dcomplex(i, 0.0) * P[i];
    QP[i - 1] = PP[i - 1];
  }
  L2:
    //initialize the root with a "small" imaginary part
    x = x0;
    v = Evaluate(n, x, Q);
    m = pow(abs(v), 2.0);
  L3:
    c = 0;
    dx = v / Evaluate(n - 1, x, QP);
    if (abs(dx) < epsilon) goto L5;
  L4:
    y = x - dx;
    v1 = Evaluate(n, y, Q);
    m1 = pow(abs(v1), 2.0);
    if (m1 < m) {
      x = y, v = v1, m = m1;
      goto L3;
    }
    c++;
    dx /= dcomplex(4.0, 0.0);
    if (c < 20) goto L4;
  L5:
    //evaluate the polynomial and its derivative
    Px = Evaluate(degreeP, x, P);
    PPx = Evaluate(degreeP - 1, x, PP);
    //apply Newton's formula
    if (PPx != dcomplex(0.0, 0.0)) {
      x -= Px / PPx;
      Px = Evaluate(degreeP, x, P);
      PPx = Evaluate(degreeP - 1, x, PP);
      //apply Newton's formula a second time
      if (PPx != dcomplex(0.0, 0.0)) x -= Px / PPx;
    }
    if (f == 0 || f == 1 && fabs(x.imag()) <= epsilon) {
      x = dcomplex(x.real(), 0.0);
      z[count++] = x;
      //deflate the polynomial by a linear factor
      A[n - 1] = Q[n];
      for (i = n - 1; i >= 1; i--)
        A[i - 1] = Q[i] + x * A[i];
      n--;
      for (i = 0; i <= n; i++) Q[i] = A[i];
    }
    else {
      z[count] = x;
      z[count + 1] = conj(x);
      count += 2;
      //deflate the polynomial by a quadratic factor
      alpha = dcomplex(2.0 * x.real(), 0.0);
      beta = dcomplex(pow(abs(x), 2.0), 0.0);
      B[n - 2] = Q[n];
      B[n - 3] = Q[n - 1] + alpha * B[n - 2];
      for (i = n - 2; i >= 2; i--)
        B[i - 2] = Q[i] + alpha * B[i - 1] - beta * B[i];
      n -= 2;
      for (i = 0; i <= n; i++) Q[i] = B[i];
    }
    //calculate the derivative of Q(x)
    for (i = 1; i <= n; i++)
      QP[i - 1] = dcomplex(i, 0.0) * Q[i];
    if (n > 0) goto L2;
}

void main(void)
{
  bool real = true;
  double im, re;
  int i, degreeP;
  dcomplex P[SIZE], roots[SIZE];

  cout << "the degree of the polynomial = ";
  cin >> degreeP;
  for (i = 0; i <= degreeP; i++) {
    cout << "P[" << i << "].real = ";
    cin >> re;
    cout << "P[" << i << "].imag = ";
    cin >> im;
    P[i] = dcomplex(re, im);
    if (real && im != 0.0) real = false;
  }
  //real = true means the coefficients of P(x) are all real
  Roots(real, 1.0e-08, degreeP, P, roots);
  cout << "z\tP(z)" << endl;
  for (i = 0; i < degreeP; i++)
    cout << roots[i] << ' ' << Evaluate(degreeP, roots[i], P) << endl;
}