/********************************************************
* Calculate the coefficients of the characteristic poly-*
* nomial P(l)=A-l*I of a square complex matrix A(i,j)   *
* ----------------------------------------------------- *
* Ref.: "Algèbre, Algorithmes et programmes en Pascal   *
*        By Jean-Louis Jardrin, DUNOD Paris, 1988"      *
*        [BIBLI 10].                                    * 
* ----------------------------------------------------- *
* SAMPLE RUN:                                           *
* (Find the characteristic polynomial of matrix:        *
*           1 -1  3                                     *
*      A = -1  i -1-2i                                  *
*           i  1 -2+i   )                               *
*                                                       *
* Input matrix size (max 20): 3                         *
*                                                       *
* Line 1                                                *
*  Column 1                                             *
*   Real  part: 1                                       *
*   Imag. part: 0                                       *
*  Column 2                                             *
*   Real  part: -1                                      *
*   Imag. part: 0                                       *
*  Column 3                                             *
*   Real  part: 3                                       *
*   Imag. part: 0                                       *
* Line 2                                                *
*  Column 1                                             *
*   Real  part: -1                                      *
*   Imag. part: 0                                       *
*  Column 2                                             *
*   Real  part: 0                                       *
*   Imag. part: 1                                       *
*  Column 3                                             *
*   Real  part: -1                                      *
*   Imag. part: -2                                      *
* Line 3                                                *
*  Column 1                                             *
*   Real  part: 0                                       *
*   Imag. part: 1                                       *
*  Column 2                                             *
*   Real  part: 1                                       *
*   Imag. part: 0                                       *
*  Column 3                                             *
*   Real  part: -2                                      *
*   Imag. part: 1                                       *
*                                                       *
* Characteristic polynomial P(l):                       *
* Coefficient for degree 3:                             *
* -1.000000 +0.000000 i                                 *
* Coefficient for degree 2:                             *
* -1.000000 +2.000000 i                                 *
* Coefficient for degree 1:                             *
* 3.000000 +1.000000 i                                  *
* Coefficient for degree 0:                             *
* 0.000000 +0.000000 i                                  *
*                                                       *
* (The characteristic polynomial of matrix A is:        *
*   P(l)=-l^3+(-1+2i)l^2+(3+i)l                         *
*                                                       *
*                 C++ version by J-P Moreau, Paris      *
*                (with statically allocated vectors).   *
*                       (www.jpmoreau.fr)               *
*                                                       *
* To be linked with ucomplex.cpp.                       *
********************************************************/
#include "ucomplex1.h" //elementary operations on complex numbers

#define  NMAX  20
#define  NMAXP 21

typedef  Complex MATC[NMAX][NMAX];
typedef  Complex VECC[NMAXP];

int k,n;

MATC  A;
VECC  P;


void Read_data(int n) {
  int i,j;
  for (i=1; i<=n; i++) {
    printf("\n Line %d\n",i);
    for (j=1; j<=n; j++) {
      printf("  Column %d\n",j);
      printf("   real  part: "); scanf("%lf",&A[i][j].r);
      printf("   imag. part: "); scanf("%lf",&A[i][j].i);
    }
  }
}

//calculate complex trace of complex C matrix
void TRC(int n, MATC C, Complex *t0) {
  int i; Complex c0;
  t0->r=0.0; t0->i=0.0;
  for (i=1; i<=n; i++) {
    c0.r=t0->r; c0.i=t0->i;
    CADD(c0,C[i][i],t0);
  }
}

/*****************************************************
* The procedure PCCS calculates the complex coeffi-  *
* cients of the characteristic polynomial P(l)=A-l*I *
* of a complex square matrix A(i,j).                 *
*                                                    *
* Note: the roots of P(l) are the complex eigenvalues*
*       of matrix A(i,j).                            *
* -------------------------------------------------- *
* INPUTS:                                            *
*          n: zize of matrix (integer)               *
*          A: complex matrix of size (n,n)           *
*             (see main program).                    *
* OUTPUT:                                            *
*          P: vector of complex coefficients of P(l) *
*****************************************************/
void PCCS(int n, MATC A, VECC P) {
  int i,j,k,l;
  Complex c0,c1,s,t0,z0,z1;
  MATC B, C; //temporary complex matrices

  z0.r=0.0; z0.i=0.0;
  z1.r=1.0; z1.i=0.0;
  if (n%2!=0) CPRO(-1.0,z1,&P[1]);
              else P[1]=z1;
  for (l=1; l<=n; l++) {
    if (l==1) {
      for (i=1; i<=n; i++)
        for (j=1; j<=n; j++)
          C[i][j]=A[i][j];
    }
    else {
      for (i=1; i<=n; i++) {
        for (j=1; j<=n; j++) {
          s.r=z0.r; s.i=z0.i;
          for (k=1; k<=n; k++) {
            c0.r=s.r;  c0.i=s.i;
            CMUL(B[i][k],A[k][j],&c1);
            CADD(c0,c1,&s);
          }
          C[i][j].r=s.r;
          C[i][j].i=s.i;
        }
      }
    }
    TRC(n,C,&c0);
    CPRO(1.0/l,c0,&t0);
    CMUL(t0,P[1],&c0);
    CPRO(-1.0,c0,&P[l+1]);
    if (l<n)
      for (i=1;i<=n;i++)
        for (j=1;j<=n;j++)
          if (j==i) CDIF(C[i][j],t0,&B[i][j]);
                    else B[i][j]=C[i][j];
  }
} //PCCS()

void main()  {
  printf("\n Input matrix size (max %d): ",NMAX); scanf("%d",&n);
  Read_data(n);
  //call routine to calculate coefficients
  PCCS(n,A,P);
  //print results
  printf("\n Characteristic polynomial P(l):\n");
  for (k=1; k<=n+1; k++) {
    printf(" Coefficient for degree %d: \n", n-k+1);
    printf("%f",P[k].r);
    if (P[k].i>=0) printf(" +");
                   else printf(" -");
    printf("%f i\n",fabs(P[k].i));               
  }
  printf("\n\n");
}

//end of file carpol1.cpp