sparse_dual.hpp

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of the realroot library.
*      Author(s): B. Mourrain, GALAAD, INRIA                 
**********************************************************************/

#ifndef realroot_sparse_dual_hpp
#define realroot_sparse_dual_hpp

#include <realroot/sparse_monomials.hpp>

//--------------------------------------------------------------------
namespace mmx {
  namespace sparse {
//----------------------------------------------------------------------
template <class C,
          class O=DegRevLex>
struct dual: public monomial_seq<C,O> {

  // Types
    
  typedef monomial_seq<C,O>                  R;
  typedef monomial_seq<C,O>                  base_t;

  typedef typename R::monom_t                monom_t;
  typedef C                                  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 monomial_seq<C,O>                  poly_t;
  typedef dual<C,O>                          self_t;

  // Constructors
  dual(): poly_t() {}
  dual(const coeff_t & c): poly_t(c) {}
  dual(const coeff_t & c, int d, int v): poly_t(c,-d,v) {}
  dual(const monom_t & m):  poly_t(m) {neg_exponent(*this);}
  //  dual(int n, const monom_t * t): poly_t(n,t) {neg_exponent(*this);}

  dual(const self_t & P): poly_t(P){}
  dual(const char * s, variables& vars):
    poly_t(s,vars) {neg_exponent(*this);}
  dual(const char * s, const variables& vars):
    poly_t(s,vars) {neg_exponent(*this);}

  template<class Q>
  dual(const dual<C,Q>& p) : poly_t(p) {}
  
  // Destructor
  // no dynamic allocation at this level.

  //  static Variable var;

  // Assignment operators
  self_t & operator= (const self_t & x)
  { 
    this->base_t::operator=(x); return *this;
  }
  self_t & operator= (const coeff_t & c) {*this=self_t(c); return *this;}
  self_t & operator= (int n)             {*this=self_t(n); return *this;}
  self_t & operator= (const monom_t & m) {*this=self_t(m); return *this;}

  C operator()(const poly_t& p);

};


template<class C, class O>
void neg_exponent(dual<C,O>& p)
{
  typedef typename dual<C,O>::iterator iterator;
  for(iterator it=p.begin();it != p.end();it++)
    {
      for(int j=0;j<(it->nvars()+1);j++)
    {
      it->set_expt(j,-(*it)[j]);
    }
    }
}
//--------------------------------------------------------------------
template<class C, class O>
C dual<C,O>::operator()(const monomial_seq<C,O>& p)
{
  monomial_seq<C,O> t= ((monomial_seq<C,O>)(*this))*p;
  typedef typename monomial_seq<C,O>::reverse_iterator iterator;
  iterator it;
  for(it=t.rbegin(); it != t.rend();it++)
    {
      bool nul=true;
      for(unsigned j=0;j<(it->nvars()+1) && nul;j++)  
    {
      if((*it)[j]!=0) nul=false;
    }
      if(nul) return it->coeff();
    }
  return C(0);
}
//--------------------------------------------------------------------
template<class C, class O> inline void
add(dual<C,O>& res, const dual<C,O>& a, const dual<C,O>& b) {
  add((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,(const monomial_seq<C,O>&)b);
} 
//--------------------------------------------------------------------
template<class C, class O> inline void
add(dual<C,O>& res, const dual<C,O>& a, const C& b) {
  add((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,b);
}
//--------------------------------------------------------------------
template<class C, class O> inline void
add(dual<C,O>& res, const C& b , const dual<C,O>& a) {
  add((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,b);
}
//--------------------------------------------------------------------
template<class C, class O> inline void
sub(dual<C,O>& res, const dual<C,O>& a, const dual<C,O>& b) {
  sub((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,(const monomial_seq<C,O>&)b);
} 
//--------------------------------------------------------------------
template<class C, class O> inline void
sub(dual<C,O>& res, const dual<C,O>& a, const C& b) {
  sub((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,b);
}
//--------------------------------------------------------------------
template<class C, class O> inline void
sub(dual<C,O>& res, const C& a, const dual<C,O>& b) {
  sub((monomial_seq<C,O>&)res,a,(const monomial_seq<C,O>&)b);
} 
//--------------------------------------------------------------------
template<class C, class O> inline void
mul(dual<C,O>& res, const dual<C,O>& a, const dual<C,O>& b) {
  mul((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,(const monomial_seq<C,O>&)b);
}
//--------------------------------------------------------------------
template<class C, class O> inline void
mul(dual<C,O>& res, const dual<C,O>& a, const C& b) {
  mul((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,b);
} 
//--------------------------------------------------------------------
template<class C, class O> inline void
mul(dual<C,O>& res, const C& b, const dual<C,O>& a) {
  mul((monomial_seq<C,O>&)res,(const monomial_seq<C,O>&)a,b);
} 
//--------------------------------------------------------------------
template<class C, class O> void
mul(dual<C,O>& res, const monomial_seq<C,O>& p, const dual<C,O>& l) {
  monomial_seq<C,O> t;
  mul(t,p,(monomial_seq<C,O>)l);
  typedef typename dual<C,O>::const_iterator iterator;
  for(iterator it=t.begin();it != t.end();it++)
    {
      bool neg=true;
      for(unsigned j=0;j<it->size() && neg;j++) {
    if((*it)[j]>0) neg=false;
      }
      if(neg) res.push_back(*it);
    }
}

// template<class C, class O> Seq<polynom< monomial_seq<C,O> >
// coefficients(const Polynomial & f, const int & v ) {
//     unsigned N = degree(f,v)+1;
//     Seq<Polynomial> r(N);
//     typename Polynomial::Monomial m;
//     for(typename Polynomial::const_iterator it=f.begin(); it != f.end(); it++) {
//       m=*it;
//       if (v<(int)m.size()) m.set_expt(v,0);
//       r[(*it)[v]]+= Polynomial(m);
//     }
//     return r;
//   }


//--------------------------------------------------------------------
template<class C, class O> void
mul(monomial_seq<C,O>& res, const dual<C,O>& l, const monomial_seq<C,O>& p) {
  monomial_seq<C,O> t;
  mul(t,(monomial_seq<C,O>)l,p);
  typedef typename dual<C,O>::const_iterator iterator;
  for(iterator it=t.begin();it != t.end();it++)
    {
      bool pos=true;
      for(unsigned j=0;j<it->size() && pos;j++) {
    if((*it)[j]<0) pos=false;
      }
      if(pos) res.push_back(*it);
    }
}

//----------------------------------------------------------------------    
template<class C, class O> inline void 
div(dual<C,O> &f, const C& c) {
  for(typename dual<C,O>::iterator it=f.begin(); it != f.end(); it++) 
    it->coeff()/=c;
}

//----------------------------------------------------------------------
template<class C, class O> void 
div(dual<C,O> &r, const dual<C,O>& p, const C & c) {
  r = p; div(r,c);
}

//--------------------------------------------------------------------
#ifndef MMX_WITH_PLUS_SIGN
#define MMX_WITH_PLUS_SIGN
  template<class X> bool with_plus_sign(const X& x) {return x>0;}
#endif  //MMX_WITH_PLUS_SIGN

template <class OSTREAM, class C> inline OSTREAM & 
print_dual (OSTREAM & os, const monom<C>& m, const variables& V) { 
  if (m.coeff() ==0 ) {
    os<<"0"; return os;
  }
  int i=m.size()-1;
  while(i> -1 && m[i]==0) i--;
  if(i<0)
    os <<m.coeff();
  else {
    if(m.coeff() != 1) {
      if(m.coeff()==-1) 
     os << "-"; 
      else {  
    os << m.coeff(); os <<"*"; 
      }
    }
  }
  bool first=true;
  for (unsigned i = 0; i < m.size(); ++i) {
    if (m[i] !=0) {
      if(first)
    first =false;
      else
    os<<'*';
      os << 'd'<<V[i];
      if (m[i] != -1)
    os << '^' <<(-m[i]);
    }
  }
  return os;
}
//--------------------------------------------------------------------

template <class OSTREAM, class C, class O> OSTREAM & 
print(OSTREAM & os, const dual<C,O> & P, const variables& V ) {
 if ( P.begin() == P.end() ) {
   os << '0';
   return os;
 }
 typename dual<C,O>::const_iterator i = P.begin();
 while ( i != P.end() ) {
   if ( with_plus_sign(i->coeff()) && i !=P.begin())
     os <<'+';
   print_dual(os,*i++,V);
 }
 return os;
}
//--------------------------------------------------------------------__
} //namespace sparse

//======================================================================
# define Polynomial sparse::dual<C,O>
# define TMPL template<class C, class O>

template<class A, class B> struct use;
struct operators_of;

TMPL struct use<operators_of,Polynomial> {

    static inline void 
    add (Polynomial &r, const Polynomial &a ) { 
      sparse::add (r, a); 
    }
    static inline void 
    add (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      sparse::add (r, a, b); 
    }
    static inline void
    add (Polynomial &r, const Polynomial &a, const C & b) {
      sparse::add (r, a, b); 
    }
    static inline void
    add (Polynomial &r, const C & a, const Polynomial &b) {
      sparse::add (r, b, a); 
    }
    static inline void
    sub (Polynomial &r, const Polynomial &a ) {
      sparse::sub (r, a ); 
    }
    static inline void
    sub (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      sparse::sub (r, a, b); 
    }
    static inline void
    sub (Polynomial &r, const C & a, const Polynomial &b) {
      sparse::sub (r, a, b); 
    }
    static inline void
    sub (Polynomial &r, const Polynomial &a, const C & b) {
      sparse::sub (r, a, b); 
    }
    static inline void
    mul (Polynomial &r, const Polynomial &a ) {
      sparse::mul (r, a ); 
    }
    static inline void
    mul (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      sparse::mul (r, a, b); 
}
    static inline void
    mul (Polynomial &r, const Polynomial &a, const C & b) {
      sparse::mul (r, a, b); }
    static inline void
    mul (Polynomial &r, const C & a, const Polynomial &b) {
      sparse::mul (r, b, a); 
    }
    static inline void
    div (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      sparse::div (r, a, b); 
    }
    static inline void
    div (Polynomial &r, const Polynomial &a, const C & b) {
      sparse::div (r, a, b); 
    }
    static inline void
    rem (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      sparse::rem (r, b, a); 
    }
};
# undef Polynomial
# undef TMPL
//====================================================================
} //namespace mmx
#endif // realroot_sparse_dual_hpp