diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp index 69cdb3364..2e219f5bd 100644 --- a/src/math/lp/monomial_bounds.cpp +++ b/src/math/lp/monomial_bounds.cpp @@ -21,8 +21,7 @@ namespace nla { bool propagated = false; for (lpvar v : c().m_to_refine) { monic const& m = c().emons()[v]; - if (propagate(m)) - propagated = true; + propagated |= propagate(m); } return propagated; } @@ -32,14 +31,12 @@ namespace nla { */ void monomial_bounds::compute_product(unsigned start, monic const& m, scoped_dep_interval& product) { scoped_dep_interval vi(dep); + unsigned power = 1; for (unsigned i = start; i < m.size(); ) { lpvar v = m.vars()[i]; - unsigned power = 1; var2interval(v, vi); ++i; - for (; i < m.size() && m.vars()[i] == v; ++i) { - ++power; - } + for (power = 1; i < m.size() && m.vars()[i] == v; ++i, ++power); dep.power(vi, power, vi); dep.mul(product, vi, product); } @@ -79,6 +76,68 @@ namespace nla { } } + /** + * val(v)^p should be in range. + * if val(v)^p > upper(range) add + * v <= root(p, upper(range)) and v >= -root(p, upper(range)) if p is even + * v <= root(p, upper(range)) if p is odd + * if val(v)^p < lower(range) add + * v >= root(p, lower(range)) or v <= -root(p, lower(range)) if p is even + * v >= root(p, lower(range)) if p is odd + */ + + bool monomial_bounds::propagate_value(dep_interval& range, lpvar v, unsigned p) { + SASSERT(p > 0); + if (p == 1) + return propagate_value(range, v); + auto val = c().val(v); + val = power(val, p); + rational r; + if (dep.is_below(range, val)) { + lp::explanation ex; + dep.get_upper_dep(range, ex); + if (p % 2 == 0 && rational(dep.upper(range)).is_neg()) { + new_lemma lemma(c(), "range requires a non-negative upper bound"); + lemma &= ex; + return true; + } + if (rational(dep.upper(range)).root(p, r)) { + { + auto le = dep.upper_is_open(range) ? llc::LT : llc::LE; + new_lemma lemma(c(), "propagate value - root case - lower bound of range is below value"); + lemma &= ex; + lemma |= ineq(v, le, r); + } + if (p % 2 == 0) { + SASSERT(!r.is_neg()); + auto ge = dep.upper_is_open(range) ? llc::GT : llc::GE; + new_lemma lemma(c(), "propagate value - root case - lower bound of range is below value"); + lemma &= ex; + lemma |= ineq(v, ge, -r); + } + return true; + } + // TBD: add bounds as long as difference to val is above some epsilon. + } + else if (dep.is_above(range, val)) { + if (rational(dep.lower(range)).root(p, r)) { + lp::explanation ex; + dep.get_lower_dep(range, ex); + auto ge = dep.lower_is_open(range) ? llc::GT : llc::GE; + auto le = dep.lower_is_open(range) ? llc::LT : llc::LE; + new_lemma lemma(c(), "propagate value - root case - lower bound of range is above value"); + lemma &= ex; + lemma |= ineq(v, ge, r); + if (p % 2 == 0) { + lemma |= ineq(v, le, -r); + } + return true; + } + // TBD: add bounds as long as difference to val is above some epsilon. + } + return false; + } + void monomial_bounds::var2interval(lpvar v, scoped_dep_interval& i) { lp::constraint_index ci; rational bound; @@ -114,14 +173,10 @@ namespace nla { unsigned num_free, power; lpvar free_var; analyze_monomial(m, num_free, free_var, power); - bool m_is_free = is_free(m.var()); - if (num_free >= 2) - return false; - if (num_free >= 1 && m_is_free) - return false; - SASSERT(num_free == 0 || !m_is_free); bool do_propagate_up = num_free == 0; - bool do_propagate_down = !m_is_free; + bool do_propagate_down = !is_free(m.var()) && num_free <= 1; + if (!do_propagate_up && !do_propagate_down) + return false; scoped_dep_interval product(dep); scoped_dep_interval vi(dep), mi(dep); scoped_dep_interval other_product(dep); @@ -130,16 +185,14 @@ namespace nla { for (unsigned i = 0; i < m.size(); ) { lpvar v = m.vars()[i]; ++i; - unsigned power = 1; - for (; i < m.size() && v == m.vars()[i]; ++i) - ++power; + for (power = 1; i < m.size() && v == m.vars()[i]; ++i, ++power); var2interval(v, vi); dep.power(vi, power, vi); - if (power == 1 && do_propagate_down && (num_free == 0 || free_var == v)) { + if (do_propagate_down && (num_free == 0 || free_var == v)) { dep.set(other_product, product); compute_product(i, m, other_product); - if (propagate_down(m, mi, v, other_product)) + if (propagate_down(m, mi, v, power, other_product)) return true; } dep.mul(product, vi, product); @@ -147,12 +200,12 @@ namespace nla { return do_propagate_up && propagate_value(product, m.var()); } - bool monomial_bounds::propagate_down(monic const& m, dep_interval& mi, lpvar v, dep_interval& product) { + bool monomial_bounds::propagate_down(monic const& m, dep_interval& mi, lpvar v, unsigned power, dep_interval& product) { if (!dep.separated_from_zero(product)) return false; scoped_dep_interval range(dep); dep.div(mi, product, range); - return propagate_value(range, v); + return propagate_value(range, v, power); } bool monomial_bounds::is_free(lpvar v) const { @@ -164,19 +217,24 @@ namespace nla { c().has_lower_bound(v) && c().has_upper_bound(v) && c().get_lower_bound(v).is_zero() && - c().get_lower_bound(v) == c().get_upper_bound(v); + c().get_upper_bound(v).is_zero(); } + /** + * Count the number of unbound (free) variables. + * Variables with no lower and no upper bound multiplied + * to an odd degree have unbound ranges when it comes to + * bounds propagation. + */ void monomial_bounds::analyze_monomial(monic const& m, unsigned& num_free, lpvar& fv, unsigned& fv_power) const { - unsigned power = 0; + unsigned power = 1; num_free = 0; fv = null_lpvar; fv_power = 0; for (unsigned i = 0; i < m.vars().size(); ) { lpvar v = m.vars()[i]; - unsigned power = 1; ++i; - for (; i < m.vars().size() && m.vars()[i] == v; ++i, ++power); + for (power = 1; i < m.vars().size() && m.vars()[i] == v; ++i, ++power); if (is_zero(v)) { num_free = 0; return; diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h index aaf39d175..e231e282c 100644 --- a/src/math/lp/monomial_bounds.h +++ b/src/math/lp/monomial_bounds.h @@ -19,9 +19,10 @@ namespace nla { void var2interval(lpvar v, scoped_dep_interval& i); bool propagate_down(monic const& m, lpvar u); bool propagate_value(dep_interval& range, lpvar v); + bool propagate_value(dep_interval& range, lpvar v, unsigned power); void compute_product(unsigned start, monic const& m, scoped_dep_interval& i); bool propagate(monic const& m); - bool propagate_down(monic const& m, dep_interval& mi, lpvar v, dep_interval& product); + bool propagate_down(monic const& m, dep_interval& mi, lpvar v, unsigned power, dep_interval& product); void analyze_monomial(monic const& m, unsigned& num_free, lpvar& free_v, unsigned& power) const; bool is_free(lpvar v) const; bool is_zero(lpvar v) const; diff --git a/src/math/lp/nla_tangent_lemmas.cpp b/src/math/lp/nla_tangent_lemmas.cpp index 42741a07e..d01dc51ca 100644 --- a/src/math/lp/nla_tangent_lemmas.cpp +++ b/src/math/lp/nla_tangent_lemmas.cpp @@ -11,53 +11,56 @@ namespace nla { -struct tangent_imp { - point m_a; - point m_b; - point m_xy; - rational m_correct_v; +class tangent_imp { + point m_a; + point m_b; + point m_xy; + rational m_correct_v; // "below" means that the incorrect value is less than the correct one, that is m_v < m_correct_v - bool m_below; - rational m_v; // the monomial value - lpvar m_j; // the monic variable - const monic& m_m; + bool m_below; + rational m_v; // the monomial value + lpvar m_j; // the monic variable + const monic& m_m; const factor& m_x; const factor& m_y; - lpvar m_jx; - lpvar m_jy; - tangents& m_tang; - bool m_is_mon; + lpvar m_jx; + lpvar m_jy; + tangents& m_tang; + bool m_is_mon; + +public: tangent_imp(point xy, const rational& v, - lpvar j, // the monic variable const monic& m, const factorization& f, tangents& tang) : m_xy(xy), m_correct_v(xy.x * xy.y), m_below(v < m_correct_v), m_v(v), - m_j(tang.var(m)), + m_j(m.var()), m_m(m), m_x(f[0]), m_y(f[1]), - m_jx(tang.var(m_x)), - m_jy(tang.var(m_y)), + m_jx(m_x.var()), + m_jy(m_y.var()), m_tang(tang), m_is_mon(f.is_mon()) { SASSERT(f.size() == 2); } - - - core & c() { return m_tang.c(); } - - void tangent_lemma_on_bf() { - get_tang_points(); - TRACE("nla_solver", tout << "tang domain = "; print_tangent_domain(tout) << std::endl;); - generate_two_tang_lines(); - generate_tang_plane(m_a); - generate_tang_plane(m_b); + + void operator()() { + get_points(); + TRACE("nla_solver", print_tangent_domain(tout << "tang domain = ") << std::endl;); + generate_line1(); + generate_line2(); + generate_plane(m_a); + generate_plane(m_b); } +private: + + core & c() { return m_tang.c(); } + void explain(new_lemma& lemma) { if (!m_is_mon) { lemma &= m_m; @@ -66,7 +69,7 @@ struct tangent_imp { } } - void generate_tang_plane(const point & pl) { + void generate_plane(const point & pl) { new_lemma lemma(c(), "generate tangent plane"); c().negate_relation(lemma, m_jx, m_x.rat_sign()*pl.x); c().negate_relation(lemma, m_jy, m_y.rat_sign()*pl.y); @@ -86,24 +89,24 @@ struct tangent_imp { lemma |= ineq(t, m_below? llc::GT : llc::LT, - pl.x*pl.y); explain(lemma); } - - void generate_two_tang_lines() { - { - new_lemma lemma(c(), "two tangent planes 1"); - // Should be v = val(m_x)*val(m_y), and val(factor) = factor.rat_sign()*var(factor.var()) - lemma |= ineq(m_jx, llc::NE, c().val(m_jx)); - lemma |= ineq(lp::lar_term(m_j, - m_y.rat_sign() * m_xy.x, m_jy), llc::EQ, 0); - explain(lemma); - } - { - new_lemma lemma(c(), "two tangent planes 2"); - lemma |= ineq(m_jy, llc::NE, c().val(m_jy)); - lemma |= ineq(lp::lar_term(m_j, - m_x.rat_sign() * m_xy.y, m_jx), llc::EQ, 0); - explain(lemma); - } + + void generate_line1() { + new_lemma lemma(c(), "tangent line 1"); + // Should be v = val(m_x)*val(m_y), and val(factor) = factor.rat_sign()*var(factor.var()) + lemma |= ineq(m_jx, llc::NE, c().val(m_jx)); + lemma |= ineq(lp::lar_term(m_j, - m_y.rat_sign() * m_xy.x, m_jy), llc::EQ, 0); + explain(lemma); } + + void generate_line2() { + new_lemma lemma(c(), "tangent line 2"); + lemma |= ineq(m_jy, llc::NE, c().val(m_jy)); + lemma |= ineq(lp::lar_term(m_j, - m_x.rat_sign() * m_xy.y, m_jx), llc::EQ, 0); + explain(lemma); + } + // Get two planes tangent to surface z = xy, one at point a, and another at point b, creating a cut - void get_initial_tang_points() { + void get_initial_points() { const rational& x = m_xy.x; const rational& y = m_xy.y; bool all_ints = m_v.is_int() && x.is_int() && y.is_int(); @@ -130,7 +133,7 @@ struct tangent_imp { } } - void push_tang_point(point & a) { + void push_point(point & a) { SASSERT(plane_is_correct_cut(a)); int steps = 10; point del = a - m_xy; @@ -139,7 +142,7 @@ struct tangent_imp { point na = m_xy + del; TRACE("nla_solver_tp", tout << "del = " << del << std::endl;); if (!plane_is_correct_cut(na)) { - TRACE("nla_solver_tp", tout << "exit";tout << std::endl;); + TRACE("nla_solver_tp", tout << "exit\n";); return; } a = na; @@ -147,25 +150,24 @@ struct tangent_imp { } rational tang_plane(const point& a) const { - return a.x * m_xy.y + a.y * m_xy.x - a.x * a.y; + return a.x * m_xy.y + a.y * m_xy.x - a.x * a.y; } - void get_tang_points() { - get_initial_tang_points(); + void get_points() { + get_initial_points(); TRACE("nla_solver", tout << "xy = " << m_xy << ", correct val = " << m_correct_v; - tout << "\ntang points:"; print_tangent_domain(tout);tout << std::endl;); - push_tang_point(m_a); - TRACE("nla_solver", tout << "pushed a = " << m_a << std::endl;); - - push_tang_point(m_b); - TRACE("nla_solver", tout << "pushed b = " << m_b << std::endl;); + print_tangent_domain(tout << "\ntang points:") << std::endl;); + push_point(m_a); + push_point(m_b); TRACE("nla_solver", - tout << "tang_plane(a) = " << tang_plane(m_a) << " , val = " << m_v << ", tang_plane(b) = " << tang_plane(m_b) << " , val = " << std::endl;); + tout << "pushed a = " << m_a << std::endl + << "pushed b = " << m_b << std::endl + << "tang_plane(a) = " << tang_plane(m_a) << " , val = " << m_a << ", " + << "tang_plane(b) = " << tang_plane(m_b) << " , val = " << m_b << std::endl;); } std::ostream& print_tangent_domain(std::ostream& out) { - out << "(" << m_a << ", " << m_b << ")"; - return out; + return out << "(" << m_a << ", " << m_b << ")"; } bool plane_is_correct_cut(const point& plane) const { @@ -173,7 +175,7 @@ struct tangent_imp { tout << "tang_plane() = " << tang_plane(plane) << ", v = " << m_v << ", correct_v = " << m_correct_v << "\n";); SASSERT((m_below && m_v < m_correct_v) || ((!m_below) && m_v > m_correct_v)); - rational sign = m_below? rational(1) : rational(-1); + rational sign = rational(m_below ? 1 : -1); rational px = tang_plane(plane); return ((m_correct_v - px)*sign).is_pos() && !((px - m_v)*sign).is_neg(); } @@ -182,21 +184,12 @@ struct tangent_imp { tangents::tangents(core * c) : common(c) {} void tangents::tangent_lemma() { - if (!c().m_nla_settings.run_tangents()) { - TRACE("nla_solver", tout << "not generating tangent lemmas\n";); - return; - } factorization bf(nullptr); - const monic* m; - if (c().find_bfc_to_refine(m, bf)) { - unsigned j = m->var(); - tangent_imp i(point(val(bf[0]), val(bf[1])), - c().val(j), - j, - *m, - bf, - *this); - i.tangent_lemma_on_bf(); + const monic* m = nullptr; + if (c().m_nla_settings.run_tangents() && c().find_bfc_to_refine(m, bf)) { + lpvar j = m->var(); + tangent_imp tangent(point(val(bf[0]), val(bf[1])), c().val(j), *m, bf, *this); + tangent(); } } diff --git a/src/smt/seq_skolem.h b/src/smt/seq_skolem.h index 4800047ad..21711ed7b 100644 --- a/src/smt/seq_skolem.h +++ b/src/smt/seq_skolem.h @@ -35,7 +35,7 @@ namespace smt { symbol m_seq_first, m_seq_last; symbol m_indexof_left, m_indexof_right; // inverse of indexof: (indexof_left s t) + s + (indexof_right s t) = t, for s in t. symbol m_aut_step; // regex unfolding state - symbol m_accept, m_reject; // regex + symbol m_accept; // regex symbol m_pre, m_post; // inverse of at: (pre s i) + (at s i) + (post s i) = s if 0 <= i < (len s) symbol m_eq; // equality atom symbol m_seq_align;