divide polynomials and attempt to build Lagrange basis

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

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}));
   }
 }