- algorithme pour calculer le polynome caracteristique d'une matrice (carree) :
Pc(M)(x) = det(M - x.I)
- l'algorithme consiste a faire un FFT a partir de certaines valeurs du polynome
que l'on sait calculer avec le determinant d'un matrice (algo. en n^3)
- trois fonctions :
- det -> calcule le determinant d'un matrice (complexe)
- polyFFT -> calcule la FFT d'un tableau de nombres complexes
- Pc -> retourne le polynome caracteristique
- voir captures d'ecran pour des exemples de matrices
Source / Exemple :
//-------------------------------------------------
// 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);
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()
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.