Fix Euclidean division

src/analysis/Polynomial.hpp

@@ -146,7 +146,8 @@ class Polynomial
   }
 
   template <size_t M>
-  PUGS_INLINE constexpr Polynomial<M + N> operator*(const Polynomial<M>& Q) const
+  PUGS_INLINE constexpr Polynomial<M + N>
+  operator*(const Polynomial<M>& Q) const
   {
     Polynomial<M + N> P;
     P.coefficients() = zero;
@@ -177,7 +178,8 @@ class Polynomial
   }
 
   PUGS_INLINE
-  constexpr Polynomial<N> operator*(const double& lambda) const
+  constexpr Polynomial<N>
+  operator*(const double& lambda) const
   {
     TinyVector<N + 1> mult_coefs = lambda * m_coefficients;
     Polynomial<N> M(mult_coefs);
@@ -266,7 +268,8 @@ class Polynomial
   }
 
   PUGS_INLINE
-  constexpr friend Polynomial<N> operator*(const double& lambda, const Polynomial<N> P)
+  constexpr friend Polynomial<N>
+  operator*(const double& lambda, const Polynomial<N> P)
   {
     return P * lambda;
   }
@@ -305,13 +308,13 @@ class Polynomial
     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 (ssize_t k = Nr - Mr; k >= 0; --k) {
+      Q.coefficients()[k] = R.coefficients()[Mr + k] / P2.coefficients()[Mr];
+      for (ssize_t j = Mr + k; j >= k; --j) {
+        R.coefficients()[j] -= Q.coefficients()[k] * P2.coefficients()[j - k];
       }
     }
-    for (size_t j = Mr; j < Nr + 1; ++j) {
+    for (size_t j = Mr; j <= Nr; ++j) {
       R.coefficients()[j] = 0;
     }
   }

tests/test_Polynomial.cpp

@@ -102,6 +102,20 @@ TEST_CASE("Polynomial", "[analysis]")
     REQUIRE(Polynomial<2>({1, 0, 0}) == S);
     REQUIRE(Polynomial<2>({0, 1, 0}) == R);
   }
+  SECTION("divide2")
+  {
+    Polynomial<7> P({1, -3, 0.6, -2.6, 7, 1.8, 4, -2});
+    Polynomial<3> Q({-1, 1, 0, 2});
+
+    Polynomial<7> R;
+    Polynomial<7> S;
+    REQUIRE(P.realDegree() == 7);
+    REQUIRE(Q.realDegree() == 3);
+
+    divide(P, Q, R, S);
+    REQUIRE(Polynomial<7>({-1, 2, 1.4, 2, -1, 0, 0, 0}) == R);
+    REQUIRE(Polynomial<7>({0, 0, 0, 0, 0, 0, 0, 0}) == S);
+  }
   SECTION("evaluation")
   {
     Polynomial<2> P({2, -3, 4});