z3-z3-4.13.0.src.smt.theory_arith_nl.h Maven / Gradle / Ivy
The newest version!
/*++
Copyright (c) 2006 Microsoft Corporation
Module Name:
theory_arith_nl.h
Abstract:
Author:
Leonardo de Moura (leonardo) 2008-12-08.
Revision History:
--*/
#pragma once
#include "ast/ast_smt2_pp.h"
namespace smt {
template
expr * theory_arith::mk_nary_mul(unsigned sz, expr * const * args, bool is_int) {
if (sz == 0)
return m_util.mk_numeral(rational(1), is_int);
if (sz == 1)
return args[0];
if (sz == 2)
return m_util.mk_mul(args[0], args[1]);
if (m_util.is_numeral(args[0]))
return m_util.mk_mul(args[0], m_util.mk_mul(sz - 1, args + 1));
return m_util.mk_mul(sz, args);
}
template
expr * theory_arith::mk_nary_add(unsigned sz, expr * const * args, bool is_int) {
if (sz == 0)
return m_util.mk_numeral(rational(0), is_int);
if (sz == 1)
return args[0];
return m_util.mk_add(sz, args);
}
template
expr * theory_arith::mk_nary_add(unsigned sz, expr * const * args) {
SASSERT(sz != 0);
return mk_nary_add(sz, args, false);
}
/**
\brief Insert v into vars and already_found if v is not already in already_found.
*/
template
void theory_arith::mark_var(theory_var v, svector & vars, var_set & already_found) {
if (already_found.contains(v))
return;
already_found.insert(v);
vars.push_back(v);
}
/**
\brief Invoke mark_var for all variables in rows that contain v.
*/
template
void theory_arith::mark_dependents(theory_var v, svector & vars, var_set & already_found, row_set & already_visited_rows) {
if (is_pure_monomial(v)) {
expr * n = var2expr(v);
SASSERT(m_util.is_mul(n));
for (expr * curr : *to_app(n)) {
if (ctx.e_internalized(curr)) {
theory_var v = expr2var(curr);
SASSERT(v != null_theory_var);
mark_var(v, vars, already_found);
}
}
}
if (is_fixed(v))
return;
column & c = m_columns[v];
for (auto& ce : c) {
if (ce.is_dead() || already_visited_rows.contains(ce.m_row_id))
continue;
TRACE("non_linear_bug", tout << "visiting row: " << ce.m_row_id << "\n";);
already_visited_rows.insert(ce.m_row_id);
row & r = m_rows[ce.m_row_id];
theory_var s = r.get_base_var();
// ignore quasi base vars... actually they should not be used if the problem is non linear...
if (is_quasi_base(s))
continue;
// If s is a base variable different from v and it is free, then this row can be ignored.
// It doesn't need to be part of the non linear cluster. For all purposes, this variable
// was eliminated by substitution.
if (s != null_theory_var && is_free(s) && s != v)
continue;
for (auto& re : r) {
if (!re.is_dead() && !is_fixed(re.m_var))
mark_var(re.m_var, vars, already_found);
CTRACE("non_linear", !re.is_dead() && is_fixed(re.m_var), tout << "skipped fixed\n");
}
}
}
/**
\brief Store in vars the variables that are in the non linear cluster of constraints,
and are not satisfied by the current assignment.
*/
template
void theory_arith::get_non_linear_cluster(svector & vars) {
if (m_nl_monomials.empty())
return;
var_set already_found;
row_set already_visited_rows;
for (theory_var v : m_nl_monomials) {
expr * n = var2expr(v);
if (ctx.is_relevant(n))
mark_var(v, vars, already_found);
}
// NB: vars changes inside of loop
for (unsigned idx = 0; idx < vars.size(); ++idx) {
theory_var v = vars[idx];
TRACE("non_linear", tout << "marking dependents of: v" << v << "\n";);
mark_dependents(v, vars, already_found, already_visited_rows);
}
TRACE("non_linear", tout << "variables in non linear cluster: ";
for (theory_var v : vars) tout << "v" << v << " "; tout << "\n";
for (theory_var v : m_nl_monomials) tout << "non-linear v" << v << " " << mk_pp(var2expr(v), m) << "\n";);
}
/**
\brief Return the number of variables that
do not have bounds associated with it.
The result is 0, 1, or 2. The value 2 means "2 or more".
The second value is the idx of the variable that does not
have bounds associated with it. It is only useful when the first value is 1.
The second value is -1 if such variable does not exist, that is, the first
value is 0.
\remark if a variables has an even number of occurrences, then
I consider that it has a bound associated with it.
Examples:
1) Assume x1, x4 have bounds:
analyze_monomial(x1 * x2 * x2 * x3 * x3 * x3 * x4)
-->
(1,2)
Explanation: x2 doesn't have bounds, but x2 has an even power.
So x2*x2 has bound [0, oo). So, there is one variable without bounds x3.
It is the third variable in the monomial, then its idx is 2.
*/
template
typename theory_arith::n_var_power_pair theory_arith::analyze_monomial(expr * m) const {
buffer vp;
decompose_monomial(m, vp);
unsigned c = 0;
var_power_pair q(nullptr, 0);
for (auto const& p : vp) {
if (p.second % 2 == 1 && is_free(p.first)) {
c++;
q = p;
if (c > 1)
break;
}
}
return std::make_pair(c, q);
}
template
rational theory_arith::decompose_monomial(expr* m, buffer& vp) const {
rational coeff(1);
vp.reset();
expr_fast_mark1 mark;
auto insert = [&](expr* arg) {
rational r;
if (m_util.is_numeral(arg, r))
coeff *= r;
else if (mark.is_marked(arg)) {
for (unsigned i = vp.size(); i-- > 0; ) {
if (vp[i].first == arg) {
vp[i].second++;
break;
}
}
}
else {
mark.mark(arg);
vp.push_back(var_power_pair(arg, 1));
}
};
while (m_util.is_mul(m)) {
unsigned sz = to_app(m)->get_num_args();
for (unsigned i = 0; i + 1 < sz; ++i) {
insert(to_app(m)->get_arg(i));
}
m = to_app(m)->get_arg(sz - 1);
}
insert(m);
return coeff;
}
/**
\brief Return an interval using the bounds for v.
*/
template
interval theory_arith::mk_interval_for(theory_var v) {
bound * l = lower(v);
bound * u = upper(v);
if (l && u) {
// optimization may introduce non-standard bounds.
if (l->get_value() == u->get_value() && !l->get_value().get_infinitesimal().to_rational().is_zero()) {
return interval(m_dep_manager);
}
return interval(m_dep_manager,
l->get_value().get_rational().to_rational(),
l->get_value().get_infinitesimal().to_rational().is_pos(),
m_dep_manager.mk_leaf(l),
u->get_value().get_rational().to_rational(),
u->get_value().get_infinitesimal().to_rational().is_neg(),
m_dep_manager.mk_leaf(u));
}
else if (l) {
return interval(m_dep_manager,
l->get_value().get_rational().to_rational(),
l->get_value().get_infinitesimal().to_rational().is_pos(),
true,
m_dep_manager.mk_leaf(l));
}
else if (u) {
return interval(m_dep_manager,
u->get_value().get_rational().to_rational(),
u->get_value().get_infinitesimal().to_rational().is_neg(),
false,
m_dep_manager.mk_leaf(u));
}
else {
return interval(m_dep_manager);
}
}
/**
\brief Return an interval for the given expression using its bounds.
*/
template
interval theory_arith::mk_interval_for(expr * n) {
if (!has_var(n))
return interval(m_dep_manager);
return mk_interval_for(expr2var(n));
}
/**
\brief target *= [lower(var), upper(var)]^power
*/
template
void theory_arith::mul_bound_of(expr * var, unsigned power, interval & target) {
theory_var v = expr2var(var);
interval i = mk_interval_for(v);
TRACE("non_linear",
display_interval(tout << "bound: ",i); tout << i << "\n";
tout << mk_pp(var, get_manager()) << "\n";
tout << "power " << power << ": " << expt(i, power) << "\n";
display_interval(tout << "target before: ", target); tout << "\n";);
i.expt(power);
target *= i;
get_manager().limit().inc((target.is_lower_open() || target.minus_infinity()) ? 1 : target.get_lower_value().bitsize());
get_manager().limit().inc((target.is_upper_open() || target.plus_infinity()) ? 1 : target.get_upper_value().bitsize());
TRACE("non_linear", display_interval(tout << "target after: ", target); tout << "\n";);
}
/**
\brief Evaluate the given expression using interval arithmetic.
- If a subexpression is internalized, then mk_interval_for is used to
compute its interval.
- Only +, *, and numerals are handled.
*/
template
interval theory_arith::evaluate_as_interval(expr * n) {
expr* arg;
rational val;
#define TR() TRACE("nl_evaluate", tout << "eval: " << mk_bounded_pp(n, get_manager(), 10) << "\n";\
display_nested_form(tout, n); tout << "\ninterval: " << r << "\n";);
if (has_var(n)) {
interval r = mk_interval_for(n);
TR();
return r;
}
else if (m_util.is_add(n)) {
interval r(m_dep_manager, rational(0));
for (expr* arg : *to_app(n)) {
r += evaluate_as_interval(arg);
}
TR();
return r;
}
else if (m_util.is_mul(n)) {
buffer vp;
rational coeff = decompose_monomial(n, vp);
interval r(m_dep_manager, coeff);
for (var_power_pair const& p : vp) {
expr * var = p.first;
unsigned power = p.second;
interval it = evaluate_as_interval(var);
it.expt(power);
r *= it;
}
TR();
return r;
}
else if (m_util.is_to_real(n, arg)) {
return evaluate_as_interval(arg);
}
else if (m_util.is_numeral(n, val)) {
interval r = interval(m_dep_manager, val);
TR();
return r;
}
else {
TRACE("nl_evaluate", tout << "is unknown\n";);
return interval(m_dep_manager);
}
}
template
void theory_arith::display_monomial(std::ostream & out, expr * n) const {
bool first = true;
buffer vp;
rational coeff = decompose_monomial(n, vp);
if (!coeff.is_one()) {
out << coeff;
first = false;
}
for (auto const& p : vp) {
SASSERT(p.first != 0);
if (first) first = false; else out << " * ";
out << mk_bounded_pp(p.first, get_manager()) << "^" << p.second;
}
}
template
void theory_arith::dependency2new_bound(v_dependency * dep, derived_bound& new_bound) {
ptr_vector bounds;
m_dep_manager.linearize(dep, bounds);
m_tmp_lit_set.reset();
m_tmp_eq_set.reset();
for (void* _b : bounds) {
bound * b = static_cast(_b);
accumulate_justification(*b, new_bound, numeral::zero(), m_tmp_lit_set, m_tmp_eq_set);
}
}
/**
\brief Create a new derived bound. The justification is stored in the object dep.
*/
template
void theory_arith::mk_derived_nl_bound(theory_var v, inf_numeral const & coeff, bound_kind k, v_dependency * dep) {
inf_numeral coeff_norm = normalize_bound(v, coeff, k);
derived_bound * new_bound = alloc(derived_bound, v, coeff_norm, k);
m_bounds_to_delete.push_back(new_bound);
m_asserted_bounds.push_back(new_bound);
// copy justification to new bound
dependency2new_bound(dep, *new_bound);
TRACE("non_linear", new_bound->display(*this, tout); tout << "\n";);
}
/**
\brief Update the bounds of v, using the interval i.
Return true if i improves the bounds of v.
*/
template
bool theory_arith::update_bounds_using_interval(theory_var v, interval const & i) {
SASSERT(v != null_theory_var);
bool r = false;
if (!i.minus_infinity()) {
inf_numeral new_lower(i.get_lower_value());
if (i.is_lower_open()) {
if (is_int(v)) {
if (new_lower.is_int()) {
new_lower += rational::one();
}
else {
new_lower = ceil(new_lower.get_rational());
}
}
else {
new_lower += get_epsilon(v);
}
}
bound * old_lower = lower(v);
if (old_lower == nullptr || new_lower > old_lower->get_value()) {
TRACE("non_linear", tout << "NEW lower bound for v" << v << " " << mk_pp(var2expr(v), get_manager())
<< " " << new_lower << "\n";
display_interval(tout, i); tout << "\n";);
mk_derived_nl_bound(v, new_lower, B_LOWER, i.get_lower_dependencies());
r = true;
}
}
if (!i.plus_infinity()) {
inf_numeral new_upper(i.get_upper_value());
if (i.is_upper_open()) {
if (is_int(v)) {
if (new_upper.is_int()) {
new_upper -= rational::one();
}
else {
new_upper = floor(new_upper.get_rational());
}
}
else {
new_upper -= get_epsilon(v);
}
}
bound * old_upper = upper(v);
if (old_upper == nullptr || new_upper < old_upper->get_value()) {
TRACE("non_linear", tout << "NEW upper bound for v" << v << " " << new_upper << "\n";
display_interval(tout, i); tout << "\n";);
mk_derived_nl_bound(v, new_upper, B_UPPER, i.get_upper_dependencies());
r = true;
}
}
return r;
}
template
bool theory_arith::update_bounds_using_interval(expr * n, interval const & i) {
SASSERT(expr2var(n) != null_theory_var);
TRACE("non_linear", tout << "NL bounds for m: " << i << "\n" << mk_pp(n, get_manager()) << "\n";);
return update_bounds_using_interval(expr2var(n), i);
}
/**
\brief Use the bounds of the variables to build a bound for m.
*/
template
bool theory_arith::propagate_nl_upward(expr * m) {
SASSERT(is_pure_monomial(m));
TRACE("nl_arith_bug", tout << "processing upward:\n" << mk_pp(m, get_manager()) << "\n";);
buffer vp;
rational coeff = decompose_monomial(m, vp);
interval new_bounds(m_dep_manager, coeff);
for (var_power_pair const& p : vp) {
expr * var = p.first;
unsigned power = p.second;
TRACE("nl_arith_bug", tout << "interval before: " << new_bounds << "\n";
theory_var v = expr2var(var);
interval i = mk_interval_for(v);
display_var(tout, v);
tout << "interval for v" << i << " " << mk_pp(var, get_manager()) << "\npower: " << power << " " << expt(i, power) << "\n";);
mul_bound_of(var, power, new_bounds);
TRACE("nl_arith_bug", tout << "interval after: " << new_bounds << "\n";);
}
return update_bounds_using_interval(m, new_bounds);
}
/**
\brief Propagate a bound to the i-th variable of the given monomial
using the bounds of m and other variables in m.
\remark We do not support roots in interval... so, if the i-th var has power != 1
the method returns without doing anything.
*/
template
bool theory_arith::propagate_nl_downward(expr * n, var_power_pair const& p) {
SASSERT(is_pure_monomial(n));
expr * v = p.first;
unsigned power = p.second;
if (power != 1)
return false; // TODO: remove, when the n-th root is implemented in interval.
buffer vp;
rational coeff = decompose_monomial(n, vp);
interval other_bounds(m_dep_manager, coeff);
// TODO: the following code can be improved it is quadratic on the degree of the monomial.
for (var_power_pair const& p : vp) {
if (p.first == v)
continue;
expr * var = p.first;
unsigned power = p.second;
mul_bound_of(var, power, other_bounds);
}
if (other_bounds.contains_zero())
return false; // interval division requires that divisor doesn't contain 0.
interval r = mk_interval_for(n);
TRACE("nl_arith_bug", tout << "m: " << mk_ismt2_pp(n, get_manager()) << "\nv: " << mk_ismt2_pp(v, get_manager()) <<
"\npower: " << power << "\n";
display_interval(tout << "monomial bounds\n", r);
display_interval(tout << "other bounds\n", other_bounds);
);
r /= other_bounds;
return update_bounds_using_interval(v, r);
}
/**
\brief The given monomial and its elements have bounds.
Propagate bounds to all of them.
Return true if some bound was propagated.
*/
template
bool theory_arith::propagate_nl_bounds(expr * m) {
TRACE("non_linear", tout << "propagate several bounds using:\n"; display_monomial(tout, m); tout << "\n";);
bool result = propagate_nl_upward(m);
buffer vp;
rational coeff = decompose_monomial(m, vp);
for (auto const& p : vp) {
if (propagate_nl_downward(m, p)) {
m_stats.m_nl_bounds++;
result = true;
}
}
return result;
}
/**
\brief Try to propagate bounds using non linear monomials.
Return true if some bound was propagated.
*/
template
bool theory_arith::propagate_nl_bounds() {
m_dep_manager.reset();
bool propagated = false;
for (unsigned i = 0; i < m_nl_monomials.size(); i++) {
theory_var v = m_nl_monomials[i];
expr * m = var2expr(v);
if (!ctx.is_relevant(m))
continue;
auto p = analyze_monomial(m);
TRACE("propagate_nl_bound", tout << "m: " << mk_ismt2_pp(m, get_manager()) << "\n" << "p: " << p.first << "\n";);
unsigned num_bad_vars = p.first;
var_power_pair q = p.second;
SASSERT(num_bad_vars != 1 || q.first != nullptr);
if (num_bad_vars >= 2)
continue;
bool is_free_m = is_free(m);
TRACE("propagate_nl_bound", tout << "is_free_m: " << is_free_m << "\n";);
if (num_bad_vars == 1 && is_free_m)
continue;
if (num_bad_vars == 0) {
if (!is_free_m) {
if (propagate_nl_bounds(m))
propagated = true;
}
else {
if (propagate_nl_upward(m)) {
m_stats.m_nl_bounds++;
propagated = true;
}
}
}
else {
SASSERT (!is_free_m);
if (propagate_nl_downward(m, q)) {
m_stats.m_nl_bounds++;
propagated = true;
}
}
}
return propagated;
}
/**
\brief Return the value of v as a rational. If computed_epsilon = false and v has an infinitesimal, then
compute_epsilon() is invoked.
*/
template
rational theory_arith::get_value(theory_var v, bool & computed_epsilon) {
inf_numeral const & val = get_value(v);
if (!val.get_infinitesimal().is_zero() && !computed_epsilon) {
compute_epsilon();
refine_epsilon();
computed_epsilon = true;
m_model_depends_on_computed_epsilon = true;
}
return val.get_rational().to_rational() + m_epsilon.to_rational() * val.get_infinitesimal().to_rational();
}
/**
\brief Return true if for the monomial x_1 * ... * x_n associated with v,
the following holds:
get_value(x_1) * ... * get_value(x_n) = get_value(v)
*/
template
bool theory_arith::check_monomial_assignment(theory_var v, bool & computed_epsilon) {
SASSERT(is_pure_monomial(var2expr(v)));
expr * m = var2expr(v);
rational val(1), v_val;
for (expr* arg : *to_app(m)) {
theory_var curr = expr2var(arg);
SASSERT(curr != null_theory_var);
v_val = get_value(curr, computed_epsilon);
TRACE("non_linear", tout << mk_pp(arg, get_manager()) << " = " << v_val << "\n";);
val *= v_val;
}
v_val = get_value(v, computed_epsilon);
TRACE("non_linear", tout << "v" << v << " := " << v_val << " == " << val << "\n";);
return v_val == val;
}
/**
\brief Return true if for every monomial x_1 * ... * x_n,
get_value(x_1) * ... * get_value(x_n) = get_value(x_1 * ... * x_n)
*/
template
bool theory_arith::check_monomial_assignments() {
bool computed_epsilon = false;
for (theory_var v : m_nl_monomials) {
TRACE("non_linear", tout << "v" << v << " is relevant: " << ctx.is_relevant(get_enode(v)) << "\n");
if (ctx.is_relevant(get_enode(v)) && !check_monomial_assignment(v, computed_epsilon)) {
TRACE("non_linear", tout << "check_monomial_assignment failed for:\n" << mk_ismt2_pp(var2expr(v), get_manager()) << "\n";
display_var(tout, v););
return false;
}
}
return true;
}
/**
\brief Try to find an integer variable for performing branching
in the non linear cluster.
The idea is select a variable in a monomial with an invalid
assignment. I give preference to variables with small ranges.
If no variable is bounded, then select a random one.
Free variables are not considered.
*/
template
theory_var theory_arith::find_nl_var_for_branching() {
TRACE("nl_branching", tout << "looking for variable to branch...\n"; display(tout););
theory_var target = null_theory_var;
bool bounded = false;
unsigned n = 0;
numeral range;
for (unsigned j = 0; j < m_nl_monomials.size(); ++j) {
theory_var v = m_nl_monomials[j];
if (is_real(v))
continue;
bool computed_epsilon = false;
bool r = check_monomial_assignment(v, computed_epsilon);
if (!r) {
expr * m = get_enode(v)->get_expr();
SASSERT(is_pure_monomial(m));
for (expr * arg : *to_app(m)) {
theory_var curr = ctx.get_enode(arg)->get_th_var(get_id());
TRACE("nl_branching", tout << "target: v" << target << ", curr: v" << curr << "\n";);
if (!is_fixed(curr) && is_int(curr)) {
if (is_bounded(curr)) {
numeral new_range;
new_range = upper_bound(curr).get_rational();
new_range -= lower_bound(curr).get_rational();
if (!bounded || new_range < range) {
target = curr;
range = new_range;
bounded = true;
}
}
else if (!bounded) {
n++;
TRACE("nl_branching", tout << "n: " << n << "\n";);
if (m_random()%n == 0)
target = curr;
SASSERT(target != null_theory_var);
}
SASSERT(target != null_theory_var);
}
TRACE("nl_branching", tout << "after target: v" << target << "\n";);
}
}
}
return target;
}
/**
\brief Branch on an integer variable. This method is invoked when v is part
of a non linear monomial that is not satisfied by the current assignment.
if v >= l, then create the case split v >= l+1
else v <= u, then create the case split v <= u-1
else create the bound v = 0 and case split on it.
*/
template
bool theory_arith::branch_nl_int_var(theory_var v) {
TRACE("non_linear", tout << "BRANCHING on v" << v << "\n";);
m_stats.m_nl_branching++;
SASSERT(is_int(v));
expr_ref bound(get_manager());
if (lower(v))
bound = m_util.mk_le(var2expr(v), m_util.mk_numeral(lower_bound(v).get_rational().to_rational(), true));
else if (upper(v))
bound = m_util.mk_ge(var2expr(v), m_util.mk_numeral(upper_bound(v).get_rational().to_rational(), true));
else
bound = m_util.mk_eq(var2expr(v), m_util.mk_numeral(rational(0), true));
TRACE("non_linear", tout << "new bound:\n" << mk_pp(bound, get_manager()) << "\n";);
ast_manager & m = get_manager();
{
std::function fn = [&]() { return m.mk_or(bound, m.mk_not(bound)); };
scoped_trace_stream _sts(*this, fn);
ctx.internalize(bound, true);
}
ctx.mark_as_relevant(bound.get());
literal l = ctx.get_literal(bound);
SASSERT(!l.sign());
ctx.set_true_first_flag(l.var()); // force the context to case split to true first, independently of the phase selection strategy.
return true;
}
/**
\brief Return true if the given monomial is linear.
*/
template
bool theory_arith::is_monomial_linear(expr * m) const {
SASSERT(is_pure_monomial(m));
unsigned num_nl_vars = 0;
for (expr* arg : *to_app(m)) {
if (!ctx.e_internalized(arg))
return false;
theory_var _var = expr2var(arg);
CTRACE("non_linear", is_fixed(_var),
tout << mk_pp(arg, get_manager()) << " is fixed: " << lower_bound(_var) << "\n";);
if (!is_fixed(_var)) {
num_nl_vars++;
}
else if (lower_bound(_var).is_zero()) {
return true;
}
}
return num_nl_vars <= 1;
}
/**
\brief Return the product of the value of the fixed variables in the
monomial m.
*/
template
typename theory_arith::numeral theory_arith::get_monomial_fixed_var_product(expr * m) const {
SASSERT(is_pure_monomial(m));
numeral r(1);
for (expr * arg : *to_app(m)) {
theory_var _var = expr2var(arg);
if (is_fixed(_var))
r *= lower_bound(_var).get_rational();
}
TRACE("arith", tout << mk_pp(m, get_manager()) << " " << r << "\n";);
return r;
}
/**
\brief Return the first non fixed variable in the given monomial.
Return 0, if the monomial does not have a non fixed variable.
*/
template
expr * theory_arith::get_monomial_non_fixed_var(expr * m) const {
SASSERT(is_pure_monomial(m));
for (expr* arg : *to_app(m)) {
if (!is_fixed(expr2var(arg)))
return arg;
}
return nullptr;
}
/**
\brief Propagate linear monomial. Check whether the give
monomial became linear and propagate.
*/
template
bool theory_arith::propagate_linear_monomial(theory_var v) {
TRACE("non_linear_verbose", tout << "checking whether v" << v << " became linear...\n";);
if (m_data[v].m_nl_propagated)
return false; // already propagated this monomial.
expr * m = var2expr(v);
if (!is_monomial_linear(m))
return false; // monomial is not linear.
m_stats.m_nl_linear++;
m_data[v].m_nl_propagated = true;
m_nl_propagated.push_back(v);
TRACE("non_linear", tout << "v" << v << " is linear " << mk_pp(m, get_manager()) << "\n";);
numeral k = get_monomial_fixed_var_product(m);
TRACE("non_linear", tout << "new linear monomial... k: " << k << "\n";);
expr * x_n = k.is_zero() ? nullptr : get_monomial_non_fixed_var(m);
TRACE("non_linear_bug", if (x_n != 0) { tout << "x_n: " << mk_bounded_pp(x_n, get_manager()) << "\nx_n: #" << x_n->get_id() << "\n"; });
derived_bound * new_lower = nullptr;
derived_bound * new_upper = nullptr;
if (x_n != nullptr) {
// All but one of the x_i variables are assigned.
// Let x_n be the unassigned variable.
// Then, we know that x_1*...*x_n = k*x_n, where k is the product of beta(x_1)*...*beta(x_{n-1})
// beta(x_i) == lower(x_i)
// Let m be (* x_1 ... x_n), then assert equality
// (= (+ (* x_1 ... x_n) (* -k x_n)) 0) when x_1 ... x_{n-1} are fixed variables.
// where k = lower(x_1)*...*lower(x_{n-1})
TRACE("non_linear", tout << "x_n: " << mk_pp(x_n, get_manager()) << "\n";);
k.neg();
expr * k_x_n = k.is_one() ? x_n : m_util.mk_mul(m_util.mk_numeral(k.to_rational(), is_int(v)), x_n);
expr * rhs = m_util.mk_add(m, k_x_n);
TRACE("non_linear_bug", tout << "rhs: " << mk_bounded_pp(rhs, get_manager(),5) << "\ninternalized: " << ctx.e_internalized(rhs) << "\n";);
if (!has_var(rhs)) {
ctx.internalize(rhs, false);
ctx.mark_as_relevant(rhs);
}
TRACE("non_linear_bug", tout << "enode: " << ctx.get_enode(rhs) << " enode_id: " << ctx.get_enode(rhs)->get_owner_id() << "\n";);
IF_VERBOSE(3,
for (auto* arg : *to_app(m))
if (is_fixed(expr2var(arg)))
verbose_stream() << mk_pp(arg, get_manager()) << " = " << -k << "\n");
theory_var new_v = expr2var(rhs);
SASSERT(new_v != null_theory_var);
new_lower = alloc(derived_bound, new_v, inf_numeral(0), B_LOWER);
new_upper = alloc(derived_bound, new_v, inf_numeral(0), B_UPPER);
}
else {
// One of the x_i variables is zero,
// or all of them are assigned.
// Assert the equality
// (= (* x_1 ... x_n) k)
TRACE("non_linear", tout << "all variables are fixed, and bound is: " << k << "\n";);
new_lower = alloc(derived_bound, v, inf_numeral(k), B_LOWER);
new_upper = alloc(derived_bound, v, inf_numeral(k), B_UPPER);
}
SASSERT(new_lower);
SASSERT(new_upper);
m_bounds_to_delete.push_back(new_lower);
m_asserted_bounds.push_back(new_lower);
m_bounds_to_delete.push_back(new_upper);
m_asserted_bounds.push_back(new_upper);
// Add the justification for new_lower and new_upper.
// The justification is the lower and upper bounds of all fixed variables.
m_tmp_lit_set.reset();
m_tmp_eq_set.reset();
SASSERT(is_pure_monomial(m));
bool found_zero = false;
for (unsigned i = 0; !found_zero && i < to_app(m)->get_num_args(); i++) {
expr * arg = to_app(m)->get_arg(i);
theory_var _var = expr2var(arg);
if (is_fixed(_var)) {
bound * l = lower(_var);
bound * u = upper(_var);
if (l->get_value().is_zero()) {
/* if zero was found, then it is the explanation */
SASSERT(k.is_zero());
found_zero = true;
m_tmp_lit_set.reset();
m_tmp_eq_set.reset();
new_lower->m_lits.reset();
new_lower->m_eqs.reset();
}
accumulate_justification(*l, *new_lower, numeral::zero(), m_tmp_lit_set, m_tmp_eq_set);
TRACE("non_linear",
for (literal l : new_lower->m_lits) {
ctx.display_detailed_literal(tout, l) << " ";
}
tout << "\n";);
accumulate_justification(*u, *new_lower, numeral::zero(), m_tmp_lit_set, m_tmp_eq_set);
TRACE("non_linear",
for (literal l : new_lower->m_lits) {
ctx.display_detailed_literal(tout, l) << " ";
}
tout << "\n";);
}
}
new_upper->m_lits.append(new_lower->m_lits);
new_upper->m_eqs.append(new_lower->m_eqs);
TRACE("non_linear",
new_lower->display(*this, tout << "lower: "); tout << "\n";
new_upper->display(*this, tout << "upper: "); tout << "\n";
for (literal lit : new_upper->m_lits) {
ctx.display_detailed_literal(tout, lit) << " ";
}
tout << "\n";);
return true;
}
/**
\brief Traverse all non linear monomials, and check the ones that became
linear and propagate. Return true if propagated.
*/
template
bool theory_arith::propagate_linear_monomials() {
if (!m_params.m_nl_arith_propagate_linear_monomials)
return false;
if (!reflection_enabled())
return false;
TRACE("non_linear_verbose", tout << "propagating linear monomials...\n";);
bool p = false;
// CMW: m_nl_monomials can grow during this loop, so
// don't use iterators.
for (unsigned i = 0; i < m_nl_monomials.size(); i++) {
if (propagate_linear_monomial(m_nl_monomials[i]))
p = true;
}
CTRACE("non_linear", p, display(tout););
return p;
}
/*
Interval arithmetic does not satisfy distributivity.
Actually, it satisfies the sub-distributivity property:
x*(y + z) \subseteq x*y + x*z
The sub-distributivity property only holds if condensation
is not used. For example:
x * (x^3 + 1) \subseteq x*x^3 + x,
but it is not the case that
x * (x^3 + 1) \subseteq x^4 + x
for example, for x = [-2,1]
x*(x^3+1) = [-7, 14]
x^4 + x = [-2, 17]
This weakness of AI is known as the "dependency problem",
which comes from the decorrelation of the multiple occurrences
of one variable during interval evaluation.
Given a polynomial:
p(x) = a_0 + a_1 * x + ... + a_n * x^n
The horner extension is:
h_p(x) = a_0 + x*(a_1 + ... + x*(a_{n-1} + a_n * x) + ...)
The horner extension of p(x) = x^4 + x^3 + 2*x is:
h_p(x) = x(2 + x^3(1 + x))
The horner extension evaluates tighter intervals when
condensation is not used.
Remark: there is no guarantee that horner extension will
provide a tighter interval than a sum of monomials when
condensation is used.
For multivariate polynomials nested (or cross nested) forms
are used. The idea is to select one variable, and pretend the
other are parameters. The horner form is computed for the selected
variable, and the computation continues for the polynomials on the
parameters.
As described above, the horner form is not optimal with respect to
to condensation. I use the following two properties to deal with
monovariate polynomials with two monomials:
p(x) = a*x^n + b*x^{n+m} for n >= m
is equivalent to
b*x^{n-m}*[(x^{m} + a/(2b))^2 - (a/2b)^2]
This polynomial provides tight bound when n and m have the same parity and:
1) a*b > 0 and (lower(x) >= 0 or upper(x)^m <= -a/b)
2) a*b < 0 and (upper(x) <= 0 or lower(x)^m >= a/b)
This polynomial also provides tight bounds when n = m,
and the polynomial is simplified to, and n and m may have arbitrary parities:
b*[(x^{n} + a/(2b))^2 - (a/2b)^2]
Example:
x > 1
x^2 - x <= 0
is unsatisfiable
If we compute the bounds for x^2 - x we obtain
(-oo, oo).
On the other hand, if we compute the bounds for
(x - 1/2)^2 - 1/4
we obtain the bounds (0, oo), and the inconsistency
is detected.
Remark: In Z3, I condensate multiple occurrences of a variable
when evaluating monomials. So, the interval for a monomial is
always tight.
Remark: M1*(M2 + M3) is more precise than M1 * M2 + M1 * M3,
if intersection(Vars(M1), union(Vars(M2), Vars(M3))) = empty-set,
Remark: A trivial consequence of Moore's theorem for interval
arithmetic. If two monomials M1 and M2 do not share variables,
then the interval for M1 + M2 is tight.
*/
/**
\brief Check whether the same variable occurs in two different monomials.
\remark Fixed variables are ignored.
\remark A trivial consequence of Moore's theorem for interval
arithmetic. If two monomials M1 and M2 do not share variables,
then the interval for M1 + M2 is tight.
*/
template
bool theory_arith::is_problematic_non_linear_row(row const & r) {
m_tmp_var_set.reset();
typename vector::const_iterator it = r.begin_entries();
typename vector::const_iterator end = r.end_entries();
for (; it != end; ++it) {
if (!it->is_dead()) {
theory_var v = it->m_var;
if (is_fixed(v))
continue;
if (is_pure_monomial(v)) {
expr * m = var2expr(v);
for (expr* arg : *to_app(m)) {
theory_var curr = expr2var(arg);
if (m_tmp_var_set.contains(curr))
return true;
}
SASSERT(m == var2expr(v));
for (expr* arg : *to_app(m)) {
theory_var curr = expr2var(arg);
if (!is_fixed(curr))
m_tmp_var_set.insert(curr);
}
}
else {
if (m_tmp_var_set.contains(v))
return true;
SASSERT(!is_fixed(v));
m_tmp_var_set.insert(v);
}
}
}
return false;
}
/**
\brief Return true if the row mixes real and integer variables.
This kind of row cannot be converted back to an expression, since
expressions in Z3 cannot have mixed sorts.
*/
template
bool theory_arith::is_mixed_real_integer(row const & r) const {
bool found_int = false;
bool found_real = false;
typename vector::const_iterator it = r.begin_entries();
typename vector::const_iterator end = r.end_entries();
for (; it != end; ++it) {
if (it->is_dead())
continue;
theory_var v = it->m_var;
// TODO: possible improvement... ignore fixed variables.
// If we implement this improvement, we are actually changing the contract of this function
// and we will also have to fix the affected functions.
if (is_int(v))
found_int = true;
if (is_real(v))
found_real = true;
if (found_int && found_real)
return true;
}
return false;
}
/**
\brief Return true if the row contains only integer variables.
*/
template
bool theory_arith::is_integer(row const & r) const {
typename vector::const_iterator it = r.begin_entries();
typename vector::const_iterator end = r.end_entries();
for (; it != end; ++it) {
if (it->is_dead())
continue;
theory_var v = it->m_var;
// TODO: possible improvement... ignore fixed variables.
if (!is_int(v))
return false;
}
return true;
}
template
void theory_arith::display_coeff_exprs(std::ostream & out, buffer const & p) const {
bool first = true;
for (coeff_expr const& ce : p) {
if (first)
first = false;
else
out << "+\n";
out << ce.first << " * " << mk_pp(ce.second, get_manager()) << "\n";
}
}
/**
\brief Traverse p and store in vars the (non-fixed) variables
that occur in more than one monomial. The number of
occurrences is also stored.
*/
template
bool theory_arith::get_polynomial_info(buffer const & p, sbuffer & varinfo) {
varinfo.reset();
m_var2num_occs.reset();
auto add_occ = [&](expr* v) {
if (has_var(v) && !is_fixed(expr2var(v))) {
TRACE("nl_info", tout << "adding occ: " << mk_bounded_pp(v, get_manager()) << "\n";);
unsigned occs = 0;
m_var2num_occs.find(v, occs);
m_var2num_occs.insert(v, 1 + occs);
}
};
for (auto const& ce : p) {
expr * m = ce.second;
if (m_util.is_numeral(m)) {
continue;
}
else if (false && m_util.is_add(m)) // introduced by #4532, disabled for #4765
return false;
else if (ctx.e_internalized(m) && !is_pure_monomial(m))
add_occ(m);
else {
buffer vp;
decompose_monomial(m, vp);
for (auto const& p : vp) {
add_occ(p.first);
}
}
}
// Update the number of occurrences in the result vector.
for (auto const& vn : m_var2num_occs) {
if (vn.m_value > 1)
varinfo.push_back(var_num_occs(vn.m_key, vn.m_value));
}
return true;
}
/**
\brief Convert p into an expression.
*/
template
expr_ref theory_arith::p2expr(buffer & p) {
SASSERT(!p.empty());
TRACE("p2expr_bug", display_coeff_exprs(tout, p););
ptr_buffer args;
rational c2;
for (coeff_expr const& ce : p) {
rational const & c = ce.first;
expr * var = ce.second;
if (m_util.is_numeral(var, c2)) {
expr* m = m_util.mk_numeral(c * c2, c.is_int() && m_util.is_int(var));
m_nl_new_exprs.push_back(m);
args.push_back(m);
}
else if (!c.is_one()) {
expr * m = m_util.mk_mul(m_util.mk_numeral(c, c.is_int() && m_util.is_int(var)), var);
m_nl_new_exprs.push_back(m);
args.push_back(m);
}
else {
args.push_back(var);
}
}
SASSERT(!args.empty());
expr_ref r(mk_nary_add(args.size(), args.data()), get_manager());
m_nl_new_exprs.push_back(r);
return r;
}
/**
\brief Return expression representing: var^power
*/
template
expr * theory_arith::power(expr * var, unsigned power) {
SASSERT(power > 0);
expr * r = var;
for (unsigned i = 1; i < power; i++)
r = m_util.mk_mul(var, r);
m_nl_new_exprs.push_back(r);
return r;
}
/**
\brief Return true if var only occurs in two monovariate monomials,
and return its power and coefficients and these monomials.
The arguments i1 and i2 contain the position in p of the two monomials.
*/
template
bool theory_arith::in_monovariate_monomials(buffer & p, expr * var,
unsigned & i1, rational & c1, unsigned & n1, unsigned & i2, rational & c2, unsigned & n2) {
int idx = 0;
auto set_result = [&](unsigned i, unsigned power, coeff_expr const& ce) {
if (idx == 0) {
c1 = ce.first;
n1 = power;
idx = 1;
i1 = i;
}
else if (idx == 1) {
c2 = ce.first;
n2 = power;
idx = 2;
i2 = i;
}
else {
idx = 3;
}
};
for (unsigned i = 0; i < p.size() && idx != 3; ++i) {
auto const& ce = p[i];
expr * m = ce.second;
if (is_pure_monomial(m)) {
buffer vp;
decompose_monomial(m, vp);
for (auto const& p : vp) {
if (p.first == var) {
if (vp.size() > 1)
return false;
set_result(i, p.second, ce); }
}
}
else if (m == var) {
set_result(i, 1, ce);
}
}
if (idx != 2)
return false;
return true;
}
template
bool theory_arith::is_pure_monomial(expr* mon) const {
if (!m_util.is_mul(mon))
return false;
app* p = to_app(mon);
for (expr* arg : *p)
if (m_util.is_numeral(arg) || m_util.is_mul(arg))
return false;
return true;
}
/**
\brief Display a nested form expression
*/
template
void theory_arith::display_nested_form(std::ostream & out, expr * p) {
if (has_var(p)) {
out << "#" << p->get_id();
}
else if (m_util.is_add(p)) {
SASSERT(!has_var(p));
out << "(";
for (unsigned i = 0; i < to_app(p)->get_num_args(); i++) {
if (i > 0) out << " + ";
display_nested_form(out, to_app(p)->get_arg(i));
}
out << ")";
}
else if (m_util.is_mul(p)) {
buffer vp;
rational c = decompose_monomial(p, vp);
bool first = true;
if (!c.is_one()) {
out << c;
first = false;
}
for (auto const& pair : vp) {
if (first) first = false; else out << "*";
expr * var = pair.first;
unsigned power = pair.second;
display_nested_form(out, var);
if (power != 1)
out << "^" << power;
}
}
else {
rational val;
if (m_util.is_numeral(p, val))
out << val;
else
out << "[unknown #" << p->get_id() << "]";
}
}
/**
\brief Return the degree of var in m.
*/
template
unsigned theory_arith::get_degree_of(expr * m, expr * var) {
if (m == var)
return 1;
if (is_pure_monomial(m)) {
buffer vp;
decompose_monomial(m, vp);
for (auto const& p : vp) {
if (p.first == var)
return p.second;
}
}
return 0;
}
/**
\brief Return the minimal degree of var in the polynomial p.
*/
template
unsigned theory_arith::get_min_degree(buffer & p, expr * var) {
SASSERT(!p.empty());
SASSERT(var != 0);
// get monomial where the degree of var is min.
unsigned d = UINT_MAX; // min. degree of var
for (coeff_expr const& ce : p) {
expr * m = ce.second;
d = std::min(d, get_degree_of(m, var));
if (d == 0)
return d;
}
SASSERT(d != UINT_MAX);
return d;
}
/**
\brief Divide m by var^d.
*/
template
expr * theory_arith::factor(expr * m, expr * var, unsigned d) {
TRACE("factor", tout << "m: " << mk_pp(m, get_manager()) << "\nvar: " << mk_pp(var, get_manager()) << "\nd: " << d << "\n";);
if (d == 0)
return m;
if (m == var) {
SASSERT(d == 1);
expr * result = m_util.mk_numeral(rational(1), m_util.is_int(var));
m_nl_new_exprs.push_back(result);
return result;
}
ptr_buffer new_args;
unsigned idx = 0;
auto insert = [&](expr* arg) {
if (idx < d && var == arg)
++idx;
else
new_args.push_back(arg);
};
while (m_util.is_mul(m) && idx < d) {
unsigned sz = to_app(m)->get_num_args();
for (unsigned i = 0; i + 1 < sz; ++i) {
insert(to_app(m)->get_arg(i));
}
m = to_app(m)->get_arg(sz - 1);
}
insert(m);
SASSERT(idx == d);
TRACE("factor_bug", tout << "new_args:\n"; for(unsigned i = 0; i < new_args.size(); i++) tout << mk_pp(new_args[i], get_manager()) << "\n";);
expr * result = mk_nary_mul(new_args.size(), new_args.data(), m_util.is_int(var));
m_nl_new_exprs.push_back(result);
TRACE("factor", tout << "result: " << mk_pp(result, get_manager()) << "\n";);
return result;
}
/**
\brief Return the horner extension of p with respect to var.
*/
template
expr_ref theory_arith::horner(unsigned depth, buffer & p, expr * var) {
SASSERT(!p.empty());
SASSERT(var != 0);
unsigned d = get_min_degree(p, var);
TRACE("horner_bug", tout << "poly:\n";
for (unsigned i = 0; i < p.size(); i++) { if (i > 0) tout << " + "; tout << p[i].first << "*" << mk_pp(p[i].second, get_manager()); } tout << "\n";
tout << "var: " << mk_pp(var, get_manager()) << "\n";
tout << "min_degree: " << d << "\n";);
buffer e; // monomials/x^d where var occurs with degree d
buffer r; // rest
for (auto const& kv : p) {
expr * m = kv.second;
expr * f = factor(m, var, d);
if (get_degree_of(m, var) == d) {
e.push_back(coeff_expr(kv.first, f));
}
else {
SASSERT(get_degree_of(m, var) > d);
r.push_back(coeff_expr(kv.first, f));
}
}
expr_ref s = cross_nested(depth + 1, e, nullptr);
if (!r.empty()) {
expr_ref q = horner(depth + 1, r, var);
// TODO: improve here
s = m_util.mk_add(q, s);
}
expr_ref result = s;
if (d != 0) {
expr * xd = power(var, d);
result = m_util.mk_mul(xd, s);
}
m_nl_new_exprs.push_back(result);
return result;
}
/**
\brief Convert the polynomial p into an equivalent cross nested
expression. The idea is to obtain an expression e where
evaluate_as_interval(e) is more precise than
evaluate_as_interval(p).
If var != 0, then it is used for performing the horner extension
*/
template
expr_ref theory_arith::cross_nested(unsigned depth, buffer & p, expr * var) {
TRACE("non_linear", tout << "p.size: " << p.size() << "\n";);
if (var == nullptr) {
sbuffer varinfo;
if (!get_polynomial_info(p, varinfo))
return p2expr(p);
if (varinfo.empty())
return p2expr(p);
unsigned max = 0;
for (auto const& kv : varinfo) {
if (kv.second >= max) {
max = kv.second;
var = kv.first;
}
}
}
if (depth > 20)
return p2expr(p);
SASSERT(var != nullptr);
unsigned i1 = UINT_MAX;
unsigned i2 = UINT_MAX;
rational a, b;
unsigned n = UINT_MAX;
unsigned nm = UINT_MAX;
if (in_monovariate_monomials(p, var, i1, a, n, i2, b, nm)) {
CTRACE("in_monovariate_monomials", n == nm,
for (unsigned i = 0; i < p.size(); i++) {
if (i > 0) tout << " + "; tout << p[i].first << "*" << mk_pp(p[i].second, get_manager());
}
tout << "\n";
tout << "var: " << mk_pp(var, get_manager()) << "\n";
tout << "i1: " << i1 << "\n";
tout << "a: " << a << "\n";
tout << "n: " << n << "\n";
tout << "i2: " << i2 << "\n";
tout << "b: " << b << "\n";
tout << "nm: " << nm << "\n";);
if (n == nm) return horner(depth, p, var);
SASSERT(n != nm);
expr * new_expr = nullptr;
if (nm < n) {
std::swap(n, nm);
std::swap(a, b);
}
SASSERT(nm > n);
unsigned m = nm - n;
if (n % 2 == m % 2 && n >= m) {
// b*x^{n-m}*[(x^{m} + a/(2b))^2 - (a/2b)^2]
// b*[(x^{m} + a/(2b))^2 - (a/2b)^2] for n == m
rational a2b = a;
expr_ref xm(power(var, m), get_manager());
a2b /= (rational(2) * b);
// we cannot create a numeral that has sort int, but it is a rational.
if (!m_util.is_int(var) || a2b.is_int()) {
rational ma2b2 = a2b * a2b;
ma2b2.neg();
expr * xm_a2b = m_util.mk_add(m_util.mk_numeral(a2b, m_util.is_int(var)), xm);
expr * xm_a2b2 = m_util.mk_mul(xm_a2b, xm_a2b);
expr * rhs = m_util.mk_add(xm_a2b2, m_util.mk_numeral(ma2b2, m_util.is_int(var)));
expr * rhs2 = nullptr;
if (n > m)
rhs2 = m_util.mk_mul(power(var, n - m), rhs);
else
rhs2 = rhs;
new_expr = b.is_one() ? rhs2 : m_util.mk_mul(m_util.mk_numeral(b, m_util.is_int(var)), rhs2);
m_nl_new_exprs.push_back(new_expr);
TRACE("non_linear", tout << "new_expr:\n"; display_nested_form(tout, new_expr); tout << "\n";);
buffer rest;
unsigned sz = p.size();
for (unsigned i = 0; i < sz; i++) {
if (i != i1 && i != i2)
rest.push_back(p[i]);
}
if (rest.empty())
return expr_ref(new_expr, get_manager());
TRACE("non_linear", tout << "rest size: " << rest.size() << ", i1: " << i1 << ", i2: " << i2 << "\n";);
expr_ref h = cross_nested(depth + 1, rest, nullptr);
expr * r = m_util.mk_add(new_expr, h);
m_nl_new_exprs.push_back(r);
return expr_ref(r, get_manager());
}
}
}
return horner(depth, p, var);
}
/**
\brief Check whether the given polynomial is consistent with respect to the known
bounds. The polynomial is converted into an equivalent cross nested form.
*/
template
bool theory_arith::is_cross_nested_consistent(buffer & p) {
sbuffer varinfo;
if (!get_polynomial_info(p, varinfo))
return true;
if (varinfo.empty())
return true;
std::stable_sort(varinfo.begin(), varinfo.end(), var_num_occs_lt());
TRACE("cross_nested", tout << "var num occs:\n";
for (auto const& kv : varinfo) {
tout << mk_bounded_pp(kv.first, get_manager()) << " -> " << kv.second << "\n";
});
for (auto const& kv : varinfo) {
m_nl_new_exprs.reset();
expr * var = kv.first;
expr_ref cn = cross_nested(0, p, var);
// Remark: cn may not be well-sorted because, since a row may contain mixed integer/real monomials.
// This is not really a problem, since evaluate_as_interval will work even if cn is not well-sorted.
if (!cn)
continue;
TRACE("cross_nested", tout << "nested form for var:\n" << mk_ismt2_pp(var, get_manager()) << "\n";
display_nested_form(tout, cn); tout << "\n";
tout << "c:\n" << mk_ismt2_pp(cn, get_manager()) << "\n";);
interval i = evaluate_as_interval(cn);
TRACE("cross_nested", tout << "interval: " << i << "\n";);
v_dependency * d = nullptr;
if (!i.minus_infinity() && (i.get_lower_value().is_pos() || (i.get_lower_value().is_zero() && i.is_lower_open())))
d = i.get_lower_dependencies();
else if (!i.plus_infinity() && (i.get_upper_value().is_neg() || (i.get_upper_value().is_zero() && i.is_upper_open())))
d = i.get_upper_dependencies();
if (d) {
TRACE("cross_nested", tout << "nested form conflict: " << i << "\n";);
set_conflict(d);
return false;
}
}
return true;
}
/**
\brief Check whether the polynomial represented by the current row is
consistent with respect to the known bound when converted into a
equivalent cross nested form.
*/
template
bool theory_arith::is_cross_nested_consistent(row const & r) {
if (!is_problematic_non_linear_row(r))
return true;
TRACE("cross_nested", tout << "is_cross_nested_consistent:\n"; display_row(tout, r, false);
display(tout);
);
/*
The method is_cross_nested converts rows back to expressions.
The conversion back to expressions may create sort incorrect expressions.
This is in some sense ok, since these expressions are temporary, but
the sort incorrect expressions may generate assertion violations.
Sort incorrect expressions may be created in the following cases:
1) mixed real int rows.
2) int rows that contain non integer coefficients.
3) int rows that when converted to cross nested form use non integer coefficients.
There are several ways to deal with this problem:
a) Ignore the assertion violations. Disadvantage: it will prevent us from running Z3 in debug mode on some benchmarks.
b) Remove the assertions. Disadvantage: these assertions helped us to find many important bugs in Z3
c) Disable the assertions temporally. This sounds like a big HACK.
d) Use a different data-structure to represent polynomials in cross-nested form. Disadvantage: code duplication, the data-structure
is essentially identical to the ASTs we are using right now.
e) Disable the test when we cannot create a well-sorted expression.
I'm temporally using this solution.
I implemented the following logic:
1) (mixed real int) Disable the test. Most benchmarks do not contain mixed int real variables.
2) (int coeffs) I multiply the row by a constant to force it to have only integer coefficients.
3) (new non-int coeffs) This only happens in an optional step in the conversion. Now, for int rows, I only apply this optional step only if non-int coeffs are not created.
*/
if (!get_manager().int_real_coercions() && is_mixed_real_integer(r))
return true; // giving up... see comment above
TRACE("cross_nested", tout << "checking problematic row...\n";);
rational c = rational::one();
if (is_integer(r))
c = r.get_denominators_lcm().to_rational();
TRACE("non_linear", tout << "check problematic row:\n"; display_row(tout, r); display_row(tout, r, false););
buffer p;
for (auto & col : r) {
if (!col.is_dead()) {
p.push_back(coeff_expr(col.m_coeff.to_rational() * c, var2expr(col.m_var)));
}
}
SASSERT(!p.empty());
CTRACE("cross_nested_bug", !c.is_one(), tout << "c: " << c << "\n"; display_row(tout, r); tout << "---> p (coeffs, exprs):\n"; display_coeff_exprs(tout, p););
return is_cross_nested_consistent(p);
}
/**
\brief Check whether an inconsistency can be found using cross nested
form in the non linear cluster.
*/
template
bool theory_arith::is_cross_nested_consistent(svector const & nl_cluster) {
for (theory_var v : nl_cluster) {
if (!is_base(v))
continue;
m_stats.m_nl_cross_nested++;
row const & r = m_rows[get_var_row(v)];
if (!is_cross_nested_consistent(r))
return false;
}
return true;
}
#define FIXED 0
#define QUOTED_FIXED 1
#define BOUNDED 2
#define QUOTED_BOUNDED 3
#define NOT_FREE 4
#define QUOTED_NOT_FREE 5
#define FREE 6
#define QUOTED_FREE 7
#define MAX_DEFAULT_WEIGHT 7
/**
\brief Initialize variable order for grobner basis computation.
Make:
"quoted free vars" > "free vars" > "quoted variables with lower or upper bounds" >
"variables with lower or upper bounds" > "quoted bounded variables" >
"bounded variables" > "quoted fixed variables" > "fixed variables"
*/
template
void theory_arith::init_grobner_var_order(svector const & nl_cluster, grobner & gb) {
// Initialize variable order
for (theory_var v : nl_cluster) {
expr * var = var2expr(v);
if (is_fixed(v)) {
gb.set_weight(var, is_pure_monomial(var) ? QUOTED_FIXED : FIXED);
}
else if (is_bounded(v)) {
gb.set_weight(var, is_pure_monomial(var) ? QUOTED_BOUNDED : BOUNDED);
}
else if (lower(v) || upper(v)) {
gb.set_weight(var, is_pure_monomial(var) ? QUOTED_NOT_FREE : NOT_FREE);
}
else {
SASSERT(is_free(v));
gb.set_weight(var, is_pure_monomial(var) ? QUOTED_FREE : FREE);
}
}
}
/**
\brief Create a new monomial using the given coeff and m. Fixed
variables in m are substituted by their values. The arg dep is
updated to store these dependencies. The set already_found is
updated with the fixed variables in m. A variable is only
added to dep if it is not already in already_found.
Return null if the monomial was simplified to 0.
*/
template
grobner::monomial * theory_arith::mk_gb_monomial(rational const & _coeff, expr * m, grobner & gb, v_dependency * & dep, var_set & already_found) {
ptr_buffer vars;
rational coeff = _coeff;
rational r;
auto proc_var = [&](expr* v) {
if (m_util.is_numeral(v, r)) {
coeff *= r;
return;
}
theory_var _var = expr2var(v);
if (is_fixed(_var)) {
if (!already_found.contains(_var)) {
already_found.insert(_var);
dep = m_dep_manager.mk_join(dep, m_dep_manager.mk_join(m_dep_manager.mk_leaf(lower(_var)), m_dep_manager.mk_leaf(upper(_var))));
}
coeff *= lower_bound(_var).get_rational().to_rational();
}
else {
vars.push_back(v);
}
};
while (m_util.is_mul(m)) {
unsigned sz = to_app(m)->get_num_args();
for (unsigned i = 0; i + 1 < sz; ++i) {
proc_var(to_app(m)->get_arg(i));
}
m = to_app(m)->get_arg(sz-1);
}
proc_var(m);
if (!coeff.is_zero())
return gb.mk_monomial(coeff, vars.size(), vars.data());
else
return nullptr;
}
/**
\brief Send the given row to the grobner basis object.
All fixed variables are substituted before sending the row to gb.
*/
template
void theory_arith::add_row_to_gb(row const & r, grobner & gb) {
TRACE("grobner", tout << "adding row to gb\n"; display_row(tout, r););
ptr_buffer monomials;
v_dependency * dep = nullptr;
m_tmp_var_set.reset();
for (auto& re : r) {
if (re.is_dead())
continue;
rational coeff = re.m_coeff.to_rational();
expr * m = var2expr(re.m_var);
grobner::monomial * new_m = mk_gb_monomial(coeff, m, gb, dep, m_tmp_var_set);
if (new_m)
monomials.push_back(new_m);
TRACE("grobner",
tout << "monomial: " << mk_pp(m, get_manager()) << "\n";
tout << "new monomial: "; if (new_m) gb.display_monomial(tout, *new_m); else tout << "null"; tout << "\n";);
}
gb.assert_eq_0(monomials.size(), monomials.data(), dep);
}
/**
\brief v must be a pure monomial. That is, v = (quote (* x_1 ... x_n))
Add the monomial (quote (* x_1 ... x_n)) = x_1 * ... * x_n.
Fixed variables are substituted.
*/
template
void theory_arith::add_monomial_def_to_gb(theory_var v, grobner & gb) {
ptr_buffer monomials;
v_dependency * dep = nullptr;
m_tmp_var_set.reset();
expr * m = var2expr(v);
SASSERT(is_pure_monomial(m));
grobner::monomial * new_m = mk_gb_monomial(rational(1), m, gb, dep, m_tmp_var_set);
if (new_m)
monomials.push_back(new_m);
rational coeff(-1);
if (is_fixed(v)) {
dep = m_dep_manager.mk_join(dep, m_dep_manager.mk_join(m_dep_manager.mk_leaf(lower(v)), m_dep_manager.mk_leaf(upper(v))));
coeff *= lower_bound(v).get_rational().to_rational();
if (!coeff.is_zero())
monomials.push_back(gb.mk_monomial(coeff, 0, nullptr));
}
else {
monomials.push_back(gb.mk_monomial(coeff, 1, &m));
}
gb.assert_eq_0(monomials.size(), monomials.data(), dep);
}
/**
Initialize grobner basis data structure using the non linear cluster.
The GB is initialized using rows and non linear monomials.
*/
template
void theory_arith::init_grobner(svector const & nl_cluster, grobner & gb) {
init_grobner_var_order(nl_cluster, gb);
for (theory_var v : nl_cluster) {
if (is_base(v)) {
row const & r = m_rows[get_var_row(v)];
add_row_to_gb(r, gb);
}
if (is_pure_monomial(v) && !m_data[v].m_nl_propagated && is_fixed(v)) {
add_monomial_def_to_gb(v, gb);
}
}
}
/**
\brief Return the interval for the given monomial
*/
template
interval theory_arith::mk_interval_for(grobner::monomial const * m) {
interval r(m_dep_manager, rational(m->get_coeff()));
expr * var = nullptr;
unsigned power = 0;
unsigned num_vars = m->get_degree();
for (unsigned i = 0; i < num_vars; i++) {
expr * curr = m->get_var(i);
if (var == nullptr) {
var = curr;
power = 1;
}
else if (curr == var) {
power++;
}
else {
mul_bound_of(var, power, r);
var = curr;
power = 1;
}
}
if (var != nullptr)
mul_bound_of(var, power, r);
return r;
}
/**
\brief Set a conflict using a dependency object.
*/
template
void theory_arith::set_conflict(v_dependency * d) {
antecedents ante(*this);
derived_bound b(null_theory_var, inf_numeral(0), B_LOWER);
dependency2new_bound(d, b);
set_conflict(b, ante, "arith_nl");
TRACE("non_linear", for (literal lit : b.m_lits) ctx.display_literal_verbose(tout, lit) << "\n"; tout << "\n";);
}
/**
\brief Return true if I.get_lower() <= - M_1 - ... - M_n <= I.get_upper() is inconsistent.
Where M_i is monomials[i] and n = num_monomials.
A conflict will also be set using the bounds of the variables occurring in the monomials M_i's.
*/
template
bool theory_arith::is_inconsistent(interval const & I, unsigned num_monomials, grobner::monomial * const * monomials, v_dependency * dep) {
interval r(I);
for (unsigned i = 0; i < num_monomials; i++) {
grobner::monomial const * m = monomials[i];
r += mk_interval_for(m);
if (r.minus_infinity() && r.plus_infinity())
return false;
}
TRACE("non_linear_bug", tout << "is_inconsistent, r: " << r << "\n";);
v_dependency * interval_deps = nullptr;
bool conflict = false;
if (!r.minus_infinity() && (r.get_lower_value().is_pos() || (r.get_lower_value().is_zero() && r.is_lower_open()))) {
interval_deps = r.get_lower_dependencies();
conflict = true;
TRACE("non_linear_bug", tout << "is inconsistent, interval_deps: " << interval_deps << "\n";);
}
else if (!r.plus_infinity() && (r.get_upper_value().is_neg() || (r.get_upper_value().is_zero() && r.is_upper_open()))) {
interval_deps = r.get_upper_dependencies();
conflict = true;
TRACE("non_linear_bug", tout << "is inconsistent, interval_deps: " << interval_deps << "\n";);
}
// interval_deps cannot be used to check if a conflict was detected, since interval_deps may be 0 even when r does not contain 0
if (conflict) {
TRACE("non_linear", tout << "conflicting interval for = 0 equation: " << r << "\n";);
set_conflict(m_dep_manager.mk_join(interval_deps, dep));
return true;
}
return false;
}
/**
\brief Return true if the equation is inconsistent,
and sign a conflict.
*/
template
bool theory_arith::is_inconsistent(grobner::equation const * eq, grobner & gb) {
interval zero(m_dep_manager, rational(0));
if (is_inconsistent(zero, eq->get_num_monomials(), eq->get_monomials(), eq->get_dependency())) {
TRACE("non_linear", tout << "found conflict\n"; gb.display_equation(tout, *eq););
return true;
}
return false;
}
/**
\brief Return true if the given monomial c*M is squared.
The square root of the c is stored in r.
*/
bool is_perfect_square(grobner::monomial const * m, rational & r) {
unsigned num_vars = m->get_degree();
if (num_vars % 2 == 1)
return false;
if (!m->get_coeff().is_perfect_square(r))
return false;
expr * var = nullptr;
unsigned power = 0;
for (unsigned i = 0; i < num_vars; i++) {
expr * curr = m->get_var(i);
if (var == nullptr) {
var = curr;
power = 1;
}
else if (var == curr) {
power++;
}
else {
if (power % 2 == 1)
return false;
var = curr;
power = 1;
}
}
return power % 2 == 0;
}
/**
\brief Return m1m2 is of the form (-2ab)*M1*M2
assuming that
m1_sq = a^2*M1*M1
m2_sq = b^2*M2*M2
*/
bool is_perfect_square(grobner::monomial const * m1_sq, rational const & a, grobner::monomial const * m2_sq, rational const & b, grobner::monomial const * m1m2) {
DEBUG_CODE({
rational a1;
rational b1;
SASSERT(is_perfect_square(m1_sq, a1) && a == a1 && is_perfect_square(m2_sq, b1) && b == b1);
});
if (m1m2->get_coeff().is_nonneg())
return false;
rational c(-2);
c *= a;
c *= b;
if (m1m2->get_coeff() != c)
return false;
unsigned num1 = m1_sq->get_degree();
unsigned num2 = m2_sq->get_degree();
unsigned num12 = m1m2->get_degree();
if (num1 + num2 != num12 * 2)
return false;
unsigned i1, i2, i12;
i1 = i2 = i12 = 0;
while (true) {
expr * v1 = nullptr;
expr * v2 = nullptr;
expr * v12 = nullptr;
if (i1 < num1)
v1 = m1_sq->get_var(i1);
if (i2 < num2)
v2 = m2_sq->get_var(i2);
if (i12 < num12)
v12 = m1m2->get_var(i12);
if (v1 == nullptr && v2 == nullptr && v12 == nullptr)
return true;
if (v12 == nullptr)
return false;
if (v1 == v12) {
SASSERT(m1_sq->get_var(i1+1) == v1);
i1 += 2;
i12 ++;
}
else if (v2 == v12) {
SASSERT(m2_sq->get_var(i2+1) == v2);
i2 += 2;
i12 ++;
}
else {
return false;
}
}
}
/**
\brief Return true if the equation is inconsistent. In this
version, perfect squares are eliminated, and replaced with the
interval [0, oo), if the interval associated with them is less
precise than [0, oo).
\remark I track only simple perfect squares of the form (M1 - M2)^2,
where M1 and M2 are arbitrary monomials.
*/
template
bool theory_arith::is_inconsistent2(grobner::equation const * eq, grobner & gb) {
// TODO: a possible improvement: create a quotation for (M1 - M2)^2
// instead of trying to find it in a specific equation.
// This approach is more precise, but more expensive
// since a new row must be created.
buffer intervals;
unsigned num = eq->get_num_monomials();
for (unsigned i = 0; i < num; i++) {
grobner::monomial const * m = eq->get_monomial(i);
intervals.push_back(mk_interval_for(m));
}
sbuffer deleted;
deleted.resize(num, false);
ptr_buffer monomials;
// try to eliminate monomials that form perfect squares of the form (M1 - M2)^2
for (unsigned i = 0; i < num; i++) {
grobner::monomial const * m1 = eq->get_monomial(i);
rational a;
if (deleted[i])
continue;
if (!is_perfect_square(m1, a)) {
monomials.push_back(const_cast(m1));
continue;
}
TRACE("non_linear", tout << "found perfect square monomial m1: "; gb.display_monomial(tout, *m1); tout << "\n";);
// try to find another perfect square
unsigned j = i + 1;
for (; j < num; j++) {
if (deleted[j])
continue;
grobner::monomial const * m2 = eq->get_monomial(j);
rational b;
if (!is_perfect_square(m2, b))
continue;
TRACE("non_linear", tout << "found perfect square monomial m2: "; gb.display_monomial(tout, *m2); tout << "\n";);
// try to find -2*root(m1)*root(m2)
// This monomial must be smaller than m1, since m2 is smaller than m1.
unsigned k = i + 1;
for (; k < num; k++) {
if (deleted[k])
continue;
grobner::monomial const * m1m2 = eq->get_monomial(k);
if (!is_perfect_square(m1, a, m2, b, m1m2))
continue;
// m1, m2, and m1m2 form a perfect square.
// check if [0, oo) provides a better lowerbound than adding the intervals of m1, m2 and m1m2;
TRACE("non_linear", tout << "found perfect square (M1-M2)^2:\n";
gb.display_monomial(tout, *m1); tout << "\n";
gb.display_monomial(tout, *m2); tout << "\n";
gb.display_monomial(tout, *m1m2); tout << "\n";);
interval I = intervals[i];
I += intervals[j];
I += intervals[k];
if (I.minus_infinity() || I.get_lower_value().is_neg()) {
TRACE("non_linear", tout << "the lower bound improved when perfect square is eliminated.\n";);
// Found improvement...
// mark these monomials as deleted
deleted[i] = true;
deleted[j] = true;
deleted[k] = true;
break;
}
}
if (k < num)
break; // found perfect square
}
if (j == num) {
// didn't find perfect square of the form (M1-M2)^2
monomials.push_back(const_cast(m1));
}
}
if (monomials.size() == num)
return false; // didn't find any perfect square.
interval ge_zero(m_dep_manager, rational(0), false, true, nullptr);
if (is_inconsistent(ge_zero, monomials.size(), monomials.data(), eq->get_dependency())) {
TRACE("non_linear", tout << "found conflict\n"; gb.display_equation(tout, *eq););
return true;
}
return false;
}
template
expr * theory_arith::monomial2expr(grobner::monomial const * m, bool is_int) {
unsigned num_vars = m->get_degree();
ptr_buffer args;
if (!m->get_coeff().is_one())
args.push_back(m_util.mk_numeral(m->get_coeff(), is_int));
for (unsigned j = 0; j < num_vars; j++)
args.push_back(m->get_var(j));
return mk_nary_mul(args.size(), args.data(), is_int);
}
/**
\brief Assert the new equation in the simplex tableau.
*/
template
bool theory_arith::internalize_gb_eq(grobner::equation const * eq) {
bool is_int = false;
unsigned num_monomials = eq->get_num_monomials();
for (unsigned i = 0; i < num_monomials; i++) {
grobner::monomial const * m = eq->get_monomial(i);
unsigned degree = m->get_degree();
if (degree > m_params.m_nl_arith_max_degree)
return false;
if (degree > 0)
is_int = m_util.is_int(m->get_var(0));
}
rational k;
ptr_buffer args;
for (unsigned i = 0; i < num_monomials; i++) {
grobner::monomial const * m = eq->get_monomial(i);
if (m->get_degree() == 0)
k -= m->get_coeff();
else
args.push_back(monomial2expr(eq->get_monomial(i), is_int));
}
th_rewriter& s = ctx.get_rewriter();
expr_ref pol(get_manager());
SASSERT(!args.empty());
pol = mk_nary_add(args.size(), args.data());
expr_ref s_pol(get_manager());
proof_ref pr(get_manager());
TRACE("gb_bug", tout << mk_ll_pp(pol, get_manager()) << "\n";);
s(pol, s_pol, pr);
if (!has_var(s_pol)) {
TRACE("spol_bug", tout << "internalizing...\n" << mk_ll_pp(s_pol, get_manager()) << "\n";);
ctx.internalize(s_pol, false);
ctx.mark_as_relevant(s_pol.get());
}
SASSERT(has_var(s_pol.get()));
// s_pol = k
theory_var v = expr2var(s_pol);
// v = k
CTRACE("spol_bug", v == null_theory_var, tout << mk_ll_pp(s_pol, get_manager()) << "\n"; display(tout););
SASSERT(v != null_theory_var);
// assert bounds for s_pol
mk_derived_nl_bound(v, inf_numeral(k), B_LOWER, eq->get_dependency());
mk_derived_nl_bound(v, inf_numeral(k), B_UPPER, eq->get_dependency());
TRACE("non_linear", tout << "inserted new equation into the tableau\n"; display_var(tout, v););
return true;
}
template
void theory_arith::update_statistics(grobner & gb) {
m_stats.m_gb_simplify += gb.m_stats.m_simplify;
m_stats.m_gb_superpose += gb.m_stats.m_superpose;
m_stats.m_gb_num_processed += gb.m_stats.m_num_processed;
m_stats.m_gb_compute_basis++;
}
template
void theory_arith::set_gb_exhausted() {
IF_VERBOSE(3, verbose_stream() << "Grobner basis computation interrupted. Increase threshold using NL_ARITH_GB_THRESHOLD=\n";);
ctx.push_trail(value_trail(m_nl_gb_exhausted));
m_nl_gb_exhausted = true;
}
// Scan the grobner basis eqs, and look for inconsistencies.
template
bool theory_arith::get_gb_eqs_and_look_for_conflict(ptr_vector& eqs, grobner& gb) {
eqs.reset();
gb.get_equations(eqs);
TRACE("grobner", tout << "after gb\n";
std::function _fn = [&](std::ostream& out, expr* v) { out << "v" << expr2var(v); };
for (grobner::equation* eq : eqs)
gb.display_equation(tout, *eq, _fn);
);
for (grobner::equation* eq : eqs) {
if (is_inconsistent(eq, gb) || is_inconsistent2(eq, gb)) {
TRACE("grobner", tout << "inconsistent: "; gb.display_equation(tout, *eq););
return true;
}
}
TRACE("grobner", tout << "not found\n"; );
return false;
}
// Try to change the variable order... in such a way the leading term is modified.
// I only consider linear equations... (HACK)
// Moreover, I do not change the weight of a variable more than once in this loop.
template
bool theory_arith::try_to_modify_eqs(ptr_vector& eqs, grobner& gb, unsigned & next_weight) {
bool modified = false;
for (grobner::equation const* eq : eqs) {
unsigned num_monomials = eq->get_num_monomials();
CTRACE("grobner_bug", num_monomials <= 0, gb.display_equation(tout, *eq););
if (num_monomials == 0)
continue; // HACK: the equation 0 = 0, should have been discarded by the GB module.
if (eq->get_monomial(0)->get_degree() != 1)
continue;
for (unsigned j = 1; j < num_monomials; j++) {
grobner::monomial const * m = eq->get_monomial(j);
if (m->get_degree() != 1)
continue;
expr * var = m->get_var(0);
if (gb.get_weight(var) > MAX_DEFAULT_WEIGHT)
continue; // variable was already updated
TRACE("non_linear", tout << "increased weight of: " << mk_bounded_pp(var, get_manager()) << "\n";);
gb.set_weight(var, next_weight);
next_weight++;
gb.update_order();
TRACE("non_linear", tout << "after updating order\n"; gb.display(tout););
modified = true;
break;
}
if (modified)
break;
}
return modified;
}
// Scan the grobner basis eqs for equations of the form x - k = 0 or x = 0 is found, and x is not fixed,
// then assert bounds for x, and continue
template
bool theory_arith::scan_for_linear(ptr_vector& eqs, grobner& gb) {
bool result = false;
if (m_params.m_nl_arith_gb_eqs) { // m_nl_arith_gb_eqs is false by default
SASSERT(false);
for (grobner::equation* eq : eqs) {
if (!eq->is_linear_combination()) {
TRACE("non_linear", tout << "processing new equality:\n"; gb.display_equation(tout, *eq););
TRACE("non_linear_bug", tout << "processing new equality:\n"; gb.display_equation(tout, *eq););
if (internalize_gb_eq(eq))
result = true;
}
}
}
return result;
}
template
bool theory_arith::compute_basis_loop(grobner & gb) {
while (gb.get_num_new_equations() < m_params.m_nl_arith_gb_threshold) {
if (ctx.get_cancel_flag()) {
return false;
}
if (gb.compute_basis_step())
return true;
}
return false;
}
template
void theory_arith::compute_basis(grobner& gb, bool& warn) {
gb.compute_basis_init();
if (!compute_basis_loop(gb) && !warn) {
set_gb_exhausted();
warn = true;
}
}
/**
\brief Compute Grobner basis, return true if a conflict or new fixed variables were detected.
*/
template
typename theory_arith::gb_result theory_arith::compute_grobner(svector const & nl_cluster) {
if (m_nl_gb_exhausted)
return GB_FAIL;
grobner gb(get_manager(), m_dep_manager);
init_grobner(nl_cluster, gb);
TRACE("non_linear", display(tout););
bool warn = false;
unsigned next_weight = MAX_DEFAULT_WEIGHT + 1; // next weight using during perturbation phase.
ptr_vector eqs;
do {
TRACE("grobner", tout << "before grobner:\n";
std::function _fn = [&](std::ostream& out, expr* v) { out << "v" << expr2var(v); };
gb.display(tout, _fn));
compute_basis(gb, warn);
update_statistics(gb);
TRACE("non_linear_gb", tout << "after:\n"; gb.display(tout););
if (ctx.get_cancel_flag())
return GB_FAIL;
if (get_gb_eqs_and_look_for_conflict(eqs, gb))
return GB_PROGRESS;
if (scan_for_linear(eqs, gb))
return GB_NEW_EQ;
}
while(m_params.m_nl_arith_gb_perturbate &&
(!m_nl_gb_exhausted) && try_to_modify_eqs(eqs, gb, next_weight));
return GB_FAIL;
}
/**
\brief Maximize/Minimize variables in non linear monomials.
*/
template
bool theory_arith::max_min_nl_vars() {
if (!m_params.m_nl_arith_optimize_bounds)
return true;
var_set already_found;
svector vars;
for (theory_var v : m_nl_monomials) {
mark_var(v, vars, already_found);
expr * n = var2expr(v);
SASSERT(is_pure_monomial(n));
for (expr * curr : *to_app(n)) {
if (ctx.e_internalized(curr)) {
theory_var v = expr2var(curr);
SASSERT(v != null_theory_var);
mark_var(v, vars, already_found);
}
}
}
return max_min(vars);
}
/**
\brief Process non linear constraints.
*/
template
final_check_status theory_arith::process_non_linear() {
m_model_depends_on_computed_epsilon = false;
if (m_nl_monomials.empty())
return FC_DONE;
if (!reflection_enabled())
return FC_GIVEUP;
if (check_monomial_assignments()) {
return FC_DONE;
}
if (!m_params.m_nl_arith) {
TRACE("non_linear", tout << "Non-linear is not enabled\n";);
return FC_GIVEUP;
}
TRACE("process_non_linear", display(tout););
if (m_nl_rounds > m_params.m_nl_arith_rounds) {
TRACE("non_linear", tout << "GIVEUP non linear problem...\n";);
IF_VERBOSE(3, verbose_stream() << "Max. non linear arithmetic rounds. Increase threshold using NL_ARITH_ROUNDS=\n";);
return FC_GIVEUP;
}
ctx.push_trail(value_trail(m_nl_rounds));
m_nl_rounds++;
elim_quasi_base_rows();
move_non_base_vars_to_bounds();
TRACE("non_linear_verbose", tout << "processing non linear constraints...\n"; ctx.display(tout););
if (!make_feasible()) {
TRACE("non_linear", tout << "failed to move variables to bounds.\n";);
failed();
return FC_CONTINUE;
}
if (!max_min_nl_vars())
return FC_CONTINUE;
if (check_monomial_assignments()) {
return m_liberal_final_check || !m_changed_assignment ? FC_DONE : FC_CONTINUE;
}
svector vars;
get_non_linear_cluster(vars);
bool progress;
unsigned old_idx = m_nl_strategy_idx;
ctx.push_trail(value_trail(m_nl_strategy_idx));
do {
progress = false;
switch (m_nl_strategy_idx) {
case 0:
if (propagate_nl_bounds()) {
propagate_core();
progress = true;
}
break;
case 1:
if (m_params.m_nl_arith_cross_nested && !is_cross_nested_consistent(vars))
progress = true;
break;
case 2:
if (m_params.m_nl_arith_gb) {
switch(compute_grobner(vars)) {
case GB_PROGRESS:
progress = true;
break;
case GB_NEW_EQ:
progress = true;
propagate_core();
break;
case GB_FAIL:
break;
}
}
break;
case 3:
if (m_params.m_nl_arith_branching) {
theory_var target = find_nl_var_for_branching();
if (target != null_theory_var && branch_nl_int_var(target))
progress = true;
}
break;
}
m_nl_strategy_idx = (m_nl_strategy_idx + 1) % 4;
if (progress)
return FC_CONTINUE;
}
while (m_nl_strategy_idx != old_idx);
if (check_monomial_assignments()) {
return m_liberal_final_check || !m_changed_assignment ? FC_DONE : FC_CONTINUE;
}
TRACE("non_linear", display(tout););
return FC_GIVEUP;
}
};