diff --git a/.ci/build.sh b/.ci/build.sh index 446ab1ad..e8832ebe 100644 --- a/.ci/build.sh +++ b/.ci/build.sh @@ -2,7 +2,7 @@ mkdir -p build || return 1 cd build/ || return 1 -cmake -D DEVELOPER=ON -D USE_BLISS=ON -D USE_COCOA=ON ../ || return 1 +cmake -D DEVELOPER=ON -D USE_BLISS=ON -D USE_COCOA=ON -D USE_LIBPOLY=ON ../ || return 1 function keep_waiting() { while true; do @@ -61,11 +61,11 @@ elif [[ ${TASK} == "documentation" ]]; then git push -f origin master || return 1 elif [[ ${TASK} == "tidy" ]]; then - /usr/bin/time make ${MAKE_PARALLEL} tidy || return 1 -elif [[ ${TASK} == "parallel" ]]; then +elif [[ ${TASK} == "all" ]]; then /usr/bin/time make ${MAKE_PARALLEL} || return 1 + else /usr/bin/time make || return 1 fi \ No newline at end of file diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d19af612..1df3a918 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,68 +12,71 @@ cache: #policy: pull-push stages: - - build-gcc - - build-clang + - build - test - documentation build-gcc12: dependencies: [] - stage: build-gcc + stage: build script: - export CC=/usr/bin/gcc-12 && export CXX=/usr/bin/g++-12 - - MAKE_PARALLEL=-j8 TASK=parallel source .ci/build.sh + - MAKE_PARALLEL=-j8 TASK=all source .ci/build.sh artifacts: name: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG" paths: - build/ - cmake/ - src/ + expire_in: 1 week -build-gcc11: +.build-gcc11: dependencies: [] - stage: build-gcc + stage: build script: - export CC=/usr/bin/gcc-11 && export CXX=/usr/bin/g++-11 - - MAKE_PARALLEL=-j8 TASK=parallel source .ci/build.sh + - MAKE_PARALLEL=-j8 TASK=all source .ci/build.sh artifacts: name: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG" paths: - build/ - cmake/ - src/ + expire_in: 1 week only: - development build-clang14: dependencies: [] - stage: build-clang + stage: build script: - export CC=/usr/bin/clang-14 && export CXX=/usr/bin/clang++-14 - - MAKE_PARALLEL=-j8 TASK=parallel source .ci/build.sh + - MAKE_PARALLEL=-j8 TASK=all source .ci/build.sh artifacts: name: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG" paths: - build/ - cmake/ - src/ + expire_in: 1 week -build-clang13: +.build-clang13: dependencies: [] - stage: build-clang + stage: build script: - export CC=/usr/bin/clang-13 && export CXX=/usr/bin/clang++-13 - - MAKE_PARALLEL=-j8 TASK=parallel source .ci/build.sh + - MAKE_PARALLEL=-j8 TASK=all source .ci/build.sh artifacts: name: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG" paths: - build/ - cmake/ - src/ + expire_in: 1 week only: - development -test-clang: +.test-clang: dependencies: [build-clang14] stage: test script: diff --git a/CMakeLists.txt b/CMakeLists.txt index 06a573a5..45e14164 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,7 +14,7 @@ include(carlmacros) set(PROJECT_FULLNAME "carl") set(PROJECT_DESCRIPTION "Computer ARithmetic Library") -set_version(24 02) +set_version(24 04) message(STATUS "Version: ${PROJECT_FULLNAME} ${PROJECT_VERSION_FULL}") # # # # # # # # # # # # # # # # # # # # # # diff --git a/cmake/FindPoly.cmake b/cmake/FindPoly.cmake index 7fd4b012..0c451da3 100644 --- a/cmake/FindPoly.cmake +++ b/cmake/FindPoly.cmake @@ -2,13 +2,13 @@ if(NOT FORCE_SHIPPED_RESOURCES) find_path(LIBPOLY_INCLUDE_DIR NAMES poly/poly.h) find_library(LIBPOLY_SHARED NAMES poly${DYNAMIC_EXT}) find_library(LIBPOLY_STATIC NAMES poly${STATIC_EXT}) - find_library(LIBPOLYXX_SHARED NAMES polyxx${DYNAMIC_EXT}) - find_library(LIBPOLYXX_STATIC NAMES polyxx${STATIC_EXT}) + # find_library(LIBPOLYXX_SHARED NAMES polyxx${DYNAMIC_EXT}) + # find_library(LIBPOLYXX_STATIC NAMES polyxx${STATIC_EXT}) endif() set(LIBPOLY_FOUND_SYSTEM FALSE) -if(LIBPOLY_INCLUDE_DIR AND LIBPOLY_SHARED AND LIBPOLY_STATIC AND LIBPOLYXX_SHARED AND LIBPOLYXX_STATIC) +if(LIBPOLY_INCLUDE_DIR AND LIBPOLY_SHARED AND LIBPOLY_STATIC) # AND LIBPOLYXX_SHARED AND LIBPOLYXX_STATIC) #message(STATUS "Found libpoly: ${LIBPOLY_INCLUDE_DIR} ${LIBPOLY_SHARED} ${LIBPOLY_STATIC} ${LIBPOLYXX_SHARED} ${LIBPOLYXX_STATIC}") #set(LIBPOLY_FOUND_SYSTEM TRUE) #TODO: Somehow this does not work, so we use the shipped version for now @@ -21,8 +21,11 @@ if(NOT LIBPOLY_FOUND_SYSTEM) ExternalProject_Add( LIBPOLY-EP GIT_REPOSITORY https://github.com/SRI-CSL/libpoly + GIT_TAG v0.1.13 PATCH_COMMAND git reset --hard - COMMAND git apply ${CMAKE_SOURCE_DIR}/cmake/patches/libpoly_variable_db.patch + # UPDATE_COMMAND "" + COMMAND git apply ${CMAKE_SOURCE_DIR}/cmake/patches/libpoly.patch + # PATCH_COMMAND git apply ${CMAKE_SOURCE_DIR}/cmake/patches/libpoly.patch CMAKE_ARGS -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -DCMAKE_INSTALL_PREFIX= -DLIBPOLY_BUILD_PYTHON_API=OFF @@ -31,13 +34,13 @@ if(NOT LIBPOLY_FOUND_SYSTEM) -DGMP_INCLUDE_DIR=${GMP_LIBRARY_DIR} -DGMP_LIBRARY=${GMP_LIB} BUILD_COMMAND ${CMAKE_COMMAND} - COMMAND ${CMAKE_MAKE_PROGRAM} poly polyxx static_poly static_polyxx static_pic_poly static_pic_polyxx + COMMAND ${CMAKE_MAKE_PROGRAM} poly static_poly static_pic_poly polyxx static_polyxx static_pic_polyxx COMMAND ${CMAKE_MAKE_PROGRAM} install BUILD_BYPRODUCTS /lib/libpoly${STATIC_EXT} /lib/libpoly${DYNAMIC_EXT} - /lib/libpolyxx${STATIC_EXT} - /lib/libpolyxx${DYNAMIC_EXT} + # /lib/libpolyxx${STATIC_EXT} + # /lib/libpolyxx${DYNAMIC_EXT} ) # ExternalProject_Add_Step( # LIBPOLY-EP cleanup @@ -53,8 +56,8 @@ if(NOT LIBPOLY_FOUND_SYSTEM) set(LIBPOLY_SHARED "${CMAKE_BINARY_DIR}/resources/lib/libpoly${DYNAMIC_EXT}") set(LIBPOLY_STATIC "${CMAKE_BINARY_DIR}/resources/lib/libpoly${STATIC_EXT}") - set(LIBPOLYXX_SHARED "${CMAKE_BINARY_DIR}/resources/lib/libpolyxx${DYNAMIC_EXT}") - set(LIBPOLYXX_STATIC "${CMAKE_BINARY_DIR}/resources/lib/libpolyxx${STATIC_EXT}") + # set(LIBPOLYXX_SHARED "${CMAKE_BINARY_DIR}/resources/lib/libpolyxx${DYNAMIC_EXT}") + # set(LIBPOLYXX_STATIC "${CMAKE_BINARY_DIR}/resources/lib/libpolyxx${STATIC_EXT}") add_dependencies(LIBPOLY-EP GMP_SHARED GMPXX_SHARED) @@ -78,34 +81,35 @@ set_target_properties(LIBPOLY_STATIC PROPERTIES INTERFACE_LINK_LIBRARIES "${GMP_STATIC}" ) -add_library(LIBPOLYXX_SHARED SHARED IMPORTED GLOBAL) -target_link_libraries(LIBPOLY_STATIC INTERFACE GMPXX_SHARED) -set_target_properties(LIBPOLYXX_SHARED PROPERTIES - IMPORTED_LOCATION ${LIBPOLYXX_SHARED} - INTERFACE_INCLUDE_DIRECTORIES ${LIBPOLY_INCLUDE_DIR} - INTERFACE_LINK_LIBRARIES "${LIBPOLY_SHARED}" -) - -add_library(LIBPOLYXX_STATIC STATIC IMPORTED GLOBAL) -target_link_libraries(LIBPOLY_STATIC INTERFACE GMPXX_STATIC) -set_target_properties(LIBPOLYXX_STATIC PROPERTIES - IMPORTED_LOCATION ${LIBPOLYXX_STATIC} - INTERFACE_INCLUDE_DIRECTORIES ${LIBPOLY_INCLUDE_DIR} - INTERFACE_LINK_LIBRARIES "${LIBPOLY_STATIC}" -) +# add_library(LIBPOLYXX_SHARED SHARED IMPORTED GLOBAL) +# target_link_libraries(LIBPOLY_STATIC INTERFACE GMPXX_SHARED) +# set_target_properties(LIBPOLYXX_SHARED PROPERTIES +# IMPORTED_LOCATION ${LIBPOLYXX_SHARED} +# INTERFACE_INCLUDE_DIRECTORIES ${LIBPOLY_INCLUDE_DIR} +# INTERFACE_LINK_LIBRARIES "${LIBPOLY_SHARED}" +# ) +# +# add_library(LIBPOLYXX_STATIC STATIC IMPORTED GLOBAL) +# target_link_libraries(LIBPOLY_STATIC INTERFACE GMPXX_STATIC) +# set_target_properties(LIBPOLYXX_STATIC PROPERTIES +# IMPORTED_LOCATION ${LIBPOLYXX_STATIC} +# INTERFACE_INCLUDE_DIRECTORIES ${LIBPOLY_INCLUDE_DIR} +# INTERFACE_LINK_LIBRARIES "${LIBPOLY_STATIC}" +# ) add_dependencies(LIBPOLY_SHARED GMP_SHARED) add_dependencies(LIBPOLY_STATIC GMP_STATIC) -add_dependencies(LIBPOLYXX_SHARED GMPXX_SHARED) -add_dependencies(LIBPOLYXX_STATIC GMPXX_STATIC) +# add_dependencies(LIBPOLYXX_SHARED GMPXX_SHARED) +# add_dependencies(LIBPOLYXX_STATIC GMPXX_STATIC) if(NOT LIBPOLY_FOUND_SYSTEM) add_dependencies(LIBPOLY_SHARED LIBPOLY-EP) add_dependencies(LIBPOLY_STATIC LIBPOLY-EP) - add_dependencies(LIBPOLYXX_SHARED LIBPOLY-EP) - add_dependencies(LIBPOLYXX_STATIC LIBPOLY-EP) +# add_dependencies(LIBPOLYXX_SHARED LIBPOLY-EP) +# add_dependencies(LIBPOLYXX_STATIC LIBPOLY-EP) endif() set(LIBPOLY_FOUND TRUE) -mark_as_advanced(LIBPOLY_INCLUDE_DIR LIBPOLY_SHARED LIBPOLY_STATIC LIBPOLYXX_SHARED LIBPOLYXX_STATIC LIBPOLY_FOUND) \ No newline at end of file +# mark_as_advanced(LIBPOLY_INCLUDE_DIR LIBPOLY_SHARED LIBPOLY_STATIC LIBPOLYXX_SHARED LIBPOLYXX_STATIC LIBPOLY_FOUND) +mark_as_advanced(LIBPOLY_INCLUDE_DIR LIBPOLY_SHARED LIBPOLY_STATIC LIBPOLY_FOUND) \ No newline at end of file diff --git a/cmake/export.cmake b/cmake/export.cmake index de28ac0e..32e6d376 100644 --- a/cmake/export.cmake +++ b/cmake/export.cmake @@ -31,8 +31,8 @@ endif() if(USE_LIBPOLY) export_target(DEP_TARGETS LIBPOLY_SHARED GMPXX_SHARED GMP_SHARED) export_target(DEP_TARGETS LIBPOLY_STATIC GMPXX_STATIC GMP_STATIC) - export_target(DEP_TARGETS LIBPOLYXX_SHARED GMPXX_SHARED GMP_SHARED) - export_target(DEP_TARGETS LIBPOLYXX_STATIC GMPXX_STATIC GMP_STATIC) + # export_target(DEP_TARGETS LIBPOLYXX_SHARED GMPXX_SHARED GMP_SHARED) + # export_target(DEP_TARGETS LIBPOLYXX_STATIC GMPXX_STATIC GMP_STATIC) endif() export_target(DEP_TARGETS GMP_SHARED) export_target(DEP_TARGETS GMP_STATIC) diff --git a/cmake/patches/libpoly.patch b/cmake/patches/libpoly.patch new file mode 100644 index 00000000..76407d2d --- /dev/null +++ b/cmake/patches/libpoly.patch @@ -0,0 +1,103 @@ +diff --git a/include/polynomial.h b/include/polynomial.h +index 8a9533f..2905513 100644 +--- a/include/polynomial.h ++++ b/include/polynomial.h +@@ -109,6 +109,8 @@ int lp_polynomial_lc_sgn(const lp_polynomial_t* A); + /** Get the context of the given polynomial */ + const lp_polynomial_context_t* lp_polynomial_get_context(const lp_polynomial_t* A); + ++void lp_polynomial_set_context(lp_polynomial_t* A, const lp_polynomial_context_t* ctx); ++ + /** Returns all the variables (it will not clear the output list vars) */ + void lp_polynomial_get_variables(const lp_polynomial_t* A, lp_variable_list_t* vars); + +@@ -337,6 +339,8 @@ lp_polynomial_t* lp_polynomial_constraint_explain_infer_bounds(const lp_polynomi + */ + int lp_polynomial_constraint_evaluate(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, const lp_assignment_t* M); + ++int lp_polynomial_constraint_evaluate_subs(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, const lp_assignment_t* M); ++ + /** + * Given a polynomial A(x1, ..., xn, y) with y being the top variable, a root index, + * a sign condition, and an assignment M that assigns x1, ..., xn, the function +diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c +index 11d948f..7468ab6 100644 +--- a/src/polynomial/polynomial.c ++++ b/src/polynomial/polynomial.c +@@ -1031,6 +1031,18 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t + lp_value_construct_none(&x_value_backup); + } + ++ lp_polynomial_t* B; ++ { ++ coefficient_t A_rat; ++ lp_integer_t multiplier; ++ integer_construct(&multiplier); ++ coefficient_construct(A->ctx, &A_rat); ++ coefficient_evaluate_rationals(A->ctx, &A->data, M, &A_rat, &multiplier); ++ B = lp_polynomial_new_from_coefficient(A->ctx, &A_rat); ++ coefficient_destruct(&A_rat); ++ integer_destruct(&multiplier); ++ } ++ + size_t i; + + lp_polynomial_t** factors = 0; +@@ -1044,11 +1056,13 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t + // Get the reduced polynomial + lp_polynomial_t A_r; + lp_polynomial_construct(&A_r, A->ctx); +- lp_polynomial_reductum_m(&A_r, A, M); +- assert(x == lp_polynomial_top_variable(A)); +- +- // Get the square-free factorization +- lp_polynomial_factor_square_free(&A_r, &factors, &multiplicities, &factors_size); ++ if (x == lp_polynomial_top_variable(B)) { ++ lp_polynomial_reductum_m(&A_r, B, M); ++ assert(x == lp_polynomial_top_variable(B)); ++ // Get the square-free factorization ++ lp_polynomial_factor_square_free(&A_r, &factors, &multiplicities, &factors_size); ++ } ++ lp_polynomial_delete(B); + + // Count the max number of roots + size_t total_degree = 0; +@@ -1943,6 +1957,22 @@ int lp_polynomial_constraint_evaluate(const lp_polynomial_t* A, lp_sign_conditio + return lp_sign_condition_consistent(sgn_condition, p_sign); + } + ++int lp_polynomial_constraint_evaluate_subs(const lp_polynomial_t* A, lp_sign_condition_t sgn_condition, const lp_assignment_t* M) { ++ coefficient_t A_rat; ++ lp_integer_t multiplier; ++ integer_construct(&multiplier); ++ coefficient_construct(A->ctx, &A_rat); ++ coefficient_evaluate_rationals(A->ctx, &A->data, M, &A_rat, &multiplier); ++ integer_destruct(&multiplier); ++ int res = -1; ++ if (A_rat.type == COEFFICIENT_NUMERIC) { ++ int sgn = integer_sgn(lp_Z, &A_rat.value.num); ++ res = lp_sign_condition_consistent(sgn_condition, sgn); ++ } ++ coefficient_destruct(&A_rat); ++ return res; ++} ++ + int lp_polynomial_root_constraint_evaluate(const lp_polynomial_t* A, size_t root_index, lp_sign_condition_t sgn_condition, const lp_assignment_t* M) { + + int eval; +diff --git a/src/variable/variable_db.c b/src/variable/variable_db.c +index 60f3df4..33866c4 100644 +--- a/src/variable/variable_db.c ++++ b/src/variable/variable_db.c +@@ -63,9 +63,11 @@ void lp_variable_db_add_variable(lp_variable_db_t* var_db, lp_variable_t var, co + assert(var_db); + while (var >= var_db->capacity) { + lp_variable_db_resize(var_db, 2*var_db->capacity); ++ var_db->size = var_db->capacity < var ? var_db->capacity : var; + } + assert(var_db->variable_names[var] == 0); + var_db->variable_names[var] = strdup(name); ++ var_db->size = var_db->size < var ? var : var_db->size; + } + + void lp_variable_db_construct(lp_variable_db_t* var_db) { diff --git a/cmake/patches/libpoly_variable_db.patch b/cmake/patches/libpoly_variable_db.patch deleted file mode 100644 index 09344ebc..00000000 --- a/cmake/patches/libpoly_variable_db.patch +++ /dev/null @@ -1,16 +0,0 @@ -diff --git a/src/variable/variable_db.c b/src/variable/variable_db.c -index 60f3df4..33866c4 100644 ---- a/src/variable/variable_db.c -+++ b/src/variable/variable_db.c -@@ -63,9 +63,11 @@ void lp_variable_db_add_variable(lp_variable_db_t* var_db, lp_variable_t var, co - assert(var_db); - while (var >= var_db->capacity) { - lp_variable_db_resize(var_db, 2*var_db->capacity); -+ var_db->size = var_db->capacity < var ? var_db->capacity : var; - } - assert(var_db->variable_names[var] == 0); - var_db->variable_names[var] = strdup(name); -+ var_db->size = var_db->size < var ? var : var_db->size; - } - - void lp_variable_db_construct(lp_variable_db_t* var_db) { diff --git a/resources/resources.cmake b/resources/resources.cmake index 4f58bcb1..2e261729 100644 --- a/resources/resources.cmake +++ b/resources/resources.cmake @@ -170,9 +170,9 @@ endif() IF(USE_LIBPOLY) set(LIBPOLY_VERSION "0.1.11") find_package(Poly ${LIBPOLY_VERSION} REQUIRED QUIET) - print_resource_info("LibPoly" LIBPOLYXX_STATIC ${LIBPOLY_VERSION}) + print_resource_info("LibPoly" LIBPOLY_STATIC ${LIBPOLY_VERSION}) add_dependencies(resources LIBPOLY_SHARED LIBPOLY_STATIC) - add_dependencies(resources LIBPOLYXX_SHARED LIBPOLYXX_STATIC) + # add_dependencies(resources LIBPOLYXX_SHARED LIBPOLYXX_STATIC) else() message(STATUS "LibPoly is disabled") endif() diff --git a/src/carl-arith/CMakeLists.txt b/src/carl-arith/CMakeLists.txt index c12f792c..ac4a8c03 100644 --- a/src/carl-arith/CMakeLists.txt +++ b/src/carl-arith/CMakeLists.txt @@ -30,9 +30,9 @@ if(USE_MPFR_FLOAT) target_link_libraries(carl-arith-static PUBLIC MPFR_STATIC) endif() if(USE_LIBPOLY) -target_include_dirs_from(carl-arith-objects SYSTEM PUBLIC LIBPOLYXX_SHARED LIBPOLYXX_STATIC LIBPOLY_SHARED LIBPOLY_STATIC) -target_link_libraries(carl-arith-shared PUBLIC LIBPOLYXX_SHARED LIBPOLY_SHARED) -target_link_libraries(carl-arith-static PUBLIC LIBPOLYXX_STATIC LIBPOLY_STATIC) +target_include_dirs_from(carl-arith-objects SYSTEM PUBLIC LIBPOLY_SHARED LIBPOLY_STATIC) # LIBPOLYXX_SHARED LIBPOLYXX_STATIC +target_link_libraries(carl-arith-shared PUBLIC LIBPOLY_SHARED) # LIBPOLYXX_SHARED +target_link_libraries(carl-arith-static PUBLIC LIBPOLY_STATIC) # LIBPOLYXX_STATIC endif() include(install) diff --git a/src/carl-arith/extended/VariableComparison.h b/src/carl-arith/extended/VariableComparison.h index f9540861..d1c528a9 100644 --- a/src/carl-arith/extended/VariableComparison.h +++ b/src/carl-arith/extended/VariableComparison.h @@ -70,6 +70,9 @@ namespace carl { const std::variant& value() const { return m_value; } + std::variant& value() { + return m_value; + } bool is_equality() const { return negated() ? relation() == Relation::NEQ : relation() == Relation::EQ; } @@ -153,15 +156,26 @@ namespace carl { } template - boost::tribool evaluate(const VariableComparison& f, const Assignment::RAN>& a) { + boost::tribool evaluate(const VariableComparison& f, const Assignment::RAN>& a, bool evaluate_non_welldef = false) { typename VariableComparison::RAN lhs; + if (a.find(f.var()) == a.end()) return boost::indeterminate; typename VariableComparison::RAN rhs = a.at(f.var()); if (std::holds_alternative::RAN>(f.value())) { lhs = std::get::RAN>(f.value()); } else { - auto res = carl::evaluate(std::get::MR>(f.value()), a); - if (!res) return boost::indeterminate; - else lhs = *res; + if (!evaluate_non_welldef) { + auto res = carl::evaluate(std::get::MR>(f.value()), a); + if (!res) return boost::indeterminate; + else lhs = *res; + } else { + auto vars = carl::variables(std::get::MR>(f.value())); + for (const auto& v : vars) { + if (a.find(v) == a.end()) return boost::indeterminate; + } + auto res = carl::evaluate(std::get::MR>(f.value()), a); + if (!res) return f.negated(); + else lhs = *res; + } } if (!f.negated()) return evaluate(rhs, f.relation(), lhs); else return !evaluate(rhs, f.relation(), lhs); diff --git a/src/carl-arith/poly/Conversion.h b/src/carl-arith/poly/Conversion.h index f8f5d8e7..015829c0 100644 --- a/src/carl-arith/poly/Conversion.h +++ b/src/carl-arith/poly/Conversion.h @@ -43,37 +43,24 @@ struct ConvertHelper> { if (denominator < 0) { denominator *= -1; } - /* - // LPPolynomial can only have integer coefficients -> so we have to store the lcm of the coefficients of every term - auto coprimeFactor = p.main_denom(); - mpz_class denominator; - if (carl::get_denom(coprimeFactor) != 1) { - // if coprime factor is not an integer - denominator = coprimeFactor > 0 ? mpz_class(1) : mpz_class(-1); - } else { - denominator = mpz_class(coprimeFactor); - } - */ + CARL_LOG_DEBUG("carl.converter", "Coprime Factor/ Denominator: " << denominator); - LPPolynomial res(context); - // iterate over terms + lp_polynomial_t* res = lp_polynomial_new(context.lp_context()); for (const auto& term : p) { - // Multiply by denominator to make the coefficient an integer assert(carl::get_denom(term.coeff() * denominator) == 1); - LPPolynomial t(context, mpz_class(term.coeff() * denominator)); - // iterate over monomial + lp_monomial_t t; + lp_monomial_construct(context.lp_context(), &t); + lp_monomial_set_coefficient(context.lp_context(), &t, (lp_integer_t*)mpz_class(term.coeff() * denominator).get_mpz_t()); if (term.monomial()) { // A monomial is a vector of pairs for (const std::pair& var_pow : (*term.monomial()).exponents()) { - // multiply 1*var^pow - t *=LPPolynomial(context, var_pow.first, mpz_class(1), (unsigned)var_pow.second); + lp_monomial_push(&t, context.lp_variable(var_pow.first), var_pow.second); } } - CARL_LOG_TRACE("carl.converter", "converted term: " << term << " -> " << t); - res += t; + lp_polynomial_add_monomial(res, &t); + lp_monomial_destruct(&t); } - CARL_LOG_DEBUG("carl.converter", "Got Polynomial: " << res); - return res; + return LPPolynomial(res, context); } }; diff --git a/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h b/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h index 4d7b9be1..347be59e 100644 --- a/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h +++ b/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h @@ -55,6 +55,9 @@ class CoCoAAdaptorLP { static CoCoA::BigInt convert(const mpz_class& n) { return CoCoA::BigIntFromMPZ(n.get_mpz_t()); } + static CoCoA::BigInt convert(const mpz_t n) { + return CoCoA::BigIntFromMPZ(n); + } mpz_class convert(const CoCoA::BigInt& n) const { return mpz_class(CoCoA::mpzref(n)); } @@ -108,9 +111,7 @@ class CoCoAAdaptorLP { } } - auto coefPol = poly::Integer(&(m->a)); - mpz_class* coef = poly::detail::cast_to_gmp(&coefPol); - *data->resPoly += CoCoA::monomial(data->mRing, convert(*(coef)), exponents); + *data->resPoly += CoCoA::monomial(data->mRing, convert(&(m->a)), exponents); }; DataLP data{&res, mSymbolThere, mSymbolBack, mQ, mRing, mContext}; @@ -119,27 +120,31 @@ class CoCoAAdaptorLP { } LPPolynomial convert(const CoCoA::RingElem& p) const { - poly::Polynomial temPoly(mContext.lp_context()); + lp_polynomial_t* res = lp_polynomial_new(mContext.lp_context()); for (CoCoA::SparsePolyIter i = CoCoA::BeginIter(p); !CoCoA::IsEnded(i); ++i) { mpq_class coeff; convert(coeff, CoCoA::coeff(i)); std::vector exponents; CoCoA::exponents(exponents, CoCoA::PP(i)); - poly::Polynomial termPoly = poly_helper::construct_poly(mContext.lp_context(), (poly::Integer)carl::get_num(coeff)); + + assert(coeff != 0); + assert(carl::get_denom(coeff) == 1); + + lp_monomial_t t; + lp_monomial_construct(mContext.lp_context(), &t); + lp_monomial_set_coefficient(mContext.lp_context(), &t, carl::get_num(coeff).get_mpz_t()); for (std::size_t i = 0; i < exponents.size(); ++i) { if (exponents[i] == 0) continue; - std::optional var = mContext.lp_variable(mSymbolBack[i]); - assert(var.has_value()); - poly::Polynomial polyVar = poly_helper::construct_poly(mContext.lp_context(), *var); - termPoly *= poly::pow(polyVar, (unsigned int)exponents[i]); + lp_variable_t var = mContext.lp_variable(mSymbolBack[i]); + lp_monomial_push(&t, var, (unsigned int)exponents[i]); } - temPoly += termPoly; + lp_polynomial_add_monomial(res, &t); + lp_monomial_destruct(&t); } - LPPolynomial res(mContext, temPoly); - return res; + return LPPolynomial(res, mContext);; } std::vector convert(const std::vector& p) const { diff --git a/src/carl-arith/poly/libpoly/Functions.cpp b/src/carl-arith/poly/libpoly/Functions.cpp new file mode 100644 index 00000000..2508a551 --- /dev/null +++ b/src/carl-arith/poly/libpoly/Functions.cpp @@ -0,0 +1,86 @@ +#include "Functions.h" + +#include +#ifdef USE_LIBPOLY + +#include "CoCoAAdaptorLP.h" + +namespace carl { + +inline Factors factorization(const LPPolynomial& p, const CoCoAAdaptorLP& adaptor, bool includeConstant) { + #ifndef USE_COCOA + CARL_LOG_FATAL("carl.poly", "factorization is not supported without libcocoa"); + #endif + return adaptor.factorize(p, includeConstant); +} + +/* + * @brief wrapper for the factorization of LPPolynomials + * @param LPPolynomial p + * @return factorization(p) + */ +Factors factorization(const LPPolynomial& p, bool includeConstant) { + CoCoAAdaptorLP adaptor = CoCoAAdaptorLP(p.context()); + return factorization(p, adaptor, includeConstant); +} + +inline std::vector irreducible_factors(const LPPolynomial& p, const CoCoAAdaptorLP& adaptor, bool includeConstants) { + #ifndef USE_COCOA + CARL_LOG_FATAL("carl.poly", "factorization is not supported without libcocoa"); + #endif + + if (is_constant(p)) { + if (includeConstants) { + return {p}; + } else { + return {}; + } + } else if (p.total_degree() == 1) { + return {p}; + } + + return adaptor.irreducible_factors(p, includeConstants); +} + +std::vector irreducible_factors(const LPPolynomial& p, bool includeConstants) { + CoCoAAdaptorLP adaptor = CoCoAAdaptorLP(p.context()); + return irreducible_factors(p, adaptor, includeConstants); +} + +std::vector square_free_factors(const LPPolynomial& p) { + lp_polynomial_t** factors = nullptr; + std::size_t* multiplicities = nullptr; + std::size_t size = 0; + lp_polynomial_factor_square_free(p.get_internal(), &factors, &multiplicities, &size); + std::vector res; + for (std::size_t i = 0; i < size; ++i) { + res.emplace_back(factors[i], p.context()); + } + free(factors); + free(multiplicities); + return res; +} + +std::vector content_free_factors(const LPPolynomial& p) { + lp_polynomial_t** factors = nullptr; + std::size_t* multiplicities = nullptr; + std::size_t size = 0; + lp_polynomial_factor_content_free(p.get_internal(), &factors, &multiplicities, &size); + std::vector res; + for (std::size_t i = 0; i < size; ++i) { + res.emplace_back(factors[i], p.context()); + } + free(factors); + free(multiplicities); + return res; +} + +std::vector groebner_basis(const std::vector& polys) { + if (polys.size() <= 1) return polys; + CoCoAAdaptorLP adaptor = CoCoAAdaptorLP(polys.at(0).context()); + return adaptor.GBasis(polys); +} + +} + +#endif \ No newline at end of file diff --git a/src/carl-arith/poly/libpoly/Functions.h b/src/carl-arith/poly/libpoly/Functions.h index be61f7e2..eb091385 100644 --- a/src/carl-arith/poly/libpoly/Functions.h +++ b/src/carl-arith/poly/libpoly/Functions.h @@ -5,8 +5,6 @@ #include "LPPolynomial.h" -#include "CoCoAAdaptorLP.h" - namespace carl { /* @@ -15,7 +13,9 @@ namespace carl { * @return gcd(p,q) */ inline LPPolynomial gcd(const LPPolynomial& p, const LPPolynomial& q) { - return LPPolynomial(p.context(), poly::gcd(p.get_polynomial(), q.get_polynomial())); + lp_polynomial_t* res = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_gcd(res, p.get_internal(), q.get_internal()); + return LPPolynomial(res, p.context()); } /* @@ -24,7 +24,9 @@ inline LPPolynomial gcd(const LPPolynomial& p, const LPPolynomial& q) { * @return lcm(p,q) */ inline LPPolynomial lcm(const LPPolynomial& p, const LPPolynomial& q) { - return LPPolynomial(p.context(), poly::lcm(p.get_polynomial(), q.get_polynomial())); + lp_polynomial_t* res = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_lcm(res, p.get_internal(), q.get_internal()); + return LPPolynomial(res, p.context()); } /* @@ -33,11 +35,15 @@ inline LPPolynomial lcm(const LPPolynomial& p, const LPPolynomial& q) { * @return content(p) */ inline LPPolynomial content(const LPPolynomial& p) { - return LPPolynomial(p.context(), poly::content(p.get_polynomial())); + lp_polynomial_t* res = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_cont(res, p.get_internal()); + return LPPolynomial(res, p.context()); } inline LPPolynomial derivative(const LPPolynomial& p) { - return LPPolynomial(p.context(), poly::derivative(p.get_polynomial())); + lp_polynomial_t* res = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_derivative(res, p.get_internal()); + return LPPolynomial(res, p.context()); } /* @@ -46,7 +52,9 @@ inline LPPolynomial derivative(const LPPolynomial& p) { * @return primitive_part(p) */ inline LPPolynomial primitive_part(const LPPolynomial& p) { - return LPPolynomial(p.context(), poly::primitive_part(p.get_polynomial())); + lp_polynomial_t* res = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_pp(res, p.get_internal()); + return LPPolynomial(res, p.context()); } /* @@ -55,7 +63,9 @@ inline LPPolynomial primitive_part(const LPPolynomial& p) { * @return resultant(p,q) */ inline LPPolynomial resultant(const LPPolynomial& p, const LPPolynomial& q) { - return LPPolynomial(p.context(), poly::resultant(p.get_polynomial(), q.get_polynomial())); + lp_polynomial_t* res = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_resultant(res, p.get_internal(), q.get_internal()); + return LPPolynomial(res, p.context()); } /* @@ -64,23 +74,19 @@ inline LPPolynomial resultant(const LPPolynomial& p, const LPPolynomial& q) { * @return discriminant(p) */ inline LPPolynomial discriminant(const LPPolynomial& p) { - assert(p.context().lp_context() == lp_polynomial_get_context(p.get_polynomial().get_internal())); - if (poly::degree(p.get_internal()) == 1) { // workaround for bug in the libpoly c++ wrapper + assert(p.context().lp_context() == lp_polynomial_get_context(p.get_internal())); + if (lp_polynomial_degree(p.get_internal()) == 1) { // workaround for bug in the libpoly c++ wrapper return LPPolynomial(p.context(), 1); } - return LPPolynomial(p.context(), poly::discriminant(p.get_polynomial())); -} - -/* - * @brief wrapper for the factorization of LPPolynomials, this is to be preferred when factorizing multiple polynomials with the same context - * @param LPPolynomial p, CoCoaAdaptorLP adaptor - * @return factorization(p) - */ -inline Factors factorization(const LPPolynomial& p, const CoCoAAdaptorLP& adaptor, bool includeConstant = true) { - #ifndef USE_COCOA - CARL_LOG_FATAL("carl.poly", "factorization is not supported without libcocoa"); - #endif - return adaptor.factorize(p, includeConstant); + // div(resultant(p, derivative(p)), leading_coefficient(p)); + lp_polynomial_t* ldcf = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_get_coefficient(ldcf, p.get_internal(), lp_polynomial_degree(p.get_internal())); + lp_polynomial_t* res = lp_polynomial_new(p.context().lp_context()); + lp_polynomial_derivative(res, p.get_internal()); + lp_polynomial_resultant(res, p.get_internal(), res); + lp_polynomial_div(res, res, ldcf); + lp_polynomial_delete(ldcf); + return LPPolynomial(res, p.context()); } /* @@ -88,67 +94,25 @@ inline Factors factorization(const LPPolynomial& p, const CoCoAAda * @param LPPolynomial p * @return factorization(p) */ -inline Factors factorization(const LPPolynomial& p, bool includeConstant = true) { - CoCoAAdaptorLP adaptor = CoCoAAdaptorLP(p.context()); - return factorization(p, adaptor, includeConstant); -} +Factors factorization(const LPPolynomial& p, bool includeConstant = true); -inline std::vector irreducible_factors(const LPPolynomial& p, const CoCoAAdaptorLP& adaptor, bool includeConstants = true) { - #ifndef USE_COCOA - CARL_LOG_FATAL("carl.poly", "factorization is not supported without libcocoa"); - #endif - - if (is_constant(p)) { - if (includeConstants) { - return {p}; - } else { - return {}; - } - } else if (p.total_degree() == 1) { - return {p}; - } - - return adaptor.irreducible_factors(p, includeConstants); -} - -inline std::vector irreducible_factors(const LPPolynomial& p, bool includeConstants = true) { - CoCoAAdaptorLP adaptor = CoCoAAdaptorLP(p.context()); - return irreducible_factors(p, adaptor, includeConstants); -} +std::vector irreducible_factors(const LPPolynomial& p, bool includeConstants = true); /* * @brief wrapper for the square_free_factorization of LPPolynomials * @param LPPolynomial p * @return square_free_factorization(p) */ -inline std::vector square_free_factors(const LPPolynomial& p) { - std::vector factors = poly::square_free_factors(p.get_polynomial()); - std::vector result; - for (const auto& factor : factors) { - result.push_back(LPPolynomial(p.context(), std::move(factor))); - } - return result; -} +std::vector square_free_factors(const LPPolynomial& p); /* * @brief wrapper for the content_free_factors of LPPolynomials * @param LPPolynomial p * @return content_free_factors(p) */ -inline std::vector content_free_factors(const LPPolynomial& p) { - std::vector factors = poly::content_free_factors(p.get_polynomial()); - std::vector result; - for (const auto& factor : factors) { - result.push_back(LPPolynomial(p.context(), std::move(factor))); - } - return result; -} +std::vector content_free_factors(const LPPolynomial& p); -inline std::vector groebner_basis(const std::vector& polys) { - if (polys.size() <= 1) return polys; - CoCoAAdaptorLP adaptor = CoCoAAdaptorLP(polys.at(0).context()); - return adaptor.GBasis(polys); -} +std::vector groebner_basis(const std::vector& polys); } // namespace carl diff --git a/src/carl-arith/poly/libpoly/LPContext.cpp b/src/carl-arith/poly/libpoly/LPContext.cpp new file mode 100644 index 00000000..1632647a --- /dev/null +++ b/src/carl-arith/poly/libpoly/LPContext.cpp @@ -0,0 +1,33 @@ +#include +#ifdef USE_LIBPOLY + +#include "LPContext.h" + +namespace carl { + +LPContext::Data::Data(const std::vector& v) : variable_order(v) { + lp_var_order = lp_variable_order_new(); + // lp_context = lp_polynomial_context_new(0, LPVariables::getInstance().lp_var_db, lp_var_order); + lp_context = (lp_polynomial_context_t*) malloc(sizeof(lp_polynomial_context_t)); + //lp_polynomial_context_construct(lp_context, 0, LPVariables::getInstance().lp_var_db, lp_var_order); + lp_context->ref_count = 0; + lp_context->var_db = LPVariables::getInstance().lp_var_db; + lp_context->K = 0; + lp_context->var_order = lp_var_order; + #define TEMP_VARIABLE_SIZE 10 + lp_context->var_tmp = (lp_variable_t*)malloc(sizeof(lp_variable_t)*TEMP_VARIABLE_SIZE); + for (size_t i = 0; i < TEMP_VARIABLE_SIZE; ++ i) { + lp_context->var_tmp[i] = LPVariables::getInstance().lp_var_tmp[i]; + } + lp_context->var_tmp_size = 0; + lp_polynomial_context_attach(lp_context); +} + +LPContext::Data::~Data() { + lp_variable_order_detach(lp_var_order); + lp_polynomial_context_detach(lp_context); +} + +} + +#endif \ No newline at end of file diff --git a/src/carl-arith/poly/libpoly/LPContext.h b/src/carl-arith/poly/libpoly/LPContext.h index 960ee733..7ca640dd 100644 --- a/src/carl-arith/poly/libpoly/LPContext.h +++ b/src/carl-arith/poly/libpoly/LPContext.h @@ -6,35 +6,27 @@ #include "../../core/Variable.h" #include +#include +#include #include #include -#include #include +#include "LPVariables.h" + namespace carl { class LPContext { struct Data { std::vector variable_order; - // mapping from carl variables to libpoly variables - std::map vars_carl_libpoly; - // mapping from libpoly variables to carl variables - std::map vars_libpoly_carl; - - // lp_variable_db_t* lp_var_db; - // lp_variable_order_t* lp_var_order; - // lp_polynomial_context_t* lp_context; - - poly::Context poly_context; - - Data(const std::vector& v) : variable_order(v) {} - // ~Data() { - // lp_variable_db_detach(lp_var_db); - // lp_variable_order_detach(lp_var_order); - // lp_polynomial_context_detach(lp_context); - // } + + lp_variable_order_t* lp_var_order; + lp_polynomial_context_t* lp_context; + + Data(const std::vector& v); + ~Data(); }; std::shared_ptr m_data; @@ -55,17 +47,15 @@ class LPContext { } std::optional carl_variable(lp_variable_t var) const { - auto it = m_data->vars_libpoly_carl.find(var); - if(it == m_data->vars_libpoly_carl.end()) return std::nullopt; - CARL_LOG_TRACE("carl.poly", "Mapping libpoly variable " << lp_variable_db_get_name(m_data->poly_context.get_variable_db(), var) << " -> " << it->second); - return it->second; + return LPVariables::getInstance().carl_variable(var); } - std::optional lp_variable(carl::Variable var) const { - auto it = m_data->vars_carl_libpoly.find(var); - if(it == m_data->vars_carl_libpoly.end()) return std::nullopt; - CARL_LOG_TRACE("carl.poly", "Mapping carl variable " << var << " -> " << lp_variable_db_get_name(m_data->poly_context.get_variable_db(), it->second)); - return it->second; + lp_variable_t lp_variable(carl::Variable var) const { + return LPVariables::getInstance().lp_variable(var); + } + + std::optional lp_variable_opt(carl::Variable var) const { + return LPVariables::getInstance().lp_variable_opt(var); } @@ -74,33 +64,18 @@ class LPContext { */ LPContext(const std::vector& carl_var_order) : m_data(std::make_shared(carl_var_order)) { CARL_LOG_FUNC("carl.poly", carl_var_order); - // m_data->lp_var_db = lp_variable_db_new(); - // m_data->lp_var_order = lp_variable_order_new(); - // m_data->poly_context = lp_polynomial_context_new(0, m_data->lp_var_db, m_data->lp_var_order); for (size_t i = 0; i < carl_var_order.size(); i++) { - std::string var_name = carl_var_order[i].name(); - lp_variable_t poly_var = lp_variable_db_new_variable(m_data->poly_context.get_variable_db(), &var_name[0]); - lp_variable_order_push(m_data->poly_context.get_variable_order(), poly_var); - CARL_LOG_TRACE("carl.poly", "Creating lp var for carl var " << carl_var_order[i] << " -> " << lp_variable_db_get_name(m_data->poly_context.get_variable_db(), poly_var)); - m_data->vars_carl_libpoly.emplace(carl_var_order[i], poly_var); - m_data->vars_libpoly_carl.emplace(poly_var, carl_var_order[i]); + lp_variable_t poly_var = LPVariables::getInstance().lp_variable(carl_var_order[i]); + lp_variable_order_push(m_data->lp_var_order, poly_var); } }; - poly::Context& poly_context() { - return m_data->poly_context; - } - - const poly::Context& poly_context() const { - return m_data->poly_context; - } - lp_polynomial_context_t* lp_context() { - return m_data->poly_context.get_polynomial_context(); + return m_data->lp_context; }; const lp_polynomial_context_t* lp_context() const { - return m_data->poly_context.get_polynomial_context(); + return m_data->lp_context; }; const std::vector& variable_ordering() const { @@ -117,6 +92,16 @@ class LPContext { inline bool operator==(const LPContext& rhs) const { return m_data == rhs.m_data; } + + bool is_extension_of(const LPContext& other) const { + auto it_a = variable_ordering().begin(); + auto it_b = other.variable_ordering().begin(); + while (it_a != variable_ordering().end() && it_b != other.variable_ordering().end() && *it_a == *it_b) { + it_a++; + it_b++; + } + return it_b == other.variable_ordering().end(); + } }; inline std::ostream& operator<<(std::ostream& os, const LPContext& ctx) { diff --git a/src/carl-arith/poly/libpoly/LPPolynomial.cpp b/src/carl-arith/poly/libpoly/LPPolynomial.cpp index 55c92335..090d96b0 100644 --- a/src/carl-arith/poly/libpoly/LPPolynomial.cpp +++ b/src/carl-arith/poly/libpoly/LPPolynomial.cpp @@ -1,5 +1,6 @@ #include "LPPolynomial.h" #include "helper.h" +#include #include #ifdef USE_LIBPOLY @@ -7,154 +8,155 @@ namespace carl { LPPolynomial::LPPolynomial(const LPPolynomial& rhs) - : m_poly(rhs.m_poly), m_context(rhs.m_context) { - assert(lp_polynomial_check_order(m_poly.get_internal())); + : m_internal(lp_polynomial_new_copy(rhs.m_internal)), m_context(rhs.m_context) { + assert(lp_polynomial_check_order(get_internal())); } LPPolynomial::LPPolynomial(LPPolynomial&& rhs) - : m_poly(std::move(rhs.m_poly)), m_context(std::move(rhs.m_context)) { - assert(lp_polynomial_check_order(m_poly.get_internal())); + : m_internal(rhs.m_internal), m_context(std::move(rhs.m_context)) { + assert(lp_polynomial_check_order(get_internal())); } LPPolynomial& LPPolynomial::operator=(const LPPolynomial& rhs) { - m_poly = rhs.m_poly; + m_internal = lp_polynomial_new_copy(rhs.m_internal); m_context = rhs.m_context; - assert(lp_polynomial_check_order(m_poly.get_internal())); + assert(lp_polynomial_check_order(get_internal())); return *this; } LPPolynomial& LPPolynomial::operator=(LPPolynomial&& rhs) { - m_poly = std::move(rhs.m_poly); + m_internal = rhs.m_internal; m_context = std::move(rhs.m_context); - assert(lp_polynomial_check_order(m_poly.get_internal())); + assert(lp_polynomial_check_order(get_internal())); return *this; } LPPolynomial::LPPolynomial(const LPContext& context) - : m_poly(context.lp_context()), m_context(context) { - //lp_polynomial_set_external(m_poly.get_internal()); - assert(lp_polynomial_check_order(m_poly.get_internal())); + : m_internal(lp_polynomial_new(context.lp_context())), m_context(context) { + //lp_polynomial_set_external(get_internal()); + assert(lp_polynomial_check_order(get_internal())); } -LPPolynomial::LPPolynomial(const LPContext& context, const poly::Polynomial& p) - : m_poly(p), m_context(context) { - //lp_polynomial_set_external(m_poly.get_internal()); - - assert(lp_polynomial_check_order(m_poly.get_internal())); - assert(context.lp_context() == lp_polynomial_get_context(get_internal())); -} -LPPolynomial::LPPolynomial(const LPContext& context, poly::Polynomial&& p) - : m_poly(std::move(p)), m_context(context) { - //lp_polynomial_set_external(m_poly.get_internal()); - - assert(lp_polynomial_check_order(m_poly.get_internal())); +LPPolynomial::LPPolynomial(lp_polynomial_t* p, const LPContext& context) + : m_internal(p), m_context(context) { + //lp_polynomial_set_external(get_internal()); + assert(lp_polynomial_check_order(get_internal())); assert(context.lp_context() == lp_polynomial_get_context(get_internal())); } LPPolynomial::LPPolynomial(const LPContext& context, long val) - : m_context(context) { - lp_polynomial_construct_simple(m_poly.get_internal(), context.lp_context(), poly::Integer(val).get_internal(), 0, 0); - //lp_polynomial_set_external(m_poly.get_internal()); - - assert(lp_polynomial_check_order(m_poly.get_internal())); + : m_internal(lp_polynomial_alloc()), m_context(context) { + lp_polynomial_construct_simple(get_internal(), context.lp_context(), mpz_class(val).get_mpz_t(), 0, 0); + //lp_polynomial_set_external(get_internal()); + assert(lp_polynomial_check_order(get_internal())); } LPPolynomial::LPPolynomial(const LPContext& context, const mpz_class& val) - : m_context(context) { - lp_polynomial_construct_simple(m_poly.get_internal(), context.lp_context(), val.get_mpz_t(), lp_variable_null, 0) ; - //lp_polynomial_set_external(m_poly.get_internal()); - - assert(lp_polynomial_check_order(m_poly.get_internal())); + : m_internal(lp_polynomial_alloc()), m_context(context) { + lp_polynomial_construct_simple(get_internal(), context.lp_context(), val.get_mpz_t(), lp_variable_null, 0) ; + //lp_polynomial_set_external(get_internal()); + assert(lp_polynomial_check_order(get_internal())); } LPPolynomial::LPPolynomial(const LPContext& context, const mpq_class& val) : LPPolynomial(context, carl::get_num(val)) {} - - LPPolynomial::LPPolynomial(const LPContext& context, const Variable& var, const mpz_class& coeff, unsigned int degree) - : m_context(context) { - lp_polynomial_construct_simple(m_poly.get_internal(), context.lp_context(), poly::Integer(coeff).get_internal(), *context.lp_variable(var), degree); - //lp_polynomial_set_external(m_poly.get_internal()); - - assert(lp_polynomial_check_order(m_poly.get_internal())); + : m_internal(lp_polynomial_alloc()), m_context(context) { + lp_polynomial_construct_simple(get_internal(), context.lp_context(), mpz_class(coeff).get_mpz_t(), context.lp_variable(var), degree); + //lp_polynomial_set_external(get_internal()); + assert(lp_polynomial_check_order(get_internal())); } LPPolynomial::LPPolynomial(const LPContext& context, const Variable& var) - : m_context(context) { - lp_polynomial_construct_simple(m_poly.get_internal(), context.lp_context(), poly::Integer(1).get_internal(), *context.lp_variable(var), 1); - //lp_polynomial_set_external(m_poly.get_internal()); - - assert(lp_polynomial_check_order(m_poly.get_internal())); + : m_internal(lp_polynomial_alloc()), m_context(context) { + lp_polynomial_construct_simple(get_internal(), context.lp_context(), mpz_class(1).get_mpz_t(), context.lp_variable(var), 1); + //lp_polynomial_set_external(get_internal()); + assert(lp_polynomial_check_order(get_internal())); } LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, const std::initializer_list& coefficients) - : m_poly(context.lp_context()), m_context(context) { + : LPPolynomial(context) { auto var = context.lp_variable(mainVar); - assert(var.has_value()); auto pow = coefficients.size(); for (const mpz_class& coeff : coefficients) { pow--; - poly::Polynomial temp; - lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(coeff).get_internal(), *var, (unsigned int)pow); - m_poly += temp; + if (is_zero(coeff)) continue; + lp_monomial_t t; + lp_monomial_construct(context.lp_context(), &t); + lp_monomial_set_coefficient(context.lp_context(), &t, mpz_class(coeff).get_mpz_t()); + if (pow>0) + lp_monomial_push(&t, var, (unsigned int)pow); + lp_polynomial_add_monomial(get_internal(), &t); + lp_monomial_destruct(&t); } - //lp_polynomial_set_external(m_poly.get_internal()); + //lp_polynomial_set_external(get_internal()); } LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, const std::vector& coefficients) - : m_poly(context.lp_context()), m_context(context) { + : LPPolynomial(context) { auto var = context.lp_variable(mainVar); - assert(var.has_value()); auto pow = coefficients.size(); for (const mpz_class& coeff : coefficients) { pow--; if (is_zero(coeff)) continue; - poly::Polynomial temp; - lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(coeff).get_internal(), *var, (unsigned int)pow); - m_poly += temp; + lp_monomial_t t; + lp_monomial_construct(context.lp_context(), &t); + lp_monomial_set_coefficient(context.lp_context(), &t, mpz_class(coeff).get_mpz_t()); + if (pow>0) + lp_monomial_push(&t, var, (unsigned int)pow); + lp_polynomial_add_monomial(get_internal(), &t); + lp_monomial_destruct(&t); } - //lp_polynomial_set_external(m_poly.get_internal()); + //lp_polynomial_set_external(get_internal()); } LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, std::vector&& coefficients) - : m_poly(context.lp_context()), m_context(context) { + : LPPolynomial(context) { auto var = context.lp_variable(mainVar); - assert(var.has_value()); auto pow = coefficients.size(); for (const mpz_class& coeff : coefficients) { pow--; - poly::Polynomial temp; - lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(std::move(coeff)).get_internal(), *var, (unsigned int)pow); - m_poly += temp; + if (is_zero(coeff)) continue; + lp_monomial_t t; + lp_monomial_construct(context.lp_context(), &t); + lp_monomial_set_coefficient(context.lp_context(), &t, mpz_class(coeff).get_mpz_t()); + if (pow>0) + lp_monomial_push(&t, var, (unsigned int)pow); + lp_polynomial_add_monomial(get_internal(), &t); + lp_monomial_destruct(&t); } - //lp_polynomial_set_external(m_poly.get_internal()); + //lp_polynomial_set_external(get_internal()); } LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, const std::map& coefficients) - : m_poly(context.lp_context()), m_context(context) { + : LPPolynomial(context) { auto var = context.lp_variable(mainVar); - assert(var.has_value()); - for (const auto& coef : coefficients) { - poly::Polynomial temp; - lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(coef.second).get_internal(), *var, coef.first); - m_poly += temp; + for (const auto& coeff : coefficients) { + if (is_zero(coeff.second)) continue; + lp_monomial_t t; + lp_monomial_construct(context.lp_context(), &t); + lp_monomial_set_coefficient(context.lp_context(), &t, mpz_class(coeff.second).get_mpz_t()); + if (coeff.first>0) + lp_monomial_push(&t, var, (unsigned int)coeff.first); + lp_polynomial_add_monomial(get_internal(), &t); + lp_monomial_destruct(&t); } - //lp_polynomial_set_external(m_poly.get_internal()); + //lp_polynomial_set_external(get_internal()); } bool LPPolynomial::has(const Variable& var) const { lp_variable_list_t varList; lp_variable_list_construct(&varList); - lp_polynomial_get_variables(m_poly.get_internal(), &varList); - auto lp_variable = context().lp_variable(var); + lp_polynomial_get_variables(get_internal(), &varList); + auto lp_variable = context().lp_variable_opt(var); if (!lp_variable) return false; bool contains = lp_variable_list_contains(&varList, *lp_variable); lp_variable_list_destruct(&varList); @@ -184,52 +186,51 @@ bool operator!=(const mpz_class& lhs, const LPPolynomial& rhs) { return !(lhs == rhs); } +inline auto cmp_util(const LPPolynomial& lhs, const mpz_class& rhs) { + lp_polynomial_t* tmp = poly_helper::construct_lp_poly(lp_polynomial_get_context(lhs.get_internal()), rhs); + auto res = lp_polynomial_cmp(lhs.get_internal(), tmp); + lp_polynomial_delete(tmp); + return res; +} + bool operator<(const LPPolynomial& lhs, const LPPolynomial& rhs) { return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) < 0; } bool operator<(const LPPolynomial& lhs, const mpz_class& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(lhs.get_internal()), poly::Integer(rhs)); - return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) < 0; + return cmp_util(lhs,rhs) < 0; } bool operator<(const mpz_class& lhs, const LPPolynomial& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(rhs.get_internal()), poly::Integer(lhs)); - return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) < 0; + return cmp_util(rhs,lhs) > 0; } bool operator<=(const LPPolynomial& lhs, const LPPolynomial& rhs) { return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) <= 0; } bool operator<=(const LPPolynomial& lhs, const mpz_class& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(lhs.get_internal()), poly::Integer(rhs)); - return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) <= 0; + return cmp_util(lhs,rhs) <= 0; } bool operator<=(const mpz_class& lhs, const LPPolynomial& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(rhs.get_internal()), poly::Integer(lhs)); - return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) <= 0; + return cmp_util(rhs,lhs) >= 0; } bool operator>(const LPPolynomial& lhs, const LPPolynomial& rhs) { return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) > 0; } bool operator>(const LPPolynomial& lhs, const mpz_class& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(lhs.get_internal()), poly::Integer(rhs)); - return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) > 0; + return cmp_util(lhs,rhs) > 0; } bool operator>(const mpz_class& lhs, const LPPolynomial& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(rhs.get_internal()), poly::Integer(lhs)); - return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) > 0; + return cmp_util(rhs,lhs) < 0; } bool operator>=(const LPPolynomial& lhs, const LPPolynomial& rhs) { return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) >= 0; } bool operator>=(const LPPolynomial& lhs, const mpz_class& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(lhs.get_internal()), poly::Integer(rhs)); - return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) >= 0; + return cmp_util(lhs,rhs) >= 0; } bool operator>=(const mpz_class& lhs, const LPPolynomial& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(rhs.get_internal()), poly::Integer(lhs)); - return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) >= 0; + return cmp_util(rhs,lhs) <= 0; } LPPolynomial operator+(const LPPolynomial& lhs, const LPPolynomial& rhs) { @@ -281,8 +282,9 @@ LPPolynomial& operator+=(LPPolynomial& lhs, const LPPolynomial& rhs) { return lhs; } LPPolynomial& operator+=(LPPolynomial& lhs, const mpz_class& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(lhs.get_internal()), poly::Integer(rhs)); - lp_polynomial_add(lhs.get_internal(), lhs.get_internal(), tmp.get_internal()); + lp_polynomial_t* tmp = poly_helper::construct_lp_poly(lp_polynomial_get_context(lhs.get_internal()), rhs); + lp_polynomial_add(lhs.get_internal(), lhs.get_internal(), tmp); + lp_polynomial_delete(tmp); return lhs; } @@ -293,8 +295,9 @@ LPPolynomial& operator-=(LPPolynomial& lhs, const LPPolynomial& rhs) { return lhs; } LPPolynomial& operator-=(LPPolynomial& lhs, const mpz_class& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(lhs.get_internal()), poly::Integer(rhs)); - lp_polynomial_sub(lhs.get_internal(), lhs.get_internal(), tmp.get_internal()); + lp_polynomial_t* tmp = poly_helper::construct_lp_poly(lp_polynomial_get_context(lhs.get_internal()), rhs); + lp_polynomial_sub(lhs.get_internal(), lhs.get_internal(), tmp); + lp_polynomial_delete(tmp); return lhs; } @@ -305,8 +308,9 @@ LPPolynomial& operator*=(LPPolynomial& lhs, const LPPolynomial& rhs) { return lhs; } LPPolynomial& operator*=(LPPolynomial& lhs, const mpz_class& rhs) { - poly::Polynomial tmp = poly_helper::construct_poly(lp_polynomial_get_context(lhs.get_internal()), poly::Integer(rhs)); - lp_polynomial_mul(lhs.get_internal(), lhs.get_internal(), tmp.get_internal()); + lp_polynomial_t* tmp = poly_helper::construct_lp_poly(lp_polynomial_get_context(lhs.get_internal()), rhs); + lp_polynomial_mul(lhs.get_internal(), lhs.get_internal(), tmp); + lp_polynomial_delete(tmp); return lhs; } @@ -341,8 +345,11 @@ mpz_class LPPolynomial::coprime_factor() const { LPPolynomial LPPolynomial::coprime_coefficients() const { mpz_class g = coprime_factor(); if (g == 1) return *this; - poly::Polynomial temp = poly_helper::construct_poly(lp_polynomial_get_context(get_internal()), poly::Integer(g)); - return LPPolynomial(context(), poly::div(get_polynomial(), temp)); + lp_polynomial_t* temp = poly_helper::construct_lp_poly(lp_polynomial_get_context(get_internal()), g); + lp_polynomial_t* res = lp_polynomial_new(context().lp_context()); + lp_polynomial_div(res, get_internal(), temp); + lp_polynomial_delete(temp); + return LPPolynomial(res, context()); } std::size_t LPPolynomial::total_degree() const { @@ -393,9 +400,7 @@ std::size_t LPPolynomial::degree(Variable::Arg var) const { }; degree_travers travers; - auto lp_var = context().lp_variable(var); - assert(lp_var.has_value()); - travers.var = *lp_var; + travers.var = context().lp_variable(var); lp_polynomial_traverse(get_internal(), getDegree, &travers); return travers.degree; @@ -448,9 +453,7 @@ std::vector LPPolynomial::monomial_degrees(Variable::Arg var) cons }; degree_travers travers; - auto lp_var = context().lp_variable(var); - assert(lp_var.has_value()); - travers.var = *lp_var; + travers.var = context().lp_variable(var); lp_polynomial_traverse(get_internal(), getDegree, &travers); return travers.degree; @@ -463,8 +466,7 @@ mpz_class LPPolynomial::unit_part() const { if(is_zero(*this)) { return 1 ; } - return poly::lc_sgn(this->get_polynomial()) ; - + return lp_polynomial_lc_sgn(get_internal()) ; } LPPolynomial LPPolynomial::normalized() const { @@ -529,9 +531,7 @@ LPPolynomial LPPolynomial::coeff(Variable::Arg var, std::size_t exp) const { }; coeff_travers travers; - auto lp_var = context().lp_variable(var); - assert(lp_var.has_value()); - travers.var = *lp_var; + travers.var = context().lp_variable(var); travers.exp = exp; travers.ctx = lp_polynomial_get_context(get_internal()); lp_polynomial_traverse(get_internal(), getCoeff, &travers); @@ -555,10 +555,23 @@ bool LPPolynomial::is_normal() const { } std::ostream& operator<<(std::ostream& os, const LPPolynomial& p) { - os << lp_polynomial_to_string(p.m_poly.get_internal()); + os << lp_polynomial_to_string(p.get_internal()); return os; } +void LPPolynomial::set_context(const LPContext& c) { + for (auto& v : variables(*this)) assert(c.has(v)); + if (context() == c) return; + + bool reorder = !(c.is_extension_of(context()) || context().is_extension_of(c)); + m_context = c; + lp_polynomial_set_context(get_internal(), m_context.lp_context()); + if (reorder) { + lp_polynomial_ensure_order(get_internal()); + } + assert(lp_polynomial_check_order(get_internal())); +} + } // namespace carl #endif \ No newline at end of file diff --git a/src/carl-arith/poly/libpoly/LPPolynomial.h b/src/carl-arith/poly/libpoly/LPPolynomial.h index 28c557b0..6328e9bf 100644 --- a/src/carl-arith/poly/libpoly/LPPolynomial.h +++ b/src/carl-arith/poly/libpoly/LPPolynomial.h @@ -8,6 +8,7 @@ #ifdef USE_LIBPOLY #include "LPContext.h" +#include #include "helper.h" #include @@ -17,7 +18,6 @@ #include #include #include -#include #include #include @@ -26,7 +26,7 @@ namespace carl { class LPPolynomial { private: /// The libpoly polynomial. - poly::Polynomial m_poly; + lp_polynomial_t* m_internal; LPContext m_context; @@ -68,18 +68,9 @@ class LPPolynomial { explicit LPPolynomial(const LPContext& context); /** - * Construct a LPPolynomial with the given libpoly polynomial. - * Also uses the given context. - * @param mainPoly Libpoly Polynomial. - */ - LPPolynomial(const LPContext& context, const poly::Polynomial& mainPoly); - - /** - * Moves a LPPolynomial with the given libpoly polynomial. - * Also uses the given context. - * @param mainPoly Libpoly Polynomial. + * Move constructor. */ - LPPolynomial(const LPContext& context, poly::Polynomial&& mainPoly); + LPPolynomial(lp_polynomial_t* p, const LPContext& context); /** * Construct a LPPolynomial with only a integer. @@ -166,31 +157,44 @@ class LPPolynomial { * @return The only variable occurring in the term. */ Variable single_variable() const { - assert(poly::is_univariate(m_poly)); - auto carl_var = context().carl_variable(poly::main_variable(m_poly).get_internal()); + assert(lp_polynomial_is_univariate(get_internal())); + auto carl_var = context().carl_variable(lp_polynomial_top_variable(get_internal())); assert(carl_var.has_value()); return *carl_var; } + LPPolynomial coeff(std::size_t k) const { + lp_polynomial_t* res = lp_polynomial_alloc(); + lp_polynomial_construct(res, m_context.lp_context()); + lp_polynomial_get_coefficient(res, get_internal(), k); + return LPPolynomial(res, m_context); + } + + /** + * Get the maximal exponent of the main variable. + * As the degree of the zero polynomial is \f$-\infty\f$, we assert that this polynomial is not zero. This must be checked by the caller before calling this method. + * @see @cite GCL92, page 38 + * @return Degree. + */ + size_t degree() const { + return lp_polynomial_degree(get_internal()); + } + /** * Returns the leading coefficient. * @return The leading coefficient. */ LPPolynomial lcoeff() const { - return LPPolynomial(m_context, poly::leading_coefficient(m_poly)); - } - - LPPolynomial coeff(std::size_t k) const { - return LPPolynomial(m_context, poly::coefficient(m_poly, k)); + return coeff(degree()); } /** Obtain all non-zero coefficients of a polynomial. */ std::vector coefficients() const { std::vector res; - for (std::size_t deg = 0; deg <= poly::degree(m_poly); ++deg) { - auto coeff = poly::coefficient(m_poly, deg); - if (lp_polynomial_is_zero(coeff.get_internal())) continue; - res.emplace_back(context(), std::move(coeff)); + for (std::size_t deg = 0; deg <= degree(); ++deg) { + auto cf = coeff(deg); + if (lp_polynomial_is_zero(cf.get_internal())) continue; + res.emplace_back(std::move(cf)); } return res; } @@ -214,25 +218,20 @@ class LPPolynomial { }; LPPolynomial_constantPart_visitor visitor; - lp_polynomial_traverse(m_poly.get_internal(), getConstantPart, &visitor); + lp_polynomial_traverse(get_internal(), getConstantPart, &visitor); return visitor.part; } - /** - * Get the maximal exponent of the main variable. - * As the degree of the zero polynomial is \f$-\infty\f$, we assert that this polynomial is not zero. This must be checked by the caller before calling this method. - * @see @cite GCL92, page 38 - * @return Degree. - */ - size_t degree() const { - return poly::degree(m_poly); - } - /** * Removes the leading term from the polynomial. */ void truncate() { - m_poly -= poly::leading_coefficient(m_poly); + lp_polynomial_t* lcoeff = lp_polynomial_alloc(); + lp_polynomial_construct(lcoeff, m_context.lp_context()); + lp_polynomial_get_coefficient(lcoeff, get_internal(), lp_polynomial_degree(get_internal())); + lp_polynomial_sub(get_internal(), get_internal(), lcoeff); + lp_polynomial_delete(lcoeff); + //*this -= lcoeff(); } /** @@ -240,8 +239,8 @@ class LPPolynomial { * @return Main variable. */ Variable main_var() const { - if (poly::is_constant(m_poly)) return carl::Variable::NO_VARIABLE; - else return *(context().carl_variable(poly::main_variable(m_poly).get_internal())); + if (lp_polynomial_is_constant(get_internal())) return carl::Variable::NO_VARIABLE; + else return *(context().carl_variable(lp_polynomial_top_variable(get_internal()))); } /** @@ -250,7 +249,7 @@ class LPPolynomial { * @return Libpoly Polynomial. */ lp_polynomial_t* get_internal() { - return m_poly.get_internal(); + return m_internal; } /** @@ -258,16 +257,7 @@ class LPPolynomial { * @return Libpoly Polynomial. */ const lp_polynomial_t* get_internal() const { - return m_poly.get_internal(); - } - - /** - * @brief Get the underlying Polynomial object - * - * @return const poly::Polynomial& - */ - const poly::Polynomial& get_polynomial() const { - return m_poly; + return m_internal; } /** @@ -288,6 +278,8 @@ class LPPolynomial { return m_context; } + void set_context(const LPContext& c); + /** * Checks if the given variable occurs in the polynomial. * @param v Variable. @@ -447,20 +439,15 @@ LPPolynomial& operator*=(LPPolynomial& lhs, const mpz_class& rhs); * @return If polynomial is zero. */ inline bool is_zero(const LPPolynomial& p) { - return poly::is_zero(p.get_polynomial()); + return lp_polynomial_is_zero(p.get_internal()); } -// bool isNegative(LPPolynomial& p) { -// // return poly::is_zero(*p.mainPoly()); -// return true; -// } - /** * Checks if the polynomial is linear or not. * @return If polynomial is linear. */ inline bool is_constant(const LPPolynomial& p) { - return poly::is_constant(p.get_polynomial()); + return lp_polynomial_is_constant(p.get_internal()); } /** @@ -471,9 +458,13 @@ inline bool is_one(const LPPolynomial& p) { if (!is_constant(p)) { return false; } - poly::Polynomial temp(lp_polynomial_get_context(p.get_internal())); - temp += poly::Integer(mpz_class(1)); - return lp_polynomial_eq(p.get_internal(), temp.get_internal()); + + lp_polynomial_t* one = lp_polynomial_alloc(); + mpz_class one_int(1); + lp_polynomial_construct_simple(one, p.context().lp_context(), one_int.get_mpz_t(), lp_variable_null, 0); + bool res = lp_polynomial_eq(p.get_internal(), one); + lp_polynomial_delete(one); + return res; } /** @@ -488,14 +479,14 @@ inline bool is_number(const LPPolynomial& p) { * @brief Check if the given polynomial is linear. */ inline bool is_linear(const LPPolynomial& p) { - return poly::is_linear(p.get_polynomial()); + return lp_polynomial_is_linear(p.get_internal()); } /** * @brief Check if the given polynomial is univariate. */ inline bool is_univariate(const LPPolynomial& p) { - return poly::is_univariate(p.get_polynomial()); + return lp_polynomial_is_univariate(p.get_internal()); } inline std::size_t level_of(const LPPolynomial& p) { diff --git a/src/carl-arith/poly/libpoly/LPVariables.cpp b/src/carl-arith/poly/libpoly/LPVariables.cpp new file mode 100644 index 00000000..cfd8c625 --- /dev/null +++ b/src/carl-arith/poly/libpoly/LPVariables.cpp @@ -0,0 +1,55 @@ +#include "LPVariables.h" + + +#include +#ifdef USE_LIBPOLY + + +namespace carl { + +LPVariables::LPVariables() { + lp_var_db = lp_variable_db_new(); + for (size_t i = 0; i < TEMP_VARIABLE_SIZE; ++ i) { + char name[10]; + sprintf(name, "#%zu", i); + lp_var_tmp[i] = lp_variable_db_new_variable(lp_var_db, name); + } +} + +LPVariables::~LPVariables() { + lp_variable_db_detach(lp_var_db); +} + + +std::optional LPVariables::carl_variable(lp_variable_t var) const { + auto it = vars_libpoly_carl.find(var); + if(it == vars_libpoly_carl.end()) return std::nullopt; + CARL_LOG_TRACE("carl.poly", "Mapping libpoly variable " << lp_variable_db_get_name(lp_var_db, var) << " (" << var << ") -> " << it->second << " (" << it->second.id() << ")"); + return it->second; +} + +std::optional LPVariables::lp_variable_opt(carl::Variable var) const { + auto it = vars_carl_libpoly.find(var); + if(it == vars_carl_libpoly.end()) return std::nullopt; + CARL_LOG_TRACE("carl.poly", "Mapping carl variable " << var << " (" << var.id() << ") -> " << lp_variable_db_get_name(lp_var_db, it->second) << " (" << it->second << ")"); + return it->second; +} + +lp_variable_t LPVariables::lp_variable(carl::Variable var) { + auto it = vars_carl_libpoly.find(var); + if(it != vars_carl_libpoly.end()) { + CARL_LOG_TRACE("carl.poly", "Mapping carl variable " << var << " (" << var.id() << ") -> " << lp_variable_db_get_name(lp_var_db, it->second) << " (" << it->second << ")"); + return it->second; + } else { + std::string var_name = var.name(); + lp_variable_t poly_var = lp_variable_db_new_variable(lp_var_db, &var_name[0]); + vars_carl_libpoly.emplace(var, poly_var); + vars_libpoly_carl.emplace(poly_var, var); + CARL_LOG_TRACE("carl.poly", "Mapping carl variable " << var << " (" << var.id() << ") -> " << lp_variable_db_get_name(lp_var_db, poly_var) << " (" << poly_var << ")"); + return poly_var; + } +} + +} // namespace carl + +#endif \ No newline at end of file diff --git a/src/carl-arith/poly/libpoly/LPVariables.h b/src/carl-arith/poly/libpoly/LPVariables.h new file mode 100644 index 00000000..be6661f1 --- /dev/null +++ b/src/carl-arith/poly/libpoly/LPVariables.h @@ -0,0 +1,40 @@ +#pragma once + +#include +#ifdef USE_LIBPOLY + +#include +#include +#include +#include "../../core/Variable.h" +#include +#include +#include + +namespace carl { + +class LPVariables : public Singleton { + friend Singleton; + + // mapping from carl variables to libpoly variables + std::map vars_carl_libpoly; + // mapping from libpoly variables to carl variables + std::map vars_libpoly_carl; + +public: + lp_variable_db_t* lp_var_db; + + #define TEMP_VARIABLE_SIZE 10 + + lp_variable_t lp_var_tmp[TEMP_VARIABLE_SIZE]; + + LPVariables(); + ~LPVariables(); + std::optional carl_variable(lp_variable_t var) const; + std::optional lp_variable_opt(carl::Variable var) const; + lp_variable_t lp_variable(carl::Variable var); +}; + +} // namespace carl + +#endif \ No newline at end of file diff --git a/src/carl-arith/poly/libpoly/helper.h b/src/carl-arith/poly/libpoly/helper.h index bf7a27d8..f37954f3 100644 --- a/src/carl-arith/poly/libpoly/helper.h +++ b/src/carl-arith/poly/libpoly/helper.h @@ -1,21 +1,23 @@ #pragma once +#include + /** * Helpers due to the shortcomings of libpoly's C++ API. * */ namespace carl::poly_helper { -inline poly::Polynomial construct_poly(const lp_polynomial_context_t* c, lp_variable_t v) { - poly::Polynomial poly; - lp_polynomial_construct_simple(poly.get_internal(), c, poly::Integer(1).get_internal(), v, 1); - return poly; +inline lp_polynomial_t* construct_lp_poly(const lp_polynomial_context_t* c, lp_variable_t v) { + lp_polynomial_t* p = lp_polynomial_alloc(); + lp_polynomial_construct_simple(p, c, mpz_class(1).get_mpz_t(), v, 1); + return p; } -inline poly::Polynomial construct_poly(const lp_polynomial_context_t* c, poly::Integer i) { - poly::Polynomial poly; - lp_polynomial_construct_simple(poly.get_internal(), c, i.get_internal(), lp_variable_null, 0) ; - return poly; +inline lp_polynomial_t* construct_lp_poly(const lp_polynomial_context_t* c, const mpz_class& i) { + lp_polynomial_t* p = lp_polynomial_alloc(); + lp_polynomial_construct_simple(p, c, i.get_mpz_t(), lp_variable_null, 0) ; + return p; } } diff --git a/src/carl-arith/ran/interval/Evaluation.h b/src/carl-arith/ran/interval/Evaluation.h index 3af4b843..c8048756 100644 --- a/src/carl-arith/ran/interval/Evaluation.h +++ b/src/carl-arith/ran/interval/Evaluation.h @@ -83,7 +83,7 @@ std::optional> evaluate(MultivariatePolynomial } CARL_LOG_TRACE("carl.ran.interval", "Compute result polynomial"); - Variable v = fresh_real_variable(); + static Variable v = fresh_real_variable(); std::vector>> algebraic_information; for (const auto& [var, ran] : m) { if (var_to_interval.find(var) == var_to_interval.end()) continue; @@ -211,7 +211,7 @@ boost::tribool evaluate(const BasicConstraint>& c } // compute the result polynomial - Variable v = fresh_real_variable(); + static Variable v = fresh_real_variable(); std::vector>> algebraic_information; for (const auto& [var, ran] : m) { if (var_to_interval.find(var) == var_to_interval.end()) continue; diff --git a/src/carl-arith/ran/libpoly/Evaluation.cpp b/src/carl-arith/ran/libpoly/Evaluation.cpp index 6786deeb..aea5be17 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -2,91 +2,66 @@ #ifdef USE_LIBPOLY +#include "LPAssignment.h" + + namespace carl { std::optional evaluate( const LPPolynomial& polynomial, const std::map& evalMap) { + lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); + auto result = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); - //Turn into poly::Assignment - poly::Assignment assignment(polynomial.context().poly_context()); - for (const auto& entry : evalMap) { - lp_value_t val; - //Turn into value - lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); - //That copies the value into the assignment - assignment.set(*(polynomial.context().lp_variable(entry.first)), poly::Value(&val)); - lp_value_destruct(&val); - } - - poly::Value result = poly::Value(lp_polynomial_evaluate(polynomial.get_internal(), assignment.get_internal())); - - CARL_LOG_DEBUG("carl.ran.libpoly", " Result: " << result); - - if(poly::is_none(result)) { + if (result->type == LP_VALUE_NONE) { + CARL_LOG_DEBUG("carl.ran.libpoly", " Result: NONE"); return std::nullopt; } - return LPRealAlgebraicNumber::create_from_value(std::move(result.get_internal())); + auto ran = LPRealAlgebraicNumber(std::move(*result)); + CARL_LOG_DEBUG("carl.ran.libpoly", " Result: " << ran); + return ran; } -boost::tribool evaluate(const BasicConstraint& constraint, const std::map& evalMap) { - - CARL_LOG_DEBUG("carl.ran.libpoly", " Evaluation constraint " << constraint << " for assignment " << evalMap); - - //Easy checks of lhs of constraint is constant - if (is_constant(constraint.lhs())) { - auto num = constraint.lhs().constant_part(); - switch (constraint.relation()) { - case Relation::EQ: - return num == 0; - case Relation::NEQ: - return num != 0; - case Relation::LESS: - return num < 0; - case Relation::LEQ: - return num <= 0; - case Relation::GREATER: - return num > 0; - case Relation::GEQ: - return num >= 0; - default: - assert(false); - return false; - } - } - - //denominator can be omitted - const poly::Polynomial& poly_pol = constraint.lhs().get_polynomial(); - - //Turn into poly::Assignment - poly::Assignment assignment(constraint.lhs().context().poly_context()); - for (const auto& entry : evalMap) { - lp_value_t val; - //Turn into value - lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); - //That copies the value into the assignment - assignment.set(*(constraint.lhs().context().lp_variable(entry.first)), poly::Value(&val)); - lp_value_destruct(&val); - } - - switch (constraint.relation()) { +inline auto lp_sign(carl::Relation rel) { + switch (rel) { case Relation::EQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::EQ); + return LP_SGN_EQ_0; case Relation::NEQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::NE); + return LP_SGN_NE_0; case Relation::LESS: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::LT); + return LP_SGN_LT_0; case Relation::LEQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::LE); + return LP_SGN_LE_0; case Relation::GREATER: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::GT); + return LP_SGN_GT_0; case Relation::GEQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::GE); + return LP_SGN_GE_0; default: assert(false); - return false; + return LP_SGN_EQ_0; + } +} + +boost::tribool evaluate(const BasicConstraint& constraint, const std::map& evalMap) { + CARL_LOG_DEBUG("carl.ran.libpoly", " Evaluation constraint " << constraint << " for assignment " << evalMap); + + if (is_constant(constraint.lhs())) { + return carl::evaluate(constraint.lhs().constant_part(), constraint.relation()); + } + + auto poly_pol = constraint.lhs().get_internal(); + lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); + + for (const auto& v : carl::variables(constraint)) { + if (evalMap.find(v) == evalMap.end()) { + int result = lp_polynomial_constraint_evaluate_subs(poly_pol, lp_sign(constraint.relation()), &assignment); + if (result == -1) return boost::indeterminate; + else return result; + } } + + return lp_polynomial_constraint_evaluate(poly_pol, lp_sign(constraint.relation()), &assignment); } } diff --git a/src/carl-arith/ran/libpoly/LPAssignment.cpp b/src/carl-arith/ran/libpoly/LPAssignment.cpp new file mode 100644 index 00000000..51a850f9 --- /dev/null +++ b/src/carl-arith/ran/libpoly/LPAssignment.cpp @@ -0,0 +1,39 @@ +#include "LPAssignment.h" + + +#include +#ifdef USE_LIBPOLY + +#include "../../poly/libpoly/LPVariables.h" + + +namespace carl { + +LPAssignment::LPAssignment() { + lp_assignment_construct(&lp_assignment, LPVariables::getInstance().lp_var_db); +} + +LPAssignment::~LPAssignment() { + lp_assignment_destruct(&lp_assignment); +} + +lp_assignment_t& LPAssignment::get(const carl::Assignment& ass) { + if (last_assignment == ass) { + return lp_assignment; + } else { + last_assignment = ass; + if (lp_assignment.values) { + for (size_t i = 0; i < lp_assignment.size; ++ i) { + lp_assignment_set_value(&lp_assignment, i, 0); + } + } + for (const auto& entry : ass) { + lp_assignment_set_value(&lp_assignment, LPVariables::getInstance().lp_variable(entry.first), entry.second.get_internal()); + } + return lp_assignment; + } +} + +} // namespace carl + +#endif \ No newline at end of file diff --git a/src/carl-arith/ran/libpoly/LPAssignment.h b/src/carl-arith/ran/libpoly/LPAssignment.h new file mode 100644 index 00000000..4ea19ccb --- /dev/null +++ b/src/carl-arith/ran/libpoly/LPAssignment.h @@ -0,0 +1,33 @@ +#pragma once + +#include +#ifdef USE_LIBPOLY + +#include +#include +#include +#include "../../core/Variable.h" +#include +#include +#include +#include +#include "LPRan.h" + + +namespace carl { + +class LPAssignment : public Singleton { + friend Singleton; + + lp_assignment_t lp_assignment; + carl::Assignment last_assignment; + +public: + LPAssignment(); + ~LPAssignment(); + lp_assignment_t& get(const carl::Assignment& ass); +}; + +} // namespace carl + +#endif \ No newline at end of file diff --git a/src/carl-arith/ran/libpoly/LPRan.cpp b/src/carl-arith/ran/libpoly/LPRan.cpp index 70d20fa1..f769ebfc 100644 --- a/src/carl-arith/ran/libpoly/LPRan.cpp +++ b/src/carl-arith/ran/libpoly/LPRan.cpp @@ -7,79 +7,71 @@ namespace carl { using NumberType = LPRealAlgebraicNumber::NumberType; LPRealAlgebraicNumber::~LPRealAlgebraicNumber() { - lp_algebraic_number_destruct(get_internal()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber() { - lp_algebraic_number_construct_zero(get_internal()); +LPRealAlgebraicNumber::LPRealAlgebraicNumber() : m_content(std::make_shared()) { + lp_value_construct_zero(get_internal()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const poly::AlgebraicNumber& num) { - lp_algebraic_number_construct_copy(get_internal(), num.get_internal()); +LPRealAlgebraicNumber::LPRealAlgebraicNumber(const lp_value_t& num) : m_content(std::make_shared()) { + lp_value_construct_copy(get_internal(), &num); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const lp_algebraic_number_t& num) { - lp_algebraic_number_construct_copy(get_internal(), &num); +LPRealAlgebraicNumber::LPRealAlgebraicNumber(lp_value_t&& num) : m_content(std::make_shared()) { + *get_internal() = num; } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(poly::AlgebraicNumber&& num) - : LPRealAlgebraicNumber() { - lp_algebraic_number_swap(get_internal(), num.get_internal()); -} - -LPRealAlgebraicNumber::LPRealAlgebraicNumber(lp_algebraic_number_t&& num) - : LPRealAlgebraicNumber() { - lp_algebraic_number_swap(get_internal(), &num); -} - -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial& p, const Interval& i) { - CARL_LOG_DEBUG("carl.ran.libpoly", " Create safe from poly: " << p << " in interval: " << i); +LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial& p, const Interval& i) : LPRealAlgebraicNumber() { + CARL_LOG_DEBUG("carl.ran.libpoly", " Create from poly: " << p << " in interval: " << i); - poly::UPolynomial upoly = to_libpoly_upolynomial(p); - poly::Interval inter_poly = to_libpoly_interval(i); + lp_upolynomial_t* upoly = to_libpoly_upolynomial(p); + lp_interval_t* inter_poly = to_libpoly_interval(i); - CARL_LOG_DEBUG("carl.ran.libpoly", " Converted to poly: " << upoly << " in interval: " << inter_poly); + lp_algebraic_number_t* roots = new lp_algebraic_number_t[lp_upolynomial_degree(upoly)]; + std::size_t roots_size; + lp_upolynomial_roots_isolate(upoly, roots, &roots_size); - //actual calculations - std::vector roots = poly::isolate_real_roots(upoly); - std::vector res; - - for (const auto& val : roots) { - // filter out roots not in interval - if (poly::contains(inter_poly, poly::Value(val))) { - CARL_LOG_DEBUG("carl.ran.libpoly", " Found Root " << val); - res.emplace_back(val); + bool found = false; + for (std::size_t i = 0; i < roots_size; ++i) { + lp_value_t tmp; + lp_value_construct(&tmp, LP_VALUE_ALGEBRAIC, &roots[i]); + if (lp_interval_contains(inter_poly, &tmp)) { + *get_internal() = tmp; + found = true; + break; } + lp_value_destruct(&tmp); + } + assert(found); + + for (std::size_t i = 0; i < roots_size; ++i) { + lp_algebraic_number_destruct(&roots[i]); } + delete[] roots; - assert(res.size() == 1); - lp_algebraic_number_construct_zero(get_internal()); - lp_algebraic_number_swap(get_internal(), res.front().get_internal()); - //We dont need to free the items in res, as they're from the C++ interface + lp_upolynomial_delete(upoly); + + lp_interval_destruct(inter_poly); + delete inter_poly; } /** * Construct from NumberType (usually mpq_class) */ -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const NumberType& num) { - poly::Rational rat = to_libpoly_rational(num); - lp_algebraic_number_construct_from_rational(get_internal(), rat.get_internal()); +LPRealAlgebraicNumber::LPRealAlgebraicNumber(const NumberType& num) : LPRealAlgebraicNumber() { + lp_value_construct(get_internal(), LP_VALUE_RATIONAL, (lp_rational_t*)num.get_mpq_t()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const LPRealAlgebraicNumber& ran) { - lp_algebraic_number_construct_copy(get_internal(), ran.get_internal()); +LPRealAlgebraicNumber::LPRealAlgebraicNumber(const LPRealAlgebraicNumber& ran) : m_content(ran.m_content) { } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(LPRealAlgebraicNumber&& ran) - : LPRealAlgebraicNumber() { - lp_algebraic_number_swap(get_internal(), ran.get_internal()); +LPRealAlgebraicNumber::LPRealAlgebraicNumber(LPRealAlgebraicNumber&& ran) : m_content(std::move(ran.m_content)) { } LPRealAlgebraicNumber& LPRealAlgebraicNumber::operator=(const LPRealAlgebraicNumber& n) { - lp_algebraic_number_destruct(get_internal()); - lp_algebraic_number_construct_copy(get_internal(), n.get_internal()); + m_content = n.m_content; return *this; } LPRealAlgebraicNumber& LPRealAlgebraicNumber::operator=(LPRealAlgebraicNumber&& n) { - lp_algebraic_number_swap(get_internal(), n.get_internal()); + m_content = std::move(n.m_content); return *this; } @@ -90,67 +82,28 @@ LPRealAlgebraicNumber LPRealAlgebraicNumber::create_safe(const carl::UnivariateP return LPRealAlgebraicNumber(p, i); } -/** - * Convert a libpoly Value into an algebraic number - * This does not free the value - * In case the value is none or infinity this conversion is not possible - * @param Libpoly Value (C interface) - * @return Copy of val as a algebraic number - */ -LPRealAlgebraicNumber LPRealAlgebraicNumber::create_from_value(const lp_value_t* val) { - CARL_LOG_DEBUG("carl.ran.libpoly", "Converting value into algebraic number"); - lp_algebraic_number_t mVal; - if (val->type == lp_value_type_t::LP_VALUE_NONE || val->type == lp_value_type_t::LP_VALUE_MINUS_INFINITY || val->type == lp_value_type_t::LP_VALUE_PLUS_INFINITY) { - //If value is not assigned at all or pm infty just return an error - assert(false && "Invalid RAN creation"); - } else if (val->type == lp_value_type_t::LP_VALUE_ALGEBRAIC) { - //val is already an algebraic number - lp_algebraic_number_construct_copy(&mVal, &val->value.a); - LPRealAlgebraicNumber ret(std::move(mVal)); - lp_algebraic_number_destruct(&mVal); - return ret; - } else if (lp_value_is_rational(val)) { - //if the value is a rational number: either an integer, dyadic rational or a rational - lp_rational_t rat; - lp_rational_construct(&rat); - lp_value_get_rational(val, &rat); - lp_algebraic_number_construct_from_rational(&mVal, &rat); - LPRealAlgebraicNumber ret(std::move(mVal)); - lp_algebraic_number_destruct(&mVal); - lp_rational_destruct(&rat); - return ret; - } else { - assert(false && "Invalid RAN creation"); - } - return LPRealAlgebraicNumber(); -} - bool is_integer(const LPRealAlgebraicNumber& n) { - return lp_algebraic_number_is_integer(n.get_internal()); + return lp_value_is_integer(n.get_internal()); } LPRealAlgebraicNumber::NumberType integer_below(const LPRealAlgebraicNumber& n) { - poly::Integer val; - lp_algebraic_number_floor(n.get_internal(), val.get_internal()); - return to_rational(val); + lp_integer_t result; + lp_integer_construct(&result); + lp_value_floor(n.get_internal(), &result); + NumberType num = to_rational(&result) ; + lp_integer_destruct(&result) ; + return num; } /** * @return true if the interval is a point or the poly has degree 1 */ bool LPRealAlgebraicNumber::is_numeric() const { - return lp_algebraic_number_is_rational(get_internal()); -} - -const poly::UPolynomial LPRealAlgebraicNumber::libpoly_polynomial() const { - return poly::UPolynomial(static_cast(get_internal()->f)); -} - -const poly::DyadicInterval& LPRealAlgebraicNumber::libpoly_interval() const { - return *reinterpret_cast(&get_internal()->I); + return lp_value_is_rational(get_internal()); } const UnivariatePolynomial LPRealAlgebraicNumber::polynomial() const { - return to_carl_univariate_polynomial(libpoly_polynomial(), auxVariable); + assert(!is_numeric()); + return to_carl_univariate_polynomial(get_internal()->value.a.f, auxVariable); } const Interval LPRealAlgebraicNumber::interval() const { @@ -161,12 +114,14 @@ const Interval LPRealAlgebraicNumber::interval() const { } } -const NumberType LPRealAlgebraicNumber::get_upper_bound() const { - return to_rational(poly::get_upper(libpoly_interval())); +const NumberType LPRealAlgebraicNumber::get_lower_bound() const { + if (is_numeric()) return value(); + else return to_rational(&get_internal()->value.a.I.a); } -const NumberType LPRealAlgebraicNumber::get_lower_bound() const { - return to_rational(poly::get_lower(libpoly_interval())); +const NumberType LPRealAlgebraicNumber::get_upper_bound() const { + if (is_numeric()) return value(); + else return to_rational(&get_internal()->value.a.I.b); } /** @@ -175,32 +130,12 @@ const NumberType LPRealAlgebraicNumber::get_lower_bound() const { */ const NumberType LPRealAlgebraicNumber::value() const { assert(is_numeric()); - lp_rational_t result ; - if(lp_dyadic_interval_is_point(&get_internal()->I)){ - // It's a point value, so we just get it - lp_rational_construct_from_dyadic(&result, lp_dyadic_interval_get_point(&get_internal()->I)); - }else{ - const lp_upolynomial_t* poly = get_internal()->f ; - assert(lp_upolynomial_degree(poly) == 1) ; - // p = ax + b = 0 => x = -b/a - const lp_integer_t* b = lp_upolynomial_const_term(poly); - const lp_integer_t* a = lp_upolynomial_lead_coeff(poly); - if(b){ - //b != 0 - //we can do this, as lp_rational_t == __mpq_struct - mpq_init(&result); - mpq_set_num(&result, b); - mpq_set_den(&result, a); - mpq_canonicalize(&result); - mpq_neg(&result, &result); - }else{ - // b = 0 - mpq_init(&result); - } - } + lp_rational_t result; + lp_rational_construct(&result); + lp_value_get_rational(get_internal(), &result); NumberType num = to_rational(&result) ; lp_rational_destruct(&result) ; - return num ; + return num; } /** @@ -208,38 +143,35 @@ const NumberType LPRealAlgebraicNumber::value() const { * @return Absolute Value of the stored RAN */ LPRealAlgebraicNumber abs(const LPRealAlgebraicNumber& n) { - if (n.is_numeric()) { CARL_LOG_DEBUG("carl.ran.libpoly", "Algebraic NumberType abs got numeric value"); return LPRealAlgebraicNumber(carl::abs(n.value())); } - int sign = lp_algebraic_number_sgn(n.get_internal()); + int sign = lp_value_sgn(n.get_internal()); CARL_LOG_DEBUG("carl.ran.libpoly", "Algebraic NumberType abs got sign : " << sign << " of " << n); if (sign >= 0) { return LPRealAlgebraicNumber(n); } else { - lp_algebraic_number_t val; - lp_algebraic_number_construct_zero(&val); - lp_algebraic_number_neg(&val, n.get_internal()); + lp_value_t val; + lp_value_construct_none(&val); + lp_value_neg(&val, n.get_internal()); auto ret = LPRealAlgebraicNumber(std::move(val)); - lp_algebraic_number_destruct(&val); - return ret ; + return ret; } } std::size_t bitsize(const LPRealAlgebraicNumber& n) { - //From Ran.h if (n.is_numeric()) { return carl::bitsize(n.get_lower_bound()) + carl::bitsize(n.get_upper_bound()); } else { - return carl::bitsize(n.get_lower_bound()) + carl::bitsize(n.get_upper_bound()) + poly::degree(n.libpoly_polynomial()); + return carl::bitsize(n.get_lower_bound()) + carl::bitsize(n.get_upper_bound()) + lp_upolynomial_degree(n.get_internal()->value.a.f); } } Sign sgn(const LPRealAlgebraicNumber& n) { - return static_cast(lp_algebraic_number_sgn(n.get_internal())); + return static_cast(lp_value_sgn(n.get_internal())); } Sign sgn(const LPRealAlgebraicNumber&, const UnivariatePolynomial&) { @@ -249,17 +181,11 @@ Sign sgn(const LPRealAlgebraicNumber&, const UnivariatePolynomial& i) { CARL_LOG_DEBUG("carl.ran.libpoly", "ran " << n << " contained in " << i); - - poly::Interval inter = to_libpoly_interval(i); - - lp_value_t val; - lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, n.get_internal()); - - bool b = lp_interval_contains(inter.get_internal(), &val); - //Destruct the value, but do not free m_content - lp_value_destruct(&val); - - return b; + lp_interval_t* inter = to_libpoly_interval(i); + bool res = lp_interval_contains(inter, n.get_internal()); + lp_interval_destruct(inter); + delete inter; + return res; } @@ -271,8 +197,10 @@ bool contained_in(const LPRealAlgebraicNumber& n, const Interval 0 && !n.is_numeric()) { - lp_algebraic_number_refine_const(n.get_internal()); + if (n.is_numeric()) return; + + while (lp_dyadic_interval_size(&n.get_internal()->value.a.I) > 0 && !n.is_numeric()) { + lp_algebraic_number_refine_const(&n.get_internal()->value.a); } CARL_LOG_DEBUG("carl.ran.libpoly", "Finished Refining Algebraic NumberType : " << n); } @@ -282,104 +210,99 @@ void LPRealAlgebraicNumber::refine() const { } NumberType branching_point(const LPRealAlgebraicNumber& n) { - //return carl::sample(n.interval_int()); - refine(n); - poly::DyadicRational res; - lp_algebraic_number_get_dyadic_midpoint(n.get_internal(), res.get_internal()); - NumberType num = to_rational(res); - return num; + if (n.is_numeric()) return n.value(); + refine(n); + lp_dyadic_rational_t res; + lp_dyadic_rational_construct(&res); + lp_algebraic_number_get_dyadic_midpoint(&n.get_internal()->value.a, &res); + auto result = to_rational(&res); + lp_dyadic_rational_destruct(&res); + return result; } NumberType sample_above(const LPRealAlgebraicNumber& n) { - //return carl::floor(n.interval_int().upper()) + 1; From Ran.h CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling above: " << n); - refine(n); + + lp_value_t inf; + lp_value_construct(&inf, LP_VALUE_PLUS_INFINITY, 0); - return carl::floor(n.get_upper_bound()) + 1; + lp_value_t res; + lp_value_construct_none(&res); + lp_value_get_value_between(n.get_internal(), true, &inf, true, &res); + auto val = to_rational(&res); + lp_value_destruct(&res); + return val; } NumberType sample_below(const LPRealAlgebraicNumber& n) { //return carl::ceil(n.interval_int().lower()) - 1; From Ran.h CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling below: " << n); - refine(n); + + lp_value_t inf; + lp_value_construct(&inf, LP_VALUE_MINUS_INFINITY, 0); - return carl::ceil(n.get_lower_bound()) - 1; + lp_value_t res; + lp_value_construct_none(&res); + lp_value_get_value_between(&inf, true, n.get_internal(), true, &res); + auto val = to_rational(&res); + lp_value_destruct(&res); + return val; } NumberType sample_between(const LPRealAlgebraicNumber& lower, const LPRealAlgebraicNumber& upper) { CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling between: " << lower << " and " << upper); - //Make sure that the intervals are disjoint - while (lower.get_upper_bound() > upper.get_lower_bound()) { - if(lower.is_numeric()){ - break ; - } - if(upper.is_numeric()){ - break ; - } - lp_algebraic_number_refine_const(lower.get_internal()); - lp_algebraic_number_refine_const(upper.get_internal()); - } - - if (lower.is_numeric()) { - CARL_LOG_DEBUG("carl.ran.libpoly", "Found sample between: " << sample_between(lower.value(), upper)); - return sample_between(lower.value(), upper); - } else if (upper.is_numeric()) { - CARL_LOG_DEBUG("carl.ran.libpoly", "Found sample between: " << sample_between(lower, upper.value())); - return sample_between(lower, upper.value()); - } else { - CARL_LOG_DEBUG("carl.ran.libpoly", "Found sample between: " << carl::sample(Interval(lower.get_upper_bound(), upper.get_lower_bound()), true)); - return carl::sample(Interval(lower.get_upper_bound(), upper.get_lower_bound()), true); - } + lp_value_t res; + lp_value_construct_none(&res); + lp_value_get_value_between(lower.get_internal(), true, upper.get_internal(), true, &res); + auto val = to_rational(&res); + lp_value_destruct(&res); + return val; } NumberType sample_between(const LPRealAlgebraicNumber& lower, const NumberType& upper) { - CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling between: " << lower << " and " << upper); - //Make sure that the intervals are disjoint - while (!(lower.get_upper_bound() < upper) && !lower.is_numeric()) { - lp_algebraic_number_refine_const(lower.get_internal()); - } + CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling between: " << lower << " and " << upper); - if (lower.is_numeric()) { - NumberType sample = sample_between(lower.value(), upper); - CARL_LOG_DEBUG("carl.ran.libpoly", "Found sample between: " << sample); - return sample; - } else { - //sample from interval - NumberType sample = carl::sample(Interval(lower.get_upper_bound(), BoundType::WEAK, upper, BoundType::STRICT), false); - CARL_LOG_DEBUG("carl.ran.libpoly", "Found sample between: " << sample); - return sample; - } -} + lp_value_t v_upper; + lp_value_construct(&v_upper, LP_VALUE_RATIONAL, upper.get_mpq_t()); + + lp_value_t res; + lp_value_construct_none(&res); + lp_value_get_value_between(lower.get_internal(), true, &v_upper, true, &res); -//Make sure that the intervals are disjoint + lp_value_destruct(&v_upper); + + auto val = to_rational(&res); + lp_value_destruct(&res); + return val; +} NumberType sample_between(const NumberType& lower, const LPRealAlgebraicNumber& upper) { - CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling between: " << lower << " and " << upper); + CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling between: " << lower << " and " << upper); - while (!(lower < upper.get_lower_bound()) && !upper.is_numeric()) { - lp_algebraic_number_refine_const(upper.get_internal()); - } + lp_value_t v_lower; + lp_value_construct(&v_lower, LP_VALUE_RATIONAL, lower.get_mpq_t()); + lp_value_t res; + lp_value_construct_none(&res); + lp_value_get_value_between(&v_lower, true, upper.get_internal(), true, &res); - if (upper.is_numeric()) { - CARL_LOG_DEBUG("carl.ran.libpoly", "Found sample between: " << sample_between(lower, upper.value())); - return sample_between(lower, upper.value()); - } else { - CARL_LOG_DEBUG("carl.ran.libpoly", "Found sample between: " << carl::sample(Interval(lower, BoundType::STRICT, upper.get_lower_bound(), BoundType::WEAK), false)); - return carl::sample(Interval(lower, BoundType::STRICT, upper.get_lower_bound(), BoundType::WEAK), false); - } + lp_value_destruct(&v_lower); + + auto val = to_rational(&res); + lp_value_destruct(&res); + return val; } NumberType floor(const LPRealAlgebraicNumber& n) { CARL_LOG_DEBUG("carl.ran.libpoly", "Floor of: " << n); - refine(n); lp_integer_t val; - lp_algebraic_number_floor(n.get_internal(), &val); - NumberType ret = to_rational(*poly::detail::cast_from(&val)); + lp_integer_construct(&val); + lp_value_floor(n.get_internal(), &val); + NumberType ret = to_rational(&val); lp_integer_destruct(&val) ; return ret ; } @@ -387,19 +310,33 @@ NumberType floor(const LPRealAlgebraicNumber& n) { NumberType ceil(const LPRealAlgebraicNumber& n) { CARL_LOG_DEBUG("carl.ran.libpoly", "Ceil of: " << n); - refine(n); lp_integer_t val; - lp_algebraic_number_ceiling(n.get_internal(), &val); - NumberType ret = to_rational(*poly::detail::cast_from(&val)); + lp_integer_construct(&val); + lp_value_ceiling(n.get_internal(), &val); + NumberType ret = to_rational(&val); lp_integer_destruct(&val) ; return ret ; } - bool compare(const LPRealAlgebraicNumber& lhs, const LPRealAlgebraicNumber& rhs, const Relation relation) { - int cmp = lp_algebraic_number_cmp(lhs.get_internal(), rhs.get_internal()); + if (lhs.m_content == rhs.m_content) { + switch (relation) { + case Relation::EQ: + case Relation::LEQ: + case Relation::GEQ: + return true; + case Relation::NEQ: + case Relation::LESS: + case Relation::GREATER: + default: + return false; + } + } + + int cmp = lp_value_cmp(lhs.get_internal(), rhs.get_internal()); switch (relation) { case Relation::EQ: + if (cmp == 0) lhs.m_content = rhs.m_content; return cmp == 0; case Relation::NEQ: return cmp != 0; @@ -422,9 +359,7 @@ const carl::Variable LPRealAlgebraicNumber::auxVariable = fresh_real_variable("_ bool compare(const LPRealAlgebraicNumber& lhs, const NumberType& rhs, const Relation relation) { - poly::Rational rat = to_libpoly_rational(rhs); - - int cmp = lp_algebraic_number_cmp_rational(lhs.get_internal(), rat.get_internal()); + int cmp = lp_value_cmp_rational(lhs.get_internal(), (lp_rational_t*)rhs.get_mpq_t()); switch (relation) { case Relation::EQ: @@ -447,7 +382,7 @@ bool compare(const LPRealAlgebraicNumber& lhs, const NumberType& rhs, const Rela std::ostream& operator<<(std::ostream& os, const LPRealAlgebraicNumber& ran) { - char* str = lp_algebraic_number_to_string(ran.get_internal()); + char* str = lp_value_to_string(ran.get_internal()); os << str; free(str); return os; diff --git a/src/carl-arith/ran/libpoly/LPRan.h b/src/carl-arith/ran/libpoly/LPRan.h index 4fc4e397..5ba37c93 100644 --- a/src/carl-arith/ran/libpoly/LPRan.h +++ b/src/carl-arith/ran/libpoly/LPRan.h @@ -25,12 +25,16 @@ class LPRealAlgebraicNumber { using NumberType = mpq_class; private: - lp_algebraic_number_t m_content; + struct Content { + lp_value_t c; + ~Content() { + lp_value_destruct(&c); + } + }; + mutable std::shared_ptr m_content; static const Variable auxVariable; - - public: ~LPRealAlgebraicNumber(); @@ -43,25 +47,14 @@ class LPRealAlgebraicNumber { * Construct from libpoly Algebraic NumberType and takes ownership * @param num, LibPoly Algebraic Number */ - LPRealAlgebraicNumber(const poly::AlgebraicNumber& num); + LPRealAlgebraicNumber(const lp_value_t& num); - /** - * Construct from libpoly Algebraic NumberType and takes ownership - * @param num, LibPoly Algebraic Number - */ - LPRealAlgebraicNumber(const lp_algebraic_number_t& num); - - /** - * Move from libpoly Algebraic NumberType (C++ Interface) - * @param num, LibPoly Algebraic Number - */ - LPRealAlgebraicNumber(poly::AlgebraicNumber&& num); /** * Move from libpoly Algebraic NumberType (C Interface) * @param num, LibPoly Algebraic Number */ - LPRealAlgebraicNumber(lp_algebraic_number_t&& num); + LPRealAlgebraicNumber(lp_value_t&& num); /** * Construct from Polynomial and Interval @@ -87,27 +80,18 @@ class LPRealAlgebraicNumber { */ static LPRealAlgebraicNumber create_safe(const carl::UnivariatePolynomial& p, const Interval& i); - /** - * Convert a libpoly Value into an algebraic number - * This does not free the value - * In case the value is none or infinity this conversion is not possible - * @param Libpoly Value (C interface) - * @return Copy of val as a algebraic number - */ - static LPRealAlgebraicNumber create_from_value(const lp_value_t* val) ; - /** * @return Pointer to the internal libpoly algebraic number (C interface) */ - inline lp_algebraic_number_t* get_internal() { - return &m_content; + inline lp_value_t* get_internal() { + return &(m_content->c); } /** * @return Pointer to the internal libpoly algebraic number (C interface) */ - inline const lp_algebraic_number_t* get_internal() const { - return &m_content; + inline const lp_value_t* get_internal() const { + return &(m_content->c); } /** @@ -115,9 +99,6 @@ class LPRealAlgebraicNumber { */ bool is_numeric() const; - const poly::UPolynomial libpoly_polynomial() const; - - const poly::DyadicInterval& libpoly_interval() const; const UnivariatePolynomial polynomial() const; @@ -164,7 +145,7 @@ void refine(const LPRealAlgebraicNumber& n); inline bool is_zero(const LPRealAlgebraicNumber& n) { refine(n); - return lp_algebraic_number_sgn(n.get_internal()) == 0; + return lp_value_sgn(n.get_internal()) == 0; } bool is_integer(const LPRealAlgebraicNumber& n); @@ -187,7 +168,7 @@ template<> struct hash { std::size_t operator()(const carl::LPRealAlgebraicNumber& n) const { //Todo test if the precisions needs to be adjusted - return lp_algebraic_number_hash_approx(n.get_internal(), 0); + return lp_value_hash_approx(n.get_internal(), 0); } }; } // namespace std diff --git a/src/carl-arith/ran/libpoly/RealRoots.cpp b/src/carl-arith/ran/libpoly/RealRoots.cpp index 5d47a8fc..169f70aa 100644 --- a/src/carl-arith/ran/libpoly/RealRoots.cpp +++ b/src/carl-arith/ran/libpoly/RealRoots.cpp @@ -2,6 +2,8 @@ #ifdef USE_LIBPOLY +#include "LPAssignment.h" + namespace carl { @@ -10,9 +12,8 @@ RealRootsResult real_roots( const Interval& interval) { CARL_LOG_DEBUG("carl.ran.libpoly", " Real roots of " << polynomial << " within " << interval); - assert(poly::is_univariate(polynomial.get_polynomial())); + assert(carl::is_univariate(polynomial)); - // Easy checks if (carl::is_zero(polynomial)) { CARL_LOG_TRACE("carl.ran.libpoly", "poly is 0 -> nullified"); return RealRootsResult::nullified_response(); @@ -21,30 +22,39 @@ RealRootsResult real_roots( return RealRootsResult::no_roots_response(); } - poly::Interval inter_poly = to_libpoly_interval(interval); + lp_upolynomial_t* upoly = lp_polynomial_to_univariate(polynomial.get_internal()); - // actual calculations - std::vector roots = poly::isolate_real_roots(poly::to_univariate(polynomial.get_polynomial())); + lp_algebraic_number_t* roots = new lp_algebraic_number_t[lp_upolynomial_degree(upoly)]; + std::size_t roots_size; + lp_upolynomial_roots_isolate(upoly, roots, &roots_size); - if (roots.empty()) { + if (roots_size == 0) { + delete[] roots; + lp_upolynomial_delete(upoly); CARL_LOG_DEBUG("carl.ran.libpoly", "Poly has no roots"); return RealRootsResult::no_roots_response(); } - // sort roots in ascending order - std::sort(roots.begin(), roots.end(), std::less()); - - // turn into RealRootsResult + lp_interval_t* inter_poly = to_libpoly_interval(interval); std::vector res; - for (const auto& val : roots) { - auto tmp = LPRealAlgebraicNumber(val); - // filter out roots not in interval - if (poly::contains(inter_poly, poly::Value(val))) { - CARL_LOG_DEBUG("carl.ran.libpoly", " Found Root " << val); - res.emplace_back(tmp); + for (std::size_t i = 0; i < roots_size; ++i) { + lp_value_t tmp; + lp_value_construct(&tmp, LP_VALUE_ALGEBRAIC, &roots[i]); + if (lp_interval_contains(inter_poly, &tmp)) { + res.emplace_back(std::move(tmp)); + } else { + lp_value_destruct(&tmp); } + lp_algebraic_number_destruct(&roots[i]); } + lp_interval_destruct(inter_poly); + delete inter_poly; + + delete[] roots; + lp_upolynomial_delete(upoly); + + std::sort(res.begin(), res.end()); return RealRootsResult::roots_response(std::move(res)); } @@ -54,7 +64,7 @@ RealRootsResult real_roots( const Interval& interval) { CARL_LOG_DEBUG("carl.ran.libpoly", polynomial << " " << m << " " << interval); - if (poly::is_univariate(polynomial.get_polynomial())) { + if (carl::is_univariate(polynomial)) { return real_roots(polynomial, interval); } @@ -67,55 +77,65 @@ RealRootsResult real_roots( return RealRootsResult::no_roots_response(); } - poly::Interval inter_poly = to_libpoly_interval(interval); + for (const auto& v : carl::variables(polynomial)) { + if (v != polynomial.main_var() && m.find(v) == m.end()) return RealRootsResult::non_univariate_response(); + } // Multivariate Polynomial // build the assignment - poly::Assignment assignment(polynomial.context().poly_context()); Variable mainVar = polynomial.main_var(); - for (Variable& var : carl::variables(polynomial)) { - if (var == mainVar) continue; - // We convert numbers to libpoly values and add to assignment so we can substitute them later using libpoly - lp_value_t val; - // Turn into value - lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, m.at(var).get_internal()); - // That copies the value into the assignment - assignment.set(*(polynomial.context().lp_variable(var)), poly::Value(&val)); - // Destroy the value, but dont free the algebraic number! - lp_value_destruct(&val); - } + auto evalMap = m; + evalMap.erase(mainVar); + lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); CARL_LOG_TRACE("carl.ran.libpoly", "Call libpoly"); - std::vector roots = poly::isolate_real_roots(polynomial.get_polynomial(), assignment); - CARL_LOG_TRACE("carl.ran.libpoly", "Libpoly returned"); + lp_value_t* roots = new lp_value_t[lp_polynomial_degree(polynomial.get_internal())]; + std::size_t roots_size; + lp_polynomial_roots_isolate(polynomial.get_internal(), &assignment, roots, &roots_size); - if (roots.empty()) { + CARL_LOG_TRACE("carl.ran.libpoly", "Libpoly returned"); + if (roots_size == 0) { + delete[] roots; CARL_LOG_DEBUG("carl.ran.libpoly", " Checking for nullification -> Evaluation at " << mainVar << "= 1"); - assignment.set(*(polynomial.context().lp_variable(mainVar)), poly::Value(long(1))); - poly::Value eval_val = poly::evaluate(polynomial.get_polynomial(), assignment); - CARL_LOG_DEBUG("carl.ran.libpoly", " Got eval_val " << eval_val); - - if (eval_val == poly::Value(long(0))) { + lp_value_t val; + // here we know that no value for mainVar has been set, so we can safely set it! + lp_value_construct_int(&val, long(1)); + lp_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(mainVar), &val); + lp_value_destruct(&val); + auto eval_val = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); + lp_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(mainVar), 0); + //CARL_LOG_DEBUG("carl.ran.libpoly", " Got eval_val " << eval_val); + + lp_value_t zero_value; + lp_value_construct_zero(&zero_value); + if (lp_value_cmp(eval_val, &zero_value) == 0) { + lp_value_destruct(&zero_value); CARL_LOG_DEBUG("carl.ran.libpoly", "poly is 0 after substituting rational assignments -> nullified"); + lp_value_delete(eval_val); return RealRootsResult::nullified_response(); } else { + lp_value_destruct(&zero_value); CARL_LOG_DEBUG("carl.ran.libpoly", "Poly has no roots"); + lp_value_delete(eval_val); return RealRootsResult::no_roots_response(); } } - std::sort(roots.begin(), roots.end(), std::less()); - - // turn result into RealRootsResult + lp_interval_t* inter_poly = to_libpoly_interval(interval); std::vector res; - for (const poly::Value& val : roots) { - auto tmp = LPRealAlgebraicNumber::create_from_value(val.get_internal()); - // filter out roots not in interval - if (poly::contains(inter_poly, val)) { - CARL_LOG_DEBUG("carl.ran.libpoly", " Found root " << val); - res.emplace_back(tmp); + for (std::size_t i = 0; i < roots_size; ++i) { + if (lp_interval_contains(inter_poly, &roots[i])) { + res.emplace_back(std::move(roots[i])); + CARL_LOG_DEBUG("carl.ran.libpoly", " Found root " << res.back()); + } else { + lp_value_destruct(&roots[i]); } } + lp_interval_destruct(inter_poly); + delete inter_poly; + delete[] roots; + + std::sort(res.begin(), res.end()); return RealRootsResult::roots_response(std::move(res)); } diff --git a/src/carl-arith/ran/libpoly/helper.h b/src/carl-arith/ran/libpoly/helper.h index 33439120..a3b91daa 100644 --- a/src/carl-arith/ran/libpoly/helper.h +++ b/src/carl-arith/ran/libpoly/helper.h @@ -1,7 +1,5 @@ #pragma once -// TODO cleanup, move to ran folder - #include #ifdef USE_LIBPOLY @@ -18,84 +16,47 @@ #include #include -#include +#include namespace carl { -inline carl::UnivariatePolynomial to_carl_univariate_polynomial(const poly::UPolynomial& p, const carl::Variable& main_var){ - CARL_LOG_FUNC("carl.converter", p << ", " << main_var); - std::vector carl_coeffs; - for (const poly::Integer& coeff : poly::coefficients(p)) { - carl_coeffs.emplace_back(*(poly::detail::cast_to_gmp(&coeff))); - } - CARL_LOG_TRACE("carl.converter", "-> coefficients: " << carl_coeffs); - return carl::UnivariatePolynomial(main_var, carl_coeffs); +inline mpz_class to_integer(const lp_integer_t* m) { + return *reinterpret_cast(m); } -template> = dummy> -inline poly::UPolynomial to_libpoly_upolynomial(const carl::UnivariatePolynomial& p) { - CARL_LOG_FUNC("carl.converter", p); - - mpz_class denominator; - - if (carl::is_constant(p)) { - CARL_LOG_TRACE("carl.converter", "Poly is constant"); - denominator = carl::get_denom(p.lcoeff()); - CARL_LOG_TRACE("carl.converter", "Coprime Factor/ Denominator: " << denominator); - return poly::UPolynomial(poly::Integer(carl::get_num(p.lcoeff()))); - } - Coeff coprimeFactor = p.coprime_factor(); - - CARL_LOG_TRACE("carl.converter", "Coprime Factor: " << coprimeFactor); - - if (carl::get_denom(coprimeFactor) != 1) { - // if coprime factor is not an integer - denominator = coprimeFactor > 0 ? mpz_class(carl::get_num(coprimeFactor)) : mpz_class(-carl::get_num(coprimeFactor)); - } else if (coprimeFactor == 0) { - // Wie kann das überhaupt sein? TODO - denominator = mpz_class(1); - } else { - denominator = mpz_class(coprimeFactor); - } - CARL_LOG_TRACE("carl.converter", "Coprime factor/ denominator: " << denominator); - - std::vector coefficients; - for (const auto& c : p.coefficients()) { - coefficients.emplace_back(mpz_class(c * denominator)); - } - CARL_LOG_TRACE("carl.converter", "-> coefficients: " << coefficients); - return poly::UPolynomial(coefficients); -} - -inline mpq_class to_rational(const poly::Integer& m) { - return *poly::detail::cast_to_gmp(&m); -} - -inline mpq_class to_rational(const poly::Rational& m) { - return *poly::detail::cast_to_gmp(&m); +inline mpq_class to_rational(const lp_integer_t* m) { + return mpq_class(to_integer(m)); } inline mpq_class to_rational(const lp_rational_t* m) { return *reinterpret_cast(m); } +inline mpq_class to_rational(const lp_dyadic_rational_t* m) { + lp_integer_t num; + lp_integer_construct(&num); + lp_dyadic_rational_get_num(m, &num); -inline mpz_class to_integer(const poly::Integer& m) { - return *poly::detail::cast_to_gmp(&m); -} + lp_integer_t den; + lp_integer_construct(&den); + lp_dyadic_rational_get_den(m, &den); + + auto res = mpq_class(to_integer(&num), to_integer(&den)); -inline mpq_class to_rational(const poly::DyadicRational& m) { - return mpq_class(to_integer(poly::numerator(m)), to_integer(poly::denominator(m))); + lp_integer_destruct(&num); + lp_integer_destruct(&den); + + return res; } -inline mpq_class to_rational(const poly::Value& m){ - switch(m.get_internal()->type){ +inline mpq_class to_rational(const lp_value_t* m){ + switch(m->type){ case lp_value_type_t::LP_VALUE_INTEGER: - return to_rational(poly::as_integer(m)); - case lp_value_type_t::LP_VALUE_DYADIC_RATIONAL: - return to_rational(poly::as_dyadic_rational(m)); + return to_rational(&m->value.z); case lp_value_type_t::LP_VALUE_RATIONAL: - return to_rational(poly::as_rational(m)); + return to_rational(&m->value.q); + case lp_value_type_t::LP_VALUE_DYADIC_RATIONAL: + return to_rational(&m->value.dy_q); default: CARL_LOG_ERROR("carl.converter", "Cannot convert libpoly value: " << m << " to rational."); assert(false); @@ -103,105 +64,107 @@ inline mpq_class to_rational(const poly::Value& m){ } } -inline poly::Rational to_libpoly_rational(const mpq_class& num) { - return poly::Rational(num); -} +inline carl::UnivariatePolynomial to_carl_univariate_polynomial(const lp_upolynomial_t* p, const carl::Variable& main_var){ + CARL_LOG_FUNC("carl.converter", p << ", " << main_var); + + lp_integer_t coeffs[lp_upolynomial_degree(p) + 1]; + for (std::size_t i = 0; i < lp_upolynomial_degree(p) + 1; ++i) { + lp_integer_construct(&coeffs[i]); + } + lp_upolynomial_unpack(p, coeffs); -//Exact for whole numbers -inline poly::DyadicRational get_libpoly_dyadic_approximation(const mpz_class& num) { - return poly::DyadicRational(poly::Integer(num)); + std::vector carl_coeffs; + for (std::size_t i = 0; i < lp_upolynomial_degree(p) + 1; ++i) { + carl_coeffs.emplace_back(to_rational(&coeffs[i])); + lp_integer_destruct(&coeffs[i]); + } + + CARL_LOG_TRACE("carl.converter", "-> coefficients: " << carl_coeffs); + return carl::UnivariatePolynomial(main_var, carl_coeffs); } -/** - * Tries to convert num = a/b into a dyadic rational of the form c/2^d - * We take d = precision * ceil(log(2,b)) and c = ceil((a * 2^d)/(b)) - * @return Approximation to num by a dyadic rational - */ -inline poly::DyadicRational get_libpoly_dyadic_approximation(const mpq_class& num, const unsigned int& precision) { - CARL_LOG_DEBUG("carl.converter", "Starting Dyadic Approximation with: " << num); - - mpz_class a = num.get_num(); - mpz_class b = num.get_den(); - mpz_class d = (precision)*carl::ceil(carl::log(b)); //unsigned long - assert(d.fits_ulong_p()); - - mpz_class c = carl::ceil((a * carl::pow(mpq_class(2), d.get_ui())) / (b)); +template> = dummy> +inline lp_upolynomial_t* to_libpoly_upolynomial(const carl::UnivariatePolynomial& p) { + CARL_LOG_FUNC("carl.converter", p); - if (!c.fits_slong_p()) { - CARL_LOG_DEBUG("carl.converter", "Dyadic Rational: Fallback for overflow of numerator"); - poly::DyadicRational tmp = get_libpoly_dyadic_approximation(c); //Construct with arbitrary precision - return poly::div_2exp(tmp, d.get_ui()); - } + if (carl::is_constant(p)) { + CARL_LOG_TRACE("carl.converter", "Poly is constant"); + assert(carl::get_denom(p.lcoeff())>0); + lp_upolynomial_t* res = lp_upolynomial_construct(lp_Z, 0, mpz_class(carl::get_num(p.lcoeff())).get_mpz_t()); + return res; + } - CARL_LOG_DEBUG("carl.converter", "Got d: " << d << " and c " << c); - CARL_LOG_DEBUG("carl.converter", "Got dyadic number: " << poly::DyadicRational(c.get_si(), d.get_ui())); + Coeff coprimeFactor = p.coprime_factor(); + CARL_LOG_TRACE("carl.converter", "Coprime Factor: " << coprimeFactor); + mpz_class denominator; + if (carl::get_denom(coprimeFactor) != 1) { + // if coprime factor is not an integer + denominator = coprimeFactor > 0 ? mpz_class(carl::get_num(coprimeFactor)) : mpz_class(-carl::get_num(coprimeFactor)); + } else if (coprimeFactor == 0) { + // Wie kann das überhaupt sein? TODO + denominator = mpz_class(1); + } else { + denominator = mpz_class(coprimeFactor); + } + CARL_LOG_TRACE("carl.converter", "Coprime factor/ denominator: " << denominator); - return poly::DyadicRational(c.get_si(), d.get_ui()); + std::vector coefficients; + std::vector coefficients_ptr; + for (const auto& c : p.coefficients()) { + coefficients.emplace_back(c * denominator); + coefficients_ptr.emplace_back(*coefficients.back().get_mpz_t()); + } + CARL_LOG_TRACE("carl.converter", "coefficients: " << coefficients); + lp_upolynomial_t* res = lp_upolynomial_construct(lp_Z, p.coefficients().size()-1, reinterpret_cast(coefficients_ptr.data())); + return res; } /** * Convert a carl interval to a libpoly interval * This keeps the bound types */ -inline poly::Interval to_libpoly_interval(const carl::Interval& inter) { - poly::Value low; +inline lp_interval_t* to_libpoly_interval(const carl::Interval& inter) { + lp_value_t low; bool low_open; switch (inter.lower_bound_type()) { case BoundType::STRICT: - low = poly::Value(poly::Rational(inter.lower())); + lp_value_construct(&low, LP_VALUE_RATIONAL, &inter.lower()); low_open = true; break; case BoundType::WEAK: - low = poly::Value(poly::Rational(inter.lower())); + lp_value_construct(&low, LP_VALUE_RATIONAL, &inter.lower()); low_open = false; break; case BoundType::INFTY: - low = poly::Value::minus_infty(); + default: + low = *lp_value_minus_infinity(); low_open = true; break; } - poly::Value up; + lp_value_t up; bool up_open; switch (inter.upper_bound_type()) { case BoundType::STRICT: - up = poly::Value(poly::Rational(inter.upper())); + lp_value_construct(&up, LP_VALUE_RATIONAL, &inter.upper()); up_open = true; break; case BoundType::WEAK: - up = poly::Value(poly::Rational(inter.upper())); + lp_value_construct(&up, LP_VALUE_RATIONAL, &inter.upper()); up_open = false; break; case BoundType::INFTY: - up = poly::Value::plus_infty(); + default: + up = *lp_value_plus_infinity(); up_open = false; break; } - return poly::Interval(low, low_open, up, up_open); -} - -/** - * Converts a libpoly dyadic interval to a carl interval - * Keeps the bound types - */ -template -inline carl::Interval to_carl_interval(const poly::DyadicInterval& inter) { - BoundType upper_bound = inter.get_internal()->a_open ? BoundType::STRICT : BoundType::WEAK ; - BoundType lower_bound = inter.get_internal()->b_open ? BoundType::STRICT : BoundType::WEAK ; - - return carl::Interval(to_rational(poly::get_lower(inter)), lower_bound, to_rational(poly::get_upper(inter)), upper_bound); -} - -/** - * Converts a libpoly interval to a carl interval - * Keeps the bound types - */ -template -inline carl::Interval to_carl_interval(const poly::Interval& inter) { - BoundType upper_bound = inter.get_internal()->a_open ? BoundType::STRICT : BoundType::WEAK ; - BoundType lower_bound = inter.get_internal()->b_open ? BoundType::STRICT : BoundType::WEAK ; - return carl::Interval(to_rational(poly::get_lower(inter)), lower_bound, to_rational(poly::get_upper(inter)), upper_bound); + lp_interval_t* res = new lp_interval_t; + lp_interval_construct(res, &low, low_open, &up, up_open); + lp_value_destruct(&low); + lp_value_destruct(&up); + return res; } } // namespace carl diff --git a/src/carl-statistics/MultiCounter.h b/src/carl-statistics/MultiCounter.h index fc9e365b..b97731e4 100644 --- a/src/carl-statistics/MultiCounter.h +++ b/src/carl-statistics/MultiCounter.h @@ -7,10 +7,12 @@ namespace carl::statistics { template class MultiCounter { boost::container::flat_map m_data; + std::size_t m_total = 0; public: void inc(const T& key, std::size_t inc) { m_data.try_emplace(key).first->second += inc; + m_total += inc; } void collect(std::map& data, const std::string& key) const { @@ -20,6 +22,7 @@ class MultiCounter { ss << "=" << v << ";"; } data.emplace(key, ss.str()); + data.emplace(key + ".total", std::to_string(m_total)); } }; diff --git a/src/carl-statistics/Statistics.h b/src/carl-statistics/Statistics.h index 04be9ef0..efa9b163 100644 --- a/src/carl-statistics/Statistics.h +++ b/src/carl-statistics/Statistics.h @@ -3,7 +3,7 @@ #include #include #include -#include +#include #include "StatisticsCollector.h" #include "Serialization.h" @@ -32,7 +32,7 @@ class Statistics { mCollected.emplace(key, value); } - void addKeyValuePair(const std::string& key, const Timer& value) { + void addKeyValuePair(const std::string& key, Timer& value) { value.collect(mCollected, key); } diff --git a/src/carl-statistics/Timing.h b/src/carl-statistics/Timing.h index 9dc8d124..7be4d06a 100644 --- a/src/carl-statistics/Timing.h +++ b/src/carl-statistics/Timing.h @@ -30,7 +30,7 @@ namespace timing { class Timer { std::size_t m_count = 0; timing::duration m_overall = timing::zero(); - timing::time_point m_current_start; + timing::time_point m_current_start = timing::time_point::min(); public: static timing::time_point start() { @@ -45,7 +45,15 @@ class Timer { } void finish() { finish(m_current_start); + m_current_start = timing::time_point::min(); } + bool check_finish() { + if (m_current_start != timing::time_point::min()) { + finish(); + return true; + } + return false; + } auto count() const { return m_count; } @@ -53,9 +61,11 @@ class Timer { return m_overall.count(); } - void collect(std::map& data, const std::string& key) const { + void collect(std::map& data, const std::string& key) { + bool active_at_timeout = check_finish(); data.emplace(key+".count", std::to_string(count())); data.emplace(key+".overall_ms", std::to_string(overall_ms())); + data.emplace(key+".active_at_timeout", active_at_timeout ? "1" : "0"); } }; diff --git a/src/carl-vs/substitute.h b/src/carl-vs/substitute.h index 68fa18e9..7f6d28ea 100644 --- a/src/carl-vs/substitute.h +++ b/src/carl-vs/substitute.h @@ -419,8 +419,8 @@ namespace carl::vs { assert(zeros.size() == 1); // calculate subVar1-subVar2 ~ 0 [term//subVar1][zero//subVar2] - carl::Variable subVar1 = carl::fresh_real_variable("__subVar1"); - carl::Variable subVar2 = carl::fresh_real_variable("__subVar2"); + static carl::Variable subVar1 = carl::fresh_real_variable("__subVar1"); + static carl::Variable subVar2 = carl::fresh_real_variable("__subVar2"); auto subRes1 = substitute(Constraint(Poly(subVar1) - subVar2, varcompRelation), subVar1, term); assert(subRes1); CaseDistinction result; @@ -452,7 +452,7 @@ namespace carl::vs { } return result; } else if(term.is_minus_infty() ) { - carl::Variable subVar1 = carl::fresh_real_variable("__subVar1"); + static carl::Variable subVar1 = carl::fresh_real_variable("__subVar1"); return substitute(Constraint(subVar1, varcompRelation), subVar1, Term::minus_infty()); } else { return std::nullopt; diff --git a/src/tests/carl-arith-poly/Test_LPPolynomialFunctions.cpp b/src/tests/carl-arith-poly/Test_LPPolynomialFunctions.cpp index ac4a623e..35884cb3 100644 --- a/src/tests/carl-arith-poly/Test_LPPolynomialFunctions.cpp +++ b/src/tests/carl-arith-poly/Test_LPPolynomialFunctions.cpp @@ -158,9 +158,8 @@ TEST(LPPOLYNOMIAL, factorization) { startTime = std::chrono::high_resolution_clock::now(); - CoCoAAdaptorLP adaptor(context); for (const auto& pol : polys) { - const auto& factors = carl::factorization(pol, adaptor); + const auto& factors = carl::factorization(pol); LPPolynomial productOfFactors = LPPolynomial(context, x, (mpz_class)1); for (const auto& factor : factors) { EXPECT_NE((unsigned)0, factor.second); diff --git a/src/tests/carl-arith-ran/Test_Libpoly.cpp b/src/tests/carl-arith-ran/Test_Libpoly.cpp index 9464db84..50e7e482 100644 --- a/src/tests/carl-arith-ran/Test_Libpoly.cpp +++ b/src/tests/carl-arith-ran/Test_Libpoly.cpp @@ -14,26 +14,26 @@ using namespace carl; -TEST(LIBPOLY, convertToLibpolyUnivariate) { - auto x = fresh_real_variable("x"); - - UnivariatePolynomial carl_pol1(x, {Rational(-6), Rational(-5), Rational(1)}); - poly::UPolynomial lp_pol1({-6, -5, 1}); - - EXPECT_EQ(lp_pol1, to_libpoly_upolynomial(carl_pol1)); - - UnivariatePolynomial carl_pol2(x, {mpq_class(2, 13)}); - to_libpoly_upolynomial(carl_pol2); -} - -TEST(LIBPOLY, convertToCarlUnivariate) { - auto x = fresh_real_variable("x"); - - UnivariatePolynomial carl_pol1(x, {Rational(-6), Rational(-5), Rational(1)}); - poly::UPolynomial lp_pol1({-6, -5, 1}); - - EXPECT_EQ(carl_pol1, to_carl_univariate_polynomial(lp_pol1, x)); -} +// TEST(LIBPOLY, convertToLibpolyUnivariate) { +// auto x = fresh_real_variable("x"); +// +// UnivariatePolynomial carl_pol1(x, {Rational(-6), Rational(-5), Rational(1)}); +// poly::UPolynomial lp_pol1({-6, -5, 1}); +// +// EXPECT_EQ(lp_pol1, to_libpoly_upolynomial(carl_pol1)); +// +// UnivariatePolynomial carl_pol2(x, {mpq_class(2, 13)}); +// to_libpoly_upolynomial(carl_pol2); +// } + +// TEST(LIBPOLY, convertToCarlUnivariate) { +// auto x = fresh_real_variable("x"); +// +// UnivariatePolynomial carl_pol1(x, {Rational(-6), Rational(-5), Rational(1)}); +// poly::UPolynomial lp_pol1({-6, -5, 1}); +// +// EXPECT_EQ(carl_pol1, to_carl_univariate_polynomial(lp_pol1, x)); +// } TEST(LIBPOLY, convertRan) { lp_trace_enable("coefficient::arith");