IsinB3sdp.hpp

Go to the documentation of this file.

#ifndef ALREADY_IsinB
#define ALREADY_IsinB
//void reinit_multiple();
template<typename T>
int Iszero(const T&);

template<typename monom>
struct predicat 
{
  int taille1;
  monom * accept;
  int taille2;
  monom * refuse;
  typedef monom monom_t;
};

template<typename predicat>
struct Base
{
  list<predicat> def;
  int nbvars;
  typedef predicat predicat_t;
  typedef typename predicat::monom_t mon;
  typedef typename list<predicat>::iterator iterator;
  typedef typename list<predicat>::const_iterator const_iterator;
  typedef typename list<predicat>::reverse_iterator reverse_iterator;
  typedef typename list<predicat>::const_reverse_iterator const_reverse_iterator;

  iterator begin () { return def.begin(); }
  const_iterator begin () const { return def.begin(); }
  reverse_iterator rbegin ()    { return def.rbegin(); }
  const_reverse_iterator rbegin () const { return def.rbegin(); }

  iterator end () { return def.end(); }
  const_iterator end () const { return def.end(); }
  reverse_iterator rend () { return def.rend(); } 
  const_reverse_iterator rend () const { return def.rend(); }

  iterator erase(iterator first, iterator last){return def.erase(first,last);};
  iterator erase(iterator first){return def.erase(first);};
  int nbvar() {return nbvars;};
  int nbvar() const {return nbvars;};
};

template<typename mon,typename predicat,typename Base>
int isinb(const mon &m, predicat p,const Base &b)
{
  int inaccept=0,inrefuse=0;
   for(int i=0;i<p.taille2;i++)
    {
      inrefuse|=isdivisible(p.refuse[i],m,b);
      //      cout<<"p.refuse[i] "<<p.refuse[i]<<" "<<m<<" "<<inrefuse<<endl;
      if(inrefuse) {
          goto jesais; break;}
    }
   if(!inrefuse)
     for(int i=0;i<p.taille1;i++)
       {
     // cout<<"p.accept[i] "<<p.accept[i]<<" "<<m<<endl;
     inaccept|=isdivisible(p.accept[i],m,b);
     if(inaccept) break;
       }
 jesais:
   return((!inrefuse)&&inaccept);
}

template<typename mon,typename Base>
int IsinB(const mon &m,const Base &b)
{
  int res=0;
  //cout<<"monom "<<m<<endl;
  for(typename Base::const_iterator iter=b.begin();iter!=b.end();iter++)
    {
      res=res||isinb(m,*iter,b);
      //cout<<"verdict isinb "<<res<<endl;
      if (res) break;
    }

  return res;
}

template<typename mon,typename predicat,typename Base>
int isinb2(const mon &m, predicat p,const Base &b,mon &div,mon &coc)
{
  int inaccept=0,inrefuse=0;
  //mon div=mon(1);
   for(int i=0;i<p.taille2;i++)
    {
      inrefuse|=isdivisible2(p.refuse[i],m,b,coc);
      //      cout<<"p.refuse[i] "<<p.refuse[i]<<" "<<m<<" "<<inrefuse<<endl;
      if(inrefuse) {
          div=p.refuse[i];
          goto jesais; break;}
    }
   if(!inrefuse)
     for(int i=0;i<p.taille1;i++)
       {
         // cout<<"p.accept[i] "<<p.accept[i]<<" "<<m<<endl;
         inaccept|=isdivisible(p.accept[i],m,b);
         if(inaccept) break;
       }
 jesais:
   return((!inrefuse)&&inaccept);
}




template<typename mon,typename Base>
int IsinB2(const mon &m,const Base &b, mon &div,mon &coc)
{
  int res=0;
  //mon coc=mon(1);
  //cout<<"monom "<<m<<endl;
  for(typename Base::const_iterator iter=b.begin();iter!=b.end();iter++)
    {
      res=res||isinb2(m,*iter,b,div,coc);
      //cout<<"verdict isinb "<<res<<endl;
      if (res) break;
    }

  return res;
}

template<typename typMk,typename typdump,typename Base>
void Dovoisins(const typename typMk::value_type &m, typMk &voisins,
           const typdump & dump,const Base &b)
{
  typedef typename typMk::value_type mon;  
  int i=0;
  typename typdump::const_iterator iter=dump.begin();
  //on determine un parent
  typMk ancien;
  for(i=0;i<b.nbvar();i++)
    if(m.GetDegree(i)>0)
      ancien.push_back(my_divmon(m,i,b));
  //cout<<"Dans Dovoisins "<<m<<endl;
  for(;iter->k!=(m.GetDegree()-1);iter++);
  for(i=0;i<iter->size;i++)
    if(member(ancien,iter->nf[i].ind)) break;
  //maintenant on l'a on prend toutes les consequences de ce parent
  for(int j=0;j<b.nbvar();j++)
    if(!member(voisins,iter->nf[i].ind*mon(j,1)))
      voisins.push_back(iter->nf[i].ind*mon(j,1));
}

#ifdef uninf
  template<typename pol> void reinit_nf();
#endif

template<typename typPk,typename typMk ,typename typdump,typename Base, 
     typename monomial_server, typename workspace>
void rajoute_voisins(typPk &Pk,const typMk & Mk,
             const typdump &dump, Base &b,workspace &w,
             monomial_server &serv)
  // l'object if de la fonction est d'eviter les 
  //oublis de voisin lors d'un rajout au quotient
{
  typedef typename typPk::value_type pol;
  typedef typename pol::monom_t mon;
  typedef typename pol::coeff_t coeff;

#ifdef unimul
  // void reinit_multiple();
  serv.reinit_multiple();
#endif
#ifdef uninf

   serv.reinit_nf();
#endif 
#ifdef BONRAJOUT
  return;
#endif

  mon * baccept=b.rbegin()->accept;
  int taille=b.rbegin()->taille1;
  typMk voisins;
  for(int i=0;i<taille;i++)
    Dovoisins(baccept[i],voisins,dump,b);
  //Dovoisins ne cree pas de monomes redondants!!!!!!!!!!!!
  for(typename typMk::iterator iter=voisins.begin();iter!=voisins.end();iter++)
    {
      for(int i=0;i<taille;i++)
    if(baccept[i].rep==iter->rep)
      {
        voisins.erase(iter--);
        break;
    
      }  
    }
  for(typename typMk::iterator iter=voisins.begin();iter!=voisins.end();iter++)
    {
      if(member(Mk,*iter)) 
    voisins.erase(iter--);
    }
  //on vire les voisins que l on a rajoute a l'accept quotient
  //et ceux que l'on a deja
  for(typename typMk::iterator iter=voisins.begin();iter!=voisins.end();iter++)
    {
      typMk ancien;
      for(int i=0;i<b.nbvar();i++)
    if(iter->GetDegree(i)>0)
      ancien.push_back(my_divmon(*iter,i,b));
      //on vient de creer les parents possibles
      //maintenant on cherche un pere possible
      
      {
    typename typdump::const_iterator iterdump;
      for(iterdump=dump.begin();iterdump!=dump.end();iterdump++)
        if(iterdump->k==(iter->GetDegree()-1)) break;
    {
      int i,j;
      for(i=0;i<iterdump->size;i++)
        if(member(ancien,iterdump->nf[i].ind)) 
          break;
      //on a un parent on calcul la nf du monome *iter maintenant
      //d'abord on calcul la bonne variable
      for(j=0;j<b.nbvar();j++)
        if(iter->rep==(iterdump->nf[i].ind*mon(j,1)).rep)
          break;
      mon * tmp=NULL;
      int tmpsize=0;
     //Tue Dec  4 16:39:25 MET 2001
      //pol tmppol;
      //densify(iterdump->nf[i],tmppol);
      //cout<<"avant mult"<<endl;
      //Tue Dec  4 16:39:25 MET 2001
      //cout<<"Le parent"<<endl<<invconv<Poly>(iterdump->nf[i])<<endl;
      //Pk.push_back(mult(j,tmppol,dump,b,&tmp,tmpsize));
      Pk.push_back(mult(j,iterdump->nf[i],dump,b,&tmp,tmpsize,w,serv));
      //cout<<"j'ai pushe ici"<<endl<<invconv<Poly>(*(Pk.rbegin()))<<endl;
      //il n devrait y avoir personne dans tmp normalement
      //int toint=0;
      //for(;Pk.back().ind.GetCoeff()>toint;toint++);
      for(int k=0;k<tmpsize-1;k++)
        {
          for(typename typPk::iterator iterpk=Pk.begin()
            ;iterpk!=Pk.end();iterpk++)
        if(iterpk->ind.rep==tmp[k].rep)
          {
            //          my_plus(Pk.back(),*iterpk,tmp[k].GetCoeff(),&w,serv);//ICIOMP
            my_plus(Pk.back(),*iterpk,tmp[k].GetCoeff(),&w);//ICIOMP
            break;
          }
        }
      MAC_REV_FREE<mon>(tmp,sizeof(mon)*(tmpsize));
    }
      }
      //cout<<"coucou de rajoutevoisins "<<Pk.back().ind<<endl;
      //  cout<<"le pol rajoute vaut"<<endl<<invconv<Poly>(Pk.back())<<endl;
      Pk.back().ind.SetCoeff((coeff)1);
      //cout<<"Je l'ai pas !!! "<<endl;
      //exit(10);
    }
  //on rajoute a B les monomes correspondants au voisin
  //on enleve en fait...
  //
  // cout<<"rajoute voisins .size "<<voisins.size()<<endl;
  b.rbegin()->refuse=(mon*)MAC_REV_REALLOC<mon>(b.rbegin()->refuse
                       ,sizeof(mon)*b.rbegin()->taille2,
                sizeof(mon)*(b.rbegin()->taille2+voisins.size()));
  typename typMk::iterator iter=voisins.begin();
  for(unsigned int i=b.rbegin()->taille2;i<b.rbegin()->taille2+voisins.size()
    ;i++,iter++)
    b.rbegin()->refuse[i]=*iter;
  b.rbegin()->taille2+=voisins.size();
} 

//AddB n'est normalement appele que lorsque la matruice n'a pas rang plein 

template<typename typmat,typename typPk, typename typMk,typename typdump,
     typename typB,typename monomial_server>
void AddB(const typmat & mat, typPk & Pk,const typMk &Mk,const typdump & dump,
      typB &b, monomial_server &serv)
{
  typedef typename typPk::value_type pol;
  typedef typename pol::monom_t mon;
  typedef typename pol::coeff_t coeff_t;
  typename typPk::iterator iterpk;
  predicat<mon> tmp;
  typename typMk::const_iterator iter=Mk.begin();
  int r=my_rank(mat),comprefuse=0,compaccept=0;
  unsigned int i;
  cout<<"DANS ADDB"<<endl;
  cout<<"rang de la matrice m "<<r<<" Mk.size "<<Mk.size()<<endl;
  cout<<"taille de la matrice "<<mat.nrow<<"x"<<mat.ncol<<endl;
  tmp.refuse=(mon*)MAC_REV_MALLOC<mon>(r*sizeof(mon));
  tmp.taille2=r;
  tmp.accept=(mon*)MAC_REV_MALLOC<mon>((Mk.size()-r)*sizeof(mon));
  tmp.taille1=Mk.size()-r;
  //for(i=0;i<(unsigned)r;i++,iter++)
  //  tmp.refuse[i]=*iter;
  iter=Mk.begin();
  for(i=0;i<Mk.size();i++,iter++)
    {
      if(Iszero(mat(i,i)))
        {
          tmp.accept[compaccept++]=*iter;
          cout<<"monome rajoute au quotient "<<*iter<<endl;
        }
      else
        {
      tmp.refuse[comprefuse++]=*iter;
        }
    }
  b.def.push_back(tmp);
  //bon maintenant que l'on vient de rajouter des 
  //monomes au quotient il faut remttre les polynomes comme il faut
  iter=Mk.begin();
  for(unsigned int i=0;i<Mk.size();i++,iter++)
    {
      if(Iszero(mat(i,i)))
    {
      int index=serv.mon2int(*iter);
      iterpk=Pk.begin();
      //on commence a faire de l'algebre lineaire
      coeff_t *vec=(coeff_t*)MAC_REV_MALLOC<coeff_t>(i*sizeof(coeff_t));
      for(unsigned int j=0;j<i;j++)
        vec[j]=mat(j,i);
      //le vecteur est fait reste plus qu'a 
      //resoudre le syteme lineaire
      for(int k=i-1;k>=0;k--)
        {
          if(Iszero(mat(k,k))) continue;
          vec[k]/=mat(k,k);
          for(int j=k-1;j>=0;j--)
        vec[j]-=vec[k]*mat(j,k);
        }
      //bon ca a l'air d'etre bon pour la resolution du systeme lineaire
      //on remet les coeff de vec dans les pol a la bonne place
      for(unsigned int j=0;j<i;j++,iterpk++)
        {
          //cout<<"tour deboucle addB "<<j<<endl;
          if((!Iszero(vec[j])))//&&(!Iszero(mat(j,j))))
          {
        int i;
        iterpk->nf=(typename mon::coeff_t*)
          MAC_REV_REALLOC<typename mon::coeff_t>(iterpk->nf
                 ,iterpk->size*sizeof(typename mon::coeff_t)
                 ,(iterpk->size+1)*sizeof(typename mon::coeff_t));
        if((index/8)>=iterpk->sizenf)
          {
            iterpk->nfind=(unsigned char*)
              MAC_REV_REALLOC<unsigned char>(iterpk->nfind
                ,iterpk->sizenf,index/8+1);
            for(int i=iterpk->sizenf;i<index/8+1;i++)
              iterpk->nfind[i]=0;
            iterpk->sizenf=index/8+1;
            
          }
        
        i=iterpk->size;
        for(int k=iterpk->sizenf-1;k>index/8;k--)
          {
            //cout<<"k "<<k<<endl;
            if(iterpk->nfind[k]!=0)
              {
            for(int j=0;j<nbbits[iterpk->nfind[k]];j++)
              {
                iterpk->nf[i]=iterpk->nf[i-1];
                i--;
              }
            }
          }
        for(int k=7;k>(index%8);k--)
          if((iterpk->nfind[index/8]>>k)&1)
            {
              iterpk->nf[i]=iterpk->nf[i-1];
              i--;
            }
        iterpk->nf[i]=vec[j];
        iterpk->nfind[index/8]|=(1<<(index%8));
        iterpk->size++;
        //cout<<"rajoute le monome "<<*iter<<" a la ligne "<<j<<" "<<vec[j]<<" i "<<i<<endl;
        //cout<<"le pol"<<endl<<invconv<Poly>(*iterpk)<<endl;
          }
        }
      MAC_REV_FREE<coeff_t>(vec,i*sizeof(coeff_t));
    }
    }
  //rajoute_voisins(Pk,Mk,dump,b);
  return;
}
      

//AddB n'est normalement appele que lorsque la matrice n'a pas rang plein 
//et que le second membre se reduit a 0
#if 1
template<typename mon,typename typB>
predicat<mon> findpred(const mon &m, const typB &b)
//findpred est sense renvoyer le predicat qui fait que le
// monome "m" se trouve sur la frontiere de b
// il faut que m soit sur la frontiere du predicat
{
  int ok=0;
  typename typB::const_iterator iter;
  for(iter=b.begin();iter!=b.end();iter++)
    {
      for(int i=0;i<iter->taille1;i++)
    if(isdivisible(iter->accept[i],m,b))
      ok=1;
      if(ok)
    {
      for(int i=0;i<iter->taille2;i++)
        if(isdivisible(iter->refuse[i],m,b))
          {
        for(int j=0;j<b.nbvar();j++)
          ok+=(m.GetDegree(i)>m.GetDegree(i));
        if(ok==1+b.nbvar())
          ok=0;
          }
    }
      if(ok)
    break;
    }
  cout<<"je cherche le pred pere de "<<m<<endl;
  cout<<"je renvoie ce pred"<<endl;
  for(int i=0;i<iter->taille1;i++)
    {
      cout<<"iter.accept[i] "<<iter->accept[i]<<endl;
    }
  for(int i=0;i<iter->taille2;i++)
    {
      cout<<"iter.refuse[i] "<<iter->refuse[i]<<endl;
    }

  return *iter;
}

template<typename typmat,typename typPk, typename typMk,typename typdump,
     typename typB, typename typP>
void AddB_dim(const typmat & mat, const typmat &L, typP &P, typPk & Pk
          ,const typMk &Mk,const typdump & dump, typB &b)
{
  
  typedef typename typPk::value_type pol;
  typedef typename pol::monom_t mon;
  typedef typename pol::coeff_t coeff_t;
  typename typPk::iterator iterpk;
  
  typename typMk::const_iterator iter=Mk.begin();
   int r=my_rank(mat),comprefuse=0,compaccept=0;
  unsigned int i;
  //cout<<"DANS ADDB"<<endl;
  //cout<<"rang de la matrice m "<<r<<" Mk.size "<<Mk.size()<<endl;
  //cout<<"taille de la matrice "<<mat.nrow<<"x"<<mat.ncol<<endl;
  //  tmp.refuse=(mon*)MAC_REV_MALLOC<mon>(r*sizeof(mon));
  //tmp.taille2=r;
  //tmp.accept=(mon*)MAC_REV_MALLOC<mon>((Mk.size()-r)*sizeof(mon));
  //tmp.taille1=Mk.size()-r;
  //for(i=0;i<(unsigned)r;i++,iter++)
  //  tmp.refuse[i]=*iter;
 
 //bon maintenant que l'on vient de rajouter des 
  //monomes au quotient il faut remttre les polynomes comme il faut
  iter=Mk.begin();  
  for(unsigned int i=0;i<Mk.size();i++,iter++)
    {
      if(Iszero(mat(i,i)))
    {
      int index=mon2int(*iter);
      iterpk=Pk.begin();
      //on commence a faire de l'algebre lineaire
      coeff_t *vec=(coeff_t*)MAC_REV_MALLOC<coeff_t>(i*sizeof(coeff_t));
      for(unsigned int j=0;j<i;j++)
        vec[j]=mat(j,i);
      //le vecteur est fait reste plus qu'a 
      //resoudre le syteme lineaire
      for(int k=i-1;k>=0;k--)
        {
          if(Iszero(mat(k,k))) continue;
          vec[k]/=mat(k,k);
          for(int j=k-1;j>=0;j--)
        vec[j]-=vec[k]*mat(j,k);
        }
      //bon ca a l'air d'etre bon pour la resolution du systeme lineaire
      //on remet les coeff de vec dans les pol a la bonne place
      for(unsigned int j=0;j<i;j++,iterpk++)
        {
          //cout<<"tour deboucle addB "<<j<<endl;
          if((!Iszero(vec[j])))//&&(!Iszero(mat(j,j))))
          {
        int i;
        iterpk->nf=(typename mon::coeff_t*)
          MAC_REV_REALLOC<typename mon::coeff_t>(iterpk->nf
                 ,iterpk->size*sizeof(typename mon::coeff_t)
                 ,(iterpk->size+1)*sizeof(typename mon::coeff_t));
        if((index/8)>=iterpk->sizenf)
          {
            iterpk->nfind=(unsigned char*)
              MAC_REV_REALLOC<unsigned char>(iterpk->nfind
                ,iterpk->sizenf,index/8+1);
            for(int i=iterpk->sizenf;i<index/8+1;i++)
              iterpk->nfind[i]=0;
            iterpk->sizenf=index/8+1;
            
          }
        
        i=iterpk->size;
        for(int k=iterpk->sizenf-1;k>index/8;k--)
          {
            //cout<<"k "<<k<<endl;
            if(iterpk->nfind[k]!=0)
              {
            for(int j=0;j<nbbits[iterpk->nfind[k]];j++)
              {
                iterpk->nf[i]=iterpk->nf[i-1];
                i--;
              }
            }
          }
        for(int k=7;k>(index%8);k--)
          if((iterpk->nfind[index/8]>>k)&1)
            {
              iterpk->nf[i]=iterpk->nf[i-1];
              i--;
            }
        iterpk->nf[i]=vec[j];
        iterpk->nfind[index/8]|=(1<<(index%8));
        iterpk->size++;
        //cout<<"rajoute le monome "<<*iter<<" a la ligne "<<j<<" "<<vec[j]<<" i "<<i<<endl;
        //cout<<"le pol"<<endl<<invconv<Poly>(*iterpk)<<endl;
          }
        }
      MAC_REV_FREE<coeff_t>(vec,i*sizeof(coeff_t));
    }
    }
  //rajoute_voisins(Pk,Mk,dump,b);
  iter=Mk.begin();
  for(i=0;i<Mk.size();i++,iter++)
    {
      if(Iszero(mat(i,i))) 
    {
      list<mon> listlcm;
      mon pere;
      predicat<mon> tmp,stockpred=findpred(*iter,b);
      int j=0,k=0;
      tmp.refuse=(mon*)MAC_REV_MALLOC<mon>
        ((stockpred.taille2)*sizeof(mon));
      tmp.taille2=(stockpred.taille2-1);
      tmp.accept=(mon*)MAC_REV_MALLOC<mon>(sizeof(mon));
      tmp.taille1=1;
      tmp.accept[0]=*iter;
      for(int j=0;j<stockpred.taille2;j++)
        if(isdivisible(stockpred.refuse[j],*iter,b))
          {
        pere.rep=stockpred.refuse[j].rep;
        break;
          }
      cout<<"le pere "<<pere<<endl;
      for(int j=0;j<stockpred.taille2;j++)
        if(stockpred.refuse[j].rep!=pere.rep)
          listlcm.push_back(lcm(pere,stockpred.refuse[j]));
      for(typename list<mon>::iterator iterlcm=listlcm.begin();
          iterlcm!=listlcm.end();iterlcm++)
        for(typename list<mon>::iterator iterlcm2=listlcm.begin();
        iterlcm2!=listlcm.end();iterlcm2++)
          if ((iterlcm2!=iterlcm) && (isdivisible(*iterlcm,*iterlcm2,b)))
        {
          cout<<"je vire "<<*iterlcm2<<" parce que divisible par "<<*iterlcm<<endl;
          listlcm.erase(iterlcm2--);
        }
      cout<<"apres tri il me reste"<<endl;
      for(typename list<mon>::iterator iterlcm=listlcm.begin();
          iterlcm!=listlcm.end();iterlcm++)
        cout<<*iterlcm<<endl;
      //On a maintenant les minimaux pour la division
      //reste a retirer les bons
      for(typename list<mon>::iterator iterlcm=listlcm.begin();
          iterlcm!=listlcm.end();iterlcm++)
        if(iterlcm->GetDegree()>=Mk.begin()->GetDegree())
          {
        for(j=0;j<stockpred.taille2;j++)
          if ((isdivisible(stockpred.refuse[j],*iterlcm,b))
              && (stockpred.refuse[j].rep != pere.rep))
            break;      
        tmp.refuse[k++]=stockpred.refuse[j];
          }
#if 0
      iterpk=Pk.begin();
      cout<<"mon rajoute "<<*iter<<endl;
      for(typename typMk::const_iterator iter2=Mk.begin()
        ;iter2!=iter;iter2++,j++,iterpk++)
        if(!Iszero(L(i,j)))
          {
        predicat<mon> tmppred=findpred(*iter2,b);
        mon tmpmon=lcm(*iter,*iter2), tmpmon2=mon(1);
        for(int l=0;l<b.nbvar();l++)
          tmpmon2*=mon(l,tmpmon.GetDegree(l)-iter2->GetDegree(l));
        tmp.refuse[k++]=tmpmon;
        cout<<"iterpk->ind*tmpmon2 "<<iterpk->ind*tmpmon2<<endl;
        cout<<invconv<Poly>(*iterpk)*tmpmon2<<endl;;
        P.push_back(invconv<Poly>(*iterpk)*tmpmon2);
          }
#endif  
      tmp.refuse=(mon*)MAC_REV_REALLOC<mon>(tmp.refuse,
                        (stockpred.taille2)*sizeof(mon)
                        ,k*sizeof(mon));
      tmp.taille2=k;
      b.def.push_back(tmp);  
      //cout<<"monome rajoute au quotient "<<*iter<<endl;
    }
    }
  cout<<"Dans Addb"<<endl;
      for(typename typB::iterator p=b.begin();p!=b.end();p++)
        {
    for(int i=0;i<p->taille1;i++)
      {
        cout<<"p.accept[i] "<<p->accept[i]<<endl;
      }
    for(int i=0;i<p->taille2;i++)
      {
        cout<<"p.refuse[i] "<<p->refuse[i]<<endl;
      }
        }
      cout<<"eeeeeeeeeeeeeeeeeeeeeeeeeee"<<endl;
return;
}

#endif

#endif //ALREADY_IsinB