corealgo3_omp.hpp

Go to the documentation of this file.

#ifndef ALREADY_corealgo3
#define ALREADY_corealgo3
#include<omp.h>

template<typename typP,typename typPk,typename Base>
void Crochetchoix(const typP & P,typPk & Pk, typPk & redspol,int k,
            Base &b);

template<typename typPk>
void stab(typPk & Pk)
{
  typedef typename typPk::value_type pol;
  for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
    {
      pol tmpres=mon(0);
      for(typename pol::iterator iterpol=iter->begin();iterpol!=iter->end();
      iterpol++)
    {
      if(!Iszero(iterpol->GetCoeff()))
        tmpres+=*iterpol;
    }
      *iter=tmpres;
    }
}

template<typename T>
void  my_merge(T &l1,const T &l2)
{
  typename T::const_iterator iter=l2.begin();
  for(;iter!=l2.end();iter++)
     l1.push_back(*iter);
}
template<typename mon>
mon lcm(const mon &m1,const mon &m2)
{
  mon m(1);
  int maxi=max(lvar(m1.rep),lvar(m2.rep));
  for(int i=0;i<=maxi;i++)
    m*=mon(i,max(m1.GetDegree(i),m2.GetDegree(i)));
  return m;
}

static int max_degppcm;

template<typename typppcm,typename Base>
int compute_ppcm(typppcm &ppcm,const Base & b)
{
  typedef typename Base::mon mon;
  //typppcm tmpppcm;
  max_degppcm=0;
  int maxdeg=-1;
  list<mon> list1,list2,tmplist2;
 

 for(typename Base::const_reverse_iterator iter=b.rbegin()
    ;iter==b.rbegin();iter++)
    {
#ifndef oldppcm
      for(int i=0;i<iter->taille2;i++)
    for(int j=i+1;j<iter->taille2;j++)
      {
        mon tmpmon=lcm(iter->refuse[i],iter->refuse[j]);
        
        if((!ppcm.count(tmpmon))&& 
           (tmpmon.GetDegree()!=
        (iter->refuse[i].GetDegree()+iter->refuse[j].GetDegree())))
          {
        ppcm[tmpmon]=1;
        if(tmpmon.GetDegree()>max_degppcm)
          max_degppcm=tmpmon.GetDegree();
          }
        
        //      cout<<"ppcm rajoute "<<tmpmon<<" "<<iter->refuse[i]<<" "<<iter->refuse[j]<<endl;
      }
#else      
      for(int i=0;i<iter->taille2;i++)
    for(int j=i+1;j<iter->taille2;j++)
      {
        mon tmpmon=lcm(iter->refuse[i],iter->refuse[j]);
    
             if((tmpmon.GetDegree()!=
        (iter->refuse[i].GetDegree()+iter->refuse[j].GetDegree())))
          {
        int ok=1;
        for(int k=0;k<iter->taille2;k++)
          {
            if((k!=i) && (k!=j) 
               && (isdivisible(iter->refuse[k],tmpmon,b)))
              {
            mon lc_i=lcm(iter->refuse[k],iter->refuse[i]);
            mon lc_j=lcm(iter->refuse[k],iter->refuse[j]);
            if((lc_i.rep!=tmpmon.rep) 
               && (isdivisible(lc_i,tmpmon,b)) &&
               (lc_j.rep!=tmpmon.rep)
               &&(isdivisible(lc_j,tmpmon,b)))
              ok=0;
              }
          }
        if(ok)
          list2.push_back(tmpmon);
          }
      }
#endif
    }
#ifdef oldppcm  
  maxdeg=0;
  //  cout<<"liste2.size "<<list2.size()<<endl;
  for(typename list<mon>::iterator iter2=list2.begin();
      iter2!=list2.end();iter2++)
    {
      // cout<<"je passe ici "<<*iter2<<endl;
      ppcm[*iter2]=1;
      if(iter2->GetDegree()>maxdeg)
    maxdeg=iter2->GetDegree();
    }
  max_degppcm=maxdeg;
#endif
  return maxdeg;
}

int GB_stop(int k)
{
  int res=0;
  if(k>max_degppcm) res=1;
  return res;
}

template<typename mon,typename Base,typename typppcm,typename typMk>
int goodspol(const mon &m1,const mon &m2, typppcm & ppcm,const Base &b
         ,const typMk &Spoldejavu)
{
  //list<mon> ppcm;
  int res=0,comp=0;//,nobplus=1;
  mon m=m1;
  for(int i=0;i<b.nbvar();i++)
    if(m1.GetDegree(i)!=m2.GetDegree(i))
      {
    comp++;
    if (m1.GetDegree(i)-m2.GetDegree(i)==1)
      res+=1;
    if (m1.GetDegree(i)-m2.GetDegree(i)==-1)
      {
        res+=1;
        m*=mon(i,1);
      }
      };
  //compute_ppcm(ppcm,b);
  //if((res==2)&&(comp==2)) 
  //  for(int i=0;i<b.nbvar();i++)
  //    if(m.GetDegree(i)>0)
  //    nobplus=nobplus && !IsinB(my_divmon(m,i,b),b);
   //cout<<"m1 "<<m1<<" m2 "<<m2<<" resultat "<<((res==2)&&(comp==2))<<endl;
   //return ((res==2)&&(comp==2)&&nobplus);
  //cout<<m1<<" < "<<m2<<" res "<<(member(ppcm,m)&&(res==2)&&(comp==2))<<endl;
  
  //  cout<<"m "<<m<<endl;
  //il y a un pbm a ne pas refaire les spol deja vu je ne comprends 
  //pas pourquoi voir example caprass
  #ifdef NO_TWICE
  return (res==2)&&(comp==2)&&ppcm[m]&&!member(Spoldejavu,m);
  #else
  return (res==2)&&(comp==2)&&ppcm.count(m);
  #endif
}

template<typename typPk,typename Base,typename typMk>
typPk Spolordsup(const typPk &Pk,const Base & b, const typMk & Spoldejavu)
{
  typedef typename typPk::value_type::monom_t mon;
  //typedef typename typ::order_t O;
  map<mon,int> ppcm;
  typPk res;
  int maxdeg;
  //cout<<"compute ppcm"<<endl;
  maxdeg=compute_ppcm(ppcm,b);
  //cout<<"test goodspol"<<endl;
  for(typename typPk::const_iterator iter1=Pk.begin();iter1!=Pk.end();iter1++)
    for(typename typPk::const_iterator iter2=iter1;iter2!=Pk.end();iter2++)
      if (goodspol(iter1->ind,iter2->ind,ppcm,b,Spoldejavu))
    {
      res.push_back(*iter1);
      res.push_back(*iter2);
    }
  return res;
}
template<typename typP>
int minmaxdegP(typP P,int &maxdeg)
{
  int res=Degree(P.front());
  maxdeg=res;
  for(typename typP::iterator iter=P.begin();iter!=P.end();iter++)
    {
      res=(res>Degree(*iter))?Degree(*iter):res;
      maxdeg=(maxdeg<Degree(*iter))?Degree(*iter):maxdeg;
    }
  return res;
}
template<typename typP>
typP Crochet(const typP &P, int k)
{
  typP res;
  for(typename typP::const_iterator iter=P.begin();iter!=P.end();iter++)
    {
      if(Degree(*iter)==k) res.push_back(*iter);
    }
  return res;
}

template<typename mon, typename T, typename monomial_server>
std::list<pol<mon,T> > Crochet(const std::list<pol<mon,T> > &P, int k, 
                   monomial_server &serv)
{
  std::list<pol<mon,T> > res;
  for(typename std::list<pol<mon,T> >::const_iterator iter=P.begin()
    ;iter!=P.end();iter++)
    {
      if(Degree(*iter,serv)==k) res.push_back(*iter);
    }
  return res;
}

template<typename typP, typename Base>
void Setnbvars(const typP &tmpPk,Base& b)
{
  int tmpvar=0;
  for(typename typP::const_iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
    for(typename typP::value_type::const_iterator iterpol=iter->begin();
    iterpol!=iter->end();iterpol++)
      if(tmpvar<lvar(iterpol->rep)) tmpvar=lvar(iterpol->rep);
  b.nbvars=tmpvar+1;
}

template<typename typP, typename typMk, typename Base, 
     typename monomial_server>
void DoBaseMk(const typP &tmpPk, typMk &Mk, Base &b, monomial_server &serv)
{
  typedef typename typMk::value_type mon;
  typename Base::predicat_t tmp;
  int i=0;
  //int tmpvar=0;
   
  for(typename typP::const_iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
    {
      mon tmp(choice(*iter,b,serv));
      int yes=1;
      //  cout<<"pol dans DoBaseMk "<<endl<<*iter<<endl<<"g choisi "<<tmp<<endl;
      tmp.SetCoeff(1);//*=mon(1/tmp.GetCoeff());
      for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
    if(*iter==tmp) {yes=0;break;};
      if(yes) Mk.push_back(tmp);
    }
  //cout<<"Mk.size()"<<Mk.size()<<endl;
  tmp.taille2=Mk.size();
  tmp.refuse=(mon*)MAC_REV_MALLOC<mon>(tmp.taille2*sizeof(mon));
  tmp.taille1=1;
  tmp.accept=(mon*)MAC_REV_MALLOC<mon>(sizeof(mon));
  // cout<<"mon(1) "<<mon(1)<<endl;
  tmp.accept[0]=mon(1);
  //cout<<"tmp.accept "<<tmp.accept[0]<<endl;
  i=0;//pour en etre bien sur
  for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++,i++)
    {
      tmp.refuse[i]=*iter;
    };
  if(!b.def.empty())
    {
      typename Base::predicat_t verytmp;
      for(typename Base::iterator iter=b.begin();iter!=b.end();iter++)
    {
      verytmp=*iter;
      MAC_REV_FREE<mon>(verytmp.refuse,verytmp.taille2*sizeof(mon));
      MAC_REV_FREE<mon>(verytmp.accept,verytmp.taille1*sizeof(mon));
    }
      b.erase(b.begin(),b.end());
    }
  b.def.push_back(tmp);
  return;
}

template <typename typPk,typename typMk,typename typmat >
typmat*
MatrixOf(const typPk & L, 
     const typMk & lm,const typmat &mat) 
{
  typedef typename typPk::value_type POL;
  typedef typename typMk::value_type monom_t;  
  typmat *res=NULL;
  map<monom_t,int> index;
  int l=1;
  if (res!=NULL) {delete res;};
  for(typename typMk::const_iterator m = lm.begin(); m!=lm.end(); m++){
    index[*m]=l;
    l++;
  }
  int c;
  res=new typmat(L.size(),lm.size());
  l=0;
  // cout<<"L.size() "<<L.size()<<" lm.size() "<<lm.size()<<endl;
  for(typename typPk::const_iterator pi = L.begin(); pi != L.end(); ++pi){
    for(typename POL::const_iterator m = pi->begin(); m!=pi->end(); m++)
      {
    
    if((index.count(*m))) 
      {
        c=index[*m];
        (*res)(l,c-1,m->GetCoeff());
      }
      }
    l++;
  }
  return res;   
}

template<typename typMat,typename typP,typename typPk,typename typMk, typename Base,typename typdump, class monomial_server, typename workspace>
void init_algo(typP &P,int &nouveaumon, int &maxdeg, int & k,
           typPk &Pk, typMk & Mk, Base & b, typdump & dump, workspace &w,
           monomial_server &serv)
{
  typedef typename typPk::value_type pol;
  typedef typename pol::coeff_t coeff;
  typedef typename pol::monom_t mon;
  const typename typP::value_type::order_t (tmpord)(const mon,const mon);
  typP tmpPk=P;
  pol debugpol;
  typMat *Mat;
  int *pr,*pc,flag=1;
  serv.setorder(my_ord);
  nouveaumon=1;
  Setnbvars(P,b);
  init_space(&w);
#ifdef unimul
  serv.init_multiple(b.nbvars);
#endif
#ifdef uninf
  serv.init_nf(b.nbvars);
#endif 
  flag=1;
while(flag)
    {
      Mk.erase(Mk.begin(),Mk.end());
      k=minmaxdegP(tmpPk,maxdeg);
      tmpPk=Crochet(tmpPk,k);
      DoBaseMk(tmpPk,Mk,b,serv);
      Mat=MatrixOf(tmpPk,Mk,typMat());
      tmpPk=Solve(*Mat,tmpPk,&pr,&pc);
      selectnozero(tmpPk);
      stab(tmpPk);
      if(rank(*Mat)==(int)tmpPk.size()) flag=0;
      Mat->destroystore();
      delete(Mat);
      free(pr);
      free(pc);
    }
      for(typename typP::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
    Pk.push_back(convert<pol>(*iter,b,serv));
      Dump(Pk,dump);
  return ;
}

template<typename pol,typename coeff, typename workspace,typename typserv>
void inline my_plus(pol &res,const pol & nf,const coeff toto,workspace *w,
            typserv &serv)
{
  typedef typename pol::monom_t mon;
  typedef typename mon::coeff_t mcoeff_t;
  //symbolic size
  int i=0,j=0,k=0,p=0,q=0;
  int slicep[omp_get_max_threads()+1];
  int sliceq[omp_get_max_threads()+1];
  int tmpsize=0,tmp;
  mcoeff_t *tmpnf,*tmpspace,zero=0;
  unsigned char *tmpnfind;
  cout<<"toto "<<toto<<endl;
  cout<<"res "<<invconv<Poly>(res,serv)<<endl;
  cout<<"nf "<<invconv<Poly>(nf,serv)<<endl;
  if(res.sizenf<=nf.sizenf)
    {
      //      cout<<"res.sizenf<nf.sizenf "<<endl;
      if(w->sizenf<nf.sizenf)
    increase_space(w);
      for(i=0;i<omp_get_max_threads()+1;i++)
    {
      slicep[i]=0;
      sliceq[i]=0;
    }
#pragma omp parallel for shared(w,toto,slicep,sliceq)
      for(i=0;i<res.sizenf;i++)
    {
      slicep[omp_get_thread_num()+1]+=nbbits[res.nfind[i]];
      sliceq[omp_get_thread_num()+1]+=nbbits[nf.nfind[i]];
    }
      for(i=1;i<omp_get_thread_num()+1;i++)
    {
      slicep[i]+=slicep[i-1];
      sliceq[i]+=sliceq[i-1];
    }
#pragma omp parallel for shared(w,toto,slicep,sliceq) firstprivate(p,q)
      for(i=0;i<res.sizenf;i++)
      {
        if(res.nfind[i]|nf.nfind[i])
          {
        w->tabnf[i]=res.nfind[i]|nf.nfind[i];
        for(int j=0;j<8;j++)
          {
            if ((res.nfind[i]>>j)&1)
              w->tabcoeff[8*i+j]=
            res.nf[p++ +slicep[omp_get_thread_num()]];
            if((nf.nfind[i]>>j)&1)
              w->tabcoeff[8*i+j]-=toto*nf.nf[q++ 
                     +sliceq[omp_get_thread_num()]];
            if(Iszero(w->tabcoeff[8*i+j]))
              w->tabnf[i]&=~(1<<j);
          }
          }
      }
      //      cout<<"slsliceq[omp_get_max_threads()] "<<sliceq[omp_get_max_threads()]<<endl;
      //int ttmp=0;
      //for(int j=0;j<res.sizenf;j++)
      //    ttmp+=nbbits[nf.nfind[j]];
      //cout<<" ttmp "<<ttmp<<endl;
      q=sliceq[omp_get_max_threads()];
//#pragma omp parallel for shared(w,toto)
      for(i=res.sizenf;i<nf.sizenf;i++)
      {
        if(nf.nfind[i])
          for(int j=0;j<8;j++)
        if ((nf.nfind[i]>>j)&1)
          if (!Iszero(nf.nf[q++]))
            w->tabcoeff[8*i+j]=((mcoeff_t)-1)*toto*(nf.nf[q-1]);
          else
            w->tabnf[i]&=~(1<<j);
        w->tabnf[i]=nf.nfind[i];
      }
      for(int i=0;i<nf.sizenf;i++)
    tmpsize+=nbbits[w->tabnf[i]];
      res.nf=(COEFF*)MAC_REV_REALLOC<COEFF>(res.nf,res.size*sizeof(COEFF)
                        ,tmpsize*sizeof(COEFF));
      res.size=tmpsize;
      res.nfind=(unsigned char*)realloc(res.nfind,nf.sizenf);
      res.sizenf=nf.sizenf;

      for(int i=0;i<nf.sizenf;i++)
    {
      if(w->tabnf[i])
        for(int j=0;j<8;j++)
          if ((w->tabnf[i]>>j)&1)
        //if (!Iszero(w->tabcoeff[8*i+j]))
        {
          res.nf[k++]=w->tabcoeff[8*i+j];
          w->tabcoeff[8*i+j]=zero;
        }
      //else
      //  w->tabnf[i]&=~(1<<j);
      res.nfind[i]=w->tabnf[i];
      w->tabnf[i]=0;
    }
      cout<<"le resultat "<<invconv<Poly>(res,serv)<<endl;
    }
  else
    {
      if(w->sizenf<res.sizenf)
    increase_space(w);
      for(i=0;i<omp_get_max_threads()+1;i++)
    {
      slicep[i]=0;
      sliceq[i]=0;
    }
      //#pragma omp parallel for shared(w,toto,slicep,sliceq)
      for(i=0;i<nf.sizenf;i++)
    {
      slicep[omp_get_thread_num()+1]+=nbbits[res.nfind[i]];
      sliceq[omp_get_thread_num()+1]+=nbbits[nf.nfind[i]];
    }
      for(i=1;i<omp_get_thread_num()+1;i++)
    {
      slicep[i]+=slicep[i-1];
      sliceq[i]+=sliceq[i-1];
    }
      //#pragma omp parallel for shared(w,toto,slicep,sliceq) firstprivate(p,q)
      for(i=0;i<nf.sizenf;i++)
#if 1
    {
        if(res.nfind[i])
          {
        unsigned char mask=1;
        w->tabnf[i]=res.nfind[i];
        if(res.nfind[i]==255)
          for(int j=0;j<8;j++)
            w->tabcoeff[8*i+j]=res.nf[p++];
        else
          for(int j=0;j<8;j++,mask<<=1)
            {
              //if ((res.nfind[i]>>j)&1) 
              if (res.nfind[i]&mask)
            w->tabcoeff[8*i+j]=res.nf[p++];
            }
          }
        if(nf.nfind[i])
        {
          unsigned char mask=1;
          w->tabnf[i]|=nf.nfind[i];
          if(nf.nfind[i]==255)
            for(int j=0;j<8;j++,mask<<=1)
              {
            w->tabcoeff[8*i+j]-=toto*nf.nf[q++];
            if(Iszero(w->tabcoeff[8*i+j]))
              w->tabnf[i]&=~mask;
              }
          else
            for(int j=0;j<8;j++,mask<<=1)
              {
            if(nf.nfind[i]&mask)
              w->tabcoeff[8*i+j]-=toto*nf.nf[q++];
            if(Iszero(w->tabcoeff[8*i+j]))
              w->tabnf[i]&=~mask;
              }
        }
      }
#endif
#if 0
        if(res.nfind[i]|nf.nfind[i])
          {
        w->tabnf[i]=res.nfind[i]|nf.nfind[i];
        for(int j=0;j<8;j++)
          {
            if ((res.nfind[i]>>j)&1)
              {
            w->tabcoeff[8*i+j]=res.nf[p++];
              }
            if((nf.nfind[i]>>j)&1)
              w->tabcoeff[8*i+j]-=toto*nf.nf[q++];
            if(Iszero(w->tabcoeff[8*i+j]))
              w->tabnf[i]&=~(1<<j);
          }
          }
#endif
        //    p=slicep[omp_get_max_threads()];
//#pragma omp parallel for shared(w,toto)
      for(i=nf.sizenf;i<res.sizenf;i++)
      {
        if(res.nfind[i])
          for(int j=0;j<8;j++)
        if ((res.nfind[i]>>j)&1) 
          if (!Iszero(res.nf[p++]))
            w->tabcoeff[8*i+j]=res.nf[p-1];
          else
            w->tabnf[i]&=~(1<<j);
        w->tabnf[i]=res.nfind[i];
      }
      for(int i=0;i<res.sizenf;i++)
    tmpsize+=nbbits[w->tabnf[i]];
      res.nf=(COEFF*)MAC_REV_REALLOC<COEFF>(res.nf,res.size*sizeof(COEFF)
                        ,tmpsize*sizeof(COEFF));
      res.size=tmpsize;
      //  res.nfind=(unsigned char*)realloc(res.nfind,nf.sizenf);
      //res.sizenf=nf.sizenf;
      for(int i=0;i<res.sizenf;i++)
    {
      if(w->tabnf[i])
        for(int j=0;j<8;j++)
          if ((w->tabnf[i]>>j)&1)
        //if (!Iszero(w->tabcoeff[8*i+j]))
        {
            res.nf[k++]=w->tabcoeff[8*i+j];
            w->tabcoeff[8*i+j]=zero;
        }
          //else
          //  w->tabnf[i]&=~(1<<j);
    
      res.nfind[i]=w->tabnf[i];
      w->tabnf[i]=0;
    }
    }

#ifdef DEB
  cout<<"k "<<k<<" tmpsize "<<tmpsize<<endl;
  int cpt=0;
  for(int i=0;i<res.sizenf;i++)
    if(res.nfind[i])
      for(int j=0;j<8;j++)
    if ((res.nfind[i]>>j)&1) 
      if (Iszero(res.nf[cpt++]))
        cout<<"j'ai fait n'importe quoi ici"<<endl;
  cout<<"les resultat "<<endl<<invconv<Poly>(res)<<endl;
#endif
  return;
}
template<typename pol,typename coeff>
void inline my_plus(pol &res,const pol & nf,const coeff toto,int removeme)
{
  typedef typename pol::monom_t mon;
  typedef typename mon::coeff_t mcoeff_t;
  //symbolic size
  int i=0,j=0,k=0,p=0,q=0;
  int tmpsize=0,tmp;
  mcoeff_t *tmpnf,*tmpspace;
  unsigned char *tmpnfind;
#ifdef DEB
  cout<<"dans myplus je calcul"<<endl<<invconv<Poly>(res,serv)<<endl<<
    "toto "<<toto<<endl<<"nf"<<endl<<invconv<Poly>(nf,serv)<<endl;
#endif
  if(res.sizenf<nf.sizenf)
    {
      tmpspace=(mcoeff_t*)
    MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*nf.sizenf*8);
      for(int i=0;i<nf.sizenf*8;i++)
    tmpspace[i]=0;
      for(int i=0;i<res.sizenf;i++)
    tmpsize+=nbbits[res.nfind[i]|nf.nfind[i]];
      for(int i=0;i<nf.sizenf;i++)
    tmpsize+=nbbits[nf.nfind[i]];
      tmpnf=(mcoeff_t*)
    MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*tmpsize);
      tmpnfind=(unsigned char*)
    MAC_REV_MALLOC<unsigned char>(nf.sizenf);
      //      memset(tmpnfind,0,nf.sizenf);

      for(int i=0;i<res.sizenf;i++)
      {
        if(res.nfind[i]|nf.nfind[i])
          for(int j=0;j<8;j++)
          {
        if ((res.nfind[i]>>j)&1)
          tmpspace[8*i+j]=res.nf[p++];
        if((nf.nfind[i]>>j)&1)
          tmpspace[8*i+j]-=toto*nf.nf[q++];
          }
        tmpnfind[i]=res.nfind[i]|nf.nfind[i];
      }

      for(int i=res.sizenf;i<nf.sizenf;i++)
      {
        if(nf.nfind[i])
          for(int j=0;j<8;j++)
        if ((nf.nfind[i]>>j)&1)
          tmpspace[8*i+j]=((mcoeff_t)-1)*toto*(nf.nf[q++]);
        tmpnfind[i]=nf.nfind[i];
      }
      for(int i=0;i<nf.sizenf;i++)
    if(tmpnfind[i])
      for(int j=0;j<8;j++)
        if ((tmpnfind[i]>>j)&1)
          if (!Iszero(tmpspace[8*i+j]))
        tmpnf[k++]=tmpspace[8*i+j];
          else
        tmpnfind[i]&=~(1<<j);
      MAC_REV_FREE<mcoeff_t>(tmpspace,nf.sizenf*8*sizeof(mcoeff_t));
      res.sizenf=nf.sizenf;
    }
  else
    {
      tmpspace=(mcoeff_t*)
    MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*res.sizenf*8);
      for(int i=0;i<res.sizenf*8;i++)
    tmpspace[i]=0;
      for(int i=0;i<nf.sizenf;i++)
    tmpsize+=nbbits[res.nfind[i]|nf.nfind[i]];
      for(int i=0;i<res.sizenf;i++)
    tmpsize+=nbbits[res.nfind[i]];
      tmpnf=(mcoeff_t*)
    MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*tmpsize);
      tmpnfind=(unsigned char*)
    MAC_REV_MALLOC<unsigned char>(res.sizenf);

      for(int i=0;i<nf.sizenf;i++)
    {
      if(res.nfind[i]|nf.nfind[i])
        for(int j=0;j<8;j++)
          {
        if((res.nfind[i]>>j)&1)
          tmpspace[8*i+j]=res.nf[p++];
        if((nf.nfind[i]>>j)&1)
          tmpspace[8*i+j]-=toto*nf.nf[q++];
          }
        tmpnfind[i]=res.nfind[i]|nf.nfind[i];
      }

      for(int i=nf.sizenf;i<res.sizenf;i++)
      {
        if(res.nfind[i])
          for(int j=0;j<8;j++)
          if((res.nfind[i]>>j)&1)
        tmpspace[8*i+j]=res.nf[p++];
        tmpnfind[i]=res.nfind[i];
      }

      for(int i=0;i<res.sizenf;i++)
    if(tmpnfind[i])
      for(int j=0;j<8;j++)
        if  ((tmpnfind[i]>>j)&1) 
          if (!Iszero(tmpspace[8*i+j]))
        tmpnf[k++]=tmpspace[8*i+j];
          else
        tmpnfind[i]&=~(1<<j);
        
      MAC_REV_FREE<mcoeff_t>(tmpspace,res.sizenf*8*sizeof(mcoeff_t));
    }
#ifdef DEB
  cout<<"k "<<k<<" tmpsize "<<tmpsize<<endl;
#endif
  MAC_REV_FREE<mcoeff_t>(res.nf,res.size*sizeof(mcoeff_t));
  //res.nf=tmpnf;
  MAC_REV_FREE<unsigned char>(res.nfind,res.sizenf);
  //res.nfind=tmpnfind;
  
  res.nf=tmpnf;
  res.nf=(mcoeff_t*)MAC_REV_REALLOC<mcoeff_t>(res.nf,res.size*sizeof(mcoeff_t)
                    ,k*sizeof(mcoeff_t));
                               
  //cout<<"tmpnfind "<<tmpnfind<<endl;
  res.nfind=tmpnfind;
  //  cout<<"res.nfind "<<res.nfind<<endl;
  res.size=k;
#ifdef DEB
  int cpt=0;
  for(int i=0;i<res.sizenf;i++)
    if(tmpnfind[i])
      for(int j=0;j<8;j++)
    if  ((tmpnfind[i]>>j)&1) 
      if (Iszero(res.nf[cpt++]))
        cout<<"j'ai fait n'importe quoi ici"<<endl;
  cout<<"les resultat "<<endl<<invconv<Poly>(res,serv)<<endl;
#endif
  return;
}

//#undef DEB
//renvoit l'indice du monome x_i*m_j
template<typename mon,typename Base>
int indice2(int i,int j,const Base & b)
{
 int res=-1;
 mon tmp;
 int2mon(j,tmp);
 tmp*=mon(i,1);//on suppose qd meme que j correspond a qqchose
 //cout<<"indice2 IsinB("<<mon(i,1)*tmp<<",b) "<<IsinB(mon(i,1)*tmp,b)<<endl;
 if(IsinB(tmp,b))
   {
     res=mon2int(tmp);
     //  cout<<"le monome "<<mon(i,1)*tmp<<" l'indice "<<res<<endl;
   }
 return res;
}

//on lui donne l'ensemble des formes normales, dump, l'indice de b(cf ci dessus)
// le numero de la variable par laquelle on multiplie
//et elle renvoie un pol(cf def du type) sur la forme normale de var*m_{indmon}
#ifndef uninf 
template<typename typpol,typename typdump,typename Base>
typpol nf(int var ,int indmon,const typdump &dump,const Base & b)
{
  //on cherche dans dump...
  //a optimiser plus tard
  typedef typename typdump::value_type dumpstruct;
  typedef typpol pol;
  typedef typename pol::monom_t mon;
  pol res;
  res.nf=NULL;
  res.size=-1;
  mon tmpmon;//=mon(1,var)*int2mon<mon>(indmon);
  int deg=0;//=tmpmon.GetDegree()+1;
  typename typdump::const_iterator iter;
  int2mon(indmon,tmpmon);
  deg=tmpmon.GetDegree()+1;
  //cout<<"mon("<<var<<",1) "<<mon(var,1)<<" et deg "<<deg<<endl;
  tmpmon*=mon(var,1);
   for(iter=dump.begin();
       (iter!=dump.end()) && (iter->k!=deg);iter++);
   //    cout<<"iter-k "<<iter->k<<endl;
   // if(iter!=dump.end()) cout<<"iter->k en sortie "<<iter->k<<" et deg "<<deg<<endl;
  //cout<<"valuer a la sortie "<<(--iter)++->k<<endl;
   // cout<<"iter!=dump.end() "<<(iter!=dump.end())
   //   <<" dump.empty() "<<dump.empty()<<endl;
  if(iter!=dump.end())//on cherche pas n'importe quoi qd meme
    {
      int i;
      for(i=0;i<iter->size;i++)
    {
      //  cout<<"dans nf iter->nf[i].ind "<<iter->nf[i].ind<<endl;
      if(iter->nf[i].ind.rep==tmpmon.rep)
        break;
    }
      if(i!=iter->size)
    //Tue Dec  4 16:39:25 MET 2001
    //  densify(iter->nf[i],res);

    res=(iter->nf[i]);
    }
  //cout<<"la gueule de res "<<res<<endl;
  //if (res!=NULL) cout<<"suite de res "<<res->ind<<" "<<res->nf<<endl;
  return res;
}
#endif

void getperm(int *perm, int sizeperm, int *indextab)
{
  int median=indextab[sizeperm/2],i,j=sizeperm-1;
  if(1)
  for(i=0;i<=j;i++)
    if(indextab[i]>median)
      {
    int tmp=indextab[j];
    indextab[j]=indextab[i];
    indextab[i]=tmp;
    j--;i--;
      }
  perm[sizeperm/2]=sizeperm-i;
  getperm(perm,sizeperm-i-1,indextab);
  getperm(perm+sizeperm-i,i,indextab+sizeperm-i);
} 

//effectue la multiplication du pol T par la variable x_i
template<typename typpol,typename typdump,typename Base,typename mon,
     typename workspace, typename monomial_server> 
typpol inline mult(int i,const typpol &pol,const typdump &dump,
           const Base &b,mon **exce,int &exsize,workspace &w,
           monomial_server &serv)
{
  typedef typename typpol::coeff_t coeff;
  static const mon monzero=mon(0);
  static COEFF zero=0;
  typpol ret;
  (*exce)=(mon*)MAC_REV_MALLOC<mon>((pol.size+1)*sizeof(mon));
  int *stockindicetab=(int*)malloc(pol.size*sizeof(int));
  int noquo=0,maxind=0,p=0,q=0;
  coeff *tmpspace;
  //  for(int i=0;i<w.sizenf;i++)
  //  p+=nbbits[w.tabnf[i]];
  for(int j=0;j<pol.sizenf;++j)
    if (pol.nfind[j])
     {
       for(int k=0;k<8;k++)
     if((pol.nfind[j]>>k)&1)
       {
         int stockindice2=mulind(8*j+k,i,b);
         if(stockindice2<0)
           stockindice2=-1*(8*j+k)-1;
         stockindicetab[p++]=stockindice2;
               
         //   cout<<"stockindice2 "<<stockindice2<<" pol.size "<<
         //pol.size<<" J "<<j<<endl;
         if (stockindice2>=maxind)
           maxind=stockindice2;
       }
     }
  //cout<<"p "<<p<<" pol.size "<<pol.size<<endl;
  for(;maxind/8+1>=w.sizenf;increase_space(&w));
  //tmpspace=(coeff*)MAC_REV_MALLOC<coeff>((maxind/8+1)*8*sizeof(coeff));
  //ret.sizenf=(maxind/8+1);
  //ret.nfind=(unsigned char*)MAC_REV_MALLOC<unsigned char>(ret.sizenf);
  //memset(ret.nfind,0,ret.sizenf);
    


  //for(int j=0;j<ret.sizenf*8;j++)
  //  tmpspace[j]=0;
#ifdef DEB 
  q=0;
  for(int j=0;j<p;j++)
    q+=nbbits[pol.nfind[j]];
  cout<<" a '"<<q<<endl;
#endif
  int j;
//#pragma omp parallel for shared(stockindicetab,w,serv)
  for(j=0;j<p;j++)
    if(stockindicetab[j]>=0)
      {
    
#ifdef DEB
    mon tmpmm;
    int2mon(stockindicetab[j],tmpmm);
    cout<<"le monome la "<<tmpmm<<" "<<pol.nf[j]<<endl;
#endif
    w.tabcoeff[stockindicetab[j]]+=pol.nf[j];
    if(!Iszero(w.tabcoeff[stockindicetab[j]]))
      w.tabnf[stockindicetab[j]/8]|=1<<(stockindicetab[j]%8);
      }
    else  
      {
    //  cout<<"ici aussi"<<endl;
    typpol tmp;
    tmp=serv.nf(i,-1*(stockindicetab[j]+1),dump,b);
          if (tmp.size!=-1)
        {
          int compteurtmp=0;
          //On fait un plus dense dans tmpspace
          for(;tmp.sizenf>=w.sizenf;increase_space(&w));
//                  if(tmp.sizenf>res.sizenf)
//                   {
//            ret.nfind=(unsigned char*)realloc(ret.nfind,tmp.sizenf);
//            for(int k=ret.sizenf;k<tmp.sizenf;k++)
//          ret.nfind[k]=0;
//            tmpspace=(coeff*)MAC_REV_REALLOC<coeff>(tmpspace,
//                   ret.sizenf*8*sizeof(coeff),
//                   tmp.sizenf*8*sizeof(coeff));
//            for(int k=ret.sizenf*8;k<tmp.sizenf*8;k++)
//          tmpspace[k]=0;
//            ret.sizenf=tmp.sizenf;
              
//            //normalement on ne devrait JAMAIS passer par la
//          }
#if 1
          int k,l;
          unsigned *tmpptr=(unsigned*)tmp.nfind;
          unsigned mask=1;
          for(k=0;k<tmp.sizenf/sizeof(unsigned);k++)
            if(tmpptr[k])
              {
            ((unsigned*)w.tabnf)[k]|=tmpptr[k];
            mask=1;
            for(l=0;l<8*sizeof(unsigned);l++,mask<<=1)
              {
                if(tmpptr[k]&mask)
                  w.tabcoeff[8*sizeof(unsigned)*k+l]
                -=pol.nf[j]*tmp.nf[compteurtmp++];
                //    if(Iszero(w.tabcoeff[8*k*sizeof(unsigned)+l]))
                //  ((unsigned*)w.tabnf)[k]&=~mask;
              }
              }
          for(k=tmp.sizenf-(tmp.sizenf%sizeof(unsigned))
            ;k<tmp.sizenf;k++)
            if(tmp.nfind[k])
              {
            w.tabnf[k]|=tmp.nfind[k];
            mask=1;
            for(l=0;l<8;l++,mask<<=1)
              {
                if(tmp.nfind[k]&mask)
                  w.tabcoeff[8*k+l]-=pol.nf[j]*tmp.nf[compteurtmp++];
                //if(Iszero(w.tabcoeff[8*k+l]))
                //  w.tabnf[k]&=~mask;
              }
              }
#endif
#if 0
          for(int k=0;k<tmp.sizenf;k++)
            {
              if(tmp.nfind[k])
            {
              if(tmp.nfind[k]==255)
                {
                  unsigned char mask=1;
                  w.tabnf[k]|=tmp.nfind[k];
                  for(int l=0;l<8;l++,mask<<=1)
                {
                  w.tabcoeff[8*k+l]-=
                    pol.nf[j]*tmp.nf[compteurtmp++];
                  if(Iszero(w.tabcoeff[8*k+l]))
                    w.tabnf[k]&=~mask;
                }
                }
              else
                {
                  unsigned char mask=1;
                  w.tabnf[k]|=tmp.nfind[k];
                  for(int l=0;l<8;l++,mask<<=1)
                if(tmp.nfind[k]&mask)
                  {
                    w.tabcoeff[8*k+l]-=
                      pol.nf[j]*tmp.nf[compteurtmp++];
                    if(Iszero(w.tabcoeff[8*k+l]))
                      w.tabnf[k]&=~mask;
                  }
                }
            }
              
            }
#endif
        }
          else  
        {
          mon tmpmon;
          (*exce)[noquo]=mon(i,1);
          int2mon(-1*(stockindicetab[j]+1),tmpmon);
          (*exce)[noquo]*=tmpmon;//int2mon<mon>(j);
          (*exce)[noquo]*=pol.nf[j];
          ++noquo;
        };
      }
  ret.ind=pol.ind;
  ret.ind*=mon(i,1);
  (*exce)[noquo++]=ret.ind;
  exsize=noquo;
  (*exce)=(mon*)MAC_REV_REALLOC<mon>((*exce),(pol.size+1)*sizeof(mon)
                   ,(noquo)*sizeof(mon));
#ifdef DEB  
  for(int j=0;j<exsize;j++)
    cout<<" "<<(*exce)[j];
  cout<<endl;
#endif
  ret.ind.SetCoeff(noquo+1);
  ret.size=0;
  p=0;
  for(int j=0;j<w.sizenf;j++)
    if(w.tabnf[j])
      {
    ret.size+=nbbits[w.tabnf[j]];
    p=j;
      }
#ifdef DEB
  cout<<" val p "<<p<<" et pol.size "<<pol.size<<endl;
#endif
  ret.nf=(coeff*)MAC_REV_MALLOC<coeff>(ret.size*sizeof(coeff));
  ret.sizenf=p+1;
  ret.nfind=(unsigned char *)malloc(ret.sizenf);
  p=0;
  for(int j=0;j<ret.sizenf;j++)
    {
      ret.nfind[j]=w.tabnf[j];
      if(w.tabnf[j])
    for(int k=0;k<8;k++)
      if((w.tabnf[j]>>k)&1)
        {
          if(!Iszero(w.tabcoeff[8*j+k]))
        {
          ret.nf[p++]=w.tabcoeff[8*j+k];
          w.tabcoeff[8*j+k]=zero;
        }
          else
        {
          ret.nfind[j]&=~(1<<k);
        }
        }
      w.tabnf[j]=0;
    }
  ret.nf=(coeff*)MAC_REV_REALLOC<coeff>(ret.nf,ret.size*sizeof(coeff),
                    p*sizeof(coeff));
  ret.size=p;
  //MAC_REV_FREE<coeff>(tmpspace,ret.sizenf*8*sizeof(coeff));
  free(stockindicetab);
#ifdef DEB
  cout<<"ret.size "<<ret.size<<endl;
  cout<<invconv<Poly>(ret)<<endl;
  cout<<"-----------------------------------"<<endl;
  int cpt=0;
  for(int i=0;i<ret.sizenf;i++)
    if(ret.nfind[i])
      for(int j=0;j<8;j++)
    if  ((ret.nfind[i]>>j)&1) 
      if (Iszero(ret.nf[cpt++]))
        cout<<"nimpici"<<endl;
#endif  
  return ret;
}

// renvoit mon/x_i .... mais suppose qu'on ne fait pas n'importe quoi
 template <typename Mon,typename Base>
Mon my_divmon(const Mon &mon,int i,const Base& b)
{
  Mon res(1);
  for(int j=0;j<b.nbvar();j++)
    if(i==j)
      {
    res*=Mon(j,mon.GetDegree(j)-1);
      }
     else 
       {
     res*=Mon(j,mon.GetDegree(j));
       }
  return res;
}

template<typename typmon,typename Base>
int GoodDir(const typmon &m,int i, const Base & b)
{
  int j,res=0;
  mon tmpm=m*mon(i,1);
  //cout<<"dans GoodDir les monomes "<<m<<"********"<<typmon(i,1)<<endl;
  //cout<<"b.nbvar() "<<b.nbvar()<<endl;
#pragma omp parallel for shared(tmpm,res)
  for(j=0;j<b.nbvar();j++)
    {
      if(tmpm.GetDegree(j)>=1)
    if(IsinB(my_divmon(tmpm,j,b),b)) {res|=1;};
    }
  //cout<<" le verdict de GoodDir "<<res<<endl;
  return res;
}

//fonction adhoc qui merge une liste de polyalp avec une liste de pol
//en cours de calcul cad liste de pol + liste excedentaire 
template<typename typPk,typename typP,typename typexceed,
     typename typsize ,typename Base,typename monomial_server>
void conv_merge(typPk &Pk,typP &P,typexceed &ex,typsize &sizeex,
        const Base &b, monomial_server &serv)
{
  typedef typename typP::value_type polalp;
  typedef typename polalp::order_t ord;
  typedef typename polalp::monom_t mon;
  typedef typename typPk::value_type pol;
  typedef typename pol::coeff_t coeff;
  coeff * tabcoeff=NULL;
  unsigned char *tabindices=NULL;
  int tabsize=0,tabnfsize=0;
  
  for(typename typP::iterator iter=P.begin();iter!=P.end();iter++)
    {
      //cout<<"dans for"<<endl;
      int comp=0;int i=0;
      pol tmp;
      int maxind=0;
      tmp.ind=mon(1);
      for(typename polalp::iterator iterpol=iter->begin();
      iterpol!=iter->end();iterpol++)
    if (IsinB(*iterpol,b) && (serv.mon2int(*iterpol)>maxind))
      maxind=serv.mon2int(*iterpol);
      tabcoeff=(coeff*)MAC_REV_MALLOC<coeff>((maxind+1)*sizeof(coeff));
      tabindices=(unsigned char*)malloc(maxind/8+1);
      
      tabsize=maxind+1;
      tabnfsize=maxind/8+1;
      memset((void*)tabindices,0,tabnfsize);
      maxind=0;
      //ici il faudrait un type template
      for(i=0;i<tabsize;i++)
    {
      tabcoeff[i]=0;
    }
      //memset((void*)tabcoeff,0,tabsize*sizeof(coeff));
      //ici aussi
      i=0;
      mon * tmpex=NULL;
      tmpex=(mon*)MAC_REV_REALLOC<mon>((void*)tmpex,0,iter->size()
                       *sizeof(mon));
        
    for(typename polalp::iterator iterpol=iter->begin();
        iterpol!=iter->end();iterpol++)
      {
#ifdef BONRAJOUT    
    cout<<"IsinB("<<*iterpol<<","<<"b)"<<IsinB(*iterpol,b)<<endl;
#endif     
           if (IsinB(*iterpol,b))
          {
#ifdef BONRAJOUT        
        cout<<"dedans"<<endl;
#endif
        int tmp=serv.mon2int(*iterpol);
        tabcoeff[tmp]=iterpol->GetCoeff();
        tabindices[tmp/8]|=1<<(tmp%8);
        i++;
        maxind++;
          }
        else
          { 
#ifdef BONRAJOUT    
        cout<<"y en a un qu'est pas dans B "<<*iterpol<<endl;
#endif
        if(!Iszero(iterpol->GetCoeff()))
           tmpex[comp++]=*iterpol;
          }
      }
    tmpex=(mon*)MAC_REV_REALLOC<mon>((void*)tmpex,iter->size()*sizeof(mon)
                     ,(comp+1)*sizeof(mon));
    tmp.size=maxind;
    tmp.nf=(coeff*)MAC_REV_MALLOC<coeff>(tmp.size*sizeof(coeff));
    tmp.nfind=tabindices;
    tmp.sizenf=tabnfsize;
    tmp.ind=mon(comp+1);
    maxind=0;
    for(int i=0;i<tabnfsize;i++)
      if(tabindices[i])
        {
          for(int j=0;j<8;j++)
        if (((tabindices[i]>>j)&1) && !Iszero(tabcoeff[8*i+j]))
            tmp.nf[maxind++]=tabcoeff[8*i+j];
        else
          tmp.nfind[i]&=~(1<<j);
        }
    tmp.nf=(coeff*)MAC_REV_REALLOC<coeff>(tmp.nf,tmp.size*sizeof(coeff)
                          ,maxind*sizeof(coeff));
    tmp.size=maxind;
    MAC_REV_FREE<coeff>(tabcoeff,tabsize*sizeof(coeff));
    Pk.push_back(tmp);
    ex.push_back(tmpex);
    sizeex.push_back(comp);
#ifdef BONRAJOUT
    cout<<"ici dans conv merge"<<endl;
    for(int i=0;i<comp;i++)
      cout<<tmpex[i]<<" ";
    cout<<endl;
#endif
    }
}
//member ne test que les exposants pas les coefficients
template<typename typMk>
int member(const typMk & tmpMk,const typename typMk::value_type mon)
{
  int res=0;
  for(typename typMk::const_iterator iter=tmpMk.begin();
      iter!=tmpMk.end();iter++)
    if(iter->rep==mon.rep) {res=1;break;}
  return res;
}

template<typename typMat,typename typPk,typename typP,
     typename typMk,typename Base,typename typdump,typename typrecall,
     typename workspace, typename monomial_server>
int newPkMkMatrixof(typMat &mat,typPk &Pk, typP &P,typMk &Mk,
            int &k,Base &b, typdump &dump,typrecall &recall
            ,workspace &w, monomial_server &serv)
{
  typedef typename typMk::value_type mon;
  typedef typename typP::value_type::order_t ord;
  typedef typename typPk::value_type pol;
  typedef typename pol::coeff_t coeff;
  typedef typename typP::value_type polyalp;
  typPk tmpPk,workvect;
  typMk tmpMk;
  list<mon *> exceed;
  list<int> sizeexceed;
  int nbproj=0,flagdoplus=1,*ncol;
  //prise en compte des pol pas encore introduits
  cerr<<"k on entering newpkmkmatrixof "<<k<<endl;
  typP ttmpP=Crochet(P,k+1);//
  typrecall verytmp=Crochet(recall,k+1,serv);
  recall.erase(recall.begin(),recall.end());
  if (!ttmpP.empty()) 
    {
      int tmpdeg;
      proj(ttmpP,dump,b,serv);
      selectnozero(ttmpP);
      nbproj=ttmpP.size();
      if (!ttmpP.empty())
    {
      tmpdeg=minmaxdegP(ttmpP,tmpdeg);
      if(tmpdeg<k+1)
        {
          flagdoplus=0;
          cout<<"un des pol de pretendu degre le bon se reduit mal"<<endl;
          //en fait a un polynome de degre inferieur
          typPk tmpconv;
          recover(tmpPk,Mk,dump,tmpdeg);
          cout<<"forget "<<tmpdeg<<endl;
          recall=verytmp;
          forget(b,dump,tmpdeg,recall);
          serv.compress(tmpPk,dump,b,tmpdeg);
          verytmp.erase(verytmp.begin(),verytmp.end());
          ttmpP=Crochet(ttmpP,tmpdeg);
          nbproj=ttmpP.size();
          Mk.erase(Mk.begin(),Mk.end());
          for (typename typPk::iterator iterpk=tmpPk.begin();
           iterpk!=tmpPk.end();iterpk++)
        {
          // cout<<"iterpk->ind "<<iterpk->ind
          //    <<"et Isinb(iterpk-ind,b) "<<IsinB(iterpk->ind,b)<<endl;
          if (IsinB(iterpk->ind,b))
            {
              for(typename Base::iterator iterb=b.begin()
                ;iterb!=b.end();iterb++)
            {
              iterb->refuse=(mon*)MAC_REV_REALLOC<mon>
                (iterb->refuse,(iterb->taille2)*sizeof(mon),
                 (iterb->taille2+1)*sizeof(mon));
              iterb->refuse[iterb->taille2]=iterpk->ind;
              iterb->refuse[iterb->taille2].SetCoeff(1);
              iterb->taille2++;
            }
            }
        }
          for(typename typPk::iterator iter=tmpPk.begin()
            ;iter!=tmpPk.end();iter++)
        {
          mon *tmp=(mon*)MAC_REV_MALLOC<mon>(sizeof(mon));
          *tmp=iter->ind;
          tmpMk.push_back(*tmp);
          exceed.push_back(tmp);
          sizeexceed.push_back(1);
          //      cout<<"pour Mk iter->ind "<<iter->ind<<endl;
        }
          //on vient de chuter de 
          //degre et on a reconstitue les bons mk et pk
          k=tmpdeg-1;
          //Dump(Pk,dump);
        }
    }
    }
  if(flagdoplus)
    {
#ifdef unimul
        serv.reinit_multiple();
#endif      
#ifdef uninf
    serv.reinit_nf();
#endif 

      //list<mon> dejavu;
      //construction des pol ord sup
//
// Problewme de gestion des accumulateurs ici   
int i;
//#pragma omp parallel for shared(tmpPk,tmpMk,exceed,sizeexceed) private(i)
for(i=0;i<b.nbvar();i++)
    {

      for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
        {
          //cout<<"avant gooddir l'ind "<<iter->ind<<" et i "<<i<<endl;
          if(GoodDir((*iter).ind,i,b))
        {
          mon tmpmon=iter->ind*mon(i,1);
          tmpmon.SetCoeff(1);
          //cout<<"tmpmon ici "<<tmpmon<<endl;
          if(!member(tmpMk,tmpmon))
            {
              //
              //ici aussi rendre les op atomic ne suffit pas 
              //il faut aussi faire l'ajout.
                      //#pragma omp atomic
              typename typrecall::iterator iterrec;
              
              mon * tmp=NULL;
              int tmpsize;
              //cout<<"avant mult"<<endl;
              pol stockmult=mult(i,*iter,dump,b,&tmp,tmpsize,w
                     ,serv);


                       
            //attention le coefficient de ind est la taille de 
            // l'exceed correspondant
#ifdef DEB           
            cout<<"apres mult "<<tmpsize<<endl;
#endif
//#pragma omp critical
            {
              tmpMk.push_back(tmpmon);
              tmpPk.push_back(stockmult);
              exceed.push_back(tmp);
              sizeexceed.push_back(tmpsize);
            }
            }
        }
        }
    }
      //      recall=verytmp;
    }
  //  cout<<"tmpPK.size "<<tmpPk.size()<<endl;
  //Constructyion de pol terminee
  //cout<<"ttmpP.size "<<ttmpP.size()<<endl;
  //cout<<"avant conv merge"<<endl;
  conv_merge(tmpPk,ttmpP,exceed,sizeexceed,b,serv);
  //cout<<"apres conv merge"<<endl;
  //cout<<"tmpMK==empty?? "<<tmpMk.empty()<<endl;
#ifdef DEB
     for(typename typPk::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
       cout<<"invocn"<<invconv<Poly>(*iter)<<endl;
#endif
       //   {
  //    cout<<"iter->ind.GetCoeff() "<<iter->ind.GetCoeff()-1
  //        <<" et ind "<<iter->ind<<endl;
  //   }
    //cout<<"tmpPK.size "<<tmpPk.size()<<endl;

  //Tue Nov 27 15:57:25 MET 2001
  /*
   for(typename typPk::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
    
   {
   mon tmp=iter->ind;
   tmp.SetCoeff(1);
   if((iter->ind.GetDegree()!=0)&&(!member(tmpMk,tmp)))
   {
   //cout<<"iter->ind.GetCoeff()"<<iter->ind.GetCoeff()-(coeff)1<<" "
   //    <<tmp<<endl;
   tmpMk.push_back(tmp);
   //tmpMk.back().SetCoeff(1);
   }
   }
  */
  //construction de la matrice maintenant
  //positionnement des pointeurs de matrices pour gagner en tps
   {
     int nnz=0;
     //-----------------
     //for(typename typMk::iterator iter=tmpMk.begin();iter!=tmpMk.end();iter++)
     // cout<<"mon "<<*iter<<endl;
     //-----------------
     //     cout<<"nb proj "<<nbproj<<endl;
     //for(list<int>::iterator iter=sizeexceed.begin();iter!=sizeexceed.end()
     //    ;iter++)
     // {
     // nnz+=*iter;
     //  cout<<"on a donc de monme pour celui la "<<
     //  (int)((iter->ind).GetCoeff()-1)<<" d'ou nnz "<<nnz<<endl;
     //       }
     //for(typename typPk::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
     //  {
     //  int tmp=0;
     //  for(;(iter->ind).GetCoeff()>(coeff)tmp;tmp++);
     //  nnz+=tmp-1;
     //  cout<<"on a donc de monme pour celui la "<<
     //  (int)((iter->ind).GetCoeff()-1)<<" d'ou nnz "<<nnz<<endl;
     //       }
     // nnz+=tmpPk.size()-nbproj;//le 10/05/2001 avant: nnz+=tmpMk.size()-nbproj;
     //le -nbproj rajoute le 23/01/2001!!
     //     #define DEBUG
#ifdef DEBUG
     typename typPk::iterator tmpiter=tmpPk.begin();
     list<int>::iterator itersize=sizeexceed.begin();
     for(typename list<mon*>::iterator iter=exceed.begin();iter!=exceed.end()
      ;tmpiter++,iter++,itersize++)
      {
    cout<<"pour un exceeed la taille attendue est "
        <<*itersize<<endl;
    for(int i=0;*itersize>i;i++)//;tmpiter->ind.GetCoeff()>i+1;i++)
       {
         cout<<"le gus que je regarde "<<(*iter)[i]<<endl;
         if(!member(tmpMk,(*iter)[i]))
         {
           cout<<(*iter)[i]<<endl;
           cout<<"dans ce pol "<<invconv<Poly>(*tmpiter)<<endl;
           cout<<"l'ind "<<tmpiter->ind<<endl;
           exit(-10);
         }
       }
     }
#endif
     // #undef DEBUG
     mat.nrow=tmpPk.size();
     mat.ncol=tmpMk.size();
     ncol=(int*)calloc(mat.ncol,sizeof(int));
     //memset((void*)ncol,0,(*mat)->ncol*sizeof(int));
     map<mon,int,ord> mapint;
     typename typMk::iterator itermk=tmpMk.begin();
     typename typPk::iterator iterpk=tmpPk.begin();
     typename list<mon*>::iterator iterex=exceed.begin();
     list<int>::iterator itersize=sizeexceed.begin();
     int line_num=0;
     for(int i=0;itermk!=tmpMk.end();i++,itermk++)
       {
     mapint[*itermk]=i;
       } 
     iterpk=tmpPk.begin();
     nnz=0;
     for(;iterpk!=tmpPk.end();iterpk++,iterex++,itersize++,line_num++)
       {
     nnz+=*itersize;
#ifdef DEB
     cout<<"itesize "<<*itersize<<endl;
#endif   
     // if(iterpk->ind.GetDegree()>0) 
     //  {
     //    mon tmpmon=iterpk->ind;
     //    tmpmon.SetCoeff(1);
     //   ncol[mapint[tmpmon]]++;
         //      cout<<"on est a la ligne "<<line_num<<" ai augmente "<<mapint[(iterpk->ind*(1/iterpk->ind.GetCoeff()))]<<" a cause de son monome dom "<<endl;
         // cout<<"le voila "<<(iterpk->ind*(1/iterpk->ind.GetCoeff()))<<endl;
     //}
     
    {
      // coeff tmpcoeff=(iterpk->ind).GetCoeff()-(coeff)1;
     for(int i=0;*itersize>i;i++)
       {
         //cout<<"iterpk icici ->GetCoeff()"<<iterpk->ind.GetCoeff()<<endl;
         if (mapint.count((*iterex)[i]))
           {
#ifdef DEB
         cout<<"iterex ici "<<(*iterex)[i]<<endl;
#endif      
         ++ncol[mapint[(*iterex)[i]]];
         // cout<<" ai augmente "<<mapint[(*iterex)[i]]<<"a cause de iterex "<<endl;
           }
#ifdef FAST
         else
           {
         int j,tmpsize;
         cout<<"OUILLE OUILLE OUILLE "<<(*iterex)[i]<<endl;
         mat.ncol++;
         mat.nrow++;
         mapint[(*iterex)[i]]=mat.ncol-1;
         ncol=(int*)realloc(ncol,mat.ncol*sizeof(int));
         ncol[mapint[(*iterex)[i]]]=1;
         //reste a calculer un proj de ce monome et le tour sera joue
         Poly tmppolyalp;
         mon tmpmon=(*iterex)[i];
         pol tmppol;
         mon *tmp;
         tmpMk.push_back(tmpmon);
         for(j=0;j<b.nbvar()&&(!tmpmon.GetDegree(j))
               &&IsinB(my_divmon(tmpmon,j,b),b);j++); 
         tmppolyalp=Poly(my_divmon(tmpmon,j,b));
         //proj(tmppolyalp,dump,b);
         tmppol=convert<pol>(tmppolyalp,b);
         tmpPk.push_back(mult(j,tmppol,dump,b,&tmp,tmpsize,w));
         exceed.push_back(tmp);
         sizeexceed.push_back(tmpsize);
         mat.nnz+=tmpsize;
           }
#endif
       }
    }
       }
     //rajoute le 23/01/2001
     //enleve le 17.05/2001
     //     for(int i=mat.nrow-nbproj;i<mat.nrow;i++)
     // ncol[i]--;
     //fin rajout c pour prendre en compte que lors de 
     //de la proj le monome de tete se trouve dans exceed
     //fin enlevement
     cout<<"nnz of matrix "<<nnz<<endl;
     //en general la matrice est rectangulaire
     void * toto=MAC_REV_MALLOC<coeff>(nnz*sizeof(typename typMat::coeff_t));
     //+nnz*sizeof(int)+((*mat)->ncol+1)*sizeof(int));
     if (mat.nzval!=NULL) MAC_REV_FREE<typename typMat::coeff_t>
                (mat.nzval,mat.nnz
                 *sizeof(typename typMat::coeff_t));
     mat.nnz=nnz;
     mat.nzval=toto;
     
     
     //for(int i=0;i<nnz;i++)
     //  ((typename typMat::coeff_t*)mat.nzval)[i]=(coeff)1000000000000.0;
     
     
     // mstore->rowind=(int*)((double*)toto+nnz);
     if(mat.rowind!=NULL) free(mat.rowind);
     mat.rowind=(int*)malloc(nnz*sizeof(int));
     //mstore->colptr=(int*)(mstore->rowind+nnz);//
     if(mat.colptr!=NULL) free(mat.colptr);
     mat.colptr=(int*)malloc((mat.ncol+1)*sizeof(int));
     cout<<"matrix of size "<<mat.nrow<<"x"<<mat.ncol<<endl;
   //matrice cree rest plus qu'a ecrire dedans ca
     typename typMat::coeff_t* nzval=(typename typMat::coeff_t*)mat.nzval;
     int *rowind=mat.rowind;
     int *colptr=mat.colptr;
    //mise  en place des colptr
    *colptr=0;
    for(int i=0;i<mat.ncol-1;i++)
      {
    colptr[i+1]=colptr[i]+ncol[i];
    //cout<<"nb element colonne "<<i<<" "<<ncol[i]<<endl;
    //cout<<"les valeurs de debut de colonne "<<8*colptr[i+1]<<endl;
      }
    colptr[mat.ncol]=mat.nnz;
    
    //cout<<"mat.nnz "<<mat.nnz<<" "<<mat.nnz*8<<endl;
    //fin mise en place
    //la on fait des manip symboliques pour savoir ou ecrire 
    //les coeffs de chaques colonnes ca evite d'avoir a les
    //recopier 50000 fois pour construire la matrice 
    //c bcp plus rapide
    itermk=tmpMk.begin();
    iterpk=tmpPk.begin();
    iterex=exceed.begin();
    itersize=sizeexceed.begin();
    line_num=0;
    //cout<<"nbproj "<<nbproj<<endl;
    for(;0&&(line_num < (int)(tmpPk.size()-nbproj));
    iterpk++,iterex++,itersize++,line_num++)
      //les polynome projette sont foutut au bout de la liste!!
      //d'ou le tmpPk.size()-nbproj
      {
    mon tmpmon=iterpk->ind;
    tmpmon.SetCoeff(1);
    //const unsigned int tmpplus=line_num+1,tmp=line_num;
    const unsigned int tmpplus
      =mapint[tmpmon]+1,tmp=tmpplus-1;
    //  cout<<"ici line_num"<<line_num<<" le mon dom est col "<<tmp<<endl;
    //ici couille car le ind n'est pas forcement le mon diag
    nzval[colptr[tmpplus]-ncol[tmp]]=1;
    //fin couille changer un truc
    //  cout<<"mat.nnz "<<mat.nnz<<endl;
    //cout<<"colptr[...]"<<colptr[tmpplus]<<endl;
    //cout<<"ncol[....]"<<ncol[tmp]<<endl;
    //cout<<"1 a la pos "<<colptr[tmpplus]-ncol[tmp]<<endl;
    rowind[colptr[tmpplus]-ncol[tmp]]=line_num;
    //cout<<"1 a la ligne "<<rowind[colptr[tmpplus]-ncol[tmp]]<<endl;
    ncol[tmp]--;
    {
      //coeff tmpcoeff=iterpk->ind.GetCoeff()-(coeff)1;
      for(int i=0;*itersize>i;i++)
      {
        const coeff stock=(*iterex)[i].GetCoeff();
        (*iterex)[i].SetCoeff(1);
        const unsigned int tmpplusex=mapint[(*iterex)[i]]+1;
        const unsigned int tmpex=mapint[(*iterex)[i]];
        nzval[colptr[tmpplusex]-ncol[tmpex]]
          =stock;
        rowind[colptr[tmpplusex]-ncol[tmpex]]=line_num; 
        ncol[tmpex]--;
        //    cout<<"valeur de ncol["<<tmpex<<"] "<<ncol[tmpex]<<endl;
      }
    }
    MAC_REV_FREE<mon>(*iterex,*itersize*sizeof(mon));
      }
    for(;(iterpk!=tmpPk.end());iterpk++,iterex++,itersize++,line_num++)
      {
    //const unsigned int tmpplus=line_num+1,tmp=line_num;
    {
      //coeff tmpcoeff=iterpk->ind.GetCoeff()-(coeff)1;
      for(int i=0;*itersize>i;i++)
      {
        const coeff stock=(*iterex)[i].GetCoeff();
        (*iterex)[i].SetCoeff(1);
        const unsigned int tmpplusex=mapint[(*iterex)[i]]+1;
        const unsigned int tmpex=mapint[(*iterex)[i]];
        nzval[colptr[tmpplusex]-ncol[tmpex]]
          =stock;
        rowind[colptr[tmpplusex]-ncol[tmpex]]=line_num; 
        ncol[tmpex]--;
      }
    }
    MAC_REV_FREE<mon>(*iterex,*itersize*sizeof(mon));
    //free(*iterex);
      }
    free(ncol);
   for(typename typPk::iterator iterpk=tmpPk.begin()
     ;iterpk!=tmpPk.end();iterpk++)
     iterpk->ind.SetCoeff(1);
     //iterpk->ind*=1/iterpk->ind.GetCoeff();
   //  cout<<"la matrice dans newpkmkmatrixof"<<endl;
   //cout<<" nombre de non zero "<<mat.nnz<<endl;
#ifdef DEB
   affdebug(mat);
#endif
   //cout<<"le colptr "<<endl;
   //for(int i=0;i<mat.ncol;i++)
   // cout<<mat.colptr[i]<<" "<<endl;
   //   cout<<"les indices de lignes et le nnz"<<endl;
   //for(int i=0;i<mat.nnz;i++)
   //  cout<<mat.rowind[i]<<" "<<((typename typMat::coeff_t*)mat.nzval)[i]<<endl;
   //cout<<"fin matrice"<<endl;
   Pk.swap(tmpPk);
   Mk.swap(tmpMk);
   }
   //  for(int i=0;i<mat.nnz;i++)
   // cout<<"coeff mat "<<((typename typMat::coeff_t*)mat.nzval)[i]<<" et i "<<i<<endl; 
   return flagdoplus;
}
//newpkmkmatrixof termine reste plus qu'a tester .... ;-(

template<typename typdump,typename mon,typename pol>
void findcorresppol(const mon & m,const typdump &dump,pol &res)
{
  int i,d=m.GetDegree();
  typename typdump::const_iterator iter=dump.begin();
  for(;(iter!=dump.end())&&(iter->k!=d);iter++);
  //localise la dump struct de bon degre
  if(iter!=dump.end())
    {
      for(i=0;(i<iter->size)&&(iter->nf[i].ind.rep!=m.rep);i++);
      //cout<<"monomes passe en revue pour chercher "<<m<<" : "<<iter->nf[i].ind<<endl;;
      //trouve dans la dump struct le pol qui marche
      if (i<iter->size)//estce que je connais?
    //oui super
    {
      //cout<<"g trouve un gus qui convient"<<endl;
      //for(int j=0;j<iter->nf[i].size;j++)
      // cout<<"et ses coeff sont "<<iter->nf[i].nf[j]<<endl;
      //Tue Dec  4 16:12:10 MET 2001
      //densify(iter->nf[i],res);}
      //Tue Dec  4 16:12:10 MET 2001
      res=iter->nf[i];}
      else {
    //non ben alors il se reecrit en lui meme
    res.ind=mon(0);//m*(1/m.GetCoeff());
    res.size=-1;//pour proj avant ct -1
    res.nf=NULL;
    res.nfind=NULL;
    res.sizenf=0;
      }
    }
  return;
}
template<typename iterator, typename mon,typename Base>
void varmult(iterator iter1,iterator iter2,mon &var1, mon &var2,const Base &b)
{
  for(int i=0;i<b.nbvar();i++)
    {
      if(iter1->ind.GetDegree(i)>iter2->ind.GetDegree(i))
    var2=mon(i,1);
      if(iter2->ind.GetDegree(i)>iter1->ind.GetDegree(i))
    var1=mon(i,1);
    }
}

template <typename mon,typename pol, typename typdump,typename Base, 
      typename monomial_server>
pol halfred(const mon &targetind, const pol & p, const typdump & dump
        , const Base & b, monomial_server & serv)
{
  typedef typename pol::coeff_t coeff;
  static const mon monzero=mon(0);
  static const typename pol::coeff_t zero(0);
  int *place=(int*)MAC_REV_MALLOC<int>(p.size*sizeof(int)),comp=0;
  pol res;
#ifdef uninf
  int varref;
#else
  mon varref;
#endif 
  int i,maxref=0,ppp=0,nbbitsun=0;
  coeff *tmpspace;
  for(i=0;i<b.nbvar();i++)
    {
      if(targetind.GetDegree(i)>p.ind.GetDegree(i))
    break;
    }
#ifdef uninf
  varref=i;
#else
  varref=mon(i,1);
#endif
  
  //  for(i=0;i<p.size;i++)
  // cout<<p.nf[i]<<" ";
  //cout<<endl;
  //cout<<"ppp "<<ppp<<" p.size "<<p.size<<endl;
  //cout<<"p en entree"<<invconv<Poly>(p)<<endl;
  //cout<<"mulind(0,1,b) "<<mulind(0,1,b)<<endl;
  //cout<<"target "<<targetind<<" et ind "<<p.ind<<" et var "<<varref<<endl;

  for(i=0;i<p.sizenf;i++)
    {
      if(p.nfind[i])
    {
#ifdef DEB  
      cout<<nbbitsun<<" p.nfind[i] "<<(int)p.nfind[i]<<endl;
#endif
      nbbitsun+=nbbits[p.nfind[i]];
      for(int j=0;j<8;j++)
        {
          if((p.nfind[i]>>j)&1)
        {
#ifdef unimul   
          
          place[comp]=serv.mulind(8*i+j,varref,b);
          if(place[comp]<0)
            place[comp]=-1*(8*i+j+1);
          if (place[comp]>maxref)
            maxref=place[comp];
          
#else
          //ATTENTION UNIQUEMENT UNIMUL SINON PAS DEBUGGE 
          int indtmpmon;
          mon tmpmon;
          int2mon(8*i+j,tmpmon);
          tmpmon*=varref;
          if(IsinB(tmpmon,b))
            {
              //int j;
              indtmpmon=mon2int(tmpmon);
              //for(j=0;j<p.size&&p.nfind[j]!=indtmpmon;j++);
              //if(j<p.size) 
              maxref++;
              place[comp]=indtmpmon;
            }
          else
            {
              place[comp]=-1*(8*i+j);
            }
#endif
          comp++;
        }
        }
    }
    }
  tmpspace=(coeff*)MAC_REV_MALLOC<coeff>((maxref/8+1)*8*sizeof(coeff));
#ifdef DEB
  cout<<"comp "<<comp<<" et p.size "<<p.size<<endl;
#endif 
  res.sizenf=(maxref/8+1);
  for(int j=0;j<8*res.sizenf;j++)
    tmpspace[j]=0;
  res.nfind=(unsigned char*)MAC_REV_MALLOC<unsigned char>(res.sizenf);
  memset(res.nfind,0,res.sizenf);
  for(int j=0;j<nbbitsun;j++)
    if(place[j]>=0)
      {
#ifdef DEB
    mon tmpmm;
    int2mon(place[j],tmpmm);
    cout<<"le monome la "<<tmpmm<<" "<<p.nf[j]<<endl;
#endif
    tmpspace[place[j]]+=p.nf[j];
    res.nfind[place[j]/8]|=1<<(place[j]%8);
      }
    else  
      {
    pol tmp;
#ifdef unimul      
    tmp=serv.nf(varref,-1*place[j]-1,dump,b);
#else
          mon tmpmon;
          int2mon(-1*place[j],tmpmon);
          tmpmon*=mon(i,1);
          tmp=nf<typpol>(tmpmon,dump,b);
#endif   

          if (tmp.size>=0)
        {
          int compteurtmp=0;
          //cout<<" je reduit prr lui"<<endl<<invconv<Poly>(tmp)<<endl;
          //On fait un plus dense dans tmpspace
          if(tmp.sizenf>res.sizenf)
            {
              res.nfind=(unsigned char*)realloc(res.nfind,tmp.sizenf);
              for(int k=res.sizenf;k<tmp.sizenf;k++)
            res.nfind[k]=0;
              tmpspace=(coeff*)MAC_REV_REALLOC<coeff>(tmpspace,
                     res.sizenf*8*sizeof(coeff),
                     tmp.sizenf*8*sizeof(coeff));
              for(int k=res.sizenf*8;k<tmp.sizenf*8;k++)
            tmpspace[k]=0;
              res.sizenf=tmp.sizenf;
              
            }
          for(int k=0;k<tmp.sizenf;k++)
            {
              if(tmp.nfind[k])
            {
              for(int l=0;l<8;l++)
              if((tmp.nfind[k]>>l)&1)
                tmpspace[8*k+l]-=p.nf[j]*tmp.nf[compteurtmp++];
            }
              res.nfind[k]|=tmp.nfind[k];
            }
        }
          else
        {
          mon tmpmon;
          serv.int2mon(-1*place[j]-1,tmpmon);
          tmpmon*=mon(varref,1);
          for(int i=0;i<serv.sizeplacemon();i++)
            cout<<i<<" "<<serv.int2mon_[i]<<endl;
          cout<<"monom "<<tmpmon<<" "<<-1*place[j]
              <<" "<<p.nf[j]<<endl;
          cerr<<"Je n'aurrai pas du passer par la "<<endl;
          exit(10);
        };
      } 
 res.ind=mon(0);
 comp=0;
  for(int j=0;j<res.sizenf;j++)
    if(res.nfind[j])
      comp+=nbbits[res.nfind[j]];
  res.nf=(coeff*)MAC_REV_MALLOC<coeff>(comp*sizeof(coeff));
  res.size=comp;
  comp=0;
  for(int j=0;j<res.sizenf;j++)
    if(res.nfind[j])
      for(int k=0;k<8;k++)
    if((res.nfind[j]>>k)&1)
      if(tmpspace[8*j+k]!=0)
        {res.nf[comp++]=tmpspace[8*j+k];}
      else
        res.nfind[j]&=~(1<<k);
  res.nf=(coeff*)MAC_REV_REALLOC<coeff>(res.nf,res.size*sizeof(coeff)
                    ,comp*sizeof(coeff));
  res.size=comp;
  MAC_REV_FREE<coeff>(tmpspace,res.sizenf*8*sizeof(coeff));
  free(place);
  //cout<<" apres reduction "<<invconv<Poly>(res)<<endl;
 return res;
}

template <typename typPk, typename typdump,typename Base,
      typename workspace ,typename monomial_server>//,typename typMk>
void ReduceSpol(typPk &S,const typdump & dump,typPk &redspol,const Base &b
        ,workspace &w, monomial_server &serv)
        //, typMk & Spoldejavu)
{
  typedef typename typPk::value_type pol;
  typedef typename pol::monom_t mon;
  typedef typename pol::coeff_t coeff;
  list<mon> leadmon;
  list<typPk> membermon;
  cout<<"en entree de RedSpol S.size "<<S.size()/2<<endl;

  for(typename typPk::iterator iterS=S.begin();iterS!=S.end();iterS++,iterS++)
    {
      typename typPk::iterator iterS2=iterS;
      typename list<mon>::iterator itermon=leadmon.begin();
      typename list<typPk>::iterator itermemb=membermon.begin();
      mon var1,var2,tmpleadmon;
      iterS2++;
      varmult(iterS,iterS2,var1,var2,b);
      tmpleadmon=iterS->ind*var1;
      // Spoldejavu.push_back(tmpleadmon);
      for(;itermon!=leadmon.end()&&itermon->rep!=tmpleadmon.rep;
      itermon++,itermemb++);
      if(itermon!=leadmon.end())
    {
      int yes1=1,yes2=1;
      for(typename typPk::iterator tmpiter=itermemb->begin();
          tmpiter!=itermemb->end();tmpiter++)
        {
          if(tmpiter->ind.rep==iterS->ind.rep) yes1=0;
          if(tmpiter->ind.rep==iterS2->ind.rep) yes2=0;
        }
      if(yes1) itermemb->push_back(*iterS);
      if(yes2) itermemb->push_back(*iterS2);
    }
      else
    {
      typPk verytmp;
      leadmon.push_back(tmpleadmon);
      //cout<<"le tmpleadmon "<<tmpleadmon<<endl;
      membermon.push_back(verytmp);
      membermon.back().push_back(*iterS); 
      membermon.back().push_back(*iterS2);
    }
    }
  //maintenant on a mis en place la structure de Redspol
  //reste plus qu'a calculer les bons et a reduire tout ca!!
  { 
    typename list<mon>::iterator itermon=leadmon.begin();
    typename list<typPk>::iterator itermemb=membermon.begin();
    for(;itermon!=leadmon.end();itermon++,itermemb++)
      {
    typename typPk::value_type ref;
    typename typPk::iterator tmpiter=itermemb->begin();
    ref=halfred(*itermon,*tmpiter,dump,b,serv);
    //  /cout<<"ref.size "<<ref.size<<" "<<ref.ind<<endl;
    //cout<<"ref en entier"<<endl<<invconv<Poly>(ref)<<endl;
    tmpiter++;
    for(;tmpiter!=itermemb->end();tmpiter++)
      {
        typename typPk::value_type redcur;
        redcur=halfred(*itermon,*tmpiter,dump,b,serv);
        my_plus(redcur,ref,(coeff)1,&w,serv);//ICIOMP
        if(!Ispolzero(redcur))
          {
        redspol.push_back(redcur);
        //cout<<"redcur "<<endl<<invconv<Poly>(redcur)<<endl;
        //cout<<"le leadmon qui donne un non zero "<<*itermon<<endl;
          }
        else
          {
        MAC_REV_FREE<coeff>(redcur.nf,redcur.size*sizeof(coeff));
        MAC_REV_FREE<unsigned char>(redcur.nfind,redcur.sizenf);
          }
      }
    MAC_REV_FREE<coeff>(ref.nf,ref.size*sizeof(coeff));
    MAC_REV_FREE<unsigned char>(ref.nfind,ref.sizenf);
      }
  }
}

template<typename typPk,typename monomial_server>
int mindeg(const typPk &l,monomial_server &serv)
{
  typedef typename typPk::value_type::monom_t mon;
  typedef typename mon::coeff_t coeff;
  int res=-1;
  //  cout<<"dans mindeg res au debut"<<res<<endl;
  for(typename typPk::const_iterator iter=l.begin() ;iter!=l.end();iter++)
   {
     int mmax=-1;
     // cout<<"tour de boucle de mindeg "<<endl;
     for(int i=0;i<iter->sizenf;i++)
       if(iter->nfind[i])
     for(int j=0;j<8;j++)
       if((iter->nfind[i]>>j)&1)
       {
         mon tmpmon;
         serv.int2mon(8*i+j,tmpmon);
         mmax=(mmax<tmpmon.GetDegree())?tmpmon.GetDegree():mmax;
       }
     // cout<<"mmax ici "<<mmax<<endl;
     if (iter->ind.GetCoeff()==(coeff)1)
       mmax=(mmax<iter->ind.GetDegree())?iter->ind.GetDegree():mmax;
     //cout<<"les gens que je regarde "<<iter->ind<<endl;
     if(mmax!=-1)
     {
       if(res==-1) res=mmax;
       res=(mmax>res)?res:mmax;
     }
   }
  //cout<<"resultat de mindeg "<<res<<endl;
  return res;
}

template<typename coeff,typename typiterpk, typename workspace,typename typserv>
void inline update_crochetchoix(typiterpk &iterPk, typiterpk &iter,int tmpind
                ,workspace &w, typserv &serv)
{
 
  if((iterPk->sizenf>tmpind/8)&&(((iterPk->nfind[tmpind/8])>>(tmpind%8))&1))
    {
      //      cout<<"g trouve un gus qui a un coeff du mon choisi "<<endl;
      //cout<<"le voila "<<invconv<Poly>(*iterPk)<<endl;
        int i=0;
        coeff tmp;
        for(int j=0;8*j<=tmpind;j++)
          i+=nbbits[iterPk->nfind[j]];
        for(int j=7;j>=(tmpind % 8);j--)
          if((iterPk->nfind[tmpind/8]>>j)&1)
            i--;
        
#ifdef DEB
        cout<<"i ici "<<i<<" et iterPk->size "<<iterPk->size<<endl;
#endif
        tmp=iterPk->nf[i];
        //  cout<<"le coeff "<<tmp<<" tmpind "<<tmpind<<endl;
        //cout<<"les voisins "<<iterPk->nf[i-1]<<" "<<iterPk->nf[i+1]<<endl;
        
        //cout<<"les voisins plus loin"<<iterPk->nf[i-2]<<" "<<iterPk->nf[i+2]<<endl;
        for(int j=i;j<iterPk->size-1;j++){
          iterPk->nf[j]=iterPk->nf[j+1];
        }
        iterPk->nfind[tmpind/8]&=~(1<< (tmpind%8));
        iterPk->nf=(coeff*)MAC_REV_REALLOC<coeff>(iterPk->nf,
            iterPk->size*sizeof(coeff),(iterPk->size-1)*sizeof(coeff));
        iterPk->size--;
        //iterPk->nf[i]=0;
        my_plus(*iterPk,*iter,tmp,&w,serv);//ICIOMP
        
        //cout<<"apres plus le coeff vaut "
        //    <<iterPk->nf[mon2int(iter->ind)]<<endl; 
          }
  //cout<<"apres mise a jour"<<endl<<invconv<Poly>(*iterPk)<<endl;
}
template<typename typP,typename typPk,typename typMk
         ,typename Base,typename typdump, typename typrecall,
     typename workspace ,typename monomial_server>
void NewCrochetchoix(typP & P,typPk & Pk,typMk & Mk, typPk & redspol
             ,int &k,typdump &dump, Base &b, typrecall & recall
             ,workspace &w, monomial_server &serv)
             //,typMk &Spoldejavu)
  //template<typename typP,typename typPk,typename typMk,typename Base,typename typdump>
//void NewCrochetchoix(const typP & P,typPk & Pk,typMk & Mk, typPk & redspol
//           ,int &k,typdump &dump, Base &b)
  {
  typedef typename typP::value_type polyalp;
  typedef typename polyalp::monom_t mon;
  typedef typename polyalp::coeff_t coeff;
  //for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
  //  P.push_back(invconv<Poly>(*iter));
  //cout<<"dans newcrochetchoix"<<endl;
  /*
    
    //cout<<"les sizes "<<iter->size<<endl;
    {
    int cpt=0;
    cout<<invconv<typename typP::value_type>(*iter)<<endl;       
      for(int i=0;i<iter->sizenf;i++)
    for(int j=0;j<8;j++)
      if (((iter->nfind[i]>>j)&1) && (iter->nf[cpt++]==0))
        {
          cout<<"ici ca couille a la place "<<8*i+j<<endl;
          exit(1);
        }
    }
  */
  typPk res=Crochet(redspol,k,serv),tmpres;
  //cout<<"redspol.size "<<redspol.size()<<endl;
  //cout<<"res.size "<<res.size()<<endl;
  for(typename typPk::iterator iter=res.begin();iter!=res.end();iter++)
    {
      for(typename typPk::iterator iterredspol=redspol.begin();
          iterredspol!=redspol.end();iterredspol++)
        if(iter->nf==iterredspol->nf)
          {
        redspol.erase(iterredspol);
        break;
          }
    }
  NewDestroyredspol(redspol,Pk);
  redspol.erase(redspol.begin(),redspol.end());
  //for(typename typPk::iterator iter=res.begin();iter!=res.end();iter++)
  //   cout<<"res->sizenf<< "<<iter->sizenf<<endl;
  for(typename typPk::iterator iter=res.begin();iter!=res.end();iter++)
    {
      //cout<<"un nouveau tour de boucle "<<endl;
      if (!Ispolzero(*iter))
    {
      int tmpind,i;
      iter->ind=choice(*iter,b,serv);
      //cout<<"le pol pas zero "<<(invconv<polyalp>(*iter))<<endl;
      //cout<<"iter -> size "<<iter->size<<" "<<iter->ind<<endl;
      //for (int i=0;i<iter->size;i++)
      // cout<<iter->nf[i]<<" "<<iter->nfind[i]<<endl;
      if(iter->ind.GetDegree()<k)
        {
         //cout<<"le pol qui fait chuter"<<(invconv<polyalp>(*iter))<<endl;
          k=iter->ind.GetDegree();
          iter->ind=mon(0);
          //cout<<"passage par bout de code en plus "<<endl;
          my_merge(res,tmpres);
          my_merge(redspol,res);
          tmpres.erase(tmpres.begin(),tmpres.end());
          res.erase(res.begin(),res.end());
          res=Crochet(redspol,k,serv);
          //cout<<"newres.size "<<res.size()<<endl;
          for(typename typPk::iterator iterres=res.begin()
            ;iterres!=res.end();iterres++)
        {
          for(typename typPk::iterator iterredspol=redspol.begin();
              iterredspol!=redspol.end();iterredspol++)
            if(iterres->nf==iterredspol->nf)
              {
            redspol.erase(iterredspol);
            break;
              }
        }
          for(typename typPk::iterator iterredspol=redspol.begin();
              iterredspol!=redspol.end();iterredspol++)
        if(Degree(*iterredspol,serv)<=k)
          {
            tmpres.push_back(*iterredspol);
            redspol.erase(iterredspol--);
          }
          NewDestroyredspol(redspol,Pk);
          redspol.erase(redspol.begin(),redspol.end());
          my_merge(redspol,tmpres);
          tmpres.erase(tmpres.begin(),tmpres.end());
          iter=res.begin();
          iter->ind=choice(*iter,b,serv);
          
          Dump(Pk,dump);//c pour detruire Pk dans forget
                  //ca ne respecte pas la structure du dump!!!
          recover(Pk,Mk,dump,k);
          // Listedejafait(Mk,Spoldejavu);
          forget(b,dump,k,recall);
          //forget(b,dump,k);
          //          cout<<"un coup de compress"<<endl;
          //ce compresss la vire le 06 10 2006
          serv.compress(Pk,redspol,dump,b,k);
          //pose un probleme de oherence pour le int2mon 
          //si on le laisse 
          //swap(*iter,*res.rbegin());//on veut ne plus 


          //iter=--res.end();//refaire d'autre tours
        }
      //cout<<"g choisi dans Crochetchoix "<<iter->ind<<endl;
      //cout<<"a t'il un mon2int ?"<<IsinB(iter->ind,b)<<endl;
      //mon tmp10,tmp11;
      //      int2mon(10,tmp10);int2mon(11,tmp11);
      //cout<<"qu'en est il de int2mon(10) et 11 "<<tmp10<<" "<<tmp11<<endl;
      //      
      //
      //
      //UGLY HACK
      remiseenforme(*iter);
      //SALE HACK
      //
      //
      //
      //
      //
      *iter/=iter->ind.GetCoeff();//pour la "nf"
      //cout<<"le coeff "<<iter->ind<<endl;
      iter->ind*=(coeff)1/iter->ind.GetCoeff();//pour l'ind
      tmpind=serv.mon2int(iter->ind);
      i=0;
      for(int j=0;8*j<=tmpind;j++)
        i+=nbbits[iter->nfind[j]];
      for(int j=7;j>=(tmpind % 8);j--)
        if((iter->nfind[tmpind/8]>>j)&1)
          i--;
      for(int j=i;j<iter->size-1;j++)
        iter->nf[j]=iter->nf[j+1];
      iter->nf=(coeff*)MAC_REV_REALLOC<coeff>(iter->nf,
                         iter->size*sizeof(coeff),
                         (iter->size-1)*sizeof(coeff));
      iter->nfind[tmpind/8]&=~(1<<(tmpind%8));
      iter->size--;
#ifdef DEB
      cout<<"apres tambouille"<<endl<<invconv<Poly>(*iter)<<endl;
#endif

      //      
      //
          //
      //UGLY HACK
      remiseenforme(*iter);
      //SALE HACK
          //
      //
          //
      //
      //
      
      tmpres.push_back(*iter);
      //on met a 0 le coeff du mon passe en ind
      //maintenant on a une regle en forme comme il faut reste a:
      //-reduire les regles deja existantes par rapport a celle ci
      //-reduire les autres redspol par rapport a celle ci
      for(typename typPk::iterator iterPk=Pk.begin();iterPk!=Pk.end();iterPk++)
        update_crochetchoix<coeff>(iterPk,iter,tmpind,w,serv);
      typename typPk::iterator iterredspol=iter;
      iterredspol++;
      for(;iterredspol!=res.end();iterredspol++)
        {
          //cout<<"redspol.size() "<<redspol.size()<<endl;
          //cout<<"iterredspol-ind "<<iterredspol->ind<<" "<<iterredspol->nf<<endl;
          update_crochetchoix<coeff>(iterredspol,iter,tmpind,w,serv);
        }
      //au tour de tmpres...
      typename typPk::iterator iterend=tmpres.end();
      iterend--;
      for(iterredspol=tmpres.begin();iterredspol!=iterend;iterredspol++)
        {
          //cout<<"redspol.size() "<<redspol.size()<<endl;
          //cout<<"iterredspol-ind "<<iterredspol->ind<<" "<<iterredspol->nf<<endl;
          update_crochetchoix<coeff>(iterredspol,iter,tmpind,w,serv);
        }
          
      //voila c fait 
      //maintenant on met a jour la Base b.
      for(typename Base::iterator iterb=b.begin();iterb!=b.end();iterb++)
        {
          int tmptaille1=0;
          mon *tmp=(mon*)MAC_REV_MALLOC<mon>(iterb->taille1*sizeof(mon));
          for(int i=0;i<iterb->taille1;i++)
        if(!isdivisible(iter->ind,iterb->accept[i],b))
          tmp[tmptaille1++]=iterb->accept[i];
          tmp=(mon*)MAC_REV_REALLOC<mon>(tmp,iterb->taille1*sizeof(mon)
                         ,tmptaille1*sizeof(mon));
          MAC_REV_FREE<mon>(iterb->accept,iterb->taille1*sizeof(mon));
          iterb->accept=tmp;
          iterb->taille1=tmptaille1;
          if(iterb->taille1==0)
        {
          MAC_REV_FREE<mon>(iterb->refuse,iterb->taille2*sizeof(mon));
          b.erase(iterb--);
          continue;
        }
          iterb->refuse=(mon*)MAC_REV_REALLOC<mon>
        (iterb->refuse,(iterb->taille2)*sizeof(mon),
         (iterb->taille2+1)*sizeof(mon));
          iterb->refuse[iterb->taille2]=iter->ind;
          // cout<<"reffuse rajoute "<<iter->ind<<endl;
          iterb->refuse[iterb->taille2].SetCoeff(1);
          iterb->taille2++;
        }
      res.erase(iter--);
    }
      //break;//comme ca on fait pas 50 tours.
    }
  for (typename typPk::iterator iterpk=Pk.begin();iterpk!=Pk.end();iterpk++)
    {
      //cout<<"iterpk->ind "<<iterpk->ind<<"et Isinb(iterpk-ind,b) "<<IsinB(iterpk->ind,b)<<endl;
      //cout<<"le pol "<<invconv<Poly>(*iterpk)<<endl;
      if (IsinB(iterpk->ind,b))
    {
      for(typename Base::iterator iterb=b.begin();iterb!=b.end();iterb++)
        {
          iterb->refuse=(mon*)MAC_REV_REALLOC<mon>
        (iterb->refuse,(iterb->taille2)*sizeof(mon),
         (iterb->taille2+1)*sizeof(mon));
          iterb->refuse[iterb->taille2]=iterpk->ind;
          iterb->refuse[iterb->taille2].SetCoeff(1);
          iterb->taille2++;
        }
    }
    }
#ifdef NO_TWICE_FALL
  for(typename typPk::iterator iter=tmpres.begin();iter!=tmpres.end();iter++)
    P.push_back(invconv<polyalp>(*iter));
#endif
#ifdef BONRAJOUT
  for(typename typPk::iterator iter=tmpres.begin();iter!=tmpres.end();iter++)
    P.push_back(invconv<polyalp>(*iter));
  cout<<"a la fin de §New"<<endl;
#endif
#ifdef DEB
  for(typename typPk::iterator iter=tmpres.begin();iter!=tmpres.end();iter++)
    cout<<"tmpres"<<endl<<invconv<polyalp>(*iter)<<endl;
  cout<<"tmpres.size "<<tmpres.size()<<endl;
#endif 
 my_merge(Pk,tmpres);//les gens dans res\tmpres ne sont pas detruits ici!!
  my_merge(redspol,res);
  
  
  /*  for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
       //cout<<"les sizes "<<iter->size<<endl;
      {
    int cpt=0;
    cout<<invconv<typename typP::value_type>(*iter)<<endl;       
    for(int i=0;i<iter->sizenf;i++)
      for(int j=0;j<8;j++)
        if (((iter->nfind[i]>>j)&1) && (iter->nf[cpt++]==0))
          {
        cout<<"ici ca couille a la place "<<8*i+j<<endl;
        exit(1);
          }
      }
  */
  return;
}
//#undef DEB
template<typename typPk>
void NewDestroyredspol(typPk &redspol,const typPk &Pk)
{
//  typename typPk::const_iterator iterpk=--Pk.end(); 
   for(typename typPk::iterator iter=redspol.begin();
      iter!=redspol.end();iter++)
    {
      //    if(iter->nf!=iterpk->nf)
      MAC_REV_FREE<typename typPk::value_type::coeff_t>(iter->nf,iter->size*sizeof(typename typPk::value_type::coeff_t));
      MAC_REV_FREE<unsigned char>(iter->nfind,iter->sizenf);
    }
}

template<typename typMk>
void Listedejafait(const typMk &Mk,typMk &Spoldejavu)
{
  Spoldejavu.erase(Spoldejavu.begin(),Spoldejavu.end());
  for(typename typMk::const_iterator iter1=Mk.begin();iter1!=Mk.end();
      iter1++)
    for(typename typMk::const_iterator iter2=iter1;iter2!=Mk.end();
      iter2++)
      {
    typename typMk::value_type tmp=lcm(*iter1,*iter2);
    if (tmp.GetDegree()==iter1->GetDegree()+1)
      Spoldejavu.push_back(tmp);
      }
}

template<typename typP,typename typPk,typename typdump,typename Base
     , template<typename  T> class tmptypmat ,typename monomial_server>
void algo(typP P,typdump &dump,Base &b, monomial_server &serv)
{
#ifdef PEQUAN
  int toto=0;
  FILE *fic;
  char name[8];
#endif
  typedef typename typP::value_type polyalp;
  typedef list<typename typPk::value_type::monom_t> typMk;
  typedef typename typMk::value_type mon;
  typedef typename typPk::value_type pol;

  typMk Mk;
  typMk Spoldejavu;
  typPk Pk,S;
  //typP recall 
  typPk recall;
  //Base b;
  workspace<typename typP::value_type::coeff_t> w;
  typedef tmptypmat<typename typP::value_type::coeff_t> typMat;
  typMat Mat,L;
  int nouveaumon,allpolused=0,maxdeg,k,*pr,*pc;
  //cout<<"avant init"<<endl;
  
  init_algo<typMat>(P,nouveaumon,maxdeg,k,Pk,Mk,b,dump,w,serv);
  //cout<<"apres init"<<endl;
#ifndef GROEBNER

  //  while((critere(b,k)==-1)||(k<=maxdeg))
  while(nouveaumon||(k<=maxdeg))
#else
    while(!GB_stop(k)||(k<=maxdeg))
#endif
 
    {
    //if (toto>=5) {exit(0);}
    //toto++;
        
    int killspol;
#ifdef BONRAJOUT
      for(typename Base::const_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;
      }
        }
#endif
#ifdef DEB

    //cout<<"mon a la fin de la boulcle de l'algo "<<Mk.size()<<endl;
      //for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
      // cout<<*iter<<endl;  
  cout<<"poly a la fin de la boulcle de l'algo "<<Pk.size()<<endl;
    for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
       //cout<<"les sizes "<<iter->size<<endl;
      {
    int cpt=0;
    cout<<invconv<typename typP::value_type>(*iter)<<endl;       
    for(int i=0;i<iter->sizenf;i++)
      for(int j=0;j<8;j++)
        if (((iter->nfind[i]>>j)&1) && (iter->nf[cpt++]==0))
          {
        cout<<"ici ca couille a la place"<<8*i+j<<endl;
        cout<<iter->sizenf<<" "<<iter->size<<endl;
        exit(1);
          }
      }
#endif   

    //cout<<"Spolordsup"<<endl;
    S=Spolordsup(Pk,b,Spoldejavu);
    //S.erase(S.begin(),S.end());
    //S=Spolordsup(Pk,b);
        cout<<endl;
    //<<endl<<"newpkmkmatrixof"<<endl;
    //cout<<"Pk.size() "<<Pk.size()<<endl;
    //killspol=!newPkMkMatrixof(Mat,Pk,P,Mk,k,b,dump,recall,Spoldejavu);
    killspol=!newPkMkMatrixof(Mat,Pk,P,Mk,k,b,dump,recall,w,serv);
#ifdef PEQUAN
    sprintf(name,"mat%d",toto);
    fic=fopen(name,"w");
    printmatfile_pequan(Mat,fic);
    fclose(fic);
    toto++;
#endif
    //    if (toto>=5) {cout<<"killredspol "<<killspol<<endl; exit(0);}
    // killspol=!newPkMkMatrixof(Mat,Pk,P,Mk,k,b,dump);
    if(killspol) S.erase(S.begin(),S.end());
    if (k+1>=maxdeg) allpolused=1;
    if(Mk.size()==0) nouveaumon=0;
    //cout<<"j'en suis a k="<<k<<endl;
    //cout<<"Ici dans algo Mat"<<endl;
#ifdef BONRAJOUT
    affdebug(Mat);
#endif
    //Solve affecte a Mat la valeur U
    //cout<<"solve"<<endl;
    //for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
    //  cout<<"avant solve les ind "<<iter->ind<<endl;
    //for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
    //  cout<<"avant solve les pol "<<invconv<Poly>(*iter)<<endl;
    //if(k==3) exit(-10);
    //for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
    //  cout<<"monomes dans Mk "<<*iter<<endl;
    
    if (Mat.ncol!=0)
      {
    //cout<<"debut du solve"<<endl;
    Solve(Mat,L,Pk,&pr,&pc,w);
    cout<<"fin du solve"<<endl;
    //for(int kk=0;kk<Mk.size();kk++)
    //  cout<<"permc "<<kk<<" "<<pc[kk]<<endl;
    inverse(pc,Mk.size());
    Dopermu(Mk,pc);
    {
      //ici on remet les ind a la bonne valuer 
      //a inclure dans une petite fonction inline apres
      int j=0;
      typename typMk::iterator itermon=Mk.begin();
      typename typPk::iterator iter;
      for(iter=Pk.begin();
          itermon!=Mk.end();iter++,itermon++,j++)
        {
          if(!Iszero(Mat(j,j)))
        {
          //          cout<<"ligne "<<j<<" je mets le monome "<<*itermon<<endl;
          iter->ind=*itermon*(iter->ind.GetCoeff());
        }
          else
        {
          iter->ind=mon(0);
        }
        }
      //for(;itermon!=Mk.end();iter++,itermon++)
      //  iter->ind=mon(0);
    }

    free(pr);
    free(pc);
      }

 //appliquer les permutation a MK!!!!!!!!
    // for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
      //  cout<<"apres solve les ind "<<iter->ind<<endl;
    //for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
    // cout<<"apres solve les pol "<<invconv<Poly>(*iter)<<endl;
    //
    //
    //APPLIQUER LES PERMUTYATIONS A MK!!
    //
    //
    //
    //
    if ((unsigned)rank(Mat)==Mk.size())
      //matrice de rang plein
      {
    typMk tmpSpoldejavu;
    typPk redspol,secd=secondmembre(Mat,Pk);;
    //secondmembre2(Mat,Pk,secd);
    //typename typPk::iterator iter=Pk.begin();
    cout<<"matrice de rang plein "<<endl;
    //for(unsigned int i=0;i<Mk.size();i++,iter++);
    //Pk.erase(iter,Pk.end());
    //cout<<"Dump "<<endl;
    Dump(Pk,dump);
    //cout<<"Reduce Spol "<<endl;
    
    ReduceSpol(S,dump,redspol,b,w,serv);//tmpSpoldejavu);
    //cout<<"secd.size() "<<secd.size()<<endl;
    my_merge(redspol,secd);
    if(redspol.size()==0)
      {
        //a voir....    Dump(Pk,dump);
        k++;
        my_merge(Spoldejavu,tmpSpoldejavu);
        L.destroystore();
        continue;
        //la matrice a rang plein 
        //ET les Spol\secondmembre de ce niveau se ,reduisent a 0
      }
    else
      {
        //on a des spol\secondmembre qui se reduisent pas a 0
        //cad un defaut de commuatation
        int tmpdeg=mindeg(redspol,serv);
        max_degppcm=k+2;
        nouveaumon=1;
        cout<<"redspol.size() "<<redspol.size()<<endl;
        //cout<<invconv<typename typP::value_type>(redspol.front())<<endl;
        if (tmpdeg==(k+1))
          my_merge(Spoldejavu,tmpSpoldejavu);
        k=tmpdeg;
        //cout<<"le mindeg ici vaut "<<k<<endl;
        //cout<<"Pk.size() avant recover "<<Pk.size()<<endl;
        recover(Pk,Mk,dump,k);

        //Listedejafait(Mk,Spoldejavu);
        //cout<<"Pk.size() apres recover "<<Pk.size()<<endl;
        forget(b,dump,k,recall);
        //forget(b,dump,k);
#ifdef DEB     
        cout<<"Pk.size() avant Crochetchoix "<<Pk.size()<<endl;
#endif     
        //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,Spoldejavu);
        NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,w,serv);
        //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b);
#ifdef DEB     
        cout<<"Pk.size apres Crochetchoix "<<Pk.size()<<" et k "<<k<<endl;
#endif     
        NewDestroyredspol(redspol,Pk);
        serv.compress(Pk,dump,b,k);
        Dump(Pk,dump);
        L.destroystore();
        //misa jour de B dans crochet choix!!
        continue;
      }
      }
    else
      //la matrice n'est pas de rang plein
      {
    cout<<"matrice de rang pas plein "<<endl;
    typPk secd;
#ifdef BONRAJOUT
    AddB_dim(Mat,L,P,Pk,Mk,dump,b);
#else
    AddB(Mat,Pk,Mk,dump,b,serv);
#endif
    //
    //secondmembre supprime des trucs dans Pk!!
    //
    secd=secondmembre(Mat,Pk);
    rajoute_voisins(Pk,Mk,dump,b,w,serv);
    //ne selectionne que les non nuls
       if(secd.size()==0)
     {
       typPk redspol;
       //for(int i=Pk.size();i>rank(Mat);i--)
       //{
       //  if(Pk.back().size!=0)
       // MAC_REV_FREE<typename typMat::coeff_t>(Pk.rbegin()->nf
       //    ,Pk.rbegin()->size*sizeof(typename typMat::coeff_t));
       //  Pk.pop_back();
       //}//on vire les 0 avant de dumper
       //      AddB(Mat,Pk,Mk,dump,b);
       Dump(Pk,dump);
       ReduceSpol(S,dump,redspol,b,w,serv);//,Spoldejavu);
#ifdef NO_TWICE    
for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
      P.push_back(invconv<polyalp>(*iter));
#endif
       if(redspol.size()==0)
         {
           // forget(b,dump,k);
           //AddB(Mat,Pk,Mk,b);
           //Dump(Pk,dump);
           //ici on modifie Pk dans AddB 
           //<=> forget(b,dump,k) AddB(..) Dump(Pk,dump);
           k++;
           L.destroystore();
           continue;
           //la matrice a pas rang plein second membre =0 
           //ET les Spol de ce niveau se ,reduisent a 0
         }
       else
         {
           //on a des spol qui se reduisent pas a 0
           //cad un defaut de commutation 
           k=mindeg(redspol,serv);
           recover(Pk,Mk,dump,k);
           //Listedejafait(Mk,Spoldejavu);
           forget(b,dump,k,recall);
           //forget(b,dump,k);
           cout<<"il ya un passage par mat rang pas"
               <<"plein et pol non zero"<<endl;
           //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,Spoldejavu);
           NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,w,serv);
           //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b);
           NewDestroyredspol(redspol,Pk);
           serv.compress(Pk,dump,b,k);
           Dump(Pk,dump); 
           L.destroystore();
           continue;
         }
     }
       else
     {
       //on a des pol qui se reduisent pas a 0
       //cad un defaut de commuatation dans les coins
       typPk redspol;
       //for(int i=Pk.size();i>rank(Mat);i--)
       //  Pk.pop_back();
       //cout<<"PK. size "<<Pk.size()<<" Mat.nbcol "<<Mat.ncol<<endl;
       //cout<<"rank(Mat) "<<rank(Mat)<<endl;
       //int j=0;
       //cout<<"avant recover et meme avant le dump "<<endl;
       //for(typename typPk::iterator iter = secd.begin();iter!=secd.end()
       //    ;iter++,j++)
       // cout<<"pol de secd num "<<j<<endl<<invconv<Poly>(*iter)<<endl;
       //      AddB(Mat,Pk,Mk,b);
       //cout<<"apres ADDB le dump "<<endl;
       //for(typename typPk::iterator iter = Pk.begin();iter!=Pk.end()
       //    ;iter++)
       // cout<<invconv<Poly>(*iter)<<endl;
       Dump(Pk,dump);
       //cout<<"secd.size()" <<secd.size()<<endl;
       ReduceSpol(S,dump,redspol,b,w,serv);//,Spoldejavu);
       my_merge(secd,redspol);
#ifdef NO_TWICE    
for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
         P.push_back(invconv<polyalp>(*iter));
#endif
           k=mindeg(secd,serv);
       recover(Pk,Mk,dump,k);
       //cout<<"k "<<k<<endl;
       //Listedejafait(Mk,Spoldejavu);
       //cout<<"apres recover"<<endl;
       //      for(typename typPk::iterator iter = Pk.begin();iter!=Pk.end()
       //    ;iter++)
       //  cout<<invconv<Poly>(*iter)<<endl;
       //cout<<"les Spol non reduits"<<endl;
       //for(typename typPk::iterator iter = secd.begin();iter!=secd.end()
       //    ;iter++)
       //  cout<<invconv<Poly>(*iter)<<endl;
       forget(b,dump,k,recall);
       
       //forget(b,dump,k);
       //cout<<"il ya un passage par mat pas rang plein et Spol non zero"<<endl;
       //for(typename typPk::iterator iter = secd.begin();iter!=secd.end()
       //    ;iter++)
       // cout<<invconv<Poly>(*iter)<<endl;
       //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,Spoldejavu);
       NewCrochetchoix(P,Pk,Mk,secd,k,dump,b,recall,w,serv);
       //NewCrochetchoix(P,Pk,Mk,secd,k,dump,b);
       NewDestroyredspol(secd,Pk);
#ifdef BONRAJOUT
       cout<<"Apres NewCrochetchoix"<<endl;

       for(typename Base::const_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;
         }
         }
#endif
       //for(typename typPk::iterator iter = Pk.begin();iter!=Pk.end()
       //    ;iter++)
       // cout<<invconv<Poly>(*iter)<<endl;
       serv.compress(Pk,dump,b,k);
       Dump(Pk,dump);
       L.destroystore();
       continue;
     }
      
      }
  }
  
  Mat.destroystore();
  //while nouveau mon etc...
  // Pk=recomposePk(dump);
  //return Pk;
  //  return;
#ifdef DEBUG
  for(typename Base::const_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;
    }
    }
#endif
  destroy_space(&w);
  cout<<"sortie d'algo flash et dimension "<<endl;//<<critere(b,k)<<endl;
  return ;//dump;
  //comme ca on conserve le tri par degre des formes normales et 
  //donc la construction des operateurs de multiplication se 
  //fait plus vite.. :-)

}

#endif //ALREADY_corealgo3