alpar@2086: /* -*- C++ -*-
alpar@2086:  *
alpar@2086:  * This file is a part of LEMON, a generic C++ optimization library
alpar@2086:  *
alpar@2553:  * Copyright (C) 2003-2008
alpar@2086:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@2086:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@2086:  *
alpar@2086:  * Permission to use, modify and distribute this software is granted
alpar@2086:  * provided that this copyright notice appears in all copies. For
alpar@2086:  * precise terms see the accompanying LICENSE file.
alpar@2086:  *
alpar@2086:  * This software is provided "AS IS" with no warranty of any kind,
alpar@2086:  * express or implied, and with no claim as to its suitability for any
alpar@2086:  * purpose.
alpar@2086:  *
alpar@2086:  */
alpar@2086: 
alpar@2086: #ifndef LEMON_BEZIER_H
alpar@2086: #define LEMON_BEZIER_H
alpar@2086: 
alpar@2086: ///\ingroup misc
alpar@2086: ///\file
alpar@2086: ///\brief A simple class implementing polynomials.
alpar@2086: ///
alpar@2086: ///\author Alpar Juttner
alpar@2086: 
alpar@2086: #include<vector>
alpar@2086: 
alpar@2086: namespace lemon {
alpar@2086: 
alpar@2086:   /// \addtogroup misc
alpar@2086:   /// @{
alpar@2086: 
alpar@2086:   ///Simple polinomial class
alpar@2086: 
alpar@2086:   ///This class implements a polynomial where the coefficients are of
alpar@2086:   ///type \c T.
alpar@2086:   ///
alpar@2086:   ///The coefficients are stored in an std::vector.
alpar@2086:   template<class T>
alpar@2086:   class Polynomial
alpar@2086:   {
alpar@2086:     std::vector<T> _coeff;
alpar@2086:   public:
alpar@2086:     ///Construct a polynomial of degree \c d.
alpar@2086:     explicit Polynomial(int d=0) : _coeff(d+1) {}
alpar@2086:     ///\e
alpar@2086:     template<class U> Polynomial(const U &u) : _coeff(1,u) {}
alpar@2086:     ///\e
alpar@2086:     template<class U> Polynomial(const Polynomial<U> &u) : _coeff(u.deg()+1)
alpar@2086:     {
deba@2386:       for(int i=0;i<int(_coeff.size());i++) _coeff[i]=u[i];
alpar@2086:     }
alpar@2086:     ///Query the degree of the polynomial.
alpar@2086:     
alpar@2086:     ///Query the degree of the polynomial.
alpar@2086:     ///\warning This number differs from real degree of the polinomial if
alpar@2086:     ///the coefficient of highest degree is 0.
alpar@2086:     int deg() const { return _coeff.size()-1; }
alpar@2086:     ///Set the degree of the polynomial.
alpar@2086: 
alpar@2086:     ///Set the degree of the polynomial. In fact it resizes the
alpar@2086:     ///coefficient vector.
alpar@2086:     void deg(int d) { _coeff.resize(d+1);}
alpar@2086: 
alpar@2086:     ///Returns (as a reference) the coefficient of degree \c d.
alpar@2086:     typename std::vector<T>::reference operator[](int d) { return _coeff[d]; }
alpar@2086:     ///Returns (as a const reference) the coefficient of degree \c d.
alpar@2086:     typename std::vector<T>::const_reference
alpar@2086:     operator[](int d) const {return _coeff[d];}
alpar@2086:     
alpar@2086:     ///Substitute the value u into the polinomial.
alpar@2086: 
alpar@2086:     ///Substitute the value u into the polinomial.
alpar@2086:     ///The calculation will be done using type \c R.
alpar@2086:     ///The following examples shows the usage of the template parameter \c R.
alpar@2086:     ///\code
alpar@2207:     ///  Polynomial<dim2::Point<double> > line(1);
alpar@2207:     ///  line[0]=dim2::Point<double>(12,25);
alpar@2207:     ///  line[1]=dim2::Point<double>(2,7);
alpar@2086:     ///  ...
alpar@2207:     ///  dim2::Point<double> d = line.subst<dim2::Point<double> >(23.2);
alpar@2086:     ///\endcode
alpar@2086:     ///
alpar@2086:     ///\code
alpar@2086:     ///  Polynomial<double> p;
alpar@2086:     ///  Polynomial<double> q;
alpar@2086:     ///  ...
alpar@2086:     ///  Polynomial<double> s = p.subst<Polynomial<double> >(q);
alpar@2086:     ///\endcode
alpar@2086:     template<class R,class U>
alpar@2086:     R subst(const U &u) const
alpar@2086:     {
alpar@2086:       typename std::vector<T>::const_reverse_iterator i=_coeff.rbegin();
alpar@2086:       R v=*i;
alpar@2086:       for(++i;i!=_coeff.rend();++i) v=v*u+*i;
alpar@2086:       return v;
alpar@2086:     }
alpar@2086:     ///Substitute the value u into the polinomial.
alpar@2086: 
alpar@2086:     ///Substitute the value u into the polinomial.
alpar@2086:     ///The calculation will be done using type \c T
alpar@2086:     ///(i.e. using the type of the coefficients.)
alpar@2086:     template<class U>
alpar@2086:     T operator()(const U &u) const 
alpar@2086:     {
alpar@2086:       return subst<T>(u);
alpar@2086:     }
alpar@2086:     
alpar@2086:     ///Derivate the polynomial (in place)
alpar@2086:     Polynomial &derivateMyself()
alpar@2086:     {
deba@2386:       for(int i=1;i<int(_coeff.size());i++) _coeff[i-1]=i*_coeff[i];
alpar@2086:       _coeff.pop_back();
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086:     
alpar@2086:     ///Return the derivate of the polynomial
alpar@2086:     Polynomial derivate() const
alpar@2086:     {
alpar@2086:       Polynomial tmp(deg()-1);
deba@2386:       for(int i=1;i<int(_coeff.size());i++) tmp[i-1]=i*_coeff[i];
alpar@2086:       return tmp;
alpar@2086:     }
alpar@2086: 
alpar@2086:     ///Integrate the polynomial (in place)
alpar@2086:     Polynomial &integrateMyself()
alpar@2086:     {
alpar@2086:       _coeff.push_back(T());
deba@2199:       for(int i=_coeff.size()-1;i>0;i--) _coeff[i]=_coeff[i-1]/i;
alpar@2086:       _coeff[0]=0;
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086:     
alpar@2086:     ///Return the integrate of the polynomial
alpar@2086:     Polynomial integrate() const
alpar@2086:     {
alpar@2086:       Polynomial tmp(deg()+1);
alpar@2086:       tmp[0]=0;
deba@2386:       for(int i=0;i<int(_coeff.size());i++) tmp[i+1]=_coeff[i]/(i+1);
alpar@2086:       return tmp;
alpar@2086:     }
alpar@2086: 
alpar@2086:     ///\e
alpar@2086:     template<class U>
alpar@2086:     Polynomial &operator+=(const Polynomial<U> &p)
alpar@2086:     {
alpar@2086:       if(p.deg()>deg()) _coeff.resize(p.deg()+1);
deba@2386:       for(int i=0;i<=int(std::min(deg(),p.deg()));i++)
alpar@2086: 	_coeff[i]+=p[i];
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086:     ///\e
alpar@2086:     template<class U>
alpar@2086:     Polynomial &operator-=(const Polynomial<U> &p)
alpar@2086:     {
alpar@2086:       if(p.deg()>deg()) _coeff.resize(p.deg()+1);
alpar@2086:       for(int i=0;i<=std::min(deg(),p.deg());i++) _coeff[i]-=p[i];
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086:     ///\e
alpar@2086:     template<class U>
alpar@2086:     Polynomial &operator+=(const U &u)
alpar@2086:     {
alpar@2086:       _coeff[0]+=u;
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086:     ///\e
alpar@2086:     template<class U>
alpar@2086:     Polynomial &operator-=(const U &u)
alpar@2086:     {
alpar@2086:       _coeff[0]+=u;
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086:     ///\e
alpar@2086:     template<class U>
alpar@2086:     Polynomial &operator*=(const U &u)
alpar@2086:     {
alpar@2086:       for(typename std::vector<T>::iterator i=_coeff.begin();i!=_coeff.end();++i)
alpar@2086: 	*i*=u;
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086:     ///\e
alpar@2086:     template<class U>
alpar@2086:     Polynomial &operator/=(const U &u)
alpar@2086:     {
alpar@2086:       for(typename std::vector<T>::iterator i=_coeff.begin();i!=_coeff.end();++i)
alpar@2086: 	*i/=u;
alpar@2086:       return *this;
alpar@2086:     }
alpar@2086: 
alpar@2086:   };
alpar@2086:   
alpar@2086:   ///Equality comparison
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///\warning Two polynomials are defined to be unequal if their degrees differ,
alpar@2086:   ///even if the non-zero coefficients are the same.
alpar@2086:   template<class U,class V>
alpar@2086:   bool operator==(const Polynomial<U> &u,const Polynomial<V> &v)
alpar@2086:   {
alpar@2086:     if(u.deg()!=v.deg()) return false;
alpar@2086:     for(int i=0;i<=u.deg();i++) if(u[i]!=v[i]) return false;
alpar@2086:     return true;
alpar@2086:   }
alpar@2086: 
alpar@2086:   ///Non-equality comparison
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///\warning Two polynomials are defined to be unequal if their degrees differ,
alpar@2086:   ///even if the non-zero coefficients are the same.
alpar@2086:   template<class U,class V>
alpar@2086:   bool operator!=(const Polynomial<U> &u,const Polynomial<V> &v)
alpar@2086:   {
alpar@2086:     return !(u==v);
alpar@2086:   }
alpar@2086: 
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator+(const Polynomial<U> &u,const Polynomial<V> &v)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp+=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086: 
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator-(const Polynomial<U> &u,const Polynomial<V> &v)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp-=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086: 
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator*(const Polynomial<U> &u,const Polynomial<V> &v)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp(u.deg()+v.deg());
alpar@2086:     for(int i=0;i<=v.deg();i++)
alpar@2086:       for(int j=0;j<=u.deg();j++)
alpar@2086: 	tmp[i+j]+=v[i]*u[j];
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator+(const Polynomial<U> &u,const V &v)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp+=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator+(const V &v,const Polynomial<U> &u)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp+=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator-(const Polynomial<U> &u,const V &v)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp-=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U>
alpar@2086:   Polynomial<U> operator-(const Polynomial<U> &u)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp(u.deg());
alpar@2086:     for(int i=0;i<=u.deg();i++) tmp[i]=-u[i];
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator-(const V &v,const Polynomial<U> &u)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=-u;
alpar@2086:     tmp+=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator*(const Polynomial<U> &u,const V &v)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp*=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator*(const V &v,const Polynomial<U> &u)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp*=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:   ///\e
alpar@2086: 
alpar@2086:   ///\relates Polynomial
alpar@2086:   ///
alpar@2086:   template<class U,class V>
alpar@2086:   Polynomial<U> operator/(const Polynomial<U> &u,const V &v)
alpar@2086:   {
alpar@2086:     Polynomial<U> tmp=u;
alpar@2086:     tmp/=v;
alpar@2086:     return tmp;
alpar@2086:   }
alpar@2086:     
alpar@2086:   /// @}
alpar@2086: 
alpar@2086: } //END OF NAMESPACE LEMON
alpar@2086: 
alpar@2086: #endif // LEMON_POLYNOMIAL_H