synaps/mpol/approximate.h

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of SYNAPS.               *
*      V.t TRAN , GALAAD, INRIA                                      *
**********************************************************************/
#ifndef SYNAPS_RESULTANT_APPROXIMATE_H
#define SYNAPS_RESULTANT_APPROXIMATE_H
//---------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/linalg.h>
#include <synaps/mpol/MPOLDST.m>
#include <synaps/mpol/matrixof.h>

__BEGIN_NAMESPACE_SYNAPS


template<class POL>
POL MonomV(const int v, const int n)
{
  typedef typename POL::monom_t MON;
  int d1,d2;
  POL p(0);
  for (int j=0;j<=v;j++) {
    POL tmp(MON(0,j));
    p = p+tmp;
  }

  for (int i=1;i<n;i++) {
    POL p1(0), p3(0);
    for (int j=1;j<=v;j++) {
      POL tmp(MON(i,j));
      p1 = p1+tmp;//creation of monome xi^j
    }

    for (typename POL::const_iterator pp=p.begin();pp!=p.end();pp++) {
      POL k1(*pp); d1= MPOLDST::degree(k1);
      POL p2(0);
      for (typename POL::const_iterator tt=p1.begin();tt!=p1.end();tt++){
                POL k2(*tt); d2 = MPOLDST::degree(k2);
                if (d1+d2<=v)//satisfied
                  p2+=k2*k1;    
      }
      p3+=p2;
    }
    p+=p3;
  }
  return p;
};

template<class LP>
LP baseLocal(const LP & f, const std::vector<int> & d, const int nv)
{
  typedef typename LP::value_type POL;
  LP bas;
  int n = f.size();
  for (int i=0;i<n;i++) {
    POL mv = MonomV<POL>(d[i],nv); //number of variables = nv
    for (int j=0;j<nv;j++) //shift the variables
      mv = MPOLDST::swap(mv, j, j+n);
    for (typename POL::const_iterator mm=mv.begin();mm!=mv.end();mm++) {
      bas.push_back(f[i]*(*mm));
    }
  }
  return bas;
};

template <class MAT, class LL>
MAT evaluation(const LL& f,
               const std::vector<int> & d,
               const std::vector<double>& u,
               const int cas=1)
{
  typedef typename LL::value_type POL;
  typedef typename POL::monom_t MON;
  int n=f.size();
  int nv = numvar(f)-n; //number of variables
  LL bas = baseLocal<LL>(f,d,nv);
  int nz = bas.size(), v=0;

  MAT res;
  POL tmp;
  matrixof::reserve(res,nz,nz,1);
  switch (cas){
  case 1://e(U) simple
        for (int i=0;i<nz;i++) {
          v = 0;
          for (typename LL::const_iterator pp = bas.begin();pp!=bas.end();pp++) {
                POL tmp = (*pp);
                for (int j=n;j<n+nv;j++)
                  tmp = MPOLDST::subs(tmp,j,u[i*nv+j-n]);//replace variables by values
                matrixof::assigncoeff(res,i,v++,tmp);
          }
        }
        break;
  case 2://e(U) & e(U')
        int vu=0;
        int i=0;
        while (i<nz) {
          v = 0;
          for (typename LL::const_iterator pp = bas.begin();pp!=bas.end();pp++) {
                tmp = (*pp);    //evaluation
                for (int j=n;j<n+nv;j++)
                  tmp = MPOLDST::subs(tmp,j,u[vu*nv+j-n]);
                matrixof::assigncoeff(res,i,v++,tmp);
          }
          i++; int k=0;

          while (i<nz && k<nv) {
                v=0;
                for (typename LL::const_iterator pp = bas.begin();pp!=bas.end();pp++) {
                  POL tmp1 = (*pp);
                  MPOLDST::diff(tmp,tmp1,n+k);//derivation with n+k
                  for (int j=n;j<n+nv;j++)
                        tmp = MPOLDST::subs(tmp,j,u[vu*nv+j-n]);
                  matrixof::assigncoeff(res,i,v++,tmp);
                }
                i++; k++;
          }
          vu++;
        }
        break;
  }

  return res;
};

template<class LL>
typename LL::value_type Elim_Approximate(const LL& f,
                                         const std::vector<int> & d,
                                         const std::vector<double> & u,
                                         const int cas=1)
{
  typedef typename LL::value_type POL;
  typedef MatrDse<POL, linalg::rep2d<POL> > MAT;

  MAT res = evaluation<MAT, LL>(f,d,u,cas);
#ifdef DEBUG
  std::cout<<"Matrix E :\n"<<res<<std::endl;
#endif
  POL D = Det(/*Berko(),*/res);
  return D;
};


__END_NAMESPACE_SYNAPS

//---------------------------------------------------------------------
#endif // SYNAPS_RESULTANT_APPROXIMATE_H