diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp index f9929b0bd79eb85d75ea412b8f9a53b7ae1644ce..a2e82f9c1bece29c400a101cdd6a33ea8d109a9d 100644 --- a/src/analysis/Polynomial.hpp +++ b/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 diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp index 1e5d6f83fbb36416f09dd0dbe6e31d5c66f3e3dc..a935d3992329836642bfdf0a2044c75638aaa153 100644 --- a/tests/test_Polynomial.cpp +++ b/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); + } }