synaps/mpol/macaulay.h

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
**********************************************************************/
#ifndef SYNAPS_RESULTANT_MACAULAY_H
#define SYNAPS_RESULTANT_MACAULAY_H
//----------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/linalg/VECTOR.m>
#include <synaps/linalg/rep2d.h>
#include <synaps/linalg/MATRIX.m>
#include <synaps/mpol/MPol.h>
#include <synaps/mpol/matrixof.h>
#include <synaps/mpol/MPOLDST.m>
#include <map>
#include <synaps/arithm/binomial.h>
__BEGIN_NAMESPACE_SYNAPS
//====================================================================
//--------------------------------------------------------------------
unsigned int SizeOfS(unsigned int n, unsigned int k) 
{
  //  double factn = 1;
  //  double factnu=1;
  //  double factnun=1;
  //  for (unsigned int i = 0; i < n;i++) factn *= (n-i);
  //  for (unsigned int i = k+n; i >1;i--) factnun *= i; 
  //  for (unsigned int i = 0; i < k;i++) factnu*=(k-i);
  //   return (int) (factnun/(factn*factnu));
  /* peut eviter un depassement de capacite dans le calcul des factorielles */
  return binomial<int>(n+k,k); 
}
//----------------------------------------------------------------------
template<class POL>
int choose(int k, int n, POL & ML){

  typedef typename POL::monom_t MON;
  
  MON m;
  POL P, Sigma(1);

  for(int i=0;i<n;i++) Sigma +=MON(i,1);
  ML = POL(1);
  for(int j=1; j<=k; j++) {
    
    ML *= Sigma;
  }
  
  int c=0;
  for (typename POL::iterator i = ML.begin(); i!=ML.end(); i++)
    {
      c++;
      i->setcoeff(1.);
    }
  return c;
}

template<class POL>
POL HomChoose(int k, int n){

  typedef typename POL::monom_t MON;
  
  MON m;
  POL ML, P, Sigma =POL(0);

  for(int i=0;i<n;i++) 
    {
      Sigma += POL(MON(i,1));
    }
  ML = POL(1);
  for(int j=0; j<k; j++) 
    ML *= Sigma;
  
  for (typename POL::iterator i = ML.begin(); i!=ML.end(); i++)
    {
      i->setcoeff(1.);
    }
  return ML;
}


//----------------------------------------------------------------------
template<class R, class L>
R macaulay(const L & f, char z='N'){

  typedef typename L::value_type POL;
  typedef typename POL::monom_t MON;
  typedef typename POL::order_t O;

  int d[f.size()], n = f.size()-1, nu = 0;
  POL F;

  // Calcul de mu et des differents degres di
  typename L::const_iterator fi;int i;
  for(i=0, fi=f.begin(); i<=n; i++, fi++)
  {
    d[i] = std::max(MPOLDST::degree(*fi),1);
    nu += d[i];
  }
  nu -= n; 

  MON* xdi = new MON[n+1];
  POL* B   = new POL[n+1];

  xdi[0]= MON(1);
  B[0]  = POL(1);
  for(int j=1; j<=n; j++){
    xdi[j] = MON(j-1,d[j]);
    B[j]=0;
    for(int i=0;i<d[j];i++) {B[j]+=POL(MON(j-1,i));}
    //cout<<"B "<< B[j]<<endl;
  }

  // Construction des bases locales 
  POL* E= new POL[n+1];
  E[0]=1;

  POL U(1), V(1);
  int mu;
    
  for(int i=n;i>0;i--){

    mu = nu-d[i];
    for (typename POL::const_iterator j = E[0].begin(); j!=E[0].end(); j++)
    {
      choose(mu- degree(*j),i,V);
      //XXX A Ameliorer
      E[i] += V*(*j);
    }
    E[0]*=B[i];
    //cout<<"E "<< E[i]<<endl;
  } 
    //cout<<"E "<< E[0]<<endl;
  //Index des monomes
  std::map<MON,int,O> index;
  int c =0, nz=0;
  MON m;
  for(int i=0; i <=n;i++){
    for (typename POL::reverse_iterator j = E[i].rbegin(); j!=E[i].rend();j++){
      m = (*j)*xdi[i];
      index[m]=c;
      //cout <<"index "<<m<<" "<<index[m]<<endl;
      c++;
    }
    nz += E[i].size()*(f[i].size());
    
  }
  
  //Construction de la matrice de Macaulay
  int taille = SizeOfS(n,nu);
  assert(taille==c);

  R Res;
  matrixof::reserve(Res,taille,taille,nz);
  // if(!Res) Error(cout<<"Not enough memory for a matrix"<<i<<"*"<<j<<endl);
  int v=0; 
  POL P;
  
  if(z =='T')
    for(i=0, fi=f.begin(); i <= n; i++, fi++){
      for (typename POL::reverse_iterator j=E[i].rbegin();j!=E[i].rend();j++){
        P = (*fi)*(*j);
        for (typename POL::const_iterator k = P.begin(); k!=P.end(); k++){
          matrixof::assigncoeff(Res,v,index[*k],k->coeff());
        }
        v++;
      }
    }
  else
    for(i=0, fi=f.begin(); i <= n; i++, fi++){
      for (typename POL::reverse_iterator j=E[i].rbegin();j!=E[i].rend();j++){
        P = (*fi)*(*j);
        for (typename POL::const_iterator k = P.begin(); k!=P.end(); k++){
          //cout <<"I    "<<*k<<" "<<index[*k]<<endl;
          matrixof::assigncoeff(Res,index[*k],v,k->coeff());
        }
        v++;
      }
    } 
  return Res;
};

//----------------------------------------------------------------------
template<class R, class L>
R macaulay(const L & f, unsigned int nu, unsigned int nv, char z='T'){

  typedef typename L::value_type POL;
  typedef typename POL::monom_t MON;
  typedef typename POL::order_t O;

  int n=f.size()-1;

  // Construction des bases locales 
  POL* E= new POL[n+1];
  int c=0,taille=0;
  for (typename L::const_iterator j = f.begin(); j !=f.end(); ++j) 
    {
      Choose(nu- Degree(*j),nv,E[c]);
      taille+=E[c].size();c++;
    }
  //XXX A Ameliorer     

  //Index des monomes
  std::map<MON,int,O> index;
  POL S;  Choose(nu,nv,S);
  
  int nbrow =0;
  for (typename POL::reverse_iterator j = S.rbegin(); j!=S.rend(); j++)
    {
      index[*j]=nbrow; nbrow++;
    }
  
  
  //Construction de la matrice de Macaulay
  R Res;
  matrixof::reserve(Res,taille,nbrow,0);
  // if(!Res) Error(cout<<"Not enough memory for a matrix"<<i<<"*"<<j<<endl);
  int v=0;
  POL P;typename L::const_iterator fi=f.begin();
  for(int i=0; i <= n; i++,fi++) {
    for (typename POL::reverse_iterator j = E[i].rbegin(); j!=E[i].rend();j++)
      {
        P = (*fi)*(*j);
        for (typename POL::const_iterator k = P.begin(); k!=P.end(); k++)
          matrixof::assigncoeff(Res,v, index[*k],k->GetCoeff());
        v++;
      }
  }
  return Res;
};
//----------------------------------------------------------------------------
#ifdef Sparse
struct compcol_mat_double;
// Sylvester-like matrix under the specialized form @MatSps<compcol_mat_double> >@
template <class POL>
MatRef<compcol_mat_double>
macaulay(const VectStd<array1d<POL> > & f, compcol_mat_double S) {

  typedef typename POL::monom_t MON;
  compcol_mat_double *Res;
  int nz=0;
  typedef typename compcol_mat_double::coeff_t C;
  C zero = C(0);

  int d[f.size()], n = (f.size()-1), nu = 0;
  POL F;

  // Calcul de mu et des differents degres di
  for(int i=0; i<=n; i++){
    d[i] = MPOLY::degree(f[i]);
    nu += d[i];
   }
  nu -= n; 

  MON* xdi = new MON[n+1];
  POL* D = new POL[n+1];

  xdi[0]= MON(1.0);
  D[0] = POL(xdi[0]);
  for(int j=1; j<=n; j++){
    xdi[j] = MON(j,d[j]);
    for(int i=0;i<d[j];i++) D[j]+= MON(j,i);
  }

  // Construction de l'ensemble des monomes de degre mu
  int taille = SizeOfS(n,nu);
 
  // Construction des bases locales 
  POL* E= new POL[n+1];
  POL U(1.),V(1.);
  int mu=nu;

  for(int i=n;i>0;i--){
    mu = nu-d[i];
    for (typename POL::const_iterator j = U.begin(); j!=U.end(); j++){
      Choose(mu-j->GetDegree(),i,V);
      //XXX A Ameliorer
      E[i] += V*(*j);
    }
    U*=D[i];
  }
  E[0] = U;

  //Index des monomes
  map<MON,int> index;
  int c =0,nbcol=0;
  MON m;
  for(int i=0; i <=n;i++){
    for (typename POL::const_iterator j = E[i].begin(); j!=E[i].end(); j++){
      m = (*j)*xdi[i];
      index[m]=c;
      c++;nbcol++;
    }
    nz+=NbCoef(f[i])*(nbcol); // Calcul du nombre de non-zero
    nbcol=0;
  }

  int * col = new int[taille+1];
  int * r = new int[nz];
  C * val = new C[nz];
  col[0] = 0;

  int u,t=0,l=0, v=0;
  POL P;
  for(int i=0; i <= n; i++){
    for (typename POL::const_iterator j = E[i].begin(); j!=E[i].end(); j++){
      P = f[i]*(*j);
      for (typename POL::const_iterator k = P.begin(); k!=P.end(); k++) {
        if (k->GetCoeff()!=zero) {
          val[t]=k->GetCoeff(); 
          r[t]=index[*k];
          if ((t!=0) && (l!=v)) {
            l++;
            col[l] = t;}
          t++;  
        }
      }
      v++;
    }
  }
  col[taille] = nz;

  Res = new compcol_mat_double(taille, taille, nz, val, r, col, 0);

  return MatRef<compcol_mat_double>(Res);
}

//--------------------------------------------------------------------------
// Transposed Sylvester-like matrix under the specialized form @MatSps<compcol_mat_double> >@
template <class POL>
MatRef<compcol_mat_double>
macaulay_T(const VectStd<array1d<POL> > & f, compcol_mat_double S) {
  typedef typename POL::monom_t MON;
  compcol_mat_double *Res;
  int nz=0;
  typedef typename compcol_mat_double::coeff_t C;
  C zero = C(0);

  int d[f.size()], n = (f.size()-1), nu = 0;
  POL F;

  // Calcul de mu et des differents degres di
  for(int i=0; i<=n; i++){
    d[i] = MPOLY::degree(f[i]);
    nu += d[i];
   }
  nu -= n; 

  MON* xdi = new MON[n+1];
  POL* D = new POL[n+1];

  xdi[0]= MON(1.0);
  D[0] = POL(xdi[0]);
  for(int j=1; j<=n; j++){
    xdi[j] = MON(j,d[j]);
    for(int i=0;i<d[j];i++) D[j]+= MON(j,i);
  }

  // Construction de l'ensemble des monomes de degre mu
  // 
  int taille = SizeOfS(n,nu);


  // Construction des bases locales 
  POL* E= new POL[n+1];
  POL U(1.),V(1.);
  int mu=nu;

  for(int i=n;i>0;i--){
    mu = nu-d[i];
    for (typename POL::const_iterator j = U.begin(); j!=U.end(); j++){
      Choose(mu-j->GetDegree(),i,V);
      //XXX A Ameliorer
      E[i] += V*(*j);
    }
    U*=D[i];

  }
  E[0] = U;

  //Index des monomes
  map<MON,int> index;
  int c =0,nbcol=0;
  MON m;
  for(int i=0; i <=n;i++){
    for (typename POL::const_iterator j = E[i].begin(); j!=E[i].end(); j++){
      m = (*j)*xdi[i];
      index[m]=c;
      c++;nbcol++;
    }
    nz+=NbCoef(f[i])*(nbcol); // Calcul du nombre de non-zero
    nbcol=0;
  }

  int * col = new int[taille+1];
  int * r = new int[nz];
  C * val = new C[nz];
  col[0] = 0;

  int t=0,l=0, v=0;
  POL P;
  for(int i=0; i <= n; i++){
    for (typename POL::const_iterator j = E[i].begin(); j!=E[i].end(); j++){
      P = f[i]*(*j);
      for (typename POL::const_iterator k = P.begin(); k!=P.end(); k++) {
        if (k->GetCoeff()!=zero) {
          val[t]=k->GetCoeff(); 
          r[t]=index[*k];
          if ((t!=0) && (l!=v)) {
            l++;
            col[l] = t;}
          t++;  
        }
      }
      v++;
    }
  }

  col[taille] = nz;

  Res = new compcol_mat_double(taille, taille, nz, val, r, col);
  MatSps<compcol_mat_double>::transpose(*Res);

  return MatRef<compcol_mat_double>(Res);
}

// Sylvester-like matrix under the specialized form @MatSps<umfpack<C> >@
template <class POL, class C>
MatRef<umfpack<C> >
macaulay(const VectStd<array1d<POL> > & f, umfpack<C> S, char tt='N') {
  typedef typename POL::monom_t MON;
  typedef double C;
  umfpack<C> *Res;
  int nz=0;
  C zero = C(0);
  int d[f.size()], n = (f.size()-1), nu = 0;
  POL F;

  // Calcul de mu et des differents degres di
  for(int i=0; i<=n; i++){
    d[i] = MPOLY::degree(f[i]);
    nu += d[i];
   }
  nu -= n; 

  MON* xdi = new MON[n+1];
  POL* D = new POL[n+1];

  xdi[0]= MON(C(1.0));
  D[0] = POL(xdi[0]);
  for(int j=1; j<=n; j++){
    xdi[j] = MON(j,d[j]);
    for(int i=0;i<d[j];i++) D[j]+= MON(j,i);
  }

  // Construction de l'ensemble des monomes de degre mu
  int taille = SizeOfS(n,nu);

  // Construction des bases locales 
  POL* E= new POL[n+1];
  POL U(C(1.)),V(C(1.));
  int mu=nu;

  for(int i=n;i>0;i--){
    mu = nu-d[i];
    for (typename POL::const_iterator j = U.begin(); j!=U.end(); j++){
      Choose(mu-j->GetDegree(),i,V);
      //XXX A Ameliorer
      E[i] += V*(*j);
    }
    U*=D[i];
  }
  E[0] = U;

  //Index des monomes
  map<MON,int> index;
  int c =0,nbcol=0;
  MON m;
  for(int i=0; i <=n;i++){
    for (typename POL::const_iterator j = E[i].begin(); j!=E[i].end(); j++){
      m = (*j)*xdi[i];
      index[m]=c;
      c++;nbcol++;
    }
    nz+=NbCoef(f[i])*(nbcol); // Calcul du nombre de non-zero
    nbcol=0;
  }

  Res = new umfpack<C>(taille, taille, nz);
  int t=0,v=0;
  POL P;
  if(tt =='N')
    for(int i=0; i <= n; i++){
      for (typename POL::const_iterator j = E[i].begin(); j!=E[i].end(); j++){
        P = f[i]*(*j);
        for (typename POL::const_iterator k = P.begin(); k!=P.end(); k++) {
          if (k->GetCoeff()!=zero) {
            Res->tab_[t]=k->GetCoeff(); 
            Res->index_[t]=index[*k];
            Res->index_[t+nz] = v;
            t++;
          }
        }
        v++;
      }
    }
 else
   for(int i=0; i <= n; i++){
     for (typename POL::const_iterator j = E[i].begin(); j!=E[i].end(); j++){
       P = f[i]*(*j);
       for (typename POL::const_iterator k = P.begin(); k!=P.end(); k++) {
         if (k->GetCoeff()!=zero) {
           Res->tab_[t]=k->GetCoeff(); 
           Res->index_[t]=v;
           Res->index_[t+nz] = index[*k];
           t++;
         }
       }
       v++;
     }
   }

  return MatRef<umfpack<C> >(Res);
}
template <class POL, class C>
MatRef<umfpack<C> >
macaulay_T(const VectStd<array1d<POL> > & f, umfpack<C> S) {
  return macaulay(f, umfpack<C>(), 'Y');
}

#endif //Sparse
//-------------------------------------------------------------------------

__END_NAMESPACE_SYNAPS

#endif // SYNAPS_RESULTANT_MACAULAY_H