diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt index 6ec8ba12d..477577345 100644 --- a/src/math/lp/CMakeLists.txt +++ b/src/math/lp/CMakeLists.txt @@ -38,6 +38,7 @@ z3_add_component(lp nla_intervals.cpp nla_monotone_lemmas.cpp nla_order_lemmas.cpp + nla_powers.cpp nla_solver.cpp nla_tangent_lemmas.cpp nra_solver.cpp diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index efd510824..b072112f0 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -29,6 +29,7 @@ core::core(lp::lar_solver& s, reslimit & lim) : m_basics(this), m_order(this), m_monotone(this), + m_powers(*this), m_intervals(this, lim), m_monomial_bounds(this), m_horner(this), @@ -120,9 +121,8 @@ bool core::canonize_sign(const monic& m) const { bool core::canonize_sign(const factorization& f) const { bool r = false; - for (const factor & a : f) { + for (const factor & a : f) r ^= canonize_sign(a); - } return r; } @@ -1477,6 +1477,10 @@ void core::check_weighted(unsigned sz, std::pair& l_vec) { + m_lemma_vec = &l_vec; + return m_powers.check(r, x, y, l_vec); +} lbool core::check(vector& l_vec) { lp_settings().stats().m_nla_calls++; diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h index 0e759cf59..1e66e4cd7 100644 --- a/src/math/lp/nla_core.h +++ b/src/math/lp/nla_core.h @@ -19,6 +19,7 @@ #include "math/lp/nla_order_lemmas.h" #include "math/lp/nla_monotone_lemmas.h" #include "math/lp/nla_grobner.h" +#include "math/lp/nla_powers.h" #include "math/lp/emonics.h" #include "math/lp/nla_settings.h" #include "math/lp/nex.h" @@ -42,104 +43,6 @@ bool try_insert(const A& elem, B& collection) { return true; } -typedef lp::constraint_index lpci; -typedef lp::lconstraint_kind llc; -typedef lp::constraint_index lpci; -typedef lp::explanation expl_set; -typedef lp::var_index lpvar; -const lpvar null_lpvar = UINT_MAX; - -inline int rat_sign(const rational& r) { return r.is_pos()? 1 : ( r.is_neg()? -1 : 0); } -inline rational rrat_sign(const rational& r) { return rational(rat_sign(r)); } -inline bool is_set(unsigned j) { return j != null_lpvar; } -inline bool is_even(unsigned k) { return (k & 1) == 0; } -class ineq { - lp::lconstraint_kind m_cmp; - lp::lar_term m_term; - rational m_rs; -public: - ineq(lp::lconstraint_kind cmp, const lp::lar_term& term, const rational& rs) : m_cmp(cmp), m_term(term), m_rs(rs) {} - ineq(const lp::lar_term& term, lp::lconstraint_kind cmp, int i) : m_cmp(cmp), m_term(term), m_rs(rational(i)) {} - ineq(const lp::lar_term& term, lp::lconstraint_kind cmp, const rational& rs) : m_cmp(cmp), m_term(term), m_rs(rs) {} - ineq(lpvar v, lp::lconstraint_kind cmp, int i): m_cmp(cmp), m_term(v), m_rs(rational(i)) {} - ineq(lpvar v, lp::lconstraint_kind cmp, rational const& r): m_cmp(cmp), m_term(v), m_rs(r) {} - bool operator==(const ineq& a) const { - return m_cmp == a.m_cmp && m_term == a.m_term && m_rs == a.m_rs; - } - const lp::lar_term& term() const { return m_term; }; - lp::lconstraint_kind cmp() const { return m_cmp; }; - const rational& rs() const { return m_rs; }; -}; - -class lemma { - vector m_ineqs; - lp::explanation m_expl; -public: - void push_back(const ineq& i) { m_ineqs.push_back(i);} - size_t size() const { return m_ineqs.size() + m_expl.size(); } - const vector& ineqs() const { return m_ineqs; } - vector& ineqs() { return m_ineqs; } - lp::explanation& expl() { return m_expl; } - const lp::explanation& expl() const { return m_expl; } - bool is_conflict() const { return m_ineqs.empty() && !m_expl.empty(); } -}; - -class core; -// -// lemmas are created in a scope. -// when the destructor of new_lemma is invoked -// all constraints are assumed added to the lemma -// correctness of the lemma can be checked at this point. -// -class new_lemma { - char const* name; - core& c; - lemma& current() const; - -public: - new_lemma(core& c, char const* name); - ~new_lemma(); - lemma& operator()() { return current(); } - std::ostream& display(std::ostream& out) const; - new_lemma& operator&=(lp::explanation const& e); - new_lemma& operator&=(const monic& m); - new_lemma& operator&=(const factor& f); - new_lemma& operator&=(const factorization& f); - new_lemma& operator&=(lpvar j); - new_lemma& operator|=(ineq const& i); - new_lemma& explain_fixed(lpvar j); - new_lemma& explain_equiv(lpvar u, lpvar v); - new_lemma& explain_var_separated_from_zero(lpvar j); - new_lemma& explain_existing_lower_bound(lpvar j); - new_lemma& explain_existing_upper_bound(lpvar j); - - lp::explanation& expl() { return current().expl(); } - - unsigned num_ineqs() const { return current().ineqs().size(); } -}; - - -inline std::ostream& operator<<(std::ostream& out, new_lemma const& l) { - return l.display(out); -} - -struct pp_fac { - core const& c; - factor const& f; - pp_fac(core const& c, factor const& f): c(c), f(f) {} -}; - -struct pp_var { - core const& c; - lpvar v; - pp_var(core const& c, lpvar v): c(c), v(v) {} -}; - -struct pp_factorization { - core const& c; - factorization const& f; - pp_factorization(core const& c, factorization const& f): c(c), f(f) {} -}; class core { friend struct common; @@ -149,6 +52,7 @@ class core { friend struct basics; friend struct tangents; friend class monotone; + friend class powers; friend struct nla_settings; friend class intervals; friend class horner; @@ -183,6 +87,7 @@ class core { basics m_basics; order m_order; monotone m_monotone; + powers m_powers; intervals m_intervals; monomial_bounds m_monomial_bounds; nla_settings m_nla_settings; @@ -463,9 +368,8 @@ public: bool conflict_found() const; - lbool check(vector& l_vec); - - void set_lemma_vec(vector& l_vec) { m_lemma_vec = &l_vec; } + lbool check(vector& l_vec); + lbool check_power(lpvar r, lpvar x, lpvar y, vector& l_vec); bool no_lemmas_hold() const; diff --git a/src/math/lp/nla_powers.cpp b/src/math/lp/nla_powers.cpp new file mode 100644 index 000000000..ee3c41ee8 --- /dev/null +++ b/src/math/lp/nla_powers.cpp @@ -0,0 +1,181 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + nla_powers.cpp + +Author: + Lev Nachmanson (levnach) + Nikolaj Bjorner (nbjorner) + +Description: + Refines bounds on powers. + + Reference: TOCL-2018, Cimatti et al. + + +Special cases: + +1. Exponentiation. x is fixed numeral. + +TOCL18 axioms: + a^y > 0 (if a > 0) + y = 0 <=> a^y = 1 (if a != 0) + y < 0 <=> a^y < 1 (if a > 1) + y > 0 <=> a^y > 1 (if a > 1) + y != 0 <=> a^y > y + 1 (if a >= 2) + y1 < y2 <=> a^y1 < a^y2 (**) + +Other special case: + + y = 1 <= a^y = a + +TOCL18 approach: Polynomial abstractions + +Taylor: a^y = sum_i ln(a)*y^i/i! + +Truncation: P(n, a) = sum_{i=0}^n ln(a)*y^i/i! + +y = 0: handled by axiom a^y = 1 +y < 0: P(2n-1, y) <= a^y <= P(2n, y), n > 0 because Taylor contribution is negative at odd powers. +y > 0: P(n, y) <= a^y <= P(n, y)*(1 - y^{n+1}/(n+1)!) + + +2. Powers. y is fixed positive integer. + +3. Other + +General case: + + For now the solver integrates just weak monotonicity lemmas: + + - x >= x0 > 0, y >= y0 => x^y >= x0^y0 + - 0 < x <= x0, y <= y0 => x^y <= x0^y0 + + +TODO: + +- Comprehensive integration for truncation polynomial approximation. +- TOCL18 approach includes refinement loop based on precision epsilon. +- accept solvability if r is within a small range of x^y, when x^y is not rational. +- integrate algebraic numbers, or even extension fields (for 'e'). +- integrate monotonicy axioms (**) by tracking exponents across instances. + +anum isn't initialized unless nra_solver is invoked. +there is no proviso for using algebraic numbers outside of the nra solver. +so either we have a rational refinement version _and_ an algebraic numeral refinement +loop or we introduce algebraic numerals outside of the nra_solver + +scoped_anum xval(am()), yval(am()), rval(am()); + +am().set(xval, am_value(x)); +am().set(yval, am_value(y)); +am().set(rval, am_value(r)); + +--*/ +#include "math/lp/nla_core.h" + +namespace nla { + + lbool powers::check(lpvar r, lpvar x, lpvar y, vector& lemmas) { + if (x == null_lpvar || y == null_lpvar || r == null_lpvar) + return l_undef; + + core& c = m_core; + + auto xval = c.val(x); + auto yval = c.val(y); + auto rval = c.val(r); + + lemmas.reset(); + + + if (xval != 0 && yval == 0 && rval != 1) { + new_lemma lemma(c, "x != 0 => x^0 = 1"); + lemma |= ineq(x, llc::EQ, rational::zero()); + lemma |= ineq(y, llc::NE, rational::zero()); + lemma |= ineq(r, llc::EQ, rational::one()); + return l_false; + } + + if (xval == 0 && yval > 0 && rval != 0) { + new_lemma lemma(c, "y != 0 => 0^y = 0"); + lemma |= ineq(x, llc::NE, rational::zero()); + lemma |= ineq(y, llc::EQ, rational::zero()); + lemma |= ineq(r, llc::EQ, rational::zero()); + return l_false; + } + + if (xval > 0 && rval < 0 && rval <= 0) { + new_lemma lemma(c, "x > 0 => x^y > 0"); + lemma |= ineq(x, llc::LE, rational::zero()); + lemma |= ineq(r, llc::GT, rational::zero()); + return l_false; + } + + if (xval > 1 && rval >= 1 && yval < 0) { + new_lemma lemma(c, "x > 1, y < 0 => x^y < 1"); + lemma |= ineq(x, llc::LE, rational::one()); + lemma |= ineq(y, llc::GE, rational::zero()); + lemma |= ineq(r, llc::LT, rational::one()); + return l_false; + } + + if (xval > 1 && rval <= 1 && yval > 0) { + new_lemma lemma(c, "x > 1, y > 0 => x^y > 1"); + lemma |= ineq(x, llc::LE, rational::one()); + lemma |= ineq(y, llc::LE, rational::zero()); + lemma |= ineq(r, llc::GT, rational::one()); + return l_false; + } + + if (xval >= 2 && yval != 0 & rval <= yval + 1) { + new_lemma lemma(c, "x >= 2, y != 0 => x^y > y + 1"); + lemma |= ineq(x, llc::LT, rational(2)); + lemma |= ineq(y, llc::EQ, rational::zero()); + lemma |= ineq(lp::lar_term(r, rational::minus_one(), y), llc::GT, rational::one()); + return l_false; + } + + if (xval > 0 && yval.is_unsigned()) { + auto r2val = power(xval, yval.get_unsigned()); + if (rval == r2val) + return l_true; + if (xval > 0 && r2val < rval) { + SASSERT(yval > 0); + new_lemma lemma(c, "x >= x0 > 0, y >= y0 > 0 => r >= x0^y0"); + lemma |= ineq(x, llc::LT, xval); + lemma |= ineq(y, llc::LT, yval); + lemma |= ineq(r, llc::GE, r2val); + return l_false; + } + if (xval > 0 && r2val < rval) { + new_lemma lemma(c, "x >= x0 > 0, y <= y0 => r <= x0^y0"); + lemma |= ineq(x, llc::LT, xval); + lemma |= ineq(y, llc::GT, yval); + lemma |= ineq(r, llc::LE, r2val); + return l_false; + } + } + if (xval > 0 && yval > 0 && !yval.is_int()) { + auto ynum = numerator(yval); + auto yden = denominator(yval); + if (!ynum.is_unsigned()) + return l_undef; + if (!yden.is_unsigned()) + return l_undef; + // r = x^{yn/yd} + // <=> + // r^yd = x^yn + auto ryd = power(rval, yden.get_unsigned()); + auto xyn = power(xval, ynum.get_unsigned()); + if (ryd == xyn) + return l_true; + } + + return l_undef; + + } + +} diff --git a/src/math/lp/nla_powers.h b/src/math/lp/nla_powers.h new file mode 100644 index 000000000..f74417ae3 --- /dev/null +++ b/src/math/lp/nla_powers.h @@ -0,0 +1,29 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + nla_powers.h + +Author: + Lev Nachmanson (levnach) + Nikolaj Bjorner (nbjorner) + +Description: + Refines bounds on powers. + +--*/ + +#include "math/lp/nla_types.h" + +namespace nla { + + class core; + + class powers { + core& m_core; + public: + powers(core& c):m_core(c) {} + lbool check(lpvar r, lpvar x, lpvar y, vector&); + }; +} diff --git a/src/math/lp/nla_solver.cpp b/src/math/lp/nla_solver.cpp index e0df39503..b5a030a92 100644 --- a/src/math/lp/nla_solver.cpp +++ b/src/math/lp/nla_solver.cpp @@ -77,104 +77,7 @@ namespace nla { // ensure r = x^y, add abstraction/refinement lemmas lbool solver::check_power(lpvar r, lpvar x, lpvar y, vector& lemmas) { - if (x == null_lpvar || y == null_lpvar || r == null_lpvar) - return l_undef; - - if (use_nra_model()) - return l_undef; - - auto xval = m_core->val(x); - auto yval = m_core->val(y); - auto rval = m_core->val(r); - - core& c = get_core(); - c.set_lemma_vec(lemmas); - lemmas.reset(); - - // x >= x0 > 0, y >= y0 > 0 => r >= x0^y0 - // x >= x0 > 0, y <= y0 => r <= x0^y0 - // x != 0, y = 0 => r = 1 - // x = 0, y != 0 => r = 0 - // - // for x fixed, it is exponentiation - // => use tangent lemmas and error tolerance. - - if (xval > 0 && yval.is_unsigned()) { - auto r2val = power(xval, yval.get_unsigned()); - if (rval == r2val) - return l_true; - if (xval != 0 && yval == 0) { - new_lemma lemma(c, "x != 0 => x^0 = 1"); - lemma |= ineq(x, llc::EQ, rational::zero()); - lemma |= ineq(y, llc::NE, rational::zero()); - lemma |= ineq(r, llc::EQ, rational::one()); - return l_false; - } - if (xval == 0 && yval > 0) { - new_lemma lemma(c, "y != 0 => 0^y = 0"); - lemma |= ineq(x, llc::NE, rational::zero()); - lemma |= ineq(y, llc::EQ, rational::zero()); - lemma |= ineq(r, llc::EQ, rational::zero()); - return l_false; - } - if (xval > 0 && r2val < rval) { - SASSERT(yval > 0); - new_lemma lemma(c, "x >= x0 > 0, y >= y0 > 0 => r >= x0^y0"); - lemma |= ineq(x, llc::LT, xval); - lemma |= ineq(y, llc::LT, yval); - lemma |= ineq(r, llc::GE, r2val); - return l_false; - } - if (xval > 0 && r2val < rval) { - new_lemma lemma(c, "x >= x0 > 0, y <= y0 => r <= x0^y0"); - lemma |= ineq(x, llc::LT, xval); - lemma |= ineq(y, llc::GT, yval); - lemma |= ineq(r, llc::LE, r2val); - return l_false; - } - } - if (xval > 0 && yval > 0 && !yval.is_int()) { - auto ynum = numerator(yval); - auto yden = denominator(yval); - if (!ynum.is_unsigned()) - return l_undef; - if (!yden.is_unsigned()) - return l_undef; - // r = x^{yn/yd} - // <=> - // r^yd = x^yn - auto ryd = power(rval, yden.get_unsigned()); - auto xyn = power(xval, ynum.get_unsigned()); - if (ryd == xyn) - return l_true; -#if 0 - // try some root approximation instead? - if (ryd > xyn) { - // todo - } - if (ryd < xyn) { - // todo - } -#endif - - } - - - return l_undef; - - // anum isn't initialized unless nra_solver is invoked. - // there is no proviso for using algebraic numbers outside of the nra solver. - // so either we have a rational refinement version _and_ an algebraic numeral refinement - // loop or we introduce algebraic numerals outside of the nra_solver - -#if 0 - scoped_anum xval(am()), yval(am()), rval(am()); - - am().set(xval, am_value(x)); - am().set(yval, am_value(y)); - am().set(rval, am_value(r)); -#endif - + return m_core->check_power(r, x, y, lemmas); } } diff --git a/src/math/lp/nla_types.h b/src/math/lp/nla_types.h new file mode 100644 index 000000000..8169266cc --- /dev/null +++ b/src/math/lp/nla_types.h @@ -0,0 +1,120 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + nla_types.h + +Author: + Lev Nachmanson (levnach) + Nikolaj Bjorner (nbjorner) + +Description: + Types used for nla solver. + +--*/ + +#pragma once + +namespace nla { + + typedef lp::constraint_index lpci; + typedef lp::lconstraint_kind llc; + typedef lp::constraint_index lpci; + typedef lp::explanation expl_set; + typedef lp::var_index lpvar; + const lpvar null_lpvar = UINT_MAX; + + inline int rat_sign(const rational& r) { return r.is_pos()? 1 : ( r.is_neg()? -1 : 0); } + inline rational rrat_sign(const rational& r) { return rational(rat_sign(r)); } + inline bool is_set(unsigned j) { return j != null_lpvar; } + inline bool is_even(unsigned k) { return (k & 1) == 0; } + class ineq { + lp::lconstraint_kind m_cmp; + lp::lar_term m_term; + rational m_rs; + public: + ineq(lp::lconstraint_kind cmp, const lp::lar_term& term, const rational& rs) : m_cmp(cmp), m_term(term), m_rs(rs) {} + ineq(const lp::lar_term& term, lp::lconstraint_kind cmp, int i) : m_cmp(cmp), m_term(term), m_rs(rational(i)) {} + ineq(const lp::lar_term& term, lp::lconstraint_kind cmp, const rational& rs) : m_cmp(cmp), m_term(term), m_rs(rs) {} + ineq(lpvar v, lp::lconstraint_kind cmp, int i): m_cmp(cmp), m_term(v), m_rs(rational(i)) {} + ineq(lpvar v, lp::lconstraint_kind cmp, rational const& r): m_cmp(cmp), m_term(v), m_rs(r) {} + bool operator==(const ineq& a) const { + return m_cmp == a.m_cmp && m_term == a.m_term && m_rs == a.m_rs; + } + const lp::lar_term& term() const { return m_term; }; + lp::lconstraint_kind cmp() const { return m_cmp; }; + const rational& rs() const { return m_rs; }; + }; + + class lemma { + vector m_ineqs; + lp::explanation m_expl; + public: + void push_back(const ineq& i) { m_ineqs.push_back(i);} + size_t size() const { return m_ineqs.size() + m_expl.size(); } + const vector& ineqs() const { return m_ineqs; } + vector& ineqs() { return m_ineqs; } + lp::explanation& expl() { return m_expl; } + const lp::explanation& expl() const { return m_expl; } + bool is_conflict() const { return m_ineqs.empty() && !m_expl.empty(); } + }; + + class core; + // + // lemmas are created in a scope. + // when the destructor of new_lemma is invoked + // all constraints are assumed added to the lemma + // correctness of the lemma can be checked at this point. + // + class new_lemma { + char const* name; + core& c; + lemma& current() const; + + public: + new_lemma(core& c, char const* name); + ~new_lemma(); + lemma& operator()() { return current(); } + std::ostream& display(std::ostream& out) const; + new_lemma& operator&=(lp::explanation const& e); + new_lemma& operator&=(const monic& m); + new_lemma& operator&=(const factor& f); + new_lemma& operator&=(const factorization& f); + new_lemma& operator&=(lpvar j); + new_lemma& operator|=(ineq const& i); + new_lemma& explain_fixed(lpvar j); + new_lemma& explain_equiv(lpvar u, lpvar v); + new_lemma& explain_var_separated_from_zero(lpvar j); + new_lemma& explain_existing_lower_bound(lpvar j); + new_lemma& explain_existing_upper_bound(lpvar j); + + lp::explanation& expl() { return current().expl(); } + + unsigned num_ineqs() const { return current().ineqs().size(); } + }; + + + inline std::ostream& operator<<(std::ostream& out, new_lemma const& l) { + return l.display(out); + } + + struct pp_fac { + core const& c; + factor const& f; + pp_fac(core const& c, factor const& f): c(c), f(f) {} + }; + + struct pp_var { + core const& c; + lpvar v; + pp_var(core const& c, lpvar v): c(c), v(v) {} + }; + + struct pp_factorization { + core const& c; + factorization const& f; + pp_factorization(core const& c, factorization const& f): c(c), f(f) {} + }; + +}