From 7c752a5ebfbb3b17d9bbd4ab6a6bb72ecd3e24d2 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 1 Mar 2024 12:25:14 +0100 Subject: [PATCH 01/29] more efficient libpoly adaption --- src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h | 5 +- src/carl-arith/poly/libpoly/LPContext.h | 70 +++++++------------- src/carl-arith/poly/libpoly/LPPolynomial.cpp | 30 +++------ src/carl-arith/poly/libpoly/LPVariables.cpp | 40 +++++++++++ src/carl-arith/poly/libpoly/LPVariables.h | 44 ++++++++++++ src/carl-arith/ran/libpoly/Evaluation.cpp | 61 ++++++++++------- src/carl-arith/ran/libpoly/RealRoots.cpp | 45 ++++++++----- 7 files changed, 188 insertions(+), 107 deletions(-) create mode 100644 src/carl-arith/poly/libpoly/LPVariables.cpp create mode 100644 src/carl-arith/poly/libpoly/LPVariables.h diff --git a/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h b/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h index 4d7b9be1..b5227269 100644 --- a/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h +++ b/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h @@ -130,9 +130,8 @@ class CoCoAAdaptorLP { 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); + lp_variable_t var = mContext.lp_variable(mSymbolBack[i]); + poly::Polynomial polyVar = poly_helper::construct_poly(mContext.lp_context(), var); termPoly *= poly::pow(polyVar, (unsigned int)exponents[i]); } temPoly += termPoly; diff --git a/src/carl-arith/poly/libpoly/LPContext.h b/src/carl-arith/poly/libpoly/LPContext.h index 960ee733..b6fde9fe 100644 --- a/src/carl-arith/poly/libpoly/LPContext.h +++ b/src/carl-arith/poly/libpoly/LPContext.h @@ -13,28 +13,25 @@ #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) : 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); + } + ~Data() { + lp_variable_order_detach(lp_var_order); + lp_polynomial_context_detach(lp_context); + } }; std::shared_ptr m_data; @@ -55,17 +52,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 +69,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 { diff --git a/src/carl-arith/poly/libpoly/LPPolynomial.cpp b/src/carl-arith/poly/libpoly/LPPolynomial.cpp index 55c92335..1991067c 100644 --- a/src/carl-arith/poly/libpoly/LPPolynomial.cpp +++ b/src/carl-arith/poly/libpoly/LPPolynomial.cpp @@ -73,7 +73,7 @@ LPPolynomial::LPPolynomial(const LPContext& context, const mpq_class& val) : LPP 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_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())); @@ -81,7 +81,7 @@ LPPolynomial::LPPolynomial(const LPContext& context, const Variable& var, const 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_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())); @@ -91,13 +91,12 @@ LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, co : m_poly(context.lp_context()), m_context(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); + lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(coeff).get_internal(), var, (unsigned int)pow); m_poly += temp; } //lp_polynomial_set_external(m_poly.get_internal()); @@ -107,14 +106,13 @@ LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, co : m_poly(context.lp_context()), m_context(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); + lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(coeff).get_internal(), var, (unsigned int)pow); m_poly += temp; } //lp_polynomial_set_external(m_poly.get_internal()); @@ -124,13 +122,12 @@ LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, st : m_poly(context.lp_context()), m_context(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); + lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(std::move(coeff)).get_internal(), var, (unsigned int)pow); m_poly += temp; } //lp_polynomial_set_external(m_poly.get_internal()); @@ -140,11 +137,10 @@ LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, co : m_poly(context.lp_context()), m_context(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); + lp_polynomial_construct_simple(temp.get_internal(), context.lp_context(), poly::Integer(coef.second).get_internal(), var, coef.first); m_poly += temp; } //lp_polynomial_set_external(m_poly.get_internal()); @@ -154,7 +150,7 @@ 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); + 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); @@ -393,9 +389,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 +442,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; @@ -529,9 +521,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); diff --git a/src/carl-arith/poly/libpoly/LPVariables.cpp b/src/carl-arith/poly/libpoly/LPVariables.cpp new file mode 100644 index 00000000..9aabaf37 --- /dev/null +++ b/src/carl-arith/poly/libpoly/LPVariables.cpp @@ -0,0 +1,40 @@ +#include "LPVariables.h" + + +#include +#ifdef USE_LIBPOLY + + +namespace carl { + +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) << " -> " << it->second); + 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 << " -> " << lp_variable_db_get_name(lp_var_db, 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()) { + 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); + return poly_var; + } + // CARL_LOG_TRACE("carl.poly", "Mapping carl variable " << var << " -> " << lp_variable_db_get_name(lp_var_db, it->second)); +} + +} // 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..7f935c08 --- /dev/null +++ b/src/carl-arith/poly/libpoly/LPVariables.h @@ -0,0 +1,44 @@ +#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; + + LPVariables() { + std::cout << "init vardb" << std::endl; + lp_var_db = lp_variable_db_new(); + } + + ~LPVariables() { + lp_variable_db_detach(lp_var_db); + } + + 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/ran/libpoly/Evaluation.cpp b/src/carl-arith/ran/libpoly/Evaluation.cpp index 6786deeb..80a29d74 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -8,26 +8,27 @@ std::optional evaluate( const LPPolynomial& polynomial, const std::map& evalMap) { - //Turn into poly::Assignment - poly::Assignment assignment(polynomial.context().poly_context()); + lp_assignment_t assignment; + lp_assignment_construct(&assignment, LPVariables::getInstance().lp_var_db); 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_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(entry.first), &val); lp_value_destruct(&val); } + auto result = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); + lp_assignment_destruct(&assignment); - 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::create_from_value(result); + CARL_LOG_DEBUG("carl.ran.libpoly", " Result: " << ran); + lp_value_destruct(result); + free(result); + return ran; } boost::tribool evaluate(const BasicConstraint& constraint, const std::map& evalMap) { @@ -59,34 +60,46 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s //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()); + lp_assignment_t assignment; + lp_assignment_construct(&assignment, LPVariables::getInstance().lp_var_db); for (const auto& entry : evalMap) { - lp_value_t val; - //Turn into value + lp_value_t val; 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_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(entry.first), &val); lp_value_destruct(&val); } + bool result = false; switch (constraint.relation()) { case Relation::EQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::EQ); + result = lp_polynomial_constraint_evaluate( + poly_pol.get_internal(), LP_SGN_EQ_0, &assignment); + break; case Relation::NEQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::NE); + result = lp_polynomial_constraint_evaluate( + poly_pol.get_internal(), LP_SGN_NE_0, &assignment); + break; case Relation::LESS: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::LT); + result = lp_polynomial_constraint_evaluate( + poly_pol.get_internal(), LP_SGN_LT_0, &assignment); + break; case Relation::LEQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::LE); + result = lp_polynomial_constraint_evaluate( + poly_pol.get_internal(), LP_SGN_LE_0, &assignment); + break; case Relation::GREATER: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::GT); + result = lp_polynomial_constraint_evaluate( + poly_pol.get_internal(), LP_SGN_GT_0, &assignment); + break; case Relation::GEQ: - return evaluate_constraint(poly_pol, assignment, poly::SignCondition::GE); + result = lp_polynomial_constraint_evaluate( + poly_pol.get_internal(), LP_SGN_GE_0, &assignment); + break; default: assert(false); - return false; } + lp_assignment_destruct(&assignment); + return result; } } diff --git a/src/carl-arith/ran/libpoly/RealRoots.cpp b/src/carl-arith/ran/libpoly/RealRoots.cpp index 5d47a8fc..38d9baf1 100644 --- a/src/carl-arith/ran/libpoly/RealRoots.cpp +++ b/src/carl-arith/ran/libpoly/RealRoots.cpp @@ -71,7 +71,8 @@ RealRootsResult real_roots( // Multivariate Polynomial // build the assignment - poly::Assignment assignment(polynomial.context().poly_context()); + lp_assignment_t assignment; + lp_assignment_construct(&assignment, LPVariables::getInstance().lp_var_db); Variable mainVar = polynomial.main_var(); for (Variable& var : carl::variables(polynomial)) { if (var == mainVar) continue; @@ -80,43 +81,57 @@ RealRootsResult real_roots( // 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)); + lp_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(var), &val); // Destroy the value, but dont free the algebraic number! lp_value_destruct(&val); } 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[poly::degree(polynomial.get_polynomial())]; + 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); + lp_value_t val; + 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); + //CARL_LOG_DEBUG("carl.ran.libpoly", " Got eval_val " << eval_val); + + lp_assignment_destruct(&assignment); - if (eval_val == poly::Value(long(0))) { + if (lp_value_cmp(eval_val, poly::Value(long(0)).get_internal()) == 0) { CARL_LOG_DEBUG("carl.ran.libpoly", "poly is 0 after substituting rational assignments -> nullified"); + lp_value_delete(eval_val); return RealRootsResult::nullified_response(); } else { 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()); + lp_assignment_destruct(&assignment); - // turn result into RealRootsResult std::vector res; - for (const poly::Value& val : roots) { - auto tmp = LPRealAlgebraicNumber::create_from_value(val.get_internal()); + for (std::size_t i = 0; i < roots_size; ++i) { + auto tmp = LPRealAlgebraicNumber::create_from_value(&roots[i]); // filter out roots not in interval - if (poly::contains(inter_poly, val)) { - CARL_LOG_DEBUG("carl.ran.libpoly", " Found root " << val); + if (lp_interval_contains(inter_poly.get_internal(), &roots[i])) { + CARL_LOG_DEBUG("carl.ran.libpoly", " Found root " << tmp); res.emplace_back(tmp); } + lp_value_destruct(&roots[i]); } + delete[] roots; + + std::sort(res.begin(), res.end()); + return RealRootsResult::roots_response(std::move(res)); } } From c700612d69ba72592157c61be32e6fcd0a1343b2 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 1 Mar 2024 20:25:46 +0100 Subject: [PATCH 02/29] remove cout --- src/carl-arith/poly/libpoly/LPVariables.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/carl-arith/poly/libpoly/LPVariables.h b/src/carl-arith/poly/libpoly/LPVariables.h index 7f935c08..7d06b2c6 100644 --- a/src/carl-arith/poly/libpoly/LPVariables.h +++ b/src/carl-arith/poly/libpoly/LPVariables.h @@ -26,7 +26,6 @@ class LPVariables : public Singleton { lp_variable_db_t* lp_var_db; LPVariables() { - std::cout << "init vardb" << std::endl; lp_var_db = lp_variable_db_new(); } From 1ed71ec6e7ddc2dc0e5ca4569941867cd543f016 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Tue, 5 Mar 2024 16:49:29 +0100 Subject: [PATCH 03/29] add feature --- src/carl-arith/extended/VariableComparison.h | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/carl-arith/extended/VariableComparison.h b/src/carl-arith/extended/VariableComparison.h index f9540861..c7823dad 100644 --- a/src/carl-arith/extended/VariableComparison.h +++ b/src/carl-arith/extended/VariableComparison.h @@ -153,15 +153,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); From 01843a24a653d50ac30d81da0fdca107c2bec914 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 6 Mar 2024 14:24:17 +0100 Subject: [PATCH 04/29] performance improvements --- src/carl-arith/poly/libpoly/LPContext.h | 15 ++++++++- src/carl-arith/poly/libpoly/LPVariables.cpp | 35 ++++++++++++++++++--- src/carl-arith/poly/libpoly/LPVariables.h | 13 ++++---- src/carl-arith/ran/interval/Evaluation.h | 4 +-- src/carl-arith/ran/libpoly/Evaluation.cpp | 13 ++++---- src/carl-arith/ran/libpoly/RealRoots.cpp | 11 ++++--- src/carl-vs/substitute.h | 6 ++-- 7 files changed, 69 insertions(+), 28 deletions(-) diff --git a/src/carl-arith/poly/libpoly/LPContext.h b/src/carl-arith/poly/libpoly/LPContext.h index b6fde9fe..69b3ba3e 100644 --- a/src/carl-arith/poly/libpoly/LPContext.h +++ b/src/carl-arith/poly/libpoly/LPContext.h @@ -26,7 +26,20 @@ class LPContext { 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_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); } ~Data() { lp_variable_order_detach(lp_var_order); diff --git a/src/carl-arith/poly/libpoly/LPVariables.cpp b/src/carl-arith/poly/libpoly/LPVariables.cpp index 9aabaf37..e2f33d9a 100644 --- a/src/carl-arith/poly/libpoly/LPVariables.cpp +++ b/src/carl-arith/poly/libpoly/LPVariables.cpp @@ -7,32 +7,59 @@ namespace carl { +LPVariables::LPVariables() { + lp_var_db = lp_variable_db_new(); + lp_assignment_construct(&lp_assignment, lp_var_db); + + 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_assignment_destruct(&lp_assignment); + 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) << " -> " << it->second); + 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 << " -> " << lp_variable_db_get_name(lp_var_db, it->second)); + 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) { +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; } - // CARL_LOG_TRACE("carl.poly", "Mapping carl variable " << var << " -> " << lp_variable_db_get_name(lp_var_db, it->second)); +} + +lp_assignment_t& LPVariables::get_assignment() { + if (lp_assignment.values) { + for (size_t i = 0; i < lp_assignment.size; ++ i) { + lp_assignment_set_value(&lp_assignment, i, 0); + } + } + return lp_assignment; } } // namespace carl diff --git a/src/carl-arith/poly/libpoly/LPVariables.h b/src/carl-arith/poly/libpoly/LPVariables.h index 7d06b2c6..e6693ba7 100644 --- a/src/carl-arith/poly/libpoly/LPVariables.h +++ b/src/carl-arith/poly/libpoly/LPVariables.h @@ -22,20 +22,21 @@ class LPVariables : public Singleton { // mapping from libpoly variables to carl variables std::map vars_libpoly_carl; + lp_assignment_t lp_assignment; + public: lp_variable_db_t* lp_var_db; - LPVariables() { - lp_var_db = lp_variable_db_new(); - } + #define TEMP_VARIABLE_SIZE 10 - ~LPVariables() { - lp_variable_db_detach(lp_var_db); - } + 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); + lp_assignment_t& get_assignment(); }; } // namespace carl 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 80a29d74..a42f565d 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -7,9 +7,7 @@ namespace carl { std::optional evaluate( const LPPolynomial& polynomial, const std::map& evalMap) { - - lp_assignment_t assignment; - lp_assignment_construct(&assignment, LPVariables::getInstance().lp_var_db); + lp_assignment_t& assignment = LPVariables::getInstance().get_assignment(); for (const auto& entry : evalMap) { lp_value_t val; lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); @@ -17,7 +15,6 @@ std::optional evaluate( lp_value_destruct(&val); } auto result = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); - lp_assignment_destruct(&assignment); if (result->type == LP_VALUE_NONE) { CARL_LOG_DEBUG("carl.ran.libpoly", " Result: NONE"); @@ -57,11 +54,14 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s } } + for (const auto& v : carl::variables(constraint)) { + if (evalMap.find(v) == evalMap.end()) return boost::indeterminate; + } + //denominator can be omitted const poly::Polynomial& poly_pol = constraint.lhs().get_polynomial(); - lp_assignment_t assignment; - lp_assignment_construct(&assignment, LPVariables::getInstance().lp_var_db); + lp_assignment_t& assignment = LPVariables::getInstance().get_assignment(); for (const auto& entry : evalMap) { lp_value_t val; lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); @@ -98,7 +98,6 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s default: assert(false); } - lp_assignment_destruct(&assignment); return result; } diff --git a/src/carl-arith/ran/libpoly/RealRoots.cpp b/src/carl-arith/ran/libpoly/RealRoots.cpp index 38d9baf1..f15d4d14 100644 --- a/src/carl-arith/ran/libpoly/RealRoots.cpp +++ b/src/carl-arith/ran/libpoly/RealRoots.cpp @@ -67,12 +67,16 @@ RealRootsResult real_roots( return RealRootsResult::no_roots_response(); } + for (const auto& v : carl::variables(polynomial)) { + if (v != polynomial.main_var() && m.find(v) == m.end()) return RealRootsResult::non_univariate_response(); + } + poly::Interval inter_poly = to_libpoly_interval(interval); // Multivariate Polynomial // build the assignment - lp_assignment_t assignment; - lp_assignment_construct(&assignment, LPVariables::getInstance().lp_var_db); + lp_assignment_t& assignment = LPVariables::getInstance().get_assignment(); + Variable mainVar = polynomial.main_var(); for (Variable& var : carl::variables(polynomial)) { if (var == mainVar) continue; @@ -102,7 +106,6 @@ RealRootsResult real_roots( auto eval_val = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); //CARL_LOG_DEBUG("carl.ran.libpoly", " Got eval_val " << eval_val); - lp_assignment_destruct(&assignment); if (lp_value_cmp(eval_val, poly::Value(long(0)).get_internal()) == 0) { CARL_LOG_DEBUG("carl.ran.libpoly", "poly is 0 after substituting rational assignments -> nullified"); @@ -115,8 +118,6 @@ RealRootsResult real_roots( } } - lp_assignment_destruct(&assignment); - std::vector res; for (std::size_t i = 0; i < roots_size; ++i) { auto tmp = LPRealAlgebraicNumber::create_from_value(&roots[i]); 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; From d51ec8111e305cd9613f59d1226784bdd23bdb59 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Thu, 7 Mar 2024 14:22:12 +0100 Subject: [PATCH 05/29] libpoly integration profiling --- src/carl-arith/poly/libpoly/LPContext.cpp | 33 ++++++++++++++++ src/carl-arith/poly/libpoly/LPContext.h | 23 +---------- src/carl-arith/poly/libpoly/LPVariables.cpp | 12 ------ src/carl-arith/poly/libpoly/LPVariables.h | 3 -- src/carl-arith/ran/libpoly/Evaluation.cpp | 19 +++------ src/carl-arith/ran/libpoly/LPAssignment.cpp | 44 +++++++++++++++++++++ src/carl-arith/ran/libpoly/LPAssignment.h | 33 ++++++++++++++++ src/carl-arith/ran/libpoly/LPRan.cpp | 44 ++++++++++++--------- src/carl-arith/ran/libpoly/LPRan.h | 14 ++++--- src/carl-arith/ran/libpoly/RealRoots.cpp | 18 +++------ 10 files changed, 157 insertions(+), 86 deletions(-) create mode 100644 src/carl-arith/poly/libpoly/LPContext.cpp create mode 100644 src/carl-arith/ran/libpoly/LPAssignment.cpp create mode 100644 src/carl-arith/ran/libpoly/LPAssignment.h 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 69b3ba3e..4d1a57f5 100644 --- a/src/carl-arith/poly/libpoly/LPContext.h +++ b/src/carl-arith/poly/libpoly/LPContext.h @@ -24,27 +24,8 @@ class LPContext { lp_variable_order_t* lp_var_order; lp_polynomial_context_t* lp_context; - 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); - } - ~Data() { - lp_variable_order_detach(lp_var_order); - lp_polynomial_context_detach(lp_context); - } + Data(const std::vector& v); + ~Data(); }; std::shared_ptr m_data; diff --git a/src/carl-arith/poly/libpoly/LPVariables.cpp b/src/carl-arith/poly/libpoly/LPVariables.cpp index e2f33d9a..cfd8c625 100644 --- a/src/carl-arith/poly/libpoly/LPVariables.cpp +++ b/src/carl-arith/poly/libpoly/LPVariables.cpp @@ -9,8 +9,6 @@ namespace carl { LPVariables::LPVariables() { lp_var_db = lp_variable_db_new(); - lp_assignment_construct(&lp_assignment, lp_var_db); - for (size_t i = 0; i < TEMP_VARIABLE_SIZE; ++ i) { char name[10]; sprintf(name, "#%zu", i); @@ -19,7 +17,6 @@ LPVariables::LPVariables() { } LPVariables::~LPVariables() { - lp_assignment_destruct(&lp_assignment); lp_variable_db_detach(lp_var_db); } @@ -52,15 +49,6 @@ lp_variable_t LPVariables::lp_variable(carl::Variable var) { return poly_var; } } - -lp_assignment_t& LPVariables::get_assignment() { - if (lp_assignment.values) { - for (size_t i = 0; i < lp_assignment.size; ++ i) { - lp_assignment_set_value(&lp_assignment, i, 0); - } - } - return lp_assignment; -} } // namespace carl diff --git a/src/carl-arith/poly/libpoly/LPVariables.h b/src/carl-arith/poly/libpoly/LPVariables.h index e6693ba7..b169d088 100644 --- a/src/carl-arith/poly/libpoly/LPVariables.h +++ b/src/carl-arith/poly/libpoly/LPVariables.h @@ -22,8 +22,6 @@ class LPVariables : public Singleton { // mapping from libpoly variables to carl variables std::map vars_libpoly_carl; - lp_assignment_t lp_assignment; - public: lp_variable_db_t* lp_var_db; @@ -36,7 +34,6 @@ class LPVariables : public Singleton { 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); - lp_assignment_t& get_assignment(); }; } // namespace carl diff --git a/src/carl-arith/ran/libpoly/Evaluation.cpp b/src/carl-arith/ran/libpoly/Evaluation.cpp index a42f565d..98b84ead 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -2,18 +2,15 @@ #ifdef USE_LIBPOLY +#include "LPAssignment.h" + + namespace carl { std::optional evaluate( const LPPolynomial& polynomial, const std::map& evalMap) { - lp_assignment_t& assignment = LPVariables::getInstance().get_assignment(); - for (const auto& entry : evalMap) { - lp_value_t val; - lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); - lp_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(entry.first), &val); - lp_value_destruct(&val); - } + lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); auto result = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); if (result->type == LP_VALUE_NONE) { @@ -61,13 +58,7 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s //denominator can be omitted const poly::Polynomial& poly_pol = constraint.lhs().get_polynomial(); - lp_assignment_t& assignment = LPVariables::getInstance().get_assignment(); - for (const auto& entry : evalMap) { - lp_value_t val; - lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); - lp_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(entry.first), &val); - lp_value_destruct(&val); - } + lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); bool result = false; switch (constraint.relation()) { diff --git a/src/carl-arith/ran/libpoly/LPAssignment.cpp b/src/carl-arith/ran/libpoly/LPAssignment.cpp new file mode 100644 index 00000000..252cf0d4 --- /dev/null +++ b/src/carl-arith/ran/libpoly/LPAssignment.cpp @@ -0,0 +1,44 @@ +#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_value_t val; + lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); + lp_assignment_set_value(&lp_assignment, LPVariables::getInstance().lp_variable(entry.first), &val); + lp_value_destruct(&val); + } + 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..7a9d9bf9 --- /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..a256fa84 100644 --- a/src/carl-arith/ran/libpoly/LPRan.cpp +++ b/src/carl-arith/ran/libpoly/LPRan.cpp @@ -7,32 +7,29 @@ namespace carl { using NumberType = LPRealAlgebraicNumber::NumberType; LPRealAlgebraicNumber::~LPRealAlgebraicNumber() { - lp_algebraic_number_destruct(get_internal()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber() { +LPRealAlgebraicNumber::LPRealAlgebraicNumber() : m_content(std::make_shared()) { lp_algebraic_number_construct_zero(get_internal()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const poly::AlgebraicNumber& num) { +LPRealAlgebraicNumber::LPRealAlgebraicNumber(const poly::AlgebraicNumber& num) : LPRealAlgebraicNumber() { lp_algebraic_number_construct_copy(get_internal(), num.get_internal()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const lp_algebraic_number_t& num) { +LPRealAlgebraicNumber::LPRealAlgebraicNumber(const lp_algebraic_number_t& num) : LPRealAlgebraicNumber() { lp_algebraic_number_construct_copy(get_internal(), &num); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(poly::AlgebraicNumber&& num) - : LPRealAlgebraicNumber() { +LPRealAlgebraicNumber::LPRealAlgebraicNumber(poly::AlgebraicNumber&& num) : LPRealAlgebraicNumber() { lp_algebraic_number_swap(get_internal(), num.get_internal()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(lp_algebraic_number_t&& num) - : LPRealAlgebraicNumber() { +LPRealAlgebraicNumber::LPRealAlgebraicNumber(lp_algebraic_number_t&& num) : LPRealAlgebraicNumber() { lp_algebraic_number_swap(get_internal(), &num); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial& p, const Interval& i) { +LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial& p, const Interval& i) : LPRealAlgebraicNumber() { CARL_LOG_DEBUG("carl.ran.libpoly", " Create safe from poly: " << p << " in interval: " << i); poly::UPolynomial upoly = to_libpoly_upolynomial(p); @@ -61,25 +58,21 @@ LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial m_content; static const Variable auxVariable; - - public: ~LPRealAlgebraicNumber(); @@ -100,14 +104,14 @@ class LPRealAlgebraicNumber { * @return Pointer to the internal libpoly algebraic number (C interface) */ inline lp_algebraic_number_t* get_internal() { - return &m_content; + 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; + return &(m_content->c); } /** diff --git a/src/carl-arith/ran/libpoly/RealRoots.cpp b/src/carl-arith/ran/libpoly/RealRoots.cpp index f15d4d14..ffccf5f7 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 { @@ -75,20 +77,10 @@ RealRootsResult real_roots( // Multivariate Polynomial // build the assignment - lp_assignment_t& assignment = LPVariables::getInstance().get_assignment(); - 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 - lp_assignment_set_value(&assignment, LPVariables::getInstance().lp_variable(var), &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"); lp_value_t* roots = new lp_value_t[poly::degree(polynomial.get_polynomial())]; From 707ae0055de4dfb0184c669df3918860da01f176 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Thu, 7 Mar 2024 15:47:18 +0100 Subject: [PATCH 06/29] context --- cmake/patches/libpoly_variable_db.patch | 13 +++++++++++++ src/carl-arith/extended/VariableComparison.h | 3 +++ src/carl-arith/poly/libpoly/LPContext.h | 10 ++++++++++ src/carl-arith/poly/libpoly/LPPolynomial.cpp | 13 +++++++++++++ src/carl-arith/poly/libpoly/LPPolynomial.h | 2 ++ 5 files changed, 41 insertions(+) diff --git a/cmake/patches/libpoly_variable_db.patch b/cmake/patches/libpoly_variable_db.patch index 09344ebc..93ebc768 100644 --- a/cmake/patches/libpoly_variable_db.patch +++ b/cmake/patches/libpoly_variable_db.patch @@ -1,3 +1,16 @@ +diff --git a/include/polynomial.h b/include/polynomial.h +index 8a9533f..12123e8 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); + diff --git a/src/variable/variable_db.c b/src/variable/variable_db.c index 60f3df4..33866c4 100644 --- a/src/variable/variable_db.c diff --git a/src/carl-arith/extended/VariableComparison.h b/src/carl-arith/extended/VariableComparison.h index c7823dad..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; } diff --git a/src/carl-arith/poly/libpoly/LPContext.h b/src/carl-arith/poly/libpoly/LPContext.h index 4d1a57f5..065629d5 100644 --- a/src/carl-arith/poly/libpoly/LPContext.h +++ b/src/carl-arith/poly/libpoly/LPContext.h @@ -91,6 +91,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 1991067c..1781a7fa 100644 --- a/src/carl-arith/poly/libpoly/LPPolynomial.cpp +++ b/src/carl-arith/poly/libpoly/LPPolynomial.cpp @@ -549,6 +549,19 @@ std::ostream& operator<<(std::ostream& os, const LPPolynomial& p) { 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..ac3c6107 100644 --- a/src/carl-arith/poly/libpoly/LPPolynomial.h +++ b/src/carl-arith/poly/libpoly/LPPolynomial.h @@ -288,6 +288,8 @@ class LPPolynomial { return m_context; } + void set_context(const LPContext& c); + /** * Checks if the given variable occurs in the polynomial. * @param v Variable. From 2d7a69b38abc36f5d2ee8e7665bb24f3151ffc6f Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Tue, 12 Mar 2024 09:43:07 +0100 Subject: [PATCH 07/29] use lp_vale_t in LPRan --- src/carl-arith/ran/libpoly/Evaluation.cpp | 6 +- src/carl-arith/ran/libpoly/LPAssignment.cpp | 5 +- src/carl-arith/ran/libpoly/LPRan.cpp | 271 +++++++------------- src/carl-arith/ran/libpoly/LPRan.h | 42 +-- src/carl-arith/ran/libpoly/RealRoots.cpp | 11 +- src/carl-arith/ran/libpoly/helper.h | 35 ++- 6 files changed, 134 insertions(+), 236 deletions(-) diff --git a/src/carl-arith/ran/libpoly/Evaluation.cpp b/src/carl-arith/ran/libpoly/Evaluation.cpp index 98b84ead..91545894 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -18,10 +18,8 @@ std::optional evaluate( return std::nullopt; } - auto ran = LPRealAlgebraicNumber::create_from_value(result); + auto ran = LPRealAlgebraicNumber(std::move(*result)); CARL_LOG_DEBUG("carl.ran.libpoly", " Result: " << ran); - lp_value_destruct(result); - free(result); return ran; } @@ -60,6 +58,8 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); + // std::cout << lp_assignment_to_string(&assignment) << std::endl; + bool result = false; switch (constraint.relation()) { case Relation::EQ: diff --git a/src/carl-arith/ran/libpoly/LPAssignment.cpp b/src/carl-arith/ran/libpoly/LPAssignment.cpp index 252cf0d4..f7cb13c1 100644 --- a/src/carl-arith/ran/libpoly/LPAssignment.cpp +++ b/src/carl-arith/ran/libpoly/LPAssignment.cpp @@ -28,10 +28,7 @@ lp_assignment_t& LPAssignment::get(const carl::Assignment } } for (const auto& entry : ass) { - lp_value_t val; - lp_value_construct(&val, lp_value_type_t::LP_VALUE_ALGEBRAIC, entry.second.get_internal()); - lp_assignment_set_value(&lp_assignment, LPVariables::getInstance().lp_variable(entry.first), &val); - lp_value_destruct(&val); + lp_assignment_set_value(&lp_assignment, LPVariables::getInstance().lp_variable(entry.first), entry.second.get_internal()); } return lp_assignment; } diff --git a/src/carl-arith/ran/libpoly/LPRan.cpp b/src/carl-arith/ran/libpoly/LPRan.cpp index a256fa84..51f6db48 100644 --- a/src/carl-arith/ran/libpoly/LPRan.cpp +++ b/src/carl-arith/ran/libpoly/LPRan.cpp @@ -10,23 +10,15 @@ LPRealAlgebraicNumber::~LPRealAlgebraicNumber() { } LPRealAlgebraicNumber::LPRealAlgebraicNumber() : m_content(std::make_shared()) { - lp_algebraic_number_construct_zero(get_internal()); + lp_value_construct_zero(get_internal()); } -LPRealAlgebraicNumber::LPRealAlgebraicNumber(const poly::AlgebraicNumber& num) : LPRealAlgebraicNumber() { - 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) : LPRealAlgebraicNumber() { - lp_algebraic_number_construct_copy(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(lp_value_t&& num) : m_content(std::make_shared()) { + *get_internal() = num; } LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial& p, const Interval& i) : LPRealAlgebraicNumber() { @@ -37,30 +29,24 @@ LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial 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); } } - 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_value_construct_copy(get_internal(), poly::Value(res.front()).get_internal()); } /** * Construct from NumberType (usually mpq_class) */ LPRealAlgebraicNumber::LPRealAlgebraicNumber(const NumberType& num) : LPRealAlgebraicNumber() { - poly::Rational rat = to_libpoly_rational(num); - lp_algebraic_number_construct_from_rational(get_internal(), rat.get_internal()); + lp_value_construct(get_internal(), LP_VALUE_RATIONAL, (lp_rational_t*)num.get_mpq_t()); } LPRealAlgebraicNumber::LPRealAlgebraicNumber(const LPRealAlgebraicNumber& ran) : m_content(ran.m_content) { @@ -83,63 +69,33 @@ 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()); + return lp_value_is_rational(get_internal()); } const poly::UPolynomial LPRealAlgebraicNumber::libpoly_polynomial() const { - return poly::UPolynomial(static_cast(get_internal()->f)); + assert(!is_numeric()); + return poly::UPolynomial(static_cast(get_internal()->value.a.f)); } const poly::DyadicInterval& LPRealAlgebraicNumber::libpoly_interval() const { - return *reinterpret_cast(&get_internal()->I); + assert(!is_numeric()); + return *reinterpret_cast(&get_internal()->value.a.I); } const UnivariatePolynomial LPRealAlgebraicNumber::polynomial() const { @@ -155,11 +111,13 @@ const Interval LPRealAlgebraicNumber::interval() const { } const NumberType LPRealAlgebraicNumber::get_upper_bound() const { - return to_rational(poly::get_upper(libpoly_interval())); + if (is_numeric()) return value(); + else return to_rational(poly::get_upper(libpoly_interval())); } const NumberType LPRealAlgebraicNumber::get_lower_bound() const { - return to_rational(poly::get_lower(libpoly_interval())); + if (is_numeric()) return value(); + else return to_rational(poly::get_lower(libpoly_interval())); } /** @@ -168,32 +126,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; } /** @@ -201,29 +139,26 @@ 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_zero(&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 { @@ -232,7 +167,7 @@ std::size_t bitsize(const LPRealAlgebraicNumber& n) { } 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&) { @@ -242,17 +177,8 @@ 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; + return lp_interval_contains(inter.get_internal(), n.get_internal()); } @@ -264,8 +190,10 @@ bool contained_in(const LPRealAlgebraicNumber& n, const Interval 0 && !n.is_numeric()) { - lp_algebraic_number_refine_const(n.get_internal()); + lp_algebraic_number_refine_const(&n.get_internal()->value.a); } CARL_LOG_DEBUG("carl.ran.libpoly", "Finished Refining Algebraic NumberType : " << n); } @@ -275,103 +203,90 @@ void LPRealAlgebraicNumber::refine() const { } NumberType branching_point(const LPRealAlgebraicNumber& n) { - //return carl::sample(n.interval_int()); - refine(n); + if (n.is_numeric()) return n.value(); + + refine(n); poly::DyadicRational res; - lp_algebraic_number_get_dyadic_midpoint(n.get_internal(), res.get_internal()); - NumberType num = to_rational(res); - return num; + lp_algebraic_number_get_dyadic_midpoint(&n.get_internal()->value.a, res.get_internal()); + return to_rational(res); } 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_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_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_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_get_value_between(lower.get_internal(), true, &v_upper, true, &res); + + lp_value_destruct(&v_upper); -//Make sure that the intervals are disjoint + 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_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); + lp_value_floor(n.get_internal(), &val); NumberType ret = to_rational(*poly::detail::cast_from(&val)); lp_integer_destruct(&val) ; return ret ; @@ -380,15 +295,13 @@ 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); + lp_value_ceiling(n.get_internal(), &val); NumberType ret = to_rational(*poly::detail::cast_from(&val)); lp_integer_destruct(&val) ; return ret ; } - bool compare(const LPRealAlgebraicNumber& lhs, const LPRealAlgebraicNumber& rhs, const Relation relation) { if (lhs.m_content == rhs.m_content) { switch (relation) { @@ -404,7 +317,7 @@ bool compare(const LPRealAlgebraicNumber& lhs, const LPRealAlgebraicNumber& rhs, } } - int cmp = lp_algebraic_number_cmp(lhs.get_internal(), rhs.get_internal()); + 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; @@ -430,9 +343,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: @@ -455,7 +366,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 8158a63c..9776daa6 100644 --- a/src/carl-arith/ran/libpoly/LPRan.h +++ b/src/carl-arith/ran/libpoly/LPRan.h @@ -26,9 +26,9 @@ class LPRealAlgebraicNumber { private: struct Content { - lp_algebraic_number_t c; + lp_value_t c; ~Content() { - lp_algebraic_number_destruct(&c); + lp_value_destruct(&c); } }; mutable std::shared_ptr m_content; @@ -47,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 @@ -91,26 +80,17 @@ 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() { + 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 { + inline const lp_value_t* get_internal() const { return &(m_content->c); } @@ -119,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; @@ -150,6 +127,9 @@ class LPRealAlgebraicNumber { friend NumberType ceil(const LPRealAlgebraicNumber& n); void refine() const; + + const poly::UPolynomial libpoly_polynomial() const; + const poly::DyadicInterval& libpoly_interval() const; }; bool compare(const LPRealAlgebraicNumber&, const LPRealAlgebraicNumber&, const Relation); @@ -168,7 +148,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); @@ -191,7 +171,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 ffccf5f7..19732100 100644 --- a/src/carl-arith/ran/libpoly/RealRoots.cpp +++ b/src/carl-arith/ran/libpoly/RealRoots.cpp @@ -39,11 +39,10 @@ RealRootsResult real_roots( // turn into RealRootsResult 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); + res.emplace_back(*poly::Value(val).get_internal()); } } @@ -98,7 +97,6 @@ RealRootsResult real_roots( auto eval_val = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); //CARL_LOG_DEBUG("carl.ran.libpoly", " Got eval_val " << eval_val); - if (lp_value_cmp(eval_val, poly::Value(long(0)).get_internal()) == 0) { CARL_LOG_DEBUG("carl.ran.libpoly", "poly is 0 after substituting rational assignments -> nullified"); lp_value_delete(eval_val); @@ -112,13 +110,10 @@ RealRootsResult real_roots( std::vector res; for (std::size_t i = 0; i < roots_size; ++i) { - auto tmp = LPRealAlgebraicNumber::create_from_value(&roots[i]); - // filter out roots not in interval if (lp_interval_contains(inter_poly.get_internal(), &roots[i])) { - CARL_LOG_DEBUG("carl.ran.libpoly", " Found root " << tmp); - res.emplace_back(tmp); + res.emplace_back(std::move(roots[i])); + CARL_LOG_DEBUG("carl.ran.libpoly", " Found root " << res.back()); } - lp_value_destruct(&roots[i]); } delete[] roots; diff --git a/src/carl-arith/ran/libpoly/helper.h b/src/carl-arith/ran/libpoly/helper.h index 33439120..775e3116 100644 --- a/src/carl-arith/ran/libpoly/helper.h +++ b/src/carl-arith/ran/libpoly/helper.h @@ -67,23 +67,42 @@ inline poly::UPolynomial to_libpoly_upolynomial(const carl::UnivariatePolynomial 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(*reinterpret_cast(m)); } inline mpq_class to_rational(const lp_rational_t* m) { return *reinterpret_cast(m); } +inline mpq_class to_rational(const lp_value_t* m){ + switch(m->type){ + case lp_value_type_t::LP_VALUE_INTEGER: + return to_rational(&m->value.z); + case lp_value_type_t::LP_VALUE_RATIONAL: + return to_rational(&m->value.q); + // case lp_value_type_t::LP_VALUE_DYADIC_RATIONAL: + // return to_rational(m->value.q); + default: + CARL_LOG_ERROR("carl.converter", "Cannot convert libpoly value: " << m << " to rational."); + assert(false); + return mpq_class(0); + } +} + inline mpz_class to_integer(const poly::Integer& m) { return *poly::detail::cast_to_gmp(&m); } +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 poly::DyadicRational& m) { return mpq_class(to_integer(poly::numerator(m)), to_integer(poly::denominator(m))); } @@ -103,10 +122,6 @@ inline mpq_class to_rational(const poly::Value& m){ } } -inline poly::Rational to_libpoly_rational(const mpq_class& num) { - return poly::Rational(num); -} - //Exact for whole numbers inline poly::DyadicRational get_libpoly_dyadic_approximation(const mpz_class& num) { return poly::DyadicRational(poly::Integer(num)); From a3e9441b6ef152a7894217b3fe1770ee2c4f9606 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Tue, 12 Mar 2024 15:10:37 +0100 Subject: [PATCH 08/29] fixes --- src/carl-arith/poly/libpoly/helper.h | 12 +-- src/carl-arith/ran/libpoly/LPRan.cpp | 48 ++++++------ src/carl-arith/ran/libpoly/LPRan.h | 3 - src/carl-arith/ran/libpoly/helper.h | 112 ++++++--------------------- 4 files changed, 52 insertions(+), 123 deletions(-) diff --git a/src/carl-arith/poly/libpoly/helper.h b/src/carl-arith/poly/libpoly/helper.h index bf7a27d8..2dd831dc 100644 --- a/src/carl-arith/poly/libpoly/helper.h +++ b/src/carl-arith/poly/libpoly/helper.h @@ -7,15 +7,15 @@ 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; + lp_polynomial_t* p = lp_polynomial_alloc(); + lp_polynomial_construct_simple(p, c, poly::Integer(1).get_internal(), v, 1); + return poly::Polynomial(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; + lp_polynomial_t* p = lp_polynomial_alloc(); + lp_polynomial_construct_simple(p, c, i.get_internal(), lp_variable_null, 0) ; + return poly::Polynomial(p); } } diff --git a/src/carl-arith/ran/libpoly/LPRan.cpp b/src/carl-arith/ran/libpoly/LPRan.cpp index 51f6db48..f6f28d98 100644 --- a/src/carl-arith/ran/libpoly/LPRan.cpp +++ b/src/carl-arith/ran/libpoly/LPRan.cpp @@ -88,18 +88,9 @@ bool LPRealAlgebraicNumber::is_numeric() const { return lp_value_is_rational(get_internal()); } -const poly::UPolynomial LPRealAlgebraicNumber::libpoly_polynomial() const { - assert(!is_numeric()); - return poly::UPolynomial(static_cast(get_internal()->value.a.f)); -} - -const poly::DyadicInterval& LPRealAlgebraicNumber::libpoly_interval() const { - assert(!is_numeric()); - return *reinterpret_cast(&get_internal()->value.a.I); -} - const UnivariatePolynomial LPRealAlgebraicNumber::polynomial() const { - return to_carl_univariate_polynomial(libpoly_polynomial(), auxVariable); + assert(!is_numeric()); + return to_carl_univariate_polynomial(poly::UPolynomial(static_cast(get_internal()->value.a.f)), auxVariable); } const Interval LPRealAlgebraicNumber::interval() const { @@ -110,14 +101,14 @@ const Interval LPRealAlgebraicNumber::interval() const { } } -const NumberType LPRealAlgebraicNumber::get_upper_bound() const { +const NumberType LPRealAlgebraicNumber::get_lower_bound() const { if (is_numeric()) return value(); - else return to_rational(poly::get_upper(libpoly_interval())); + else return to_rational(&get_internal()->value.a.I.a); } -const NumberType LPRealAlgebraicNumber::get_lower_bound() const { +const NumberType LPRealAlgebraicNumber::get_upper_bound() const { if (is_numeric()) return value(); - else return to_rational(poly::get_lower(libpoly_interval())); + else return to_rational(&get_internal()->value.a.I.b); } /** @@ -151,7 +142,7 @@ LPRealAlgebraicNumber abs(const LPRealAlgebraicNumber& n) { return LPRealAlgebraicNumber(n); } else { lp_value_t val; - lp_value_construct_zero(&val); + lp_value_construct_none(&val); lp_value_neg(&val, n.get_internal()); auto ret = LPRealAlgebraicNumber(std::move(val)); return ret; @@ -162,7 +153,7 @@ std::size_t bitsize(const LPRealAlgebraicNumber& n) { 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); } } @@ -192,7 +183,7 @@ void refine(const LPRealAlgebraicNumber& n) { if (n.is_numeric()) return; - while (poly::log_size(n.libpoly_interval()) > 0 && !n.is_numeric()) { + 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); @@ -204,11 +195,13 @@ void LPRealAlgebraicNumber::refine() const { NumberType branching_point(const LPRealAlgebraicNumber& n) { if (n.is_numeric()) return n.value(); - refine(n); - poly::DyadicRational res; - lp_algebraic_number_get_dyadic_midpoint(&n.get_internal()->value.a, res.get_internal()); - return to_rational(res); + 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; } @@ -219,6 +212,7 @@ NumberType sample_above(const LPRealAlgebraicNumber& n) { lp_value_construct(&inf, LP_VALUE_PLUS_INFINITY, 0); 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); @@ -234,6 +228,7 @@ NumberType sample_below(const LPRealAlgebraicNumber& n) { lp_value_construct(&inf, LP_VALUE_MINUS_INFINITY, 0); 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); @@ -245,6 +240,7 @@ NumberType sample_between(const LPRealAlgebraicNumber& lower, const LPRealAlgebr CARL_LOG_DEBUG("carl.ran.libpoly", "Sampling between: " << lower << " and " << upper); 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); @@ -258,6 +254,7 @@ NumberType sample_between(const LPRealAlgebraicNumber& lower, const NumberType& 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); lp_value_destruct(&v_upper); @@ -274,6 +271,7 @@ NumberType sample_between(const NumberType& lower, const LPRealAlgebraicNumber& 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); lp_value_destruct(&v_lower); @@ -286,8 +284,9 @@ NumberType sample_between(const NumberType& lower, const LPRealAlgebraicNumber& NumberType floor(const LPRealAlgebraicNumber& n) { CARL_LOG_DEBUG("carl.ran.libpoly", "Floor of: " << n); lp_integer_t val; + lp_integer_construct(&val); lp_value_floor(n.get_internal(), &val); - NumberType ret = to_rational(*poly::detail::cast_from(&val)); + NumberType ret = to_rational(&val); lp_integer_destruct(&val) ; return ret ; } @@ -296,8 +295,9 @@ NumberType floor(const LPRealAlgebraicNumber& n) { NumberType ceil(const LPRealAlgebraicNumber& n) { CARL_LOG_DEBUG("carl.ran.libpoly", "Ceil of: " << n); lp_integer_t val; + lp_integer_construct(&val); lp_value_ceiling(n.get_internal(), &val); - NumberType ret = to_rational(*poly::detail::cast_from(&val)); + NumberType ret = to_rational(&val); lp_integer_destruct(&val) ; return ret ; } diff --git a/src/carl-arith/ran/libpoly/LPRan.h b/src/carl-arith/ran/libpoly/LPRan.h index 9776daa6..5ba37c93 100644 --- a/src/carl-arith/ran/libpoly/LPRan.h +++ b/src/carl-arith/ran/libpoly/LPRan.h @@ -127,9 +127,6 @@ class LPRealAlgebraicNumber { friend NumberType ceil(const LPRealAlgebraicNumber& n); void refine() const; - - const poly::UPolynomial libpoly_polynomial() const; - const poly::DyadicInterval& libpoly_interval() const; }; bool compare(const LPRealAlgebraicNumber&, const LPRealAlgebraicNumber&, const Relation); diff --git a/src/carl-arith/ran/libpoly/helper.h b/src/carl-arith/ran/libpoly/helper.h index 775e3116..5e394d96 100644 --- a/src/carl-arith/ran/libpoly/helper.h +++ b/src/carl-arith/ran/libpoly/helper.h @@ -67,54 +67,43 @@ inline poly::UPolynomial to_libpoly_upolynomial(const carl::UnivariatePolynomial return poly::UPolynomial(coefficients); } +inline mpz_class to_integer(const lp_integer_t* m) { + return *reinterpret_cast(m); +} + inline mpq_class to_rational(const lp_integer_t* m) { - return mpq_class(*reinterpret_cast(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_value_t* m){ - switch(m->type){ - case lp_value_type_t::LP_VALUE_INTEGER: - return to_rational(&m->value.z); - case lp_value_type_t::LP_VALUE_RATIONAL: - return to_rational(&m->value.q); - // case lp_value_type_t::LP_VALUE_DYADIC_RATIONAL: - // return to_rational(m->value.q); - default: - CARL_LOG_ERROR("carl.converter", "Cannot convert libpoly value: " << m << " to rational."); - assert(false); - return mpq_class(0); - } -} - +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); -inline mpq_class to_rational(const poly::Integer& m) { - return *poly::detail::cast_to_gmp(&m); -} + auto res = mpq_class(to_integer(&num), to_integer(&den)); -inline mpq_class to_rational(const poly::Rational& m) { - return *poly::detail::cast_to_gmp(&m); -} + lp_integer_destruct(&num); + lp_integer_destruct(&den); -inline mpq_class to_rational(const poly::DyadicRational& m) { - return mpq_class(to_integer(poly::numerator(m)), to_integer(poly::denominator(m))); + 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); @@ -122,39 +111,6 @@ inline mpq_class to_rational(const poly::Value& m){ } } -//Exact for whole numbers -inline poly::DyadicRational get_libpoly_dyadic_approximation(const mpz_class& num) { - return poly::DyadicRational(poly::Integer(num)); -} - -/** - * 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)); - - 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()); - } - - 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())); - - return poly::DyadicRational(c.get_si(), d.get_ui()); -} - /** * Convert a carl interval to a libpoly interval * This keeps the bound types @@ -195,30 +151,6 @@ inline poly::Interval to_libpoly_interval(const carl::Interval& inter 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); -} - } // namespace carl #endif \ No newline at end of file From c52ac1a103fa8509d95380cacc53b378b5e21b09 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Tue, 12 Mar 2024 16:05:25 +0100 Subject: [PATCH 09/29] more efficient conversions --- src/carl-arith/poly/Conversion.h | 31 ++---- src/carl-arith/ran/libpoly/Evaluation.cpp | 20 ++-- src/carl-arith/ran/libpoly/LPRan.cpp | 22 ++-- src/carl-arith/ran/libpoly/RealRoots.cpp | 25 +++-- src/carl-arith/ran/libpoly/helper.h | 127 ++++++++++++---------- src/tests/carl-arith-ran/Test_Libpoly.cpp | 16 +-- 6 files changed, 125 insertions(+), 116 deletions(-) diff --git a/src/carl-arith/poly/Conversion.h b/src/carl-arith/poly/Conversion.h index f8f5d8e7..47146eec 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(context, poly::Polynomial(res)); } }; diff --git a/src/carl-arith/ran/libpoly/Evaluation.cpp b/src/carl-arith/ran/libpoly/Evaluation.cpp index 91545894..fdf68bfd 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -54,7 +54,7 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s } //denominator can be omitted - const poly::Polynomial& poly_pol = constraint.lhs().get_polynomial(); + auto poly_pol = constraint.lhs().get_polynomial().get_internal(); lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); @@ -63,28 +63,22 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s bool result = false; switch (constraint.relation()) { case Relation::EQ: - result = lp_polynomial_constraint_evaluate( - poly_pol.get_internal(), LP_SGN_EQ_0, &assignment); + result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_EQ_0, &assignment); break; case Relation::NEQ: - result = lp_polynomial_constraint_evaluate( - poly_pol.get_internal(), LP_SGN_NE_0, &assignment); + result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_NE_0, &assignment); break; case Relation::LESS: - result = lp_polynomial_constraint_evaluate( - poly_pol.get_internal(), LP_SGN_LT_0, &assignment); + result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_LT_0, &assignment); break; case Relation::LEQ: - result = lp_polynomial_constraint_evaluate( - poly_pol.get_internal(), LP_SGN_LE_0, &assignment); + result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_LE_0, &assignment); break; case Relation::GREATER: - result = lp_polynomial_constraint_evaluate( - poly_pol.get_internal(), LP_SGN_GT_0, &assignment); + result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_GT_0, &assignment); break; case Relation::GEQ: - result = lp_polynomial_constraint_evaluate( - poly_pol.get_internal(), LP_SGN_GE_0, &assignment); + result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_GE_0, &assignment); break; default: assert(false); diff --git a/src/carl-arith/ran/libpoly/LPRan.cpp b/src/carl-arith/ran/libpoly/LPRan.cpp index f6f28d98..91f650c2 100644 --- a/src/carl-arith/ran/libpoly/LPRan.cpp +++ b/src/carl-arith/ran/libpoly/LPRan.cpp @@ -25,21 +25,26 @@ LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial roots = poly::isolate_real_roots(upoly); std::vector res; for (const auto& val : roots) { - if (poly::contains(inter_poly, poly::Value(val))) { + lp_value_t tmp; + lp_value_construct(&tmp, LP_VALUE_ALGEBRAIC, val.get_internal()); + if (lp_interval_contains(inter_poly, &tmp)) { CARL_LOG_DEBUG("carl.ran.libpoly", " Found Root " << val); res.emplace_back(val); } + lp_value_destruct(&tmp); } assert(res.size() == 1); - lp_value_construct_copy(get_internal(), poly::Value(res.front()).get_internal()); + lp_interval_destruct(inter_poly); + delete inter_poly; + lp_value_construct(get_internal(), LP_VALUE_ALGEBRAIC, res.front().get_internal()); } /** @@ -90,7 +95,7 @@ bool LPRealAlgebraicNumber::is_numeric() const { const UnivariatePolynomial LPRealAlgebraicNumber::polynomial() const { assert(!is_numeric()); - return to_carl_univariate_polynomial(poly::UPolynomial(static_cast(get_internal()->value.a.f)), auxVariable); + return to_carl_univariate_polynomial(get_internal()->value.a.f, auxVariable); } const Interval LPRealAlgebraicNumber::interval() const { @@ -168,8 +173,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); - return lp_interval_contains(inter.get_internal(), n.get_internal()); + 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; } diff --git a/src/carl-arith/ran/libpoly/RealRoots.cpp b/src/carl-arith/ran/libpoly/RealRoots.cpp index 19732100..f787f2fc 100644 --- a/src/carl-arith/ran/libpoly/RealRoots.cpp +++ b/src/carl-arith/ran/libpoly/RealRoots.cpp @@ -12,7 +12,7 @@ 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)) { @@ -23,8 +23,6 @@ RealRootsResult real_roots( return RealRootsResult::no_roots_response(); } - poly::Interval inter_poly = to_libpoly_interval(interval); - // actual calculations std::vector roots = poly::isolate_real_roots(poly::to_univariate(polynomial.get_polynomial())); @@ -37,15 +35,20 @@ RealRootsResult real_roots( 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) { - // 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(*poly::Value(val).get_internal()); + lp_value_t tmp; + lp_value_construct(&tmp, LP_VALUE_ALGEBRAIC, val.get_internal()); + if (lp_interval_contains(inter_poly, &tmp)) { + res.emplace_back(std::move(tmp)); } } + lp_interval_destruct(inter_poly); + delete inter_poly; + return RealRootsResult::roots_response(std::move(res)); } @@ -72,8 +75,6 @@ RealRootsResult real_roots( if (v != polynomial.main_var() && m.find(v) == m.end()) return RealRootsResult::non_univariate_response(); } - poly::Interval inter_poly = to_libpoly_interval(interval); - // Multivariate Polynomial // build the assignment Variable mainVar = polynomial.main_var(); @@ -108,14 +109,16 @@ RealRootsResult real_roots( } } + lp_interval_t* inter_poly = to_libpoly_interval(interval); std::vector res; for (std::size_t i = 0; i < roots_size; ++i) { - if (lp_interval_contains(inter_poly.get_internal(), &roots[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()); } } - + lp_interval_destruct(inter_poly); + delete inter_poly; delete[] roots; std::sort(res.begin(), res.end()); diff --git a/src/carl-arith/ran/libpoly/helper.h b/src/carl-arith/ran/libpoly/helper.h index 5e394d96..f2fd138c 100644 --- a/src/carl-arith/ran/libpoly/helper.h +++ b/src/carl-arith/ran/libpoly/helper.h @@ -22,51 +22,6 @@ 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); -} - -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 mpz_class to_integer(const lp_integer_t* m) { return *reinterpret_cast(m); } @@ -111,44 +66,106 @@ inline mpq_class to_rational(const lp_value_t* m){ } } +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); + + 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); +} + + +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); +} + /** * 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); + + 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/tests/carl-arith-ran/Test_Libpoly.cpp b/src/tests/carl-arith-ran/Test_Libpoly.cpp index 9464db84..f3ff54f4 100644 --- a/src/tests/carl-arith-ran/Test_Libpoly.cpp +++ b/src/tests/carl-arith-ran/Test_Libpoly.cpp @@ -26,14 +26,14 @@ TEST(LIBPOLY, convertToLibpolyUnivariate) { 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, 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"); From c36d2923cc6f9ff810ad740e822e9bda0c26a174 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Tue, 12 Mar 2024 16:49:05 +0100 Subject: [PATCH 10/29] more efficient conversion --- src/carl-arith/ran/libpoly/LPRan.cpp | 30 +++++++++++------- src/carl-arith/ran/libpoly/RealRoots.cpp | 40 +++++++++++++++--------- src/carl-arith/ran/libpoly/helper.h | 26 ++++++++------- 3 files changed, 59 insertions(+), 37 deletions(-) diff --git a/src/carl-arith/ran/libpoly/LPRan.cpp b/src/carl-arith/ran/libpoly/LPRan.cpp index 91f650c2..f769ebfc 100644 --- a/src/carl-arith/ran/libpoly/LPRan.cpp +++ b/src/carl-arith/ran/libpoly/LPRan.cpp @@ -22,29 +22,37 @@ LPRealAlgebraicNumber::LPRealAlgebraicNumber(lp_value_t&& num) : m_content(std:: } LPRealAlgebraicNumber::LPRealAlgebraicNumber(const carl::UnivariatePolynomial& p, const Interval& i) : LPRealAlgebraicNumber() { - CARL_LOG_DEBUG("carl.ran.libpoly", " Create safe from poly: " << p << " in interval: " << i); + CARL_LOG_DEBUG("carl.ran.libpoly", " Create from poly: " << p << " in interval: " << i); - poly::UPolynomial upoly = to_libpoly_upolynomial(p); + 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: "); + 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); - std::vector roots = poly::isolate_real_roots(upoly); - std::vector res; - for (const auto& val : roots) { + bool found = false; + for (std::size_t i = 0; i < roots_size; ++i) { lp_value_t tmp; - lp_value_construct(&tmp, LP_VALUE_ALGEBRAIC, val.get_internal()); + lp_value_construct(&tmp, LP_VALUE_ALGEBRAIC, &roots[i]); if (lp_interval_contains(inter_poly, &tmp)) { - CARL_LOG_DEBUG("carl.ran.libpoly", " Found Root " << val); - res.emplace_back(val); + *get_internal() = tmp; + found = true; + break; } lp_value_destruct(&tmp); } - assert(res.size() == 1); + assert(found); + + for (std::size_t i = 0; i < roots_size; ++i) { + lp_algebraic_number_destruct(&roots[i]); + } + delete[] roots; + + lp_upolynomial_delete(upoly); lp_interval_destruct(inter_poly); delete inter_poly; - lp_value_construct(get_internal(), LP_VALUE_ALGEBRAIC, res.front().get_internal()); } /** diff --git a/src/carl-arith/ran/libpoly/RealRoots.cpp b/src/carl-arith/ran/libpoly/RealRoots.cpp index f787f2fc..323cacb0 100644 --- a/src/carl-arith/ran/libpoly/RealRoots.cpp +++ b/src/carl-arith/ran/libpoly/RealRoots.cpp @@ -14,7 +14,6 @@ RealRootsResult real_roots( 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(); @@ -23,32 +22,39 @@ RealRootsResult real_roots( return RealRootsResult::no_roots_response(); } - // actual calculations - std::vector roots = poly::isolate_real_roots(poly::to_univariate(polynomial.get_polynomial())); + lp_upolynomial_t* upoly = lp_polynomial_to_univariate(polynomial.get_internal()); - if (roots.empty()) { + 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_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) { + for (std::size_t i = 0; i < roots_size; ++i) { lp_value_t tmp; - lp_value_construct(&tmp, LP_VALUE_ALGEBRAIC, val.get_internal()); + 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)); } @@ -58,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); } @@ -83,7 +89,7 @@ RealRootsResult real_roots( lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); CARL_LOG_TRACE("carl.ran.libpoly", "Call libpoly"); - lp_value_t* roots = new lp_value_t[poly::degree(polynomial.get_polynomial())]; + 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); @@ -98,11 +104,15 @@ RealRootsResult real_roots( auto eval_val = lp_polynomial_evaluate(polynomial.get_internal(), &assignment); //CARL_LOG_DEBUG("carl.ran.libpoly", " Got eval_val " << eval_val); - if (lp_value_cmp(eval_val, poly::Value(long(0)).get_internal()) == 0) { + 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(); @@ -115,6 +125,8 @@ RealRootsResult real_roots( 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); diff --git a/src/carl-arith/ran/libpoly/helper.h b/src/carl-arith/ran/libpoly/helper.h index f2fd138c..f57aedf9 100644 --- a/src/carl-arith/ran/libpoly/helper.h +++ b/src/carl-arith/ran/libpoly/helper.h @@ -87,21 +87,20 @@ inline carl::UnivariatePolynomial to_carl_univariate_polynomial(const template> = dummy> -inline poly::UPolynomial to_libpoly_upolynomial(const carl::UnivariatePolynomial& p) { +inline lp_upolynomial_t* 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()))); + assert(carl::get_denom(p.lcoeff())>0); + lp_upolynomial_t* res; + lp_upolynomial_construct(lp_Z, 0, (lp_integer_t*)mpz_class(carl::get_num(p.lcoeff())).get_mpz_t()); + return res; } - Coeff coprimeFactor = p.coprime_factor(); + 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)); @@ -113,12 +112,15 @@ inline poly::UPolynomial to_libpoly_upolynomial(const carl::UnivariatePolynomial } CARL_LOG_TRACE("carl.converter", "Coprime factor/ denominator: " << denominator); - std::vector coefficients; + std::vector coefficients; + std::vector coefficients_ptr; for (const auto& c : p.coefficients()) { - coefficients.emplace_back(mpz_class(c * denominator)); + coefficients.emplace_back(c * denominator); + coefficients_ptr.emplace_back((lp_integer_t*)coefficients.back().get_mpz_t()); } - CARL_LOG_TRACE("carl.converter", "-> coefficients: " << coefficients); - return poly::UPolynomial(coefficients); + lp_upolynomial_t* res; + lp_upolynomial_construct(lp_Z, p.coefficients().size()-1, reinterpret_cast(coefficients_ptr.data())); + return res; } /** From d382927cb745e1f4275ecbc0afeb6fb70f0d8a89 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 13 Mar 2024 10:37:07 +0100 Subject: [PATCH 11/29] bugfixes and further refactoring --- src/carl-arith/poly/Conversion.h | 2 +- src/carl-arith/poly/libpoly/Functions.h | 2 +- src/carl-arith/poly/libpoly/LPPolynomial.cpp | 8 ++ src/carl-arith/poly/libpoly/LPPolynomial.h | 79 ++++++++++++-------- src/carl-arith/ran/libpoly/Evaluation.cpp | 4 +- src/carl-arith/ran/libpoly/helper.h | 11 ++- src/tests/carl-arith-ran/Test_Libpoly.cpp | 22 +++--- 7 files changed, 75 insertions(+), 53 deletions(-) diff --git a/src/carl-arith/poly/Conversion.h b/src/carl-arith/poly/Conversion.h index 47146eec..015829c0 100644 --- a/src/carl-arith/poly/Conversion.h +++ b/src/carl-arith/poly/Conversion.h @@ -60,7 +60,7 @@ struct ConvertHelper> { lp_polynomial_add_monomial(res, &t); lp_monomial_destruct(&t); } - return LPPolynomial(context, poly::Polynomial(res)); + return LPPolynomial(res, context); } }; diff --git a/src/carl-arith/poly/libpoly/Functions.h b/src/carl-arith/poly/libpoly/Functions.h index be61f7e2..cb9612a9 100644 --- a/src/carl-arith/poly/libpoly/Functions.h +++ b/src/carl-arith/poly/libpoly/Functions.h @@ -65,7 +65,7 @@ inline LPPolynomial resultant(const LPPolynomial& p, const LPPolynomial& q) { */ 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 + 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())); diff --git a/src/carl-arith/poly/libpoly/LPPolynomial.cpp b/src/carl-arith/poly/libpoly/LPPolynomial.cpp index 1781a7fa..1bf521f2 100644 --- a/src/carl-arith/poly/libpoly/LPPolynomial.cpp +++ b/src/carl-arith/poly/libpoly/LPPolynomial.cpp @@ -51,6 +51,14 @@ LPPolynomial::LPPolynomial(const LPContext& context, poly::Polynomial&& p) assert(context.lp_context() == lp_polynomial_get_context(get_internal())); } +LPPolynomial::LPPolynomial(lp_polynomial_t* p, const LPContext& context) + : 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, long val) : m_context(context) { lp_polynomial_construct_simple(m_poly.get_internal(), context.lp_context(), poly::Integer(val).get_internal(), 0, 0); diff --git a/src/carl-arith/poly/libpoly/LPPolynomial.h b/src/carl-arith/poly/libpoly/LPPolynomial.h index ac3c6107..f918a985 100644 --- a/src/carl-arith/poly/libpoly/LPPolynomial.h +++ b/src/carl-arith/poly/libpoly/LPPolynomial.h @@ -81,6 +81,11 @@ class LPPolynomial { */ LPPolynomial(const LPContext& context, poly::Polynomial&& mainPoly); + /** + * Move constructor. + */ + LPPolynomial(lp_polynomial_t* p, const LPContext& context); + /** * Construct a LPPolynomial with only a integer. * @param mainPoly Libpoly Polynomial. @@ -166,31 +171,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; } @@ -218,21 +236,16 @@ class LPPolynomial { 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 +253,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()))); } /** @@ -449,7 +462,7 @@ 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) { @@ -462,7 +475,7 @@ inline bool is_zero(const LPPolynomial& p) { * @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()); } /** @@ -473,9 +486,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; } /** @@ -490,14 +507,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/ran/libpoly/Evaluation.cpp b/src/carl-arith/ran/libpoly/Evaluation.cpp index fdf68bfd..e10a769a 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -54,12 +54,10 @@ boost::tribool evaluate(const BasicConstraint& constraint, const s } //denominator can be omitted - auto poly_pol = constraint.lhs().get_polynomial().get_internal(); + auto poly_pol = constraint.lhs().get_internal(); lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); - // std::cout << lp_assignment_to_string(&assignment) << std::endl; - bool result = false; switch (constraint.relation()) { case Relation::EQ: diff --git a/src/carl-arith/ran/libpoly/helper.h b/src/carl-arith/ran/libpoly/helper.h index f57aedf9..e9e326d5 100644 --- a/src/carl-arith/ran/libpoly/helper.h +++ b/src/carl-arith/ran/libpoly/helper.h @@ -93,8 +93,7 @@ inline lp_upolynomial_t* to_libpoly_upolynomial(const carl::UnivariatePolynomial 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, (lp_integer_t*)mpz_class(carl::get_num(p.lcoeff())).get_mpz_t()); + lp_upolynomial_t* res = lp_upolynomial_construct(lp_Z, 0, mpz_class(carl::get_num(p.lcoeff())).get_mpz_t()); return res; } @@ -113,13 +112,13 @@ inline lp_upolynomial_t* to_libpoly_upolynomial(const carl::UnivariatePolynomial CARL_LOG_TRACE("carl.converter", "Coprime factor/ denominator: " << denominator); std::vector coefficients; - std::vector coefficients_ptr; + std::vector coefficients_ptr; for (const auto& c : p.coefficients()) { coefficients.emplace_back(c * denominator); - coefficients_ptr.emplace_back((lp_integer_t*)coefficients.back().get_mpz_t()); + coefficients_ptr.emplace_back(*coefficients.back().get_mpz_t()); } - lp_upolynomial_t* res; - lp_upolynomial_construct(lp_Z, p.coefficients().size()-1, reinterpret_cast(coefficients_ptr.data())); + 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; } diff --git a/src/tests/carl-arith-ran/Test_Libpoly.cpp b/src/tests/carl-arith-ran/Test_Libpoly.cpp index f3ff54f4..50e7e482 100644 --- a/src/tests/carl-arith-ran/Test_Libpoly.cpp +++ b/src/tests/carl-arith-ran/Test_Libpoly.cpp @@ -14,17 +14,17 @@ 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, 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"); From 443098ce48eb0c58f62ed68abf853e24cee3475b Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 13 Mar 2024 13:16:18 +0100 Subject: [PATCH 12/29] remove usages of libpolyxx --- cmake/FindPoly.cmake | 57 +++--- cmake/export.cmake | 4 +- resources/resources.cmake | 4 +- src/carl-arith/CMakeLists.txt | 6 +- src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h | 26 ++- src/carl-arith/poly/libpoly/Functions.cpp | 86 ++++++++ src/carl-arith/poly/libpoly/Functions.h | 102 +++------- src/carl-arith/poly/libpoly/LPContext.h | 2 +- src/carl-arith/poly/libpoly/LPPolynomial.cpp | 204 ++++++++++--------- src/carl-arith/poly/libpoly/LPPolynomial.h | 38 +--- src/carl-arith/poly/libpoly/LPVariables.h | 3 +- src/carl-arith/poly/libpoly/helper.h | 14 +- src/carl-arith/ran/libpoly/LPAssignment.h | 2 +- src/carl-arith/ran/libpoly/helper.h | 4 +- 14 files changed, 291 insertions(+), 261 deletions(-) create mode 100644 src/carl-arith/poly/libpoly/Functions.cpp diff --git a/cmake/FindPoly.cmake b/cmake/FindPoly.cmake index 7fd4b012..08c91f57 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 @@ -31,13 +31,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 +53,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 +78,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/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/poly/libpoly/CoCoAAdaptorLP.h b/src/carl-arith/poly/libpoly/CoCoAAdaptorLP.h index b5227269..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,26 +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; lp_variable_t var = mContext.lp_variable(mSymbolBack[i]); - poly::Polynomial polyVar = poly_helper::construct_poly(mContext.lp_context(), var); - termPoly *= poly::pow(polyVar, (unsigned int)exponents[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 cb9612a9..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())); + 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.h b/src/carl-arith/poly/libpoly/LPContext.h index 065629d5..fd5b2920 100644 --- a/src/carl-arith/poly/libpoly/LPContext.h +++ b/src/carl-arith/poly/libpoly/LPContext.h @@ -6,10 +6,10 @@ #include "../../core/Variable.h" #include +#include #include #include -#include #include diff --git a/src/carl-arith/poly/libpoly/LPPolynomial.cpp b/src/carl-arith/poly/libpoly/LPPolynomial.cpp index 1bf521f2..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,111 +8,94 @@ 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())); -} - -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())); - assert(context.lp_context() == lp_polynomial_get_context(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(lp_polynomial_t* p, const LPContext& context) - : m_poly(p), m_context(context) { - //lp_polynomial_set_external(m_poly.get_internal()); - - assert(lp_polynomial_check_order(m_poly.get_internal())); + : 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); 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); auto pow = coefficients.size(); @@ -119,45 +103,59 @@ LPPolynomial::LPPolynomial(const LPContext& context, const Variable& mainVar, co 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); 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); - 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); + 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); @@ -188,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) { @@ -285,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; } @@ -297,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; } @@ -309,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; } @@ -345,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 { @@ -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 { @@ -553,7 +555,7 @@ 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; } diff --git a/src/carl-arith/poly/libpoly/LPPolynomial.h b/src/carl-arith/poly/libpoly/LPPolynomial.h index f918a985..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; @@ -67,20 +67,6 @@ 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. - */ - LPPolynomial(const LPContext& context, poly::Polynomial&& mainPoly); - /** * Move constructor. */ @@ -232,7 +218,7 @@ 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; } @@ -263,7 +249,7 @@ class LPPolynomial { * @return Libpoly Polynomial. */ lp_polynomial_t* get_internal() { - return m_poly.get_internal(); + return m_internal; } /** @@ -271,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; } /** @@ -465,11 +442,6 @@ inline bool is_zero(const LPPolynomial& p) { 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. diff --git a/src/carl-arith/poly/libpoly/LPVariables.h b/src/carl-arith/poly/libpoly/LPVariables.h index b169d088..be6661f1 100644 --- a/src/carl-arith/poly/libpoly/LPVariables.h +++ b/src/carl-arith/poly/libpoly/LPVariables.h @@ -8,9 +8,8 @@ #include #include "../../core/Variable.h" #include +#include #include -#include - namespace carl { diff --git a/src/carl-arith/poly/libpoly/helper.h b/src/carl-arith/poly/libpoly/helper.h index 2dd831dc..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) { +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, poly::Integer(1).get_internal(), v, 1); - return poly::Polynomial(p); + 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) { +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_internal(), lp_variable_null, 0) ; - return poly::Polynomial(p); + lp_polynomial_construct_simple(p, c, i.get_mpz_t(), lp_variable_null, 0) ; + return p; } } diff --git a/src/carl-arith/ran/libpoly/LPAssignment.h b/src/carl-arith/ran/libpoly/LPAssignment.h index 7a9d9bf9..4ea19ccb 100644 --- a/src/carl-arith/ran/libpoly/LPAssignment.h +++ b/src/carl-arith/ran/libpoly/LPAssignment.h @@ -9,7 +9,7 @@ #include "../../core/Variable.h" #include #include -#include +#include #include #include "LPRan.h" diff --git a/src/carl-arith/ran/libpoly/helper.h b/src/carl-arith/ran/libpoly/helper.h index e9e326d5..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,7 +16,7 @@ #include #include -#include +#include namespace carl { From 11af43c3a44ffe9675dd83e24dc0c069741c318c Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 13 Mar 2024 13:18:31 +0100 Subject: [PATCH 13/29] include --- src/carl-arith/poly/libpoly/LPContext.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/carl-arith/poly/libpoly/LPContext.h b/src/carl-arith/poly/libpoly/LPContext.h index fd5b2920..7ca640dd 100644 --- a/src/carl-arith/poly/libpoly/LPContext.h +++ b/src/carl-arith/poly/libpoly/LPContext.h @@ -7,6 +7,7 @@ #include #include +#include #include #include From 497e6a169d51108b8f51ba3736297b221e937045 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 15 Mar 2024 09:40:29 +0100 Subject: [PATCH 14/29] libpoly version --- cmake/FindPoly.cmake | 1 + cmake/patches/libpoly_variable_db.patch | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cmake/FindPoly.cmake b/cmake/FindPoly.cmake index 08c91f57..93242a8b 100644 --- a/cmake/FindPoly.cmake +++ b/cmake/FindPoly.cmake @@ -21,6 +21,7 @@ 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 CMAKE_ARGS -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} diff --git a/cmake/patches/libpoly_variable_db.patch b/cmake/patches/libpoly_variable_db.patch index 93ebc768..66f4e367 100644 --- a/cmake/patches/libpoly_variable_db.patch +++ b/cmake/patches/libpoly_variable_db.patch @@ -1,8 +1,8 @@ diff --git a/include/polynomial.h b/include/polynomial.h -index 8a9533f..12123e8 100644 +index 734bfa9..6a2cf75 100644 --- a/include/polynomial.h +++ b/include/polynomial.h -@@ -109,6 +109,8 @@ int lp_polynomial_lc_sgn(const lp_polynomial_t* A); +@@ -109,6 +109,8 @@ void lp_polynomial_lc_constant(const lp_polynomial_t* A, lp_integer_t *out); /** Get the context of the given polynomial */ const lp_polynomial_context_t* lp_polynomial_get_context(const lp_polynomial_t* A); From a35b22345340739b134a808f5f386252e1527fce Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 15 Mar 2024 10:03:00 +0100 Subject: [PATCH 15/29] set UPDATE_COMMAND to "" --- cmake/FindPoly.cmake | 1 + 1 file changed, 1 insertion(+) diff --git a/cmake/FindPoly.cmake b/cmake/FindPoly.cmake index 93242a8b..1459fc5b 100644 --- a/cmake/FindPoly.cmake +++ b/cmake/FindPoly.cmake @@ -23,6 +23,7 @@ if(NOT LIBPOLY_FOUND_SYSTEM) GIT_REPOSITORY https://github.com/SRI-CSL/libpoly GIT_TAG v0.1.13 PATCH_COMMAND git reset --hard + UPDATE_COMMAND "" COMMAND git apply ${CMAKE_SOURCE_DIR}/cmake/patches/libpoly_variable_db.patch CMAKE_ARGS -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -DCMAKE_INSTALL_PREFIX= From 36b616fd78226dba337c5e4810146afcbe908078 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 15 Mar 2024 12:52:47 +0100 Subject: [PATCH 16/29] path libpoly for faster rri --- cmake/FindPoly.cmake | 9 +++---- ...ibpoly_variable_db.patch => libpoly.patch} | 24 +++++++++++++++++-- 2 files changed, 27 insertions(+), 6 deletions(-) rename cmake/patches/{libpoly_variable_db.patch => libpoly.patch} (60%) diff --git a/cmake/FindPoly.cmake b/cmake/FindPoly.cmake index 1459fc5b..e424d1d7 100644 --- a/cmake/FindPoly.cmake +++ b/cmake/FindPoly.cmake @@ -22,9 +22,10 @@ if(NOT LIBPOLY_FOUND_SYSTEM) LIBPOLY-EP GIT_REPOSITORY https://github.com/SRI-CSL/libpoly GIT_TAG v0.1.13 - PATCH_COMMAND git reset --hard - UPDATE_COMMAND "" - COMMAND git apply ${CMAKE_SOURCE_DIR}/cmake/patches/libpoly_variable_db.patch + # PATCH_COMMAND git reset --hard + # 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 @@ -33,7 +34,7 @@ if(NOT LIBPOLY_FOUND_SYSTEM) -DGMP_INCLUDE_DIR=${GMP_LIBRARY_DIR} -DGMP_LIBRARY=${GMP_LIB} BUILD_COMMAND ${CMAKE_COMMAND} - COMMAND ${CMAKE_MAKE_PROGRAM} poly static_poly static_pic_poly # polyxx static_polyxx 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} diff --git a/cmake/patches/libpoly_variable_db.patch b/cmake/patches/libpoly.patch similarity index 60% rename from cmake/patches/libpoly_variable_db.patch rename to cmake/patches/libpoly.patch index 66f4e367..00e5470a 100644 --- a/cmake/patches/libpoly_variable_db.patch +++ b/cmake/patches/libpoly.patch @@ -1,8 +1,8 @@ diff --git a/include/polynomial.h b/include/polynomial.h -index 734bfa9..6a2cf75 100644 +index 8a9533f..12123e8 100644 --- a/include/polynomial.h +++ b/include/polynomial.h -@@ -109,6 +109,8 @@ void lp_polynomial_lc_constant(const lp_polynomial_t* A, lp_integer_t *out); +@@ -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); @@ -11,6 +11,26 @@ index 734bfa9..6a2cf75 100644 /** 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); +diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c +index 11d948f..68f4a8d 100644 +--- a/src/polynomial/polynomial.c ++++ b/src/polynomial/polynomial.c +@@ -1031,6 +1031,15 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t + lp_value_construct_none(&x_value_backup); + } + ++ { ++ 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); ++ A = lp_polynomial_new_from_coefficient(A->ctx, &A_rat); ++ } ++ + size_t i; + + lp_polynomial_t** factors = 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 From b36adbac65f08efb890f5db6cdf5ec157bf420e2 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 15 Mar 2024 12:59:45 +0100 Subject: [PATCH 17/29] fix --- cmake/patches/libpoly.patch | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/cmake/patches/libpoly.patch b/cmake/patches/libpoly.patch index 00e5470a..63900bd2 100644 --- a/cmake/patches/libpoly.patch +++ b/cmake/patches/libpoly.patch @@ -12,25 +12,40 @@ index 8a9533f..12123e8 100644 void lp_polynomial_get_variables(const lp_polynomial_t* A, lp_variable_list_t* vars); diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c -index 11d948f..68f4a8d 100644 +index 11d948f..9869451 100644 --- a/src/polynomial/polynomial.c +++ b/src/polynomial/polynomial.c -@@ -1031,6 +1031,15 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t +@@ -1031,6 +1031,17 @@ 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); -+ A = lp_polynomial_new_from_coefficient(A->ctx, &A_rat); ++ B = lp_polynomial_new_from_coefficient(A->ctx, &A_rat); ++ coefficient_destruct(&A_rat); + } + size_t i; lp_polynomial_t** factors = 0; +@@ -1044,8 +1055,10 @@ 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)); ++ lp_polynomial_reductum_m(&A_r, B, M); ++ assert(x == lp_polynomial_top_variable(B)); ++ ++ lp_polynomial_delete(B); + + // Get the square-free factorization + lp_polynomial_factor_square_free(&A_r, &factors, &multiplicities, &factors_size); diff --git a/src/variable/variable_db.c b/src/variable/variable_db.c index 60f3df4..33866c4 100644 --- a/src/variable/variable_db.c From a104ab0af818bf67cea257d3c8679dba3cec3013 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 15 Mar 2024 13:53:02 +0100 Subject: [PATCH 18/29] fix FindPoly --- cmake/FindPoly.cmake | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cmake/FindPoly.cmake b/cmake/FindPoly.cmake index e424d1d7..0c451da3 100644 --- a/cmake/FindPoly.cmake +++ b/cmake/FindPoly.cmake @@ -22,10 +22,10 @@ if(NOT LIBPOLY_FOUND_SYSTEM) LIBPOLY-EP GIT_REPOSITORY https://github.com/SRI-CSL/libpoly GIT_TAG v0.1.13 - # PATCH_COMMAND git reset --hard + PATCH_COMMAND git reset --hard # UPDATE_COMMAND "" - # COMMAND git apply ${CMAKE_SOURCE_DIR}/cmake/patches/libpoly.patch - PATCH_COMMAND git apply ${CMAKE_SOURCE_DIR}/cmake/patches/libpoly.patch + 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 From 7f937c55aa8c19f916720afed4c8f4f94e09fe3c Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 15 Mar 2024 14:33:15 +0100 Subject: [PATCH 19/29] fix --- cmake/patches/libpoly.patch | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/cmake/patches/libpoly.patch b/cmake/patches/libpoly.patch index 63900bd2..4966c117 100644 --- a/cmake/patches/libpoly.patch +++ b/cmake/patches/libpoly.patch @@ -12,7 +12,7 @@ index 8a9533f..12123e8 100644 void lp_polynomial_get_variables(const lp_polynomial_t* A, lp_variable_list_t* vars); diff --git a/src/polynomial/polynomial.c b/src/polynomial/polynomial.c -index 11d948f..9869451 100644 +index 11d948f..e507f52 100644 --- a/src/polynomial/polynomial.c +++ b/src/polynomial/polynomial.c @@ -1031,6 +1031,17 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t @@ -33,19 +33,25 @@ index 11d948f..9869451 100644 size_t i; lp_polynomial_t** factors = 0; -@@ -1044,8 +1055,10 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t +@@ -1044,11 +1055,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)); -+ 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); ++ 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); - // Get the square-free factorization - lp_polynomial_factor_square_free(&A_r, &factors, &multiplicities, &factors_size); + // Count the max number of roots + size_t total_degree = 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 From 574a621ba596e976766e10323f05106733931dff Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Mon, 18 Mar 2024 10:56:49 +0100 Subject: [PATCH 20/29] better evaluation --- cmake/patches/libpoly.patch | 41 +++++++++-- src/carl-arith/ran/libpoly/Evaluation.cpp | 78 ++++++++------------- src/carl-arith/ran/libpoly/LPAssignment.cpp | 2 - src/carl-arith/ran/libpoly/RealRoots.cpp | 2 + 4 files changed, 68 insertions(+), 55 deletions(-) diff --git a/cmake/patches/libpoly.patch b/cmake/patches/libpoly.patch index 4966c117..76407d2d 100644 --- a/cmake/patches/libpoly.patch +++ b/cmake/patches/libpoly.patch @@ -1,5 +1,5 @@ diff --git a/include/polynomial.h b/include/polynomial.h -index 8a9533f..12123e8 100644 +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); @@ -11,11 +11,20 @@ index 8a9533f..12123e8 100644 /** 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..e507f52 100644 +index 11d948f..7468ab6 100644 --- a/src/polynomial/polynomial.c +++ b/src/polynomial/polynomial.c -@@ -1031,6 +1031,17 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t +@@ -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); } @@ -28,12 +37,13 @@ index 11d948f..e507f52 100644 + 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 +1055,13 @@ void lp_polynomial_roots_isolate(const lp_polynomial_t* A, const lp_assignment_t +@@ -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); @@ -52,6 +62,29 @@ index 11d948f..e507f52 100644 // 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 diff --git a/src/carl-arith/ran/libpoly/Evaluation.cpp b/src/carl-arith/ran/libpoly/Evaluation.cpp index e10a769a..aea5be17 100644 --- a/src/carl-arith/ran/libpoly/Evaluation.cpp +++ b/src/carl-arith/ran/libpoly/Evaluation.cpp @@ -23,65 +23,45 @@ std::optional evaluate( return ran; } -boost::tribool evaluate(const BasicConstraint& constraint, const std::map& evalMap) { +inline auto lp_sign(carl::Relation rel) { + switch (rel) { + case Relation::EQ: + return LP_SGN_EQ_0; + case Relation::NEQ: + return LP_SGN_NE_0; + case Relation::LESS: + return LP_SGN_LT_0; + case Relation::LEQ: + return LP_SGN_LE_0; + case Relation::GREATER: + return LP_SGN_GT_0; + case Relation::GEQ: + return LP_SGN_GE_0; + default: + assert(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); - //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; - } + return carl::evaluate(constraint.lhs().constant_part(), constraint.relation()); } - for (const auto& v : carl::variables(constraint)) { - if (evalMap.find(v) == evalMap.end()) return boost::indeterminate; - } - - //denominator can be omitted auto poly_pol = constraint.lhs().get_internal(); - lp_assignment_t& assignment = LPAssignment::getInstance().get(evalMap); - bool result = false; - switch (constraint.relation()) { - case Relation::EQ: - result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_EQ_0, &assignment); - break; - case Relation::NEQ: - result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_NE_0, &assignment); - break; - case Relation::LESS: - result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_LT_0, &assignment); - break; - case Relation::LEQ: - result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_LE_0, &assignment); - break; - case Relation::GREATER: - result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_GT_0, &assignment); - break; - case Relation::GEQ: - result = lp_polynomial_constraint_evaluate(poly_pol, LP_SGN_GE_0, &assignment); - break; - default: - assert(false); + 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 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 index f7cb13c1..51a850f9 100644 --- a/src/carl-arith/ran/libpoly/LPAssignment.cpp +++ b/src/carl-arith/ran/libpoly/LPAssignment.cpp @@ -32,8 +32,6 @@ lp_assignment_t& LPAssignment::get(const carl::Assignment } return lp_assignment; } - - } } // namespace carl diff --git a/src/carl-arith/ran/libpoly/RealRoots.cpp b/src/carl-arith/ran/libpoly/RealRoots.cpp index 323cacb0..169f70aa 100644 --- a/src/carl-arith/ran/libpoly/RealRoots.cpp +++ b/src/carl-arith/ran/libpoly/RealRoots.cpp @@ -98,10 +98,12 @@ RealRootsResult real_roots( delete[] roots; CARL_LOG_DEBUG("carl.ran.libpoly", " Checking for nullification -> Evaluation at " << mainVar << "= 1"); 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; From 9c4ea0eca09a4b14fe460a9d1d8b604be3e31211 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Mon, 18 Mar 2024 13:52:53 +0100 Subject: [PATCH 21/29] cassert include --- src/carl-statistics/Statistics.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/carl-statistics/Statistics.h b/src/carl-statistics/Statistics.h index 04be9ef0..0b1386b0 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" From 69e138356528d83ecdfc8428ae9ac87eb2951b41 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Tue, 26 Mar 2024 15:28:07 +0100 Subject: [PATCH 22/29] better timing --- src/carl-statistics/Statistics.h | 2 +- src/carl-statistics/Timing.h | 14 ++++++++++++-- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/carl-statistics/Statistics.h b/src/carl-statistics/Statistics.h index 0b1386b0..efa9b163 100644 --- a/src/carl-statistics/Statistics.h +++ b/src/carl-statistics/Statistics.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"); } }; From f1c33f77383c500df4b7fbb77758f410c7961021 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Thu, 28 Mar 2024 17:17:20 +0100 Subject: [PATCH 23/29] multicounter total --- src/carl-statistics/MultiCounter.h | 3 +++ 1 file changed, 3 insertions(+) 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)); } }; From 5075ae06ff6e7f2470046cb76b2162feae5cb67d Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 3 Apr 2024 09:08:48 +0200 Subject: [PATCH 24/29] fix build --- src/tests/carl-arith-poly/Test_LPPolynomialFunctions.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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); From 6efba4a72bf57c671fba6f56e153e8afa943b9ae Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 3 Apr 2024 11:19:19 +0200 Subject: [PATCH 25/29] update gitlab ci --- .gitlab-ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d19af612..defa2652 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -30,7 +30,7 @@ build-gcc12: - cmake/ - src/ -build-gcc11: +.build-gcc11: dependencies: [] stage: build-gcc script: @@ -58,7 +58,7 @@ build-clang14: - cmake/ - src/ -build-clang13: +.build-clang13: dependencies: [] stage: build-clang script: @@ -73,7 +73,7 @@ build-clang13: only: - development -test-clang: +.test-clang: dependencies: [build-clang14] stage: test script: From 624af89759ab8545d6d9b61e755e1a202e20e3c0 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 3 Apr 2024 13:15:04 +0200 Subject: [PATCH 26/29] gitlab ci --- .gitlab-ci.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index defa2652..5d3da32e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -29,6 +29,7 @@ build-gcc12: - build/ - cmake/ - src/ + expire_in: 1 week .build-gcc11: dependencies: [] @@ -42,6 +43,7 @@ build-gcc12: - build/ - cmake/ - src/ + expire_in: 1 week only: - development @@ -57,6 +59,7 @@ build-clang14: - build/ - cmake/ - src/ + expire_in: 1 week .build-clang13: dependencies: [] @@ -70,6 +73,7 @@ build-clang14: - build/ - cmake/ - src/ + expire_in: 1 week only: - development From 70b99793f9cbec85860fc18410e2ab95afa2f301 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 3 Apr 2024 16:25:36 +0200 Subject: [PATCH 27/29] gitlab ci --- .ci/build.sh | 6 +++--- .gitlab-ci.yml | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) 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 5d3da32e..5955eeed 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -22,7 +22,7 @@ build-gcc12: stage: build-gcc 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: @@ -36,7 +36,7 @@ build-gcc12: stage: build-gcc 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: @@ -52,7 +52,7 @@ build-clang14: stage: build-clang 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: @@ -66,7 +66,7 @@ build-clang14: stage: build-clang 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: From 88eecc9a8dd1fc982c230ae1e30252e7d7787fc2 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Wed, 3 Apr 2024 17:44:15 +0200 Subject: [PATCH 28/29] gitlab ci --- .gitlab-ci.yml | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5955eeed..1df3a918 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,14 +12,13 @@ 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=all source .ci/build.sh @@ -33,7 +32,7 @@ build-gcc12: .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=all source .ci/build.sh @@ -49,7 +48,7 @@ build-gcc12: 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=all source .ci/build.sh @@ -63,7 +62,7 @@ build-clang14: .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=all source .ci/build.sh From 45a2501b4a3ecf8c2be351c8dc8d88ba2aeb4f86 Mon Sep 17 00:00:00 2001 From: Jasper Nalbach Date: Fri, 26 Apr 2024 10:55:48 +0200 Subject: [PATCH 29/29] bump version --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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}") # # # # # # # # # # # # # # # # # # # # # #