From 87e0224bcd9306a33d102523727d28f9623fced0 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Wed, 22 Jun 2022 21:53:06 +0200
Subject: [PATCH] Add several features to multivariate polynomials and tests

---
 src/analysis/PolynomialP.hpp | 257 +++++++++--------------------------
 tests/CMakeLists.txt         |   1 +
 tests/test_Polynomial.cpp    |   4 +
 tests/test_PolynomialP.cpp   |  39 +++---
 4 files changed, 88 insertions(+), 213 deletions(-)

diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index b3521136b..5712e5ef5 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -90,194 +90,76 @@ class PolynomialP
     return P;
   }
 
-  PUGS_INLINE constexpr PolynomialP() noexcept = default;
-  ~PolynomialP()                               = default;
-};
-
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<std::max(M, N)>
-//   operator-(const PolynomialP<M>& Q) const
-//   {
-//     PolynomialP<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 PolynomialP<N>&
-//   operator=(const PolynomialP<M>& Q)
-//   {
-//     coefficients() = zero;
-//     for (size_t i = N + 1; i <= M; ++i) {
-//       Assert(Q.coefficients()[i] == 0, "degree of polynomialP to small in assignation");
-//     }
-//     //    static_assert(N >= M, "degree of polynomialP to small in assignation");
-//     for (size_t i = 0; i <= std::min(M, N); ++i) {
-//       coefficients()[i] = Q.coefficients()[i];
-//     }
-//     return *this;
-//   }
-
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<N>&
-//   operator+=(const PolynomialP<M>& Q)
-//   {
-//     static_assert(N >= M, "PolynomialP degree to small in affectation addition");
-//     for (size_t i = 0; i <= M; ++i) {
-//       coefficients()[i] += Q.coefficients()[i];
-//     }
-//     return *this;
-//   }
-
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<N>&
-//   operator-=(const PolynomialP<M>& Q)
-//   {
-//     static_assert(N >= M, "PolynomialP degree to small in affectation addition");
-//     for (size_t i = 0; i <= M; ++i) {
-//       coefficients()[i] -= Q.coefficients()[i];
-//     }
-//     return *this;
-//   }
-
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<M + N>
-//   operator*(const PolynomialP<M>& Q) const
-//   {
-//     PolynomialP<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;
-//   }
-
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<N>&
-//   operator*=(const PolynomialP<M>& Q)
-//   {
-//     static_assert(N >= M, "Degree to small in affectation product");
-//     for (size_t i = N - M + 1; i <= N; ++i) {
-//       Assert(coefficients()[i] == 0, "Degree of affectation product greater than the degree of input polynomialP");
-//     }
-//     PolynomialP<N> P(zero);
-//     for (size_t i = 0; i <= N - M; ++i) {
-//       for (size_t j = 0; j <= M; ++j) {
-//         P.coefficients()[i + j] += coefficients()[i] * Q.coefficients()[j];
-//       }
-//     }
-//     coefficients() = P.coefficients();
-//     return *this;
-//   }
-
-//   PUGS_INLINE
-//   constexpr PolynomialP<N>
-//   operator*(const double& lambda) const
-//   {
-//     TinyVector<N + 1> mult_coefs = lambda * m_coefficients;
-//     PolynomialP<N> M(mult_coefs);
-//     return M;
-//   }
+  PUGS_INLINE constexpr PolynomialP
+  operator-(const PolynomialP Q) const
+  {
+    PolynomialP<N, Dimension> P(m_coefficients);
+    P = P + (-Q);
+    return P;
+  }
 
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<M * N>
-//   compose(const PolynomialP<M>& Q) const
-//   {
-//     PolynomialP<M * N> P;
-//     P.coefficients() = zero;
-//     PolynomialP<M * N> R;
-//     R.coefficients() = zero;
-
-//     for (size_t i = 0; i <= N; ++i) {
-//       R = Q.template pow<N>(i) * coefficients()[i];
-//       P += R;   // R;
-//     }
-//     return P;
-//   }
+  template <size_t M, size_t Dim>
+  PUGS_INLINE constexpr PolynomialP&
+  operator=(const PolynomialP<M, Dim>& Q)
+  {
+    coefficients() = zero;
+    for (size_t i = 0; i < size_coef; ++i) {
+      coefficients()[i] = Q.coefficients()[i];
+    }
+    return *this;
+  }
+  PUGS_INLINE constexpr PolynomialP&
+  operator+=(const PolynomialP& Q)
+  {
+    m_coefficients += Q.coefficients();
+    return *this;
+  }
 
-//   template <size_t M, size_t I>
-//   PUGS_INLINE constexpr PolynomialP<M * N>
-//   power(const PolynomialP<M>& Q) const
-//   {
-//     return coefficients()[I] * Q.template pow<N>(I);
-//   }
+  template <size_t M>
+  PUGS_INLINE constexpr PolynomialP&
+  operator-=(const PolynomialP& Q)
+  {
+    m_coefficients -= Q.coefficients();
+    return *this;
+  }
 
-//   template <size_t M, size_t... I>
-//   PUGS_INLINE constexpr PolynomialP<M * N>
-//   compose2(const PolynomialP<M>& Q, std::index_sequence<I...>) const
-//   {
-//     PolynomialP<M * N> P;
-//     P.coefficients() = zero;
-//     P                = (power<M, I>(Q) + ...);
-//     return P;
-//   }
+  PUGS_INLINE
+  constexpr PolynomialP
+  operator*(const double& lambda) const
+  {
+    TinyVector<size_coef> mult_coefs = lambda * m_coefficients;
+    PolynomialP<N, Dimension> M(mult_coefs);
+    return M;
+  }
 
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<M * N>
-//   compose2(const PolynomialP<M>& Q) const
-//   {
-//     PolynomialP<M * N> P;
-//     P.coefficients()    = zero;
-//     using IndexSequence = std::make_index_sequence<N + 1>;
-//     return compose2<M>(Q, IndexSequence{});
-//     //    for (size_t i = 0; i <= N; ++i) {
-//     //      P += Q.template pow<N>(i) * coefficients()[i];
-//     //    }
-//     return P;
-//   }
+  PUGS_INLINE
+  constexpr friend PolynomialP<N, Dimension>
+  operator*(const double& lambda, const PolynomialP<N, Dimension> P)
+  {
+    return P * lambda;
+  }
 
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<M * N>
-//   operator()(const PolynomialP<M>& Q) const
-//   {
-//     PolynomialP<M * N> P;
-//     P.coefficients() = zero;
-//     PolynomialP<M * N> R;
-//     R.coefficients() = zero;
-
-//     for (size_t i = 0; i <= N; ++i) {
-//       R = Q.template pow<N>(i) * coefficients()[i];
-//       P += R;   // R;
-//     }
-//     return P;
-//   }
+  PUGS_INLINE
+  constexpr double
+  coef(const TinyVector<Dimension, size_t> relative_pos) const
+  {
+    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;
+    return absolute_coefs[absolute_position];
+  }
 
-//   template <size_t M>
-//   PUGS_INLINE constexpr PolynomialP<M * N>
-//   pow(size_t power) const
-//   {
-//     Assert(power <= M, "You declared a polynomialP of degree too small for return of the pow function");
-//     PolynomialP<M * N> R;
-//     R.coefficients() = zero;
-//     if (power == 0) {
-//       R.coefficients()[0] = 1;
-//     } else {
-//       R = *this;
-//       for (size_t i = 1; i < power; ++i) {
-//         R = R * *this;
-//       }
-//     }
-//     return R;
-//   }
+  PUGS_INLINE constexpr PolynomialP() noexcept = default;
+  ~PolynomialP()                               = default;
+};
 
-//   PUGS_INLINE
-//   constexpr friend PolynomialP<N>
-//   operator*(const double& lambda, const PolynomialP<N> P)
-//   {
-//     return P * lambda;
-//   }
 //   // evaluation using Horner's method https://en.wikipedia.org/wiki/Horner's_method
 //   PUGS_INLINE
 //   constexpr double
@@ -292,19 +174,6 @@ class PolynomialP
 //     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;
-//   }
-
 //   template <size_t M>
 //   PUGS_INLINE constexpr friend void
 //   divide(const PolynomialP<N>& P1, const PolynomialP<M>& P2, PolynomialP<N>& Q, PolynomialP<N>& R)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index bf3e97872..a23fd187b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -117,6 +117,7 @@ add_executable (unit_tests
   test_ParseError.cpp
   test_PETScUtils.cpp
   test_Polynomial.cpp
+  test_PolynomialP.cpp
   test_PolynomialBasis.cpp
   test_PrimalToDiamondDualConnectivityDataMapper.cpp
   test_PrimalToDual1DConnectivityDataMapper.cpp
diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp
index e9f92dbea..2b44b74de 100644
--- a/tests/test_Polynomial.cpp
+++ b/tests/test_Polynomial.cpp
@@ -7,6 +7,7 @@
 
 #include <algebra/TinyMatrix.hpp>
 #include <analysis/Polynomial.hpp>
+#include <analysis/PolynomialP.hpp>
 
 // Instantiate to ensure full coverage is performed
 template class Polynomial<0>;
@@ -15,6 +16,9 @@ template class Polynomial<2>;
 template class Polynomial<3>;
 template class Polynomial<4>;
 template class Polynomial<5>;
+template class PolynomialP<0, 2>;
+template class PolynomialP<1, 2>;
+template class PolynomialP<3, 2>;
 
 // clazy:excludeall=non-pod-global-static
 
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
index 1f1f2e496..1a9b6dbc6 100644
--- a/tests/test_PolynomialP.cpp
+++ b/tests/test_PolynomialP.cpp
@@ -65,26 +65,27 @@ TEST_CASE("PolynomialP", "[analysis]")
     REQUIRE(-P == PolynomialP<2, 2>(coef2));
   }
 
-  // SECTION("difference")
-  // {
-  //   Polynomial<2> P(2, 3, 4);
-  //   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});
-  //   R -= P;
-  //   REQUIRE(R == Polynomial<3>(0, 0, 0, 1));
-  // }
+  SECTION("difference")
+  {
+    TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+    TinyVector<6> coef2(1, 2, 3, -2, -1, -3);
+    TinyVector<6> coef3(0, 0, 0, 6, 6, 9);
+    PolynomialP<2, 2> P(coef);
+    PolynomialP<2, 2> Q(coef2);
+    PolynomialP<2, 2> R(coef3);
+    R = P - Q;
+    REQUIRE(R == PolynomialP<2, 2>(coef3));
+  }
 
-  // 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("product_by_scalar")
+  {
+    TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+    TinyVector<6> coef2(2, 4, 6, 8, 10, 12);
+    PolynomialP<2, 2> P(coef);
+    PolynomialP<2, 2> Q(coef2);
+    REQUIRE(Q == (P * 2));
+    REQUIRE(Q == (2 * P));
+  }
 
   // SECTION("product")
   // {
-- 
GitLab