MPoly.hpp

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of borderbasix           *
*      (C) B. Mourrain, Ph. Trebuchet                                *
**********************************************************************/
#ifndef _MPoly_hpp_
#define _MPoly_hpp_

#include "MPOLYNOMIAL.hpp"
#include <stdlib.h>
#include <string.h>
#include <string>
#include <cassert>
#include "parser.hpp"
#include <borderbasix/arithm/has_sign.hpp>

template <class R, class O> class MPoly {
public:
    // Types
    typedef typename R::value_type monom_t;
    typedef typename monom_t::coeff_t coeff_t;
    typedef O order_t;
    typedef typename R::iterator iterator;
    typedef typename R::const_iterator const_iterator;
    typedef typename R::reverse_iterator reverse_iterator;
    typedef typename R::const_reverse_iterator const_reverse_iterator;

    typedef MPoly<R,O> self_t;

    // Constructors
    MPoly():rep(){}
    MPoly(const monom_t & m):rep() { rep.push_back( m); }
    MPoly(const coeff_t & c):rep() { rep.push_back(monom_t(c));}
    MPoly(int, const monom_t *);
    MPoly(const self_t & P) {
        MPOLYNOMIAL__copy(rep, P.rep);
    }
    MPoly(const int & i): rep(){ rep.push_back(monom_t(i));}
    MPoly(const char *);
    MPoly(char *);

    //  void make();
    // Destructor
    // no dynamic allocation at this level.

    // Assignment operators
    self_t & operator=(const self_t &);
    self_t & operator+=(const monom_t &);
    self_t & operator+=(const self_t &);
    self_t & operator-=(const monom_t &);
    self_t & operator-=(const self_t &);
    self_t & operator*=(const coeff_t &);
    self_t & operator*=(const monom_t &);
    self_t & operator*=(const self_t &);
    self_t & operator/=(const monom_t &);
    self_t & operator/=(const self_t & Q);

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

    //   friend bool operator == <>(const self_t &, const self_t &);
    //   friend bool operator != <>(const self_t &, const self_t &);
    //   friend bool operator < <>(const self_t &, const self_t &);
    //   friend bool operator <= <> (const self_t &, const self_t &);
    //   friend bool operator > <>(const self_t &, const self_t &);
    //   friend bool operator >= <> (const self_t &, const self_t &);

    // Comput operators and functions
    self_t operator-() const;
    //   friend self_t operator + <>(const self_t &, const self_t &);
    //   friend self_t & operator - <>(const self_t &, const self_t &);
    //   friend self_t & operator * <>(const self_t &, const monom_t &);
    //   friend self_t & operator * <>(const self_t &, const self_t &);
    //   friend self_t & operator / <>(const self_t &, const self_t &);

    // Input/Output operator
    //friend template <class R, class O>
    //istream & operator >> (istream & is, MPoly<R,O> & P);
    //istream & operator >> <>(istream &, self_t &);
    //friend template <class R, class O>
    //ostream & operator << (ostream & is, MPoly<R,O> & P);
    //ostream & operator << <>(ostream &, const self_t &);

    // Iterator functions
    iterator begin () { return rep.begin(); }
    const_iterator begin () const { return rep.begin(); }
    reverse_iterator rbegin ()    { return rep.rbegin(); }
    const_reverse_iterator rbegin () const { return rep.rbegin(); }

    iterator end () { return rep.end(); }
    const_iterator end () const { return rep.end(); }
    reverse_iterator rend () { return rep.rend(); }
    const_reverse_iterator rend () const { return rep.rend(); }

    int size() const {return rep.size();}

    // Data
    R rep;

    // private:
    //   // private member functions
    //   iterator add (const self_t &, iterator);
    //   iterator sub (const self_t &, iterator);
};

//----------------------------------------------------------------------
//char *myinput;
//char *myinputptr; 
//char *myinputlim;
//char *yylval;

//extern "C" int yywrap (void);
//extern "C" int yylex(void);

//----------------------------------------------------------------------
// Constructors
//----------------------------------------------------------------------
template <class R, class O>
inline MPoly<R,O>::MPoly (int s, const monom_t *t) {
    for (int i = 0; i<s; i++) *this += t[i];
}

template<class R, class O>
MPoly<R,O>::MPoly(const char* s) {
    int n = strlen(s);
    if(s[n-1]!=';'){
        char* str = new char[n+2];
        memcpy(str,s,n);
        str[n]=';'; str[n+1]='\0';
        yy_scan_string(str);
    } else {
        yy_scan_string(s);
    }
    parser<coeff_t,monom_t,self_t> p;
    p.yyparse();
    if(p.yypol!=0)
        MPOLYNOMIAL__copy(rep,p.yypol->rep);
}

template<class R, class O>
MPoly<R,O>::MPoly(char* s) {
    yy_scan_string(const_cast<const char*>(s));
    parser<coeff_t,monom_t,self_t> p;
    p.yyparse();
    if(p.yypol!=0)
        MPOLYNOMIAL__copy(rep,p.yypol->rep);
}

//----------------------------------------------------------------------
// Assignment operators
//----------------------------------------------------------------------
template <class R, class O>
inline MPoly<R,O>& MPoly<R,O>::operator = (const MPoly<R,O> & P)
{
    //  cout <<"  =[p] "<<P<<endl;
    MPOLYNOMIAL__copy(rep,P.rep);
    return *this;
}

template <class R, class O>
inline MPoly<R,O>& MPoly<R,O>::operator += (const typename MPoly<R,O>::monom_t & m){
    typedef typename MPoly<R,O>::monom_t mon;
    typedef typename mon::coeff_t C;
    if (m.GetCoeff()!=(C)0)
        MPOLYNOMIAL__insert<R,O>(rep, rep.begin(), m);
    return *this;
}

template <class R, class O>
inline MPoly<R,O>& MPoly<R,O>::operator += (const MPoly<R,O> & P)
{
    if (P==0) return *this;
    if (*this==0) {*this = P; return *this;}

    MPOLYNOMIAL__plus<R,O>(rep,P.rep,begin());
    return *this;
}

template <class R, class O>
inline MPoly<R,O>&  MPoly<R,O>::operator -= (const typename MPoly<R,O>::monom_t & m){
    if (m.GetCoeff()!=(typename R::value_type::coeff_t)0){
        MPOLYNOMIAL__insert<R,O>(rep,rep.begin(),-m);
    }
    return *this;
}

template <class R, class O>
inline MPoly<R,O>&  MPoly<R,O>::operator -= (const MPoly<R,O> & P){
    //cout<<"*this"<<*this<<endl;
    //cout<<"P "<<P<<endl;
    if (P==0) return *this;
    if (*this==0) {*this = -P; return *this;}
    MPOLYNOMIAL__sub<R,O>(rep,P.rep,begin());
    return *this;
}

template <class R, class O>
inline MPoly<R,O>& MPoly<R,O>::operator /= (const typename MPoly<R,O>::monom_t & m){
    if ( (*this!=0)  && (m.GetCoeff() !=(typename R::value_type::coeff_t)0) ) {
        MPoly<R,O> P(*this);
        for (iterator i = P.begin(); i != P.end(); i++) {
            int n=m.Lvar();
            if (n<=(*i).Lvar())
                while (((*i).GetDegree(n) >= m.GetDegree(n) ) && ( n >= 0)) n--;
            if (n<0)
                *i=div(*i,m);
            else return *this;
        }
        MPOLYNOMIAL__copy(this->rep,P.rep);
    }
    return *this;
}
template <class R, class O>
inline MPoly<R,O>& MPoly<R,O>::operator *= (const typename MPoly<R,O>::coeff_t & m){
    //xxx IF m<> 0
    MPoly<R,O> P;
    for (iterator i = begin(); i != end(); i++) {
        (*i)*= m;
        //monom_t mon((*i) * m);
        //MPOLYNOMIAL__insert<monom_t,O>(P.rep,P.begin(),mon);
    }
    //xxx Polypack::sort(this->rep,O,IFN<O,MONOMIALORDER>::VAL());
    return *this;
}
template <class R, class O>
inline MPoly<R,O>& MPoly<R,O>::operator *= (const typename MPoly<R,O>::monom_t & m){
    //xxx IF m<> 0
    MPoly<R,O> P;
    for (iterator i = begin(); i != end(); i++) {
        (*i)*= m;
        //monom_t mon((*i) * m);
        //MPOLYNOMIAL__insert<monom_t,O>(P.rep,P.begin(),mon);
    }
    //xxx Polypack::sort(this->rep,O,IFN<O,MONOMIALORDER>::VAL());
    return *this;
}

template <class R, class O>
inline  MPoly<R,O>& MPoly<R,O>::operator *= (const MPoly<R,O> & Q){
    //   cout << "*="<<endl;
    if ( (*this !=0) && (Q !=0) ) {
        bool copy;
        const MPoly<R,O> *NQ;
        copy = (&rep == &Q.rep);
        if (copy)
            NQ = new MPoly<R,O>(Q);
        else
            NQ = &Q;

        iterator i, j;
        const_iterator iP,iQ = NQ->begin();
        monom_t  *m = new monom_t;

        MPoly<R,O> P(*this);
        //    swap(P.rep,this->rep);
        //    cout <<P<<" "<<*this<<endl;
        for(i=begin();i != end();++i) (*i)*=(*iQ);
        //xxx should not imply a copy of m;
        //  this->rep.push_back(m); //xxx use output iterator

        (iQ++);
        i=begin();
        for (; iQ != NQ->end();iQ++) {
            j = begin();
            for (iP = P.begin(); iP != P.end(); iP++) {
                *m = (*iP) * (*iQ);
                j = MPOLYNOMIAL__insert<R,O>(this->rep,j,*m);
            }
        }
        delete m;
        if (copy) delete(NQ);
    }
    else *this=MPoly<R,O>();
    return *this;
}

template <class R, class O>
inline  MPoly<R,O>& MPoly<R,O>::operator /= (const MPoly<R,O> & Q){
    if (*this!=0 && Q!=0){
        bool copy;
        if (Degree(Q) > 0) {
            const MPoly<R,O> *NQ;
            copy = (&rep == &Q.rep);
            if (copy)
                NQ = new MPoly<R,O>(Q);
            else
                NQ = &Q;
            MPoly<R,O> Quo;// XP((*this));
            while (*this!=0) {
                int n=Q.begin()->Lvar();
                if (n<=(*this).begin()->Lvar())
                    while (((*this).begin()->GetDegree(n) >= Q.begin()->GetDegree(n) ) && ( n >= 0)) n--;
                if (n<0) {
                    typename MPoly<R,O>::monom_t m=(div(*((*this).begin()),*(Q.begin())));
                    Quo+=m;
                    (*this)-=Q*m;
                }
                else
                    (*this)-=(*((*this).begin()));
            }
            //Rem = XP-Q*Quo;
            (*this)=Quo;
        }
        else (*this)/= *(Q.begin());
    }
    else
        (*this)=MPoly<R,O>(0);
    return *this;
}

//----------------------------------------------------------------------
// Boolean operators and functions
//----------------------------------------------------------------------
template <class R, class O> inline
bool MPoly<R,O>::operator==(int n) const
{
    if(n==0)
        if (this->rep.empty())
            return true;
        else
            return false;
    else
        if (this->rep.empty())
            return false;
        else
            return( this->rep.size()==1 && (*this - MPoly<R,O>(n)) ==0);
}
template <class R, class O> inline
bool operator == (const MPoly<R,O> & P, const MPoly<R,O> & Q)
{
    return P.rep == Q.rep;
}

template <class R, class O> inline
bool operator != (const MPoly<R,O> & P, const MPoly<R,O> & Q)
{
    return !(P.rep == Q.rep);
}

template <class R, class O> inline
bool operator < (const MPoly<R,O> & P, const MPoly<R,O> & Q) {
    typename MPoly<R,O>::const_iterator iP = P.begin();
    typename MPoly<R,O>::const_iterator iQ = Q.begin();
    while (iP != P.end() && iQ != Q.end() && IsComparable(*iP,*iQ)) {
        ++iP;
        ++iQ;
    }
    if (iQ == Q.end())
        return false;
    if (iP == P.end())
        return true;
    return O::less(*iP,*iQ);
}

template <class R, class O> inline
bool operator > (const MPoly<R,O> & P, const MPoly<R,O> & Q){
    return Q < P;
}

template <class R, class O> inline
bool operator >= (const MPoly<R,O> & P, const MPoly<R,O> & Q){
    return !(P < Q);
}

template <class R, class O> inline
bool operator <= (const MPoly<R,O> & P, const MPoly<R,O> & Q){
    return !(Q < P);
}

//----------------------------------------------------------------------
// Comput operators and functions
//----------------------------------------------------------------------
template <class R, class O> inline
MPoly<R,O> MPoly<R,O>::operator - () const {
    MPoly<R,O> r(*this);
    for (iterator i = r.begin(); i != r.end(); i++) *i = -(*i);
    return (r);
}
template <class R, class O> inline
MPoly<R,O> operator / (const MPoly<R,O> & P, const MPoly<R,O> & Q)
{
    return (MPoly<R,O>(P)/= Q);
}

template <class R, class O> inline
MPoly<R,O>  operator + (const MPoly<R,O> & P, const MPoly<R,O> & Q)
{
    return (MPoly<R,O>(P)+= Q);
}

template <class R, class O> inline
MPoly<R,O>  operator - (const MPoly<R,O> & P, const MPoly<R,O> & Q)
{
    return (MPoly<R,O>(P)-= Q);
}

template <class R, class O> inline
MPoly<R,O> operator * (const MPoly<R,O> & P,
                       const typename MPoly<R,O>::monom_t & m)
{
    return (MPoly<R,O>(P) *= m);
}

template <class R, class O> inline
MPoly<R,O> operator * (const MPoly<R,O> & P, const MPoly<R,O> & Q)
{
    return (MPoly<R,O>(P)*= Q);
}
//----------------------------------------------------------------------
template <class R, class O> inline
int leading_var (const MPoly<R,O> & P) {
    return MPOLYNOMIAL__lvar(P.rep);
}

template <class R, class O> inline
int Degree (const MPoly<R,O> & P) {
    return MPOLYNOMIAL__degree(P.rep);
}

//----------------------------------------------------------------------
// Input/putput operators
//----------------------------------------------------------------------
template <class R, class O> inline
std::istream & operator >> (std::istream & is, MPoly<R,O> & P) {
    std::string s,ss = "";
    int p;
    char *buf;

    for (;;) {
        is >> s;
        p = s.find(';');
        if (p == -1)
            ss += s;
        else {
            ss += s.substr(0,p+1);
            break;
        }
    }
    buf = (char *)malloc(sizeof(char)*ss.length() +1);
    strcpy(buf,ss.c_str());
    P = MPoly<R,O>(buf);
    free(buf);
    return is;
}

template <class R, class O> inline
std::ostream & operator << (std::ostream & os, const MPoly<R,O> & P){
    typename MPoly<R,O>::const_iterator i = P.begin();
    if ( i == P.end() )
        os << '0';
    while ( i != P.end() )
    {
        os << (*i++);
        if ( i != P.end() && has_positive_sign(i->GetCoeff()) )
            os << '+';
    }
    return os;
}

template <class R, class O, class T>
MPoly<R,O> pow(const MPoly<R,O> & P, T n)
{
    MPoly<R,O> result(typename MPoly<R,O>::coeff_t(1));
    MPoly<R,O> temp(P);
    assert(n>=0);
    while(n > 0 )
    {
        if (n & 1) result *= temp;
        temp *= temp;
        n /= 2;
    }
    return result;
}


// Substitute the i-th variable of the polynomial P by the polynomial Q
template <class R, class O>
MPoly<R,O> subs(MPoly<R,O> P, int i, MPoly<R,O> Q)
{
    typename MPoly<R,O>::iterator it;
    MPoly<R,O> result;
    for(it=P.begin(); it!=P.end(); it++) {
        MPoly<R,O> temp_pol;
        typename MPoly<R,O>::monom_t temp(*it);
        typename MPoly<R,O>::monom_t::exponent_t t;
        t = temp[i];
        if(t > 0) {
            temp.SetDegree(i,0);
            temp_pol= pow(Q,t) * temp;
        } else {
            temp_pol = MPoly<R,O>(temp);
        }
        result += temp_pol;
    }

    return result;
}

#endif //