1 #ifndef utilitaire_algebre_h
2 #define utilitaire_algebre_h
4 /* rpoly.cpp -- Jenkins-Traub real polynomial root finder.
6 * (C) 2000, C. Bond. All rights reserved.
8 * Translation of TOMS493 from FORTRAN to C. This
9 * implementation of Jenkins-Traub partially adapts
10 * the original code to a C environment by restruction
11 * many of the 'goto' controls to better fit a block
12 * structured form. It also eliminates the global memory
13 * allocation in favor of local, dynamic memory management.
15 * The calling conventions are slightly modified to return
16 * the number of roots found as the function value.
19 * op - double precision vector of coefficients in order of
21 * degree - integer degree of polynomial
24 * zeror,zeroi - output double precision vectors of the
25 * real and imaginary parts of the zeros.
28 * returnval: -1 if leading coefficient is zero, otherwise
29 * number of roots found.
32 /* A roots_polynoms class */
34 /*! \class roots_polynoms utilitaire_algebre.hxx "utilitaire_algebre.hxx"
35 * \brief Computes the roots of a given polynomial
36 * \details Translation of TOMS493 from FORTRAN to C. This
37 * implementation of Jenkins-Traub partially adapts
38 * the original code to a C environment by restruction
39 * many of the 'goto' controls to better fit a block
40 * structured form. It also eliminates the global memory
41 * allocation in favor of local, dynamic memory management.
42 * The calling conventions are slightly modified to return
43 * the number of roots found as the function value.
59 * \brief calcule les racines d'un polynome de degré 2
60 * \Details Calculate the zeros of the quadratic a*x² + b*x + c.
61 * The quadratic formula, modified to avoid overflow, is used
62 * to find the larger zero if the zeros are real and both
63 * are complex. The smaller real zero is found directly from
64 * the product of the zeros c/a
65 * @param a double correspond au coefficient de x²
66 * @param b double correspond au coefficient de x
67 * @param c double correspond au coefficient de 1
68 * @param *sr double correspond à la pratie réelle de la premiere racine
69 * @param *si double correspond à la pratie imaginaire de la premiere racine
70 * @param *lr double correspond à la pratie réelle de la deuxième racine
71 * @param *li double correspond à la pratie imaginaire de la deuxième racine
72 * @return calcule les racine du polynome dans sr+i*si et lr+i*li
74 void quad(double a,double b,double c,double *sr,double *si, double *lr,double *li);
78 * \Details Computes up to L2 fixed shift k-polynomials,
79 * testing for convergence in the linear or quadratic
80 * case. Initiates one of the variable shift
81 * iterations and returns with the number of zeros found
85 void fxshfr(int l2, int *nz);
89 * \Details /* Variable-shift k-polynomial iteration for a
90 * quadratic factor converges only if the zeros are
91 * equimodular or nearly so.
92 * @param *uu, *vv double coefficients of starting quadratic
93 * @param *nz entier number of zeros found
95 void quadit(double *uu,double *vv,int *nz);
98 * \brief Variable-shift H polynomial iteration for a real zero
99 * @param *sss starting iterate (double)
100 * @param *nz number of zeros found (int)
101 * @param *iflag flag to indicate a pair of zeros near real axis (int)
103 void realit(double *sss, int *nz, int *iflag);
106 * \Details This routine calculates scalar quantities used to
107 * compute the next k polynomial and new estimates of
108 * the quadratic coefficients.
109 * @param *type integer variable set here
110 * indicating how the calculations are normalized to avoid overflow
112 void calcsc(int *type);
115 * \brief Computes the next k polynomials using scalars computed in calcsc().
116 * @param *type integer variable set here
117 * indicating how the calculations are normalized to avoid overflow
119 void nextk(int *type);
122 * \brief Compute new estimates of the quadratic coefficients
123 * using the scalars computed in calcsc ().
125 * @param *uu, *vv double coefficients of starting quadratic
127 void newest(int type,double *uu,double *vv);
130 * \brief Divides p by the quadratic 1,u,v placing the quotient
131 * in q and the remainder in a,b.
133 void quadsd(int n,double *u,double *v,double *p,double *q,
134 double *a,double *b);
143 int rpoly(double *op, int degree, double *zeror, double *zeroi);
150 double *p,*qp,*k,*qk,*svk;
151 double sr,si,u,v,a,b,c,d,a1,a2;
152 double a3,a6,a7,e,f,g,h,szr,szi,lzr,lzi;
161 static void abs_par_interp_directe(int n, vector< complex< double > > vp, double * A, int sizeA, double epsilon, double *result,vector< complex< double > > y);
163 /** \fn polynome_caracteristique
164 * \brief Calcule le polynome caracteristique
166 static vector< double > polynome_caracteristique(double alpha_v, double alpha_l, double Uv, double Ul, double rhov, double rhol, double cv_2, double cl_2, double dPiv, double dPil);//c= inverse vitesse du son!
167 /** \fn polynome_caracteristique
184 static vector< double > polynome_caracteristique(double alpha1, double alpha2, double u1, double u2, double rho1, double rho2, double invc1sq, double invc2sq, double dpi1, double dpi2, double g2press, double g2alpha, double g2, double epsilon);
187 * \brief calcule le module d'un nombre complexe
188 * @param z est un nombre complexe
189 * @return calcule le module de z*/
190 static double module(complex< double > z);
191 /** \fn modulec calcule le module² de z */
192 static double modulec(complex< double > z);
193 /** \fn abs_generalise
194 * \brief calcule la valeur absolue d'un nombre complexe
195 * \Details calcule la valeur absolue d'un nombre complexe en prenant en compte
196 * que la partie réelle c-a-d si z= a+ib abs_generalize() renvoie |a|+ib
198 * @return si z = a+ib elle renvoie |a|+ib */
199 static complex< double > abs_generalise(complex< double > z);
201 static int new_tri_selectif(vector< complex< double > > &L, int n, double epsilon);
209 template< typename T >
210 static int tri_selectif(T& L, int n, double epsilon);
212 /** \fn matrixProduct
213 * \brief calcule le produit de 2 matrices .
214 * @param *Matrix1 la 1ere matrice
215 * @param R1,C1 la taille de la matrice1 Nbr de ligne et le nbr de colonne
216 * @param Matrix2 la 2ieme matrice
217 * @param R2,C2 la taille de la matrice2 Nbr de ligne et le nbr de colonne
218 * @param out le resultat de Matrix1*Matrix2
220 static void matrixProduct
222 const double *Matrix1,
225 const double *Matrix2,
231 /** \fn matrixProdVec
232 * \brief Calcule le produit Matrice vecteur
233 * @param *Matrix correspond à la matrice du produit
234 * @param R1,C1 la taille de la matrice( Nbr de ligne et le nbr de colonne)
235 * @param Vector correspond au vecteur du produit
236 * @param out le resultat de Matrix*Vector
238 static void matrixProdVec
239 ( const double *Matrix,
242 const double *Vector,
246 static bool belongTo(complex< double > a , vector< complex <double > > v, double epsilon);
251 * \brief Calcule la somme de 2vecteurs ( U=U+V)
252 * @param *U,*V les 2 vecteurs
253 * @param N taille des 2 vecteurs
254 * @return la somme de U et V dans U*/
257 double *U, //vecteur auquel on ajoute
258 int N, //taille du vecteur
259 const double *V //ajout
262 * \brief Calcule le tenseur de 2 vecteurs
263 * @param *a le 1er vecteur (*double)
264 * @param N la taille du premier vecteur (int)
265 * @param *b le 2ieme vecteur (*double)
266 * @param M la taille du 2ieme vecteur (int)
267 * @param *ab le resultat
271 ( const double *a, //vecteur gauche
272 int N, //taille du vecteur gauche
273 const double *b, //vecteur droit
274 int M, //taille du vecteur droit
275 double *ab //conteneur de sortie
278 /** \fn shift_diagonal
279 * \brief Calcule la trace d'une matrice carrée
280 * @param *Matrix correspond à la matrice
281 * @param N taille de la matrice
282 * @param shift résultat
283 * @return renvoie la trace de la matrice Matrix dans shift */
284 static void shift_diagonal( double *Matrix, const int N, double shift);
285 //qques outils d'algebre des polynomes
287 /** \fn remplir_trinome
288 * \brief remplie un trinome
290 static void remplir_trinome(double coeff, double u, double alpha, vector< double >& trinome);
293 * \brief Additionne 2 vecteurs de tailles differentes
294 * @param P,Q les 2 vecteurs à additionner
295 * @param n,m les 2 tailles respectives des vecteurs P et Q
296 * @return un vecteur qui correspond a P+Q*/
297 static vector< double > additionne(const vector< double >& P, const vector< double >& Q, int n, int m);
300 * \brief Calcule le produit scalaire de 2 vecteurs de tailles differentes
301 * @param P,Q les 2 vecteurs
302 * @param n,m les 2 tailles respectives des 2 vecteurs
303 * @return un vecteur correspond à (P,Q) */
304 static vector< double > multiplie(const vector< double >& P, const vector< double >& Q, int n, int m);
306 //Pour le tri des valeurs propres
308 static void ordre_croissant_abs(vector< complex< double > > &L, int n);
310 //Calcul des coefficients du polynome d'interpolation x->y dans la base de Newton
320 static vector< complex< double > > dif_div(int n, const vector< complex< double > >& x, const vector< T >& y, double epsilon);//n=nb valeurs à interpoler
322 //attention, Les vp complexes conjugu�es doivent se suivre dans la liste x
323 static void appliquer_dif_div(int n, const vector< complex< double > >& dif_div, const vector< complex< double > >& x, const double* A, const int sizeA, double epsilon, double *p);//p=result
325 static double avr(double a, double b);