include/realroot/vector_monomials.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
*
*      Modification: et@di.uoa.gr
*                    Add lcoeff, tcoeff functions
*                    Add diff function that returns the derivative
**********************************************************************/
#ifndef realroot_univariate_hpp
#define realroot_univariate_hpp
//--------------------------------------------------------------------
#include <iostream>
#include <math.h>
#include <realroot/assign.hpp>
#include <realroot/array.hpp>
# define TMPL  template<class C>
# define TMPLX template<class C, class X>
# define Polynomial monomials<C>

//====================================================================
namespace mmx {
//--------------------------------------------------------------------
#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
//--------------------------------------------------------------------
  struct AsSize;

  //====================================================================

namespace univariate {

  template <class C>
  struct monomials {

    typedef C            value_type;
    typedef unsigned int size_type;
    typedef C*           iterator;
    typedef iterator     const_iterator;
    typedef std::reverse_iterator<const_iterator>  const_reverse_iterator;
    typedef std::reverse_iterator<iterator>        reverse_iterator;
    typedef C coeff_t;

    C * tab_;
    size_type size_;
    int degree_;

    monomials(): tab_(0), size_(0), degree_(-1) {}
    monomials(const C& c);
    monomials(const C& c, size_type d, int v=0);
    monomials(const size_type& s, C *t);
    monomials(const size_type& s, const char* t);
    monomials(C* b, C *e);
    //    {degree_=this->size_-1;checkdegree(*this);}
    monomials(const monomials<C> & p);
    //    {}//{cout <<"copy rep"<<endl;}
    ~monomials () { delete[] tab_;}

    iterator        begin()       {return iterator(tab_); }
    const_iterator  begin() const {return const_iterator(tab_); }
    iterator        end()         {return iterator(this->tab_+degree_+1); }
    const_iterator  end()   const {return const_iterator(this->tab_+degree_+1); }
    
    reverse_iterator       rend()         {return reverse_iterator(tab_); }
    const_reverse_iterator rend()   const {return reverse_iterator(tab_); }
    
    reverse_iterator       rbegin() {
      return reverse_iterator(this->tab_+degree_+1); 
    }
    const_reverse_iterator rbegin() const {
      return const_reverse_iterator(this->tab_+degree_+1); 
    }

    unsigned size()   const { return size_;}
    unsigned degree() const { return degree_;}

    monomials& operator= (const monomials<C> & v);

    inline       C & operator[] (size_type i)       {return tab_[i];}
    inline const C & operator[] (size_type i) const {return tab_[i];}
    // No test {if(i<size_) return tab_[i]; else return C(0);}

    bool operator==( const C& c ) const;

    void resize(const size_type& n);

  };

  TMPL
  Polynomial::monomials(const Polynomial & v) : size_(v.size_), degree_(v.degree_) {
    tab_ = new C[v.size_]; 
    for(size_type i=0;i<v.size_;i++) {tab_[i] = v[i];}
  }

  TMPL
  Polynomial::monomials(const C& c): size_(1), degree_(0) {
    tab_ = new C[1];
    tab_[0]=c;
  }
  TMPL
  Polynomial::monomials(const C& c, size_type l, int v): size_(l+1), degree_(l) {
    tab_ = new C[l+1];
    for (size_type i = 0; i < l; i++) tab_[i]=C(0);
    tab_[l]=c;
  }
  
  TMPL
  Polynomial::monomials(const size_type& i, C* t): size_(i), degree_(i-1) {
    tab_ = new C[i];
    for(size_type l=0;l<i;l++) tab_[l] = t[l];
    check_degree(*this);
  }

  TMPL 
  Polynomial::monomials(C* b, C* e) {
    size_=0; C* p;
    for(p=b;p!=e;p++,size_++) ;
    tab_ = new C[size_];
    degree_=0;
    for(p=b; p!=e; p++,degree_++) {tab_[degree_]=(*p);}
    check_degree(*this);
  }


  TMPL
  Polynomial::monomials(const size_type& l, const char * nm) {
    std::istringstream ins(nm);
    size_=l;
    tab_ = new C[size_];
    for(size_type i=0; i< size_; i++) ins >>tab_[i];
    degree_=l-1;
    check_degree(*this);
  }
  //----------------------------------------------------------------------
  TMPL void 
  Polynomial::resize(const size_type & i) {
    if(size_!=i){
      if (size_ != 0) {
        C* tmp = tab_;
        tab_ = new C[i];
        for(size_type j=0;j<std::min(i,size_);j++) tab_[j]=tmp[j];
        for(size_type j=std::min(i,size_);j<i;j++) tab_[j]=value_type();
        size_=i;
        delete [] tmp;      
      } else {
        size_=i; tab_= new C[i];
        for(size_type j=0;j<size_;j++) tab_[j]=value_type();
      }
    }
    degree_=i-1;
  }

  TMPL Polynomial& 
  Polynomial::operator=(const Polynomial & v) {
    if (size_ != v.size_) {
      if (size_ !=0) delete [] tab_;
      tab_ = new C[v.size_];
      size_ = v.size_;
    }
    degree_=v.degree_;
    for(size_type i=0;i<size_;i++) tab_[i] = v[i];
    return *this;
  }

  TMPL bool 
  Polynomial::operator== (const C& c) const {
    if(c==0)
      return (this->degree()<0);
    else
      return (this->degree()==0 && tab_[0]== c);
  }
  //--------------------------------------------------------------------
  template <class R>
  int degree(const R & p) {
    // std::cout <<"(Using default degree) "<<p<<std::endl;
    int i = p.size()-1;
    while(i>-1 && p[i]==(typename R::value_type)0) i--;
    return i;
  }

  TMPL inline int degree(const Polynomial& a) { return a.degree_;}


  //-------------------------------------------------------------------

  TMPL inline int 
  check_degree(Polynomial& a) {
    int l = a.degree_;
    while(l >=0 && (a.tab_[l]==0)) {l--;}
    a.degree_=l;
    return a.degree_;
  }
    
  //----------------------------------------------------------------------
  TMPL void inline
  add(Polynomial &r, const Polynomial& a, const Polynomial & b) {
    array::add(r,a,b);
    check_degree(r);
  }

  //----------------------------------------------------------------------
  TMPL void inline
  add(Polynomial &r, const Polynomial& a, const C& c) {
    if(degree(a)>-1) {
      r=a; r[0]+=c;
    }
    else
      r= Polynomial(c,0);
    check_degree(r);
  }

  //----------------------------------------------------------------------
  TMPL void inline
  add(Polynomial &r, const C& c) {
    if(degree(r)>-1) {
     r[0]+=c;
    }
    else
      r= Polynomial(c,0);
    check_degree(r);
  }
  //----------------------------------------------------------------------
  TMPL void inline
  add(Polynomial &r, const Polynomial& a) {
    array::add(r,a);
    check_degree(r);
  }
  //----------------------------------------------------------------------
  TMPL void inline
  sub(Polynomial &r, const Polynomial& a, const Polynomial & b) {
    array::sub(r,a,b);
    check_degree(r);
  }
  //----------------------------------------------------------------------
  TMPL void inline
  sub(Polynomial &r, const Polynomial& a) {
    array::sub(r,a);
    check_degree(r);
  }
  //----------------------------------------------------------------------
  TMPL void inline
  sub(Polynomial &r, const Polynomial& a, const C& c) {
    if(degree(a)>-1) {
      r=a; r[0]-=c;
    }
    else
      r= Polynomial(-c,0);
    check_degree(r);
  }
  //----------------------------------------------------------------------
  TMPLX void inline
  sub(Polynomial &r, const Polynomial& a, const X& x) {
    sub(r,a,(C)x);
  }
  //--------------------------------------------------------------------
  TMPL void inline
  sub(Polynomial &r, const C& c) {
    if(degree(r)>-1) {
     r[0]-=c;
    }
    else
      r= Polynomial(-c,0);
    check_degree(r);
  }
  //----------------------------------------------------------------------
  TMPL void inline
  mul(Polynomial &r, const Polynomial& a, const Polynomial & b) {
    assert(&r!=&a);
    
    r.resize(std::max(degree(a)+degree(b)+1,0));
    array::set_cst(r,(C)0);
    mul_index(r,a,b);
    check_degree(r);
  }

  //--------------------------------------------------------------------
  TMPL void inline 
  mul(Polynomial &r, const C & c) {
    array::mul_ext(r,c);
  }

  //--------------------------------------------------------------------
  TMPL void inline 
  mul(Polynomial &a, const Polynomial & b) {
    Polynomial r; mul(r,a,b); a=r;
  }
  //----------------------------------------------------------------------
  TMPL void inline
  mul(Polynomial &r, const Polynomial& p, const C & c) {
    r=p; mul(r,c); check_degree(r);
  }


  //----------------------------------------------------------------------
  TMPL void inline
  shift(Polynomial &r, const Polynomial& p, int d, int v=0) {
    if(d>=0) {
      r.resize(degree(p)+d+1);
      for(int i=0; i<d;i++) 
        r[i]= (typename Polynomial::coeff_t)0;
      for(int i=0; i<=degree(p);i++) 
        r[i+d]=p[i];
    }
  }

  //----------------------------------------------------------------------
  TMPL void inline
  div(Polynomial &r, const Polynomial& a, const Polynomial & b) {
    assert(&r!=&a);
    Polynomial tmp(a);
    r = Polynomial(0);
    div_rem(r,tmp,b);
    check_degree(r);
  }

  //----------------------------------------------------------------------
  TMPL void inline
  div(Polynomial &r, const Polynomial& b) {
    Polynomial tmp(r);
    r = Polynomial(0);
    div_rem(r,tmp,b);
    check_degree(r);
  }

  //--------------------------------------------------------------------
  TMPL void inline 
  div(Polynomial &r, const C & c) {
    array::div_ext(r,c);
    check_degree(r);
  }

  //-------------------------------------------------------------------
  template< class R > inline typename R::value_type
  lcoeff(const R& p) {
    return p[degree(p)];
  }
  //-------------------------------------------------------------------
  template< class R > inline typename R::value_type    
  tcoeff(const R& p)     { return p[0]; }

  //-------------------------------------------------------------------
  template <class OSTREAM, class C> OSTREAM& 
  print_as_coeff(OSTREAM & os, const C & c, bool plus)
  {
    if(plus) 
      os <<"+";
    os<<"("<<c<<")";
    return os;
  }
  
  template <class OSTREAM, class C, class VAR>  inline OSTREAM& 
  print(OSTREAM& os, const Polynomial& p, const VAR& var)
  {
    typedef typename Polynomial::value_type coeff_t;
    bool plus=false, not_one;
    if ( degree(p)<0) 
      os << coeff_t(0) ;
    else {
      for(int i= 0; i<=degree(p); i++)
      {
        if(p[i]!= (coeff_t)0) 
        {
          not_one = (p[i] != (coeff_t)1);
          if(not_one || i==0)
            print_as_coeff(os,p[i],plus);
          plus=true;
          //      if(i != degree(p) && with_plus_sign(p[i])) os<<"+"; 
          if(i>0)
            {
              //              if(p[i] == (coeff_t)(-1)) os<<"-";
              //              else 
              if(not_one) 
                os<<"*";
              os<<"x";
              if(i !=1) 
                {os<<'^'; os<<i;}
            }
//        else
//            os<<p[0];
        }
      }
      if(!plus) os << 0 ;
    }
    return os;
  }

  template <class OSTREAM, class C>  inline OSTREAM& 
  print(OSTREAM& os, const Polynomial& p) {
    return print(os,p, "x");
  }
 
  //--------------------------------------------------------------------
  template <class R, class C>
  inline void set_monomial(R & x, const C & c, unsigned n)
  {

    x.resize(n+1);
    array::set_cst(x,(typename R::value_type)0);
    x[n]=c;
    checkdegree(x);

  }
  //--------------------------------------------------------------------
  template <class R, class S>
  inline void add_cst(R & r, const S & c)  { r[0]+=c; }

  template <class R, class S>
  inline void sub_cst(R & r, const S & c)  { r[0]-=c;}
  //-------------------------------------------------------------------
  template <class R>
  inline void mul_index(R & r, const R & a, const R & b)
  {
    typedef typename R::coeff_t Scalar;
    int da=degree(a), db=degree(b),i,j;
    for (i=0;i<=db;i++) {
      if  (b[i] == 0) continue;
      for (j=0;j<=da;j++) r[i+j] += a[j]*b[i];
    }
  }
  //-------------------------------------------------------------------
  template <class R> inline
  void mul_index_it(R & r, const R & a, const R & b) 
  {
    typename R::const_iterator ia, ib;
    typename R::iterator       ir, it;
    ir=r.begin();
    for (ib=b.begin();ib !=b.end()&& ir != r.end();++ib) {
      typename R::value_type c = *ib;
      it=ir;
      if  (c==0) {++ir;continue;}
      for (ia=a.begin();ia!=a.end();++ia,++it) {
        (*it)+=(*ia)*c;
      }
      ++ir;
    }
  }
  //-------------------------------------------------------------------
  template <class R> inline
  void mul(R & a, const R & b)
  {
    if(degree(a)>=0 && degree(b)>=0)
      {
        R ta(a);
        a.resize(degree(a)+degree(b)+1);
        array::set_cst(a,0);
        mul_index(a,ta,b);
      }
    else
      array::set_cst(a,0);
  }
  //--------------------------------------------------------------------
  template <class R> inline
  void mul_index(R & a, const R & b) 
  {
    R ta(a);
    a.resize(degree(a)+degree(b)+1);
    array::set_cst(a,0);
    mul_index(a,ta,b);
  }
  //----------------------------------------------------------------------
  template <class C, class R>
  inline C eval_horner(const R & p, const C & c)
  {
    using namespace let;
    C r,cf;
    if(degree(p)>0) {
      typedef typename R::const_reverse_iterator const_reverse_iterator;
      const_reverse_iterator it=p.rbegin();
        assign(r,*it); it++;
        //      for(; it !=p.rend(); ++it) {r*=c; assign(cf,*it); r+=cf;}
        for(; it !=p.rend(); ++it) {r=r*c; assign(cf,*it); r=r+cf;}
        return r;
    } else if(degree(p)>-1)
      {
        assign(r,p[0]); return r;
      }
    else
      return C(0);
  }
  //----------------------------------------------------------------------
  template <class C, class R>
  inline C eval_horner_idx(const R & p, const C & c)
  {
    if(degree(p)>0) {
      int d=degree(p);
      C r=p[d];
      for(int i=d-1; i>-1; --i) {r*=c; r+= p[i];}
      return r;
    } else if(degree(p)>-1)
      return p[0];
    else
      return C(0);
  }
  //--------------------------------------------------------------------
  template <class C, class R>
  inline C eval(const R & p, const C & c)
  {
    return eval_horner(p,c);
  }
  //----------------------------------------------------------------------
  template <class C, class R>
  inline C eval_homogeneous(const R & p, const C & a, const C& b)
  {
    using namespace let;
    C r,y=1,cf;
    if(degree(p)>0) {
      typedef typename R::const_reverse_iterator const_reverse_iterator;
      const_reverse_iterator it=p.rbegin();
      assign(r,*it); it++;
      //        for(; it !=p.rend(); ++it) {r*=c; assign(cf,*it); r+=cf;}
      for(; it !=p.rend(); ++it) {r*=a; y *=b; assign(cf,*it); r=r+cf*y;}
      return r;
    } else if(degree(p)==0)
      {
        assign(r,p[0]); return r;
      }
    else
      return C(0);
  }

  //--------------------------------------------------------------------
  template<typename POL, typename X> inline
  int sign_at(const POL & p, const X& x)
  {
    X v= eval_horner(p,x);
    return sign(v);
  }

  template<class R>
  void div_rem(R & q, R & a, const R & b) 
  {
    typedef typename R::iterator iterator;
    int da = degree(a), db=degree(b);
    typename R::value_type lb=b[db];
    R t;
    q=(R)0;
    while (da>=db) {
      while(a[da]==0 && da>=db) da--;
      if(da>=db) {
        R m( a[da]/lb, da-db );
        mul(t,b,m);
        t[da]=a[da];
        sub(a,t);
        add(q,m);
      }
    }
  } 

//----------------------------------------------------------------------
template<class R> inline void checkdegree(R & p) { 
  int d = p.size()-1; 
  if ( d < 0 ) return;
  while ( p[d] == 0 ) d--;
  p.resize(d+1);
}
//----------------------------------------------------------------------
template <class R> inline  R diff(const R & p);
//----------------------------------------------------------------------
template<class R> void reciprocal(R & w,const R & p);
//----------------------------------------------------------------------
template<class R> void reverse(R & p, typename R::size_type n)
{
  int l=n-1;
  for(typename R::size_type i=0;i<n/2;i++) 
    std::swap(p[i],p[l-i]);
}
//----------------------------------------------------------------------
template<class R>
typename R::value_type derive(const R p,typename R::value_type x);
//----------------------------------------------------------------------
template<class R> void reduce(R & p, const typename R::size_type & e);
//----------------------------------------------------------------------
template<class R, class C> 
  void scale(R & t,const R & p,const C & l);
//----------------------------------------------------------------------
template<class R> inline
void diff(R & r, const R & p) // version + simple
{
  int n,i;
  r.resize(n=degree(p));
  for( i = 1; i < n+1; ++i) r[i-1] =p[i]*i;
}
  // ancienne version
  // template <class R> inline
  // void diff(R & r, const R & p)
  // {
  //   if(degree(p)>0)//  assert(degree(p)>0);
  //    {
  //      for(unsigned i=1;i<p.size();++i)
  //        r[i-1] =p[i]*i;
  //      r.resize(degree(p));
  //    } else r.resize(0);
  // }

//---------------------------------------------------------------------
template <class R> inline
R diff(const R & p)
{
  R r(p); 
  diff(r,p);
  return r;
}
// ancienne version
// template <class R> inline 
// R diff(const R & p)
// {
//   R r(p);
//   if(degree(p)>0) {
//      for(unsigned i = 1; i < p.size(); ++i)
//        r[i-1] = p[i] * i;
//      r.resize(degree(p));
//   } else r.resize(0);
//   return r;
// }
//---------------------------------------------------------------------
template <class R>
void reciprocal(R & w, const R & p)
{
  typedef typename R::size_type size_type;
  typedef typename R::value_type C;

  const size_type deg = univariate::degree(p)+1;
  R w0(deg),v(deg),xp(p);
  size_type K = (size_type) (std::log(p.degree()+1)/std::log(2)),j=1;
  C p0 = xp[deg-1];

  xp /= p0;
  w0 =xp/xp[0];
  v = w0;
  v *= C(-1.0);
  v[0]+=C(2);
  w=v/xp[0];
  w0=w;
  while (j <= K) {
    univariate::reduce(w0,pow(2,j+1));
    w=xp*w0;
    v =w;
    v*=C(-1.0);
    v[0]+=C(2);
    univariate::reduce(v,pow(2,j+1));
    w=w0*v;
    w0=w;
    j++;
  }
   univariate::reduce(w,deg);
   w/=p0;
}
//-------------------------------------------------------------------------
template <class T> 
inline void reduce(T & p, const typename T::size_type & e)
{
  // il me semble que tout ca c'est la meme chose que p.resize(e) 
  T temp(e);
  for (typename T::size_type i = 0; i < e; i++) temp[i]=p[i];
  p.resize(e);
  //  p.degree()=e-1;
  for (typename T::size_type i = 0; i < e; i++) p[i]=temp[i];
}
//-------------------------------------------------------------------------
template <class T> // version + simple
inline void reverse(T & p, int n)
{ for ( int i = 0; i < n/2; i++ ) std::swap(p[i],p[n-i]) ; }

// template <class T> // ancienne version
// inline void reverse(T & p,typename T::size_type n)
// {
//   for (typename T::size_type i = 0, j=n-1; i<j; i++, j--)
//     {
//       swap(p[i],p[j]);
//     }
// }

//----------------------------------------------------------------------------
template<class O, class R, class I> inline
void eval( O& p, O& dp, const R& Pol, const I& x )
{
  int n = p.size()-1;
  p  = Pol[n];
  dp = O(0);
  for ( int j = n-1; j >= 0; dp = dp*x + p, p =p*x + Pol[j], j-- ) ; 
};

template <class R> inline // wrap sur eval
typename R::value_type derive( const R& Pol, const typename R::value_type& x )
{ 
  typename R::value_type p,res;
  eval(p,res,Pol,x); 
  return res;
};

// template <class R> // ancienne version
// inline typename R::value_type
// derive(const R p, typename R::value_type x)
// {
//   typename R::size_type degree = p.degree();
//   typename R::value_type  deriv(0);
//   for (typename R::size_type j=1; j<=p.degree();j++)
//     deriv+=p[j]*j*pow(x,j-1);
//   return (typename R::value_type) deriv;
// }

//--------------------------------------------------------------------
template< class R,class C> // version + simple
void shift( R& p, const C& c )
{
  int j,k,s;
  for ( s = p.size(), j = 0; j <= s-2; j++ )
    for( k = s-2; k >= j; p[k] += c*p[k+1], k-- ) ;
};
//--------------------------------------------------------------------
template<class R, class C> 
void shift(R & r,const R & p,const C& x0)
{
  r=p; shift(r,x0);
}

// ancienne version
// {
//   R dp(p);
//   r[0]=eval_horner(dp,x0);
//   for(unsigned i=1; i< p.size();i++){
//     diff(dp,dp);
//     array::div(dp,i);
//     r[i]=eval_horner(dp,x0);
//   }
// }
//--------------------------------------------------------------------

template<class R, class C>
void scale(R & r, const R & p, const C & l)
{
  r=p;
  C s(l);
  for(unsigned i=1; i< p.size();i++){
    r[i]*=s;
    s *=l;
  }
}
//--------------------------------------------------------------------

template<class R, class C>
void inv_scale(R & r, const R & p, const C & l)
{
  int sz=p.size();
  r.resize(sz);
  sz--;
  r[sz]=p[sz];
  typename R::value_type s=l;
  for(int i=sz-1; i>=0;i--){
    r[i]=p[i]*s;
    s *=l;
  }
}

//--------------------------------------------------------------------
template<class T,class P, class C> void  
convertm2b(T& bz, const P& p, unsigned size, const C& a, const C& b) {
  typedef typename T::value_type coeff_t;
  T pps(p.size());
  array::assign(pps,p);
  
  if(a !=0)
    shift(pps,pps,a);
  
  scale(pps, pps, C(b-a));
  
  array::set_cst(bz,pps[0]);
  int dn=1, nm=1;
  
  for(unsigned j=1; j<size; j++) {
    dn *= (size-j); dn/= j;
    nm=1;
    bz[j] += pps[j]/dn;
    for(unsigned i=j+1; i<size; i++) {
      nm *= i; nm/= (i-j);
      bz[i] += pps[j]*nm/dn;
    }
  }
}
//--------------------------------------------------------------------

template<class R>
void coeff_modulo( R & r, const typename R::value_type & x ) {
  for ( typename R::iterator i = r.begin(); i != r.end(); *i %= x, ++i ) ;
};
//--------------------------------------------------------------------
template<class S, class R>
S numer(const R & f) {
  S p(0,degree(f));
  typename S::value_type d=1;
  for(int i=0;i<degree(f)+1;i++) {
    if (f[i] !=0) d=lcm(d,denominator(f[i]));
  }
  for(int i=0;i<degree(f)+1;i++) 
    p[i]=numerator(f[i])* (d/denominator(f[i]));
  return p;
}
//--------------------------------------------------------------------
} //namespace univariate
//======================================================================
# undef  Polynomial
# define Polynomial univariate::monomials<C>
# define TMPL template<class C>

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 ) { 
      univariate::add (r, a); 
    }
    static inline void 
    add (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      univariate::add (r, a, b); 
    }
    static inline void
    add (Polynomial &r, const Polynomial &a, const C & b) {
      univariate::add (r, a, b); 
    }
    static inline void
    add (Polynomial &r, const C & a, const Polynomial &b) {
      univariate::add (r, b, a); 
    }
    static inline void
    sub (Polynomial &r, const Polynomial &a ) {
      univariate::sub (r, a ); 
    }
    static inline void
    sub (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      univariate::sub (r, a, b); 
    }
    static inline void
    sub (Polynomial &r, const C & a, const Polynomial &b) {
      univariate::sub (r, a, b); 
    }
    static inline void
    sub (Polynomial &r, const Polynomial &a, const C & b) {
      univariate::sub (r, a, b); 
    }
    static inline void
    mul (Polynomial &r, const Polynomial &a ) {
      univariate::mul (r, a ); 
    }
    static inline void
    mul (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      univariate::mul (r, a, b); 
}
    static inline void
    mul (Polynomial &r, const Polynomial &a, const C & b) {
      univariate::mul (r, a, b); }
    static inline void
    mul (Polynomial &r, const C & a, const Polynomial &b) {
      univariate::mul (r, b, a); 
    }
    static inline void
    div (Polynomial &r, const Polynomial &a, const Polynomial &b) {
      univariate::div (r, a, b); 
    }
    static inline void
    div (Polynomial &r, const Polynomial &a, const C & b) {
      univariate::div (r, a, b); 
    } 
    static inline void
    div (Polynomial &a, const C & c) {
      univariate::div (a, c); 
    }
//     static inline void
//     rem (Polynomial &r, const Polynomial &a, const Polynomial &b) {
//       univariate::rem (r, b, a); 
//     }
};
# undef Polynomial
# undef TMPL
//====================================================================
} //namespace mmx 
//======================================================================
# undef TMPL
# undef TMPLX
# undef Polynomial
#endif // realroot_univariate_hpp

---