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

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of realroot kernel.        *
*      Author(s): B. Mourrain, GALAAD, INRIA                         *
**********************************************************************
$Id: monomial.hpp,v 1.1 2005/07/11 11:17:56 mourrain Exp $
**********************************************************************/
#ifndef realroot_monomial_hpp
#define realroot_monomial_hpp
//======================================================================
//======================================================================
#include <assert.h>
#include <map>
#include <string>
#include <cstring>
#include <realroot/shared_object.hpp>
#include <realroot/parser.hpp>
#include <realroot/array.hpp>
#include <realroot/dynamicexp.hpp>
#include <realroot/variables.hpp>
#include <realroot/print.hpp>
//======================================================================

#define TMPL template<class C, class E>
#define Monomial monom<C,E>

namespace mmx {

template <class C, class E=dynamic_exp<int> >
struct monom {

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

  static variables var;

  // Type definitions
  typedef int                    index_t;
  typedef E                      container_t;
  typedef typename E::exponent_t exponent_t; 
  //typename E::coeff_t exponent_t;
  typedef typename E::size_type  size_type;
  typedef C                      coeff_t;
  typedef C                      Scalar;
  typedef Monomial               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);
  } 
  monom(int i):_coef(i), _expt() {
    rep().resize(0);
  }
  monom(const std::string & s, int d, variables& vars = variables::default_ );  
  monom(const C & c, unsigned s, const AsSize & ):_coef(c),_expt(s){}

  monom(const C & c, int d, unsigned i); 
  monom(const C & c, int s, int *t):_coef(c), _expt(s){
    for(int i=0;i<s;i++) rep()[i]= t[i];
  }
  monom(const char *, variables& vars =  monom::var);

  //Copy constructor
  monom(const monom & m):_coef(m._coef),_expt(m._expt) {}

  // Destructor
  // no dynamic allocation at this level

  // Assignement
  self_t & operator = (const self_t & m)
    {
      _coef = m.coeff(); 
      _expt = m._expt; 
      return *this;
    };

  // Accessors
  const C & coeff () const     { return _coef; }
  C&        coeff ()           { return _coef; }
  void  set_coeff (const C & c){ _coef = c; }

  int*     begin()             { return rep().begin(); }
  int*     begin() const       { return rep().begin(); }

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

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

  self_t & set_expt (size_type i, exponent_t d)
    {rep().set_expt(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 C,class E> variables Monomial::var;
//----------------------------------------------------------------------//
// template<class R,class O> unsigned Monomial::_number_var=1;
// template<class R,class O> std::map<std::string,int> Monomial::_index_of_var;
// template<class R,class O> std::map<int,std::string> Monomial::_var_of_index;
//----------------------------------------------------------------------//
namespace let 
{
  TMPL inline void 
  assign( C & r, const Monomial& m ) { r = m.coeff(); };
 //  TMPL inline const C&  
//   convert( const Monomial & m, const meta::As<C>& ) { return m.coeff(); };
};
//----------------------------------------------------------------------//
//#define MONOM  Monomial
//  declare_binary_operator(TMPL,MONOM,MONOM,meta::_mul_,operator*);
//  declare_binary_operator(TMPL,MONOM,MONOM,meta::_div_,operator/);
//#undef MONOM
TMPL inline Monomial 
operator* (const Monomial& m1, const Monomial& m2) {
  Monomial r; mul(r,m1,m2); return r;
}

TMPL inline Monomial 
operator/ (const Monomial& m1, const Monomial& m2) {
  Monomial r; div(r,m1,m2); return r;
}

TMPL inline bool
divide(const Monomial& m1, const Monomial& m2, Monomial& r) {
  div(r,m1,m2); 
  for (int i = 0; i <= r.nvars(); ++i)
    if(r[i]<0) return false;
  return true;
}
//----------------------------------------------------------------------//
TMPL inline void 
mul(Monomial & a, const Monomial& m1, const Monomial& m2) {
   a.coeff() = (m1.coeff()*m2.coeff());
   if (a.coeff() !=(typename Monomial::Scalar)0)
     {
       a.rep().reserve(std::max(m1.rep().size(),m2.rep().size()));
       array::add(a.rep(),m1.rep(),m2.rep());
     }
   else
     {
       erase(a.rep());
     }
}
//----------------------------------------------------------------------//
TMPL inline void 
div(Monomial & a, const Monomial& m1, const Monomial& m2) {
  a.coeff() = (m1.coeff()/m2.coeff());
  if (a.coeff() !=0)
    {
      a.rep().reserve(std::max(m1.rep().size(),m2.rep().size()));
      array::sub(a.rep(),m1.rep(),m2.rep());
    }
  else
    {
      erase(a.rep());
    }
}
//----------------------------------------------------------------------
TMPL inline Monomial 
div(const Monomial & m1, const Monomial & m2) {
  Monomial res(1); res*=(m1.coeff()/m2.coeff());
  for (int i = 0; i <= m1.nvars(); ++i)
    res.set_expt(i,m1[i] - m2[i]);
  return res;
}
//----------------------------------------------------------------------
TMPL inline Monomial 
div_quo(const Monomial & m1, const Monomial & m2) {
  Monomial res(1); res*=(m1.coeff()/m2.coeff());
  for (int i = 0; i <= m1.nvars(); ++i)
    res.set_expt(i,m1[i] - m2[i]);
  return res;
}
//----------------------------------------------------------------------
template <class C, class E> inline
typename E::degree_t degree(const Monomial & m) {
  return degree(m.rep());
}
//----------------------------------------------------------------------//
TMPL Monomial
MGcd(const Monomial & m1, const Monomial & m2) {
  int d = Gcd(m1.coeff(),m2.coeff());
  Monomial res(d);
  int v = Min(m1.nvars(),m2.nvars());
  for (int i = 0; i <= v; ++i)
    res.set_expt(i, Min(m1[i],m2[i]));
  return res;
}


TMPL inline Monomial 
pow(const Monomial& m1, int n) {
  Monomial res(m1);
  for (int i = 0; i <= m1.nvars(); ++i)
    res.set_expt(i,m1[i]*n);
  return res;
}
//----------------------------------------------------------------------
template <class C, class E>
Monomial::monom(const std::string & s, int d, 
                   variables& vars  ):_coef(1) {
  if(d == 0)
    rep().resize(0);
 else
     {
       int i=vars[s];
       rep().reserve(i+1);
       for(int j=0;j<i;j++) rep().set_expt(j,0);
       rep().set_expt(i,d);
     }
}
//----------------------------------------------------------------------
template <class C, class E>
Monomial::monom(const C & c, int d, unsigned v):_coef(c) {
  if(d == 0)
    rep().resize(0);
  else
    {
      rep().reserve(v+1);
      for(unsigned j=0;j<v;j++) rep().set_expt(j,0);
      rep().set_expt(v,d);
    }
}
//----------------------------------------------------------------------
extern "C" char *synaps_inputptr;
template <class C, class E>
Monomial::monom(const char * s, variables&  vars ) {
  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;
  *this = Monomial((C)1);
  int Ask;

  for (;;) {
    Ask = yylex();
    if (Ask == COEF) { /* monome constant */
      C c;
      let::assign(c,yylval);
      free(yylval);
      *this *= c;
    }
    else if (Ask == XID) {
      char* p = strchr(yylval,'^');
      if (p == NULL)
        {
          *this *= Monomial(vars[std::string(yylval)],1);
        }
      else {
        p++;
        *this *= Monomial(vars[std::string(yylval,p-1)], atoi(p));
      }
      free(yylval);
    }
    else if (Ask == TERMINATE) {
      break;
    }
  }
  delete[] synaps_input;
}
//----------------------------------------------------------------------
template <class C, class E> inline bool 
Monomial::operator== (int n) const {
  return (coeff()==n);
}
//----------------------------------------------------------------------
// arithemtic operators
//----------------------------------------------------------------------
template <class C, class E> inline Monomial & 
Monomial::operator += (const Monomial & m) {
  coeff() += m.coeff();
  if (coeff()==0) erase(rep());
  return *this;
}

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

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

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

template <class C, class E> inline Monomial & 
Monomial::operator *= (const Monomial & m) {
  coeff() *= m.coeff();
  if (coeff()!=0)
    array::add(rep(),m.rep());
  else
    erase(rep());
  return *this;
}

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

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

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

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


//----------------------------------------------------------------------
// operateurs d'entrees/sorties
//----------------------------------------------------------------------
template <class OSTREAM, class C, class E> inline OSTREAM& 
print_pp (OSTREAM & os, const Monomial & m, const variables& V) { 
    bool first=true;
    for (typename E::size_type i = 0; i < m.rep().size(); ++i)
      {
        if (m[i] !=0)
          {
            if(first) {
              first =false;
            } else {
              os<<'*';
            }
            os << V[i];
            if (m[i] != 1) {
              os << '^' <<m[i];
            }
          }
      }
    return os;
  }
//--------------------------------------------------------------------

 inline bool is_polynomial(const std::ostringstream& sm){
    for(unsigned i=1; i<sm.str().length();i++)
      if(sm.str()[i]=='+' || (sm.str()[i]=='-' && (i>0 && sm.str()[i-1]=='e'))) return true;
    return false;
  }


template <class OSTREAM, class C, class E> inline OSTREAM & 
print (OSTREAM & os, const Monomial & m, 
       const variables& V = variables::default_) {
  if (m.coeff() == C(0) )
    os<<"0";
  else {
    int i=m.nvars();
    while(i> -1 && m[i]==0) i--;
    if(i<0)
      print(os,m.coeff());
    else {
      if(m.coeff() != (C)1) {
        if(m.coeff()==(C)-1) 
          { os << "-"; }
        else 
          {  
            std::ostringstream sm;
            print(sm,m.coeff()); 
            bool is_pol = is_polynomial(sm);
            if (is_pol) os <<"+(";
            os<<sm.str();
            if (is_pol) os <<")";
            os <<"*"; 
          }
      }
      print_pp(os,m,V);
    }
  }
  return os;
}



template <class OSTREAM, class C, class E> inline OSTREAM & 
print_as_double (OSTREAM & os, const Monomial & m, 
       const variables& V = variables::default_) {
  if (m.coeff() == C(0) )
    os<<"0";
  else {
    int i=m.nvars();
    while(i> -1 && m[i]==0) i--;
    if(i<0)
      print(os,as_double(m.coeff()));
    else {
      if(m.coeff() != (C)1) {
        if(m.coeff()==(C)-1) 
          { os << "-"; }
        else 
          {  
            std::ostringstream sm;
            print(sm,as_double(m.coeff())); 
            bool is_pol = is_polynomial(sm);
            if (is_pol) os <<"+(";
            os<<sm.str();
            if (is_pol) os <<")";
            os <<"*"; 
          }
      }
      print_pp(os,m,V);
    }
  }
  return os;
}




template <class C, class E> inline
std::ostream & operator << (std::ostream & os, const monom <C,E> & m)
{
  if (m.coeff() ==0 )    os<<"(0)";
  else 
    print(os,m);
  
  return os;
}

//--------------------------------------------------------------------
} //  namespace mmx 
//======================================================================
#undef TMPL
#undef Monomial
#endif // realroot_monom_hpp