#ifndef POLYNOMIAL_1D_HPP #define POLYNOMIAL_1D_HPP #include <utils/Array.hpp> #include <utils/Exceptions.hpp> #include <algorithm> class [[nodiscard]] Polynomial1D { private: std::vector<double> m_coefficients; PUGS_INLINE size_t _getRealDegree() const { size_t real_degree = this->degree(); while (real_degree > 0 and std::abs(m_coefficients[real_degree]) < 1E-14) { --real_degree; } return real_degree; } PUGS_INLINE friend Polynomial1D _simplify(Polynomial1D && p) { // return std::move(p); size_t real_degree = p._getRealDegree(); if (real_degree != p.degree()) { Polynomial1D q(real_degree); for (size_t i = 0; i <= real_degree; ++i) { q.coefficient(i) = p.coefficient(i); } return q; } else { return std::move(p); } } public: friend std::ostream& operator<<(std::ostream& os, const Polynomial1D& p) { bool has_written = false; for (size_t i = 0; i < p.m_coefficients.size(); ++i) { const double& coef = p.m_coefficients[i]; // if (coef != 0) { if (coef > 0) { if (has_written) { os << " + "; } os << coef; } else { if (has_written) { os << " - " << -coef; } else { os << coef; } } if (i > 0) { os << "*x"; if (i > 1) { os << '^' << i; } } has_written = true; // } } if (not has_written) { os << 0; } return os; } PUGS_INLINE Polynomial1D operator-() const { Polynomial1D opposite(this->degree()); for (size_t i = 0; i < m_coefficients.size(); ++i) { opposite.m_coefficients[i] = -m_coefficients[i]; } return _simplify(std::move(opposite)); } PUGS_INLINE friend Polynomial1D operator*(const double a, const Polynomial1D& p) { Polynomial1D product(p.degree()); for (size_t i = 0; i < p.m_coefficients.size(); ++i) { product.m_coefficients[i] = a * p.m_coefficients[i]; } return _simplify(std::move(product)); } PUGS_INLINE Polynomial1D operator*(const Polynomial1D& p) const { Polynomial1D product(this->degree() + p.degree()); std::fill(begin(product.m_coefficients), end(product.m_coefficients), 0); for (size_t i = 0; i < this->m_coefficients.size(); ++i) { for (size_t j = 0; j < p.m_coefficients.size(); ++j) { product.m_coefficients[i + j] += this->m_coefficients[i] * p.m_coefficients[j]; } } return _simplify(std::move(product)); } PUGS_INLINE Polynomial1D operator+(const Polynomial1D& p) const { auto sum_left_is_bigger = [](const Polynomial1D& greater, const Polynomial1D& smaller) { Polynomial1D result(greater.degree()); for (size_t i = 0; i < greater.m_coefficients.size(); ++i) { result.m_coefficients[i] = greater.m_coefficients[i]; } for (size_t i = 0; i < smaller.m_coefficients.size(); ++i) { result.m_coefficients[i] += smaller.m_coefficients[i]; } return _simplify(std::move(result)); }; if (m_coefficients.size() >= p.m_coefficients.size()) { return sum_left_is_bigger(*this, p); } else { return sum_left_is_bigger(p, *this); } } PUGS_INLINE Polynomial1D operator-(const Polynomial1D& p) const { if (m_coefficients.size() >= p.m_coefficients.size()) { Polynomial1D result(this->degree()); for (size_t i = 0; i < p.m_coefficients.size(); ++i) { result.m_coefficients[i] = m_coefficients[i] - p.m_coefficients[i]; } for (size_t i = p.m_coefficients.size(); i < m_coefficients.size(); ++i) { result.m_coefficients[i] = m_coefficients[i]; } return _simplify(std::move(result)); } else { Polynomial1D result(p.degree()); for (size_t i = 0; i < m_coefficients.size(); ++i) { result.m_coefficients[i] = m_coefficients[i] - p.m_coefficients[i]; } for (size_t i = m_coefficients.size(); i < p.m_coefficients.size(); ++i) { result.m_coefficients[i] = -p.m_coefficients[i]; } return _simplify(std::move(result)); } } PUGS_INLINE Polynomial1D operator%(const Polynomial1D& q) const { const Polynomial1D& p = *this; Polynomial1D ratio(this->degree()); std::fill(begin(ratio.m_coefficients), end(ratio.m_coefficients), 0); Polynomial1D remaining(this->degree()); for (size_t i = 0; i < m_coefficients.size(); ++i) { remaining.m_coefficients[i] = m_coefficients[i]; } const size_t p_degree = p.degree(); const size_t q_degree = q._getRealDegree(); for (ssize_t i = p_degree - q_degree; i >= 0; --i) { ratio.m_coefficients[i] = remaining.m_coefficients[q_degree + i] / q.m_coefficients[q_degree]; for (ssize_t j = q_degree + i; j >= i; --j) { remaining.m_coefficients[j] -= ratio.m_coefficients[i] * q.m_coefficients[j - i]; } } if (q_degree == 0) { remaining.m_coefficients.resize(1); remaining.m_coefficients[0] = 0; } else { remaining.m_coefficients.resize(std::max(q_degree - 1, 1ul)); } // for (size_t j = q_degree; j <= remaining.degree(); ++j) { // remaining.m_coefficients[j] = 0; // } return _simplify(std::move(remaining)); } PUGS_INLINE Polynomial1D operator/(const Polynomial1D& q) const { const Polynomial1D& p = *this; Polynomial1D ratio(this->degree()); std::fill(begin(ratio.m_coefficients), end(ratio.m_coefficients), 0); Polynomial1D remaining(this->degree()); for (size_t i = 0; i < m_coefficients.size(); ++i) { remaining.m_coefficients[i] = m_coefficients[i]; } const size_t p_degree = p.degree(); const size_t q_degree = q._getRealDegree(); for (ssize_t i = p_degree - q_degree; i >= 0; --i) { ratio.m_coefficients[i] = remaining.m_coefficients[q_degree + i] / q.m_coefficients[q_degree]; for (ssize_t j = q_degree + i; j >= i; --j) { remaining.m_coefficients[j] -= ratio.m_coefficients[i] * q.m_coefficients[j - i]; } } return _simplify(std::move(ratio)); } PUGS_INLINE size_t degree() const { Assert(m_coefficients.size() > 0); return m_coefficients.size() - 1; } PUGS_INLINE double& coefficient(const size_t i) { Assert(i < m_coefficients.size(), "invalid coefficient number"); return m_coefficients[i]; } PUGS_INLINE const double& coefficient(const size_t i) const { Assert(i < m_coefficients.size(), "invalid coefficient number"); return m_coefficients[i]; } PUGS_INLINE friend Polynomial1D derive(const Polynomial1D& p) { if (p.degree() > 0) { Polynomial1D derivative{p.degree() - 1}; for (size_t i = 0; i <= derivative.degree(); ++i) { derivative.m_coefficients[i] = (i + 1) * p.m_coefficients[i + 1]; } return derivative; } else { return Polynomial1D(std::vector<double>{0}); } } PUGS_INLINE friend Polynomial1D primitive(const Polynomial1D& p) { if (p.degree() > 0) { Polynomial1D primitive{p.degree() + 1}; primitive.m_coefficients[0] = 0; for (size_t i = 0; i <= p.degree(); ++i) { primitive.m_coefficients[i + 1] = p.m_coefficients[i] / (i + 1); } return _simplify(std::move(primitive)); } else { return Polynomial1D(std::vector<double>{0}); } } PUGS_INLINE double operator()(const double x) const { double value = m_coefficients[this->degree()]; for (ssize_t i = this->degree() - 1; i >= 0; --i) { value *= x; value += m_coefficients[i]; } return value; } Polynomial1D& operator=(Polynomial1D&&) = default; Polynomial1D& operator=(const Polynomial1D&) = default; PUGS_INLINE explicit Polynomial1D(const std::vector<double>& coeffients) : m_coefficients(coeffients.size()) { Assert(coeffients.size() > 0, "empty list is not allowed"); for (size_t i = 0; i < coeffients.size(); ++i) { m_coefficients[i] = coeffients[i]; } size_t real_degree = this->_getRealDegree(); if (real_degree + 1 < m_coefficients.size()) { std::vector<double> real_coefficients(real_degree + 1); for (size_t i = 0; i <= real_degree; ++i) { real_coefficients[i] = m_coefficients[i]; } m_coefficients = std::move(real_coefficients); } } PUGS_INLINE explicit Polynomial1D(const size_t degree) : m_coefficients(degree + 1) { std::fill(m_coefficients.begin(), m_coefficients.end(), 0); } Polynomial1D(const Polynomial1D&) = default; Polynomial1D(Polynomial1D &&) = default; ~Polynomial1D() = default; }; #endif // POLYNOMIAL_1D_HPP