diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt
index 1f32dd6f..ecbb97c9 100644
--- a/include/CMakeLists.txt
+++ b/include/CMakeLists.txt
@@ -1,10 +1,10 @@
 # Configure the version
 configure_file(version.h.in ${CMAKE_CURRENT_SOURCE_DIR}/version.h)
 
-# Find all the base headers and add them to install terget
+# Find all the base headers and add them to install target
 file(GLOB headers *.h)
 install(FILES ${headers} DESTINATION include/poly)
 
-# Find all the cxx headers and add them to install terget
+# Find all the cxx headers and add them to install target
 file(GLOB headers_cxx polyxx/*.h)
 install(FILES ${headers_cxx} DESTINATION include/poly/polyxx)
\ No newline at end of file
diff --git a/include/polynomial.h b/include/polynomial.h
index 8a9533fc..48dfa20d 100644
--- a/include/polynomial.h
+++ b/include/polynomial.h
@@ -88,9 +88,6 @@ void lp_polynomial_swap(lp_polynomial_t* A1, lp_polynomial_t* A2);
 /** Assign the polynomial a given polynomial. */
 void lp_polynomial_assign(lp_polynomial_t* A, const lp_polynomial_t* from);
 
-/** Returns the context of the polynomial */
-const lp_polynomial_context_t* lp_polynomial_context(const lp_polynomial_t* A);
-
 /** Returns the degree of the polynomial (in it's top variable) */
 size_t lp_polynomial_degree(const lp_polynomial_t* A);
 
@@ -214,14 +211,14 @@ void lp_polynomial_add_mul(lp_polynomial_t* S, const lp_polynomial_t* A1, const
 void lp_polynomial_sub_mul(lp_polynomial_t* S, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
 
 /**
- * Reduce the polynomial A in Z[y,x] using B in Z[y,x] so that
+ * Reduce the polynomial A in Z[x,y] using B in Z[x,y] so that
  *
  *   P*A = Q*B + R
  *
  * and
  *
- *   P in Z[y]
- *   Q, R in Z[y,x]
+ *   P in Z[x]
+ *   Q, R in Z[x,y]
  *
  * with
  *
@@ -239,9 +236,15 @@ void lp_polynomial_rem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_p
 /** Compute a*A1 = D*A2 + R (pseudo remainder). */
 void lp_polynomial_prem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
 
+/** Compute a*A1 = D*A2 + R (pseudo remainder). */
+void lp_polynomial_pdivrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
+
 /** Compute a*A1 = D*A2 + R (sparse pseudo remainder). */
 void lp_polynomial_sprem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
 
+/** Compute a*A1 = D*A2 + R (sparse pseudo remainder). */
+void lp_polynomial_spdivrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
+
 /** Compute A1 = D*A2 + R (assumes that exact division). */
 void lp_polynomial_divrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
 
diff --git a/python/polypyPolynomial2.c b/python/polypyPolynomial2.c
index 6903da2b..da42fa68 100644
--- a/python/polypyPolynomial2.c
+++ b/python/polypyPolynomial2.c
@@ -84,6 +84,12 @@ Polynomial_prem(PyObject* self, PyObject* args);
 static PyObject*
 Polynomial_sprem(PyObject* self, PyObject* args);
 
+static PyObject*
+Polynomial_pdivrem(PyObject* self, PyObject* args);
+
+static PyObject*
+Polynomial_spdivrem(PyObject* self, PyObject* args);
+
 static PyObject*
 Polynomial_gcd(PyObject* self, PyObject* args);
 
@@ -183,6 +189,8 @@ PyMethodDef Polynomial_methods[] = {
     {"rem", (PyCFunction)Polynomial_rem, METH_VARARGS, "Returns the remainder of current and given polynomial"},
     {"prem", (PyCFunction)Polynomial_prem, METH_VARARGS, "Returns the pseudo remainder of current and given polynomial"},
     {"sprem", (PyCFunction)Polynomial_sprem, METH_VARARGS, "Returns the sparse pseudo remainder of current and given polynomial"},
+    {"pdivrem", (PyCFunction)Polynomial_pdivrem, METH_VARARGS, "Returns the pseudo quotient and pseudo remainder of current and given polynomial"},
+    {"spdivprem", (PyCFunction)Polynomial_spdivrem, METH_VARARGS, "Returns the sparse pseudo quotient and sparse pseudo remainder of current and given polynomial"},
     {"gcd", (PyCFunction)Polynomial_gcd, METH_VARARGS, "Returns the gcd of current and given polynomial"},
     {"lcm", (PyCFunction)Polynomial_lcm, METH_VARARGS, "Returns the lcm of current and given polynomial"},
     {"extended_gcd", (PyCFunction)Polynomial_extended_gcd, METH_VARARGS, "Returns the extended gcd, i.e. (gcd, u, v), of current and given polynomial"},
@@ -749,8 +757,14 @@ Polynomial_rem_operator(PyObject* self, PyObject* other) {
   return Polynomial_create(rem);
 }
 
+enum rem_type {
+  REM_EXACT,
+  REM_PSEUDO,
+  REM_SPARSE_PSEUDO
+};
+
 static PyObject*
-Polynomial_divmod(PyObject* self, PyObject* other) {
+Polynomial_divmod_general(PyObject* self, PyObject* args, enum rem_type type) {
 
   int dec_other = 0;
 
@@ -763,7 +777,18 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
   Polynomial* p1 = (Polynomial*) self;
   const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
 
-  // Make sure other is a polynomial
+  // check that there is only one other
+  PyObject *other;
+  if (PyTuple_Check(args)) {
+    if (PyTuple_Size(args) != 1) {
+      Py_INCREF(Py_NotImplemented);
+      return Py_NotImplemented;
+    }
+    other = PyTuple_GetItem(args, 0);
+  } else {
+    other = args;
+  }
+  // make sure other is a polynomial
   if (!PyPolynomial_CHECK(other)) {
     if (PyVariable_CHECK(other)) {
       other = PyPolynomial_FromVariable(other, p1_ctx);
@@ -785,10 +810,20 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
     return Py_NotImplemented;
   }
 
-  // Multiply the polynomials
+  // Divide the polynomials
   lp_polynomial_t* rem = lp_polynomial_new(p1_ctx);
   lp_polynomial_t* div = lp_polynomial_new(p1_ctx);
-  lp_polynomial_divrem(div, rem, p1->p, p2->p);
+  switch (type) {
+    case REM_EXACT:
+      lp_polynomial_divrem(div, rem, p1->p, p2->p);
+      break;
+    case REM_PSEUDO:
+      lp_polynomial_pdivrem(div, rem, p1->p, p2->p);
+      break;
+    case REM_SPARSE_PSEUDO:
+      lp_polynomial_spdivrem(div, rem, p1->p, p2->p);
+      break;
+  }
 
   if (dec_other) {
     Py_DECREF(other);
@@ -805,6 +840,21 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
   return pair;
 }
 
+static PyObject*
+Polynomial_divmod(PyObject* self, PyObject* other) {
+  return Polynomial_divmod_general(self, other, REM_EXACT);
+}
+
+static PyObject*
+Polynomial_pdivrem(PyObject* self, PyObject* other) {
+  return Polynomial_divmod_general(self, other, REM_PSEUDO);
+}
+
+static PyObject*
+Polynomial_spdivrem(PyObject* self, PyObject* other) {
+  return Polynomial_divmod_general(self, other, REM_SPARSE_PSEUDO);
+}
+
 static int
 Polynomial_nonzero(PyObject* self) {
   // Get arguments
@@ -813,12 +863,6 @@ Polynomial_nonzero(PyObject* self) {
   return !lp_polynomial_is_zero(p->p);
 }
 
-enum rem_type {
-  REM_EXACT,
-  REM_PSEUDO,
-  REM_SPARSE_PSEUDO
-};
-
 static PyObject*
 Polynomial_rem_general(PyObject* self, PyObject* args, enum rem_type type) {
   int dec_other = 0;
diff --git a/python/polypyPolynomial3.c b/python/polypyPolynomial3.c
index 33c10b86..fdad3bae 100644
--- a/python/polypyPolynomial3.c
+++ b/python/polypyPolynomial3.c
@@ -81,6 +81,12 @@ Polynomial_prem(PyObject* self, PyObject* args);
 static PyObject*
 Polynomial_sprem(PyObject* self, PyObject* args);
 
+static PyObject*
+Polynomial_pdivrem(PyObject* self, PyObject* args);
+
+static PyObject*
+Polynomial_spdivrem(PyObject* self, PyObject* args);
+
 static PyObject*
 Polynomial_gcd(PyObject* self, PyObject* args);
 
@@ -180,6 +186,8 @@ PyMethodDef Polynomial_methods[] = {
     {"rem", (PyCFunction)Polynomial_rem, METH_VARARGS, "Returns the remainder of current and given polynomial"},
     {"prem", (PyCFunction)Polynomial_prem, METH_VARARGS, "Returns the pseudo remainder of current and given polynomial"},
     {"sprem", (PyCFunction)Polynomial_sprem, METH_VARARGS, "Returns the sparse pseudo remainder of current and given polynomial"},
+    {"pdivrem", (PyCFunction)Polynomial_pdivrem, METH_VARARGS, "Returns the pseudo quotient and pseudo remainder of current and given polynomial"},
+    {"spdivprem", (PyCFunction)Polynomial_spdivrem, METH_VARARGS, "Returns the sparse pseudo quotient and sparse pseudo remainder of current and given polynomial"},
     {"gcd", (PyCFunction)Polynomial_gcd, METH_VARARGS, "Returns the gcd of current and given polynomial"},
     {"lcm", (PyCFunction)Polynomial_lcm, METH_VARARGS, "Returns the lcm of current and given polynomial"},
     {"extended_gcd", (PyCFunction)Polynomial_extended_gcd, METH_VARARGS, "Returns the extended gcd, i.e. (gcd, u, v), of current and given polynomial"},
@@ -730,8 +738,14 @@ Polynomial_rem_operator(PyObject* self, PyObject* other) {
   return Polynomial_create(rem);
 }
 
+enum rem_type {
+  REM_EXACT,
+  REM_PSEUDO,
+  REM_SPARSE_PSEUDO
+};
+
 static PyObject*
-Polynomial_divmod(PyObject* self, PyObject* other) {
+Polynomial_divmod_general(PyObject* self, PyObject* args, enum rem_type type) {
 
   int dec_other = 0;
 
@@ -741,10 +755,21 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
   }
 
   // self is always a polynomial
-  Polynomial* p1 = (Polynomial*) self;
-  const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
+  Polynomial *p1 = (Polynomial *) self;
+  const lp_polynomial_context_t *p1_ctx = lp_polynomial_get_context(p1->p);
 
-  // Make sure other is a polynomial
+  // check that there is only one other
+  PyObject *other;
+  if (PyTuple_Check(args)) {
+    if (PyTuple_Size(args) != 1) {
+      Py_INCREF(Py_NotImplemented);
+      return Py_NotImplemented;
+    }
+    other = PyTuple_GetItem(args, 0);
+  } else {
+    other = args;
+  }
+  // make sure other is a polynomial
   if (!PyPolynomial_CHECK(other)) {
     if (PyVariable_CHECK(other)) {
       other = PyPolynomial_FromVariable(other, p1_ctx);
@@ -759,26 +784,36 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
   }
 
   // other can be a variable or a number
-  Polynomial* p2 = (Polynomial*) other;
-  const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
+  Polynomial *p2 = (Polynomial *) other;
+  const lp_polynomial_context_t *p2_ctx = lp_polynomial_get_context(p2->p);
   if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
     Py_INCREF(Py_NotImplemented);
     return Py_NotImplemented;
   }
 
-  // Multiply the polynomials
-  lp_polynomial_t* rem = lp_polynomial_new(p1_ctx);
-  lp_polynomial_t* div = lp_polynomial_new(p1_ctx);
-  lp_polynomial_divrem(div, rem, p1->p, p2->p);
+  // Divide the polynomials
+  lp_polynomial_t *rem = lp_polynomial_new(p1_ctx);
+  lp_polynomial_t *div = lp_polynomial_new(p1_ctx);
+  switch (type) {
+    case REM_EXACT:
+      lp_polynomial_divrem(div, rem, p1->p, p2->p);
+      break;
+    case REM_PSEUDO:
+      lp_polynomial_pdivrem(div, rem, p1->p, p2->p);
+      break;
+    case REM_SPARSE_PSEUDO:
+      lp_polynomial_spdivrem(div, rem, p1->p, p2->p);
+      break;
+  }
 
   if (dec_other) {
     Py_DECREF(other);
   }
 
   // Return the result
-  PyObject* pair = PyTuple_New(2);
-  PyObject* divObj = Polynomial_create(div);
-  PyObject* remObj = Polynomial_create(rem);
+  PyObject *pair = PyTuple_New(2);
+  PyObject *divObj = Polynomial_create(div);
+  PyObject *remObj = Polynomial_create(rem);
   Py_INCREF(divObj);
   Py_INCREF(remObj);
   PyTuple_SetItem(pair, 0, divObj);
@@ -786,6 +821,21 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
   return pair;
 }
 
+static PyObject*
+Polynomial_divmod(PyObject* self, PyObject* other) {
+  return Polynomial_divmod_general(self, other, REM_EXACT);
+}
+
+static PyObject*
+Polynomial_pdivrem(PyObject* self, PyObject* other) {
+  return Polynomial_divmod_general(self, other, REM_PSEUDO);
+}
+
+static PyObject*
+Polynomial_spdivrem(PyObject* self, PyObject* other) {
+  return Polynomial_divmod_general(self, other, REM_SPARSE_PSEUDO);
+}
+
 static int
 Polynomial_nonzero(PyObject* self) {
   // Get arguments
@@ -794,16 +844,15 @@ Polynomial_nonzero(PyObject* self) {
   return !lp_polynomial_is_zero(p->p);
 }
 
-enum rem_type {
-  REM_EXACT,
-  REM_PSEUDO,
-  REM_SPARSE_PSEUDO
-};
-
 static PyObject*
 Polynomial_rem_general(PyObject* self, PyObject* args, enum rem_type type) {
   int dec_other = 0;
 
+  if (!PyPolynomial_CHECK(self)) {
+    Py_INCREF(Py_NotImplemented);
+    return Py_NotImplemented;
+  }
+
   // self is always a polynomial
   Polynomial* p1 = (Polynomial*) self;
   const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
diff --git a/src/number/value.c b/src/number/value.c
index a964edc7..cfd4d5f0 100644
--- a/src/number/value.c
+++ b/src/number/value.c
@@ -276,7 +276,7 @@ int lp_value_sgn(const lp_value_t* v) {
 int lp_value_cmp(const lp_value_t* v1, const lp_value_t* v2) {
 
   if (trace_is_enabled("value::cmp")) {
-    tracef("lp_value_cmp()\n")
+    tracef("lp_value_cmp()\n");
     tracef("v1 = "); lp_value_print(v1, trace_out); tracef("\n");
     tracef("v2 = "); lp_value_print(v2, trace_out); tracef("\n");
   }
@@ -612,7 +612,7 @@ char* lp_value_to_string(const lp_value_t* v) {
 void lp_value_get_value_between(const lp_value_t* a, int a_strict, const lp_value_t* b, int b_strict, lp_value_t* v) {
 
   if (trace_is_enabled("value::get_value_between")) {
-    tracef("lp_value_get_value_between()\n")
+    tracef("lp_value_get_value_between()\n");
     tracef("a = "); lp_value_print(a, trace_out); tracef(", a_strict = %d\n", a_strict);
     tracef("b = "); lp_value_print(b, trace_out); tracef(", b_strict = %d\n", b_strict);
   }
diff --git a/src/polynomial/coefficient.c b/src/polynomial/coefficient.c
index c8cc6433..447140e0 100644
--- a/src/polynomial/coefficient.c
+++ b/src/polynomial/coefficient.c
@@ -404,7 +404,7 @@ void coefficient_reductum(const lp_polynomial_context_t* ctx, coefficient_t* R,
 
   assert(C->type == COEFFICIENT_POLYNOMIAL);
 
-  // Locate the first non-zero ceofficient past the top one
+  // Locate the first non-zero coefficient past the top one
   int i = SIZE(C) - 2;
   while (i >= 0 && coefficient_is_zero(ctx, COEFF(C, i))) {
     -- i;
@@ -436,7 +436,7 @@ void coefficient_reductum_m(const lp_polynomial_context_t* ctx, coefficient_t* R
 
   assert(C->type == COEFFICIENT_POLYNOMIAL);
 
-  // Locate the first non-zero ceofficient (normal reductum is the next nonzero)
+  // Locate the first non-zero coefficient (normal reductum is the next nonzero)
   int i = SIZE(C) - 1;
   while (i >= 0 && coefficient_sgn(ctx, COEFF(C, i), m) == 0) {
     if (assumptions != 0 && !coefficient_is_constant(COEFF(C, i))) {
@@ -1136,7 +1136,7 @@ int coefficient_cmp_general(const lp_polynomial_context_t* ctx, const coefficien
     }
   }
 
-  TRACE("coefficien::internal", "coefficient_cmp() => %d\n", cmp);
+  TRACE("coefficient::internal", "coefficient_cmp() => %d\n", cmp);
   return cmp;
 }
 
@@ -2109,6 +2109,8 @@ void coefficient_div(const lp_polynomial_context_t* ctx, coefficient_t* D, const
     coefficient_div_constant(ctx, D, &C2->value.num);
     return;
   }
+  // A polynomial does not divide a constant
+  assert(!coefficient_is_constant(C1));
 
   // If different variables
   if (VAR(C1) != VAR(C2)) {
@@ -2275,6 +2277,88 @@ void coefficient_prem(const lp_polynomial_context_t* ctx, coefficient_t* R, cons
   assert(coefficient_is_normalized(ctx, R));
 }
 
+STAT_DECLARE(int, coefficient, pdivrem)
+
+void coefficient_pdivrem(const lp_polynomial_context_t* ctx, coefficient_t* D, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2) {
+    TRACE("coefficient", "coefficient_pdivrem()\n");
+    STAT_INCR(coefficient, pdivrem)
+
+    if (trace_is_enabled("coefficient")) {
+        tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
+        tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
+    }
+
+    assert(!coefficient_is_zero(ctx, C2));
+
+    int cmp_type = coefficient_cmp_type(ctx, C1, C2);
+
+    assert(cmp_type >= 0);
+
+    if (cmp_type == 0 && C1->type == COEFFICIENT_NUMERIC) {
+        assert(C2->type == COEFFICIENT_NUMERIC);
+        if (R->type == COEFFICIENT_POLYNOMIAL) {
+            coefficient_destruct(R);
+            coefficient_construct(ctx, R);
+        }
+        if (D->type == COEFFICIENT_POLYNOMIAL) {
+            coefficient_destruct(D);
+            coefficient_construct(ctx, D);
+        }
+        integer_div_rem_Z(&D->value.num, &R->value.num, &C1->value.num, &C2->value.num);
+    } else {
+        coefficient_reduce(ctx, C1, C2, 0, D, R, REMAINDERING_PSEUDO_DENSE);
+    }
+
+    if (trace_is_enabled("coefficient")) {
+        tracef("coefficient_spdivrem() => \n");
+        tracef("D = "); coefficient_print(ctx, D, trace_out); tracef("\n");
+        tracef("R = "); coefficient_print(ctx, R, trace_out); tracef("\n");
+    }
+
+    assert(coefficient_is_normalized(ctx, D));
+    assert(coefficient_is_normalized(ctx, R));
+}
+
+STAT_DECLARE(int, coefficient, spdivrem)
+
+void coefficient_spdivrem(const lp_polynomial_context_t* ctx, coefficient_t* D, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2) {
+    TRACE("coefficient", "coefficient_spdivrem()\n");
+    STAT_INCR(coefficient, spdivrem)
+
+    if (trace_is_enabled("coefficient")) {
+        tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
+        tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
+    }
+
+    assert(!coefficient_is_zero(ctx, C2));
+
+    int cmp_type = coefficient_cmp_type(ctx, C1, C2);
+
+    assert(cmp_type >= 0);
+
+    if (cmp_type == 0 && C1->type == COEFFICIENT_NUMERIC) {
+        assert(C2->type == COEFFICIENT_NUMERIC);
+        if (R->type == COEFFICIENT_POLYNOMIAL) {
+            coefficient_destruct(R);
+            coefficient_construct(ctx, R);
+        }
+        if (D->type == COEFFICIENT_POLYNOMIAL) {
+            coefficient_destruct(D);
+            coefficient_construct(ctx, D);
+        }
+        integer_div_rem_Z(&D->value.num, &R->value.num, &C1->value.num, &C2->value.num);
+    } else {
+        coefficient_reduce(ctx, C1, C2, 0, D, R, REMAINDERING_PSEUDO_SPARSE);
+    }
+
+    if (trace_is_enabled("coefficient")) {
+        tracef("coefficient_spdivrem() => \n");
+        tracef("D = "); coefficient_print(ctx, D, trace_out); tracef("\n");
+        tracef("R = "); coefficient_print(ctx, R, trace_out); tracef("\n");
+    }
+
+    assert(coefficient_is_normalized(ctx, R));
+}
 
 STAT_DECLARE(int, coefficient, divrem)
 
@@ -2301,7 +2385,11 @@ void coefficient_divrem(const lp_polynomial_context_t* ctx, coefficient_t* D, co
         coefficient_destruct(R);
         coefficient_construct(ctx, R);
       }
-      integer_rem_Z(&R->value.num, &C1->value.num, &C2->value.num);
+      if (D->type == COEFFICIENT_POLYNOMIAL) {
+        coefficient_destruct(D);
+        coefficient_construct(ctx, D);
+      }
+      integer_div_rem_Z(&D->value.num, &R->value.num, &C1->value.num, &C2->value.num);
       break;
     case COEFFICIENT_POLYNOMIAL:
     {
@@ -2320,7 +2408,7 @@ void coefficient_divrem(const lp_polynomial_context_t* ctx, coefficient_t* D, co
   if (trace_is_enabled("coefficient")) {
     tracef("coefficient_divrem() => \n");
     tracef("D = "); coefficient_print(ctx, D, trace_out); tracef("\n");
-    tracef("D = "); coefficient_print(ctx, R, trace_out); tracef("\n");
+    tracef("R = "); coefficient_print(ctx, R, trace_out); tracef("\n");
   }
 
   assert(coefficient_is_normalized(ctx, R));
@@ -3217,10 +3305,10 @@ lp_value_t* coefficient_evaluate(const lp_polynomial_context_t* ctx, const coeff
 
       if (trace_is_enabled("coefficient")) {
         tracef("coefficient_evaluate(): filtered roots:\n");
-        size_t i = 0;
-        for (i = 0; i < roots_size; ++ i) {
-          tracef("%zu: ", i);
-          lp_value_print(roots + i, trace_out);
+        size_t j = 0;
+        for (j = 0; j < roots_size; ++ j) {
+          tracef("%zu: ", j);
+          lp_value_print(roots + j, trace_out);
           tracef("\n");
         }
       }
diff --git a/src/polynomial/coefficient.h b/src/polynomial/coefficient.h
index 31835b63..3f31950d 100644
--- a/src/polynomial/coefficient.h
+++ b/src/polynomial/coefficient.h
@@ -261,14 +261,14 @@ void coefficient_add_mul(const lp_polynomial_context_t* ctx, coefficient_t* S, c
 void coefficient_sub_mul(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C1, const coefficient_t* C2);
 
 /**
- * Reduce the coefficient A in Z[y,x] using B in Z[y,x] so that
+ * Reduce the coefficient A in Z[x,y] using B in Z[x,y] so that
  *
  *   P*A = Q*B + R
  *
  * and
  *
- *   P in Z[y]
- *   Q, R in Z[y,x]
+ *   P in Z[x]
+ *   Q, R in Z[x,y]
  *
  * with
  *
@@ -293,9 +293,15 @@ void coefficient_rem(const lp_polynomial_context_t* ctx, coefficient_t* R, const
 /** Compute the pseudo remainder */
 void coefficient_prem(const lp_polynomial_context_t* ctx, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2);
 
+/** Compute a*C1 = D*C2 + R, in the given ring (using pseudo division). */
+void coefficient_pdivrem(const lp_polynomial_context_t* ctx, coefficient_t* D, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2);
+
 /** Compute the sparse pseudo remainder */
 void coefficient_sprem(const lp_polynomial_context_t* ctx, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2);
 
+/** Compute a*C1 = D*C2 + R, in the given ring (using sparse pseudo division). */
+void coefficient_spdivrem(const lp_polynomial_context_t* ctx, coefficient_t* D, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2);
+
 /** Compute C1 = D*C2 + R, in the given ring (assumes division is exact). */
 void coefficient_divrem(const lp_polynomial_context_t* ctx, coefficient_t* D, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2);
 
diff --git a/src/polynomial/gcd.c b/src/polynomial/gcd.c
index fd10bb9e..7aa033f3 100644
--- a/src/polynomial/gcd.c
+++ b/src/polynomial/gcd.c
@@ -183,7 +183,7 @@ void coefficient_gcd_pp_euclid(const lp_polynomial_context_t* ctx, coefficient_t
   STAT_INCR(coefficient, gcd_pp_euclid)
 
   if (trace_is_enabled("coefficient::gcd")) {
-    tracef("gcd\n")
+    tracef("gcd\n");
     tracef("P = "); coefficient_print(ctx, P, trace_out); tracef("\n");
     tracef("Q = "); coefficient_print(ctx, Q, trace_out); tracef("\n");
   }
@@ -257,7 +257,7 @@ void coefficient_gcd_pp_subresultant(const lp_polynomial_context_t* ctx, coeffic
   STAT_INCR(coefficient, gcd_pp_subresultant)
 
   if (trace_is_enabled("coefficient::gcd")) {
-    tracef("gcd\n")
+    tracef("gcd\n");
     tracef("P = "); coefficient_print(ctx, P, trace_out); tracef("\n");
     tracef("Q = "); coefficient_print(ctx, Q, trace_out); tracef("\n");
   }
@@ -761,7 +761,7 @@ lp_polynomial_vector_t* coefficient_mgcd_pp_subresultant(const lp_polynomial_con
   coefficient_construct(ctx, &cont);
 
   if (trace_is_enabled("coefficient::mgcd")) {
-    tracef("mgcd\n")
+    tracef("mgcd\n");
     tracef("P = "); coefficient_print(ctx, &P, trace_out); tracef("\n");
     tracef("Q = "); coefficient_print(ctx, &Q, trace_out); tracef("\n");
   }
diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c
index 11d948f0..2403f30e 100644
--- a/src/polynomial/polynomial.c
+++ b/src/polynomial/polynomial.c
@@ -23,6 +23,7 @@
 #include <variable_list.h>
 
 #include "polynomial/polynomial.h"
+#include "polynomial/coefficient.h"
 
 #include "polynomial/gcd.h"
 #include "polynomial/factorization.h"
@@ -40,7 +41,7 @@
 #include <stdlib.h>
 #include <string.h>
 
-#define SWAP(type, x, y) { type tmp = x; x = y; y = tmp; }
+#define SWAP(type, x, y) ({ type tmp = x; x = y; y = tmp; })
 
 static
 void check_polynomial_assignment(const lp_polynomial_t* A, const lp_assignment_t* M, lp_variable_t x) {
@@ -60,7 +61,7 @@ void check_polynomial_assignment(const lp_polynomial_t* A, const lp_assignment_t
     lp_variable_t y = vars.list[i];
     if (x != y && lp_assignment_get_value(M, y)->type == LP_VALUE_NONE) {
       lp_polynomial_print(A, trace_out);
-      tracef("\n")
+      tracef("\n");
       assert(0);
     }
   }
@@ -172,8 +173,6 @@ void lp_polynomial_set_external(lp_polynomial_t* A) {
   }
 }
 
-#define SWAP(type, x, y) { type tmp = x; x = y; y = tmp; }
-
 void lp_polynomial_swap(lp_polynomial_t* A1, lp_polynomial_t* A2) {
   // Swap everything, but keep the external flags
   lp_polynomial_t tmp = *A1; *A1 = *A2; *A2 = tmp;
@@ -626,6 +625,31 @@ void lp_polynomial_prem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_
   }
 }
 
+void lp_polynomial_pdivrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2) {
+
+    if (trace_is_enabled("polynomial")) {
+        tracef("polynomial_pdirvrem("); lp_polynomial_print(D, trace_out); tracef(", "); lp_polynomial_print(R, trace_out); tracef(", "); lp_polynomial_print(A1, trace_out); tracef(", "); lp_polynomial_print(A2, trace_out); tracef(")\n");
+        lp_variable_order_print(
+                A1->ctx->var_order, A1->ctx->var_db,
+                trace_out);
+        tracef("\n");
+    }
+
+    assert(lp_polynomial_context_equal(A1->ctx, A2->ctx));
+
+    lp_polynomial_external_clean(A1);
+    lp_polynomial_external_clean(A2);
+
+    lp_polynomial_set_context(D, A1->ctx);
+    lp_polynomial_set_context(R, A1->ctx);
+
+    coefficient_pdivrem(D->ctx, &D->data, &R->data, &A1->data, &A2->data);
+
+    if (trace_is_enabled("polynomial")) {
+        tracef("polynomial_pdirvrem() => ("); lp_polynomial_print(D, trace_out); tracef(", "); lp_polynomial_print(R, trace_out); tracef(")\n");
+    }
+}
+
 void lp_polynomial_sprem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2) {
 
   if (trace_is_enabled("polynomial")) {
@@ -650,6 +674,31 @@ void lp_polynomial_sprem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp
   }
 }
 
+void lp_polynomial_spdivrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2) {
+
+    if (trace_is_enabled("polynomial")) {
+        tracef("lp_polynomial_spdirvrem("); lp_polynomial_print(D, trace_out); tracef(", "); lp_polynomial_print(R, trace_out); tracef(", "); lp_polynomial_print(A1, trace_out); tracef(", "); lp_polynomial_print(A2, trace_out); tracef(")\n");
+        lp_variable_order_print(
+                A1->ctx->var_order, A1->ctx->var_db,
+                trace_out);
+        tracef("\n");
+    }
+
+    assert(lp_polynomial_context_equal(A1->ctx, A2->ctx));
+
+    lp_polynomial_external_clean(A1);
+    lp_polynomial_external_clean(A2);
+
+    lp_polynomial_set_context(D, A1->ctx);
+    lp_polynomial_set_context(R, A1->ctx);
+
+    coefficient_spdivrem(D->ctx, &D->data, &R->data, &A1->data, &A2->data);
+
+    if (trace_is_enabled("polynomial")) {
+        tracef("lp_polynomial_spdirvrem() => ("); lp_polynomial_print(D, trace_out); tracef(", "); lp_polynomial_print(R, trace_out); tracef(")\n");
+    }
+}
+
 void lp_polynomial_divrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2) {
 
   if (trace_is_enabled("polynomial")) {
@@ -668,7 +717,7 @@ void lp_polynomial_divrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polyn
   lp_polynomial_set_context(D, A1->ctx);
   lp_polynomial_set_context(R, A1->ctx);
 
-  coefficient_divrem(D->ctx, &R->data, &R->data, &A1->data, &A2->data);
+  coefficient_divrem(D->ctx, &D->data, &R->data, &A1->data, &A2->data);
 
   if (trace_is_enabled("polynomial")) {
     tracef("polynomial_rem() => ("); lp_polynomial_print(D, trace_out); tracef(", "); lp_polynomial_print(R, trace_out); tracef(")\n");
@@ -754,10 +803,10 @@ void lp_polynomial_reduce(
   lp_polynomial_set_context(Q, ctx);
   lp_polynomial_set_context(R, ctx);
 
-  coefficient_reduce(ctx, &A->data, &B->data, &P->data, &Q->data, &R->data, 1);
+  coefficient_reduce(ctx, &A->data, &B->data, &P->data, &Q->data, &R->data, REMAINDERING_PSEUDO_DENSE);
 
   if (trace_is_enabled("polynomial")) {
-    tracef("polynomial_derivative() =>\n");
+    tracef("polynomial_reduce() =>\n");
     tracef("\t P = "); lp_polynomial_print(P, trace_out); tracef("\n");
     tracef("\t Q = "); lp_polynomial_print(Q, trace_out); tracef("\n");
     tracef("\t R = "); lp_polynomial_print(R, trace_out); tracef("\n");
@@ -1095,7 +1144,7 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t
   }
 
   if (trace_is_enabled("polynomial")) {
-    tracef("polynomial_root_isolate("); lp_polynomial_print(A, trace_out); tracef("): unsorted roots\n")
+    tracef("polynomial_root_isolate("); lp_polynomial_print(A, trace_out); tracef("): unsorted roots\n");
     for (i = 0; i < roots_tmp_size; ++ i) {
       tracef("%zu :", i); lp_value_print(roots_tmp + i, trace_out); tracef("\n");
     }
@@ -1106,7 +1155,7 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t
     qsort(roots_tmp, roots_tmp_size, sizeof(lp_value_t), lp_value_cmp_void);
 
     if (trace_is_enabled("polynomial")) {
-      tracef("polynomial_root_isolate("); lp_polynomial_print(A, trace_out); tracef("): sorted roots\n")
+      tracef("polynomial_root_isolate("); lp_polynomial_print(A, trace_out); tracef("): sorted roots\n");
       for (i = 0; i < roots_tmp_size; ++ i) {
         tracef("%zu :", i); lp_value_print(roots_tmp + i, trace_out); tracef("\n");
       }
diff --git a/src/polynomial/psc.c b/src/polynomial/psc.c
index 97f9d950..a068d70b 100644
--- a/src/polynomial/psc.c
+++ b/src/polynomial/psc.c
@@ -373,7 +373,7 @@ void coefficient_psc_optimized(const lp_polynomial_context_t* ctx, coefficient_t
   coefficient_pow(ctx, &s, coefficient_lc(Q), P_deg - Q_deg);
 
   if (trace_is_enabled("coefficient::resultant")) {
-    tracef("s = "); coefficient_print(ctx, &s, trace_out); tracef("\n")
+    tracef("s = "); coefficient_print(ctx, &s, trace_out); tracef("\n");
   }
 
   // Set the final position
diff --git a/src/polyxx/polynomial.cpp b/src/polyxx/polynomial.cpp
index 0624a4d6..45decf4d 100644
--- a/src/polyxx/polynomial.cpp
+++ b/src/polyxx/polynomial.cpp
@@ -296,6 +296,22 @@ namespace poly {
                         rhs.get_internal());
     return res;
   }
+  std::pair<Polynomial, Polynomial> p_div_rem(const Polynomial& lhs,
+                                              const Polynomial& rhs) {
+    Polynomial d(detail::context(lhs, rhs));
+    Polynomial r(detail::context(lhs, rhs));
+    lp_polynomial_pdivrem(d.get_internal(), r.get_internal(), lhs.get_internal(),
+                         rhs.get_internal());
+    return {d, r};
+  }
+  std::pair<Polynomial, Polynomial> sp_div_rem(const Polynomial& lhs,
+                                               const Polynomial& rhs) {
+    Polynomial d(detail::context(lhs, rhs));
+    Polynomial r(detail::context(lhs, rhs));
+    lp_polynomial_spdivrem(d.get_internal(), r.get_internal(), lhs.get_internal(),
+                          rhs.get_internal());
+    return {d, r};
+  }
   std::pair<Polynomial, Polynomial> div_rem(const Polynomial& lhs,
                                             const Polynomial& rhs) {
     Polynomial d(detail::context(lhs, rhs));
diff --git a/src/utils/debug_trace.h b/src/utils/debug_trace.h
index 5c0d3d39..cce99f26 100644
--- a/src/utils/debug_trace.h
+++ b/src/utils/debug_trace.h
@@ -28,7 +28,7 @@ FILE* trace_out_real;
 #define trace_out (trace_out_real ? trace_out_real : stderr)
 
 /** Print to the debug trace printf style */
-#define tracef(...) fprintf(trace_out, __VA_ARGS__);
+#define tracef(...) fprintf(trace_out, __VA_ARGS__)
 
 /** Set the output file for tracing */
 void trace_set_output(FILE* file);
@@ -41,17 +41,17 @@ void trace_disable(const char* tag);
 
 int trace_is_enabled(const char* tag);
 
-#define TRACE(tag, ...) { \
+#define TRACE(tag, ...) do { \
   if (trace_is_enabled(tag)) { \
     tracef(__VA_ARGS__); \
   } \
-}
+} while(0)
 
-#define TRACE_CMD(tag, cmd) { \
+#define TRACE_CMD(tag, cmd) do { \
   if (trace_is_enabled(tag)) { \
     cmd; \
   } \
-} \
+} while(0)
 
 #else
 
@@ -61,6 +61,6 @@ int trace_is_enabled(const char* tag) {
   return 0;
 }
 
-#define TRACE(tag, ...)
-#define TRACE_CMD(tag, cmd)
+#define TRACE(tag, ...) ((void)0)
+#define TRACE_CMD(tag, cmd) ((void)0)
 #endif
diff --git a/test/python/check.py b/test/python/check.py
index 174d85c6..5cf7d50c 100755
--- a/test/python/check.py
+++ b/test/python/check.py
@@ -14,7 +14,6 @@ def forkexec(test, env):
 
 
 def main():
-
     parser = argparse.ArgumentParser(description='Process some integers.')
     parser.add_argument('--sympy', action="store_true")
     parser.add_argument('--stats', action="store_true")
@@ -36,15 +35,14 @@ def main():
              "tests/polynomial_feasibility.py", 
              "tests/value.py"]
 
-    if (args.sympy):
+    if args.sympy:
         print("Sympy checking enabled")
         polypy_test.sympy_checker.enabled = True
     else:
         print("Sympy checking disabled")
         polypy_test.sympy_checker.enabled = False
 
-
-
+    failed = 0
     for test in tests:
         print("Running {0}:".format(test))
         context = dict()
@@ -52,11 +50,14 @@ def main():
         module = context["polypy_test"]
         print("PASS: {0}".format(module.PASS))
         print("FAIL: {0}".format(module.FAIL))
+        failed += module.FAIL
 
-    if (args.stats):
+    if args.stats:
         print("Statistics:")
         polypy.stats_print()
 
+    return 1 if failed > 0 else 0
+
 
 if __name__ == '__main__':
     sys.exit(main())
diff --git a/test/python/tests/algebraic_number.py b/test/python/tests/algebraic_number.py
index f6437d7d..ebad1647 100755
--- a/test/python/tests/algebraic_number.py
+++ b/test/python/tests/algebraic_number.py
@@ -159,4 +159,4 @@ def check_comparison(a1, a2, cmp, result, expected):
     p = functools.reduce(lambda x, y: x*y, sample, one)
     random.shuffle(sample)
     p = functools.reduce(lambda x, y: x/y, sample, p)
-    polypy_test.check(p == 1)
+    polypy_test.check(p == one)
diff --git a/test/python/tests/polynomial_arithmetic.py b/test/python/tests/polynomial_arithmetic.py
index da62b7ea..e04c53a4 100755
--- a/test/python/tests/polynomial_arithmetic.py
+++ b/test/python/tests/polynomial_arithmetic.py
@@ -3,32 +3,46 @@
 import polypy
 import polypy_test
 
+
 def check_binary(op, op_name, p, q, expected):
     result = op(p, q)
     ok = result == expected
-    if (not ok):
+    if not ok:
+        print("p = {0}".format(p))
+        print("q = {0}".format(q))
+        print("{0} = {1}".format(op_name, result))
+        print("expected = {0}".format(expected))
+    polypy_test.check(ok)
+
+
+def check_binary_tuple(op, op_name, p, q, n, expected):
+    result = (op(p, q))[n]
+    ok = result == expected
+    if not ok:
         print("p = {0}".format(p))
         print("q = {0}".format(q))
         print("{0} = {1}".format(op_name, result))
         print("expected = {0}".format(expected))
     polypy_test.check(ok)
 
+
 def check_unary(op, op_name, p, expected):
     result = op(p)
     ok = result == expected
-    if (not ok):
+    if not ok:
         print("p = {0}".format(p))
         print("{0} = {1}".format(op_name, result))
         print("expected = {0}".format(expected))
     polypy_test.check(ok)
 
+
 polypy_test.init()
 
 polypy_test.start("Addition")
 
-x = polypy.Variable("x");
-y = polypy.Variable("y");
-z = polypy.Variable("z");
+x = polypy.Variable("x")
+y = polypy.Variable("y")
+z = polypy.Variable("z")
 
 polypy.variable_order.set([z, y, x])
 
@@ -164,56 +178,91 @@ def poly_div(p, q):
 def poly_rem(p, q):
     return p % q
 
+def poly_prem(p, q):
+    return p.prem(q)
+
+def poly_divmod(p, q):
+    return divmod(p, q)
+
+def poly_pdivrem(p, q):
+    return p.pdivrem(q)
+
+
 p = 2*x + 2
 q = 2
 expected = x + 1
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 p = x**2 - 1
 q = x - 1
 expected = x + 1
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 p = x**2 - 1
 q = x - 1
 expected = 0
 check_binary(poly_rem, "rem", p, q, expected)
+check_binary(poly_prem, "rem", p, q, expected)
+check_binary_tuple(poly_divmod, "rem", p, q, 1, expected)
+check_binary_tuple(poly_pdivrem, "rem", p, q, 1, expected)
 
 p = (1*z**11 - 1*z**5)
 q = (z**5)
 expected = z**6 - 1
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 p = x*y*z
 q = (x*y)
 expected = z
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 p = 6*y*(x + 1)
 q = (2*y)
 expected = 3*(x + 1)
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 p = 6*y*(x + 1)
 q = 3*(x+1)
 expected = 2*y
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 q = x + 1
 expected = 6*y
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 p = 2*x**2 + 4*x + 6
 q = 2
 expected = x**2 + 2*x + 3
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
 
 p = 2*x**2 + 4*x + 6
 q = (x + 2) - x
 expected = x**2 + 2*x + 3
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
+
+expected = 0
+check_binary(poly_rem, "rem", p, q, expected)
+check_binary(poly_prem, "rem", p, q, expected)
+check_binary_tuple(poly_divmod, "rem", p, q, 1, expected)
+check_binary_tuple(poly_pdivrem, "rem", p, q, 1, expected)
 
 p = (x + 6) - x
 q = (x + 3) - x
 expected = 2
 check_binary(poly_div, "div", p, q, expected)
+check_binary_tuple(poly_divmod, "div", p, q, 0, expected)
+
+p = 2*x**2 + 4*x + 6
+q = (x + 2)
+expected = 6
+check_binary(poly_prem, "rem", p, q, expected)
+check_binary_tuple(poly_pdivrem, "rem", p, q, 1, expected)