Merge remote-tracking branch 'origin/feature/polynomials' into hyperplastic

src/analysis/Polynomial.hpp

@@ -288,15 +288,15 @@ class Polynomial
     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];
+      Polynomial<N - M> Q1;
+      Q1.coefficients()        = zero;
+      Q1.coefficients()[k - 1] = Q.coefficients()[k - 1];
+      R -= Q1 * P2;
     }
-    }
-    for (size_t j = Mr; j < Nr + 1; ++j) {
+    for (size_t j = Mr + 1; j < Nr + 1; ++j) {
       R.coefficients()[j] = 0;
     }
   }
-
   PUGS_INLINE
   constexpr friend Polynomial<N + 1>
   primitive(const Polynomial& P)

tests/test_Polynomial.cpp

@@ -113,6 +113,15 @@ TEST_CASE("Polynomial", "[analysis]")
     divide(P, Q1, R, S);
     REQUIRE(Polynomial<2>{1, 0, 0} == S);
     REQUIRE(Polynomial<2>{0, 1, 0} == R);
+    Polynomial<3> P2(5, 1, 3, 2);
+    Polynomial<2> Q2(1, 0, 1);
+
+    Polynomial<3> R2;
+    Polynomial<3> S2;
+
+    divide(P2, Q2, R2, S2);
+    REQUIRE(Polynomial<3>{3, 2, 0, 0} == R2);
+    REQUIRE(Polynomial<3>{2, -1, 0, 0} == S2);
   }
 
   SECTION("evaluation")