From 00eea8e2d35fb2b6db845863b3b0169cd9cb1d2b Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Thu, 30 Jun 2022 18:28:14 +0200
Subject: [PATCH] Normal error fixed
---
src/analysis/PolynomialP.hpp | 154 ++++++++++++++++++++++++++---------
tests/test_PolynomialP.cpp | 2 +-
2 files changed, 116 insertions(+), 40 deletions(-)
diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index b76fe6559..bcb89042b 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -109,6 +109,7 @@ class PolynomialP
}
return *this;
}
+
PUGS_INLINE constexpr PolynomialP&
operator+=(const PolynomialP& Q)
{
@@ -170,6 +171,65 @@ class PolynomialP
return absolute_coefs[absolute_position];
}
+ PUGS_INLINE
+ constexpr double
+ operator[](const TinyVector<Dimension, size_t> relative_pos)
+ {
+ size_t total_degree = 0;
+ for (size_t i = 0; i < Dimension; ++i) {
+ Assert((relative_pos[i] <= N),
+ "You are looking for a coefficient of order greater than the degree of the polynomial");
+ total_degree += relative_pos[i];
+ }
+ Assert((total_degree <= N), "The sum of the degrees of the coefficient you are looking for is greater than the "
+ "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
+ absolute_position(const TinyVector<Dimension, size_t> relative_pos)
+ {
+ size_t total_degree = 0;
+ for (size_t i = 0; i < Dimension; ++i) {
+ Assert((relative_pos[i] <= N),
+ "You are looking for a coefficient of order greater than the degree of the polynomial");
+ total_degree += relative_pos[i];
+ }
+ Assert((total_degree <= N), "The sum of the degrees of the coefficient you are looking for is greater than the "
+ "degree of the polynomial");
+ size_t abs_pos = 0;
+ if constexpr (Dimension == 1) {
+ abs_pos = relative_pos[0];
+ } else if constexpr (Dimension == 2) {
+ size_t total_degree = relative_pos[0] + relative_pos[1];
+ abs_pos = 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 abs_pos;
+ }
+
PUGS_INLINE
constexpr double
operator()(const TinyVector<Dimension> x) const
@@ -187,23 +247,15 @@ class PolynomialP
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);
+ relative_pos = TinyVector<Dimension, 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);
+ relative_pos = TinyVector<Dimension, 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");
}
@@ -211,24 +263,63 @@ class PolynomialP
return value;
}
+ PUGS_INLINE size_t
+ find_size_coef(const size_t degree)
+ {
+ if constexpr (Dimension == 1) {
+ return degree + 1;
+ } else if constexpr (Dimension == 2) {
+ return (degree + 1) * (degree + 2) / 2;
+ } else {
+ static_assert(Dimension == 3);
+ return (degree + 1) * (degree + 2) * (degree + 3) / 6;
+ }
+ }
+
+ // PUGS_INLINE constexpr auto
+ // derivative(const PolynomialP& P, const size_t var)
+ // {
+ // TinyVector<size_coef> coefs(zero);
+ // PolynomialP<N, Dimension> Q(coefs);
+ // if constexpr (N == 0) {
+ // return Q;
+ // } else {
+ // Assert(var < Dimension,
+ // "You can not derive a polynomial with respect to a variable of rank greater than the dimension");
+ // if constexpr (Dimension == 1) {
+ // for (size_t i = 0; i < size_coef; ++i) {
+ // coefs[i] = double(i + 1) * P.coefficients()[i + 1];
+ // }
+ // return Q;
+ // } else if constexpr (Dimension == 2) {
+ // if (var == 0) {
+ // for (size_t i = 0; i < N; ++i) {
+ // for (size_t j = 0; j < N; ++i) {
+ // TinyVector<Dimension, size_t> relative_pos(i, j);
+ // TinyVector<Dimension, size_t> relative_posp(i + 1, j);
+ // Q[relative_pos] = double(i + 1) * P[relative_posp];
+ // }
+ // }
+ // } else {
+ // for (size_t i = 0; i < N; ++i) {
+ // for (size_t j = 0; j < N; ++i) {
+ // TinyVector<Dimension, size_t> relative_pos(i, j);
+ // TinyVector<Dimension, size_t> relative_posp(i, j + 1);
+ // Q[relative_pos] = double(j + 1) * P[relative_posp];
+ // }
+ // }
+ // }
+ // return Q;
+ // } else {
+ // throw UnexpectedError("Not yet Available in 3D");
+ // }
+ // }
+ // }
+
PUGS_INLINE constexpr PolynomialP() noexcept = default;
~PolynomialP() = default;
};
-// // 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;
-// }
-
// template <size_t M>
// PUGS_INLINE constexpr friend void
// divide(const PolynomialP<N>& P1, const PolynomialP<M>& P2, PolynomialP<N>& Q, PolynomialP<N>& R)
@@ -378,21 +469,6 @@ class PolynomialP
// return integral;
// }
-// template <size_t N>
-// PUGS_INLINE constexpr auto
-// derivative(const PolynomialP<N>& P)
-// {
-// if constexpr (N == 0) {
-// return PolynomialP<0>(0);
-// } else {
-// TinyVector<N> coefs;
-// for (size_t i = 0; i < N; ++i) {
-// coefs[i] = double(i + 1) * P.coefficients()[i + 1];
-// }
-// return PolynomialP<N - 1>(coefs);
-// }
-// }
-
// template <size_t N>
// PUGS_INLINE constexpr PolynomialP<N>
// lagrangePolynomialP(const TinyVector<N + 1> zeros, const size_t k)
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
index d6323a37b..1855c7b67 100644
--- a/tests/test_PolynomialP.cpp
+++ b/tests/test_PolynomialP.cpp
@@ -116,7 +116,7 @@ TEST_CASE("PolynomialP", "[analysis]")
REQUIRE(Py(pos) == 11);
REQUIRE(Px(pos2) == 10);
REQUIRE(P(pos2) == 62);
- REQUIRE(Q(pos) == -20);
+ REQUIRE(Q(pos) == -24);
}
// // SECTION("product")
--
GitLab