/***************************************************
*   Find all roots of a complex polynomial using   *
*   Newton's iterative formulation.                *
* ------------------------------------------------ *
* SAMPLE RUN:                                      *
* (Find roots of complex polynomial:               *
*  (2.5 + I) X3 + (7.5 - 12 I) X2 + (-3.75 + 0.4 I *
*  ) X + (4.25 - I)  )                             *
*                                                  *
*  Error code = 0                                  *
*                                                  *
*  Complex roots are:                              *
*   0.25996740  -0.39966139                        *
*   -0.11184574  0.66800696                        *
*   -1.07915614  4.90406822                        *
*                                                  *
* ------------------------------------------------ *
* Ref.: From Numath Library By Tuan Dang Trong in  *
*       Fortran 77 [BIBLI 18].                     *
*                                                  *
*                C++ Release By J-P Moreau, Paris. *
*                      (www.jpmoreau.fr)           *
***************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  50


     int I,IERR, N;        //IERR: error code, N: order of polynomial
     double A[NMAX],B[NMAX],RR[NMAX],RI[NMAX];    
		          /* A: real parts of complex coefficients
                             B: imaginary parts of complex coefficients
                                (stored in decreasing powers)
                             RR: real parts of roots
                             RI: imaginary parts of roots */
     double W1[NMAX],W2[NMAX],W3[NMAX],W4[NMAX];    
	                       //working zones



    void CNEWTON (int N, double *A, double *B, double *X, double *Y, 
		 int *IER, double *G, double *H, double *T, double *F) {
/*--------------------------------------------------------------------
!     CALCULATE ALL THE ROOTS OF A COMPLEX POLYNOMIAL USING
!     NEWTON'S ITERATIVE FORMULATION
!
!     INPUTS:
!     N        ORDER OF POLYNOMIAL
!     A,B      TABLES OF DIMENSION N+1, RESPECTIVELY REAL,
!              IMAGINARY PART OF POLYNOMIAL COMPLEX COEFFICIENTS 
!              C = A+I*B, STORED IN DECREASING POWERS.
!     OUTPUTS:
!     X,Y      TABLES OF DIMENSION N, RESPECTIVELY REAL,
!              IMAGINARY PART OF COMPLEX ROOTS FOUND,  Z = X+I*Y
!     IER      ERROR CODE = 0, CONVERGENCE OK
!                         = I, NO CONVERGENCE FOR I-TH ROOT
!     WORKING ZONES:
!     G,H,T,F  TABLES OF DIMENSION N
!----------------------------------------------------------------------
!     REFERENCE:
!     E.DURAND -SOLUTIONS NUMERIQUES DES EQUATIONS ALGEBRIQUES,
!               TOME 1, MASSON & CIE, 1960, PP.260-266
!---------------------------------------------------------------------*/
//  Labels: e5,e7,e10,e15,e30,e40,e45,e55,e60,e70,e75
      int I,IT1,IT2,K,K1,K2,L,M;
      double EPS,EPS1,P,Q,R,TOL,X0,X1,Y0,Y1;
      double CABS,DX,DY,GD,GX,GY,R1,R2,R3,X2,Y2;

      IT1=25; IT2=100; TOL=1e-2; EPS=0.0;
      IER = 0;
      if (EPS != 0.0)  goto e7;
      EPS = 1.0;
e5:   EPS = 0.50*EPS;
      EPS1 = EPS+1.0;
      if (EPS1 > 1.0) goto e5;
e7:   L = N;
      X0 = 1.0;
      Y0 = 1.0;
      G[1] = A[1];
      T[1] = A[1];
      H[1] = B[1];
      F[1] = B[1];
e10:  if (L == 1) goto e55;
      M = L+1;
      I = N-L+1;
      K1 = 0;
      K2 = 0;
      X1 = X0;
      Y1 = Y0;
e15:  for (K=2; K<=M; K++) {
        G[K] = A[K]+X1*G[K-1]-Y1*H[K-1];
        H[K] = B[K]+Y1*G[K-1]+X1*H[K-1];
      }
      for (K=2; K<=L; K++) {
        T[K] = G[K]+X1*T[K-1]-Y1*F[K-1];
        F[K] = H[K]+Y1*T[K-1]+X1*F[K-1];
      }
      GD = T[L]*T[L]+F[L]*F[L];
      GX = G[M]*T[L]+H[M]*F[L];
      GY = H[M]*T[L]-G[M]*F[L];
      DX = -GX/GD;
      DY = -GY/GD;
      X2 = X1+DX;
      Y2 = Y1+DY;
      R1 = (fabs(DX)+fabs(DY))/(fabs(X2)+fabs(Y2));
      K2 = K2+1;
      X1 = X2;
      Y1 = Y2;
      if (R1 < TOL) goto e30;
      if (K2 <= IT2) goto e15;
      *IER = I;
      goto e60;
e30:  K1++;
      if (K1 <= IT1) goto e15;
      X[I] = X2;
      Y[I] = Y2;
      if (fabs(X2) <= EPS)  X2 = EPS;
      R2 = fabs(Y2/X2);
      if (R2 < 0.1) goto e40;
      if (R2*EPS > 1.0)  X[I] = 0.0;
      X0 = X2;
      Y0 = Y2;
      goto e45;
e40:  if (R2 <= EPS)  Y[I] = 0.0;
      X0 = X2+Y2;
      Y0 = X2-Y2;
e45:  for (K=1; K<=L; K++) {
        A[K] = G[K];
        B[K] = H[K];
      }
      L--;
      goto e10;
e55:  R3 = G[1]*G[1]+H[1]*H[1];
      X[N] = -(G[1]*G[2]+H[1]*H[2])/R3;
      Y[N] =  (H[1]*G[2]-G[1]*H[2])/R3;
      if (fabs(X[N]) <= EPS)   X[N] = EPS;
      R2 = fabs(Y[N]/X[N]);
      if (R2*EPS > 1.0)   X[N] = 0.0;
      if (R2 <= EPS)   Y[N] = 0.0;
e60:  for (L=1; L<=N; L++) {
        P=X[L];
        Q=Y[L];
        R=P*P+Q*Q;
        if (L == 1) goto e70;
        for (I=L; I>1; I--) {
          CABS=X[I-1]*X[I-1]+Y[I-1]*Y[I-1];
          if (R >= CABS) goto e75;
          X[I]=X[I-1];
          Y[I]=Y[I-1];
        }
e70:    I=1;
e75:    X[I]=P;
        Y[I]=Q;
      }
    } //CNEWTON

void main()  {

  N = 3;   //order of polynomial
 
  //define vectors A, B of coefficients
  A[1]= 2.5;  B[1]=1.0;
  A[2]= 7.5;  B[2]=-12.0;
  A[3]=-3.75; B[3]=0.4;
  A[4]= 4.25; B[4]=-1.0;

/* Example #2 (real roots are -3, 1, 2)
  A[1]= 1.;   B[1]=0.0;
  A[2]= 0.;   B[2]=0.0;
  A[3]=-7.;   B[3]=0.0;
  A[4]= 6.;   B[4]=0.0;  */

  //call main function
  CNEWTON(N,A,B,RR,RI,&IERR,W1,W2,W3,W4);

  //print results
  printf("\n Error code = %d\n", IERR);
  printf("\n Complex roots are:\n");

  for (I=1; I<=N; I++)
    printf("  %10.8f  %10.8f\n", RR[I], RI[I]);
  printf("\n");

}

// end of file tnewton.cpp