diff --git a/CMakeLists.txt b/CMakeLists.txt index c22590e..32b5534 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/include/feasibility_set_int.h b/include/feasibility_set_int.h new file mode 100644 index 0000000..5858b15 --- /dev/null +++ b/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 . + */ + +#pragma once + +#include "poly.h" +#include "integer.h" +#include + +#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 diff --git a/include/poly.h b/include/poly.h index 99b43cb..557a959 100644 --- a/include/poly.h +++ b/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; diff --git a/include/polynomial.h b/include/polynomial.h index 734bfa9..58604ec 100644 --- a/include/polynomial.h +++ b/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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7c790e1..eb9a7e5 100755 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/polynomial/feasibility_set.c b/src/polynomial/feasibility_set.c index af7bea5..79d919e 100644 --- a/src/polynomial/feasibility_set.c +++ b/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; diff --git a/src/polynomial/feasibility_set.h b/src/polynomial/feasibility_set.h index 3ca631e..ea90efa 100644 --- a/src/polynomial/feasibility_set.h +++ b/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); diff --git a/src/polynomial/feasibility_set_int.c b/src/polynomial/feasibility_set_int.c new file mode 100644 index 0000000..dbd5ae0 --- /dev/null +++ b/src/polynomial/feasibility_set_int.c @@ -0,0 +1,589 @@ +/** + * 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 . + */ + +#include + +#include +#include + +#include "utils/u_memstream.h" + +static +int int_cmp(const void *i1, const void *i2) { + return lp_integer_cmp(lp_Z, i1, i2); +} + +static +void lp_feasibility_set_int_construct(lp_int_ring_t *K, lp_feasibility_set_int_t* s, size_t size, bool inverted) { + s->K = K; + lp_int_ring_attach(s->K); + s->size = size; + s->inverted = inverted; + s->elements = size ? malloc(s->size * sizeof(lp_integer_t)) : NULL; +} + +lp_feasibility_set_int_t* lp_feasibility_set_int_new_empty(lp_int_ring_t *K) { + lp_feasibility_set_int_t *s = malloc(sizeof(lp_feasibility_set_int_t)); + lp_feasibility_set_int_construct(K, s, 0, false); + return s; +} + +lp_feasibility_set_int_t* lp_feasibility_set_int_new_full(lp_int_ring_t *K) { + lp_feasibility_set_int_t *s = malloc(sizeof(lp_feasibility_set_int_t)); + lp_feasibility_set_int_construct(K, s, 0, true); + return s; +} + +static +size_t unique_lp_integer(lp_integer_t *b, size_t size) { + if (size <= 1) return size; + + size_t new_size = 1; + lp_integer_t *first = b, *last = b + size; + lp_integer_t *result = first; + while (++first != last) { + if (lp_integer_cmp(lp_Z, result, first) != 0) { + lp_integer_swap(++result, first); + ++new_size; + } + } + while (++result != last) { + lp_integer_destruct(result); + } + assert(new_size <= size); + return new_size; +} + +static +void lp_feasibility_set_int_construct_from_integer(lp_feasibility_set_int_t *set, lp_int_ring_t *K, const lp_integer_t *elements, size_t size, bool inverted) { + lp_feasibility_set_int_construct(K, set, size, inverted); + for (size_t i = 0; i < size; ++ i) { + lp_integer_construct_copy(K, set->elements + i, elements + i); + } + qsort(set->elements, size, sizeof(lp_integer_t), int_cmp); + size_t new_size = unique_lp_integer(set->elements, size); + if (new_size < size) { + set->elements = realloc(set->elements, new_size * sizeof(lp_integer_t)); + } + set->size = new_size; +} + +lp_feasibility_set_int_t* lp_feasibility_set_int_new_from_integer(lp_int_ring_t *K, const lp_integer_t* elements, size_t elements_size, bool inverted) { + lp_feasibility_set_int_t *result = malloc(sizeof(lp_feasibility_set_int_t)); + lp_feasibility_set_int_construct_from_integer(result, K, elements, elements_size, inverted); + return result; +} + +static +void lp_feasibility_set_int_construct_copy(lp_feasibility_set_int_t *set, const lp_feasibility_set_int_t* from) { + // we assume that from is valid, thus no sorting + lp_feasibility_set_int_construct(from->K, set, from->size, from->inverted); + for (size_t i = 0; i < from->size; ++ i) { + lp_integer_construct_copy(lp_Z, set->elements + i, from->elements + i); + } + set->size = from->size; +} + +lp_feasibility_set_int_t* lp_feasibility_set_int_new_copy(const lp_feasibility_set_int_t* set) { + lp_feasibility_set_int_t *result = malloc(sizeof(lp_feasibility_set_int_t)); + lp_feasibility_set_int_construct_copy(result, set); + return result; +} + +static +void lp_feasibility_set_int_destruct(lp_feasibility_set_int_t* set) { + lp_int_ring_detach(set->K); + free(set->elements); +} + +void lp_feasibility_set_int_delete(lp_feasibility_set_int_t* set) { + assert(set); + lp_feasibility_set_int_destruct(set); + free(set); +} + +void lp_feasibility_set_int_assign(lp_feasibility_set_int_t* set, const lp_feasibility_set_int_t* from) { + if (set != from) { + lp_feasibility_set_int_destruct(set); + lp_feasibility_set_int_construct_copy(set, from); + } +} + +void lp_feasibility_set_int_swap(lp_feasibility_set_int_t* s1, lp_feasibility_set_int_t* s2) { + lp_feasibility_set_int_t tmp = *s1; + *s1 = *s2; + *s2 = tmp; +} + +int lp_feasibility_set_int_is_empty(const lp_feasibility_set_int_t* set) { + if (!set->inverted && set->size == 0) return 1; + if (set->inverted && lp_integer_cmp_int(lp_Z, &set->K->M, set->size) == 0) return 1; + return 0; +} + +int lp_feasibility_set_int_is_full(const lp_feasibility_set_int_t* set) { + if (set->inverted && set->size == 0) return 1; + if (!set->inverted && lp_integer_cmp_int(lp_Z, &set->K->M, set->size) == 0) return 1; + return 0; +} + +void lp_feasibility_set_int_size(const lp_feasibility_set_int_t *set, lp_integer_t *out) { + lp_integer_assign_int(lp_Z, out, set->size); + if (set->inverted) { + lp_integer_sub(lp_Z, out, &set->K->M, out); + } +} + +size_t lp_feasibility_set_int_size_approx(const lp_feasibility_set_int_t *set) { + if (!set->inverted) { + return set->size; + } else { + if (mpz_fits_ulong_p(&set->K->M)) { + return lp_integer_to_int(&set->K->M) - set->size; + } else { + // size can't be big enough for an actual size of 0 or 1 + return ULONG_MAX; + } + } +} + +int lp_feasibility_set_int_is_point(const lp_feasibility_set_int_t* set) { + assert(lp_integer_cmp_int(lp_Z, &set->K->M, 2) > 0); + if (!set->inverted && set->size == 1) return 1; + if (set->inverted && lp_integer_cmp_int(lp_Z, &set->K->M, set->size + 1) == 0) return 1; + return 0; +} + +/** returns true if value is in elements. Assumes that the elements are sorted and that value is normalized wrt set->K */ +static +bool lp_feasibility_set_int_find(const lp_feasibility_set_int_t *set, const lp_integer_t *value) { + if (set->size == 0) { + return false; + } + assert(set->elements); + // including l and r + long l = 0, r = set->size - 1; + while (l <= r) { + long p = (r + l) / 2; + int cmp = lp_integer_cmp(lp_Z, set->elements + p, value); + if (cmp > 0) { + // continue left + r = p - 1; + } else if (cmp < 0) { + // continue right + l = p + 1; + } else { + // found + return true; + } + } + return false; +} + +int lp_feasibility_set_int_contains(const lp_feasibility_set_int_t* set, const lp_integer_t* value) { + lp_integer_t value_normalized; + // normalize value before check + lp_integer_construct_copy(set->K, &value_normalized, value); + bool found = lp_feasibility_set_int_find(set, &value_normalized); + lp_integer_destruct(&value_normalized); + return found != set->inverted; +} + +void lp_feasibility_set_int_pick_value(const lp_feasibility_set_int_t* set, lp_integer_t* value) { + assert(!lp_feasibility_set_int_is_empty(set)); + if (!set->inverted) { + size_t pos = random() % set->size; + lp_integer_assign(lp_Z, value, set->elements + pos); + } else { + lp_integer_construct_from_int(lp_Z, value, 0); + // check 0 + if (!lp_feasibility_set_int_find(set, value)) { + return; + } + while (true) { + lp_integer_inc(lp_Z, value); + assert(lp_integer_in_ring(set->K, value)); + if (!lp_feasibility_set_int_find(set, value)) { + return; + } + lp_integer_neg(lp_Z, value, value); + if (!lp_feasibility_set_int_find(set, value)) { + return; + } + lp_integer_neg(lp_Z, value, value); + } + // TODO proper implementation + // get random element between (incl.) 0 and (|K| - size) + // find and add number of <= element + assert(false); + } +} + +typedef enum { + NONE = 0, + S1 = 1, + S2 = 2, + BOTH = 3, +} set_status_internal_t; + +/** Calculates i1 \cup i2 */ +static +set_status_internal_t ordered_integer_set_union( + lp_integer_t **result, size_t *result_size, + const lp_integer_t *i1, size_t i1_size, + const lp_integer_t *i2, size_t i2_size) { + + size_t max_size = i1_size + i2_size; + if (i1_size == 0 && i2_size == 0) { + if (result_size) { *result_size = 0; } + return BOTH; + } + + if (result) { *result = malloc(max_size * sizeof(lp_integer_t)); } + + bool just_i1 = true, just_i2 = true; + size_t p1 = 0, p2 = 0, pr = 0; + while (p1 < i1_size && p2 < i2_size) { + int cmp = lp_integer_cmp(lp_Z, i1 + p1, i2 + p2); + if (cmp < 0) { + if (result) { lp_integer_construct_copy(lp_Z, *result + pr, i1 + p1); } + just_i2 = false; + p1 ++; + } else if (cmp > 0) { + if (result) { lp_integer_construct_copy(lp_Z, *result + pr, i2 + p2); } + just_i1 = false; + p2 ++; + } else { + if (result) { lp_integer_construct_copy(lp_Z, *result + pr, i1 + p1); } + p1 ++; p2 ++; + } + pr ++; + } + if (result) { + while (p1 < i1_size) { + lp_integer_construct_copy(lp_Z, *result + pr, i1 + p1); + p1++; pr++; + just_i2 = false; + } + while (p2 < i2_size) { + lp_integer_construct_copy(lp_Z, *result + pr, i2 + p2); + p2++; pr++; + just_i1 = false; + } + *result = realloc(*result, pr * sizeof(lp_integer_t)); + } else { + if (p1 < i1_size) { just_i2 = false; pr += (i1_size - p1); } + if (p2 < i2_size) { just_i1 = false; pr += (i2_size - p2); } + } + + if (result_size) { *result_size = pr; } + return (just_i1 ? S1 : NONE) | (just_i2 ? S2 : NONE); +} + +/** Calculates i1 \cap i2 */ +static +set_status_internal_t ordered_integer_set_intersect( + lp_integer_t **result, size_t *result_size, + const lp_integer_t *i1, size_t i1_size, + const lp_integer_t *i2, size_t i2_size) { + + size_t max_size = i1_size < i2_size ? i1_size : i2_size; + if (i1_size == 0 || i2_size == 0) { + if (result_size) { *result_size = 0; } + return (i1_size == 0 ? S1 : NONE) | (i2_size == 0 ? S2 : NONE); + } + + if (result) { *result = malloc(max_size * sizeof(lp_integer_t)); } + + size_t p1 = 0, p2 = 0, pr = 0; + bool all_i1 = true, all_i2 = true; + while (p1 < i1_size && p2 < i2_size) { + int cmp = lp_integer_cmp(lp_Z, i1 + p1, i2 + p2); + if (cmp > 0) { + p2 ++; + all_i2 = false; + } else if (cmp < 0) { + p1 ++; + all_i1 = false; + } else { + if (result) { lp_integer_construct_copy(lp_Z, *result + pr, i1 + p1); } + p1 ++; p2 ++; pr ++; + } + } + if (p1 < i1_size) { all_i1 = false; } + if (p2 < i2_size) { all_i2 = false; } + + if (result) { *result = realloc(*result, pr * sizeof(lp_integer_t)); } + if (result_size) { *result_size = pr; } + return (all_i1 ? S1 : NONE) | (all_i2 ? S2 : NONE); +} + +/** Calculates i1 \setminus i2 */ +static +set_status_internal_t ordered_integer_set_minus( + lp_integer_t **result, size_t *result_size, + const lp_integer_t *i1, size_t i1_size, + const lp_integer_t *i2, size_t i2_size) { + + size_t max_size = i1_size; + if (result) { *result = malloc(max_size * sizeof(lp_integer_t)); } + + size_t p1 = 0, p2 = 0, pr = 0; + while (p1 < i1_size) { + bool found = false; + while(p2 < i2_size) { + // as long i2's value is smaller than current i1, we continue + int cmp = lp_integer_cmp(lp_Z, i1 + p1, i2 + p2); + if (cmp == 0) { + found = true; + } + if (cmp <= 0) { + break; + } + p2 ++; + } + if (!found) { + if (result) { lp_integer_construct_copy(lp_Z, *result + pr, i1 + p1); } + pr ++; + } + p1 ++; + } + + if (result) { *result = realloc(*result, pr * sizeof(lp_integer_t)); } + if (result_size) { *result_size = pr; } + return (pr == i1_size ? S1 : NONE); +} + +static +void lp_feasibility_set_int_invert(lp_feasibility_set_int_t *set) { + assert(set->K != lp_Z); + // don't invert big fields + assert(lp_integer_cmp_int(lp_Z, &set->K->M, 10000) < 0); + + size_t cnt = lp_integer_to_int(&set->K->M); + assert(set->size <= cnt); + cnt -= set->size; + + lp_integer_t *new = malloc(cnt * sizeof(lp_integer_t)); + lp_integer_t *old = set->elements; + size_t pos_new = 0; + size_t pos_old = 0; + + long lb = lp_integer_to_int(&set->K->lb); + long ub = lp_integer_to_int(&set->K->ub); + for (long val = lb; val <= ub; ++val) { + assert(pos_old >= set->size || lp_integer_cmp_int(lp_Z, old + pos_old, val) >= 0); + if (pos_old < set->size && lp_integer_cmp_int(lp_Z, old + pos_old, val) == 0) { + ++ pos_old; + } else { + assert(pos_new < cnt); + lp_integer_construct_from_int(lp_Z, new + pos_new, val); + ++ pos_new; + } + } + + free(old); + set->elements = new; + set->size = cnt; + set->inverted = !set->inverted; +} + +static inline +set_status_internal_t invert_i1_i2(set_status_internal_t status) { + return (status & S1 ? S2 : NONE) | (status & S2 ? S1: NONE); +} + +static inline +lp_feasibility_set_int_status_t status_internal_to_external(set_status_internal_t is) { + switch (is) { + case BOTH: + case S1: + return LP_FEASIBILITY_SET_INT_S1; + case S2: + return LP_FEASIBILITY_SET_INT_S2; + default: + case NONE: + return LP_FEASIBILITY_SET_INT_NEW; + } +} + +static +lp_feasibility_set_int_t* lp_feasibility_set_int_intersect_internal(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2, set_status_internal_t *status) { + assert(lp_int_ring_equal(s1->K, s2->K)); + + lp_feasibility_set_int_t *result; + + if (s1->inverted && s2->inverted) { + result = lp_feasibility_set_int_new_empty(s1->K); + *status = ordered_integer_set_union(&result->elements, &result->size, + s1->elements, s1->size, s2->elements, + s2->size); + result->inverted = true; + } else if (!s1->inverted && !s2->inverted) { + result = lp_feasibility_set_int_new_empty(s1->K); + *status = ordered_integer_set_intersect(&result->elements, &result->size, + s1->elements, s1->size, s2->elements, + s2->size); + result->inverted = false; + } else if (s1->inverted && !s2->inverted) { + result = lp_feasibility_set_int_intersect_internal(s2, s1, status); + *status = invert_i1_i2(*status); + } else { + assert(!s1->inverted && s2->inverted); + if (s1->size > lp_feasibility_set_int_size_approx(s2)) { + // s2 could be the smaller set, LP_FEASIBILITY_SET_INT_S2 is possible + // TODO this effort could be saved in case we don't care about the status + lp_feasibility_set_int_t *tmp = lp_feasibility_set_int_new_copy(s2); + lp_feasibility_set_int_invert(tmp); + result = lp_feasibility_set_int_intersect_internal(s1, tmp, status); + lp_feasibility_set_int_delete(tmp); + } else { + // s1 is the smaller set, LP_FEASIBILITY_SET_INT_S2 is not possible + result = lp_feasibility_set_int_new_empty(s1->K); + *status = ordered_integer_set_minus(&result->elements, &result->size, + s1->elements, s1->size, + s2->elements, s2->size); + result->inverted = false; + } + } + + return result; +} + +lp_feasibility_set_int_t* lp_feasibility_set_int_intersect(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2) { + set_status_internal_t is; + return lp_feasibility_set_int_intersect_internal(s1, s2, &is); +} + +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) { + set_status_internal_t is; + lp_feasibility_set_int_t *result = lp_feasibility_set_int_intersect_internal(s1, s2, &is); + *status = lp_feasibility_set_int_is_empty(result) ? LP_FEASIBILITY_SET_INT_EMPTY : status_internal_to_external(is); + return result; +} + +static +lp_feasibility_set_int_t* lp_feasibility_set_int_union_internal(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2, set_status_internal_t * status) { + assert(lp_int_ring_equal(s1->K, s2->K)); + + lp_feasibility_set_int_t *result; + + if (s1->inverted && s2->inverted) { + result = lp_feasibility_set_int_new_empty(s1->K); + *status = ordered_integer_set_intersect(&result->elements, &result->size, + s1->elements, s1->size, + s2->elements, s2->size); + result->inverted = true; + } else if (!s1->inverted && !s2->inverted) { + result = lp_feasibility_set_int_new_empty(s1->K); + *status = ordered_integer_set_union(&result->elements, &result->size, + s1->elements, s1->size, + s2->elements, s2->size); + result->inverted = false; + } else if (!s1->inverted && s2->inverted) { + result = lp_feasibility_set_int_union_internal(s2, s1, status); + *status = invert_i1_i2(*status); + } else { + assert (s1->inverted && !s2->inverted); + if (s2->size > lp_feasibility_set_int_size_approx(s1)) { + // s2 is be the bigger set, LP_FEASIBILITY_SET_INT_S2 is possible + // TODO this effort could be saved in case we don't care about the status + lp_feasibility_set_int_t *tmp = lp_feasibility_set_int_new_copy(s1); + lp_feasibility_set_int_invert(tmp); + result = lp_feasibility_set_int_union_internal(tmp, s2, status); + lp_feasibility_set_int_delete(tmp); + } else { + result = lp_feasibility_set_int_new_empty(s1->K); + // s1 is the bigger set, LP_FEASIBILITY_SET_INT_S2 is not possible + *status = ordered_integer_set_minus(&result->elements, &result->size, + s1->elements, s1->size, + s2->elements, s2->size); + result->inverted = true; + } + } + + return result; +} + +lp_feasibility_set_int_t* lp_feasibility_set_int_union(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2) { + set_status_internal_t status; + return lp_feasibility_set_int_union_internal(s1, s2, &status); +} + +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) { + set_status_internal_t is; + lp_feasibility_set_int_t *result = lp_feasibility_set_int_union_internal(s1, s2, &is); + *status = lp_feasibility_set_int_is_empty(result) ? LP_FEASIBILITY_SET_INT_EMPTY : status_internal_to_external(is); + return result; +} + +void lp_feasibility_set_int_add(lp_feasibility_set_int_t* s, const lp_feasibility_set_int_t* from) { + lp_feasibility_set_int_t *tmp = lp_feasibility_set_int_union(s, from); + lp_feasibility_set_int_swap(tmp, s); + lp_feasibility_set_int_delete(tmp); +} + +bool lp_feasibility_set_int_eq(const lp_feasibility_set_int_t* s1, const lp_feasibility_set_int_t* s2) { + if (s1->inverted == s2->inverted) { + if (s1->size != s2->size) { + return false; + } + for (size_t i = 0; i < s1->size; ++i) { + if (lp_integer_cmp(lp_Z, s1->elements + i, s2->elements + i) != 0) { + return false; + } + } + return true; + } else { + if (lp_feasibility_set_int_size_approx(s1) != lp_feasibility_set_int_size_approx(s2)) { + return false; + } + size_t count; + ordered_integer_set_intersect(NULL, &count, s1->elements, s1->size, s2->elements, s2->size); + return count == 0; + } +} + +int lp_feasibility_set_int_print(const lp_feasibility_set_int_t* set, FILE* out) { + int ret = 0; + if (set->inverted) { + ret += fprintf(out, "F \\"); + } + ret += fprintf(out, "{ "); + for(size_t i = 0; i < set->size; ++ i) { + if (i) { + ret += fprintf(out, ", "); + } + ret += lp_integer_print(&set->elements[i], out); + } + ret += fprintf(out, " } in "); + ret += lp_int_ring_print(set->K, out); + return ret; +} + +char* lp_feasibility_set_int_to_string(const lp_feasibility_set_int_t* set) { + struct u_memstream mem; + char* str = 0; + size_t size = 0; + u_memstream_open(&mem, &str, &size); + FILE* f = u_memstream_get(&mem); + lp_feasibility_set_int_print(set, f); + u_memstream_close(&mem); + return str; +} diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c index 80c07c2..5c3474f 100644 --- a/src/polynomial/polynomial.c +++ b/src/polynomial/polynomial.c @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -2036,6 +2037,50 @@ lp_feasibility_set_t* lp_polynomial_root_constraint_get_feasible_set(const lp_po } +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) { + const lp_polynomial_context_t *ctx = lp_polynomial_get_context(A); + assert(ctx->K != lp_Z); + assert(lp_sign_condition_Zp_valid(sgn_condition)); + + // The ring + lp_int_ring_t *K = ctx->K; + + // Negate the constraint if negated + if (negated) { + sgn_condition = lp_sign_condition_negate(sgn_condition); + } + + assert(lp_polynomial_is_univariate_m(A, M)); + + lp_upolynomial_t *upoly = lp_polynomial_to_univariate_m(A, M); + if (lp_upolynomial_degree(upoly) == 0) { + if (lp_upolynomial_is_zero(upoly)) { + return sgn_condition == LP_SGN_EQ_0 ? lp_feasibility_set_int_new_full(K) : lp_feasibility_set_int_new_empty(K); + } else { + return sgn_condition == LP_SGN_EQ_0 ? lp_feasibility_set_int_new_empty(K) : lp_feasibility_set_int_new_full(K); + } + } + + lp_feasibility_set_int_t *result = sgn_condition == LP_SGN_EQ_0 ? lp_feasibility_set_int_new_empty(K) : lp_feasibility_set_int_new_full(K); + + lp_upolynomial_roots_find_Zp(upoly, &result->elements, &result->size); + +#ifndef NDEBUG + { + lp_integer_t tmp; + lp_integer_construct(&tmp); + for (size_t i = 0; i < result->size; ++i) { + lp_upolynomial_evaluate_at_integer(upoly, &result->elements[i], &tmp); + assert(lp_integer_is_zero(lp_upolynomial_ring(upoly) , &tmp)); + } + lp_integer_destruct(&tmp); + } +#endif + + lp_upolynomial_delete(upoly); + return result; +} + void lp_polynomial_get_variables(const lp_polynomial_t* A, lp_variable_list_t* vars) { coefficient_get_variables(&A->data, vars); } diff --git a/test/poly/CMakeLists.txt b/test/poly/CMakeLists.txt new file mode 100644 index 0000000..49bcf87 --- /dev/null +++ b/test/poly/CMakeLists.txt @@ -0,0 +1,14 @@ +set(poly_tests + test_feasible_int_set +) + +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -Wextra -std=c++11") + +include_directories(${GMP_INCLUDE_DIR}) + +foreach(file ${poly_tests}) + add_executable(${file} ${file}.cpp) + target_include_directories(${file} PUBLIC ${CMAKE_SOURCE_DIR}/include) + add_test(NAME ${file} COMMAND ${file}) + target_link_libraries(${file} poly polyxx) +endforeach() diff --git a/test/poly/test_feasible_int_set.cpp b/test/poly/test_feasible_int_set.cpp new file mode 100644 index 0000000..0108dd2 --- /dev/null +++ b/test/poly/test_feasible_int_set.cpp @@ -0,0 +1,718 @@ +#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN +#include "../polyxx/doctest.h" + +#include +#include + +#include +using namespace poly; +using namespace std; + +/** + * Struct to handle a list of lp_integer_t (and proper deallocate it) + */ +struct integer_list { +private: + vector *v; +public: + integer_list() : v(new vector()) {} + integer_list(initializer_list values) : integer_list() { + for (int x : values) { + lp_integer_t tmp; + lp_integer_construct_from_int(lp_Z, &tmp, x); + this->v->push_back(tmp); + } + } + integer_list(const integer_list &other) : integer_list() { + for (const lp_integer_t x : *other.v) { + lp_integer_t tmp; + lp_integer_construct_copy(lp_Z, &tmp, &x); + this->v->push_back(tmp); + } + } + integer_list(integer_list &&other) noexcept : v(other.v) { + other.v = nullptr; + } + ~integer_list() { + if (this->v == nullptr) return; + for (lp_integer_t i : *this->v) { + lp_integer_destruct(&i); + } + delete this->v; + } + lp_feasibility_set_int_t* gen_set(lp_int_ring_t *K, bool inverted) const { + return lp_feasibility_set_int_new_from_integer(K, this->v->data(), this->v->size(), inverted); + } + lp_feasibility_set_int_t* gen_set(IntegerRing &K, bool inverted) const { + return this->gen_set(K.get_internal(), inverted); + } + const lp_integer_t* data() const { return this->v->data(); } + size_t size() const { return this->v->size(); } +}; + + +TEST_CASE("feasibility_set_int::construct") { + lp_integer_t prime; + lp_integer_construct_from_int(lp_Z, &prime, 13); + lp_int_ring_t *K = lp_int_ring_create(&prime, 1); + + SUBCASE("feasibility_set_int::construct::full") { + lp_feasibility_set_int_t *set = lp_feasibility_set_int_new_full(K); + CHECK(!lp_feasibility_set_int_is_empty(set)); + CHECK(lp_feasibility_set_int_is_full(set)); + CHECK(!lp_feasibility_set_int_is_point(set)); + lp_feasibility_set_int_delete(set); + } + SUBCASE("feasibility_set_int::construct::empty") { + lp_feasibility_set_int_t *set = lp_feasibility_set_int_new_empty(K); + CHECK(lp_feasibility_set_int_is_empty(set)); + CHECK(!lp_feasibility_set_int_is_full(set)); + CHECK(!lp_feasibility_set_int_is_point(set)); + lp_feasibility_set_int_delete(set); + } + SUBCASE("feasibility_set_int::construct::integer") { + integer_list l({2,1,4,3}); + lp_feasibility_set_int_t *set = lp_feasibility_set_int_new_from_integer(K, l.data(), l.size(), false); + CHECK_FALSE(lp_feasibility_set_int_is_empty(set)); + CHECK_FALSE(lp_feasibility_set_int_is_full(set)); + CHECK(lp_feasibility_set_int_size_approx(set) == 4); + lp_feasibility_set_int_delete(set); + } + SUBCASE("feasibility_set_int::construct::point") { + integer_list l({2}); + lp_feasibility_set_int_t *set = lp_feasibility_set_int_new_from_integer(K, l.data(), l.size(), false); + CHECK(!lp_feasibility_set_int_is_empty(set)); + CHECK(!lp_feasibility_set_int_is_full(set)); + CHECK(lp_feasibility_set_int_is_point(set)); + lp_feasibility_set_int_delete(set); + } + SUBCASE("feasibility_set_int::construct::point_inv") { + integer_list l({0,1,/*2,*/3,4,5,6,7,8,9,10,11,12}); + lp_feasibility_set_int_t *set = lp_feasibility_set_int_new_from_integer(K, l.data(), l.size(), true); + CHECK(!lp_feasibility_set_int_is_empty(set)); + CHECK(!lp_feasibility_set_int_is_full(set)); + CHECK(lp_feasibility_set_int_is_point(set)); + lp_feasibility_set_int_delete(set); + } + SUBCASE("feasibility_set_int::construct::integer::invert") { + integer_list l({2,1,4,3}); + lp_feasibility_set_int_t *set = lp_feasibility_set_int_new_from_integer(K, l.data(), l.size(), true); + CHECK_FALSE(lp_feasibility_set_int_is_empty(set)); + CHECK_FALSE(lp_feasibility_set_int_is_full(set)); + CHECK(lp_feasibility_set_int_size_approx(set) == 9); + CHECK(!lp_feasibility_set_int_contains(set, Integer(3).get_internal())); + CHECK(!lp_feasibility_set_int_contains(set, Integer(1).get_internal())); + CHECK(lp_feasibility_set_int_contains(set, Integer(-1).get_internal())); + CHECK(lp_feasibility_set_int_contains(set, Integer(13).get_internal())); + CHECK(!lp_feasibility_set_int_contains(set, Integer(14).get_internal())); + lp_feasibility_set_int_delete(set); + } + SUBCASE("feasibility_set_int::construct::integer::duplicates") { + integer_list l({1,2,2,14}); + lp_feasibility_set_int_t *set = lp_feasibility_set_int_new_from_integer(K, l.data(), l.size(), false); + CHECK_FALSE(lp_feasibility_set_int_is_empty(set)); + CHECK_FALSE(lp_feasibility_set_int_is_full(set)); + CHECK(lp_feasibility_set_int_size_approx(set) == 2); + lp_feasibility_set_int_delete(set); + } + + CHECK_EQ(K->ref_count, 1); + lp_int_ring_detach(K); + lp_integer_destruct(&prime); +} + +TEST_CASE("feasibility_set_int::intersect::inv") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set2 = integer_list({2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *set2_inv = integer_list({0,1,4,5,6}).gen_set(K, true); + + SUBCASE("feasibility_set_int::intersect::inv1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::inv2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::inv3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::inv4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(set2); + lp_feasibility_set_int_delete(set2_inv); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::intersect::inv_dual") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set2 = integer_list({1,2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set1 = integer_list({2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set2_inv = integer_list({0,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *set1_inv = integer_list({0,1,4,5,6}).gen_set(K, true); + + SUBCASE("feasibility_set_int::intersect::inv_dual1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::inv_dual2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::inv_dual") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::inv4_dual") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(set2); + lp_feasibility_set_int_delete(set2_inv); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::union::inv") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set2 = integer_list({2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *set2_inv = integer_list({0,1,4,5,6}).gen_set(K, true); + + SUBCASE("feasibility_set_int::union::inv1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union::inv2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union::inv3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union::inv4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(set2); + lp_feasibility_set_int_delete(set2_inv); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::union::inv_dual") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set2 = integer_list({1,2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set1 = integer_list({2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set2_inv = integer_list({0,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *set1_inv = integer_list({0,1,4,5,6}).gen_set(K, true); + + SUBCASE("feasibility_set_int::union::inv_dual1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union::inv_dual2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union::inv_dual3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union::inv_dual4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set2)); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(set2); + lp_feasibility_set_int_delete(set2_inv); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::union") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set2 = integer_list({2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *set2_inv = integer_list({0,1,4,5,6}).gen_set(K, true); + + SUBCASE("feasibility_set_int::union1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 3); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(set2); + lp_feasibility_set_int_delete(set2_inv); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::intersect") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set2 = integer_list({2,3}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *set2_inv = integer_list({0,1,4,5,6}).gen_set(K, true); + + SUBCASE("feasibility_set_int::intersect1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_NEW); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(set2); + lp_feasibility_set_int_delete(set2_inv); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::intersect_wrap_around") { + Integer prime(3); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,-1}).gen_set(K, false); + lp_feasibility_set_int_t *set2 = integer_list({-1}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0}).gen_set(K, true); + lp_feasibility_set_int_t *set2_inv = integer_list({0,1}).gen_set(K, true); + + SUBCASE("feasibility_set_int::intersect_wrap_around1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_contains(s, Integer(K, -1).get_internal())); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect_wrap_around2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_contains(s, Integer(K, -1).get_internal())); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect_wrap_around3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_contains(s, Integer(K, -1).get_internal())); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect_wrap_around2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, set2_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 1); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_contains(s, Integer(K, -1).get_internal())); + lp_feasibility_set_int_delete(s); + } + + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(set2); + lp_feasibility_set_int_delete(set2_inv); + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::intersect_empty") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *empty = lp_feasibility_set_int_new_empty(K.get_internal()); + + SUBCASE("feasibility_set_int::intersect_empty1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, empty, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 0); + CHECK(status == LP_FEASIBILITY_SET_INT_EMPTY); + CHECK(lp_feasibility_set_int_is_empty(s)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect_empty2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(empty, set1, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 0); + CHECK(status == LP_FEASIBILITY_SET_INT_EMPTY); + CHECK(lp_feasibility_set_int_is_empty(s)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect_empty3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, empty, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 0); + CHECK(status == LP_FEASIBILITY_SET_INT_EMPTY); + CHECK(lp_feasibility_set_int_is_empty(s)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect_empty4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(empty, set1_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 0); + CHECK(status == LP_FEASIBILITY_SET_INT_EMPTY); + CHECK(lp_feasibility_set_int_is_empty(s)); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(empty); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::intersect_full") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *full = lp_feasibility_set_int_new_full(K.get_internal()); + + SUBCASE("feasibility_set_int::intersect::intersect_full1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1, full, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::intersect_full2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(full, set1, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::intersect_full3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(set1_inv, full, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::intersect::intersect_full4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_intersect_with_status(full, set1_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(full); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::union_full") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *full = lp_feasibility_set_int_new_full(K.get_internal()); + + SUBCASE("feasibility_set_int::union_full1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, full, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 7); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_is_full(s)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union_full2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(full, set1, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 7); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_is_full(s)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union_full3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, full, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 7); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_is_full(s)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union_full4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(full, set1_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 7); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_is_full(s)); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(full); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::union_empty") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *empty = lp_feasibility_set_int_new_empty(K.get_internal()); + + SUBCASE("feasibility_set_int::union_empty1") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1, empty, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union_empty2") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(empty, set1, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union_empty3") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(set1_inv, empty, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S1); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + SUBCASE("feasibility_set_int::union_empty4") { + lp_feasibility_set_int_status_t status; + lp_feasibility_set_int_t *s = lp_feasibility_set_int_union_with_status(empty, set1_inv, &status); + CHECK(lp_feasibility_set_int_size_approx(s) == 2); + CHECK(status == LP_FEASIBILITY_SET_INT_S2); + CHECK(lp_feasibility_set_int_eq(s, set1)); + lp_feasibility_set_int_delete(s); + } + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(empty); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::size") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + lp_feasibility_set_int_t *empty = lp_feasibility_set_int_new_empty(K.get_internal()); + lp_feasibility_set_int_t *full = lp_feasibility_set_int_new_full(K.get_internal()); + + SUBCASE("feasibility_set_int::size::1") { + Integer size; + lp_feasibility_set_int_size(set1, size.get_internal()); + size_t size_approx = lp_feasibility_set_int_size_approx(set1); + CHECK(size == 2); + CHECK(size_approx == 2); + } + SUBCASE("feasibility_set_int::size::2") { + Integer size; + lp_feasibility_set_int_size(set1_inv, size.get_internal()); + size_t size_approx = lp_feasibility_set_int_size_approx(set1_inv); + CHECK(size == 2); + CHECK(size_approx == 2); + } + SUBCASE("feasibility_set_int::size::full") { + Integer size; + lp_feasibility_set_int_size(full, size.get_internal()); + size_t size_approx = lp_feasibility_set_int_size_approx(full); + CHECK(size == prime); + CHECK(size_approx == prime); + } + SUBCASE("feasibility_set_int::size::empty") { + Integer size; + lp_feasibility_set_int_size(empty, size.get_internal()); + size_t size_approx = lp_feasibility_set_int_size_approx(empty); + CHECK(size == 0); + CHECK(size_approx == 0); + } + + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + lp_feasibility_set_int_delete(full); + lp_feasibility_set_int_delete(empty); + + CHECK_EQ(K.get_internal()->ref_count, 1); +} + +TEST_CASE("feasibility_set_int::pick") { + Integer prime(7); + IntegerRing K(prime, true); + + lp_feasibility_set_int_t *set1 = integer_list({1,2}).gen_set(K, false); + lp_feasibility_set_int_t *set1_inv = integer_list({0,3,4,5,6}).gen_set(K, true); + + SUBCASE("feasibility_set_int::pick") { + Integer pick; + lp_feasibility_set_int_pick_value(set1, pick.get_internal()); + + CHECK((pick == 1 || pick == 2)); + } + SUBCASE("feasibility_set_int::pick_inv") { + Integer pick; + lp_feasibility_set_int_pick_value(set1_inv, pick.get_internal()); + CHECK((pick == 1 || pick == 2)); + } + + lp_feasibility_set_int_delete(set1); + lp_feasibility_set_int_delete(set1_inv); + + CHECK_EQ(K.get_internal()->ref_count, 1); +}