From 47464abf08fa748f1facbe46c133a8a61449b4ae Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Fri, 20 Sep 2024 19:20:26 +0200 Subject: [PATCH 1/8] code readability update --- src/polynomial/feasibility_set.c | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/polynomial/feasibility_set.c b/src/polynomial/feasibility_set.c index 79d919e..adb3873 100644 --- a/src/polynomial/feasibility_set.c +++ b/src/polynomial/feasibility_set.c @@ -457,7 +457,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; } From 84e85a387684fa987204776629fb40825d0f9ee6 Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Fri, 20 Sep 2024 20:42:57 +0200 Subject: [PATCH 2/8] Code cleanup --- include/interval.h | 4 ++-- src/interval/interval.c | 6 +++--- src/number/value.c | 3 +-- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/include/interval.h b/include/interval.h index f77c783..728b633 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). */ @@ -113,7 +113,7 @@ 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); /** diff --git a/src/interval/interval.c b/src/interval/interval.c index 85594f7..5ada085 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; } diff --git a/src/number/value.c b/src/number/value.c index cfd4d5f..62f45fc 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; From 61438d1ed0f8751e863b729bdd1d4d365bbea60f Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Sat, 21 Sep 2024 14:11:25 +0200 Subject: [PATCH 3/8] Implemented has_integer for interval and FeasibleSet --- include/feasibility_set.h | 5 ++++ include/interval.h | 3 ++ include/polyxx/interval.h | 3 ++ python/polypyFeasibilitySet2.c | 45 ++++++++++++++++++++++++++-- python/polypyFeasibilitySet3.c | 51 ++++++++++++++++++++++++++++---- python/polypyInterval2.c | 14 +++++++-- python/polypyInterval3.c | 17 ++++++++--- src/interval/interval.c | 28 ++++++++++++++++++ src/polynomial/feasibility_set.c | 9 ++++++ src/polyxx/interval.cpp | 5 ++++ test/polyxx/test_interval.cpp | 11 +++++++ 11 files changed, 177 insertions(+), 14 deletions(-) diff --git a/include/feasibility_set.h b/include/feasibility_set.h index c0d9d76..d4a6271 100644 --- a/include/feasibility_set.h +++ b/include/feasibility_set.h @@ -98,6 +98,11 @@ int lp_feasibility_set_is_point(const lp_feasibility_set_t* 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/interval.h b/include/interval.h index 728b633..3c1a828 100644 --- a/include/interval.h +++ b/include/interval.h @@ -83,6 +83,9 @@ 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); + /** Returns an approximation of the log interval size */ int lp_interval_size_approx(const lp_interval_t* I); diff --git a/include/polyxx/interval.h b/include/polyxx/interval.h index 432cbdc..2b78815 100644 --- a/include/polyxx/interval.h +++ b/include/polyxx/interval.h @@ -81,6 +81,9 @@ 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); + /** 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 1958901..3705d8c 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 e5ecf39..979d4f1 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 aa6d90e..8dbdca0 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 260bcc8..adf34c1 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 5ada085..9d9e622 100644 --- a/src/interval/interval.c +++ b/src/interval/interval.c @@ -686,6 +686,34 @@ int lp_interval_contains(const lp_interval_t* I, const lp_value_t* v) { return 1; } +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; +} 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/polynomial/feasibility_set.c b/src/polynomial/feasibility_set.c index adb3873..206670c 100644 --- a/src/polynomial/feasibility_set.c +++ b/src/polynomial/feasibility_set.c @@ -179,6 +179,15 @@ int lp_feasibility_set_contains(const lp_feasibility_set_t* set, const lp_value_ 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_int(set->intervals + i)) { + return 1; + } + } + return 0; +} + // We get smallest integer < rational < algebraic // If same we get one with largest interval size static inline diff --git a/src/polyxx/interval.cpp b/src/polyxx/interval.cpp index f5dd72b..6eaa4fe 100644 --- a/src/polyxx/interval.cpp +++ b/src/polyxx/interval.cpp @@ -68,6 +68,11 @@ 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()); + } + int log_size(const Interval& i) { return lp_interval_size_approx(i.get_internal()); } diff --git a/test/polyxx/test_interval.cpp b/test/polyxx/test_interval.cpp index 8fbdc54..b7b5e5c 100644 --- a/test/polyxx/test_interval.cpp +++ b/test/polyxx/test_interval.cpp @@ -53,6 +53,17 @@ 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::pick_value") { Interval i(1,3); CHECK(contains(i, pick_value(i))); From ece4f77704aa8bace6f8896ccde1ad3972539ee0 Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Mon, 23 Sep 2024 09:01:56 +0200 Subject: [PATCH 4/8] Created lp_integer_fits_int --- include/integer.h | 5 +++++ src/number/integer.c | 4 ++++ src/number/integer.h | 5 +++++ src/upolynomial/factorization.c | 10 +++++----- src/upolynomial/root_finding.c | 2 +- 5 files changed, 20 insertions(+), 6 deletions(-) diff --git a/include/integer.h b/include/integer.h index affcca5..efc49df 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/src/number/integer.c b/src/number/integer.c index f0f8788..e89ecef 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 dc0d664..d8e2dc8 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/upolynomial/factorization.c b/src/upolynomial/factorization.c index 44d9efd..e165d6d 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 c64d1b5..e02ba9c 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); From 13bc940658ef343a09ed24e70de1886e43fd5047 Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Mon, 23 Sep 2024 09:09:24 +0200 Subject: [PATCH 5/8] Created lp_feasibility_set_count_int --- include/feasibility_set.h | 10 ++++++++ include/interval.h | 3 +++ include/polyxx/interval.h | 3 ++- src/interval/interval.c | 40 ++++++++++++++++++++++++++++++++ src/polynomial/feasibility_set.c | 30 ++++++++++++++++++++++-- src/polyxx/interval.cpp | 4 ++++ test/polyxx/test_interval.cpp | 15 ++++++++++++ 7 files changed, 102 insertions(+), 3 deletions(-) diff --git a/include/feasibility_set.h b/include/feasibility_set.h index d4a6271..6afa343 100644 --- a/include/feasibility_set.h +++ b/include/feasibility_set.h @@ -93,6 +93,16 @@ 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. */ diff --git a/include/interval.h b/include/interval.h index 3c1a828..e1d4642 100644 --- a/include/interval.h +++ b/include/interval.h @@ -86,6 +86,9 @@ 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); diff --git a/include/polyxx/interval.h b/include/polyxx/interval.h index 2b78815..66ddb4b 100644 --- a/include/polyxx/interval.h +++ b/include/polyxx/interval.h @@ -83,7 +83,8 @@ namespace poly { 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/src/interval/interval.c b/src/interval/interval.c index 9d9e622..2009f11 100644 --- a/src/interval/interval.c +++ b/src/interval/interval.c @@ -715,6 +715,46 @@ int lp_interval_contains_int(const lp_interval_t* I) { 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); if (I->is_point) { diff --git a/src/polynomial/feasibility_set.c b/src/polynomial/feasibility_set.c index 206670c..cf1d9bf 100644 --- a/src/polynomial/feasibility_set.c +++ b/src/polynomial/feasibility_set.c @@ -144,11 +144,37 @@ 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); + 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, ", "); } diff --git a/src/polyxx/interval.cpp b/src/polyxx/interval.cpp index 6eaa4fe..867815e 100644 --- a/src/polyxx/interval.cpp +++ b/src/polyxx/interval.cpp @@ -73,6 +73,10 @@ namespace poly { 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/test/polyxx/test_interval.cpp b/test/polyxx/test_interval.cpp index b7b5e5c..016bc59 100644 --- a/test/polyxx/test_interval.cpp +++ b/test/polyxx/test_interval.cpp @@ -64,6 +64,21 @@ TEST_CASE("interval::contains_int") { 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))); From 43aa5f5f3abd04775a04e6067d570e797e8854be Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Mon, 23 Sep 2024 15:35:50 +0200 Subject: [PATCH 6/8] Implement binary search for feasibility_set::contains --- include/interval.h | 3 +++ src/interval/interval.c | 18 +++++++++++------- src/polynomial/feasibility_set.c | 15 ++++++++++++--- 3 files changed, 26 insertions(+), 10 deletions(-) diff --git a/include/interval.h b/include/interval.h index e1d4642..f843b5b 100644 --- a/include/interval.h +++ b/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. */ diff --git a/src/interval/interval.c b/src/interval/interval.c index 2009f11..dce080e 100644 --- a/src/interval/interval.c +++ b/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) { diff --git a/src/polynomial/feasibility_set.c b/src/polynomial/feasibility_set.c index cf1d9bf..a3099b4 100644 --- a/src/polynomial/feasibility_set.c +++ b/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; } } From 31839cb9167b5652a80d1b852fc738b922a73526 Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Mon, 7 Oct 2024 15:50:02 +0200 Subject: [PATCH 7/8] fixed whitespaces --- test/polyxx/test_interval.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/polyxx/test_interval.cpp b/test/polyxx/test_interval.cpp index 016bc59..d871099 100644 --- a/test/polyxx/test_interval.cpp +++ b/test/polyxx/test_interval.cpp @@ -70,13 +70,13 @@ TEST_CASE("interval::count_int") { 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(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); + CHECK(count_int(Interval(LONG_MIN,LONG_MAX)) == LONG_MAX); } TEST_CASE("interval::pick_value") { From db09ece0677a3c252ae527d729278608bace67fb Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Wed, 9 Oct 2024 15:55:39 +0200 Subject: [PATCH 8/8] Added comment --- src/polynomial/feasibility_set.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/polynomial/feasibility_set.c b/src/polynomial/feasibility_set.c index a3099b4..7e59ab0 100644 --- a/src/polynomial/feasibility_set.c +++ b/src/polynomial/feasibility_set.c @@ -149,6 +149,7 @@ int lp_feasibility_set_is_point_int(const lp_feasibility_set_t* set) { 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; }