Support for subresultant generation and features for Zp

include/assignment.h

@@ -65,9 +65,19 @@ const lp_value_t* lp_assignment_get_value(const lp_assignment_t* m, lp_variable_
 /** Get an approximate value of the variable */
 void lp_assignment_get_value_approx(const lp_assignment_t* m, lp_variable_t x, lp_rational_interval_t* approx);
 
+/** Returns 1 if the variable is set, 0 otherwise. */
+int lp_assignment_is_set(const lp_assignment_t* m, lp_variable_t x);
+
+/** Returns 1 if the assignment assigns integer values only */
+int lp_assignment_is_integer(const lp_assignment_t *m);
+
+/** Returns 1 if the assignment is empty */
+int lp_assignment_is_empty(const lp_assignment_t *m);
+
 /** Get the sign of the polynomial in the model */
 int lp_assignment_sgn(const lp_assignment_t* m, const lp_polynomial_t* A);
 
+
 struct lp_interval_assignment_struct {
   /** Size of the map */
   size_t size;

include/polynomial.h

@@ -103,6 +103,9 @@ int lp_polynomial_lc_is_constant(const lp_polynomial_t* A);
 /** In case lc is constant, this returns the sign */
 int lp_polynomial_lc_sgn(const lp_polynomial_t* A);
 
+/** Returns the constant part of the leading coefficient */
+void lp_polynomial_lc_constant(const lp_polynomial_t* A, lp_integer_t *out);
+
 /** Get the context of the given polynomial */
 const lp_polynomial_context_t* lp_polynomial_get_context(const lp_polynomial_t* A);
 
@@ -136,6 +139,18 @@ int lp_polynomial_is_univariate_m(const lp_polynomial_t* A, const lp_assignment_
 /** Returns the univariate polynomial (if univariate, or 0 otherwise) */
 lp_upolynomial_t* lp_polynomial_to_univariate(const lp_polynomial_t* A);
 
+/**
+ * Returns the univariate polynomial of A partially evaluated with m.
+ * m must assign all variables in A (but the top variable) to an integer value.
+ */
+lp_upolynomial_t* lp_polynomial_to_univariate_m(const lp_polynomial_t* A, const lp_assignment_t* m);
+
+/** Returns true if the polynomial is a monomial */
+int lp_polynomial_is_monomial(const lp_polynomial_t* A);
+
+/** Returns A as a monomial. Assumes that A is a monomial */
+void lp_polynomial_to_monomial(const lp_polynomial_t* A, lp_monomial_t* out);
+
 /** Returns true if all of polynomial's variables are assigned */
 int lp_polynomial_is_assigned(const lp_polynomial_t* A, const lp_assignment_t* m);
 
@@ -145,9 +160,11 @@ int lp_polynomial_sgn(const lp_polynomial_t* A, const lp_assignment_t* m);
 /** Returns the interval approximation of the polynomial value */
 void lp_polynomial_interval_value(const lp_polynomial_t* A, const lp_interval_assignment_t* m, lp_interval_t* result);
 
-/** returns the sign of the polynomial in the model */
+/** returns the value of the polynomial in the model */
 lp_value_t* lp_polynomial_evaluate(const lp_polynomial_t* A, const lp_assignment_t* m);
 
+/** copies the integer value of a polynomial to out; assignment m must only assign integer values */
+void lp_polynomial_evaluate_integer(const lp_polynomial_t* A, const lp_assignment_t* m, lp_integer_t* out);
 
 /**
  * Compare the two polynomials in the ring. Not necessarily +/- 1, could be
@@ -273,12 +290,19 @@ void lp_polynomial_pp_cont(lp_polynomial_t* pp, lp_polynomial_t* cont, const lp_
 void lp_polynomial_resultant(lp_polynomial_t* res, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
 
 /**
- * Compute the principal subresultant coefficients (psc) of A1 and A1. Bot A1
+ * Compute the principal subresultant coefficients (psc) of A1 and A2. Both A1
  * and A2 must be (non-trivial) polynomials over the same variable. If
- *  deg(A1) = m, deg(A2) = n, and output will be of size min(m, n).
+ * deg(A1) = m, deg(A2) = n, and output will be of size min(m, n)+1.
  */
 void lp_polynomial_psc(lp_polynomial_t** psc, const lp_polynomial_t* A1, const lp_polynomial_t* A2);
 
+/**
+ * Compute the subresultant regular sub-chain (subres) of A1 and A2. Bot A1
+ * and A2 must be (non-trivial) polynomials over the same variable. If
+ * deg(A1) = m, deg(A2) = n, and output will be of size min(m, n)+1.
+ */
+void lp_polynomial_subres(lp_polynomial_t** subres, const lp_polynomial_t* A, const lp_polynomial_t* B);
+
 /**
  * Compute the model-based GCD of the two polynomials. Adds the assumptions to
  * the vector.
@@ -340,6 +364,12 @@ lp_polynomial_t* lp_polynomial_constraint_explain_infer_bounds(const lp_polynomi
  */
 int lp_polynomial_constraint_evaluate(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, const lp_assignment_t* M);
 
+/**
+ * Given a polynomial constraint over an integer ring Zp, evaluate its truth value.
+ * sgn_condition must either be (== 0) or (!= 0) and M must assign to Zp.
+ */
+int lp_polynomial_constraint_evaluate_Zp(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, const lp_assignment_t* M);
+
 /**
  * Given a polynomial A(x1, ..., xn, y) with y being the top variable, a root index,
  * a sign condition, and an assignment M that assigns x1, ..., xn, the function
@@ -383,7 +413,11 @@ int lp_polynomial_constraint_resolve_fm(
     lp_polynomial_t* R, lp_sign_condition_t* R_sgn,
     lp_polynomial_vector_t* assumptions);
 
-
+/**
+ * Reduces the degree of the polynomials wrt. to the polynomial's int_ring_t K.
+ * E.g. for K = x mod 5, x = x^5 for every value in K
+ */
+void lp_polynomial_reduce_degree_Zp(lp_polynomial_t *R, const lp_polynomial_t *A);
 
 #ifdef __cplusplus
 } /* close extern "C" { */

include/polyxx/polynomial.h

@@ -115,6 +115,10 @@ namespace poly {
   /** Check if the given polynomial is completely assigned over a given
    * assignment. */
   bool is_assigned_over_assignment(const Polynomial& p, const Assignment& a);
+  /** Evaluates p over a given assignment and returns an univariate polynomial.
+   * Assumes that a assigns all variable in p but the top variable.
+   * Assumes that a assigns to integer only. */
+  UPolynomial to_univariate(const Polynomial& p, const Assignment& a);
   /** Compute the sign of a polynomial over an assignment. */
   int sgn(const Polynomial& p, const Assignment& a);
   /** Evaluate a polynomial over an assignment. */

include/sign_condition.h

@@ -53,7 +53,7 @@ typedef enum lp_sign_condition_enum lp_sign_condition_t;
 int lp_sign_condition_consistent(lp_sign_condition_t sgn_condition, int sign);
 
 /**
- * CHeck if the sign condition is consistent with the given interval.
+ * Check if the sign condition is consistent with the given interval.
  */
 int lp_sign_condition_consistent_interval(lp_sign_condition_t sgn_condition, const lp_interval_t* I);
 
@@ -62,6 +62,11 @@ int lp_sign_condition_consistent_interval(lp_sign_condition_t sgn_condition, con
  */
 lp_sign_condition_t lp_sign_condition_negate(lp_sign_condition_t sgn_condition);
 
+/**
+ * Check if the sign condition is valid for a integer ring with a non-zero characteristic.
+ */
+int lp_sign_condition_Zp_valid(lp_sign_condition_t sgn_condition);
+
 /**
  * Print the sign condition.
  */

python/polypyPolynomial3.c

@@ -150,6 +150,9 @@ Polynomial_resultant(PyObject* self, PyObject* args);
 static PyObject*
 Polynomial_psc(PyObject* self, PyObject* args);
 
+static PyObject*
+Polynomial_subres(PyObject* self, PyObject* args);
+
 static PyObject*
 Polynomial_mgcd(PyObject* self, PyObject* args);
 
@@ -199,6 +202,7 @@ PyMethodDef Polynomial_methods[] = {
     {"derivative", (PyCFunction)Polynomial_derivative, METH_NOARGS, "Returns the derivative of the polynomial"},
     {"resultant", (PyCFunction)Polynomial_resultant, METH_VARARGS, "Returns the resultant of the current and given polynomial"},
     {"psc", (PyCFunction)Polynomial_psc, METH_VARARGS, "Returns the principal subresultant coefficients of the current and given polynomial"},
+    {"subres", (PyCFunction) Polynomial_subres, METH_VARARGS, "Returns the subresultant chain of the current and given polynomial"},
     {"mgcd", (PyCFunction)Polynomial_mgcd, METH_VARARGS, "Returns assumptions that the GCD of two polynomials is of same degree"},
     {"factor_square_free", (PyCFunction)Polynomial_factor_square_free, METH_NOARGS, "Returns the square-free factorization of the polynomial"},
     {"evaluate", (PyCFunction)Polynomial_evaluate, METH_VARARGS, "Returns the value of the polynomial in the given assignment (or null if it doesn't fully evaluate"},
@@ -1023,8 +1027,13 @@ Polynomial_lcm(PyObject* self, PyObject* args) {
   return Polynomial_create(lcm);
 }
 
+enum subres_type {
+  SUBRES,
+  SUBRES_PSC
+};
+
 static PyObject*
-Polynomial_psc(PyObject* self, PyObject* args) {
+Polynomial_subres_impl(PyObject* self, PyObject* args, enum subres_type type) {
 
   // self is always a polynomial
   Polynomial* p1 = (Polynomial*) self;
@@ -1053,7 +1062,7 @@ Polynomial_psc(PyObject* self, PyObject* args) {
     }
   }
 
-  // Othe polynomial
+  // Other polynomial
   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)) {
@@ -1077,21 +1086,30 @@ Polynomial_psc(PyObject* self, PyObject* args) {
   size_t p2_deg = lp_polynomial_degree(p2->p);
   int size = p1_deg > p2_deg ? p2_deg + 1 : p1_deg + 1;
 
-  lp_polynomial_t** psc = malloc(sizeof(lp_polynomial_t*)*size);
+  lp_polynomial_t** S = malloc(sizeof(lp_polynomial_t*)*size);
   int i;
   for (i = 0; i < size; ++ i) {
-    psc[i] = lp_polynomial_new(p1_ctx);
+    S[i] = lp_polynomial_new(p1_ctx);
   }
 
-  // Compute the psc
-  lp_polynomial_psc(psc, p1->p, p2->p);
+  switch (type) {
+    case SUBRES:
+      // Compute the full subres
+      lp_polynomial_subres(S, p1->p, p2->p);
+      break;
+    case SUBRES_PSC:
+      // Compute the psc
+      lp_polynomial_psc(S, p1->p, p2->p);
+      break;
+  }
 
   // Copy the polynomials into a list
   PyObject* list = PyList_New(size);
   for (i = 0; i < size; ++i) {
-    PyObject* p = Polynomial_create(psc[i]);
+    PyObject* p = Polynomial_create(S[i]);
     PyList_SetItem(list, i, p);
   }
+  free(S);
 
   if (dec_other) {
     Py_DECREF(other);
@@ -1101,6 +1119,16 @@ Polynomial_psc(PyObject* self, PyObject* args) {
   return list;
 }
 
+static PyObject*
+Polynomial_psc(PyObject* self, PyObject* args) {
+  return Polynomial_subres_impl(self, args, SUBRES_PSC);
+}
+
+static PyObject*
+Polynomial_subres(PyObject* self, PyObject* args) {
+  return Polynomial_subres_impl(self, args, SUBRES);
+}
+
 static PyObject*
 Polynomial_mgcd(PyObject* self, PyObject* args) {
 
@@ -1208,7 +1236,7 @@ Polynomial_resultant(PyObject* self, PyObject* args) {
     }
   }
 
-  // Othe polynomial
+  // Other polynomial
   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)) {
@@ -1241,7 +1269,6 @@ Polynomial_resultant(PyObject* self, PyObject* args) {
   return Polynomial_create(resultant);
 }
 
-
 static PyObject*
 Polynomial_extended_gcd(PyObject* self, PyObject* args) {
   return 0;

src/CMakeLists.txt

@@ -28,7 +28,7 @@ set(poly_SOURCES
   polynomial/coefficient.c
   polynomial/output.c
   polynomial/gcd.c
-  polynomial/psc.c
+  polynomial/subres.c
   polynomial/factorization.c
   polynomial/polynomial.c
   polynomial/polynomial_context.c