diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp
index 9014679bfbe861078cb35802a839f5cf1f2972e9..cf716a246a4bfa81d2de207c0101f23ea4aa7c09 100644
--- a/src/analysis/Polynomial.hpp
+++ b/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)
@@ -391,7 +391,7 @@ class Polynomial
   }
 
   PUGS_INLINE constexpr Polynomial& operator=(const Polynomial&) = default;
-  PUGS_INLINE constexpr Polynomial& operator=(Polynomial&&) = default;
+  PUGS_INLINE constexpr Polynomial& operator=(Polynomial&&)      = default;
 
   PUGS_INLINE
   constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients(coefficients) {}
diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp
index 2b44b74ded4e90b0aa88f88404ddc6bd9010de31..4089bf27d42b6c704e1f707ccb01387a79121565 100644
--- a/tests/test_Polynomial.cpp
+++ b/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")