Rabin root finding (#78)

* extended assertion in berlekamp factorization

* some code cleanup (made functions static, added assertions)

* Added lp_upolynomial_roots_find_Zp and

* Added comment with reference

* Added todos and comments

* Added C++ code for find_roots_Zp

* Added python lib and test for root_finding_Zp

* Revert changes to libpoly0.symbols

* Removed (already fixed) TODO

* Fixed typo

* Fixed comment

* added missing comments

---------

Co-authored-by: Thomas Hader <[email protected]>

include/integer.h

@@ -168,7 +168,7 @@ int lp_integer_sgn(const lp_int_ring_t* K, const lp_integer_t* c);
 int lp_integer_cmp(const lp_int_ring_t* K, const lp_integer_t* c, const lp_integer_t* to);
 
 /**
- * Compare the two integers in the ring. Same as for sgn.
+ * Compare the two integers in the ring. Same as for cmp.
  */
 int lp_integer_cmp_int(const lp_int_ring_t* K, const lp_integer_t* c, long to);
 

include/polyxx/upolynomial.h

@@ -235,4 +235,9 @@ namespace poly {
    */
   std::vector<AlgebraicNumber> isolate_real_roots(const UPolynomial& p);
 
+  /**
+   * Finds the roots for a polynomial mod p using rabin root finding.
+   */
+  std::vector<Integer> find_roots_Zp(const UPolynomial& p);
+
 }  // namespace poly

include/upolynomial.h

@@ -191,10 +191,20 @@ void lp_upolynomial_subst_x_pow_in_place(lp_upolynomial_t* f, size_t n);
 lp_upolynomial_t* lp_upolynomial_neg(const lp_upolynomial_t* p);
 
 /**
- * Negates p in place.
+ * Negates p in-place.
  */
 void lp_upolynomial_neg_in_place(lp_upolynomial_t* p);
 
+/**
+ * Makes the polynomial monic.
+ */
+lp_upolynomial_t* lp_upolynomial_make_monic(const lp_upolynomial_t* p);
+
+/**
+ * Makes the polynomial monic in-place.
+ */
+void lp_upolynomial_make_monic_in_place(lp_upolynomial_t* p);
+
 /**
  * Add two polynomials (all operations in the same ring).
  */
@@ -365,7 +375,14 @@ int lp_upolynomial_roots_count(const lp_upolynomial_t* p, const lp_rational_inte
 void lp_upolynomial_roots_isolate(const lp_upolynomial_t* p, lp_algebraic_number_t* roots, size_t* roots_size);
 
 /**
- * Reverses the coefficient of p in place. The result polynomial is
+ * Finds the root for an univariate polynomial in Zp. The array roots will be allocated,
+ * and the user should de-allocate it. The size parameter will be updated with the size
+ * of the array.
+ */
+void lp_upolynomial_roots_find_Zp(const lp_upolynomial_t* f, lp_integer_t** roots, size_t* roots_size);
+
+/**
+ * Reverses the coefficient of p in-place. The result polynomial is
  * p'(x) = x^n * p(1/x) = a_n + ... + a_0 * x^n
  */
 void lp_upolynomial_reverse_in_place(lp_upolynomial_t* p);

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/CMakeLists.txt

@@ -24,6 +24,7 @@ set(poly_SOURCES
   upolynomial/factors.c
   upolynomial/factorization.c
   upolynomial/root_finding.c
+  upolynomial/upolynomial_vector.c
   polynomial/monomial.c
   polynomial/coefficient.c
   polynomial/output.c

src/number/integer.h

@@ -184,6 +184,16 @@ int integer_is_zero(const lp_int_ring_t* K, const lp_integer_t* c) {
   }
 }
 
+static inline
+int integer_is_odd(const lp_integer_t* c) {
+  return mpz_odd_p(c);
+}
+
+static inline
+int integer_is_even(const lp_integer_t* c) {
+  return mpz_even_p(c);
+}
+
 static inline
 int integer_sgn(const lp_int_ring_t* K, const lp_integer_t* c) {
   if (K) {

src/polynomial/polynomial.c

@@ -261,6 +261,7 @@ int lp_polynomial_is_univariate(const lp_polynomial_t* A) {
   return coefficient_is_univariate(&A->data);
 }
 
+// TODO is a constant polynomial univariate?
 int lp_polynomial_is_univariate_m(const lp_polynomial_t* A, const lp_assignment_t* m) {
   if (lp_polynomial_is_constant(A)) {
     return 0;
@@ -286,6 +287,7 @@ int lp_polynomial_is_univariate_m(const lp_polynomial_t* A, const lp_assignment_
 }
 
 lp_upolynomial_t* lp_polynomial_to_univariate(const lp_polynomial_t* A) {
+  lp_polynomial_external_clean(A);
   if (!(coefficient_is_constant(&A->data) || coefficient_is_univariate(&A->data))) {
     return NULL;
   } else {

src/polyxx/upolynomial.cpp

@@ -342,5 +342,19 @@ namespace poly {
     return res;
   }
 
+  std::vector<Integer> find_roots_Zp(const UPolynomial& p) {
+    lp_integer_t *roots;
+    std::size_t roots_size;
+    lp_upolynomial_roots_find_Zp(p.get_internal(), &roots, &roots_size);
+    std::vector<Integer> res;
+    for (std::size_t i = 0; i < roots_size; ++i) {
+      res.emplace_back(&roots[i]);
+    }
+    for (std::size_t i = 0; i < roots_size; ++i) {
+      lp_integer_destruct(&roots[i]);
+    }
+    free(roots);
+    return res;
+  }
 
 }  // namespace poly