diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp index b76fe65591afdd1c5c85e802c0eb0c628b340547..bcb89042b4b124a99e4e0ccbb4e0758690f26653 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 d6323a37bb5ece3fa9acc9f8e934658306b4e257..1855c7b671d4201f55e471f2655bd86e6a4f9c5e 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")