From 0e5133223d7aded0d81d22c3330807d203c0a3aa Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Thu, 30 Jun 2022 15:58:58 +0200
Subject: [PATCH] dd evaluation feature with multi-D Horner method
---
src/analysis/PolynomialP.hpp | 57 +++++++++++++++++++++++++++++++++++-
tests/test_PolynomialP.cpp | 42 +++++++++++++++++++++-----
2 files changed, 90 insertions(+), 9 deletions(-)
diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index 5712e5ef5..b76fe6559 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -2,6 +2,7 @@
#define POLYNOMIALP_HPP
#include <algebra/TinyVector.hpp>
+#include <utils/Exceptions.hpp>
template <size_t N, size_t Dimension>
class PolynomialP
@@ -141,7 +142,7 @@ class PolynomialP
PUGS_INLINE
constexpr double
- coef(const TinyVector<Dimension, size_t> relative_pos) const
+ operator[](const TinyVector<Dimension, size_t> relative_pos) const
{
size_t total_degree = 0;
for (size_t i = 0; i < Dimension; ++i) {
@@ -153,9 +154,63 @@ class PolynomialP
"degree of the polynomial");
TinyVector<size_coef> absolute_coefs = this->coefficients();
size_t absolute_position = 0;
+ if constexpr (Dimension == 1) {
+ absolute_position = relative_pos[0];
+ } else if constexpr (Dimension == 2) {
+ size_t total_degree = relative_pos[0] + relative_pos[1];
+ absolute_position = total_degree * (total_degree + 1) / 2 + relative_pos[1];
+ } else {
+ throw UnexpectedError("Not yet Available in 3D");
+ // static_assert(Dimension == 3);
+ // size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
+ // return (total_degree + 1) * (total_degree + 2) * (total_degree + 3) / 6 + relative_pos[1];
+ // return (N + 1) * (N + 2) * (N + 3) / 6;
+ }
+
return absolute_coefs[absolute_position];
}
+ PUGS_INLINE
+ constexpr double
+ operator()(const TinyVector<Dimension> x) const
+ {
+ const auto& P = *this;
+ double value = 0.;
+ if constexpr (Dimension == 1) {
+ value = m_coefficients[N];
+ for (size_t i = N; i > 0; --i) {
+ value *= x;
+ value += m_coefficients[i - 1];
+ }
+ } else if constexpr (Dimension == 2) {
+ TinyVector<Dimension, size_t> relative_pos(0, N);
+ value = P[relative_pos];
+ for (size_t i = N; i > 0; --i) {
+ value *= x[1];
+ relative_pos = TinyVector<2, size_t>(N - i + 1, i - 1);
+ double valuex = P[relative_pos];
+ for (size_t j = N - i + 1; j > 0; --j) {
+ valuex *= x[0];
+ relative_pos = TinyVector<2, size_t>(j - 1, i - 1);
+ valuex += P[relative_pos];
+ }
+ value += valuex;
+ }
+ // relative_pos = TinyVector<2, size_t>(N, 0);
+ // double valuex = P[relative_pos];
+ // for (size_t j = N; j > 0; --j) {
+ // valuex *= x[0];
+ // relative_pos = TinyVector<2, size_t>(j - 1, 0);
+ // valuex += P[relative_pos];
+ // }
+ // value += valuex;
+ } else {
+ throw UnexpectedError("Not yet Available in 3D");
+ }
+
+ return value;
+ }
+
PUGS_INLINE constexpr PolynomialP() noexcept = default;
~PolynomialP() = default;
};
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
index 1a9b6dbc6..d6323a37b 100644
--- a/tests/test_PolynomialP.cpp
+++ b/tests/test_PolynomialP.cpp
@@ -87,11 +87,43 @@ TEST_CASE("PolynomialP", "[analysis]")
REQUIRE(Q == (2 * P));
}
- // SECTION("product")
+ SECTION("access_coef")
+ {
+ TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+ TinyVector<10> coef2(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+ PolynomialP<2, 2> P(coef);
+ PolynomialP<3, 2> Q(coef2);
+ TinyVector<2, size_t> relative_coef(1, 1);
+ TinyVector<2, size_t> relative_coef2(1, 2);
+ REQUIRE(P[relative_coef] == 2);
+ REQUIRE(Q[relative_coef2] == 1);
+ REQUIRE(Q[relative_coef] == 3);
+ }
+
+ SECTION("evaluation")
+ {
+ TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+ TinyVector<10> coef2(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+ PolynomialP<2, 2> P(coef);
+ PolynomialP<3, 2> Q(coef2);
+ TinyVector<6> coefx(1, -2, 0, 7, 0, 0);
+ TinyVector<6> coefy(2, 0, -2, 0, 0, 7);
+ TinyVector<2> pos(1, -1);
+ TinyVector<2> pos2(-1, 2);
+ PolynomialP<2, 2> Px(coefx);
+ PolynomialP<2, 2> Py(coefy);
+ REQUIRE(Px(pos) == 6);
+ REQUIRE(Py(pos) == 11);
+ REQUIRE(Px(pos2) == 10);
+ REQUIRE(P(pos2) == 62);
+ REQUIRE(Q(pos) == -20);
+ }
+
+ // // SECTION("product")
// {
// Polynomial<2> P(2, 3, 4);
// Polynomial<3> Q(1, 2, -1, 1);
- // Polynomial<4> R;
+ // Polynomial<4 > R;
// Polynomial<5> S;
// R = P;
// S = P;
@@ -118,12 +150,6 @@ TEST_CASE("PolynomialP", "[analysis]")
// REQUIRE(Polynomial<2>{0, 1, 0} == R);
// }
- // SECTION("evaluation")
- // {
- // Polynomial<2> P(2, -3, 4);
- // REQUIRE(P(3) == 29);
- // }
-
// SECTION("primitive")
// {
// Polynomial<2> P(2, -3, 4);
--
GitLab