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 <basix/basix.hpp>
# 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_iterator(r,a.begin(),a.end(),b);
    else
        mul_iterator(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_iterator(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);
    //std::cout<<"mul iter "<<n<<std::endl;
    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, tmp1; //((C)0), tmp1((C)0);
        mul_iterator(tmp0,b,med,p);
        mul_iterator(tmp1,med,e,p);
        add(r,tmp0,tmp1);
    }
    //std::cout<<"mul iter end "<<n <<std::endl;
}
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 j=0; j< (unsigned)x.size() && j< (unsigned)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 C, class O, class MONOM, class REP, class VARIABLES>
std::string
to_string(const Polynomial & P, const VARIABLES & V) {
    //     std::cout<<"print"<<std::endl;
    //     std::string s;
    if ( P.begin() == P.end() )
    {
        return std::string("0");
    }
    std::string os;
    for(typename Polynomial::const_iterator i = P.begin(); i != P.end();i++ ) {
        std::string s = to_string(*i, V);
        if(s[0] != '-' && s[0] != '+' && i != P.begin())
            os += "+";
        os += s;
        //print(s,*i,V);
    }
    //      std::cout<<"print end"<<std::endl;
    return os;
}

//--------------------------------------------------------------------
template <class OS, class C, class O, class MONOM, class REP, class VARIABLES> OS&
print(OS & os, const Polynomial & P, const VARIABLES & V) {
    //     std::cout<<"print"<<std::endl;
    //     std::string s;
    if ( P.begin() == P.end() )
    {
        os<<"0";
        return os;
    }
    for(typename Polynomial::const_iterator i = P.begin(); i != P.end();i++ ) {
        std::string s = to_string(*i, V);
        if(s[0] != '-' && s[0] != '+' && i != P.begin())
            os<< "+";
        os<<s;
        //print(s,*i,V);
    }
    //      std::cout<<"print end"<<std::endl;
    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().c_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());
}


//----------------------------------------------------------------------
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((int)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!! " << std::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(coeff_t(1),(*it)[ind],ind);
                    m*= mon;
                }
            }
        }
        if (m.nvars() == -1)
        {
            C += m.coeff();
        }
        else if (m != 0)
        {
            Result+= m;
        }
    }
    if (C != 0)
        Result += monom_t(C);

    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