diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp
index 4d9941996e8379d2aed73d368d8d75a4a55f7b13..252db7b8deab35c3aa5991c03a74e6aa28469a0e 100644
--- a/src/analysis/Polynomial.hpp
+++ b/src/analysis/Polynomial.hpp
@@ -9,6 +9,7 @@ class Polynomial
private:
using Coefficients = TinyVector<N + 1, double>;
Coefficients m_coefficients;
+ static_assert((N >= 0), "Polynomial degree must be non-negative");
public:
PUGS_INLINE
@@ -25,10 +26,12 @@ class Polynomial
return m_coefficients;
}
- PUGS_INLINE
- constexpr bool
- operator==(const Polynomial<N>& q) const
+ template <size_t M>
+ PUGS_INLINE constexpr bool
+ operator==(const Polynomial<M>& q) const
{
+ if (M != N)
+ return false;
if (m_coefficients != q.coefficients()) {
return false;
}
@@ -42,22 +45,55 @@ class Polynomial
return not this->operator==(q);
}
- PUGS_INLINE
- constexpr Polynomial<N>
- operator+(const Polynomial<N>& q) const
+ template <size_t M>
+ PUGS_INLINE constexpr Polynomial<std::max(M, N)>
+ operator+(const Polynomial<M>& Q) const
{
- TinyVector<N + 1> sum_coefs = m_coefficients + q.coefficients();
- Polynomial<N> S(sum_coefs);
- return S;
+ Polynomial<std::max(M, N)> P;
+ if constexpr (M > N) {
+ P.coefficients() = Q.coefficients();
+ for (size_t i = 0; i <= N; ++i) {
+ P.coefficients()[i] += coefficients()[i];
+ }
+ } else {
+ P.coefficients() = coefficients();
+ for (size_t i = 0; i <= M; ++i) {
+ P.coefficients()[i] += Q.coefficients()[i];
+ }
+ }
+ return P;
}
- PUGS_INLINE
- constexpr Polynomial<N>
- operator-(const Polynomial<N>& q) const
+ template <size_t M>
+ PUGS_INLINE constexpr Polynomial<std::max(M, N)>
+ operator-(const Polynomial<M>& Q) const
+ {
+ Polynomial<std::max(M, N)> P;
+ if constexpr (M > N) {
+ P.coefficients() = -Q.coefficients();
+ for (size_t i = 0; i <= N; ++i) {
+ P.coefficients()[i] += coefficients()[i];
+ }
+ } else {
+ P.coefficients() = coefficients();
+ for (size_t i = 0; i <= M; ++i) {
+ P.coefficients()[i] -= Q.coefficients()[i];
+ }
+ }
+ return P;
+ }
+
+ template <size_t M>
+ PUGS_INLINE constexpr Polynomial<M + N> operator*(const Polynomial<M>& Q) const
{
- TinyVector<N + 1> diff_coefs = m_coefficients - q.coefficients();
- Polynomial<N> D(diff_coefs);
- return D;
+ Polynomial<M + N> P;
+ P.coefficients() = zero;
+ for (size_t i = 0; i <= N; ++i) {
+ for (size_t j = 0; j <= M; ++j) {
+ P.coefficients()[i + j] += coefficients()[i] * Q.coefficients()[j];
+ }
+ }
+ return P;
}
PUGS_INLINE
@@ -112,47 +148,89 @@ class Polynomial
return Polynomial<N + 1>{coefs};
}
+ PUGS_INLINE
+ constexpr friend double
+ integrate(const Polynomial<N>& P, const double& xinf, const double& xsup)
+ {
+ Polynomial<N + 1> Q = primitive(P);
+ return (Q(xsup) - Q(xinf));
+ }
+
+ PUGS_INLINE
+ constexpr friend auto
+ derivative(const Polynomial<N>& P)
+ {
+ if constexpr (N == 0) {
+ return Polynomial<0>(0);
+ } else {
+ TinyVector<N> coefs;
+ for (size_t i = 0; i < N; ++i) {
+ coefs[i] = double(i + 1) * P.coefficients()[i + 1];
+ }
+ 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 (N == 0) {
+ os << P.coefficients()[0];
+ return os;
+ }
+ if (N != 1) {
+ if (P.coefficients()[N] != 0.) {
+ if (P.coefficients()[N] < 0.) {
+ os << "- ";
+ }
+ if (P.coefficients()[N] != 1 && P.coefficients()[N] != -1) {
+ os << std::abs(P.coefficients()[N]);
+ }
+ os << "x^" << N;
+ all_coef_zero = false;
+ }
+ }
+ for (size_t i = N - 1; i > 1; --i) {
if (P.coefficients()[i] != 0.) {
+ if (P.coefficients()[i] > 0.) {
+ os << " + ";
+ } else if (P.coefficients()[i] < 0.) {
+ os << " - ";
+ }
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] > 0. && N != 1) {
+ os << " + ";
+ } else if (P.coefficients()[1] < 0.) {
+ os << " - ";
+ }
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)
+ if (P.coefficients()[0] != 0. || all_coef_zero) {
+ if (P.coefficients()[0] > 0.) {
+ os << " + ";
+ } else if (P.coefficients()[0] < 0.) {
+ os << " - ";
+ }
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(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {}
PUGS_INLINE
constexpr Polynomial(TinyVector<N + 1>&& coefficients) noexcept : m_coefficients{coefficients} {}
diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp
index 214be6ca85555730783975514172c0bfa7edc1a5..b5570d9dee27d3d182003f6120e17db013bc3acc 100644
--- a/tests/test_Polynomial.cpp
+++ b/tests/test_Polynomial.cpp
@@ -9,10 +9,12 @@
#include <analysis/Polynomial.hpp>
// Instantiate to ensure full coverage is performed
+template class Polynomial<0>;
template class Polynomial<1>;
template class Polynomial<2>;
template class Polynomial<3>;
template class Polynomial<4>;
+template class Polynomial<5>;
// clazy:excludeall=non-pod-global-static
@@ -26,15 +28,20 @@ TEST_CASE("Polynomial", "[analysis]")
{
Polynomial<2> P({2, 3, 4});
Polynomial<2> Q({2, 3, 4});
+ Polynomial<2> S({2, 3, 5});
REQUIRE(P == Q);
+ REQUIRE(P != S);
}
SECTION("addition")
{
Polynomial<2> P({2, 3, 4});
Polynomial<2> Q({-1, -3, 2});
Polynomial<2> S({1, 0, 6});
+ Polynomial<3> T({0, 3, 1, -2});
+ Polynomial<3> U({2, 6, 5, -2});
REQUIRE(S == (P + Q));
+ REQUIRE((T + P) == U);
}
SECTION("difference")
{
@@ -42,6 +49,9 @@ TEST_CASE("Polynomial", "[analysis]")
Polynomial<2> Q({3, 4, 5});
Polynomial<2> D({-1, -1, -1});
REQUIRE(D == (P - Q));
+ Polynomial<3> R({2, 3, 4, 1});
+ REQUIRE(D == (P - Q));
+ REQUIRE((P - R) == Polynomial<3>({0, 0, 0, -1}));
}
SECTION("product_by_scalar")
{
@@ -50,6 +60,12 @@ TEST_CASE("Polynomial", "[analysis]")
REQUIRE(M == (P * 3));
REQUIRE(M == (3 * P));
}
+ SECTION("product")
+ {
+ Polynomial<2> P({2, 3, 4});
+ Polynomial<3> Q({1, 2, -1, 1});
+ REQUIRE(Polynomial<5>({2, 7, 8, 7, -1, 4}) == (P * Q));
+ }
SECTION("evaluation")
{
Polynomial<2> P({2, -3, 4});
@@ -67,4 +83,23 @@ TEST_CASE("Polynomial", "[analysis]")
Polynomial<3> R({0, 2, -3. / 2, 4. / 3});
REQUIRE(Q == R);
}
+ SECTION("integrate")
+ {
+ Polynomial<2> P({2, -3, 3});
+ double xinf = -1;
+ double xsup = 1;
+ double result = integrate(P, xinf, xsup);
+ REQUIRE(result == 6);
+ }
+ SECTION("derivative")
+ {
+ Polynomial<2> P({2, -3, 3});
+ Polynomial<1> Q = derivative(P);
+ REQUIRE(Q == Polynomial<1>({-3, 6}));
+
+ Polynomial<0> P2(3);
+
+ Polynomial<0> R(0);
+ REQUIRE(derivative(P2) == R);
+ }
}