synaps/mpol/MPOLDST.m

/*********************************************************************
*   This file is part of the source code of the SYNAPS.       *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
**********************************************************************
$Id: MPOLDST.h,v 1.2 2005/09/28 06:40:50 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_MPOL_MPOLDST_H
#define SYNAPS_MPOL_MPOLDST_H

#include <list>
#include <vector>
#include <algorithm>
#include <synaps/init.h>
#include <synaps/linalg/VECTOR.m>
#include <synaps/arithm/let.h>
#include <synaps/arithm/binomial.h>
#include <synaps/arithm/sign.h>
#include <synaps/mpol/Variables.h>
#include <synaps/mpol/binomial.h>
#include <synaps/mpol/Monom.h>
#include <synaps/mpol/MONOMIAL.m>
#include <synaps/arithm/pow.h>


__BEGIN_NAMESPACE_SYNAPS

namespace MPOLDST
{
  template<class O, class R>
  void add(R & result, const R& p1, const R& p2)
  {
    assert(&result != &p1); assert(&result != &p2);

    typedef typename R::const_iterator  const_iterator;
    typedef typename R::iterator iterator;
    typedef typename R::value_type::coeff_t 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(O::less(*b2,*b1))
        *i++=*b1++;
      else if(O::less(*b1,*b2))
        *i++=*b2++;
      else {
          coeff_type c=b1->coeff()+ b2->coeff();
          if (c !=0) {
                *i =*b1;
                i->setcoeff(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);
  }

  template<class R, class O>
  void sub(R & result, const R& p1, const R& p2)
  {
    assert(&result != &p1); assert(&result != &p2);
    typedef typename R::const_iterator  const_iterator;
    typedef typename R::iterator iterator;
    typedef typename R::value_type::coeff_t 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(O::less(*b2,*b1))
        *i++=*b1++;
      else if(O::less(*b1,*b2))
        {
          *i=*b2++; i->coeff()*= -1; i++;
        }
      else {
          coeff_type c=b1->coeff()- b2->coeff();
          if (c !=0) {
                *i =*b1;
                i->setcoeff(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);
  }

  template<class O, class R>
  void add(R & p1, const R& p2)
  {
    assert(&p1 != &p2);
    typedef typename R::const_iterator  const_iterator;
    typedef typename R::iterator iterator;
    typedef typename R::value_type::coeff_t 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(O::less(*b2,*b1))
        b1++;
       else if(O::less(*b1,*b2))
        {
          b1 = p1.insert(b1,*b2); b1++; b2++;
        }
      else
        {
          coeff_type c=b1->coeff() + b2->coeff();
          if (c !=0) {
            b1->setcoeff(c); ++b1;
          }else{
            b1 = p1.erase(b1);
          }
          ++b2;
        }
    }
    p1.insert(b1,b2,p2.end());
  }
//----------------------------------------------------------------------
 template<class R, class I, class J, class M, class O>
 I add(R & p1, I b1, J e1, const M & m2, const O & o)
  {
    typedef typename R::value_type::coeff_t coeff_type;

    while (b1!=e1) {
      if(O::less(m2,*b1))
        {
          b1++; //COUNT<M>('i');
        }
      else if(O::less(*b1,m2))
        {
          b1=p1.insert(b1,m2);
          return b1;
        }
      else
        {
          coeff_type c=b1->coeff() + m2.coeff();
          if (c !=0) {
                b1->setcoeff(c);
          }else
            b1 = p1.erase(b1);
          return(b1);
        }
    }
    b1 = p1.insert(b1,m2);
    return b1;
  }
//----------------------------------------------------------------------
  template <class R, class S, class C> inline
  void mul_ext(R &result, const S& a, const C & m)
  {
    VECTOR::apply_mult(result, a, typename R::value_type(m));
  }
//----------------------------------------------------------------------
  template <class R,class C> inline
  void mul_ext_e(std::list<R> &result, const std::list<R>& a, const C & m)
  {
    typedef typename std::list<R>::const_iterator  const_iterator;
    typedef typename std::list<R>::iterator iterator;

    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))++;
      }
  }
//----------------------------------------------------------------------
  template<class R, class O> inline
  void mul(R &r, const R& a, const R& b, const O & o)
  {
    if (a.size()<b.size())
      MPOLDST::mul<R,O>(r,a.begin(),a.end(),b);
    else
      MPOLDST::mul<R,O>(r,b.begin(),b.end(),a);
  }

  template <class M, class O> inline
  void mul(std::list<M> &r, const std::list<M>& a, const std::list<M>& b, const O & o)
  {
    typedef typename std::list<M>::const_iterator const_iterator;
    typedef typename std::list<M>::iterator       iterator;

    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 = MPOLDST::add(r,ip,r.end(),m,O()); }
            //MPOLDST::mul_ext(tmp,a,*i);
            //MPOLDST::add<list<M>,vector<M>,O>(r,tmp);
          }
        else
          MPOLDST::mul_ext(r,a,*i);
      }
  }

//----------------------------------------------------------------------
template <class R, class O> inline
void mul(R &r,
          typename R::const_iterator b,
          typename R::const_iterator e, const R& p)
{
    typedef typename R::const_iterator  const_iterator;
    long int n=std::distance(b,e);

    assert(n>=0);
    if (n==0) r.resize(0);
    if (n==1) MPOLDST::mul_ext(r,p,*b);
    else {

      const_iterator med=b+(n/2);
      R tmp(0), tmp1(0);
      MPOLDST::mul<R,O>(tmp,b,med,p);
      MPOLDST::mul<R,O>(tmp1,med,e,p);
      MPOLDST::add<R,O>(r,tmp,tmp1);
    }

}

//----------------------------------------------------------------------
template <class R>  int lvar (const R & p)
{
  unsigned v = 0;
  for (typename R::const_iterator it = p.begin(); it != p.end(); ++it)
    v = std::max(v,it->size());
  return (v-1);
}
//--------------------------------------------------------------------
template <class R>  unsigned nbvar (const R & p)
{
  return lvar(p)+1;
}
//----------------------------------------------------------------------
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 R>
void div_rem(R & q, R & a, const R & b)
{

  typedef typename R::monom_t monom_t;

  q= R(0);
  if(a!=0)
    {
      typename R::monom_t mb = (*b.begin());
      typename R::monom_t m = (*a.begin())/mb;
      R t;
      while( a !=0 && m != 0) {
        t=b*m;
        a-=t;
        q+=m;
        if(a !=0) m =(*a.begin())/mb;
      }
    }
}
//----------------------------------------------------------------------
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 R>
void div_rem_x(R & q, R & a, const R & b0)
{
  //  cout <<"-------------------------------------------------------------"<<endl;
  typedef typename R::monom_t       monom_t;
  typedef typename monom_t::coeff_t coeff_t;
  q= R(0);
  if(a==0) return;

  R b(b0), a0(a);
  coeff_t cb =b0.begin()->coeff();
  //  b/=cb;
  monom_t mb = (*b.begin()),
          m  = (*a.begin())/mb,
          null =  monom_t(0);
  R t;
  //  cout <<">>> "<<a<<"\n||| "<<b<<endl;
  while(m != null && a !=0)
    {
      typename R::const_iterator ia =a.begin();
      while(ia !=a.end() && (m =(*ia/mb))== null ) ia++;
      t=b*m;
      a-= t;
      if ( (m) != 0 )
        q+= (m);
  }
//   cout <<"============================================================="<<endl;
//   cout <<q<<endl;
//   cout <<"============================================================="<<endl;
}

//----------------------------------------------------------------------
template<class R>
void diff(R & r, const R & p, int i)
{
    assert(&r != &p);
    typedef typename R::monom_t monom_t;
    r= R(0);
    for (typename R::const_iterator it = p.begin(); it != p.end(); ++it)
      {
        if((*it)[i]>1){
          monom_t m(*it);
          m*=(typename R::coeff_t)(*it);
          m.setexpt(i,(*it)[i]-1);
          r+=m;
        }
        else if ((*it)[i]==1)
          {
            monom_t m(it->coeff());
            for( int ind=0; ind<(it->nvars()+1); ++ind)
              {
                if( ((*it)[ind] !=0) && (ind != i) )
                  m*=monom_t(ind,(*it)[ind]);
              }
            r+=m;
          }
      }
}

//----------------------------------------------------------------------
template <class R>
void copy(R& r, const R& a)
{
  //to nbe modified when bug fixed in gcc
  //  using VECTOR::reserve;
  using namespace VECTOR;
  reserve(r,a.size());
  std::copy(a.begin(),a.end(),r.begin());
}
//----------------------------------------------------------------------
template <class R>
typename R::coeff_t eval(const R& p,
                         const typename R::coeff_t &x,
                         const typename R::coeff_t &y)
{
  typename R::coeff_t r(0);
  for(typename R::const_iterator it=p.begin(); it!=p.end(); ++it)
    {
      typename R::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 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);
      int puissy=(*i);
      int d=(*i);
      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 R, class VARIABLES>
OS & print(OS & os, const R & P, const VARIABLES & V)
{
  if ( P.begin() == P.end() )
    {
      os << '0';
      return os;
    }
 typename R::const_iterator i = P.begin();
 while ( i != P.end() )
   {
     if ( with_plus_sign(i->coeff()) && i !=P.begin())
       os <<'+';
     MONOMIAL::print(os,*i++,V);
   }
 return os;
}
//--------------------------------------------------------------------
template <class OS,class R>
OS & print(OS & os, const R & P)
{
  return print(os,P,Variables::default_);
}
//----------------------------------------------------------------------
template <class OS,class R>
OS & print_verbatim(OS & os, const R & P)
{
  if ( P.begin() == P.end() )
    {
      os << '0';
      return os;
    }
 for( typename R::const_iterator i=P.begin(); i != P.end();++i)
   {
     os<<" "<<i->coeff()<<" [";
     for(int l=0;l<i->size();l++) os<<(*i)<<" ";
     os<<"]";
   }
 return os;
}
//----------------------------------------------------------------------
template<class MP>
MP Rotate_yz(const MP & P,typename MP::coeff_t t)
{
  typedef typename MP::const_iterator  const_iterator;
  typedef typename MP::coeff_t coeff_t;

  coeff_t cos_teta = (1-pow<coeff_t>(t,2))/(1+pow<coeff_t>(t,2));
  coeff_t sin_teta = 2*t/(1+pow<coeff_t>(t,2));
  coeff_t coeff;

  MP pi1,pi2,pi3,P5;

  for(const_iterator i = P.begin(); i != P.end(); ++i)
    {
      int puiss1=(*i);
      int puiss2=(*i);
      int puiss3=(*i);

      coeff = i->coeff();
      MP P1,P2,P3;

      if (puiss1)
        {
          pi1 = Monom<coeff_t,dynamicexp<unsigned> >(1,0,puiss1);
          if ( Degree(pi1) > 0 )
            P1+=pi1;
        }

      for(int j=0;j<=puiss2;j++)
        {
          int puiss_j = puiss2-j;
          coeff_t coeffB=binomial<coeff_t>(puiss2,j);
          coeff_t k=pow<coeff_t>(cos_teta,j);
          pi2 = Monom<coeff_t,dynamicexp<unsigned> >(coeffB*k,1,j);
          k=pow<coeff_t>(-sin_teta,puiss_j);
          pi3 = Monom<coeff_t,dynamicexp<unsigned> >(k,2,puiss_j);
          if ( Degree(pi2*pi3) > 0 )
            P2 += (pi2*pi3);
        }

      for(int j=0;j<=puiss3;j++)
        {
          int puiss_j = puiss3-j;
          coeff_t coeffB=binomial<coeff_t>(puiss3,j);
          coeff_t k=pow<coeff_t>(sin_teta,j);
          pi2 = Monom<coeff_t,dynamicexp<unsigned> >(coeffB*k,1,j);
          k=pow<coeff_t>(cos_teta,puiss_j);
          pi3 = Monom<coeff_t,dynamicexp<unsigned> >(k,2,puiss_j);
          if ( Degree(pi2*pi3) > 0)
            P3 += (pi2*pi3);
        }

      MP P4("0;");
      if (coeff != 0)
        P4 += coeff;
      if (puiss1)
        P4 *= P1;
      if (puiss2)
        P4 *= P2;
      if (puiss3)
        P4 *= P3;

      if (P4 != MP("0;"))
        P5 += P4;
    }
  return P5;
}

//----------------------------------------------------------------------
template<class MP>
MP Rotate_xz(const MP & P,typename MP::coeff_t t)
{
  typedef typename MP::const_iterator  const_iterator;
  typedef typename MP::coeff_t coeff_t;

  coeff_t cos_teta = (1-pow<coeff_t>(t,2))/(1+pow<coeff_t>(t,2));
  coeff_t sin_teta = 2*t/(1+pow<coeff_t>(t,2));
  coeff_t coeff;

  MP pi1,pi2,pi3,P5;

  for(const_iterator i = P.begin(); i != P.end(); ++i)
    {
      int puiss1=(*i);
      int puiss2=(*i);
      int puiss3=(*i);

      coeff = i->coeff();
      MP P1,P2,P3;

      for(int j=0;j<=puiss1;j++)
        {
          int puiss_j = puiss1-j;
          coeff_t coeffB=binomial<coeff_t>(puiss1,j);
          coeff_t k=pow<coeff_t>(cos_teta,j);
          pi1 = Monom<coeff_t,dynamicexp<unsigned> >(coeffB*k,0,j);
          k=pow<coeff_t>(-sin_teta,puiss_j);
          pi3 = Monom<coeff_t,dynamicexp<unsigned> >(k,2,puiss_j);
          if ( Degree(pi1*pi3) > 0)
            P1 += (pi1*pi3);
        }

      if (puiss2)
        {
          pi2 = Monom<coeff_t,dynamicexp<unsigned> >(1,1,puiss2);
          if ( Degree(pi2) > 0 )
            P2+=pi2;
        }

      for(int j=0;j<=puiss3;j++)
        {
          int puiss_j = puiss3-j;
          coeff_t coeffB=binomial<coeff_t>(puiss3,j);
          coeff_t k=pow<coeff_t>(sin_teta,j);
          pi1 = Monom<coeff_t,dynamicexp<unsigned> >(coeffB*k,0,j);
          k=pow<coeff_t>(cos_teta,puiss_j);
          pi3 = Monom<coeff_t,dynamicexp<unsigned> >(k,2,puiss_j);
          if ( Degree(pi1*pi3) > 0)
            P3 += (pi1*pi3);
        }

      MP P4("0;");
      if (coeff != 0)
        P4 += coeff;
      if (puiss1)
        P4 *= P1;
      if (puiss2)
        P4 *= P2;
      if (puiss3)
        P4 *= P3;

      if (P4 != MP("0;"))
        P5 += P4;
    }
  return P5;
}

//----------------------------------------------------------------------
template<class MP>
MP Rotate_xy(const MP & P,typename MP::coeff_t t)
{
  typedef typename MP::const_iterator  const_iterator;
  typedef typename MP::coeff_t coeff_t;

  coeff_t cos_teta = (1-pow<coeff_t>(t,2))/(1+pow<coeff_t>(t,2));
  coeff_t sin_teta = 2*t/(1+pow<coeff_t>(t,2));
  coeff_t coeff;

  MP pi1,pi2,pi3,P5;

  for(const_iterator i = P.begin(); i != P.end(); ++i)
    {
      int puiss1=(*i);
      int puiss2=(*i);
      int puiss3=(*i);

      coeff = i->coeff();
      MP P1,P2,P3;

      for(int j=0;j<=puiss1;j++)
        {
          int puiss_j = puiss1-j;
          coeff_t coeffB=binomial<coeff_t>(puiss1,j);
          coeff_t k=pow<coeff_t>(cos_teta,j);
          pi1 = Monom<coeff_t,dynamicexp<unsigned> >(coeffB*k,0,j);
          k=pow<coeff_t>(-sin_teta,puiss_j);
          pi2 = Monom<coeff_t,dynamicexp<unsigned> >(k,1,puiss_j);
          if ( Degree(pi1*pi2) > 0)
            P1 += (pi1*pi2);
        }

      for(int j=0;j<=puiss2;j++)
        {
          int puiss_j = puiss2-j;
          coeff_t coeffB=binomial<coeff_t>(puiss2,j);
          coeff_t k=pow<coeff_t>(sin_teta,j);
          pi1 = Monom<coeff_t,dynamicexp<unsigned> >(coeffB*k,0,j);
          k=pow<coeff_t>(cos_teta,puiss_j);
          pi2 = Monom<coeff_t,dynamicexp<unsigned> >(k,1,puiss_j);
          if ( Degree(pi1*pi2) > 0)
            P2 += (pi1*pi2);
        }

      if (puiss3)
        {
          pi3 = Monom<coeff_t,dynamicexp<unsigned> >(1,2,puiss3);
          if ( Degree(pi3) > 0 )
            P3+=pi3;
        }

      MP P4("0;");
      if (coeff != 0)
        P4 += coeff;
      if (puiss1)
        P4 *= P1;
      if (puiss2)
        P4 *= P2;
      if (puiss3)
        P4 *= P3;

      if (P4 != MP("0;"))
        P5 += P4;
    }
  return P5;
}

//----------------------------------------------------------------------
template<class MP>
MP Rotate_xyz(const MP & P,typename MP::coeff_t t)
{
  typedef typename MP::const_iterator  const_iterator;
  typedef typename MP::coeff_t coeff_t;

  MP P1,P2,P3;

  P1 = Rotate_xy(P,t);
  P2 = Rotate_xz(P1,t);
  P3 = Rotate_yz(P2,t);
  return P3;
}

//----------------------------------------------------------------------
template <class POL, class C>
POL shift(typename POL::const_iterator monom, C a, int i)
{

  POL polynom;
  // Definition du type monome.
  typedef typename POL::monom_t Monom;
  // Monome temporaire.
  Monom tmpMonom;

  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;
      expt = 0;
      expt = monom->operator(i1);
      expt = monom->operator(i2);

      if ( ( expt[i0] == 0 ) && ( expt[i1] == 0 ) && ( expt[i2] == 0 ) )
        tmpMonom = Monom( monom->coeff(), 0, 0 );
      else   tmpMonom = Monom( monom->coeff(), 3, expt );

      polynom *= tmpMonom;
    }
  // Si a = 0 alors le monome reste inchange.
  else polynom = ( *monom );

  return polynom;

}
//---------------------------------------------------------------------
template < class POL, class C>
POL shift(const POL &polynom, C x0, C y0, C z0 ) {

  // Definition du type monome.
  typedef typename POL::monom_t Monom;

  // Les polynomes temporaires.
  POL tmpPolynom1, tmpPolynom2, tmpPolynom3;

  // L'iterateur va pointer sur le monome courrant.
  typename POL::const_iterator currentMonom = polynom.begin();

  if ( currentMonom == polynom.end() ) {
    std::cerr << "Error : empty polynom ... in Shift ";
    exit(1);
  }

  // Pour la premiere variable.
  while ( currentMonom != polynom.end() ) {
    // Changement de variable sur le monome courant.
    tmpPolynom1 += MPOLDST::shift<POL,C>(currentMonom, x0, 0);

    // Le prochain Monome.
    ( *currentMonom++ );
  }

  if ( tmpPolynom1 == 0 ) tmpPolynom1= polynom ;

  currentMonom = tmpPolynom1.begin();

  // Pour la seconde variable.
  while ( currentMonom != tmpPolynom1.end() ) {
    tmpPolynom2 += shift<POL,C>(currentMonom, y0, 1);
    ( *currentMonom++ );
  }

  if ( tmpPolynom2 == 0 ) tmpPolynom2 =polynom;
  currentMonom = tmpPolynom2.begin();

  // Pour la 3ieme variable.
  while ( currentMonom != tmpPolynom2.end() ) {
    tmpPolynom3 += shift<POL,C>(currentMonom, z0, 2);
    ( *currentMonom++ );
  }

  return  tmpPolynom3;
}

//----------------------------------------------------------------------
template < class POL, class C >
POL scale( const POL &polynom, C a, C b, C c ) {

  typedef typename POL::monom_t Monom;

  C newCoeff;
  POL  newPolynom;

  for ( typename POL::const_iterator i = polynom.begin();
        i != polynom.end(); ++i ) {

    Monom tmpMonom;

    newCoeff = i->coeff()
      * pow( a, int( i->operator[]( 0 ) ) )
      * pow( b, int( i->operator[]( 1 ) ) )
      * pow( c, int( i->operator[]( 2 ) ) );

    int expt;
    expt = i->operator(0);
    expt = i->operator(1);
    expt = i->operator(2);

    if ( ( expt[0] == 0 ) && ( expt[1] == 0 ) && ( expt[2] == 0 ) )
      tmpMonom = Monom( i->coeff(), 0, 0 );
    else   tmpMonom = Monom( newCoeff, 3, expt );

    newPolynom += tmpMonom;
  }
  return newPolynom;

}

//----------------------------------------------------------------------
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::monom_t 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() = 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::monom_t monom_t;
  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);
  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 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)
        for(unsigned  k=0;k<(*it);k++)
          {
            c*=v[i];
          }
      r+=c;
    }
  return r;
}

//======================================================================
} //CLOSE MPOLDST
__END_NAMESPACE_SYNAPS

#endif // SYNAPS_MPOL_MPOLDST_H