synaps/upol/resultant.h

Go to the documentation of this file.
00001 /*********************************************************************
00002 *      This file is part of the source code of SYNAPS kernel.        *
00003 *   Author(s): B. Mourrain, GALAAD, INRIA                            *
00004 *********************************************************************/
00005 #ifndef SYNAPS_RESULTANT_UNIRESLT_H
00006 #define SYNAPS_RESULTANT_UNIRESLT_H
00007 //====================================================================
00008 #include <complex>
00009 #include <synaps/init.h>
00010 #include <synaps/upol/UPolDse.h>
00011 #include <synaps/linalg/MatrDse.h>
00012 #include <synaps/linalg/Det.h>
00013 #include <synaps/linalg/hankel.h>
00014 __BEGIN_NAMESPACE_SYNAPS
00015 //====================================================================
00018 //====================================================================
00020 
00025 template <class MATR, class R>
00026 MATR bezout(const UPolDse<R> & P, const UPolDse<R> & Q) 
00027 {
00028   unsigned int l = std::max(degree(P),degree(Q));
00029   MATR M(l,l);
00030   if(degree(P)<degree(Q)){
00031     for(int i=0; i<degree(P)+1; i++){
00032       for( int j=0; j<degree(Q)-i; j++)
00033         for( int l=0; l<=j; l++)
00034           M(i+j-l,i+l) += P[i]*Q[j+i+1];
00035       for( int j=0; j<i; j++)
00036         for( int l=0; l<=j; l++)
00037           M(i-1-l,i-1-j+l) -= P[i]*Q[i-1-j];
00038     }
00039   } else {
00040     for( int i=0; i<degree(Q)+1; i++){
00041       for( int j=0; j<degree(P)-i; j++)
00042         for( int l=0; l<=j; l++)
00043           M(i+j-l,i+l) -= Q[i]*P[j+i+1];
00044       for( int j=0; j<i; j++)
00045         for( int l=0; l<=j; l++)
00046           M(i-1-l,i-1-j+l) += Q[i]*P[i-1-j];
00047     }
00048   }
00049   return (M);
00050 } 
00051 //------------------------------------------------------------------------
00053 template <class R>
00054 MatrStr<linalg::hankel<typename R::value_type> > bezout(const R & P)
00055 {
00056   typedef typename R::value_type C;
00057   typedef MatrStr<linalg::hankel<C> >  Mat;
00058   Mat H(degree(P),degree(P));
00059   for(int i=0; i<degree(P); i++) {
00060     (H.rep())[i]=P[i+1];
00061     if (i!=degree(P)-1) (H.rep())[degree(P)+i]=(C)0;
00062   }
00063   return H;
00064 }
00065 //-------------------------------------------------------------------------
00067 
00072 template <class MATR, class R>
00073 MATR sylvester(const UPolDse<R> & P, const UPolDse<R> & Q) 
00074 {
00075   int l = degree(P)+degree(Q); 
00076   MATR M(l,l);
00077   for( int i=0; i<l; i++)
00078     for( int j=0; j<l; j++) M(i,j)=0;
00079   for( int i=0; i< degree(Q); i++)
00080     for( int j=0; j< degree(P)+1; j++)
00081       M(j+i,i)=P[j];
00082   for( int i=0; i< degree(P); i++)
00083     for( int j=0; j< degree(Q)+1; j++)
00084       M(j+i,degree(Q)+i)=Q[j];
00085 
00086   return (M);
00087 }
00088 //-------------------------------------------------------------------------
00090 template <class C, class R, class MATR>
00091 MATR sylvester(const UPolDse<C,R> & P, const UPolDse<C,R> & Q, char c)
00092 {
00093   int dp = degree(P), dq = degree(Q);
00094   int l = dp+dq;
00095   MATR M(l,l);
00096   for( int i=0; i<l; i++)
00097     for( int j=0; j<l; j++) M(i,j)=0;
00098   for( int i=0; i< dq; i++)
00099     for( int j=0; j< dp+1; j++)
00100       M(j+i,i)=P[dp-j];
00101   for( int i=0; i< dp; i++)
00102     for( int j=0; j< dq+1; j++)
00103       M(j+i,dq+i)=Q[dq-j];
00104 
00105   return (M);
00106 }
00107 
00108 //-------------------------------------------------------------------------
00116 #if 0
00117 template<class UPOL, class MPOL>
00118 MatrDse<UPOL> u_sylvester(const MPOL & p1, const MPOL & p2)
00119 {
00120    int d1 = degree(p1,1), d2=degree(p2,1);
00121    MatrDse<UPOL> A(d1+d2,d1+d2);
00122 
00123    for(typename MPOL::const_iterator it=p1.begin();it!=p1.end();it++)
00124    {
00125       UPOL m(it->coeff(),(*it)[0]);
00126       for(int k=0;k<d2;k++) A((*it)[1]+k,k) +=m;
00127    }
00128    for(typename MPOL::const_iterator it=p2.begin();it!=p2.end();it++)
00129    {
00130       UPOL m(it->coeff(),(*it)[0]);
00131       for(int k=0;k<d1;k++) A((*it)[1]+k,k+d2) +=m;
00132    }
00133 
00134    return A;
00135 }
00136 #endif
00137 //#if 0
00138 template<class UPOL, class MPOL>
00139 MatrDse<UPOL> u_sylvester(const MPOL & p1, const MPOL & p2,
00140         typename MPOL::monom_t::index_t i = 1)
00141 {
00142   typename MPOL::monom_t::index_t j = (i == 0 ? 1 : 0);
00143 
00144   int d1 = degree(p1, i), d2 = degree(p2, i);
00145   MatrDse<UPOL> A(d1+d2, d1+d2);
00146 
00147   for(typename MPOL::const_iterator it=p1.begin();it!=p1.end();it++)
00148   {
00149         UPOL m(it->coeff(),(*it)[j]);
00150 //      cout << *it << " " << m << endl;
00151         //for(int k=0;k<d2;k++) A((*it)[0]+k,k) +=m;
00152         for(int k=0;k<d2;k++) A(k, k+ (*it)[i]) +=m;
00153   }
00154 //  cout << "---" << endl;
00155   for(typename MPOL::const_iterator it=p2.begin();it!=p2.end();it++)
00156   {
00157         UPOL m(it->coeff(),(*it)[j]);
00158 //      cout << m << endl;
00159         //for(int k=0;k<d1;k++) A((*it)[0]+k,k+d2) +=m;
00160         for(int k=0;k<d1;k++) A(d2 + k, k+ (*it)[i]) +=m;
00161   }
00162 
00163   return A;
00164 }
00165 //#endif
00166 //-------------------------------------------------------------------------
00168 template <class MATR, class R>
00169 MATR companion(const UPolDse<R> & P)
00170 {
00171   typedef typename UPolDse<R>::size_type size_type;
00172   size_type n = degree(P);
00173   MATR M(n,n);
00174   for(size_type i=0; i<n; i++)
00175     for(size_type j=0; j<n-1; j++)
00176       if(i==j+1) M(i,j)=1; else M(i,j)=0;
00177   for(size_type i=0; i<n; i++) M(i,n-1)=-(P[i]/P[n]);
00178 
00179   return M;
00180 }
00181 //--------------------------------------------------------------------
00187 template<class MPOL>
00188 MatrDse<MPOL>  u_sylvester(const MPOL & p1, const MPOL & p2,int I)
00189 {
00190   typedef typename MPOL::monom_t monom_t;
00191   int d1 = degree(p1,I), d2=degree(p2,I);
00192   MatrDse<MPOL> A(d1+d2,d1+d2);
00193   for(typename MPOL::const_iterator it=p1.begin();it!=p1.end();it++)
00194     {
00195       int d=(*it).nvars();
00196       MPOL P("1;");
00197       if(d<0)
00198         {
00199           P=(*it);
00200         }
00201       else
00202         {
00203           for(int j=0;j<=d;j++)
00204             {
00205               MPOL P1;
00206               if(j!=I)
00207                 {
00208                   int power=(*it)[j];
00209                   if(power!=0)
00210                     {
00211                       P1= monom_t(1,j,power);
00212                       P*=P1; 
00213                     }
00214                 }
00215               else if(j==I && (*it)[I]>0)
00216               {
00217                 P1=monom_t(1,0,0);
00218                 P*=P1;
00219               }
00220             }
00221           P*=(it->coeff());
00222         }
00223       for(int k=0;k<d2;k++)
00224         {
00225           A((*it)[I]+k,k) +=P;
00226         }
00227     }
00228   for(typename MPOL::const_iterator it=p2.begin();it!=p2.end();it++)
00229     {
00230       int d=(*it).nvars();
00231       MPOL P("1;");
00232       if(d<0)
00233         {
00234           P=(*it);
00235         }
00236       else
00237         {
00238           for(int j=0;j<=d;j++)
00239             {
00240               MPOL P1;
00241               if(j!=I)
00242                 {
00243                   int power=(*it)[j];
00244                   if(power!=0)
00245                     { 
00246                       P1 = monom_t(1,j,power);
00247                       P*=P1;
00248                     }
00249                 }
00250               else if(j==I && (*it)[I]>0)
00251               {
00252                 P1=monom_t(1,0,0);
00253                 P*=P1;
00254               }
00255             }
00256           P*=(it->coeff());
00257         }
00258       for(int k=0;k<d1;k++) 
00259         {
00260           A((*it)[I]+k,k+d2) = A((*it)[I]+k,k+d2)+P;
00261         }
00262     }
00263   return A;
00264 }
00265 
00266 //--------------------------------------------------------------------
00268 
00269 //pour l'instant elle ne calcule le resultant que losque p1 et p2 sont trois variable (a generaliser)
00270 
00271 template<class MPOL>
00272 MPOL  Resultant(const MPOL & p1, const MPOL & p2,int I)
00273 {
00274   typedef typename MPOL::monom_t monom_t;
00275   MatrDse<MPOL> M=u_sylvester(p1,p2,I);
00276   return Det(M);
00277 }
00278 
00279 //--------------------------------------------------------------------
00280 template<class C, class O, class R> struct MPol;
00281 //--------------------------------------------------------------------
00284 template<class C, class O, class R>
00285 C  Subresultant(const MPol<C,O,R> & p1, const MPol<C,O,R> & p2,int I, int Ord,int Deg, C x,C y) 
00286 {
00287   int d1= degree(p1,I), d2=degree(p2,I);
00288   int taille = d1+d2-2*Ord;
00289   if (taille <1)
00290     {
00291       std::cerr<<"unable computing Subresultant"<<std::endl; //a voir proprement
00292       return 0.;
00293     }
00294   else
00295     {
00296       MatrDse<MPol<C,O,R> > A(d1+d2-2*Ord,d1+d2-Ord);
00297       typedef typename MPol<C,O,R>::monom_t monom_t;
00298       for(typename MPol<C,O,R>::const_iterator it=p1.begin();it!=p1.end();it++)
00299         {
00300           int d=(*it).nvars();
00301           MPol<C,O,R> P;
00302           if(d<0)
00303             {
00304               P=(*it);
00305             }
00306           for(int j=0;j<=d;j++)
00307             {
00308               MPol<C,O,R> P1;
00309               if(j!=I)
00310                 {
00311                   int power=(*it)[j];
00312                   if(power>0)
00313                     {
00314                       P1= monom_t(1,j,power);
00315                       P+=P1; 
00316                     }
00317                 }
00318               else if(j==I && (*it)[I]>0)
00319                 {
00320                   P1=monom_t(1,0,0);
00321                   P+=P1;
00322                 }
00323               P*=(it->coeff());
00324             }
00325           for(int k=0;k<d2-Ord;k++)
00326             {
00327               A(k,(*it)[I]+k) +=P;
00328             }
00329         }
00330       for(typename MPol<C,O,R>::const_iterator it=p2.begin();it!=p2.end();it++)
00331         {
00332           int d=(*it).nvars();
00333           MPol<C,O,R> P;
00334           if(d<0)
00335             {
00336               P=(*it);
00337             }
00338           else
00339             {
00340               for(int j=0;j<=d;j++)
00341                 {
00342                   MPol<C,O,R> P1;
00343                   if(j!=I)
00344                     {
00345                       int power=(*it)[j];
00346                       if(power!=0)
00347                         {
00348                           P1 = monom_t(1,j,power);
00349                           P+=P1;
00350                         }
00351                     }
00352                   else if(j==I && (*it)[I]>0)
00353                     {
00354                       P1=monom_t(1,0,0);
00355                       P+=P1;
00356                     }
00357                 }
00358               P*=(it->coeff());
00359             }
00360           for(int k=d2-Ord;k<d1+d2-2*Ord;k++)
00361             {
00362               A(k,k+(*it)[I]-d2+Ord) +=P;
00363             }
00364         }
00365       MatrDse<MPol<C,O,R> > Minor(d1+d2-2*Ord,d1+d2-2*Ord);
00366       for(int j=0;j<d1+d2-2*Ord-1;j++)
00367         {
00368           for(int i=0;i<d1+d2-2*Ord;i++)
00369             {
00370               Minor(i,j)=A(i,j);
00371             }
00372         }
00373       for(int i=0;i<d1+d2-2*Ord;i++)
00374         {
00375           Minor(i,d1+d2-2*Ord-1)=A(i,d1+d2-Ord-Deg-1);
00376         }
00377       MPol<C,O,R> subres=(Det(Berko(),Minor));
00378       C z=0.;
00379       int var1=0;
00380       for(int j=0;j<=2;j++)
00381         {
00382           if (j!=I)
00383             {
00384               var1=j;break;
00385             }
00386         }
00387    int var2=0;
00388       for(int j=0;j<=2;j++)
00389         {
00390           if (j!=I && j!=var1)
00391             {
00392               var2=j;break;
00393             }
00394         }
00395       for(typename MPol<C,O,R>::const_iterator it=subres.begin();it!=subres.end();it++)
00396         {
00397           int dx=(*it)[var1];
00398           int dy=(*it)[var2];
00399           C z1=pow(x,dx)*pow(y,dy);
00400           z1*=(*it).coeff();
00401           z+=z1;
00402         }
00403       return(z);
00404     }
00405 }
00406 //----------------------------------------------------------------------
00407 
00408 
00409 __END_NAMESPACE_SYNAPS
00410 
00411 #endif // SYNAPS_RESULTANT_UNIRESLT_H
00412 

SYNAPS DOCUMENTATION
logo