borderbasix: include/borderbasix/pol2ter.cc Source File

00001 #define polter pol
00002 #include"ugly.cc"
00003 
00004 int nbbits[256]={0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8};
00005 
00006 template<typename T> void *MAC_REV_MALLOC(int size);
00007 template<typename T> void *MAC_REV_REALLOC(void*ptr,int,int);
00008 template<typename T> void MAC_REV_FREE(void*ptr,int size);
00009 
00010 template<typename mon,typename T>
00011 struct polter  
00012 {
00013   typedef mon monom_t;
00014   typedef T coeff_t;
00015   mon ind;
00016   T * nf;
00017   unsigned *nfind;
00018   int size,sizenf;
00019   // pol operator*(const mon&t);//multi par un monome
00020   polter operator*(const T &t);
00021   polter &operator-=(const polter &Q);
00022   polter &operator=(const polter &Q);
00023   polter &operator/=(const T &t);
00024 };
00025 
00026 //
00027 //l'autre par un coeff constant
00028 //
00029 template<typename mon,typename T>
00030 polter<mon,T> 
00031 polter<mon,T>::operator*(const typename polter<mon,T>::coeff_t & toto)
00032 {
00033   int tmpsize=this->size;
00034   for(int i=0;i<this->size;i++)
00035     this->nf[i]*=toto;
00036   return *this;
00037 }
00038 
00039 template<typename mon,typename T>
00040 polter<mon,T>& polter<mon,T>::operator-=(const polter<mon,T> &Q)
00041 {
00042   int sizeres=0;
00043   polter<mon,T> res;
00044   if(sizenf<Q.sizenf)
00045     {
00046       T *ptrnf=nf,*ptrresnf,*ptrQnf;
00047       for(int i=0;i<sizenf;i++)
00048         sizeres+=nbbits[nfind[i]&Q.nfind[i]];
00049  
00050       for(int i=sizenf;i<Q.sizenf;i++)
00051         sizeres+=nbbits[Q.nfind[i]];
00052       res.size=sizeres;
00053       //res.nfind=malloc(Q.sizenf);
00054       res.nf=(T*)MAC_REV_MALLOC<T>(res.size*sizeof(T));
00055       ptrresnf.nf;
00056       ptrQnf=Q.nf;
00057       for(int i=0;i<sizenf;i++)
00058         {
00059           for(int j=0;j<8;j++)
00060             {
00061               if((nfind[i]>>j)&1)
00062                 *(ptrresnf++)=*(ptrnf++);
00063               if ((Q.nfind[i]>>j)&1)
00064                 *(ptrresnf-1)-=*(ptrQnf++);
00065             }
00066         }
00067       for(int i=sizenf;i<Q.sizenf;i++)
00068         {
00069           for(int j=0;j<8;j++)
00070             if((Q.nfind[i]>>j)&1)
00071               *(ptrresnf++)=((T)-1)*(*(ptrQnf++));
00072         }
00073     }
00074   else
00075     {
00076       T *ptrnf=nf,*ptrresnf,*ptrQnf;
00077       for(int i=0;i<Q.nizenf;i++)
00078         sizeres+=nbbits[nfind[i]&Q.nfind[i]];
00079       
00080       for(int i=Q.sizenf;i<sizenf;i++)
00081         sizeres+=nbbits[nfind[i]];
00082       
00083       res.size=sizeres;
00084       // res.nfind=malloc(sizenf);
00085       res.nf=(T*)MAC_REV_MALLOC<T>(res.size*sizeof(T));
00086       ptrresnf.nf;
00087       ptrQnf=Q.nf;
00088       for(int i=0;i<Q.sizenf;i++)
00089         {
00090           for(int j=0;j<8;j++)
00091             {
00092               if((nfind[i]>>j)&1)
00093                 *(ptrresnf++)=*(ptrnf++);
00094               if((Q.nfind[i]>>j)&1)
00095                 *(ptrresnf-1)-=*(ptrQnf++);
00096             }
00097         }
00098       for(int i=Q.sizenf;i<sizenf;i++)
00099         {
00100           for(int j=0;j<8;j++)
00101             if((nfind[i]>>j)&1)
00102               *(ptrresnf++)=*(ptrnf++);
00103         }
00104 
00105     }
00106   if(sizenf<Q.sizenf)
00107     nfind=(unsigned char*)realloc(nfind,Q.sizenf);
00108   for(int i=0;i<Q.sizenf;i++)
00109     nfind[i]|=Q.nfind[i];
00110   MAC_REV_FREE<T>(sizenf*sizeof(T));
00111   free(nfind);
00112   nfind=res.nfind;
00113   nf=res.nf;
00114   sizenf=res.sizenf;
00115   size=res.size;
00116   return *this;
00117 }
00118 //attention l'operateur -= n'agit que sur les nf
00119 
00120 //
00121 //operateur d'affectation sans recopie de nf!!!!
00122 //
00123 template<typename mon,typename T>
00124 polter<mon,T> & polter<mon,T>::operator= (const polter<mon,T>&Q)
00125 {
00126   nf=Q.nf;
00127   nfind=Q.nfind;
00128   size=Q.size;
00129   sizenf=Q.sizenf;
00130   ind=Q.ind;
00131   return *this;
00132 }
00133 
00134 template<typename mon,typename T>
00135 polter<mon,T> & polter<mon,T>::operator/= (const T&t)
00136 {
00137  for(int i=0;i<this->size;i++)
00138     this->nf[i]/=t;
00139   return *this;
00140 }
00141 
00142 template<typename polyalp, typename mon, typename T> 
00143 polyalp invconv(const polter<mon,T> & p)
00144 {
00145   typedef typename polyalp::order_t ord;
00146   polyalp res(0);
00147   T *tmpptr=p.nf;
00148   int counter=0,decalage=0;
00149   mon tmp;
00150   for(int i=0;i<p.sizenf;i++)
00151     {
00152       unsigned char motif=p.nfind[i]&255;
00153       unsigned nb_repete=p.nfind[i]>>8;
00154       decalage=0;
00155       for(int j=0;j<8;j++)
00156         if ((p.nfind[i]>>j)&1)
00157           {
00158             for(int k=0;k<nb_repete;k++)
00159               {
00160                 int2mon(counter+8*k+j,tmp);
00161                 res+=tmp*(*(tmpptr+decalage+k*nbbits[motif]));
00162               }
00163             decalage++;
00164           }
00165       tmpptr+=nbbits[motif]*nb_repete;
00166     }
00167   if (!Iszero(p.ind.GetCoeff()))
00168     res+=p.ind;
00169   //cout<<"res dans invconv"<<endl<<res<<endl;
00170   return res;
00171 }
00172 
00173 void compress_ind();
00174 
00175 // renvoit mon/x_i .... mais suppose qu'on ne fait pas n'importe quoi
00176 //necessaire pour convert
00177 
00178 template <typename Mon,typename Base>
00179 Mon divmon(const Mon &mon,int i,const Base &b)
00180 {
00181   Mon res(1);
00182   for(int j=0;j<b.nvar();j++)
00183     if(i==j)
00184       {
00185         res*=Mon(j,mon.GetDegree(j)-1);
00186       }
00187   else 
00188     {
00189       res*=Mon(j,mon.GetDegree(j));
00190     }
00191   return res;
00192 }
00193 
00194 //converti un polyalp en polter
00195 template<typename polter,typename polyalp,typename Base>
00196 polter convert(const polyalp &P,const Base &b)
00197 {
00198   typedef typename polyalp::order_t ord;
00199   typedef typename polyalp::monom_t mon;
00200   typedef typename mon::coeff_t coeff;
00201   coeff *tmp;
00202   polter res;
00203   int i,maxind=-1,j,p=0;
00204   i=0;
00205   res.ind=mon(0);
00206   for(typename polyalp::const_iterator iter=P.begin();iter!=P.end();iter++)
00207     if ((IsinB(*iter,b))&&(mon2int(*iter)>maxind))
00208       maxind=mon2int(*iter);
00209   maxind++;
00210   res.sizenf=(maxind/8)+1;
00211   res.nfind=(unsigned char*)malloc((maxind/8+1));
00212   memset(res.nfind,0,res.sizenf);
00213   tmp=(coeff*)MAC_REV_MALLOC<coeff>((maxind)*sizeof(coeff));
00214   for(int i=0;i<maxind;i++)
00215     tmp[i]=0;
00216   for(typename polyalp::const_iterator iter=P.begin();iter!=P.end();iter++)
00217     {
00218       if(IsinB(*iter,b))
00219         {
00220           if(!Iszero(iter->GetCoeff()))
00221             {
00222               int stockind=mon2int(*iter);
00223               tmp[stockind]=iter->GetCoeff();
00224               res.nfind[stockind/8]|=(unsigned char)(1<<(stockind%8));
00225               p++;
00226             }
00227         }
00228       else
00229         res.ind=*iter;
00230     }
00231   res.size=p;//l'ind n'est pas comptabilise dans la nf
00232   res.nf=(coeff*)MAC_REV_MALLOC<coeff>(res.size*sizeof(coeff));
00233    for( i=0;i<res.size;i++)
00234     res.nf[i]=0;
00235  
00236   j=0;
00237   for(int i=0;i<maxind;i++)
00238     if(!Iszero(coeff(tmp[i])))
00239       res.nf[j++]=tmp[i];
00240   MAC_REV_FREE<coeff>(tmp,(maxind)*sizeof(coeff));
00241   
00242   p=0;
00243   for(int i=0;i<res.sizenf;i++)
00244     p+=nbbits[res.nfind[i]];
00245 #ifdef DEB
00246   cout<<"pol.size "<<P.size()<<" kkk "<<p<<endl; 
00247 #endif
00248   return res;
00249 }
00250 
00251 template<typename mon,typename T>
00252 int Degree(const polter<mon,T> &p)
00253 {
00254   int mmax=-1;
00255   for(int i=0;i<p.sizenf;i++)
00256     if(p.nfind[i])
00257       {
00258         mon tmpmon;
00259         for(int j=0;j<8;j++)
00260           if((p.nfind[i]>>j)&1)
00261             {
00262               int2mon(8*i+j,tmpmon);
00263               mmax=(mmax<tmpmon.GetDegree())?tmpmon.GetDegree():mmax;
00264             }
00265       }
00266   if(!Iszero(p.ind.GetCoeff()))
00267      mmax=(mmax<p.ind.GetDegree())?p.ind.GetDegree():mmax;
00268   return mmax;//p.ind.GetDegree();
00269 }
00270 
00271 template<typename mon,typename T>
00272 int hasmonom(const polter<mon,T>& p,const mon & m)
00273 {
00274   int res=0,i,tmpind=mon2int(m);
00275   return (tmpind/8<p.sizenf) && ((p.nfind[tmpind/8]>>(tmpind%8))&1);
00276 }
00277