Polynomial.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */



#ifndef SLIP_POLYNOMIAL_HPP
#define SLIP_POLYNOMIAL_HPP
#include <iostream>
#include <iterator>
#include <cassert>
#include <numeric>
#include <algorithm>
#include <complex>
#include <cstddef>
#include <string>
#include "Array.hpp"
#include "Vector.hpp"
#include "Matrix.hpp"
#include "convolution.hpp"
#include "statistics.hpp"
#include "linear_least_squares.hpp"
#include "linear_algebra_eigen.hpp"
#include "polynomial_algo.hpp"

#include <boost/serialization/access.hpp>
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/string.hpp>
#include <boost/serialization/complex.hpp>
#include <boost/serialization/version.hpp>

namespace slip
{

template <typename T>
class Polynomial;

template <typename T>
std::ostream& operator<<(std::ostream & out, 
                                                 const Polynomial<T>& v);

template<typename T>
bool operator==(const Polynomial<T>& x, 
                                const Polynomial<T>& y);

template<typename T>
bool operator!=(const Polynomial<T>& x, 
                                const Polynomial<T>& y);

template<typename T>
bool operator<(const Polynomial<T>& x, 
                       const Polynomial<T>& y);

template<typename T>
bool operator>(const Polynomial<T>& x, 
                       const Polynomial<T>& y);

template<typename T>
bool operator<=(const Polynomial<T>& x, 
                                const Polynomial<T>& y);

template<typename T>
bool operator>=(const Polynomial<T>& x, 
                                const Polynomial<T>& y);

template <typename T>
class Polynomial
{
public :

  typedef T value_type;
  typedef Polynomial<T> self;
  typedef const Polynomial<T> const_self;

  typedef value_type* pointer;
  typedef const value_type* const_pointer;
  typedef value_type& reference;
  typedef const value_type& const_reference;

  typedef ptrdiff_t difference_type;
  typedef std::size_t size_type;

  typedef pointer iterator;
  typedef const_pointer const_iterator;

  typedef std::reverse_iterator<iterator> reverse_iterator;
  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;


   //default iterator_category of the container
  typedef std::random_access_iterator_tag iterator_category;

  
  Polynomial();
 
  Polynomial(const size_type n);
  
  Polynomial(const size_type n,
                     const T& val);
 
  Polynomial(const size_type n,
                     const T* val);
  
  template<typename InputIterator>
  Polynomial(InputIterator first,
                     InputIterator last)
  {
    assert(last > first);
    array_->resize(last - first);
    std::copy(first,last, this->begin());
  }

  Polynomial(const self& rhs);

  ~Polynomial();
  

 
  iterator begin();

  const_iterator begin() const;

  iterator end();
 
  const_iterator end() const;
  
  reverse_iterator rbegin();
  
  const_reverse_iterator rbegin() const;

  reverse_iterator rend();

  const_reverse_iterator rend() const;
 
  friend std::ostream& operator<< <>(std::ostream & out, 
                                                                     const self& v);
  self& operator=(const self & rhs);

  self& operator=(const T& val);

  void fill(const T& value)
  {
    if(value != T(0))
      {
                std::fill_n(this->begin(),this->size(),value);
      }
    else
      {
                array_->resize(1);
                (*array_)[0] = value;
      }      
  } 
  
  void fill(const T* value)
  {
    std::copy(value,value + this->size(), this->begin());
    this->kill_zeros();
  } 

  template<typename InputIterator>
  void fill(InputIterator first,
                    InputIterator last)
  {
    std::copy(first,last, this->begin());
    this->kill_zeros();
  }
  friend bool operator== <>(const Polynomial<T>& x, 
                                                    const Polynomial<T>& y);

  friend bool operator!= <>(const Polynomial<T>& x, 
                                                    const Polynomial<T>& y);

  friend bool operator< <>(const Polynomial<T>& x, 
                                                   const Polynomial<T>& y);

  friend bool operator> <>(const Polynomial<T>& x, 
                                                   const Polynomial<T>& y);

  friend bool operator<= <>(const Polynomial<T>& x, 
                                                    const Polynomial<T>& y);

  friend bool operator>= <>(const Polynomial<T>& x, 
                                                    const Polynomial<T>& y);


  reference operator[](const size_type i);

  const_reference operator[](const size_type i) const;

  reference operator()(const size_type i);

  const_reference operator()(const size_type i) const;


  self& operator+=(const T& val);
  self& operator-=(const T& val);
  self& operator*=(const T& val);
  self& operator/=(const T& val);

  self  operator-() const;


  self& operator+=(const self& rhs);
  self& operator-=(const self& rhs);
  self& operator*=(const self& rhs);
  //  self& operator/=(const self& rhs);
    


  self& derivative();


  self& integral();


  double fit(iterator xi_first,
                     iterator xi_last,
                     iterator yi_first,
                     const std::size_t order);
 

  slip::Matrix<T> companion() const;
  
  slip::Vector<std::complex<T> > roots() const;
  
  T evaluate(const T& x) const;

  template<typename InputIterator, typename OutputIterator>
  void multi_evaluate(InputIterator first, InputIterator last,
                                      OutputIterator result) const
  {
    for(;first != last; first++, result++)
      {
                *result = this->evaluate(*first);
      }

  }


  
  

  std::string name() const;


  size_type size() const;

  size_type degree() const;

  size_type order() const;

  size_type max_size() const;

  bool empty()const;

  void swap(self& M);



  
private:
  slip::Array<T>* array_;


  void kill_zeros(); 
  
  std::size_t count_zeros(const std::size_t size);

  friend class boost::serialization::access;
    template<class Archive>
    void save(Archive & ar, const unsigned int version) const
    {
      ar & this->array_ ;
    }
    template<class Archive>
    void load(Archive & ar, const unsigned int version)
    {
      ar & this->array_ ;
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER();

};


template<typename T>
Polynomial<T> operator+(const Polynomial<T>& V1, 
                                                const Polynomial<T>& V2);

template<typename T>
Polynomial<T> operator+(const Polynomial<T>& V, 
                                    const T& val);

template<typename T>
Polynomial<T> operator+(const T& val, 
                                    const Polynomial<T>& V);


template<typename T>
Polynomial<T> operator-(const Polynomial<T>& V1, 
                                    const Polynomial<T>& V2);

template<typename T>
Polynomial<T> operator-(const Polynomial<T>& V, 
                                    const T& val);

template<typename T>
Polynomial<T> operator-(const T& val, 
                                    const Polynomial<T>& V);

template<typename T>
Polynomial<T> operator*(const Polynomial<T>& V1, 
                                    const Polynomial<T>& V2);

template<typename T>
Polynomial<T> operator*(const Polynomial<T>& V, 
                                    const T& val);

template<typename T>
Polynomial<T> operator*(const T& val, 
                                    const Polynomial<T>& V);



  template<typename T>
  Polynomial<T> operator/(const Polynomial<T>& V, 
                                      const T& val);

}//slip::

namespace slip
{
  template<typename T>
  inline
  Polynomial<T>::Polynomial():
    array_(new slip::Array<T>())
  {}

  template<typename T>
  inline
  Polynomial<T>::Polynomial(const typename Polynomial<T>::size_type n):
    array_(new slip::Array<T>(n+1))
  {}
  
  template<typename T>
  inline
  Polynomial<T>::Polynomial(const typename Polynomial<T>::size_type n,
                                                    const T& val):
    array_(new slip::Array<T>(n+1,val))
  {}

  template<typename T>
  inline
  Polynomial<T>::Polynomial(const typename Polynomial<T>::size_type n,
                                                    const T* val):
    array_(new slip::Array<T>(n+1,val))   
  {
    this->kill_zeros();
  }
  
  template<typename T>
  inline
  Polynomial<T>::Polynomial(const Polynomial<T>& rhs):
    array_(new slip::Array<T>((*rhs.array_)))
  {}
  
  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator=(const Polynomial<T> & rhs)
  {
    if(this != &rhs)
      {
                *array_ = *(rhs.array_);
      }
    
    return *this;
  }

  template<typename T>
  inline
  Polynomial<T>::~Polynomial()
  {
    delete array_;
  }

  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator=(const T& val)
  {
    array_->operator=(val);
    return *this;
  }
  

  template<typename T>
  inline
  typename Polynomial<T>::iterator Polynomial<T>::begin()
  {
    return array_->begin();
  }

  template<typename T>
  inline
  typename Polynomial<T>::iterator Polynomial<T>::end()
  {
    return array_->end();
  }
  
  template<typename T>
  inline
  typename Polynomial<T>::const_iterator Polynomial<T>::begin() const
  {
    return array_->begin();
  }

  template<typename T>
  inline
   typename Polynomial<T>::const_iterator Polynomial<T>::end() const
  {
    return array_->end();
  }
  

  template<typename T>
  inline
  typename Polynomial<T>::reverse_iterator Polynomial<T>::rbegin()
  {
    return array_->rbegin();
  }

  template<typename T>
  inline
   typename Polynomial<T>::reverse_iterator Polynomial<T>::rend()
  {
    return array_->rend();
  }

  template<typename T>
  inline
   typename Polynomial<T>::const_reverse_iterator Polynomial<T>::rbegin() const
  {
    return array_->rbegin();
  }

  template<typename T>
  inline
  typename Polynomial<T>::const_reverse_iterator Polynomial<T>::rend() const
  {
    return array_->rend();
  }


  template <typename T>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const Polynomial<T>& v)
  {
    if(v.size() != 0)
      {
                out<<v[0];
                for(typename Polynomial<T>::size_type i = 1; i < v.size(); ++i)
                  {
                    out<<" + "<<v[i]<<" x^"<<i;
                  }
                out<<std::endl;
      }
    else
      {
                out<<""<<std::endl;
      }
                return out;
  }

  template<typename T>
  inline
  typename Polynomial<T>::reference Polynomial<T>::operator[](const typename Polynomial<T>::size_type i)
  {
    return (*array_)[i];
  }
  
  template<typename T>
  inline
  typename Polynomial<T>::const_reference Polynomial<T>::operator[](const typename Polynomial<T>::size_type i) const 
  {
    return (*array_)[i];
  }
  
  template<typename T>
  inline
  typename Polynomial<T>::reference Polynomial<T>::operator()(const typename Polynomial<T>::size_type i)
  {
    return (*array_)(i);
  }
  
  template<typename T>
  inline
  typename Polynomial<T>::const_reference Polynomial<T>::operator()(const typename Polynomial<T>::size_type i) const
  {
    return (*array_)(i);
  }

  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::derivative() 
  {
      
    
   if(this->degree()!=0)
      {
                Polynomial<T> Ptmp(*this);
                array_->resize(this->degree());
    for(typename Polynomial<T>::size_type i = 1; i <= Ptmp.degree(); ++i)
      {
                (*this)[i-1] = T(i)*Ptmp[i]; 
      }
    this->kill_zeros();
      }
     else
       {
                 Polynomial<T> Ptmp(0,static_cast<T>(0));
                 (*this) = Ptmp;
       }
    return *this;
  }

  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::integral() 
  {
    Polynomial<T> Ptmp(*this);
    array_->resize(this->size()+1);
    (*this)[0] = T();
    for(typename Polynomial<T>::size_type i = 1; i < array_->size(); ++i)
      {
                (*this)[i] = Ptmp[i-1] / i; 
      }
    this->kill_zeros();
    return *this;
  }

  template<typename T>
  inline
  double Polynomial<T>::fit(typename Polynomial<T>::iterator xi_first,
                                                    typename Polynomial<T>::iterator xi_last,
                                                    typename Polynomial<T>::iterator yi_first,
                                                    const std::size_t order)
  {
     typedef typename Polynomial<T>::iterator iterator;         
                
    //input data
    std::size_t coeff_size = order + 1;
       
    //output data
    double chi_square = 0.0;
    array_->resize(coeff_size);
    slip::EvalPowerBasis<T,iterator> power_basis;
    chi_square = slip::svd_least_square<iterator,iterator,iterator>(xi_first,xi_last,yi_first,this->begin(),this->end(),power_basis);
    return chi_square;
  }


  template<typename T>
  inline
  slip::Matrix<T> Polynomial<T>::companion() const
   {
     
     std::size_t degree = this->degree();
     assert(degree >= 2);
     slip::Matrix<T> tmp(degree,degree,T(0));
       
     T norm = (*array_)[degree];
     for(std::size_t k = 0, j = (tmp.cols() - 1); k < degree; j--,++k)
       {
                 tmp[0][j] = -(*array_)[k] / norm;
       }
     
     std::size_t i = 1;
     std::size_t j = 0;
     
     for(std::size_t k = 1; k < degree; ++k)
       {
                 tmp[i][j] = 1;
                 
                 i += 1;
                 j += 1;
       }
     return tmp;
   }

  template<typename T>
  inline
  slip::Vector<std::complex<T> > Polynomial<T>::roots() const
   {
     slip::Matrix<T> companion = this->companion();
     slip::Vector<std::complex<T> > tmp(companion.rows());
     slip::eigen(companion,tmp,false);

     return tmp;
   }

  template<typename T>
  inline
  T Polynomial<T>::evaluate(const T& x) const
  {
    return slip::eval_horner(this->begin(),this->end(),x);
  }

 
  template<typename T>
  inline
  std::string 
  Polynomial<T>::name() const {return "Polynomial";} 
  
  template<typename T>
  inline
  typename Polynomial<T>::size_type Polynomial<T>::size() const {return array_->size();}


  template<typename T>
  inline
  typename Polynomial<T>::size_type Polynomial<T>::degree() const {
    if(this->size() != 0)
      return (this->size() - 1);
    else
      return 0;
  }

  template<typename T>
  inline
  typename Polynomial<T>::size_type Polynomial<T>::order() const 
  {
    return this->degree();
  }

  template<typename T>
  inline
  typename Polynomial<T>::size_type Polynomial<T>::max_size() const 
  {
    return typename Polynomial<T>::size_type(-1)/sizeof(T);
  }
  
  template<typename T>
  inline
  bool Polynomial<T>::empty()const {return array_->empty();}

  template<typename T>
  inline
  void Polynomial<T>::swap(Polynomial<T>& V)
  {
    array_->swap(*(V.array_));
  }


  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator+=(const T& val)
  {
    //    std::transform(array_->begin(),array_->end(),array_->begin(),std::bind2nd(std::plus<T>(),val));
    //this->kill_zeros();
    *array_->begin()+=val;
    return *this;
  }

  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator-=(const T& val)
  {
   //  std::transform(array_->begin(),array_->end(),array_->begin(),std::bind2nd(std::minus<T>(),val));
//     this->kill_zeros();
    *array_->begin()-=val;
    return *this;
  }

  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator*=(const T& val)
  {
    std::transform(array_->begin(),array_->end(),array_->begin(),std::bind2nd(std::multiplies<T>(),val));
    this->kill_zeros();
    return *this;
  }
  
  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator/=(const T& val)
  {
    std::transform(array_->begin(),array_->end(),array_->begin(),std::bind2nd(std::multiplies<T>(),T(1)/val));
    this->kill_zeros();
    return *this;
  }
  

  template<typename T>
  inline
  Polynomial<T> Polynomial<T>::operator-() const
  {
    Polynomial<T> tmp(*this);
    std::transform(array_->begin(),array_->end(),tmp.begin(),std::negate<T>());
    return tmp;
  }

  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator+=(const Polynomial<T>& rhs)
  {
    std::size_t max_size = std::max(this->size(),rhs.size());
    Polynomial<T> tmp(*this);
    array_->resize(max_size);
    
    if(max_size == tmp.size())
      {
                std::transform((*(rhs.array_)).begin(),(*(rhs.array_)).end(),tmp.begin(),this->begin(),std::plus<T>());
                std::copy(tmp.begin() + rhs.size(),tmp.end(),this->begin()+rhs.size());
      }
    else
      {
                std::transform(tmp.begin(),tmp.end(),(*(rhs.array_)).begin(),this->begin(),std::plus<T>());
                std::copy(rhs.begin() + tmp.size(),rhs.end(),this->begin()+tmp.size());
    
      }
    this->kill_zeros();
    return *this;
  }
  
  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator-=(const Polynomial<T>& rhs)
  {
    std::size_t max_size = std::max(this->size(),rhs.size());
    Polynomial<T> tmp(*this);
    array_->resize(max_size);
    if(max_size == tmp.size())
      {
                std::transform(tmp.begin(),tmp.begin()+rhs.size(),rhs.begin(),this->begin(),std::minus<T>());
                std::copy(tmp.begin() + rhs.size(),tmp.end(),this->begin()+rhs.size());
      }
    else
      {
                std::transform(tmp.begin(),tmp.end(),(*(rhs.array_)).begin(),this->begin(),std::minus<T>());
                std::copy(rhs.begin() + tmp.size(),rhs.end(),this->begin()+tmp.size());
                std::transform(this->begin()+tmp.size(),this->end(),this->begin()+tmp.size(),std::negate<T>());
    
      }
    this->kill_zeros();
    return *this;
  }
  

  template<typename T>
  inline
  Polynomial<T>& Polynomial<T>::operator*=(const Polynomial<T>& rhs)
  {
    assert(rhs.size() > 0);
    assert(this->size() > 0);
    Polynomial Ptmp(*this);
    
    array_->resize((rhs.degree() + Ptmp.degree()) + 1);
    if(Ptmp.size() > rhs.size())
      {
                slip::full_convolution(Ptmp.begin(),Ptmp.end(),rhs.begin(),rhs.end(),this->begin());
      }
    else
      {
                slip::full_convolution(rhs.begin(),rhs.end(),Ptmp.begin(),Ptmp.end(),this->begin());
      }
    this->kill_zeros();
    return *this;
  }


  template<typename T>
  inline
  std::size_t Polynomial<T>::count_zeros(std::size_t size)
  {
    if(size != 0)
      {
                std::size_t new_size  = size;
                if((*this)[new_size-1] == T(0))
                  {
                    return count_zeros(new_size-1);
                  }
                else
                  {
                    return new_size;
                  }
      }
    else
      {
                return 1;
      }
   }  
  
  template<typename T>
  inline
  void Polynomial<T>::kill_zeros()
  {
    std::size_t new_size = this->count_zeros(this->size());
    if(new_size != this->size())
      {
                slip::Polynomial<T> tmp(*this);
                (this->array_)->resize(this->count_zeros(new_size));
                std::copy(tmp.begin(),tmp.begin()+new_size,this->begin());
      }
                                                  
  }
  
  /* @{ */
  template<typename T>
  inline
  Polynomial<T> operator+(const Polynomial<T>& V1, 
                                                  const Polynomial<T>& V2)
  {
    Polynomial<T> tmp(V1);
    tmp+=V2;
    return tmp;
  }

  template<typename T>
  inline
  Polynomial<T> operator+(const Polynomial<T>& V1, 
                                                  const T& val)
  {
    Polynomial<T> tmp(V1);
    tmp+=val;
    return tmp;
  }

  template<typename T>
  inline
  Polynomial<T> operator+(const T& val,
                                                  const Polynomial<T>& V1)
  {
    return V1 + val;
  }
  
  template<typename T>
  inline
  Polynomial<T> operator-(const Polynomial<T>& V1, 
                                                  const Polynomial<T>& V2)
  {
    Polynomial<T> tmp(V1);
    tmp-=V2;
    return tmp;
  }
  
  template<typename T>
  inline
  Polynomial<T> operator-(const Polynomial<T>& V1, 
                                                  const T& val)
  {
    Polynomial<T> tmp(V1);
    tmp-=val;
    return tmp;
  }

  template<typename T>
  inline
  Polynomial<T> operator-(const T& val,
                                                  const Polynomial<T>& V1)
  {
    return -(V1 - val);
  }

  template<typename T>
  inline
  Polynomial<T> operator*(const Polynomial<T>& V1, 
                                                  const Polynomial<T>& V2)
  {
    Polynomial<T> tmp(V1);
    tmp*=V2;
    return tmp;
   //  Polynomial<T> Ptmp(V1.degree() + V2.degree());

//     if(V1.size() > V2.size())
//       {
//              slip::full_convolution(V1.begin(),V1.end(),V2.end(),0,V2.size()-1,Ptmp.begin());
//       }
//     else
//       {
//              slip::full_convolution(V2.begin(),V2.end(),V1.end(),0,V1.size()-1,Ptmp.begin());
//       }
    
//     return Ptmp;
  }
  
  template<typename T>
  inline
  Polynomial<T> operator*(const Polynomial<T>& V1, 
                                                  const T& val)
{
    Polynomial<T> tmp(V1);
    tmp*=val;
    return tmp;
}

template<typename T>
inline
Polynomial<T> operator*(const T& val,
                                                const Polynomial<T>& V1)
{
  return V1 * val;
}



template<typename T>
inline
Polynomial<T> operator/(const Polynomial<T>& V1, 
                                    const T& val)
{
    Polynomial<T> tmp(V1);
    tmp/=val;
    return tmp;
}
/* @} */
  /* @{ */
template<typename T>
inline
bool operator==(const Polynomial<T>& x, 
                                const Polynomial<T>& y)
{
  return *(x.array_) == *(y.array_); 
}

template<typename T>
inline
bool operator!=(const Polynomial<T>& x, 
                                const Polynomial<T>& y)
{
  return !(x == y);
}
/* @} */
  /* @{ */
template<typename T>
inline
bool operator<(const Polynomial<T>& x, 
                       const Polynomial<T>& y)
{
  return  *(x.array_) < *(y.array_);
}

template<typename T>
inline
bool operator>(const Polynomial<T>& x, 
                       const Polynomial<T>& y)
{
  return (y < x);
}

template<typename T>
inline
bool operator<=(const Polynomial<T>& x, 
                                const Polynomial<T>& y)
{
  return !(y < x);
}

template<typename T>
inline
bool operator>=(const Polynomial<T>& x, 
                                const Polynomial<T>& y)
{
  return !(x < y);
} 
/* @} */
}//slip::

#endif //SLIP_POLYNOMIAL_HPP