From 03385bf78d79a6eb884ab2d29a5ad685b89037f2 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Fri, 12 Aug 2022 10:20:43 -0400
Subject: [PATCH 01/16] improve quantifier elimination for arithmetic
This update changes the handling of mod and adds support for nested div terms.
Simple use cases that are handled using small results are given below.
```
(declare-const x Int)
(declare-const y Int)
(declare-const z Int)
(assert (exists ((x Int)) (and (<= y (* 10 x)) (<= (* 10 x) z))))
(apply qe2)
(reset)
(declare-const y Int)
(assert (exists ((x Int)) (and (> x 0) (= (div x 41) y))))
(apply qe2)
(reset)
(declare-const y Int)
(assert (exists ((x Int)) (= (mod x 41) y)))
(apply qe2)
(reset)
```
The main idea is to introduce definition rows for mod/div terms.
Elimination of variables under mod/div is defined by rewriting the variable to multiples of the mod/divisior and remainder.
The functionality is disabled in this push.
---
src/math/simplex/model_based_opt.cpp | 2864 ++++++++++++++------------
src/math/simplex/model_based_opt.h | 38 +-
src/qe/mbp/mbp_arith.cpp | 325 +--
3 files changed, 1795 insertions(+), 1432 deletions(-)
diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index 0223e8bff..be2b4fe2b 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -1,1300 +1,1564 @@
-/*++
-Copyright (c) 2016 Microsoft Corporation
-
-Module Name:
-
- model_based_opt.cpp
-
-Abstract:
-
- Model-based optimization and projection for linear real, integer arithmetic.
-
-Author:
-
- Nikolaj Bjorner (nbjorner) 2016-27-4
-
-Revision History:
-
-
---*/
-
-#include "math/simplex/model_based_opt.h"
-#include "util/uint_set.h"
-#include "util/z3_exception.h"
-
-std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) {
- switch (ie) {
- case opt::t_eq: return out << " = ";
- case opt::t_lt: return out << " < ";
- case opt::t_le: return out << " <= ";
- case opt::t_mod: return out << " mod ";
- }
- return out;
-}
-
-
-namespace opt {
-
- /**
- * Convert a row ax + coeffs + coeff = value into a definition for x
- * x = (value - coeffs - coeff)/a
- * as backdrop we have existing assignments to x and other variables that
- * satisfy the equality with value, and such that value satisfies
- * the row constraint ( = , <= , < , mod)
- */
- model_based_opt::def::def(row const& r, unsigned x) {
- for (var const & v : r.m_vars) {
- if (v.m_id != x) {
- m_vars.push_back(v);
- }
- else {
- m_div = -v.m_coeff;
- }
- }
- m_coeff = r.m_coeff;
- switch (r.m_type) {
- case opt::t_lt:
- m_coeff += m_div;
- break;
- case opt::t_le:
- // for: ax >= t, then x := (t + a - 1) div a
- if (m_div.is_pos()) {
- m_coeff += m_div;
- m_coeff -= rational::one();
- }
- break;
- default:
- break;
- }
- normalize();
- SASSERT(m_div.is_pos());
- }
-
- model_based_opt::def model_based_opt::def::operator+(def const& other) const {
- def result;
- vector<var> const& vs1 = m_vars;
- vector<var> const& vs2 = other.m_vars;
- vector<var> & vs = result.m_vars;
- rational c1(1), c2(1);
- if (m_div != other.m_div) {
- c1 = other.m_div;
- c2 = m_div;
- }
- 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 == v2) {
- vs.push_back(vs1[i]);
- vs.back().m_coeff *= c1;
- vs.back().m_coeff += c2 * 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 *= c1;
- }
- else {
- vs.push_back(vs2[j]);
- vs.back().m_coeff *= c2;
- }
- }
- result.m_div = c1*m_div;
- result.m_coeff = (m_coeff*c1) + (other.m_coeff*c2);
- result.normalize();
- return result;
- }
-
- model_based_opt::def model_based_opt::def::operator/(rational const& r) const {
- def result(*this);
- result.m_div *= r;
- result.normalize();
- return result;
- }
-
- model_based_opt::def model_based_opt::def::operator*(rational const& n) const {
- def result(*this);
- for (var& v : result.m_vars) {
- v.m_coeff *= n;
- }
- result.m_coeff *= n;
- result.normalize();
- return result;
- }
-
- model_based_opt::def model_based_opt::def::operator+(rational const& n) const {
- def result(*this);
- result.m_coeff += n * result.m_div;
- result.normalize();
- return result;
- }
-
- void model_based_opt::def::normalize() {
- if (!m_div.is_int()) {
- rational den = denominator(m_div);
- SASSERT(den > 1);
- for (var& v : m_vars)
- v.m_coeff *= den;
- m_coeff *= den;
- m_div *= den;
-
- }
- if (m_div.is_neg()) {
- for (var& v : m_vars)
- v.m_coeff.neg();
- m_coeff.neg();
- m_div.neg();
- }
- if (m_div.is_one())
- return;
- rational g(m_div);
- if (!m_coeff.is_int())
- return;
- g = gcd(g, m_coeff);
- for (var const& v : m_vars) {
- if (!v.m_coeff.is_int())
- return;
- g = gcd(g, abs(v.m_coeff));
- if (g.is_one())
- break;
- }
- if (!g.is_one()) {
- for (var& v : m_vars)
- v.m_coeff /= g;
- m_coeff /= g;
- m_div /= g;
- }
- }
-
- model_based_opt::model_based_opt() {
- m_rows.push_back(row());
- }
-
- bool model_based_opt::invariant() {
- for (unsigned i = 0; i < m_rows.size(); ++i) {
- if (!invariant(i, m_rows[i])) {
- return false;
- }
- }
- return true;
- }
-
-#define PASSERT(_e_) { CTRACE("qe", !(_e_), display(tout, r); display(tout);); SASSERT(_e_); }
-
- bool model_based_opt::invariant(unsigned index, row const& r) {
- vector<var> const& vars = r.m_vars;
- for (unsigned i = 0; i < vars.size(); ++i) {
- // variables in each row are sorted and have non-zero coefficients
- PASSERT(i + 1 == vars.size() || vars[i].m_id < vars[i+1].m_id);
- PASSERT(!vars[i].m_coeff.is_zero());
- PASSERT(index == 0 || m_var2row_ids[vars[i].m_id].contains(index));
- }
-
- PASSERT(r.m_value == eval(r));
- PASSERT(r.m_type != t_eq || r.m_value.is_zero());
- // values satisfy constraints
- PASSERT(index == 0 || r.m_type != t_lt || r.m_value.is_neg());
- PASSERT(index == 0 || r.m_type != t_le || !r.m_value.is_pos());
- PASSERT(index == 0 || r.m_type != t_mod || (mod(r.m_value, r.m_mod).is_zero()));
- return true;
- }
-
- // a1*x + obj
- // a2*x + t2 <= 0
- // a3*x + t3 <= 0
- // a4*x + t4 <= 0
- // a1 > 0, a2 > 0, a3 > 0, a4 < 0
- // x <= -t2/a2
- // x <= -t2/a3
- // determine lub among these.
- // then resolve lub with others
- // e.g., -t2/a2 <= -t3/a3, then
- // replace inequality a3*x + t3 <= 0 by -t2/a2 + t3/a3 <= 0
- // mark a4 as invalid.
- //
-
- // a1 < 0, a2 < 0, a3 < 0, a4 > 0
- // x >= t2/a2
- // x >= t3/a3
- // determine glb among these
- // the resolve glb with others.
- // e.g. t2/a2 >= t3/a3
- // then replace a3*x + t3 by t3/a3 - t2/a2 <= 0
- //
- inf_eps model_based_opt::maximize() {
- SASSERT(invariant());
- unsigned_vector bound_trail, bound_vars;
- TRACE("opt", display(tout << "tableau\n"););
- while (!objective().m_vars.empty()) {
- var v = objective().m_vars.back();
- unsigned x = v.m_id;
- rational const& coeff = v.m_coeff;
- unsigned bound_row_index;
- rational bound_coeff;
- if (find_bound(x, bound_row_index, bound_coeff, coeff.is_pos())) {
- SASSERT(!bound_coeff.is_zero());
- TRACE("opt", display(tout << "update: " << v << " ", objective());
- for (unsigned above : m_above) {
- display(tout << "resolve: ", m_rows[above]);
- });
- for (unsigned above : m_above) {
- resolve(bound_row_index, bound_coeff, above, x);
- }
- for (unsigned below : m_below) {
- resolve(bound_row_index, bound_coeff, below, x);
- }
- // coeff*x + objective <= ub
- // a2*x + t2 <= 0
- // => coeff*x <= -t2*coeff/a2
- // objective + t2*coeff/a2 <= ub
-
- mul_add(false, m_objective_id, - coeff/bound_coeff, bound_row_index);
- retire_row(bound_row_index);
- bound_trail.push_back(bound_row_index);
- bound_vars.push_back(x);
- }
- else {
- TRACE("opt", display(tout << "unbound: " << v << " ", objective()););
- update_values(bound_vars, bound_trail);
- return inf_eps::infinity();
- }
- }
-
- //
- // update the evaluation of variables to satisfy the bound.
- //
-
- update_values(bound_vars, bound_trail);
-
- rational value = objective().m_value;
- if (objective().m_type == t_lt) {
- return inf_eps(inf_rational(value, rational(-1)));
- }
- else {
- return inf_eps(inf_rational(value));
- }
- }
-
-
- void model_based_opt::update_value(unsigned x, rational const& val) {
- rational old_val = m_var2value[x];
- m_var2value[x] = val;
- SASSERT(val.is_int() || !is_int(x));
- unsigned_vector const& row_ids = m_var2row_ids[x];
- for (unsigned row_id : row_ids) {
- rational coeff = get_coefficient(row_id, x);
- if (coeff.is_zero()) {
- continue;
- }
- row & r = m_rows[row_id];
- rational delta = coeff * (val - old_val);
- r.m_value += delta;
- SASSERT(invariant(row_id, r));
- }
- }
-
-
- void model_based_opt::update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail) {
- for (unsigned i = bound_trail.size(); i-- > 0; ) {
- unsigned x = bound_vars[i];
- row& r = m_rows[bound_trail[i]];
- rational val = r.m_coeff;
- rational old_x_val = m_var2value[x];
- rational new_x_val;
- rational x_coeff, eps(0);
- vector<var> const& vars = r.m_vars;
- for (var const& v : vars) {
- if (x == v.m_id) {
- x_coeff = v.m_coeff;
- }
- else {
- val += m_var2value[v.m_id]*v.m_coeff;
- }
- }
- SASSERT(!x_coeff.is_zero());
- new_x_val = -val/x_coeff;
-
- if (r.m_type == t_lt) {
- eps = abs(old_x_val - new_x_val)/rational(2);
- eps = std::min(rational::one(), eps);
- SASSERT(!eps.is_zero());
-
- //
- // ax + t < 0
- // <=> x < -t/a
- // <=> x := -t/a - epsilon
- //
- if (x_coeff.is_pos()) {
- new_x_val -= eps;
- }
- //
- // -ax + t < 0
- // <=> -ax < -t
- // <=> -x < -t/a
- // <=> x > t/a
- // <=> x := t/a + epsilon
- //
- else {
- new_x_val += eps;
- }
- }
- TRACE("opt", display(tout << "v" << x
- << " coeff_x: " << x_coeff
- << " old_x_val: " << old_x_val
- << " new_x_val: " << new_x_val
- << " eps: " << eps << " ", r); );
- m_var2value[x] = new_x_val;
-
- r.m_value = eval(r);
- SASSERT(invariant(bound_trail[i], r));
- }
-
- // update and check bounds for all other affected rows.
- for (unsigned i = bound_trail.size(); i-- > 0; ) {
- unsigned x = bound_vars[i];
- unsigned_vector const& row_ids = m_var2row_ids[x];
- for (unsigned row_id : row_ids) {
- row & r = m_rows[row_id];
- r.m_value = eval(r);
- SASSERT(invariant(row_id, r));
- }
- }
- SASSERT(invariant());
- }
-
- bool model_based_opt::find_bound(unsigned x, unsigned& bound_row_index, rational& bound_coeff, bool is_pos) {
- bound_row_index = UINT_MAX;
- rational lub_val;
- rational const& x_val = m_var2value[x];
- unsigned_vector const& row_ids = m_var2row_ids[x];
- uint_set visited;
- m_above.reset();
- m_below.reset();
- for (unsigned row_id : row_ids) {
- SASSERT(row_id != m_objective_id);
- if (visited.contains(row_id)) {
- continue;
- }
- visited.insert(row_id);
- row& r = m_rows[row_id];
- if (r.m_alive) {
- rational a = get_coefficient(row_id, x);
- if (a.is_zero()) {
- // skip
- }
- else if (a.is_pos() == is_pos || r.m_type == t_eq) {
- rational value = x_val - (r.m_value/a);
- if (bound_row_index == UINT_MAX) {
- lub_val = value;
- bound_row_index = row_id;
- bound_coeff = a;
- }
- else if ((value == lub_val && r.m_type == opt::t_lt) ||
- (is_pos && value < lub_val) ||
-
- (!is_pos && value > lub_val)) {
- m_above.push_back(bound_row_index);
- lub_val = value;
- bound_row_index = row_id;
- bound_coeff = a;
- }
- else {
- m_above.push_back(row_id);
- }
- }
- else {
- m_below.push_back(row_id);
- }
- }
- }
- return bound_row_index != UINT_MAX;
- }
-
- void model_based_opt::retire_row(unsigned row_id) {
- m_rows[row_id].m_alive = false;
- m_retired_rows.push_back(row_id);
- }
-
- rational model_based_opt::eval(unsigned x) const {
- return m_var2value[x];
- }
-
- rational model_based_opt::eval(def const& d) const {
- vector<var> const& vars = d.m_vars;
- rational val = d.m_coeff;
- for (var const& v : vars) {
- val += v.m_coeff * eval(v.m_id);
- }
- val /= d.m_div;
- return val;
- }
-
- rational model_based_opt::eval(row const& r) const {
- vector<var> const& vars = r.m_vars;
- rational val = r.m_coeff;
- for (var const& v : vars) {
- val += v.m_coeff * eval(v.m_id);
- }
- return val;
- }
-
- rational model_based_opt::eval(vector<var> const& coeffs) const {
- rational val(0);
- for (var const& v : coeffs)
- val += v.m_coeff * eval(v.m_id);
- return val;
- }
-
- rational model_based_opt::get_coefficient(unsigned row_id, unsigned var_id) const {
- return m_rows[row_id].get_coefficient(var_id);
- }
-
- rational model_based_opt::row::get_coefficient(unsigned var_id) const {
- if (m_vars.empty()) {
- return rational::zero();
- }
- unsigned lo = 0, hi = m_vars.size();
- while (lo < hi) {
- unsigned mid = lo + (hi - lo)/2;
- SASSERT(mid < hi);
- unsigned id = m_vars[mid].m_id;
- if (id == var_id) {
- lo = mid;
- break;
- }
- if (id < var_id) {
- lo = mid + 1;
- }
- else {
- hi = mid;
- }
- }
- if (lo == m_vars.size()) {
- return rational::zero();
- }
- unsigned id = m_vars[lo].m_id;
- if (id == var_id) {
- return m_vars[lo].m_coeff;
- }
- else {
- return rational::zero();
- }
- }
-
- model_based_opt::row& model_based_opt::row::normalize() {
-#if 0
- if (m_type == t_mod)
- return *this;
- rational D(denominator(abs(m_coeff)));
- if (D == 0)
- D = 1;
- for (auto const& [id, coeff] : m_vars)
- if (coeff != 0)
- D = lcm(D, denominator(abs(coeff)));
- if (D == 1)
- return *this;
- SASSERT(D > 0);
- for (auto & [id, coeff] : m_vars)
- coeff *= D;
- m_coeff *= D;
-#endif
- return *this;
- }
-
- //
- // Let
- // row1: t1 + a1*x <= 0
- // row2: t2 + a2*x <= 0
- //
- // assume a1, a2 have the same signs:
- // (t2 + a2*x) <= (t1 + a1*x)*a2/a1
- // <=> t2*a1/a2 - t1 <= 0
- // <=> t2 - t1*a2/a1 <= 0
- //
- // assume a1 > 0, -a2 < 0:
- // t1 + a1*x <= 0, t2 - a2*x <= 0
- // t2/a2 <= -t1/a1
- // t2 + t1*a2/a1 <= 0
- // assume -a1 < 0, a2 > 0:
- // t1 - a1*x <= 0, t2 + a2*x <= 0
- // t1/a1 <= -t2/a2
- // t2 + t1*a2/a1 <= 0
- //
- // the resolvent is the same in all cases (simpler proof should exist)
- //
-
- void model_based_opt::resolve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) {
-
- SASSERT(a1 == get_coefficient(row_src, x));
- SASSERT(!a1.is_zero());
- SASSERT(row_src != row_dst);
-
- if (m_rows[row_dst].m_alive) {
- rational a2 = get_coefficient(row_dst, x);
- if (is_int(x)) {
- TRACE("opt",
- tout << a1 << " " << a2 << ": ";
- display(tout, m_rows[row_dst]);
- display(tout, m_rows[row_src]););
- if (a1.is_pos() != a2.is_pos() || m_rows[row_src].m_type == opt::t_eq) {
- mul_add(x, a1, row_src, a2, row_dst);
- }
- else {
- mul(row_dst, abs(a1));
- mul_add(false, row_dst, -abs(a2), row_src);
- }
- TRACE("opt", display(tout, m_rows[row_dst]););
- normalize(row_dst);
- }
- else {
- mul_add(row_dst != m_objective_id && a1.is_pos() == a2.is_pos(), row_dst, -a2/a1, row_src);
- }
- }
- }
-
- /**
- * a1 > 0
- * a1*x + r1 = value
- * a2*x + r2 <= 0
- * ------------------
- * a1*r2 - a2*r1 <= value
- */
- void model_based_opt::solve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) {
- SASSERT(a1 == get_coefficient(row_src, x));
- SASSERT(a1.is_pos());
- SASSERT(row_src != row_dst);
- if (!m_rows[row_dst].m_alive) return;
- rational a2 = get_coefficient(row_dst, x);
- mul(row_dst, a1);
- mul_add(false, row_dst, -a2, row_src);
- SASSERT(get_coefficient(row_dst, x).is_zero());
- }
-
- // resolution for integer rows.
- void model_based_opt::mul_add(
- unsigned x, rational const& src_c, unsigned row_src, rational const& dst_c, unsigned row_dst) {
- row& dst = m_rows[row_dst];
- row const& src = m_rows[row_src];
- SASSERT(is_int(x));
- SASSERT(t_le == dst.m_type && t_le == src.m_type);
- SASSERT(src_c.is_int());
- SASSERT(dst_c.is_int());
- SASSERT(m_var2value[x].is_int());
-
- rational abs_src_c = abs(src_c);
- rational abs_dst_c = abs(dst_c);
- rational x_val = m_var2value[x];
- rational slack = (abs_src_c - rational::one()) * (abs_dst_c - rational::one());
- rational dst_val = dst.m_value - x_val*dst_c;
- rational src_val = src.m_value - x_val*src_c;
- rational distance = abs_src_c * dst_val + abs_dst_c * src_val + slack;
- bool use_case1 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one();
-
-#if 0
- if (distance.is_nonpos() && !abs_src_c.is_one() && !abs_dst_c.is_one()) {
- unsigned r = copy_row(row_src);
- mul_add(false, r, rational::one(), row_dst);
- del_var(r, x);
- add(r, slack);
- TRACE("qe", tout << m_rows[r];);
- SASSERT(!m_rows[r].m_value.is_pos());
- }
-#endif
-
- if (use_case1) {
- TRACE("opt", tout << "slack: " << slack << " " << src_c << " " << dst_val << " " << dst_c << " " << src_val << "\n";);
- // dst <- abs_src_c*dst + abs_dst_c*src + slack
- mul(row_dst, abs_src_c);
- add(row_dst, slack);
- mul_add(false, row_dst, abs_dst_c, row_src);
- return;
- }
-
- //
- // create finite disjunction for |b|.
- // exists x, z in [0 .. |b|-2] . b*x + s + z = 0 && ax + t <= 0 && bx + s <= 0
- // <=>
- // exists x, z in [0 .. |b|-2] . b*x = -z - s && ax + t <= 0 && bx + s <= 0
- // <=>
- // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && bx + s <= 0
- // <=>
- // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z - s + s <= 0
- // <=>
- // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z <= 0
- // <=>
- // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0
- // <=>
- // exists x, z in [0 .. |b|-2] . b*x = -z - s && a*n_sign(b)(s + z) + |b|t <= 0
- // <=>
- // exists z in [0 .. |b|-2] . |b| | (z + s) && a*n_sign(b)(s + z) + |b|t <= 0
- //
-
- TRACE("qe", tout << "finite disjunction " << distance << " " << src_c << " " << dst_c << "\n";);
- vector<var> coeffs;
- if (abs_dst_c <= abs_src_c) {
- rational z = mod(dst_val, abs_dst_c);
- if (!z.is_zero()) z = abs_dst_c - z;
- mk_coeffs_without(coeffs, dst.m_vars, x);
- add_divides(coeffs, dst.m_coeff + z, abs_dst_c);
- add(row_dst, z);
- mul(row_dst, src_c * n_sign(dst_c));
- mul_add(false, row_dst, abs_dst_c, row_src);
- }
- else {
- // z := b - (s + bx) mod b
- // := b - s mod b
- // b | s + z <=> b | s + b - s mod b <=> b | s - s mod b
- rational z = mod(src_val, abs_src_c);
- if (!z.is_zero()) z = abs_src_c - z;
- mk_coeffs_without(coeffs, src.m_vars, x);
- add_divides(coeffs, src.m_coeff + z, abs_src_c);
- mul(row_dst, abs_src_c);
- add(row_dst, z * dst_c * n_sign(src_c));
- mul_add(false, row_dst, dst_c * n_sign(src_c), row_src);
- }
- }
-
- void model_based_opt::mk_coeffs_without(vector<var>& dst, vector<var> const& src, unsigned x) {
- for (var const & v : src) {
- if (v.m_id != x) dst.push_back(v);
- }
- }
-
- rational model_based_opt::n_sign(rational const& b) const {
- return rational(b.is_pos()?-1:1);
- }
-
- void model_based_opt::mul(unsigned dst, rational const& c) {
- if (c.is_one()) return;
- row& r = m_rows[dst];
- for (auto & v : r.m_vars) {
- v.m_coeff *= c;
- }
- r.m_coeff *= c;
- r.m_value *= c;
- }
-
- void model_based_opt::add(unsigned dst, rational const& c) {
- row& r = m_rows[dst];
- r.m_coeff += c;
- r.m_value += c;
- }
-
- void model_based_opt::sub(unsigned dst, rational const& c) {
- row& r = m_rows[dst];
- r.m_coeff -= c;
- r.m_value -= c;
- }
-
- void model_based_opt::del_var(unsigned dst, unsigned x) {
- row& r = m_rows[dst];
- unsigned j = 0;
- for (var & v : r.m_vars) {
- if (v.m_id == x) {
- r.m_value -= eval(x)*r.m_coeff;
- }
- else {
- r.m_vars[j++] = v;
- }
- }
- r.m_vars.shrink(j);
- }
-
-
- void model_based_opt::normalize(unsigned row_id) {
- row& r = m_rows[row_id];
- if (r.m_vars.empty()) {
- retire_row(row_id);
- return;
- }
- if (r.m_type == t_mod) return;
- rational g(abs(r.m_vars[0].m_coeff));
- bool all_int = g.is_int();
- for (unsigned i = 1; all_int && !g.is_one() && i < r.m_vars.size(); ++i) {
- rational const& coeff = r.m_vars[i].m_coeff;
- if (coeff.is_int()) {
- g = gcd(g, abs(coeff));
- }
- else {
- all_int = false;
- }
- }
- if (all_int && !r.m_coeff.is_zero()) {
- if (r.m_coeff.is_int()) {
- g = gcd(g, abs(r.m_coeff));
- }
- else {
- all_int = false;
- }
- }
- if (all_int && !g.is_one()) {
- SASSERT(!g.is_zero());
- mul(row_id, rational::one()/g);
- }
- }
-
- //
- // set row1 <- row1 + c*row2
- //
- void model_based_opt::mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2) {
- if (c.is_zero()) {
- return;
- }
-
- m_new_vars.reset();
- row& r1 = m_rows[row_id1];
- row const& r2 = m_rows[row_id2];
- unsigned i = 0, j = 0;
- while (i < r1.m_vars.size() || j < r2.m_vars.size()) {
- if (j == r2.m_vars.size()) {
- m_new_vars.append(r1.m_vars.size() - i, r1.m_vars.data() + i);
- break;
- }
- if (i == r1.m_vars.size()) {
- for (; j < r2.m_vars.size(); ++j) {
- m_new_vars.push_back(r2.m_vars[j]);
- m_new_vars.back().m_coeff *= c;
- if (row_id1 != m_objective_id) {
- m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
- }
- }
- break;
- }
-
- unsigned v1 = r1.m_vars[i].m_id;
- unsigned v2 = r2.m_vars[j].m_id;
- if (v1 == v2) {
- m_new_vars.push_back(r1.m_vars[i]);
- m_new_vars.back().m_coeff += c*r2.m_vars[j].m_coeff;
- ++i;
- ++j;
- if (m_new_vars.back().m_coeff.is_zero()) {
- m_new_vars.pop_back();
- }
- }
- else if (v1 < v2) {
- m_new_vars.push_back(r1.m_vars[i]);
- ++i;
- }
- else {
- m_new_vars.push_back(r2.m_vars[j]);
- m_new_vars.back().m_coeff *= c;
- if (row_id1 != m_objective_id) {
- m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
- }
- ++j;
- }
- }
- r1.m_coeff += c*r2.m_coeff;
- r1.m_vars.swap(m_new_vars);
- r1.m_value += c*r2.m_value;
-
- if (!same_sign && r2.m_type == t_lt) {
- r1.m_type = t_lt;
- }
- else if (same_sign && r1.m_type == t_lt && r2.m_type == t_lt) {
- r1.m_type = t_le;
- }
- SASSERT(invariant(row_id1, r1));
- }
-
- void model_based_opt::display(std::ostream& out) const {
- for (auto const& r : m_rows) {
- display(out, r);
- }
- for (unsigned i = 0; i < m_var2row_ids.size(); ++i) {
- unsigned_vector const& rows = m_var2row_ids[i];
- out << i << ": ";
- for (auto const& r : rows) {
- out << r << " ";
- }
- out << "\n";
- }
- }
-
- void model_based_opt::display(std::ostream& out, vector<var> const& vars, rational const& coeff) {
- unsigned i = 0;
- for (var const& v : vars) {
- if (i > 0 && v.m_coeff.is_pos()) {
- out << "+ ";
- }
- ++i;
- if (v.m_coeff.is_one()) {
- out << "v" << v.m_id << " ";
- }
- else {
- out << v.m_coeff << "*v" << v.m_id << " ";
- }
- }
- if (coeff.is_pos()) {
- out << " + " << coeff << " ";
- }
- else if (coeff.is_neg()) {
- out << coeff << " ";
- }
- }
-
- std::ostream& model_based_opt::display(std::ostream& out, row const& r) {
- out << (r.m_alive?"a":"d") << " ";
- display(out, r.m_vars, r.m_coeff);
- if (r.m_type == opt::t_mod) {
- out << r.m_type << " " << r.m_mod << " = 0; value: " << r.m_value << "\n";
- }
- else {
- out << r.m_type << " 0; value: " << r.m_value << "\n";
- }
- return out;
- }
-
- std::ostream& model_based_opt::display(std::ostream& out, def const& r) {
- display(out, r.m_vars, r.m_coeff);
- if (!r.m_div.is_one()) {
- out << " / " << r.m_div;
- }
- return out;
- }
-
- unsigned model_based_opt::add_var(rational const& value, bool is_int) {
- unsigned v = m_var2value.size();
- m_var2value.push_back(value);
- m_var2is_int.push_back(is_int);
- SASSERT(value.is_int() || !is_int);
- m_var2row_ids.push_back(unsigned_vector());
- return v;
- }
-
- rational model_based_opt::get_value(unsigned var) {
- return m_var2value[var];
- }
-
- void model_based_opt::set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel) {
- row& r = m_rows[row_id];
- rational val(c);
- SASSERT(r.m_vars.empty());
- r.m_vars.append(coeffs.size(), coeffs.data());
- bool is_int_row = !coeffs.empty();
- std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
- for (auto const& c : coeffs) {
- val += m_var2value[c.m_id] * c.m_coeff;
- SASSERT(!is_int(c.m_id) || c.m_coeff.is_int());
- is_int_row &= is_int(c.m_id);
- }
- r.m_alive = true;
- r.m_coeff = c;
- r.m_value = val;
- r.m_type = rel;
- r.m_mod = m;
- if (is_int_row && rel == t_lt) {
- r.m_type = t_le;
- r.m_coeff += rational::one();
- r.m_value += rational::one();
- }
- }
-
- unsigned model_based_opt::new_row() {
- unsigned row_id = 0;
- if (m_retired_rows.empty()) {
- row_id = m_rows.size();
- m_rows.push_back(row());
- }
- else {
- row_id = m_retired_rows.back();
- m_retired_rows.pop_back();
- m_rows[row_id].reset();
- m_rows[row_id].m_alive = true;
- }
- return row_id;
- }
-
- unsigned model_based_opt::copy_row(unsigned src) {
- unsigned dst = new_row();
- row const& r = m_rows[src];
- set_row(dst, r.m_vars, r.m_coeff, r.m_mod, r.m_type);
- for (auto const& v : r.m_vars) {
- m_var2row_ids[v.m_id].push_back(dst);
- }
- SASSERT(invariant(dst, m_rows[dst]));
- return dst;
- }
-
- void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, ineq_type rel) {
- add_constraint(coeffs, c, rational::zero(), rel);
- }
-
- void model_based_opt::add_divides(vector<var> const& coeffs, rational const& c, rational const& m) {
- add_constraint(coeffs, c, m, t_mod);
- }
-
- void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel) {
- unsigned row_id = new_row();
- set_row(row_id, coeffs, c, m, rel);
- for (var const& coeff : coeffs) {
- m_var2row_ids[coeff.m_id].push_back(row_id);
- }
- SASSERT(invariant(row_id, m_rows[row_id]));
- }
-
- void model_based_opt::set_objective(vector<var> const& coeffs, rational const& c) {
- set_row(m_objective_id, coeffs, c, rational::zero(), t_le);
- }
-
- void model_based_opt::get_live_rows(vector<row>& rows) {
- for (row & r : m_rows) {
- if (r.m_alive) {
- rows.push_back(r.normalize());
- }
- }
- }
-
- //
- // pick glb and lub representative.
- // The representative is picked such that it
- // represents the fewest inequalities.
- // The constraints that enforce a glb or lub are not forced.
- // The constraints that separate the glb from ub or the lub from lb
- // are not forced.
- // In other words, suppose there are
- // . N inequalities of the form t <= x
- // . M inequalities of the form s >= x
- // . t0 is glb among N under valuation.
- // . s0 is lub among M under valuation.
- // If N < M
- // create the inequalities:
- // t <= t0 for each t other than t0 (N-1 inequalities).
- // t0 <= s for each s (M inequalities).
- // If N >= M the construction is symmetric.
- //
- model_based_opt::def model_based_opt::project(unsigned x, bool compute_def) {
- unsigned_vector& lub_rows = m_lub;
- unsigned_vector& glb_rows = m_glb;
- unsigned_vector& mod_rows = m_mod;
- unsigned lub_index = UINT_MAX, glb_index = UINT_MAX;
- bool lub_strict = false, glb_strict = false;
- rational lub_val, glb_val;
- rational const& x_val = m_var2value[x];
- unsigned_vector const& row_ids = m_var2row_ids[x];
- uint_set visited;
- lub_rows.reset();
- glb_rows.reset();
- mod_rows.reset();
- bool lub_is_unit = true, glb_is_unit = true;
- unsigned eq_row = UINT_MAX;
- // select the lub and glb.
- for (unsigned row_id : row_ids) {
- if (visited.contains(row_id)) {
- continue;
- }
- visited.insert(row_id);
- row& r = m_rows[row_id];
- if (!r.m_alive) {
- continue;
- }
- rational a = get_coefficient(row_id, x);
- if (a.is_zero()) {
- continue;
- }
- if (r.m_type == t_eq) {
- eq_row = row_id;
- continue;
- }
- if (r.m_type == t_mod) {
- mod_rows.push_back(row_id);
- }
- else if (a.is_pos()) {
- rational lub_value = x_val - (r.m_value/a);
- if (lub_rows.empty() ||
- lub_value < lub_val ||
- (lub_value == lub_val && r.m_type == t_lt && !lub_strict)) {
- lub_val = lub_value;
- lub_index = row_id;
- lub_strict = r.m_type == t_lt;
- }
- lub_rows.push_back(row_id);
- lub_is_unit &= a.is_one();
- }
- else {
- SASSERT(a.is_neg());
- rational glb_value = x_val - (r.m_value/a);
- if (glb_rows.empty() ||
- glb_value > glb_val ||
- (glb_value == glb_val && r.m_type == t_lt && !glb_strict)) {
- glb_val = glb_value;
- glb_index = row_id;
- glb_strict = r.m_type == t_lt;
- }
- glb_rows.push_back(row_id);
- glb_is_unit &= a.is_minus_one();
- }
- }
-
- if (!mod_rows.empty()) {
- return solve_mod(x, mod_rows, compute_def);
- }
-
- if (eq_row != UINT_MAX) {
- return solve_for(eq_row, x, compute_def);
- }
-
- def result;
- unsigned lub_size = lub_rows.size();
- unsigned glb_size = glb_rows.size();
- unsigned row_index = (lub_size <= glb_size) ? lub_index : glb_index;
-
- // There are only upper or only lower bounds.
- if (row_index == UINT_MAX) {
- if (compute_def) {
- if (lub_index != UINT_MAX) {
- result = solve_for(lub_index, x, true);
- }
- else if (glb_index != UINT_MAX) {
- result = solve_for(glb_index, x, true);
- }
- else {
- result = def() + m_var2value[x];
- }
- SASSERT(eval(result) == eval(x));
- }
- else {
- for (unsigned row_id : lub_rows) retire_row(row_id);
- for (unsigned row_id : glb_rows) retire_row(row_id);
- }
- return result;
- }
-
- SASSERT(lub_index != UINT_MAX);
- SASSERT(glb_index != UINT_MAX);
- if (compute_def) {
- if (lub_size <= glb_size) {
- result = def(m_rows[lub_index], x);
- }
- else {
- result = def(m_rows[glb_index], x);
- }
- }
-
- // The number of matching lower and upper bounds is small.
- if ((lub_size <= 2 || glb_size <= 2) &&
- (lub_size <= 3 && glb_size <= 3) &&
- (!is_int(x) || lub_is_unit || glb_is_unit)) {
- for (unsigned i = 0; i < lub_size; ++i) {
- unsigned row_id1 = lub_rows[i];
- bool last = i + 1 == lub_size;
- rational coeff = get_coefficient(row_id1, x);
- for (unsigned row_id2 : glb_rows) {
- if (last) {
- resolve(row_id1, coeff, row_id2, x);
- }
- else {
- unsigned row_id3 = copy_row(row_id2);
- resolve(row_id1, coeff, row_id3, x);
- }
- }
- }
- for (unsigned row_id : lub_rows) retire_row(row_id);
-
- return result;
- }
-
- // General case.
- rational coeff = get_coefficient(row_index, x);
- for (unsigned row_id : lub_rows) {
- if (row_id != row_index) {
- resolve(row_index, coeff, row_id, x);
- }
- }
- for (unsigned row_id : glb_rows) {
- if (row_id != row_index) {
- resolve(row_index, coeff, row_id, x);
- }
- }
- retire_row(row_index);
- return result;
- }
-
- //
- // compute D and u.
- //
- // D = lcm(d1, d2)
- // u = eval(x) mod D
- //
- // d1 | (a1x + t1) & d2 | (a2x + t2)
- // =
- // d1 | (a1(D*x' + u) + t1) & d2 | (a2(D*x' + u) + t2)
- // =
- // d1 | (a1*u + t1) & d2 | (a2*u + t2)
- //
- // x := D*x' + u
- //
-
- model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def) {
- SASSERT(!mod_rows.empty());
- rational D(1);
- for (unsigned idx : mod_rows) {
- D = lcm(D, m_rows[idx].m_mod);
- }
- if (D.is_zero()) {
- throw default_exception("modulo 0 is not defined");
- }
- if (D.is_neg()) D = abs(D);
- TRACE("opt1", display(tout << "lcm: " << D << " x: v" << x << " tableau\n"););
- rational val_x = m_var2value[x];
- rational u = mod(val_x, D);
- SASSERT(u.is_nonneg() && u < D);
- for (unsigned idx : mod_rows) {
- replace_var(idx, x, u);
- SASSERT(invariant(idx, m_rows[idx]));
- normalize(idx);
- }
- TRACE("opt1", display(tout << "tableau after replace x under mod\n"););
- //
- // update inequalities such that u is added to t and
- // D is multiplied to coefficient of x.
- // the interpretation of the new version of x is (x-u)/D
- //
- // a*x + t <= 0
- // a*(D*x' + u) + t <= 0
- // a*D*x' + a*u + t <= 0
- //
- rational new_val = (val_x - u) / D;
- SASSERT(new_val.is_int());
- unsigned y = add_var(new_val, true);
- unsigned_vector const& row_ids = m_var2row_ids[x];
- uint_set visited;
- for (unsigned row_id : row_ids) {
- if (!visited.contains(row_id)) {
- // x |-> D*y + u
- replace_var(row_id, x, D, y, u);
- visited.insert(row_id);
- normalize(row_id);
- }
- }
- TRACE("opt1", display(tout << "tableau after replace x by y := v" << y << "\n"););
- def result = project(y, compute_def);
- if (compute_def) {
- result = (result * D) + u;
- m_var2value[x] = eval(result);
- }
- TRACE("opt1", display(tout << "tableau after project y" << y << "\n"););
-
- return result;
- }
-
- // update row with: x |-> C
- void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& C) {
- row& r = m_rows[row_id];
- SASSERT(!get_coefficient(row_id, x).is_zero());
- unsigned sz = r.m_vars.size();
- unsigned i = 0, j = 0;
- rational coeff(0);
- for (; i < sz; ++i) {
- if (r.m_vars[i].m_id == x) {
- coeff = r.m_vars[i].m_coeff;
- }
- else {
- if (i != j) {
- r.m_vars[j] = r.m_vars[i];
- }
- ++j;
- }
- }
- if (j != sz) {
- r.m_vars.shrink(j);
- }
- r.m_coeff += coeff*C;
- r.m_value += coeff*(C - m_var2value[x]);
- }
-
- // update row with: x |-> A*y + B
- void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B) {
- row& r = m_rows[row_id];
- rational coeff = get_coefficient(row_id, x);
- if (coeff.is_zero()) return;
- if (!r.m_alive) return;
- replace_var(row_id, x, B);
- r.m_vars.push_back(var(y, coeff*A));
- r.m_value += coeff*A*m_var2value[y];
- if (!r.m_vars.empty() && r.m_vars.back().m_id > y) {
- std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
- }
- m_var2row_ids[y].push_back(row_id);
- SASSERT(invariant(row_id, r));
- }
-
- // 3x + t = 0 & 7 | (c*x + s) & ax <= u
- // 3 | -t & 21 | (-ct + 3s) & a-t <= 3u
-
- model_based_opt::def model_based_opt::solve_for(unsigned row_id1, unsigned x, bool compute_def) {
- TRACE("opt", tout << "v" << x << " := " << eval(x) << "\n" << m_rows[row_id1] << "\n";);
- rational a = get_coefficient(row_id1, x), b;
- row& r1 = m_rows[row_id1];
- ineq_type ty = r1.m_type;
- SASSERT(!a.is_zero());
- SASSERT(r1.m_alive);
- if (a.is_neg()) {
- a.neg();
- r1.neg();
- }
- SASSERT(a.is_pos());
- if (ty == t_lt) {
- SASSERT(compute_def);
- r1.m_coeff -= r1.m_value;
- r1.m_type = t_le;
- r1.m_value = 0;
- }
-
- if (m_var2is_int[x] && !a.is_one()) {
- r1.m_coeff -= r1.m_value;
- r1.m_value = 0;
- vector<var> coeffs;
- mk_coeffs_without(coeffs, r1.m_vars, x);
- rational c = mod(-eval(coeffs), a);
- add_divides(coeffs, c, a);
- }
- unsigned_vector const& row_ids = m_var2row_ids[x];
- uint_set visited;
- visited.insert(row_id1);
- for (unsigned row_id2 : row_ids) {
- if (!visited.contains(row_id2)) {
- visited.insert(row_id2);
- b = get_coefficient(row_id2, x);
- if (b.is_zero())
- continue;
- row& dst = m_rows[row_id2];
- switch (dst.m_type) {
- case t_eq:
- case t_lt:
- case t_le:
- solve(row_id1, a, row_id2, x);
- break;
- case t_mod:
- // mod reduction already done.
- UNREACHABLE();
- break;
- }
- }
- }
- def result;
- if (compute_def) {
- result = def(m_rows[row_id1], x);
- m_var2value[x] = eval(result);
- TRACE("opt1", tout << "updated eval " << x << " := " << eval(x) << "\n";);
- }
- retire_row(row_id1);
- return result;
- }
-
- vector<model_based_opt::def> model_based_opt::project(unsigned num_vars, unsigned const* vars, bool compute_def) {
- vector<def> result;
- for (unsigned i = 0; i < num_vars; ++i) {
- result.push_back(project(vars[i], compute_def));
- TRACE("opt", display(tout << "After projecting: v" << vars[i] << "\n"););
- }
- return result;
- }
-
-}
-
+/*++
+Copyright (c) 2016 Microsoft Corporation
+
+Module Name:
+
+ model_based_opt.cpp
+
+Abstract:
+
+ Model-based optimization and projection for linear real, integer arithmetic.
+
+Author:
+
+ Nikolaj Bjorner (nbjorner) 2016-27-4
+
+Revision History:
+
+
+--*/
+
+#include "math/simplex/model_based_opt.h"
+#include "util/uint_set.h"
+#include "util/z3_exception.h"
+
+std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) {
+ switch (ie) {
+ case opt::t_eq: return out << " = ";
+ case opt::t_lt: return out << " < ";
+ case opt::t_le: return out << " <= ";
+ case opt::t_divides: return out << " divides ";
+ case opt::t_mod: return out << " mod ";
+ case opt::t_div: return out << " div ";
+ }
+ return out;
+}
+
+
+namespace opt {
+
+ /**
+ * Convert a row ax + coeffs + coeff = value into a definition for x
+ * x = (value - coeffs - coeff)/a
+ * as backdrop we have existing assignments to x and other variables that
+ * satisfy the equality with value, and such that value satisfies
+ * the row constraint ( = , <= , < , mod)
+ */
+ model_based_opt::def::def(row const& r, unsigned x) {
+ for (var const & v : r.m_vars) {
+ if (v.m_id != x) {
+ m_vars.push_back(v);
+ }
+ else {
+ m_div = -v.m_coeff;
+ }
+ }
+ m_coeff = r.m_coeff;
+ switch (r.m_type) {
+ case opt::t_lt:
+ m_coeff += m_div;
+ break;
+ case opt::t_le:
+ // for: ax >= t, then x := (t + a - 1) div a
+ if (m_div.is_pos()) {
+ m_coeff += m_div;
+ m_coeff -= rational::one();
+ }
+ break;
+ default:
+ break;
+ }
+ normalize();
+ SASSERT(m_div.is_pos());
+ }
+
+ model_based_opt::def model_based_opt::def::operator+(def const& other) const {
+ def result;
+ vector<var> const& vs1 = m_vars;
+ vector<var> const& vs2 = other.m_vars;
+ vector<var> & vs = result.m_vars;
+ rational c1(1), c2(1);
+ if (m_div != other.m_div) {
+ c1 = other.m_div;
+ c2 = m_div;
+ }
+ 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 == v2) {
+ vs.push_back(vs1[i]);
+ vs.back().m_coeff *= c1;
+ vs.back().m_coeff += c2 * 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 *= c1;
+ }
+ else {
+ vs.push_back(vs2[j]);
+ vs.back().m_coeff *= c2;
+ }
+ }
+ result.m_div = c1*m_div;
+ result.m_coeff = (m_coeff*c1) + (other.m_coeff*c2);
+ result.normalize();
+ return result;
+ }
+
+ model_based_opt::def model_based_opt::def::operator/(rational const& r) const {
+ def result(*this);
+ result.m_div *= r;
+ result.normalize();
+ return result;
+ }
+
+ model_based_opt::def model_based_opt::def::operator*(rational const& n) const {
+ def result(*this);
+ for (var& v : result.m_vars) {
+ v.m_coeff *= n;
+ }
+ result.m_coeff *= n;
+ result.normalize();
+ return result;
+ }
+
+ model_based_opt::def model_based_opt::def::operator+(rational const& n) const {
+ def result(*this);
+ result.m_coeff += n * result.m_div;
+ result.normalize();
+ return result;
+ }
+
+ void model_based_opt::def::normalize() {
+ if (!m_div.is_int()) {
+ rational den = denominator(m_div);
+ SASSERT(den > 1);
+ for (var& v : m_vars)
+ v.m_coeff *= den;
+ m_coeff *= den;
+ m_div *= den;
+
+ }
+ if (m_div.is_neg()) {
+ for (var& v : m_vars)
+ v.m_coeff.neg();
+ m_coeff.neg();
+ m_div.neg();
+ }
+ if (m_div.is_one())
+ return;
+ rational g(m_div);
+ if (!m_coeff.is_int())
+ return;
+ g = gcd(g, m_coeff);
+ for (var const& v : m_vars) {
+ if (!v.m_coeff.is_int())
+ return;
+ g = gcd(g, abs(v.m_coeff));
+ if (g.is_one())
+ break;
+ }
+ if (!g.is_one()) {
+ for (var& v : m_vars)
+ v.m_coeff /= g;
+ m_coeff /= g;
+ m_div /= g;
+ }
+ }
+
+ model_based_opt::model_based_opt() {
+ m_rows.push_back(row());
+ }
+
+ bool model_based_opt::invariant() {
+ for (unsigned i = 0; i < m_rows.size(); ++i) {
+ if (!invariant(i, m_rows[i])) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+#define PASSERT(_e_) { CTRACE("qe", !(_e_), display(tout, r); display(tout);); SASSERT(_e_); }
+
+ bool model_based_opt::invariant(unsigned index, row const& r) {
+ vector<var> const& vars = r.m_vars;
+ for (unsigned i = 0; i < vars.size(); ++i) {
+ // variables in each row are sorted and have non-zero coefficients
+ PASSERT(i + 1 == vars.size() || vars[i].m_id < vars[i+1].m_id);
+ PASSERT(!vars[i].m_coeff.is_zero());
+ PASSERT(index == 0 || m_var2row_ids[vars[i].m_id].contains(index));
+ }
+
+ PASSERT(r.m_value == eval(r));
+ PASSERT(r.m_type != t_eq || r.m_value.is_zero());
+ // values satisfy constraints
+ PASSERT(index == 0 || r.m_type != t_lt || r.m_value.is_neg());
+ PASSERT(index == 0 || r.m_type != t_le || !r.m_value.is_pos());
+ PASSERT(index == 0 || r.m_type != t_divides || (mod(r.m_value, r.m_mod).is_zero()));
+ return true;
+ }
+
+ // a1*x + obj
+ // a2*x + t2 <= 0
+ // a3*x + t3 <= 0
+ // a4*x + t4 <= 0
+ // a1 > 0, a2 > 0, a3 > 0, a4 < 0
+ // x <= -t2/a2
+ // x <= -t2/a3
+ // determine lub among these.
+ // then resolve lub with others
+ // e.g., -t2/a2 <= -t3/a3, then
+ // replace inequality a3*x + t3 <= 0 by -t2/a2 + t3/a3 <= 0
+ // mark a4 as invalid.
+ //
+
+ // a1 < 0, a2 < 0, a3 < 0, a4 > 0
+ // x >= t2/a2
+ // x >= t3/a3
+ // determine glb among these
+ // the resolve glb with others.
+ // e.g. t2/a2 >= t3/a3
+ // then replace a3*x + t3 by t3/a3 - t2/a2 <= 0
+ //
+ inf_eps model_based_opt::maximize() {
+ SASSERT(invariant());
+ unsigned_vector bound_trail, bound_vars;
+ TRACE("opt", display(tout << "tableau\n"););
+ while (!objective().m_vars.empty()) {
+ var v = objective().m_vars.back();
+ unsigned x = v.m_id;
+ rational const& coeff = v.m_coeff;
+ unsigned bound_row_index;
+ rational bound_coeff;
+ if (find_bound(x, bound_row_index, bound_coeff, coeff.is_pos())) {
+ SASSERT(!bound_coeff.is_zero());
+ TRACE("opt", display(tout << "update: " << v << " ", objective());
+ for (unsigned above : m_above) {
+ display(tout << "resolve: ", m_rows[above]);
+ });
+ for (unsigned above : m_above) {
+ resolve(bound_row_index, bound_coeff, above, x);
+ }
+ for (unsigned below : m_below) {
+ resolve(bound_row_index, bound_coeff, below, x);
+ }
+ // coeff*x + objective <= ub
+ // a2*x + t2 <= 0
+ // => coeff*x <= -t2*coeff/a2
+ // objective + t2*coeff/a2 <= ub
+
+ mul_add(false, m_objective_id, - coeff/bound_coeff, bound_row_index);
+ retire_row(bound_row_index);
+ bound_trail.push_back(bound_row_index);
+ bound_vars.push_back(x);
+ }
+ else {
+ TRACE("opt", display(tout << "unbound: " << v << " ", objective()););
+ update_values(bound_vars, bound_trail);
+ return inf_eps::infinity();
+ }
+ }
+
+ //
+ // update the evaluation of variables to satisfy the bound.
+ //
+
+ update_values(bound_vars, bound_trail);
+
+ rational value = objective().m_value;
+ if (objective().m_type == t_lt) {
+ return inf_eps(inf_rational(value, rational(-1)));
+ }
+ else {
+ return inf_eps(inf_rational(value));
+ }
+ }
+
+
+ void model_based_opt::update_value(unsigned x, rational const& val) {
+ rational old_val = m_var2value[x];
+ m_var2value[x] = val;
+ SASSERT(val.is_int() || !is_int(x));
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ for (unsigned row_id : row_ids) {
+ rational coeff = get_coefficient(row_id, x);
+ if (coeff.is_zero()) {
+ continue;
+ }
+ row & r = m_rows[row_id];
+ rational delta = coeff * (val - old_val);
+ r.m_value += delta;
+ SASSERT(invariant(row_id, r));
+ }
+ }
+
+
+ void model_based_opt::update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail) {
+ for (unsigned i = bound_trail.size(); i-- > 0; ) {
+ unsigned x = bound_vars[i];
+ row& r = m_rows[bound_trail[i]];
+ rational val = r.m_coeff;
+ rational old_x_val = m_var2value[x];
+ rational new_x_val;
+ rational x_coeff, eps(0);
+ vector<var> const& vars = r.m_vars;
+ for (var const& v : vars) {
+ if (x == v.m_id) {
+ x_coeff = v.m_coeff;
+ }
+ else {
+ val += m_var2value[v.m_id]*v.m_coeff;
+ }
+ }
+ SASSERT(!x_coeff.is_zero());
+ new_x_val = -val/x_coeff;
+
+ if (r.m_type == t_lt) {
+ eps = abs(old_x_val - new_x_val)/rational(2);
+ eps = std::min(rational::one(), eps);
+ SASSERT(!eps.is_zero());
+
+ //
+ // ax + t < 0
+ // <=> x < -t/a
+ // <=> x := -t/a - epsilon
+ //
+ if (x_coeff.is_pos()) {
+ new_x_val -= eps;
+ }
+ //
+ // -ax + t < 0
+ // <=> -ax < -t
+ // <=> -x < -t/a
+ // <=> x > t/a
+ // <=> x := t/a + epsilon
+ //
+ else {
+ new_x_val += eps;
+ }
+ }
+ TRACE("opt", display(tout << "v" << x
+ << " coeff_x: " << x_coeff
+ << " old_x_val: " << old_x_val
+ << " new_x_val: " << new_x_val
+ << " eps: " << eps << " ", r); );
+ m_var2value[x] = new_x_val;
+
+ r.m_value = eval(r);
+ SASSERT(invariant(bound_trail[i], r));
+ }
+
+ // update and check bounds for all other affected rows.
+ for (unsigned i = bound_trail.size(); i-- > 0; ) {
+ unsigned x = bound_vars[i];
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ for (unsigned row_id : row_ids) {
+ row & r = m_rows[row_id];
+ r.m_value = eval(r);
+ SASSERT(invariant(row_id, r));
+ }
+ }
+ SASSERT(invariant());
+ }
+
+ bool model_based_opt::find_bound(unsigned x, unsigned& bound_row_index, rational& bound_coeff, bool is_pos) {
+ bound_row_index = UINT_MAX;
+ rational lub_val;
+ rational const& x_val = m_var2value[x];
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ uint_set visited;
+ m_above.reset();
+ m_below.reset();
+ for (unsigned row_id : row_ids) {
+ SASSERT(row_id != m_objective_id);
+ if (visited.contains(row_id)) {
+ continue;
+ }
+ visited.insert(row_id);
+ row& r = m_rows[row_id];
+ if (r.m_alive) {
+ rational a = get_coefficient(row_id, x);
+ if (a.is_zero()) {
+ // skip
+ }
+ else if (a.is_pos() == is_pos || r.m_type == t_eq) {
+ rational value = x_val - (r.m_value/a);
+ if (bound_row_index == UINT_MAX) {
+ lub_val = value;
+ bound_row_index = row_id;
+ bound_coeff = a;
+ }
+ else if ((value == lub_val && r.m_type == opt::t_lt) ||
+ (is_pos && value < lub_val) ||
+
+ (!is_pos && value > lub_val)) {
+ m_above.push_back(bound_row_index);
+ lub_val = value;
+ bound_row_index = row_id;
+ bound_coeff = a;
+ }
+ else {
+ m_above.push_back(row_id);
+ }
+ }
+ else {
+ m_below.push_back(row_id);
+ }
+ }
+ }
+ return bound_row_index != UINT_MAX;
+ }
+
+ void model_based_opt::retire_row(unsigned row_id) {
+ m_rows[row_id].m_alive = false;
+ m_retired_rows.push_back(row_id);
+ }
+
+ rational model_based_opt::eval(unsigned x) const {
+ return m_var2value[x];
+ }
+
+ rational model_based_opt::eval(def const& d) const {
+ vector<var> const& vars = d.m_vars;
+ rational val = d.m_coeff;
+ for (var const& v : vars) {
+ val += v.m_coeff * eval(v.m_id);
+ }
+ val /= d.m_div;
+ return val;
+ }
+
+ rational model_based_opt::eval(row const& r) const {
+ vector<var> const& vars = r.m_vars;
+ rational val = r.m_coeff;
+ for (var const& v : vars) {
+ val += v.m_coeff * eval(v.m_id);
+ }
+ return val;
+ }
+
+ rational model_based_opt::eval(vector<var> const& coeffs) const {
+ rational val(0);
+ for (var const& v : coeffs)
+ val += v.m_coeff * eval(v.m_id);
+ return val;
+ }
+
+ rational model_based_opt::get_coefficient(unsigned row_id, unsigned var_id) const {
+ return m_rows[row_id].get_coefficient(var_id);
+ }
+
+ rational model_based_opt::row::get_coefficient(unsigned var_id) const {
+ if (m_vars.empty())
+ return rational::zero();
+ unsigned lo = 0, hi = m_vars.size();
+ while (lo < hi) {
+ unsigned mid = lo + (hi - lo)/2;
+ SASSERT(mid < hi);
+ unsigned id = m_vars[mid].m_id;
+ if (id == var_id) {
+ lo = mid;
+ break;
+ }
+ if (id < var_id)
+ lo = mid + 1;
+ else
+ hi = mid;
+ }
+ if (lo == m_vars.size())
+ return rational::zero();
+ unsigned id = m_vars[lo].m_id;
+ if (id == var_id)
+ return m_vars[lo].m_coeff;
+ else
+ return rational::zero();
+ }
+
+ model_based_opt::row& model_based_opt::row::normalize() {
+#if 0
+ if (m_type == t_divides)
+ return *this;
+ rational D(denominator(abs(m_coeff)));
+ if (D == 0)
+ D = 1;
+ for (auto const& [id, coeff] : m_vars)
+ if (coeff != 0)
+ D = lcm(D, denominator(abs(coeff)));
+ if (D == 1)
+ return *this;
+ SASSERT(D > 0);
+ for (auto & [id, coeff] : m_vars)
+ coeff *= D;
+ m_coeff *= D;
+#endif
+ return *this;
+ }
+
+ //
+ // Let
+ // row1: t1 + a1*x <= 0
+ // row2: t2 + a2*x <= 0
+ //
+ // assume a1, a2 have the same signs:
+ // (t2 + a2*x) <= (t1 + a1*x)*a2/a1
+ // <=> t2*a1/a2 - t1 <= 0
+ // <=> t2 - t1*a2/a1 <= 0
+ //
+ // assume a1 > 0, -a2 < 0:
+ // t1 + a1*x <= 0, t2 - a2*x <= 0
+ // t2/a2 <= -t1/a1
+ // t2 + t1*a2/a1 <= 0
+ // assume -a1 < 0, a2 > 0:
+ // t1 - a1*x <= 0, t2 + a2*x <= 0
+ // t1/a1 <= -t2/a2
+ // t2 + t1*a2/a1 <= 0
+ //
+ // the resolvent is the same in all cases (simpler proof should exist)
+ //
+ // assume a1 < 0, -a1 = a2:
+ // t1 <= a2*div(t2, a2)
+ //
+
+ void model_based_opt::resolve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) {
+
+ SASSERT(a1 == get_coefficient(row_src, x));
+ SASSERT(!a1.is_zero());
+ SASSERT(row_src != row_dst);
+
+ if (m_rows[row_dst].m_alive) {
+ rational a2 = get_coefficient(row_dst, x);
+ if (is_int(x)) {
+ TRACE("opt",
+ tout << a1 << " " << a2 << ": ";
+ display(tout, m_rows[row_dst]);
+ display(tout, m_rows[row_src]););
+ if (a1.is_pos() != a2.is_pos() || m_rows[row_src].m_type == opt::t_eq) {
+ mul_add(x, a1, row_src, a2, row_dst);
+ }
+ else {
+ mul(row_dst, abs(a1));
+ mul_add(false, row_dst, -abs(a2), row_src);
+ }
+ TRACE("opt", display(tout, m_rows[row_dst]););
+ normalize(row_dst);
+ }
+ else {
+ mul_add(row_dst != m_objective_id && a1.is_pos() == a2.is_pos(), row_dst, -a2/a1, row_src);
+ }
+ }
+ }
+
+ /**
+ * a1 > 0
+ * a1*x + r1 = value
+ * a2*x + r2 <= 0
+ * ------------------
+ * a1*r2 - a2*r1 <= value
+ */
+ void model_based_opt::solve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) {
+ SASSERT(a1 == get_coefficient(row_src, x));
+ SASSERT(a1.is_pos());
+ SASSERT(row_src != row_dst);
+ if (!m_rows[row_dst].m_alive) return;
+ rational a2 = get_coefficient(row_dst, x);
+ mul(row_dst, a1);
+ mul_add(false, row_dst, -a2, row_src);
+ SASSERT(get_coefficient(row_dst, x).is_zero());
+ }
+
+ // resolution for integer rows.
+ void model_based_opt::mul_add(
+ unsigned x, rational const& src_c, unsigned row_src, rational const& dst_c, unsigned row_dst) {
+ row& dst = m_rows[row_dst];
+ row const& src = m_rows[row_src];
+ SASSERT(is_int(x));
+ SASSERT(t_le == dst.m_type && t_le == src.m_type);
+ SASSERT(src_c.is_int());
+ SASSERT(dst_c.is_int());
+ SASSERT(m_var2value[x].is_int());
+
+ rational abs_src_c = abs(src_c);
+ rational abs_dst_c = abs(dst_c);
+ rational x_val = m_var2value[x];
+ rational slack = (abs_src_c - rational::one()) * (abs_dst_c - rational::one());
+ rational dst_val = dst.m_value - x_val*dst_c;
+ rational src_val = src.m_value - x_val*src_c;
+ rational distance = abs_src_c * dst_val + abs_dst_c * src_val + slack;
+ bool use_case1 = abs_src_c == abs_dst_c && src_c.is_pos() != dst_c.is_pos() && !abs_src_c.is_one() && t_le == dst.m_type && t_le == src.m_type;
+ bool use_case2 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one();
+
+ if (use_case1 && false) {
+ //
+ // x*src_c + s <= 0
+ // -x*src_c + t <= 0
+ //
+ // -src_c*div(-s, src_c) + t <= 0
+ //
+ // Example:
+ // t <= 100*x <= s
+ // Then t <= 100*div(s, 100)
+ //
+
+ if (src_c < 0)
+ std::swap(row_src, row_dst);
+
+ vector<var> src_coeffs, dst_coeffs;
+ rational src_coeff = m_rows[row_src].m_coeff;
+ rational dst_coeff = m_rows[row_dst].m_coeff;
+ for (auto const& v : m_rows[row_src].m_vars)
+ if (v.m_id != x)
+ src_coeffs.push_back(var(v.m_id, -v.m_coeff));
+ for (auto const& v : m_rows[row_dst].m_vars)
+ if (v.m_id != x)
+ dst_coeffs.push_back(v);
+ unsigned v = add_div(src_coeffs, -src_coeff, abs_src_c);
+ dst_coeffs.push_back(var(v, -abs_src_c));
+ add_constraint(dst_coeffs, dst_coeff, t_le);
+ retire_row(row_dst);
+ return;
+ }
+
+ if (use_case2) {
+ TRACE("opt", tout << "slack: " << slack << " " << src_c << " " << dst_val << " " << dst_c << " " << src_val << "\n";);
+ // dst <- abs_src_c*dst + abs_dst_c*src + slack
+ mul(row_dst, abs_src_c);
+ add(row_dst, slack);
+ mul_add(false, row_dst, abs_dst_c, row_src);
+ return;
+ }
+
+ //
+ // create finite disjunction for |b|.
+ // exists x, z in [0 .. |b|-2] . b*x + s + z = 0 && ax + t <= 0 && bx + s <= 0
+ // <=>
+ // exists x, z in [0 .. |b|-2] . b*x = -z - s && ax + t <= 0 && bx + s <= 0
+ // <=>
+ // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && bx + s <= 0
+ // <=>
+ // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z - s + s <= 0
+ // <=>
+ // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z <= 0
+ // <=>
+ // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0
+ // <=>
+ // exists x, z in [0 .. |b|-2] . b*x = -z - s && a*n_sign(b)(s + z) + |b|t <= 0
+ // <=>
+ // exists z in [0 .. |b|-2] . |b| | (z + s) && a*n_sign(b)(s + z) + |b|t <= 0
+ //
+
+ TRACE("qe", tout << "finite disjunction " << distance << " " << src_c << " " << dst_c << "\n";);
+ vector<var> coeffs;
+ if (abs_dst_c <= abs_src_c) {
+ rational z = mod(dst_val, abs_dst_c);
+ if (!z.is_zero()) z = abs_dst_c - z;
+ mk_coeffs_without(coeffs, dst.m_vars, x);
+ add_divides(coeffs, dst.m_coeff + z, abs_dst_c);
+ add(row_dst, z);
+ mul(row_dst, src_c * n_sign(dst_c));
+ mul_add(false, row_dst, abs_dst_c, row_src);
+ }
+ else {
+ // z := b - (s + bx) mod b
+ // := b - s mod b
+ // b | s + z <=> b | s + b - s mod b <=> b | s - s mod b
+ rational z = mod(src_val, abs_src_c);
+ if (!z.is_zero()) z = abs_src_c - z;
+ mk_coeffs_without(coeffs, src.m_vars, x);
+ add_divides(coeffs, src.m_coeff + z, abs_src_c);
+ mul(row_dst, abs_src_c);
+ add(row_dst, z * dst_c * n_sign(src_c));
+ mul_add(false, row_dst, dst_c * n_sign(src_c), row_src);
+ }
+ }
+
+ void model_based_opt::mk_coeffs_without(vector<var>& dst, vector<var> const& src, unsigned x) {
+ for (var const & v : src) {
+ if (v.m_id != x) dst.push_back(v);
+ }
+ }
+
+ rational model_based_opt::n_sign(rational const& b) const {
+ return rational(b.is_pos()?-1:1);
+ }
+
+ void model_based_opt::mul(unsigned dst, rational const& c) {
+ if (c.is_one())
+ return;
+ row& r = m_rows[dst];
+ for (auto & v : r.m_vars)
+ v.m_coeff *= c;
+ r.m_coeff *= c;
+ r.m_value *= c;
+ }
+
+ void model_based_opt::add(unsigned dst, rational const& c) {
+ row& r = m_rows[dst];
+ r.m_coeff += c;
+ r.m_value += c;
+ }
+
+ void model_based_opt::sub(unsigned dst, rational const& c) {
+ row& r = m_rows[dst];
+ r.m_coeff -= c;
+ r.m_value -= c;
+ }
+
+ void model_based_opt::normalize(unsigned row_id) {
+ row& r = m_rows[row_id];
+ if (r.m_vars.empty()) {
+ retire_row(row_id);
+ return;
+ }
+ if (r.m_type == t_divides) return;
+ rational g(abs(r.m_vars[0].m_coeff));
+ bool all_int = g.is_int();
+ for (unsigned i = 1; all_int && !g.is_one() && i < r.m_vars.size(); ++i) {
+ rational const& coeff = r.m_vars[i].m_coeff;
+ if (coeff.is_int()) {
+ g = gcd(g, abs(coeff));
+ }
+ else {
+ all_int = false;
+ }
+ }
+ if (all_int && !r.m_coeff.is_zero()) {
+ if (r.m_coeff.is_int()) {
+ g = gcd(g, abs(r.m_coeff));
+ }
+ else {
+ all_int = false;
+ }
+ }
+ if (all_int && !g.is_one()) {
+ SASSERT(!g.is_zero());
+ mul(row_id, rational::one()/g);
+ }
+ }
+
+ //
+ // set row1 <- row1 + c*row2
+ //
+ void model_based_opt::mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2) {
+ if (c.is_zero()) {
+ return;
+ }
+
+ m_new_vars.reset();
+ row& r1 = m_rows[row_id1];
+ row const& r2 = m_rows[row_id2];
+ unsigned i = 0, j = 0;
+ while (i < r1.m_vars.size() || j < r2.m_vars.size()) {
+ if (j == r2.m_vars.size()) {
+ m_new_vars.append(r1.m_vars.size() - i, r1.m_vars.data() + i);
+ break;
+ }
+ if (i == r1.m_vars.size()) {
+ for (; j < r2.m_vars.size(); ++j) {
+ m_new_vars.push_back(r2.m_vars[j]);
+ m_new_vars.back().m_coeff *= c;
+ if (row_id1 != m_objective_id) {
+ m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
+ }
+ }
+ break;
+ }
+
+ unsigned v1 = r1.m_vars[i].m_id;
+ unsigned v2 = r2.m_vars[j].m_id;
+ if (v1 == v2) {
+ m_new_vars.push_back(r1.m_vars[i]);
+ m_new_vars.back().m_coeff += c*r2.m_vars[j].m_coeff;
+ ++i;
+ ++j;
+ if (m_new_vars.back().m_coeff.is_zero()) {
+ m_new_vars.pop_back();
+ }
+ }
+ else if (v1 < v2) {
+ m_new_vars.push_back(r1.m_vars[i]);
+ ++i;
+ }
+ else {
+ m_new_vars.push_back(r2.m_vars[j]);
+ m_new_vars.back().m_coeff *= c;
+ if (row_id1 != m_objective_id) {
+ m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
+ }
+ ++j;
+ }
+ }
+ r1.m_coeff += c*r2.m_coeff;
+ r1.m_vars.swap(m_new_vars);
+ r1.m_value += c*r2.m_value;
+
+ if (!same_sign && r2.m_type == t_lt) {
+ r1.m_type = t_lt;
+ }
+ else if (same_sign && r1.m_type == t_lt && r2.m_type == t_lt) {
+ r1.m_type = t_le;
+ }
+ SASSERT(invariant(row_id1, r1));
+ }
+
+ void model_based_opt::display(std::ostream& out) const {
+ for (auto const& r : m_rows) {
+ display(out, r);
+ }
+ for (unsigned i = 0; i < m_var2row_ids.size(); ++i) {
+ unsigned_vector const& rows = m_var2row_ids[i];
+ out << i << ": ";
+ for (auto const& r : rows) {
+ out << r << " ";
+ }
+ out << "\n";
+ }
+ }
+
+ void model_based_opt::display(std::ostream& out, vector<var> const& vars, rational const& coeff) {
+ unsigned i = 0;
+ for (var const& v : vars) {
+ if (i > 0 && v.m_coeff.is_pos()) {
+ out << "+ ";
+ }
+ ++i;
+ if (v.m_coeff.is_one()) {
+ out << "v" << v.m_id << " ";
+ }
+ else {
+ out << v.m_coeff << "*v" << v.m_id << " ";
+ }
+ }
+ if (coeff.is_pos())
+ out << " + " << coeff << " ";
+ else if (coeff.is_neg())
+ out << coeff << " ";
+ }
+
+ std::ostream& model_based_opt::display(std::ostream& out, row const& r) {
+ out << (r.m_alive?"a":"d") << " ";
+ display(out, r.m_vars, r.m_coeff);
+ switch (r.m_type) {
+ case opt::t_divides:
+ out << r.m_type << " " << r.m_mod << " = 0; value: " << r.m_value << "\n";
+ break;
+ case opt::t_mod:
+ out << r.m_type << " " << r.m_mod << " = v" << r.m_id << " ; mod: " << mod(r.m_value, r.m_mod) << "\n";
+ break;
+ case opt::t_div:
+ out << r.m_type << " " << r.m_mod << " = v" << r.m_id << " ; div: " << div(r.m_value, r.m_mod) << "\n";
+ break;
+ default:
+ out << r.m_type << " 0; value: " << r.m_value << "\n";
+ break;
+ }
+ return out;
+ }
+
+ std::ostream& model_based_opt::display(std::ostream& out, def const& r) {
+ display(out, r.m_vars, r.m_coeff);
+ if (!r.m_div.is_one()) {
+ out << " / " << r.m_div;
+ }
+ return out;
+ }
+
+ unsigned model_based_opt::add_var(rational const& value, bool is_int) {
+ unsigned v = m_var2value.size();
+ m_var2value.push_back(value);
+ m_var2is_int.push_back(is_int);
+ SASSERT(value.is_int() || !is_int);
+ m_var2row_ids.push_back(unsigned_vector());
+ return v;
+ }
+
+ rational model_based_opt::get_value(unsigned var) {
+ return m_var2value[var];
+ }
+
+ void model_based_opt::set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel) {
+ row& r = m_rows[row_id];
+ rational val(c);
+ SASSERT(r.m_vars.empty());
+ r.m_vars.append(coeffs.size(), coeffs.data());
+ bool is_int_row = !coeffs.empty();
+ std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
+ for (auto const& c : coeffs) {
+ val += m_var2value[c.m_id] * c.m_coeff;
+ SASSERT(!is_int(c.m_id) || c.m_coeff.is_int());
+ is_int_row &= is_int(c.m_id);
+ }
+ r.m_alive = true;
+ r.m_coeff = c;
+ r.m_value = val;
+ r.m_type = rel;
+ r.m_mod = m;
+ if (is_int_row && rel == t_lt) {
+ r.m_type = t_le;
+ r.m_coeff += rational::one();
+ r.m_value += rational::one();
+ }
+ }
+
+ unsigned model_based_opt::new_row() {
+ unsigned row_id = 0;
+ if (m_retired_rows.empty()) {
+ row_id = m_rows.size();
+ m_rows.push_back(row());
+ }
+ else {
+ row_id = m_retired_rows.back();
+ m_retired_rows.pop_back();
+ m_rows[row_id].reset();
+ m_rows[row_id].m_alive = true;
+ }
+ return row_id;
+ }
+
+ unsigned model_based_opt::copy_row(unsigned src, unsigned excl) {
+ unsigned dst = new_row();
+ row const& r = m_rows[src];
+ set_row(dst, r.m_vars, r.m_coeff, r.m_mod, r.m_type);
+ for (auto const& v : r.m_vars) {
+ if (v.m_id != excl)
+ m_var2row_ids[v.m_id].push_back(dst);
+ }
+ SASSERT(invariant(dst, m_rows[dst]));
+ return dst;
+ }
+
+ 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);
+ }
+
+ void model_based_opt::add_upper_bound(unsigned x, rational const& hi) {
+ vector<var> coeffs;
+ coeffs.push_back(var(x, rational::one()));
+ add_constraint(coeffs, -hi, t_le);
+ }
+
+ void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, ineq_type rel) {
+ add_constraint(coeffs, c, rational::zero(), rel);
+ }
+
+ void model_based_opt::add_divides(vector<var> const& coeffs, rational const& c, rational const& m) {
+ add_constraint(coeffs, c, m, t_divides);
+ }
+
+ unsigned model_based_opt::add_mod(vector<var> const& coeffs, rational const& c, rational const& m) {
+ add_constraint(coeffs, c, m, t_mod);
+ unsigned row_id = m_rows.size() - 1;
+ auto& r = m_rows[row_id];
+ unsigned v = add_var(mod(r.m_value, m), true);
+ r.m_id = v;
+ m_var2row_ids[v].push_back(row_id);
+ return v;
+ }
+
+ unsigned model_based_opt::add_div(vector<var> const& coeffs, rational const& c, rational const& m) {
+ add_constraint(coeffs, c, m, t_div);
+ unsigned row_id = m_rows.size() - 1;
+ auto& r = m_rows[row_id];
+ unsigned v = add_var(div(r.m_value, m), true);
+ r.m_id = v;
+ m_var2row_ids[v].push_back(row_id);
+ return v;
+ }
+
+ void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel) {
+ unsigned row_id = new_row();
+ set_row(row_id, coeffs, c, m, rel);
+ for (var const& coeff : coeffs)
+ m_var2row_ids[coeff.m_id].push_back(row_id);
+ SASSERT(invariant(row_id, m_rows[row_id]));
+ }
+
+ void model_based_opt::set_objective(vector<var> const& coeffs, rational const& c) {
+ set_row(m_objective_id, coeffs, c, rational::zero(), t_le);
+ }
+
+ void model_based_opt::get_live_rows(vector<row>& rows) {
+ for (row & r : m_rows) {
+ if (r.m_alive) {
+ rows.push_back(r.normalize());
+ }
+ }
+ }
+
+ //
+ // pick glb and lub representative.
+ // The representative is picked such that it
+ // represents the fewest inequalities.
+ // The constraints that enforce a glb or lub are not forced.
+ // The constraints that separate the glb from ub or the lub from lb
+ // are not forced.
+ // In other words, suppose there are
+ // . N inequalities of the form t <= x
+ // . M inequalities of the form s >= x
+ // . t0 is glb among N under valuation.
+ // . s0 is lub among M under valuation.
+ // If N < M
+ // create the inequalities:
+ // t <= t0 for each t other than t0 (N-1 inequalities).
+ // t0 <= s for each s (M inequalities).
+ // If N >= M the construction is symmetric.
+ //
+ model_based_opt::def model_based_opt::project(unsigned x, bool compute_def) {
+ unsigned_vector& lub_rows = m_lub;
+ unsigned_vector& glb_rows = m_glb;
+ unsigned_vector& divide_rows = m_divides;
+ unsigned_vector& mod_rows = m_mod;
+ unsigned_vector& div_rows = m_div;
+ unsigned lub_index = UINT_MAX, glb_index = UINT_MAX;
+ bool lub_strict = false, glb_strict = false;
+ rational lub_val, glb_val;
+ rational const& x_val = m_var2value[x];
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ uint_set visited;
+ lub_rows.reset();
+ glb_rows.reset();
+ divide_rows.reset();
+ mod_rows.reset();
+ div_rows.reset();
+ bool lub_is_unit = true, glb_is_unit = true;
+ unsigned eq_row = UINT_MAX;
+ // select the lub and glb.
+ for (unsigned row_id : row_ids) {
+ if (visited.contains(row_id)) {
+ continue;
+ }
+ visited.insert(row_id);
+ row& r = m_rows[row_id];
+ if (!r.m_alive) {
+ continue;
+ }
+ rational a = get_coefficient(row_id, x);
+ if (a.is_zero()) {
+ continue;
+ }
+ if (r.m_type == t_eq) {
+ eq_row = row_id;
+ continue;
+ }
+ if (r.m_type == t_mod)
+ mod_rows.push_back(row_id);
+ else if (r.m_type == t_div)
+ div_rows.push_back(row_id);
+ else if (r.m_type == t_divides)
+ divide_rows.push_back(row_id);
+ else if (a.is_pos()) {
+ rational lub_value = x_val - (r.m_value/a);
+ if (lub_rows.empty() ||
+ lub_value < lub_val ||
+ (lub_value == lub_val && r.m_type == t_lt && !lub_strict)) {
+ lub_val = lub_value;
+ lub_index = row_id;
+ lub_strict = r.m_type == t_lt;
+ }
+ lub_rows.push_back(row_id);
+ lub_is_unit &= a.is_one();
+ }
+ else {
+ SASSERT(a.is_neg());
+ rational glb_value = x_val - (r.m_value/a);
+ if (glb_rows.empty() ||
+ glb_value > glb_val ||
+ (glb_value == glb_val && r.m_type == t_lt && !glb_strict)) {
+ glb_val = glb_value;
+ glb_index = row_id;
+ glb_strict = r.m_type == t_lt;
+ }
+ glb_rows.push_back(row_id);
+ glb_is_unit &= a.is_minus_one();
+ }
+ }
+
+ if (!mod_rows.empty())
+ return solve_mod(x, mod_rows, compute_def);
+
+ if (!div_rows.empty())
+ return solve_div(x, div_rows, compute_def);
+
+ if (!divide_rows.empty())
+ return solve_divides(x, divide_rows, compute_def);
+
+ if (eq_row != UINT_MAX)
+ return solve_for(eq_row, x, compute_def);
+
+ def result;
+ unsigned lub_size = lub_rows.size();
+ unsigned glb_size = glb_rows.size();
+ unsigned row_index = (lub_size <= glb_size) ? lub_index : glb_index;
+
+ // There are only upper or only lower bounds.
+ if (row_index == UINT_MAX) {
+ if (compute_def) {
+ if (lub_index != UINT_MAX) {
+ result = solve_for(lub_index, x, true);
+ }
+ else if (glb_index != UINT_MAX) {
+ result = solve_for(glb_index, x, true);
+ }
+ else {
+ result = def() + m_var2value[x];
+ }
+ SASSERT(eval(result) == eval(x));
+ }
+ else {
+ for (unsigned row_id : lub_rows) retire_row(row_id);
+ for (unsigned row_id : glb_rows) retire_row(row_id);
+ }
+ return result;
+ }
+
+ SASSERT(lub_index != UINT_MAX);
+ SASSERT(glb_index != UINT_MAX);
+ if (compute_def) {
+ if (lub_size <= glb_size) {
+ result = def(m_rows[lub_index], x);
+ }
+ else {
+ result = def(m_rows[glb_index], x);
+ }
+ }
+
+ // The number of matching lower and upper bounds is small.
+ if ((lub_size <= 2 || glb_size <= 2) &&
+ (lub_size <= 3 && glb_size <= 3) &&
+ (!is_int(x) || lub_is_unit || glb_is_unit)) {
+ for (unsigned i = 0; i < lub_size; ++i) {
+ unsigned row_id1 = lub_rows[i];
+ bool last = i + 1 == lub_size;
+ rational coeff = get_coefficient(row_id1, x);
+ for (unsigned row_id2 : glb_rows) {
+ if (last) {
+ resolve(row_id1, coeff, row_id2, x);
+ }
+ else {
+ unsigned row_id3 = copy_row(row_id2);
+ resolve(row_id1, coeff, row_id3, x);
+ }
+ }
+ }
+ for (unsigned row_id : lub_rows) retire_row(row_id);
+
+ return result;
+ }
+
+ // General case.
+ rational coeff = get_coefficient(row_index, x);
+
+ for (unsigned row_id : lub_rows)
+ if (row_id != row_index)
+ resolve(row_index, coeff, row_id, x);
+
+ for (unsigned row_id : glb_rows)
+ if (row_id != row_index)
+ resolve(row_index, coeff, row_id, x);
+ retire_row(row_index);
+ return result;
+ }
+
+
+ //
+ // Given v = a*x + b mod m
+ //
+ // - remove v = a*x + b mod m
+ //
+ // case a = 1:
+ // - add w = b mod m
+ // - x |-> m*y + z, 0 <= z < m
+ // - if z.value + w.value < m:
+ // add z + w - v = 0
+ // - if z.value + w.value >= m:
+ // add z + w - v - m = 0
+ //
+ // case a != 1, gcd(a, m) = 1
+ // - x |-> x*y + a^-1*z, 0 <= z < m
+ // - add w = b mod m
+ // if z.value + w.value < m
+ // add z + w - v = 0
+ // if z.value + w.value >= m
+ // add z + w - v - m = 0
+ //
+ // case a != 1, gcd(a,m) = g != 1
+ // - x |-> x*y + a^-1*z, 0 <= z < m
+ // a*x + b mod m = v is now
+ // g*z + b mod m = v
+ // - add w = b mod m
+ // - 0 <= g*z.value + w.value < m*(g+1)
+ // - add g*z + w - v - k*m = 0 for suitable k from 0 .. g based on model
+ //
+
+ model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def) {
+ def result;
+ SASSERT(!mod_rows.empty());
+ unsigned row_index = mod_rows.back();
+ rational a = get_coefficient(row_index, x);
+ rational m = m_rows[row_index].m_mod;
+ unsigned v = m_rows[row_index].m_id;
+ rational x_value = m_var2value[x];
+ // disable row
+ m_rows[row_index].m_alive = false;
+ replace_var(row_index, x, rational::zero());
+
+
+ // compute a_inv
+ rational a_inv, m_inv;
+ rational g = gcd(a, m, a_inv, m_inv);
+ if (a_inv.is_neg())
+ a_inv = mod(a_inv, m);
+ SASSERT(mod(a_inv * a, m) == g);
+
+ // solve for x_value = m*y_value + a^-1*z_value, 0 <= z_value < m.
+ rational z_value = mod(x_value, m);
+ rational y_value = div(x_value, m) - div(a_inv*z_value, m);
+ SASSERT(x_value == m*y_value + a_inv*z_value);
+ SASSERT(0 <= z_value && z_value < m);
+
+ // add new variables
+ unsigned y = add_var(y_value, true);
+ unsigned z = add_var(z_value, true);
+ // TODO: we could recycle x by either y or z.
+
+ // replace x by m*y + a^-1*z in other rows.
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ uint_set visited;
+ visited.insert(row_index);
+ for (unsigned row_id : row_ids) {
+ if (visited.contains(row_id))
+ continue;
+ replace_var(row_id, x, m, y, a_inv, z);
+ visited.insert(row_id);
+ normalize(row_id);
+ }
+
+ // add bounds for z
+ add_lower_bound(z, rational::zero());
+ add_upper_bound(z, m - 1);
+
+
+ // add w = b mod m
+ vector<var> coeffs = m_rows[row_index].m_vars;
+ rational coeff = m_rows[row_index].m_coeff;
+
+ unsigned w = add_mod(coeffs, coeff, m);
+
+ rational w_value = m_var2value[w];
+
+ // add g*z + w - v - k*m = 0, for k = (g*z_value + w_value) div m
+ rational km = div(g*z_value + w_value, m)*m;
+ vector<var> mod_coeffs;
+ mod_coeffs.push_back(var(z, g));
+ mod_coeffs.push_back(var(w, rational::one()));
+ mod_coeffs.push_back(var(v, rational::minus_one()));
+ add_constraint(mod_coeffs, km, t_eq);
+ add_lower_bound(v, rational::zero());
+ add_upper_bound(v, m - 1);
+
+ // allow to recycle row.
+ m_retired_rows.push_back(row_index);
+
+ project(v, false);
+ def y_def = project(y, compute_def);
+ def z_def = project(z, compute_def);
+
+ if (compute_def) {
+ result = (y_def * m) + z_def;
+ m_var2value[x] = eval(result);
+ }
+
+ return result;
+ }
+
+ //
+ // Given v = a*x + b div m
+ // Replace x |-> m*y + a_inv*z
+ // - w = b div m
+ // - v = ((m*y + g*z) + b) div m
+ // = a*y + (a_inv*z + b) div m
+ // = a*y + b div m + (b mod m + g*z) div m
+ // = a*y + b div m + k
+ // where k := (b.value mod m + g*z.value) div m
+ //
+ model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& div_rows, bool compute_def) {
+ def result;
+ SASSERT(!div_rows.empty());
+ unsigned row_index = div_rows.back();
+ rational a = get_coefficient(row_index, x);
+ rational m = m_rows[row_index].m_mod;
+ unsigned v = m_rows[row_index].m_id;
+ rational x_value = m_var2value[x];
+ // disable row
+ m_rows[row_index].m_alive = false;
+ replace_var(row_index, x, rational::zero());
+ rational b_value = m_rows[row_index].m_value;
+
+ // compute a_inv
+ rational a_inv, m_inv;
+ rational g = gcd(a, m, a_inv, m_inv);
+ if (a_inv.is_neg())
+ a_inv = mod(a_inv, m);
+ SASSERT(mod(a_inv * a, m) == g);
+
+ // solve for x_value = m*y_value + a^-1*z_value, 0 <= z_value < m.
+ rational z_value = mod(x_value, m);
+ rational y_value = div(x_value, m) - div(a_inv*z_value, m);
+ SASSERT(x_value == m*y_value + a_inv*z_value);
+ SASSERT(0 <= z_value && z_value < m);
+
+ // add new variables
+ unsigned y = add_var(y_value, true);
+ unsigned z = add_var(z_value, true);
+ // TODO: we could recycle x by either y or z.
+
+ // replace x by m*y + a^-1*z in other rows.
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ uint_set visited;
+ visited.insert(row_index);
+ for (unsigned row_id : row_ids) {
+ if (visited.contains(row_id))
+ continue;
+ replace_var(row_id, x, m, y, a_inv, z);
+ visited.insert(row_id);
+ normalize(row_id);
+ }
+
+ // add bounds for z
+ add_lower_bound(z, rational::zero());
+ add_upper_bound(z, m - 1);
+
+ // add w = b div m
+ vector<var> coeffs = m_rows[row_index].m_vars;
+ rational coeff = m_rows[row_index].m_coeff;
+
+ unsigned w = add_div(coeffs, coeff, m);
+ rational k = div(g*z_value + mod(b_value, m), m);
+ vector<var> div_coeffs;
+ div_coeffs.push_back(var(v, rational::minus_one()));
+ div_coeffs.push_back(var(y, a));
+ div_coeffs.push_back(var(w, rational::one()));
+ add_constraint(div_coeffs, k, t_eq);
+
+ // allow to recycle row.
+ m_retired_rows.push_back(row_index);
+ project(v, false);
+ def y_def = project(y, compute_def);
+ def z_def = project(z, compute_def);
+
+ if (compute_def) {
+ result = (y_def * m) + z_def;
+ m_var2value[x] = eval(result);
+ }
+
+ return result;
+
+ }
+
+ //
+ // compute D and u.
+ //
+ // D = lcm(d1, d2)
+ // u = eval(x) mod D
+ //
+ // d1 | (a1x + t1) & d2 | (a2x + t2)
+ // =
+ // d1 | (a1(D*x' + u) + t1) & d2 | (a2(D*x' + u) + t2)
+ // =
+ // d1 | (a1*u + t1) & d2 | (a2*u + t2)
+ //
+ // x := D*x' + u
+ //
+
+ model_based_opt::def model_based_opt::solve_divides(unsigned x, unsigned_vector const& divide_rows, bool compute_def) {
+ SASSERT(!divide_rows.empty());
+ rational D(1);
+ for (unsigned idx : divide_rows) {
+ D = lcm(D, m_rows[idx].m_mod);
+ }
+ if (D.is_zero()) {
+ throw default_exception("modulo 0 is not defined");
+ }
+ if (D.is_neg()) D = abs(D);
+ TRACE("opt1", display(tout << "lcm: " << D << " x: v" << x << " tableau\n"););
+ rational val_x = m_var2value[x];
+ rational u = mod(val_x, D);
+ SASSERT(u.is_nonneg() && u < D);
+ for (unsigned idx : divide_rows) {
+ replace_var(idx, x, u);
+ SASSERT(invariant(idx, m_rows[idx]));
+ normalize(idx);
+ }
+ TRACE("opt1", display(tout << "tableau after replace x under mod\n"););
+ //
+ // update inequalities such that u is added to t and
+ // D is multiplied to coefficient of x.
+ // the interpretation of the new version of x is (x-u)/D
+ //
+ // a*x + t <= 0
+ // a*(D*x' + u) + t <= 0
+ // a*D*x' + a*u + t <= 0
+ //
+ rational new_val = (val_x - u) / D;
+ SASSERT(new_val.is_int());
+ unsigned y = add_var(new_val, true);
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ uint_set visited;
+ for (unsigned row_id : row_ids) {
+ if (!visited.contains(row_id)) {
+ // x |-> D*y + u
+ replace_var(row_id, x, D, y, u);
+ visited.insert(row_id);
+ normalize(row_id);
+ }
+ }
+ TRACE("opt1", display(tout << "tableau after replace x by y := v" << y << "\n"););
+ def result = project(y, compute_def);
+ if (compute_def) {
+ result = (result * D) + u;
+ m_var2value[x] = eval(result);
+ }
+ TRACE("opt1", display(tout << "tableau after project y" << y << "\n"););
+
+ return result;
+ }
+
+ // update row with: x |-> C
+ void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& C) {
+ row& r = m_rows[row_id];
+ SASSERT(!get_coefficient(row_id, x).is_zero());
+ unsigned sz = r.m_vars.size();
+ unsigned i = 0, j = 0;
+ rational coeff(0);
+ for (; i < sz; ++i) {
+ if (r.m_vars[i].m_id == x) {
+ coeff = r.m_vars[i].m_coeff;
+ }
+ else {
+ if (i != j) {
+ r.m_vars[j] = r.m_vars[i];
+ }
+ ++j;
+ }
+ }
+ if (j != sz) {
+ r.m_vars.shrink(j);
+ }
+ r.m_coeff += coeff*C;
+ r.m_value += coeff*(C - m_var2value[x]);
+ }
+
+ // update row with: x |-> A*y + B
+ void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B) {
+ row& r = m_rows[row_id];
+ rational coeff = get_coefficient(row_id, x);
+ if (coeff.is_zero()) return;
+ if (!r.m_alive) return;
+ replace_var(row_id, x, B);
+ r.m_vars.push_back(var(y, coeff*A));
+ r.m_value += coeff*A*m_var2value[y];
+ if (!r.m_vars.empty() && r.m_vars.back().m_id > y)
+ std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
+ m_var2row_ids[y].push_back(row_id);
+ SASSERT(invariant(row_id, r));
+ }
+
+ // update row with: x |-> A*y + B*z
+ void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B, unsigned z) {
+ row& r = m_rows[row_id];
+ rational coeff = get_coefficient(row_id, x);
+ if (coeff.is_zero() || !r.m_alive)
+ return;
+ replace_var(row_id, x, rational::zero());
+ r.m_vars.push_back(var(y, coeff*A));
+ r.m_vars.push_back(var(z, coeff*B));
+ r.m_value += coeff*A*m_var2value[y];
+ r.m_value += coeff*B*m_var2value[z];
+ std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
+ m_var2row_ids[y].push_back(row_id);
+ m_var2row_ids[z].push_back(row_id);
+ SASSERT(invariant(row_id, r));
+ }
+
+ // 3x + t = 0 & 7 | (c*x + s) & ax <= u
+ // 3 | -t & 21 | (-ct + 3s) & a-t <= 3u
+
+ model_based_opt::def model_based_opt::solve_for(unsigned row_id1, unsigned x, bool compute_def) {
+ TRACE("opt", tout << "v" << x << " := " << eval(x) << "\n" << m_rows[row_id1] << "\n";);
+ rational a = get_coefficient(row_id1, x), b;
+ row& r1 = m_rows[row_id1];
+ ineq_type ty = r1.m_type;
+ SASSERT(!a.is_zero());
+ SASSERT(r1.m_alive);
+ if (a.is_neg()) {
+ a.neg();
+ r1.neg();
+ }
+ SASSERT(a.is_pos());
+ if (ty == t_lt) {
+ SASSERT(compute_def);
+ r1.m_coeff -= r1.m_value;
+ r1.m_type = t_le;
+ r1.m_value = 0;
+ }
+
+ if (m_var2is_int[x] && !a.is_one()) {
+ r1.m_coeff -= r1.m_value;
+ r1.m_value = 0;
+ vector<var> coeffs;
+ mk_coeffs_without(coeffs, r1.m_vars, x);
+ rational c = mod(-eval(coeffs), a);
+ add_divides(coeffs, c, a);
+ }
+ unsigned_vector const& row_ids = m_var2row_ids[x];
+ uint_set visited;
+ visited.insert(row_id1);
+ for (unsigned row_id2 : row_ids) {
+ if (!visited.contains(row_id2)) {
+ visited.insert(row_id2);
+ b = get_coefficient(row_id2, x);
+ if (b.is_zero())
+ continue;
+ row& dst = m_rows[row_id2];
+ switch (dst.m_type) {
+ case t_eq:
+ case t_lt:
+ case t_le:
+ solve(row_id1, a, row_id2, x);
+ break;
+ case t_divides:
+ // mod reduction already done.
+ UNREACHABLE();
+ break;
+ }
+ }
+ }
+ def result;
+ if (compute_def) {
+ result = def(m_rows[row_id1], x);
+ m_var2value[x] = eval(result);
+ TRACE("opt1", tout << "updated eval " << x << " := " << eval(x) << "\n";);
+ }
+ retire_row(row_id1);
+ return result;
+ }
+
+ vector<model_based_opt::def> model_based_opt::project(unsigned num_vars, unsigned const* vars, bool compute_def) {
+ vector<def> result;
+ for (unsigned i = 0; i < num_vars; ++i) {
+ result.push_back(project(vars[i], compute_def));
+ TRACE("opt", display(tout << "After projecting: v" << vars[i] << "\n"););
+ }
+ return result;
+ }
+
+}
+
diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h
index cb8c18542..9e1a67721 100644
--- a/src/math/simplex/model_based_opt.h
+++ b/src/math/simplex/model_based_opt.h
@@ -30,7 +30,9 @@ namespace opt {
t_eq,
t_lt,
t_le,
- t_mod
+ t_divides,
+ t_mod,
+ t_div
};
@@ -57,6 +59,7 @@ namespace opt {
ineq_type m_type; // inequality type
rational m_value; // value of m_vars + m_coeff under interpretation of m_var2value.
bool m_alive; // rows can be marked dead if they have been processed.
+ unsigned m_id; // variable defined by row (used for mod_t and div_t)
void reset() { m_vars.reset(); m_coeff.reset(); m_value.reset(); }
row& normalize();
@@ -85,9 +88,9 @@ namespace opt {
static const unsigned m_objective_id = 0;
vector<unsigned_vector> m_var2row_ids;
vector<rational> m_var2value;
- bool_vector m_var2is_int;
+ bool_vector m_var2is_int;
vector<var> m_new_vars;
- unsigned_vector m_lub, m_glb, m_mod;
+ unsigned_vector m_lub, m_glb, m_divides, m_mod, m_div;
unsigned_vector m_above, m_below;
unsigned_vector m_retired_rows;
@@ -122,14 +125,18 @@ namespace opt {
void sub(unsigned dst, rational const& c);
- void del_var(unsigned dst, unsigned x);
+ void set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel);
- void set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel);
+ void add_lower_bound(unsigned x, rational const& lo);
+
+ void add_upper_bound(unsigned x, rational const& hi);
void add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type r);
void replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B);
+ void replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B, unsigned z);
+
void replace_var(unsigned row_id, unsigned x, rational const& C);
void normalize(unsigned row_id);
@@ -138,7 +145,7 @@ namespace opt {
unsigned new_row();
- unsigned copy_row(unsigned row_id);
+ unsigned copy_row(unsigned row_id, unsigned excl = UINT_MAX);
rational n_sign(rational const& b) const;
@@ -150,8 +157,12 @@ namespace opt {
def solve_for(unsigned row_id, unsigned x, bool compute_def);
- def solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def);
-
+ def solve_divides(unsigned x, unsigned_vector const& divide_rows, bool compute_def);
+
+ def solve_mod(unsigned x, unsigned_vector const& divide_rows, bool compute_def);
+
+ def solve_div(unsigned x, unsigned_vector const& divide_rows, bool compute_def);
+
bool is_int(unsigned x) const { return m_var2is_int[x]; }
void retire_row(unsigned row_id);
@@ -173,6 +184,17 @@ namespace opt {
// add a divisibility constraint. The row should divide m.
void add_divides(vector<var> const& coeffs, rational const& c, rational const& m);
+
+ // add sub-expression for modulus
+ // v = add_mod(coeffs, m) means
+ // v = coeffs mod m
+ unsigned add_mod(vector<var> const& coeffs, rational const& c, rational const& m);
+
+ // add sub-expression for div
+ // v = add_div(coeffs, m) means
+ // v = coeffs div m
+ unsigned add_div(vector<var> const& coeffs, rational const& c, rational const& m);
+
// Set the objective function (linear).
void set_objective(vector<var> const& coeffs, rational const& c);
diff --git a/src/qe/mbp/mbp_arith.cpp b/src/qe/mbp/mbp_arith.cpp
index 377b62ac0..7674a1577 100644
--- a/src/qe/mbp/mbp_arith.cpp
+++ b/src/qe/mbp/mbp_arith.cpp
@@ -23,7 +23,6 @@ Revision History:
#include "ast/ast_util.h"
#include "ast/arith_decl_plugin.h"
#include "ast/ast_pp.h"
-#include "ast/rewriter/th_rewriter.h"
#include "ast/expr_functors.h"
#include "ast/rewriter/expr_safe_replace.h"
#include "math/simplex/model_based_opt.h"
@@ -32,13 +31,13 @@ Revision History:
#include "model/model_v2_pp.h"
namespace mbp {
-
+
struct arith_project_plugin::imp {
- ast_manager& m;
+ ast_manager& m;
arith_util a;
- bool m_check_purified { true }; // check that variables are properly pure
- bool m_apply_projection { false };
+ bool m_check_purified = true; // check that variables are properly pure
+ bool m_apply_projection = false;
imp(ast_manager& m) :
@@ -48,10 +47,10 @@ namespace mbp {
void insert_mul(expr* x, rational const& v, obj_map<expr, rational>& ts) {
rational w;
- if (ts.find(x, w))
- ts.insert(x, w + v);
- else
- ts.insert(x, v);
+ if (ts.find(x, w))
+ ts.insert(x, w + v);
+ else
+ ts.insert(x, v);
}
@@ -65,15 +64,15 @@ namespace mbp {
rational c(0), mul(1);
expr_ref t(m);
opt::ineq_type ty = opt::t_le;
- expr* e1, *e2;
- DEBUG_CODE(expr_ref val(m);
- eval(lit, val);
- CTRACE("qe", !m.is_true(val), tout << mk_pp(lit, m) << " := " << val << "\n";);
- SASSERT(m.limit().is_canceled() || !m.is_false(val)););
+ expr* e1, * e2;
+ DEBUG_CODE(expr_ref val(m);
+ eval(lit, val);
+ CTRACE("qe", !m.is_true(val), tout << mk_pp(lit, m) << " := " << val << "\n";);
+ SASSERT(m.limit().is_canceled() || !m.is_false(val)););
- if (!m.inc())
+ if (!m.inc())
return false;
-
+
TRACE("opt", tout << mk_pp(lit, m) << " " << a.is_lt(lit) << " " << a.is_gt(lit) << "\n";);
bool is_not = m.is_not(lit, lit);
if (is_not) {
@@ -86,37 +85,35 @@ namespace mbp {
ty = is_not ? opt::t_lt : opt::t_le;
}
else if ((a.is_lt(lit, e1, e2) || a.is_gt(lit, e2, e1))) {
- linearize(mbo, eval, mul, e1, c, fmls, ts, tids);
+ linearize(mbo, eval, mul, e1, c, fmls, ts, tids);
linearize(mbo, eval, -mul, e2, c, fmls, ts, tids);
- ty = is_not ? opt::t_le: opt::t_lt;
+ ty = is_not ? opt::t_le : opt::t_lt;
}
else if (m.is_eq(lit, e1, e2) && !is_not && is_arith(e1)) {
- linearize(mbo, eval, mul, e1, c, fmls, ts, tids);
+ linearize(mbo, eval, mul, e1, c, fmls, ts, tids);
linearize(mbo, eval, -mul, e2, c, fmls, ts, tids);
ty = opt::t_eq;
- }
+ }
else if (m.is_eq(lit, e1, e2) && is_not && is_arith(e1)) {
-
+
rational r1, r2;
- expr_ref val1 = eval(e1);
+ expr_ref val1 = eval(e1);
expr_ref val2 = eval(e2);
- //TRACE("qe", tout << mk_pp(e1, m) << " " << val1 << "\n";);
- //TRACE("qe", tout << mk_pp(e2, m) << " " << val2 << "\n";);
if (!a.is_numeral(val1, r1)) return false;
if (!a.is_numeral(val2, r2)) return false;
SASSERT(r1 != r2);
if (r1 < r2) {
std::swap(e1, e2);
- }
+ }
ty = opt::t_lt;
- linearize(mbo, eval, mul, e1, c, fmls, ts, tids);
- linearize(mbo, eval, -mul, e2, c, fmls, ts, tids);
- }
+ linearize(mbo, eval, mul, e1, c, fmls, ts, tids);
+ linearize(mbo, eval, -mul, e2, c, fmls, ts, tids);
+ }
else if (m.is_distinct(lit) && !is_not && is_arith(to_app(lit)->get_arg(0))) {
expr_ref val(m);
rational r;
app* alit = to_app(lit);
- vector<std::pair<expr*,rational> > nums;
+ vector<std::pair<expr*, rational> > nums;
for (expr* arg : *alit) {
val = eval(arg);
TRACE("qe", tout << mk_pp(arg, m) << " " << val << "\n";);
@@ -125,8 +122,8 @@ namespace mbp {
}
std::sort(nums.begin(), nums.end(), compare_second());
for (unsigned i = 0; i + 1 < nums.size(); ++i) {
- SASSERT(nums[i].second < nums[i+1].second);
- expr_ref fml(a.mk_lt(nums[i].first, nums[i+1].first), m);
+ SASSERT(nums[i].second < nums[i + 1].second);
+ expr_ref fml(a.mk_lt(nums[i].first, nums[i + 1].first), m);
if (!linearize(mbo, eval, fml, fmls, tids)) {
return false;
}
@@ -139,14 +136,14 @@ namespace mbp {
map<rational, expr*, rational::hash_proc, rational::eq_proc> values;
bool found_eq = false;
for (unsigned i = 0; !found_eq && i < to_app(lit)->get_num_args(); ++i) {
- expr* arg1 = to_app(lit)->get_arg(i), *arg2 = nullptr;
+ expr* arg1 = to_app(lit)->get_arg(i), * arg2 = nullptr;
rational r;
expr_ref val = eval(arg1);
TRACE("qe", tout << mk_pp(arg1, m) << " " << val << "\n";);
if (!a.is_numeral(val, r)) return false;
if (values.find(r, arg2)) {
ty = opt::t_eq;
- linearize(mbo, eval, mul, arg1, c, fmls, ts, tids);
+ linearize(mbo, eval, mul, arg1, c, fmls, ts, tids);
linearize(mbo, eval, -mul, arg2, c, fmls, ts, tids);
found_eq = true;
}
@@ -169,27 +166,37 @@ namespace mbp {
//
// convert linear arithmetic term into an inequality for mbo.
//
- void linearize(opt::model_based_opt& mbo, model_evaluator& eval, rational const& mul, expr* t, rational& c,
- expr_ref_vector& fmls, obj_map<expr, rational>& ts, obj_map<expr, unsigned>& tids) {
- expr* t1, *t2, *t3;
+ void linearize(opt::model_based_opt& mbo, model_evaluator& eval, rational const& mul, expr* t, rational& c,
+ expr_ref_vector& fmls, obj_map<expr, rational>& ts, obj_map<expr, unsigned>& tids) {
+ expr* t1, * t2, * t3;
rational mul1;
expr_ref val(m);
- if (a.is_mul(t, t1, t2) && is_numeral(t1, mul1))
- linearize(mbo, eval, mul* mul1, t2, c, fmls, ts, tids);
- else if (a.is_mul(t, t1, t2) && is_numeral(t2, mul1))
- linearize(mbo, eval, mul* mul1, t1, c, fmls, ts, tids);
+
+ auto add_def = [&](expr* t1, rational const& m, vars& coeffs) {
+ obj_map<expr, rational> ts0;
+ rational mul0(1), c0(0);
+ linearize(mbo, eval, mul0, t1, c0, fmls, ts0, tids);
+ extract_coefficients(mbo, eval, ts0, tids, coeffs);
+ insert_mul(t, mul, ts);
+ return c0;
+ };
+
+ if (a.is_mul(t, t1, t2) && is_numeral(t1, mul1))
+ linearize(mbo, eval, mul * mul1, t2, c, fmls, ts, tids);
+ else if (a.is_mul(t, t1, t2) && is_numeral(t2, mul1))
+ linearize(mbo, eval, mul * mul1, t1, c, fmls, ts, tids);
else if (a.is_uminus(t, t1))
linearize(mbo, eval, -mul, t1, c, fmls, ts, tids);
- else if (a.is_numeral(t, mul1))
- c += mul * mul1;
+ else if (a.is_numeral(t, mul1))
+ c += mul * mul1;
else if (a.is_add(t)) {
- for (expr* arg : *to_app(t))
- linearize(mbo, eval, mul, arg, c, fmls, ts, tids);
+ for (expr* arg : *to_app(t))
+ linearize(mbo, eval, mul, arg, c, fmls, ts, tids);
}
else if (a.is_sub(t, t1, t2)) {
- linearize(mbo, eval, mul, t1, c, fmls, ts, tids);
+ linearize(mbo, eval, mul, t1, c, fmls, ts, tids);
linearize(mbo, eval, -mul, t2, c, fmls, ts, tids);
- }
+ }
else if (m.is_ite(t, t1, t2, t3)) {
val = eval(t1);
@@ -210,21 +217,30 @@ namespace mbp {
else if (a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && !mul1.is_zero()) {
rational r;
val = eval(t);
- if (!a.is_numeral(val, r)) {
+ if (!a.is_numeral(val, r))
throw default_exception("mbp evaluation didn't produce an integer");
- }
- c += mul*r;
+ c += mul * r;
// t1 mod mul1 == r
- rational c0(-r), mul0(1);
- obj_map<expr, rational> ts0;
- linearize(mbo, eval, mul0, t1, c0, fmls, ts0, tids);
vars coeffs;
- extract_coefficients(mbo, eval, ts0, tids, coeffs);
- mbo.add_divides(coeffs, c0, mul1);
+ rational c0 = add_def(t1, mul1, coeffs);
+ mbo.add_divides(coeffs, c0 - r, mul1);
}
- else {
+ else if (false && a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && !mul1.is_zero()) {
+ // v = t1 mod mul1
+ vars coeffs;
+ rational c0 = add_def(t1, mul1, coeffs);
+ unsigned v = mbo.add_mod(coeffs, c0, mul1);
+ tids.insert(t, v);
+ }
+ else if (false && a.is_idiv(t, t1, t2) && is_numeral(t2, mul1) && mul1 > 0) {
+ // v = t1 div mul1
+ vars coeffs;
+ rational c0 = add_def(t1, mul1, coeffs);
+ unsigned v = mbo.add_div(coeffs, c0, mul1);
+ tids.insert(t, v);
+ }
+ else
insert_mul(t, mul, ts);
- }
}
bool is_numeral(expr* t, rational& r) {
@@ -232,8 +248,8 @@ namespace mbp {
}
struct compare_second {
- bool operator()(std::pair<expr*, rational> const& a,
- std::pair<expr*, rational> const& b) const {
+ bool operator()(std::pair<expr*, rational> const& a,
+ std::pair<expr*, rational> const& b) const {
return a.second < b.second;
}
};
@@ -243,14 +259,14 @@ namespace mbp {
}
rational n_sign(rational const& b) {
- return rational(b.is_pos()?-1:1);
+ return rational(b.is_pos() ? -1 : 1);
}
bool operator()(model& model, app* v, app_ref_vector& vars, expr_ref_vector& lits) {
app_ref_vector vs(m);
vs.push_back(v);
vector<def> defs;
- return project(model, vs, lits, defs, false) && vs.empty();
+ return project(model, vs, lits, defs, false) && vs.empty();
}
typedef opt::model_based_opt::var var;
@@ -267,10 +283,10 @@ namespace mbp {
bool project(model& model, app_ref_vector& vars, expr_ref_vector& fmls, vector<def>& result, bool compute_def) {
bool has_arith = false;
- for (expr* v : vars)
- has_arith |= is_arith(v);
- if (!has_arith)
- return true;
+ for (expr* v : vars)
+ has_arith |= is_arith(v);
+ if (!has_arith)
+ return true;
model_evaluator eval(model);
TRACE("qe", tout << model;);
eval.set_model_completion(true);
@@ -281,7 +297,7 @@ namespace mbp {
expr_ref_vector pinned(m);
unsigned j = 0;
TRACE("qe", tout << "vars: " << vars << "\n";
- for (expr* f : fmls) tout << mk_pp(f, m) << " := " << model(f) << "\n";);
+ for (expr* f : fmls) tout << mk_pp(f, m) << " := " << model(f) << "\n";);
for (unsigned i = 0; i < fmls.size(); ++i) {
expr* fml = fmls.get(i);
if (!linearize(mbo, eval, fml, fmls, tids)) {
@@ -294,8 +310,8 @@ namespace mbp {
}
fmls.shrink(j);
TRACE("qe", tout << "formulas\n" << fmls << "\n";
- for (auto const& [e, id] : tids)
- tout << mk_pp(e, m) << " -> " << id << "\n";);
+ for (auto const& [e, id] : tids)
+ tout << mk_pp(e, m) << " -> " << id << "\n";);
// fmls holds residue,
// mbo holds linear inequalities that are in scope
@@ -306,7 +322,7 @@ namespace mbp {
// return those to fmls.
expr_mark var_mark, fmls_mark;
- for (app * v : vars) {
+ for (app* v : vars) {
var_mark.mark(v);
if (is_arith(v) && !tids.contains(v)) {
rational r;
@@ -321,60 +337,70 @@ namespace mbp {
}
// bail on variables in non-linear sub-terms
+ auto is_pure = [&](expr* e) {
+ expr* x, * y;
+ rational r;
+ if (a.is_mod(e, x, y) && a.is_numeral(y))
+ return true;
+ if (a.is_idiv(e, x, y) && a.is_numeral(y, r) && r > 0)
+ return true;
+ return false;
+ };
+
for (auto& kv : tids) {
expr* e = kv.m_key;
- if (is_arith(e) && !var_mark.is_marked(e))
- mark_rec(fmls_mark, e);
+ if (is_arith(e) && !is_pure(e) && !var_mark.is_marked(e))
+ mark_rec(fmls_mark, e);
}
if (m_check_purified) {
- for (expr* fml : fmls)
- mark_rec(fmls_mark, fml);
+ for (expr* fml : fmls)
+ mark_rec(fmls_mark, fml);
for (auto& kv : tids) {
expr* e = kv.m_key;
- if (!var_mark.is_marked(e))
- mark_rec(fmls_mark, e);
+ if (!var_mark.is_marked(e) && !is_pure(e))
+ mark_rec(fmls_mark, e);
}
}
ptr_vector<expr> index2expr;
- for (auto& kv : tids)
- index2expr.setx(kv.m_value, kv.m_key, nullptr);
+ for (auto& kv : tids)
+ index2expr.setx(kv.m_value, kv.m_key, nullptr);
j = 0;
unsigned_vector real_vars;
for (app* v : vars) {
- if (is_arith(v) && !fmls_mark.is_marked(v))
- real_vars.push_back(tids.find(v));
- else
- vars[j++] = v;
+ if (is_arith(v) && !fmls_mark.is_marked(v))
+ real_vars.push_back(tids.find(v));
+ else
+ vars[j++] = v;
}
vars.shrink(j);
-
- TRACE("qe", tout << "remaining vars: " << vars << "\n";
- for (unsigned v : real_vars) tout << "v" << v << " " << mk_pp(index2expr[v], m) << "\n";
- mbo.display(tout););
+
+ TRACE("qe", tout << "remaining vars: " << vars << "\n";
+ for (unsigned v : real_vars) tout << "v" << v << " " << mk_pp(index2expr[v], m) << "\n";
+ mbo.display(tout););
vector<opt::model_based_opt::def> defs = mbo.project(real_vars.size(), real_vars.data(), compute_def);
vector<row> rows;
mbo.get_live_rows(rows);
rows2fmls(rows, index2expr, fmls);
TRACE("qe", mbo.display(tout << "mbo result\n");
- for (auto const& d : defs) tout << "def: " << d << "\n";
- tout << fmls << "\n";);
-
- if (compute_def)
- optdefs2mbpdef(defs, index2expr, real_vars, result);
+ for (auto const& d : defs) tout << "def: " << d << "\n";
+ tout << fmls << "\n";);
+
+ if (compute_def)
+ optdefs2mbpdef(defs, index2expr, real_vars, result);
if (m_apply_projection && !apply_projection(eval, result, fmls))
return false;
TRACE("qe",
for (auto const& [v, t] : result)
tout << v << " := " << t << "\n";
- for (auto* f : fmls)
- tout << mk_pp(f, m) << " := " << eval(f) << "\n";
- tout << "fmls:" << fmls << "\n";);
+ for (auto* f : fmls)
+ tout << mk_pp(f, m) << " := " << eval(f) << "\n";
+ tout << "fmls:" << fmls << "\n";);
return true;
- }
+ }
void optdefs2mbpdef(vector<opt::model_based_opt::def> const& defs, ptr_vector<expr> const& index2expr, unsigned_vector const& real_vars, vector<def>& result) {
SASSERT(defs.size() == real_vars.size());
@@ -384,31 +410,92 @@ namespace mbp {
bool is_int = a.is_int(x);
expr_ref_vector ts(m);
expr_ref t(m);
- for (var const& v : d.m_vars)
- ts.push_back(var2expr(index2expr, v));
+ for (var const& v : d.m_vars)
+ ts.push_back(var2expr(index2expr, v));
if (!d.m_coeff.is_zero())
ts.push_back(a.mk_numeral(d.m_coeff, is_int));
if (ts.empty())
ts.push_back(a.mk_numeral(rational(0), is_int));
t = mk_add(ts);
- if (!d.m_div.is_one() && is_int)
- t = a.mk_idiv(t, a.mk_numeral(d.m_div, is_int));
- else if (!d.m_div.is_one() && !is_int)
- t = a.mk_div(t, a.mk_numeral(d.m_div, is_int));
+ if (!d.m_div.is_one() && is_int)
+ t = a.mk_idiv(t, a.mk_numeral(d.m_div, is_int));
+ else if (!d.m_div.is_one() && !is_int)
+ t = a.mk_div(t, a.mk_numeral(d.m_div, is_int));
result.push_back(def(expr_ref(x, m), t));
}
}
- void rows2fmls(vector<row> const& rows, ptr_vector<expr> const& index2expr, expr_ref_vector& fmls) {
- for (row const& r : rows) {
- expr_ref_vector ts(m);
- expr_ref t(m), s(m), val(m);
- if (r.m_vars.empty()) {
+ expr_ref id2expr(u_map<row> const& def_vars, ptr_vector<expr> const& index2expr, unsigned id) {
+ row r;
+ if (def_vars.find(id, r))
+ return row2expr(def_vars, index2expr, r);
+ return expr_ref(index2expr[id], m);
+ }
+
+ expr_ref row2expr(u_map<row> const& def_vars, ptr_vector<expr> const& index2expr, row const& r) {
+ expr_ref_vector ts(m);
+ expr_ref t(m);
+ rational n;
+ for (var const& v : r.m_vars) {
+ t = id2expr(def_vars, index2expr, v.m_id);
+ if (a.is_numeral(t, n) && n == 0)
continue;
+ else if (a.is_numeral(t, n))
+ t = a.mk_numeral(v.m_coeff * n, a.is_int(t));
+ else if (!v.m_coeff.is_one())
+ t = a.mk_mul(a.mk_numeral(v.m_coeff, a.is_int(t)), t);
+ ts.push_back(t);
+ }
+ switch (r.m_type) {
+ case opt::t_mod:
+ if (ts.empty()) {
+ t = a.mk_int(mod(r.m_coeff, r.m_mod));
+ return t;
}
- if (r.m_vars.size() == 1 && r.m_vars[0].m_coeff.is_neg() && r.m_type != opt::t_mod) {
+ ts.push_back(a.mk_int(r.m_coeff));
+ t = mk_add(ts);
+ t = a.mk_mod(t, a.mk_int(r.m_mod));
+ return t;
+ case opt::t_div:
+ if (ts.empty()) {
+ t = a.mk_int(div(r.m_coeff, r.m_mod));
+ return t;
+ }
+ ts.push_back(a.mk_int(r.m_coeff));
+ t = mk_add(ts);
+ t = a.mk_idiv(t, a.mk_int(r.m_mod));
+ return t;
+ case opt::t_divides:
+ ts.push_back(a.mk_int(r.m_coeff));
+ return mk_add(ts);
+ default:
+ return mk_add(ts);
+ }
+ }
+
+ void rows2fmls(vector<row> const& rows, ptr_vector<expr> const& index2expr, expr_ref_vector& fmls) {
+
+ u_map<row> def_vars;
+ for (row const& r : rows) {
+ if (r.m_type == opt::t_mod)
+ def_vars.insert(r.m_id, r);
+ else if (r.m_type == opt::t_div)
+ def_vars.insert(r.m_id, r);
+ }
+
+ for (row const& r : rows) {
+ expr_ref t(m), s(m), val(m);
+
+ if (r.m_vars.empty())
+ continue;
+ if (r.m_type == opt::t_mod)
+ continue;
+ if (r.m_type == opt::t_div)
+ continue;
+
+ if (r.m_vars.size() == 1 && r.m_vars[0].m_coeff.is_neg() && r.m_type != opt::t_divides) {
var const& v = r.m_vars[0];
- t = index2expr[v.m_id];
+ t = id2expr(def_vars, index2expr, v.m_id);
if (!v.m_coeff.is_minus_one()) {
t = a.mk_mul(a.mk_numeral(-v.m_coeff, a.is_int(t)), t);
}
@@ -422,24 +509,14 @@ namespace mbp {
fmls.push_back(t);
continue;
}
- for (var const& v : r.m_vars) {
- t = index2expr[v.m_id];
- if (!v.m_coeff.is_one()) {
- t = a.mk_mul(a.mk_numeral(v.m_coeff, a.is_int(t)), t);
- }
- ts.push_back(t);
- }
- t = mk_add(ts);
+ t = row2expr(def_vars, index2expr, r);
s = a.mk_numeral(-r.m_coeff, r.m_coeff.is_int() && a.is_int(t));
switch (r.m_type) {
case opt::t_lt: t = a.mk_lt(t, s); break;
case opt::t_le: t = a.mk_le(t, s); break;
case opt::t_eq: t = a.mk_eq(t, s); break;
- case opt::t_mod:
- if (!r.m_coeff.is_zero()) {
- t = a.mk_sub(t, s);
- }
- t = a.mk_eq(a.mk_mod(t, a.mk_numeral(r.m_mod, true)), a.mk_int(0));
+ case opt::t_divides:
+ t = a.mk_eq(a.mk_mod(t, a.mk_int(r.m_mod)), a.mk_int(0));
break;
}
fmls.push_back(t);
@@ -468,11 +545,11 @@ namespace mbp {
SASSERT(validate_model(eval, fmls0));
// extract linear constraints
-
- for (expr * fml : fmls) {
+
+ for (expr* fml : fmls) {
linearize(mbo, eval, fml, fmls, tids);
}
-
+
// find optimal value
value = mbo.maximize();
@@ -543,7 +620,7 @@ namespace mbp {
}
CTRACE("qe", kv.m_value.is_zero(), tout << mk_pp(v, m) << " has coefficeint 0\n";);
if (!kv.m_value.is_zero()) {
- coeffs.push_back(var(id, kv.m_value));
+ coeffs.push_back(var(id, kv.m_value));
}
}
}
@@ -571,7 +648,7 @@ namespace mbp {
};
- arith_project_plugin::arith_project_plugin(ast_manager& m):project_plugin(m) {
+ arith_project_plugin::arith_project_plugin(ast_manager& m) :project_plugin(m) {
m_imp = alloc(imp, m);
}
From b6d71fccd86f09a0bafb6092feb16c152eef23fc Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Fri, 12 Aug 2022 10:22:22 -0400
Subject: [PATCH 02/16] fix #6265
---
src/smt/theory_lra.cpp | 21 ++++++++-------------
1 file changed, 8 insertions(+), 13 deletions(-)
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index b8001174e..7a46e1a1e 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -3507,16 +3507,15 @@ public:
case lp::lp_status::OPTIMAL: {
init_variable_values();
TRACE("arith", display(tout << st << " v" << v << " vi: " << vi << "\n"););
- inf_rational val = get_value(v);
- // inf_rational val(term_max.x, term_max.y);
+ auto val = value(v);
blocker = mk_gt(v);
- return inf_eps(rational::zero(), val);
+ return val;
}
case lp::lp_status::FEASIBLE: {
- inf_rational val = get_value(v);
+ auto val = value(v);
TRACE("arith", display(tout << st << " v" << v << " vi: " << vi << "\n"););
blocker = mk_gt(v);
- return inf_eps(rational::zero(), val);
+ return val;
}
default:
SASSERT(st == lp::lp_status::UNBOUNDED);
@@ -3533,23 +3532,19 @@ public:
rational r = val.x;
expr_ref e(m);
if (a.is_int(obj->get_sort())) {
- if (r.is_int()) {
+ if (r.is_int())
r += rational::one();
- }
- else {
+ else
r = ceil(r);
- }
e = a.mk_numeral(r, obj->get_sort());
e = a.mk_ge(obj, e);
}
else {
e = a.mk_numeral(r, obj->get_sort());
- if (val.y.is_neg()) {
+ if (val.y.is_neg())
e = a.mk_ge(obj, e);
- }
- else {
+ else
e = a.mk_gt(obj, e);
- }
}
TRACE("opt", tout << "v" << v << " " << val << " " << r << " " << e << "\n";);
return e;
From f989521a8cba520cffdea018843b3a367a88abe9 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Fri, 12 Aug 2022 11:54:10 -0400
Subject: [PATCH 03/16] add initial skeleton for xor-solver
---
src/sat/smt/xor_solver.cpp | 79 ++++++++++++++++++++++++++++++++++++++
src/sat/smt/xor_solver.h | 48 +++++++++++++++++++++++
2 files changed, 127 insertions(+)
create mode 100644 src/sat/smt/xor_solver.cpp
create mode 100644 src/sat/smt/xor_solver.h
diff --git a/src/sat/smt/xor_solver.cpp b/src/sat/smt/xor_solver.cpp
new file mode 100644
index 000000000..27354eac2
--- /dev/null
+++ b/src/sat/smt/xor_solver.cpp
@@ -0,0 +1,79 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+ xor_solver.h
+
+Abstract:
+
+ XOR solver.
+ Interface outline.
+
+--*/
+
+
+#include "sat/smt/xor_solver.h"
+#include "sat/sat_simplifier_params.hpp"
+#include "sat/sat_xor_finder.h"
+
+namespace xr {
+
+ solver::solver(euf::solver& ctx):
+ th_solver(ctx.get_manager(), symbol("xor-solver"), ctx.get_manager().get_family_id("xor-solver"))
+ {}
+
+ euf::th_solver* solver::clone(euf::solver& ctx) {
+ // and relevant copy internal state
+ return alloc(solver, ctx);
+ }
+
+ void solver::asserted(sat::literal l) {
+
+ }
+
+ bool solver::unit_propagate() {
+ return false;
+ }
+
+ void solver::get_antecedents(sat::literal l, sat::ext_justification_idx idx,
+ sat::literal_vector & r, bool probing) {
+
+ }
+
+ sat::check_result solver::check() {
+ return sat::check_result::CR_DONE;
+ }
+
+ void solver::push() {
+ }
+
+ void solver::pop(unsigned n) {
+ }
+
+ // inprocessing
+ // pre_simplify: decompile all xor constraints to allow other in-processing.
+ // simplify: recompile clauses to xor constraints
+ // literals that get added to xor constraints are tagged with the theory.
+ void solver::pre_simplify() {
+
+ }
+
+ void solver::simplify() {
+
+ }
+
+ std::ostream& solver::display(std::ostream& out) const {
+ return out;
+ }
+
+ std::ostream& solver::display_justification(std::ostream& out, sat::ext_justification_idx idx) const {
+ return out;
+ }
+
+ std::ostream& solver::display_constraint(std::ostream& out, sat::ext_constraint_idx idx) const {
+ return out;
+ }
+
+}
+
diff --git a/src/sat/smt/xor_solver.h b/src/sat/smt/xor_solver.h
new file mode 100644
index 000000000..3da30c580
--- /dev/null
+++ b/src/sat/smt/xor_solver.h
@@ -0,0 +1,48 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+ xor_solver.h
+
+Abstract:
+
+ XOR solver.
+ Interface outline.
+
+--*/
+
+#pragma once
+
+#include "sat/smt/euf_solver.h"
+
+namespace xr {
+ class solver : public euf::th_solver {
+ public:
+ solver(euf::solver& ctx);
+
+ th_solver* clone(euf::solver& ctx) override;
+
+ sat::literal internalize(expr* e, bool sign, bool root, bool redundant) override { UNREACHABLE(); return sat::null_literal; }
+
+ void internalize(expr* e, bool redundant) override { UNREACHABLE(); }
+
+
+ void asserted(sat::literal l) override;
+ bool unit_propagate() override;
+ void get_antecedents(sat::literal l, sat::ext_justification_idx idx, sat::literal_vector & r, bool probing) override;
+
+ void pre_simplify() override;
+ void simplify() override;
+
+ sat::check_result check() override;
+ void push() override;
+ void pop(unsigned n) override;
+
+ std::ostream& display(std::ostream& out) const override;
+ std::ostream& display_justification(std::ostream& out, sat::ext_justification_idx idx) const override;
+ std::ostream& display_constraint(std::ostream& out, sat::ext_constraint_idx idx) const override;
+
+ };
+
+}
From d272becade506ef4a247a26a427222c7f2826839 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Fri, 12 Aug 2022 11:54:26 -0400
Subject: [PATCH 04/16] fixes for division
---
src/math/simplex/model_based_opt.cpp | 70 +++++++++++++---------------
1 file changed, 32 insertions(+), 38 deletions(-)
diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index be2b4fe2b..577ce4633 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -1101,15 +1101,12 @@ namespace opt {
// There are only upper or only lower bounds.
if (row_index == UINT_MAX) {
if (compute_def) {
- if (lub_index != UINT_MAX) {
- result = solve_for(lub_index, x, true);
- }
- else if (glb_index != UINT_MAX) {
- result = solve_for(glb_index, x, true);
- }
- else {
- result = def() + m_var2value[x];
- }
+ if (lub_index != UINT_MAX)
+ result = solve_for(lub_index, x, true);
+ else if (glb_index != UINT_MAX)
+ result = solve_for(glb_index, x, true);
+ else
+ result = def() + m_var2value[x];
SASSERT(eval(result) == eval(x));
}
else {
@@ -1122,12 +1119,10 @@ namespace opt {
SASSERT(lub_index != UINT_MAX);
SASSERT(glb_index != UINT_MAX);
if (compute_def) {
- if (lub_size <= glb_size) {
- result = def(m_rows[lub_index], x);
- }
- else {
- result = def(m_rows[glb_index], x);
- }
+ if (lub_size <= glb_size)
+ result = def(m_rows[lub_index], x);
+ else
+ result = def(m_rows[glb_index], x);
}
// The number of matching lower and upper bounds is small.
@@ -1148,7 +1143,8 @@ namespace opt {
}
}
}
- for (unsigned row_id : lub_rows) retire_row(row_id);
+ for (unsigned row_id : lub_rows)
+ retire_row(row_id);
return result;
}
@@ -1281,13 +1277,14 @@ namespace opt {
//
// Given v = a*x + b div m
- // Replace x |-> m*y + a_inv*z
+ // Replace x |-> m*y + z
// - w = b div m
- // - v = ((m*y + g*z) + b) div m
- // = a*y + (a_inv*z + b) div m
- // = a*y + b div m + (b mod m + g*z) div m
+ // - v = ((a*m*y + a*z) + b) div m
+ // = a*y + (a*z + b) div m
+ // = a*y + b div m + (b mod m + a*z) div m
// = a*y + b div m + k
- // where k := (b.value mod m + g*z.value) div m
+ // where k := (b.value mod m + a*z.value) div m
+ // k is between 0 and a
//
model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& div_rows, bool compute_def) {
def result;
@@ -1302,32 +1299,24 @@ namespace opt {
replace_var(row_index, x, rational::zero());
rational b_value = m_rows[row_index].m_value;
- // compute a_inv
- rational a_inv, m_inv;
- rational g = gcd(a, m, a_inv, m_inv);
- if (a_inv.is_neg())
- a_inv = mod(a_inv, m);
- SASSERT(mod(a_inv * a, m) == g);
-
- // solve for x_value = m*y_value + a^-1*z_value, 0 <= z_value < m.
+ // solve for x_value = m*y_value + z_value, 0 <= z_value < m.
rational z_value = mod(x_value, m);
- rational y_value = div(x_value, m) - div(a_inv*z_value, m);
- SASSERT(x_value == m*y_value + a_inv*z_value);
+ rational y_value = div(x_value, m);
+ SASSERT(x_value == m*y_value + z_value);
SASSERT(0 <= z_value && z_value < m);
// add new variables
unsigned y = add_var(y_value, true);
unsigned z = add_var(z_value, true);
- // TODO: we could recycle x by either y or z.
- // replace x by m*y + a^-1*z in other rows.
+ // replace x by m*y + z in other rows.
unsigned_vector const& row_ids = m_var2row_ids[x];
uint_set visited;
visited.insert(row_index);
for (unsigned row_id : row_ids) {
if (visited.contains(row_id))
continue;
- replace_var(row_id, x, m, y, a_inv, z);
+ replace_var(row_id, x, m, y, rational::one(), z);
visited.insert(row_id);
normalize(row_id);
}
@@ -1339,9 +1328,16 @@ namespace opt {
// add w = b div m
vector<var> coeffs = m_rows[row_index].m_vars;
rational coeff = m_rows[row_index].m_coeff;
-
unsigned w = add_div(coeffs, coeff, m);
- rational k = div(g*z_value + mod(b_value, m), m);
+
+
+ //
+ // w = b div m
+ // v = a*y + w + k
+ // k = (a*z_value + (b_value mod m)) div m
+ //
+
+ rational k = div(a*z_value + mod(b_value, m), m);
vector<var> div_coeffs;
div_coeffs.push_back(var(v, rational::minus_one()));
div_coeffs.push_back(var(y, a));
@@ -1358,9 +1354,7 @@ namespace opt {
result = (y_def * m) + z_def;
m_var2value[x] = eval(result);
}
-
return result;
-
}
//
From 550d6914b1748c4bde23a7c661e0e0b32f32c4ba Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Fri, 12 Aug 2022 14:39:33 -0400
Subject: [PATCH 05/16] updates to div/mod handling in quantifier projection
note: the new code remains disabled at this point.
---
src/math/simplex/model_based_opt.cpp | 55 ++++++++++++++++++++++------
src/qe/mbp/mbp_arith.cpp | 8 ++--
2 files changed, 47 insertions(+), 16 deletions(-)
diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index 577ce4633..d0b4a03ea 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -1206,7 +1206,6 @@ namespace opt {
m_rows[row_index].m_alive = false;
replace_var(row_index, x, rational::zero());
-
// compute a_inv
rational a_inv, m_inv;
rational g = gcd(a, m, a_inv, m_inv);
@@ -1268,7 +1267,7 @@ namespace opt {
def z_def = project(z, compute_def);
if (compute_def) {
- result = (y_def * m) + z_def;
+ result = (y_def * m) + z_def * a_inv;
m_var2value[x] = eval(result);
}
@@ -1286,6 +1285,16 @@ namespace opt {
// where k := (b.value mod m + a*z.value) div m
// k is between 0 and a
//
+ // - k*m <= b mod m + a*z < (k+1)*m
+ //
+ // A better version using a^-1
+ // - v = (a*m*y + a^-1*a*z + b) div m
+ // = a*y + ((m*A + g)*z + b) div m where we write a*a^-1 = m*A + g
+ // = a*y + A + (g*z + b) div m
+ // - k*m <= b mod m + gz < (k+1)*m
+ // where k is between 0 and g
+ // when gcd(a, m) = 1, then there are only two cases.
+ //
model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& div_rows, bool compute_def) {
def result;
SASSERT(!div_rows.empty());
@@ -1299,24 +1308,31 @@ namespace opt {
replace_var(row_index, x, rational::zero());
rational b_value = m_rows[row_index].m_value;
- // solve for x_value = m*y_value + z_value, 0 <= z_value < m.
+ // compute a_inv
+ rational a_inv, m_inv;
+ rational g = gcd(a, m, a_inv, m_inv);
+ if (a_inv.is_neg())
+ a_inv = mod(a_inv, m);
+ SASSERT(mod(a_inv * a, m) == g);
+
+ // solve for x_value = m*y_value + a_inv*z_value, 0 <= z_value < m.
rational z_value = mod(x_value, m);
- rational y_value = div(x_value, m);
- SASSERT(x_value == m*y_value + z_value);
+ rational y_value = div(x_value, m) - div(z_value*a_inv, m);
+ SASSERT(x_value == m*y_value + a_inv*z_value);
SASSERT(0 <= z_value && z_value < m);
// add new variables
unsigned y = add_var(y_value, true);
unsigned z = add_var(z_value, true);
- // replace x by m*y + z in other rows.
+ // replace x by m*y + a^-1*z in other rows.
unsigned_vector const& row_ids = m_var2row_ids[x];
uint_set visited;
visited.insert(row_index);
for (unsigned row_id : row_ids) {
if (visited.contains(row_id))
continue;
- replace_var(row_id, x, m, y, rational::one(), z);
+ replace_var(row_id, x, m, y, a_inv, z);
visited.insert(row_id);
normalize(row_id);
}
@@ -1330,28 +1346,45 @@ namespace opt {
rational coeff = m_rows[row_index].m_coeff;
unsigned w = add_div(coeffs, coeff, m);
-
//
// w = b div m
- // v = a*y + w + k
- // k = (a*z_value + (b_value mod m)) div m
+ // v = a*y + w + (a*a_inv div m) + k
+ // k = (g*z_value + (b_value mod m)) div m
+ // k*m <= g*z + b mod m < (k+1)*m
//
rational k = div(a*z_value + mod(b_value, m), m);
+ rational n = div(a_inv * a, m);
vector<var> div_coeffs;
div_coeffs.push_back(var(v, rational::minus_one()));
div_coeffs.push_back(var(y, a));
div_coeffs.push_back(var(w, rational::one()));
+ if (n != 0) div_coeffs.push_back(var(z, n));
add_constraint(div_coeffs, k, t_eq);
+ unsigned u = add_mod(coeffs, coeff, m);
+
+ // add g*z + (b mod m) < (k + 1)*m
+ vector<var> bound_coeffs;
+ bound_coeffs.push_back(var(z, g));
+ bound_coeffs.push_back(var(u, rational::one()));
+ add_constraint(bound_coeffs, 1 - m * (k + 1), t_le);
+
+ // add k*m <= gz + (b mod m)
+ for (auto& c : bound_coeffs)
+ c.m_coeff.neg();
+ add_constraint(bound_coeffs, k * m, t_le);
+
// allow to recycle row.
m_retired_rows.push_back(row_index);
+
+ // project internal variables.
project(v, false);
def y_def = project(y, compute_def);
def z_def = project(z, compute_def);
if (compute_def) {
- result = (y_def * m) + z_def;
+ result = (y_def * m) + (z_def * a_inv);
m_var2value[x] = eval(result);
}
return result;
diff --git a/src/qe/mbp/mbp_arith.cpp b/src/qe/mbp/mbp_arith.cpp
index 7674a1577..11db720ef 100644
--- a/src/qe/mbp/mbp_arith.cpp
+++ b/src/qe/mbp/mbp_arith.cpp
@@ -225,19 +225,17 @@ namespace mbp {
rational c0 = add_def(t1, mul1, coeffs);
mbo.add_divides(coeffs, c0 - r, mul1);
}
- else if (false && a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && !mul1.is_zero()) {
+ else if (false && a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && mul1 > 0) {
// v = t1 mod mul1
vars coeffs;
rational c0 = add_def(t1, mul1, coeffs);
- unsigned v = mbo.add_mod(coeffs, c0, mul1);
- tids.insert(t, v);
+ tids.insert(t, mbo.add_mod(coeffs, c0, mul1));
}
else if (false && a.is_idiv(t, t1, t2) && is_numeral(t2, mul1) && mul1 > 0) {
// v = t1 div mul1
vars coeffs;
rational c0 = add_def(t1, mul1, coeffs);
- unsigned v = mbo.add_div(coeffs, c0, mul1);
- tids.insert(t, v);
+ tids.insert(t, mbo.add_div(coeffs, c0, mul1));
}
else
insert_mul(t, mul, ts);
From 72f4ee923016cf3931f471e020dce6c317c6f87f Mon Sep 17 00:00:00 2001
From: Bruce Mitchener <bruce.mitchener@gmail.com>
Date: Fri, 12 Aug 2022 11:47:56 +0700
Subject: [PATCH 06/16] api: Correctly map OP_BSREM0 to Z3_BSREM0.
---
src/api/api_ast.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/api/api_ast.cpp b/src/api/api_ast.cpp
index 647cc8631..df3059d4b 100644
--- a/src/api/api_ast.cpp
+++ b/src/api/api_ast.cpp
@@ -1178,7 +1178,7 @@ extern "C" {
case OP_BSMOD: return Z3_OP_BSMOD;
case OP_BSDIV0: return Z3_OP_BSDIV0;
case OP_BUDIV0: return Z3_OP_BUDIV0;
- case OP_BSREM0: return Z3_OP_BUREM0;
+ case OP_BSREM0: return Z3_OP_BSREM0;
case OP_BUREM0: return Z3_OP_BUREM0;
case OP_BSMOD0: return Z3_OP_BSMOD0;
case OP_ULEQ: return Z3_OP_ULEQ;
From 88b6c4a30d08f8a1d3d483f9c07833a5d9fa723c Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Fri, 12 Aug 2022 13:44:59 -0700
Subject: [PATCH 07/16] pdate decl collection to include functions under arrays
Signedoff-by: Nikolaj Bjorner <nbjorner@microsoft.com>
---
src/ast/decl_collector.cpp | 17 ++++++++++++++---
src/ast/decl_collector.h | 2 ++
2 files changed, 16 insertions(+), 3 deletions(-)
diff --git a/src/ast/decl_collector.cpp b/src/ast/decl_collector.cpp
index ecad96521..5b634abbd 100644
--- a/src/ast/decl_collector.cpp
+++ b/src/ast/decl_collector.cpp
@@ -48,6 +48,7 @@ bool decl_collector::is_bool(sort * s) {
}
void decl_collector::visit_func(func_decl * n) {
+ func_decl* g;
if (!m_visited.is_marked(n)) {
family_id fid = n->get_family_id();
if (fid == null_family_id)
@@ -57,6 +58,8 @@ void decl_collector::visit_func(func_decl * n) {
recfun::util u(m());
m_todo.push_back(u.get_def(n).get_rhs());
}
+ else if (m_ar_util.is_as_array(n, g))
+ m_todo.push_back(g);
m_visited.mark(n, true);
m_trail.push_back(n);
}
@@ -65,7 +68,8 @@ void decl_collector::visit_func(func_decl * n) {
decl_collector::decl_collector(ast_manager & m):
m_manager(m),
m_trail(m),
- m_dt_util(m) {
+ m_dt_util(m),
+ m_ar_util(m) {
m_basic_fid = m_manager.get_basic_family_id();
m_dt_fid = m_dt_util.get_family_id();
recfun::util rec_util(m);
@@ -156,8 +160,15 @@ void decl_collector::collect_deps(sort* s, sort_set& set) {
for (unsigned i = 0; i < num_cnstr; i++) {
func_decl * cnstr = m_dt_util.get_datatype_constructors(s)->get(i);
set.insert(cnstr->get_range());
- for (unsigned j = 0; j < cnstr->get_arity(); ++j)
- set.insert(cnstr->get_domain(j));
+ for (unsigned j = 0; j < cnstr->get_arity(); ++j) {
+ sort* n = cnstr->get_domain(j);
+ set.insert(n);
+ for (unsigned i = n->get_num_parameters(); i-- > 0; ) {
+ parameter const& p = n->get_parameter(i);
+ if (p.is_ast() && is_sort(p.get_ast()))
+ set.insert(to_sort(p.get_ast()));
+ }
+ }
}
}
diff --git a/src/ast/decl_collector.h b/src/ast/decl_collector.h
index 876b97188..31cabe796 100644
--- a/src/ast/decl_collector.h
+++ b/src/ast/decl_collector.h
@@ -23,6 +23,7 @@ Revision History:
#include "util/lim_vector.h"
#include "ast/ast.h"
#include "ast/datatype_decl_plugin.h"
+#include "ast/array_decl_plugin.h"
class decl_collector {
ast_manager & m_manager;
@@ -35,6 +36,7 @@ class decl_collector {
family_id m_basic_fid;
family_id m_dt_fid;
datatype_util m_dt_util;
+ array_util m_ar_util;
family_id m_rec_fid;
ptr_vector<ast> m_todo;
From 5669cf65bc99863f21db6fbe0d49df613c11474c Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Sat, 13 Aug 2022 06:18:13 -0700
Subject: [PATCH 08/16] bug fixes to mod/div quantifier elimination features
---
.gitignore | 3 ++-
src/math/simplex/model_based_opt.cpp | 23 ++++++++++++++---------
2 files changed, 16 insertions(+), 10 deletions(-)
diff --git a/.gitignore b/.gitignore
index 3fe3a3110..ffc50c1ba 100644
--- a/.gitignore
+++ b/.gitignore
@@ -91,4 +91,5 @@ examples/**/obj
CMakeSettings.json
# Editor temp files
*.swp
-.DS_Store
\ No newline at end of file
+.DS_Store
+dbg/**
diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index d0b4a03ea..7e654a3f6 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -871,6 +871,7 @@ namespace opt {
unsigned model_based_opt::add_var(rational const& value, bool is_int) {
unsigned v = m_var2value.size();
+ verbose_stream() << "add var " << v << "\n";
m_var2value.push_back(value);
m_var2is_int.push_back(is_int);
SASSERT(value.is_int() || !is_int);
@@ -1208,7 +1209,7 @@ namespace opt {
// compute a_inv
rational a_inv, m_inv;
- rational g = gcd(a, m, a_inv, m_inv);
+ rational g = mod(gcd(a, m, a_inv, m_inv), m);
if (a_inv.is_neg())
a_inv = mod(a_inv, m);
SASSERT(mod(a_inv * a, m) == g);
@@ -1252,7 +1253,7 @@ namespace opt {
// add g*z + w - v - k*m = 0, for k = (g*z_value + w_value) div m
rational km = div(g*z_value + w_value, m)*m;
vector<var> mod_coeffs;
- mod_coeffs.push_back(var(z, g));
+ if (g != 0) mod_coeffs.push_back(var(z, g));
mod_coeffs.push_back(var(w, rational::one()));
mod_coeffs.push_back(var(v, rational::minus_one()));
add_constraint(mod_coeffs, km, t_eq);
@@ -1270,7 +1271,7 @@ namespace opt {
result = (y_def * m) + z_def * a_inv;
m_var2value[x] = eval(result);
}
-
+ TRACE("opt", display(tout << "solve_mod\n"));
return result;
}
@@ -1308,11 +1309,14 @@ namespace opt {
replace_var(row_index, x, rational::zero());
rational b_value = m_rows[row_index].m_value;
+ TRACE("opt", display(tout << "solve_div\n"));
+
// compute a_inv
rational a_inv, m_inv;
- rational g = gcd(a, m, a_inv, m_inv);
+ rational g = mod(gcd(a, m, a_inv, m_inv), m);
if (a_inv.is_neg())
a_inv = mod(a_inv, m);
+
SASSERT(mod(a_inv * a, m) == g);
// solve for x_value = m*y_value + a_inv*z_value, 0 <= z_value < m.
@@ -1366,7 +1370,7 @@ namespace opt {
// add g*z + (b mod m) < (k + 1)*m
vector<var> bound_coeffs;
- bound_coeffs.push_back(var(z, g));
+ if (g != 0) bound_coeffs.push_back(var(z, g));
bound_coeffs.push_back(var(u, rational::one()));
add_constraint(bound_coeffs, 1 - m * (k + 1), t_le);
@@ -1387,6 +1391,7 @@ namespace opt {
result = (y_def * m) + (z_def * a_inv);
m_var2value[x] = eval(result);
}
+ TRACE("opt", display(tout << "solve_div\n"));
return result;
}
@@ -1505,13 +1510,13 @@ namespace opt {
if (coeff.is_zero() || !r.m_alive)
return;
replace_var(row_id, x, rational::zero());
- r.m_vars.push_back(var(y, coeff*A));
- r.m_vars.push_back(var(z, coeff*B));
+ if (A != 0) r.m_vars.push_back(var(y, coeff*A));
+ if (B != 0) r.m_vars.push_back(var(z, coeff*B));
r.m_value += coeff*A*m_var2value[y];
r.m_value += coeff*B*m_var2value[z];
std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
- m_var2row_ids[y].push_back(row_id);
- m_var2row_ids[z].push_back(row_id);
+ if (A != 0) m_var2row_ids[y].push_back(row_id);
+ if (B != 0) m_var2row_ids[z].push_back(row_id);
SASSERT(invariant(row_id, r));
}
From fa91a644d3eb90d9a75a0ed35e244f71182932cd Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Sat, 13 Aug 2022 07:07:14 -0700
Subject: [PATCH 09/16] make extensionality commutative
---
src/ast/array_decl_plugin.cpp | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/src/ast/array_decl_plugin.cpp b/src/ast/array_decl_plugin.cpp
index b336dd2e5..975eaa07a 100644
--- a/src/ast/array_decl_plugin.cpp
+++ b/src/ast/array_decl_plugin.cpp
@@ -326,7 +326,9 @@ func_decl * array_decl_plugin::mk_array_ext(unsigned arity, sort * const * domai
}
sort * r = to_sort(s->get_parameter(i).get_ast());
parameter param(i);
- return m_manager->mk_func_decl(m_array_ext_sym, arity, domain, r, func_decl_info(m_family_id, OP_ARRAY_EXT, 1, ¶m));
+ func_decl_info info(func_decl_info(m_family_id, OP_ARRAY_EXT, 1, ¶m));
+ info.set_commutative(true);
+ return m_manager->mk_func_decl(m_array_ext_sym, arity, domain, r, info);
}
From d80e2fb61d7d60790d652e830fa6df0622c911ac Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Sat, 13 Aug 2022 08:49:07 -0700
Subject: [PATCH 10/16] fix build
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
---
src/math/simplex/model_based_opt.cpp | 96 ++++++++++++++++++----------
src/math/simplex/model_based_opt.h | 12 +++-
2 files changed, 71 insertions(+), 37 deletions(-)
diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index 7e654a3f6..98ff63f35 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -202,6 +202,8 @@ namespace opt {
PASSERT(index == 0 || r.m_type != t_lt || r.m_value.is_neg());
PASSERT(index == 0 || r.m_type != t_le || !r.m_value.is_pos());
PASSERT(index == 0 || r.m_type != t_divides || (mod(r.m_value, r.m_mod).is_zero()));
+ PASSERT(index == 0 || r.m_type != t_mod || r.m_id < m_var2value.size());
+ PASSERT(index == 0 || r.m_type != t_div || r.m_id < m_var2value.size());
return true;
}
@@ -536,7 +538,7 @@ namespace opt {
rational a2 = get_coefficient(row_dst, x);
if (is_int(x)) {
TRACE("opt",
- tout << a1 << " " << a2 << ": ";
+ tout << x << ": " << a1 << " " << a2 << ": ";
display(tout, m_rows[row_dst]);
display(tout, m_rows[row_src]););
if (a1.is_pos() != a2.is_pos() || m_rows[row_src].m_type == opt::t_eq) {
@@ -546,7 +548,7 @@ namespace opt {
mul(row_dst, abs(a1));
mul_add(false, row_dst, -abs(a2), row_src);
}
- TRACE("opt", display(tout, m_rows[row_dst]););
+ TRACE("opt", display(tout << "result ", m_rows[row_dst]););
normalize(row_dst);
}
else {
@@ -594,7 +596,7 @@ namespace opt {
bool use_case1 = abs_src_c == abs_dst_c && src_c.is_pos() != dst_c.is_pos() && !abs_src_c.is_one() && t_le == dst.m_type && t_le == src.m_type;
bool use_case2 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one();
- if (use_case1 && false) {
+ if (use_case1) {
//
// x*src_c + s <= 0
// -x*src_c + t <= 0
@@ -620,8 +622,10 @@ namespace opt {
dst_coeffs.push_back(v);
unsigned v = add_div(src_coeffs, -src_coeff, abs_src_c);
dst_coeffs.push_back(var(v, -abs_src_c));
- add_constraint(dst_coeffs, dst_coeff, t_le);
+ if (src_c < 0)
+ std::swap(row_src, row_dst);
retire_row(row_dst);
+ add_constraint(dst_coeffs, dst_coeff, t_le);
return;
}
@@ -871,7 +875,6 @@ namespace opt {
unsigned model_based_opt::add_var(rational const& value, bool is_int) {
unsigned v = m_var2value.size();
- verbose_stream() << "add var " << v << "\n";
m_var2value.push_back(value);
m_var2is_int.push_back(is_int);
SASSERT(value.is_int() || !is_int);
@@ -946,40 +949,43 @@ namespace opt {
add_constraint(coeffs, -hi, t_le);
}
- void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, ineq_type rel) {
- add_constraint(coeffs, c, rational::zero(), rel);
+ unsigned model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, ineq_type rel) {
+ return add_constraint(coeffs, c, rational::zero(), rel, 0);
}
void model_based_opt::add_divides(vector<var> const& coeffs, rational const& c, rational const& m) {
- add_constraint(coeffs, c, m, t_divides);
+ add_constraint(coeffs, c, m, t_divides, 0);
}
unsigned model_based_opt::add_mod(vector<var> const& coeffs, rational const& c, rational const& m) {
- add_constraint(coeffs, c, m, t_mod);
- unsigned row_id = m_rows.size() - 1;
- auto& r = m_rows[row_id];
- unsigned v = add_var(mod(r.m_value, m), true);
- r.m_id = v;
- m_var2row_ids[v].push_back(row_id);
+ rational value = c;
+ for (auto const& var : coeffs)
+ value += var.m_coeff * m_var2value[var.m_id];
+ unsigned v = add_var(mod(value, m), true);
+ add_constraint(coeffs, c, m, t_mod, v);
return v;
}
- unsigned model_based_opt::add_div(vector<var> const& coeffs, rational const& c, rational const& m) {
- add_constraint(coeffs, c, m, t_div);
- unsigned row_id = m_rows.size() - 1;
- auto& r = m_rows[row_id];
- unsigned v = add_var(div(r.m_value, m), true);
- r.m_id = v;
- m_var2row_ids[v].push_back(row_id);
+ unsigned model_based_opt::add_div(vector<var> const& coeffs, rational const& c, rational const& m) {
+ rational value = c;
+ for (auto const& var : coeffs)
+ value += var.m_coeff * m_var2value[var.m_id];
+ unsigned v = add_var(div(value, m), true);
+ add_constraint(coeffs, c, m, t_div, v);
return v;
}
- void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel) {
+ unsigned model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel, unsigned id) {
+ auto const& r = m_rows.back();
+ if (r.m_vars == coeffs && r.m_coeff == c && r.m_mod == m && r.m_type == rel && r.m_id == id && r.m_alive)
+ return m_rows.size() - 1;
unsigned row_id = new_row();
set_row(row_id, coeffs, c, m, rel);
+ m_rows[row_id].m_id = id;
for (var const& coeff : coeffs)
m_var2row_ids[coeff.m_id].push_back(row_id);
SASSERT(invariant(row_id, m_rows[row_id]));
+ return row_id;
}
void model_based_opt::set_objective(vector<var> const& coeffs, rational const& c) {
@@ -1246,15 +1252,20 @@ namespace opt {
vector<var> coeffs = m_rows[row_index].m_vars;
rational coeff = m_rows[row_index].m_coeff;
- unsigned w = add_mod(coeffs, coeff, m);
+ unsigned w = UINT_MAX;
+ rational offset(0);
+ if (coeffs.empty())
+ offset = mod(coeff, m);
+ else
+ w = add_mod(coeffs, coeff, m);
- rational w_value = m_var2value[w];
+ rational w_value = w == UINT_MAX ? offset : m_var2value[w];
// add g*z + w - v - k*m = 0, for k = (g*z_value + w_value) div m
- rational km = div(g*z_value + w_value, m)*m;
+ rational km = div(g*z_value + w_value, m)*m + offset;
vector<var> mod_coeffs;
if (g != 0) mod_coeffs.push_back(var(z, g));
- mod_coeffs.push_back(var(w, rational::one()));
+ if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one()));
mod_coeffs.push_back(var(v, rational::minus_one()));
add_constraint(mod_coeffs, km, t_eq);
add_lower_bound(v, rational::zero());
@@ -1304,6 +1315,9 @@ namespace opt {
rational m = m_rows[row_index].m_mod;
unsigned v = m_rows[row_index].m_id;
rational x_value = m_var2value[x];
+
+ // display(verbose_stream(), m_rows[row_index]);
+
// disable row
m_rows[row_index].m_alive = false;
replace_var(row_index, x, rational::zero());
@@ -1348,7 +1362,12 @@ namespace opt {
// add w = b div m
vector<var> coeffs = m_rows[row_index].m_vars;
rational coeff = m_rows[row_index].m_coeff;
- unsigned w = add_div(coeffs, coeff, m);
+ unsigned w = UINT_MAX;
+ rational offset(0);
+ if (coeffs.empty())
+ offset = div(coeff, m);
+ else
+ w = add_div(coeffs, coeff, m);
//
// w = b div m
@@ -1357,28 +1376,32 @@ namespace opt {
// k*m <= g*z + b mod m < (k+1)*m
//
- rational k = div(a*z_value + mod(b_value, m), m);
+ rational k = div(a*z_value + mod(b_value, m), m) + offset;
rational n = div(a_inv * a, m);
vector<var> div_coeffs;
div_coeffs.push_back(var(v, rational::minus_one()));
div_coeffs.push_back(var(y, a));
- div_coeffs.push_back(var(w, rational::one()));
+ if (w != UINT_MAX) div_coeffs.push_back(var(w, rational::one()));
if (n != 0) div_coeffs.push_back(var(z, n));
add_constraint(div_coeffs, k, t_eq);
- unsigned u = add_mod(coeffs, coeff, m);
+ unsigned u = UINT_MAX;
+ offset = 0;
+ if (coeffs.empty())
+ offset = mod(coeff, m);
+ else
+ u = add_mod(coeffs, coeff, m);
// add g*z + (b mod m) < (k + 1)*m
vector<var> bound_coeffs;
if (g != 0) bound_coeffs.push_back(var(z, g));
- bound_coeffs.push_back(var(u, rational::one()));
- add_constraint(bound_coeffs, 1 - m * (k + 1), t_le);
+ if (u != UINT_MAX) bound_coeffs.push_back(var(u, rational::one()));
+ add_constraint(bound_coeffs, 1 - m * (k + 1) + offset, t_le);
// add k*m <= gz + (b mod m)
for (auto& c : bound_coeffs)
c.m_coeff.neg();
- add_constraint(bound_coeffs, k * m, t_le);
-
+ add_constraint(bound_coeffs, k * m - offset, t_le);
// allow to recycle row.
m_retired_rows.push_back(row_index);
@@ -1524,7 +1547,8 @@ namespace opt {
// 3 | -t & 21 | (-ct + 3s) & a-t <= 3u
model_based_opt::def model_based_opt::solve_for(unsigned row_id1, unsigned x, bool compute_def) {
- TRACE("opt", tout << "v" << x << " := " << eval(x) << "\n" << m_rows[row_id1] << "\n";);
+ TRACE("opt", tout << "v" << x << " := " << eval(x) << "\n" << m_rows[row_id1] << "\n";
+ display(tout));
rational a = get_coefficient(row_id1, x), b;
row& r1 = m_rows[row_id1];
ineq_type ty = r1.m_type;
@@ -1548,6 +1572,7 @@ namespace opt {
vector<var> coeffs;
mk_coeffs_without(coeffs, r1.m_vars, x);
rational c = mod(-eval(coeffs), a);
+ TRACE("opt", tout << "add divides " << eval(coeffs) << " " << a << " " << c << "\n");
add_divides(coeffs, c, a);
}
unsigned_vector const& row_ids = m_var2row_ids[x];
@@ -1580,6 +1605,7 @@ namespace opt {
TRACE("opt1", tout << "updated eval " << x << " := " << eval(x) << "\n";);
}
retire_row(row_id1);
+ TRACE("opt", display(tout << "solved v" << x << "\n"));
return result;
}
diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h
index 9e1a67721..40abf1333 100644
--- a/src/math/simplex/model_based_opt.h
+++ b/src/math/simplex/model_based_opt.h
@@ -50,6 +50,14 @@ namespace opt {
return x.m_id < y.m_id;
}
};
+
+ bool operator==(var const& other) const {
+ return m_id == other.m_id && m_coeff == other.m_coeff;
+ }
+
+ bool operator!=(var const& other) const {
+ return !(*this == other);
+ }
};
struct row {
row(): m_type(t_le), m_value(0), m_alive(false) {}
@@ -131,7 +139,7 @@ namespace opt {
void add_upper_bound(unsigned x, rational const& hi);
- void add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type r);
+ unsigned add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type r, unsigned id);
void replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B);
@@ -179,7 +187,7 @@ namespace opt {
// add a constraint. We assume that the constraint is
// satisfied under the values provided to the variables.
- void add_constraint(vector<var> const& coeffs, rational const& c, ineq_type r);
+ unsigned add_constraint(vector<var> const& coeffs, rational const& c, ineq_type r);
// add a divisibility constraint. The row should divide m.
void add_divides(vector<var> const& coeffs, rational const& c, rational const& m);
From f014e30d46a80fa3c58ae33dd97cbaf57b957fb2 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Sat, 13 Aug 2022 08:53:19 -0700
Subject: [PATCH 11/16] disable case1
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
---
src/math/simplex/model_based_opt.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index 98ff63f35..46ab55167 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -596,7 +596,7 @@ namespace opt {
bool use_case1 = abs_src_c == abs_dst_c && src_c.is_pos() != dst_c.is_pos() && !abs_src_c.is_one() && t_le == dst.m_type && t_le == src.m_type;
bool use_case2 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one();
- if (use_case1) {
+ if (use_case1 && false) {
//
// x*src_c + s <= 0
// -x*src_c + t <= 0
From 1d87592b130e6c2a11187edce2fc3590beb7e731 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Sun, 14 Aug 2022 11:34:03 -0700
Subject: [PATCH 12/16] fixes to mod/div elimination
elimination of mod/div should be applied to all occurrences of x under mod/div at the same time. It affects performance and termination to perform elimination on each occurrence since substituting in two new variables for eliminated x doubles the number of variables under other occurrences.
Also generalize inequality resolution to use div.
The new features are still disabled.
---
src/math/simplex/model_based_opt.cpp | 517 ++++++++++++++-------------
src/math/simplex/model_based_opt.h | 6 +-
src/qe/qsat.cpp | 50 ++-
3 files changed, 299 insertions(+), 274 deletions(-)
diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index 46ab55167..55f2616b8 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -485,7 +485,7 @@ namespace opt {
model_based_opt::row& model_based_opt::row::normalize() {
#if 0
- if (m_type == t_divides)
+ if (m_type == t_divides || m_type == t_mod || m_type == t_div)
return *this;
rational D(denominator(abs(m_coeff)));
if (D == 0)
@@ -572,12 +572,13 @@ namespace opt {
rational a2 = get_coefficient(row_dst, x);
mul(row_dst, a1);
mul_add(false, row_dst, -a2, row_src);
+ normalize(row_dst);
SASSERT(get_coefficient(row_dst, x).is_zero());
}
// resolution for integer rows.
void model_based_opt::mul_add(
- unsigned x, rational const& src_c, unsigned row_src, rational const& dst_c, unsigned row_dst) {
+ unsigned x, rational src_c, unsigned row_src, rational dst_c, unsigned row_dst) {
row& dst = m_rows[row_dst];
row const& src = m_rows[row_src];
SASSERT(is_int(x));
@@ -593,11 +594,22 @@ namespace opt {
rational dst_val = dst.m_value - x_val*dst_c;
rational src_val = src.m_value - x_val*src_c;
rational distance = abs_src_c * dst_val + abs_dst_c * src_val + slack;
- bool use_case1 = abs_src_c == abs_dst_c && src_c.is_pos() != dst_c.is_pos() && !abs_src_c.is_one() && t_le == dst.m_type && t_le == src.m_type;
- bool use_case2 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one();
+ bool use_case1 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one();
+ bool use_case2 = false && abs_src_c == abs_dst_c && src_c.is_pos() != dst_c.is_pos() && !abs_src_c.is_one() && t_le == dst.m_type && t_le == src.m_type;
+ bool use_case3 = false && src_c.is_pos() != dst_c.is_pos() && t_le == dst.m_type && t_le == src.m_type;
- if (use_case1 && false) {
- //
+
+ if (use_case1) {
+ TRACE("opt", tout << "slack: " << slack << " " << src_c << " " << dst_val << " " << dst_c << " " << src_val << "\n";);
+ // dst <- abs_src_c*dst + abs_dst_c*src + slack
+ mul(row_dst, abs_src_c);
+ add(row_dst, slack);
+ mul_add(false, row_dst, abs_dst_c, row_src);
+ return;
+ }
+
+ if (use_case2 || use_case3) {
+ // case2:
// x*src_c + s <= 0
// -x*src_c + t <= 0
//
@@ -607,10 +619,21 @@ namespace opt {
// t <= 100*x <= s
// Then t <= 100*div(s, 100)
//
-
- if (src_c < 0)
- std::swap(row_src, row_dst);
+ // case3:
+ // x*src_c + s <= 0
+ // -x*dst_c + t <= 0
+ // t <= x*dst_c, x*src_c <= -s ->
+ // t <= dst_c*div(-s, src_c) ->
+ // -dst_c*div(-s,src_c) + t <= 0
+ //
+ bool swapped = false;
+ if (src_c < 0) {
+ std::swap(row_src, row_dst);
+ std::swap(src_c, dst_c);
+ std::swap(abs_src_c, abs_dst_c);
+ swapped = true;
+ }
vector<var> src_coeffs, dst_coeffs;
rational src_coeff = m_rows[row_src].m_coeff;
rational dst_coeff = m_rows[row_dst].m_coeff;
@@ -620,23 +643,18 @@ namespace opt {
for (auto const& v : m_rows[row_dst].m_vars)
if (v.m_id != x)
dst_coeffs.push_back(v);
- unsigned v = add_div(src_coeffs, -src_coeff, abs_src_c);
- dst_coeffs.push_back(var(v, -abs_src_c));
- if (src_c < 0)
+ unsigned v = UINT_MAX;
+ if (src_coeffs.empty())
+ dst_coeff -= abs_dst_c*div(-src_coeff, abs_src_c);
+ else
+ v = add_div(src_coeffs, -src_coeff, abs_src_c);
+ if (v != UINT_MAX) dst_coeffs.push_back(var(v, -abs_dst_c));
+ if (swapped)
std::swap(row_src, row_dst);
retire_row(row_dst);
add_constraint(dst_coeffs, dst_coeff, t_le);
return;
- }
-
- if (use_case2) {
- TRACE("opt", tout << "slack: " << slack << " " << src_c << " " << dst_val << " " << dst_c << " " << src_val << "\n";);
- // dst <- abs_src_c*dst + abs_dst_c*src + slack
- mul(row_dst, abs_src_c);
- add(row_dst, slack);
- mul_add(false, row_dst, abs_dst_c, row_src);
- return;
- }
+ }
//
// create finite disjunction for |b|.
@@ -698,8 +716,10 @@ namespace opt {
row& r = m_rows[dst];
for (auto & v : r.m_vars)
v.m_coeff *= c;
+ r.m_mod *= c;
r.m_coeff *= c;
- r.m_value *= c;
+ if (r.m_type != t_div && r.m_type != t_mod)
+ r.m_value *= c;
}
void model_based_opt::add(unsigned dst, rational const& c) {
@@ -720,7 +740,12 @@ namespace opt {
retire_row(row_id);
return;
}
- if (r.m_type == t_divides) return;
+ if (r.m_type == t_divides)
+ return;
+ if (r.m_type == t_mod)
+ return;
+ if (r.m_type == t_div)
+ return;
rational g(abs(r.m_vars[0].m_coeff));
bool all_int = g.is_int();
for (unsigned i = 1; all_int && !g.is_one() && i < r.m_vars.size(); ++i) {
@@ -750,9 +775,9 @@ namespace opt {
// set row1 <- row1 + c*row2
//
void model_based_opt::mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2) {
- if (c.is_zero()) {
+ if (c.is_zero())
return;
- }
+
m_new_vars.reset();
row& r1 = m_rows[row_id1];
@@ -767,9 +792,8 @@ namespace opt {
for (; j < r2.m_vars.size(); ++j) {
m_new_vars.push_back(r2.m_vars[j]);
m_new_vars.back().m_coeff *= c;
- if (row_id1 != m_objective_id) {
- m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
- }
+ if (row_id1 != m_objective_id)
+ m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
}
break;
}
@@ -781,9 +805,8 @@ namespace opt {
m_new_vars.back().m_coeff += c*r2.m_vars[j].m_coeff;
++i;
++j;
- if (m_new_vars.back().m_coeff.is_zero()) {
- m_new_vars.pop_back();
- }
+ if (m_new_vars.back().m_coeff.is_zero())
+ m_new_vars.pop_back();
}
else if (v1 < v2) {
m_new_vars.push_back(r1.m_vars[i]);
@@ -792,9 +815,8 @@ namespace opt {
else {
m_new_vars.push_back(r2.m_vars[j]);
m_new_vars.back().m_coeff *= c;
- if (row_id1 != m_objective_id) {
- m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
- }
+ if (row_id1 != m_objective_id)
+ m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1);
++j;
}
}
@@ -802,25 +824,21 @@ namespace opt {
r1.m_vars.swap(m_new_vars);
r1.m_value += c*r2.m_value;
- if (!same_sign && r2.m_type == t_lt) {
- r1.m_type = t_lt;
- }
- else if (same_sign && r1.m_type == t_lt && r2.m_type == t_lt) {
- r1.m_type = t_le;
- }
+ if (!same_sign && r2.m_type == t_lt)
+ r1.m_type = t_lt;
+ else if (same_sign && r1.m_type == t_lt && r2.m_type == t_lt)
+ r1.m_type = t_le;
SASSERT(invariant(row_id1, r1));
}
void model_based_opt::display(std::ostream& out) const {
- for (auto const& r : m_rows) {
- display(out, r);
- }
+ for (auto const& r : m_rows)
+ display(out, r);
for (unsigned i = 0; i < m_var2row_ids.size(); ++i) {
unsigned_vector const& rows = m_var2row_ids[i];
out << i << ": ";
- for (auto const& r : rows) {
- out << r << " ";
- }
+ for (auto const& r : rows)
+ out << r << " ";
out << "\n";
}
}
@@ -828,16 +846,13 @@ namespace opt {
void model_based_opt::display(std::ostream& out, vector<var> const& vars, rational const& coeff) {
unsigned i = 0;
for (var const& v : vars) {
- if (i > 0 && v.m_coeff.is_pos()) {
- out << "+ ";
- }
+ if (i > 0 && v.m_coeff.is_pos())
+ out << "+ ";
++i;
- if (v.m_coeff.is_one()) {
- out << "v" << v.m_id << " ";
- }
- else {
- out << v.m_coeff << "*v" << v.m_id << " ";
- }
+ if (v.m_coeff.is_one())
+ out << "v" << v.m_id << " ";
+ else
+ out << v.m_coeff << "*v" << v.m_id << " ";
}
if (coeff.is_pos())
out << " + " << coeff << " ";
@@ -949,11 +964,16 @@ namespace opt {
add_constraint(coeffs, -hi, t_le);
}
- unsigned model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, ineq_type rel) {
- return add_constraint(coeffs, c, rational::zero(), rel, 0);
+ void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, ineq_type rel) {
+ add_constraint(coeffs, c, rational::zero(), rel, 0);
}
void model_based_opt::add_divides(vector<var> const& coeffs, rational const& c, rational const& m) {
+ rational g(c);
+ for (auto const& [v, coeff] : coeffs)
+ g = gcd(coeff, g);
+ if ((g/m).is_int())
+ return;
add_constraint(coeffs, c, m, t_divides, 0);
}
@@ -975,17 +995,17 @@ namespace opt {
return v;
}
- unsigned model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel, unsigned id) {
+ void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel, unsigned id) {
auto const& r = m_rows.back();
if (r.m_vars == coeffs && r.m_coeff == c && r.m_mod == m && r.m_type == rel && r.m_id == id && r.m_alive)
- return m_rows.size() - 1;
+ return;
unsigned row_id = new_row();
set_row(row_id, coeffs, c, m, rel);
m_rows[row_id].m_id = id;
for (var const& coeff : coeffs)
m_var2row_ids[coeff.m_id].push_back(row_id);
SASSERT(invariant(row_id, m_rows[row_id]));
- return row_id;
+ normalize(row_id);
}
void model_based_opt::set_objective(vector<var> const& coeffs, rational const& c) {
@@ -993,11 +1013,9 @@ namespace opt {
}
void model_based_opt::get_live_rows(vector<row>& rows) {
- for (row & r : m_rows) {
- if (r.m_alive) {
- rows.push_back(r.normalize());
- }
- }
+ for (row & r : m_rows)
+ if (r.m_alive)
+ rows.push_back(r.normalize());
}
//
@@ -1172,114 +1190,108 @@ namespace opt {
//
- // Given v = a*x + b mod m
+ // Given v = a*x + b mod K
//
- // - remove v = a*x + b mod m
+ // - remove v = a*x + b mod K
//
// case a = 1:
- // - add w = b mod m
- // - x |-> m*y + z, 0 <= z < m
- // - if z.value + w.value < m:
+ // - add w = b mod K
+ // - x |-> K*y + z, 0 <= z < K
+ // - if z.value + w.value < K:
// add z + w - v = 0
- // - if z.value + w.value >= m:
- // add z + w - v - m = 0
+ // - if z.value + w.value >= K:
+ // add z + w - v - K = 0
//
- // case a != 1, gcd(a, m) = 1
- // - x |-> x*y + a^-1*z, 0 <= z < m
- // - add w = b mod m
- // if z.value + w.value < m
+ // case a != 1, gcd(a, K) = 1
+ // - x |-> x*y + a^-1*z, 0 <= z < K
+ // - add w = b mod K
+ // if z.value + w.value < K
// add z + w - v = 0
- // if z.value + w.value >= m
- // add z + w - v - m = 0
+ // if z.value + w.value >= K
+ // add z + w - v - K = 0
//
- // case a != 1, gcd(a,m) = g != 1
- // - x |-> x*y + a^-1*z, 0 <= z < m
- // a*x + b mod m = v is now
- // g*z + b mod m = v
- // - add w = b mod m
- // - 0 <= g*z.value + w.value < m*(g+1)
- // - add g*z + w - v - k*m = 0 for suitable k from 0 .. g based on model
+ // case a != 1, gcd(a,K) = g != 1
+ // - x |-> x*y + a^-1*z, 0 <= z < K
+ // a*x + b mod K = v is now
+ // g*z + b mod K = v
+ // - add w = b mod K
+ // - 0 <= g*z.value + w.value < K*(g+1)
+ // - add g*z + w - v - k*K = 0 for suitable k from 0 .. g based on model
//
- model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def) {
+ model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& _mod_rows, bool compute_def) {
def result;
- SASSERT(!mod_rows.empty());
- unsigned row_index = mod_rows.back();
- rational a = get_coefficient(row_index, x);
- rational m = m_rows[row_index].m_mod;
- unsigned v = m_rows[row_index].m_id;
+ unsigned_vector mod_rows(_mod_rows);
+ rational K(1);
+ for (unsigned ri : mod_rows)
+ K = lcm(K, m_rows[ri].m_mod);
+
rational x_value = m_var2value[x];
- // disable row
- m_rows[row_index].m_alive = false;
- replace_var(row_index, x, rational::zero());
-
- // compute a_inv
- rational a_inv, m_inv;
- rational g = mod(gcd(a, m, a_inv, m_inv), m);
- if (a_inv.is_neg())
- a_inv = mod(a_inv, m);
- SASSERT(mod(a_inv * a, m) == g);
-
- // solve for x_value = m*y_value + a^-1*z_value, 0 <= z_value < m.
- rational z_value = mod(x_value, m);
- rational y_value = div(x_value, m) - div(a_inv*z_value, m);
- SASSERT(x_value == m*y_value + a_inv*z_value);
- SASSERT(0 <= z_value && z_value < m);
-
+ rational y_value = div(x_value, K);
+ rational z_value = mod(x_value, K);
+ SASSERT(x_value == K * y_value + z_value);
+ SASSERT(0 <= z_value && z_value < K);
// add new variables
- unsigned y = add_var(y_value, true);
unsigned z = add_var(z_value, true);
- // TODO: we could recycle x by either y or z.
+ unsigned y = add_var(y_value, true);
- // replace x by m*y + a^-1*z in other rows.
- unsigned_vector const& row_ids = m_var2row_ids[x];
uint_set visited;
- visited.insert(row_index);
- for (unsigned row_id : row_ids) {
- if (visited.contains(row_id))
+ for (unsigned ri : mod_rows) {
+ m_rows[ri].m_alive = false;
+ visited.insert(ri);
+ }
+
+ // replace x by K*y + z in other rows.
+ for (unsigned ri : m_var2row_ids[x]) {
+ if (visited.contains(ri))
continue;
- replace_var(row_id, x, m, y, a_inv, z);
- visited.insert(row_id);
- normalize(row_id);
+ replace_var(ri, x, K, y, rational::one(), z);
+ visited.insert(ri);
+ normalize(ri);
}
// add bounds for z
add_lower_bound(z, rational::zero());
- add_upper_bound(z, m - 1);
+ add_upper_bound(z, K - 1);
-
- // add w = b mod m
- vector<var> coeffs = m_rows[row_index].m_vars;
- rational coeff = m_rows[row_index].m_coeff;
-
- unsigned w = UINT_MAX;
- rational offset(0);
- if (coeffs.empty())
- offset = mod(coeff, m);
- else
- w = add_mod(coeffs, coeff, m);
+ for (unsigned ri : mod_rows) {
+ rational a = get_coefficient(ri, x);
+ // add w = b mod K
+ vector<var> coeffs = m_rows[ri].m_vars;
+ rational coeff = m_rows[ri].m_coeff;
+ unsigned v = m_rows[ri].m_id;
- rational w_value = w == UINT_MAX ? offset : m_var2value[w];
+ unsigned w = UINT_MAX;
+ rational offset(0);
+ if (coeffs.empty())
+ offset = mod(coeff, K);
+ else
+ w = add_mod(coeffs, coeff, K);
- // add g*z + w - v - k*m = 0, for k = (g*z_value + w_value) div m
- rational km = div(g*z_value + w_value, m)*m + offset;
- vector<var> mod_coeffs;
- if (g != 0) mod_coeffs.push_back(var(z, g));
- if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one()));
- mod_coeffs.push_back(var(v, rational::minus_one()));
- add_constraint(mod_coeffs, km, t_eq);
- add_lower_bound(v, rational::zero());
- add_upper_bound(v, m - 1);
+ rational w_value = w == UINT_MAX ? offset : m_var2value[w];
- // allow to recycle row.
- m_retired_rows.push_back(row_index);
+ // add v = a*z + w - k*K = 0, 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
+ rational kK = div(a * z_value + w_value, K) * K;
+ vector<var> mod_coeffs;
+ 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, kK + offset, t_eq);
+ add_lower_bound(v, rational::zero());
+ add_upper_bound(v, K - 1);
+
+ // allow to recycle row.
+ m_retired_rows.push_back(ri);
+
+ project(v, false);
+ }
- project(v, false);
def y_def = project(y, compute_def);
def z_def = project(z, compute_def);
if (compute_def) {
- result = (y_def * m) + z_def * a_inv;
+ result = (y_def * K) + z_def;
m_var2value[x] = eval(result);
}
TRACE("opt", display(tout << "solve_mod\n"));
@@ -1287,134 +1299,152 @@ namespace opt {
}
//
- // Given v = a*x + b div m
- // Replace x |-> m*y + z
- // - w = b div m
- // - v = ((a*m*y + a*z) + b) div m
- // = a*y + (a*z + b) div m
- // = a*y + b div m + (b mod m + a*z) div m
- // = a*y + b div m + k
- // where k := (b.value mod m + a*z.value) div m
+ // Given v = a*x + b div K
+ // Replace x |-> K*y + z
+ // - w = b div K
+ // - v = ((a*K*y + a*z) + b) div K
+ // = a*y + (a*z + b) div K
+ // = a*y + b div K + (b mod K + a*z) div K
+ // = a*y + b div K + k
+ // where k := (b.value mod K + a*z.value) div K
// k is between 0 and a
//
- // - k*m <= b mod m + a*z < (k+1)*m
+ // - k*K <= b mod K + a*z < (k+1)*K
//
// A better version using a^-1
- // - v = (a*m*y + a^-1*a*z + b) div m
- // = a*y + ((m*A + g)*z + b) div m where we write a*a^-1 = m*A + g
- // = a*y + A + (g*z + b) div m
- // - k*m <= b mod m + gz < (k+1)*m
+ // - v = (a*K*y + a^-1*a*z + b) div K
+ // = a*y + ((K*A + g)*z + b) div K where we write a*a^-1 = K*A + g
+ // = a*y + A + (g*z + b) div K
+ // - k*K <= b Kod m + gz < (k+1)*K
// where k is between 0 and g
- // when gcd(a, m) = 1, then there are only two cases.
+ // when gcd(a, K) = 1, then there are only two cases.
//
- model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& div_rows, bool compute_def) {
+ model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& _div_rows, bool compute_def) {
def result;
+ unsigned_vector div_rows(_div_rows);
SASSERT(!div_rows.empty());
- unsigned row_index = div_rows.back();
- rational a = get_coefficient(row_index, x);
- rational m = m_rows[row_index].m_mod;
- unsigned v = m_rows[row_index].m_id;
+
+ rational K(1);
+ for (unsigned ri : div_rows)
+ K = lcm(K, m_rows[ri].m_mod);
+
rational x_value = m_var2value[x];
-
- // display(verbose_stream(), m_rows[row_index]);
-
- // disable row
- m_rows[row_index].m_alive = false;
- replace_var(row_index, x, rational::zero());
- rational b_value = m_rows[row_index].m_value;
-
- TRACE("opt", display(tout << "solve_div\n"));
-
- // compute a_inv
- rational a_inv, m_inv;
- rational g = mod(gcd(a, m, a_inv, m_inv), m);
- if (a_inv.is_neg())
- a_inv = mod(a_inv, m);
-
- SASSERT(mod(a_inv * a, m) == g);
-
- // solve for x_value = m*y_value + a_inv*z_value, 0 <= z_value < m.
- rational z_value = mod(x_value, m);
- rational y_value = div(x_value, m) - div(z_value*a_inv, m);
- SASSERT(x_value == m*y_value + a_inv*z_value);
- SASSERT(0 <= z_value && z_value < m);
-
+ rational z_value = mod(x_value, K);
+ rational y_value = div(x_value, K);
+ SASSERT(x_value == K * y_value + z_value);
+ SASSERT(0 <= z_value && z_value < K);
// add new variables
- unsigned y = add_var(y_value, true);
unsigned z = add_var(z_value, true);
+ unsigned y = add_var(y_value, true);
- // replace x by m*y + a^-1*z in other rows.
- unsigned_vector const& row_ids = m_var2row_ids[x];
uint_set visited;
- visited.insert(row_index);
- for (unsigned row_id : row_ids) {
- if (visited.contains(row_id))
+ for (unsigned ri : div_rows) {
+ row& r = m_rows[ri];
+ mul(ri, K / r.m_mod);
+ r.m_alive = false;
+ visited.insert(ri);
+ }
+
+ // replace x by K*y + z in other rows.
+ for (unsigned ri : m_var2row_ids[x]) {
+ if (visited.contains(ri))
continue;
- replace_var(row_id, x, m, y, a_inv, z);
- visited.insert(row_id);
- normalize(row_id);
+ replace_var(ri, x, K, y, rational::one(), z);
+ visited.insert(ri);
+ normalize(ri);
}
// add bounds for z
add_lower_bound(z, rational::zero());
- add_upper_bound(z, m - 1);
-
- // add w = b div m
- vector<var> coeffs = m_rows[row_index].m_vars;
- rational coeff = m_rows[row_index].m_coeff;
- unsigned w = UINT_MAX;
- rational offset(0);
- if (coeffs.empty())
- offset = div(coeff, m);
- else
- w = add_div(coeffs, coeff, m);
-
- //
- // w = b div m
- // v = a*y + w + (a*a_inv div m) + k
- // k = (g*z_value + (b_value mod m)) div m
- // k*m <= g*z + b mod m < (k+1)*m
- //
+ add_upper_bound(z, K - 1);
- rational k = div(a*z_value + mod(b_value, m), m) + offset;
- rational n = div(a_inv * a, m);
- vector<var> div_coeffs;
- div_coeffs.push_back(var(v, rational::minus_one()));
- div_coeffs.push_back(var(y, a));
- if (w != UINT_MAX) div_coeffs.push_back(var(w, rational::one()));
- if (n != 0) div_coeffs.push_back(var(z, n));
- add_constraint(div_coeffs, k, t_eq);
+ TRACE("opt", display(tout));
- unsigned u = UINT_MAX;
- offset = 0;
- if (coeffs.empty())
- offset = mod(coeff, m);
- else
- u = add_mod(coeffs, coeff, m);
+ // solve for x_value = K*y_value + z_value, 0 <= z_value < K.
- // add g*z + (b mod m) < (k + 1)*m
- vector<var> bound_coeffs;
- if (g != 0) bound_coeffs.push_back(var(z, g));
- if (u != UINT_MAX) bound_coeffs.push_back(var(u, rational::one()));
- add_constraint(bound_coeffs, 1 - m * (k + 1) + offset, t_le);
+ for (unsigned ri : div_rows) {
- // add k*m <= gz + (b mod m)
- for (auto& c : bound_coeffs)
- c.m_coeff.neg();
- add_constraint(bound_coeffs, k * m - offset, t_le);
- // allow to recycle row.
- m_retired_rows.push_back(row_index);
+ rational a = get_coefficient(ri, x);
+ replace_var(ri, x, rational::zero());
+ // add w = b div m
+ vector<var> coeffs = m_rows[ri].m_vars;
+ rational coeff = m_rows[ri].m_coeff;
+ unsigned w = UINT_MAX;
+ rational offset(0);
+ if (coeffs.empty())
+ offset = div(coeff, K);
+ else
+ w = add_div(coeffs, coeff, K);
+
+ //
+ // w = b div K
+ // v = a*y + w + k
+ // k = (a*z_value + (b_value mod K)) div K
+ // k*K <= a*z + b mod K < (k+1)*K
+ //
+ /**
+ * It is based on the following claim (tested for select values of a, K)
+ * (define-const K Int 13)
+ * (declare-const b Int)
+ * (define-const a Int -11)
+ * (declare-const y Int)
+ * (declare-const z Int)
+ * (define-const w Int (div b K))
+ * (define-const k1 Int (+ (* a z) (mod b K)))
+ * (define-const k Int (div k1 K))
+ * (define-const x Int (+ (* K y) z))
+ * (define-const u Int (+ (* a x) b))
+ * (define-const v Int (+ (* a y) w k))
+ * (assert (<= 0 z))
+ * (assert (< z K))
+ * (assert (<= (* K k) k1))
+ * (assert (< k1 (* K (+ k 1))))
+ * (assert (not (= (div u K) v)))
+ * (check-sat)
+ */
+ unsigned v = m_rows[ri].m_id;
+ rational b_value = eval(coeffs) + coeff;
+ rational k = div(a * z_value + mod(b_value, K), K);
+ vector<var> div_coeffs;
+ div_coeffs.push_back(var(v, rational::minus_one()));
+ div_coeffs.push_back(var(y, a));
+ if (w != UINT_MAX) div_coeffs.push_back(var(w, rational::one()));
+ add_constraint(div_coeffs, k + offset, t_eq);
+
+ unsigned u = UINT_MAX;
+ offset = 0;
+ if (coeffs.empty())
+ offset = mod(coeff, K);
+ else
+ u = add_mod(coeffs, coeff, K);
+
+ // add a*z + (b mod K) < (k + 1)*K
+ vector<var> bound_coeffs;
+ bound_coeffs.push_back(var(z, a));
+ if (u != UINT_MAX) bound_coeffs.push_back(var(u, rational::one()));
+ add_constraint(bound_coeffs, 1 - K * (k + 1) + offset, t_le);
+
+ // add k*K <= az + (b mod K)
+ for (auto& c : bound_coeffs)
+ c.m_coeff.neg();
+ add_constraint(bound_coeffs, k * K - offset, t_le);
+ // allow to recycle row.
+ m_retired_rows.push_back(ri);
+ project(v, false);
+ }
+
+ TRACE("opt", display(tout << "solve_div reduced " << y << " " << z << "\n"));
// project internal variables.
- project(v, false);
+
def y_def = project(y, compute_def);
def z_def = project(z, compute_def);
if (compute_def) {
- result = (y_def * m) + (z_def * a_inv);
+ result = (y_def * K) + z_def;
m_var2value[x] = eval(result);
}
- TRACE("opt", display(tout << "solve_div\n"));
+ TRACE("opt", display(tout << "solve_div done v" << x << "\n"));
return result;
}
@@ -1572,7 +1602,6 @@ namespace opt {
vector<var> coeffs;
mk_coeffs_without(coeffs, r1.m_vars, x);
rational c = mod(-eval(coeffs), a);
- TRACE("opt", tout << "add divides " << eval(coeffs) << " " << a << " " << c << "\n");
add_divides(coeffs, c, a);
}
unsigned_vector const& row_ids = m_var2row_ids[x];
diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h
index 40abf1333..a62c265a9 100644
--- a/src/math/simplex/model_based_opt.h
+++ b/src/math/simplex/model_based_opt.h
@@ -125,7 +125,7 @@ namespace opt {
void mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2);
- void mul_add(unsigned x, rational const& a1, unsigned row_src, rational const& a2, unsigned row_dst);
+ void mul_add(unsigned x, rational a1, unsigned row_src, rational a2, unsigned row_dst);
void mul(unsigned dst, rational const& c);
@@ -139,7 +139,7 @@ namespace opt {
void add_upper_bound(unsigned x, rational const& hi);
- unsigned add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type r, unsigned id);
+ void add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type r, unsigned id);
void replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B);
@@ -187,7 +187,7 @@ namespace opt {
// add a constraint. We assume that the constraint is
// satisfied under the values provided to the variables.
- unsigned add_constraint(vector<var> const& coeffs, rational const& c, ineq_type r);
+ void add_constraint(vector<var> const& coeffs, rational const& c, ineq_type r);
// add a divisibility constraint. The row should divide m.
void add_divides(vector<var> const& coeffs, rational const& c, rational const& m);
diff --git a/src/qe/qsat.cpp b/src/qe/qsat.cpp
index 767b96e5d..6023273c0 100644
--- a/src/qe/qsat.cpp
+++ b/src/qe/qsat.cpp
@@ -241,9 +241,8 @@ namespace qe {
while (sz0 != todo.size()) {
app* a = to_app(todo.back());
todo.pop_back();
- if (mark.is_marked(a)) {
+ if (mark.is_marked(a))
continue;
- }
mark.mark(a);
if (m_lit2pred.find(a, p)) {
@@ -284,9 +283,8 @@ namespace qe {
m_elevel.insert(r, l);
eq = m.mk_eq(r, a);
defs.push_back(eq);
- if (!is_predicate(a, l.max())) {
+ if (!is_predicate(a, l.max()))
insert(r, l);
- }
level.merge(l);
}
}
@@ -637,57 +635,55 @@ namespace qe {
check_cancel();
expr_ref_vector asms(m_asms);
m_pred_abs.get_assumptions(m_model.get(), asms);
- if (m_model.get()) {
+ if (m_model.get())
validate_assumptions(*m_model.get(), asms);
- }
TRACE("qe", tout << asms << "\n";);
solver& s = get_kernel(m_level).s();
lbool res = s.check_sat(asms);
switch (res) {
case l_true:
s.get_model(m_model);
+ CTRACE("qe", !m_model, tout << "no model\n");
if (!m_model)
return l_undef;
SASSERT(validate_defs("check_sat"));
SASSERT(!m_model.get() || validate_assumptions(*m_model.get(), asms));
SASSERT(validate_model(asms));
TRACE("qe", s.display(tout); display(tout << "\n", *m_model.get()); display(tout, asms); );
- if (m_level == 0) {
+ if (m_level == 0)
m_model_save = m_model;
- }
push();
- if (m_level == 1 && m_mode == qsat_maximize) {
+ if (m_level == 1 && m_mode == qsat_maximize)
maximize_model();
- }
break;
case l_false:
switch (m_level) {
case 0:
return l_false;
case 1:
- if (m_mode == qsat_sat) {
+ if (m_mode == qsat_sat)
return l_true;
- }
if (m_model.get()) {
SASSERT(validate_assumptions(*m_model.get(), asms));
- if (!project_qe(asms)) return l_undef;
+ if (!project_qe(asms))
+ return l_undef;
}
- else {
+ else
pop(1);
- }
break;
default:
if (m_model.get()) {
- if (!project(asms)) return l_undef;
+ if (!project(asms))
+ return l_undef;
}
- else {
+ else
pop(1);
- }
break;
}
break;
case l_undef:
+ TRACE("qe", tout << "check-sat is undef\n");
return res;
}
}
@@ -833,11 +829,10 @@ namespace qe {
}
}
- bool get_core(expr_ref_vector& core, unsigned level) {
+ void get_core(expr_ref_vector& core, unsigned level) {
SASSERT(validate_defs("get_core"));
get_kernel(level).get_core(core);
m_pred_abs.pred2lit(core);
- return true;
}
bool minimize_core(expr_ref_vector& core, unsigned level) {
@@ -905,9 +900,7 @@ namespace qe {
SASSERT(m_level == 1);
expr_ref fml(m);
model& mdl = *m_model.get();
- if (!get_core(core, m_level)) {
- return false;
- }
+ get_core(core, m_level);
SASSERT(validate_core(mdl, core));
get_vars(m_level);
SASSERT(validate_assumptions(mdl, core));
@@ -927,7 +920,7 @@ namespace qe {
}
bool project(expr_ref_vector& core) {
- if (!get_core(core, m_level)) return false;
+ get_core(core, m_level);
TRACE("qe", display(tout); display(tout << "core\n", core););
SASSERT(m_level >= 2);
expr_ref fml(m);
@@ -950,14 +943,17 @@ namespace qe {
if (level.max() == UINT_MAX) {
num_scopes = 2*(m_level/2);
}
+ else if (level.max() + 2 > m_level) {
+ // fishy - this can happen.
+ TRACE("qe", tout << "max-level: " << level.max() << " level: " << m_level << "\n");
+ return false;
+ }
else {
- if (level.max() + 2 > m_level) return false;
SASSERT(level.max() + 2 <= m_level);
num_scopes = m_level - level.max();
SASSERT(num_scopes >= 2);
- if ((num_scopes % 2) != 0) {
+ if ((num_scopes % 2) != 0)
--num_scopes;
- }
}
pop(num_scopes);
From 138f0d269c3a1bb8ce9612165b44e390d78eac6c Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Sun, 14 Aug 2022 12:26:33 -0700
Subject: [PATCH 13/16] fix regression found by fuzzers fix #6271
---
src/qe/mbp/mbp_arith.cpp | 11 +++++++----
1 file changed, 7 insertions(+), 4 deletions(-)
diff --git a/src/qe/mbp/mbp_arith.cpp b/src/qe/mbp/mbp_arith.cpp
index 11db720ef..713a63c14 100644
--- a/src/qe/mbp/mbp_arith.cpp
+++ b/src/qe/mbp/mbp_arith.cpp
@@ -220,10 +220,13 @@ namespace mbp {
if (!a.is_numeral(val, r))
throw default_exception("mbp evaluation didn't produce an integer");
c += mul * r;
- // t1 mod mul1 == r
- vars coeffs;
- rational c0 = add_def(t1, mul1, coeffs);
- mbo.add_divides(coeffs, c0 - r, mul1);
+
+ rational c0(-r), mul0(1);
+ obj_map<expr, rational> ts0;
+ linearize(mbo, eval, mul0, t1, c0, fmls, ts0, tids);
+ vars coeffs;
+ extract_coefficients(mbo, eval, ts0, tids, coeffs);
+ mbo.add_divides(coeffs, c0, mul1);
}
else if (false && a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && mul1 > 0) {
// v = t1 mod mul1
From a0d4a8c21c12a5ee6e8e0c338d9a1389683d8f47 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Mon, 15 Aug 2022 00:12:44 -0700
Subject: [PATCH 14/16] update diagnostics
---
src/smt/qi_queue.cpp | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/src/smt/qi_queue.cpp b/src/smt/qi_queue.cpp
index de52cdb01..f8fbe739e 100644
--- a/src/smt/qi_queue.cpp
+++ b/src/smt/qi_queue.cpp
@@ -23,6 +23,7 @@ Revision History:
#include "ast/rewriter/var_subst.h"
#include "smt/smt_context.h"
#include "smt/qi_queue.h"
+#include <iostream>
namespace smt {
@@ -228,15 +229,12 @@ namespace smt {
TRACE("qi_queue", tout << "new instance:\n" << mk_pp(instance, m) << "\n";);
- TRACE("qi_queue_instance", tout << "new instance:\n" << mk_pp(instance, m) << "\n";);
expr_ref s_instance(m);
proof_ref pr(m);
m_context.get_rewriter()(instance, s_instance, pr);
TRACE("qi_queue_bug", tout << "new instance after simplification:\n" << s_instance << "\n";);
if (m.is_true(s_instance)) {
- TRACE("checker", tout << "reduced to true, before:\n" << mk_ll_pp(instance, m););
-
STRACE("instance", tout << "Instance reduced to true\n";);
stat -> inc_num_instances_simplify_true();
if (m.has_trace_stream()) {
@@ -250,9 +248,11 @@ namespace smt {
std::cout << "instantiate\n";
enode_vector _bindings(num_bindings, bindings);
for (auto * b : _bindings)
- std::cout << enode_pp(b, m_context) << " ";
+ std::cout << mk_pp(b->get_expr(), m) << " ";
std::cout << "\n";
std::cout << mk_pp(q, m) << "\n";
+ std::cout << "instance\n";
+ std::cout << instance << "\n";
#endif
TRACE("qi_queue", tout << "simplified instance:\n" << s_instance << "\n";);
From e0aa32e6c59047109f25dc8c0273a68b5ab54a60 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Mon, 15 Aug 2022 00:13:32 -0700
Subject: [PATCH 15/16] fix #6270
MBQI asserts auxiliary function definitions to handle models of arrays. This is unsound if the definition contains a model value.
---
src/smt/smt_model_checker.cpp | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/src/smt/smt_model_checker.cpp b/src/smt/smt_model_checker.cpp
index 039416f73..0ffe8108d 100644
--- a/src/smt/smt_model_checker.cpp
+++ b/src/smt/smt_model_checker.cpp
@@ -218,7 +218,7 @@ namespace smt {
TRACE("model_checker", tout << "Got some value " << sk_value << "\n";);
if (use_inv) {
- unsigned sk_term_gen = 0;
+ unsigned sk_term_gen = 0;
expr * sk_term = m_model_finder.get_inv(q, i, sk_value, sk_term_gen);
if (sk_term != nullptr) {
TRACE("model_checker", tout << "Found inverse " << mk_pp(sk_term, m) << "\n";);
@@ -243,6 +243,8 @@ namespace smt {
func_decl * f = nullptr;
if (autil.is_as_array(sk_value, f) && cex->get_func_interp(f) && cex->get_func_interp(f)->get_interp()) {
expr_ref body(cex->get_func_interp(f)->get_interp(), m);
+ if (contains_model_value(body))
+ return false;
ptr_vector<sort> sorts(f->get_arity(), f->get_domain());
svector<symbol> names;
for (unsigned i = 0; i < f->get_arity(); ++i) {
From 88f4664c652c3523a248239ba90a256181815f2b Mon Sep 17 00:00:00 2001
From: jofleish <jofleish@microsoft.com>
Date: Mon, 15 Aug 2022 09:11:20 -0400
Subject: [PATCH 16/16] Standardize ubutu-latest vmImage
---
azure-pipelines.yml | 6 +++---
scripts/jsdoctest.yml | 2 +-
scripts/mk_nuget_task.py | 3 ++-
scripts/nightly.yaml | 6 +++---
scripts/release.yml | 2 +-
5 files changed, 10 insertions(+), 9 deletions(-)
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index d29a3dd51..f338a5d98 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -21,7 +21,7 @@ jobs:
- job: "LinuxPythonDebug"
displayName: "Ubuntu build - python make - debug"
pool:
- vmImage: "Ubuntu-latest"
+ vmImage: "ubuntu-latest"
strategy:
matrix:
MT:
@@ -102,7 +102,7 @@ jobs:
displayName: "Ubuntu build - cmake"
condition: eq(0,1)
pool:
- vmImage: "Ubuntu-latest"
+ vmImage: "ubuntu-latest"
strategy:
matrix:
msanClang:
@@ -135,7 +135,7 @@ jobs:
- job: "UbuntuCMake"
displayName: "Ubuntu build - cmake"
pool:
- vmImage: "Ubuntu-latest"
+ vmImage: "ubuntu-latest"
strategy:
matrix:
releaseClang:
diff --git a/scripts/jsdoctest.yml b/scripts/jsdoctest.yml
index 1a3a5286f..0556a6cf5 100644
--- a/scripts/jsdoctest.yml
+++ b/scripts/jsdoctest.yml
@@ -12,7 +12,7 @@ stages:
- job: UbuntuDoc
displayName: "Ubuntu Doc build"
pool:
- vmImage: "Ubuntu-latest"
+ vmImage: "ubuntu-latest"
steps:
# TODO setup emscripten with no-install, then run
- script: npm --prefix=src/api/js ci
diff --git a/scripts/mk_nuget_task.py b/scripts/mk_nuget_task.py
index 75d26a1c5..dd81349df 100644
--- a/scripts/mk_nuget_task.py
+++ b/scripts/mk_nuget_task.py
@@ -21,7 +21,8 @@ def mk_dir(d):
if not os.path.exists(d):
os.makedirs(d)
-os_info = { 'ubuntu-18' : ('so', 'linux-x64'),
+os_info = { 'ubuntu-latest' : ('so', 'linux-x64'),
+ 'ubuntu-18' : ('so', 'linux-x64'),
'ubuntu-20' : ('so', 'linux-x64'),
'glibc-2.31' : ('so', 'linux-x64'),
'x64-win' : ('dll', 'win-x64'),
diff --git a/scripts/nightly.yaml b/scripts/nightly.yaml
index 33bed5b55..98a378e59 100644
--- a/scripts/nightly.yaml
+++ b/scripts/nightly.yaml
@@ -55,7 +55,7 @@ stages:
- job: UbuntuDoc
displayName: "Ubuntu Doc build"
pool:
- vmImage: "Ubuntu-latest"
+ vmImage: "ubuntu-latest"
steps:
- script: sudo apt-get install ocaml opam libgmp-dev
- script: opam init -y
@@ -92,7 +92,7 @@ stages:
name: ManyLinux
displayName: "ManyLinux build"
pool:
- vmImage: "Ubuntu-18.04"
+ vmImage: "ubuntu-latest"
container: "quay.io/pypa/manylinux2010_x86_64:latest"
steps:
- script: $(python) scripts/mk_unix_dist.py --nodotnet --nojava
@@ -111,7 +111,7 @@ stages:
# name: MuslLinux
# displayName: "MuslLinux build"
# pool:
-# vmImage: "Ubuntu-18.04"
+# vmImage: "ubuntu-latest"
# container: "quay.io/pypa/musllinux_1_1_x86_64:latest"
# steps:
# - script: $(python) scripts/mk_unix_dist.py --nodotnet --nojava
diff --git a/scripts/release.yml b/scripts/release.yml
index f305e7c74..969e3697e 100644
--- a/scripts/release.yml
+++ b/scripts/release.yml
@@ -91,7 +91,7 @@ stages:
- job: UbuntuDoc
displayName: "Ubuntu Doc build"
pool:
- vmImage: "Ubuntu-latest"
+ vmImage: "ubuntu-latest"
steps:
- script: sudo apt-get install ocaml opam libgmp-dev
- script: opam init -y