3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-05 17:14:07 +00:00

extend monomial bounds to handle powers

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2020-05-14 19:13:17 -07:00
parent 73fa5995d4
commit b43ed70874
4 changed files with 151 additions and 99 deletions

View file

@ -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<dep_intervals::with_deps>(vi, power, vi);
dep.mul<dep_intervals::with_deps>(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<dep_intervals::with_deps>(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<dep_intervals::with_deps>(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<dep_intervals::with_deps>(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<dep_intervals::with_deps>(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;

View file

@ -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;

View file

@ -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();
}
}

View file

@ -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;