synaps/mpol/Monom.h

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*      Author(s): B. Mourrain, GALAAD, INRIA                         *
**********************************************************************
$Id: Monom.h,v 1.1 2005/07/11 11:17:56 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_MPOL_MONOM_H
#define SYNAPS_MPOL_MONOM_H
//--------------------------------------------------------------------
#include <assert.h>
#include <map>
#include <string>
#include <synaps/init.h>
#include <synaps/base/COUNT.h>
#include <synaps/base/shared_object.h>
#include <synaps/base/type.h>
#include <synaps/linalg/VECTOR.m>
#include <synaps/arithm/REF.h>
#include <synaps/mpol/MONOMIAL.m>
#include <synaps/mpol/dynamicexp.h>
#include <synaps/mpol/Variables.h>

__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class C, class R=dynamicexp<unsigned> >
struct Monom {

  //data
  C                 _coef;
#ifdef SHARED_MONOM
  shared_object<R>  _expt;
  R &       rep()       {return (*_expt);}
  const R & rep() const {return (*_expt);}
  typedef  shared_object<R> exp_t;
#else
  R  _expt;
  R &       rep()       {return (_expt);}
  const R & rep() const {return (_expt);}
  typedef  R exp_t;
#endif

  static Variables var;
  //   //  unsigned                  _number_var;
  //   //  static std::map<std::string,int> _index_of_var;
  //   //  static std::map<int,std::string> _var_of_index;
  //   static  int index_of_var(const std::string & s) { return var[s]; }
  //   static std::string var_of_index(int i)          { return var[i]; }

  // Type definitions
  typedef int                    index_t;
  typedef R                      container_t;
  typedef typename R::exponent_t exponent_t; //typename R::coeff_t exponent_t;
  typedef typename R::size_type  size_type;
  typedef C                      coeff_t;
  typedef Monom<C,R>             self_t;

  // Constructors
  Monom(): _coef(),_expt()
    {
      rep().resize(0); //COUNT<self_t>(' ');
    }
  Monom(const C & c, const exp_t & r):_coef(c), _expt(r){}
  Monom(const C & c): _coef(c),_expt()
    {
      rep().resize(0);
    } // monome constant
  Monom(int i):_coef(i), _expt()
    {
      rep().resize(0);
    }
  Monom(int i, int d);                 // monome x[i]**d
  Monom(const std::string & s, int d, Variables& vars = Variables::default_ );                // monome x[i]**d
  Monom(const C & c, unsigned s, const type::AsSize & S):_coef(c),_expt(s){}
  Monom(const C & c, int i, int d =1); // monome c*x[i]**d
  Monom(const C & c, int s, int *t):_coef(c), _expt(s){
    for(int i=0;i<s;i++) rep()[i]= t[i];
  }

  //Copy constructor
  Monom(const Monom & m):_coef(m._coef),_expt(m._expt)
    {
      // COUNT<self_t>('c');
    }

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

  // Destructor
  // not dynamic allocation at this level

  // Assignement
  self_t & operator = (const self_t & m)
    {
      //      COUNT<self_t>('=');
      _coef = m.coeff(); 
      _expt = m._expt; 
      return *this;
    };
  template<class S>
  self_t & operator = (const VAL<S> & M)
    {
      assign(*this,M.rep()); return *this;
    }

  // Accessors
  const C & coeff () const { return _coef; }
  C& coeff ()       { return _coef; }

  void setcoeff (const C & c) { _coef = c; }

  unsigned size() const     { return rep().size(); }
  int      nvars() const    { return lvar(rep()); }
  int      nbvar() const    { return lvar(rep())+1; }

  exponent_t operator[] (size_type i) const { return rep()[i]; }

  self_t & setexpt (size_type i, exponent_t d)
    {rep().setExponent(i,d);return *this; }

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


  friend bool operator == (const self_t & m1, const self_t & m2) {
    return m1._coef == m2._coef && m1.rep() == m2.rep();
  }

  friend bool  operator != (const self_t & m1, const self_t & m2) {
    return !(m1 == m2);
  }

  friend bool IsComparable (const self_t & m1, const self_t & m2){
    return m1.rep() == m2.rep();
  }

  // Arithmetic operators
  self_t   operator  - () const {return self_t(-_coef,rep());}
  self_t & operator += (const self_t &);
  self_t & operator -= (const self_t &);
  self_t & operator *= (const self_t &);
  self_t & operator *= (const C &);
  self_t & operator /= (const C &);

};
//======================================================================
template<class R,class O> Variables Monom<R,O>::var;
//----------------------------------------------------------------------//
// template<class R,class O> unsigned Monom<R,O>::_number_var=1;
// template<class R,class O> std::map<std::string,int> Monom<R,O>::_index_of_var;
// template<class R,class O> std::map<int,std::string> Monom<R,O>::_var_of_index;
//----------------------------------------------------------------------//
template <class R,class O>
inline VAL<OP<'*',Monom<R,O>,Monom<R,O> > >
  operator*(const Monom<R,O> & a, const Monom<R,O> & b)
{
  return VAL<OP<'*',Monom<R,O>,Monom<R,O> > >
   (OP<'*',Monom<R,O>,Monom<R,O> >(a,b));
}
//----------------------------------------------------------------------//
template <class R,class O>
inline VAL<OP<'/',Monom<R,O>,Monom<R,O> > >
  operator/(const Monom<R,O> & a, const Monom<R,O> & b)
{
  return VAL<OP<'/',Monom<R,O>,Monom<R,O> > >
   (OP<'/',Monom<R,O>,Monom<R,O> >(a,b));
}

//----------------------------------------------------------------------//
namespace let 
{
  template<class C, class R> inline
  void assign( C & r, const Monom<C,R>& m )
  {
    r = m.coeff();
  };
  template<class C, class R> inline
  const C&  convert( const Monom<C,R> & m, const type::As<C>& ) { return m.coeff(); };
};
//----------------------------------------------------------------------//
template <class C,class R> inline
void assign(Monom<C,R> & a, const OP<'*',Monom<C,R>,Monom<C,R> > & M)
{
  a.coeff() = (M.op1.coeff()*M.op2.coeff());
  if (a.coeff() !=0)
    {
      a.rep().reserve(std::max(M.op1.rep().size(),M.op2.rep().size()));
      using namespace VECTOR; add(a.rep(),M.op1.rep(),M.op2.rep());
    }
  else
    {
      erase(a.rep());
    }
}
//----------------------------------------------------------------------//
template <class R,class O> inline
void assign(Monom<R,O> & res, const OP<'/',Monom<R,O>,Monom<R,O> > & M)
{
  typedef typename Monom<R,O>::exponent_t exponent_t;
  if(M.op1.nvars()< M.op2.nvars())
    res= Monom<R,O>(0);
  else {
    res= Monom<R,O>(M.op1);
    res.coeff() /=M.op2.coeff();
    for (int i = 0; i <= M.op1.nvars() && res.coeff() !=0 ; ++i)
      if(M.op1[i] <M.op2[i])
        {
          res= Monom<R,O>(0);
        }
      else
        res.rep().setExponent(i,exponent_t(M.op1[i] - M.op2[i]));
  }
}
//----------------------------------------------------------------------
template <class C, class R> inline
typename R::degree_t degree(const Monom<C,R> & m)
{
  return degree(m.rep());
}
//----------------------------------------------------------------------//
template<class M>
M MGcd(const M & m1, const M & m2) {
  int d = Gcd(m1.coeff(),m2.coeff());
  M res(d);
  int v = Min(m1.nvars(),m2.nvars());
  for (int i = 0; i <= v; ++i)
    res.rep().setExponent(i, Min(m1[i],m2[i]));
  return res;
}

template <class M>
M div(const M & m1, const M & m2) {
  M res(1); res*=(m1.coeff()/m2.coeff());
  for (int i = 0; i <= m1.nvars(); ++i)
    res.rep().setExponent(i,m1[i] - m2[i]);
  return res;
}

template <class M>
M div_quo(const M & m1, const M & m2) {
  M res(1); res*=(m1.coeff()/m2.coeff());
  for (int i = 0; i <= m1.nvars(); ++i)
    res.rep().setExponent(i,m1[i] - m2[i]);
  return res;
}



template <class C, class R>
Monom<C,R>::Monom(int i, int d):_coef(1)
{
  if(d == 0)
    rep().resize(0);
  else
    {
      rep().reserve(i+1);
      for(int j=0;j<i;j++) rep().setExponent(j,0);
      rep().setExponent(i,d);
    }
}
//----------------------------------------------------------------------

template <class C, class R>
Monom<C,R>::Monom(const std::string & s, int d, 
                  Variables& vars                 ):_coef(1)
{
  if(d == 0)
    rep().resize(0);
 else
     {
       int i=vars[s];
       //  std::cout <<":: "<<s<<" "<<i<<std::endl;
       rep().reserve(i+1);
       for(int j=0;j<i;j++) rep().setExponent(j,0);
       rep().setExponent(i,d);
     }
 }

//----------------------------------------------------------------------
template <class C, class R>
Monom<C,R>::Monom(const C & c, int i, int d ):_coef(c)
{
  if(d == 0)
    rep().resize(0);
  else
    {
      rep().reserve(i+1);
      for(int j=0;j<i;j++) rep().setExponent(j,0);
      rep().setExponent(i,d);
    }
}
//----------------------------------------------------------------------
// arithemtic operators
//----------------------------------------------------------------------
template <class C, class R> inline
bool Monom<C,R>::operator== (int n) const
{
  return (coeff()==n);
}
//----------------------------------------------------------------------
// template <class C, class R> inline
// int Monom<C,R>::index_of_var(const std::string & s)
// {
//   //cout <<">> "<<s<<endl;
//   if(s[0]=='x' && s.size()>1 )
//     {
//       const char* var=s.data()+1;
//       return atoi(var);
//     }
//   else
//     {
//       int i=_index_of_var[s]-1;
//       if(i<0)
//      {
//        while(_var_of_index[_number_var]!="")
//          {
//            _number_var++;
//          }
//        i=_number_var-1;
//        _index_of_var[s]=_number_var;
//        _var_of_index[_number_var]=s;
//      }
//       return i;
//     }
// }
// //----------------------------------------------------------------------
// template <class C, class R> inline
// std::string Monom<C,R>::var_of_index(int i)
// {
//   if(_var_of_index[i+1]!="")
//     return _var_of_index[i+1];
//   else
//     {
//       int sz=1,m=1;
//       while(m<i) {m*=10; sz++;}
//       char str[sz];
//       sprintf(str,"%d",i);
//       return std::string("x")+str;
//     }
// }
//----------------------------------------------------------------------
// arithemtic operators
//----------------------------------------------------------------------
template <class C, class R> inline
Monom<C,R> & Monom<C,R>::operator += (const Monom<C,R> & m)
{
  coeff() += m.coeff();
  if (coeff()==0) erase(rep());
  return *this;
}

template <class C, class R> inline
Monom<C,R> & Monom<C,R>::operator -= (const Monom<C,R> & m)
{
  coeff() -= m.coeff();
  if (coeff()==0) erase(rep());
  return *this;
}

template <class C, class R> inline
Monom<C,R>  operator +  (const Monom<C,R> & m1, const Monom<C,R> & m2)
{
  Monom<C,R> res(m1);
  res.coeff() += m2.coeff();
  return (*res);
}

template <class C, class R> inline
Monom<C,R> operator - (const Monom<C,R> & m1, const Monom<C,R> & m2)
{
        Monom<C,R> res(m1);
        res.coeff() -= m2.coeff();
        return (res);
}

template <class C, class R> inline
Monom<C,R> & Monom<C,R>::operator *= (const Monom<C,R> & m)
{
  using namespace VECTOR;
  coeff() *= m.coeff();
  if (coeff()!=0)
    add(rep(),m.rep());
  else
    erase(rep());
  return *this;
}

template <class C, class R> inline
Monom<C,R> & Monom<C,R>::operator *= (const C & c)
{
  coeff() *= c;
  if (coeff()==0) erase(rep());
  return *this;
}

template <class C, class R> inline
Monom<C,R> & Monom<C,R>::operator /= (const C & c)
{
  assert(c !=0);
  coeff() /= c;
  if (coeff()==0) erase(rep());
  return *this;
}

template <class C, class R> inline
Monom<C,R> operator * (const Monom<C,R> & m1, const C & c2)
{
  Monom<C,R> res(m1);
  res.coeff() *= c2;
  if (res.coeff() ==0) erase(res.rep());
  return (res);
}
template <class C, class R> inline
Monom<C,R> operator * (const C & c2, const Monom<C,R> & m1)
{
  return(m1*c2);
}

template <class C, class R> inline
Monom<C,R> operator / (const Monom<C,R> & m1, const C & c2)
{
  Monom<C,R> res(m1);
  res.coeff() /= c2;
  if (res.coeff() ==0) erase(res.rep());
  return (res);
}


//----------------------------------------------------------------------
// operateurs d'entrees/sorties
//----------------------------------------------------------------------
template <class C, class R> inline
std::ostream & operator << (std::ostream & os, const Monom <C,R> & m)
{
  if (m.coeff() ==0 )    os<<"(0)";
  else {
    MONOMIAL::print(os,m);
  }
  return os;
}
//----------------------------------------------------------------------
__END_NAMESPACE_SYNAPS
//----------------------------------------------------------------------
#endif // SYNAPS_MPOL_MONOM_H