implement a polynomial class with elementary opertaions

src/analysis/Polynomial.hpp

@@ -9,6 +9,7 @@ class Polynomial
  private:
   using Coefficients = TinyVector<N + 1, double>;
   Coefficients m_coefficients;
+  static_assert((N >= 0), "Polynomial degree must be non-negative");
 
  public:
   PUGS_INLINE
@@ -25,10 +26,12 @@ class Polynomial
     return m_coefficients;
   }
 
-  PUGS_INLINE
-  constexpr bool
-  operator==(const Polynomial<N>& q) const
+  template <size_t M>
+  PUGS_INLINE constexpr bool
+  operator==(const Polynomial<M>& q) const
   {
+    if (M != N)
+      return false;
     if (m_coefficients != q.coefficients()) {
       return false;
     }
@@ -42,22 +45,55 @@ class Polynomial
     return not this->operator==(q);
   }
 
-  PUGS_INLINE
-  constexpr Polynomial<N>
-  operator+(const Polynomial<N>& q) const
+  template <size_t M>
+  PUGS_INLINE constexpr Polynomial<std::max(M, N)>
+  operator+(const Polynomial<M>& Q) const
   {
-    TinyVector<N + 1> sum_coefs = m_coefficients + q.coefficients();
-    Polynomial<N> S(sum_coefs);
-    return S;
+    Polynomial<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;
   }
 
-  PUGS_INLINE
-  constexpr Polynomial<N>
-  operator-(const Polynomial<N>& q) const
+  template <size_t M>
+  PUGS_INLINE constexpr Polynomial<std::max(M, N)>
+  operator-(const Polynomial<M>& Q) const
+  {
+    Polynomial<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 Polynomial<M + N> operator*(const Polynomial<M>& Q) const
   {
-    TinyVector<N + 1> diff_coefs = m_coefficients - q.coefficients();
-    Polynomial<N> D(diff_coefs);
-    return D;
+    Polynomial<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;
   }
 
   PUGS_INLINE
@@ -112,47 +148,89 @@ class Polynomial
     return Polynomial<N + 1>{coefs};
   }
 
+  PUGS_INLINE
+  constexpr friend double
+  integrate(const Polynomial<N>& P, const double& xinf, const double& xsup)
+  {
+    Polynomial<N + 1> Q = primitive(P);
+    return (Q(xsup) - Q(xinf));
+  }
+
+  PUGS_INLINE
+  constexpr friend auto
+  derivative(const Polynomial<N>& P)
+  {
+    if constexpr (N == 0) {
+      return Polynomial<0>(0);
+    } else {
+      TinyVector<N> coefs;
+      for (size_t i = 0; i < N; ++i) {
+        coefs[i] = double(i + 1) * P.coefficients()[i + 1];
+      }
+      return Polynomial<N - 1>(coefs);
+    }
+  }
+
   PUGS_INLINE
   constexpr friend std::ostream&
   operator<<(std::ostream& os, const Polynomial<N>& P)
   {
     //    os << "P(x) = ";
     bool all_coef_zero = true;
-    for (size_t i = N; i > 1; --i) {
+    if (N == 0) {
+      os << P.coefficients()[0];
+      return os;
+    }
+    if (N != 1) {
+      if (P.coefficients()[N] != 0.) {
+        if (P.coefficients()[N] < 0.) {
+          os << "- ";
+        }
+        if (P.coefficients()[N] != 1 && P.coefficients()[N] != -1) {
+          os << std::abs(P.coefficients()[N]);
+        }
+        os << "x^" << N;
+        all_coef_zero = false;
+      }
+    }
+    for (size_t i = N - 1; i > 1; --i) {
       if (P.coefficients()[i] != 0.) {
+        if (P.coefficients()[i] > 0.) {
+          os << " + ";
+        } else if (P.coefficients()[i] < 0.) {
+          os << " - ";
+        }
         if (P.coefficients()[i] != 1 && P.coefficients()[i] != -1) {
           os << std::abs(P.coefficients()[i]);
         }
         os << "x^" << i;
-        if (P.coefficients()[i - 1] > 0.)
-          os << " + ";
-        else
-          os << " - ";
         all_coef_zero = false;
-      } else {
-        os << " ";
       }
     }
     if (P.coefficients()[1] != 0.) {
+      if (P.coefficients()[1] > 0. && N != 1) {
+        os << " + ";
+      } else if (P.coefficients()[1] < 0.) {
+        os << " - ";
+      }
       if (P.coefficients()[1] != 1 && P.coefficients()[1] != -1) {
         os << std::abs(P.coefficients()[1]);
       }
       os << "x";
-      if (P.coefficients()[0] > 0.)
+      all_coef_zero = false;
+    }
+    if (P.coefficients()[0] != 0. || all_coef_zero) {
+      if (P.coefficients()[0] > 0.) {
         os << " + ";
-      else if (P.coefficients()[0] < 0.)
+      } else if (P.coefficients()[0] < 0.) {
         os << " - ";
-      all_coef_zero = false;
-    } else {
-      os << " ";
       }
-    if (P.coefficients()[0] != 0. || all_coef_zero)
       os << std::abs(P.coefficients()[0]);
+    }
     return os;
   }
 
-  PUGS_INLINE
-  constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {}
+  PUGS_INLINE constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {}
 
   PUGS_INLINE
   constexpr Polynomial(TinyVector<N + 1>&& coefficients) noexcept : m_coefficients{coefficients} {}

tests/test_Polynomial.cpp

@@ -9,10 +9,12 @@
 #include <analysis/Polynomial.hpp>
 
 // Instantiate to ensure full coverage is performed
+template class Polynomial<0>;
 template class Polynomial<1>;
 template class Polynomial<2>;
 template class Polynomial<3>;
 template class Polynomial<4>;
+template class Polynomial<5>;
 
 // clazy:excludeall=non-pod-global-static
 
@@ -26,15 +28,20 @@ TEST_CASE("Polynomial", "[analysis]")
   {
     Polynomial<2> P({2, 3, 4});
     Polynomial<2> Q({2, 3, 4});
+    Polynomial<2> S({2, 3, 5});
 
     REQUIRE(P == Q);
+    REQUIRE(P != S);
   }
   SECTION("addition")
   {
     Polynomial<2> P({2, 3, 4});
     Polynomial<2> Q({-1, -3, 2});
     Polynomial<2> S({1, 0, 6});
+    Polynomial<3> T({0, 3, 1, -2});
+    Polynomial<3> U({2, 6, 5, -2});
     REQUIRE(S == (P + Q));
+    REQUIRE((T + P) == U);
   }
   SECTION("difference")
   {
@@ -42,6 +49,9 @@ TEST_CASE("Polynomial", "[analysis]")
     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}));
   }
   SECTION("product_by_scalar")
   {
@@ -50,6 +60,12 @@ TEST_CASE("Polynomial", "[analysis]")
     REQUIRE(M == (P * 3));
     REQUIRE(M == (3 * P));
   }
+  SECTION("product")
+  {
+    Polynomial<2> P({2, 3, 4});
+    Polynomial<3> Q({1, 2, -1, 1});
+    REQUIRE(Polynomial<5>({2, 7, 8, 7, -1, 4}) == (P * Q));
+  }
   SECTION("evaluation")
   {
     Polynomial<2> P({2, -3, 4});
@@ -67,4 +83,23 @@ TEST_CASE("Polynomial", "[analysis]")
     Polynomial<3> R({0, 2, -3. / 2, 4. / 3});
     REQUIRE(Q == R);
   }
+  SECTION("integrate")
+  {
+    Polynomial<2> P({2, -3, 3});
+    double xinf   = -1;
+    double xsup   = 1;
+    double result = integrate(P, xinf, xsup);
+    REQUIRE(result == 6);
+  }
+  SECTION("derivative")
+  {
+    Polynomial<2> P({2, -3, 3});
+    Polynomial<1> Q = derivative(P);
+    REQUIRE(Q == Polynomial<1>({-3, 6}));
+
+    Polynomial<0> P2(3);
+
+    Polynomial<0> R(0);
+    REQUIRE(derivative(P2) == R);
+  }
 }