symbo2.hpp

Go to the documentation of this file.

#ifdef COEFF
#define tutu(a) a
#define toto(a) tutu(a)
#define OLD_COEFF toto(COEFF)
//#undef COEFF
#endif
OLD_COEFF machin;
#define COEFF mpq_t 
#define PREFIX mpq
#include"ptpol.c"

template<>
mpq_t  &harewell<mpq_t>::operator() (const int &i,const int &j,const mpq_t &d)
{
  cout<<"Nein Nein!!!!!!!!!!!"<<endl;
  static mpq_t toto;
  return toto;
}

template<>
ptpol_t  &harewell<ptpol_t>::operator() (const int &i,const int &j,const ptpol_t &d)
{
  cout<<"Nein Nein!!!!!!!!!!!"<<endl;
  static ptpol_t toto;
  return toto;
}

template<>
mpq_t & harewell<mpq_t>::operator() (const int  &i,const int &j) 
{
  static mpq_t d;
  cout<<"Nein!!!!!!!!!!!"<<endl;
  //cout<<*(int*)NULL<<endl;
  int colbegin=colptr[j]
    ,colend=(j+1<this->ncol)?colptr[j+1]:nnz;
  int t;
  //for(t=colbegin;(t<colend)&&(i>rowind[t]);t++) ;
  //le 17/07/2001 idee pouvoir permutter les ligne plus facilement 
  for(t=colbegin;(t<colend)&&(i!=rowind[t]);t++);
  if (t>=colend) return d;
  return ((i==(rowind[t]))?(((mpq_t*)(nzval))[t]):d);
}

template<>
ptpol_t & harewell<ptpol_t>::operator() (const int  &i,const int &j) 
{
  static ptpol_t d;
  cout<<"Nein!!!!!!!!!!!"<<endl;
  //  cout<<*(int*)NULL<<endl;
  int colbegin=colptr[j]
    ,colend=(j+1<this->ncol)?colptr[j+1]:nnz;
  int t;
  //for(t=colbegin;(t<colend)&&(i>rowind[t]);t++) ;
  //le 17/07/2001 idee pouvoir permutter les ligne plus facilement 
  for(t=colbegin;(t<colend)&&(i!=rowind[t]);t++);
  if (t>=colend) return d;
  return ((i==(rowind[t]))?(((ptpol_t*)(nzval))[t]):d);
}

template<> 
harewell<mpq_t>::harewell(const int nrow,const int ncol)
{
  //Pour eviter le undefined behavior on met tous lespointeurs non initialises a NULL
  this->nrow=nrow;
  this->ncol=ncol;
  this->rowind=(int*)malloc(sizeof(int)*(this->ncol));
  this->colptr=(int*)malloc(sizeof(int)*(1+this->ncol));
  this->nzval=malloc(sizeof(mpq_t)*(this->ncol));
  this->nnz=this->ncol;
  for(int i=0;i<this->ncol;i++)
    {
      this->rowind[i]=i;
      this->colptr[i]=i;
      mpq_init(((mpq_t*)this->nzval)[i]);
    }
  this->colptr[this->ncol]=this->nnz;
}


template<> 
harewell<ptpol_t>::harewell(const int nrow,const int ncol)
{
  //Pour eviter le undefined behavior on met tous lespointeurs non initialises a NULL
  this->nrow=nrow;
  this->ncol=ncol;
  this->rowind=(int*)malloc(sizeof(int)*(this->ncol));
  this->colptr=(int*)malloc(sizeof(int)*(1+this->ncol));
  this->nzval=malloc(sizeof(ptpol_t)*(this->ncol));
  this->nnz=this->ncol;
  for(int i=0;i<this->ncol;i++)
    {
      this->rowind[i]=i;
      this->colptr[i]=i;
      ptpol_init(((ptpol_t*)this->nzval)[i]);
    }
  this->colptr[this->ncol]=this->nnz;
}

void harewell_set_coeff(harewell<mpq_t> &mat,int i,int j, mpq_t coeff)
{
  mpq_t zero;
  mpq_init(zero);
  int colbegin=mat.colptr[j],colend=(j+1<mat.ncol)?mat.colptr[j+1]:mat.nnz;
  int t,carry=0;;
  for(t=colbegin;(t<colend)&&(i!=(mat.rowind[t]));t++);
  if (t>=colend) {t--;carry=1;};
  if (mat.rowind[t]==i) // je suis tombe sur un coeff qui etait deja pas nul
    {
      if (mpq_cmp(coeff,zero))//super on ne vire rien on change un valeur!!!!
    {
      mpq_set(((mpq_t *)(mat.nzval))[t],coeff);
    }
      else//aie on vire un truc
    {
      if ((t==colbegin)&&(colend==colbegin+1))
        mpq_set(((mpq_t *)(mat.nzval))[t],coeff);//on ne cree pas de colonnes de zeros!!!
      else
        {
          //on ne cree jamais de colonnes de 0
          //le formatNC ne les gerent pas!!!
          memmove(((mpq_t*)(mat.nzval)+t)
              ,((mpq_t*)(mat.nzval))+t+1,(mat.nnz-t-1)*sizeof(mpq_t));
          memmove(((int*)(mat.rowind)+t),((int*)(mat.rowind))+t+1
              ,(mat.nnz-t-1)*sizeof(int));
          mat.nnz--;
          //reglage des colptr
          for(int k=j+1;k<mat.ncol;k++)
        mat.colptr[k]--;
          //reglage des nnz
          mat.colptr[mat.ncol]=mat.nnz;
        }
    }
    }
  else//il faut creer la place pour le nouveau coeff
    // il n'y a pas l'equivalent realloc en C++ !!!!!
    if (!mpq_cmp(coeff,zero)) {}
    else
      {
    mat.nnz++;
    t+=carry;
    //
    //les champs dynamiques sont supposes alloues par malloc
    //et pas par new!!!!!! c pour connexion avec super lu 
    // qui bricole sa memoire avec malloc
    //
    mpq_t* tmpnzval=(mpq_t*)malloc(sizeof(mpq_t)*(mat.nnz));    
    int *tmprowind=(int*)malloc(sizeof(int)*mat.nnz);
    memcpy(((mpq_t*)(tmpnzval)),((mpq_t*)(mat.nzval)),(t)*sizeof(mpq_t));
    memcpy(tmprowind,(mat.rowind),(t)*sizeof(int));
    mpq_init(tmpnzval[t]);
    mpq_set(tmpnzval[t],coeff);
    tmprowind[t]=i;
    if (t==colbegin) mat.colptr[j]=t;//on met a jour pour la nvlle valeur
    for (int i=j+1;i<mat.ncol;i++)
      mat.colptr[i]++;//il faut decaler les indices de 1
    memcpy(((mpq_t*)(tmpnzval))+t+1,
           ((mpq_t*)(mat.nzval))+t,(mat.nnz-t-1)*sizeof(mpq_t));
    memcpy(tmprowind+t+1,mat.rowind+t,(mat.nnz-t-1)*sizeof(int));
    //for (k=t;k<nnz-1;k++)
    //{
    //  tmpnzval[k+1]=((T*)(nzval))[k];
    //  tmprowind[k+1]=rowind[k];
    //};
    free(mat.nzval);
    mat.nzval=(void*)tmpnzval;
    free(mat.rowind);
    mat.rowind=tmprowind;
    mat.colptr[mat.ncol]=mat.nnz;//le colptr[ncol] c nnz recopie pour 
    // usage interne a la lib apres.
      }
}

void harewell_set_coeff(harewell<ptpol_t> &mat,int i,int j, ptpol_t coeff)
{
  int colbegin=mat.colptr[j],colend=(j+1<mat.ncol)?mat.colptr[j+1]:mat.nnz;
  int t,carry=0;;
  for(t=colbegin;(t<colend)&&(i!=(mat.rowind[t]));t++);
  if (t>=colend) {t--;carry=1;};
  if (mat.rowind[t]==i) // je suis tombe sur un coeff qui etait deja pas nul
    {
      //cout<<"je passe par la"<<endl;
      //ptpol_out_str(stdout,coeff);
      //cout<<endl<<"la comparaison "<<ptpol_cmp_ui(coeff,0)<<endl;
      if (ptpol_cmp_ui(coeff,0))//super on ne vire rien on change un valeur!!!!
    {
      //cout<<"je fait un set"<<endl;
      ptpol_set(((ptpol_t *)(mat.nzval))[t],coeff);
    }
      else//aie on vire un truc
    {
      if ((t==colbegin)&&(colend==colbegin+1))
        ptpol_set(((ptpol_t *)(mat.nzval))[t],coeff);//on ne cree pas de colonnes de zeros!!!
      else
        {
          //on ne cree jamais de colonnes de 0
          //le formatNC ne les gerent pas!!!
          memmove(((ptpol_t*)(mat.nzval)+t)
              ,((ptpol_t*)(mat.nzval))+t+1,(mat.nnz-t-1)*sizeof(ptpol_t));
          memmove(((int*)(mat.rowind)+t),((int*)(mat.rowind))+t+1
              ,(mat.nnz-t-1)*sizeof(int));
          mat.nnz--;
          //reglage des colptr
          for(int k=j+1;k<mat.ncol;k++)
        mat.colptr[k]--;
          //reglage des nnz
          mat.colptr[mat.ncol]=mat.nnz;
        }
    }
    }
  else//il faut creer la place pour le nouveau coeff
    // il n'y a pas l'equivalent realloc en C++ !!!!!
    if (!ptpol_cmp_ui(coeff,0)) {}
    else
      {
    mat.nnz++;
    t+=carry;
    //
    //les champs dynamiques sont supposes alloues par malloc
    //et pas par new!!!!!! c pour connexion avec super lu 
    // qui bricole sa memoire avec malloc
    //
    ptpol_t* tmpnzval=(ptpol_t*)malloc(sizeof(ptpol_t)*(mat.nnz));
    int *tmprowind=(int*)malloc(sizeof(int)*mat.nnz);
    memcpy(((ptpol_t*)(tmpnzval)),((ptpol_t*)(mat.nzval)),(t)*sizeof(ptpol_t));
    memcpy(tmprowind,(mat.rowind),(t)*sizeof(int));
    ptpol_init(tmpnzval[t]);
    ptpol_set(tmpnzval[t],coeff);
    tmprowind[t]=i;
    if (t==colbegin) mat.colptr[j]=t;//on met a jour pour la nvlle valeur
    for (int i=j+1;i<mat.ncol;i++)
      mat.colptr[i]++;//il faut decaler les indices de 1
    memcpy(((ptpol_t*)(tmpnzval))+t+1,
           ((ptpol_t*)(mat.nzval))+t,(mat.nnz-t-1)*sizeof(ptpol_t));
    memcpy(tmprowind+t+1,mat.rowind+t,(mat.nnz-t-1)*sizeof(int));
    //for (k=t;k<nnz-1;k++)
    //{
    //  tmpnzval[k+1]=((T*)(nzval))[k];
    //  tmprowind[k+1]=rowind[k];
    //};
    free(mat.nzval);
    mat.nzval=(void*)tmpnzval;
    free(mat.rowind);
    mat.rowind=tmprowind;
    mat.colptr[mat.ncol]=mat.nnz;//le colptr[ncol] c nnz recopie pour 
    // usage interne a la lib apres.
      }
}

void harwell_get_coeff(ptpol_t & res,const harewell<ptpol_t> & mat,
               int i, int j)
{
  int colbegin=mat.colptr[j]
    ,colend=(j+1<mat.ncol)?mat.colptr[j+1]:mat.nnz;
  int t;
  //  cout<<"colbeg "<<colbegin<<endl;
  //cout<<"colend "<<colend<<endl;
  for(t=colbegin;
      (t<colend)
    &&(i!=mat.rowind[t]);t++);
    //cout<<mat.rowind[t]<<endl;
  if (t>=colend) {ptpol_set_d(res,0);return;}
  if (i==(mat.rowind[t]))
    ptpol_set(res,(((ptpol_t*)(mat.nzval))[t]));
  else
    ptpol_set_d(res,0);
}

void harwell_get_coeff(mpq_t & res,const harewell<mpq_t> & mat,
               int i, int j)
{
  int colbegin=mat.colptr[j]
    ,colend=(j+1<mat.ncol)?mat.colptr[j+1]:mat.nnz;
  int t;
  for(t=colbegin;(t<colend)&&(i!=mat.rowind[t]);t++);
  if (t>=colend) 
    {
      mpq_set_d(res,0);
      return;
    }
  
  if (i==(mat.rowind[t]))
    mpq_set(res,(((mpq_t*)(mat.nzval))[t]));
  else
    mpq_set_d(res,0);
}
void harwell_get_coeff(__mpq_struct *res,const harewell<__mpq_struct*> & mat,
               int i, int j)
{
  int colbegin=mat.colptr[j]
    ,colend=(j+1<mat.ncol)?mat.colptr[j+1]:mat.nnz;
  int t;
  for(t=colbegin;(t<colend)&&(i!=mat.rowind[t]);t++);
  if (t>=colend) mpq_set_d(res,0);
  if (i==(mat.rowind[t]))
    mpq_set(res,(((mpq_t*)(mat.nzval))[t]));
  else
    mpq_set_d(res,0);
}

int Iszerop(const ptpol_t c)
{
  return ptpol_isNull(c);
}

void swapcol(harewell<ptpol_t> &m,int k, int l)
{
  ptpol_t tmp,tmp2;
  ptpol_init(tmp);  ptpol_init(tmp2);
  for(int i=0;i<m.nrow;i++)
    {
      harwell_get_coeff(tmp,m,i,k);
      //cout<<"coeff lu"<<endl;
      //ptpol_out_str(stdout,tmp);
      //cout<<endl;
      harwell_get_coeff(tmp2,m,i,l);    
      //cout<<"coeff lu2"<<endl;
      //ptpol_out_str(stdout,tmp2);
      //cout<<endl;
      harewell_set_coeff(m,i,k,tmp2);
      //cout<<"-------------"<<endl;
      harewell_set_coeff(m,i,l,tmp);
    
      harwell_get_coeff(tmp,m,i,k);    
      //cout<<"coeff relu"<<endl;
      //ptpol_out_str(stdout,tmp);
      //cout<<endl;
      //harwell_get_coeff(tmp2,m,i,l);    
      //cout<<"coeff relu2"<<endl;
      //ptpol_out_str(stdout,tmp2);
      //cout<<endl;
     
}
  ptpol_clear(tmp);  ptpol_clear(tmp2);
}

void swaprow(harewell<ptpol_t> &m,int k, int l)
{
  ptpol_t tmp,tmp2;
  ptpol_init(tmp);  ptpol_init(tmp2);
  for(int i=0;i<m.ncol;i++)
    {
      harwell_get_coeff(tmp,m,k,i);
      harwell_get_coeff(tmp2,m,l,i);    
      harewell_set_coeff(m,k,i,tmp2);
      harewell_set_coeff(m,l,i,tmp);
    }
  ptpol_clear(tmp);  ptpol_clear(tmp2);
}

template <typename typmat>
void Triang(typmat & M)
{
  ptpol_t d,un,tmp,tmp2;
  ptpol_init(d);
  ptpol_init(tmp);  
  ptpol_init(tmp2);
  ptpol_init(un);
  ptpol_set_d(d,1);
  ptpol_set_d(un,1);
  int l,tmpl;
  for(int k=0; k<M.nrow;k++)
    {
      //cout<<"la matrice ici"<<endl;
      //printmatpol(M);
      ptpol_set_d(tmp,0);
      ptpol_set_d(tmp2,0);
      for(l=k;(l<M.nrow);l++)
      //for(l=k;(l>=0) && (Iszerop(tmp));l--)
    {
      harwell_get_coeff(tmp,M,l,k);
      if (Iszerop(tmp2) 
          || ((!Iszerop(tmp)) && (tmp->degree<tmp2->degree)))
        {
          tmpl=l+1;
          ptpol_set(tmp2,tmp);
        }
      //cout<<"valeur de tmp"<<endl;
      //ptpol_out_str(stdout,tmp);
      //cout<<" "<<Iszerop(tmp)<<endl;
      
    }
      //l++;
      l=tmpl;
      ptpol_set(tmp,tmp2);
      l--;
      cout<<"mon pivot "<<l<<" et "<<k<<endl;
      ptpol_out_str(stdout,tmp);cout<<endl;
      swaprow(M,k,l);
      //cout<<"la matrice apres swap"<<endl;
      //printmatpol(M);
      //ptpol_out_str(stdout,tmp);
      //cout<<endl;

      if(l!=M.nrow){
      //if(l!=-1){
    if(k>0) harwell_get_coeff(d,M,k-1,k-1); //d=M(k-1,k-1); 
    //if(k<M.col-1) harwell_get_coeff(d,M,k-1,k-1); //d=M(k-1,k-1); 
    for(int i=k+1; i<M.nrow; i++)
      //for(int i=k-1; i>=0; i--)
      {
        ptpol_t c;
        ptpol_init(c);
        harwell_get_coeff(c,M,i,k);
        for(int j=k; j<M.ncol; j++) 
        //for(int j=k; j<M.nrow; j++) 
          {
        //cout<<"i "<<i<<" j "<<j<<endl;
        harwell_get_coeff(tmp,M,i,j);
        harwell_get_coeff(tmp2,M,k,k);
        ptpol_karatsuba_mul(tmp2,tmp,tmp2);
        harwell_get_coeff(tmp,M,k,j);
        ptpol_karatsuba_mul(tmp,tmp,c);
        ptpol_sub(tmp2,tmp2,tmp);
        
              
        //M(i,j)-=c*M(k,j);
        if (!Iszerop(tmp2) && (!Iszerop(d)) && (ptpol_cmp(d, un)))  
          ptpol_div(tmp2,tmp2,d);
        harewell_set_coeff(M,i,j,tmp2);
          }
        ptpol_clear(c);

        //M(i,k)=C(0);   
      }
      }
    }
  ptpol_clear(tmp);     
  ptpol_clear(tmp2);
  ptpol_clear(d);
  ptpol_clear(un);
}

template <typename typmat>
void Triangsave(typmat & M)
{
  ptpol_t d,un,tmp,tmp2;
  ptpol_init(d);
  ptpol_init(tmp);  
  ptpol_init(tmp2);
  ptpol_init(un);
  ptpol_set_d(d,1);
  ptpol_set_d(un,1);
  int l;
  for(int k=0; k<M.ncol;k++)
  //for(int k=M.ncol-1; k>=0;k--)
    {
      //cout<<"la matrice ici"<<endl;
      //printmatpol(M);
      ptpol_set_d(tmp,0);
      for(l=k;(l<M.ncol) && (Iszerop(tmp));l++)
      //for(l=k;(l>=0) && (Iszerop(tmp));l--)
    {
      harwell_get_coeff(tmp,M,k,l);
      //cout<<"valeur de tmp"<<endl;
      //ptpol_out_str(stdout,tmp);
      //cout<<" "<<Iszerop(tmp)<<endl;
      
    }
      //l++;
      l--;
      //cout<<"mon pivot "<<l<<" et "<<k<<endl;
      swapcol(M,k,l);
      //cout<<"la matrice apres swap"<<endl;
      //printmatpol(M);

      ptpol_out_str(stdout,tmp);
      cout<<endl;
      if(l!=M.ncol){
      //if(l!=-1){
    if(k>0) harwell_get_coeff(d,M,k-1,k-1); //d=M(k-1,k-1); 
    //if(k<M.col-1) harwell_get_coeff(d,M,k-1,k-1); //d=M(k-1,k-1); 
    for(int i=k+1; i<M.ncol; i++)
      //for(int i=k-1; i>=0; i--)
      {
        ptpol_t c;
        ptpol_init(c);
        harwell_get_coeff(c,M,k,i);
        for(int j=k; j<M.nrow; j++) 
        //for(int j=k; j<M.nrow; j++) 
          {
        //cout<<"i "<<i<<" j "<<j<<endl;
        harwell_get_coeff(tmp,M,j,i);
        harwell_get_coeff(tmp2,M,k,k);
        ptpol_karatsuba_mul(tmp2,tmp,tmp2);
        harwell_get_coeff(tmp,M,j,k);
        ptpol_karatsuba_mul(tmp,tmp,c);
        ptpol_sub(tmp2,tmp2,tmp);
        
              
        //M(i,j)-=c*M(k,j);
        if (!Iszerop(tmp2) && (!Iszerop(d)) && (ptpol_cmp(d, un)))  
          ptpol_div(tmp2,tmp2,d);
        harewell_set_coeff(M,j,i,tmp2);
          }
        ptpol_clear(c);

        //M(i,k)=C(0);   
      }
      }
    }
  ptpol_clear(tmp);     
  ptpol_clear(tmp2);
  ptpol_clear(d);
  ptpol_clear(un);
}

template<typename typmat, typename typsec>
void SolveUff_restricted(const typmat &mat, typsec &sec)
{
  ptpol_t accu;
  ptpol_init(accu);
  ptpol_set_d(accu,1);
  cout<<mat.ncol<<" "<<mat.nrow<<endl;
  for(int i=mat.nrow-2;i>=0;i--)
    {
      int j;
      
      cout<<"je tourne la "<<i<<endl;
      //la matrice est permutee par rapport a la forme 
      //triangulaire sup
      ptpol_t tmp;
      ptpol_init(tmp);
      for(j=0;(j<(mat.ncol-1)) && Iszerop(tmp);j++)
    {
      harwell_get_coeff(tmp,mat,i,j);
      //cout<<"valeur de tmp"<<endl;
      //ptpol_out_str(stdout,tmp);
      //cout<<" "<<Iszerop(tmp)<<endl;
    }
      j--;
      //cout<<"j vaut "<<j<<" "<<mat.ncol<<endl;
      if(j==mat.ncol-1)
    {
      cout<<"Merde a la ligne "<<i<<endl;
      exit(10);
    }
      for(int k=0;k<mat.nrow-1;k++)
    {
      //cout<<"ici"<<endl;
      //harwell_get_coeff(tmp,mat,k,j);
      if(k!=i)
        ptpol_karatsuba_mul(sec[k],sec[k],tmp);
#if 0   
      ptpol_div_r(accu,sec[k],sec[mat.nrow-1]);
      cout<<endl<<"valeur de accu"<<endl;
      ptpol_out_str(stdout,accu);
      cout<<" "<<endl;
      cout<<"valeur de scxe[k]"<<endl;
      ptpol_out_str(stdout,sec[k]);
      cout<<" "<<endl;
      cout<<"valeur de charpoly"<<endl;
      ptpol_out_str(stdout,sec[mat.nrow-1]);
      cout<<" "<<endl;
#endif  
      ptpol_div_r(sec[k],sec[k],sec[mat.nrow-1]);
      //cout<<"apres mul"<<endl;
      //if(!Iszerop(accu))
      //  ptpol_div(sec[k],sec[k],accu);
    }
      //harwell_get_coeff(accu,mat,i,j);
      for(int k=i-1;k>=0;k--)
    {
      //  cout<<"la"<<endl;
      harwell_get_coeff(tmp,mat,k,j);
      ptpol_karatsuba_mul(tmp,sec[i],tmp);
      //ptpol_div_r(sec[k],sec[k],sec[mat.nrow-1]);
      ptpol_sub(sec[k],sec[k],tmp);//sec[k]-=sec[i]*mat(k,j);
    }
      ptpol_clear(tmp);
    }
}

template<typename typmat, typename typsec>
void SolveUff_restrictedsave(const typmat &mat, typsec &sec)
{
  ptpol_t accu;
  ptpol_init(accu);
  ptpol_set_d(accu,1);
  for(int i=1;i<mat.nrow;i++)
    {
      int j;
      
      cout<<"je tourne la"<<endl;
      //la matrice est permutee par rapport a la forme 
      //triangulaire sup
      ptpol_t tmp;
      ptpol_init(tmp);
      for(j=mat.ncol-1;j>=0 && Iszerop(tmp);j--)
    {
      harwell_get_coeff(tmp,mat,i,j);
      //cout<<"valeur de tmp"<<endl;
      //  ptpol_out_str(stdout,tmp);
      //  cout<<" "<<Iszerop(tmp)<<endl;
    }
      if(j==-1)
    {
      cout<<"Merde a la ligne "<<i<<endl;
      exit(10);
    }
      for(int k=i;k<mat.nrow;k++)
    {
      //cout<<"ici"<<endl;
      //harwell_get_coeff(tmp,mat,k,j);
      ptpol_karatsuba_mul(sec[k],sec[k],tmp);
      //cout<<"apres mul"<<endl;
      if(!Iszerop(accu))
        ptpol_div(sec[k],sec[k],accu);
    }
      harwell_get_coeff(accu,mat,i,j);
      for(int k=i+1;k<mat.nrow;k++)
    {
      //  cout<<"la"<<endl;
      harwell_get_coeff(tmp,mat,k,j);
      ptpol_karatsuba_mul(tmp,sec[i-1],tmp);
      ptpol_sub(sec[k-1],sec[k-1],tmp);//sec[k]-=sec[i]*mat(k,j);

    }
      ptpol_clear(tmp);
    }
}

void printmatpol(harewell<ptpol_t> &m)
{
  int i;
  ptpol_t tmp;
  ptpol_init(tmp);
  cout<<"nrow "<<m.nrow<<" ncol "<<m.ncol<<endl;
  cout<<"[";
  for(i=0;i<m.nrow-1;i++)
    {
      int j;
      cout<<"[";
      for(j=0;j<m.ncol-1;j++)
    {
      harwell_get_coeff(tmp,m,i,j);   ptpol_out_str(stdout,tmp);
      cout<<",";
    }
      harwell_get_coeff(tmp,m,i,j);
      ptpol_out_str(stdout,tmp);
      cout<<"],"<<endl;
    }
  {
    int j;
    cout<<"[";
    for(j=0;j<m.ncol-1;j++)
      {
    harwell_get_coeff(tmp,m,i,j);
    ptpol_out_str(stdout,tmp);
    cout<<",";
      }
    harwell_get_coeff(tmp,m,i,j);
    ptpol_out_str(stdout,tmp);
    cout<<"]"<<endl;
  }
  cout<<"]"<<endl;
  ptpol_clear(tmp);
}

template<typename typmat, typename typres, typename typmon
     , typename typdump, typename Base>
void symbo(typres &res,  const list<typmat> &lmat, const typmon &stockmon
       ,const typdump &dump, const Base &b)
{
  typedef typename typdump::value_type::pol_t mpoly;
  typedef typename mpoly::coeff_t coeff;
  ptpol_t pol_x;
  ptpol_t *secmembre;
  int stockcol;
  int coldim=lmat.front().nrow,rowdim=lmat.front().ncol;
  harewell<ptpol_t> combilin(coldim,rowdim),Xid(coldim,rowdim);
  ptpol_init(pol_x);
  ptpol_set_d(pol_x,1);
  ptpol_mulvar(pol_x,pol_x,1);
  cout<<"debut symbo"<<endl;
  for(typename list<typmat>::const_iterator itercomb=lmat.begin()
    ;itercomb!=lmat.end();itercomb++)
    {
      COEFF tmpcoeff,tmpcoeff2;
      ptpol_t tmp,tmp2;
      gen_init(tmpcoeff); gen_init(tmpcoeff2);
      gen_set_d(tmpcoeff2,rand()%32051);
      ptpol_init(tmp);ptpol_init(tmp2);
      for(int i=0;i<combilin.nrow;i++)
    for(int j=0;j<combilin.ncol;j++)
      {
        harwell_get_coeff(tmpcoeff,*itercomb,i,j);
        gen_mul(tmpcoeff,tmpcoeff,tmpcoeff2);

        harwell_get_coeff(tmp,combilin,i,j);
        ptpol_set_coeff(tmp2,tmpcoeff);

        ptpol_add(tmp,tmp,tmp2);
        harewell_set_coeff(combilin,i,j,tmp);
        //    combilin(i,j,combilin(i,j)
        //   +UPolyDense<upar<coeff> >(coeff(rand()%13)
        //               *((*itercomb)(i,j)),0));
      }
      gen_clear(tmpcoeff);gen_clear(tmpcoeff2);
      ptpol_clear(tmp);
      ptpol_clear(tmp2);
    }
  for(int i=0;i<Xid.nrow;i++)
    {
      ptpol_t tmp;
      ptpol_init(tmp);
      harwell_get_coeff(tmp,combilin,i,i);
      ptpol_sub(tmp,tmp,pol_x);
      harewell_set_coeff(combilin,i,i,tmp);
      ptpol_clear(tmp);
    }
  //cout<<combilin.nnz<<endl;
  //  for(int i=0;i<combilin.ncol;i++)
  //cout<<"colptr "<<i<<" "<<combilin.colptr[i]<<endl;
  cout<<"combilin"<<endl;
  printmatpol(combilin);
  Triang(combilin);
  cout<<"apres triang"<<endl;
  printmatpol(combilin);
  //for(int i=0;i<combilin.ncol;i++)
  //
  //  for(int j=i;j<combilin.ncol;j++)
  //{
  //  harwell_get_coeff(pol_x,combilin,i,j);
  //  ptpol_out_str(stdout,pol_x);
  //  printf("\n");
  //    }
  harwell_get_coeff(res[b.nbvar()],combilin,combilin.nrow-1,combilin.ncol-1);
  secmembre=(ptpol_t*)malloc((combilin.nrow)*sizeof(ptpol_t));
  for(int i=0;i<combilin.nrow;i++)
    ptpol_init(secmembre[i]);
  for(stockcol=0;Iszerop(secmembre[0]);stockcol++)
    harwell_get_coeff(secmembre[0],combilin,combilin.nrow-1,stockcol);
  stockcol--;
  cout<<"stockcol "<<stockcol<<endl;
  for(int i=0;i<combilin.nrow;i++)
    {
      harwell_get_coeff(secmembre[i],combilin,i,stockcol);
    }
#if 0
  printf("\navant------------------\n");
  for(int j=0;j<combilin.nrow-1;j++)
    {
      ptpol_out_str(stdout,secmembre[j]);
      printf("\n------------------\n");
    }
#endif
  //combilin.ncol--;  
  //combilin.nrow--;
  SolveUff_restricted(combilin,secmembre);
#if 0
  printf("\napres------------------\n");
  for(int j=0;j<combilin.nrow-1;j++)
    {
      ptpol_out_str(stdout,secmembre[j]);
      printf("\n------------------\n");
    }
#endif
  ptpol_out_str(stdout,secmembre[combilin.nrow-1]);
  printf("\n------------------\n");
 for(int i=0;i<b.nbvar();i++)
    {
      mpoly tmp;
      tmp.size=-1;
      int j=0;
      findcorresppol(mon(i,1),dump,tmp);
      for(typename typmon::const_iterator iter=stockmon.begin()
        ;*iter!=mon(i,1) && iter!=stockmon.end();j++,iter++);
      if (tmp.size == -1)
    {
      *res[i]=*secmembre[j];
      cout<<"je trouve x"<<i<<" a la position "<<j<<endl;
    }
      else
    {
      cout<<"je ne trouve pas x"<<i<<" a la position "<<j<<endl;
      *res[i]=*pol_x;
    }
    }
 *res[b.nbvar()]=*secmembre[0];
 ptpol_out_str(stdout, res[b.nbvar()]);
   printf("ttttttttttttttttttttt\n");
}
#ifdef OLD_COEFF

//#undef COEFF
#define COEFF toto(OLD_COEFF)
#endif