/* -*- C++ -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library

 *

 * Copyright (C) 2003-2006

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_BEZIER_H

#define LEMON_BEZIER_H

///\ingroup misc

///\file

///\brief A simple class implementing polynomials.

///

///\author Alpar Juttner

#include<vector>

namespace lemon {

  /// \addtogroup misc

  /// @{

  ///Simple polinomial class

  ///This class implements a polynomial where the coefficients are of

  ///type \c T.

  ///

  ///The coefficients are stored in an std::vector.

  template<class T>

  class Polynomial

  {

    std::vector<T> _coeff;

  public:

    ///Construct a polynomial of degree \c d.

    explicit Polynomial(int d=0) : _coeff(d+1) {}

    ///\e

    template<class U> Polynomial(const U &u) : _coeff(1,u) {}

    ///\e

    template<class U> Polynomial(const Polynomial<U> &u) : _coeff(u.deg()+1)

    {

      for(int i=0;i<(int)_coeff.size();i++) _coeff[i]=u[i];

    }

    ///Query the degree of the polynomial.

    ///Query the degree of the polynomial.

    ///\warning This number differs from real degree of the polinomial if

    ///the coefficient of highest degree is 0.

    int deg() const { return _coeff.size()-1; }

    ///Set the degree of the polynomial.

    ///Set the degree of the polynomial. In fact it resizes the

    ///coefficient vector.

    void deg(int d) { _coeff.resize(d+1);}

    ///Returns (as a reference) the coefficient of degree \c d.

    typename std::vector<T>::reference operator[](int d) { return _coeff[d]; }

    ///Returns (as a const reference) the coefficient of degree \c d.

    typename std::vector<T>::const_reference

    operator[](int d) const {return _coeff[d];}

    ///Substitute the value u into the polinomial.

    ///Substitute the value u into the polinomial.

    ///The calculation will be done using type \c R.

    ///The following examples shows the usage of the template parameter \c R.

    ///\code

    ///  Polynomial<dim2::Point<double> > line(1);

    ///  line[0]=dim2::Point<double>(12,25);

    ///  line[1]=dim2::Point<double>(2,7);

    ///  ...

    ///  dim2::Point<double> d = line.subst<dim2::Point<double> >(23.2);

    ///\endcode

    ///

    ///\code

    ///  Polynomial<double> p;

    ///  Polynomial<double> q;

    ///  ...

    ///  Polynomial<double> s = p.subst<Polynomial<double> >(q);

    ///\endcode

    template<class R,class U>

    R subst(const U &u) const

    {

      typename std::vector<T>::const_reverse_iterator i=_coeff.rbegin();

      R v=*i;

      for(++i;i!=_coeff.rend();++i) v=v*u+*i;

      return v;

    }

    ///Substitute the value u into the polinomial.

    ///Substitute the value u into the polinomial.

    ///The calculation will be done using type \c T

    ///(i.e. using the type of the coefficients.)

    template<class U>

    T operator()(const U &u) const 

    {

      return subst<T>(u);

    }

    ///Derivate the polynomial (in place)

    Polynomial &derivateMyself()

    {

      for(int i=1;i<(int)_coeff.size();i++) _coeff[i-1]=i*_coeff[i];

      _coeff.pop_back();

      return *this;

    }

    ///Return the derivate of the polynomial

    Polynomial derivate() const

    {

      Polynomial tmp(deg()-1);

      for(int i=1;i<(int)_coeff.size();i++) tmp[i-1]=i*_coeff[i];

      return tmp;

    }

    ///Integrate the polynomial (in place)

    Polynomial &integrateMyself()

    {

      _coeff.push_back(T());

      for(int i=_coeff.size()-1;i>0;i--) _coeff[i]=_coeff[i-1]/i;

      _coeff[0]=0;

      return *this;

    }

    ///Return the integrate of the polynomial

    Polynomial integrate() const

    {

      Polynomial tmp(deg()+1);

      tmp[0]=0;

      for(int i=0;i<(int)_coeff.size();i++) tmp[i+1]=_coeff[i]/(i+1);

      return tmp;

    }

    ///\e

    template<class U>

    Polynomial &operator+=(const Polynomial<U> &p)

    {

      if(p.deg()>deg()) _coeff.resize(p.deg()+1);

      for(int i=0;i<=(int)std::min(deg(),p.deg());i++)

	_coeff[i]+=p[i];

      return *this;

    }

    ///\e

    template<class U>

    Polynomial &operator-=(const Polynomial<U> &p)

    {

      if(p.deg()>deg()) _coeff.resize(p.deg()+1);

      for(int i=0;i<=std::min(deg(),p.deg());i++) _coeff[i]-=p[i];

      return *this;

    }

    ///\e

    template<class U>

    Polynomial &operator+=(const U &u)

    {

      _coeff[0]+=u;

      return *this;

    }

    ///\e

    template<class U>

    Polynomial &operator-=(const U &u)

    {

      _coeff[0]+=u;

      return *this;

    }

    ///\e

    template<class U>

    Polynomial &operator*=(const U &u)

    {

      for(typename std::vector<T>::iterator i=_coeff.begin();i!=_coeff.end();++i)

	*i*=u;

      return *this;

    }

    ///\e

    template<class U>

    Polynomial &operator/=(const U &u)

    {

      for(typename std::vector<T>::iterator i=_coeff.begin();i!=_coeff.end();++i)

	*i/=u;

      return *this;

    }

  };

  ///Equality comparison

  ///\relates Polynomial

  ///\warning Two polynomials are defined to be unequal if their degrees differ,

  ///even if the non-zero coefficients are the same.

  template<class U,class V>

  bool operator==(const Polynomial<U> &u,const Polynomial<V> &v)

  {

    if(u.deg()!=v.deg()) return false;

    for(int i=0;i<=u.deg();i++) if(u[i]!=v[i]) return false;

    return true;

  }

  ///Non-equality comparison

  ///\relates Polynomial

  ///\warning Two polynomials are defined to be unequal if their degrees differ,

  ///even if the non-zero coefficients are the same.

  template<class U,class V>

  bool operator!=(const Polynomial<U> &u,const Polynomial<V> &v)

  {

    return !(u==v);

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator+(const Polynomial<U> &u,const Polynomial<V> &v)

  {

    Polynomial<U> tmp=u;

    tmp+=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator-(const Polynomial<U> &u,const Polynomial<V> &v)

  {

    Polynomial<U> tmp=u;

    tmp-=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator*(const Polynomial<U> &u,const Polynomial<V> &v)

  {

    Polynomial<U> tmp(u.deg()+v.deg());

    for(int i=0;i<=v.deg();i++)

      for(int j=0;j<=u.deg();j++)

	tmp[i+j]+=v[i]*u[j];

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator+(const Polynomial<U> &u,const V &v)

  {

    Polynomial<U> tmp=u;

    tmp+=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator+(const V &v,const Polynomial<U> &u)

  {

    Polynomial<U> tmp=u;

    tmp+=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator-(const Polynomial<U> &u,const V &v)

  {

    Polynomial<U> tmp=u;

    tmp-=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U>

  Polynomial<U> operator-(const Polynomial<U> &u)

  {

    Polynomial<U> tmp(u.deg());

    for(int i=0;i<=u.deg();i++) tmp[i]=-u[i];

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator-(const V &v,const Polynomial<U> &u)

  {

    Polynomial<U> tmp=-u;

    tmp+=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator*(const Polynomial<U> &u,const V &v)

  {

    Polynomial<U> tmp=u;

    tmp*=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator*(const V &v,const Polynomial<U> &u)

  {

    Polynomial<U> tmp=u;

    tmp*=v;

    return tmp;

  }

  ///\e

  ///\relates Polynomial

  ///

  template<class U,class V>

  Polynomial<U> operator/(const Polynomial<U> &u,const V &v)

  {

    Polynomial<U> tmp=u;

    tmp/=v;

    return tmp;

  }

  /// @}

} //END OF NAMESPACE LEMON

#endif // LEMON_POLYNOMIAL_H