diff --git a/include/integer.h b/include/integer.h
index 0e18e155..affcca51 100644
--- a/include/integer.h
+++ b/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);
 
diff --git a/include/polyxx/upolynomial.h b/include/polyxx/upolynomial.h
index 850a7f17..f0a5ac56 100644
--- a/include/polyxx/upolynomial.h
+++ b/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
diff --git a/include/upolynomial.h b/include/upolynomial.h
index 51ad71e1..04744347 100644
--- a/include/upolynomial.h
+++ b/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);
diff --git a/python/polypyUPolynomial3.c b/python/polypyUPolynomial3.c
index 38fb2f05..e3b65056 100644
--- a/python/polypyUPolynomial3.c
+++ b/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) {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 2fc95935..7c790e18 100755
--- a/src/CMakeLists.txt
+++ b/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
diff --git a/src/number/integer.h b/src/number/integer.h
index ab0efbfe..dc0d6645 100644
--- a/src/number/integer.h
+++ b/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) {
diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c
index 91976367..80c07c21 100644
--- a/src/polynomial/polynomial.c
+++ b/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 {
diff --git a/src/polyxx/upolynomial.cpp b/src/polyxx/upolynomial.cpp
index 234c3242..90c21c1f 100644
--- a/src/polyxx/upolynomial.cpp
+++ b/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
diff --git a/src/upolynomial/factorization.c b/src/upolynomial/factorization.c
index d2673c63..44d9efd0 100644
--- a/src/upolynomial/factorization.c
+++ b/src/upolynomial/factorization.c
@@ -31,6 +31,10 @@
 #include <assert.h>
 #include <stdlib.h>
 
+/* Limits on the prime field order p when performing operations with O(p) time and
+ * space complexity, respectively. Assertions will fail when the limits are exceeded. */
+#define FIELD_ORDER_LIMIT 1000
+#define FIELD_ORDER_LIMIT_BERLEKAMP 250
 
 STAT_DECLARE(int, upolynomial, factor_square_free)
 STAT_DECLARE(int, upolynomial, factor_distinct_degree)
@@ -77,7 +81,7 @@ STAT_DECLARE(int, upolynomial, factor_berlekamp_square_free)
  * The loop ends when L = 1 and then we know that P = \prod_{p \div k} f_k^k
  * which we know how to special case (if P != 1).
  */
-lp_upolynomial_factors_t* lp_upolynomial_factor_square_free_primitive(const lp_upolynomial_t* f) {
+lp_upolynomial_factors_t* upolynomial_factor_square_free_primitive(const lp_upolynomial_t* f) {
 
   if (trace_is_enabled("factorization")) {
     tracef("upolynomial_factor_square_free("); lp_upolynomial_print(f, trace_out); tracef(")\n");
@@ -107,11 +111,11 @@ lp_upolynomial_factors_t* lp_upolynomial_factor_square_free_primitive(const lp_u
     // f' is zero for a non-zero polynomial => f has to be of the form
     // f = \sum a_k x^(p*d_k) = f_p(x^p) where f_p = \sum a_k x^d_k
     // we factor f_p and then return f_p(x^p)=(f_p)^p
-    int p = integer_to_int(&f->K->M);
+    assert(mpz_fits_slong_p(&f->K->M));
+    long p = integer_to_int(&f->K->M);
     lp_upolynomial_t* f_p = lp_upolynomial_div_degrees(f, p);
-    factors = lp_upolynomial_factor_square_free_primitive(f_p);
-    size_t i;
-    for (i = 0; i < factors->size; ++ i) {
+    factors = upolynomial_factor_square_free_primitive(f_p);
+    for (size_t i = 0; i < factors->size; ++ i) {
       factors->multiplicities[i] *= p;
     }
     lp_upolynomial_delete(f_p);
@@ -167,9 +171,10 @@ lp_upolynomial_factors_t* lp_upolynomial_factor_square_free_primitive(const lp_u
 
     // If P has content, it is a power of p
     if (lp_upolynomial_degree(P) > 0) {
-      int p = integer_to_int(&f->K->M);
+      assert(mpz_fits_slong_p(&f->K->M));
+      long p = integer_to_int(&f->K->M);
       lp_upolynomial_t* P_p = lp_upolynomial_div_degrees(P, p);
-      lp_upolynomial_factors_t* sub_factors = lp_upolynomial_factor_square_free_primitive(P_p);
+      lp_upolynomial_factors_t* sub_factors = upolynomial_factor_square_free_primitive(P_p);
       size_t i;
       for (i = 0; i < sub_factors->size; ++ i) {
         lp_upolynomial_factors_add(factors, sub_factors->factors[i], sub_factors->multiplicities[i] * p);
@@ -220,7 +225,7 @@ lp_upolynomial_factors_t* lp_upolynomial_factor_square_free(const lp_upolynomial
   // Take out the power of x^k
   if (lp_upolynomial_const_term(f_pp)) {
     // Get a square-free decomposition of f
-    sq_free_factors = lp_upolynomial_factor_square_free_primitive(f_pp);
+    sq_free_factors = upolynomial_factor_square_free_primitive(f_pp);
   } else {
     // Get a copy without the power of x
     lp_upolynomial_t* f_pp_nonzero = lp_upolynomial_construct_copy(f_pp);
@@ -230,7 +235,7 @@ lp_upolynomial_factors_t* lp_upolynomial_factor_square_free(const lp_upolynomial
       f_pp_nonzero->monomials[i].degree -= x_degree;
     }
     // Get a square-free decomposition of f
-    sq_free_factors = lp_upolynomial_factor_square_free_primitive(f_pp_nonzero);
+    sq_free_factors = upolynomial_factor_square_free_primitive(f_pp_nonzero);
     // Add x^k to the factorization
     lp_upolynomial_t* x_upoly = lp_upolynomial_construct_power(f->K, 1, 1);
     lp_upolynomial_factors_add(sq_free_factors, x_upoly, x_degree);
@@ -283,8 +288,9 @@ lp_upolynomial_factors_t* upolynomial_factor_distinct_degree(const lp_upolynomia
   assert(lp_upolynomial_is_monic(f));
 
   // The prime
-  int p = integer_to_int(&K->M);
-  assert(p < 100);
+  assert(mpz_fits_slong_p(&K->M));
+  long p = integer_to_int(&K->M);
+  assert(p < FIELD_ORDER_LIMIT);
 
   lp_upolynomial_factors_t* factors = lp_upolynomial_factors_construct();
 
@@ -360,11 +366,13 @@ lp_upolynomial_factors_t* upolynomial_factor_distinct_degree(const lp_upolynomia
   return factors;
 }
 
-
-static void Q_construct(lp_integer_t* Q, size_t size, const lp_upolynomial_t* u) {
+static
+void Q_construct(lp_integer_t* Q, size_t size, const lp_upolynomial_t* u) {
 
   const lp_int_ring_t* K = lp_upolynomial_ring(u);
+  assert(mpz_fits_slong_p(&K->M));
   size_t p = (size_t) integer_to_int(&K->M);
+  assert(p < FIELD_ORDER_LIMIT);
 
   size_t k;
 
@@ -393,7 +401,8 @@ static void Q_construct(lp_integer_t* Q, size_t size, const lp_upolynomial_t* u)
   integer_destruct(&tmp); integer_destruct(&one);
 }
 
-static void Q_column_multiply(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, int j, const lp_integer_t* m) {
+static
+void Q_column_multiply(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, int j, const lp_integer_t* m) {
 
   if (trace_is_enabled("nullspace")) {
     tracef("Multiplying column %d of Q with ", j); integer_print(m, trace_out); tracef("\n");
@@ -417,7 +426,8 @@ static void Q_column_multiply(const lp_int_ring_t* K, lp_integer_t* Q, size_t si
 }
 
 // add column j into i multiplied by m
-static void Q_column_add(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, int i, const lp_integer_t* m, int j) {
+static
+void Q_column_add(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, int i, const lp_integer_t* m, int j) {
 
   assert(i != j);
 
@@ -453,7 +463,8 @@ static void Q_column_add(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, i
   }
 }
 
-static void Q_null_space(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, lp_integer_t** v, size_t* v_size) {
+static
+void Q_null_space(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, lp_integer_t** v, size_t* v_size) {
 
   *v_size = 0;
 
@@ -520,9 +531,9 @@ static void Q_null_space(const lp_int_ring_t* K, lp_integer_t* Q, size_t size, l
 
 }
 
-static void Q_destruct(lp_integer_t* Q, size_t size) {
-  size_t i;
-  for (i = 0; i < size*size; ++ i) {
+static
+void Q_destruct(lp_integer_t* Q, size_t size) {
+  for (size_t i = 0; i < size*size; ++ i) {
     integer_destruct(Q + i);
   }
 }
@@ -584,6 +595,7 @@ static void Q_destruct(lp_integer_t* Q, size_t size) {
  * does not filter the factorization, so we start at 2. If the size of the
  * basis is 1, then the polynomial is already irreducible.
  */
+static
 lp_upolynomial_factors_t* upolynomial_factor_berlekamp_square_free(const lp_upolynomial_t* f) {
 
   if (trace_is_enabled("berlekamp")) {
@@ -603,7 +615,7 @@ lp_upolynomial_factors_t* upolynomial_factor_berlekamp_square_free(const lp_upol
   } else {
 
     // The prime (should be small)
-    assert(integer_cmp_int(lp_Z, &K->M, 100) < 0);
+    assert(integer_cmp_int(lp_Z, &K->M, FIELD_ORDER_LIMIT_BERLEKAMP) < 0);
     int p = integer_to_int(&K->M);
 
     // Construct the Q matrix
@@ -630,16 +642,14 @@ lp_upolynomial_factors_t* upolynomial_factor_berlekamp_square_free(const lp_upol
     lp_upolynomial_factors_add(factors, lp_upolynomial_construct_copy(f), 1);
 
     // Filter through the null-basis
-    size_t i;
     int done = 0;
-    for (i = 1; !done && i < v_size; ++i) {
+    for (size_t i = 1; !done && i < v_size; ++i) {
 
       // v-s polynomial
       lp_upolynomial_t* v_poly[p];
 
       // Construct the v[i](x) - s polynomials
-      int s;
-      for (s = 0; s < p; ++ s) {
+      for (int s = 0; s < p; ++ s) {
         v_poly[s] = lp_upolynomial_construct(K, deg_f - 1, v[i]);
         integer_dec(K, &v[i][0]);
       }
@@ -656,8 +666,7 @@ lp_upolynomial_factors_t* upolynomial_factor_berlekamp_square_free(const lp_upol
         }
 
         // Filter the current factorization with v_i
-        int s;
-        for (s = 0; !done && s < p; ++ s) {
+        for (int s = 0; !done && s < p; ++ s) {
           // Compute the gcd
           lp_upolynomial_t* gcd = lp_upolynomial_gcd(factor, v_poly[s]);
           // We only split if gcd splits factor
@@ -677,15 +686,14 @@ lp_upolynomial_factors_t* upolynomial_factor_berlekamp_square_free(const lp_upol
       }
 
       // Free  the v[i](x) -s polynomials
-      for (s = 0; s < p; ++ s) {
+      for (int s = 0; s < p; ++ s) {
         lp_upolynomial_delete(v_poly[s]);
       }
     }
 
     // Deallocate the null vectors
-    size_t j;
-    for (i = 0; i < v_size; ++i) {
-      for (j = 0; j < deg_f; ++j) {
+    for (size_t i = 0; i < v_size; ++i) {
+      for (size_t j = 0; j < deg_f; ++j) {
         integer_destruct(v[i] + j);
       }
       free(v[i]);
@@ -717,7 +725,7 @@ lp_upolynomial_factors_t* upolynomial_factor_berlekamp_square_free(const lp_upol
  */
 lp_upolynomial_factors_t* upolynomial_factor_Zp(const lp_upolynomial_t* f) {
 
-  if (trace_is_enabled("berlekamp")) {
+  if (trace_is_enabled("factorization")) {
     tracef("upolynomial_factor_Zp("); lp_upolynomial_print(f, trace_out); tracef(")\n");
   }
 
@@ -749,7 +757,9 @@ lp_upolynomial_factors_t* upolynomial_factor_Zp(const lp_upolynomial_t* f) {
     size_t f_i_multiplicity = sq_free_factors->multiplicities[i];
 
     // Extract linear factors
-    int x_int, p = integer_to_int(&K->M);
+    assert(mpz_fits_slong_p(&K->M));
+    long x_int, p = integer_to_int(&K->M);
+    assert(p < FIELD_ORDER_LIMIT);
     lp_upolynomial_t* linear_factors_product = 0;
     for (x_int = 0; x_int < p; ++ x_int) {
       // Values
@@ -813,7 +823,6 @@ lp_upolynomial_factors_t* upolynomial_factor_Zp(const lp_upolynomial_t* f) {
   return result;
 }
 
-
 /**
  * All polynomials here are in Zp for some prime p.
  *
@@ -844,6 +853,7 @@ lp_upolynomial_factors_t* upolynomial_factor_Zp(const lp_upolynomial_t* f) {
  *
  * [6] (P2/A1)*U2 + ... + (Pr/A1)*Ur = V1
  */
+static
 void hensel_lift_initialize(const lp_upolynomial_factors_t* A, lp_upolynomial_factors_t* U) {
 
   int i;
@@ -931,6 +941,7 @@ void hensel_lift_initialize(const lp_upolynomial_factors_t* A, lp_upolynomial_fa
  *
  * All products are in Z[x]
  */
+static
 void hensel_lift_compute_products(const lp_upolynomial_factors_t* A, lp_upolynomial_t** P_z) {
 
   int k;
@@ -1008,6 +1019,7 @@ void hensel_lift_compute_products(const lp_upolynomial_factors_t* A, lp_upolynom
  *  Same as input for (q^2, F, Bk, Vk)
  *
  */
+static
 void hensel_lift_quadratic(const lp_upolynomial_t* F,
     const lp_upolynomial_factors_t* A, const lp_upolynomial_factors_t* U,
     lp_upolynomial_factors_t* B, lp_upolynomial_factors_t* V)
@@ -1252,6 +1264,7 @@ void hensel_lift_quadratic(const lp_upolynomial_t* F,
  * the result in factors. Basically try combinations of factors and see if they
  * divide f.
  */
+static
 void factorization_recombination(const lp_upolynomial_t* f, const lp_upolynomial_factors_t* factors_p, lp_upolynomial_factors_t* factors) {
 
   // Our copy of f to work with
@@ -1428,6 +1441,7 @@ static size_t primes_count = sizeof(primes)/sizeof(long);
 /**
  * Factors a square-free f in Z.
  */
+static
 lp_upolynomial_factors_t* upolynomial_factor_Z_square_free(const lp_upolynomial_t* f) {
 
   if (trace_is_enabled("factorization")) {
@@ -1577,7 +1591,7 @@ lp_upolynomial_factors_t* upolynomial_factor_Z(const lp_upolynomial_t* f) {
   lp_upolynomial_t* f_pp = lp_upolynomial_primitive_part_Z(f);
 
   // Get a square-free decomposition of f
-  lp_upolynomial_factors_t* sq_free_factors = lp_upolynomial_factor_square_free_primitive(f_pp);
+  lp_upolynomial_factors_t* sq_free_factors = upolynomial_factor_square_free_primitive(f_pp);
   assert(integer_cmp_int(lp_Z, &sq_free_factors->constant, 1) == 0);
 
   // Factor individuals
diff --git a/src/upolynomial/factorization.h b/src/upolynomial/factorization.h
index afa992a0..0f2b0296 100644
--- a/src/upolynomial/factorization.h
+++ b/src/upolynomial/factorization.h
@@ -48,4 +48,4 @@ lp_upolynomial_factors_t* upolynomial_factor_Z(const lp_upolynomial_t* f);
  * Factors a primitive polynomial into its square-free factors. If x is a factor,
  * it is factored into a separate factor.
  */
-lp_upolynomial_factors_t* lp_upolynomial_factor_square_free_primitive(const lp_upolynomial_t* f);
+lp_upolynomial_factors_t* upolynomial_factor_square_free_primitive(const lp_upolynomial_t* f);
diff --git a/src/upolynomial/root_finding.c b/src/upolynomial/root_finding.c
index f295e9ff..c64d1b57 100644
--- a/src/upolynomial/root_finding.c
+++ b/src/upolynomial/root_finding.c
@@ -24,15 +24,18 @@
 #include "upolynomial/factors.h"
 #include "upolynomial/root_finding.h"
 #include "upolynomial/factorization.h"
+#include "upolynomial/upolynomial_vector.h"
 #include "upolynomial/upolynomial_dense.h"
 #include "upolynomial/output.h"
 
 #include <assert.h>
 #include <stdlib.h>
-#include <stdlib.h>
 
 #include "utils/debug_trace.h"
 
+/** Field order up to we're doing brute force. */
+#define FIELD_ORDER_LIMIT 1000
+
 /** Negative infinity */
 #define INF_N (void*) 0
 /** Positive infinity */
@@ -102,6 +105,7 @@ void upolynomial_compute_sturm_sequence(const lp_upolynomial_t* f, upolynomial_d
  * Compute the number sgn_changes(a) given Sturm sequence and a. If a is 0 or
  * 1 as a pointer, we evaluate at -inf, +inf respectively.
  */
+static
 int sturm_seqence_count_sign_changes(
     const upolynomial_dense_t* sturm_sequence, int sturm_sequence_size,
     const lp_rational_t* a, int max_changes)
@@ -128,6 +132,7 @@ int sturm_seqence_count_sign_changes(
  * Compute the number sgn_changes(a) given Sturm sequence and a. If a is 0 or
  * 1 as a pointer, we evaluate at -inf, +inf respectively.
  */
+static
 int sturm_seqence_count_sign_changes_dyadic(
     const upolynomial_dense_t* sturm_sequence, int sturm_sequence_size,
     const lp_dyadic_rational_t* a, int max_changes)
@@ -150,6 +155,7 @@ int sturm_seqence_count_sign_changes_dyadic(
 }
 
 /** Count roots in the given interval */
+static
 int sturm_seqence_count_roots(
     const upolynomial_dense_t* sturm_sequence, int sturm_sequence_size,
     const lp_rational_interval_t* interval)
@@ -178,6 +184,7 @@ int sturm_seqence_count_roots(
 }
 
 /** Count roots in the given interval */
+static
 int sturm_seqence_count_roots_dyadic(
     const upolynomial_dense_t* sturm_sequence, int sturm_sequence_size,
     const lp_dyadic_interval_t* interval)
@@ -263,6 +270,7 @@ int upolynomial_roots_count_sturm(const lp_upolynomial_t* f, const lp_rational_i
 /**
  * Recursive root isolation on an interval (a, b].
  */
+static
 void sturm_seqence_isolate_roots(
     const upolynomial_dense_t* S, size_t S_size,
     lp_algebraic_number_t* roots, size_t* roots_size,
@@ -432,8 +440,7 @@ void upolynomial_roots_isolate_sturm(const lp_upolynomial_t* f, lp_algebraic_num
 
     // Destroy the temporaries
     lp_dyadic_interval_destruct(&interval_all);
-    size_t i;
-    for (i = 0; i < sturm_sequence_size; ++i) {
+    for (size_t i = 0; i < sturm_sequence_size; ++i) {
       upolynomial_dense_destruct(sturm_sequence + i);
     }
     free(sturm_sequence);
@@ -451,3 +458,206 @@ void upolynomial_roots_isolate_sturm(const lp_upolynomial_t* f, lp_algebraic_num
   // Destroy the factors
   lp_upolynomial_factors_destruct(square_free_factors, 1);
 }
+
+/**
+ * helper function for rabin root finding
+ * returns x^p mod f
+ */
+static
+lp_upolynomial_t* upolynomial_power_mod(const lp_upolynomial_t *b, const lp_integer_t *e, const lp_upolynomial_t *m) {
+  lp_integer_t ee;
+  lp_integer_construct_copy(lp_Z, &ee, e);
+  lp_upolynomial_t *acc = lp_upolynomial_construct_power(m->K, 0, 1);
+  lp_upolynomial_t *power = lp_upolynomial_construct_copy(b);
+
+  while (!lp_integer_is_zero(lp_Z, &ee)) {
+    if (integer_is_odd(&ee)) {
+      upolynomial_op_inplace(lp_upolynomial_mul, &acc, power);
+      upolynomial_op_inplace(lp_upolynomial_rem_exact, &acc, m);
+    }
+    upolynomial_op_inplace(lp_upolynomial_mul, &power, power);
+    upolynomial_op_inplace(lp_upolynomial_rem_exact, &power, m);
+
+    integer_div_floor_pow2(&ee, &ee, 1);
+  }
+
+  lp_integer_destruct(&ee);
+  lp_upolynomial_delete(power);
+
+  return acc;
+}
+
+/**
+ * helper function for rabin root finding
+ * returns gcd(f, x^p - x)
+ */
+static
+lp_upolynomial_t* upolynomial_reduce_to_linear(const lp_upolynomial_t *f ) {
+  lp_upolynomial_t *x = lp_upolynomial_construct_power(f->K, 1, 1);
+  lp_upolynomial_t *pm = upolynomial_power_mod(x, &f->K->M, f);
+  lp_upolynomial_t *fp = lp_upolynomial_sub(pm, x);
+  lp_upolynomial_t *linear = lp_upolynomial_gcd(f, fp);
+
+  lp_upolynomial_delete(x);
+  lp_upolynomial_delete(pm);
+  lp_upolynomial_delete(fp);
+
+  return linear;
+}
+
+/**
+ * Implements Rabin root-finding.
+ * Rabin, Michael O. "Probabilistic algorithms in finite fields." SIAM Journal on computing 9.2 (1980): 273-280.
+ */
+static
+void upolynomial_roots_find_rabin(const lp_upolynomial_t* f, lp_integer_t* roots, size_t* roots_size) {
+  const lp_int_ring_t* K = f->K;
+
+  *roots_size = 0;
+
+  gmp_randstate_t state;
+  gmp_randinit_default(state);
+
+  lp_integer_t s;
+  lp_integer_construct_copy(lp_Z, &s, &K->M);
+  // right shift by p1
+  integer_div_floor_pow2(&s, &s, 1);
+
+  // construct helper polynomials (1 and x)
+  lp_upolynomial_t *p1 = lp_upolynomial_construct_power(K, 0, 1);
+  lp_upolynomial_t *px = lp_upolynomial_construct_power(K, 1, 1);
+
+  lp_upolynomial_vector_t *to_factor = lp_upolynomial_vector_construct();
+  lp_upolynomial_vector_move_back(to_factor, upolynomial_reduce_to_linear(f));
+
+  while (lp_upolynomial_vector_size(to_factor)) {
+    lp_upolynomial_t *p = lp_upolynomial_vector_pop(to_factor);
+    size_t d = lp_upolynomial_degree(p);
+    bool const_is_zero = !lp_upolynomial_const_term(p) || lp_integer_is_zero(lp_Z, lp_upolynomial_const_term(p));
+    if (d == 0) {
+      // it's dead
+    } else if (const_is_zero) {
+      // it has a zero root
+      assert(lp_upolynomial_degree(f) > *roots_size);
+      lp_integer_construct_from_int(lp_Z, &roots[(*roots_size)++], 0);
+      lp_upolynomial_vector_move_back(to_factor, lp_upolynomial_div_exact(p, px));
+    } else if (d == 1) {
+      // push the zero
+      lp_upolynomial_make_monic_in_place(p);
+      lp_upolynomial_neg_in_place(p);
+      assert(lp_upolynomial_const_term(p));
+      assert(lp_upolynomial_degree(f) > *roots_size);
+      lp_integer_construct_copy(lp_Z, &roots[(*roots_size)++], lp_upolynomial_const_term(p));
+    } else {
+      // gen (px - delta)
+      lp_upolynomial_t *xd = lp_upolynomial_construct_from_int(K, 1, (int[2]){1, 1});
+      lp_integer_t *delta = (lp_integer_t*) lp_upolynomial_const_term(xd);
+      assert(delta);
+
+      bool done = false;
+      while (!done) {
+        // set constant term delta in xd to random value
+        mpz_urandomm(delta, state, &K->M);
+        integer_ring_normalize(K, delta);
+        integer_neg(K, delta, delta);
+        // calculate gcd
+        lp_upolynomial_t *h = upolynomial_power_mod(xd, &s, p);
+        upolynomial_op_inplace(lp_upolynomial_sub, &h, p1);
+        upolynomial_op_inplace(lp_upolynomial_gcd, &h, p);
+        if (0 < lp_upolynomial_degree(h) && lp_upolynomial_degree(h) < d) {
+          lp_upolynomial_vector_push_back(to_factor, h);
+          lp_upolynomial_vector_move_back(to_factor, lp_upolynomial_div_exact(p, h));
+          done = true;
+        }
+        lp_upolynomial_delete(h);
+      }
+      lp_upolynomial_delete(xd);
+    }
+    lp_upolynomial_delete(p);
+  }
+
+  lp_integer_destruct(&s);
+  lp_upolynomial_delete(p1);
+  lp_upolynomial_delete(px);
+  lp_upolynomial_vector_delete(to_factor);
+  gmp_randclear(state);
+}
+
+/**
+ * iterates over all p possible values of Zp and check each
+ */
+static
+void upolynomial_roots_find_brute_force(const lp_upolynomial_t* f, lp_integer_t* roots, size_t* roots_size) {
+  const lp_int_ring_t* K = f->K;
+
+  size_t d = lp_upolynomial_degree(f);
+  assert(d > 0);
+
+  assert(mpz_fits_slong_p(&K->M));
+  long p = integer_to_int(&K->M);
+  assert(p < FIELD_ORDER_LIMIT);
+
+  *roots_size = 0;
+
+  // Values
+  lp_integer_t x, value;
+  integer_construct_copy(lp_Z, &x, &K->lb);
+  integer_construct(&value);
+
+  for (long x_int = 0; x_int < p; ++x_int) {
+    assert(integer_in_ring(K, &x));
+
+    // Evaluate the polynomial
+    lp_upolynomial_evaluate_at_integer(f, &x, &value);
+
+    // If zero we found one
+    if (integer_sgn(K, &value) == 0) {
+      lp_integer_construct_copy(K, roots + *roots_size, &x);
+      (*roots_size) ++;
+
+      // check if we found all
+      if (*roots_size == d) {
+        break;
+      }
+    }
+
+    // next element
+    integer_inc(lp_Z, &x);
+  }
+  assert(*roots_size == d || lp_integer_cmp(lp_Z, &x, &K->ub));
+
+  integer_destruct(&value);
+  integer_destruct(&x);
+}
+
+void upolynomial_roots_find_Zp(const lp_upolynomial_t* f, lp_integer_t** roots, size_t* roots_size) {
+  if (trace_is_enabled("roots")) {
+    tracef("upolynomial_roots_find_Zp("); lp_upolynomial_print(f, trace_out); tracef(")\n");
+  }
+
+  const lp_int_ring_t* K = f->K;
+  assert(K && K->is_prime);
+
+  size_t d = lp_upolynomial_degree(f);
+
+  *roots = malloc(sizeof(lp_integer_t) * d);
+
+  // depending on the finite field size, choose appropriate function
+  if (integer_cmp_int(lp_Z, &K->M, FIELD_ORDER_LIMIT) < 0) {
+    upolynomial_roots_find_brute_force(f, *roots, roots_size);
+  } else {
+    upolynomial_roots_find_rabin(f, *roots, roots_size);
+  }
+
+  // realloc has undefined behavior in case size is zero
+  if (*roots_size == 0) {
+    free(*roots);
+    *roots = NULL;
+  } else if (*roots_size < d) {
+    *roots = realloc(*roots, sizeof(lp_integer_t) * (*roots_size));
+  }
+
+  if (trace_is_enabled("roots")) {
+    tracef("upolynomial_roots_find_Zp() found %zu roots\n", *roots_size);
+  }
+}
diff --git a/src/upolynomial/root_finding.h b/src/upolynomial/root_finding.h
index 00fc9c97..ff2a3386 100644
--- a/src/upolynomial/root_finding.h
+++ b/src/upolynomial/root_finding.h
@@ -52,3 +52,9 @@ int upolynomial_roots_count_sturm(const lp_upolynomial_t* f, const lp_rational_i
  * will be updated to the number of roots.
  */
 void upolynomial_roots_isolate_sturm(const lp_upolynomial_t* f, lp_algebraic_number_t* roots, size_t* roots_size);
+
+/**
+* Finds the roots for a polynomial over Zp. Uses brute force or rabin root
+* 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);
diff --git a/src/upolynomial/upolynomial.c b/src/upolynomial/upolynomial.c
index cbf048ed..c01f2b3a 100644
--- a/src/upolynomial/upolynomial.c
+++ b/src/upolynomial/upolynomial.c
@@ -251,7 +251,7 @@ const lp_integer_t* lp_upolynomial_const_term(const lp_upolynomial_t* p) {
   if (p->monomials[0].degree == 0) {
     return &p->monomials[0].coefficient;
   } else {
-    return 0;
+    return NULL;
   }
 }
 
@@ -774,7 +774,7 @@ void lp_upolynomial_div_rem_exact(const lp_upolynomial_t* p, const lp_upolynomia
     upolynomial_dense_t rem_buffer;
     upolynomial_dense_t div_buffer;
     lp_upolynomial_div_general(p, q, &div_buffer, &rem_buffer, /** exact */ 1);
-    *div= upolynomial_dense_to_upolynomial(&div_buffer, K);
+    *div = upolynomial_dense_to_upolynomial(&div_buffer, K);
     *rem = upolynomial_dense_to_upolynomial(&rem_buffer, K);
     upolynomial_dense_destruct(&div_buffer);
     upolynomial_dense_destruct(&rem_buffer);
@@ -1091,6 +1091,7 @@ void lp_upolynomial_solve_bezout(const lp_upolynomial_t* p, const lp_upolynomial
 }
 
 lp_upolynomial_factors_t* lp_upolynomial_factor(const lp_upolynomial_t* p) {
+  assert(p);
 
   if (trace_is_enabled("factorization")) {
     tracef("upolynomial_factor("); lp_upolynomial_print(p, trace_out); tracef(")\n");
@@ -1156,6 +1157,15 @@ void lp_upolynomial_sturm_sequence(const lp_upolynomial_t* f, lp_upolynomial_t**
   free(S_dense);
 }
 
+void lp_upolynomial_roots_find_Zp(const lp_upolynomial_t* f, lp_integer_t** roots, size_t* roots_size) {
+  if (trace_is_enabled("roots")) {
+    tracef("upolynomial_roots_find_Zp("); lp_upolynomial_print(f, trace_out); tracef(")\n");
+  }
+  assert(f->K != lp_Z);
+
+  upolynomial_roots_find_Zp(f, roots, roots_size);
+}
+
 lp_upolynomial_t* lp_upolynomial_neg(const lp_upolynomial_t* p) {
   lp_upolynomial_t* neg = lp_upolynomial_construct_copy(p);
   lp_upolynomial_neg_in_place(neg);
@@ -1169,6 +1179,27 @@ void lp_upolynomial_neg_in_place(lp_upolynomial_t* p) {
   }
 }
 
+lp_upolynomial_t* lp_upolynomial_make_monic(const lp_upolynomial_t* p) {
+  lp_upolynomial_t* monic = lp_upolynomial_construct_copy(p);
+  lp_upolynomial_make_monic_in_place(monic);
+  return monic;
+}
+
+void lp_upolynomial_make_monic_in_place(lp_upolynomial_t* p) {
+  if (lp_upolynomial_is_zero(p)) {
+    return;
+  }
+  const lp_int_ring_t* K = p->K;
+
+  const lp_integer_t *c = lp_upolynomial_lead_coeff(p);
+  assert(integer_cmp_int(K, c, 0));
+  for (size_t i = 0; i < p->size; ++i) {
+    assert(integer_divides(K, c, &p->monomials[i].coefficient));
+    integer_div_exact(K, &p->monomials[i].coefficient, &p->monomials[i].coefficient, c);
+  }
+  assert(lp_upolynomial_is_monic(p));
+}
+
 void lp_upolynomial_reverse_in_place(lp_upolynomial_t* p) {
   size_t p_deg;
   ulp_monomial_t tmp, *m_i, *m_j;
diff --git a/src/upolynomial/upolynomial.h b/src/upolynomial/upolynomial.h
index db2c62d0..0eeb9d46 100644
--- a/src/upolynomial/upolynomial.h
+++ b/src/upolynomial/upolynomial.h
@@ -20,6 +20,7 @@
 #pragma once
 
 #include <integer.h>
+#include <upolynomial.h>
 #include "upolynomial/umonomial.h"
 
 /**
@@ -33,3 +34,12 @@ struct lp_upolynomial_struct {
   /** The monomials */
   ulp_monomial_t monomials[];
 };
+
+typedef lp_upolynomial_t* (*upolynomial_op)(const lp_upolynomial_t*, const lp_upolynomial_t*);
+
+static inline
+void upolynomial_op_inplace(upolynomial_op op, lp_upolynomial_t **a, const lp_upolynomial_t *b) {
+  lp_upolynomial_t *r = op(*a, b);
+  lp_upolynomial_delete(*a);
+  *a = r;
+}
diff --git a/src/upolynomial/upolynomial_vector.c b/src/upolynomial/upolynomial_vector.c
new file mode 100644
index 00000000..715305eb
--- /dev/null
+++ b/src/upolynomial/upolynomial_vector.c
@@ -0,0 +1,82 @@
+/**
+ * Copyright 2015, SRI International.
+ *
+ * This file is part of LibPoly.
+ *
+ * LibPoly is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LibPoly is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with LibPoly.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <stdlib.h>
+
+#include "upolynomial_vector.h"
+
+#include <upolynomial.h>
+#include "upolynomial/upolynomial.h"
+
+lp_upolynomial_vector_t *lp_upolynomial_vector_construct(void) {
+  lp_upolynomial_vector_t*v = malloc(sizeof(lp_upolynomial_vector_t));
+  v->size = 0;
+  v->capacity = 10;
+  v->data = malloc(v->capacity * sizeof(lp_upolynomial_t*));
+  return v;
+}
+
+void lp_upolynomial_vector_swap(lp_upolynomial_vector_t *v1, lp_upolynomial_vector_t *v2) {
+  lp_upolynomial_vector_t tmp = *v1;
+  *v1 = *v2;
+  *v2 = tmp;
+}
+
+void lp_upolynomial_vector_clear(lp_upolynomial_vector_t *v) {
+  for (size_t i = 0; i < v->size; ++i) {
+    lp_upolynomial_delete(v->data[i]);
+  }
+  v->size = 0;
+}
+
+void lp_upolynomial_vector_delete(lp_upolynomial_vector_t *v) {
+  lp_upolynomial_vector_clear(v);
+  free(v->data);
+  free(v);
+}
+
+size_t lp_upolynomial_vector_size(const lp_upolynomial_vector_t *v) {
+  return v->size;
+}
+
+lp_upolynomial_t *lp_upolynomial_vector_at(lp_upolynomial_vector_t *v, size_t i) {
+  return lp_upolynomial_construct_copy(v->data[i]);
+}
+
+void lp_upolynomial_vector_move_back(lp_upolynomial_vector_t *v, lp_upolynomial_t *p) {
+  if (v->size == v->capacity) {
+    v->capacity *= 2;
+    v->data = realloc(v->data, v->capacity * sizeof(lp_upolynomial_t*));
+  }
+  v->data[v->size] = p;
+  v->size ++;
+}
+
+void lp_upolynomial_vector_push_back(lp_upolynomial_vector_t *v, const lp_upolynomial_t *p) {
+  lp_upolynomial_vector_move_back(v, lp_upolynomial_construct_copy(p));
+}
+
+lp_upolynomial_t *lp_upolynomial_vector_pop(lp_upolynomial_vector_t *v) {
+  if (v->size == 0) {
+    return NULL;
+  } else {
+    v->size --;
+    return v->data[v->size];
+  }
+}
diff --git a/src/upolynomial/upolynomial_vector.h b/src/upolynomial/upolynomial_vector.h
new file mode 100644
index 00000000..7862e09d
--- /dev/null
+++ b/src/upolynomial/upolynomial_vector.h
@@ -0,0 +1,78 @@
+/**
+ * Copyright 2015, SRI International.
+ *
+ * This file is part of LibPoly.
+ *
+ * LibPoly is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LibPoly is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with LibPoly.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <poly.h>
+
+struct lp_upolynomial_vector_struct {
+  /** Number of elements in the vector */
+  size_t size;
+  /** Size of the array */
+  size_t capacity;
+  /** The elements */
+  lp_upolynomial_t**data;
+};
+
+typedef struct lp_upolynomial_vector_struct lp_upolynomial_vector_t;
+
+/**
+ * Construct a new upolynomial vector
+ */
+lp_upolynomial_vector_t* lp_upolynomial_vector_construct(void);
+
+/**
+ * Swaps two polynomial vectors
+ */
+void lp_upolynomial_vector_swap(lp_upolynomial_vector_t *v1, lp_upolynomial_vector_t *v2);
+
+/**
+ * Clears the upolynomial vector
+ */
+void lp_upolynomial_vector_clear(lp_upolynomial_vector_t *v);
+
+/**
+ * Deletes the vector.
+ */
+void lp_upolynomial_vector_delete(lp_upolynomial_vector_t *v);
+
+/**
+ * Returns the size of the vector.
+ */
+size_t lp_upolynomial_vector_size(const lp_upolynomial_vector_t *v);
+
+/**
+ * Returns a copy of the polynomial at position i
+ */
+lp_upolynomial_t* lp_upolynomial_vector_at(lp_upolynomial_vector_t *v, size_t i);
+
+/**
+ * Moves a polynomial to the vector. p must not be used afterwards.
+ */
+void lp_upolynomial_vector_move_back(lp_upolynomial_vector_t *v, lp_upolynomial_t* p);
+
+/**
+ * Makes a copy of p and adds it to the vector.
+ */
+void lp_upolynomial_vector_push_back(lp_upolynomial_vector_t *v, const lp_upolynomial_t* p);
+
+/**
+ * Removes and returns the last element of the vector.
+ */
+lp_upolynomial_t* lp_upolynomial_vector_pop(lp_upolynomial_vector_t *v);
diff --git a/test/polyxx/test_upolynomial.cpp b/test/polyxx/test_upolynomial.cpp
index 5b6d7c52..6d561524 100644
--- a/test/polyxx/test_upolynomial.cpp
+++ b/test/polyxx/test_upolynomial.cpp
@@ -274,3 +274,45 @@ TEST_CASE("upolynomial::operator<<") {
   out << p;
   CHECK(out.str() == "5*x^4 + 4*x^3 + 3*x^2 + 2*x + 1");
 }
+
+TEST_CASE("upolynomial::find_roots_Zp") {
+  // test brute force
+  {
+    Integer mod(13);
+    IntegerRing ring(mod, true);
+    UPolynomial p(ring, {1, 2, 3, 4, 5});
+    std::vector<Integer> roots = find_roots_Zp(p);
+    CHECK(roots.size() == 1);
+    CHECK(compare(ring, roots[0], -5) == 0);
+    CHECK(compare(ring, roots[0], 8) == 0);
+  }
+  {
+    Integer mod(211);
+    IntegerRing ring(mod, true);
+    UPolynomial p(ring, {1, 2, 3, 4, 5});
+    std::vector<Integer> roots = find_roots_Zp(p);
+    CHECK(roots.size() == 2);
+    CHECK(compare(ring, roots[0], -31) == 0);
+    CHECK(compare(ring, roots[0], 180) == 0);
+    CHECK(compare(ring, roots[1], 51) == 0);
+  }
+  // test rabin root finding
+  {
+    Integer mod(1000003);
+    IntegerRing ring(mod, true);
+    UPolynomial p(ring, {1, 2, 3, 4, 5});
+    std::vector<Integer> roots = find_roots_Zp(p);
+    CHECK(roots.size() == 2);
+    CHECK(compare(ring, roots[0], 682693) == 0);
+    CHECK(compare(ring, roots[1], 939713) == 0);
+  }
+  {
+    Integer mod("2425967623052370772757633156976982469681", 10);
+    IntegerRing ring(mod, true);
+    UPolynomial p(ring, {1, 2, 3, 4, 5});
+    std::vector<Integer> roots = find_roots_Zp(p);
+    CHECK(roots.size() == 2);
+    CHECK(compare(ring, roots[0], Integer("-1812258784923425884426588992727582130833", 10)) == 0);
+    CHECK(compare(ring, roots[1], Integer("-2083615688050021089599246736669243933656", 10)) == 0);
+  }
+}
diff --git a/test/python/tests/upolynomial_roots.py b/test/python/tests/upolynomial_roots.py
index a4db4867..fc475f1c 100755
--- a/test/python/tests/upolynomial_roots.py
+++ b/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])