pol2bis.hpp

Go to the documentation of this file.

#ifndef polbis
#define polbis pol
#include"ugly.hpp"

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};

template<typename T> void *MAC_REV_MALLOC(int size);
template<typename T> void *MAC_REV_REALLOC(void*ptr,int,int);
template<typename T> void MAC_REV_FREE(void*ptr,int size);

template<typename mon,typename T>
struct polbis  
{
  typedef mon monom_t;
  typedef T coeff_t;
  mon ind;
  T * nf;
  unsigned char *nfind;
  int size,sizenf;
  //size correspond a nf et sizenf a nfind
  // pol operator*(const mon&t);//multi par un monome
  polbis operator*(const T &t);
  polbis &operator-=(const polbis &Q);
  polbis &operator=(const polbis &Q);
  polbis &operator/=(const T &t);
};

//
//l'autre par un coeff constant
//
template<typename mon,typename T>
polbis<mon,T> 
polbis<mon,T>::operator*(const typename polbis<mon,T>::coeff_t & toto)
{
  //int tmpsize=this->size;
  for(int i=0;i<this->size;i++)
    this->nf[i]*=toto;
  return *this;
}

template<typename mon,typename T>
polbis<mon,T>& polbis<mon,T>::operator-=(const polbis<mon,T> &Q)
{
  int sizeres=0;
  polbis<mon,T> res;
  //  cout<<"on m'appelle"<<endl;
  if(sizenf<Q.sizenf)
    {
      T *ptrnf=nf,*ptrresnf,*ptrQnf;
      for(int i=0;i<sizenf;i++)
    sizeres+=nbbits[nfind[i]&Q.nfind[i]];
 
      for(int i=sizenf;i<Q.sizenf;i++)
    sizeres+=nbbits[Q.nfind[i]];
      res.size=sizeres;
      //res.nfind=malloc(Q.sizenf);
      res.nf=(T*)MAC_REV_MALLOC<T>(res.size*sizeof(T));
      ptrresnf=res.nf;
      ptrQnf=Q.nf;
      for(int i=0;i<sizenf;i++)
    {
      for(int j=0;j<8;j++)
        {
          if((nfind[i]>>j)&1)
        *(ptrresnf++)=*(ptrnf++);
          if ((Q.nfind[i]>>j)&1)
        *(ptrresnf-1)-=*(ptrQnf++);
        }
    }
      for(int i=sizenf;i<Q.sizenf;i++)
    {
      for(int j=0;j<8;j++)
        if((Q.nfind[i]>>j)&1)
          *(ptrresnf++)=((T)-1)*(*(ptrQnf++));
    }
    }
  else
    {
      T *ptrnf=nf,*ptrresnf,*ptrQnf;
      for(int i=0;i<Q.nizenf;i++)
    sizeres+=nbbits[nfind[i]&Q.nfind[i]];
      
      for(int i=Q.sizenf;i<sizenf;i++)
    sizeres+=nbbits[nfind[i]];
      
      res.size=sizeres;
      // res.nfind=malloc(sizenf);
      res.nf=(T*)MAC_REV_MALLOC<T>(res.size*sizeof(T));
      ptrresnf=res.nf;
      ptrQnf=Q.nf;
      for(int i=0;i<Q.sizenf;i++)
    {
      for(int j=0;j<8;j++)
        {
          if((nfind[i]>>j)&1)
        *(ptrresnf++)=*(ptrnf++);
          if((Q.nfind[i]>>j)&1)
        *(ptrresnf-1)-=*(ptrQnf++);
        }
    }
      for(int i=Q.sizenf;i<sizenf;i++)
    {
      for(int j=0;j<8;j++)
        if((nfind[i]>>j)&1)
          *(ptrresnf++)=*(ptrnf++);
    }

    }
  if(sizenf<Q.sizenf)
    nfind=(unsigned char*)realloc(nfind,Q.sizenf);
  for(int i=0;i<Q.sizenf;i++)
    nfind[i]|=Q.nfind[i];
  MAC_REV_FREE<T>(sizenf*sizeof(T));
  free(nfind);
  nfind=res.nfind;
  nf=res.nf;
  sizenf=res.sizenf;
  size=res.size;
  return *this;
}
//attention l'operateur -= n'agit que sur les nf

//
//operateur d'affectation sans recopie de nf!!!!
//
template<typename mon,typename T>
polbis<mon,T> & polbis<mon,T>::operator= (const polbis<mon,T>&Q)
{
  nf=Q.nf;
  nfind=Q.nfind;
  size=Q.size;
  sizenf=Q.sizenf;
  ind=Q.ind;
  return *this;
}

template<typename mon,typename T>
polbis<mon,T> & polbis<mon,T>::operator/= (const T&t)
{
 for(int i=0;i<this->size;i++)
    this->nf[i]/=t;
  return *this;
}

template<typename polyalp, typename mon, typename T, typename monomial_server> 
polyalp invconv(const polbis<mon,T> & p, monomial_server &serv)
{
  typedef typename polyalp::order_t ord;
  polyalp res(0);
  T *tmpptr=p.nf;
  for(int i=0;i<p.sizenf;i++)
    {
      //cout<<"i dans invconv "<<i<<" sizeplace "<<sizeplacemon()<<endl;
      for(int j=0;j<8;j++)
    if ((p.nfind[i]>>j)&1)
      {
        mon tmp;
        serv.int2mon(8*i+j,tmp);
        //      cout<<"8*i+j vaut"<<8*i+j<<endl;
        //cout<<"j'aoute le monome "<<tmp<<endl;
        //cout<<"LE COEFF EST "<<p.nf[i]<<endl;  
        //cout<<"res dans invconv"<<endl<<res<<endl;
        res+=tmp*(*(tmpptr++));
      }
    }
  if (!Iszero(p.ind.GetCoeff()))
    res+=p.ind;
  //cout<<"res dans invconv"<<endl<<res<<endl;
  return res;
}


// renvoit mon/x_i .... mais suppose qu'on ne fait pas n'importe quoi
//necessaire pour convert

template <typename Mon,typename Base>
Mon divmon(const Mon &mon,int i,const Base &b)
{
  Mon res(1);
  for(int j=0;j<b.nvar();j++)
    if(i==j)
      {
    res*=Mon(j,mon.GetDegree(j)-1);
      }
  else 
    {
      res*=Mon(j,mon.GetDegree(j));
    }
  return res;
}

//converti un polyalp en polbis
template<typename polbis,typename polyalp,typename Base, 
     typename monomial_server>
polbis convert(const polyalp &P,const Base &b, monomial_server &serv)
{
  typedef typename polyalp::order_t ord;
  typedef typename polyalp::monom_t mon;
  typedef typename mon::coeff_t coeff;
  coeff *tmp;
  polbis res;
  int i,maxind=-1,j,p=0;
  i=0;
  res.ind=mon(0);
  for(typename polyalp::const_iterator iter=P.begin();iter!=P.end();iter++)
    if ((IsinB(*iter,b))&&(serv.mon2int(*iter)>maxind))
      maxind=serv.mon2int(*iter);
  maxind++;
  res.sizenf=(maxind/8)+1;
  res.nfind=(unsigned char*)malloc((maxind/8+1));
  memset(res.nfind,0,res.sizenf);
  tmp=(coeff*)MAC_REV_MALLOC<coeff>((maxind)*sizeof(coeff));
  for(int i=0;i<maxind;i++)
    tmp[i]=0;
  for(typename polyalp::const_iterator iter=P.begin();iter!=P.end();iter++)
    {
      if(IsinB(*iter,b))
    {
      if(!Iszero(iter->GetCoeff()))
        {
          int stockind=serv.mon2int(*iter);
          tmp[stockind]=iter->GetCoeff();
          res.nfind[stockind/8]|=(unsigned char)(1<<(stockind%8));
          p++;
        }
    }
      else
    res.ind=*iter;
    }
  res.size=p;//l'ind n'est pas comptabilise dans la nf
  res.nf=(coeff*)MAC_REV_MALLOC<coeff>(res.size*sizeof(coeff));
   for( i=0;i<res.size;i++)
    res.nf[i]=0;
 
  j=0;
  for(int i=0;i<maxind;i++)
    if(!Iszero(coeff(tmp[i])))
      res.nf[j++]=tmp[i];
  MAC_REV_FREE<coeff>(tmp,(maxind)*sizeof(coeff));
  
  p=0;
  for(int i=0;i<res.sizenf;i++)
    p+=nbbits[res.nfind[i]];
#ifdef DEB
  cout<<"pol.size "<<P.size()<<" kkk "<<p<<endl; 
#endif
  return res;
}

template<typename mon,typename T, typename monomial_server>
int Degree(const polbis<mon,T> &p,monomial_server &serv)
{
  int mmax=-1;
  for(int i=0;i<p.sizenf;i++)
    if(p.nfind[i])
      {
    mon tmpmon;
    for(int j=0;j<8;j++)
      if((p.nfind[i]>>j)&1)
        {
          serv.int2mon(8*i+j,tmpmon);
          mmax=(mmax<tmpmon.GetDegree())?tmpmon.GetDegree():mmax;
        }
      }
  if(!Iszero(p.ind.GetCoeff()))
     mmax=(mmax<p.ind.GetDegree())?p.ind.GetDegree():mmax;
  return mmax;//p.ind.GetDegree();
}

template<typename mon,typename T>
int hasmonom(const polbis<mon,T>& p,const mon & m)
{
  int tmpind=mon2int(m);
  return (tmpind/8<p.sizenf) && ((p.nfind[tmpind/8]>>(tmpind%8))&1);
}

template<typename mon, typename T>
mon monofindex(int j, pol<mon,T> &p)
{
  mon tmp;
  int cpt=0,i,tmpcpt;
  //for(i=0;(i<p.sizenf) ;i++)
  // cpt+=nbbits[p.nfind[i]];
  //cout<<endl<<"nbbits "<<cpt<<endl<<" et j "<<j<<endl;
  //cpt=0;
  for(i=0;(cpt<j) && (i<p.sizenf) ;i++)
    cpt+=nbbits[p.nfind[i]];
  tmpcpt=8*i;
  (i==p.sizenf?i--:i);
  for(int k=128;cpt>j;k>>=1)
    if(p.nfind[i]&&k)
      {tmpcpt--;cpt--;}
  int2mon(tmpcpt,tmp);
  return tmp;
}
#endif