mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-10-06 07:53:59 +00:00
Code Activity

Merge branch 'master' into polysat

Browse source
This commit is contained in:
Jakob Rath 2023-02-01 16:28:57 +01:00
commit 20b5455d08
669 changed files with 26145 additions and 20652 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/math/interval/dep_intervals.h
View file

@ -172,6 +172,7 @@ public:
void set_upper_is_inf(interval& a, bool inf) const { m_config.set_upper_is_inf(a, inf); }
void set_lower_dep(interval& a, u_dependency* d) const { m_config.set_lower_dep(a, d); }
void set_upper_dep(interval& a, u_dependency* d) const { m_config.set_upper_dep(a, d); }
void reset(interval& a) const { set_lower_is_inf(a, true); set_upper_is_inf(a, true); }
void set_value(interval& a, rational const& n) const {
set_lower(a, n);
set_upper(a, n);
@ -331,6 +332,7 @@ public:
}
mpq const& lower(interval const& a) const { return m_config.lower(a); }
mpq const& upper(interval const& a) const { return m_config.upper(a); }
bool is_empty(interval const& a) const;
void set_interval_for_scalar(interval&, const rational&);
template <typename T>

308
src/math/interval/mod_interval.h
View file

@ -1,5 +1,5 @@
/*++
Copyright (c) 2014 Microsoft Corporation
Copyright (c) 2017 Microsoft Corporation
Module Name:
@ -7,105 +7,235 @@ Module Name:
Abstract:
Intervals over fixed precision modular arithmetic
Modular interval for bit-vector comparisons
Author:
Nikolaj Bjorner (nbjorner) 2021-03-19
Jakob Rath 2021-04-6
Nikolaj and Nuno
--*/
#pragma once
#include "util/rational.h"
namespace bv {
template<typename T, typename Base>
struct interval_tpl : public Base {
T l, h;
unsigned sz = 0;
bool tight = true;
interval_tpl(T const& l, T const& h, unsigned sz, bool tight = false): l(l), h(h), sz(sz), tight(tight) {}
interval_tpl() {}
bool invariant() const {
return
0 <= l && (l <= Base::bound(sz)) &&
0 <= h && (h <= Base::bound(sz)) &&
(!is_wrapped() || l != h + 1);
}
bool is_full() const {
return l == 0 && h == Base::bound(sz);
}
bool is_wrapped() const { return l > h; }
bool is_singleton() const { return l == h; }
bool operator==(const interval_tpl<T, Base>& b) const {
SASSERT(sz == b.sz);
return l == b.l && h == b.h && tight == b.tight;
}
bool operator!=(const interval_tpl& b) const { return !(*this == b); }
bool implies(const interval_tpl<T, Base>& b) const {
if (b.is_full())
return true;
else if (is_full())
return false;
else if (is_wrapped())
// l >= b.l >= b.h >= h
return b.is_wrapped() && h <= b.h && l >= b.l;
else if (b.is_wrapped())
// b.l > b.h >= h >= l
// h >= l >= b.l > b.h
return h <= b.h || l >= b.l;
else
return l >= b.l && h <= b.h;
}
/// return false if intersection is unsat
bool intersect(const interval_tpl<T, Base>& b, interval_tpl& result) const {
if (is_full() || *this == b) {
result = b;
return true;
}
if (b.is_full()) {
result = *this;
return true;
}
if (is_wrapped()) {
if (b.is_wrapped()) {
if (h >= b.l)
result = b;
else if (b.h >= l)
result = *this;
else
result = interval_tpl(std::max(l, b.l), std::min(h, b.h), sz);
}
else
return b.intersect(*this, result);
}
else if (b.is_wrapped()) {
// ... b.h ... l ... h ... b.l ..
if (h < b.l && l > b.h)
return false;
// ... l ... b.l ... h ...
if (h >= b.l && l <= b.h)
result = b;
else if (h >= b.l)
result = interval_tpl(b.l, h, sz);
else {
// ... l .. b.h .. h .. b.l ...
SASSERT(l <= b.h);
result = interval_tpl(l, std::min(h, b.h), sz);
}
} else {
if (l > b.h || h < b.l)
return false;
// 0 .. l.. l' ... h ... h'
result = interval_tpl(std::max(l, b.l), std::min(h, b.h), sz, tight && b.tight);
}
return true;
}
/// return false if negation is empty
bool negate(interval_tpl<T, Base>& result) const {
if (!tight)
result = interval_tpl(Base::zero(), Base::bound(sz), sz, true);
else if (is_full())
return false;
else if (l == 0 && Base::bound(sz) == h)
result = interval_tpl(Base::zero(), Base::bound(sz), sz);
else if (l == 0)
result = interval_tpl(h + 1, Base::bound(sz), sz);
else if (Base::bound(sz) == h)
result = interval_tpl(Base::zero(), l - 1, sz);
else
result = interval_tpl(h + 1, l - 1, sz);
return true;
}
template<typename Numeral>
struct pp {
Numeral n;
pp(Numeral const& n):n(n) {}
};
};
struct rinterval_base {
static rational bound(unsigned sz) {
return rational::power_of_two(sz) - 1;
}
static rational zero() { return rational::zero(); }
};
struct rinterval : public interval_tpl<rational, rinterval_base> {
rinterval(rational const& l, rational const& h, unsigned sz, bool tight = false) {
this->l = l; this->h = h; this->sz = sz; this->tight = tight;
}
rinterval() { l = 0; h = 0; tight = true; }
};
struct iinterval_base {
static uint64_t uMaxInt(unsigned sz) {
SASSERT(sz <= 64);
return ULLONG_MAX >> (64u - sz);
}
static uint64_t bound(unsigned sz) { return uMaxInt(sz); }
static uint64_t zero() { return 0; }
};
struct iinterval : public interval_tpl<uint64_t, iinterval_base> {
iinterval(uint64_t l, uint64_t h, unsigned sz, bool tight = false) {
this->l = l; this->h = h; this->sz = sz; this->tight = tight;
}
iinterval() { l = 0; h = 0; sz = 0; tight = true; }
};
struct interval {
bool is_small = true;
iinterval i;
rinterval r;
interval() {}
interval(rational const& l, rational const& h, unsigned sz, bool tight = false) {
if (sz <= 64) {
is_small = true;
i.l = l.get_uint64();
i.h = h.get_uint64();
i.tight = tight;
i.sz = sz;
}
else {
is_small = false;
r.l = l;
r.h = h;
r.tight = tight;
r.sz = sz;
}
}
unsigned size() const {
return is_small ? i.sz : r.sz;
}
bool negate(interval& result) const {
result.is_small = is_small;
if (is_small)
return i.negate(result.i);
else
return r.negate(result.r);
}
bool intersect(interval const& b, interval & result) const {
result.is_small = is_small;
SASSERT(b.is_small == is_small);
if (is_small)
return i.intersect(b.i, result.i);
else
return r.intersect(b.r, result.r);
}
bool operator==(interval const& other) const {
SASSERT(is_small == other.is_small);
return is_small ? i == other.i : r == other.r;
}
bool operator!=(interval const& other) const {
return !(*this == other);
}
bool is_singleton() const { return is_small ? i.is_singleton() : r.is_singleton(); }
bool is_full() const { return is_small ? i.is_full() : r.is_full(); }
bool tight() const { return is_small ? i.tight : r.tight; }
bool implies(const interval& b) const {
SASSERT(is_small == b.is_small);
return is_small ? i.implies(b.i) : r.implies(b.r);
}
rational lo() const { return is_small ? rational(i.l, rational::ui64()) : r.l; }
rational hi() const { return is_small ? rational(i.h, rational::ui64()) : r.h; }
};
inline std::ostream& operator<<(std::ostream& out, pp<uint8_t> const& p) {
return out << (unsigned)p.n;
}
template<typename Numeral>
inline std::ostream& operator<<(std::ostream& out, pp<Numeral> const& p) {
if ((Numeral)(0 - p.n) < p.n)
return out << "-" << (Numeral)(0 - p.n);
return out << p.n;
}
inline std::ostream& operator<<(std::ostream& out, pp<rational> const& p) {
return out << p.n;
}
template<typename Numeral>
class mod_interval {
bool emp = false;
public:
Numeral lo { 0 };
Numeral hi { 0 };
mod_interval() {}
mod_interval(Numeral const& l, Numeral const& h): lo(l), hi(h) {}
virtual ~mod_interval() {}
static mod_interval free() { return mod_interval(0, 0); }
static mod_interval empty() { mod_interval i(0, 0); i.emp = true; return i; }
bool is_free() const { return !emp && lo == hi; }
bool is_empty() const { return emp; }
bool is_singleton() const { return !is_empty() && (lo + 1 == hi || (hi == 0 && is_max(lo))); }
bool contains(Numeral const& n) const;
bool contains(mod_interval const& other) const;
virtual bool is_max(Numeral const& n) const { return (Numeral)(n + 1) == 0; }
Numeral max() const;
Numeral min() const;
void set_free() { lo = hi = 0; emp = false; }
void set_bounds(Numeral const& l, Numeral const& h) { lo = l; hi = h; }
void set_empty() { emp = true; }
mod_interval& intersect_ule(Numeral const& h);
mod_interval& intersect_uge(Numeral const& l);
mod_interval& intersect_ult(Numeral const& h);
mod_interval& intersect_ugt(Numeral const& l);
mod_interval& intersect_fixed(Numeral const& n);
mod_interval& intersect_diff(Numeral const& n);
mod_interval& update_lo(Numeral const& new_lo);
mod_interval& update_hi(Numeral const& new_hi);
mod_interval operator&(mod_interval const& other) const;
mod_interval operator+(mod_interval const& other) const;
mod_interval operator-(mod_interval const& other) const;
mod_interval operator*(mod_interval const& other) const;
mod_interval operator-() const;
mod_interval operator*(Numeral const& n) const;
mod_interval operator+(Numeral const& n) const { return mod_interval(lo + n, hi + n); }
mod_interval operator-(Numeral const& n) const { return mod_interval(lo - n, hi - n); }
mod_interval& operator+=(mod_interval const& other) { *this = *this + other; return *this; }
std::ostream& display(std::ostream& out) const {
if (is_empty()) return out << "empty";
if (is_free()) return out << "free";
return out << "[" << pp(lo) << ", " << pp(hi) << "[";
inline std::ostream& operator<<(std::ostream& o, const interval& I) {
if (I.is_small)
return o << "[" << I.i.l << ", " << I.i.h << "]";
else
return o << "[" << I.r.l << ", " << I.r.h << "]";
}
Numeral closest_value(Numeral const& n) const;
bool operator==(mod_interval const& other) const {
if (is_empty())
return other.is_empty();
if (is_free())
return other.is_free();
return lo == other.lo && hi == other.hi;
}
bool operator!=(mod_interval const& other) const {
return !(*this == other);
}
};
template<typename Numeral>
inline std::ostream& operator<<(std::ostream& out, mod_interval<Numeral> const& i) {
return i.display(out);
}

2
src/math/lp/CMakeLists.txt
View file

@ -34,10 +34,12 @@ z3_add_component(lp
nla_basics_lemmas.cpp
nla_common.cpp
nla_core.cpp
nla_divisions.cpp
nla_grobner.cpp
nla_intervals.cpp
nla_monotone_lemmas.cpp
nla_order_lemmas.cpp
nla_powers.cpp
nla_solver.cpp
nla_tangent_lemmas.cpp
nra_solver.cpp

101
src/math/lp/emonics.cpp
View file

@ -40,43 +40,39 @@ void emonics::push() {
TRACE("nla_solver_mons", display(tout << "push\n"););
SASSERT(invariant());
m_u_f_stack.push_scope();
m_lim.push_back(m_monics.size());
m_region.push_scope();
m_ve.push();
SASSERT(monics_are_canonized());
SASSERT(invariant());
}
void emonics::pop_monic() {
m_ve.pop(1);
monic& m = m_monics.back();
TRACE("nla_solver_mons", display(tout << m << "\n"););
remove_cg_mon(m);
m_var2index[m.var()] = UINT_MAX;
do_canonize(m);
// variables in vs are in the same state as they were when add was called
lpvar last_var = UINT_MAX;
for (lpvar v : m.rvars()) {
if (v != last_var) {
remove_cell(m_use_lists[v]);
last_var = v;
}
}
m_ve.pop(1);
m_monics.pop_back();
}
void emonics::pop(unsigned n) {
TRACE("nla_solver_mons", tout << "pop: " << n << "\n";);
SASSERT(invariant());
for (unsigned j = 0; j < n; ++j) {
unsigned old_sz = m_lim[m_lim.size() - 1];
for (unsigned i = m_monics.size(); i-- > old_sz; ) {
m_ve.pop(1);
monic & m = m_monics[i];
TRACE("nla_solver_mons", display(tout << m << "\n"););
remove_cg_mon(m);
m_var2index[m.var()] = UINT_MAX;
do_canonize(m);
// variables in vs are in the same state as they were when add was called
lpvar last_var = UINT_MAX;
for (lpvar v : m.rvars()) {
if (v != last_var) {
remove_cell(m_use_lists[v]);
last_var = v;
}
}
m_ve.pop(1);
}
m_ve.pop(1);
m_monics.shrink(old_sz);
m_region.pop_scope(1);
m_lim.pop_back();
for (unsigned i = 0; i < n; ++i) {
m_u_f_stack.pop_scope(1);
SASSERT(invariant());
SASSERT(monics_are_canonized());
m_ve.pop(1);
}
SASSERT(invariant());
SASSERT(monics_are_canonized());
}
void emonics::remove_cell(head_tail& v) {
@ -96,7 +92,7 @@ void emonics::remove_cell(head_tail& v) {
void emonics::insert_cell(head_tail& v, unsigned mIndex) {
cell*& cur_head = v.m_head;
cell*& cur_tail = v.m_tail;
cell* new_head = new (m_region) cell(mIndex, cur_head);
cell* new_head = new (m_u_f_stack.get_region()) cell(mIndex, cur_head);
cur_head = new_head;
if (!cur_tail) cur_tail = new_head;
cur_tail->m_next = new_head;
@ -331,6 +327,14 @@ void emonics::add(lpvar v, unsigned sz, lpvar const* vs) {
m_monics.push_back(monic(v, sz, vs, idx));
do_canonize(m_monics.back());
class pop_mon : public trail {
emonics& p;
public:
pop_mon(emonics& p) :p(p) {}
void undo() override { p.pop_monic(); }
};
m_u_f_stack.push(pop_mon(*this));
// variables in m_vs are canonical and sorted,
// so use last_var to skip duplicates,
// while updating use-lists
@ -351,9 +355,8 @@ void emonics::add(lpvar v, unsigned sz, lpvar const* vs) {
void emonics::do_canonize(monic & m) const {
TRACE("nla_solver_mons", tout << m << "\n";);
m.reset_rfields();
for (lpvar v : m.vars()) {
m.push_rvar(m_ve.find(v));
}
for (lpvar v : m.vars())
m.push_rvar(m_ve.find(v));
m.sort_rvars();
TRACE("nla_solver_mons", tout << m << "\n";);
}
@ -365,40 +368,34 @@ bool emonics::is_canonized(const monic & m) const {
}
void emonics::ensure_canonized() {
for (auto & m : m_monics) {
do_canonize(m);
}
for (auto & m : m_monics)
do_canonize(m);
}
bool emonics::monics_are_canonized() const {
for (auto & m: m_monics) {
if (!is_canonized(m)) {
return false;
}
}
for (auto & m: m_monics)
if (!is_canonized(m))
return false;
return true;
}
bool emonics::canonize_divides(monic& m, monic & n) const {
if (m.size() > n.size()) return false;
if (m.size() > n.size())
return false;
unsigned ms = m.size(), ns = n.size();
unsigned i = 0, j = 0;
while (true) {
if (i == ms) {
return true;
}
else if (j == ns) {
return false;
}
if (i == ms)
return true;
else if (j == ns)
return false;
else if (m.rvars()[i] == n.rvars()[j]) {
++i; ++j;
}
else if (m.rvars()[i] < n.rvars()[j]) {
return false;
}
else {
++j;
}
else if (m.rvars()[i] < n.rvars()[j])
return false;
else
++j;
}
}

6
src/math/lp/emonics.h
View file

@ -87,15 +87,15 @@ class emonics {
var_eqs<emonics>& m_ve;
mutable vector<monic> m_monics; // set of monics
mutable unsigned_vector m_var2index; // var_mIndex -> mIndex
unsigned_vector m_lim; // backtracking point
mutable unsigned m_visited; // timestamp of visited monics during pf_iterator
region m_region; // region for allocating linked lists
mutable svector<head_tail> m_use_lists; // use list of monics where variables occur.
hash_canonical m_cg_hash;
eq_canonical m_cg_eq;
map<lpvar, unsigned_vector, hash_canonical, eq_canonical> m_cg_table; // congruence (canonical) table.
void pop_monic();
void inc_visited() const;
void remove_cell(head_tail& v);
@ -115,6 +115,8 @@ class emonics {
std::ostream& display_use(std::ostream& out) const;
std::ostream& display_uf(std::ostream& out) const;
std::ostream& display(std::ostream& out, cell* c) const;
public:
unsigned number_of_monics() const { return m_monics.size(); }
/**

2
src/math/lp/factorization.cpp
View file

@ -23,7 +23,7 @@ bool const_iterator_mon::get_factors(factor& k, factor& j, rational& sign) const
std::sort(k_vars.begin(), k_vars.end());
std::sort(j_vars.begin(), j_vars.end());
if (false && m_num_failures > 10) {
if (m_num_failures > 1000) {
for (bool& m : m_mask) m = true;
m_mask[0] = false;
m_full_factorization_returned = true;

44
src/math/lp/lar_solver.h
View file

@ -295,70 +295,86 @@ class lar_solver : public column_namer {
public:
// this function just looks at the status
bool is_feasible() const;
const map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>& fixed_var_table_int() const {
return m_fixed_var_table_int;
}
map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>& fixed_var_table_int() {
return m_fixed_var_table_int;
}
const map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>& fixed_var_table_real() const {
return m_fixed_var_table_real;
}
map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>& fixed_var_table_real() {
return m_fixed_var_table_real;
}
bool find_in_fixed_tables(const rational& mpq, bool is_int, unsigned & j) const {
return is_int? fixed_var_table_int().find(mpq, j) :
fixed_var_table_real().find(mpq, j);
return is_int? fixed_var_table_int().find(mpq, j) : fixed_var_table_real().find(mpq, j);
}
template <typename T> void remove_non_fixed_from_table(T&);
unsigned external_to_column_index(unsigned) const;
bool inside_bounds(lpvar, const impq&) const;
inline void set_column_value(unsigned j, const impq& v) {
m_mpq_lar_core_solver.m_r_solver.update_x(j, v);
}
inline void set_column_value_test(unsigned j, const impq& v) {
set_column_value(j, v);
}
unsigned get_total_iterations() const;
var_index add_named_var(unsigned ext_j, bool is_integer, const std::string&);
lp_status maximize_term(unsigned j_or_term, impq &term_max);
inline
core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) const { return
core_solver_pretty_printer<lp::mpq, lp::impq>(m_mpq_lar_core_solver.m_r_solver, out); }
inline core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) const {
return core_solver_pretty_printer<lp::mpq, lp::impq>(m_mpq_lar_core_solver.m_r_solver, out);
}
void get_infeasibility_explanation(explanation &) const;
inline void backup_x() { m_backup_x = m_mpq_lar_core_solver.m_r_x; }
inline void restore_x() { m_mpq_lar_core_solver.m_r_x = m_backup_x; }
template <typename T>
void explain_implied_bound(const implied_bound & ib, lp_bound_propagator<T> & bp) {
unsigned i = ib.m_row_or_term_index;
int bound_sign = ib.m_is_lower_bound? 1: -1;
int j_sign = (ib.m_coeff_before_j_is_pos ? 1 :-1) * bound_sign;
int bound_sign = (ib.m_is_lower_bound ? 1 : -1);
int j_sign = (ib.m_coeff_before_j_is_pos ? 1 : -1) * bound_sign;
unsigned bound_j = ib.m_j;
if (tv::is_term(bound_j)) {
if (tv::is_term(bound_j))
bound_j = m_var_register.external_to_local(bound_j);
}
for (auto const& r : A_r().m_rows[i]) {
for (auto const& r : get_row(i)) {
unsigned j = r.var();
if (j == bound_j) continue;
if (j == bound_j)
continue;
mpq const& a = r.coeff();
int a_sign = is_pos(a)? 1: -1;
int a_sign = is_pos(a) ? 1 : -1;
int sign = j_sign * a_sign;
const ul_pair & ul = m_columns_to_ul_pairs[j];
auto witness = sign > 0? ul.upper_bound_witness(): ul.lower_bound_witness();
auto witness = sign > 0 ? ul.upper_bound_witness() : ul.lower_bound_witness();
lp_assert(is_valid(witness));
bp.consume(a, witness);
}
}
void set_value_for_nbasic_column(unsigned j, const impq& new_val);
inline unsigned get_base_column_in_row(unsigned row_index) const {
return m_mpq_lar_core_solver.m_r_solver.get_base_column_in_row(row_index);
}
// lp_assert(implied_bound_is_correctly_explained(ib, explanation)); }
constraint_index mk_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side);
void activate_check_on_equal(constraint_index, var_index&);

2
src/math/lp/lp_primal_core_solver_def.h
View file

@ -480,7 +480,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_and_t(unsigned entering, X & t) {
if (this->m_settings.use_breakpoints_in_feasibility_search && !this->current_x_is_feasible())
return find_leaving_and_t_with_breakpoints(entering, t);
X theta;
X theta = zero_of_type<X>();
bool unlimited = get_harris_theta(theta);
lp_assert(unlimited || theta >= zero_of_type<X>());
if (try_jump_to_another_bound_on_entering(entering, theta, t, unlimited)) return entering;

4
src/math/lp/mps_reader.h
View file

@ -277,8 +277,8 @@ class mps_reader {
} else {
fail:
set_m_ok_to_false();
*m_message_stream << "cannot understand this line" << std::endl;
*m_message_stream << "line = " << m_line << ", line number is " << m_line_number << std::endl;
*m_message_stream << "cannot understand this line\n"
"line = " << m_line << ", line number is " << m_line_number << std::endl;
return;
}
}

26
src/math/lp/nla_core.cpp
View file

@ -29,6 +29,8 @@ core::core(lp::lar_solver& s, reslimit & lim) :
m_basics(this),
m_order(this),
m_monotone(this),
m_powers(*this),
m_divisions(*this),
m_intervals(this, lim),
m_monomial_bounds(this),
m_horner(this),
@ -120,9 +122,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;
}
@ -158,7 +159,13 @@ rational core::product_value(const monic& m) const {
}
// return true iff the monic value is equal to the product of the values of the factors
// or if the variable associated with the monomial is not relevant.
bool core::check_monic(const monic& m) const {
#if 0
// TODO test this
if (!is_relevant(m.var()))
return true;
#endif
SASSERT((!m_lar_solver.column_is_int(m.var())) || m_lar_solver.get_column_value(m.var()).is_int());
bool ret = product_value(m) == m_lar_solver.get_column_value(m.var()).x;
CTRACE("nla_solver_check_monic", !ret, print_monic(m, tout) << '\n';);
@ -978,6 +985,9 @@ bool core::rm_check(const monic& rm) const {
return check_monic(m_emons[rm.var()]);
}
bool core::has_relevant_monomial() const {
return any_of(emons(), [&](auto const& m) { return is_relevant(m.var()); });
}
bool core::find_bfc_to_refine_on_monic(const monic& m, factorization & bf) {
for (auto f : factorization_factory_imp(m, *this)) {
@ -1477,6 +1487,15 @@ void core::check_weighted(unsigned sz, std::pair<unsigned, std::function<void(vo
}
}
lbool core::check_power(lpvar r, lpvar x, lpvar y, vector<lemma>& l_vec) {
m_lemma_vec = &l_vec;
return m_powers.check(r, x, y, l_vec);
}
void core::check_bounded_divisions(vector<lemma>& l_vec) {
m_lemma_vec = &l_vec;
m_divisions.check_bounded_divisions();
}
lbool core::check(vector<lemma>& l_vec) {
lp_settings().stats().m_nla_calls++;
@ -1515,6 +1534,9 @@ lbool core::check(vector<lemma>& l_vec) {
if (l_vec.empty() && !done())
m_basics.basic_lemma(false);
if (l_vec.empty() && !done())
m_divisions.check();
#if 0
if (l_vec.empty() && !done() && !run_horner)
m_horner.horner_lemmas();

130
src/math/lp/nla_core.h
View file

@ -19,6 +19,8 @@
#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/nla_divisions.h"
#include "math/lp/emonics.h"
#include "math/lp/nla_settings.h"
#include "math/lp/nex.h"
@ -42,104 +44,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<ineq> 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<ineq>& ineqs() const { return m_ineqs; }
vector<ineq>& 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 +53,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;
@ -177,12 +82,15 @@ class core {
lp::lar_solver& m_lar_solver;
reslimit& m_reslim;
std::function<bool(lpvar)> m_relevant;
vector<lemma> * m_lemma_vec;
lp::u_set m_to_refine;
tangents m_tangents;
basics m_basics;
order m_order;
monotone m_monotone;
powers m_powers;
divisions m_divisions;
intervals m_intervals;
monomial_bounds m_monomial_bounds;
nla_settings m_nla_settings;
@ -204,6 +112,9 @@ class core {
void check_weighted(unsigned sz, std::pair<unsigned, std::function<void(void)>>* checks);
public:
// constructor
core(lp::lar_solver& s, reslimit&);
void insert_to_refine(lpvar j);
void erase_from_to_refine(lpvar j);
@ -212,9 +123,7 @@ public:
void insert_to_active_var_set(unsigned j) const { m_active_var_set.insert(j); }
void clear_active_var_set() const {
m_active_var_set.clear();
}
void clear_active_var_set() const { m_active_var_set.clear(); }
void clear_and_resize_active_var_set() const {
m_active_var_set.clear();
@ -226,9 +135,9 @@ public:
reslimit& reslim() { return m_reslim; }
emonics& emons() { return m_emons; }
const emonics& emons() const { return m_emons; }
// constructor
core(lp::lar_solver& s, reslimit &);
bool has_relevant_monomial() const;
bool compare_holds(const rational& ls, llc cmp, const rational& rs) const;
rational value(const lp::lar_term& r) const;
@ -294,9 +203,18 @@ public:
void deregister_monic_from_tables(const monic & m, unsigned i);
void add_monic(lpvar v, unsigned sz, lpvar const* vs);
void add_idivision(lpvar q, lpvar x, lpvar y) { m_divisions.add_idivision(q, x, y); }
void add_rdivision(lpvar q, lpvar x, lpvar y) { m_divisions.add_rdivision(q, x, y); }
void add_bounded_division(lpvar q, lpvar x, lpvar y) { m_divisions.add_bounded_division(q, x, y); }
void set_relevant(std::function<bool(lpvar)>& is_relevant) { m_relevant = is_relevant; }
bool is_relevant(lpvar v) const { return !m_relevant || m_relevant(v); }
void push();
void pop(unsigned n);
trail_stack& trail() { return m_emons.get_trail_stack(); }
rational mon_value_by_vars(unsigned i) const;
rational product_value(const monic & m) const;
@ -463,7 +381,9 @@ public:
bool conflict_found() const;
lbool check(vector<lemma>& l_vec);
lbool check(vector<lemma>& l_vec);
lbool check_power(lpvar r, lpvar x, lpvar y, vector<lemma>& l_vec);
void check_bounded_divisions(vector<lemma>&);
bool no_lemmas_hold() const;

208
src/math/lp/nla_divisions.cpp Normal file
View file

@ -0,0 +1,208 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
nla_divisions.cpp
Author:
Lev Nachmanson (levnach)
Nikolaj Bjorner (nbjorner)
Description:
Check divisions
--*/
#include "math/lp/nla_core.h"
namespace nla {
void divisions::add_idivision(lpvar q, lpvar x, lpvar y) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
return;
if (lp::tv::is_term(x) || lp::tv::is_term(y) || lp::tv::is_term(q))
return;
m_idivisions.push_back({q, x, y});
m_core.trail().push(push_back_vector(m_idivisions));
}
void divisions::add_rdivision(lpvar q, lpvar x, lpvar y) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
return;
if (lp::tv::is_term(x) || lp::tv::is_term(y) || lp::tv::is_term(q))
return;
m_rdivisions.push_back({ q, x, y });
m_core.trail().push(push_back_vector(m_rdivisions));
}
void divisions::add_bounded_division(lpvar q, lpvar x, lpvar y) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
return;
if (lp::tv::is_term(x) || lp::tv::is_term(y) || lp::tv::is_term(q))
return;
m_bounded_divisions.push_back({ q, x, y });
m_core.trail().push(push_back_vector(m_bounded_divisions));
}
typedef lp::lar_term term;
// y1 >= y2 > 0 & x1 <= x2 => x1/y1 <= x2/y2
// y2 <= y1 < 0 & x1 >= x2 >= 0 => x1/y1 <= x2/y2
// y2 <= y1 < 0 & x1 <= x2 <= 0 => x1/y1 >= x2/y2
void divisions::check() {
core& c = m_core;
if (c.use_nra_model())
return;
auto monotonicity1 = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& q1, auto& q1val,
auto x2, auto& x2val, auto y2, auto& y2val, auto& q2, auto& q2val) {
if (y1val >= y2val && y2val > 0 && x1val <= x2val && q1val > q2val) {
new_lemma lemma(c, "y1 >= y2 > 0 & x1 <= x2 => x1/y1 <= x2/y2");
lemma |= ineq(term(y1, rational(-1), y2), llc::LT, 0);
lemma |= ineq(y2, llc::LE, 0);
lemma |= ineq(term(x1, rational(-1), x2), llc::GT, 0);
lemma |= ineq(term(q1, rational(-1), q2), llc::LE, 0);
return true;
}
return false;
};
auto monotonicity2 = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& q1, auto& q1val,
auto x2, auto& x2val, auto y2, auto& y2val, auto& q2, auto& q2val) {
if (y2val <= y1val && y1val < 0 && x1val >= x2val && x2val >= 0 && q1val > q2val) {
new_lemma lemma(c, "y2 <= y1 < 0 & x1 >= x2 >= 0 => x1/y1 <= x2/y2");
lemma |= ineq(term(y1, rational(-1), y2), llc::LT, 0);
lemma |= ineq(y1, llc::GE, 0);
lemma |= ineq(term(x1, rational(-1), x2), llc::LT, 0);
lemma |= ineq(x2, llc::LT, 0);
lemma |= ineq(term(q1, rational(-1), q2), llc::LE, 0);
return true;
}
return false;
};
auto monotonicity3 = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& q1, auto& q1val,
auto x2, auto& x2val, auto y2, auto& y2val, auto& q2, auto& q2val) {
if (y2val <= y1val && y1val < 0 && x1val <= x2val && x2val <= 0 && q1val < q2val) {
new_lemma lemma(c, "y2 <= y1 < 0 & x1 <= x2 <= 0 => x1/y1 >= x2/y2");
lemma |= ineq(term(y1, rational(-1), y2), llc::LT, 0);
lemma |= ineq(y1, llc::GE, 0);
lemma |= ineq(term(x1, rational(-1), x2), llc::GT, 0);
lemma |= ineq(x2, llc::GT, 0);
lemma |= ineq(term(q1, rational(-1), q2), llc::GE, 0);
return true;
}
return false;
};
auto monotonicity = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& q1, auto& q1val,
auto x2, auto& x2val, auto y2, auto& y2val, auto& q2, auto& q2val) {
if (monotonicity1(x1, x1val, y1, y1val, q1, q1val, x2, x2val, y2, y2val, q2, q2val))
return true;
if (monotonicity1(x2, x2val, y2, y2val, q2, q2val, x1, x1val, y1, y1val, q1, q1val))
return true;
if (monotonicity2(x1, x1val, y1, y1val, q1, q1val, x2, x2val, y2, y2val, q2, q2val))
return true;
if (monotonicity2(x2, x2val, y2, y2val, q2, q2val, x1, x1val, y1, y1val, q1, q1val))
return true;
if (monotonicity3(x1, x1val, y1, y1val, q1, q1val, x2, x2val, y2, y2val, q2, q2val))
return true;
if (monotonicity3(x2, x2val, y2, y2val, q2, q2val, x1, x1val, y1, y1val, q1, q1val))
return true;
return false;
};
for (auto const & [r, x, y] : m_idivisions) {
if (!c.is_relevant(r))
continue;
auto xval = c.val(x);
auto yval = c.val(y);
auto rval = c.val(r);
// idiv semantics
if (!xval.is_int() || !yval.is_int() || yval == 0 || rval == div(xval, yval))
continue;
for (auto const& [q2, x2, y2] : m_idivisions) {
if (q2 == r)
continue;
if (!c.is_relevant(q2))
continue;
auto x2val = c.val(x2);
auto y2val = c.val(y2);
auto q2val = c.val(q2);
if (monotonicity(x, xval, y, yval, r, rval, x2, x2val, y2, y2val, q2, q2val))
return;
}
}
for (auto const& [r, x, y] : m_rdivisions) {
if (!c.is_relevant(r))
continue;
auto xval = c.val(x);
auto yval = c.val(y);
auto rval = c.val(r);
// / semantics
if (yval == 0 || rval == xval / yval)
continue;
for (auto const& [q2, x2, y2] : m_rdivisions) {
if (q2 == r)
continue;
if (!c.is_relevant(q2))
continue;
auto x2val = c.val(x2);
auto y2val = c.val(y2);
auto q2val = c.val(q2);
if (monotonicity(x, xval, y, yval, r, rval, x2, x2val, y2, y2val, q2, q2val))
return;
}
}
}
// if p is bounded, q a value, r = eval(p):
// p <= q * div(r, q) + q - 1 => div(p, q) <= div(r, q)
// p >= q * div(r, q) => div(r, q) <= div(p, q)
void divisions::check_bounded_divisions() {
core& c = m_core;
unsigned offset = c.random(), sz = m_bounded_divisions.size();
for (unsigned j = 0; j < sz; ++j) {
unsigned i = (offset + j) % sz;
auto [q, x, y] = m_bounded_divisions[i];
if (!c.is_relevant(q))
continue;
auto xv = c.val(x);
auto yv = c.val(y);
auto qv = c.val(q);
if (xv < 0 || !xv.is_int())
continue;
if (yv <= 0 || !yv.is_int())
continue;
if (qv == div(xv, yv))
continue;
rational div_v = div(xv, yv);
// y = yv & x <= yv * div(xv, yv) + yv - 1 => div(x, y) <= div(xv, yv)
// y = yv & x >= y * div(xv, yv) => div(xv, yv) <= div(x, y)
rational mul(1);
rational hi = yv * div_v + yv - 1;
rational lo = yv * div_v;
if (xv > hi) {
new_lemma lemma(c, "y = yv & x <= yv * div(xv, yv) + yv - 1 => div(p, y) <= div(xv, yv)");
lemma |= ineq(y, llc::NE, yv);
lemma |= ineq(x, llc::GT, hi);
lemma |= ineq(q, llc::LE, div_v);
return;
}
if (xv < lo) {
new_lemma lemma(c, "y = yv & x >= yv * div(xv, yv) => div(xv, yv) <= div(x, y)");
lemma |= ineq(y, llc::NE, yv);
lemma |= ineq(x, llc::LT, lo);
lemma |= ineq(q, llc::GE, div_v);
return;
}
}
}
}

37
src/math/lp/nla_divisions.h Normal file
View file

@ -0,0 +1,37 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
nla_divisions.h
Author:
Lev Nachmanson (levnach)
Nikolaj Bjorner (nbjorner)
Description:
Check division constraints.
--*/
#include "math/lp/nla_types.h"
namespace nla {
class core;
class divisions {
core& m_core;
vector<std::tuple<lpvar, lpvar, lpvar>> m_idivisions;
vector<std::tuple<lpvar, lpvar, lpvar>> m_rdivisions;
vector<std::tuple<lpvar, lpvar, lpvar>> m_bounded_divisions;
public:
divisions(core& c):m_core(c) {}
void add_idivision(lpvar q, lpvar x, lpvar y);
void add_rdivision(lpvar q, lpvar x, lpvar y);
void add_bounded_division(lpvar q, lpvar x, lpvar y);
void check();
void check_bounded_divisions();
};
}

182
src/math/lp/nla_powers.cpp Normal file
View file

@ -0,0 +1,182 @@
/*++
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 a.
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! = 1 + ln(a)*y + ln(a)^2*y^2/2 +
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<lemma>& lemmas) {
if (x == null_lpvar || y == null_lpvar || r == null_lpvar)
return l_undef;
core& c = m_core;
if (c.use_nra_model())
return l_undef;
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) {
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 && yval < 0 && rval >= 1) {
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 && yval > 0 && rval <= 1) {
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 >= 3 && yval != 0 & rval <= yval + 1) {
new_lemma lemma(c, "x >= 3, y != 0 => x^y > ln(x)y + 1");
lemma |= ineq(x, llc::LT, rational(3));
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;
}
}

29
src/math/lp/nla_powers.h Normal file
View file

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

130
src/math/lp/nla_solver.cpp
View file

@ -13,63 +13,91 @@
#include "math/lp/factorization.h"
#include "math/lp/nla_solver.h"
#include "math/lp/nla_core.h"
#include "math/polynomial/algebraic_numbers.h"
namespace nla {
nla_settings& solver::settings() { return m_core->m_nla_settings; }
nla_settings& solver::settings() { return m_core->m_nla_settings; }
void solver::add_monic(lpvar v, unsigned sz, lpvar const* vs) {
m_core->add_monic(v, sz, vs);
}
bool solver::is_monic_var(lpvar v) const {
return m_core->is_monic_var(v);
}
bool solver::need_check() { return true; }
lbool solver::check(vector<lemma>& l) {
return m_core->check(l);
}
void solver::push(){
m_core->push();
}
void solver::pop(unsigned n) {
m_core->pop(n);
}
solver::solver(lp::lar_solver& s, reslimit& limit):
m_core(alloc(core, s, limit)) {
}
bool solver::influences_nl_var(lpvar j) const {
return m_core->influences_nl_var(j);
}
solver::~solver() {
dealloc(m_core);
}
std::ostream& solver::display(std::ostream& out) const {
m_core->print_monics(out);
if( use_nra_model()) {
m_core->m_nra.display(out);
void solver::add_monic(lpvar v, unsigned sz, lpvar const* vs) {
m_core->add_monic(v, sz, vs);
}
return out;
}
bool solver::use_nra_model() const { return m_core->use_nra_model(); }
core& solver::get_core() { return *m_core; }
nlsat::anum_manager& solver::am() { return m_core->m_nra.am(); }
nlsat::anum const& solver::am_value(lp::var_index v) const {
SASSERT(use_nra_model());
return m_core->m_nra.value(v);
}
void solver::add_idivision(lpvar q, lpvar x, lpvar y) {
m_core->add_idivision(q, x, y);
}
void solver::collect_statistics(::statistics & st) {
m_core->collect_statistics(st);
}
void solver::add_rdivision(lpvar q, lpvar x, lpvar y) {
m_core->add_rdivision(q, x, y);
}
void solver::add_bounded_division(lpvar q, lpvar x, lpvar y) {
m_core->add_bounded_division(q, x, y);
}
void solver::set_relevant(std::function<bool(lpvar)>& is_relevant) {
m_core->set_relevant(is_relevant);
}
bool solver::is_monic_var(lpvar v) const {
return m_core->is_monic_var(v);
}
bool solver::need_check() { return m_core->has_relevant_monomial(); }
lbool solver::check(vector<lemma>& l) {
return m_core->check(l);
}
void solver::push(){
m_core->push();
}
void solver::pop(unsigned n) {
m_core->pop(n);
}
solver::solver(lp::lar_solver& s, reslimit& limit):
m_core(alloc(core, s, limit)) {
}
bool solver::influences_nl_var(lpvar j) const {
return m_core->influences_nl_var(j);
}
solver::~solver() {
dealloc(m_core);
}
std::ostream& solver::display(std::ostream& out) const {
m_core->print_monics(out);
if (use_nra_model())
m_core->m_nra.display(out);
return out;
}
bool solver::use_nra_model() const { return m_core->use_nra_model(); }
core& solver::get_core() { return *m_core; }
nlsat::anum_manager& solver::am() { return m_core->m_nra.am(); }
nlsat::anum const& solver::am_value(lp::var_index v) const {
SASSERT(use_nra_model());
return m_core->m_nra.value(v);
}
void solver::collect_statistics(::statistics & st) {
m_core->collect_statistics(st);
}
// ensure r = x^y, add abstraction/refinement lemmas
lbool solver::check_power(lpvar r, lpvar x, lpvar y, vector<lemma>& lemmas) {
return m_core->check_power(r, x, y, lemmas);
}
void solver::check_bounded_divisions(vector<lemma>& lemmas) {
m_core->check_bounded_divisions(lemmas);
}
}

54
src/math/lp/nla_solver.h
View file

@ -16,29 +16,37 @@ Author:
#include "math/lp/nla_settings.h"
#include "math/lp/nla_core.h"
namespace nra {
class solver;
class solver;
}
namespace nla {
class core;
// nonlinear integer incremental linear solver
class solver {
core* m_core;
public:
void add_monic(lpvar v, unsigned sz, lpvar const* vs);
solver(lp::lar_solver& s, reslimit& limit);
~solver();
nla_settings& settings();
void push();
void pop(unsigned scopes);
bool need_check();
lbool check(vector<lemma>&);
bool is_monic_var(lpvar) const;
bool influences_nl_var(lpvar) const;
std::ostream& display(std::ostream& out) const;
bool use_nra_model() const;
core& get_core();
nlsat::anum_manager& am();
nlsat::anum const& am_value(lp::var_index v) const;
void collect_statistics(::statistics & st);
};
class core;
// nonlinear integer incremental linear solver
class solver {
core* m_core;
public:
solver(lp::lar_solver& s, reslimit& limit);
~solver();
void add_monic(lpvar v, unsigned sz, lpvar const* vs);
void add_idivision(lpvar q, lpvar x, lpvar y);
void add_rdivision(lpvar q, lpvar x, lpvar y);
void add_bounded_division(lpvar q, lpvar x, lpvar y);
void check_bounded_divisions(vector<lemma>&);
void set_relevant(std::function<bool(lpvar)>& is_relevant);
nla_settings& settings();
void push();
void pop(unsigned scopes);
bool need_check();
lbool check(vector<lemma>&);
lbool check_power(lpvar r, lpvar x, lpvar y, vector<lemma>&);
bool is_monic_var(lpvar) const;
bool influences_nl_var(lpvar) const;
std::ostream& display(std::ostream& out) const;
bool use_nra_model() const;
core& get_core();
nlsat::anum_manager& am();
nlsat::anum const& am_value(lp::var_index v) const;
void collect_statistics(::statistics & st);
};
}

120
src/math/lp/nla_types.h Normal file
View file

@ -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<ineq> 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<ineq>& ineqs() const { return m_ineqs; }
vector<ineq>& 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) {}
};
}

3
src/math/lp/stacked_vector.h
View file

@ -55,6 +55,9 @@ public:
bool operator==(B const& other) const {
return m_vec.m_vector[m_i] == other;
}
bool operator!=(B const& other) const {
return m_vec.m_vector[m_i] != other;
}
B& operator+=(B const &delta) {
// not tracking the change here!
return m_vec.m_vector[m_i] += delta;

3
src/math/lp/var_eqs.h
View file

@ -68,8 +68,7 @@ class var_eqs {
T* m_merge_handler;
union_find<var_eqs> m_uf;
lp::incremental_vector<std::pair<signed_var, signed_var>>
m_trail;
lp::incremental_vector<std::pair<signed_var, signed_var>> m_trail;
vector<svector<eq_edge>> m_eqs; // signed_var.index() -> the edges adjacent to signed_var.index()
trail_stack m_stack;

2
src/math/polynomial/upolynomial_factorization.cpp
View file

@ -27,7 +27,7 @@ Notes:
#include "math/polynomial/upolynomial_factorization_int.h"
#include "util/prime_generator.h"
using namespace std;
using std::endl;
namespace upolynomial {

33
src/math/realclosure/realclosure.cpp
View file

@ -506,8 +506,7 @@ namespace realclosure {
m_bqim(lim, m_bqm),
m_plus_inf_approx(m_bqm),
m_minus_inf_approx(m_bqm) {
mpq one(1);
m_one = mk_rational(one);
m_one = mk_rational(mpq(1));
inc_ref(m_one);
m_pi = nullptr;
m_e = nullptr;
@ -2557,13 +2556,10 @@ namespace realclosure {
return new (allocator()) rational_value();
}
/**
\brief Make a rational and swap its value with v
*/
rational_value * mk_rational_and_swap(mpq & v) {
rational_value * mk_rational(mpq && v) {
SASSERT(!qm().is_zero(v));
rational_value * r = mk_rational();
::swap(r->m_value, v);
r->m_value = std::move(v);
return r;
}
@ -2585,7 +2581,7 @@ namespace realclosure {
SASSERT(!bqm().is_zero(v));
scoped_mpq v_q(qm()); // v as a rational
::to_mpq(qm(), v, v_q);
return mk_rational(v_q);
return mk_rational(std::move(v_q));
}
void reset_interval(value * a) {
@ -3270,7 +3266,7 @@ namespace realclosure {
scoped_mpq num_z(qm());
qm().div(lcm_z, to_mpq(dens[i]), num_z);
SASSERT(qm().is_int(num_z));
m = mk_rational_and_swap(num_z);
m = mk_rational(std::move(num_z));
is_z = true;
}
bool found_lt_eq = false;
@ -3432,7 +3428,7 @@ namespace realclosure {
scoped_mpq r(qm());
SASSERT(qm().is_int(to_mpq(a)));
qm().div(to_mpq(a), b, r);
a = mk_rational_and_swap(r);
a = mk_rational(std::move(r));
}
else {
rational_function_value * rf = to_rational_function(a);
@ -3592,9 +3588,8 @@ namespace realclosure {
r.reset();
if (sz > 1) {
for (unsigned i = 1; i < sz; i++) {
mpq i_mpq(i);
value_ref a_i(*this);
a_i = mk_rational_and_swap(i_mpq);
a_i = mk_rational(mpq(i));
mul(a_i, p[i], a_i);
r.push_back(a_i);
}
@ -3821,7 +3816,7 @@ namespace realclosure {
scoped_mpz mpz_twok(qm());
qm().mul2k(mpz(1), b.k(), mpz_twok);
value_ref twok(*this), twok_i(*this);
twok = mk_rational(mpz_twok);
twok = mk_rational(std::move(mpz_twok));
twok_i = twok;
value_ref c(*this);
c = mk_rational(b.numerator());
@ -5061,7 +5056,7 @@ namespace realclosure {
if (qm().is_zero(v))
r = nullptr;
else
r = mk_rational_and_swap(v);
r = mk_rational(std::move(v));
}
else {
INC_DEPTH();
@ -5090,7 +5085,7 @@ namespace realclosure {
if (qm().is_zero(v))
r = nullptr;
else
r = mk_rational_and_swap(v);
r = mk_rational(std::move(v));
}
else {
value_ref neg_b(*this);
@ -5124,7 +5119,7 @@ namespace realclosure {
scoped_mpq v(qm());
qm().set(v, to_mpq(a));
qm().neg(v);
r = mk_rational_and_swap(v);
r = mk_rational(std::move(v));
}
else {
neg_rf(to_rational_function(a), r);
@ -5269,7 +5264,7 @@ namespace realclosure {
else if (is_nz_rational(a) && is_nz_rational(b)) {
scoped_mpq v(qm());
qm().mul(to_mpq(a), to_mpq(b), v);
r = mk_rational_and_swap(v);
r = mk_rational(std::move(v));
}
else {
INC_DEPTH();
@ -5304,7 +5299,7 @@ namespace realclosure {
else if (is_nz_rational(a) && is_nz_rational(b)) {
scoped_mpq v(qm());
qm().div(to_mpq(a), to_mpq(b), v);
r = mk_rational_and_swap(v);
r = mk_rational(std::move(v));
}
else {
value_ref inv_b(*this);
@ -5557,7 +5552,7 @@ namespace realclosure {
if (is_nz_rational(a)) {
scoped_mpq v(qm());
qm().inv(to_mpq(a), v);
r = mk_rational_and_swap(v);
r = mk_rational(std::move(v));
}
else {
inv_rf(to_rational_function(a), r);

146
src/math/simplex/model_based_opt.cpp
View file

@ -98,11 +98,13 @@ namespace opt {
}
else if (v1 < v2) {
vs.push_back(vs1[i]);
vs.back().m_coeff *= c1;
vs.back().m_coeff *= c1;
++i;
}
else {
vs.push_back(vs2[j]);
vs.back().m_coeff *= c2;
vs.back().m_coeff *= c2;
++j;
}
}
result.m_div = c1*m_div;
@ -111,6 +113,65 @@ namespace opt {
return result;
}
/**
a1*x1 + a2*x2 + a3*x3 + coeff1 / c1
x2 |-> b1*x1 + b4*x4 + ceoff2 / c2
------------------------------------------------------------------------
(a1*x1 + a2*((b1*x1 + b4*x4 + coeff2) / c2) + a3*x3 + coeff1) / c1
------------------------------------------------------------------------
(c2*a1*x1 + a2*b1*x1 + a2*b4*x4 + c2*a3*x3 + c2*coeff1 + coeff2) / c1*c2
*/
void model_based_opt::def::substitute(unsigned v, def const& other) {
vector<var> const& vs1 = m_vars;
rational coeff(0);
for (auto const& [id, c] : vs1) {
if (id == v) {
coeff = c;
break;
}
}
if (coeff == 0)
return;
rational c1 = m_div;
rational c2 = other.m_div;
vector<var> const& vs2 = other.m_vars;
vector<var> vs;
unsigned i = 0, j = 0;
while (i < vs1.size() || j < vs2.size()) {
unsigned v1 = UINT_MAX, v2 = UINT_MAX;
if (i < vs1.size()) v1 = vs1[i].m_id;
if (j < vs2.size()) v2 = vs2[j].m_id;
if (v1 == v)
++i;
else if (v1 == v2) {
vs.push_back(vs1[i]);
vs.back().m_coeff *= c2;
vs.back().m_coeff += coeff * vs2[j].m_coeff;
++i; ++j;
if (vs.back().m_coeff.is_zero())
vs.pop_back();
}
else if (v1 < v2) {
vs.push_back(vs1[i]);
vs.back().m_coeff *= c2;
++i;
}
else {
vs.push_back(vs2[j]);
vs.back().m_coeff *= coeff;
++j;
}
}
m_div *= other.m_div;
m_coeff *= c2;
m_coeff += coeff*other.m_coeff;
m_vars.reset();
m_vars.append(vs);
normalize();
}
model_based_opt::def model_based_opt::def::operator/(rational const& r) const {
def result(*this);
result.m_div *= r;
@ -953,12 +1014,14 @@ namespace opt {
return dst;
}
// -x + lo <= 0
void model_based_opt::add_lower_bound(unsigned x, rational const& lo) {
vector<var> coeffs;
coeffs.push_back(var(x, rational::minus_one()));
add_constraint(coeffs, lo, t_le);
}
// x - hi <= 0
void model_based_opt::add_upper_bound(unsigned x, rational const& hi) {
vector<var> coeffs;
coeffs.push_back(var(x, rational::one()));
@ -1238,7 +1301,7 @@ namespace opt {
def result;
unsigned_vector div_rows(_div_rows), mod_rows(_mod_rows);
SASSERT(!div_rows.empty() || !mod_rows.empty());
TRACE("opt", display(tout << "solve_div " << x << "\n"));
TRACE("opt", display(tout << "solve_div v" << x << "\n"));
rational K(1);
for (unsigned ri : div_rows)
@ -1381,8 +1444,9 @@ namespace opt {
for (unsigned ri : mod_rows) {
rational a = get_coefficient(ri, x);
replace_var(ri, x, rational::zero());
rational rMod = m_rows[ri].m_mod;
// add w = b mod K
// add w = b mod rMod
vector<var> coeffs = m_rows[ri].m_vars;
rational coeff = m_rows[ri].m_coeff;
unsigned v = m_rows[ri].m_id;
@ -1390,16 +1454,46 @@ namespace opt {
unsigned w = UINT_MAX;
rational offset(0);
if (coeffs.empty() || K == 1)
offset = mod(coeff, K);
if (coeffs.empty() || rMod == 1)
offset = mod(coeff, rMod);
else
w = add_mod(coeffs, coeff, K);
w = add_mod(coeffs, coeff, rMod);
rational w_value = w == UINT_MAX ? offset : m_var2value[w];
// add v = a*z + w - V, for k = (a*z_value + w_value) div K
// claim: (= (mod x K) (- x (* K (div x K)))))) is a theorem for every x, K != 0
#if 0
// V := (a * z_value + w_value) div rMod
// V*rMod <= a*z + w < (V+1)*rMod
// v = a*z + w - V*rMod
SASSERT(a > 0);
SASSERT(z_value >= 0);
SASSERT(w_value >= 0);
SASSERT(a * z_value + w_value >= 0);
rational V = div(a * z_value + w_value, rMod);
vector<var> mod_coeffs;
SASSERT(V >= 0);
SASSERT(a * z_value + w_value >= V*rMod);
SASSERT((V+1)*rMod > a*z_value + w_value);
// -a*z - w + V*rMod <= 0
mod_coeffs.push_back(var(z, -a));
if (w != UINT_MAX) mod_coeffs.push_back(var(w, -rational::one()));
add_constraint(mod_coeffs, V*rMod - offset, t_le);
mod_coeffs.reset();
// a*z + w - (V+1)*rMod + 1 <= 0
mod_coeffs.push_back(var(z, a));
if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one()));
add_constraint(mod_coeffs, -(V+1)*rMod + offset + 1, t_le);
mod_coeffs.reset();
// -v + a*z + w - V*rMod = 0
mod_coeffs.push_back(var(v, rational::minus_one()));
mod_coeffs.push_back(var(z, a));
if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one()));
add_constraint(mod_coeffs, offset - V*rMod, t_eq);
#else
// add v = a*z + w - V, for V = v_value - a * z_value - w_value
// claim: (= (mod x rMod) (- x (* rMod (div x rMod)))))) is a theorem for every x, rMod != 0
rational V = v_value - a * z_value - w_value;
vector<var> mod_coeffs;
mod_coeffs.push_back(var(v, rational::minus_one()));
@ -1407,24 +1501,34 @@ namespace opt {
if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one()));
add_constraint(mod_coeffs, V + offset, t_eq);
add_lower_bound(v, rational::zero());
add_upper_bound(v, K - 1);
add_upper_bound(v, rMod - 1);
#endif
retire_row(ri);
vs.push_back(v);
}
for (unsigned v : vs)
project(v, false);
for (unsigned v : vs) {
def v_def = project(v, compute_def);
if (compute_def)
eliminate(v, v_def);
}
// project internal variables.
def y_def = project(y, compute_def);
def z_def = project(z, compute_def);
def y_def = project(y, compute_def); // may depend on z
if (compute_def) {
z_def.substitute(y, y_def);
eliminate(y, y_def);
eliminate(z, z_def);
result = (y_def * K) + z_def;
m_var2value[x] = eval(result);
TRACE("opt", tout << y << " := " << y_def << "\n";
tout << z << " := " << z_def << "\n";
tout << x << " := " << result << "\n");
}
TRACE("opt", display(tout << "solve_div done v" << x << "\n"));
return result;
@ -1624,14 +1728,20 @@ namespace opt {
TRACE("opt", display(tout << "solved v" << x << "\n"));
return result;
}
void model_based_opt::eliminate(unsigned v, def const& new_def) {
for (auto & d : m_result)
d.substitute(v, new_def);
}
vector<model_based_opt::def> model_based_opt::project(unsigned num_vars, unsigned const* vars, bool compute_def) {
vector<def> result;
m_result.reset();
for (unsigned i = 0; i < num_vars; ++i) {
result.push_back(project(vars[i], compute_def));
m_result.push_back(project(vars[i], compute_def));
eliminate(vars[i], m_result.back());
TRACE("opt", display(tout << "After projecting: v" << vars[i] << "\n"););
}
return result;
return m_result;
}
}

4
src/math/simplex/model_based_opt.h
View file

@ -86,6 +86,7 @@ namespace opt {
def operator/(rational const& n) const;
def operator*(rational const& n) const;
def operator+(rational const& n) const;
void substitute(unsigned v, def const& other);
void normalize();
};
@ -100,6 +101,9 @@ namespace opt {
unsigned_vector m_lub, m_glb, m_divides, m_mod, m_div;
unsigned_vector m_above, m_below;
unsigned_vector m_retired_rows;
vector<model_based_opt::def> m_result;
void eliminate(unsigned v, def const& d);
bool invariant();
bool invariant(unsigned index, row const& r);

2
src/math/subpaving/tactic/subpaving_tactic.h
View file

@ -15,6 +15,8 @@ Author:
Revision History:
## Tactic subpaving
--*/
#pragma once