synaps/mpol/elim.h

00001 /*********************************************************************
00002 *      This file is part of the source code of SYNAPS.               *
00003 *      V.t TRAN , GALAAD, INRIA                                      *
00004 **********************************************************************
00005 $Id: Elim.h,v 1.1 2005/07/11 08:15:51 mourrain Exp $
00006 \file 
00007 **********************************************************************/
00008 #ifndef SYNAPS_RESULTANT_ELIM_H
00009 #define SYNAPS_RESULTANT_ELIM_H
00010 
00011 #include <synaps/init.h>
00012 #include <synaps/linalg.h>
00013 #include <synaps/linalg/SousMatMax.h>
00014 #include <synaps/mpol/macaulay.h>
00015 #include <synaps/mpol/bezoutian.h>
00016 
00017 __BEGIN_NAMESPACE_SYNAPS
00018 
00019 
00023 template<class LL>
00024 int numvar(const LL l)
00025 {
00026   int nv=0, n=l.size();
00027   for (int i=0;i<n;i++) {
00028     int lv = MPOLDST::lvar(l[i]);
00029     nv = (nv > lv) ? nv : lv;
00030   }
00031   return ++nv;
00032 };
00033 
00040 template<class LL>
00041 LL transformType(LL h) {
00042   typedef typename LL::value_type POL;
00043   typedef typename POL::monom_t MON;
00044 
00045   LL f;
00046   int i, n = h.size()-1;
00047   char xi[20];// = new char(10);
00048 
00049   //number of variables
00050   int nv=numvar(h);
00051 
00052   //recharging values of the variables in x
00053   /*  for (i=0;i<n;i++) {
00054     sprintf(xi,"x%d",i); 
00055     strcat(xi,itoa(i));
00056     POL::var(i+1,xi);
00057   }
00058   */
00059   //  delete xi;
00060 
00061   for (i=0;i<=n;i++)
00062     for (int j=0;j<nv;j++)
00063       h[i] = MPOLDST::swap(h[i],j,j+n);
00064 
00065   POL h0 = h[0];
00066   for (i=0; i<n; i++) {
00067     POL p = h0*MON(i,1) - h[i+1];//pn = h0*xi - hi
00068     f.push_back(p);
00069   }
00070 
00071   return f;
00072 };
00073 
00074 
00083 template <class LL>
00084 LL hiddenVar(const LL & f,  const int  l0, const int l1,  LL& f2)
00085 {
00086   typedef typename LL::value_type POL;
00087   typedef typename POL::monom_t MON;
00088   typedef typename POL::coeff_t TT;
00089  
00090   LL f1, f3, f4;
00091   int N=1; 
00092   int i,j, flag;
00093   typename LL::const_iterator ll,ll1;
00094   typename POL::const_iterator pp,pp1;
00095   
00096   for (ll=f.begin();ll!=f.end();ll++) {
00097     POL tmp1 = (*ll);
00098     for (j=0;j<l0;j++)
00099       tmp1 = MPOLDST::subs(tmp1,j,1); //Replace the @l0@ first variables by 1
00100 
00101     for (pp=(*ll).begin();pp!=(*ll).end();pp++) {
00102       POL ct(*pp);
00103       if (MPOLDST::degree(ct)==0) {
00104                 int perd =1;
00105                 for (pp1=tmp1.begin();pp1!=tmp1.end();pp1++) {
00106                   POL rc(*pp1);
00107                   if (MPOLDST::degree(rc)==0){
00108                         perd = 0;
00109                         break;
00110                   }
00111                 }
00112                 if (perd == 1)
00113                   tmp1+=POL("1");
00114                 break;
00115       }
00116     }
00117     
00118     POL tmp2("0");
00119     for (pp=tmp1.begin(); pp!=tmp1.end();pp++) {
00120       MON mm = (*pp);
00121       mm.setcoeff(TT(1));  //set coefficients to 1
00122       tmp2 +=mm;
00123     }
00124     f1.push_back(tmp2);
00125   }
00126  
00127   f2.push_back(POL("0"));//value NUL of Macaulay matrix if it exists
00128 
00129   for (ll=f.begin(),ll1=f1.begin();ll!=f.end();ll++,ll1++) {
00130     POL pol("0");
00131     for (pp=(*ll1).begin(); pp!=(*ll1).end();pp++) {
00132       POL ini("0");//Init
00133       for (typename POL::const_iterator kk=(*ll).begin(); kk!=(*ll).end();kk++) {
00134         POL tmp(*kk), div(*pp);
00135         tmp = tmp/div;
00136         if (tmp!=0) {//divisible
00137           flag=1;
00138           for (j=l0;j<l0+l1;j++) {
00139             POL con(MON(j,1)); 
00140             con = tmp/con;
00141             if (con != 0) 
00142               flag=0;
00143           }
00144           if (flag==1)//yes
00145             ini = ini + tmp;//insert
00146         }
00147       }
00148       f2.push_back(ini); 
00149       
00150       MON mm =(*pp);
00151       mm.setcoeff(TT(N));//Identify the coefficient of monomial mm by an integer N
00152       pol +=mm;
00153       N++;
00154     }
00155     f3.push_back(pol);
00156   }
00157  
00158   //Shift the variables has x0,..,xl1
00159   for (ll=f3.begin();ll!=f3.end();ll++) {
00160     POL tmp = (*ll);
00161     for (j=0;j<l1;j++)
00162       tmp=MPOLDST::swap(tmp,j,j+l0);
00163     f4.push_back(tmp);
00164   }
00165 
00166   return f4;
00167 };
00168 
00173 template<class MAT, class LL>
00174 MAT Macaulay(const LL & f, const int  l0, const int l1)
00175 {
00176   typedef typename LL::value_type POL;
00177   LL f1,f2;
00178   f1 = hiddenVar<LL>(f,l0,l1,f2); //Hide the l0 firsts variables in f2
00179  
00180   //Compute the Macaulay matrix with f1.
00181   MAT M=macaulay<MAT>(f1);
00182   
00183   //Find the coefficients according to the hidden variables stored in f2
00184   int mc;
00185   for (int i=0;i<M.nbrow(); i++)
00186     for (int j=0;j<M.nbcol();j++) {
00187       POL pm = M(i,j);
00188       if (pm==0)
00189                 mc =0;
00190       else
00191                 mc = int((pm.begin()->coeff()).get_d());
00192       matrixof::assigncoeff(M,i,j,f2[mc]);
00193     }
00194   
00195   return M;  
00196 };
00197 
00201 template<class MAT, class LL>
00202 MAT Minor(const MAT& mat, const LL & f)
00203 {
00204   typedef typename LL::value_type POL;
00205   typedef typename POL::monom_t MON;
00206   typedef typename POL::order_t O;
00207   int d[f.size()], n = f.size()-1, nu = 0;
00208   POL F;
00209 
00210   typename LL::const_iterator fi;int i;
00211   for(i=0, fi=f.begin(); i<=n; i++, fi++) {
00212     d[i] = std::max(MPOLDST::degree(*fi),1);
00213     nu += d[i];
00214   }
00215   nu -= n; 
00216 
00217   MON* xdi = new MON[n+1];
00218   POL* B   = new POL[n+1];
00219 
00220   xdi[0]= MON(n+1,d[0]);
00221   B[0]  = POL(1);
00222   for(int j=1; j<=n; j++){
00223     xdi[j] = MON(j-1,d[j]);
00224     B[j]=0;
00225     for(int i=0;i<d[j];i++) 
00226       B[j]+=POL(MON(j-1,i));
00227   }
00228 
00229   POL* E= new POL[n+1];
00230   E[0]=1;
00231 
00232   POL U(1), V(1);
00233   int mu;
00234     
00235   for(int i=n;i>0;i--) {
00236     mu = nu-d[i];
00237     for (typename POL::const_iterator j = E[0].begin(); j!=E[0].end(); j++) {
00238       choose(mu- degree(*j),i,V);
00239       E[i] += V*(*j);
00240     }
00241     E[0]*=B[i];
00242   } 
00243   for (int i=1;i<=n;i++)
00244     E[i] *= xdi[i];
00245 
00246   //Index of the monomial
00247   std::map<MON,int,O> index;
00248   int c=0;
00249   MON m;
00250   for(int i=0; i<=n; i++)
00251     for (typename POL::reverse_iterator j = E[i].rbegin(); j!=E[i].rend();j++){
00252       m = (*j);
00253       index[m]=c;
00254       c++;
00255     }
00256 
00257   //Construction of the dodus monomials
00258   POL DO("0"); 
00259   for (int i=0;i<n;i++) {
00260     for (int j=i+1;j<=n;j++) {
00261       POL tmp = xdi[i];
00262       tmp = tmp*xdi[j];
00263       DO += tmp;//xdi[i]*xdi[j];
00264     }
00265   }
00266 
00267   int v=0; 
00268   std::vector<int> ind;
00269 
00270   for (int i=1;i<=n;i++) //for each E[i]
00271     for (typename POL::const_iterator ee=E[i].begin();ee!=E[i].end();ee++) {//for each monomial of E[i]
00272       POL tmp1(*ee);// tmp1*=xdi[i];
00273       tmp1 *= MON(n+1,nu-degree(*ee));
00274       for (typename POL::const_iterator pp=DO.begin();pp!=DO.end();pp++) {
00275                 POL tmp2(*pp), tmp3 = tmp1/tmp2;
00276                 if (tmp3!=0) {//pp in tmp1
00277                   ind.push_back(index[(*ee)]);
00278                   v++;
00279                   break;
00280                 }
00281       }
00282     }
00283 
00284  //Construction of maximum minor of Macaulay matrix
00285   MAT Res;
00286   if (v<2) {//Not exists
00287     POL p("1");
00288     matrixof::reserve(Res,1,1,1);
00289     matrixof::assigncoeff(Res,0,0,p);
00290   }
00291   else {
00292     matrixof::reserve(Res,v,v,1);
00293     for (int i=0;i<v;i++)
00294       for (int j=0;j<v;j++)
00295                 matrixof::assigncoeff(Res,i,j,mat(ind[i],ind[j]));
00296   }
00297 
00298   return Res;
00299 };
00300 
00301 
00307 //Applicable with QQ number
00308 template<class LL>
00309 typename LL::value_type Elim_Bezout(const LL & f) 
00310 {
00311   typedef typename LL::value_type POL;
00312   typedef typename POL::coeff_t TT;
00313   typedef MatrDse<POL, linalg::rep2d<POL> > MAT;
00314 
00315   int n = f.size();
00316   
00317   MAT B = mbezout<MAT>(f,n);
00318 #ifdef DEBUG
00319   std::cout<<"\nBezout matrix B :\n"<<B;
00320 #endif 
00321   VectDse<TT> v(n); 
00322   for (int i=0; i<n;i++)
00323     v[i] = rand();
00324 
00325   MatrDse<TT, linalg::rep2d<TT> > BS(B.nbrow(), B.nbcol());
00326   for (int i=0; i<B.nbrow(); i++) 
00327     for (int j=0; j<B.nbcol(); j++) {
00328       POL p = B(i,j);
00329       TT tmp = p(v.rep());
00330       BS(i,j) = tmp;
00331     }
00332 
00333   MAT SM = SousMatMax(BS, B);
00334 #ifdef DEBUG
00335   std::cout<<"\nMinor of B : \n"<<SM<<std::endl;
00336 #endif
00337   POL result = Det(SM);//Berko(),SM);
00338 
00339   return result;
00340 };
00341 
00347 template<class LL>
00348 typename LL::value_type Elim_Macaulay(const LL & f)
00349 {
00350   typedef typename LL::value_type POL;
00351   typedef MatrDse<POL, linalg::rep2d<POL> > MAT;
00352 
00353   int nv =numvar(f)-f.size()+1;
00354   //  LL f = transformType(h);
00355 
00356   int n = f.size();
00357   MAT M= Macaulay<MAT,LL>(f,n,nv);
00358 #ifdef DEBUG
00359   std::cout<<"Macaulay matrix M :\n"<<M<<std::endl;
00360 #endif
00361   LL f1,f2;
00362   f1 = hiddenVar(f,n,nv,f2);
00363   MAT m = Minor<MAT>(M,f1);
00364 #ifdef DEBUG
00365   std::cout<<"Minor matrix m :\n"<<m<<std::endl;
00366 #endif
00367   POL d = Det(M);
00368   POL dm = Det(m);
00369 #ifdef DEBUG
00370   std::cout<<"Determinant of M : "<<d<<"\nDeterminant of m : "<<dm<<std::endl;
00371 #endif
00372   if (dm!=0) d /=dm;
00373 
00374   return d;
00375 };
00376 
00377 __END_NAMESPACE_SYNAPS
00378 
00379 #endif // SYNAPS_RESULTANT_ELIM_H
00380 

SYNAPS DOCUMENTATION
logo