From 2aa7c90c490d406bc6bf3a194514c2f0d5b88382 Mon Sep 17 00:00:00 2001
From: LABOURASSE Emmanuel <labourassee@gmail.com>
Date: Fri, 21 Aug 2020 19:00:45 +0200
Subject: [PATCH] divide polynomials and attempt to build Lagrange basis
---
src/analysis/Polynomial.hpp | 77 +++++++++++++++++++++++++++++++++++++
tests/test_Polynomial.cpp | 36 +++++++++++++++++
2 files changed, 113 insertions(+)
diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp
index a2e82f9c1..c768d2d67 100644
--- a/src/analysis/Polynomial.hpp
+++ b/src/analysis/Polynomial.hpp
@@ -12,6 +12,23 @@ class Polynomial
static_assert((N >= 0), "Polynomial degree must be non-negative");
public:
+ PUGS_INLINE
+ constexpr size_t
+ degree() const
+ {
+ return N;
+ }
+ PUGS_INLINE
+ constexpr size_t
+ realDegree() const
+ {
+ for (size_t j = N; j > 0; j--) {
+ if (std::abs(coefficients()[j]) > 1.e-14) {
+ return j;
+ }
+ }
+ return 0;
+ }
PUGS_INLINE
constexpr const TinyVector<N + 1>&
coefficients() const
@@ -64,6 +81,14 @@ class Polynomial
return P;
}
+ PUGS_INLINE constexpr Polynomial<N>
+ operator-() const
+ {
+ Polynomial<N> P;
+ P.coefficients() = -coefficients();
+ return P;
+ }
+
template <size_t M>
PUGS_INLINE constexpr Polynomial<std::max(M, N)>
operator-(const Polynomial<M>& Q) const
@@ -109,6 +134,17 @@ class Polynomial
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
{
@@ -230,6 +266,25 @@ class Polynomial
return bcoef;
}
+ template <size_t M>
+ PUGS_INLINE constexpr friend void
+ divide(const Polynomial<N>& P1, const Polynomial<M>& P2, Polynomial<N>& Q, Polynomial<N>& R)
+ {
+ const size_t Nr = P1.realDegree();
+ const size_t Mr = P2.realDegree();
+ R.coefficients() = P1.coefficients();
+ Q.coefficients() = zero;
+ for (size_t k = Nr - Mr + 1; k > 0; --k) {
+ Q.coefficients()[k - 1] = R.coefficients()[Mr + k - 1] / P2.coefficients()[Mr];
+ for (size_t j = Mr + k - 1; j > (k - 1); --j) {
+ R.coefficients()[j] -= Q.coefficients()[k - 1] * P2.coefficients()[j - k];
+ }
+ }
+ for (size_t j = Mr; j < Nr + 1; ++j) {
+ R.coefficients()[j] = 0;
+ }
+ }
+
PUGS_INLINE
constexpr friend Polynomial<N + 1>
primitive(const Polynomial<N>& P)
@@ -342,6 +397,28 @@ class Polynomial
}
}
+ PUGS_INLINE
+ constexpr friend void
+ lagrangeBasis(const TinyVector<N + 1> zeros, TinyVector<N + 1, Polynomial<N>>& basis)
+ {
+ Polynomial<N> lj;
+ for (size_t j = 0; j < N + 1; ++j) {
+ lagrangePolynomial(zeros, j, basis[j]);
+ }
+ }
+
+ PUGS_INLINE
+ constexpr Polynomial<N>
+ lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, TinyVector<N + 1, Polynomial<N>>& basis)
+ {
+ Polynomial<N> P(zero);
+ // lagrangeBasis({0, 0, 0}, basis);
+ for (size_t j = 0; j < N + 1; ++j) {
+ P += basis[j] * lagrange_coefs[j];
+ }
+ return P;
+ }
+
PUGS_INLINE constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {}
PUGS_INLINE
diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp
index a935d3992..7c8951c90 100644
--- a/tests/test_Polynomial.cpp
+++ b/tests/test_Polynomial.cpp
@@ -24,6 +24,11 @@ TEST_CASE("Polynomial", "[analysis]")
{
REQUIRE_NOTHROW(Polynomial<2>{{2, 3, 4}});
}
+ SECTION("degree")
+ {
+ Polynomial<2> P({2, 3, 4});
+ REQUIRE(P.degree() == 2);
+ }
SECTION("equality")
{
Polynomial<2> P({2, 3, 4});
@@ -43,6 +48,12 @@ TEST_CASE("Polynomial", "[analysis]")
REQUIRE(S == (P + Q));
REQUIRE((T + P) == U);
}
+ SECTION("opposed")
+ {
+ Polynomial<2> P({2, 3, 4});
+ Polynomial<2> Q = -P;
+ REQUIRE(Q == Polynomial<2>({-2, -3, -4}));
+ }
SECTION("difference")
{
Polynomial<2> P({2, 3, 4});
@@ -52,6 +63,8 @@ TEST_CASE("Polynomial", "[analysis]")
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("product_by_scalar")
{
@@ -73,6 +86,22 @@ TEST_CASE("Polynomial", "[analysis]")
REQUIRE(Polynomial<5>({2, 7, 8, 7, -1, 4}) == S);
// REQUIRE_THROWS_AS(R *= Q, AssertError);
}
+ SECTION("divide")
+ {
+ Polynomial<2> P({1, 0, 1});
+ Polynomial<1> Q({0, 1});
+ Polynomial<1> Q1({0, 1});
+
+ Polynomial<2> R;
+ Polynomial<2> S;
+ REQUIRE(P.realDegree() == 2);
+ REQUIRE(Q.realDegree() == 1);
+ REQUIRE(Q1.realDegree() == 1);
+
+ divide(P, Q1, R, S);
+ REQUIRE(Polynomial<2>({1, 0, 0}) == S);
+ REQUIRE(Polynomial<2>({0, 1, 0}) == R);
+ }
SECTION("evaluation")
{
Polynomial<2> P({2, -3, 4});
@@ -158,5 +187,12 @@ TEST_CASE("Polynomial", "[analysis]")
Polynomial<2> R;
lagrangePolynomial({-1, 0, 1}, 0, R);
REQUIRE(R == P);
+ TinyVector<2, Polynomial<2>> basis;
+ // lagrangeBasis<2>({-1,0,1},basis);
+ basis[0] = Polynomial<2>({0, -0.5, 0.5});
+ basis[1] = Polynomial<2>({1, 0, -1});
+ basis[2] = Polynomial<2>({0, 0.5, 0.5});
+ // Q = lagrangeToCanonical({1,0,1},basis);
+ // REQUIRE(lagrangeToCanonical<2>({1, 0, 1}, basis) == Polynomial<2>({0, 0, 1}));
}
}
--
GitLab