realroot: include/realroot/monomial.hpp Source File

00001 /*********************************************************************
00002 *      This file is part of the source code of realroot kernel.        *
00003 *      Author(s): B. Mourrain, GALAAD, INRIA                         *
00004 **********************************************************************
00005 $Id: monomial.hpp,v 1.1 2005/07/11 11:17:56 mourrain Exp $
00006 **********************************************************************/
00007 #ifndef realroot_monomial_hpp
00008 #define realroot_monomial_hpp
00009 //======================================================================
00042 //======================================================================
00043 #include <assert.h>
00044 #include <map>
00045 #include <string>
00046 #include <cstring>
00047 #include <realroot/shared_object.hpp>
00048 #include <realroot/parser.hpp>
00049 #include <realroot/array.hpp>
00050 #include <realroot/dynamicexp.hpp>
00051 #include <realroot/variables.hpp>
00052 #include <realroot/print.hpp>
00053 //======================================================================
00054 
00055 #define TMPL template<class C, class E>
00056 #define Monomial monom<C,E>
00057 
00058 namespace mmx {
00059 
00061 template <class C, class E=dynamic_exp<int> >
00062 struct monom {
00063 
00064   //data
00065   C                 _coef;
00066 #ifdef SHARED_MONOM
00067   shared_object<E>  _expt;
00068   E &       rep()       {return (*_expt);}
00069   const E & rep() const {return (*_expt);}
00070   typedef  shared_object<E> exp_t;
00071 #else
00072   E  _expt;
00073   E &       rep()       {return (_expt);}
00074   const E & rep() const {return (_expt);}
00075   typedef  E exp_t;
00076 #endif
00077 
00078   static variables var;
00079 
00080   // Type definitions
00081   typedef int                    index_t;
00082   typedef E                      container_t;
00083   typedef typename E::exponent_t exponent_t; 
00084   //typename E::coeff_t exponent_t;
00085   typedef typename E::size_type  size_type;
00086   typedef C                      coeff_t;
00087   typedef C                      Scalar;
00088   typedef Monomial               self_t;
00089 
00090   // Constructors
00091   monom(): _coef(),_expt() {
00092     rep().resize(0); //COUNT<self_t>(' ');
00093   }
00094 
00095   monom(const C & c, const exp_t & r):_coef(c), _expt(r){}
00096   monom(const C & c): _coef(c),_expt() {
00097     rep().resize(0);
00098   } 
00099   monom(int i):_coef(i), _expt() {
00100     rep().resize(0);
00101   }
00103   monom(const std::string & s, int d, variables& vars = variables::default_ );  
00104   monom(const C & c, unsigned s, const AsSize & ):_coef(c),_expt(s){}
00105 
00107   monom(const C & c, int d, unsigned i); 
00108   monom(const C & c, int s, int *t):_coef(c), _expt(s){
00109     for(int i=0;i<s;i++) rep()[i]= t[i];
00110   }
00111   monom(const char *, variables& vars =  monom::var);
00112 
00113   //Copy constructor
00114   monom(const monom & m):_coef(m._coef),_expt(m._expt) {}
00115 
00116   // Destructor
00117   // no dynamic allocation at this level
00118 
00119   // Assignement
00120   self_t & operator = (const self_t & m)
00121     {
00122       _coef = m.coeff(); 
00123       _expt = m._expt; 
00124       return *this;
00125     };
00126 
00127   // Accessors
00128   const C & coeff () const     { return _coef; }
00129   C&        coeff ()           { return _coef; }
00130   void  set_coeff (const C & c){ _coef = c; }
00131 
00132   int*     begin()             { return rep().begin(); }
00133   int*     begin() const       { return rep().begin(); }
00134 
00135   unsigned size()  const       { return rep().size(); }
00136   int      nvars() const       { return lvar(rep()); }
00137   int      l_variable() const  { return lvar(rep()); }
00138   unsigned nbvar() const       { return lvar(rep())+1; }
00139 
00140   exponent_t operator[] (size_type i) const { return rep()[i]; }
00141   exponent_t expt       (size_type i) const { return rep()[i]; }
00142   exponent_t exponent   (size_type i) const { return rep()[i]; }
00143 
00144   self_t & set_expt (size_type i, exponent_t d)
00145     {rep().set_expt(i,d);return *this; }
00146 
00147   // boolean operators and functions
00148   bool operator==(int n) const;
00149   bool operator!=(int n) const {return !(*this==n);}
00150 
00151 
00152   friend bool operator == (const self_t & m1, const self_t & m2) 
00153   {
00154     return m1._coef == m2._coef && m1.rep() == m2.rep();
00155   }
00156 
00157   friend bool  operator != (const self_t & m1, const self_t & m2) 
00158   {
00159     return !(m1 == m2);
00160   }
00161 
00162   friend bool IsComparable (const self_t & m1, const self_t & m2)
00163   {
00164     return m1.rep() == m2.rep();
00165   }
00166 
00167   // Arithmetic operators
00168   self_t   operator  - () const {return self_t(-_coef,rep());}
00169   self_t & operator += (const self_t &);
00170   self_t & operator -= (const self_t &);
00171   self_t & operator *= (const self_t &);
00172   self_t & operator *= (const C &);
00173   self_t & operator /= (const C &);
00174 
00175 };
00176 
00177 //======================================================================
00178   template<class C,class E> variables Monomial::var;
00179 //----------------------------------------------------------------------//
00180 // template<class R,class O> unsigned Monomial::_number_var=1;
00181 // template<class R,class O> std::map<std::string,int> Monomial::_index_of_var;
00182 // template<class R,class O> std::map<int,std::string> Monomial::_var_of_index;
00183 //----------------------------------------------------------------------//
00184 namespace let 
00185 {
00186   TMPL inline void 
00187   assign( C & r, const Monomial& m ) { r = m.coeff(); };
00188  //  TMPL inline const C&  
00189 //   convert( const Monomial & m, const meta::As<C>& ) { return m.coeff(); };
00190 };
00191 //----------------------------------------------------------------------//
00192 //#define MONOM  Monomial
00193 //  declare_binary_operator(TMPL,MONOM,MONOM,meta::_mul_,operator*);
00194 //  declare_binary_operator(TMPL,MONOM,MONOM,meta::_div_,operator/);
00195 //#undef MONOM
00196 TMPL inline Monomial 
00197 operator* (const Monomial& m1, const Monomial& m2) {
00198   Monomial r; mul(r,m1,m2); return r;
00199 }
00200 
00201 TMPL inline Monomial 
00202 operator/ (const Monomial& m1, const Monomial& m2) {
00203   Monomial r; div(r,m1,m2); return r;
00204 }
00205 
00206 TMPL inline bool
00207 divide(const Monomial& m1, const Monomial& m2, Monomial& r) {
00208   div(r,m1,m2); 
00209   for (int i = 0; i <= r.nvars(); ++i)
00210     if(r[i]<0) return false;
00211   return true;
00212 }
00213 //----------------------------------------------------------------------//
00214 TMPL inline void 
00215 mul(Monomial & a, const Monomial& m1, const Monomial& m2) {
00216    a.coeff() = (m1.coeff()*m2.coeff());
00217    if (a.coeff() !=(typename Monomial::Scalar)0)
00218      {
00219        a.rep().reserve(std::max(m1.rep().size(),m2.rep().size()));
00220        array::add(a.rep(),m1.rep(),m2.rep());
00221      }
00222    else
00223      {
00224        erase(a.rep());
00225      }
00226 }
00227 //----------------------------------------------------------------------//
00228 TMPL inline void 
00229 div(Monomial & a, const Monomial& m1, const Monomial& m2) {
00230   a.coeff() = (m1.coeff()/m2.coeff());
00231   if (a.coeff() !=0)
00232     {
00233       a.rep().reserve(std::max(m1.rep().size(),m2.rep().size()));
00234       array::sub(a.rep(),m1.rep(),m2.rep());
00235     }
00236   else
00237     {
00238       erase(a.rep());
00239     }
00240 }
00241 //----------------------------------------------------------------------
00242 TMPL inline Monomial 
00243 div(const Monomial & m1, const Monomial & m2) {
00244   Monomial res(1); res*=(m1.coeff()/m2.coeff());
00245   for (int i = 0; i <= m1.nvars(); ++i)
00246     res.set_expt(i,m1[i] - m2[i]);
00247   return res;
00248 }
00249 //----------------------------------------------------------------------
00250 TMPL inline Monomial 
00251 div_quo(const Monomial & m1, const Monomial & m2) {
00252   Monomial res(1); res*=(m1.coeff()/m2.coeff());
00253   for (int i = 0; i <= m1.nvars(); ++i)
00254     res.set_expt(i,m1[i] - m2[i]);
00255   return res;
00256 }
00257 //----------------------------------------------------------------------
00258 template <class C, class E> inline
00259 typename E::degree_t degree(const Monomial & m) {
00260   return degree(m.rep());
00261 }
00262 //----------------------------------------------------------------------//
00263 TMPL Monomial
00264 MGcd(const Monomial & m1, const Monomial & m2) {
00265   int d = Gcd(m1.coeff(),m2.coeff());
00266   Monomial res(d);
00267   int v = Min(m1.nvars(),m2.nvars());
00268   for (int i = 0; i <= v; ++i)
00269     res.set_expt(i, Min(m1[i],m2[i]));
00270   return res;
00271 }
00272 
00273 
00274 TMPL inline Monomial 
00275 pow(const Monomial& m1, int n) {
00276   Monomial res(m1);
00277   for (int i = 0; i <= m1.nvars(); ++i)
00278     res.set_expt(i,m1[i]*n);
00279   return res;
00280 }
00281 //----------------------------------------------------------------------
00283 template <class C, class E>
00284 Monomial::monom(const std::string & s, int d, 
00285                    variables& vars  ):_coef(1) {
00286   if(d == 0)
00287     rep().resize(0);
00288  else
00289      {
00290        int i=vars[s];
00291        rep().reserve(i+1);
00292        for(int j=0;j<i;j++) rep().set_expt(j,0);
00293        rep().set_expt(i,d);
00294      }
00295 }
00296 //----------------------------------------------------------------------
00298 template <class C, class E>
00299 Monomial::monom(const C & c, int d, unsigned v):_coef(c) {
00300   if(d == 0)
00301     rep().resize(0);
00302   else
00303     {
00304       rep().reserve(v+1);
00305       for(unsigned j=0;j<v;j++) rep().set_expt(j,0);
00306       rep().set_expt(v,d);
00307     }
00308 }
00309 //----------------------------------------------------------------------
00310 extern "C" char *synaps_inputptr;
00311 template <class C, class E>
00312 Monomial::monom(const char * s, variables&  vars ) {
00313   int n = strlen(s);
00314   if (s[n-1]=='\n') {n--;}
00315   if(s[n-1]==';')
00316     {
00317       synaps_input = new char[n+1];
00318       memcpy(synaps_input,s,n); synaps_input[n]='\0';
00319       synaps_inputlim = synaps_input + n;
00320     }
00321   else
00322     {
00323       synaps_input = new char[n+2];
00324       memcpy(synaps_input,s,n); 
00325       synaps_input[n]=';'; synaps_input[n+1]='\0';
00326       synaps_inputlim = synaps_input + n+1;
00327     }
00328   synaps_inputptr=synaps_input;
00329   *this = Monomial((C)1);
00330   int Ask;
00331 
00332   for (;;) {
00333     Ask = yylex();
00334     if (Ask == COEF) { /* monome constant */
00335       C c;
00336       let::assign(c,yylval);
00337       free(yylval);
00338       *this *= c;
00339     }
00340     else if (Ask == XID) {
00341       char* p = strchr(yylval,'^');
00342       if (p == NULL)
00343         {
00344           *this *= Monomial(vars[std::string(yylval)],1);
00345         }
00346       else {
00347         p++;
00348         *this *= Monomial(vars[std::string(yylval,p-1)], atoi(p));
00349       }
00350       free(yylval);
00351     }
00352     else if (Ask == TERMINATE) {
00353       break;
00354     }
00355   }
00356   delete[] synaps_input;
00357 }
00358 //----------------------------------------------------------------------
00359 template <class C, class E> inline bool 
00360 Monomial::operator== (int n) const {
00361   return (coeff()==n);
00362 }
00363 //----------------------------------------------------------------------
00364 // arithemtic operators
00365 //----------------------------------------------------------------------
00366 template <class C, class E> inline Monomial & 
00367 Monomial::operator += (const Monomial & m) {
00368   coeff() += m.coeff();
00369   if (coeff()==0) erase(rep());
00370   return *this;
00371 }
00372 
00373 template <class C, class E> inline Monomial & 
00374 Monomial::operator -= (const Monomial & m) {
00375   coeff() -= m.coeff();
00376   if (coeff()==0) erase(rep());
00377   return *this;
00378 }
00379 
00380 template <class C, class E> inline Monomial  
00381 operator +  (const Monomial & m1, const Monomial & m2) {
00382   Monomial res(m1);
00383   res.coeff() += m2.coeff();
00384   return (*res);
00385 }
00386 
00387 template <class C, class E> inline Monomial 
00388 operator - (const Monomial & m1, const Monomial & m2) {
00389         Monomial res(m1);
00390         res.coeff() -= m2.coeff();
00391         return (res);
00392 }
00393 
00394 template <class C, class E> inline Monomial & 
00395 Monomial::operator *= (const Monomial & m) {
00396   coeff() *= m.coeff();
00397   if (coeff()!=0)
00398     array::add(rep(),m.rep());
00399   else
00400     erase(rep());
00401   return *this;
00402 }
00403 
00404 template <class C, class E> inline Monomial & 
00405 Monomial::operator *= (const C & c) {
00406   coeff() *= c;
00407   if (coeff()==0) erase(rep());
00408   return *this;
00409 }
00410 
00411 template <class C, class E> inline Monomial & 
00412 Monomial::operator /= (const C & c) {
00413   assert(c !=0);
00414   coeff() /= c;
00415   if (coeff()==0) erase(rep());
00416   return *this;
00417 }
00418 
00419 template <class C, class E> inline Monomial 
00420 operator * (const Monomial & m1, const C & c2) {
00421   Monomial res(m1);
00422   res.coeff() *= c2;
00423   if (res.coeff() ==0) erase(res.rep());
00424   return (res);
00425 }
00426 template <class C, class E> inline Monomial 
00427 operator * (const C & c2, const Monomial & m1) {
00428   return(m1*c2);
00429 }
00430 
00431 template <class C, class E> inline Monomial 
00432 operator / (const Monomial & m1, const C & c2) {
00433   Monomial res(m1);
00434   res.coeff() /= c2;
00435   if (res.coeff() ==0) erase(res.rep());
00436   return (res);
00437 }
00438 
00439 
00440 //----------------------------------------------------------------------
00441 // operateurs d'entrees/sorties
00442 //----------------------------------------------------------------------
00443 template <class OSTREAM, class C, class E> inline OSTREAM& 
00444 print_pp (OSTREAM & os, const Monomial & m, const variables& V) { 
00445     bool first=true;
00446     for (typename E::size_type i = 0; i < m.rep().size(); ++i)
00447       {
00448         if (m[i] !=0)
00449           {
00450             if(first) {
00451               first =false;
00452             } else {
00453               os<<'*';
00454             }
00455             os << V[i];
00456             if (m[i] != 1) {
00457               os << '^' <<m[i];
00458             }
00459           }
00460       }
00461     return os;
00462   }
00463 //--------------------------------------------------------------------
00464 
00465  inline bool is_polynomial(const std::ostringstream& sm){
00466     for(unsigned i=1; i<sm.str().length();i++)
00467       if(sm.str()[i]=='+' || (sm.str()[i]=='-' && (i>0 && sm.str()[i-1]=='e'))) return true;
00468     return false;
00469   }
00470 
00471 
00472 template <class OSTREAM, class C, class E> inline OSTREAM & 
00473 print (OSTREAM & os, const Monomial & m, 
00474        const variables& V = variables::default_) {
00475   if (m.coeff() == C(0) )
00476     os<<"0";
00477   else {
00478     int i=m.nvars();
00479     while(i> -1 && m[i]==0) i--;
00480     if(i<0)
00481       print(os,m.coeff());
00482     else {
00483       if(m.coeff() != (C)1) {
00484         if(m.coeff()==(C)-1) 
00485           { os << "-"; }
00486         else 
00487           {  
00488             std::ostringstream sm;
00489             print(sm,m.coeff()); 
00490             bool is_pol = is_polynomial(sm);
00491             if (is_pol) os <<"+(";
00492             os<<sm.str();
00493             if (is_pol) os <<")";
00494             os <<"*"; 
00495           }
00496       }
00497       print_pp(os,m,V);
00498     }
00499   }
00500   return os;
00501 }
00502 
00503 
00504 
00505 template <class OSTREAM, class C, class E> inline OSTREAM & 
00506 print_as_double (OSTREAM & os, const Monomial & m, 
00507        const variables& V = variables::default_) {
00508   if (m.coeff() == C(0) )
00509     os<<"0";
00510   else {
00511     int i=m.nvars();
00512     while(i> -1 && m[i]==0) i--;
00513     if(i<0)
00514       print(os,as_double(m.coeff()));
00515     else {
00516       if(m.coeff() != (C)1) {
00517         if(m.coeff()==(C)-1) 
00518           { os << "-"; }
00519         else 
00520           {  
00521             std::ostringstream sm;
00522             print(sm,as_double(m.coeff())); 
00523             bool is_pol = is_polynomial(sm);
00524             if (is_pol) os <<"+(";
00525             os<<sm.str();
00526             if (is_pol) os <<")";
00527             os <<"*"; 
00528           }
00529       }
00530       print_pp(os,m,V);
00531     }
00532   }
00533   return os;
00534 }
00535 
00536 
00537 
00538 
00540 template <class C, class E> inline
00541 std::ostream & operator << (std::ostream & os, const monom <C,E> & m)
00542 {
00543   if (m.coeff() ==0 )    os<<"(0)";
00544   else 
00545     print(os,m);
00546   
00547   return os;
00548 }
00549 
00550 //--------------------------------------------------------------------
00551 } //  namespace mmx 
00552 //======================================================================
00553 #undef TMPL
00554 #undef Monomial
00555 #endif // realroot_monom_hpp
00556