POLYNOME CARACTERISTIQUE

//-------------------------------------------------
// DET
//-------------------------------------------------
// calcul d'un determinant avec le pivot de Gauss
// la complexite est en n^3
COMPLEX *Det(COMPLEX **_mat,int n,COMPLEX *res)
{
COMPLEX **mat;
COMPLEX   det;
int       x,y;
int       i,j;

mat = CreateMat(n);
for(y=0;y<n;y++)
  {
  for(x=0;x<n;x++)
    {
    mat[y][x].re = _mat[y][x].re;
    mat[y][x].im = _mat[y][x].im;
    }
  }

det.re = 1.;
det.im = 0.;

for(j=0;j<n-1;j++)
  {

  int       rankMax,rank;
  COMPLEX   coeffMax;

  // ( etape 1 )
  rankMax = j;
  for(rank=j+1;rank<n;rank++)
    {
    if(AbsComplex(&mat[rankMax][j]) < AbsComplex(&mat[rank][j]))
      {
      rankMax = rank;
      }
    }

  coeffMax.re = mat[rankMax][j].re;
  coeffMax.im = mat[rankMax][j].im;
  if(AbsComplex(&coeffMax) <= MATH_EPSILON)
    {
    det.re = 0.;
    det.im = 0.;
    goto label_end;
    }
  // ( etape 2 )
  if(rankMax != j)
    {
    for(i=j;i<n;i++)
      {
      COMPLEX tmp;
      tmp             = mat[j][i];
      mat[j][i]       = mat[rankMax][i];
      mat[rankMax][i] = tmp;
      }
    det.re *= -1.;
    det.im *= -1.;
    }

  MulComplex(&det,&det,&coeffMax);
  // ( etape 3 )
  for(rank=j+1;rank<n;rank++)
    {
    COMPLEX coeff;
    DivComplex(&coeff,&mat[rank][j],&coeffMax);
    for(i=j;i<n;i++)
      {
      COMPLEX mul;
      MulComplex(&mul,&coeff,&mat[j][i]);
      SubComplex(&mat[rank][i],&mat[rank][i],&mul);
      }
    }
    
  }

MulComplex(&det,&det,&mat[n-1][n-1]);

label_end:
DeleteMat(mat,n);
*res = det;
return res;
} // Det()


//-------------------------------------------------
// POLY_FFT
//-------------------------------------------------
// on decompose le polynome P = A0 + B0 avec les coefficients 
// de degres pairs (en A0) et impairs (en B0)
// alors on a :
// P(w) = A0(w)   +   B0(w)
//      = A(w^2)  + w.B(w^2)
// or w^2 est une racine nieme, donc on se ramene a un cas
// de taille deux fois moins grand
void polyFFT(
              P_COMPLEX     res,
              P_COMPLEX     arrayCoeff,
              int           stepi,
              int           n,
              P_COMPLEX     w
            )
{
  
if(1 == n)
  {
  // P(z)=constante
  res->re = arrayCoeff->re;
  res->im = arrayCoeff->im;
  }
else
  {
  int     k;
  COMPLEX w2,wPowk;


  // strategie diviser pour regner

  // calcul du carre de <w>
  MulComplex(&w2,w,w);
  // appel recursif
  polyFFT(res      ,arrayCoeff         ,2*stepi,n/2,&w2); // coefficients pairs
  polyFFT(res + n/2,arrayCoeff + stepi ,2*stepi,n/2,&w2); // coefficients impairs

  // on calcule la FFT a partir de ce que l'on vient de calculer
  wPowk.re = 1.;
  wPowk.im = 0.;
  for(k=0;k<n/2;k++)
    {
    COMPLEX A,B;

    A = res[k];
    B = res[k + n/2];
    MulComplex(&B,&B,&wPowk);

    AddComplex(res + k      ,&A,&B);
    SubComplex(res + k + n/2,&A,&B);

    MulComplex(&wPowk,&wPowk,w);
    }
  }
} // polyFFT()


//-------------------------------------------------
// PC
//-------------------------------------------------
// 
// calcul du polynome caracteristique
// Pc(M)(x) = det(M - x.I)
// 
P_POLYNOMIAL Pc(COMPLEX **_mat,int n)
{
int             len;        // aligne sur une puissance de 2
P_POLYNOMIAL    res;        // le polynome caracteristique
COMPLEX        *arrayDet;   // valeurs des determinants
COMPLEX       **mat;
int             k,x,y;
COMPLEX         conjw;

len       = alignedPower2(1 + n);
res       = CreatePolynomial(len - 1);
arrayDet  = Malloc(COMPLEX,len);
mat       = CreateMat(n);

// on copie la matrice
for(y=0;y<n;y++)
  {
  for(x=0;x<n;x++)
    {
    mat[y][x].re = _mat[y][x].re;
    mat[y][x].im = _mat[y][x].im;
    }
  }

// on evalue le polynome aux racines <len>-eme de l'unite
for(k=0;k<len;k++)
  {
  COMPLEX w;
  int     i;

  GetNthRootUnitComplex(&w,k,len);
  for(i=0;i<n;i++)
    {
    SubComplex(&mat[i][i],&_mat[i][i],&w);
    }
  // calcul P(w) dans <arrayDet + k>
  Det(mat,n,arrayDet + k);
  }

// on effectue une FFT pour retrouver les coefficients du polynome
// a partir de ses valeurs en certains points
GetNthRootUnitComplex(&conjw,-1,len);
polyFFT(res->coeff,arrayDet,1,len,&conjw);
CalcDeg(res);
for(k=0;k<=res->deg;k++)
  {
  res->coeff[k].re /= len;
  res->coeff[k].im /= len;
  }

Free(arrayDet);
DeleteMat(mat,n);
return res;
} // Pc()