#ifndef __JMATH__JLEGENDRE__
#define __JMATH__JLEGENDRE__

#include <istream>
#include <ostream>
#include <vector>
#include <limits>

#include"JMath/JZero.hh"
#include"JMath/JMathSupportkit.hh"


/**
 * \author mdejong
 */

namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  /**
   * Base class for Legendre polynome.
   */
  template<class JOrdinate_t>
  struct JLegendre_t :
    public std::vector<JOrdinate_t>
  {
    typedef typename std::vector<JOrdinate_t>::const_iterator           const_iterator;
    typedef typename std::vector<JOrdinate_t>::iterator                 iterator;
    typedef typename std::vector<JOrdinate_t>::const_reverse_iterator   const_reverse_iterator;
    typedef typename std::vector<JOrdinate_t>::reverse_iterator         reverse_iterator;


    /**
     * Default constructor.
     */
    JLegendre_t() :
      xmin(-1.0),
      xmax(+1.0)
    {}


    /**
     * Constructor.
     *
     * \param  xmin             minimal abscissa value
     * \param  xmax             maximal abscissa value
     */
    JLegendre_t(const double xmin,
		const double xmax) :
      xmin(xmin),
      xmax(xmax)
    {}


    /**
     * Virtual destructor.
     */
    virtual ~JLegendre_t()
    {}


    /**
     * Function value.
     *
     * \param  x                abscissa value
     * \return                  function value
     */
    virtual JOrdinate_t getValue(const double x) const = 0;


    /**
     * Function value.
     *
     * \param  x                abscissa value
     * \return                  function value
     */
    JOrdinate_t operator()(const double x) const
    {
      return this->getValue(x);
    }


    /**
     * Get minimal abscissa value.
     *
     * \return                  minimal abscissa value
     */
    double getXmin() const 
    {
      return xmin;
    }


    /**
     * Get maximal abscissa value.
     *
     * \return                  maximal abscissa value
     */
    double getXmax() const 
    {
      return xmax;
    }


    /**
     * Get normalised abscissa value.
     *
     * \param  x                abscissa value
     * \return                  normalised abscissa value
     */
    double getX(const double x) const 
    {
      return 2.0 * (x - xmin) / (xmax - xmin)  -  1.0;
    }


    /**
     * Reset.
     */
    void reset()
    {
      for (iterator i = this->begin(); i != this->end(); ++i) {
	*i = zero;
      }
    }


    /**
     * Set abscissa values.
     *
     * \param  xmin             minimal abscissa value
     * \param  xmax             maximal abscissa value
     */
    void set(const double xmin,
	     const double xmax)
    {
      this->xmin = xmin;
      this->xmax = xmax;
    }


    /**
     * Set abscissa values.
     *
     * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following data member:
     *   - double        %first;       //
     * which conforms with std::pair.
     *
     * \param  __begin          begin of data
     * \param  __end            end   of data
     */
    template<class T>
    void set(T __begin, T __end)
    {
      this->xmin = std::numeric_limits<double>::max();
      this->xmax = std::numeric_limits<double>::lowest();

      for (T i = __begin; i != __end; ++i) {
	if (i->first < this->xmin) { this->xmin = i->first; }
	if (i->first > this->xmax) { this->xmax = i->first; }
      }
    }


    /**
     * Read Legendre polynome from input.
     *
     * \param  in               input stream
     * \param  object           Legendre polynome
     * \return                  input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JLegendre_t& object)
    {
      object.clear();

      for (JOrdinate_t x; in >> x; ) {
        object.push_back(x);
      }

      return in;
    }


    /**
     * Write Legendre polynome to output.
     *
     * \param  out              output stream
     * \param  object           Legendre polynome
     * \return                  output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JLegendre_t& object)
    {
      for (typename JLegendre_t::const_iterator i = object.begin(); i != object.end(); ++i) {
        out << ' ' << *i;
      }

      return out;
    }


  protected:
    double xmin;                                      //!< minimal abscissa
    double xmax;                                      //!< maximal abscissa
  };


  /**
   * Template definition for function evaluation of Legendre polynome.
   */
  template<class JOrdinate_t, size_t N = (size_t) -1>
  struct JLegendre;
  

  /**
   * Template specialisation for function evaluation of of Legendre polynome for undefined number of degrees.
   */
  template<class JOrdinate_t>
  struct JLegendre<JOrdinate_t, (size_t) -1> :
    public JLegendre_t<JOrdinate_t>
  {
    /**
     * Default constructor.
     */
    JLegendre()
    {}


    /**
     * Constructor.
     *
     * \param  xmin             minimal abscissa value
     * \param  xmax             maximal abscissa value
     */
    JLegendre(const double xmin,
	      const double xmax)
    {
      this->set(xmin, xmax);
    }


    /**
     * Function value.
     *
     * \param  x                abscissa value
     * \return                  function value
     */
    virtual JOrdinate_t getValue(const double x) const override 
    {
      const double z = this->getX(x);
      JOrdinate_t  y = zero;

      for (size_t n = 0; n != this->size(); ++n) {
        y += (*this)[n] * legendre(n,z);
      }

      return y;
    }
  };


  /**
   * Template specialisation for function evaluation of of Legendre polynome for defined number of degrees.
   */
  template<class JOrdinate_t, size_t N>
  struct JLegendre :
    JLegendre<JOrdinate_t>
  {
    using JLegendre_t<JOrdinate_t>::set;


    /**
     * Default constructor.
     */
    JLegendre()
    {
      this->JLegendre<JOrdinate_t>::resize(N+1, getZero<JOrdinate_t>());
    }


    /**
     * Constructor.
     *
     * \param  xmin             minimal abscissa value
     * \param  xmax             maximal abscissa value
     */
    JLegendre(const double xmin,
	      const double xmax)
    {
      this->JLegendre<JOrdinate_t>::resize(N+1, getZero<JOrdinate_t>());
      this->set(xmin, xmax);
    }


    /**
     * Constructor.
     *
     * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following data members:
     *   - double      %first;       //
     *   - JOrdinate_t %second;      //
     * which conforms with std::pair.
     *
     * \param  __begin          begin of data
     * \param  __end            end   of data
     */
    template<class T>
    JLegendre(T __begin, T __end)
    {
      this->JLegendre<JOrdinate_t>::resize(N+1, getZero<JOrdinate_t>());
      this->set(__begin, __end);
    }


    /**
     * Set Legendre polynome.
     *
     * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following data members:
     *   - double      %first;       //
     *   - JOrdinate_t %second;      //
     * which conforms with std::pair.
     *
     * \param  __begin          begin of data
     * \param  __end            end   of data
     */
    template<class T>
    void set(T __begin, T __end)
    {
      this->reset();
      this->JLegendre<JOrdinate_t>::set(__begin, __end);

      for (size_t n = 0; n != N + 1; ++n) {

	JOrdinate_t V = zero;
	double      W = 0.0;

	for (T i = __begin; i != __end; ++i) {

	  const double       x = i->first;
	  const JOrdinate_t& y = i->second;
	  const double       z = this->getX(x);
	  const double       w = legendre(n, z);

	  V += w * (y - (*this)(x));
	  W += w * (w);
	}

	(*this)[n] = V / W;
      }
    }
    

    /**
     * Read Legendre polynome from input.
     *
     * \param  in               input stream
     * \param  object           Legendre polynome
     * \return                  input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JLegendre& object)
    {
      for (typename JLegendre::iterator i = object.begin(); i != object.end(); ++i) {
        in >> *i;
      }

      return in;
    }

  private:
    void clear();
    void resize();
    void erase();
    void push_back();
    void pop_back();
  };
}

#endif