Feasibility set int (#80)

* WIP

WIP: intersection/union done

* intersection and printing done

* removed capacity

* added lp_polynomial_constraint_get_feasible_set_Zp

* added doctest tests for c files

* Update doc

* adds further test and fixes

* further tests

* Removed unnecessary header file

* Fixed assertion location

---------

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

CMakeLists.txt

@@ -82,6 +82,7 @@ endif()
 if (BUILD_TESTING)
   # Configure the C++ tests
   enable_testing()
+  add_subdirectory(test/poly)
   add_subdirectory(test/polyxx)
 endif()
 

include/feasibility_set_int.h

@@ -0,0 +1,178 @@
+/**
+ * Copyright 2024, 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"
+#include "integer.h"
+#include <stdbool.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * Set of disjoint intervals representing an algebraic set, ordered from
+ * left to right (-inf to +inf).
+ */
+struct lp_feasibility_set_int_struct {
+  /** The ring in which the feasibility set lives. */
+  lp_int_ring_t *K;
+
+  /** If true, set represents K \setminus elements */
+  bool inverted;
+
+  /** Number of elements */
+  size_t size;
+
+  /**
+   * Vector feasibility elements
+   * Values are normalized wrt to K, kept sorted and unique
+   */
+  lp_integer_t* elements;
+};
+
+/**
+ * Create a new feasibility set K.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_new_full(lp_int_ring_t *K);
+
+/**
+ * Create a new feasibility set {}.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_new_empty(lp_int_ring_t *K);
+
+/**
+ * Construct a copy.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_new_copy(const lp_feasibility_set_int_t* set);
+
+/**
+ * Construct from integers.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_new_from_integer(lp_int_ring_t *K, const lp_integer_t* integers, size_t integers_size, bool inverted);
+
+/**
+ * Delete the given feasibility set.
+ */
+void lp_feasibility_set_int_delete(lp_feasibility_set_int_t* set);
+
+/**
+ * Assignment.
+ */
+void lp_feasibility_set_int_assign(lp_feasibility_set_int_t* set, const lp_feasibility_set_int_t* from);
+
+/**
+ * Swap.
+ */
+void lp_feasibility_set_int_swap(lp_feasibility_set_int_t* s1, lp_feasibility_set_int_t* s2);
+
+/**
+ * Check if the given set is empty.
+ */
+int lp_feasibility_set_int_is_empty(const lp_feasibility_set_int_t* set);
+
+/**
+ * Check if the given set is full.
+ */
+int lp_feasibility_set_int_is_full(const lp_feasibility_set_int_t* set);
+
+/**
+ * Check if the set is a point {a}.
+ */
+int lp_feasibility_set_int_is_point(const lp_feasibility_set_int_t* set);
+
+/**
+ * assigns the size of the set to out
+ */
+void lp_feasibility_set_int_size(const lp_feasibility_set_int_t *set, lp_integer_t *out);
+
+/**
+ * returns the size; guaranteed to be correct if it fits in long
+ */
+size_t lp_feasibility_set_int_size_approx(const lp_feasibility_set_int_t *set);
+
+/**
+ * Check if the given value belongs to the set.
+ */
+int lp_feasibility_set_int_contains(const lp_feasibility_set_int_t* set, const lp_integer_t* value);
+
+/**
+ * Pick a value from the feasible set (must be non-empty).
+ */
+void lp_feasibility_set_int_pick_value(const lp_feasibility_set_int_t* set, lp_integer_t* value);
+
+/**
+ * Get intersection of the two sets.
+ * s1 and s2 must be over the same ring K.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_intersect(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2);
+
+/**
+ * Get union of the two sets.
+ * s1 and s2 must be over the same ring K.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_union(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2);
+
+typedef enum {
+  LP_FEASIBILITY_SET_INT_S1,
+  LP_FEASIBILITY_SET_INT_S2,
+  LP_FEASIBILITY_SET_INT_NEW,
+  LP_FEASIBILITY_SET_INT_EMPTY
+} lp_feasibility_set_int_status_t;
+
+/**
+ * Get intersection of the two sets, returns the status in the given variable.
+ * The set s1 is given precedence so LP_FEASIBILITY_SET_S2 is the
+ * status only if the intersect is not s1.
+ * s1 and s2 must be over the same ring K.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_intersect_with_status(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2, lp_feasibility_set_int_status_t * status);
+
+/**
+ * Get union of the two sets, returns the status in the given variable.
+ * The set s1 is given precedence so LP_FEASIBILITY_SET_S2 is the
+ * status only if the union is not s1.
+ * s1 and s2 must be over the same ring K.
+ */
+lp_feasibility_set_int_t* lp_feasibility_set_int_union_with_status(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2, lp_feasibility_set_int_status_t* status);
+
+/**
+ * Add one set to another, i.e. s = s \cup from.
+ */
+void lp_feasibility_set_int_add(lp_feasibility_set_int_t* s, const lp_feasibility_set_int_t* from);
+
+/**
+ * Returns true if both sets are equal
+ */
+bool lp_feasibility_set_int_eq(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2);
+
+/**
+ * Print the set.
+ */
+int lp_feasibility_set_int_print(const lp_feasibility_set_int_t* set, FILE* out);
+
+/**
+ * Return the string representation of the set.
+ */
+char* lp_feasibility_set_int_to_string(const lp_feasibility_set_int_t* set);
+
+#ifdef __cplusplus
+} /* close extern "C" { */
+#endif

include/poly.h

@@ -66,6 +66,7 @@ typedef struct lp_dyadic_interval_struct lp_dyadic_interval_t;
 typedef struct lp_interval_struct lp_interval_t;
 
 typedef struct lp_feasibility_set_struct lp_feasibility_set_t;
+typedef struct lp_feasibility_set_int_struct lp_feasibility_set_int_t;
 typedef struct lp_polynomial_hash_set_struct lp_polynomial_hash_set_t;
 typedef struct lp_polynomial_vector_struct lp_polynomial_vector_t;
 typedef struct lp_polynomial_heap_struct lp_polynomial_heap_t;

include/polynomial.h

@@ -343,6 +343,17 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t
  */
 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);
 
+/**
+ * Given a polynomial A(x1, ..., xn, y) with y being the top variable, a sign
+ * condition, and an assignment M that assigns x1, ..., xn, the function returns
+ * a subset of Zp where
+ *
+ *   sgn(A(M(x1), ..., M(xn), y)) = sgn_condition .
+ *
+ * If negated is true, the constraint is considered negated.
+ */
+lp_feasibility_set_int_t* lp_polynomial_constraint_get_feasible_set_Zp(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, int negated, const lp_assignment_t* M);
+
 /**
  * Given a polynomial A(x1, ..., xn) and a sign condition,  the function returns
  * tries to infer bounds on the variables and stores them into the given interval

src/CMakeLists.txt

@@ -34,6 +34,7 @@ set(poly_SOURCES
   polynomial/polynomial.c
   polynomial/polynomial_context.c
   polynomial/feasibility_set.c
+  polynomial/feasibility_set_int.c
   polynomial/polynomial_hash_set.c
   polynomial/polynomial_heap.c
   polynomial/polynomial_vector.c

src/polynomial/feasibility_set.c

@@ -32,17 +32,18 @@
 #include "polynomial/feasibility_set.h"
 
 static
-void lp_feasibility_set_ensure_capacity(lp_feasibility_set_t* s, size_t size) {
-  if (size && size > s->capacity) {
-    s->capacity = size;
+void lp_feasibility_set_ensure_capacity(lp_feasibility_set_t* s, size_t capacity) {
+  if (capacity && capacity > s->capacity) {
+    s->capacity = capacity;
     s->intervals = realloc(s->intervals, s->capacity * sizeof(lp_interval_t));
   }
 }
 
+static
 void lp_feasibility_set_construct(lp_feasibility_set_t* s, size_t size) {
   s->size = 0;
   s->capacity = 0;
-  s->intervals = 0;
+  s->intervals = NULL;
   lp_feasibility_set_ensure_capacity(s, size);
 }
 
@@ -86,6 +87,7 @@ void lp_feasibility_set_delete(lp_feasibility_set_t* set) {
   free(set);
 }
 
+static
 void lp_feasibility_set_construct_copy(lp_feasibility_set_t* set, const lp_feasibility_set_t* from) {
   lp_feasibility_set_construct(set, from->size);
   size_t i;

src/polynomial/feasibility_set.h

@@ -19,6 +19,4 @@
 
 #pragma once
 
-void lp_feasibility_set_construct(lp_feasibility_set_t* s, size_t size);
-
 lp_feasibility_set_t* lp_feasibility_set_new_internal(size_t size);