synaps/mpol/bezoutian.h

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of SYNAPS.               *
*      B. Mourrain, GALAAD, INRIA                                    *
*********************************************************************/
#ifndef SYNAPS_RESULTANT_MBEZOUT_H
#define SYNAPS_RESULTANT_MBEZOUT_H
//----------------------------------------------------------------------
#include <list>
#include <map>
#include <synaps/init.h>
#include <synaps/linalg/MatrDse.h>
#include <synaps/linalg/Det.h>
#include <synaps/mpol/MPol.h>
#include <synaps/mpol/MPOLDST.m>
//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class T>
typename T::value_type Theta(const T & l, typename T::value_type::monom_t::index_t c=0)
{
  typedef typename T::value_type     POL;
  typedef typename POL::monom_t      monom_t;
  typedef typename monom_t::coeff_t  C;
  typedef typename monom_t::index_t  index_t;

  index_t n= l.size()-1;

  std::list<POL> ldp;
  POL* dx = new POL[n];
  for(int i=0; i < n; i++)
  {
        dx[i] = POL(monom_t(C(1),c+i,1))- POL(monom_t(C(1),c+i+n,1));
  }

  MatrDse<POL, linalg::rep2d<POL> > M(n+1,n+1);
  int k=0;

  for(typename T::const_iterator p = l.begin(); p != l.end(); ++p){
        M(k,0)=  (*p);
        POL dp = (*p);
        for(typename POL::const_iterator m =(*p).begin(); m!=(*p).end(); m++){
          index_t s = m->nvars();
          // if the monomial is not zero and the coefficient is not in the ground ring
          if(degree(*m) > 0 && s >= c){
                // t will be the monomial in the new variables and has the same coefficient
                monom_t t(m->coeff(),s+n);
                for(index_t j=0; j < c; j++){
                  // expoments for the variables in the ground ring are the same
                  t.setexpt(j,(*m)[j]);
                  // the variables in y_i = x_{i+n} corresponding to a variable in the ground field has exponents zero
                  t.setexpt(j+n,0);
                }
                for(index_t j=c; j < s;j++) {
                  // exponent of y_i = x_{i+n} in t is the exponent of x_i in *m and the exponent of x_i is zero in t
                  t.setexpt(j,0);
                  t.setexpt(j+n,(*m)[j]);
                }
                //We substract this monomial to the initial polynomial
                dp-=POL(t);
          }else{
                //the monmial must be considered as a constant and removed from the initial polynomial
                dp-=(*m);
          }
        }
        // the divided differences are computed by successive division by x_i - y_i = x_i - y_{i+n}, the coefficient of theta are the successives quotients
        for (int i = 0; i<n;i++)
        {
          MPOLDST::div_rem_x(M(k,i+1), dp, dx[i]);
        }
        k++;
  }
  delete [] dx;
  // The result is the determinant of the matrix
  return Det(M);
}


//-------------------------------------------------------------------------
template <class UPOL, class MPOL>
MatrDse<UPOL> UBezout(const MPOL & PP,const MPOL QQ)
{
  typedef typename MPOL::monom_t     monom_t;
  typedef typename MPOL::order_t     order_t;
  typedef typename monom_t::coeff_t coeff_t;

  std::list<MPOL> L;
  L.push_back(PP);
  L.push_back(QQ);
  MPOL P=Theta(L,1);

  typename MPOL::const_iterator i = P.begin();

  int taille= std::max(degree(PP,1),degree(QQ,1));
  MatrDse<UPOL> res(taille,taille);

  for(typename MPOL::const_iterator i = P.begin(); i != P.end(); ++i)
    {
      int d0=(*i)[0];int d1=(*i)[1];int d2=(*i)[2];
      coeff_t k=i->coeff();
      UPOL P1("x0");
      P1=P1^d0; P1*=k;
      res(d1,d2) += P1;
    }
  return (res);
}
//----------------------------------------------------------------------
template <class T, class I> inline
void add_from_iterator(T & v, const I & it, unsigned l0)
{
  v+=it->coeff();
}
//----------------------------------------------------------------------
template<class C, class R> struct UPolDse;
template<class C,class R,class I> inline
void add_from_iterator(UPolDse<C,R> & v,
                       const I & it,
                       unsigned l0)
{
  v+=UPolDse<C,R>(it->coeff(),(*it)[0]);
}
//----------------------------------------------------------------------
template<class C, class R, class O> struct MPol;
template<class C, class R, class O> inline
void add_from_iterator(MPol<C,R,O> & v,
                       const typename MPol<C,R,O>::const_iterator & it,
                       int l0)
{
  typename MPol<C,R,O>::monom_t cf(it->coeff(),l0,type::AsSize());
  for(int k=0;k<l0;k++) cf.setexpt(k,(*it)[k]);
  v+=cf;
}

//======================================================================
template <class MAT, class POL>
MAT DeltaOf(const POL & P, int l0, int l1, int l2)
{
  typedef typename POL::monom_t     monom_t;
  typedef typename POL::order_t     order_t;
  typedef typename monom_t::coeff_t coeff_t;
  POL p1(0),p2(0);

  monom_t m1(1,l0+l1,type::AsSize()), m2(1,l2,type::AsSize());

  // We construct to polynomials p1 and p2
  // The set of monomials of p1 are the monomials of P on the first set of variables
  // The set of monomials of p2 are the monomials of P on the second set of variables
  for(typename POL::const_iterator i = P.begin(); i != P.end(); ++i){
    for(int k=l0;k<l0+l1;k++){
      m1.setexpt(k,(*i)[k]);
    }
    p1+=POL(m1);
    for(int k=0;k<l2;k++){
      m2.setexpt(k,(*i)[l0+l1+k]);
    }
    p2+=POL(m2);
  }


  unsigned int c1=0, c2=0;
  std::map<monom_t,int,order_t> index1,index2;

  for(typename POL::iterator i = p1.begin(); i != p1.end(); ++i){
    i->setcoeff(coeff_t(1));
    index1[*i]=c1;++c1;
  }
  for(typename POL::iterator i = p2.begin(); i != p2.end(); ++i){
    i->setcoeff(coeff_t(1));
    index2[*i]=c2;++c2;
  }

  MAT res(c1,c2);
  //  int j=-1;
  for(typename POL::const_iterator i = P.begin(); i != P.end(); ++i){
    for(int k=l0;k<l0+l1;k++){
      m1.setexpt(k,(*i)[k]);
    }
    for(int k=0;k<l2;k++){
      m2.setexpt(k,(*i)[l0+l1+k]);
    }
    //    if (index2[m2] !=j) j=index2[m2];

    add_from_iterator(res(index1[m1],index2[m2]),i,l0);

//      monom_t cf(i->coeff(),l0+1,AsSize());
//     for(int k=l0  ;k>=0;--k) cf.setexpt(k,(*i)[k]);
//     (res)(index1[m1],index2[m2])+=POL(cf);
  }
  return (res);
}
//----------------------------------------------------------------------
template <class MAT, class T> inline
MAT mbezout(const T & l, typename T::value_type::monom_t::index_t c=0)
{
  typename T::value_type::monom_t::index_t n= l.size()-1;
  return DeltaOf<MAT>(Theta(l,c),c,n,n);
}


__END_NAMESPACE_SYNAPS

//----------------------------------------------------------------------
#endif // SYNAPS_RESULTANT_MBEZOUT_H