diff --git a/include/polynomial.h b/include/polynomial.h
index 58604ec..36f1a9f 100644
--- a/include/polynomial.h
+++ b/include/polynomial.h
@@ -109,6 +109,12 @@ 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);
 
+/**
+ * Sets the context of a polynomial.
+ * Warning: variable DB indexes are not checked. Use with caution.
+ */
+void lp_polynomial_set_context(lp_polynomial_t* A, const lp_polynomial_context_t* ctx);
+
 /** Returns all the variables (it will not clear the output list vars) */
 void lp_polynomial_get_variables(const lp_polynomial_t* A, lp_variable_list_t* vars);
 
@@ -375,6 +381,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);
 
+/**
+ * Substitutes the polynomial with M and checks the sign of the resulting real number.
+ * If M does not fully evaluate A to a real, -1 is returned.
+ */
+int lp_polynomial_constraint_evaluate_subs(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.
diff --git a/include/polyxx/polynomial.h b/include/polyxx/polynomial.h
index d18ee57..3018bd5 100644
--- a/include/polyxx/polynomial.h
+++ b/include/polyxx/polynomial.h
@@ -117,8 +117,8 @@ namespace poly {
    * 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. */
+   * 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);
@@ -127,6 +127,11 @@ namespace poly {
   /** Evaluate a polynomial constraint over an assignment. */
   bool evaluate_constraint(const Polynomial& p, const Assignment& a,
                            SignCondition sc);
+  /** Checks if a polynomial constraint is fully evaluated over assignment.
+   *  Returns -1 if 'a' does not fully evaluate the constraint.
+   *  Otherwise, returns true if the constraint holds under 'a'. */
+  int evaluate_constraint_subs(const Polynomial& p, const Assignment& a,
+                               SignCondition sc);
   /** Evaluate a polynomial over an interval assignment. The result is only an
    * approximation. */
   Interval evaluate(const Polynomial& p, const IntervalAssignment& a);
diff --git a/include/polyxx/value.h b/include/polyxx/value.h
index 4534b51..1737684 100644
--- a/include/polyxx/value.h
+++ b/include/polyxx/value.h
@@ -29,6 +29,8 @@ namespace poly {
     /** Construct a none value. */
     Value();
     /** Construct from a native integer. */
+    Value(int i);
+    /** Construct from a native integer. */
     Value(long i);
     /** Copy from the given Value. */
     Value(const Value& val);
diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c
index 4e65889..497c914 100644
--- a/src/polynomial/polynomial.c
+++ b/src/polynomial/polynomial.c
@@ -1177,6 +1177,14 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t
     lp_value_construct_none(&x_value_backup);
   }
 
+  lp_polynomial_t B;
+  lp_polynomial_construct(&B, A->ctx);
+
+  lp_integer_t multiplier;
+  integer_construct(&multiplier);
+  coefficient_evaluate_rationals(A->ctx, &A->data, M, &B.data, &multiplier);
+  integer_destruct(&multiplier);
+
   size_t i;
 
   lp_polynomial_t** factors = 0;
@@ -1190,11 +1198,12 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t
   // Get the reduced polynomial
   lp_polynomial_t A_r;
   lp_polynomial_construct(&A_r, A->ctx);
-  lp_polynomial_reductum_m(&A_r, A, M);
-  assert(x == lp_polynomial_top_variable(A));
-
-  // Get the square-free factorization
-  lp_polynomial_factor_square_free(&A_r, &factors, &multiplicities, &factors_size);
+  if (x == lp_polynomial_top_variable(&B)) {
+    lp_polynomial_reductum_m(&A_r, &B, M);
+    assert(x == lp_polynomial_top_variable(&B));
+    // Get the square-free factorization
+    lp_polynomial_factor_square_free(&A_r, &factors, &multiplicities, &factors_size);
+  }
 
   // Count the max number of roots
   size_t total_degree = 0;
@@ -1297,6 +1306,7 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t
   free(roots_tmp);
   lp_value_destruct(&x_value_backup);
   lp_polynomial_destruct(&A_r);
+  lp_polynomial_destruct(&B);
 }
 
 lp_feasibility_set_t* lp_polynomial_constraint_get_feasible_set(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, int negated, const lp_assignment_t* M) {
@@ -2134,6 +2144,29 @@ int lp_polynomial_constraint_evaluate(const lp_polynomial_t* A, lp_sign_conditio
   return lp_sign_condition_consistent(sgn_condition, p_sign);
 }
 
+int lp_polynomial_constraint_evaluate_subs(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, const lp_assignment_t* M) {
+  coefficient_t A_rat;
+  lp_integer_t multiplier;
+
+  lp_polynomial_external_clean(A);
+  assert(A->ctx->K == lp_Z);
+
+  integer_construct(&multiplier);
+  coefficient_construct(A->ctx, &A_rat);
+  coefficient_evaluate_rationals(A->ctx, &A->data, M, &A_rat, &multiplier);
+  integer_destruct(&multiplier);
+
+  int res = -1;
+  if (A_rat.type == COEFFICIENT_NUMERIC) {
+    int sgn = integer_sgn(lp_Z, &A_rat.value.num);
+    res = lp_sign_condition_consistent(sgn_condition, sgn);
+    assert(res >= 0);
+  }
+
+  coefficient_destruct(&A_rat);
+  return res;
+}
+
 int lp_polynomial_constraint_evaluate_Zp(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, const lp_assignment_t* m) {
 
   lp_polynomial_external_clean(A);
diff --git a/src/polyxx/polynomial.cpp b/src/polyxx/polynomial.cpp
index 3dc4e40..3c7f8c9 100644
--- a/src/polyxx/polynomial.cpp
+++ b/src/polyxx/polynomial.cpp
@@ -165,6 +165,11 @@ namespace poly {
     return lp_polynomial_constraint_evaluate(
         p.get_internal(), to_sign_condition(sc), a.get_internal());
   }
+  int evaluate_constraint_subs(const Polynomial& p, const Assignment& a,
+                               SignCondition sc) {
+    return lp_polynomial_constraint_evaluate_subs(
+      p.get_internal(), to_sign_condition(sc), a.get_internal());
+  }
   Interval evaluate(const Polynomial& p, const IntervalAssignment& a) {
     Interval res;
     lp_polynomial_interval_value(p.get_internal(), a.get_internal(),
diff --git a/src/polyxx/value.cpp b/src/polyxx/value.cpp
index a4ac1d4..e56a71e 100644
--- a/src/polyxx/value.cpp
+++ b/src/polyxx/value.cpp
@@ -13,6 +13,7 @@ namespace poly {
   }
 
   Value::Value() { lp_value_construct_none(&mValue); }
+  Value::Value(int i) { lp_value_construct_int(get_internal(), i); }
   Value::Value(long i) { lp_value_construct_int(get_internal(), i); }
   Value::Value(const Value& val) {
     lp_value_construct_copy(get_internal(), val.get_internal());
diff --git a/src/variable/variable_db.c b/src/variable/variable_db.c
index eeb8232..8c942f3 100644
--- a/src/variable/variable_db.c
+++ b/src/variable/variable_db.c
@@ -63,9 +63,11 @@ void lp_variable_db_add_variable(lp_variable_db_t* var_db, lp_variable_t var, co
   assert(var_db);
   while (var >= var_db->capacity) {
     lp_variable_db_resize(var_db, 2*var_db->capacity);
+    var_db->size = var_db->capacity < var ? var_db->capacity : var;
   }
   assert(var_db->variable_names[var] == 0);
   var_db->variable_names[var] = strdup(name);
+  var_db->size = var_db->size < var ? var : var_db->size;
 }
 
 void lp_variable_db_construct(lp_variable_db_t* var_db) {
diff --git a/test/polyxx/test_polynomial.cpp b/test/polyxx/test_polynomial.cpp
index b116f44..c7aec51 100644
--- a/test/polyxx/test_polynomial.cpp
+++ b/test/polyxx/test_polynomial.cpp
@@ -100,3 +100,49 @@ TEST_CASE("polynomial::constructor") {
   CHECK(p_cp == p_cp_a);
   CHECK(ptr == p_mv_a.get_internal());
 }
+
+TEST_CASE("polynomial::evaluate") {
+  Variable x("x");
+  Variable y("y");
+  Polynomial p1 = 1 * pow(x, 6) + 2 * pow(x, 5) + 3 * y - 1;
+  Polynomial p2 = 1 * y * pow(x + 1, 4);
+
+  // full assignment
+  Assignment a;
+  a.set(x, -1);
+  a.set(y, 2);
+  CHECK(a.has(x));
+  CHECK(a.has(y));
+
+  CHECK(evaluate(p1, a) == Value(4));
+  CHECK(evaluate_constraint(p1, a, SignCondition::GE) == true);
+  CHECK(evaluate_constraint(p1, a, SignCondition::LT) == false);
+  CHECK(evaluate_constraint(p1, a, SignCondition::EQ) == false);
+  CHECK(evaluate_constraint_subs(p1, a, SignCondition::GE) == 1);
+  CHECK(evaluate_constraint_subs(p1, a, SignCondition::LT) == 0);
+  CHECK(evaluate_constraint_subs(p1, a, SignCondition::EQ) == 0);
+
+  CHECK(evaluate(p2, a) == Value(0));
+  CHECK(evaluate_constraint(p2, a, SignCondition::GE) == true);
+  CHECK(evaluate_constraint(p2, a, SignCondition::LT) == false);
+  CHECK(evaluate_constraint(p2, a, SignCondition::EQ) == true);
+  CHECK(evaluate_constraint_subs(p2, a, SignCondition::GE) == 1);
+  CHECK(evaluate_constraint_subs(p2, a, SignCondition::LT) == 0);
+  CHECK(evaluate_constraint_subs(p2, a, SignCondition::EQ) == 1);
+
+  // partial assignment
+  a.unset(y);
+  CHECK(a.has(x));
+  CHECK(!a.has(y));
+
+  CHECK(evaluate_constraint_subs(p1, a, SignCondition::GE) == -1);
+  CHECK(evaluate_constraint_subs(p1, a, SignCondition::LT) == -1);
+  CHECK(evaluate_constraint_subs(p1, a, SignCondition::EQ) == -1);
+
+  CHECK(evaluate_constraint(p2, a, SignCondition::GE) == true);
+  CHECK(evaluate_constraint(p2, a, SignCondition::LT) == false);
+  CHECK(evaluate_constraint(p2, a, SignCondition::EQ) == true);
+  CHECK(evaluate_constraint_subs(p2, a, SignCondition::GE) == 1);
+  CHECK(evaluate_constraint_subs(p2, a, SignCondition::LT) == 0);
+  CHECK(evaluate_constraint_subs(p2, a, SignCondition::EQ) == 1);
+}