1 #include "utilitaire_algebre.h"
5 roots_polynoms::roots_polynoms(): smalno(1.2e-38), infin(3.4e38),
10 int roots_polynoms::rpoly(double *op, int degree, double *zeror, double *zeroi)
12 double t,aa,bb,cc,*temp,factor,rot;
14 double lo,max,min,xx,yy,cosr,sinr,xxx,x,sc,bnd;
16 int cnt,nz,i,j,jj,l,nm1,zerok;
18 /* The following statements set machine constants. */
22 /* Initialization of constants for shift rotation. */
30 /* Algorithm fails of the leading coefficient is zero. */
31 if (fabs(op[0]) <= smalno) return -1;
32 /* Remove the zeros at the origin, if any. */
33 while (op[n] == 0.0) {
42 * Allocate memory here
44 temp = new double [degree+1];
45 pt = new double [degree+1];
46 p = new double [degree+1];
47 qp = new double [degree+1];
48 k = new double [degree+1];
49 qk = new double [degree+1];
50 svk = new double [degree+1];
51 /* Make a copy of the coefficients. */
54 /* Start the algorithm for one zero. */
57 zeror[degree-1] = -p[1]/p[0];
58 zeroi[degree-1] = 0.0;
62 /* Calculate the final zero or pair of zeros. */
64 quad(p[0],p[1],p[2],&zeror[degree-2],&zeroi[degree-2],
65 &zeror[degree-1],&zeroi[degree-1]);
69 /* Find largest and smallest moduli of coefficients. */
75 if (x != 0.0 && x < min){ min = x;
76 //cerr << " x= " << x << " smalno " << smalno << endl;
79 /* Scale if there are large or very small coefficients.
80 * Computes a scale factor to multiply the coefficients of the
81 * polynomial. The scaling si done to avoid overflow and to
82 * avoid undetected underflow interfering with the convergence
83 * criterion. The factor is a power of the base.
86 if (sc > 1.0 && infin/sc < max) goto _110;
88 if (max < 10.0) goto _110;
92 l = (int)(log(sc)/log(base) + 0.5);
93 factor = pow(base*1.0,l);
96 p[i] = factor*p[i]; /* Scale polynomial. */
99 /* Compute lower bound on moduli of roots. */
101 pt[i] = (fabs(p[i]));
104 /* Compute upper estimate of bound. */
105 x = exp((log(-pt[n])-log(pt[0])) / (double)n);
106 /* If Newton step at the origin is better, use it. */
107 //if (pt[n-1] != 0.0) {
108 if (fabs(pt[n-1] ) > smalno*fabs(pt[n])) {
112 /* Chop the interval (0,x) until ff <= 0 */
118 if (ff <= 0.0) break;
122 /* Do Newton interation until x converges to two
125 while (fabs(dx/x) > 0.005) {
137 /* Compute the derivative as the initial k polynomial
138 * and do 5 steps with no shift.
142 k[i] = (double)(n-i)*p[i]/(double)n;
146 zerok = (fabs(k[n-1])< smalno);
147 for(jj=0;jj<5;jj++) {
149 if (!zerok) {//cerr << "cc= " << cc << endl;
150 /* Use a scaled form of recurrence if value of k at 0 is nonzero. */
152 for (i=0;i<nm1;i++) {
154 k[j] = t*k[j-1]+p[j];
157 zerok = (fabs(k[n-1]) <= fabs(bb)*eta*10.0 || (fabs(k[n-1])< smalno));
160 /* Use unscaled form of recurrence. */
161 for (i=0;i<nm1;i++) {
166 zerok = (fabs(k[n-1])< smalno);
169 /* Save k for restarts with new shifts. */
172 /* Loop to select the quadratic corresponding to each new shift. */
173 for (cnt = 0;cnt < 20;cnt++) {
174 /* Quadratic corresponds to a double shift to a
175 * non-real point and its complex conjugate. The point
176 * has modulus bnd and amplitude rotated by 94 degrees
177 * from the previous shift.
179 xxx = cosr*xx - sinr*yy;
180 yy = sinr*xx + cosr*yy;
186 fxshfr(20*(cnt+1),&nz);
188 /* The second stage jumps directly to one of the third
189 * stage iterations and returns here if successful.
190 * Deflate the polynomial, store the zero or zeros and
191 * return to the main algorithm.
205 /* If the iteration is unsuccessful another quadratic
206 * is chosen after restoring k.
212 /* Return with failure if no convergence with 20 shifts. */
222 // cerr << " rpoly roots " << endl;
223 // for (int ii=0; ii<degree;ii++) {
224 // cerr << " (" << zeror[ii] << "," << zeroi[ii] << ")";
227 // cerr << " rpoly roots point " << &zeror[0] << " " << &zeroi[0] << endl;
228 // cerr << " rpoly roots point 1 " << &zeror[1] << " " << &zeroi[1] << endl;
234 /* Computes up to L2 fixed shift k-polynomials,
235 * testing for convergence in the linear or quadratic
236 * case. Initiates one of the variable shift
237 * iterations and returns with the number of zeros
240 void roots_polynoms::fxshfr(int l2,int *nz)
242 double svu,svv,ui,vi,s;
243 double betas,betav,oss,ovv,ss,vv,ts,tv;
244 double ots,otv,tvv,tss;
245 int type, i,j,iflag,vpass,spass,vtry,stry;
252 /* Evaluate polynomial by synthetic division. */
253 quadsd(n,&u,&v,p,qp,&a,&b);
256 /* Calculate next k polynomial and estimate v. */
259 newest(type,&ui,&vi);
263 if (k[n-1] != 0.0) ss = -p[n]/k[n-1];
266 if (j == 0 || type == 3) goto _70;
267 /* Compute relative measures of convergence of s and v sequences. */
268 if (vv != 0.0) tv = fabs((vv-ovv)/vv);
269 if (ss != 0.0) ts = fabs((ss-oss)/ss);
270 /* If decreasing, multiply two most recent convergence measures. */
272 if (tv < otv) tvv = tv*otv;
274 if (ts < ots) tss = ts*ots;
275 /* Compare with convergence criteria. */
276 vpass = (tvv < betav);
277 spass = (tss < betas);
278 if (!(spass || vpass)) goto _70;
279 /* At least one sequence has passed the convergence test.
280 * Store variables before iterating.
288 /* Choose iteration according to the fastest converging
293 if (spass && (!vpass) || tss < tvv) goto _40;
297 /* Quadratic iteration has failed. Flag that it has
298 * been tried and decrease the convergence criterion.
302 /* Try linear iteration if it has not been tried and
303 * the S sequence is converging.
305 if (stry || !spass) goto _50;
310 realit(&s,nz,&iflag);
312 /* Linear iteration has failed. Flag that it has been
313 * tried and decrease the convergence criterion.
317 if (iflag == 0) goto _50;
318 /* If linear iteration signals an almost double real
319 * zero attempt quadratic iteration.
324 /* Restore variables. */
331 /* Try quadratic iteration if it has not been tried
332 * and the V sequence is convergin.
334 if (vpass && !vtry) goto _20;
335 /* Recompute QP and scalar values to continue the
338 quadsd(n,&u,&v,p,qp,&a,&b);
347 /* Variable-shift k-polynomial iteration for a
348 * quadratic factor converges only if the zeros are
349 * equimodular or nearly so.
350 * uu, vv - coefficients of starting quadratic.
351 * nz - number of zeros found.
353 void roots_polynoms::quadit(double *uu,double *vv,int *nz)
356 double mp,omp,ee,relstp,t,zm;
366 quad(1.0,u,v,&szr,&szi,&lzr,&lzi);
367 /* Return if roots of the quadratic are real and not
368 * close to multiple or nearly equal and of opposite
371 if (fabs(fabs(szr)-fabs(lzr)) > 0.01 * fabs(lzr)) return;
372 /* Evaluate polynomial by quadratic synthetic division. */
373 quadsd(n,&u,&v,p,qp,&a,&b);
374 mp = fabs(a-szr*b) + fabs(szi*b);
375 /* Compute a rigorous bound on the rounding error in
379 ee = 2.0*fabs(qp[0]);
382 ee = ee*zm + fabs(qp[i]);
384 ee = ee*zm + fabs(a+t);
385 ee *= (5.0 *mre + 4.0*are);
386 ee = ee - (5.0*mre+2.0*are)*(fabs(a+t)+fabs(b)*zm)+2.0*are*fabs(t);
387 /* Iteration has converged sufficiently if the
388 * polynomial value is less than 20 times this bound.
395 /* Stop iteration after 20 steps. */
398 if (relstp > 0.01 || mp < omp || tried) goto _50;
399 /* A cluster appears to be stalling the convergence.
400 * Five fixed shift steps are taken with a u,v close
403 if (relstp < eta) relstp = eta;
404 relstp = sqrt(relstp);
407 quadsd(n,&u,&v,p,qp,&a,&b);
416 /* Calculate next k polynomial and new u and v. */
420 newest(type,&ui,&vi);
421 /* If vi is zero the iteration is not converging. */
422 if (vi == 0.0) return;
423 relstp = fabs((vi-v)/vi);
428 /* Variable-shift H polynomial iteration for a real zero.
429 * sss - starting iterate
430 * nz - number of zeros found
431 * iflag - flag to indicate a pair of zeros near real axis.
433 void roots_polynoms::realit(double *sss, int *nz, int *iflag)
446 /* Evaluate p at s. */
453 /* Compute a rigorous bound on the error in evaluating p. */
455 ee = (mre/(are+mre))*fabs(qp[0]);
457 ee = ee*ms + fabs(qp[i]);
459 /* Iteration has converged sufficiently if the polynomial
460 * value is less than 20 times this bound.
462 if (mp <= 20.0*((are+mre)*ee-mre*mp)) {
469 /* Stop iteration after 10 steps. */
472 if (fabs(t) > 0.001*fabs(s-t) || mp < omp) goto _50;
473 /* A cluster of zeros near the real axis has been
474 * encountered. Return with iflag set to initiate a
475 * quadratic iteration.
480 /* Return if the polynomial value has increased significantly. */
483 /* Compute t, the next polynomial, and the new iterate. */
490 if (fabs(kv) <= fabs(k[n-1])*10.0*eta) {
491 /* Use unscaled form. */
498 /* Use the scaled form of the recurrence if the value
499 * of k at s is nonzero.
504 k[i] = t*qk[i-1] + qp[i];
512 if (fabs(kv) > fabs(k[n-1]*10.0*eta)) t = -pv/kv;
517 /* This routine calculates scalar quantities used to
518 * compute the next k polynomial and new estimates of
519 * the quadratic coefficients.
520 * type - integer variable set here indicating how the
521 * calculations are normalized to avoid overflow.
523 void roots_polynoms::calcsc(int *type)
525 /* Synthetic division of k by the quadratic 1,u,v */
526 quadsd(n-1,&u,&v,k,qk,&c,&d);
527 if (fabs(c) > fabs(k[n-1]*100.0*eta)) goto _10;
528 if (fabs(d) > fabs(k[n-2]*100.0*eta)) goto _10;
530 /* Type=3 indicates the quadratic is almost a factor of k. */
533 if (fabs(d) < fabs(c)) {
535 /* Type=1 indicates that all formulas are divided by c. */
540 a3 = a*e + (h/c+g)*b;
546 /* Type=2 indicates that all formulas are divided by d. */
551 a3 = (a+g)*e + h*(b/d);
555 /* Computes the next k polynomials using scalars
556 * computed in calcsc.
558 void roots_polynoms::nextk(int *type)
564 /* Use unscaled form of the recurrence if type is 3. */
573 if (*type == 1) temp = b;
574 if (fabs(a1) <= fabs(temp)*eta*10.0) {
575 /* If a1 is nearly zero then use a special form of the
581 k[i] = a3*qk[i-2] - a7*qp[i-1];
585 /* Use scaled form of the recurrence. */
589 k[1] = qp[1] - a7*qp[0];
591 k[i] = a3*qk[i-2] - a7*qp[i-1] + qp[i];
594 /* Compute new estimates of the quadratic coefficients
595 * using the scalars computed in calcsc.
597 void roots_polynoms::newest(int type,double *uu,double *vv)
599 double a4,a5,b1,b2,c1,c2,c3,c4,temp;
600 // double smalno = 1.2e-38;
601 // double smalno = 0.0;
603 /* Use formulas appropriate to setting of type. */
605 /* If type=3 the quadratic is zeroed. */
618 /* Evaluate new quadratic coefficients. */
620 b2 = -(k[n-2]+b1*p[n-1])/p[n];
625 temp = a5 + b1*a4 - c4;
626 if (fabs(temp) <= smalno) {
631 *uu = u - (u*(c3+c2)+v*(b1*a1+b2*a7))/temp;
632 *vv = v*(1.0+c4/temp);
636 /* Divides p by the quadratic 1,u,v placing the quotient
637 * in q and the remainder in a,b.
639 void roots_polynoms::quadsd(int my_n,double *my_u,double *my_v,
640 double *my_p,double *my_q,
641 double *my_a,double *my_b)
647 *my_a = my_p[1] - (*my_b)*(*my_u);
649 for (i=2;i<=my_n;i++) {
650 my_c = my_p[i] - (*my_a)*(*my_u) - (*my_b)*(*my_v);
656 /* Calculate the zeros of the quadratic a*z^2 + b1*z + c.
657 * The quadratic formula, modified to avoid overflow, is used
658 * to find the larger zero if the zeros are real and both
659 * are complex. The smaller real zero is found directly from
660 * the product of the zeros c/a.
662 void roots_polynoms::quad(double my_a,double my_b1,double my_c,double *my_sr,double *my_si,
663 double *my_lr,double *my_li)
665 double my_b,my_d,my_e;
666 // double smalno = 1.2e-38;
668 if (fabs(my_a) <= smalno) { /* less than two roots */
669 if (fabs(my_b1) >= smalno)
670 *my_sr = -my_c/my_b1;
678 if (fabs(my_c) <= smalno) { /* one real root, one zero root */
680 *my_lr = -my_b1/my_a;
685 /* Compute discriminant avoiding overflow. */
687 if (fabs(my_b) < fabs(my_c)) {
692 my_e = my_b*(my_b/fabs(my_c)) - my_e;
693 my_d = sqrt(fabs(my_e))*sqrt(fabs(my_c));
696 my_e = 1.0 - (my_a/my_b)*(my_c/my_b);
697 my_d = sqrt(fabs(my_e))*fabs(my_b);
699 if (my_e < 0.0) { /* complex conjugate zeros */
702 *my_si = fabs(my_d/my_a);
706 if (my_b >= 0.0) /* real zeros. */
708 *my_lr = (-my_b+my_d)/my_a;
711 *my_sr = (my_c/ *my_lr)/my_a;
717 //qques outils d'algebre lineaire
720 * \brief Calcule la somme de 2vecteurs ( U=U+V)
721 * @param *U,*V les 2 vecteurs
722 * @param N taille des 2 vecteurs
723 * @return la somme de U et V dans U*/
726 double *U, //vecteur auquel on ajoute
727 int N, //taille du vecteur
728 const double *V //ajout
731 for(int i=0;i<N;i++) U[i] += V[i];
735 * \brief Calcule le tenseur de 2 vecteurs
736 * @param *a le 1er vecteur (*double)
737 * @param N la taille du premier vecteur (int)
738 * @param *b le 2ieme vecteur (*double)
739 * @param M la taille du 2ieme vecteur (int)
740 * @param *ab le resultat
743 void Polynoms::tensor
744 ( const double *a, //vecteur gauche
745 int N, //taille du vecteur gauche
746 const double *b, //vecteur droit
747 int M, //taille du vecteur droit
748 double *ab //conteneur de sortie
753 ab[i*M+j] = a[i]*b[j];
756 /** \fn matrixProduct
757 * \brief calcule le produit de 2 matrices .
758 * @param *Matrix1 la 1ere matrice
759 * @param R1,C1 la taille de la matrice1 Nbr de ligne et le nbr de colonne
760 * @param Matrix2 la 2ieme matrice
761 * @param R2,C2 la taille de la matrice2 Nbr de ligne et le nbr de colonne
762 * @param out le resultat de Matrix1*Matrix2
764 void Polynoms::matrixProduct
766 const double *Matrix1,
769 const double *Matrix2,
775 for(int i=0; i<R1; i++)
777 for(int j=0; j<C2; j++)
780 for(int k=0; k<C1; k++)
782 out[i*C2+j] += Matrix1[i*C1+k]*Matrix2[k*C2+j];
787 /** \fn matrixProdVec
788 * \brief Calcule le produit Matrice vecteur
789 * @param *Matrix correspond à la matrice du produit
790 * @param R1,C1 la taille de la matrice( Nbr de ligne et le nbr de colonne)
791 * @param Vector correspond au vecteur du produit
792 * @param out le resultat de Matrix*Vector
794 void Polynoms::matrixProdVec
795 ( const double *Matrix,
798 const double *Vector,
802 for(int i=0; i<R1; i++)
805 for(int j=0; j<C1; j++)
807 out[i] += Matrix[i*C1+j]*Vector[j];
812 /** \fn shift_diagonal
813 * \brief Calcule la trace d'une matrice carrée
814 * @param *Matrix correspond à la matrice
815 * @param N taille de la matrice
816 * @param shift résultat
817 * @return renvoie la trace de la matrice Matrix dans shift */
818 void Polynoms::shift_diagonal( double *Matrix, const int N, double shift)
820 for (int i=0; i<N; i++)
821 Matrix[i*N + i] += shift;
823 //qques outils d'algebre des polynomes
825 /** \fn remplir_trinome
826 * \brief remplie un trinome
828 void Polynoms::remplir_trinome(double coeff, double u, double alpha, vector< double >& trinome)
832 trinome[1] = -2*coeff;
834 trinome[0] = coeff - alpha;
839 * \brief Additionne 2 vecteurs de tailles differentes
840 * @param P,Q les 2 vecteurs à additionner
841 * @param n,m les 2 tailles respectives des vecteurs P et Q
842 * @return un vecteur qui correspond a P+Q*/
843 vector< double > Polynoms::additionne(const vector< double >& P, const vector< double >& Q, int n, int m)
846 vector< double > P_plus_Q = P;
847 for(int j = 0; j<m+1; j++){
853 vector< double > P_plus_Q = Q;
854 for(int j = 0; j<n+1; j++){
863 * \brief Calcule le produit scalaire de 2 vecteurs de tailles differentes
864 * @param P,Q les 2 vecteurs
865 * @param n,m les 2 tailles respectives des 2 vecteurs
866 * @return un vecteur correspond à (P,Q) */
867 vector< double > Polynoms::multiplie(const vector< double >& P, const vector< double >& Q, int n, int m)
869 vector< double > P_fois_Q (n+m+1,0.);
870 for(int i = 0; i<n+1; i++){
871 for(int j = 0; j<m+1; j++){
872 P_fois_Q[i+j] += P[i]*Q[j];
879 /** \fn polynome_caracteristique
880 * \brief Calcule le polynome caracteristique
882 vector< double > Polynoms::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!
884 vector< double > trinome_u(3);
885 vector< double > trinome_c(3);
886 vector< double > pol_char;
888 remplir_trinome(alpha_v*cv_2,Uv,alpha_v,trinome_c);
889 remplir_trinome(rhol,Ul,dPil,trinome_u);
891 pol_char = multiplie(trinome_u, trinome_c, 2, 2);
893 remplir_trinome(alpha_l*cl_2,Ul,alpha_l,trinome_c);
894 remplir_trinome(rhov,Uv,dPiv,trinome_u);
896 pol_char = additionne(pol_char,multiplie(trinome_u, trinome_c, 2, 2),4,4);
901 /** \fn polynome_caracteristique
918 vector< double > Polynoms::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)
920 vector< double > trinome1(3);
921 vector< double > trinome2(3);
922 vector< double > pol_char;
925 /* (x-u1)^2(x-u2)^2-K1(x-u2)^2-K2(x-u1)^2+K3 */
926 K1=alpha1*rho2*g2press + dpi1*alpha2*invc2sq*g2alpha;
927 K2=alpha2*rho1*g2press + dpi2*alpha1*invc1sq*g2alpha;
928 K3=(alpha1*dpi2+alpha2*dpi1)*g2alpha*g2press/g2;
929 //cout<<"K1= " <<K1<<", K2= " <<K2<<", K3= " <<K3<<endl;
931 if(fabs(K1)>epsilon && fabs(K2)>epsilon )
933 remplir_trinome(1/K1,u1,1,trinome1);
934 //cout<<"coeff constant trinome1= " << trinome1[0]<<endl;
935 remplir_trinome(1/K2,u2,1,trinome2);
936 //cout<<"coeff constant trinome2= " << trinome2[0]<<endl;
938 pol_char = multiplie(trinome1, trinome2, 2, 2);
939 //cout<<"coeff constant produit= " << pol_char[0]<<endl;
941 pol_char[0]+=K3/K2/K1-1;
942 //cout<<"coeff constant apres soustraction= " << pol_char[0]<<endl;
946 remplir_trinome(1,u1,K1,trinome1);
947 remplir_trinome(1,u2,K2,trinome2);
949 pol_char = multiplie(trinome1, trinome2, 2, 2);
951 pol_char[0]+=K3-K2*K1;
954 /* for(int ct =0; ct<pol_char.size()-1; ct++) */
955 /* if(fabs(pol_char[ct])< epsilon) */
956 /* pol_char[ct]=0.; */
961 //Pour le tri des valeurs propres
964 * \brief calcule le module d'un nombre complexe
965 * @param z est un nombre complexe
966 * @return calcule le module de z*/
967 double Polynoms::module(complex< double > z)
972 /** \fn modulec calcule le module² de z */
973 double Polynoms::modulec(complex< double > z)
978 /** \fn abs_generalise
979 * \brief calcule la valeur absolue d'un nombre complexe
980 * \Details calcule la valeur absolue d'un nombre complexe en prenant en compte
981 * que la partie réelle c-a-d si z= a+ib abs_generalize() renvoie |a|+ib
983 * @return si z = a+ib elle renvoie |a|+ib */
984 complex< double > Polynoms::abs_generalise(complex< double > z)
991 void Polynoms::ordre_croissant_abs(vector< complex< double > > &L, int n)
993 vector< complex< double > > copieL = L;
998 for(int j=0; j<n; j++){
1000 Lmaxc = copieL[0].real();
1001 for (int i=1; i<n; i++){
1002 Lmax_iterc = copieL[i].real();
1003 if( Lmax_iterc > Lmaxc ){
1008 L[n-1-j] = copieL[i_max];
1009 copieL[i_max] = -INFINITY;
1013 int Polynoms::new_tri_selectif(vector< complex< double > > &L, int n, double epsilon)
1017 cout<<"avant ordre croissant"<<endl;
1018 for(int i=0;i<n;i++)
1019 cout<<L[i]<<", "<<endl;
1021 ordre_croissant_abs(L,n);
1023 cout<<"après ordre croissant"<<endl;
1024 for(int i=0;i<n;i++)
1025 cout<<L[i]<<", "<<endl;
1027 vector< complex< double > >result(1,L[0]);
1028 for(int i=1;i<n;i++)
1029 if(fabs(L[i].real()-result[size-1].real())>epsilon || fabs(L[i].imag())>epsilon)
1031 result.push_back(L[i]);
1034 cout<<"tri selectif step "<< i<<endl;
1035 for(int j=0;j<size;j++)
1036 cout<<result[j]<<", "<<endl;
1043 /** \fn tri_selectif
1049 template< typename T >
1050 int Polynoms::tri_selectif(T& L, int n, double epsilon)
1056 complex< double > Lmin;
1058 for(i = 0 ; i<size_vp ; i++){//On travaille sur le sous vecteur des indices de i � size_vp-1
1062 for(j = i+1 ; j < size_vp ; j++)
1063 if( Lmin.real() > L[j].real() )
1064 {Lmin = L[j]; imin = j;}
1066 //on met la vp � zero si elle est trop petite
1067 if( fabs(Lmin.real()) < epsilon) Lmin.real()=0;
1068 if( fabs(Lmin.imag()) < epsilon) Lmin.imag()=0;
1079 if ( module(L[i-1] - Lmin) > epsilon ) {
1086 L[imin]=L[size_vp - 1];
1087 L[size_vp - 1]=Lmin;
1094 if ( (module(L[i-1] - Lmin) > epsilon) && (module(L[i-2] - Lmin) > epsilon) ) {
1101 L[imin]=L[size_vp - 1];
1102 L[size_vp - 1]=Lmin;
1113 //Calcul des coefficients du polynome d'interpolation x->y dans la base de Newton
1123 vector< complex< double > > Polynoms::dif_div(int n, const vector< complex< double > >& x, const vector< T >& y, double epsilon)//n=nb valeurs à interpoler
1125 complex< double > tab_dif_div[n*n];//La matrice de taille n*n
1126 vector< complex< double > > dif_div (n);
1128 complex< double > delta_x;
1130 for( i=0 ; i<n ; i++) {
1131 tab_dif_div[i*n]=y[i];
1132 // cerr << tab_dif_div[i*n+0] << endl;
1134 for(j=1; j<n ; j++){
1135 for(i=0;i < n-j;i++) {
1136 delta_x = x[i+j]-x[i];
1137 if(module(delta_x) < epsilon) {
1138 cout << " attention, dif_div, delta_x <epsilon, " << " delta_x= " << delta_x <<" epsilon= "<<epsilon<<endl;
1139 cout<<" vp= " <<endl;
1140 for(int k=0; k<n; k++)
1141 cout << x[k]<<", " << endl;
1143 tab_dif_div[i*n+j]=(tab_dif_div[(i+1)*n+j-1]-tab_dif_div[i*n+j-1])/delta_x;
1146 for(i=0 ; i<n ; i++)
1147 dif_div[i] = tab_dif_div[i*n+n-1-i];
1150 //attention, Les vp complexes conjugu�es doivent se suivre dans la liste x
1151 void Polynoms::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
1154 double Id[sizeA*sizeA], aux[sizeA*sizeA], matProd[sizeA*sizeA];
1156 for(int k=0; k<sizeA*sizeA; k++)
1158 for(int k=0; k<sizeA; k++)
1160 /* for(int k=0; k<sizeA*sizeA; k++) */
1161 /* aux[k] = A[k] ; */
1162 //cerr << x << endl;
1163 if( fabs(x[0].imag()) > epsilon){
1164 for(int k=0; k<sizeA*sizeA; k++)
1165 p[k] = A[k] * dif_div[0].real() ;
1166 shift_diagonal(p, sizeA,(dif_div[1] - dif_div[0]*x[1]).real());
1169 for(int k=0; k<sizeA*sizeA; k++)
1170 p[k] = Id[k] * dif_div[0].real() ;
1172 for(i=debut ; i<n ; i++){
1173 for(int k=0; k<sizeA*sizeA; k++)
1175 //cerr << " on traite la racine reele " << x[i] << endl;
1176 if ( fabs(x[i].imag()) < epsilon ) {
1177 shift_diagonal(aux, sizeA, (- x[i].real()));
1178 matrixProduct( p, sizeA, sizeA, aux, sizeA, sizeA, matProd);
1179 for(int k=0; k<sizeA*sizeA; k++)
1181 shift_diagonal(p, sizeA, dif_div[i].real());
1183 else{ //couple de conjugees
1184 if(fabs(x[i].imag() + x[i+1].imag()) < epsilon){
1185 matrixProduct( aux, sizeA, sizeA, aux, sizeA, sizeA, matProd);
1186 for(int k=0; k<sizeA*sizeA; k++)
1187 aux[k] = matProd[k];
1188 for(int k=0; k<sizeA*sizeA; k++)
1189 aux[k]+=(-2*x[i].real()) * A[k];
1190 shift_diagonal(aux, sizeA, modulec(x[i]));
1192 matrixProduct( p, sizeA, sizeA, aux, sizeA, sizeA, matProd);
1193 for(int k=0; k<sizeA*sizeA; k++)
1195 for(int k=0; k<sizeA*sizeA; k++)
1196 p[k]+= dif_div[i].real() * A[k];
1197 shift_diagonal(p, sizeA, (dif_div[i+1] - dif_div[i]*x[i+1]).real());
1201 cout << " appliquer_dif_div : les racines complexes conjuguees ne se suivent pas " << endl;
1202 for(int k=0; k<n; k++)
1203 cout << x[k]<<", " << endl;
1204 throw(" appliquer_dif_div : les racines complexes conjuguees ne se suivent pas ");
1210 void Polynoms::abs_par_interp_directe(int n, vector< complex< double > > vp, double * A, int sizeA, double epsilon, double *result,vector< complex< double > > y)
1212 vector< complex< double > > differences_divisees = dif_div(n, vp, y, epsilon);
1214 cout<<"valeurs propres"<<endl;
1215 for( int i=0 ; i<n ; i++)
1216 cout<<vp[i] <<", "<<endl;
1217 cout<<"valeurs à atteindre"<<endl;
1218 for( int i=0 ; i<n ; i++)
1219 cout<<y[i] <<", "<<endl;
1220 cout<<"differences_divisees= "<<endl;
1221 for(int i=0 ; i<n ; i++)
1222 cout<<differences_divisees[i]<<", ";
1225 appliquer_dif_div(n, differences_divisees, vp, A, sizeA, epsilon, result);
1228 bool Polynoms::belongTo(complex< double > a , vector< complex <double > > v, double epsilon)
1230 bool result = false;
1231 int size = v.size();
1233 double epsilon2 = epsilon*epsilon;
1235 while(i<size && !result)
1237 result = modulec(a-v[i]) < epsilon2;
1242 double Polynoms::avr(double a, double b)