diff --git a/src/ast/rewriter/bv_rewriter.cpp b/src/ast/rewriter/bv_rewriter.cpp index 516248e24..8333295e5 100644 --- a/src/ast/rewriter/bv_rewriter.cpp +++ b/src/ast/rewriter/bv_rewriter.cpp @@ -37,6 +37,7 @@ void bv_rewriter::updt_local_params(params_ref const & _p) { m_ite2id = p.bv_ite2id(); m_le_extra = p.bv_le_extra(); m_le2extract = p.bv_le2extract(); + m_le2extract = false; // set_sort_sums(p.bv_sort_ac()); } @@ -1724,6 +1725,9 @@ br_status bv_rewriter::mk_bv_or(unsigned num, expr * const * args, expr_ref & re } } + if (!m_le2extract) + return BR_FAILED; + if (!v1.is_zero() && new_args.size() == 1) { v1 = m_util.norm(v1, sz); #ifdef _TRACE diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp index 42d473657..eb2c8b3fd 100644 --- a/src/math/dd/dd_pdd.cpp +++ b/src/math/dd/dd_pdd.cpp @@ -1208,6 +1208,11 @@ namespace dd { return true; } + /** Return true iff p contains no variables other than v. */ + bool pdd_manager::is_univariate_in(PDD p, unsigned v) { + return (is_val(p) || var(p) == v) && is_univariate(p); + } + /** * Push coefficients of univariate polynomial in order of ascending degree. * Example: a*x^2 + b*x + c ==> [ c, b, a ] @@ -1719,7 +1724,7 @@ namespace dd { unsigned pow; if (val.is_power_of_two(pow) && pow > 10) return out << "2^" << pow; - for (int offset : {-1, 1}) + for (int offset : {-2, -1, 1, 2}) if (val < m.max_value() && (val - offset).is_power_of_two(pow) && pow > 10) return out << lparen() << "2^" << pow << (offset >= 0 ? "+" : "") << offset << rparen(); rational neg_val = mod(-val, m.two_to_N()); diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h index de29ed600..12900b9c6 100644 --- a/src/math/dd/dd_pdd.h +++ b/src/math/dd/dd_pdd.h @@ -364,6 +364,7 @@ namespace dd { bool is_monomial(PDD p); bool is_univariate(PDD p); + bool is_univariate_in(PDD p, unsigned v); void get_univariate_coefficients(PDD p, vector& coeff); // create an spoly r if leading monomials of a and b overlap @@ -430,6 +431,7 @@ namespace dd { bool is_binary() const { return m.is_binary(root); } bool is_monomial() const { return m.is_monomial(root); } bool is_univariate() const { return m.is_univariate(root); } + bool is_univariate_in(unsigned v) const { return m.is_univariate_in(root, v); } void get_univariate_coefficients(vector& coeff) const { m.get_univariate_coefficients(root, coeff); } vector get_univariate_coefficients() const { vector coeff; m.get_univariate_coefficients(root, coeff); return coeff; } bool is_never_zero() const { return m.is_never_zero(root); } diff --git a/src/math/polysat/assignment.h b/src/math/polysat/assignment.h index 2f8f65474..fb678934b 100644 --- a/src/math/polysat/assignment.h +++ b/src/math/polysat/assignment.h @@ -53,13 +53,13 @@ namespace polysat { public: substitution(dd::pdd_manager& m); - substitution add(pvar var, rational const& value) const; - pdd apply_to(pdd const& p) const; + [[nodiscard]] substitution add(pvar var, rational const& value) const; + [[nodiscard]] pdd apply_to(pdd const& p) const; - bool contains(pvar var) const; - bool value(pvar var, rational& out_value) const; + [[nodiscard]] bool contains(pvar var) const; + [[nodiscard]] bool value(pvar var, rational& out_value) const; - bool empty() const { return m_subst.is_one(); } + [[nodiscard]] bool empty() const { return m_subst.is_one(); } pdd const& to_pdd() const { return m_subst; } unsigned bit_width() const { return to_pdd().power_of_2(); } @@ -77,7 +77,7 @@ namespace polysat { /** Full variable assignment, may include variables of varying bit widths. */ class assignment { solver* m_solver; - vector m_pairs; // TODO: do we still need this? + vector m_pairs; mutable scoped_ptr_vector m_subst; vector m_subst_trail; diff --git a/src/math/polysat/boolean.cpp b/src/math/polysat/boolean.cpp index 999297bd1..d48a1fa7b 100644 --- a/src/math/polysat/boolean.cpp +++ b/src/math/polysat/boolean.cpp @@ -72,7 +72,7 @@ namespace polysat { } void bool_var_manager::eval(sat::literal lit, unsigned lvl) { - LOG_V("Evaluate " << lit << " @ " << lvl); + LOG_V(10, "Evaluate " << lit << " @ " << lvl); assign(kind_t::evaluation, lit, lvl, nullptr); SASSERT(is_evaluation(lit)); } diff --git a/src/math/polysat/clause_builder.cpp b/src/math/polysat/clause_builder.cpp index 511e341b4..e3b504c77 100644 --- a/src/math/polysat/clause_builder.cpp +++ b/src/math/polysat/clause_builder.cpp @@ -75,15 +75,23 @@ namespace polysat { m_literals.push_back(c.blit()); } - void clause_builder::insert_eval(sat::literal lit, bool status) { - insert_eval(m_solver->lit2cnstr(lit), status); + void clause_builder::insert_eval(sat::literal lit) { + insert_eval(m_solver->lit2cnstr(lit)); } - void clause_builder::insert_eval(signed_constraint c, bool status) { - if (c.bvalue(*m_solver) == l_undef) { - sat::literal lit = c.blit(); - m_solver->assign_eval(status ? lit : ~lit); - } + void clause_builder::insert_eval(signed_constraint c) { + if (c.bvalue(*m_solver) == l_undef) + m_solver->assign_eval(~c.blit()); + insert(c); + } + + void clause_builder::insert_try_eval(sat::literal lit) { + insert_eval(m_solver->lit2cnstr(lit)); + } + + void clause_builder::insert_try_eval(signed_constraint c) { + if (c.bvalue(*m_solver) == l_undef && c.is_currently_false(*m_solver)) + m_solver->assign_eval(~c.blit()); insert(c); } } diff --git a/src/math/polysat/clause_builder.h b/src/math/polysat/clause_builder.h index 4c6294138..c9e4ef8e8 100644 --- a/src/math/polysat/clause_builder.h +++ b/src/math/polysat/clause_builder.h @@ -51,19 +51,20 @@ namespace polysat { void insert(signed_constraint c); void insert(inequality const& i) { insert(i.as_signed_constraint()); } - /// Insert constraints into the clause. - /// If they are not yet on the search stack, add them as evaluated to the given status. - /// \pre Constraint must evaluate to true or false according to the given status in the current assignment. - void insert_eval(sat::literal lit, bool status); - void insert_eval(signed_constraint c, bool status); - /// Insert constraints into the clause. /// If they are not yet on the search stack, add them as evaluated to \b false. /// \pre Constraint must be \b false in the current assignment. - void insert_eval(sat::literal lit) { insert_eval(lit, false); } - void insert_eval(signed_constraint c) { insert_eval(c, false); } + void insert_eval(sat::literal lit); + void insert_eval(signed_constraint c); void insert_eval(inequality const& i) { insert_eval(i.as_signed_constraint()); } + /// Insert constraints into the clause. + /// If possible, evaluate them as in insert_eval. + /// Evaluation might not be possible if not all variables in the constraint are assigned. + /// \pre Constraint must be \b non-true in the current assignment. + void insert_try_eval(sat::literal lit); + void insert_try_eval(signed_constraint lit); + using const_iterator = decltype(m_literals)::const_iterator; const_iterator begin() const { return m_literals.begin(); } const_iterator end() const { return m_literals.end(); } diff --git a/src/math/polysat/conflict.cpp b/src/math/polysat/conflict.cpp index d605b9ca8..dc2c548d0 100644 --- a/src/math/polysat/conflict.cpp +++ b/src/math/polysat/conflict.cpp @@ -186,6 +186,8 @@ namespace polysat { SASSERT(s.at_base_level()); m_level = s.m_level; SASSERT(!empty()); + // TODO: logger().begin_conflict??? + // TODO: check uses of logger().begin_conflict(). sometimes we call it before adding constraints, sometimes after... } void conflict::init(signed_constraint c) { @@ -193,17 +195,6 @@ namespace polysat { SASSERT(empty()); m_level = s.m_level; m_narrow_queue.push_back(c.blit()); // if the conflict is only due to a missed propagation of c - set_impl(c); - logger().begin_conflict(); - } - - void conflict::set(signed_constraint c) { - SASSERT(!empty()); - remove_all(); - set_impl(c); - } - - void conflict::set_impl(signed_constraint c) { if (c.bvalue(s) == l_false) { // boolean conflict // This case should not happen: @@ -219,6 +210,7 @@ namespace polysat { insert_vars(c); } SASSERT(!empty()); + logger().begin_conflict(); } void conflict::init(clause const& cl) { @@ -234,31 +226,24 @@ namespace polysat { logger().begin_conflict(); } - void conflict::init(pvar v, bool by_viable_fallback) { - LOG("Conflict: viable v" << v); + void conflict::init_by_viable_interval(pvar v) { + LOG("Conflict: viable_interval v" << v); SASSERT(empty()); + SASSERT(!s.is_assigned(v)); m_level = s.m_level; - if (by_viable_fallback) { - logger().begin_conflict(header_with_var("unsat core from viable fallback for v", v)); - // Conflict detected by viable fallback: - // initial conflict is the unsat core of the univariate solver, - // and the assignment (under which the constraints are univariate in v) - // TODO: - // - currently we add variables directly, which is sound: - // e.g.: v^2 + w^2 == 0; w := 1 - // - but we could use side constraints on the coefficients instead (coefficients when viewed as polynomial over v): - // e.g.: v^2 + w^2 == 0; w^2 == 1 - signed_constraints unsat_core = s.m_viable_fallback.unsat_core(v); - for (auto c : unsat_core) { - insert(c); - insert_vars(c); - } - SASSERT(!m_vars.contains(v)); - } - else { - logger().begin_conflict(header_with_var("forbidden interval lemma for v", v)); - VERIFY(s.m_viable.resolve(v, *this)); - } + logger().begin_conflict(header_with_var("viable_interval v", v)); + VERIFY(s.m_viable.resolve_interval(v, *this)); + SASSERT(!empty()); + revert_pvar(v); // at this point, v is not assigned + } + + void conflict::init_by_viable_fallback(pvar v, univariate_solver& us) { + LOG("Conflict: viable_fallback v" << v); + SASSERT(empty()); + SASSERT(!s.is_assigned(v)); + m_level = s.m_level; + logger().begin_conflict(header_with_var("viable_fallback v", v)); + VERIFY(s.m_viable.resolve_fallback(v, us, *this)); SASSERT(!empty()); revert_pvar(v); // at this point, v is not assigned } @@ -318,8 +303,14 @@ namespace polysat { } void conflict::add_lemma(char const* name, clause_ref lemma) { + + for (auto lit : *lemma) + if (s.m_bvars.is_true(lit)) + verbose_stream() << "REDUNDANT lemma " << lit << " : " << show_deref(lemma) << "\n"; + LOG_H3("Lemma " << (name ? name : "") << ": " << show_deref(lemma)); SASSERT(lemma); + s.m_simplify_clause.apply(*lemma); lemma->set_redundant(true); for (sat::literal lit : *lemma) { LOG(lit_pp(s, lit)); diff --git a/src/math/polysat/conflict.h b/src/math/polysat/conflict.h index ac529b382..affdb4ba7 100644 --- a/src/math/polysat/conflict.h +++ b/src/math/polysat/conflict.h @@ -85,7 +85,7 @@ namespace polysat { scoped_ptr m_resolver; // current conflict core consists of m_literals and m_vars - indexed_uint_set m_literals; // set of boolean literals in the conflict + indexed_uint_set m_literals; // set of boolean literals in the conflict; TODO: why not sat::literal_set uint_set m_vars; // variable assignments used as premises, shorthand for literals (x := v) unsigned_vector m_var_occurrences; // for each variable, the number of constraints in m_literals that contain it @@ -102,8 +102,6 @@ namespace polysat { // Level at which the conflict was discovered unsigned m_level = UINT_MAX; - void set_impl(signed_constraint c); - public: conflict(solver& s); ~conflict(); @@ -133,11 +131,10 @@ namespace polysat { void init(signed_constraint c); /** boolean conflict with the given clause */ void init(clause const& cl); - /** conflict because there is no viable value for the variable v */ - void init(pvar v, bool by_viable_fallback); - - /** replace the current conflict by a single constraint */ - void set(signed_constraint c); + /** conflict because there is no viable value for the variable v, by interval reasoning */ + void init_by_viable_interval(pvar v); + /** conflict because there is no viable value for the variable v, by fallback solver */ + void init_by_viable_fallback(pvar v, univariate_solver& us); bool contains(signed_constraint c) const { SASSERT(c); return contains(c.blit()); } bool contains(sat::literal lit) const; diff --git a/src/math/polysat/constraint.h b/src/math/polysat/constraint.h index dc388eb95..68758729e 100644 --- a/src/math/polysat/constraint.h +++ b/src/math/polysat/constraint.h @@ -51,7 +51,7 @@ namespace polysat { /** The boolean variable associated to this constraint */ sat::bool_var m_bvar = sat::null_bool_var; - constraint(constraint_manager& m, ckind_t k): m_kind(k) {} + constraint(ckind_t k): m_kind(k) {} public: virtual ~constraint() {} @@ -115,9 +115,11 @@ namespace polysat { bool is_pwatched() const { return m_is_pwatched; } void set_pwatched(bool f) { m_is_pwatched = f; } - /// Assuming the constraint is univariate under the current assignment of 's', - /// adds the constraint to the univariate solver 'us'. - virtual void add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const = 0; + /** + * If the constraint is univariate in variable 'v' under the current assignment of 's', + * add the constraint to the univariate solver 'us'. + */ + virtual void add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const = 0; }; inline std::ostream& operator<<(std::ostream& out, constraint const& c) { return c.display(out); } @@ -162,7 +164,7 @@ namespace polysat { void narrow(solver& s, bool first) { get()->narrow(s, is_positive(), first); } clause_ref produce_lemma(solver& s, assignment const& a) { return get()->produce_lemma(s, a, is_positive()); } - void add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep) const { get()->add_to_univariate_solver(s, us, dep, is_positive()); } + void add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep) const { get()->add_to_univariate_solver(v, s, us, dep, is_positive()); } unsigned_vector const& vars() const { return m_constraint->vars(); } unsigned var(unsigned idx) const { return m_constraint->var(idx); } diff --git a/src/math/polysat/constraint_manager.cpp b/src/math/polysat/constraint_manager.cpp index 0e182ffda..ae3e383ee 100644 --- a/src/math/polysat/constraint_manager.cpp +++ b/src/math/polysat/constraint_manager.cpp @@ -56,7 +56,7 @@ namespace polysat { /** Add constraint to per-level storage */ void constraint_manager::store(constraint* c) { - LOG_V("Store constraint: " << show_deref(c)); + LOG_V(20, "Store constraint: " << show_deref(c)); m_constraints.push_back(c); } @@ -231,7 +231,7 @@ namespace polysat { } /** Look up constraint among stored constraints. */ - constraint* constraint_manager::dedup(constraint* c1) { + constraint* constraint_manager::dedup_store(constraint* c1) { constraint* c2 = nullptr; if (m_dedup.constraints.find(c1, c2)) { dealloc(c1); @@ -246,6 +246,15 @@ namespace polysat { } } + /** Find stored constraint */ + constraint* constraint_manager::dedup_find(constraint* c1) const { + constraint* c = nullptr; + if (!m_dedup.constraints.find(c1, c)) { + SASSERT(c == nullptr); + } + return c; + } + void constraint_manager::gc() { LOG_H1("gc"); gc_clauses(); @@ -294,7 +303,7 @@ namespace polysat { pdd lhs = a; pdd rhs = b; ule_constraint::simplify(is_positive, lhs, rhs); - return { dedup(alloc(ule_constraint, *this, lhs, rhs)), is_positive }; + return { dedup_store(alloc(ule_constraint, lhs, rhs)), is_positive }; } signed_constraint constraint_manager::eq(pdd const& p) { @@ -305,6 +314,19 @@ namespace polysat { return ~ule(b, a); } + signed_constraint constraint_manager::find_eq(pdd const& p) const { + return find_ule(p, p.manager().zero()); + } + + signed_constraint constraint_manager::find_ule(pdd const& a, pdd const& b) const { + bool is_positive = true; + pdd lhs = a; + pdd rhs = b; + ule_constraint::simplify(is_positive, lhs, rhs); + ule_constraint tmp(lhs, rhs); // TODO: this still allocates ule_constraint::m_vars + return { dedup_find(&tmp), is_positive }; + } + /** * encode that the i'th bit of p is 1. * It holds if p << (K - i - 1) >= 2^{K-1}, where K is the bit-width. @@ -317,19 +339,19 @@ namespace polysat { } signed_constraint constraint_manager::umul_ovfl(pdd const& a, pdd const& b) { - return { dedup(alloc(umul_ovfl_constraint, *this, a, b)), true }; + return { dedup_store(alloc(umul_ovfl_constraint, a, b)), true }; } signed_constraint constraint_manager::smul_ovfl(pdd const& a, pdd const& b) { - return { dedup(alloc(smul_fl_constraint, *this, a, b, true)), true }; + return { dedup_store(alloc(smul_fl_constraint, a, b, true)), true }; } signed_constraint constraint_manager::smul_udfl(pdd const& a, pdd const& b) { - return { dedup(alloc(smul_fl_constraint, *this, a, b, false)), true }; + return { dedup_store(alloc(smul_fl_constraint, a, b, false)), true }; } signed_constraint constraint_manager::mk_op_constraint(op_constraint::code op, pdd const& p, pdd const& q, pdd const& r) { - return { dedup(alloc(op_constraint, *this, op, p, q, r)), true }; + return { dedup_store(alloc(op_constraint, op, p, q, r)), true }; } // To do signed comparison of bitvectors, flip the msb and do unsigned comparison: @@ -405,6 +427,8 @@ namespace polysat { // addition does not overflow in (b*q) + r; for now expressed as: r <= bq+r // b ≠ 0 ==> r < b // b = 0 ==> q = -1 + // TODO: when a,b become evaluable, can we actually propagate q,r? doesn't seem like it. + // Maybe we need something like an op_constraint for better propagation. s.add_clause(eq(b * q + r - a), false); s.add_clause(~umul_ovfl(b, q), false); // r <= b*q+r diff --git a/src/math/polysat/constraint_manager.h b/src/math/polysat/constraint_manager.h index 2ef2a3ca1..2f8e8541d 100644 --- a/src/math/polysat/constraint_manager.h +++ b/src/math/polysat/constraint_manager.h @@ -67,9 +67,9 @@ namespace polysat { constraint* get_bv2c(sat::bool_var bv) const; void store(constraint* c); - void erase(constraint* c); - constraint* dedup(constraint* c); + constraint* dedup_store(constraint* c); + constraint* dedup_find(constraint* c) const; void gc_constraints(); void gc_clauses(); @@ -102,6 +102,10 @@ namespace polysat { constraint* lookup(sat::bool_var var) const; signed_constraint lookup(sat::literal lit) const; + /** Find constraint p == 0; returns null if it doesn't exist yet */ + signed_constraint find_eq(pdd const& p) const; + signed_constraint find_ule(pdd const& a, pdd const& b) const; + signed_constraint eq(pdd const& p); signed_constraint ule(pdd const& a, pdd const& b); signed_constraint ult(pdd const& a, pdd const& b); @@ -126,9 +130,65 @@ namespace polysat { constraint* const* begin() const { return m_constraints.data(); } constraint* const* end() const { return m_constraints.data() + m_constraints.size(); } - using clause_iterator = decltype(m_clauses)::const_iterator; - clause_iterator clauses_begin() const { return m_clauses.begin(); } - clause_iterator clauses_end() const { return m_clauses.end(); } + class clause_iterator { + friend class constraint_manager; + using storage_t = decltype(constraint_manager::m_clauses); + + using outer_iterator = storage_t::const_iterator; + outer_iterator outer_it; + outer_iterator outer_end; + + using inner_iterator = storage_t::data_t::const_iterator; + inner_iterator inner_it; + inner_iterator inner_end; + + clause_iterator(storage_t const& cls, bool begin) { + if (begin) { + outer_it = cls.begin(); + outer_end = cls.end(); + } + else { + outer_it = cls.end(); + outer_end = cls.end(); + } + begin_inner(); + } + + void begin_inner() { + if (outer_it == outer_end) { + inner_it = nullptr; + inner_end = nullptr; + } + else { + inner_it = outer_it->begin(); + inner_end = outer_it->end(); + } + } + + public: + clause const& operator*() const { + return *inner_it->get(); + } + + clause_iterator& operator++() { + ++inner_it; + while (inner_it == inner_end && outer_it != outer_end) { + ++outer_it; + begin_inner(); + } + return *this; + } + + bool operator==(clause_iterator const& other) const { + return outer_it == other.outer_it && inner_it == other.inner_it; + } + + bool operator!=(clause_iterator const& other) const { + return !operator==(other); + } + }; + clause_iterator clauses_begin() const { return clause_iterator(m_clauses, true); } + clause_iterator clauses_end() const { return clause_iterator(m_clauses, false); } class clauses_t { constraint_manager const* m_cm; diff --git a/src/math/polysat/forbidden_intervals.cpp b/src/math/polysat/forbidden_intervals.cpp index 8e8131651..22acc2b12 100644 --- a/src/math/polysat/forbidden_intervals.cpp +++ b/src/math/polysat/forbidden_intervals.cpp @@ -9,7 +9,7 @@ Module Name: Author: - Jakob Rath 2021-04-6 + Jakob Rath 2021-04-06 Nikolaj Bjorner (nbjorner) 2021-03-19 --*/ @@ -21,7 +21,7 @@ Author: namespace polysat { - /** + /** * * \param[in] c Original constraint * \param[in] v Variable that is bounded by constraint @@ -38,7 +38,7 @@ namespace polysat { bool forbidden_intervals::get_interval_umul_ovfl(signed_constraint const& c, pvar v, fi_record& fi) { - backtrack _backtrack(fi.side_cond); + backtrack _backtrack(fi.side_cond); fi.coeff = 1; fi.src = c; @@ -56,7 +56,7 @@ namespace polysat { std::swap(a1, a2); std::swap(e1, e2); std::swap(b1, b2); - std::swap(ok1, ok2); + std::swap(ok1, ok2); } if (ok1 && !ok2 && a1.is_one() && b1.is_zero()) { if (c.is_positive()) { @@ -69,8 +69,8 @@ namespace polysat { return true; } } - - if (!ok1 || !ok2) + + if (!ok1 || !ok2) return false; @@ -129,10 +129,10 @@ namespace polysat { } static char const* _last_function = ""; - + bool forbidden_intervals::get_interval_ule(signed_constraint const& c, pvar v, fi_record& fi) { - - backtrack _backtrack(fi.side_cond); + + backtrack _backtrack(fi.side_cond); fi.coeff = 1; fi.src = c; @@ -166,11 +166,11 @@ namespace polysat { _backtrack.released = true; // v > q - if (ok1 && !ok2 && match_non_zero(c, a1, b1, e1, fi)) + if (false && ok1 && !ok2 && match_non_zero(c, a1, b1, e1, c->to_ule().rhs(), fi)) return true; // p > v - if (!ok1 && ok2 && match_non_max(c, a2, b2, e2, fi)) + if (false && !ok1 && ok2 && match_non_max(c, c->to_ule().lhs(), a2, b2, e2, fi)) return true; if (!ok1 || !ok2 || (a1.is_zero() && a2.is_zero())) { @@ -185,6 +185,11 @@ namespace polysat { if (match_zero(c, a1, b1, e1, a2, b2, e2, fi)) return true; + // -1 <= a*v + b, a odd + // -1 > a*v + b, a odd + if (match_max(c, a1, b1, e1, a2, b2, e2, fi)) + return true; + if (match_linear1(c, a1, b1, e1, a2, b2, e2, fi)) return true; if (match_linear2(c, a1, b1, e1, a2, b2, e2, fi)) @@ -231,7 +236,7 @@ namespace polysat { out_side_cond.push_back(s.eq(q, r)); q = r; } - auto b = s.subst(e); + auto b = s.subst(e); return std::tuple(b.is_val(), q.val(), e, b); }; @@ -239,7 +244,7 @@ namespace polysat { signed_constraint const& c, bool is_trivial, rational & coeff, rational & lo_val, pdd & lo, rational & hi_val, pdd & hi) { - + dd::pdd_manager& m = lo.manager(); if (is_trivial) { @@ -256,7 +261,7 @@ namespace polysat { rational pow2 = m.max_value() + 1; if (coeff > pow2/2) { - + coeff = pow2 - coeff; SASSERT(coeff > 0); // Transform according to: y \in [l;u[ <=> -y \in [1-u;1-l[ @@ -282,10 +287,10 @@ namespace polysat { /** * Match e1 + t <= e2, with t = a1*y - * condition for empty/full: e2 == -1 + * condition for empty/full: e2 == -1 */ bool forbidden_intervals::match_linear1(signed_constraint const& c, - rational const & a1, pdd const& b1, pdd const& e1, + rational const & a1, pdd const& b1, pdd const& e1, rational const & a2, pdd const& b2, pdd const& e2, fi_record& fi) { _last_function = __func__; @@ -382,10 +387,10 @@ namespace polysat { /** * a*v <= 0, a odd - * forbidden interval for v is [1,0[ + * forbidden interval for v is [1;0[ * * a*v + b <= 0, a odd - * forbidden interval for v is [n+1,n[ where n = -b * a^-1 + * forbidden interval for v is [n+1;n[ where n = -b * a^-1 * * TODO: extend to * 2^k*a*v <= 0, a odd @@ -403,7 +408,7 @@ namespace polysat { rational a1_inv; VERIFY(a1.mult_inverse(m.power_of_2(), a1_inv)); - // interval for a*v + b > 0 is [n,n+1[ where n = -b * a^-1 + // interval for a*v + b > 0 is [n;n+1[ where n = -b * a^-1 rational lo_val = mod(-b1.val() * a1_inv, mod_value); pdd lo = -e1 * a1_inv; rational hi_val = mod(lo_val + 1, mod_value); @@ -425,34 +430,83 @@ namespace polysat { return false; } - /** + /** + * -1 <= a*v + b, a odd + * forbidden interval for v is [n+1;n[ where n = (-b-1) * a^-1 + */ + bool forbidden_intervals::match_max( + signed_constraint const& c, + rational const & a1, pdd const& b1, pdd const& e1, + rational const & a2, pdd const& b2, pdd const& e2, + fi_record& fi) { + _last_function = __func__; + if (a1.is_zero() && b1.is_max() && a2.is_odd()) { + auto& m = e2.manager(); + rational const& mod_value = m.two_to_N(); + rational a2_inv; + VERIFY(a2.mult_inverse(m.power_of_2(), a2_inv)); + + // interval for -1 > a*v + b is [n;n+1[ where n = (-b-1) * a^-1 + rational lo_val = mod((-1 - b2.val()) * a2_inv, mod_value); + pdd lo = (-1 - e2) * a2_inv; + rational hi_val = mod(lo_val + 1, mod_value); + pdd hi = lo + 1; + + // interval for -1 <= a*v + b is the complement + if (c.is_positive()) { + std::swap(lo_val, hi_val); + std::swap(lo, hi); + } + + fi.coeff = 1; + fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val); + // LHS == -1 is a precondition because we can only multiply with a^-1 in equations, not inequalities + if (b1 != e1) + fi.side_cond.push_back(s.eq(b1, e1)); + return true; + } + return false; + } + + /** * v > q * forbidden interval for v is [0,1[ - * + * * v - k > q * forbidden interval for v is [k,k+1[ - * + * + * v > q + * forbidden interval for v is [0;q+1[ but at least [0;1[ + * + * The following cases are implemented, and subsume the simple ones above. + * + * v - k > q + * forbidden interval for v is [k;k+q+1[ but at least [k;k+1[ + * * a*v - k > q, a odd * forbidden interval for v is [a^-1*k, a^-1*k + 1[ */ bool forbidden_intervals::match_non_zero( signed_constraint const& c, - rational const & a1, pdd const& b1, pdd const& e1, + rational const& a1, pdd const& b1, pdd const& e1, + pdd const& q, fi_record& fi) { _last_function = __func__; - if (a1.is_odd() && b1.is_zero() && c.is_negative()) { + SASSERT(b1.is_val()); + if (a1.is_one() && c.is_negative()) { + // v - k > q auto& m = e1.manager(); - rational lo_val(0); - auto lo = m.zero(); - rational hi_val(1); - auto hi = m.one(); + rational const& mod_value = m.two_to_N(); + rational lo_val = (-b1).val(); + pdd lo = -e1; + rational hi_val = mod(lo_val + 1, mod_value); + pdd hi = lo + q + 1; fi.coeff = 1; fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val); - if (b1 != e1) - fi.side_cond.push_back(s.eq(b1, e1)); return true; } - if (a1.is_odd() && b1.is_val() && c.is_negative()) { + if (a1.is_odd() && c.is_negative()) { + // a*v - k > q, a odd auto& m = e1.manager(); rational const& mod_value = m.two_to_N(); rational a1_inv; @@ -463,8 +517,6 @@ namespace polysat { auto hi = lo + 1; fi.coeff = 1; fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val); - if (b1 != e1) - fi.side_cond.push_back(s.eq(b1, e1)); return true; } return false; @@ -472,26 +524,45 @@ namespace polysat { /** * p > v - * forbidden interval for v is [-1,0[ + * forbidden interval for v is [p;0[ but at least [-1,0[ + * * p > v + k - * forbidden interval for v is [-k-1,-k[ + * forbidden interval for v is [p-k;-k[ but at least [-1-k,-k[ + * + * p > a*v + k, a odd + * forbidden interval for v is [ a^-1*(-1-k) ; a^-1*(-1-k) + 1 [ */ bool forbidden_intervals::match_non_max( signed_constraint const& c, - rational const & a2, pdd const& b2, pdd const& e2, + pdd const& p, + rational const& a2, pdd const& b2, pdd const& e2, fi_record& fi) { _last_function = __func__; - if (a2.is_one() && b2.is_val() && c.is_negative()) { + SASSERT(b2.is_val()); + if (a2.is_one() && c.is_negative()) { + // p > v + k auto& m = e2.manager(); rational const& mod_value = m.two_to_N(); - rational lo_val(mod(-b2.val() - 1, mod_value)); - auto lo = -e2 - 1; - rational hi_val(mod(lo_val + 1, mod_value)); - auto hi = -e2; + rational hi_val = (-b2).val(); + pdd hi = -e2; + rational lo_val = mod(hi_val - 1, mod_value); + pdd lo = p - e2; + fi.coeff = 1; + fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val); + return true; + } + if (a2.is_odd() && c.is_negative()) { + // p > a*v + k, a odd + auto& m = e2.manager(); + rational const& mod_value = m.two_to_N(); + rational a2_inv; + VERIFY(a2.mult_inverse(m.power_of_2(), a2_inv)); + rational lo_val = mod(a2_inv * (-1 - b2.val()), mod_value); + pdd lo = a2_inv * (-1 - e2); + rational hi_val = mod(lo_val + 1, mod_value); + pdd hi = lo + 1; fi.coeff = 1; fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val); - if (b2 != e2) - fi.side_cond.push_back(s.eq(b2, e2)); return true; } return false; diff --git a/src/math/polysat/forbidden_intervals.h b/src/math/polysat/forbidden_intervals.h index 9d5a6b393..4419a084a 100644 --- a/src/math/polysat/forbidden_intervals.h +++ b/src/math/polysat/forbidden_intervals.h @@ -73,12 +73,19 @@ namespace polysat { rational const & a2, pdd const& b2, pdd const& e2, fi_record& fi); - bool match_non_zero(signed_constraint const& c, + bool match_max(signed_constraint const& c, rational const & a1, pdd const& b1, pdd const& e1, + rational const & a2, pdd const& b2, pdd const& e2, + fi_record& fi); + + bool match_non_zero(signed_constraint const& c, + rational const& a1, pdd const& b1, pdd const& e1, + pdd const& q, fi_record& fi); bool match_non_max(signed_constraint const& c, - rational const & a2, pdd const& b2, pdd const& e2, + pdd const& p, + rational const& a2, pdd const& b2, pdd const& e2, fi_record& fi); bool get_interval_ule(signed_constraint const& c, pvar v, fi_record& fi); diff --git a/src/math/polysat/log.cpp b/src/math/polysat/log.cpp index 5e8642ca3..970ed5583 100644 --- a/src/math/polysat/log.cpp +++ b/src/math/polysat/log.cpp @@ -50,9 +50,6 @@ static LogLevel get_max_log_level(std::string const& fn, std::string const& pret (void)fn; (void)pretty_fn; - if (fn == "push_cjust") - return LogLevel::Verbose; - // if (fn == "pop_levels") // return LogLevel::Default; @@ -65,9 +62,11 @@ static LogLevel get_max_log_level(std::string const& fn, std::string const& pret } /// Filter log messages -bool polysat_should_log(LogLevel msg_level, std::string fn, std::string pretty_fn) { +bool polysat_should_log(unsigned verbose_lvl, LogLevel msg_level, std::string fn, std::string pretty_fn) { if (!g_log_enabled) return false; + if (get_verbosity_level() < verbose_lvl) + return false; LogLevel max_log_level = get_max_log_level(fn, pretty_fn); return msg_level <= max_log_level; } diff --git a/src/math/polysat/log.h b/src/math/polysat/log.h index 439fd3165..98e1e53db 100644 --- a/src/math/polysat/log.h +++ b/src/math/polysat/log.h @@ -56,12 +56,11 @@ enum class LogLevel : int { Heading2 = 2, Heading3 = 3, Default = 4, - Verbose = 5, }; /// Filter log messages bool -polysat_should_log(LogLevel msg_level, std::string fn, std::string pretty_fn); +polysat_should_log(unsigned verbose_lvl, LogLevel msg_level, std::string fn, std::string pretty_fn); std::pair polysat_log(LogLevel msg_level, std::string fn, std::string pretty_fn); @@ -70,31 +69,36 @@ polysat_log(LogLevel msg_level, std::string fn, std::string pretty_fn); #define __PRETTY_FUNCTION__ __FUNCSIG__ #endif -#define LOG_(lvl, x) \ - do { \ - if (polysat_should_log(lvl, __func__, __PRETTY_FUNCTION__)) { \ - auto pair = polysat_log(lvl, __func__, __PRETTY_FUNCTION__); \ - std::ostream& os = pair.first; \ - bool should_reset = pair.second; \ - os << x; \ - if (should_reset) \ - os << color_reset(); \ - os << std::endl; \ - } \ +#define LOG_(verbose_lvl, log_lvl, x) \ + do { \ + if (polysat_should_log(verbose_lvl, log_lvl, __func__, __PRETTY_FUNCTION__)) { \ + auto pair = polysat_log(log_lvl, __func__, __PRETTY_FUNCTION__); \ + std::ostream& os = pair.first; \ + bool should_reset = pair.second; \ + os << x; \ + if (should_reset) \ + os << color_reset(); \ + os << std::endl; \ + } \ } while (false) #define LOG_CONCAT_HELPER(a,b) a ## b #define LOG_CONCAT(a,b) LOG_CONCAT_HELPER(a,b) -#define LOG_INDENT(lvl, x) \ - LOG_(lvl, x); \ +#define LOG_INDENT(verbose_lvl, log_lvl, x) \ + LOG_(verbose_lvl, log_lvl, x); \ polysat_log_indent LOG_CONCAT(polysat_log_indent_obj_, __LINE__) (4); -#define LOG_H1(x) LOG_INDENT(LogLevel::Heading1, x) -#define LOG_H2(x) LOG_INDENT(LogLevel::Heading2, x) -#define LOG_H3(x) LOG_INDENT(LogLevel::Heading3, x) -#define LOG(x) LOG_(LogLevel::Default , x) -#define LOG_V(x) LOG_(LogLevel::Verbose , x) +#define LOG_H1(x) LOG_INDENT(0, LogLevel::Heading1, x) +#define LOG_H2(x) LOG_INDENT(0, LogLevel::Heading2, x) +#define LOG_H3(x) LOG_INDENT(0, LogLevel::Heading3, x) +#define LOG(x) LOG_(0, LogLevel::Default , x) + +#define LOG_H1_V(verbose_lvl, x) LOG_INDENT(verbose_lvl, LogLevel::Heading1, x) +#define LOG_H2_V(verbose_lvl, x) LOG_INDENT(verbose_lvl, LogLevel::Heading2, x) +#define LOG_H3_V(verbose_lvl, x) LOG_INDENT(verbose_lvl, LogLevel::Heading3, x) +#define LOG_V(verbose_lvl, x) LOG_(verbose_lvl, LogLevel::Default , x) + #define COND_LOG(c, x) if (c) LOG(x) #define LOGE(x) LOG(#x << " = " << (x)) @@ -110,18 +114,23 @@ inline void set_log_enabled(bool) {} inline bool get_log_enabled() { return false; } class scoped_set_log_enabled {}; -#define LOG_(lvl, x) \ - do { \ - /* do nothing */ \ +#define LOG_(vlvl, lvl, x) \ + do { \ + /* do nothing */ \ } while (false) -#define LOG_H1(x) LOG_(0, x) -#define LOG_H2(x) LOG_(0, x) -#define LOG_H3(x) LOG_(0, x) -#define LOG(x) LOG_(0, x) -#define LOG_V(x) LOG_(0, x) -#define COND_LOG(c, x) LOG_(c, x) -#define LOGE(x) LOG_(0, x) +#define LOG_H1(x) LOG_(0, 0, x) +#define LOG_H2(x) LOG_(0, 0, x) +#define LOG_H3(x) LOG_(0, 0, x) +#define LOG(x) LOG_(0, 0, x) + +#define LOG_H1_V(v, x) LOG_(v, 0, x) +#define LOG_H2_V(v, x) LOG_(v, 0, x) +#define LOG_H3_V(v, x) LOG_(v, 0, x) +#define LOG_V(v, x) LOG_(v, 0, x) + +#define COND_LOG(c, x) LOG_(0, c, x) +#define LOGE(x) LOG_(0, 0, x) #define IF_LOGGING(x) \ do { \ diff --git a/src/math/polysat/number.h b/src/math/polysat/number.h new file mode 100644 index 000000000..4a7af6024 --- /dev/null +++ b/src/math/polysat/number.h @@ -0,0 +1,30 @@ +/*++ +Copyright (c) 2021 Microsoft Corporation + +Module Name: + + polysat numbers + +Author: + + Nikolaj Bjorner (nbjorner) 2021-03-19 + Jakob Rath 2021-04-06 + +--*/ +#pragma once +#include "math/polysat/types.h" + +namespace polysat { + + inline unsigned get_parity(rational const& val, unsigned num_bits) { + if (val.is_zero()) + return num_bits; + return val.trailing_zeros(); + }; + + /** Return val with the lower k bits set to zero. */ + inline rational clear_lower_bits(rational const& val, unsigned k) { + return val - mod(val, rational::power_of_two(k)); + } + +} diff --git a/src/math/polysat/op_constraint.cpp b/src/math/polysat/op_constraint.cpp index a60eee3b0..0cf0c92ba 100644 --- a/src/math/polysat/op_constraint.cpp +++ b/src/math/polysat/op_constraint.cpp @@ -25,8 +25,8 @@ Additional possible functionality on constraints: namespace polysat { - op_constraint::op_constraint(constraint_manager& m, code c, pdd const& p, pdd const& q, pdd const& r) : - constraint(m, ckind_t::op_t), m_op(c), m_p(p), m_q(q), m_r(r) { + op_constraint::op_constraint(code c, pdd const& p, pdd const& q, pdd const& r) : + constraint(ckind_t::op_t), m_op(c), m_p(p), m_q(q), m_r(r) { m_vars.append(p.free_vars()); for (auto v : q.free_vars()) if (!m_vars.contains(v)) @@ -570,20 +570,26 @@ namespace polysat { } return true; } - - void op_constraint::add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { - auto p_coeff = s.subst(p()).get_univariate_coefficients(); - auto q_coeff = s.subst(q()).get_univariate_coefficients(); - auto r_coeff = s.subst(r()).get_univariate_coefficients(); + + void op_constraint::add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { + pdd pv = s.subst(p()); + if (!pv.is_univariate_in(v)) + return; + pdd qv = s.subst(q()); + if (!qv.is_univariate_in(v)) + return; + pdd rv = s.subst(r()); + if (!rv.is_univariate_in(v)) + return; switch (m_op) { case code::lshr_op: - us.add_lshr(p_coeff, q_coeff, r_coeff, !is_positive, dep); + us.add_lshr(pv.get_univariate_coefficients(), qv.get_univariate_coefficients(), rv.get_univariate_coefficients(), !is_positive, dep); break; case code::shl_op: - us.add_shl(p_coeff, q_coeff, r_coeff, !is_positive, dep); + us.add_shl(pv.get_univariate_coefficients(), qv.get_univariate_coefficients(), rv.get_univariate_coefficients(), !is_positive, dep); break; case code::and_op: - us.add_and(p_coeff, q_coeff, r_coeff, !is_positive, dep); + us.add_and(pv.get_univariate_coefficients(), qv.get_univariate_coefficients(), rv.get_univariate_coefficients(), !is_positive, dep); break; default: NOT_IMPLEMENTED_YET(); diff --git a/src/math/polysat/op_constraint.h b/src/math/polysat/op_constraint.h index 85166a562..8414436a6 100644 --- a/src/math/polysat/op_constraint.h +++ b/src/math/polysat/op_constraint.h @@ -35,7 +35,7 @@ namespace polysat { pdd m_q; // operand2 pdd m_r; // result - op_constraint(constraint_manager& m, code c, pdd const& p, pdd const& q, pdd const& r); + op_constraint(code c, pdd const& p, pdd const& q, pdd const& r); lbool eval(pdd const& p, pdd const& q, pdd const& r) const; clause_ref produce_lemma(solver& s, assignment const& a); @@ -73,7 +73,7 @@ namespace polysat { bool operator==(constraint const& other) const override; bool is_eq() const override { return false; } - void add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; + void add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; }; struct op_constraint_args { diff --git a/src/math/polysat/saturation.cpp b/src/math/polysat/saturation.cpp index 8d6d7229a..c9cf8bfd2 100644 --- a/src/math/polysat/saturation.cpp +++ b/src/math/polysat/saturation.cpp @@ -32,41 +32,153 @@ namespace polysat { saturation::saturation(solver& s) : s(s), m_lemma(s) {} + void saturation::log_lemma(pvar v, conflict& core) { + IF_VERBOSE(1, auto const& cl = core.lemmas().back(); + verbose_stream() << m_rule << " v" << v << " "; + for (auto lit : *cl) verbose_stream() << s.lit2cnstr(lit) << " "; + verbose_stream() << " " << *cl << "\n"); + } + void saturation::perform(pvar v, conflict& core) { + IF_VERBOSE(2, verbose_stream() << "v" << v << " " << core << "\n"); for (auto c : core) - if (perform(v, c, core)) { - IF_VERBOSE(0, auto const& cl = core.lemmas().back(); - verbose_stream() << m_rule << " v" << v << " "; - for (auto lit : *cl) verbose_stream() << s.lit2cnstr(lit) << " "; - verbose_stream() << "\n"); + if (perform(v, c, core)) return; - } } bool saturation::perform(pvar v, signed_constraint const& c, conflict& core) { - if (!c->is_ule()) - return false; if (c.is_currently_true(s)) return false; - auto i = inequality::from_ule(c); + + if (c->is_ule()) { + auto i = inequality::from_ule(c); + return try_inequality(v, i, core); + } + if (c->is_umul_ovfl()) + return try_umul_ovfl(v, c, core); + return false; + } + + bool saturation::try_inequality(pvar v, inequality const& i, conflict& core) { + bool prop = false; if (try_mul_bounds(v, core, i)) - return true; + prop = true; if (try_parity(v, core, i)) - return true; + prop = true; if (try_parity_diseq(v, core, i)) - return true; - if (try_factor_equality(v, core, i)) - return true; + prop = true; + if (try_transitivity(v, core, i)) + prop = true; + if (try_factor_equality2(v, core, i)) + prop = true; + if (try_infer_equality(v, core, i)) + prop = true; + if (try_add_overflow_bound(v, core, i)) + prop = true; + if (try_add_mul_bound(v, core, i)) + prop = true; + if (try_mul_eq_bound(v, core, i)) + prop = true; if (try_ugt_x(v, core, i)) - return true; + prop = true; if (try_ugt_y(v, core, i)) - return true; + prop = true; if (try_ugt_z(v, core, i)) - return true; - if (try_y_l_ax_and_x_l_z(v, core, i)) - return true; + prop = true; + if (try_y_l_ax_and_x_l_z(v, core, i)) + prop = true; if (false && try_tangent(v, core, i)) - return true; + prop = true; + return prop; + } + + bool saturation::try_umul_ovfl(pvar v, signed_constraint const& c, conflict& core) { +#if 1 + return false; +#else + SASSERT(c->is_umul_ovfl()); + bool prop = false; + if (c.is_positive()) { + prop = try_umul_ovfl_bounds(v, c, core); + } + else { + prop = try_umul_noovfl_bounds(v, c, core); + if (false && try_umul_noovfl_lo(v, c, core)) + prop = true; + } + return prop; +#endif + } + + bool saturation::try_umul_noovfl_lo(pvar v, signed_constraint const& c, conflict& core) { + set_rule("[x] ~ovfl(x, y) => y = 0 or x <= x * y"); + SASSERT(c->is_umul_ovfl()); + if (!c.is_negative()) + return false; + auto const& ovfl = c->to_umul_ovfl(); + auto V = s.var(v); + auto p = ovfl.p(), q = ovfl.q(); + // TODO could relax condition to be that V occurs in p + if (q == V) + std::swap(p, q); + signed_constraint q_eq_0; + if (p == V && is_forced_diseq(q, 0, q_eq_0)) { + // ~ovfl(V,q) => q = 0 or V <= V*q + m_lemma.reset(); + m_lemma.insert_eval(q_eq_0); + if (propagate(v, core, c, s.ule(p, p * q))) + return true; + } + return false; + } + + /** + * ~ovfl(p, q) & p >= k => q < 2^N/k + * TODO: subusmed by narrowing inferences? + */ + bool saturation::try_umul_noovfl_bounds(pvar x, signed_constraint const& c, conflict& core) { + set_rule("[x] ~ovfl(x, q) & x >= k => q <= (2^N-1)/k"); + SASSERT(c->is_umul_ovfl()); + SASSERT(c.is_negative()); + auto const& ovfl = c->to_umul_ovfl(); + auto p = ovfl.p(), q = ovfl.q(); + auto X = s.var(x); + auto& m = p.manager(); + rational p_val, q_val; + if (q == X) + std::swap(p, q); + if (p == X) { + vector x_ge_bound; + if (!s.try_eval(q, q_val)) + return false; + if (!has_lower_bound(x, core, p_val, x_ge_bound)) + return false; + if (p_val * q_val <= m.max_value()) + return false; + m_lemma.reset(); + m_lemma.insert_eval(~s.uge(X, p_val)); + signed_constraint conseq = s.ule(q, floor(m.max_value()/p_val)); + return propagate(x, core, c, conseq); + } + if (!s.try_eval(p, p_val) || !s.try_eval(q, q_val)) + return false; + if (p_val * q_val <= m.max_value()) + return false; + m_lemma.reset(); + m_lemma.insert_eval(~s.uge(q, q_val)); + signed_constraint conseq = s.ule(p, floor(m.max_value()/q_val)); + return propagate(x, core, c, conseq); + } + + /** + * ovfl(p, q) & p <= k => q > 2^N/k + */ + bool saturation::try_umul_ovfl_bounds(pvar v, signed_constraint const& c, conflict& core) { + SASSERT(c->is_umul_ovfl()); + SASSERT(c.is_positive()); + auto const& ovfl = c->to_umul_ovfl(); + auto p = ovfl.p(), q = ovfl.q(); + rational p_val, q_val; return false; } @@ -77,7 +189,11 @@ namespace polysat { return s.ule(lhs, rhs); } - bool saturation::propagate(conflict& core, inequality const& crit, signed_constraint c) { + bool saturation::propagate(pvar v, conflict& core, inequality const& crit, signed_constraint c) { + return propagate(v, core, crit.as_signed_constraint(), c); + } + + bool saturation::propagate(pvar v, conflict& core, signed_constraint const& crit, signed_constraint c) { if (is_forced_true(c)) return false; @@ -95,25 +211,21 @@ namespace polysat { // The current assumptions on how conflict lemmas are used do not accomodate propagation it seems. // - m_lemma.insert(~crit.as_signed_constraint()); - - IF_VERBOSE(10, verbose_stream() << "propagate " << m_rule << " "; - for (auto lit : m_lemma) verbose_stream() << s.lit2cnstr(lit) << " "; - verbose_stream() << c << "\n"; - ); + m_lemma.insert(~crit); SASSERT(all_of(m_lemma, [this](sat::literal lit) { return is_forced_false(s.lit2cnstr(lit)); })); m_lemma.insert(c); core.add_lemma(m_rule, m_lemma.build()); + log_lemma(v, core); return true; } - bool saturation::add_conflict(conflict& core, inequality const& crit1, signed_constraint c) { - return add_conflict(core, crit1, crit1, c); + bool saturation::add_conflict(pvar v, conflict& core, inequality const& crit1, signed_constraint c) { + return add_conflict(v, core, crit1, crit1, c); } - bool saturation::add_conflict(conflict& core, inequality const& _crit1, inequality const& _crit2, signed_constraint const c) { + bool saturation::add_conflict(pvar v, conflict& core, inequality const& _crit1, inequality const& _crit2, signed_constraint const c) { auto crit1 = _crit1.as_signed_constraint(); auto crit2 = _crit2.as_signed_constraint(); m_lemma.insert(~crit1); @@ -135,6 +247,7 @@ namespace polysat { m_lemma.insert_eval(c); core.add_lemma(m_rule, m_lemma.build()); + log_lemma(v, core); return true; } @@ -192,6 +305,14 @@ namespace polysat { return is_g_v(x, i); } + /* + * Match [x] y <= x + */ + bool saturation::is_Y_l_x(pvar x, inequality const& i, pdd& y) { + y = i.lhs(); + return is_l_v(x, i); + } + /* * Match [x] y <= a*x */ @@ -221,8 +342,7 @@ namespace polysat { */ bool saturation::is_AxB_l_Y(pvar x, inequality const& i, pdd& a, pdd& b, pdd& y) { y = i.rhs(); - pdd aa = a, bb = b; - return i.lhs().degree(x) == 1 && (i.lhs().factor(x, 1, aa, bb), aa == a && bb == b); + return i.lhs().degree(x) == 1 && (i.lhs().factor(x, 1, a, b), true); } bool saturation::verify_AxB_l_Y(pvar x, inequality const& i, pdd const& a, pdd const& b, pdd const& y) { @@ -232,8 +352,7 @@ namespace polysat { bool saturation::is_Y_l_AxB(pvar x, inequality const& i, pdd& y, pdd& a, pdd& b) { y = i.lhs(); - pdd aa = a, bb = b; - return i.rhs().degree(x) == 1 && (i.rhs().factor(x, 1, aa, bb), aa == a && bb == b); + return i.rhs().degree(x) == 1 && (i.rhs().factor(x, 1, a, b), true); } bool saturation::verify_Y_l_AxB(pvar x, inequality const& i, pdd const& y, pdd const& a, pdd& b) { @@ -386,7 +505,7 @@ namespace polysat { m_lemma.insert_eval(~non_ovfl); if (!xy_l_xz.is_strict()) m_lemma.insert_eval(s.eq(x)); - return add_conflict(core, xy_l_xz, ineq(xy_l_xz.is_strict(), y, z)); + return add_conflict(v, core, xy_l_xz, ineq(xy_l_xz.is_strict(), y, z)); } /** @@ -427,7 +546,7 @@ namespace polysat { pdd const& z_prime = l_y.lhs(); m_lemma.reset(); m_lemma.insert_eval(~non_ovfl); - return add_conflict(core, l_y, yx_l_zx, ineq(yx_l_zx.is_strict() || l_y.is_strict(), z_prime * x, z * x)); + return add_conflict(v, core, l_y, yx_l_zx, ineq(yx_l_zx.is_strict() || l_y.is_strict(), z_prime * x, z * x)); } /** @@ -467,7 +586,7 @@ namespace polysat { return false; m_lemma.reset(); m_lemma.insert_eval(~non_ovfl); - return add_conflict(core, yx_l_zx, z_l_y, ineq(z_l_y.is_strict() || yx_l_zx.is_strict(), y * x, y_prime * x)); + return add_conflict(z, core, yx_l_zx, z_l_y, ineq(z_l_y.is_strict() || yx_l_zx.is_strict(), y * x, y_prime * x)); } /** @@ -508,7 +627,7 @@ namespace polysat { return false; m_lemma.reset(); m_lemma.insert_eval(~non_ovfl); - return add_conflict(core, y_l_ax, x_l_z, ineq(x_l_z.is_strict() || y_l_ax.is_strict(), y, a * z)); + return add_conflict(x, core, y_l_ax, x_l_z, ineq(x_l_z.is_strict() || y_l_ax.is_strict(), y, a * z)); } /** @@ -546,7 +665,7 @@ namespace polysat { m_lemma.insert_eval(~s.eq(y)); m_lemma.insert_eval(x_eq_0); m_lemma.insert_eval(a_eq_0); - return propagate(core, axb_l_y, c); + return propagate(x, core, axb_l_y, c); }; auto prop2 = [&](signed_constraint ante, signed_constraint c) { @@ -556,7 +675,7 @@ namespace polysat { m_lemma.insert_eval(x_eq_0); m_lemma.insert_eval(a_eq_0); m_lemma.insert_eval(~ante); - return propagate(core, axb_l_y, c); + return propagate(x, core, axb_l_y, c); }; pdd minus_a = -a; @@ -612,6 +731,71 @@ namespace polysat { } + // bench 5 + // fairly ad-hoc rule derived from bench 5. + // The clause could also be added whenever narrowing the literal 2^k*x = 2^k*y + // It can be expected to be relatively common because these equalities come from + // bit-masking. + // + // A bigger hammer for detecting such propagations may be through LIA or a variant + // + // a*x - a*y + b*z = 0 0 <= x < b/a, 0 <= y < b/a => z = 0 + // and then => x = y + // + // the general lemma is that the linear term a*p = 0 is such that a*p does not overflow + // and therefore p = 0 + // + // TBD: encode the general lemma instead of this special case. + // + bool saturation::try_mul_eq_bound(pvar x, conflict& core, inequality const& axb_l_y) { + set_rule("[x] 2^k*x = 2^k*y & x < 2^N-k => y = x or y >= 2^{N-k}"); + auto& m = s.var2pdd(x); + unsigned N = m.power_of_2(); + pdd y = m.zero(); + pdd a = y, b = y, a2 = y; + pdd X = s.var(x); + if (!is_AxB_eq_0(x, axb_l_y, a, b, y)) + return false; + if (!a.is_val()) + return false; + if (!a.val().is_power_of_two()) + return false; + unsigned k = a.val().trailing_zeros(); + if (k == 0) + return false; + b = -b; + if (b.leading_coefficient() != a.val()) + return false; + for (auto c : core) { + if (!c->is_ule()) + continue; + auto i = inequality::from_ule(c); + if (!is_x_l_Y(x, i, a2)) + continue; + if (!a2.is_val()) + continue; + // x*2^k = b, x <= a2 < 2^{N-k} + rational bound = rational::power_of_two(N - k); + if (i.is_strict() && a2.val() >= bound) + continue; + if (!i.is_strict() && a2.val() > bound) + continue; + pdd Y = b.div(b.leading_coefficient()); + rational Y_val; + if (!s.try_eval(Y, Y_val) || Y_val >= bound) + continue; + signed_constraint le = s.ule(Y, bound - 1); + m_lemma.reset(); + m_lemma.insert_eval(~le); + m_lemma.insert_eval(~s.eq(y)); + m_lemma.insert(~c); + if (propagate(x, core, axb_l_y, s.eq(X, Y))) + return true; + } + + return false; + } + /* * x*y = 1 & ~ovfl(x,y) => x = 1 * x*y = -1 & ~ovfl(-x,y) => -x = 1 @@ -634,9 +818,9 @@ namespace polysat { m_lemma.insert_eval(~s.eq(b, rational(-1))); m_lemma.insert_eval(~s.eq(y)); m_lemma.insert_eval(~non_ovfl); - if (propagate(core, axb_l_y, s.eq(X, 1))) + if (propagate(x, core, axb_l_y, s.eq(X, 1))) return true; - if (propagate(core, axb_l_y, s.eq(a, 1))) + if (propagate(x, core, axb_l_y, s.eq(a, 1))) return true; return false; } @@ -651,26 +835,79 @@ namespace polysat { * * odd(x) & even(y) => x + y != 0 * + * Special case rule: a*x + y = 0 => (odd(b) <=> odd(a) & odd(x)) + * * General rule: * * a*x + y = 0 => min(K, parity(a) + parity(x)) = parity(y) - * - * Currently implemented special case: a*x + y = 0 => (odd(b) <=> odd(a) & odd(x)) - * - * general rule can be obtained by adding an - * 'is_forced_parity(x, p, x_has_parity_p)' * - * Should we also check 'is_forced_parity(a*x, p, ax_has_parity_p)' - * if a*x has a parity but not a, x? + * using inequalities: + * + * parity(x) <= i, parity(a) <= j => parity(b) <= i + j + * parity(x) >= i, parity(a) >= j => parity(b) >= i + j + * parity(x) <= i, parity(b) >= j => parity(a) >= j - i + * parity(x) >= i, parity(b) <= j => parity(a) <= j - i + * symmetric rules for swapping x, a + * + * min_parity(x) = number of trailing bits of x if x is a value + * min_parity(x) = k if 2^{N-k}*x == 0 is forced for max k + * min_parity(x1*x2) = min_parity(x1) + min_parity(x2) + * min_parity(x) = 0, otherwise + * + * max_parity(x) = number of trailing bits of x + * max_parity(x) = k if 2^{N-k-1}*x != 0 for min k + * max_parity(x1*x2) = max_parity(x1) + max_parity(x2) + * max_parity(x) = N, otherwise * */ + + unsigned saturation::min_parity(pdd const& p) { + rational val; + auto& m = p.manager(); + unsigned N = m.power_of_2(); + if (s.try_eval(p, val)) + return val == 0 ? N : val.trailing_zeros(); + +#if 0 + // TBD: factor p + auto coeff = p.leading_coefficient(); + unsigned offset = coeff.trailing_zeros(); + verbose_stream() << "COEFF " << coeff << "\n"; +#endif +#if 0 + unsigned j = 0; + while (j < N && is_forced_true(s.parity(p, j + 1))) + ++j; + return j; +#else + for (unsigned j = N; j > 0; --j) + if (is_forced_true(s.parity(p, j))) + return j; + return 0; +#endif + } + + unsigned saturation::max_parity(pdd const& p) { + auto& m = p.manager(); + unsigned N = m.power_of_2(); + rational val; + if (s.try_eval(p, val)) + return val == 0 ? N : val.trailing_zeros(); + + // TBD: factor p + + for (unsigned j = 0; j < N; ++j) + if (is_forced_true(s.parity_at_most(p, j))) + return j; + return N; + } + bool saturation::try_parity(pvar x, conflict& core, inequality const& axb_l_y) { set_rule("[x] a*x + b = 0 => (odd(a) & odd(x) <=> odd(b))"); auto& m = s.var2pdd(x); unsigned N = m.power_of_2(); pdd y = m.zero(); - pdd a = m.zero(); - pdd b = m.zero(); + pdd a = y, b = y; pdd X = s.var(x); if (!is_AxB_eq_0(x, axb_l_y, a, b, y)) return false; @@ -680,93 +917,76 @@ namespace polysat { return false; if (a.is_one()) return false; - signed_constraint b_is_odd = s.odd(b); - signed_constraint a_is_odd = s.odd(a); - signed_constraint x_is_odd = s.odd(X); + if (a.is_val() && b.is_zero()) + return false; auto propagate1 = [&](signed_constraint premise, signed_constraint conseq) { + if (is_forced_false(premise)) + return false; + IF_VERBOSE(1, verbose_stream() << "propagate " << axb_l_y << " " << premise << " => " << conseq << "\n"); m_lemma.reset(); m_lemma.insert_eval(~s.eq(y)); m_lemma.insert_eval(~premise); - return propagate(core, axb_l_y, conseq); + return propagate(x, core, axb_l_y, conseq); }; auto propagate2 = [&](signed_constraint premise1, signed_constraint premise2, signed_constraint conseq) { + if (is_forced_false(premise1)) + return false; + if (is_forced_false(premise2)) + return false; + IF_VERBOSE(1, verbose_stream() << "propagate " << axb_l_y << " " << premise1 << " " << premise2 << " => " << conseq << "\n"); m_lemma.reset(); m_lemma.insert_eval(~s.eq(y)); m_lemma.insert_eval(~premise1); m_lemma.insert_eval(~premise2); - return propagate(core, axb_l_y, conseq); + return propagate(x, core, axb_l_y, conseq); }; -#if 0 - LOG_H1("try_parity: " << X << " on: " << lit_pp(s, axb_l_y.as_signed_constraint())); - LOG("y: " << y << " a: " << a << " b: " << b); - LOG("b_is_odd: " << lit_pp(s, b_is_odd)); - LOG("a_is_odd: " << lit_pp(s, a_is_odd)); - LOG("x_is_odd: " << lit_pp(s, x_is_odd)); -#endif - if (a_is_odd.is_currently_true(s) && - x_is_odd.is_currently_true(s) && - propagate2(a_is_odd, x_is_odd, b_is_odd)) + + unsigned min_x = min_parity(X), max_x = max_parity(X); + unsigned min_b = min_parity(b), max_b = max_parity(b); + unsigned min_a = min_parity(a), max_a = max_parity(a); + SASSERT(min_x <= max_x && max_x <= N); + SASSERT(min_a <= max_a && max_a <= N); + SASSERT(min_b <= max_b && max_b <= N); + + IF_VERBOSE(2, + verbose_stream() << "try parity v" << x << " " << axb_l_y << "\n"; + verbose_stream() << X << " " << min_x << " " << max_x << "\n"; + verbose_stream() << a << " " << min_a << " " << max_a << "\n"; + verbose_stream() << b << " " << min_b << " " << max_b << "\n"); + + if (min_x >= N || min_a >= N) + return false; + + auto at_most = [&](pdd const& p, unsigned k) { + VERIFY(k < N); + return s.parity_at_most(p, k); + }; + + auto at_least = [&](pdd const& p, unsigned k) { + VERIFY(k != 0); + return s.parity(p, k); + }; + + + if (!b.is_val() && max_b > max_a + max_x && propagate2(at_most(a, max_a), at_most(X, max_x), at_most(b, max_x + max_a))) + return true; + if (!b.is_val() && min_x > min_b && propagate1(at_least(X, min_x), at_least(b, min_x))) + return true; + if (!b.is_val() && min_a > min_b && propagate1(at_least(a, min_a), at_least(b, min_a))) + return true; + if (!b.is_val() && min_x > 0 && min_a > 0 && min_x + min_a > min_b && propagate2(at_least(a, min_a), at_least(X, min_x), at_least(b, min_a + min_x))) + return true; + if (!a.is_val() && max_x <= min_b && min_a < min_b - max_x && propagate2(at_most(X, max_x), at_least(b, min_b), at_least(a, min_b - max_x))) + return true; + if (max_a <= min_b && min_x < min_b - max_a && propagate2(at_most(a, max_a), at_least(b, min_b), at_least(X, min_b - max_a))) + return true; + if (max_b < N && !a.is_val() && min_x > 0 && min_x <= max_b && max_a > max_b - min_x && propagate2(at_least(X, min_x), at_most(b, max_b), at_most(a, max_b - min_x))) + return true; + if (max_b < N && min_a > 0 && min_a <= max_b && max_x > max_b - min_a && propagate2(at_least(a, min_a), at_most(b, max_b), at_most(X, max_b - min_a))) return true; - - if (b_is_odd.is_currently_true(s)) { - if (propagate1(b_is_odd, a_is_odd)) - return true; - if (propagate1(b_is_odd, x_is_odd)) - return true; - } - - // a is divisibly by 4, - // max divisor of x is k - // -> b has parity k + 4 - unsigned a_parity = a_is_odd.is_currently_false(s) ? 1 : 0; - unsigned x_parity = x_is_odd.is_currently_false(s) ? 1 : 0; - - if ((a_parity > 0 || x_parity > 0) && !is_forced_eq(a, 0) && !is_forced_eq(X, 0)) { - while (a_parity < N && s.parity(a, a_parity+1).is_currently_true(s)) - ++a_parity; - while (x_parity < N && s.parity(X, x_parity+1).is_currently_true(s)) - ++x_parity; - unsigned b_parity = std::min(N, a_parity + x_parity); - if (a_parity > 0 && x_parity > 0 && propagate2(s.parity(a, a_parity), s.parity(X, x_parity), s.parity(b, b_parity))) - return true; - if (a_parity > 0 && x_parity == 0 && propagate1(s.parity(a, a_parity), s.parity(b, b_parity))) - return true; - if (a_parity == 0 && x_parity > 0 && propagate1(s.parity(X, x_parity), s.parity(b, b_parity))) - return true; - } - - // - // if b has at most b_parity, then a*x has at most b_parity - // - else if (!is_forced_eq(b, 0)) { - unsigned b_parity = 1; - bool found = false; - for (; b_parity < N; ++b_parity) { - if (s.parity(b, b_parity).is_currently_false(s)) { - found = true; - break; - } - } - if (found) { - if (propagate1(~s.parity(b, b_parity), ~s.parity(a, b_parity))) - return true; - if (propagate1(~s.parity(b, b_parity), ~s.parity(a, b_parity))) - return true; - - for (unsigned i = 1; i < N; ++i) { - if (s.parity(a, i).is_currently_true(s) && - propagate2(~s.parity(b, b_parity), s.parity(a, i), ~s.parity(X, b_parity - i))) - return true; - - if (s.parity(X, i).is_currently_true(s) && - propagate2(~s.parity(b, b_parity), s.parity(X, i), ~s.parity(a, b_parity - i))) - return true; - } - } - } return false; } @@ -780,7 +1000,8 @@ namespace polysat { auto& m = s.var2pdd(x); unsigned N = m.power_of_2(); pdd y = m.zero(); - pdd a = y, b = y, X = y; + pdd a = y, b = y; + pdd X = s.var(x); if (!is_AxB_diseq_0(x, axb_l_y, a, b, y)) return false; if (!is_forced_eq(b, 0)) @@ -793,7 +1014,10 @@ namespace polysat { m_lemma.reset(); m_lemma.insert_eval(~s.eq(y)); m_lemma.insert_eval(~s.eq(b)); - return propagate(core, axb_l_y, ~s.parity(X, N - k)); + if (propagate(x, core, axb_l_y, ~s.parity(X, N - k))) + return true; + + return false; } /** @@ -818,16 +1042,78 @@ namespace polysat { m_lemma.insert_eval(s.eq(y)); m_lemma.insert_eval(~s.eq(b)); m_lemma.insert_eval(a_eq_0); - if (propagate(core, axb_l_y, s.even(X))) + if (propagate(x, core, axb_l_y, s.even(X))) return true; if (!is_forced_diseq(X, 0, x_eq_0)) return false; m_lemma.insert_eval(x_eq_0); - if (propagate(core, axb_l_y, s.even(a))) + if (propagate(x, core, axb_l_y, s.even(a))) return true; return false; } + /** + * TODO If both inequalities are strict, then the implied inequality has a gap of 2 + * a < b, b < c => a + 1 < c & a + 1 != 0 + */ + bool saturation::try_transitivity(pvar x, conflict& core, inequality const& a_l_b) { + set_rule("[x] q < x & x <= p => q < p"); + auto& m = s.var2pdd(x); + pdd p = m.zero(); + pdd a = p, b = p, q = p; + // x <= p + if (!is_Ax_l_Y(x, a_l_b, a, p)) + return false; + if (!is_forced_eq(a, 1)) + return false; + for (auto c : core) { + if (!c->is_ule()) + continue; + auto i = inequality::from_ule(c); + if (c == a_l_b.as_signed_constraint()) + continue; + if (!is_Y_l_Ax(x, i, b, q)) + continue; + if (!is_forced_eq(b, 1)) + continue; + m_lemma.reset(); + m_lemma.insert_eval(~s.eq(a, 1)); + m_lemma.insert_eval(~s.eq(b, 1)); + m_lemma.insert(~c); + auto ineq = i.is_strict() || a_l_b.is_strict() ? (p.is_val() ? s.ule(q, p - 1) : s.ult(q, p)) : s.ule(q, p); + if (propagate(x, core, a_l_b, ineq)) + return true; + } + + return false; + } + + /** + * p <= q, q <= p => p - q = 0 + */ + bool saturation::try_infer_equality(pvar x, conflict& core, inequality const& a_l_b) { + set_rule("[x] p <= q, q <= p => p - q = 0"); + if (a_l_b.is_strict()) + return false; + if (a_l_b.lhs().degree(x) == 0 && a_l_b.rhs().degree(x) == 0) + return false; + for (auto c : core) { + if (!c->is_ule()) + continue; + auto i = inequality::from_ule(c); + if (i.lhs() == a_l_b.rhs() && i.rhs() == a_l_b.lhs() && !i.is_strict()) { + m_lemma.reset(); + m_lemma.insert(~c); + if (propagate(x, core, a_l_b, s.eq(i.lhs() - i.rhs()))) { + verbose_stream() << "infer equality " << s.eq(i.lhs() - i.rhs()) << "\n"; + return true; + } + } + } + return false; + } + + /** * [x] ax + p <= q, ax + r = 0 => -r + p <= q * [x] p <= ax + q, ax + r = 0 => p <= -r + q @@ -836,7 +1122,8 @@ namespace polysat { * [x] p <= abx + q, ax + r = 0 => p <= -rb + q */ - bool saturation::try_factor_equality(pvar x, conflict& core, inequality const& a_l_b) { + bool saturation::try_factor_equality1(pvar x, conflict& core, inequality const& a_l_b) { + set_rule("[x] ax + b = 0 & C[x] => C[-inv(a)*b]"); auto& m = s.var2pdd(x); unsigned N = m.power_of_2(); pdd y1 = m.zero(); @@ -846,14 +1133,10 @@ namespace polysat { bool is_axb_l_y = is_AxB_l_Y(x, a_l_b, a1, b1, y1); bool is_y_l_axb = is_Y_l_AxB(x, a_l_b, y2, a2, b2); - if (!a_l_b.is_strict() && a_l_b.rhs().is_zero()) - return false; - if (!is_axb_l_y && !is_y_l_axb) return false; - if (a1.is_val() && a2.is_val()) - return false; + bool factored = false; for (auto c : core) { if (!c->is_ule()) @@ -865,32 +1148,319 @@ namespace polysat { continue; if (c == a_l_b.as_signed_constraint()) continue; - pdd lhs = i.lhs(); - pdd rhs = i.rhs(); + pdd lhs = a_l_b.lhs(); + pdd rhs = a_l_b.rhs(); bool change = false; if (is_axb_l_y && a1 == a3) { change = true; - lhs = b3 - b1; + lhs = b1 - b3; + } + else if (is_axb_l_y && a1 == -a3) { + change = true; + lhs = b1 + b3; + } + else if (is_axb_l_y && a3.is_val() && a3.val().is_odd()) { + // a3*x + b3 == 0 + // a3 is odd => x = inverse(a3)*-b3 + change = true; + rational a3_inv; + VERIFY(a3.val().mult_inverse(m.power_of_2(), a3_inv)); + lhs = b1 - a1*(b3 * a3_inv); } if (is_y_l_axb && a2 == a3) { change = true; - rhs = b3 - b2; + rhs = b2 - b3; + } + else if (is_y_l_axb && a2 == -a3) { + change = true; + rhs = b2 + b3; + } + else if (is_y_l_axb && a3.is_val() && a3.val().is_odd()) { + change = true; + rational a3_inv; + VERIFY(a3.val().mult_inverse(m.power_of_2(), a3_inv)); + rhs = b2 - a2*(b3 * a3_inv); } if (!change) { - IF_VERBOSE(0, verbose_stream() << "missed factor equality " << c << " " << a_l_b << "\n"); + IF_VERBOSE(2, verbose_stream() << "missed factor equality? " << c << " " << a_l_b << "\n"); continue; } - signed_constraint conseq = i.is_strict() ? s.ult(lhs, rhs) : s.ule(lhs, rhs); + signed_constraint conseq = a_l_b.is_strict() ? s.ult(lhs, rhs) : s.ule(lhs, rhs); m_lemma.reset(); - m_lemma.insert(~s.eq(y3)); + m_lemma.insert_eval(~s.eq(y3)); m_lemma.insert(~c); - IF_VERBOSE(0, verbose_stream() << "factor equality " << a_l_b << "\n"); - if (propagate(core, a_l_b, conseq)) + if (propagate(x, core, a_l_b, conseq)) + factored = true; + } + return factored; + } + + bool saturation::try_factor_equality2(pvar x, conflict& core, inequality const& a_l_b) { + set_rule("[x] ax + b = 0 & C[x] => C[-inv(a)*b]"); + auto& m = s.var2pdd(x); + pdd y = m.zero(); + pdd a = y, b = y, a1 = y, b1 = y; + if (!is_AxB_eq_0(x, a_l_b, a, b, y)) + return false; + bool is_invertible = a.is_val() && a.val().is_odd(); + if (is_invertible) { + rational a_inv; + VERIFY(a.val().mult_inverse(m.power_of_2(), a_inv)); + b = -b*a_inv; + } + bool change = false; + bool prop = false; + auto replace = [&](pdd p) { + unsigned p_degree = p.degree(x); + if (p_degree == 0) + return p; + if (is_invertible) { + change = true; + return p.subst_pdd(x, b); + } + if (p_degree == 1) { + p.factor(x, 1, a1, b1); + if (a1 == a) { + change = true; + return b1 - b; + } + if (a1 == -a) { + change = true; + return b1 + b; + } + } + return p; + }; + + for (auto c : core) { + change = false; + if (c == a_l_b.as_signed_constraint()) + continue; + if (c->is_ule()) { + auto const& ule = c->to_ule(); + auto p = replace(ule.lhs()); + auto q = replace(ule.rhs()); + if (!change) + continue; + m_lemma.reset(); + m_lemma.insert(~c); + m_lemma.insert_eval(~s.eq(y)); + if (propagate(x, core, a_l_b, c.is_positive() ? s.ule(p, q) : ~s.ule(p, q))) + prop = true; + } + else if (c->is_umul_ovfl()) { + auto const& ovf = c->to_umul_ovfl(); + auto p = replace(ovf.p()); + auto q = replace(ovf.q()); + if (!change) + continue; + m_lemma.reset(); + m_lemma.insert(~c); + m_lemma.insert_eval(~s.eq(y)); + if (propagate(x, core, a_l_b, c.is_positive() ? s.umul_ovfl(p, q) : ~s.umul_ovfl(p, q))) + prop = true; + } + } + return prop; + } + + + /** + * x >= x + y & x <= n => y = 0 or y >= N - n + * x > x + y & x <= n => y >= N - n + * -x <= -x - y & x <= n => y = 0 or y >= N - n + * -x < -x - y & x <= n => y >= N - n + */ + bool saturation::try_add_overflow_bound(pvar x, conflict& core, inequality const& axb_l_y) { + set_rule("[x] x >= x + y & x <= n => y = 0 or y >= 2^N - n"); + signed_constraint y_eq_0; + vector x_ge_bound; + auto& m = s.var2pdd(x); + pdd y = m.zero(); + if (!is_add_overflow(x, axb_l_y, y)) + return false; + if (!axb_l_y.is_strict() && !is_forced_diseq(y, 0, y_eq_0)) + return false; + rational bound; + if (!has_upper_bound(x, core, bound, x_ge_bound)) + return false; + SASSERT(bound != 0); + m_lemma.reset(); + if (!axb_l_y.is_strict()) + m_lemma.insert_eval(y_eq_0); + for (auto c : x_ge_bound) + m_lemma.insert_eval(~c); + return propagate(x, core, axb_l_y, s.uge(y, m.two_to_N() - bound)); + } + + /** + * Match one of the patterns: + * x >= x + y + * x > x + y + * -x <= -x - y + * -x < -x - y + */ + bool saturation::is_add_overflow(pvar x, inequality const& i, pdd& y) { + auto& m = s.var2pdd(x); + pdd X = s.var(x); + pdd a = X; + if (i.lhs().degree(x) != 1 || i.rhs().degree(x) != 1) + return false; + if (i.rhs() == X) { + i.lhs().factor(x, 1, a, y); + if (a.is_one()) return true; } + if (i.lhs() == -X) { + i.rhs().factor(x, 1, a, y); + if ((-a).is_one()) { + y = -y; + return true; + } + } return false; } + + bool saturation::has_upper_bound(pvar x, conflict& core, rational& bound, vector& x_le_bound) { + return s.m_viable.has_upper_bound(x, bound, x_le_bound); + } + + bool saturation::has_lower_bound(pvar x, conflict& core, rational& bound, vector& x_ge_bound) { + return s.m_viable.has_lower_bound(x, bound, x_ge_bound); + } + + /* + * Bounds propagation for base case ax <= y + * where + * & y <= u_y + * & 0 < x <= u_x + * + * Special case for interval including -1 + * Compute max n, such that a \not\in [-n,0[ is implied, then + * + * => a < -n + * Claim n = floor(- u_y / u_x), + * - provided n != 0 + * + * Justification: for a1 in [1,n[ : + * 0 < a1*x <= floor (- u_y / u_x) * u_x + * <= - u_y + * and therefore for a1 in [-n-1:0[ : + * u_y < a1*x < 0 + * + * + * Bounds case for additive case ax - b <= y + * where + * & y <= u_y + * & b <= u_b + * & u_y + u_b does not overflow (implies y + b >= y) + * => ax - b + b <= y + b + * => ax <= y + b + * => ax <= u_y + u_b + * + * Base case for additive case ax + b <= y + * where + * & y <= u_y + * & b >= l_b + * & ax + b >= b + * + * => ax + b - b <= y - b + * => ax <= y - b + * => ax <= u_y - l_b + * + * TODO - deal with side condition ax + b >= b? + * It requires that ax + b does not overflow + * If the literal is already assigned, we are fine, otherwise? + * + * + * Example (bench25) + * -1*v85*v33 + v81 <= 2^128-2 + * v33 <= 2^128-1 + * v81 := -1 + * v85 := 12 + * + * Example (bench25) + * -1489: v25 > -1*v85*v25 + v81 + * 2397: v85 + 1 <= 328022915686448145675838484443875093068753497636375535522542730900603766685 + * -1195: v85 + 1 > 2^128+1 + * v25 := 353 + * v81 := -1 + * + * -1*v85*v25 + v81 < v25 + * a -1*v25 := -315 b v81 := -1 y v25 := 315 + * & v25 <= 315 + * & -v81 <= 1 + * + * The example illustrates that fixing y_val produces a weaker bound. + * The result should be a forbidden interval around v25 based on bounds for + * v85 and v81. + * + * The example also illustrates that presumably just a combination of forbidden intervals for v85 used in the conflict + * are sufficient for narrowing the bounds on v81. Querying for upper/lower bounds of v85 doesn't produce the strongest + * assumption. 2397 and -1195 come from unit intervals with fixed lo/hi. + * + * On the other hand, the bound v25 > -1*v85*v25 + v81 was obtained by evaluating v25 and v81 + * and the quantifier elimination based on viable::resolve_lemma didn't produce the most general lemma. + * Instead it relied on the evaluation at 353 and -1, respectively. + * + */ + + bool saturation::try_add_mul_bound(pvar x, conflict& core, inequality const& a_l_b) { + set_rule("[x] ax + b <= y, ... => a >= u_a"); + auto& m = s.var2pdd(x); + pdd const X = s.var(x); + pdd a = s.var(x); + pdd b = a, c = a, y = a; + rational a_val, b_val, c_val, y_val, x_bound; + vector x_le_bound, x_ge_bound; + signed_constraint b_bound; + if (is_AxB_l_Y(x, a_l_b, a, b, y) && !a.is_val() && s.try_eval(y, y_val) && s.try_eval(b, b_val) && s.try_eval(a, a_val) && !y_val.is_zero()) { + IF_VERBOSE(2, verbose_stream() << "v" << x << ": " << a_l_b << " a " << a << " := " << dd::val_pp(m, a_val, false) << " b " << b << " := " << dd::val_pp(m, b_val, false) << " y " << y << " := " << dd::val_pp(m, y_val, false) << "\n"); + SASSERT(!a.is_zero()); + + // define c := -b + c = -b; + VERIFY(s.try_eval(c, c_val)); + + if (has_upper_bound(x, core, x_bound, x_le_bound) && !x_le_bound.contains(a_l_b.as_signed_constraint())) { + // ax - c <= y + // ==> ax <= y + c if int(y) + int(c) <= 2^N, y <= int(y), c <= int(c) + // ==> a not in [-floor(-int(y+c) / int(x), 0[ + // ==> -a >= floor(-int(y+c) / int(x) + if (c_val + y_val <= m.max_value()) { + auto bound = floor((m.two_to_N() - y_val - c_val) / x_bound); + m_lemma.reset(); + for (auto c : x_le_bound) + m_lemma.insert_eval(~c); // x <= x_bound + m_lemma.insert_eval(~s.ule(c, c_val)); // c <= c_val + m_lemma.insert_eval(~s.ule(y, y_val)); // y <= y_val + auto conclusion = s.uge(-a, bound); // ==> -a >= bound + IF_VERBOSE(2, + verbose_stream() << core << "\n"; + verbose_stream() << "& " << X << " <= " << dd::val_pp(m, x_bound, false) << " := " << x_le_bound << "\n"; + verbose_stream() << "& " << s.ule(c, c_val) << "\n"; + verbose_stream() << "& " << s.ule(y, y_val) << "\n"; + verbose_stream() << "==> " << -a << " >= " << dd::val_pp(m, bound, false) << "\n"); + if (propagate(x, core, a_l_b, conclusion)) + return true; + } + // verbose_stream() << "TODO bound 1 " << a_l_b << " " << x_ge_bound << " " << b << " " << b_val << " " << y << "\n"; + } +#if 0 + if (has_lower_bound(x, core, x_bound, x_le_bound) && !x_le_bound.contains(a_l_b.as_signed_constraint())) { + + // verbose_stream() << "TODO bound 2 " << a_l_b << " " << x_le_bound << "\n"; + } +#endif + } + if (is_Y_l_AxB(x, a_l_b, y, a, b) && y.is_val() && s.try_eval(b, b_val)) { + // verbose_stream() << "TODO bound 3 " << a_l_b << "\n"; + } + return false; + } + + /* * TODO * @@ -949,7 +1519,7 @@ namespace polysat { return false; m_lemma.insert_eval(~d); auto conseq = s.ult(r_val, c.rhs()); - return add_conflict(core, c, conseq); + return add_conflict(v, core, c, conseq); } else { auto d = s.ule(c.rhs(), r_val); @@ -957,7 +1527,7 @@ namespace polysat { return false; m_lemma.insert_eval(~d); auto conseq = s.ule(c.lhs(), r_val); - return add_conflict(core, c, conseq); + return add_conflict(v, core, c, conseq); } } diff --git a/src/math/polysat/saturation.h b/src/math/polysat/saturation.h index 267581517..7f235a050 100644 --- a/src/math/polysat/saturation.h +++ b/src/math/polysat/saturation.h @@ -31,9 +31,11 @@ namespace polysat { bool is_non_overflow(pdd const& x, pdd const& y, signed_constraint& c); signed_constraint ineq(bool strict, pdd const& lhs, pdd const& rhs); - bool propagate(conflict& core, inequality const& crit1, signed_constraint c); - bool add_conflict(conflict& core, inequality const& crit1, signed_constraint c); - bool add_conflict(conflict& core, inequality const& crit1, inequality const& crit2, signed_constraint c); + void log_lemma(pvar v, conflict& core); + bool propagate(pvar v, conflict& core, signed_constraint const& crit1, signed_constraint c); + bool propagate(pvar v, conflict& core, inequality const& crit1, signed_constraint c); + bool add_conflict(pvar v, conflict& core, inequality const& crit1, signed_constraint c); + bool add_conflict(pvar v, conflict& core, inequality const& crit1, inequality const& crit2, signed_constraint c); bool try_ugt_x(pvar v, conflict& core, inequality const& c); @@ -49,10 +51,16 @@ namespace polysat { bool try_parity(pvar x, conflict& core, inequality const& axb_l_y); bool try_parity_diseq(pvar x, conflict& core, inequality const& axb_l_y); bool try_mul_bounds(pvar x, conflict& core, inequality const& axb_l_y); - bool try_factor_equality(pvar x, conflict& core, inequality const& a_l_b); + bool try_factor_equality1(pvar x, conflict& core, inequality const& a_l_b); + bool try_factor_equality2(pvar x, conflict& core, inequality const& a_l_b); + bool try_infer_equality(pvar x, conflict& core, inequality const& a_l_b); bool try_mul_eq_1(pvar x, conflict& core, inequality const& axb_l_y); bool try_mul_odd(pvar x, conflict& core, inequality const& axb_l_y); + bool try_mul_eq_bound(pvar x, conflict& core, inequality const& axb_l_y); + bool try_transitivity(pvar x, conflict& core, inequality const& axb_l_y); bool try_tangent(pvar v, conflict& core, inequality const& c); + bool try_add_overflow_bound(pvar x, conflict& core, inequality const& axb_l_y); + bool try_add_mul_bound(pvar x, conflict& core, inequality const& axb_l_y); // c := lhs ~ v // where ~ is < or <= @@ -62,7 +70,10 @@ namespace polysat { bool is_g_v(pvar v, inequality const& c); // c := x ~ Y - bool is_x_l_Y(pvar x, inequality const& c, pdd& y); + bool is_x_l_Y(pvar x, inequality const& i, pdd& y); + + // c := Y ~ x + bool is_Y_l_x(pvar x, inequality const& i, pdd& y); // c := X*y ~ X*Z bool is_Xy_l_XZ(pvar y, inequality const& c, pdd& x, pdd& z); @@ -107,6 +118,17 @@ namespace polysat { // p := coeff*x*y where coeff_x = coeff*x, x a variable bool is_coeffxY(pdd const& coeff_x, pdd const& p, pdd& y); + // i := x + y >= x or x + y > x + bool is_add_overflow(pvar x, inequality const& i, pdd& y); + + bool has_upper_bound(pvar x, conflict& core, rational& bound, vector& x_ge_bound); + + bool has_lower_bound(pvar x, conflict& core, rational& bound, vector& x_le_bound); + + // determine min/max parity of polynomial + unsigned min_parity(pdd const& p); + unsigned max_parity(pdd const& p); + bool is_forced_eq(pdd const& p, rational const& val); bool is_forced_eq(pdd const& p, int i) { return is_forced_eq(p, rational(i)); } @@ -118,6 +140,13 @@ namespace polysat { bool is_forced_true(signed_constraint const& sc); + bool try_inequality(pvar v, inequality const& i, conflict& core); + + bool try_umul_ovfl(pvar v, signed_constraint const& c, conflict& core); + bool try_umul_noovfl_lo(pvar v, signed_constraint const& c, conflict& core); + bool try_umul_noovfl_bounds(pvar v, signed_constraint const& c, conflict& core); + bool try_umul_ovfl_bounds(pvar v, signed_constraint const& c, conflict& core); + public: saturation(solver& s); void perform(pvar v, conflict& core); diff --git a/src/math/polysat/simplify_clause.cpp b/src/math/polysat/simplify_clause.cpp index b7feac925..e23db0788 100644 --- a/src/math/polysat/simplify_clause.cpp +++ b/src/math/polysat/simplify_clause.cpp @@ -77,13 +77,85 @@ namespace polysat { bool simplify_clause::apply(clause& cl) { LOG_H1("Simplifying clause: " << cl); + bool simplified = false; + if (try_remove_equations(cl)) + simplified = true; #if 0 if (try_recognize_bailout(cl)) - return true; + simplified = true; #endif if (try_equal_body_subsumptions(cl)) - return true; - return false; + simplified = true; + return simplified; + } + + /** + * If we have: + * p <= q + * p - q == 0 + * Then remove the equality. + * + * If we have: + * p < q + * p - q == 0 + * Then merge into p <= q. + */ + bool simplify_clause::try_remove_equations(clause& cl) { + LOG_H2("Remove superfluous equations from: " << cl); + bool const has_eqn = any_of(cl, [&](sat::literal lit) { return s.lit2cnstr(lit).is_eq(); }); + if (!has_eqn) + return false; + bool any_removed = false; + bool_vector& should_remove = m_bools; + should_remove.fill(cl.size(), false); + for (unsigned i = cl.size(); i-- > 0; ) { + sat::literal const lit = cl[i]; + signed_constraint const c = s.lit2cnstr(lit); + if (!c->is_ule()) + continue; + if (c->is_eq()) + continue; +#if 1 + // Disable the case pto_ule().lhs(); + pdd const q = c->to_ule().rhs(); + signed_constraint const eq = s.m_constraints.find_eq(p - q); + if (!eq) + continue; + auto const eq_it = std::find(cl.begin(), cl.end(), eq.blit()); + if (eq_it == cl.end()) + continue; + unsigned const eq_idx = std::distance(cl.begin(), eq_it); + any_removed = true; + should_remove[eq_idx] = true; + if (c.is_positive()) { + // c: p <= q + // eq: p == q + LOG("Removing " << eq.blit() << ": " << eq << " because it subsumes " << cl[i] << ": " << s.lit2cnstr(cl[i])); + } + else { + // c: p > q + // eq: p == q + cl[i] = s.ule(q, p).blit(); + LOG("Merge " << eq.blit() << ": " << eq << " and " << lit << ": " << c << " to obtain " << cl[i] << ": " << s.lit2cnstr(cl[i])); + } + } + // Remove superfluous equations + if (!any_removed) + return false; + unsigned j = 0; + for (unsigned i = 0; i < cl.size(); ++i) + if (!should_remove[i]) + cl[j++] = cl[i]; + cl.m_literals.shrink(j); + return true; } // If x != k appears among the new literals, all others are superfluous. @@ -94,7 +166,7 @@ namespace polysat { sat::literal eq = sat::null_literal; rational k; for (sat::literal lit : cl) { - LOG_V("Examine " << lit_pp(s, lit)); + LOG_V(10, "Examine " << lit_pp(s, lit)); lbool status = s.m_bvars.value(lit); // skip premise literals if (status == l_false) @@ -237,7 +309,7 @@ namespace polysat { for (unsigned i = 0; i < cl.size(); ++i) { subs_entry& entry = m_entries[i]; sat::literal lit = cl[i]; - LOG("Literal " << lit_pp(s, lit)); + LOG_V(10, "Literal " << lit_pp(s, lit)); signed_constraint c = s.lit2cnstr(lit); prepare_subs_entry(entry, c); } @@ -254,7 +326,7 @@ namespace polysat { continue; if (e.interval.currently_contains(f.interval)) { // f subset of e ==> f.src subsumed by e.src - LOG("Removing " << s.lit2cnstr(cl[i]) << " because it subsumes " << s.lit2cnstr(cl[j])); + LOG("Removing " << cl[i] << ": " << s.lit2cnstr(cl[i]) << " because it subsumes " << cl[j] << ": " << s.lit2cnstr(cl[j])); e.subsuming = true; any_subsumed = true; break; diff --git a/src/math/polysat/simplify_clause.h b/src/math/polysat/simplify_clause.h index 51963083b..75f17afca 100644 --- a/src/math/polysat/simplify_clause.h +++ b/src/math/polysat/simplify_clause.h @@ -28,7 +28,9 @@ namespace polysat { solver& s; vector m_entries; + bool_vector m_bools; + bool try_remove_equations(clause& cl); bool try_recognize_bailout(clause& cl); bool try_equal_body_subsumptions(clause& cl); diff --git a/src/math/polysat/smul_fl_constraint.cpp b/src/math/polysat/smul_fl_constraint.cpp index 98f76fac1..92cbb5f63 100644 --- a/src/math/polysat/smul_fl_constraint.cpp +++ b/src/math/polysat/smul_fl_constraint.cpp @@ -15,8 +15,8 @@ Author: namespace polysat { - smul_fl_constraint::smul_fl_constraint(constraint_manager& m, pdd const& p, pdd const& q, bool is_overflow): - constraint(m, ckind_t::smul_fl_t), m_is_overflow(is_overflow), m_p(p), m_q(q) { + smul_fl_constraint::smul_fl_constraint(pdd const& p, pdd const& q, bool is_overflow): + constraint(ckind_t::smul_fl_t), m_is_overflow(is_overflow), m_p(p), m_q(q) { simplify(); m_vars.append(m_p.free_vars()); for (auto v : m_q.free_vars()) @@ -119,13 +119,17 @@ namespace polysat { && q() == other.to_smul_fl().q(); } - void smul_fl_constraint::add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { - auto p_coeff = s.subst(p()).get_univariate_coefficients(); - auto q_coeff = s.subst(q()).get_univariate_coefficients(); + void smul_fl_constraint::add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { + auto p1 = s.subst(p()); + if (!p1.is_univariate_in(v)) + return; + auto q1 = s.subst(q()); + if (!q1.is_univariate_in(v)) + return; if (is_overflow()) - us.add_smul_ovfl(p_coeff, q_coeff, !is_positive, dep); + us.add_smul_ovfl(p1.get_univariate_coefficients(), q1.get_univariate_coefficients(), !is_positive, dep); else - us.add_smul_udfl(p_coeff, q_coeff, !is_positive, dep); + us.add_smul_udfl(p1.get_univariate_coefficients(), q1.get_univariate_coefficients(), !is_positive, dep); } } diff --git a/src/math/polysat/smul_fl_constraint.h b/src/math/polysat/smul_fl_constraint.h index a8d2c2814..07b740ea0 100644 --- a/src/math/polysat/smul_fl_constraint.h +++ b/src/math/polysat/smul_fl_constraint.h @@ -25,7 +25,7 @@ namespace polysat { pdd m_q; void simplify(); - smul_fl_constraint(constraint_manager& m, pdd const& p, pdd const& q, bool is_overflow); + smul_fl_constraint(pdd const& p, pdd const& q, bool is_overflow); public: ~smul_fl_constraint() override {} @@ -42,7 +42,7 @@ namespace polysat { bool operator==(constraint const& other) const override; bool is_eq() const override { return false; } - void add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; + void add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; }; } diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp index 9d67e94c5..5260fa11c 100644 --- a/src/math/polysat/solver.cpp +++ b/src/math/polysat/solver.cpp @@ -53,6 +53,7 @@ namespace polysat { m_config.m_max_conflicts = m_params.get_uint("max_conflicts", UINT_MAX); m_config.m_max_decisions = m_params.get_uint("max_decisions", UINT_MAX); m_config.m_log_conflicts = pp.log_conflicts(); + m_rand.set_seed(m_params.get_uint("random_seed", 0)); } bool solver::should_search() { @@ -75,6 +76,7 @@ namespace polysat { LOG("Assignment: " << assignments_pp(*this)); if (is_conflict()) LOG("Conflict: " << m_conflict); IF_LOGGING(m_viable.log()); + SASSERT(var_queue_invariant()); if (is_conflict() && at_base_level()) { LOG_H2("UNSAT"); return l_false; } else if (is_conflict()) resolve_conflict(); else if (should_add_pwatch()) add_pwatch(); @@ -241,9 +243,11 @@ namespace polysat { for (; i < sz && !is_conflict(); ++i) if (!propagate(v, wlist[i])) wlist[j++] = wlist[i]; - for (; i < sz; ++i) + for (; i < sz; ++i) wlist[j++] = wlist[i]; wlist.shrink(j); + if (is_conflict()) + shuffle(wlist.size(), wlist.data(), m_rand); DEBUG_CODE(m_locked_wlist = std::nullopt;); } @@ -394,7 +398,7 @@ namespace polysat { void solver::add_pwatch(constraint* c, pvar v) { SASSERT(m_locked_wlist != v); // the propagate loop will not discover the new size - LOG_V("Watching v" << v << " in constraint " << show_deref(c)); + LOG_V(20, "Watching v" << v << " in constraint " << show_deref(c)); m_pwatch[v].push_back(c); } @@ -495,7 +499,7 @@ namespace polysat { } case trail_instr_t::assign_i: { auto v = m_search.back().var(); - LOG_V("Undo assign_i: v" << v); + LOG_V(20, "Undo assign_i: v" << v); unsigned active_level = get_level(v); if (active_level <= target_level) { @@ -511,7 +515,7 @@ namespace polysat { } case trail_instr_t::assign_bool_i: { sat::literal lit = m_search.back().lit(); - LOG_V("Undo assign_bool_i: " << lit); + LOG_V(20, "Undo assign_bool_i: " << lit); unsigned active_level = m_bvars.level(lit); if (active_level <= target_level) @@ -588,15 +592,13 @@ namespace polysat { // Simple hack to try deciding the boolean skeleton first if (!can_bdecide()) { // enqueue all not-yet-true clauses - for (auto const& cls : m_constraints.clauses()) { - for (auto const& cl : cls) { - bool is_true = any_of(*cl, [&](sat::literal lit) { return m_bvars.is_true(lit); }); - if (is_true) - continue; - size_t undefs = count_if(*cl, [&](sat::literal lit) { return !m_bvars.is_assigned(lit); }); - if (undefs >= 2) - m_lemmas.push_back(cl.get()); - } + for (clause const& cl : m_constraints.clauses()) { + bool is_true = any_of(cl, [&](sat::literal lit) { return m_bvars.is_true(lit); }); + if (is_true) + continue; + size_t undefs = count_if(cl, [&](sat::literal lit) { return !m_bvars.is_assigned(lit); }); + if (undefs >= 2) + m_lemmas.push_back(&cl); } } #endif @@ -608,7 +610,7 @@ namespace polysat { } void solver::bdecide() { - clause& lemma = *m_lemmas[m_lemmas_qhead++]; + clause const& lemma = *m_lemmas[m_lemmas_qhead++]; on_scope_exit update_trail = [this]() { // must be done after push_level, but also if we return early. m_trail.push_back(trail_instr_t::lemma_qhead_i); @@ -654,7 +656,7 @@ namespace polysat { // (fail here in debug mode so we notice if we miss some) DEBUG_CODE( UNREACHABLE(); ); m_free_pvars.unassign_var_eh(v); - set_conflict(v, false); + SASSERT(is_conflict()); return; case find_t::singleton: // NOTE: this case may happen legitimately if all other possibilities were excluded by brute force search @@ -734,9 +736,16 @@ namespace polysat { SASSERT(!m_search.assignment().contains(v)); if (lemma) { add_clause(*lemma); - SASSERT(!is_conflict()); // if we have a conflict here, we could have produced this lemma already earlier - if (can_propagate()) + if (is_conflict()) { + // If we have a conflict here, we should have produced this lemma already earlier + LOG("Conflict after constraint::produce_lemma: TODO: should have found this lemma earlier"); + m_free_pvars.unassign_var_eh(v); return; + } + if (can_propagate()) { + m_free_pvars.unassign_var_eh(v); + return; + } } if (c) { LOG_H2("Chosen assignment " << assignment_pp(*this, v, val) << " is not actually viable!"); @@ -753,10 +762,11 @@ namespace polysat { case find_t::empty: LOG("Fallback solver: unsat"); m_free_pvars.unassign_var_eh(v); - set_conflict(v, true); + SASSERT(is_conflict()); return; case find_t::resource_out: UNREACHABLE(); // TODO: abort solving + m_free_pvars.unassign_var_eh(v); return; } } @@ -938,7 +948,6 @@ namespace polysat { clause* best_lemma = nullptr; auto appraise_lemma = [&](clause* lemma) { - m_simplify_clause.apply(*lemma); auto score = compute_lemma_score(*lemma); if (score) LOG(" score: " << *score); @@ -958,14 +967,14 @@ namespace polysat { appraise_lemma(lemmas.back()); } SASSERT(best_score < lemma_score::max()); - SASSERT(best_lemma); + VERIFY(best_lemma); unsigned const jump_level = std::max(best_score.jump_level(), base_level()); SASSERT(jump_level <= max_jump_level); LOG("best_score: " << best_score); LOG("best_lemma: " << *best_lemma); - + m_conflict.reset(); backjump(jump_level); @@ -1054,7 +1063,11 @@ namespace polysat { void solver::assign_eval(sat::literal lit) { signed_constraint const c = lit2cnstr(lit); + LOG_V(10, "Evaluate: " << lit_pp(*this ,lit)); + // assertion is false + if (!c.is_currently_true(*this)) IF_VERBOSE(0, verbose_stream() << c << " is not currently true\n"); SASSERT(c.is_currently_true(*this)); + VERIFY(c.is_currently_true(*this)); unsigned level = 0; // NOTE: constraint may be evaluated even if some variables are still unassigned (e.g., 0*x doesn't depend on x). for (auto v : c->vars()) @@ -1309,12 +1322,10 @@ namespace polysat { for (auto c : m_constraints) out << "\t" << c->bvar2string() << ": " << *c << "\n"; out << "Clauses:\n"; - for (auto const& cls : m_constraints.clauses()) { - for (auto const& cl : cls) { - out << "\t" << *cl << "\n"; - for (auto lit : *cl) - out << "\t\t" << lit << ": " << lit2cnstr(lit) << "\n"; - } + for (clause const& cl : m_constraints.clauses()) { + out << "\t" << cl << "\n"; + for (sat::literal lit : cl) + out << "\t\t" << lit << ": " << lit2cnstr(lit) << "\n"; } return out; } @@ -1430,23 +1441,20 @@ namespace polysat { // Check for missed boolean propagations: // - no clause should have exactly one unassigned literal, unless it is already true. // - no clause should be false - for (auto const& cls : m_constraints.clauses()) { - for (auto const& clref : cls) { - clause const& cl = *clref; - bool const is_true = any_of(cl, [&](auto lit) { return m_bvars.is_true(lit); }); - if (is_true) - continue; - size_t const undefs = count_if(cl, [&](auto lit) { return !m_bvars.is_assigned(lit); }); - if (undefs == 1) { - LOG("Missed boolean propagation of clause: " << cl); - for (sat::literal lit : cl) { - LOG(" " << lit_pp(*this, lit)); - } + for (clause const& cl : m_constraints.clauses()) { + bool const is_true = any_of(cl, [&](auto lit) { return m_bvars.is_true(lit); }); + if (is_true) + continue; + size_t const undefs = count_if(cl, [&](auto lit) { return !m_bvars.is_assigned(lit); }); + if (undefs == 1) { + LOG("Missed boolean propagation of clause: " << cl); + for (sat::literal lit : cl) { + LOG(" " << lit_pp(*this, lit)); } - SASSERT(undefs != 1); - bool const is_false = all_of(cl, [&](auto lit) { return m_bvars.is_false(lit); }); - SASSERT(!is_false); } + SASSERT(undefs != 1); + bool const is_false = all_of(cl, [&](auto lit) { return m_bvars.is_false(lit); }); + SASSERT(!is_false); } return true; } @@ -1456,7 +1464,7 @@ namespace polysat { } /** Check that boolean assignment and constraint evaluation are consistent */ - bool solver::assignment_invariant() { + bool solver::assignment_invariant() const { if (is_conflict()) return true; bool ok = true; @@ -1475,6 +1483,25 @@ namespace polysat { return ok; } + /** Check that each variable is either assigned or queued for decisions */ + bool solver::var_queue_invariant() const { + if (is_conflict()) + return true; + uint_set active; + for (pvar v : m_free_pvars) + active.insert(v); + for (auto const& [v, val] : assignment()) + active.insert(v); + bool ok = true; + for (pvar v = 0; v < num_vars(); ++v) { + if (!active.contains(v)) { + ok = false; + LOG("Lost variable v" << v << " (it is neither assigned nor free)"); + } + } + return ok; + } + /// Check that all constraints on the stack are satisfied by the current model. bool solver::verify_sat() { LOG_H1("Checking current model..."); @@ -1487,19 +1514,17 @@ namespace polysat { all_ok = all_ok && ok; } } - for (auto clauses : m_constraints.clauses()) { - for (auto cl : clauses) { - bool clause_ok = false; - for (sat::literal lit : *cl) { - bool ok = lit2cnstr(lit).is_currently_true(*this); - if (ok) { - clause_ok = true; - break; - } + for (clause const& cl : m_constraints.clauses()) { + bool clause_ok = false; + for (sat::literal lit : cl) { + bool ok = lit2cnstr(lit).is_currently_true(*this); + if (ok) { + clause_ok = true; + break; } - LOG((clause_ok ? "PASS" : "FAIL") << ": " << show_deref(cl) << (cl->is_redundant() ? " (redundant)" : "")); - all_ok = all_ok && clause_ok; } + LOG((clause_ok ? "PASS" : "FAIL") << ": " << cl << (cl.is_redundant() ? " (redundant)" : "")); + all_ok = all_ok && clause_ok; } if (all_ok) LOG("All good!"); return all_ok; diff --git a/src/math/polysat/solver.h b/src/math/polysat/solver.h index e74e900f8..ff964020e 100644 --- a/src/math/polysat/solver.h +++ b/src/math/polysat/solver.h @@ -157,6 +157,7 @@ namespace polysat { bool_var_manager m_bvars; // Map boolean variables to constraints var_queue m_free_pvars; // free poly vars stats m_stats; + random_gen m_rand; config m_config; // Per constraint state @@ -188,7 +189,7 @@ namespace polysat { constraints m_pwatch_trail; #endif - ptr_vector m_lemmas; ///< the non-asserting lemmas + ptr_vector m_lemmas; ///< the non-asserting lemmas unsigned m_lemmas_qhead = 0; unsigned_vector m_base_levels; // External clients can push/pop scope. @@ -215,6 +216,8 @@ namespace polysat { dd::pdd_manager& sz2pdd(unsigned sz) const; dd::pdd_manager& var2pdd(pvar v) const; + pvar num_vars() const { return m_value.size(); } + assignment_t const& assignment() const { return m_search.assignment(); } void push_level(); @@ -251,7 +254,8 @@ namespace polysat { void set_conflict_at_base_level() { m_conflict.init_at_base_level(); } void set_conflict(signed_constraint c) { m_conflict.init(c); } void set_conflict(clause& cl) { m_conflict.init(cl); } - void set_conflict(pvar v, bool by_viable_fallback) { m_conflict.init(v, by_viable_fallback); } + void set_conflict_by_viable_interval(pvar v) { m_conflict.init_by_viable_interval(v); } + void set_conflict_by_viable_fallback(pvar v, univariate_solver& us) { m_conflict.init_by_viable_fallback(v, us); } bool can_decide() const; bool can_bdecide() const; @@ -309,7 +313,8 @@ namespace polysat { static bool invariant(signed_constraints const& cs); bool wlist_invariant() const; bool bool_watch_invariant() const; - bool assignment_invariant(); + bool assignment_invariant() const; + bool var_queue_invariant() const; bool verify_sat(); bool can_propagate(); @@ -417,15 +422,33 @@ namespace polysat { signed_constraint eq(pdd const& p, unsigned q) { return eq(p - q); } signed_constraint odd(pdd const& p) { return ~even(p); } signed_constraint even(pdd const& p) { return parity(p, 1); } - /** parity(p) >= k (<=> p * 2^(K-k) == 0) */ - signed_constraint parity(pdd const& p, unsigned k) { + /** parity(p) >= k */ + signed_constraint parity(pdd const& p, unsigned k) { // TODO: rename to parity_at_least? unsigned N = p.manager().power_of_2(); + // parity(p) >= k + // <=> p * 2^(N - k) == 0 if (k >= N) return eq(p); - else if (k == 0) - return odd(p); + else if (k == 0) { + // parity(p) >= 0 is a tautology + verbose_stream() << "REDUNDANT parity constraint: parity(" << p << ", " << k << ")\n"; + return eq(p.manager().zero()); + } else - return eq(p*rational::power_of_two(N - k)); + return eq(p * rational::power_of_two(N - k)); + } + /** parity(p) <= k */ + signed_constraint parity_at_most(pdd const& p, unsigned k) { + unsigned N = p.manager().power_of_2(); + // parity(p) <= k + // <=> ~(parity(p) >= k+1) + if (k >= N) { + // parity(p) <= N is a tautology + verbose_stream() << "REDUNDANT parity constraint: parity(" << p << ", " << k << ")\n"; + return eq(p.manager().zero()); + } + else + return ~parity(p, k + 1); } signed_constraint diseq(pdd const& p, rational const& q) { return diseq(p - q); } signed_constraint diseq(pdd const& p, unsigned q) { return diseq(p - q); } diff --git a/src/math/polysat/ule_constraint.cpp b/src/math/polysat/ule_constraint.cpp index 59d2675ea..e3ba7854d 100644 --- a/src/math/polysat/ule_constraint.cpp +++ b/src/math/polysat/ule_constraint.cpp @@ -135,8 +135,8 @@ namespace { namespace polysat { - ule_constraint::ule_constraint(constraint_manager& m, pdd const& l, pdd const& r) : - constraint(m, ckind_t::ule_t), m_lhs(l), m_rhs(r) { + ule_constraint::ule_constraint(pdd const& l, pdd const& r) : + constraint(ckind_t::ule_t), m_lhs(l), m_rhs(r) { m_vars.append(m_lhs.free_vars()); for (auto v : m_rhs.free_vars()) @@ -185,9 +185,9 @@ namespace polysat { signed_constraint sc(this, is_positive); LOG_H3("Narrowing " << sc); - LOG_V("Assignment: " << assignments_pp(s)); - LOG_V("Substituted LHS: " << lhs() << " := " << p); - LOG_V("Substituted RHS: " << rhs() << " := " << q); + LOG_V(10, "Assignment: " << assignments_pp(s)); + LOG_V(10, "Substituted LHS: " << lhs() << " := " << p); + LOG_V(10, "Substituted RHS: " << rhs() << " := " << q); if (is_always_false(is_positive, p, q)) { s.set_conflict(sc); @@ -353,9 +353,16 @@ namespace polysat { return other.is_ule() && lhs() == other.to_ule().lhs() && rhs() == other.to_ule().rhs(); } - void ule_constraint::add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { - auto p_coeff = s.subst(lhs()).get_univariate_coefficients(); - auto q_coeff = s.subst(rhs()).get_univariate_coefficients(); - us.add_ule(p_coeff, q_coeff, !is_positive, dep); + void ule_constraint::add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { + pdd p = s.subst(lhs()); + pdd q = s.subst(rhs()); + bool p_ok = p.is_univariate_in(v); + bool q_ok = q.is_univariate_in(v); + if (!is_positive && !q_ok) // add p > 0 + us.add_ugt(p.get_univariate_coefficients(), rational::zero(), false, dep); + if (!is_positive && !p_ok) // add -1 > q <==> q+1 > 0 + us.add_ugt((q + 1).get_univariate_coefficients(), rational::zero(), false, dep); + if (p_ok && q_ok) + us.add_ule(p.get_univariate_coefficients(), q.get_univariate_coefficients(), !is_positive, dep); } } diff --git a/src/math/polysat/ule_constraint.h b/src/math/polysat/ule_constraint.h index a798cae47..10b5a502d 100644 --- a/src/math/polysat/ule_constraint.h +++ b/src/math/polysat/ule_constraint.h @@ -25,7 +25,7 @@ namespace polysat { pdd m_lhs; pdd m_rhs; - ule_constraint(constraint_manager& m, pdd const& l, pdd const& r); + ule_constraint(pdd const& l, pdd const& r); static void simplify(bool& is_positive, pdd& lhs, pdd& rhs); static bool is_always_true(bool is_positive, pdd const& lhs, pdd const& rhs) { return eval(lhs, rhs) == to_lbool(is_positive); } static bool is_always_false(bool is_positive, pdd const& lhs, pdd const& rhs) { return is_always_true(!is_positive, lhs, rhs); } @@ -45,7 +45,7 @@ namespace polysat { unsigned hash() const override; bool operator==(constraint const& other) const override; bool is_eq() const override { return m_rhs.is_zero(); } - void add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; + void add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; }; struct ule_pp { diff --git a/src/math/polysat/umul_ovfl_constraint.cpp b/src/math/polysat/umul_ovfl_constraint.cpp index f1f980e20..64d4253ea 100644 --- a/src/math/polysat/umul_ovfl_constraint.cpp +++ b/src/math/polysat/umul_ovfl_constraint.cpp @@ -15,8 +15,8 @@ Author: namespace polysat { - umul_ovfl_constraint::umul_ovfl_constraint(constraint_manager& m, pdd const& p, pdd const& q): - constraint(m, ckind_t::umul_ovfl_t), m_p(p), m_q(q) { + umul_ovfl_constraint::umul_ovfl_constraint(pdd const& p, pdd const& q): + constraint(ckind_t::umul_ovfl_t), m_p(p), m_q(q) { simplify(); m_vars.append(m_p.free_vars()); for (auto v : m_q.free_vars()) @@ -25,8 +25,7 @@ namespace polysat { } void umul_ovfl_constraint::simplify() { - if (m_p.is_zero() || m_q.is_zero() || - m_p.is_one() || m_q.is_one()) { + if (m_p.is_zero() || m_q.is_zero() || m_p.is_one() || m_q.is_one()) { m_q = 0; m_p = 0; return; @@ -99,8 +98,8 @@ namespace polysat { signed_constraint sc(this, is_positive); // ¬Omega(p, q) ==> q = 0 \/ p <= p*q // ¬Omega(p, q) ==> p = 0 \/ q <= p*q - s.add_clause(~sc, /* s.eq(p()), */ s.eq(q()), s.ule(p(), p()*q()), false); - s.add_clause(~sc, s.eq(p()), /* s.eq(q()), */ s.ule(q(), p()*q()), false); + s.add_clause(~sc, s.eq(q()), s.ule(p(), p()*q()), false); + s.add_clause(~sc, s.eq(p()), s.ule(q(), p()*q()), false); } } @@ -158,9 +157,13 @@ namespace polysat { return other.is_umul_ovfl() && p() == other.to_umul_ovfl().p() && q() == other.to_umul_ovfl().q(); } - void umul_ovfl_constraint::add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { - auto p_coeff = s.subst(p()).get_univariate_coefficients(); - auto q_coeff = s.subst(q()).get_univariate_coefficients(); - us.add_umul_ovfl(p_coeff, q_coeff, !is_positive, dep); + void umul_ovfl_constraint::add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const { + pdd p1 = s.subst(p()); + if (!p1.is_univariate_in(v)) + return; + pdd q1 = s.subst(q()); + if (!q1.is_univariate_in(v)) + return; + us.add_umul_ovfl(p1.get_univariate_coefficients(), q1.get_univariate_coefficients(), !is_positive, dep); } } diff --git a/src/math/polysat/umul_ovfl_constraint.h b/src/math/polysat/umul_ovfl_constraint.h index 1595f0bd6..bb270d93a 100644 --- a/src/math/polysat/umul_ovfl_constraint.h +++ b/src/math/polysat/umul_ovfl_constraint.h @@ -23,7 +23,7 @@ namespace polysat { pdd m_p; pdd m_q; - umul_ovfl_constraint(constraint_manager& m, pdd const& p, pdd const& q); + umul_ovfl_constraint(pdd const& p, pdd const& q); void simplify(); static bool is_always_true(bool is_positive, pdd const& p, pdd const& q) { return eval(p, q) == to_lbool(is_positive); } static bool is_always_false(bool is_positive, pdd const& p, pdd const& q) { return is_always_true(!is_positive, p, q); } @@ -46,7 +46,7 @@ namespace polysat { bool operator==(constraint const& other) const override; bool is_eq() const override { return false; } - void add_to_univariate_solver(solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; + void add_to_univariate_solver(pvar v, solver& s, univariate_solver& us, unsigned dep, bool is_positive) const override; }; } diff --git a/src/math/polysat/univariate/univariate_solver.cpp b/src/math/polysat/univariate/univariate_solver.cpp index e66a521c3..73a2702ce 100644 --- a/src/math/polysat/univariate/univariate_solver.cpp +++ b/src/math/polysat/univariate/univariate_solver.cpp @@ -26,6 +26,12 @@ Author: namespace polysat { + univariate_solver::dep_vector univariate_solver::unsat_core() { + dep_vector deps; + unsat_core(deps); + return deps; + } + class univariate_bitblast_solver : public univariate_solver { // TODO: does it make sense to share m and bv between different solver instances? // TODO: consider pooling solvers to save setup overhead, see if solver/solver_pool.h can be used @@ -33,6 +39,7 @@ namespace polysat { scoped_ptr bv; scoped_ptr s; unsigned bit_width; + unsigned m_scope_level = 0; func_decl_ref x_decl; expr_ref x; vector model_cache; @@ -46,6 +53,7 @@ namespace polysat { bv = alloc(bv_util, m); params_ref p; p.set_bool("bv.polysat", false); + // p.set_bool("smt", true); s = mk_solver(m, p, false, true, true, symbol::null); x_decl = m.mk_const_decl("x", bv->mk_sort(bit_width)); x = m.mk_const(x_decl); @@ -59,7 +67,8 @@ namespace polysat { } void push_cache() { - model_cache.push_back(model_cache.back()); + rational v = model_cache.back(); + model_cache.push_back(v); } void pop_cache() { @@ -67,19 +76,37 @@ namespace polysat { } void push() override { + m_scope_level++; push_cache(); s->push(); } void pop(unsigned n) override { + SASSERT(scope_level() >= n); + m_scope_level -= n; pop_cache(); s->pop(n); } + unsigned scope_level() override { + return m_scope_level; + } + expr* mk_numeral(rational const& r) const { return bv->mk_numeral(r, bit_width); } + bool is_zero(univariate const& p) const { + for (auto n : p) + if (n != 0) + return false; + return true; + } + + bool is_zero(rational const& p) const { + return p.is_zero(); + } + #if 0 // [d,c,b,a] --> ((a*x + b)*x + c)*x + d expr* mk_poly(univariate const& p) const { @@ -97,35 +124,43 @@ namespace polysat { } } #else - // TODO: shouldn't the simplification step of the underlying solver already support this transformation? how to enable? // 2^k*x --> x << k // n*x --> n * x expr* mk_poly_term(rational const& coeff, expr* xpow) const { unsigned pow; + SASSERT(!coeff.is_zero()); + if (coeff.is_one()) + return xpow; if (coeff.is_power_of_two(pow)) return bv->mk_bv_shl(xpow, mk_numeral(rational(pow))); - else - return bv->mk_bv_mul(mk_numeral(coeff), xpow); + return bv->mk_bv_mul(mk_numeral(coeff), xpow); } // [d,c,b,a] --> d + c*x + b*(x*x) + a*(x*x*x) - expr* mk_poly(univariate const& p) const { - if (p.empty()) { - return mk_numeral(rational::zero()); - } + expr_ref mk_poly(univariate const& p) { + expr_ref e(m); + if (p.empty()) + e = mk_numeral(rational::zero()); else { - expr* e = mk_numeral(p[0]); - expr* xpow = x; + if (!p[0].is_zero()) + e = mk_numeral(p[0]); + expr_ref xpow = x; for (unsigned i = 1; i < p.size(); ++i) { if (!p[i].is_zero()) { expr* t = mk_poly_term(p[i], xpow); - e = bv->mk_bv_add(e, t); + e = e ? bv->mk_bv_add(e, t) : t; } if (i + 1 < p.size()) xpow = bv->mk_bv_mul(xpow, x); } - return e; + if (!e) + e = mk_numeral(p[0]); } + return e; + } + + expr_ref mk_poly(rational const& p) { + return {mk_numeral(p), m}; } #endif @@ -133,15 +168,29 @@ namespace polysat { reset_cache(); if (sign) e = m.mk_not(e); - expr* a = m.mk_const(m.mk_const_decl(symbol(dep), m.mk_bool_sort())); - s->assert_expr(e, a); - IF_VERBOSE(10, verbose_stream() << "(assert (! " << expr_ref(e, m) << " :named " << expr_ref(a, m) << "))\n"); + if (dep == null_dep) { + s->assert_expr(e); + IF_VERBOSE(10, verbose_stream() << "(assert " << expr_ref(e, m) << ")\n"); + } + else { + expr* a = m.mk_const(m.mk_const_decl(symbol(dep), m.mk_bool_sort())); + s->assert_expr(e, a); + IF_VERBOSE(10, verbose_stream() << "(assert (! " << expr_ref(e, m) << " :named " << expr_ref(a, m) << "))\n"); + } } - void add_ule(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) override { - add(bv->mk_ule(mk_poly(lhs), mk_poly(rhs)), sign, dep); + template + void add_ule_impl(lhs_t const& lhs, rhs_t const& rhs, bool sign, dep_t dep) { + if (is_zero(rhs)) + add(m.mk_eq(mk_poly(lhs), mk_poly(rhs)), sign, dep); + else + add(bv->mk_ule(mk_poly(lhs), mk_poly(rhs)), sign, dep); } + void add_ule(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) override { add_ule_impl(lhs, rhs, sign, dep); } + void add_ule(univariate const& lhs, rational const& rhs, bool sign, dep_t dep) override { add_ule_impl(lhs, rhs, sign, dep); } + void add_ule(rational const& lhs, univariate const& rhs, bool sign, dep_t dep) override { add_ule_impl(lhs, rhs, sign, dep); } + void add_umul_ovfl(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) override { add(bv->mk_bvumul_no_ovfl(mk_poly(lhs), mk_poly(rhs)), !sign, dep); } @@ -183,11 +232,14 @@ namespace polysat { } void add_ule_const(rational const& val, bool sign, dep_t dep) override { - add(bv->mk_ule(x, mk_numeral(val)), sign, dep); + if (val == 0) + add(m.mk_eq(x, mk_poly(val)), sign, dep); + else + add(bv->mk_ule(x, mk_poly(val)), sign, dep); } void add_uge_const(rational const& val, bool sign, dep_t dep) override { - add(bv->mk_ule(mk_numeral(val), x), sign, dep); + add(bv->mk_ule(mk_poly(val), x), sign, dep); } void add_bit(unsigned idx, bool sign, dep_t dep) override { @@ -198,16 +250,16 @@ namespace polysat { return s->check_sat(); } - dep_vector unsat_core() override { - dep_vector deps; + void unsat_core(dep_vector& deps) override { + deps.reset(); expr_ref_vector core(m); s->get_unsat_core(core); for (expr* a : core) { unsigned dep = to_app(a)->get_decl()->get_name().get_num(); deps.push_back(dep); } + IF_VERBOSE(10, verbose_stream() << "core " << deps << "\n"); SASSERT(deps.size() > 0); - return deps; } rational model() override { @@ -232,7 +284,7 @@ namespace polysat { // try reducing val by setting bits to 0, starting at the msb. for (unsigned k = bit_width; k-- > 0; ) { if (!val.get_bit(k)) { - add_bit0(k, 0); + add_bit0(k, null_dep); continue; } // try decreasing k-th bit @@ -245,9 +297,9 @@ namespace polysat { } pop(1); if (result == l_true) - add_bit0(k, 0); + add_bit0(k, null_dep); else if (result == l_false) - add_bit1(k, 0); + add_bit1(k, null_dep); else return false; } @@ -274,9 +326,9 @@ namespace polysat { } pop(1); if (result == l_true) - add_bit1(k, 0); + add_bit1(k, null_dep); else if (result == l_false) - add_bit0(k, 0); + add_bit0(k, null_dep); else return false; } @@ -284,6 +336,27 @@ namespace polysat { return true; } + bool find_two(rational& out1, rational& out2) override { + out1 = model(); + bool ok = true; + push(); + add(m.mk_eq(mk_numeral(out1), x), true, null_dep); + switch (check()) { + case l_true: + out2 = model(); + break; + case l_false: + out2 = out1; + break; + default: + ok = false; + break; + } + pop(1); + IF_VERBOSE(10, verbose_stream() << "viable " << out1 << " " << out2 << "\n"); + return ok; + } + std::ostream& display(std::ostream& out) const override { return out << *s; } diff --git a/src/math/polysat/univariate/univariate_solver.h b/src/math/polysat/univariate/univariate_solver.h index 904342a79..a90ebe825 100644 --- a/src/math/polysat/univariate/univariate_solver.h +++ b/src/math/polysat/univariate/univariate_solver.h @@ -35,17 +35,21 @@ namespace polysat { /// e.g., the vector [ c, b, a ] represents a*x^2 + b*x + c. using univariate = vector; + const dep_t null_dep = UINT_MAX; + virtual ~univariate_solver() = default; virtual void push() = 0; virtual void pop(unsigned n) = 0; + virtual unsigned scope_level() = 0; virtual lbool check() = 0; /** * Precondition: check() returned l_false */ - virtual dep_vector unsat_core() = 0; + dep_vector unsat_core(); + virtual void unsat_core(dep_vector& out_deps) = 0; /** * Precondition: check() returned l_true @@ -68,7 +72,31 @@ namespace polysat { */ virtual bool find_max(rational& out_max) = 0; + /** + * Find up to two viable values. + * + * Precondition: check() returned l_true + * returns: true on success, false on resource out + */ + virtual bool find_two(rational& out1, rational& out2) = 0; + + /** lhs <= rhs */ virtual void add_ule(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) = 0; + virtual void add_ule(univariate const& lhs, rational const& rhs, bool sign, dep_t dep) = 0; + virtual void add_ule(rational const& lhs, univariate const& rhs, bool sign, dep_t dep) = 0; + /** lhs >= rhs */ + void add_uge(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) { add_ule(rhs, lhs, sign, dep); } + void add_uge(univariate const& lhs, rational const& rhs, bool sign, dep_t dep) { add_ule(rhs, lhs, sign, dep); } + void add_uge(rational const& lhs, univariate const& rhs, bool sign, dep_t dep) { add_ule(rhs, lhs, sign, dep); } + /** lhs < rhs */ + void add_ult(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) { add_ule(rhs, lhs, !sign, dep); } + void add_ult(univariate const& lhs, rational const& rhs, bool sign, dep_t dep) { add_ule(rhs, lhs, !sign, dep); } + void add_ult(rational const& lhs, univariate const& rhs, bool sign, dep_t dep) { add_ule(rhs, lhs, !sign, dep); } + /** lhs > rhs */ + void add_ugt(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) { add_ule(lhs, rhs, !sign, dep); } + void add_ugt(univariate const& lhs, rational const& rhs, bool sign, dep_t dep) { add_ule(lhs, rhs, !sign, dep); } + void add_ugt(rational const& lhs, univariate const& rhs, bool sign, dep_t dep) { add_ule(lhs, rhs, !sign, dep); } + virtual void add_umul_ovfl(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) = 0; virtual void add_smul_ovfl(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) = 0; virtual void add_smul_udfl(univariate const& lhs, univariate const& rhs, bool sign, dep_t dep) = 0; diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp index 3b832b7cf..1212d803a 100644 --- a/src/math/polysat/viable.cpp +++ b/src/math/polysat/viable.cpp @@ -23,16 +23,26 @@ TODO: improve management of the fallback univariate solvers: - incrementally add/remove constraints - set resource limit of univariate solver +TODO: plan to fix the FI "pumping": + 1. simple looping detection and bitblasting fallback. -- done + 2. intervals at multiple bit widths + - for equations, this will give us exact solutions for all coefficients + - for inequalities, a coefficient 2^k*a means that intervals are periodic because the upper k bits of x are irrelevant; + storing the interval for x[K-k:0] would take care of this. + --*/ #include "util/debug.h" #include "math/polysat/viable.h" #include "math/polysat/solver.h" +#include "math/polysat/number.h" #include "math/polysat/univariate/univariate_solver.h" namespace polysat { + using namespace viable_query; + struct inf_fi : public inference { viable& v; pvar var; @@ -137,7 +147,7 @@ namespace polysat { prop = true; break; case find_t::empty: - s.set_conflict(v, false); + SASSERT(s.is_conflict()); return true; default: break; @@ -159,6 +169,21 @@ namespace polysat { // - side conditions // - i.lo() == i.lo_val() for each unit interval i // - i.hi() == i.hi_val() for each unit interval i + + // NSB review: + // the bounds added by x < p and p < x in forbidden_intervals + // match_non_max, match_non_zero + // use values that are approximations. Then the propagations in + // try_assign_eval are incorrect. + // For example, x > p means x has forbidden interval [0, p + 1[, + // the numeric interval is [0, 1[, but p + 1 == 1 is not ensured + // even p may have free variables. + // the proper side condition on p + 1 is -1 > p or -2 >= p or p + 1 != 0 + // I am disabling match_non_max and match_non_zero from forbidden_interval + // The narrowing rules in ule_constraint already handle the bounds propagaitons + // as it propagates p != -1 and 0 != q (p < -1, q > 0), + // + for (auto const& c : get_constraints(v)) { s.try_assign_eval(c); } @@ -299,8 +324,14 @@ namespace polysat { if (!e) return true; entry const* first = e; - rational const& max_value = s.var2pdd(v).max_value(); - rational mod_value = max_value + 1; + auto& m = s.var2pdd(v); + unsigned const N = m.power_of_2(); + rational const& max_value = m.max_value(); + rational const& mod_value = m.two_to_N(); + + // Rotate the 'first' entry, to prevent getting stuck in a refinement loop + // with an early entry when a later entry could give a better interval. + m_equal_lin[v] = m_equal_lin[v]->next(); auto delta_l = [&](rational const& coeff_val) { return floor((coeff_val - e->interval.lo_val()) / e->coeff); @@ -338,22 +369,79 @@ namespace polysat { } }; do { - LOG("refine-equal-lin for src: " << e->src); rational coeff_val = mod(e->coeff * val, mod_value); if (e->interval.currently_contains(coeff_val)) { + LOG("refine-equal-lin for src: " << lit_pp(s, e->src)); + LOG("forbidden interval v" << v << " " << num_pp(s, v, val) << " " << num_pp(s, v, e->coeff, true) << " * " << e->interval); - if (e->interval.lo_val().is_one() && e->interval.hi_val().is_zero() && e->coeff.is_odd()) { - rational lo(1); - rational hi(0); - LOG("refine-equal-lin: " << " [" << lo << ", " << hi << "["); - pdd lop = s.var2pdd(v).mk_val(lo); - pdd hip = s.var2pdd(v).mk_val(hi); + if (mod(e->interval.hi_val() + 1, mod_value) == e->interval.lo_val()) { + // We have an equation: a * v == b + rational const a = e->coeff; + rational const b = e->interval.hi_val(); + LOG("refine-equal-lin: equation detected: " << dd::val_pp(m, a, true) << " * v" << v << " == " << dd::val_pp(m, b, false)); + unsigned const parity_a = get_parity(a, N); + unsigned const parity_b = get_parity(b, N); + // LOG("a " << a << " parity " << parity_a); + // LOG("b " << b << " parity " << parity_b); + if (parity_a > parity_b) { + // No solution + LOG("refined: no solution due to parity"); + entry* ne = alloc_entry(); + ne->refined = e; + ne->src = e->src; + ne->side_cond = e->side_cond; + ne->coeff = 1; + ne->interval = eval_interval::full(); + intersect(v, ne); + return false; + } + if (parity_a == 0) { + // "fast path" for odd a + rational a_inv; + VERIFY(a.mult_inverse(N, a_inv)); + rational const hi = mod(a_inv * b, mod_value); + rational const lo = mod(hi + 1, mod_value); + LOG("refined to [" << num_pp(s, v, lo) << ", " << num_pp(s, v, hi) << "["); + SASSERT_EQ(mod(a * hi, mod_value), b); // hi is the solution + entry* ne = alloc_entry(); + ne->refined = e; + ne->src = e->src; + ne->side_cond = e->side_cond; + ne->coeff = 1; + ne->interval = eval_interval::proper(m.mk_val(lo), lo, m.mk_val(hi), hi); + SASSERT(ne->interval.currently_contains(val)); + intersect(v, ne); + return false; + } + // 2^k * v == a_inv * b + // 2^k solutions because only the lower N-k bits of v are fixed. + // + // Smallest solution is v0 == a_inv * (b >> k) + // Solutions are of the form v_i = v0 + 2^(N-k) * i for i in { 0, 1, ..., 2^k - 1 }. + // Forbidden intervals: [v_i + 1; v_{i+1}[ == [ v_i + 1; v_i + 2^(N-k) [ + // We need the interval that covers val: + // v_i + 1 <= val < v_i + 2^(N-k) + // + // TODO: create one interval for v[N-k:] instead of 2^k intervals for v. + unsigned const k = parity_a; + rational const a_inv = a.pseudo_inverse(N); + unsigned const N_minus_k = N - k; + rational const two_to_N_minus_k = rational::power_of_two(N_minus_k); + rational const v0 = mod(a_inv * machine_div2k(b, k), two_to_N_minus_k); + SASSERT(mod(val, two_to_N_minus_k) != v0); // val is not a solution + rational const vi = v0 + clear_lower_bits(mod(val - v0, mod_value), N_minus_k); + rational const lo = mod(vi + 1, mod_value); + rational const hi = mod(vi + two_to_N_minus_k, mod_value); + LOG("refined to [" << num_pp(s, v, lo) << ", " << num_pp(s, v, hi) << "["); + SASSERT_EQ(mod(a * (lo - 1), mod_value), b); // lo-1 is a solution + SASSERT_EQ(mod(a * hi, mod_value), b); // hi is a solution entry* ne = alloc_entry(); ne->refined = e; ne->src = e->src; ne->side_cond = e->side_cond; ne->coeff = 1; - ne->interval = eval_interval::proper(lop, lo, hip, hi); + ne->interval = eval_interval::proper(m.mk_val(lo), lo, m.mk_val(hi), hi); + SASSERT(ne->interval.currently_contains(val)); intersect(v, ne); return false; } @@ -378,13 +466,11 @@ namespace polysat { lo = val - lambda_l; increase_hi(hi); } - LOG("forbidden interval v" << v << " " << num_pp(s, v, val) << " " << num_pp(s, v, e->coeff, true) << " * " << e->interval << " [" << num_pp(s, v, lo) << ", " << num_pp(s, v, hi) << "["); + LOG("refined to [" << num_pp(s, v, lo) << ", " << num_pp(s, v, hi) << "["); SASSERT(hi <= mod_value); bool full = (lo == 0 && hi == mod_value); if (hi == mod_value) hi = 0; - pdd lop = s.var2pdd(v).mk_val(lo); - pdd hip = s.var2pdd(v).mk_val(hi); entry* ne = alloc_entry(); ne->refined = e; ne->src = e->src; @@ -393,7 +479,7 @@ namespace polysat { if (full) ne->interval = eval_interval::full(); else - ne->interval = eval_interval::proper(lop, lo, hip, hi); + ne->interval = eval_interval::proper(m.mk_val(lo), lo, m.mk_val(hi), hi); intersect(v, ne); return false; } @@ -412,6 +498,10 @@ namespace polysat { rational const& max_value = s.var2pdd(v).max_value(); rational const mod_value = max_value + 1; + // Rotate the 'first' entry, to prevent getting stuck in a refinement loop + // with an early entry when a later entry could give a better interval. + m_diseq_lin[v] = m_diseq_lin[v]->next(); + do { LOG("refine-disequal-lin for src: " << e->src); // We compute an interval if the concrete value 'val' violates the constraint: @@ -543,88 +633,146 @@ namespace polysat { return refine_viable(v, val); } - - rational viable::min_viable(pvar v) { - refined: - rational lo(0); - auto* e = m_units[v]; - if (!e && !refine_viable(v, lo)) - goto refined; - if (!e) - return lo; - entry* first = e; - entry* last = first->prev(); - if (last->interval.currently_contains(lo)) - lo = last->interval.hi_val(); - do { - if (!e->interval.currently_contains(lo)) - break; - lo = e->interval.hi_val(); - e = e->next(); + find_t viable::find_viable(pvar v, rational& lo) { + rational hi; + switch (find_viable(v, lo, hi)) { + case l_true: + return (lo == hi) ? find_t::singleton : find_t::multiple; + case l_false: + return find_t::empty; + default: + return find_t::resource_out; } - while (e != first); - if (!refine_viable(v, lo)) - goto refined; - SASSERT(is_viable(v, lo)); - return lo; } - rational viable::max_viable(pvar v) { - refined: - rational hi = s.var2pdd(v).max_value(); - auto* e = m_units[v]; - if (!e && !refine_viable(v, hi)) - goto refined; - if (!e) - return hi; - entry* last = e->prev(); - e = last; - do { - if (!e->interval.currently_contains(hi)) - break; - hi = e->interval.lo_val() - 1; - e = e->prev(); - } - while (e != last); - if (!refine_viable(v, hi)) - goto refined; - SASSERT(is_viable(v, hi)); - return hi; + lbool viable::find_viable(pvar v, rational& lo, rational& hi) { + std::pair args{lo, hi}; + return query(v, args); } - // template - find_t viable::query(query_t mode, pvar v, rational& lo, rational& hi) { - SASSERT(mode == query_t::find_viable); // other modes are TODO + lbool viable::min_viable(pvar v, rational& lo) { + return query(v, lo); + } - auto const& max_value = s.var2pdd(v).max_value(); + lbool viable::max_viable(pvar v, rational& hi) { + return query(v, hi); + } + bool viable::has_upper_bound(pvar v, rational& out_hi, vector& out_c) { + entry const* first = m_units[v]; + entry const* e = first; + bool found = false; + out_c.reset(); + do { + found = false; + do { + if (!e->refined) { + auto const& lo = e->interval.lo(); + auto const& hi = e->interval.hi(); + if (lo.is_val() && hi.is_val()) { + if (out_c.empty() && lo.val() > hi.val()) { + out_c.push_back(e->src); + out_hi = lo.val() - 1; + found = true; + } + else if (!out_c.empty() && lo.val() <= out_hi && out_hi < hi.val()) { + out_c.push_back(e->src); + out_hi = lo.val() - 1; + found = true; + } + } + } + e = e->next(); + } + while (e != first); + } + while (found); + return !out_c.empty(); + } + + bool viable::has_lower_bound(pvar v, rational& out_lo, vector& out_c) { + entry const* first = m_units[v]; + entry const* e = first; + bool found = false; + out_c.reset(); + do { + found = false; + do { + if (!e->refined) { + auto const& lo = e->interval.lo(); + auto const& hi = e->interval.hi(); + if (lo.is_val() && hi.is_val()) { + if (out_c.empty() && hi.val() != 0 && (lo.val() == 0 || lo.val() > hi.val())) { + out_c.push_back(e->src); + out_lo = hi.val(); + found = true; + } + else if (!out_c.empty() && lo.val() <= out_lo && out_lo < hi.val()) { + out_c.push_back(e->src); + out_lo = hi.val(); + found = true; + } + } + } + e = e->next(); + } + while (e != first); + } + while (found); + return !out_c.empty(); + } + + template + lbool viable::query(pvar v, typename query_result::result_t& result) { // max number of interval refinements before falling back to the univariate solver unsigned const refinement_budget = 1000; unsigned refinements = refinement_budget; - refined: + while (refinements--) { + lbool res = l_undef; - if (!refinements) { - LOG("Refinement budget exhausted! Fall back to univariate solver."); - return find_t::resource_out; + if constexpr (mode == query_t::find_viable) + res = query_find(v, result.first, result.second); + else if constexpr (mode == query_t::min_viable) + res = query_min(v, result); + else if constexpr (mode == query_t::max_viable) + res = query_max(v, result); + else if constexpr (mode == query_t::has_viable) { + NOT_IMPLEMENTED_YET(); + } + else { + UNREACHABLE(); + } + + if (res != l_undef) + return res; } - refinements--; + LOG("Refinement budget exhausted! Fall back to univariate solver."); + return query_fallback(v, result); + } + + lbool viable::query_find(pvar v, rational& lo, rational& hi) { + auto const& max_value = s.var2pdd(v).max_value(); + lbool const refined = l_undef; // After a refinement, any of the existing entries may have been replaced // (if it is subsumed by the new entry created during refinement). // For this reason, we start chasing the intervals from the start again. lo = 0; + hi = max_value; auto* e = m_units[v]; if (!e && !refine_viable(v, lo)) - goto refined; - if (!e && !refine_viable(v, rational::one())) - goto refined; + return refined; + if (!e && !refine_viable(v, hi)) + return refined; if (!e) - return find_t::multiple; - if (e->interval.is_full()) - return find_t::empty; + return l_true; + if (e->interval.is_full()) { + s.set_conflict_by_viable_interval(v); + return l_false; + } entry* first = e; entry* last = first->prev(); @@ -635,10 +783,10 @@ namespace polysat { last->interval.hi_val() < max_value) { lo = last->interval.hi_val(); if (!refine_viable(v, lo)) - goto refined; + return refined; if (!refine_viable(v, max_value)) - goto refined; - return find_t::multiple; + return refined; + return l_true; } // find lower bound @@ -652,8 +800,10 @@ namespace polysat { } while (e != first); - if (e->interval.currently_contains(lo)) - return find_t::empty; + if (e->interval.currently_contains(lo)) { + s.set_conflict_by_viable_interval(v); + return l_false; + } // find upper bound hi = max_value; @@ -666,51 +816,22 @@ namespace polysat { } while (e != last); if (!refine_viable(v, lo)) - goto refined; + return refined; if (!refine_viable(v, hi)) - goto refined; - - if (lo == hi) - return find_t::singleton; - else - return find_t::multiple; + return refined; + return l_true; } - find_t viable::find_viable(pvar v, rational& lo) { -#if 1 - rational hi; - // return query(v, lo, hi); - return query(query_t::find_viable, v, lo, hi); -#else - refined: + lbool viable::query_min(pvar v, rational& lo) { + // TODO: should be able to deal with UNSAT case; since also min_viable has to deal with it due to fallback solver lo = 0; - auto* e = m_units[v]; + entry* e = m_units[v]; if (!e && !refine_viable(v, lo)) - goto refined; - if (!e && !refine_viable(v, rational::one())) - goto refined; + return l_undef; if (!e) - return find_t::multiple; - if (e->interval.is_full()) - return find_t::empty; - + return l_true; entry* first = e; entry* last = first->prev(); - - // quick check: last interval does not wrap around - // and has space for 2 unassigned values. - auto& max_value = s.var2pdd(v).max_value(); - if (last->interval.lo_val() < last->interval.hi_val() && - last->interval.hi_val() < max_value) { - lo = last->interval.hi_val(); - if (!refine_viable(v, lo)) - goto refined; - if (!refine_viable(v, max_value)) - goto refined; - return find_t::multiple; - } - - // find lower bound if (last->interval.currently_contains(lo)) lo = last->interval.hi_val(); do { @@ -720,12 +841,21 @@ namespace polysat { e = e->next(); } while (e != first); + if (!refine_viable(v, lo)) + return l_undef; + SASSERT(is_viable(v, lo)); + return l_true; + } - if (e->interval.currently_contains(lo)) - return find_t::empty; - - // find upper bound - rational hi = max_value; + lbool viable::query_max(pvar v, rational& hi) { + // TODO: should be able to deal with UNSAT case; since also max_viable has to deal with it due to fallback solver + hi = s.var2pdd(v).max_value(); + auto* e = m_units[v]; + if (!e && !refine_viable(v, hi)) + return l_undef; + if (!e) + return l_true; + entry* last = e->prev(); e = last; do { if (!e->interval.currently_contains(hi)) @@ -734,18 +864,126 @@ namespace polysat { e = e->prev(); } while (e != last); - if (!refine_viable(v, lo)) - goto refined; if (!refine_viable(v, hi)) - goto refined; - if (lo == hi) - return find_t::singleton; - else - return find_t::multiple; -#endif + return l_undef; + SASSERT(is_viable(v, hi)); + return l_true; } - bool viable::resolve(pvar v, conflict& core) { + template + lbool viable::query_fallback(pvar v, typename query_result::result_t& result) { + unsigned const bit_width = s.size(v); + univariate_solver* us = s.m_viable_fallback.usolver(bit_width); + sat::literal_set added; + + // First step: only query the looping constraints and see if they alone are already UNSAT. + // The constraints which caused the refinement loop will be reached from m_units. + LOG_H3("Checking looping univariate constraints for v" << v << "..."); + LOG("Assignment: " << assignments_pp(s)); + entry const* first = m_units[v]; + entry const* e = first; + do { + entry const* origin = e; + while (origin->refined) + origin = origin->refined; + signed_constraint const c = origin->src; + sat::literal const lit = c.blit(); + if (!added.contains(lit)) { + added.insert(lit); + LOG("Adding " << lit_pp(s, lit)); + c.add_to_univariate_solver(v, s, *us, lit.to_uint()); + } + e = e->next(); + } + while (e != first); + + switch (us->check()) { + case l_false: + s.set_conflict_by_viable_fallback(v, *us); + return l_false; + case l_true: + // At this point we don't know much because we did not add all relevant constraints + break; + default: + // resource limit + return l_undef; + } + + // Second step: looping constraints aren't UNSAT, so add the remaining relevant constraints + LOG_H3("Checking all univariate constraints for v" << v << "..."); + auto const& cs = s.m_viable_fallback.m_constraints[v]; + for (unsigned i = cs.size(); i-- > 0; ) { + sat::literal const lit = cs[i].blit(); + if (added.contains(lit)) + continue; + LOG("Adding " << lit_pp(s, lit)); + added.insert(lit); + cs[i].add_to_univariate_solver(v, s, *us, lit.to_uint()); + } + + switch (us->check()) { + case l_false: + s.set_conflict_by_viable_fallback(v, *us); + return l_false; + case l_true: + // pass solver to mode-specific query + break; + default: + // resource limit + return l_undef; + } + + if constexpr (mode == query_t::find_viable) + return query_find_fallback(v, *us, result.first, result.second); + + if constexpr (mode == query_t::min_viable) + return query_min_fallback(v, *us, result); + + if constexpr (mode == query_t::max_viable) + return query_max_fallback(v, *us, result); + + if constexpr (mode == query_t::has_viable) { + NOT_IMPLEMENTED_YET(); + return l_undef; + } + + UNREACHABLE(); + return l_undef; + } + + lbool viable::query_find_fallback(pvar v, univariate_solver& us, rational& lo, rational& hi) { + return us.find_two(lo, hi) ? l_true : l_undef; + } + + lbool viable::query_min_fallback(pvar v, univariate_solver& us, rational& lo) { + return us.find_min(lo) ? l_true : l_undef; + } + + lbool viable::query_max_fallback(pvar v, univariate_solver& us, rational& hi) { + return us.find_max(hi) ? l_true : l_undef; + } + + bool viable::resolve_fallback(pvar v, univariate_solver& us, conflict& core) { + // The conflict is the unsat core of the univariate solver, + // and the current assignment (under which the constraints are univariate in v) + // TODO: + // - currently we add variables directly, which is sound: + // e.g.: v^2 + w^2 == 0; w := 1 + // - but we could use side constraints on the coefficients instead (coefficients when viewed as polynomial over v): + // e.g.: v^2 + w^2 == 0; w^2 == 1 + for (unsigned dep : us.unsat_core()) { + sat::literal lit = sat::to_literal(dep); + signed_constraint c = s.lit2cnstr(lit); + core.insert(c); + core.insert_vars(c); + } + SASSERT(!core.vars().contains(v)); + core.add_lemma("viable unsat core", core.build_lemma()); + verbose_stream() << "unsat core " << core << "\n"; + return true; + } + + bool viable::resolve_interval(pvar v, conflict& core) { DEBUG_CODE( log(v); ); if (has_viable(v)) return false; @@ -810,7 +1048,7 @@ namespace polysat { auto lhs = hi - next_lo; auto rhs = next_hi - next_lo; signed_constraint c = s.m_constraints.ult(lhs, rhs); - lemma.insert_eval(~c); + lemma.insert_try_eval(~c); // "try" because linking constraint may contain unassigned variables, see test_polysat::test_bench23_fi_lemma for an example. } for (auto sc : e->side_cond) lemma.insert_eval(~sc); @@ -821,8 +1059,12 @@ namespace polysat { } while (e != first); - SASSERT(all_of(lemma, [this](sat::literal lit) { return s.m_bvars.value(lit) == l_false || s.lit2cnstr(lit).is_currently_false(s); })); + // Doesn't hold anymore: we may get new constraints with unassigned variables, see test_polysat::test_bench23_fi_lemma. + // SASSERT(all_of(lemma, [this](sat::literal lit) { return s.m_bvars.value(lit) == l_false || s.lit2cnstr(lit).is_currently_false(s); })); + // NSB review: bench23 exposes a scenario where s.m_bvars.value(lit) == l_true. So the viable lemma is mute, but the literal in the premise + // is a conflict. + // SASSERT(all_of(lemma, [this](sat::literal lit) { return s.m_bvars.value(lit) != l_true; })); core.add_lemma("viable", lemma.build()); core.logger().log(inf_fi(*this, v)); return true; @@ -947,27 +1189,36 @@ namespace polysat { return {}; } - find_t viable_fallback::find_viable(pvar v, rational& out_val) { - unsigned bit_width = s.m_size[v]; - + univariate_solver* viable_fallback::usolver(unsigned bit_width) { univariate_solver* us; + auto it = m_usolver.find_iterator(bit_width); if (it != m_usolver.end()) { us = it->m_value.get(); us->pop(1); - } else { + } + else { auto& mk_solver = *m_usolver_factory; m_usolver.insert(bit_width, mk_solver(bit_width)); us = m_usolver[bit_width].get(); } + SASSERT_EQ(us->scope_level(), 0); // push once on the empty solver so we can reset it before the next use us->push(); + return us; + } + + find_t viable_fallback::find_viable(pvar v, rational& out_val) { + unsigned const bit_width = s.m_size[v]; + univariate_solver* us = usolver(bit_width); + auto const& cs = m_constraints[v]; for (unsigned i = cs.size(); i-- > 0; ) { - LOG("Univariate constraint: " << cs[i]); - cs[i].add_to_univariate_solver(s, *us, i); + signed_constraint const c = cs[i]; + LOG("Univariate constraint: " << c); + c.add_to_univariate_solver(v, s, *us, c.blit().to_uint()); } switch (us->check()) { @@ -976,22 +1227,13 @@ namespace polysat { // we don't know whether the SMT instance has a unique solution return find_t::multiple; case l_false: + s.set_conflict_by_viable_fallback(v, *us); return find_t::empty; default: return find_t::resource_out; } } - signed_constraints viable_fallback::unsat_core(pvar v) { - unsigned bit_width = s.m_size[v]; - SASSERT(m_usolver[bit_width]); - signed_constraints cs; - for (unsigned dep : m_usolver[bit_width]->unsat_core()) { - cs.push_back(m_constraints[v][dep]); - } - return cs; - } - std::ostream& operator<<(std::ostream& out, find_t x) { switch (x) { case find_t::empty: diff --git a/src/math/polysat/viable.h b/src/math/polysat/viable.h index 6588dcc18..e677bb550 100644 --- a/src/math/polysat/viable.h +++ b/src/math/polysat/viable.h @@ -25,6 +25,7 @@ Author: #include "math/polysat/conflict.h" #include "math/polysat/constraint.h" #include "math/polysat/forbidden_intervals.h" +#include namespace polysat { @@ -39,6 +40,34 @@ namespace polysat { resource_out, }; + namespace viable_query { + enum class query_t { + has_viable, // currently only used internally in resolve_viable + min_viable, // currently unused + max_viable, // currently unused + find_viable, + }; + + template + struct query_result { + }; + + template <> + struct query_result { + using result_t = rational; + }; + + template <> + struct query_result { + using result_t = rational; + }; + + template <> + struct query_result { + using result_t = std::pair; + }; + } + std::ostream& operator<<(std::ostream& out, find_t x); class viable { @@ -76,15 +105,36 @@ namespace polysat { void propagate(pvar v, rational const& val); - enum class query_t { - has_viable, // currently only used internally in resolve_viable - min_viable, // currently unused - max_viable, // currently unused - find_viable, - }; + /** + * Interval-based queries + * @return l_true on success, l_false on conflict, l_undef on refinement + */ + lbool query_min(pvar v, rational& out_lo); + lbool query_max(pvar v, rational& out_hi); + lbool query_find(pvar v, rational& out_lo, rational& out_hi); - // template - find_t query(query_t mode, pvar v, rational& out_lo, rational& out_hi); + /** + * Bitblasting-based queries. + * The univariate solver has already been filled with all relevant constraints and check() returned l_true. + * @return l_true on success, l_false on conflict, l_undef on resource limit + */ + lbool query_min_fallback(pvar v, univariate_solver& us, rational& out_lo); + lbool query_max_fallback(pvar v, univariate_solver& us, rational& out_hi); + lbool query_find_fallback(pvar v, univariate_solver& us, rational& out_lo, rational& out_hi); + + /** + * Interval-based query with bounded refinement and fallback to bitblasting. + * @return l_true on success, l_false on conflict, l_undef on resource limit + */ + template + lbool query(pvar v, typename viable_query::query_result::result_t& out_result); + + /** + * Bitblasting-based query. + * @return l_true on success, l_false on conflict, l_undef on resource limit + */ + template + lbool query_fallback(pvar v, typename viable_query::query_result::result_t& out_result); public: viable(solver& s); @@ -122,23 +172,54 @@ namespace polysat { */ bool is_viable(pvar v, rational const& val); - /* - * Extract min and max viable values for v + /** + * Extract min viable value for v. + * @return l_true on success, l_false on conflict, l_undef on resource limit */ - rational min_viable(pvar v); - rational max_viable(pvar v); + lbool min_viable(pvar v, rational& out_lo); + + /** + * Extract max viable value for v. + * @return l_true on success, l_false on conflict, l_undef on resource limit + */ + lbool max_viable(pvar v, rational& out_hi); + + + /** + * Query for an upper bound literal for v together with justification. + * @return true if a non-trivial upper bound is found, return justifying constraint. + */ + bool has_upper_bound(pvar v, rational& out_hi, vector& out_c); + + /** + * Query for an lower bound literal for v together with justification. + * @return true if a non-trivial lower bound is found, return justifying constraint. + */ + bool has_lower_bound(pvar v, rational& out_lo, vector& out_c); /** * Find a next viable value for variable. */ find_t find_viable(pvar v, rational& out_val); + /** + * Find a next viable value for variable. Attempts to find two different values, to distinguish propagation/decision. + * @return l_true on success, l_false on conflict, l_undef on resource limit + */ + lbool find_viable(pvar v, rational& out_lo, rational& out_hi); + /** * Retrieve the unsat core for v, * and add the forbidden interval lemma for v (which eliminates v from the unsat core). - * \pre there are no viable values for v + * \pre there are no viable values for v (determined by interval reasoning) */ - bool resolve(pvar v, conflict& core); + bool resolve_interval(pvar v, conflict& core); + + /** + * Retrieve the unsat core for v. + * \pre there are no viable values for v (determined by fallback solver) + */ + bool resolve_fallback(pvar v, univariate_solver& us, conflict& core); /** Log all viable values for the given variable. * (Inefficient, but useful for debugging small instances.) @@ -251,6 +332,8 @@ namespace polysat { } class viable_fallback { + friend class viable; + solver& s; scoped_ptr m_usolver_factory; @@ -258,6 +341,8 @@ namespace polysat { vector m_constraints; svector m_constraints_trail; + univariate_solver* usolver(unsigned bit_width); + public: viable_fallback(solver& s); @@ -275,7 +360,6 @@ namespace polysat { signed_constraint find_violated_constraint(assignment const& a, pvar v); find_t find_viable(pvar v, rational& out_val); - signed_constraints unsat_core(pvar v); }; } diff --git a/src/test/polysat.cpp b/src/test/polysat.cpp index d36a60b9f..366679895 100644 --- a/src/test/polysat.cpp +++ b/src/test/polysat.cpp @@ -626,8 +626,46 @@ namespace polysat { cb.insert(s.ule(u, v)); auto cl = cb.build(); simp.apply(*cl); - std::cout << *cl << "\n"; - SASSERT(cl->size() == 2); + VERIFY_EQ(cl->size(), 2); + } + + // p <= q + // p == q (should be removed) + static void test_clause_simplify2() { + scoped_solver s(__func__); + simplify_clause simp(s); + clause_builder cb(s); + auto u = s.var(s.add_var(32)); + auto v = s.var(s.add_var(32)); + auto w = s.var(s.add_var(32)); + auto p = 2*u*v; + auto q = 7*v*w; + cb.insert(s.ule(p, q)); + cb.insert(s.eq(p, q)); + auto cl = cb.build(); + simp.apply(*cl); + VERIFY_EQ(cl->size(), 1); + VERIFY_EQ((*cl)[0], s.ule(p, q).blit()); + } + + // p < q + // p == q + // should be merged into p <= q + static void test_clause_simplify3() { + scoped_solver s(__func__); + simplify_clause simp(s); + clause_builder cb(s); + auto u = s.var(s.add_var(32)); + auto v = s.var(s.add_var(32)); + auto w = s.var(s.add_var(32)); + auto p = 2*u*v; + auto q = 7*v*w; + cb.insert(s.ult(p, q)); + cb.insert(s.eq(p, q)); + auto cl = cb.build(); + simp.apply(*cl); + VERIFY_EQ(cl->size(), 1); + VERIFY_EQ((*cl)[0], s.ule(p, q).blit()); } // 8 * x + 3 == 0 or 8 * x + 5 == 0 is unsat @@ -1099,6 +1137,7 @@ namespace polysat { if (lemmas.size() == cnt) return; LOG_H1("FAIL: Expected " << cnt << " learned lemmas; got " << lemmas.size() << "!"); + SASSERT(false); if (!collect_test_records) VERIFY(false); } @@ -1115,6 +1154,7 @@ namespace polysat { } } LOG_H1("FAIL: Lemma " << c1 << " not deduced!"); + SASSERT(false); if (!collect_test_records) VERIFY(false); } @@ -1361,7 +1401,7 @@ namespace polysat { static void test_ineq_non_axiom1(unsigned bw, unsigned i) { auto const bound = rational::power_of_two(bw - 1); - scoped_solver s(std::string(__func__) + " perm=" + std::to_string(i)); + scoped_solver s(concat(__func__, " bw=", bw, " perm=", i)); auto x = s.var(s.add_var(bw)); auto y = s.var(s.add_var(bw)); auto z = s.var(s.add_var(bw)); @@ -1383,7 +1423,7 @@ namespace polysat { static void test_ineq_axiom2(unsigned bw = 32) { auto const bound = rational::power_of_two(bw/2); for (unsigned i = 0; i < 6; ++i) { - scoped_solver s(std::string(__func__) + " perm=" + std::to_string(i)); + scoped_solver s(concat(__func__, " bw=", bw, " perm=", i)); auto x = s.var(s.add_var(bw)); auto y = s.var(s.add_var(bw)); auto z = s.var(s.add_var(bw)); @@ -1399,30 +1439,33 @@ namespace polysat { } // xy < b & a <= x & !Omega(x*y) => a*y < b - static void test_ineq_axiom3(unsigned bw = 32) { + static void test_ineq_axiom3(unsigned bw, unsigned i) { auto const bound = rational::power_of_two(bw/2); - for (unsigned i = 0; i < 24; ++i) { - scoped_solver s(std::string(__func__) + " perm=" + std::to_string(i)); - auto x = s.var(s.add_var(bw)); - auto y = s.var(s.add_var(bw)); - auto a = s.var(s.add_var(bw)); - auto b = s.var(s.add_var(bw)); - permute_args(i, x, y, a, b); - s.add_ult(x * y, b); - s.add_ule(a, x); - s.add_ult(x, bound); - s.add_ult(y, bound); - s.add_ule(b, a * y); - s.check(); - s.expect_unsat(); - } + scoped_solver s(concat(__func__, " bw=", bw, " perm=", i)); + auto x = s.var(s.add_var(bw)); + auto y = s.var(s.add_var(bw)); + auto a = s.var(s.add_var(bw)); + auto b = s.var(s.add_var(bw)); + permute_args(i, x, y, a, b); + s.add_ult(x * y, b); + s.add_ule(a, x); + s.add_ult(x, bound); + s.add_ult(y, bound); + s.add_ule(b, a * y); + s.check(); + s.expect_unsat(); + } + + static void test_ineq_axiom3(unsigned bw = 32) { + for (unsigned i = 0; i < 24; ++i) + test_ineq_axiom3(bw, i); } // x*y <= b & a <= x & !Omega(x*y) => a*y <= b static void test_ineq_axiom4(unsigned bw = 32) { auto const bound = rational::power_of_two(bw/2); for (unsigned i = 0; i < 24; ++i) { - scoped_solver s(std::string(__func__) + " perm=" + std::to_string(i)); + scoped_solver s(concat(__func__, " bw=", bw, " perm=", i)); auto x = s.var(s.add_var(bw)); auto y = s.var(s.add_var(bw)); auto a = s.var(s.add_var(bw)); @@ -1465,7 +1508,7 @@ namespace polysat { // a < xy & x <= b & !Omega(b*y) => a < b*y static void test_ineq_axiom5(unsigned bw, unsigned i) { auto const bound = rational::power_of_two(bw/2); - scoped_solver s(std::string(__func__) + " perm=" + std::to_string(i)); + scoped_solver s(concat(__func__, " bw=", bw, " perm=", i)); auto x = s.var(s.add_var(bw)); auto y = s.var(s.add_var(bw)); auto a = s.var(s.add_var(bw)); @@ -1488,7 +1531,7 @@ namespace polysat { // a <= xy & x <= b & !Omega(b*y) => a <= b*y static void test_ineq_axiom6(unsigned bw, unsigned i) { auto const bound = rational::power_of_two(bw/2); - scoped_solver s(std::string(__func__) + " perm=" + std::to_string(i)); + scoped_solver s(concat(__func__, " bw=", bw, " perm=", i)); auto x = s.var(s.add_var(bw)); auto y = s.var(s.add_var(bw)); auto a = s.var(s.add_var(bw)); @@ -1757,6 +1800,29 @@ namespace polysat { s.expect_sat({{a, 5}}); } + /** + * Motivated by weak forbidden intervals lemma in bench23: + * v66 > v41 && v67 == v66 ==> v67 != 0 + * + * Cause by omitting v41 from the symbolic bounds of v66 > v41. + * + * The stronger lemma should be: + * v66 > v41 && v67 == v66 ==> v41 + 1 <= v67 + * + * The conclusion could also be v67 > v41 (note that -1 > v41 follows from premise v66 > v41, so both conclusions are equivalent in this lemma). + */ + static void test_bench23_fi_lemma() { + scoped_solver s(__func__); + auto x = s.var(s.add_var(256)); // v41 + auto y = s.var(s.add_var(256)); // v66 + auto z = s.var(s.add_var(256)); // v67 + s.add_ult(x, y); // v66 > v41 + s.add_eq(y, z); // v66 == v67 + s.add_eq(z, 0); // v67 == 0 + s.check(); + s.expect_unsat(); + } + static void test_bench6_viable() { scoped_solver s(__func__); rational coeff("12737129816104781496592808350955669863859698313220462044340334240870444260393"); @@ -1768,6 +1834,66 @@ namespace polysat { s.check(); } + static void test_bench27_viable(const char* coeff) { + scoped_solver s(__func__); + rational a(coeff); + rational b = rational::power_of_two(128) - 1; + auto x = s.var(s.add_var(256)); + s.add_ult(0, x); + s.add_ule(a * x + b, b); + s.check(); + } + + static void test_bench27_viable1() { + test_bench27_viable("2787252867449406954228501409473613420036128"); // 2^5 * a' + } + + static void test_bench27_viable2() { + test_bench27_viable("5574846017265734846920466193554658608284160"); // 2^9 * a' + } + + // -1*v12 <= -1*v12 + v8 + v7*v1 + 2^128*v7 + 1 + // v8 := 0 v12 := 1 v4 := 0 v9 := 1 v17 := 0 v1 := 5574846017265734846920466193554658608283648 + // + // Substitute assignment: + // -1 <= (5574846017265734846920466193554658608283648 + 2^128) * v7 + // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ == 2^142 + // this is actually an equation + // + // 2^142 * v7 == -1, unsatisfiable due to parity + static void test_bench27_viable3() { + scoped_solver s(__func__); + rational const a("5574846017265734846920466193554658608283648"); + rational const b = rational::power_of_two(128); + auto const v1 = s.var(s.add_var(256)); + auto const v7 = s.var(s.add_var(256)); + auto const v8 = s.var(s.add_var(256)); + auto const v12 = s.var(s.add_var(256)); + s.add_ule(-v12, -v12 + v8 + v7*v1 + b*v7 + 1); + s.add_eq(v8, 0); + s.add_eq(v12, 1); + s.add_eq(v1, a); + s.check(); + s.expect_unsat(); + } + + // similar as test_bench27_viable3, but satisfiable (to test the other branch) + static void test_bench27_viable3_sat() { + scoped_solver s(__func__); + rational const a("5574846017265734846920466193554658608283648"); + rational const b = rational::power_of_two(128) - 1; + auto const v1 = s.var(s.add_var(256)); + auto const v7 = s.var(s.add_var(256)); + auto const v8 = s.var(s.add_var(256)); + auto const v12 = s.var(s.add_var(256)); + s.add_ule(-v12, -v12 + v8 + v7*v1 + b*v7 + 1); + s.add_eq(v8, 0); + s.add_eq(v12, 1); + s.add_eq(v1, a); + s.check(); + s.expect_sat(); + } + }; // class test_polysat @@ -1898,23 +2024,15 @@ static void STD_CALL polysat_on_ctrl_c(int) { void tst_polysat() { using namespace polysat; - polysat::test_polysat::test_band6(4); #if 0 // Enable this block to run a single unit test with detailed output. collect_test_records = false; test_max_conflicts = 50; - // test_polysat::test_parity1(); - // test_polysat::test_parity1b(); - // test_polysat::test_parity2(); - // test_polysat::test_parity3(); - // test_polysat::test_parity4(); - // test_polysat::test_parity4(32); - test_polysat::test_parity4(256); - // test_polysat::test_ineq2(); - // test_polysat::test_ineq_non_axiom4(32, 3); // stuck in viable refinement loop - // test_polysat::test_ineq_non_axiom4(32, 4); // stuck in viable refinement loop - // test_polysat::test_band1(); - // test_polysat::test_bench6_viable(); + test_polysat::test_ineq_axiom3(32, 3); // TODO: assertion + // test_polysat::test_ineq_axiom6(32, 0); // TODO: assertion + // test_polysat::test_band5(); // TODO: assertion when clause simplification (merging p>q and p=q) is enabled + // test_polysat::test_bench27_viable1(); // TODO: refinement + // test_polysat::test_bench27_viable2(); // TODO: refinement return; #endif @@ -1930,6 +2048,16 @@ void tst_polysat() { set_log_enabled(false); } +#if 0 + RUN(test_polysat::test_elim1()); + RUN(test_polysat::test_elim2()); + RUN(test_polysat::test_elim3()); + RUN(test_polysat::test_elim4()); + RUN(test_polysat::test_elim5()); + RUN(test_polysat::test_elim6()); + RUN(test_polysat::test_elim7()); +#endif + RUN(test_polysat::test_parity1()); RUN(test_polysat::test_parity1b()); RUN(test_polysat::test_parity2()); @@ -1937,6 +2065,9 @@ void tst_polysat() { RUN(test_polysat::test_parity4()); RUN(test_polysat::test_clause_simplify1()); + RUN(test_polysat::test_clause_simplify2()); + // RUN(test_polysat::test_clause_simplify3()); // TODO: corresponding simplification is disabled at the moment + RUN(test_polysat::test_bench23_fi_lemma()); RUN(test_polysat::test_add_conflicts()); RUN(test_polysat::test_wlist()); @@ -1966,14 +2097,6 @@ void tst_polysat() { RUN(test_polysat::test_var_minimize()); // works but var_minimize isn't used (UNSAT before lemma is created) - RUN(test_polysat::test_elim1()); - RUN(test_polysat::test_elim2()); - RUN(test_polysat::test_elim3()); - RUN(test_polysat::test_elim4()); - RUN(test_polysat::test_elim5()); - RUN(test_polysat::test_elim6()); - RUN(test_polysat::test_elim7()); - RUN(test_polysat::test_ineq1()); RUN(test_polysat::test_ineq2()); RUN(test_polysat::test_monot()); @@ -2002,6 +2125,10 @@ void tst_polysat() { RUN(test_polysat::test_fi_nonmax()); RUN(test_polysat::test_fi_disequal_mild()); RUN(test_polysat::test_bench6_viable()); + // RUN(test_polysat::test_bench27_viable1()); // stuck in refinement loop + fallback solver + // RUN(test_polysat::test_bench27_viable2()); // stuck in refinement loop + fallback solver + RUN(test_polysat::test_bench27_viable3()); + RUN(test_polysat::test_bench27_viable3_sat()); RUN(test_polysat::test_ineq_axiom1()); RUN(test_polysat::test_ineq_axiom2()); @@ -2010,7 +2137,7 @@ void tst_polysat() { RUN(test_polysat::test_ineq_axiom5()); RUN(test_polysat::test_ineq_axiom6()); RUN(test_polysat::test_ineq_non_axiom1()); - //RUN(test_polysat::test_ineq_non_axiom4()); + RUN(test_polysat::test_ineq_non_axiom4()); // test_fi::exhaustive(); // test_fi::randomized(); diff --git a/src/test/viable.cpp b/src/test/viable.cpp index c2614854d..ba68e377b 100644 --- a/src/test/viable.cpp +++ b/src/test/viable.cpp @@ -38,16 +38,19 @@ namespace polysat { auto c = s.ule(x + 3, x + 5); s.v.intersect(xv, c); std::cout << s.v << "\n"; - std::cout << "min-max " << s.v.min_viable(xv) << " " << s.v.max_viable(xv) << "\n"; + std::cout << "min " << s.v.min_viable(xv, val) << " " << val << "\n"; + std::cout << "max " << s.v.max_viable(xv, val) << " " << val << "\n"; s.v.intersect(xv, s.ule(x, 2)); std::cout << s.v << "\n"; s.v.intersect(xv, s.ule(1, x)); std::cout << s.v << "\n"; - std::cout << "min-max " << s.v.min_viable(xv) << " " << s.v.max_viable(xv) << "\n"; + std::cout << "min " << s.v.min_viable(xv, val) << " " << val << "\n"; + std::cout << "max " << s.v.max_viable(xv, val) << " " << val << "\n"; s.v.intersect(xv, s.ule(x, 3)); std::cout << s.v << "\n"; std::cout << s.v.find_viable(xv, val) << " " << val << "\n"; - std::cout << "min-max " << s.v.min_viable(xv) << " " << s.v.max_viable(xv) << "\n"; + std::cout << "min " << s.v.min_viable(xv, val) << " " << val << "\n"; + std::cout << "max " << s.v.max_viable(xv, val) << " " << val << "\n"; s.v.intersect(xv, s.ule(3, x)); std::cout << s.v << "\n"; std::cout << s.v.find_viable(xv, val) << " " << val << "\n"; diff --git a/src/util/mpq.h b/src/util/mpq.h index f1afcbe50..7d5a6795a 100644 --- a/src/util/mpq.h +++ b/src/util/mpq.h @@ -490,6 +490,8 @@ public: void machine_div_rem(mpz const & a, mpz const & b, mpz & c, mpz & d) { mpz_manager::machine_div_rem(a, b, c, d); } + void machine_div2k(mpz const & a, unsigned k, mpz & c) { mpz_manager::machine_div2k(a, k, c); } + void div(mpz const & a, mpz const & b, mpz & c) { mpz_manager::div(a, b, c); } void rat_div(mpz const & a, mpz const & b, mpq & c) { @@ -516,6 +518,12 @@ public: machine_div(a.m_num, b.m_num, c); } + void machine_idiv2k(mpq const & a, unsigned k, mpq & c) { + SASSERT(is_int(a)); + machine_div2k(a.m_num, k, c.m_num); + reset_denominator(c); + } + void idiv(mpq const & a, mpq const & b, mpq & c) { SASSERT(is_int(a) && is_int(b)); div(a.m_num, b.m_num, c.m_num); diff --git a/src/util/rational.cpp b/src/util/rational.cpp index af3c89ced..531669dd4 100644 --- a/src/util/rational.cpp +++ b/src/util/rational.cpp @@ -153,3 +153,21 @@ bool rational::mult_inverse(unsigned num_bits, rational & result) const { return true; } +/** + * Compute multiplicative pseudo-inverse modulo 2^num_bits: + * + * mod(n * n.pseudo_inverse(bits), 2^bits) == 2^k, + * where k is maximal such that 2^k divides n. + * + * Precondition: number is non-zero. + */ +rational rational::pseudo_inverse(unsigned num_bits) const { + rational result; + rational const& n = *this; + SASSERT(!n.is_zero()); // TODO: or we define pseudo-inverse of 0 as 0. + unsigned const k = n.trailing_zeros(); + rational const odd = machine_div2k(n, k); + VERIFY(odd.mult_inverse(num_bits, result)); + SASSERT_EQ(mod(n * result, rational::power_of_two(num_bits)), rational::power_of_two(k)); + return result; +} diff --git a/src/util/rational.h b/src/util/rational.h index 2ec822bb4..e45fbd4e3 100644 --- a/src/util/rational.h +++ b/src/util/rational.h @@ -229,6 +229,12 @@ public: return r; } + friend inline rational machine_div2k(rational const & r1, unsigned k) { + rational r; + rational::m().machine_idiv2k(r1.m_val, k, r.m_val); + return r; + } + friend inline rational mod(rational const & r1, rational const & r2) { rational r; rational::m().mod(r1.m_val, r2.m_val, r.m_val); @@ -355,6 +361,7 @@ public: } bool mult_inverse(unsigned num_bits, rational & result) const; + rational pseudo_inverse(unsigned num_bits) const; static rational const & zero() { return m_zero; diff --git a/src/util/var_queue.h b/src/util/var_queue.h index 62df77784..7245153ca 100644 --- a/src/util/var_queue.h +++ b/src/util/var_queue.h @@ -89,6 +89,10 @@ public: } return out; } + + using const_iterator = decltype(m_queue)::const_iterator; + const_iterator begin() const { return m_queue.begin(); } + const_iterator end() const { return m_queue.end(); } }; inline std::ostream& operator<<(std::ostream& out, var_queue const& queue) {