From 24233ee85ef7e395278f67e1d4131997f9b3e8d8 Mon Sep 17 00:00:00 2001
From: LABOURASSE Emmanuel <labourassee@gmail.com>
Date: Wed, 22 Jul 2020 18:19:22 +0200
Subject: [PATCH] debut du codage de la classe polynome
---
src/analysis/Polynomial.hpp | 165 ++++++++++++++++++++++++++++++++++++
tests/test_Polynomial.cpp | 70 +++++++++++++++
2 files changed, 235 insertions(+)
create mode 100644 src/analysis/Polynomial.hpp
create mode 100644 tests/test_Polynomial.cpp
diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp
new file mode 100644
index 000000000..4d9941996
--- /dev/null
+++ b/src/analysis/Polynomial.hpp
@@ -0,0 +1,165 @@
+#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
diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp
new file mode 100644
index 000000000..214be6ca8
--- /dev/null
+++ b/tests/test_Polynomial.cpp
@@ -0,0 +1,70 @@
+#include <catch2/catch.hpp>
+
+#include <Kokkos_Core.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/Types.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <analysis/Polynomial.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class Polynomial<1>;
+template class Polynomial<2>;
+template class Polynomial<3>;
+template class Polynomial<4>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Polynomial", "[analysis]")
+{
+ SECTION("construction")
+ {
+ REQUIRE_NOTHROW(Polynomial<2>{{2, 3, 4}});
+ }
+ SECTION("equality")
+ {
+ Polynomial<2> P({2, 3, 4});
+ Polynomial<2> Q({2, 3, 4});
+
+ REQUIRE(P == Q);
+ }
+ SECTION("addition")
+ {
+ Polynomial<2> P({2, 3, 4});
+ Polynomial<2> Q({-1, -3, 2});
+ Polynomial<2> S({1, 0, 6});
+ REQUIRE(S == (P + Q));
+ }
+ SECTION("difference")
+ {
+ Polynomial<2> P({2, 3, 4});
+ Polynomial<2> Q({3, 4, 5});
+ Polynomial<2> D({-1, -1, -1});
+ REQUIRE(D == (P - Q));
+ }
+ SECTION("product_by_scalar")
+ {
+ Polynomial<2> P({2, 3, 4});
+ Polynomial<2> M({6, 9, 12});
+ REQUIRE(M == (P * 3));
+ REQUIRE(M == (3 * P));
+ }
+ SECTION("evaluation")
+ {
+ Polynomial<2> P({2, -3, 4});
+ double result = P.evaluate(3);
+ REQUIRE(result == 29);
+ result = P(3);
+ REQUIRE(result == 29);
+ }
+ SECTION("primitive")
+ {
+ Polynomial<2> P({2, -3, 4});
+ TinyVector<4> coefs = zero;
+ Polynomial<3> Q(coefs);
+ Q = primitive(P);
+ Polynomial<3> R({0, 2, -3. / 2, 4. / 3});
+ REQUIRE(Q == R);
+ }
+}
--
GitLab