diff --git a/src/math/lp/nla_params.pyg b/src/math/lp/nla_params.pyg index b141c8cb1..fd7f652a4 100644 --- a/src/math/lp/nla_params.pyg +++ b/src/math/lp/nla_params.pyg @@ -13,7 +13,7 @@ def_module_params('nla', ('grobner_expr_degree_growth', UINT, 2, 'grobner\'s maximum expr degree growth'), ('grobner_max_simplified', UINT, 10000, 'grobner\'s maximum number of simplifications'), ('grobner_cnfl_to_report', UINT, 1, 'grobner\'s maximum number of conflicts to report'), - ('gr_q', UINT, 8, 'grobner\'s quota'), + ('gr_q', UINT, 10, 'grobner\'s quota'), ('grobner_subs_fixed', UINT, 2, '0 - no subs, 1 - substitute, 2 - substitute fixed zeros only') )) diff --git a/src/math/lp/nla_tangent_lemmas.cpp b/src/math/lp/nla_tangent_lemmas.cpp index 9606b609f..af39205fa 100644 --- a/src/math/lp/nla_tangent_lemmas.cpp +++ b/src/math/lp/nla_tangent_lemmas.cpp @@ -21,14 +21,200 @@ #include "math/lp/nla_core.h" namespace nla { -template rational tangents::val(T const& t) const { return m_core->val(t); } -tangents::tangents(core * c) : common(c, nullptr) {} +struct imp { + point m_a; + point m_b; + point m_xy; + bool m_a_is_ok; + bool m_b_is_ok; + 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; + const factor& m_x; + const factor& m_y; + lpvar m_jx; + lpvar m_jy; + tangents& m_tang; + imp(point xy, + const rational& v, + lpvar j, // the monic variable + const monic& m, + const factor& x, + const factor& y, + 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_m(m), + m_x(x), + m_y(y), + m_jx(tang.var(x)), + m_jy(tang.var(y)), + m_tang(tang) {} -std::ostream& tangents::print_tangent_domain(const point &a, const point &b, std::ostream& out) const { - return out << "(" << a << ", " << b << ")"; -} + core & c() { return m_tang.c(); } + + void generate_explanations_of_tang_lemma(lp::explanation& exp) { + // here we repeat the same explanation for each lemma + c().explain(m_m, exp); + c().explain(m_x, exp); + c().explain(m_y, exp); + } + void generate_simple_tangent_lemma(const monic& m, const factorization&); + void tangent_lemma_on_bf() { + get_tang_points(); + TRACE("nla_solver", tout << "tang domain = "; print_tangent_domain(tout) << std::endl;); + generate_two_tang_lines(); + if (m_a_is_ok) + generate_tang_plane(m_a); + if (m_b_is_ok) + generate_tang_plane(m_b); + } + + + void generate_tang_plane(const point & pl) { + c().add_empty_lemma(); + c().negate_relation(m_jx, pl.x); + c().negate_relation(m_jy, pl.y); +#if Z3DEBUG + int mult_sign = nla::rat_sign(pl.x - c().val(m_jx))*nla::rat_sign(pl.y - c().val(m_jy)); + SASSERT((mult_sign == 1) == m_below); + // If "mult_sign is 1" then (a - x)(b-y) > 0 and ab - bx - ay + xy > 0 + // or -ab + bx + ay < xy or -ay - bx + xy > -ab + // val(j) stands for xy. So, finally we have -ay - bx + j > - ab +#endif + + lp::lar_term t; + t.add_monomial(- pl.x, m_jy); + t.add_monomial(- pl.y, m_jx); + t.add_var(m_j); + c().mk_ineq(t, m_below? llc::GT : llc::LT, - pl.x*pl.y); + } + + void generate_two_tang_lines() { + m_tang.add_empty_lemma(); + c().mk_ineq(m_jx, llc::NE, m_xy.x); + c().mk_ineq(m_j, - m_xy.x, m_jy, llc::EQ); + + m_tang.add_empty_lemma(); + c().mk_ineq(m_jy, llc::NE, m_xy.y); + c().mk_ineq(m_j, - m_xy.y, m_jx, llc::EQ); + } + // 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() { + const rational& x = m_xy.x; + const rational& y = m_xy.y; + if (!m_below){ + m_a = point(x - rational(1), y + rational(1)); + m_b = point(x + rational(1), y - rational(1)); + } + else { + // denote x = xy.x and y = xy.y, and vx, vy - the value of x and y. + // we have val(xy) < vx*y + vy*x - vx*vy = pl(x, y); + // The plane with delta (1, 1) is (vx + 1)y + (vy + 1)x - (vx + 1)(vy + 1) = + // vx*y + vy*x - vx*vy + y + x - xv*vy - vx - vy - 1 = pl(x, y) - 1 + // For integers the last expression is greater than or equal to val(xy) when x = vx and y = vy. + // If x <= vx+1 and y <= vy+1 then (vx+1-x)*(vy+1-y) > 0, that creates a cut + // - (vx + 1)y - (vy + 1)x + xy > - (vx+1)*(vx+1) + m_a = point(x - rational(1), y - rational(1)); + m_b = point(x + rational(1), y + rational(1)); + } + } + + void push_tang_point(point & a) { + int steps = 10; + point del = a - m_xy; + while (steps--) { + del *= rational(2); + 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;); + return; + } + a = na; + } + } + + bool pull_tang_point(point & a ) { + if (plane_is_correct_cut(a)) + return true; + point del = a - m_xy; + unsigned steps = 10; + while (steps--) { + del /= rational(2); + point na = m_xy + del; + TRACE("nla_solver_tp", tout << "del = " << del << std::endl;); + if (plane_is_correct_cut(na)) { + a = na; + TRACE("nla_solver_tp", tout << "exit";tout << std::endl;); + return true; + } + } + return false; + } + + rational tang_plane(const point& a) const { + return a.x * m_xy.y + a.y * m_xy.x - a.x * a.y; + } + + void get_tang_points() { + get_initial_tang_points(); + TRACE("nla_solver", tout << "xy = " << m_xy << ", correct val = " << m_correct_v; + tout << "\ntang points:"; print_tangent_domain(tout);tout << std::endl;); + bool all_ints = m_v.is_int() && m_xy.x.is_int() && m_xy.y.is_int(); + if (!all_ints) { + m_a_is_ok = pull_tang_point(m_a); + m_b_is_ok = pull_tang_point(m_b); + } else { + m_a_is_ok = m_b_is_ok = true; + } + if (m_a_is_ok) { + push_tang_point(m_a); + TRACE("nla_solver", tout << "pushed a = " << m_a << std::endl;); + } + if (m_b_is_ok) { + push_tang_point(m_b); + TRACE("nla_solver", tout << "pushed b = " << m_b << std::endl;); + } + TRACE("nla_solver", + if (m_a_is_ok) { tout << "tang_plane(a) = " << tang_plane(m_a) << " , val = " << m_v; } + if (m_b_is_ok) { tout << "\ntang_plane(b) = " << tang_plane(m_b) << " , val = " << m_v << std::endl;}); + } + + std::ostream& print_tangent_domain(std::ostream& out) { + if (m_a_is_ok && m_b_is_ok) { + out << "(" << m_a << ", " << m_b << ")"; + } else if (m_a_is_ok) { + out << m_a; + } + else if (m_b_is_ok) { + out << m_b; + } else { + out << "no a, no b"; + } + return out; + } + + bool plane_is_correct_cut(const point& plane) const { + TRACE("nla_solver", tout << "plane = " << plane << "\n"; + 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 px = tang_plane(plane); + return ((m_correct_v - px)*sign).is_pos() && !((px - m_v)*sign).is_neg(); + } +}; + +tangents::tangents(core * c) : common(c, nullptr) {} + void tangents::tangent_lemma() { if (!c().m_nla_settings.run_tangents()) { TRACE("nla_solver", tout << "not generating tangent lemmas\n";); @@ -37,7 +223,28 @@ void tangents::tangent_lemma() { factorization bf(nullptr); const monic* m; if (c().find_bfc_to_refine(m, bf)) { - tangent_lemma_bf(*m, bf); + unsigned lemmas_size_was = c().m_lemma_vec->size(); + unsigned j = m->var(); + imp i(point(val(bf[0]), val(bf[1])), + c().val(j), + j, + *m, + bf[0], + bf[1], + *this); + i.tangent_lemma_on_bf(); + if (!bf.is_mon()) { + lp::explanation expl; + generate_explanations_of_tang_lemma(*m, bf, expl); + for (unsigned i = lemmas_size_was; i < c().m_lemma_vec->size(); i++) { + auto &l = ((*c().m_lemma_vec)[i]); + l.expl().add(expl); + } + } + TRACE("nla_solver", + for (unsigned i = lemmas_size_was; i < c().m_lemma_vec->size(); i++) + c().print_specific_lemma((*c().m_lemma_vec)[i], tout); ); + } } @@ -48,134 +255,4 @@ void tangents::generate_explanations_of_tang_lemma(const monic& rm, const factor c().explain(bf[1], exp); } -void tangents::generate_tang_plane(const rational & a, const rational& b, const factor& x, const factor& y, bool below, lpvar j) { - lpvar jx = var(x); - lpvar jy = var(y); - add_empty_lemma(); - c().negate_relation(jx, a); - c().negate_relation(jy, b); -#if Z3DEBUG - int mult_sign = nla::rat_sign(a - val(jx))*nla::rat_sign(b - val(jy)); - SASSERT((mult_sign == 1) == below); - // If "mult_sign is 1" then (a - x)(b-y) > 0 and ab - bx - ay + xy > 0 - // or -ab + bx + ay < xy or -ay - bx + xy > -ab - // val(j) stands for xy. So, finally we have -ay - bx + j > - ab -#endif - - lp::lar_term t; - t.add_monomial(-a, jy); - t.add_monomial(-b, jx); - t.add_var(j); - c().mk_ineq(t, below? llc::GT : llc::LT, - a*b); -} - -void tangents::tangent_lemma_bf(const monic& m, const factorization& bf){ - point a, b; - point xy (val(bf[0]), val(bf[1])); - rational correct_mult_val = xy.x * xy.y; - lpvar j =m.var(); - SASSERT(canonize_sign(bf) == canonize_sign(m)); - rational v = val(j); - bool below = v < correct_mult_val; - TRACE("nla_solver", tout << "below = " << below << std::endl; ); - get_tang_points(a, b, below, v, xy); - TRACE("nla_solver", tout << "tang domain = "; print_tangent_domain(a, b, tout); tout << std::endl;); - unsigned lemmas_size_was = c().m_lemma_vec->size(); - rational sign(1); - generate_two_tang_lines(bf, xy, j); - generate_tang_plane(a.x, a.y, bf[0], bf[1], below, j); - generate_tang_plane(b.x, b.y, bf[0], bf[1], below, j); - - if (!bf.is_mon()) { - lp::explanation expl; - generate_explanations_of_tang_lemma(m, bf, expl); - for (unsigned i = lemmas_size_was; i < c().m_lemma_vec->size(); i++) { - auto &l = ((*c().m_lemma_vec)[i]); - l.expl().add(expl); - } - } - TRACE("nla_solver", - for (unsigned i = lemmas_size_was; i < c().m_lemma_vec->size(); i++) - c().print_specific_lemma((*c().m_lemma_vec)[i], tout); ); -} - - -void tangents::generate_two_tang_lines(const factorization & bf, const point& xy, lpvar j) { - add_empty_lemma(); - c().mk_ineq(var(bf[0]), llc::NE, xy.x); - c().mk_ineq(j, - xy.x, var(bf[1]), llc::EQ); - - add_empty_lemma(); - c().mk_ineq(var(bf[1]), llc::NE, xy.y); - c().mk_ineq(j, - xy.y, var(bf[0]), llc::EQ); -} - -// Get two planes tangent to surface z = xy, one at point a, and another at point b. -// One can show that these planes still create a cut. -void tangents::get_initial_tang_points(point &a, point &b, const point& xy, - bool below) const { - const rational& x = xy.x; - const rational& y = xy.y; - if (!below){ - a = point(x - rational(1), y + rational(1)); - b = point(x + rational(1), y - rational(1)); - } - else { - a = point(x - rational(1), y - rational(1)); - b = point(x + rational(1), y + rational(1)); - } -} - -void tangents::push_tang_point(point &a, const point& xy, bool below, const rational& correct_val, const rational& val) const { - SASSERT(correct_val == xy.x * xy.y); - int steps = 10; - point del = a - xy; - while (steps--) { - del *= rational(2); - point na = xy + del; - TRACE("nla_solver_tp", tout << "del = " << del << std::endl;); - if (!plane_is_correct_cut(na, xy, correct_val, val, below)) { - TRACE("nla_solver_tp", tout << "exit";tout << std::endl;); - return; - } - a = na; - } -} - -void tangents::push_tang_points(point &a, point &b, const point& xy, bool below, const rational& correct_val, const rational& val) const { - push_tang_point(a, xy, below, correct_val, val); - push_tang_point(b, xy, below, correct_val, val); -} - -rational tangents::tang_plane(const point& a, const point& x) const { - return a.x * x.y + a.y * x.x - a.x * a.y; -} - -bool tangents:: plane_is_correct_cut(const point& plane, - const point& xy, - const rational & correct_val, - const rational & val, - bool below) const { - SASSERT(correct_val == xy.x * xy.y); - if (below && val > correct_val) return false; - rational sign = below? rational(1) : rational(-1); - rational px = tang_plane(plane, xy); - return ((correct_val - px)*sign).is_pos() && !((px - val)*sign).is_neg(); - -} - -// "below" means that the val is below the surface xy -void tangents::get_tang_points(point &a, point &b, bool below, const rational& val, - const point& xy) const { - get_initial_tang_points(a, b, xy, below); - auto correct_val = xy.x * xy.y; - TRACE("nla_solver", tout << "xy = " << xy << ", correct val = " << xy.x * xy.y; - tout << "\ntang points:"; print_tangent_domain(a, b, tout);tout << std::endl;); - TRACE("nla_solver", tout << "tang_plane(a, xy) = " << tang_plane(a, xy) << " , val = " << val; - tout << "\ntang_plane(b, xy) = " << tang_plane(b, xy); tout << std::endl;); - SASSERT(plane_is_correct_cut(a, xy, correct_val, val, below)); - SASSERT(plane_is_correct_cut(b, xy, correct_val, val, below)); - push_tang_points(a, b, xy, below, correct_val, val); - TRACE("nla_solver", tout << "pushed a = " << a << "\npushed b = " << b << std::endl;); -} } diff --git a/src/math/lp/nla_tangent_lemmas.h b/src/math/lp/nla_tangent_lemmas.h index 24dda0e13..95c771150 100644 --- a/src/math/lp/nla_tangent_lemmas.h +++ b/src/math/lp/nla_tangent_lemmas.h @@ -35,6 +35,11 @@ struct point { y *= a; return *this; } + inline point& operator/=(rational a) { + x /= a; + y /= a; + return *this; + } inline point operator+(const point& b) const { return point(x + b.x, y + b.y); } @@ -47,37 +52,9 @@ struct point { inline std::ostream& operator<<(std::ostream& out, point const& a) { return out << "(" << a.x << ", " << a.y << ")"; } -class tangents : common { -public: +struct tangents : common { tangents(core *core); - void tangent_lemma(); -private: - lpvar find_binomial_to_refine(); - void generate_explanations_of_tang_lemma(const monic& m, const factorization& bf, lp::explanation& exp); - void generate_simple_tangent_lemma(const monic& m, const factorization&); - void tangent_lemma_bf(const monic& m,const factorization& bf); - void generate_tang_plane(const rational & a, const rational& b, const factor& x, const factor& y, bool below, lpvar j); - - void generate_two_tang_lines(const factorization & bf, const point& xy, lpvar j); - // Get two planes tangent to surface z = xy, one at point a, and another at point b. - // One can show that these planes still create a cut. - void get_initial_tang_points(point &a, point &b, const point& xy, bool below) const; - - void push_tang_point(point &a, const point& xy, bool below, const rational& correct_val, const rational& val) const; - - void push_tang_points(point &a, point &b, const point& xy, bool below, const rational& correct_val, const rational& val) const; - - rational tang_plane(const point& a, const point& x) const; - void get_tang_points(point &a, point &b, bool below, const rational& val, const point& xy) const; - std::ostream& print_point(const point &a, std::ostream& out) const; - std::ostream& print_tangent_domain(const point &a, const point &b, std::ostream& out) const; - // "below" means that the val is below the surface xy - bool plane_is_correct_cut(const point& plane, - const point& xy, - const rational & correct_val, - const rational & val, - bool below) const; - template rational val(T const& t) const; - template lpvar var(T const& t) const { return t.var(); } + void tangent_lemma(); + void generate_explanations_of_tang_lemma(const monic& rm, const factorization& bf, lp::explanation& exp); }; // end of tangents -} +} // end of namespace diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 7a0914bec..7373d6a03 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -1735,19 +1735,36 @@ public: bool is_int = offset.is_int(); u_map coeffs; term2coeffs(term, coeffs); - TRACE("arith", - lp().print_term(term, tout << "term: ") << "\n"; - for (auto const& kv : coeffs) { - tout << "v" << kv.m_key << " * " << kv.m_value << "\n"; + TRACE("arith", + { + bool all_ints = true; + lp().print_term(term, tout << "term: ") << "\n"; + for (auto const& kv : coeffs) { + if (kv.m_value.is_int() == false) + all_ints = false; + tout << "v" << kv.m_key << " * " << kv.m_value << "\n"; + } + tout << offset << "\n"; + if (all_ints) { + rational g(0); + for (auto const& kv : coeffs) { + g = gcd(g, kv.m_value); + } + tout << "gcd: " << g << "\n"; + } } - tout << offset << "\n"; - rational g(0); - for (auto const& kv : coeffs) { - g = gcd(g, kv.m_value); - } - tout << "gcd: " << g << "\n"; ); + bool all_ints = true; if (is_int) { + for (auto const& kv : coeffs) { + if (kv.m_value.is_int() == false) { + all_ints = false; + break; + } + } + } + + if (is_int && all_ints) { // 3x + 6y >= 5 -> x + 3y >= 5/3, then x + 3y >= 2 // 3x + 6y <= 5 -> x + 3y <= 1 diff --git a/src/test/lp/nla_solver_test.cpp b/src/test/lp/nla_solver_test.cpp index e7d36bbd3..f31793d21 100644 --- a/src/test/lp/nla_solver_test.cpp +++ b/src/test/lp/nla_solver_test.cpp @@ -717,6 +717,34 @@ void test_monotone_lemma() { */ } +void test_tangent_lemma_rat() { + enable_trace("nla_solver"); + lp::lar_solver s; + unsigned a = s.number_of_vars(); + unsigned b = a + 1; + unsigned ab = b + 1; + + lpvar lp_a = s.add_named_var(a, true, "a"); + lpvar lp_b = s.add_named_var(b, false, "b"); + lpvar lp_ab = s.add_named_var(ab, false, "ab"); + s_set_column_value(s, lp_a, rational(3)); + s_set_column_value(s, lp_b, rational(4)); + rational v = rational(12) + rational (1)/rational(7); + s_set_column_value(s, lp_ab, v); + reslimit l; + params_ref p; + solver nla(s); + // create monomial ab + vector vec; + vec.push_back(lp_a); + vec.push_back(lp_b); + nla.add_monic(lp_ab, vec.size(), vec.begin()); + + vector lemma; + SASSERT(nla.get_core()->test_check(lemma) == l_false); + nla.get_core()->print_lemma(std::cout); +} + void test_tangent_lemma_reg() { enable_trace("nla_solver"); lp::lar_solver s; @@ -793,7 +821,8 @@ void test_tangent_lemma_equiv() { void test_tangent_lemma() { - test_tangent_lemma_reg(); + test_tangent_lemma_rat(); + test_tangent_lemma_reg(); test_tangent_lemma_equiv(); }