From a79b13b6979d10c89b525e237a4bd71da643a91a Mon Sep 17 00:00:00 2001
From: Thomas Hader <Ovascos@users.noreply.github.com>
Date: Wed, 9 Oct 2024 17:37:02 +0200
Subject: [PATCH] Integer handling for feasible sets (#83)

* code readability update

* Code cleanup

* Implemented has_integer for interval and FeasibleSet

* Created lp_integer_fits_int

* Created lp_feasibility_set_count_int

* Implement binary search for feasibility_set::contains

* fixed whitespaces

* Added comment

---------

Co-authored-by: Thomas Hader <tom@hader.email>
---
 include/feasibility_set.h        | 15 ++++++
 include/integer.h                |  5 ++
 include/interval.h               | 13 ++++-
 include/polyxx/interval.h        |  4 ++
 python/polypyFeasibilitySet2.c   | 45 +++++++++++++++-
 python/polypyFeasibilitySet3.c   | 51 +++++++++++++++---
 python/polypyInterval2.c         | 14 ++++-
 python/polypyInterval3.c         | 17 ++++--
 src/interval/interval.c          | 92 ++++++++++++++++++++++++++++----
 src/number/integer.c             |  4 ++
 src/number/integer.h             |  5 ++
 src/number/value.c               |  3 +-
 src/polynomial/feasibility_set.c | 59 ++++++++++++++++++--
 src/polyxx/interval.cpp          |  9 ++++
 src/upolynomial/factorization.c  | 10 ++--
 src/upolynomial/root_finding.c   |  2 +-
 test/polyxx/test_interval.cpp    | 26 +++++++++
 17 files changed, 335 insertions(+), 39 deletions(-)

diff --git a/include/feasibility_set.h b/include/feasibility_set.h
index c0d9d765..6afa3435 100644
--- a/include/feasibility_set.h
+++ b/include/feasibility_set.h
@@ -93,11 +93,26 @@ int lp_feasibility_set_is_full(const lp_feasibility_set_t* set);
  */
 int lp_feasibility_set_is_point(const lp_feasibility_set_t* set);
 
+/**
+ * Check if the set contains only one integer value.
+ */
+int lp_feasibility_set_is_point_int(const lp_feasibility_set_t* set);
+
+/**
+ * Returns the number of integers in set (up to LONG_MAX).
+ */
+long lp_feasibility_set_count_int(const lp_feasibility_set_t* set);
+
 /**
  * Check if the given value belongs to the set.
  */
 int lp_feasibility_set_contains(const lp_feasibility_set_t* set, const lp_value_t* value);
 
+/**
+ * Checks if the set contains an integer value.
+ */
+int lp_feasibility_set_contains_int(const lp_feasibility_set_t* set);
+
 /**
  * Pick a value from the feasible set (must be non-empty). If an integer value
  * is available it will be picked.
diff --git a/include/integer.h b/include/integer.h
index affcca51..efc49dff 100644
--- a/include/integer.h
+++ b/include/integer.h
@@ -129,6 +129,11 @@ int lp_integer_print_matrix(const lp_integer_t* c_array, size_t m, size_t n, FIL
  */
 char* lp_integer_to_string(const lp_integer_t* c);
 
+/**
+ * Returns true if the integer fits in long.
+ */
+int lp_integer_fits_int(const lp_integer_t* c);
+
 /**
  * Returns the int representation of the integer (truncated but keeps sign).
  */
diff --git a/include/interval.h b/include/interval.h
index f77c7830..f843b5b9 100644
--- a/include/interval.h
+++ b/include/interval.h
@@ -27,7 +27,7 @@ extern "C" {
 #endif
 
 /**
- * An interval (a, b) with both point being values. This side is open is _open
+ * An interval (a, b) with both point being values. A side is open if _open
  * is true. If interval is a point [a,a], then the value b is not used (it is
  * not constructed).
  */
@@ -83,6 +83,12 @@ void lp_interval_swap(lp_interval_t* I1, lp_interval_t* I2);
 /** Check if the value is contained in the interval */
 int lp_interval_contains(const lp_interval_t* I, const lp_value_t* v);
 
+/** Check if the interval contains an integer value */
+int lp_interval_contains_int(const lp_interval_t* I);
+
+/** Counts the number of integers in the interval (up to LONG_MAX) */
+long lp_interval_count_int(const lp_interval_t* I);
+
 /** Returns an approximation of the log interval size */
 int lp_interval_size_approx(const lp_interval_t* I);
 
@@ -113,9 +119,12 @@ void lp_interval_pick_value(const lp_interval_t* I, lp_value_t* v);
 /** Compares the lower bounds of the intervals */
 int lp_interval_cmp_lower_bounds(const lp_interval_t* I1, const lp_interval_t* I2);
 
-/** Compares the uppoer bounds of the intervals */
+/** Compares the upper bounds of the intervals */
 int lp_interval_cmp_upper_bounds(const lp_interval_t* I1, const lp_interval_t* I2);
 
+/** Compares an interval I with a value v. Returns 0 if v in I, 1 if v is below I, -1 if v is above I. */
+int lp_interval_cmp_value(const lp_interval_t* I, const lp_value_t* v);
+
 /**
  * Comparison of intervals based on upper bounds, with additional intersect info.
  */
diff --git a/include/polyxx/interval.h b/include/polyxx/interval.h
index 432cbdc0..66ddb4b7 100644
--- a/include/polyxx/interval.h
+++ b/include/polyxx/interval.h
@@ -81,6 +81,10 @@ namespace poly {
 
   /** Check whether an interval contains some value. */
   bool contains(const Interval& i, const Value& v);
+  /** Check whether an interval contains an integer value. */
+  bool contains_int(const Interval& i);
+  /** Counts the number of integers in the interval up to LONG_MAX. */
+  long count_int(const Interval& i);
   /** Get an approximation of the log interval size. */
   int log_size(const Interval& i);
 
diff --git a/python/polypyFeasibilitySet2.c b/python/polypyFeasibilitySet2.c
index 19589012..3705d8c1 100644
--- a/python/polypyFeasibilitySet2.c
+++ b/python/polypyFeasibilitySet2.c
@@ -40,9 +40,17 @@ FeasibilitySet_pick_value(PyObject* self);
 static PyObject*
 FeasibilitySet_intersect(PyObject* self, PyObject* args);
 
+static PyObject*
+FeasibilitySet_contains_value(PyObject* self, PyObject* args);
+
+static PyObject*
+FeasibilitySet_contains_int(PyObject* self);
+
 PyMethodDef FeasibilitySet_methods[] = {
     {"pick_value", (PyCFunction)FeasibilitySet_pick_value, METH_NOARGS, "Returns a value from the interval."},
-    {"intersect", (PyCFunction)FeasibilitySet_intersect, METH_VARARGS, "Returns the intersection of the interval with another one."},
+    {"intersect", (PyCFunction)FeasibilitySet_intersect, METH_VARARGS, "Returns the intersection of the set with another one."},
+    {"contains", (PyCFunction)FeasibilitySet_contains_value, METH_VARARGS, "Returns true if the value is in the feasibility set."},
+    {"contains_int", (PyCFunction)FeasibilitySet_contains_int, METH_NOARGS, "Returns true if the set contains an integer value."},
     {NULL}  /* Sentinel */
 };
 
@@ -159,9 +167,42 @@ FeasibilitySet_intersect(PyObject* self, PyObject* args) {
   lp_feasibility_set_t* S1 = ((FeasibilitySet*) self)->S;
   lp_feasibility_set_t* S2 = ((FeasibilitySet*) feasibility_set_obj)->S;
 
-  // The intersect
+  // The intersection
   lp_feasibility_set_t* P = lp_feasibility_set_intersect(S1, S2);
   PyObject* P_obj = PyFeasibilitySet_create(P);
 
   return P_obj;
 }
+
+static PyObject*
+FeasibilitySet_contains_value(PyObject* self, PyObject* args) {
+
+  if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
+    Py_INCREF(Py_NotImplemented);
+    return Py_NotImplemented;
+  }
+
+  PyObject* value_obj = PyTuple_GetItem(args, 0);
+
+  if (!PyValue_CHECK(value_obj)) {
+    Py_INCREF(Py_NotImplemented);
+    return Py_NotImplemented;
+  }
+
+  lp_feasibility_set_t* S = ((FeasibilitySet*) self)->S;
+  lp_value_t* v = &((Value*) value_obj)->v;
+
+  int result = lp_feasibility_set_contains(S, v);
+
+  PyObject* result_object = result ? Py_True : Py_False;
+  Py_INCREF(result_object);
+  return result_object;
+}
+
+static PyObject*
+FeasibilitySet_contains_int(PyObject* self) {
+  lp_feasibility_set_t* S = ((FeasibilitySet*) self)->S;
+  PyObject* result_object = lp_feasibility_set_contains_int(S) ? Py_True : Py_False;
+  Py_INCREF(result_object);
+  return result_object;
+}
diff --git a/python/polypyFeasibilitySet3.c b/python/polypyFeasibilitySet3.c
index e5ecf39f..979d4f11 100644
--- a/python/polypyFeasibilitySet3.c
+++ b/python/polypyFeasibilitySet3.c
@@ -40,9 +40,17 @@ FeasibilitySet_pick_value(PyObject* self);
 static PyObject*
 FeasibilitySet_intersect(PyObject* self, PyObject* args);
 
+static PyObject*
+FeasibilitySet_contains_value(PyObject* self, PyObject* args);
+
+static PyObject*
+FeasibilitySet_contains_int(PyObject* self);
+
 PyMethodDef FeasibilitySet_methods[] = {
     {"pick_value", (PyCFunction)FeasibilitySet_pick_value, METH_NOARGS, "Returns a value from the interval."},
-    {"intersect", (PyCFunction)FeasibilitySet_intersect, METH_VARARGS, "Returns the intersection of the interval with another one."},
+    {"intersect", (PyCFunction)FeasibilitySet_intersect, METH_VARARGS, "Returns the intersection of the set with another one."},
+    {"contains", (PyCFunction)FeasibilitySet_contains_value, METH_VARARGS, "Returns true if the value is in the feasibility set."},
+    {"contains_int", (PyCFunction)FeasibilitySet_contains_int, METH_NOARGS, "Returns true if the set contains an integer value."},
     {NULL}  /* Sentinel */
 };
 
@@ -98,8 +106,7 @@ PyTypeObject FeasibilitySetType = {
 };
 
 static void
-FeasibilitySet_dealloc(FeasibilitySet* self)
-{
+FeasibilitySet_dealloc(FeasibilitySet* self) {
   lp_feasibility_set_delete(self->S);
   ((PyObject*)self)->ob_type->tp_free((PyObject*)self);
 }
@@ -120,8 +127,7 @@ FeasibilitySet_new(PyTypeObject *type, PyObject *args, PyObject *kwds) {
 }
 
 static int
-FeasibilitySet_init(FeasibilitySet* self, PyObject* args)
-{
+FeasibilitySet_init(FeasibilitySet* self, PyObject* args) {
   if (PyTuple_Check(args) && PyTuple_Size(args) == 0) {
     // Defaults to (-inf, +inf)
     self->S = lp_feasibility_set_new_full();
@@ -168,9 +174,42 @@ FeasibilitySet_intersect(PyObject* self, PyObject* args) {
   lp_feasibility_set_t* S1 = ((FeasibilitySet*) self)->S;
   lp_feasibility_set_t* S2 = ((FeasibilitySet*) feasibility_set_obj)->S;
 
-  // The intersect
+  // The intersection
   lp_feasibility_set_t* P = lp_feasibility_set_intersect(S1, S2);
   PyObject* P_obj = PyFeasibilitySet_create(P);
 
   return P_obj;
 }
+
+static PyObject*
+FeasibilitySet_contains_value(PyObject* self, PyObject* args) {
+
+  if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
+    Py_INCREF(Py_NotImplemented);
+    return Py_NotImplemented;
+  }
+
+  PyObject* value_obj = PyTuple_GetItem(args, 0);
+
+  if (!PyValue_CHECK(value_obj)) {
+    Py_INCREF(Py_NotImplemented);
+    return Py_NotImplemented;
+  }
+
+  lp_feasibility_set_t* S = ((FeasibilitySet*) self)->S;
+  lp_value_t* v = &((Value*) value_obj)->v;
+
+  int result = lp_feasibility_set_contains(S, v);
+
+  PyObject* result_object = result ? Py_True : Py_False;
+  Py_INCREF(result_object);
+  return result_object;
+}
+
+static PyObject*
+FeasibilitySet_contains_int(PyObject* self) {
+  lp_feasibility_set_t* S = ((FeasibilitySet*) self)->S;
+  PyObject* result_object = lp_feasibility_set_contains_int(S) ? Py_True : Py_False;
+  Py_INCREF(result_object);
+  return result_object;
+}
diff --git a/python/polypyInterval2.c b/python/polypyInterval2.c
index aa6d90e3..8dbdca0b 100644
--- a/python/polypyInterval2.c
+++ b/python/polypyInterval2.c
@@ -40,9 +40,13 @@ Interval_pick_value(PyObject* self);
 static PyObject*
 Interval_contains_value(PyObject* self, PyObject* args);
 
+static PyObject*
+Interval_contains_int(PyObject* self);
+
 PyMethodDef Interval_methods[] = {
     {"pick_value", (PyCFunction)Interval_pick_value, METH_NOARGS, "Returns a value from the interval."},
     {"contains", (PyCFunction)Interval_contains_value, METH_VARARGS, "Returns true if the value is in the interval."},
+    {"contains_int", (PyCFunction)Interval_contains_int, METH_NOARGS, "Returns true if the interval contains an integer value."},
     {NULL}  /* Sentinel */
 };
 
@@ -89,8 +93,7 @@ PyTypeObject IntervalType = {
 };
 
 static void
-Interval_dealloc(Interval* self)
-{
+Interval_dealloc(Interval* self) {
   lp_interval_destruct(&self->I);
   self->ob_type->tp_free((PyObject*)self);
 }
@@ -163,6 +166,13 @@ Interval_contains_value(PyObject* self, PyObject* args) {
 
   PyObject* result_object = result ? Py_True : Py_False;
   Py_INCREF(result_object);
+  return result_object;
+}
 
+static PyObject*
+Interval_contains_int(PyObject* self) {
+  Interval* I = (Interval*) self;
+  PyObject* result_object = lp_interval_contains_int(&I->I) ? Py_True : Py_False;
+  Py_INCREF(result_object);
   return result_object;
 }
diff --git a/python/polypyInterval3.c b/python/polypyInterval3.c
index 260bcc86..adf34c11 100644
--- a/python/polypyInterval3.c
+++ b/python/polypyInterval3.c
@@ -40,9 +40,13 @@ Interval_pick_value(PyObject* self);
 static PyObject*
 Interval_contains_value(PyObject* self, PyObject* args);
 
+static PyObject*
+Interval_contains_int(PyObject* self);
+
 PyMethodDef Interval_methods[] = {
     {"pick_value", (PyCFunction)Interval_pick_value, METH_NOARGS, "Returns a value from the interval."},
     {"contains", (PyCFunction)Interval_contains_value, METH_VARARGS, "Returns true if the value is in the interval."},
+    {"contains_int", (PyCFunction)Interval_contains_int, METH_NOARGS, "Returns true if the interval contains an integer value."},
     {NULL}  /* Sentinel */
 };
 
@@ -98,8 +102,7 @@ PyTypeObject IntervalType = {
 };
 
 static void
-Interval_dealloc(Interval* self)
-{
+Interval_dealloc(Interval* self) {
   lp_interval_destruct(&self->I);
   ((PyObject*)self)->ob_type->tp_free((PyObject*)self);
 }
@@ -120,8 +123,7 @@ Interval_new(PyTypeObject *type, PyObject *args, PyObject *kwds) {
 }
 
 static int
-Interval_init(Interval* self, PyObject* args)
-{
+Interval_init(Interval* self, PyObject* args) {
   if (PyTuple_Check(args) && PyTuple_Size(args) == 0) {
     // Defaults to (-inf, +inf)
     lp_interval_construct_full(&self->I);
@@ -172,6 +174,13 @@ Interval_contains_value(PyObject* self, PyObject* args) {
 
   PyObject* result_object = result ? Py_True : Py_False;
   Py_INCREF(result_object);
+  return result_object;
+}
 
+static PyObject*
+Interval_contains_int(PyObject* self) {
+  Interval* I = (Interval*) self;
+  PyObject* result_object = lp_interval_contains_int(&I->I) ? Py_True : Py_False;
+  Py_INCREF(result_object);
   return result_object;
 }
diff --git a/src/interval/interval.c b/src/interval/interval.c
index 85594f71..dce080e3 100644
--- a/src/interval/interval.c
+++ b/src/interval/interval.c
@@ -86,8 +86,7 @@ void lp_interval_construct(lp_interval_t* I,
   }
 }
 
-void lp_rational_interval_construct_point(lp_rational_interval_t* I, const lp_rational_t* a)
-{
+void lp_rational_interval_construct_point(lp_rational_interval_t* I, const lp_rational_t* a) {
   rational_construct_copy(&I->a, a);
   I->a_open = 0;
   I->b_open = 0;
@@ -110,7 +109,8 @@ void lp_interval_construct_point(lp_interval_t* I, const lp_value_t* q) {
 
 void lp_rational_interval_construct_zero(lp_rational_interval_t* I) {
   rational_construct(&I->a);
-  I->a_open = I->b_open = 0;
+  I->a_open = 0;
+  I->b_open = 0;
   I->is_point = 1;
 }
 
@@ -673,19 +673,91 @@ int lp_rational_interval_contains_rational(const lp_rational_interval_t* I, cons
   return 1;
 }
 
-int lp_interval_contains(const lp_interval_t* I, const lp_value_t* v) {
+int lp_interval_cmp_value(const lp_interval_t* I, const lp_value_t* v) {
   int cmp_a_v = lp_value_cmp(&I->a, v);
   if (I->is_point) {
-    return cmp_a_v == 0;
+    return cmp_a_v;
   }
-  if (I->a_open && cmp_a_v >= 0) return 0;
-  if (!I->a_open && cmp_a_v > 0) return 0;
+  if (I->a_open && cmp_a_v >= 0) return 1;
+  if (!I->a_open && cmp_a_v > 0) return 1;
   int cmp_v_b = lp_value_cmp(v, &I->b);
-  if (I->b_open && cmp_v_b >= 0) return 0;
-  if (!I->b_open && cmp_v_b > 0) return 0;
-  return 1;
+  if (I->b_open && cmp_v_b >= 0) return -1;
+  if (!I->b_open && cmp_v_b > 0) return -1;
+  return 0;
 }
 
+int lp_interval_contains(const lp_interval_t* I, const lp_value_t* v) {
+  return lp_interval_cmp_value(I, v) == 0;
+}
+
+int lp_interval_contains_int(const lp_interval_t* I) {
+  if (lp_value_is_infinity(&I->a)) return 1;
+  int a_int = lp_value_is_integer(&I->a);
+  if (I->is_point) {
+    return a_int;
+  }
+  if (!I->a_open && a_int) return 1;
+  if (lp_value_is_infinity(&I->b)) return 1;
+  int b_int = lp_value_is_integer(&I->b);
+  if (!I->b_open && b_int) return 1;
+
+  lp_integer_t m, n;
+  lp_integer_construct(&m);
+  lp_integer_construct(&n);
+
+  lp_value_ceiling(&I->a, &m);
+  lp_value_floor(&I->b, &n);
+
+  if (a_int) lp_integer_inc(lp_Z, &m);
+  if (b_int) lp_integer_dec(lp_Z, &n);
+
+  int result = lp_integer_cmp(lp_Z, &n, &m) >= 0;
+
+  lp_integer_destruct(&m);
+  lp_integer_destruct(&n);
+
+  return result;
+}
+
+long lp_interval_count_int(const lp_interval_t* I) {
+  if (lp_value_is_infinity(&I->a)) return LONG_MAX;
+  int a_int = lp_value_is_integer(&I->a);
+  if (I->is_point) {
+    return a_int ? 1 : 0;
+  }
+  if (lp_value_is_infinity(&I->b)) return LONG_MAX;
+
+  long result = 0;
+  int b_int = lp_value_is_integer(&I->b);
+  if (!I->a_open && a_int) result++;
+  if (!I->b_open && b_int) result++;
+
+  lp_integer_t m, n;
+  lp_integer_construct(&m);
+  lp_integer_construct(&n);
+
+  lp_value_ceiling(&I->a, &m);
+  lp_value_floor(&I->b, &n);
+
+  if (a_int) lp_integer_inc(lp_Z, &m);
+  if (b_int) lp_integer_dec(lp_Z, &n);
+
+  lp_integer_sub(lp_Z, &n, &n, &m);
+  if (lp_integer_sgn(lp_Z, &n) >= 0) {
+    if (lp_integer_fits_int(&n)) {
+      result += lp_integer_to_int(&n) + 1;
+      // check for overflow
+      if (result < 0) result = LONG_MAX;
+    } else {
+      result = LONG_MAX;
+    }
+  }
+
+  lp_integer_destruct(&m);
+  lp_integer_destruct(&n);
+
+  return result;
+}
 
 int lp_dyadic_interval_contains_dyadic_rational(const lp_dyadic_interval_t* I, const lp_dyadic_rational_t* q) {
   int cmp_a_q = dyadic_rational_cmp(&I->a, q);
diff --git a/src/number/integer.c b/src/number/integer.c
index f0f87881..e89ecef3 100644
--- a/src/number/integer.c
+++ b/src/number/integer.c
@@ -153,6 +153,10 @@ char* lp_integer_to_string(const lp_integer_t* c) {
   return integer_to_string(c);
 }
 
+int lp_integer_fits_int(const lp_integer_t* c) {
+  return integer_fits_int(c);
+}
+
 long lp_integer_to_int(const lp_integer_t* c) {
   return integer_to_int(c);
 }
diff --git a/src/number/integer.h b/src/number/integer.h
index dc0d6645..d8e2dc8f 100644
--- a/src/number/integer.h
+++ b/src/number/integer.h
@@ -156,6 +156,11 @@ size_t integer_bits(const lp_integer_t* c) {
   return mpz_sizeinbase(c, 2);
 }
 
+static inline
+int integer_fits_int(const lp_integer_t* c) {
+  return mpz_fits_slong_p(c);
+}
+
 static inline
 long integer_to_int(const lp_integer_t* c) {
   return mpz_get_si(c);
diff --git a/src/number/value.c b/src/number/value.c
index cfd4d5f0..62f45fcd 100644
--- a/src/number/value.c
+++ b/src/number/value.c
@@ -418,7 +418,6 @@ int lp_value_is_integer(const lp_value_t* v) {
     return lp_dyadic_rational_is_integer(&v->value.dy_q);
   case LP_VALUE_RATIONAL:
     return lp_rational_is_integer(&v->value.q);
-    break;
   case LP_VALUE_ALGEBRAIC:
     return lp_algebraic_number_is_integer(&v->value.a);
   default:
@@ -670,7 +669,7 @@ void lp_value_get_value_between(const lp_value_t* a, int a_strict, const lp_valu
     } else {
       // Get the upper bound of the interval as a_ub
       rational_construct_from_dyadic(&a_ub, &a->value.a.I.b);
-      // Algebaic bound is strict so a_ub can be picked
+      // Algebraic bound is strict so a_ub can be picked
       a_ub_strict = 0;
     }
     break;
diff --git a/src/polynomial/feasibility_set.c b/src/polynomial/feasibility_set.c
index 79d919ec..7e59ab03 100644
--- a/src/polynomial/feasibility_set.c
+++ b/src/polynomial/feasibility_set.c
@@ -144,11 +144,38 @@ int lp_feasibility_set_is_point(const lp_feasibility_set_t* set) {
   return set->size == 1 && set->intervals->is_point;
 }
 
+int lp_feasibility_set_is_point_int(const lp_feasibility_set_t* set) {
+  long cnt = 0;
+  for (size_t i = 0; i < set->size; ++i) {
+    long tmp = lp_interval_count_int(set->intervals + i);
+    assert(tmp >= 0);
+    // checking tmp and the sum independently to avoid overflows
+    if (tmp > 1 || tmp + cnt > 1) {
+      return 0;
+    }
+    cnt += tmp;
+  }
+  return cnt == 1;
+}
+
+long lp_feasibility_set_count_int(const lp_feasibility_set_t* set) {
+  long cnt = 0;
+  for (size_t i = 0; i < set->size; ++i) {
+    long tmp = lp_interval_count_int(set->intervals + i);
+    assert(tmp >= 0);
+    // check for overflow
+    if (tmp >= LONG_MAX - cnt) {
+      return LONG_MAX;
+    }
+    cnt += tmp;
+  }
+  return cnt;
+}
+
 int lp_feasibility_set_print(const lp_feasibility_set_t* set, FILE* out) {
   int ret = 0;
-  size_t i;
   ret += fprintf(out, "{ ");
-  for(i = 0; i < set->size; ++ i) {
+  for(size_t i = 0; i < set->size; ++ i) {
     if (i) {
       ret += fprintf(out, ", ");
     }
@@ -170,9 +197,27 @@ char* lp_feasibility_set_to_string(const lp_feasibility_set_t* set) {
 }
 
 int lp_feasibility_set_contains(const lp_feasibility_set_t* set, const lp_value_t* value) {
-  // TODO: binary search
+  size_t l = 0, r = set->size;
+  while(r > l) {
+    size_t m = l + (r - l) / 2;
+    int cmp = lp_interval_cmp_value(set->intervals + m, value);
+    if (cmp > 0) {
+      // v below I
+      r = m;
+    } else if (cmp < 0) {
+      // v above I
+      l = m + 1;
+    } else {
+      // v in I
+      return 1;
+    }
+  }
+  return 0;
+}
+
+int lp_feasibility_set_contains_int(const lp_feasibility_set_t* set) {
   for (size_t i = 0; i < set->size; ++ i) {
-    if (lp_interval_contains(set->intervals + i, value)) {
+    if (lp_interval_contains_int(set->intervals + i)) {
       return 1;
     }
   }
@@ -457,7 +502,11 @@ int interval_sort_for_union(const void *I1_void, const void* I2_void) {
 
 void lp_feasibility_set_add(lp_feasibility_set_t* s, const lp_feasibility_set_t* from) {
 
-  if (from->size == 0) {
+  if (lp_feasibility_set_is_empty(from)) {
+    return;
+  }
+
+  if (lp_feasibility_set_is_full(s)) {
     return;
   }
 
diff --git a/src/polyxx/interval.cpp b/src/polyxx/interval.cpp
index f5dd72b8..867815ee 100644
--- a/src/polyxx/interval.cpp
+++ b/src/polyxx/interval.cpp
@@ -68,6 +68,15 @@ namespace poly {
   bool contains(const Interval& i, const Value& v) {
     return lp_interval_contains(i.get_internal(), v.get_internal());
   }
+
+  bool contains_int(const Interval& i) {
+    return lp_interval_contains_int(i.get_internal());
+  }
+
+  long count_int(const Interval& i) {
+    return lp_interval_count_int(i.get_internal());
+  }
+
   int log_size(const Interval& i) {
     return lp_interval_size_approx(i.get_internal());
   }
diff --git a/src/upolynomial/factorization.c b/src/upolynomial/factorization.c
index 44d9efd0..e165d6d2 100644
--- a/src/upolynomial/factorization.c
+++ b/src/upolynomial/factorization.c
@@ -111,7 +111,7 @@ lp_upolynomial_factors_t* upolynomial_factor_square_free_primitive(const lp_upol
     // 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
-    assert(mpz_fits_slong_p(&f->K->M));
+    assert(lp_integer_fits_int(&f->K->M));
     long p = integer_to_int(&f->K->M);
     lp_upolynomial_t* f_p = lp_upolynomial_div_degrees(f, p);
     factors = upolynomial_factor_square_free_primitive(f_p);
@@ -171,7 +171,7 @@ lp_upolynomial_factors_t* upolynomial_factor_square_free_primitive(const lp_upol
 
     // If P has content, it is a power of p
     if (lp_upolynomial_degree(P) > 0) {
-      assert(mpz_fits_slong_p(&f->K->M));
+      assert(lp_integer_fits_int(&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 = upolynomial_factor_square_free_primitive(P_p);
@@ -288,7 +288,7 @@ lp_upolynomial_factors_t* upolynomial_factor_distinct_degree(const lp_upolynomia
   assert(lp_upolynomial_is_monic(f));
 
   // The prime
-  assert(mpz_fits_slong_p(&K->M));
+  assert(lp_integer_fits_int(&K->M));
   long p = integer_to_int(&K->M);
   assert(p < FIELD_ORDER_LIMIT);
 
@@ -370,7 +370,7 @@ 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));
+  assert(lp_integer_fits_int(&K->M));
   size_t p = (size_t) integer_to_int(&K->M);
   assert(p < FIELD_ORDER_LIMIT);
 
@@ -757,7 +757,7 @@ 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
-    assert(mpz_fits_slong_p(&K->M));
+    assert(lp_integer_fits_int(&K->M));
     long x_int, p = integer_to_int(&K->M);
     assert(p < FIELD_ORDER_LIMIT);
     lp_upolynomial_t* linear_factors_product = 0;
diff --git a/src/upolynomial/root_finding.c b/src/upolynomial/root_finding.c
index c64d1b57..e02ba9c8 100644
--- a/src/upolynomial/root_finding.c
+++ b/src/upolynomial/root_finding.c
@@ -593,7 +593,7 @@ void upolynomial_roots_find_brute_force(const lp_upolynomial_t* f, lp_integer_t*
   size_t d = lp_upolynomial_degree(f);
   assert(d > 0);
 
-  assert(mpz_fits_slong_p(&K->M));
+  assert(lp_integer_fits_int(&K->M));
   long p = integer_to_int(&K->M);
   assert(p < FIELD_ORDER_LIMIT);
 
diff --git a/test/polyxx/test_interval.cpp b/test/polyxx/test_interval.cpp
index 8fbdc542..d8710997 100644
--- a/test/polyxx/test_interval.cpp
+++ b/test/polyxx/test_interval.cpp
@@ -53,6 +53,32 @@ TEST_CASE("interval::contains") {
   CHECK(contains(i, Value(Rational(5,2))));
 }
 
+TEST_CASE("interval::contains_int") {
+  CHECK_FALSE(contains_int(Interval(1,2)));
+  CHECK(contains_int(Interval(1,3)));
+  CHECK(contains_int(Interval()));
+  CHECK(contains_int(Interval(1)));
+  CHECK_FALSE(contains_int(Interval(Value(Rational(3,2)))));
+  CHECK(contains_int(Interval(1, false, 2, true)));
+  CHECK(contains_int(Interval(1, true, 2, false)));
+  CHECK_FALSE(contains_int(Interval(Rational(1,4), Rational(3,4))));
+}
+
+TEST_CASE("interval::count_int") {
+  CHECK(count_int(Interval(0,1)) == 0);
+  CHECK(count_int(Interval(0,2)) == 1);
+  CHECK(count_int(Interval()) == 1);
+  CHECK(count_int(Interval(Rational(1,2))) == 0);
+  CHECK(count_int(Interval(Rational(3,4), Rational(5,4))) == 1);
+  CHECK(count_int(Interval(0,false,2, false)) == 3);
+  CHECK(count_int(Interval(0,LONG_MAX)) == LONG_MAX-1);
+  CHECK(count_int(Interval(0,false,LONG_MAX,true)) == LONG_MAX);
+  CHECK(count_int(Interval(0,false,LONG_MAX,false)) == LONG_MAX);
+  CHECK(count_int(Interval(-1,LONG_MAX)) == LONG_MAX);
+  CHECK(count_int(Interval(-2,LONG_MAX)) == LONG_MAX);
+  CHECK(count_int(Interval(LONG_MIN,LONG_MAX)) == LONG_MAX);
+}
+
 TEST_CASE("interval::pick_value") {
   Interval i(1,3);
   CHECK(contains(i, pick_value(i)));