#ifndef __JMATH__JPOLYNOME__
#define __JMATH__JPOLYNOME__

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

#include "JLang/JEquals.hh"
#include "JLang/JManip.hh"
#include "JMath/JMath.hh"

/**
 * \author mdejong
 */

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

namespace JMATH {

  using JLANG::JEquals;

  /**
   * Polynome model.
   */
  struct JPolynome_t : 
    public std::vector<double>,
    public JMath<JPolynome_t>,
    public JEquals<JPolynome_t>
  {
    /**
     * Default constructor.
     */
    JPolynome_t()
    {}

    
    /**
     * Equality.
     *
     * \param  P             polynome
     * \param  eps           numerical precision
     * \return               true if polynomes identical; else false
     */
    bool equals(const JPolynome_t& P,
		const double       eps = std::numeric_limits<double>::min()) const
    {
      if (this->size() == P.size()) {

	for (const_iterator p = this->begin(), q = P.begin(); p != this->end(); ++p, ++q) {
	  if (fabs(*p - *q) > eps) {
	    return false;
	  }
	}

	return true;
	
      } else {

	return false;
      }
    }


    /**
     * Add polynome.
     *
     * \param  polynome      polynome
     * \return               this polynome
     */
    JPolynome_t& add(const JPolynome_t& polynome)
    {
      while (this->size() < polynome.size()) {
	this->push_back(0.0);
      }

      for (size_t i = 0; i != this->size(); ++i) {
	(*this)[i] += polynome[i];
      }

      return *this;
    }


    /**
     * Subtract polynome.
     *
     * \param  polynome      polynome
     * \return               this polynome
     */
    JPolynome_t& sub(const JPolynome_t& polynome)
    {
      while (this->size() < polynome.size()) {
	this->push_back(0.0);
      }

      for (size_t i = 0; i != this->size(); ++i) {
	(*this)[i] -= polynome[i];
      }

      return *this;
    }


    /**
     * Scale polynome.
     *
     * \param  factor        multiplication factor
     * \return               this polynome
     */
    JPolynome_t& mul(const double factor)
    {
      for (iterator i = begin(); i != end(); ++i) {
	(*i) *= factor;
      }

      return *this;
    }

    
    /**
     * Read polynome from input.
     *
     * \param  in            input stream
     * \param  object        polynome
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JPolynome_t& object)
    {
      for (double x; in >> x; ) {
	object.push_back(x);
      }

      return in;
    }


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

      return out;
    }
  };


  /**
   * Polynome function object.
   *
   * Evaluates function, derivative, integral and gradient values.
   */
  struct JPolynome : 
    public JPolynome_t
  {
    /**
     * Default constructor.
     */
    JPolynome()
    {}


    /**
     * Copy constructor.
     *
     * \param  polynome      polynome
     */
    JPolynome(const JPolynome_t& polynome) :
      JPolynome_t(polynome)
    {}


    /**
     * Constructor.
     *
     * \param  __begin       begin of data
     * \param  __end         end   of data
     */
    template<class T>
    JPolynome(T __begin,
	      T __end)
    {
      for (T i = __begin; i != __end; ++i) {
	push_back(*i);
      }
    }


    /**
     * Initialise constructor.
     *
     * \param  args          values
     */
    template<class ...Args>
    JPolynome(const Args& ...args)
    {
      set(args...);
    }


    /**
     * Set values.
     *
     * \param  args          values
     */
    template<class ...Args>
    JPolynome& set(const Args& ...args)
    {
      clear();

      __set__(args...);

      return *this;
    }


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      double y = 0.0;
      double z = 1.0;

      for (const_iterator i = begin(); i != end(); ++i, z *= x) {
	y += (*i) * z;
      }

      return y;
    }


    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      double y = 0.0;

      if (!empty()) {

	double z = 1.0;
	int    n = 1;

	for (const_iterator i = begin(); ++i != end(); ++n, z *= x) {
	  y += (*i) * z * n;
	}
      }

      return y;
    }


    /**
     * Integral value.
     *
     * \param  x             abscissa value
     * \return               integral value
     */
    double getIntegral(const double x) const
    {
      double y = 0.0;
      double z = x;
      int    n = 1;

      for (const_iterator i = begin(); i != end(); ++i, ++n, z *= x) {
	y += (*i) * z / n;
      }

      return y;
    }


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


    /**
     * Get derivative function.
     *
     * \return               derivative function
     */
    JPolynome getDerivative() const
    {
      JPolynome buffer;

      if        (size() == 0u) {
      } else if (size() == 1u) {

	buffer.push_back(0.0);

      } else {

	int n = 1;

	for (const_iterator i = begin(); ++i != end(); ++n) {
	  buffer.push_back(n * (*i));
	}
      }

      return buffer;
    }


    /**
     * get integral function.
     *
     * \return               integral function
     */
    JPolynome getIntegral() const
    {
      JPolynome buffer(0.0);

      int n = 1;

      for (const_iterator i = begin(); i != end(); ++i, ++n) {
	buffer.push_back((*i) / (double) n);
      }

      return buffer;
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JPolynome_t getGradient(const double x) const
    {
      JPolynome_t buffer;

      double z = 1.0;

      for (const_iterator i = begin(); i != end(); ++i, z *= x) {
	buffer.push_back(z);                                        // d(f)/d(x_i)
      }

      return buffer;
    }

  private:
    /**
     * Recursive method for setting values.
     *
     * \param  x             next value
     * \param  args          remaining values
     */
    template<class ...Args>
    void __set__(const double x, const Args& ...args)
    {
      this->push_back(x);

      __set__(args...);
    }


    /**
     * Termination method for setting values.
     */
    void __set__() const
    {}    
  };
}

#endif