synaps/mpol/MPol.h

/*********************************************************************
*      This file is part of the source code of the SYNAPS library.
*      Author(s): B. Mourrain, GALAAD, INRIA                 
**********************************************************************
History:
    03/03/00 by bm, update for the new version
    02/03/99 by dbondy,
                  - modification of the operator *=. "p*=p;" fonctionne.
    24/02/99 by dbondy,
        - add a constructor for the MPol class form a char *.
        - add a typedef for the arithmetic coefficient.
    12/12/98 by bm, change p in rep, add size method.
    1/12/98  by dbondy modified the structure of the file for a better
    documentation.
$Id: MPol.H,v 1.1 2005/07/11 11:17:56 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_MPOL_MPOL_H
#define SYNAPS_MPOL_MPOL_H

#include <stdlib.h>
#include <string>
#include <list>
#include <synaps/init.h>
#include "synaps/base/shared_object.h"
#include "synaps/base/parser.h"
#include "synaps/base/COUNT.h"
#include "synaps/arithm/REF.h"
#include "synaps/arithm/let.h"
#include "synaps/upol/UPolDse.h"
#include "synaps/mpol/Variables.h"
#include "synaps/mpol/dynamicexp.h"
#include "synaps/mpol/DegLex.h"
#include "synaps/mpol/Monom.h"
#include "synaps/mpol/MPOLDST.m"

//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//----------------------------------------------------------------------
namespace linalg { template<class C> struct rep1d; }
//----------------------------------------------------------------------

template <class C,
          class O=DegLex,
          class R=std::list<Monom<C,dynamicexp<unsigned> > > >
struct MPol {

  shared_object<R>  data;
  R &       rep()       {return (*data);}
  const R & rep() const {return (*data);}

  // Types
  typedef typename R::value_type             monom_t;
  typedef C                                  coeff_t;
  typedef O                                  order_t;
  typedef typename R::iterator               iterator;
  typedef typename R::const_iterator         const_iterator;
  typedef typename R::reverse_iterator       reverse_iterator;
  typedef typename R::const_reverse_iterator const_reverse_iterator;
  typedef MPol<C,O,R>                        self_t;

  // Constructors
  MPol():data() {rep().resize(0);}
  MPol(const monom_t & m): data()
    {
      rep().resize(0); rep().insert(rep().begin(),m);
    }

  MPol(const coeff_t & c):data()
    {
      if(c !=0) {
        rep().resize(0); rep().insert(rep().begin(),monom_t(c));
      }else
        rep().resize(0);
    }

  MPol(const int & i): data()
    {
      if(i !=0)
        { rep().resize(0);rep().insert(rep().begin(),monom_t(i));}
      else
        { rep().resize(0);}
    }

  MPol(int, const monom_t *);
  MPol(const self_t & P): data(P.data) {}
  MPol(const char *, Variables& vars = Variables::default_);
  MPol(const char *, const Variables& vars);

  template<class S>
  MPol(const VAL<S>& M) {using let::assign; assign(*this,M.rep());}


  template<class Q>
  MPol(const MPol<C,Q,R>& p) : data()
  {
    for(const_iterator it=p.begin();it != p.end(); it++)
      (*this) += *it;
  }
  
  // Destructor
  // no dynamic allocation at this level.

  //  static Variable var;
  static int index(const std::string & s)     { return monom_t::var[s]; }
  static std::string var(int i)               { return monom_t::var[i]; }
  static void var(int i, const std::string& s){ monom_t::var.let(i,s); }
  static int  var(char* s)                    { monom_t::var=s; return monom_t::var[s]; }


  // Assignment operators
  self_t & operator =(const self_t & x)  {data=x.data; return *this;}
  self_t & operator =(const coeff_t & c) {*this=self_t(c); return *this;}
  self_t & operator =(int n)             {*this=self_t(n); return *this;}
  self_t & operator =(const monom_t & m) {*this=self_t(m); return *this;}
  template<class S>
  self_t & operator =(const VAL<S> & M)  {assign(*this,M.rep());return *this;}

  self_t & operator*=(const coeff_t &);
  self_t & operator/=(const coeff_t &);

  self_t & operator*=(const monom_t &);
  self_t & operator/=(const monom_t &);
  //  self_t & operator+=( const monom_t & m );
  
  self_t & operator+=(const self_t &);
  self_t & operator-=(const self_t &);
  self_t & operator*=(const self_t &);
  self_t & operator/=(const self_t & Q);

  // Boolean operators and functions
  bool operator==(int n) const;
  bool operator!=(int n) const {return !(*this==n);}

  // Iterator functions
  iterator                begin ()       {return rep().begin(); }
  const_iterator          begin () const {return rep().begin(); }
  reverse_iterator       rbegin ()       {return rep().rbegin();}
  const_reverse_iterator rbegin () const {return rep().rbegin();}

  iterator                  end ()       {return rep().end(); }
  const_iterator            end () const {return rep().end(); }
  reverse_iterator         rend ()       {return rep().rend();}
  const_reverse_iterator   rend () const {return rep().rend();}

  int size()  const {return rep().size();}
  int nbvar() const {return MPOLDST::nbvar(*this);}

  template<class T>
  coeff_t operator()(const T & p) const;

};
//======================================================================
// template<class C, class O, class R> Variable MPol<C,O,R>::var;
//======================================================================
// Input/Output operator
template<class C, class O, class R>
std::istream & operator >>(std::istream &, MPol<C,O,R> &);
template<class C, class O, class R>
std::ostream & operator <<(std::ostream &, const MPol<C,O,R> &);
template<class C, class O, class R>
bool operator == (const MPol<C,O,R> &, const MPol<C,O,R> &);
template<class C, class O, class R>
bool operator != (const MPol<C,O,R> &, const MPol<C,O,R> &);
template<class C, class O, class R>
bool operator < (const MPol<C,O,R> &, const MPol<C,O,R> &);
template<class C, class O, class R>
bool operator <= (const MPol<C,O,R> &, const MPol<C,O,R> &);
template<class C, class O, class R>
bool operator > (const MPol<C,O,R> &, const MPol<C,O,R> &);
template<class C, class O, class R>
bool operator >= (const MPol<C,O,R> &, const MPol<C,O,R> &);
//======================================================================
// ARITHMETIC OPERATIONS.
//======================================================================
#ifdef WITH_TEMPLATE_EXP
        #include <synaps/MPol_with_texp.h>
#else
//======================================================================
template <class C, class O, class R>
MPol<C,O,R> operator+(const MPol<C,O,R> & a, const MPol<C,O,R> & b);
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator+(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::coeff_t & c);
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator+(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::monom_t & c);
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator+(const typename MPol<C,O,R>::coeff_t & c,const MPol<C,O,R> & a);
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator-(const MPol<C,O,R> & a, const MPol<C,O,R> & b);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const MPol<C,O,R> & a,const typename MPol<C,O,R>::coeff_t & c);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const typename MPol<C,O,R>::coeff_t & c,const MPol<C,O,R> & a);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const MPol<C,O,R> & a);
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator*(const MPol<C,O,R> & a, const MPol<C,O,R> & b);
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator*(const MPol<C,O,R> & a, const C & c);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(const typename MPol<C,O,R>::coeff_t & c,
                      const MPol<C,O,R> & a);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(int c,const MPol<C,O,R> & a);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::monom_t & m);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(const typename MPol<C,O,R>::monom_t & m,
                      const MPol<C,O,R> & a);
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator/(const MPol<C,O,R> & a, const MPol<C,O,R> & b) ;
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator%(const MPol<C,O,R> & a, const MPol<C,O,R> & b) ;
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator^(const MPol<C,O,R> & a, unsigned int n);
#endif
//----------------------------------------------------------------------
template <class C, class O, class R>
int NbVar(const MPol<C,O,R> & p)
{
  return MPOLDST::nbvar(p);
}

template <class C, class O, class R>
int LastVar(const MPol<C,O,R> & p)
{
  return MPOLDST::lvar(p);
}


//====================================================================
template <class C, class O, class R> inline
int degree (const MPol<C,O,R> & p)
{
  return MPOLDST::degree(p.rep());
}
template <class C, class O, class R> inline
int degree (const MPol<C,O,R> & p, unsigned int i)
{
  return MPOLDST::degree(p.rep(),i);
}
//====================================================================
template <class C, class O, class R> inline
MPol<C,O,R> diff(const MPol<C,O,R> & p, unsigned i)
{
  MPol<C,O,R> r;
  MPOLDST::diff(r,p,i);
  return r;
}
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> diff(const MPol<C,O,R> & p, char* x, Variables& vars = Variables::default_ )
{
  typedef typename MPol<C,O,R>::monom_t monom_t;
  int i = vars[x];
  if(i<0)
    return MPol<C,O,R>(0);
  else
    {
      MPol<C,O,R> r;
      MPOLDST::diff(r,p,i);
      return r;
    }
}
//====================================================================
namespace let 
{
  template<class P,class S, class C, class O, class R>
  void assign(MPol<double,P,S>& r, const MPol<C,O,R> & p)
  {
    typedef typename MPol<double,P,S>::monom_t monom_t;
    double c;
    for(typename MPol<C,O,R>::const_iterator it=p.begin(); it !=p.end(); ++it)
      {
        let::assign<double,C>(c,it->coeff());
        monom_t m(c,it->rep());
        r+=m;
      }
  }

  template<class Ca, class Ra,
           class Cb, class Ob, class Rb >
  void convert( UPolDse<Ca,Ra> & p, const MPol<Cb,Ob,Rb> & mp, int v )
  {
    int dv;
    p.resize( (dv = degree(mp,v))+1 );
    std::fill( p.begin(), p.end(), 0);
        
    typedef typename MPol<Cb,Ob,Rb>::const_iterator pmit;
    typedef typename MPol<Cb,Ob,Rb>::monom_t monom;
    using namespace let;
    for ( pmit it = mp.begin(); it != mp.end(); it ++ )
      {
        monom m(*it);
        if ( m.size() > v ) m.rep()[v] = 0;
        p[ (*it)[v] ] += convert(m,type::As<Ca>());
      };
  };
  
  template<class Ca, class Ra,
           class Cb, class Ob, class Rb,
           class Rp, class Ru>
  void convert( UPolDse< UPolDse<Ca,Ra>, Rp >& p, 
                const UPolDse< MPol<Cb,Ob,Rb>, Ru >&  ump, int v )
  {
    p.resize( ump.size() );
    for ( int i = 0; i < ump.size(); i++ ) convert(p[i],ump[i],v);
  };
  
  template<class Ca, class Ra, 
           class Rp, 
           class Cb, class Ob, class Rb >

  void convert( UPolDse< UPolDse<Ca,Ra>, Rp >& p, 
                const MPol<Cb,Ob,Rb> & mp, int a, int b )
  {
    UPolDse< MPol<Cb,Ob,Rb> > tmp;
    convert( tmp, mp, a );
    convert( p,  tmp, b );
  };
  
  template<class C, class O, class Ra, class Rb, class Rc>
  void convert( UPolDse< UPolDse<C,Ra>, Rb >& p, 
                const MPol<C,O,Rc> & mp, 
                const char * a, const char * b, 
                const Variables & vars = Variables::default_  )
  {
    convert( p, mp, vars[a], vars[b] );
  };

  template<class C, class O, class Rm, class Ru> inline
  void convert( MPol<C,O,Rm> & result, const UPolDse< MPol<C,O,Rm> , Ru >  & p, int v )  
  {
    int nvar = NbVar(p);
    int tmp[ nvar + 1 ];
//     MPol<C,O,Rm> result;
    for ( int i = 0; i < p.size(); i ++ ) 
      for ( typename MPol<C,O,Rm>::coeff_t::const_iterator itc = p[i].begin(); itc != p[i].end(); itc++ )
        {
          typename MPol<C,O,Rm>::monom_t tmp = *itc;
          if ( tmp.rep.size() <= v ) tmp.rep().resize( v+1 );
          tmp.rep()[v] = i;
          result += tmp;
        };
  };
}
//====================================================================
namespace UPOLDAR 
{
  template <class C, class O, class R>
  std::ostream & print_as_coeff(std::ostream & os, const MPol<C,O,R> & p)
  {
    return os<<"+("<<p<<")";;
  }
}
//----------------------------------------------------------------------
// Constructors
//----------------------------------------------------------------------
template <class C,class O, class R>
inline MPol<C,O,R>::MPol (int s, const monom_t *t)
{
        for (int i = 0; i<s; i++) *this += t[i];
}

extern "C" char *synaps_inputptr;
template <class C,class O, class R>
inline MPol<C,O,R>::MPol (const char * s, Variables&  vars ): data() {
  rep().resize(0);
  int n = strlen(s);
  if (s[n-1]=='\n') {n--;}
  if(s[n-1]==';')
    {
      synaps_input = new char[n+1];
      memcpy(synaps_input,s,n); synaps_input[n]='\0';
      synaps_inputlim = synaps_input + n;
    }
  else
    {
      synaps_input = new char[n+2];
      memcpy(synaps_input,s,n); 
      synaps_input[n]=';'; synaps_input[n+1]='\0';
      synaps_inputlim = synaps_input + n+1;
    }
  synaps_inputptr=synaps_input;
  monom_t m(1);
  bool sign = true;
  bool first = true;
  int Ask;

  for (;;) {
    Ask = yylex();
    if (Ask == COEF) { /* Monome constant */
      coeff_t c;
      let::assign(c,yylval);
      free(yylval);
      m *= c;
    }
    else if (Ask == XID) {
      char* p = strchr(yylval,'^');
      if (p == NULL)
        {
          m *= monom_t(vars[std::string(yylval)],1);
        }
      else {
        p++;
        m *= monom_t(vars[std::string(yylval,p-1)], atoi(p));
      }
      free(yylval);
    }
    else if (Ask == MULT);
    else if (Ask == SUM) {
      if (sign) {
        if(m!=0) *this += m;
        m = monom_t(1);
      } else {
        if(m!=0) *this -= m;
        sign = true;
        m = monom_t(1);
      }
    }
    else if (Ask == MINUS) {
      if (first)
        sign = false;
      else if (sign) {
        if(m!=0) *this += m;
        sign = false;
        m = monom_t(1);
      }
      else {
        if(m!=0) *this -= m;
        m = monom_t(1);
      }
    }
    else if (Ask == TERMINATOR) {
      if (sign){
        if(m!=0) *this += m;
      }else{
        if(m!=0) *this -= m;
      }
      break;
    }
    first = false;
  }
  delete[] synaps_input;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
  inline MPol<C,O,R>::MPol (const char * s, const Variables&  vars )
{
  Variables V(vars);
  *this= MPol<C,O,R>(s,V);
}
//----------------------------------------------------------------------
// template <class C,class O, class R>
// inline MPol<C,O,R> & MPol<C,O,R>::operator += (const typename MPol<C,O,R>::monom_t & m){
//   if (m.coeff()!=0)
//     MPOLDST::insert<C,O,R>(rep(), rep().begin(), m);
//   return *this;
// }
//----------------------------------------------------------------------
template <class C,class O, class R>
inline MPol<C,O,R> & MPol<C,O,R>::operator+= (const MPol<C,O,R> & P)
{
  if (P==0) return *this;
  if (*this==0)
    {*this = P; return *this;}
  else
    MPOLDST::add<O,R>(rep(),P.rep());
  return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
inline MPol<C,O,R> & MPol<C,O,R>::operator-= (const MPol<C,O,R> & P){
  if (P==0) return *this;
  if (*this==0) {*this = -P; return *this;}
  *this += (-P);
  return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
inline MPol<C,O,R> & MPol<C,O,R>::operator*= (const typename MPol<C,O,R>::coeff_t & c){
  MPol<C,O,R> P;
  for (iterator i = begin(); i != end(); i++) {
    (*i)*= c;
  }
  return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
inline MPol<C,O,R> & MPol<C,O,R>::operator/= (const typename MPol<C,O,R>::coeff_t & c){
  MPol<C,O,R> P;
  for (iterator i = begin(); i != end(); i++) { (*i)/= c;
  }
  return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
inline MPol<C,O,R> & MPol<C,O,R>::operator/= (const typename MPol<C,O,R>::monom_t & m){

  if ( (*this!=0)  && (m.coeff() !=0) ) {
    MPol<C,O,R> P(*this);
    for (iterator i = P.begin(); i != P.end(); i++) {
      int n=m.nvars();
      if (n<=(*i).nvars())
        while (((*i)[n] >= m[n] ) && ( n >= 0)) n--;
      if (n<0)
        *i=div(*i,m);
      else return *this;
    }
    (rep()=P.rep());
  }
  return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
inline MPol<C,O,R> & MPol<C,O,R>::operator*= (const typename MPol<C,O,R>::monom_t & m){
 if(m !=0)
   {
     MPol<C,O,R> P;
     for (iterator i = begin(); i != end(); i++) {
       (*i)*= m;
     }
   }
 else
   *this=0;
 return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
inline  MPol<C,O,R> & MPol<C,O,R>::operator*= (const MPol<C,O,R> & b)
{
  assert(&rep() != &b.rep());
  if(*this !=0)
    if(b!=0)
      {
        const self_t a(*this);
        *this=0;
        MPOLDST::mul(rep(),a.rep(),b.rep(),O());
      }else
        (*this)= self_t(0);
  return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
inline  MPol<C,O,R> & MPol<C,O,R>::operator/= (const MPol<C,O,R> & Q){
  assert(Q != 0);
  MPol<C,O,R> P(*this), quo;
  MPOLDST::div_rem_x(quo,P,Q);
  *this=quo;
  return *this;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
template<class T>
typename MPol<C,O,R>::coeff_t
MPol<C,O,R>::operator()(const T & p) const
{
  coeff_t r(0);
  for(const_iterator it=begin(); it!=end(); ++it)
    {
      coeff_t s(it->coeff());
      for(unsigned int i=0;i<(*it).rep().size();++i)
        s*=pow(p[i],(int)(*it)[i]);
      r+=s;
    }
  return r;
}
//----------------------------------------------------------------------
// Boolean operators and functions
//----------------------------------------------------------------------
template <class C,class O, class R> inline
bool MPol<C,O,R>::operator ==(int n) const
{
  if(n==0)
    if (rep().size()==0)
      return true;
    else
      return false;
  else
    if (rep().size() ==0)
      return false;
    else
      return(rep().size()==1 && (rep().begin()->coeff() ==n));
 /* if ( n==0 && rep().size() == 0 ) 
    return true;
  if (rep().size() == 0 )   return false;
  return(rep().size()==1 && (rep().begin()->coeff() ==n));
 */
}
//----------------------------------------------------------------------
template <class C,class O, class R>
bool operator== (const MPol<C,O,R> & P, const MPol<C,O,R> & Q)
{
        return P.rep() == Q.rep();
}
//----------------------------------------------------------------------
template <class C,class O, class R>
bool operator!= (const MPol<C,O,R> & P, const MPol<C,O,R> & Q)
{
        return !(P.rep() == Q.rep());
}
//----------------------------------------------------------------------
template <class C,class O, class R>
bool operator< (const MPol<C,O,R> & P, const MPol<C,O,R> & Q)
{
        typename MPol<C,O,R>::const_iterator iP = P.begin();
        typename MPol<C,O,R>::const_iterator iQ = Q.begin();
        while (iP != P.end() && iQ != Q.end() && IsComparable(*iP,*iQ)) {
                ++iP;
                ++iQ;
        }
        if (iQ == Q.end())
                return false;
        if (iP == P.end())
                return true;
        return O::less(*iP,*iQ);
}
//----------------------------------------------------------------------
template <class C,class O, class R>
bool operator> (const MPol<C,O,R> & P, const MPol<C,O,R> & Q){
        return Q < P;
}
//----------------------------------------------------------------------
template <class C,class O, class R>
bool operator>= (const MPol<C,O,R> & P, const MPol<C,O,R> & Q){
        return !(P < Q);
}
//----------------------------------------------------------------------
template <class C,class O, class R>
bool operator<= (const MPol<C,O,R> & P, const MPol<C,O,R> & Q){
        return !(Q < P);
}
//======================================================================
// Input/putput operators
//======================================================================
template <class C,class O, class R>
std::istream & operator>> (std::istream & is, MPol<C,O,R> & P)
{
        std::string s,ss = "";
        int p;
        for (;;) {
                is >> s;
                p = s.find(';');
                if (p== -1) s.find(':');
                if (p== -1) s.find(',');
                if (p== -1)
                        ss += s;
                else {
                        ss += s.substr(0,p+1);
                        break;
                }
        }
        P = MPol<C,O,R>(ss.c_str());
        return is;
}
//----------------------------------------------------------------------
template <class C, class O, class R>
std::ostream & operator<< (std::ostream & os, const MPol<C,O,R> & P)
{
  MPOLDST::print(os,P); return os;
}

template <class C, class O, class R>
std::ostream & print (std::ostream & os, const MPol<C,O,R> & P, const Variables& V=Variables::default_)
{
  MPOLDST::print(os,P,V); return os;
}
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator+(const MPol<C,O,R> & a, const MPol<C,O,R> & b)
{
  MPol<C,O,R> p;
  if(&p == &a && &p == &b)
    MPOLDST::add<O,R>(p.rep(),MPol<C,O,R>(b).rep());
  else if(&p == &a)
    MPOLDST::add<O,R>(p.rep(),b.rep());
  else if(&p == &b)
    MPOLDST::add<O,R>(p.rep(),a.rep());
  else
    MPOLDST::add<O,R>(p.rep(),a.rep(),b.rep());
  return p;
}
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator+(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::coeff_t & c)
{
  return a + MPol<C,O,R>(c);
}
//--------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator+(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::monom_t & m)
{
  return a + MPol<C,O,R>(m);
}
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator+(const typename MPol<C,O,R>::coeff_t & c,
                      const MPol<C,O,R> & a)
{
  return a + MPol<C,O,R>(c);
}
//----------------------------------------------------------------------
// operator -
template <class C, class O, class R>
MPol<C,O,R> operator-(const MPol<C,O,R> & a, const MPol<C,O,R> & b)
{
  return (a + (-b));
}
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::coeff_t & c)
{
  return a - MPol<C,O,R>(c);
}
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const typename MPol<C,O,R>::coeff_t & c,
                      const MPol<C,O,R> & a)
{
  return MPol<C,O,R>(c) - a;
}
//--------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::monom_t & m)
{
  return a - MPol<C,O,R>(m);
}
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const typename MPol<C,O,R>::monom_t & m,
                      const MPol<C,O,R> & a)
{
  return MPol<C,O,R>(m) - a;
}
//----------------------------------------------------------------------
// Unary operator -
template <class C, class O, class R> inline
MPol<C,O,R> operator-(const MPol<C,O,R> & a)
{
  return a*((typename MPol<C,O,R>::coeff_t)-1);
}
//----------------------------------------------------------------------
// Product of polynomials
template <class C, class O, class R>
MPol<C,O,R> operator*(const MPol<C,O,R> & a, const MPol<C,O,R> & b)
{
  MPol<C,O,R> p;
  if(a !=0 && b !=0)
    MPOLDST::mul(p.rep(),a.rep(),b.rep(),O());
  else
    p= MPol<C,O,R>(0);
  return p;
}
//----------------------------------------------------------------------
// Product of a polynomial by a "scalar"
template <class C, class O, class R>
MPol<C,O,R> operator*(const MPol<C,O,R> & a, const C & c)
{
  MPol<C,O,R> p;
  if(a!=0 && c!=0)
    MPOLDST::mul_ext(p.rep(),a.rep(),c);
 else
    p= MPol<C,O,R>(0);
  return p;
}
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(const typename MPol<C,O,R>::coeff_t & c,
                     const MPol<C,O,R> & a)
{
  return a*c;
}
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(int c,const MPol<C,O,R> & a)
{
  return a*MPol<C,O,R>::coeff_t(c);
}
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(const MPol<C,O,R> & a,
                      const typename MPol<C,O,R>::monom_t & m)
{
  MPol<C,O,R> p;
  if(a!=0)
    MPOLDST::mul_ext(p.rep(),a.rep(),m);
  else
    p= MPol<C,O,R>(0);
  return p;
}
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator*(const typename MPol<C,O,R>::monom_t & m,
                     const MPol<C,O,R> & a)
{
  return a*m;
}
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator/(const MPol<C,O,R> & A, const MPol<C,O,R> & B)
{
  MPol<C,O,R> x=A; x/=B; return x;
}
//----------------------------------------------------------------------
template <class C, class O, class R> inline
MPol<C,O,R> operator%(const MPol<C,O,R> & a, const MPol<C,O,R> & b)
{
  assert(b != 0);
  MPol<C,O,R> r(a), quo;
  MPOLDST::div_rem_x(quo,r,b);
  return r;
}
//----------------------------------------------------------------------
template <class C, class O, class R>
MPol<C,O,R> operator^(const MPol<C,O,R> & a, unsigned int n)
{
  MPol<C,O,R> p;
  if(n ==0) p = 1;
  else {
    p=1;
    unsigned N = n;
    while(N>0) {p*=a; N--;}
  }
  return p;
}


// //--------------------------------------------------------------------
// /** This function let the variable of index \c i appear as the string
//     \c s. This association is also used, when a polynomial is parsed.
//     In other words, if the variable \c s is parsed, then its index is
//     \c i.
// **/
// template<class T> void Variable(unsigned i, const std::string & s)
// {
//   T::monom_t::_var_of_index[i+1]=s;
//   T::monom_t::_index_of_var[s]=i+1;
// }
// //----------------------------------------------------------------------
// /**  Return the variable of index \c i **/
// template<class T> std::string Variable(unsigned i)
// {
//   return T::monom_t::var_of_index(i);
// }
// //----------------------------------------------------------------------
// /**  Return the index of a variable \c s **/
// template<class T> unsigned IndexOf(const std::string & s)
// {
//   return T::monom_t::index_of_var(s);
// }
//--------------------------------------------------------------------
template<class C, class O, class R, class X>
MPol<C,O,R> subs( int variable, const X & value, class MPol<C,O,R> & polynom )
{
  //  using namespace MPOLDST;
  return MPOLDST::subs( variable, value, polynom );
};
__END_NAMESPACE_SYNAPS
//======================================================================
#endif // SYNAPS_MPOL_MPOL_H