Implement binary search for feasibility_set::contains

include/interval.h

@@ -122,6 +122,9 @@ int lp_interval_cmp_lower_bounds(const lp_interval_t* I1, const lp_interval_t* I
 /** 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.
  */

src/interval/interval.c

@@ -673,17 +673,21 @@ 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) {

src/polynomial/feasibility_set.c

@@ -196,9 +196,18 @@ 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
-  for (size_t i = 0; i < set->size; ++ i) {
-    if (lp_interval_contains(set->intervals + i, value)) {
+  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;
     }
   }