construction of Lagrange polynomial

src/analysis/Polynomial.hpp

@@ -122,6 +122,24 @@ class Polynomial
     return P;
   }
 
+  template <size_t M>
+  PUGS_INLINE constexpr Polynomial<N>&
+  operator*=(const Polynomial<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 polynomial");
+    }
+    Polynomial<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 Polynomial<N> operator*(const double& lambda) const
   {
@@ -306,6 +324,24 @@ class Polynomial
     return os;
   }
 
+  PUGS_INLINE
+  constexpr friend void
+  lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k, Polynomial<N>& lk)
+  {
+    for (size_t i = 0; i < N; ++i) {
+      Assert(zeros[i] < zeros[i + 1], "Interpolation values must be strictly increasing in Lagrange polynomials");
+    }
+    lk.coefficients()    = zero;
+    lk.coefficients()[0] = 1;
+    for (size_t i = 0; i < N + 1; ++i) {
+      if (k == i)
+        continue;
+      double factor = 1. / (zeros[k] - zeros[i]);
+      Polynomial<1> P({-zeros[i] * factor, factor});
+      lk *= P;
+    }
+  }
+
   PUGS_INLINE constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {}
 
   PUGS_INLINE
@@ -315,5 +351,4 @@ class Polynomial
   constexpr Polynomial() noexcept = default;
   ~Polynomial()                   = default;
 };
-
 #endif   // POLYNOMIAL_HPP

tests/test_Polynomial.cpp

@@ -64,7 +64,14 @@ TEST_CASE("Polynomial", "[analysis]")
   {
     Polynomial<2> P({2, 3, 4});
     Polynomial<3> Q({1, 2, -1, 1});
+    Polynomial<4> R;
+    Polynomial<5> S;
+    R = P;
+    S = P;
+    S *= Q;
     REQUIRE(Polynomial<5>({2, 7, 8, 7, -1, 4}) == (P * Q));
+    REQUIRE(Polynomial<5>({2, 7, 8, 7, -1, 4}) == S);
+    // REQUIRE_THROWS_AS(R *= Q, AssertError);
   }
   SECTION("evaluation")
   {
@@ -140,4 +147,16 @@ TEST_CASE("Polynomial", "[analysis]")
     REQUIRE(P.compose(Q) == Polynomial<2>({2, -6, 12}));
     REQUIRE(R(S) == Polynomial<4>({4, 10, 22, 24, 12}));
   }
+
+  SECTION("Lagrange polynomial")
+  {
+    Polynomial<1> S({0.5, -0.5});
+    Polynomial<1> Q;
+    lagrangePolynomial({-1, 1}, 0, Q);
+    REQUIRE(S == Q);
+    Polynomial<2> P({0, -0.5, 0.5});
+    Polynomial<2> R;
+    lagrangePolynomial({-1, 0, 1}, 0, R);
+    REQUIRE(R == P);
+  }
 }