/Users/mourrain/Devel/mmx/realroot/include/realroot/sparse_monomials.hpp

Go to the documentation of this file.

#ifndef  realroot_sparse_monomials_hpp
#define  realroot_sparse_monomials_hpp
# include <list>
# include <vector>
# include <algorithm>
# include <realroot/binomials.hpp>
# include <realroot/Seq.hpp>
# include <realroot/monomial.hpp>
# include <realroot/monomial_ordering.hpp>
# define TMPL  template<class C, class O, class MONOM, class REP>
# define TMPLX template<class C, class O, class MONOM, class REP, class X>
# define Polynomial monomial_seq<C,O,MONOM,REP>
# define CLASS      monomial_seq
# define sparse_monomials  sparse::monomial_seq
//======================================================================
namespace mmx {


//--------------------------------------------------------------------

namespace sparse {

  template<class C, class O = DegRevLex, 
           class MONOM=monom<C>, 
           class REP=std::list<MONOM> >
  struct CLASS : public REP {

    typedef C                                       coeff_t;
    typedef C                                       Scalar;
    typedef C                                       value_type;
    typedef O                                       Order;
    typedef REP                                     base_t;
    typedef typename base_t::value_type             monom_t;
    typedef typename base_t::value_type             Monomial;
    typedef typename base_t::iterator               iterator;
    typedef typename base_t::const_iterator         const_iterator;
    typedef typename base_t::reverse_iterator       reverse_iterator;
    typedef typename base_t::const_reverse_iterator const_reverse_iterator;
    typedef Polynomial                              self_t;
    
    Polynomial(): base_t() {}
    Polynomial(const C& c):base_t() {
      if (c!=(coeff_t)0)
        this ->push_back(monom_t(c));
    }


    Polynomial(const C& c, int d, unsigned v);
    Polynomial(const monom_t& m): base_t() { this->push_back(m);};

    bool operator==( const Polynomial& p ) const;
    bool operator!=( const Polynomial& p ) const { return !(*this==p); }
    bool operator==( const C& c ) const;
    bool operator!=( const C& c ) const { return !(*this==c); }

    int nbvar() const;

    static MonomialOrdering * m_order;
    static MonomialOrdering * order() {return m_order;}
    static bool 
    less (const monom_t & m1, const monom_t & m2) {
      //      std::cout<<"less "<< m1<< " "<< m2<<std::endl;
      return m_order->less(m1.begin(),m1.size(), 
                           m2.begin(),m2.size());
    }
  }; 
}//namespace sparse

// #ifndef MMX_WITH_PLUS_SIGN
// #define MMX_WITH_PLUS_SIGN
//   template<class X> bool with_plus_sign(const X& x) {return true;}//x>0;}
// #endif  //MMX_WITH_PLUS_SIGN
//   TMPL bool with_plus_sign(const SparsePolynomial& x) {return true;}


  namespace sparse {
  
  TMPL MonomialOrdering* Polynomial::m_order = (MonomialOrdering*)new O; //DegreeLex;


  TMPL inline
  Polynomial::CLASS(const C& c, int d, unsigned v):base_t() {
      this ->push_back(monom_t(c,d,v));}

  TMPL int 
  Polynomial::nbvar() const{
    int v = 0;
    for (const_iterator it = this->begin(); it != this->end(); ++it)
      v = std::max(v,it->l_variable()+1);
    return v;
  }


  TMPL bool 
  Polynomial::operator== (const Polynomial& p) const {
    return ((base_t)*this==(base_t)p);
  }

  TMPL bool 
  Polynomial::operator== (const C& c) const {
    if(c==0)
      return (this->size()==0);
    else
      return(this->size()==1 && (this->begin()->coeff() == c) && degree(*this->begin())== 0);
  }

  
  //====================================================================  
  TMPL void 
  add(Polynomial & result, const Polynomial& p1, const Polynomial& p2)
  {
    assert(&result != &p1); assert(&result != &p2);
    
    typedef typename Polynomial::const_iterator  const_iterator;
    typedef typename Polynomial::iterator        iterator;
    typedef typename Polynomial::value_type      coeff_type;
    typedef typename Polynomial::Scalar          Scalar;

    const_iterator
      b1=p1.begin(),
      e1=p1.end(),
      b2=p2.begin(),
      e2=p2.end();
    result.resize(p1.size()+p2.size());
    iterator i=result.begin();
    while (b1!=e1 && b2!=e2) {
      if(Polynomial::less(*b2,*b1))
        *i++=*b1++;
      else if(Polynomial::less(*b1,*b2))
        *i++=*b2++;
      else {
        coeff_type c=b1->coeff()+ b2->coeff();
        if (c !=(Scalar)0) {
          *i =*b1;
          i->set_coeff(c);
          ++i;
        }
        ++b1;++b2;
      }
    }
    while (b1!=e1) *i++=*b1++;
    while (b2!=e2) *i++=*b2++;
    
    int m=std::distance(result.begin(),i);
    assert(m>=0);
    result.resize(m);
  }
  //----------------------------------------------------------------------
  TMPL void 
  add(Polynomial & p1, const Polynomial& p2) {
    
    assert(&p1 != &p2);
    typedef typename Polynomial::const_iterator  const_iterator;
    typedef typename Polynomial::iterator        iterator;
    typedef typename Polynomial::value_type      coeff_type;
    
    //    reserve(p1,p1.size()+p2.size());
    iterator        b1=p1.begin();
    const_iterator  b2=p2.begin();
    while (b1!=p1.end() && b2!=p2.end()) {
      if(Polynomial::less(*b2,*b1))
        b1++;
      else if(Polynomial::less(*b1,*b2))
        {
          b1 = p1.insert(b1,*b2); b1++; b2++;
        }
      else
        {
          coeff_type c=b1->coeff() + b2->coeff();
          if (c != (typename Polynomial::Scalar)0) {
            b1->set_coeff(c); ++b1;
          }else{
            b1 = p1.erase(b1);
          }
          ++b2;
        }
    }
    p1.insert(b1,b2,p2.end());
  }
  //----------------------------------------------------------------------
    template<class C, class O, class MONOM, class REP, class I, class J, class M>
  I add(Polynomial & p1, I b1, J e1, const M & m2) {
    typedef typename Polynomial::value_type::coeff_t coeff_type;
    
    while (b1!=e1) {
      if(Polynomial::less(m2,*b1))
        {
          b1++; //COUNT<M>('i');
        }
      else if(Polynomial::less(*b1,m2))
        {
          b1=p1.insert(b1,m2);
          return b1;
        }
      else
        {
          coeff_type c=b1->coeff() + m2.coeff();
          if (c !=0) {
            b1->set_coeff(c);
          }else
            b1 = p1.erase(b1);
          return(b1);
        }
    }
    b1 = p1.insert(b1,m2);
    return b1;
  }
  //----------------------------------------------------------------------
  TMPL void
  add(Polynomial & p, const typename Polynomial::monom_t& m) {
    Polynomial q(m); add(p,q);
  }

  TMPL void
  add(Polynomial & p, const C& c) {
    Polynomial m(c); add(p,m);
  }
  //----------------------------------------------------------------------
  TMPLX void
  add(Polynomial & r, const Polynomial& q, const X& c) {
    Polynomial m(c); add(r,q,m);
  } 
  //----------------------------------------------------------------------
  TMPL void
  add(Polynomial & r, const C& c, const Polynomial& q) {
    Polynomial m(c); add(r,q,m);
  }
  //----------------------------------------------------------------------
  TMPL void 
  sub(Polynomial & result, const Polynomial& p1, const Polynomial& p2)
  {
    assert(&result != &p1); assert(&result != &p2);
    typedef typename Polynomial::const_iterator  const_iterator;
    typedef typename Polynomial::iterator        iterator;
    typedef typename Polynomial::value_type      coeff_type;
    const_iterator
      b1=p1.begin(),
      e1=p1.end(),
      b2=p2.begin(),
      e2=p2.end();
    result.resize(p1.size()+p2.size());
    iterator i=result.begin();
    while (b1!=e1 && b2!=e2) {
      if(Polynomial::less(*b2,*b1))
        *i++=*b1++;
      else if(Polynomial::less(*b1,*b2))
        {
          *i=*b2++; i->coeff()*= -1; i++;
        }
      else {
        coeff_type c=b1->coeff()- b2->coeff();
        if (c !=0) {
          *i =*b1;
          i->set_coeff(c);
          ++i;
        }
        ++b1;++b2;
      }
    }
    while (b1!=e1) *i++=*b1++;
    while (b2!=e2) {*i=*b2++;i->coeff()*=-1;i++;}
    
    int m=std::distance(result.begin(),i);
    assert(m>=0);
    result.resize(m);
  }
  //--------------------------------------------------------------------  
  TMPL void
  sub(Polynomial & p, const Polynomial& q) {
    Polynomial m(q); mul(m,(C)-1);
    add(p,m);
  }
  //--------------------------------------------------------------------
  TMPLX void
  sub(Polynomial & r, const Polynomial& q, const X& c) {
    Polynomial m(C(-c)); add(r,q,m);
  }
  //--------------------------------------------------------------------    
  TMPLX void
  sub(Polynomial & r, const X& c, const Polynomial& q) {
    Polynomial m(c); sub(r,m,q);
  }
  //----------------------------------------------------------------------
  TMPL void 
  mul(Polynomial &r, const C & c) {
    for(typename Polynomial::iterator it=r.begin(); it!=r.end(); ++it) {
      it->coeff()*=c;
    }
  }
  //----------------------------------------------------------------------
  TMPL void 
  mul(Polynomial &r, const Polynomial& p, const C & c) {
    r = p; mul(r,c);
  }

  //----------------------------------------------------------------------
  TMPLX void 
  mul(Polynomial &r, const X & c, const Polynomial& p) {
    r = p; mul(r,c);
  }
  //----------------------------------------------------------------------
  template <class C, class O, class MONOM, class REP, class M> inline  void 
  mul_ext(Polynomial &r, const Polynomial& a, const M & m) {
    r.resize(a.size());
    typename Polynomial::const_iterator b=a.begin(); 
    typename Polynomial::iterator i=r.begin();
    for(; b!=a.end(); ++i, ++b) *i = (*b) * m; 
     //    VECTOR::apply_mult(r, a, m);
  }
  //----------------------------------------------------------------------
  TMPL inline  void 
  mul(Polynomial &r, const Polynomial& a, const typename Polynomial::monom_t & m) {
    mul_ext(r,a,m);
  }
  //----------------------------------------------------------------------
  TMPL inline  void 
  mul(Polynomial &r, const typename Polynomial::monom_t & m) {
    typename Polynomial::iterator i=r.begin();
    for(; i!=r.end(); ++i) *i *= m; 
  }
  //----------------------------------------------------------------------
  template<class C, class O, class MONOM, class REP, class M> inline  void 
  mul_ext_e(Polynomial &result, const Polynomial& a, const M & m)
  {
    typedef typename Polynomial::const_iterator  const_iterator;
    typedef typename Polynomial::iterator        iterator;
    typedef typename Polynomial::monom_t         R;

    result.resize(0);
    iterator i=result.begin();
    R tmp;
    for(const_iterator b=a.begin();b!=a.end(); ++b)
      {
        tmp = (*b) * m;
        (i=result.insert(i,tmp))++;
      }
  }
//----------------------------------------------------------------------
  TMPL inline void 
  mul(Polynomial &r, const Polynomial& a, const Polynomial& b)
  {
    if (a.size()==0) return;
    if (b.size()==0) return;
    if (a.size()<b.size())
      mul(r,a.begin(),a.end(),b);
    else
      mul(r,b.begin(),b.end(),a);
  }

  TMPLX inline void 
  mul(Polynomial &r, 
      const Polynomial& a, 
      const Polynomial& b, const X & o)
  {
    typedef typename Polynomial::const_iterator const_iterator;
    typedef typename Polynomial::iterator       iterator;
    typedef typename Polynomial::monomt_t       M;

    r.resize(0);
    M m;
    for(const_iterator i=b.begin();i !=b.end();i++)
      {
        if(r.size())
          {
            iterator ip = r.begin();
            for(const_iterator j=a.begin();j != a.end();j++) { 
              m  = *j * *i;
              ip = add(r,ip,r.end(),m);//,O()); 
            }
            //MPOLDST::mul_ext(tmp,a,*i);
            //MPOLDST::add<list<M>,vector<M>,O>(r,tmp);
          }
        else
          mul_ext(r,a,*i);
      }
  }

  //----------------------------------------------------------------------
  TMPL inline void 
  mul(Polynomial &r,
      typename Polynomial::const_iterator b,
      typename Polynomial::const_iterator e, 
      const Polynomial& p) {

    
    typedef typename Polynomial::const_iterator  const_iterator;
    int n=std::distance(b,e);
    assert(n>=0);
    if (n==0) r.resize(0);
    if (n==1) 
      mul_ext(r,p,*b);
    else {
      const_iterator med=b; // b+(n/2);
      for(int i=0; i<n/2 && med !=e;i++,med++) {}

      Polynomial tmp0((C)0), tmp1((C)0);
      mul(tmp0,b,med,p);
      mul(tmp1,med,e,p);
      add(r,tmp0,tmp1);
    }

  }
  TMPL inline void 
  mul(Polynomial &r, const Polynomial& p)
  {
    Polynomial s(r); mul(r,s,p); 
  }
  //----------------------------------------------------------------------
  TMPL inline void 
  div(Polynomial &q, const Polynomial& a, const Polynomial& b)
  {
    Polynomial r(a);
    div_rem(q,r,b);
  }

  TMPL inline void 
  div(Polynomial &r, const Polynomial& b)
  {
    Polynomial a(r);
    div(r,a,b);
  }

  TMPL inline void 
  div(Polynomial &f, const typename Polynomial::Scalar& c)
  {
    for(typename Polynomial::iterator it=f.begin(); it != f.end(); it++) 
      it->coeff()/=c;
  }

  //----------------------------------------------------------------------
  TMPL void 
  div(Polynomial &r, const Polynomial& p, const C & c) {
    r = p; div(r,c);
  }

  //----------------------------------------------------------------------
  TMPL inline void 
  rem(Polynomial &r, const Polynomial& a, const Polynomial& b)
  {
    Polynomial q; r=a;
    div_rem(q,r,b);
  }
 
  //----------------------------------------------------------------------
    template<class C, class O, class MONOM, class REP, class U> void
  coefficients(Seq<U>& r,const Polynomial& f,int v) {
    unsigned N = degree(f,v)+1;
    r=Seq<U>(N);
    typename Polynomial::monom_t m;
    for(typename Polynomial::const_iterator it=f.begin(); it != f.end(); it++) {
      m=*it;
      if (v<(int)m.size()) m.set_expt(v,0);
      r[(*it)[v]]+= U(m);
    }
  }

  TMPL void
  coefficients(Seq<C>& r, const Polynomial& f) {
    for(typename Polynomial::const_iterator it=f.begin(); it != f.end(); it++) {
      r<<it->coeff();
    }
  }
  //----------------------------------------------------------------------
  template <class R>  int lvar (const R & p) {
  int v = -1;
  for (typename R::const_iterator it = p.begin(); it != p.end(); ++it)
    v = std::max(v,it->l_variable());
  return v;
  }
  //--------------------------------------------------------------------
  template <class R>  unsigned nbvar (const R & p) {
    return p.nbvar();
  }
  //----------------------------------------------------------------------
  template <class R>  int 
  degree (const R & p) {
    if(!p.size())
      return -1;
    else
      {
        int d = 0;
        for (typename R::const_iterator it = p.begin(); it != p.end(); ++it)
          d = std::max(d,degree(*it));
        return d;
      }
  }

  template <class R>  int 
  degree (const R & p, int i) {
    if(!p.size())
      return -1;
    else
      {
        int d=0;
        for (typename R::const_iterator it = p.begin(); it != p.end(); ++it)
          {
            if(i<=it->nvars()) d = std::max(d,(int)(*it)[i]);
          }
        return d;
      }
  }
  //----------------------------------------------------------------------
  template<class R> typename R::coeff_t& 
  leadingcoeff(R & a)       {return a.begin()->coeff();}
  //----------------------------------------------------------------------
  template<class R> typename R::coeff_t   
  leadingcoeff(const R & a) {return a.begin()->coeff();}
  //----------------------------------------------------------------------
  template<class POL> typename POL::coeff_t 
  coeffof(const POL & p,
          const typename POL::monom_t & mono)
  {
    typedef typename POL::coeff_t coeff_t;
    typename POL::const_iterator m;
    for(m = p.begin();  m!= p.end() && !IsComparable((*m),mono); m++) {}
    if(m != p.end())
      {
        return m->coeff();
      }
    else
      {
        return  coeff_t(0);
      }
  }

  template<class POL> typename POL::const_iterator 
  last_term(const POL& p) {
    assert(p.size()>0);
    typename POL::const_iterator it=p.end();
    it--;
    return it;
  }


  //----------------------------------------------------------------------
  template<class R> void 
  div_rem(R & q, R & a, const R & b0) {
    //    std::cout <<"---- Begin div-rem sparse -----"<<std::endl;
    typedef typename R::monom_t       monom_t;
    typedef typename monom_t::coeff_t coeff_t;
    q= R((coeff_t)0);

    if(a==0) return;
    
    R b(b0);
    //coeff_t cb =b0.begin()->coeff();
    //  b/=cb;
    monom_t mb = (*b.begin()), m;
    R t;
    while( a != (coeff_t)0 && divide(*a.begin(), mb, m) )
      {
        mul(t,b,m);
        sub(a,t);
        if ( m != 0 )  add(q,m);
        //      print(std::cout,a); std::cout<<std::endl;
      }
    //    std::cout <<"---- End div-rem sparse -----"<<std::endl;
  }
  
  //----------------------------------------------------------------------
  TMPL void 
  diff(Polynomial & r, const Polynomial & p, int i)
  {
    typedef typename Polynomial::Monomial Monomial;
    typedef typename Polynomial::Scalar   Scalar;

    assert(&r != &p);
    r= Polynomial((Scalar)0);
    for (typename Polynomial::const_iterator it = p.begin(); it != p.end(); ++it) {
      if((*it)[i]>0){
        Monomial m(*it);
        m*=(typename Polynomial::coeff_t)(*it)[i];
        m.set_expt(i,(*it)[i]-1);
        add(r,m);
      }
    }
  }

  //----------------------------------------------------------------------
  TMPL void 
  shift(Polynomial & r, const Polynomial & p, const typename Polynomial::monom_t& m)
  {
    assert(&r != &p);
    typedef typename Polynomial::monom_t Monomial;
    r= Polynomial(0);
    for (typename Polynomial::const_iterator it = p.begin(); it != p.end(); ++it) {
      add(r,*it*m);
    }
  }

  //----------------------------------------------------------------------
  TMPL void 
  copy(Polynomial& r, const Polynomial& a)
  {
    //to be modified when bug fixed in gcc
    //    using dense::reserve;
    reserve(r,a.size());
    std::copy(a.begin(),a.end(),r.begin());
  }
  //----------------------------------------------------------------------
  TMPL typename Polynomial::coeff_t 
  eval(const Polynomial& p,
       const typename Polynomial::coeff_t &x,
       const typename Polynomial::coeff_t &y) {
    typename Polynomial::coeff_t r(0);
    for(typename Polynomial::const_iterator it=p.begin(); it!=p.end(); ++it)
      {
        typename Polynomial::coeff_t c(it->coeff());
        //for(unsigned int i=0;i<2;++i)
        for(int k=0;k<(int)(*it)[0];k++)
          c*=x;
        for(int k=0;k<(int)(*it)[1];k++)
          c*=y;
        r+=c;
      }
    return r;
  }

  //--------------------------------------------------------------------
  template<class C,class O, class MONOM, class REP, class R, class VCT>  void
  eval_at(R& r,
          const Polynomial& p,
          const VCT &x) {
    r=R(0);
    for(typename Polynomial::const_iterator it=p.begin(); it!=p.end(); ++it)
      {
        R c(it->coeff());
        for(unsigned int j=0;j<x.size()&& j<it->nvars();++j)
          for(int k=0;k<(int)(*it)[j];k++)
            c*=x[j];
        r+=c;
      }
    return r;
  }

  //--------------------------------------------------------------------
    template<class R, class C, class O, class MONOM, class REP, class X> inline void
  eval(R& r, const Polynomial& p, const X& x) {
    r=(R)0;
    for(typename Polynomial::const_iterator it=p.begin(); it!=p.end(); ++it)
      {
        R c(it->coeff());
        for(int k=0;k<(int)(*it)[0];k++)
          c*=x;
        r+=(R)c;
      }
  }

  //----------------------------------------------------------------------
  TMPL inline void
  homogenize(Polynomial& r, const Polynomial& p, const Polynomial& v) {
    Polynomial t;
    int d = sparse::degree(p);
    for(typename Polynomial::const_iterator it=p.begin();it != p.end();it++) {
      t=*it;
      for (int i=0;i<d-degree(*it);i++) mul(t,v);
      add(r,t);
    }
  }

  //--------------------------------------------------------------------
  TMPL inline void
  homogenize(Polynomial& r, const Polynomial& p, int n, const Polynomial& v) {
    Polynomial t;
    int d = degree(p,n);
    for(typename Polynomial::const_iterator it=p.begin();it != p.end();it++) {
      t=*it;
      for (int i=0;i<d-(*it)[n];i++) mul(t,v);
      add(r,t);
    }
  }
  //----------------------------------------------------------------------
  
  template <class MP> MP 
  convert(const MP & P,typename MP::coeff_t x,typename MP::coeff_t y,int ind) {
    typedef typename MP::const_iterator  const_iterator;
    typedef typename MP::coeff_t coeff_t;
    typedef typename MP::monom_t monom_t;
    MP  res;
    int ind1=0;
    int ind2=0;
    for(int i=0;i<=2;i++)
      {
        if(ind!=i) {ind1=i;break;}
      }
    for(int i=0;i<=2;i++)
      {
        if(ind!=i && ind1!=i) {ind2=i;break;}
      }
    for(const_iterator i = P.begin(); i != P.end(); ++i)
      {
        coeff_t  k=i->coeff();
        int puissx=(*i)[ind1];
        int puissy=(*i)[ind2];
        int d=(*i)[ind];
        coeff_t coeff1=pow(x,puissx);
        coeff_t coeff2=pow(y,puissy);
        k*=coeff1;k*=coeff2;
        MP P1=monom_t(k,ind,d);
        res=res+P1;
      }
    return res;
  }
  
  //--------------------------------------------------------------------
    template <class OS, class C, class O, class MONOM, class REP, class VARIABLES> OS& 
  print(OS & os, const Polynomial & P, const VARIABLES & V) {
    if ( P.begin() == P.end() )
      {
        os << '0';
        return os;
      }
    for(typename Polynomial::const_iterator i = P.begin(); i != P.end();i++ ) {
      std::ostringstream sm;
      print(sm,*i,V);
      if(sm.str()[0] != '-' && sm.str()[0] != '+' && i != P.begin())
        os <<'+';
      os<<sm.str();
    }
    return os;
  }


    template <class OS, class C, class O, class MONOM, class REP, class VARIABLES> OS& 
  print_as_double(OS & os, const Polynomial & P, const VARIABLES & V) {
    if ( P.begin() == P.end() )
      {
        os << '0';
        return os;
      }
    for(typename Polynomial::const_iterator i = P.begin(); i != P.end();i++ ) {
      std::ostringstream sm;
      print_as_double(sm,*i,V);
      if(sm.str()[0] != '-' && sm.str()[0] != '+' && i != P.begin())
        os <<'+';
      os<<sm.str();
    }
    return os;
  }
  //--------------------------------------------------------------------
  template <class OS, class C, class O, class MONOM, class REP> OS& 
  print(OS & os, const Polynomial & P)
  {
    return print(os,P,variables::default_);
  }

 
  //----------------------------------------------------------------------
  template <class OS, class C, class O, class MONOM, class REP> OS& 
  print_verbatim(OS & os, const Polynomial & P)
  {
    if ( P.begin() == P.end() )
      {
        os << '0';
        return os;
      }
    for( typename Polynomial::const_iterator i=P.begin(); i != P.end();++i)
      {
        os<<" "<<i->coeff()<<" [";
        for(int l=0;l<i->size();l++) os<<(*i)[l]<<" ";
        os<<"]";
      }
    return os;
  }
  //----------------------------------------------------------------------
  template <class POL, class C> POL 
  shift(typename POL::const_iterator monom, const C& a, int i) {
    
    POL polynom;
    // Definition du type monome.
    typedef typename POL::monom_t monomial;
    // monomiale temporaire.
    monomial tmpmonomial;
    
    int i0 = i;
    int i1 = (i + 1) % 3;
    int i2 = (i + 2) % 3;
    
    // Si la variable est presente dans le monome.
    if ( (a!=0) && ( monom->operator[](i0) !=0 ) )
      {
        polynom = binomial< POL >(C(1.), i0, monom->operator[](i0), a);
        int expt[3];
        expt[i0] = 0;
        expt[i1] = monom->operator[](i1);
        expt[i2] = monom->operator[](i2);
        
        if ( ( expt[i0] == 0 ) && ( expt[i1] == 0 ) && ( expt[i2] == 0 ) )
          tmpmonomial = monomial( monom->coeff(), 0, 0 );
        else   tmpmonomial = monomial( monom->coeff(), 3, expt );
        
        polynom *= tmpmonomial;
      }
    // Si a = 0 alors le monome reste inchange.
    else polynom = ( *monom );
    
    return polynom;
    
  }
  //----------------------------------------------------------------------
  template < class POL, class C > POL 
  scale( const POL &p, C a, int v) {
    
    typedef typename POL::monom_t monomial;
    
    POL  np(p);
    
    for ( typename POL::const_iterator i = np.begin(); i != np.end(); ++i ) {
      i->coeff()*= pow( a, int( i->operator[]( v ) ) );
    }
    return np;
    
  }
  //----------------------------------------------------------------------
  template<class T, class MP, class V> T 
  eval(const MP& p, const V& v) {
    T r(0);
    for(typename MP::const_iterator it=p.begin(); it!=p.end(); ++it)
      {
        T c; let::assign(c,it->coeff());
        for(unsigned i=0;i< it->size()&&i<v.size();++i)
          for(int k=0;k<(*it)[i];k++)
            {
              c*=v[i];
            }
        r+=c;
      }
    return r;
  }


  //----------------------------------------------------------------------
  template<class R, class MP, class V> void 
  eval(R& r, const MP& p, const V& v, unsigned n) {
    r=R(0);
    for(typename MP::const_iterator it=p.begin(); it!=p.end(); ++it)
      {
        R c; 
        let::assign(c,it->coeff());
        for(unsigned i=0;i< it->size()&&i<n;++i)
          for(int k=0;k<(*it)[i];k++) {
            c*=v[i];
          }
        r+=c;
      }
  }

  //----------------------------------------------------------------------
  template<class MP,class X> MP 
  subs( unsigned  var, const X& val, const MP & P)
  {
    MP Result;//("0");
    typedef typename MP::coeff_t coeff_t;
    typedef typename MP::Monomial monom_t;
    
    for( typename MP::const_iterator it = P.begin(); it != P.end(); ++it )
      {
        if ( it->size() > var )
          {
            monom_t m (*it);
            X tmp = pow(val,m[var]);
            m.rep()[var] = 0; 
            Result += tmp*m;
          }
        else 
          Result += *it;
        
      }
    return Result;
  }
  //----------------------------------------------------------------------
  template<class MP> MP 
  subs(const MP & P, int var, typename MP::coeff_t val)
  {
    MP Result;//("0");
    typedef typename MP::coeff_t coeff_t;
    typedef typename MP::Monomial monom_t;
        //  std::cout<<"Begin subs!!"<<endl;
    coeff_t C=0;
    for ( typename MP::const_iterator it = P.begin();
          it != P.end(); ++it )
      {
        monom_t m(it->coeff());
        for(int ind=0; ind<(lvar(P)+1); ++ind)
          {
            if (ind == var)
              {
                for(int k=0;k<(int)(*it)[ind];k++)
                  m*=val;
              }
            else
              {
                if ((*it)[ind] != 0)
                  {
                    monom_t mon(ind,(*it)[ind]);
                    m*= mon;
                  }
              }
          }
        if (m.nvars() == -1)
          {
            C += m.coeff();
          }
        else if (m != 0)
          Result+= m;
      }
    if (C != 0)
      Result += monom_t(C);
          //std::cout<<"End subs!!"<<endl;
    return Result;
  }

  //----------------------------------------------------------------------
  template<class MP> MP 
  subs(const MP & P, char* x, typename MP::coeff_t val) {
    int xi = MP::var(x);
    if (xi<0)
      return P;
    else
      return subs(P,xi,val);
  }
  
  //----------------------------------------------------------------------
  template<class T> void print(const T & x) {std::cout<<x<<std::endl;}
  //----------------------------------------------------------------------
  template<class MP> MP 
  swap(const MP & P, int var_i, int var_j)
  {
    MP Result;//("0");
    typedef typename MP::monom_t monom_t;
    for ( typename MP::const_iterator it = P.begin();
          it != P.end(); ++it )
      {
        monom_t m(it->coeff());
        for(int ind=0; ind<(lvar(P)+1); ++ind)
          {
            if((*it)[ind] !=0 ) {
              if (ind == var_i) {
                m*=monom_t(var_j,(*it)[ind]);
              } else { 
                if (ind == var_j)
                  m*=monom_t(var_i,(*it)[ind]);
                else
                  m*= monom_t(ind,(*it)[ind]);
              }
            }
          }
        Result+=m;
      }
    return Result;
  }

  //----------------------------------------------------------------------
  template<class MP> MP 
  swap(const MP & P, char* x_i, char* x_j)
  {
    typedef typename MP::monom_t monom_t;
    int xi = monom_t::index_of_var(x_i);
    int xj = monom_t::index_of_var(x_j);
    return swap(P,xi,xj);
  }
   
  //--------------------------------------------------------------------
  template <class C, class O, class MONOM, class REP> C
  content(const Polynomial & P)
  {
    C d=0;
    for(typename Polynomial::const_iterator i = P.begin(); i != P.end();i++ ) {
      d= gcd(d,i->coeff());
    }
    return d;
  }

} //namespace sparse

//====================================================================
namespace let {

  TMPL void
  assign(sparse::monomial_seq<C,O,MONOM,REP>& p, const C&c) {
    p = sparse::monomial_seq<C,O,MONOM,REP>(); p.push_back(c);
  }

  template<class U, class V, class O, class UMONOM, class UREP, class VMONOM, class VREP> void
  assign(sparse::monomial_seq<U,O,UMONOM,UREP>& p, 
         const sparse::monomial_seq<V,O,VMONOM,VREP>& q) {

    typedef typename sparse::monomial_seq<U,O,UMONOM,UREP>::Monomial MonomialU;
    typedef typename sparse::monomial_seq<V,O,VMONOM,VREP>::Monomial MonomialV;

    for ( typename  sparse::monomial_seq<V,O,VMONOM,VREP>::const_iterator it = q.begin(); 
          it != q.end(); ++it ) {
      U c; assign(c,it->coeff());
      MonomialU m(c,it->rep());
      p.push_back(m);
    }
  }

}
//======================================================================
# define SparsePolynomial  sparse:: monomial_seq<C,O,MONOM,REP>

template<class A, class B> struct use;
struct operators_of;

TMPL struct use<operators_of,SparsePolynomial> {

    static inline void 
    add (SparsePolynomial &r, const SparsePolynomial &a ) { 
      sparse::add (r, a); 
    }
    static inline void 
    add (SparsePolynomial &r, const SparsePolynomial &a, const SparsePolynomial &b) {
      sparse::add (r, a, b); 
    }
    static inline void
    add (SparsePolynomial &r, const SparsePolynomial &a, const C & b) {
      sparse::add (r, a, b); 
    }
    static inline void
    add (SparsePolynomial &r, const C & a, const SparsePolynomial &b) {
      sparse::add (r, b, a); 
    }
    static inline void
    sub (SparsePolynomial &r, const SparsePolynomial &a ) {
      sparse::sub (r, a ); 
    }
    static inline void
    sub (SparsePolynomial &r, const SparsePolynomial &a, const SparsePolynomial &b) {
      sparse::sub (r, a, b); 
    }
    static inline void
    sub (SparsePolynomial &r, const C & a, const SparsePolynomial &b) {
      sparse::sub (r, a, b); 
    }
    static inline void
    sub (SparsePolynomial &r, const SparsePolynomial &a, const C & b) {
      sparse::sub (r, a, b); 
    }
    static inline void
    mul (SparsePolynomial &r, const SparsePolynomial &a ) {
      sparse::mul (r, a ); 
    }
    static inline void
    mul (SparsePolynomial &r, const SparsePolynomial &a, const SparsePolynomial &b) {
      sparse::mul (r, a, b); 
    }
    static inline void
    mul (SparsePolynomial &r, const SparsePolynomial &a, const C & b) {
      sparse::mul (r, a, b); }
    static inline void
    mul (SparsePolynomial &r, const C & a, const SparsePolynomial &b) {
      sparse::mul (r, b, a); 
    }
    static inline void
    div (SparsePolynomial &r, const SparsePolynomial &a, const SparsePolynomial &b) {
      sparse::div (r, a, b); 
    }
    static inline void
    div (SparsePolynomial &r, const SparsePolynomial &a, const C & b) {
      sparse::div (r, a, b); 
    } 
    static inline void
    div (SparsePolynomial &a, const SparsePolynomial& b) {
      sparse::div (a, b); 
    } 
    static inline void
    div (SparsePolynomial &a, const C & c) {
      sparse::div (a, c); 
    }
    static inline void
    rem (SparsePolynomial &r, const SparsePolynomial &a, const SparsePolynomial &b) {
      sparse::rem (r, b, a); 
    }
};
#undef SparsePolynomial
//====================================================================
}//  namespace mmx 
#undef TMPL
#undef TMPLX
#undef CLASS
#undef Polynomial
#endif // realroot_sparse_monomial_hpp