Added python lib and test for root_finding_Zp

python/polypyUPolynomial3.c

@@ -72,6 +72,9 @@ UPolynomial_roots_isolate(PyObject* self);
 static PyObject*
 UPolynomial_sturm_sequence(PyObject* self);
 
+static PyObject*
+UPolynomial_roots_find_Zp(PyObject* self);
+
 static PyObject*
 UPolynomial_evaluate(PyObject* self, PyObject* args);
 
@@ -120,6 +123,7 @@ PyMethodDef UPolynomial_methods[] = {
     {"roots_count", (PyCFunction)UPolynomial_roots_count, METH_VARARGS, "Returns the number of real roots in the given interval"},
     {"roots_isolate", (PyCFunction)UPolynomial_roots_isolate, METH_NOARGS, "Returns the list of real roots"},
     {"sturm_sequence", (PyCFunction)UPolynomial_sturm_sequence, METH_NOARGS, "Returns the Sturm sequence"},
+    {"roots_find_Zp", (PyCFunction)UPolynomial_roots_find_Zp, METH_NOARGS, "Returns the roots of the polynomial in Zp"},
     {"derivative", (PyCFunction)UPolynomial_derivative, METH_NOARGS, "Returns the derivative of the polynomial"},
     {"evaluate", (PyCFunction)UPolynomial_evaluate, METH_VARARGS, "Returns the value of the polynomial at the given point"},
     {NULL}  /* Sentinel */
@@ -863,6 +867,36 @@ UPolynomial_sturm_sequence(PyObject* self) {
   return result;
 }
 
+static PyObject* integer_list_to_PyList(lp_integer_t *list, size_t size) {
+  // Construct the result
+  PyObject* pylist = PyList_New(size);
+
+  // Copy over the integers
+  size_t i;
+  for (i = 0; i < size; ++ i) {
+    PyObject* pylist_i = integer_to_PyLong(&list[i]);
+    Py_INCREF(pylist_i);
+    PyList_SetItem(pylist, i, pylist_i);
+  }
+
+  // Return the list
+  return pylist;
+}
+
+static PyObject*
+UPolynomial_roots_find_Zp(PyObject* self) {
+  lp_upolynomial_t* p = ((UPolynomialObject*) self)->p;
+  lp_integer_t* roots;
+  size_t roots_size;
+  lp_upolynomial_roots_find_Zp(p, &roots, &roots_size);
+  PyObject* result = integer_list_to_PyList(roots, roots_size);
+  for (size_t i = 0; i < roots_size; ++i) {
+    lp_integer_destruct(&roots[i]);
+  }
+  free(roots);
+  return result;
+}
+
 static PyObject*
 UPolynomial_evaluate(PyObject* self, PyObject* args) {
   if (PyTuple_Check(args) && PyTuple_Size(args) == 1) {

src/upolynomial/root_finding.h

@@ -55,6 +55,6 @@ void upolynomial_roots_isolate_sturm(const lp_upolynomial_t* f, lp_algebraic_num
 
 /**
 * Finds the roots for a polynomial over Zp. Uses brute force or rabin root
-* finding, depending on p.
+* finding, depending on p. Degree of p must be positive.
 */
 void upolynomial_roots_find_Zp(const lp_upolynomial_t* f, lp_integer_t** roots, size_t* roots_size);

test/python/tests/upolynomial_roots.py

@@ -113,3 +113,41 @@ def check_roots_count(p, expected, lb = None, ub = None):
 check_roots_count(p, 2, -2, 2)
 check_roots_count(p, 4, -2.5, 2.5)
 check_roots_count(p, 4, -3, 3)
+
+
+# checking root finding mod p
+polypy_test.start("Root finding in Zp")
+
+
+def check_roots_Zp(p, expected):
+    zeros = p.roots_find_Zp()
+    if set(zeros) != set(expected):
+        polypy_test.check(False)
+        print("p = {0}".format(p))
+        print("got = {0}".format(','.join(map(str, zeros))))
+        print("expected = {0}".format(','.join(map(str, expected))))
+    else:
+        polypy_test.check(True)
+
+
+primes = [13, 10007, 1230127]
+P = lambda x: [
+    (x - 1),
+    (x - 1) * (x + 1),
+    (x - 1) * (x + 1) * (x - 2),
+    (x - 1) * (x + 1) * (x - 2) * (x + 2),
+    (x - 1) ** 3 * (x + 1) ** 2 * (x - 2) ** 2 * (x + 2)
+]
+expected = [
+    [1],
+    [1, -1],
+    [1, -1, 2],
+    [1, -1, 2, -2],
+    [1, -1, 2, -2]
+]
+
+for prime in primes:
+    K = polypy.CoefficientRing(prime)
+    x = polypy.x.to_ring(K)
+    for i, p in enumerate(P(x)):
+        check_roots_Zp(p, expected[i])