Some more feature for polynomials (power,compostition)

src/analysis/Polynomial.hpp

@@ -83,6 +83,32 @@ class Polynomial
     return P;
   }
 
+  template <size_t M>
+  PUGS_INLINE constexpr Polynomial<N>&
+  operator=(const Polynomial<M>& Q)
+  {
+    coefficients() = zero;
+    for (size_t i = N + 1; i <= M; ++i) {
+      Assert(Q.coefficients()[i] == 0, "degree of polynomial to small in assignation");
+    }
+    //    static_assert(N >= M, "degree of polynomial 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 Polynomial<N>&
+  operator+=(const Polynomial<M>& Q)
+  {
+    static_assert(N >= M, "Polynomial 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 Polynomial<M + N> operator*(const Polynomial<M>& Q) const
   {
@@ -104,6 +130,56 @@ class Polynomial
     return M;
   }
 
+  template <size_t M>
+  PUGS_INLINE constexpr Polynomial<M * N>
+  compose(const Polynomial<M>& Q) const
+  {
+    Polynomial<M * N> P;
+    P.coefficients() = zero;
+    Polynomial<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>
+  PUGS_INLINE constexpr Polynomial<M * N>
+  operator()(const Polynomial<M>& Q) const
+  {
+    Polynomial<M * N> P;
+    P.coefficients() = zero;
+    Polynomial<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>
+  PUGS_INLINE constexpr Polynomial<M * N>
+  pow(size_t power) const
+  {
+    Assert(power <= M, "You declared a polynomial of degree too small for return of the pow function");
+    Polynomial<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 friend Polynomial<N> operator*(const double& lambda, const Polynomial<N> P)
   {

tests/CMakeLists.txt

@@ -68,6 +68,7 @@ add_executable (unit_tests
   test_NameProcessor.cpp
   test_OStreamProcessor.cpp
   test_PCG.cpp
+  test_Polynomial.cpp
   test_PugsFunctionAdapter.cpp
   test_PugsAssert.cpp
   test_RevisionInfo.cpp

tests/test_Polynomial.cpp

@@ -102,4 +102,42 @@ TEST_CASE("Polynomial", "[analysis]")
     Polynomial<0> R(0);
     REQUIRE(derivative(P2) == R);
   }
+
+  SECTION("affectation")
+  {
+    Polynomial<2> Q({2, -3, 3});
+    Polynomial<4> R({2, -3, 3, 0, 0});
+    Polynomial<4> P({0, 1, 2, 3, 3});
+    P = Q;
+    REQUIRE(P == R);
+  }
+
+  SECTION("affectation addition")
+  {
+    Polynomial<2> Q({2, -3, 3});
+    Polynomial<4> R({2, -2, 5, 3, 3});
+    Polynomial<4> P({0, 1, 2, 3, 3});
+    P += Q;
+    REQUIRE(P == R);
+  }
+
+  SECTION("power")
+  {
+    Polynomial<2> P({2, -3, 3});
+    Polynomial<4> R({4, -12, 21, -18, 9});
+    Polynomial<1> Q({0, 2});
+    Polynomial<2> S = Q.pow<2>(2);
+    REQUIRE(P.pow<2>(2) == R);
+    REQUIRE(S == Polynomial<2>({0, 0, 4}));
+  }
+
+  SECTION("composition")
+  {
+    Polynomial<2> P({2, -3, 3});
+    Polynomial<1> Q({0, 2});
+    Polynomial<2> R({2, -1, 3});
+    Polynomial<2> S({1, 2, 2});
+    REQUIRE(P.compose(Q) == Polynomial<2>({2, -6, 12}));
+    REQUIRE(R(S) == Polynomial<4>({4, 10, 22, 24, 12}));
+  }
 }