#ifndef POLYNOMIAL_HPP #define POLYNOMIAL_HPP #include <algebra/TinyVector.hpp> template <size_t N> class Polynomial { private: using Coefficients = TinyVector<N + 1, double>; Coefficients m_coefficients; public: PUGS_INLINE constexpr const TinyVector<N + 1>& coefficients() const { return m_coefficients; } PUGS_INLINE constexpr TinyVector<N + 1>& coefficients() { return m_coefficients; } PUGS_INLINE constexpr bool operator==(const Polynomial<N>& q) const { if (m_coefficients != q.coefficients()) { return false; } return true; } PUGS_INLINE constexpr bool operator!=(const Polynomial<N>& q) const { return not this->operator==(q); } PUGS_INLINE constexpr Polynomial<N> operator+(const Polynomial<N>& q) const { TinyVector<N + 1> sum_coefs = m_coefficients + q.coefficients(); Polynomial<N> S(sum_coefs); return S; } PUGS_INLINE constexpr Polynomial<N> operator-(const Polynomial<N>& q) const { TinyVector<N + 1> diff_coefs = m_coefficients - q.coefficients(); Polynomial<N> D(diff_coefs); return D; } PUGS_INLINE constexpr Polynomial<N> operator*(const double& lambda) const { TinyVector<N + 1> mult_coefs = lambda * m_coefficients; Polynomial<N> M(mult_coefs); return M; } PUGS_INLINE constexpr friend Polynomial<N> operator*(const double& lambda, const Polynomial<N> P) { return P * lambda; } // evaluation using Horner's method https://en.wikipedia.org/wiki/Horner's_method PUGS_INLINE constexpr double evaluate(const double& x) const { TinyVector<N + 1> coefs = this->coefficients(); double bcoef = coefs[N]; for (size_t i = N; i > 0; --i) { bcoef *= x; bcoef += coefs[i - 1]; } return bcoef; } PUGS_INLINE constexpr double operator()(const double x) const { TinyVector<N + 1> coefs = this->coefficients(); double bcoef = coefs[N]; for (size_t i = N; i > 0; --i) { bcoef *= x; bcoef += coefs[i - 1]; } return bcoef; } PUGS_INLINE constexpr friend Polynomial<N + 1> primitive(const Polynomial<N>& P) { TinyVector<N + 2> coefs; for (size_t i = 0; i < N + 1; ++i) { coefs[i + 1] = P.coefficients()[i] / double(i + 1); } coefs[0] = 0; return Polynomial<N + 1>{coefs}; } PUGS_INLINE constexpr friend std::ostream& operator<<(std::ostream& os, const Polynomial<N>& P) { // os << "P(x) = "; bool all_coef_zero = true; for (size_t i = N; i > 1; --i) { if (P.coefficients()[i] != 0.) { if (P.coefficients()[i] != 1 && P.coefficients()[i] != -1) { os << std::abs(P.coefficients()[i]); } os << "x^" << i; if (P.coefficients()[i - 1] > 0.) os << " + "; else os << " - "; all_coef_zero = false; } else { os << " "; } } if (P.coefficients()[1] != 0.) { if (P.coefficients()[1] != 1 && P.coefficients()[1] != -1) { os << std::abs(P.coefficients()[1]); } os << "x"; if (P.coefficients()[0] > 0.) os << " + "; else if (P.coefficients()[0] < 0.) os << " - "; all_coef_zero = false; } else { os << " "; } if (P.coefficients()[0] != 0. || all_coef_zero) os << std::abs(P.coefficients()[0]); return os; } PUGS_INLINE constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {} PUGS_INLINE constexpr Polynomial(TinyVector<N + 1>&& coefficients) noexcept : m_coefficients{coefficients} {} PUGS_INLINE constexpr Polynomial() noexcept = default; ~Polynomial() = default; }; #endif // POLYNOMIAL_HPP