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, const variables& vars);
    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, 
                const 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
//----------------------------------------------------------------------
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) {
    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 C, class E> std::string
to_string_pp (const Monomial & m, const variables& V) {
    std::string s;
    bool first=true;
    for (unsigned i = 0; i < (unsigned)m.rep().size(); i++) {
        if (m[i] !=0) {
            if(first) {
                first =false;
            } else {
                s+="*";
            }
            //s += "x"; s+= to_string(i);
            //std::cout<<":: "<<i<<" "<< V[i]<<std::endl;
            s += V[i];

            if (m[i] != 1) {
                s+= "^";
                s+= to_string(m[i]);
            }
        }
    }
    return s;
}
//--------------------------------------------------------------------
template <class C, class E> inline std::string
to_string (const Monomial& m, const variables& V) {
    std::string s;
    if (m.coeff() == C(0) )
        s+="0";
    else {
        int i=m.nvars();
        while(i> -1 && m[i]==0) i--;
        if(i<0)
            s += to_string(m.coeff());
        else {
            if(m.coeff() != (C)1) {
                if(m.coeff()==(C)-1)
                { s += "-"; }
                else
                {
                    std::string sm;
                    sm = to_string(m.coeff());
                    //bool is_pol = is_polynomial(sm);
                    //if (is_pol)
                        //s +="+(";
                    s+=sm;
                    //if (is_pol)
                        //s +=")";
                    s +="*";
                }
            }
            s += to_string_pp(m,V);
        }
    }
    return s;
}


template <class OSTREAM, class C, class E> inline OSTREAM & 
print_as_double (OSTREAM & os, const Monomial & m, 
                 const variables& V) {
    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,variables());

    return os;
}

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