Fixed divrem bugs, added pseudo divrem (#72)

* Added assertion in coefficient_div to warn about problem of dividing a constant by a polynomial. * Minor fixes in polynomial.c. * Added pseudo div and rem functions. * fixed various typos and minor code issues * Updated trace macros * removed function prototype with missing implementation from header file * Bump version * fixed typo * removed duplicate macro definition * updated description on lp_polynomial_reduce and coefficient_reduce * Undo patch version bump * added new functions to python bindings and updated tests * fixed algebraic number tests and added return code on failed tests * minor format --------- Co-authored-by: Thomas Hader <[email protected]>
SRI-CSL · Nov 3, 2023 · 862d7cb · 862d7cb
1 parent 5a85023
commit 862d7cb
Show file tree

Hide file tree

Showing 15 changed files with 387 additions and 82 deletions.
diff --git 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)
diff --git 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
@@ -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
@@ -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,33 +784,58 @@ 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);
   PyTuple_SetItem(pair, 1, remObj);
   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
@@ -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);
   }