#ifndef __JMATH__JPOLYNOME__ #define __JMATH__JPOLYNOME__ #include #include #include #include #include #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, public JMath, public JEquals { /** * 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::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 JPolynome(T __begin, T __end) { for (T i = __begin; i != __end; ++i) { push_back(*i); } } /** * Initialise constructor. * * \param args values */ template JPolynome(const Args& ...args) { set(args...); } /** * Set values. * * \param args values */ template 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 void __set__(const double x, const Args& ...args) { this->push_back(x); __set__(args...); } /** * Termination method for setting values. */ void __set__() const {} }; } #endif