mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-28 11:25:51 +00:00
Code Activity

Merge branch 'master' of https://github.com/Z3Prover/z3

Browse source
This commit is contained in:
Nikolaj Bjorner 2016-02-04 08:12:56 -08:00
commit 768bb84798
13 changed files with 1632 additions and 1625 deletions
Download patch file Download diff file Expand all files Collapse all files

61
src/ast/expr2polynomial.cpp
View file

@ -8,8 +8,8 @@ Module Name:
Abstract:
Translator from Z3 expressions into multivariate polynomials (and back).
Author:
Leonardo (leonardo) 2011-12-23
@ -23,6 +23,7 @@ Notes:
#include"ast_smt2_pp.h"
#include"z3_exception.h"
#include"cooperate.h"
#include"common_msgs.h"
struct expr2polynomial::imp {
struct frame {
@ -31,7 +32,7 @@ struct expr2polynomial::imp {
frame():m_curr(0), m_idx(0) {}
frame(app * t):m_curr(t), m_idx(0) {}
};
expr2polynomial & m_wrapper;
ast_manager & m_am;
arith_util m_autil;
@ -39,7 +40,7 @@ struct expr2polynomial::imp {
expr2var * m_expr2var;
bool m_expr2var_owner;
expr_ref_vector m_var2expr;
obj_map<expr, unsigned> m_cache;
expr_ref_vector m_cached_domain;
polynomial::polynomial_ref_vector m_cached_polynomials;
@ -52,7 +53,7 @@ struct expr2polynomial::imp {
bool m_use_var_idxs;
volatile bool m_cancel;
imp(expr2polynomial & w, ast_manager & am, polynomial::manager & pm, expr2var * e2v, bool use_var_idxs):
m_wrapper(w),
m_am(am),
@ -94,7 +95,7 @@ struct expr2polynomial::imp {
void checkpoint() {
if (m_cancel)
throw default_exception("canceled");
throw default_exception(Z3_CANCELED_MSG);
cooperate("expr2polynomial");
}
@ -114,7 +115,7 @@ struct expr2polynomial::imp {
SASSERT(!m_cache.contains(t));
SASSERT(m_cached_denominators.size() == m_cached_polynomials.size());
SASSERT(m_cached_denominators.size() == m_cached_domain.size());
if (t->get_ref_count() <= 1)
if (t->get_ref_count() <= 1)
return;
unsigned idx = m_cached_polynomials.size();
m_cache.insert(t, idx);
@ -176,7 +177,7 @@ struct expr2polynomial::imp {
case OP_NUM:
store_const_poly(t);
return true;
case OP_ADD: case OP_SUB: case OP_MUL: case OP_UMINUS: case OP_TO_REAL:
case OP_ADD: case OP_SUB: case OP_MUL: case OP_UMINUS: case OP_TO_REAL:
push_frame(t);
return false;
case OP_POWER: {
@ -215,7 +216,7 @@ struct expr2polynomial::imp {
store_var_poly(t);
return true;
}
SASSERT(is_app(t));
if (!m_autil.is_arith_expr(t)) {
if (m_use_var_idxs)
@ -223,7 +224,7 @@ struct expr2polynomial::imp {
store_var_poly(t);
return true;
}
return visit_arith_app(to_app(t));
}
@ -237,7 +238,7 @@ struct expr2polynomial::imp {
polynomial::polynomial * const * polynomial_args(unsigned num_args) {
SASSERT(m_presult_stack.size() >= num_args);
return m_presult_stack.c_ptr() + m_presult_stack.size() - num_args;
}
}
polynomial::numeral const * denominator_args(unsigned num_args) {
SASSERT(m_dresult_stack.size() >= num_args);
@ -330,34 +331,34 @@ struct expr2polynomial::imp {
// do nothing
cache_result(t);
}
void process_app(app * t) {
SASSERT(m_presult_stack.size() == m_dresult_stack.size());
switch (t->get_decl_kind()) {
case OP_ADD:
case OP_ADD:
process_add(t);
return;
case OP_SUB:
process_sub(t);
return;
case OP_MUL:
case OP_MUL:
process_mul(t);
return;
case OP_POWER:
process_power(t);
return;
case OP_UMINUS:
case OP_UMINUS:
process_uminus(t);
return;
case OP_TO_REAL:
case OP_TO_REAL:
process_to_real(t);
return;
default:
UNREACHABLE();
}
}
bool to_polynomial(expr * t, polynomial::polynomial_ref & p, polynomial::scoped_numeral & d) {
if (!is_int_real(t))
return false;
@ -373,7 +374,7 @@ struct expr2polynomial::imp {
while (fr.m_idx < num_args) {
expr * arg = a->get_arg(fr.m_idx);
fr.m_idx++;
if (!visit(arg))
if (!visit(arg))
goto begin_loop;
}
process_app(a);
@ -470,12 +471,12 @@ expr2polynomial::~expr2polynomial() {
dealloc(m_imp);
}
ast_manager & expr2polynomial::m() const {
return m_imp->m_am;
ast_manager & expr2polynomial::m() const {
return m_imp->m_am;
}
polynomial::manager & expr2polynomial::pm() const {
return m_imp->m_pm;
polynomial::manager & expr2polynomial::pm() const {
return m_imp->m_pm;
}
bool expr2polynomial::to_polynomial(expr * t, polynomial::polynomial_ref & p, polynomial::scoped_numeral & d) {
@ -486,14 +487,14 @@ void expr2polynomial::to_expr(polynomial::polynomial_ref const & p, bool use_pow
m_imp->to_expr(p, use_power, r);
}
bool expr2polynomial::is_var(expr * t) const {
bool expr2polynomial::is_var(expr * t) const {
SASSERT(!m_imp->m_use_var_idxs);
return m_imp->m_expr2var->is_var(t);
return m_imp->m_expr2var->is_var(t);
}
expr2var const & expr2polynomial::get_mapping() const {
expr2var const & expr2polynomial::get_mapping() const {
SASSERT(!m_imp->m_use_var_idxs);
return *(m_imp->m_expr2var);
return *(m_imp->m_expr2var);
}
void expr2polynomial::set_cancel(bool f) {
@ -511,8 +512,8 @@ bool default_expr2polynomial::is_int(polynomial::var x) const {
return m_is_int[x];
}
polynomial::var default_expr2polynomial::mk_var(bool is_int) {
polynomial::var x = pm().mk_var();
polynomial::var default_expr2polynomial::mk_var(bool is_int) {
polynomial::var x = pm().mk_var();
m_is_int.reserve(x+1, false);
m_is_int[x] = is_int;
return x;

209
src/math/interval/interval_def.h
View file

@ -24,14 +24,15 @@ Revision History:
#include"trace.h"
#include"scoped_numeral.h"
#include"cooperate.h"
#include"common_msgs.h"
#define DEFAULT_PI_PRECISION 2
// #define TRACE_NTH_ROOT
template<typename C>
interval_manager<C>::interval_manager(reslimit& lim, C const & c): m_limit(lim), m_c(c) {
m().set(m_minus_one, -1);
interval_manager<C>::interval_manager(reslimit& lim, C const & c): m_limit(lim), m_c(c) {
m().set(m_minus_one, -1);
m().set(m_one, 1);
m_pi_n = 0;
}
@ -63,7 +64,7 @@ void interval_manager<C>::del(interval & a) {
template<typename C>
void interval_manager<C>::checkpoint() {
if (!m_limit.inc())
throw default_exception("canceled");
throw default_exception(Z3_CANCELED_MSG);
cooperate("interval");
}
@ -82,7 +83,7 @@ void interval_manager<C>::nth_root_slow(numeral const & a, unsigned n, numeral c
static unsigned counter = 0;
static unsigned loop_counter = 0;
counter++;
if (counter % 1000 == 0)
if (counter % 1000 == 0)
std::cerr << "[nth-root] " << counter << " " << loop_counter << " " << ((double)loop_counter)/((double)counter) << std::endl;
#endif
@ -163,7 +164,7 @@ void interval_manager<C>::nth_root_slow(numeral const & a, unsigned n, numeral c
/**
\brief Store in o a rough approximation of a^1/n.
It uses 2^Floor[Floor(Log2(a))/n]
\pre is_pos(a)
@ -179,7 +180,7 @@ void interval_manager<C>::rough_approx_nth_root(numeral const & a, unsigned n, n
}
/*
Compute the n-th root of \c a with (suggested) precision p.
Compute the n-th root of \c a with (suggested) precision p.
The only guarantee provided by this method is that a^(1/n) is in [lo, hi].
If n is even, then a is assumed to be nonnegative.
@ -202,8 +203,8 @@ void interval_manager<C>::nth_root(numeral const & a, unsigned n, numeral const
m().abs(A);
nth_root_pos(A, n, p, lo, hi);
STRACE("nth_root_trace",
tout << "[nth-root] ("; m().display(tout, A); tout << ")^(1/" << n << ") >= "; m().display(tout, lo); tout << "\n";
STRACE("nth_root_trace",
tout << "[nth-root] ("; m().display(tout, A); tout << ")^(1/" << n << ") >= "; m().display(tout, lo); tout << "\n";
tout << "[nth-root] ("; m().display(tout, A); tout << ")^(1/" << n << ") <= "; m().display(tout, hi); tout << "\n";);
if (is_neg) {
m().swap(lo, hi);
@ -214,7 +215,7 @@ void interval_manager<C>::nth_root(numeral const & a, unsigned n, numeral const
/**
r <- A/(x^n)
If to_plus_inf, then r >= A/(x^n)
If not to_plus_inf, then r <= A/(x^n)
@ -251,7 +252,7 @@ void interval_manager<C>::A_div_x_n(numeral const & A, numeral const & x, unsign
The computation stops when the difference between current and new x is less than p.
The procedure may not terminate if m() is not precise and p is very small.
*/
template<typename C>
void interval_manager<C>::approx_nth_root(numeral const & A, unsigned n, numeral const & p, numeral & x) {
@ -261,18 +262,18 @@ void interval_manager<C>::approx_nth_root(numeral const & A, unsigned n, numeral
static unsigned counter = 0;
static unsigned loop_counter = 0;
counter++;
if (counter % 1000 == 0)
if (counter % 1000 == 0)
std::cerr << "[nth-root] " << counter << " " << loop_counter << " " << ((double)loop_counter)/((double)counter) << std::endl;
#endif
_scoped_numeral<numeral_manager> x_prime(m()), d(m());
m().set(d, 1);
if (m().lt(A, d))
if (m().lt(A, d))
m().set(x, A);
else
rough_approx_nth_root(A, n, x);
round_to_minus_inf();
if (n == 2) {
@ -297,8 +298,8 @@ void interval_manager<C>::approx_nth_root(numeral const & A, unsigned n, numeral
_scoped_numeral<numeral_manager> _n(m()), _n_1(m());
m().set(_n, n); // _n contains n
m().set(_n_1, n);
m().dec(_n_1); // _n_1 contains n-1
m().dec(_n_1); // _n_1 contains n-1
while (true) {
checkpoint();
#ifdef TRACE_NTH_ROOT
@ -311,7 +312,7 @@ void interval_manager<C>::approx_nth_root(numeral const & A, unsigned n, numeral
m().div(x_prime, _n, x_prime);
m().sub(x_prime, x, d);
m().abs(d);
TRACE("nth_root",
TRACE("nth_root",
tout << "A: "; m().display(tout, A); tout << "\n";
tout << "x: "; m().display(tout, x); tout << "\n";
tout << "x_prime: "; m().display(tout, x_prime); tout << "\n";
@ -340,7 +341,7 @@ void interval_manager<C>::nth_root_pos(numeral const & A, unsigned n, numeral co
else {
// Check if hi is really a upper bound for A^(n-1)
A_div_x_n(A, hi, n-1, true /* lo will be greater than the actual lower bound */, lo);
TRACE("nth_root_bug",
TRACE("nth_root_bug",
tout << "Assuming upper\n";
tout << "A: "; m().display(tout, A); tout << "\n";
tout << "hi: "; m().display(tout, hi); tout << "\n";
@ -397,12 +398,12 @@ template<typename C>
void interval_manager<C>::sine_series(numeral const & a, unsigned k, bool upper, numeral & o) {
SASSERT(k % 2 == 1);
// Compute sine using taylor series up to k
// x - x^3/3! + x^5/5! - x^7/7! + ...
// x - x^3/3! + x^5/5! - x^7/7! + ...
// The result should be greater than or equal to the actual value if upper == true
// Otherwise it must be less than or equal to the actual value.
// The argument upper only matter if the numeral_manager is not precise.
// Taylor series up to k with rounding to
// Taylor series up to k with rounding to
_scoped_numeral<numeral_manager> f(m());
_scoped_numeral<numeral_manager> aux(m());
m().set(o, a);
@ -412,7 +413,7 @@ void interval_manager<C>::sine_series(numeral const & a, unsigned k, bool upper,
TRACE("sine_bug", tout << "[begin-loop] o: " << m().to_rational_string(o) << "\ni: " << i << "\n";
tout << "upper: " << upper << ", upper_factor: " << upper_factor << "\n";
tout << "o (default): " << m().to_string(o) << "\n";);
set_rounding(upper_factor);
set_rounding(upper_factor);
m().power(a, i, f);
TRACE("sine_bug", tout << "a^i " << m().to_rational_string(f) << "\n";);
set_rounding(!upper_factor);
@ -441,15 +442,15 @@ void interval_manager<C>::sine(numeral const & a, unsigned k, numeral & lo, nume
m().reset(hi);
return;
}
// Compute sine using taylor series
// x - x^3/3! + x^5/5! - x^7/7! + ...
//
// Note that, the coefficient of even terms is 0.
// So, we force k to be odd to make sure the error is minimized.
if (k % 2 == 0)
k++;
k++;
// Taylor series error = |x|^(k+1)/(k+1)!
_scoped_numeral<numeral_manager> error(m());
_scoped_numeral<numeral_manager> aux(m());
@ -466,7 +467,7 @@ void interval_manager<C>::sine(numeral const & a, unsigned k, numeral & lo, nume
m().div(error, aux, error);
TRACE("sine", tout << "error: " << m().to_rational_string(error) << "\n";);
// Taylor series up to k with rounding to -oo
// Taylor series up to k with rounding to -oo
sine_series(a, k, false, lo);
if (m().precise()) {
@ -508,7 +509,7 @@ void interval_manager<C>::cosine_series(numeral const & a, unsigned k, bool uppe
// The argument upper only matter if the numeral_manager is not precise.
// Taylor series up to k with rounding to -oo
// Taylor series up to k with rounding to -oo
_scoped_numeral<numeral_manager> f(m());
_scoped_numeral<numeral_manager> aux(m());
m().set(o, 1);
@ -527,7 +528,7 @@ void interval_manager<C>::cosine_series(numeral const & a, unsigned k, bool uppe
else
m().add(o, f, o);
sign = !sign;
upper_factor = !upper_factor;
upper_factor = !upper_factor;
}
}
@ -540,15 +541,15 @@ void interval_manager<C>::cosine(numeral const & a, unsigned k, numeral & lo, nu
m().set(hi, 1);
return;
}
// Compute cosine using taylor series
// 1 - x^2/2! + x^4/4! - x^6/6! + ...
//
// Note that, the coefficient of odd terms is 0.
// So, we force k to be even to make sure the error is minimized.
if (k % 2 == 1)
k++;
k++;
// Taylor series error = |x|^(k+1)/(k+1)!
_scoped_numeral<numeral_manager> error(m());
_scoped_numeral<numeral_manager> aux(m());
@ -563,9 +564,9 @@ void interval_manager<C>::cosine(numeral const & a, unsigned k, numeral & lo, nu
m().div(error, aux, error);
TRACE("sine", tout << "error: "; m().display_decimal(tout, error, 32); tout << "\n";);
// Taylor series up to k with rounding to -oo
// Taylor series up to k with rounding to -oo
cosine_series(a, k, false, lo);
if (m().precise()) {
m().set(hi, lo);
m().sub(lo, error, lo);
@ -614,12 +615,12 @@ void interval_manager<C>::reset(interval & a) {
template<typename C>
bool interval_manager<C>::contains_zero(interval const & n) const {
return
return
(lower_is_neg(n) || (lower_is_zero(n) && !lower_is_open(n))) &&
(upper_is_pos(n) || (upper_is_zero(n) && !upper_is_open(n)));
}
template<typename C>
bool interval_manager<C>::contains(interval const & n, numeral const & v) const {
if (!lower_is_inf(n)) {
@ -688,7 +689,7 @@ void interval_manager<C>::set(interval & t, interval const & s) {
template<typename C>
bool interval_manager<C>::eq(interval const & a, interval const & b) const {
return
return
::eq(m(), lower(a), lower_kind(a), lower(b), lower_kind(b)) &&
::eq(m(), upper(a), upper_kind(a), upper(b), upper_kind(b)) &&
lower_is_open(a) == lower_is_open(b) &&
@ -792,7 +793,7 @@ void interval_manager<C>::add(interval const & a, interval const & b, interval &
add_jst(a, b, c_deps);
add(a, b, c);
}
template<typename C>
void interval_manager<C>::add(interval const & a, interval const & b, interval & c) {
ext_numeral_kind new_l_kind, new_u_kind;
@ -869,7 +870,7 @@ void interval_manager<C>::div_mul(numeral const & k, interval const & a, interva
round_to_minus_inf();
m().inv(k, m_inv_k);
::mul(m(), l, l_k, m_inv_k, EN_NUMERAL, new_l_val, new_l_kind);
round_to_plus_inf();
m().inv(k, m_inv_k);
::mul(m(), u, u_k, m_inv_k, EN_NUMERAL, new_u_val, new_u_kind);
@ -975,7 +976,7 @@ void interval_manager<C>::mul_jst(interval const & i1, interval const & i2, inte
}
}
else {
SASSERT(is_P(i1));
SASSERT(is_P(i1));
if (is_N(i2)) {
// 0 <= a <= x <= b, c <= y <= d <= 0 --> x*y <= b*c (uses the fact that x is pos (a is not neg) or y is neg (d is not pos))
// 0 <= a <= x, y <= d <= 0 --> a*d <= x*y
@ -990,7 +991,7 @@ void interval_manager<C>::mul_jst(interval const & i1, interval const & i2, inte
}
else {
SASSERT(is_P(i2));
// 0 <= a <= x, 0 <= c <= y --> a*c <= x*y
// 0 <= a <= x, 0 <= c <= y --> a*c <= x*y
// x <= b, y <= d --> x*y <= b*d (uses the fact that x is pos (a is not negative) or y is pos (c is not negative))
r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2;
r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_LOWER2
@ -1075,11 +1076,11 @@ void interval_manager<C>::mul(interval const & i1, interval const & i2, interval
// x <= b <= 0, 0 <= c <= y --> x*y <= b*c
TRACE("interval_bug", tout << "(N, P) #" << call_id << "\n";);
SASSERT(is_P(i2));
// must update upper_is_open first, since value of is_N0(i1) and is_P0(i2) may be affected by update
// must update upper_is_open first, since value of is_N0(i1) and is_P0(i2) may be affected by update
set_upper_is_open(r, (is_N0(i1) || is_P0(i2)) ? false : (b_o || c_o));
set_lower_is_open(r, a_o || d_o);
round_to_minus_inf();
::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
round_to_plus_inf();
@ -1093,7 +1094,7 @@ void interval_manager<C>::mul(interval const & i1, interval const & i2, interval
TRACE("interval_bug", tout << "(M, N) #" << call_id << "\n";);
set_lower_is_open(r, b_o || c_o);
set_upper_is_open(r, a_o || c_o);
set_upper_is_open(r, a_o || c_o);
round_to_minus_inf();
::mul(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
@ -1110,7 +1111,7 @@ void interval_manager<C>::mul(interval const & i1, interval const & i2, interval
bool bc_o = b_o || c_o;
bool ac_o = a_o || c_o;
bool bd_o = b_o || d_o;
round_to_minus_inf();
::mul(m(), a, a_k, d, d_k, ad, ad_k);
::mul(m(), b, b_k, c, c_k, bc, bc_k);
@ -1129,7 +1130,7 @@ void interval_manager<C>::mul(interval const & i1, interval const & i2, interval
set_lower_is_open(r, bc_o);
}
if (::gt(m(), ac, ac_k, bd, bd_k) || (::eq(m(), ac, ac_k, bd, bd_k) && !ac_o && bd_o)) {
m().swap(new_u_val, ac);
new_u_kind = ac_k;
@ -1148,7 +1149,7 @@ void interval_manager<C>::mul(interval const & i1, interval const & i2, interval
SASSERT(is_P(i2));
set_lower_is_open(r, a_o || d_o);
set_upper_is_open(r, b_o || d_o);
set_upper_is_open(r, b_o || d_o);
round_to_minus_inf();
::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
@ -1157,13 +1158,13 @@ void interval_manager<C>::mul(interval const & i1, interval const & i2, interval
}
}
else {
SASSERT(is_P(i1));
SASSERT(is_P(i1));
if (is_N(i2)) {
// 0 <= a <= x <= b, c <= y <= d <= 0 --> x*y <= b*c (uses the fact that x is pos (a is not neg) or y is neg (d is not pos))
// 0 <= a <= x, y <= d <= 0 --> a*d <= x*y
TRACE("interval_bug", tout << "(P, N) #" << call_id << "\n";);
// must update upper_is_open first, since value of is_P0(i1) and is_N0(i2) may be affected by update
// must update upper_is_open first, since value of is_P0(i1) and is_N0(i2) may be affected by update
set_upper_is_open(r, (is_P0(i1) || is_N0(i2)) ? false : a_o || d_o);
set_lower_is_open(r, b_o || c_o);
@ -1187,7 +1188,7 @@ void interval_manager<C>::mul(interval const & i1, interval const & i2, interval
}
else {
SASSERT(is_P(i2));
// 0 <= a <= x, 0 <= c <= y --> a*c <= x*y
// 0 <= a <= x, 0 <= c <= y --> a*c <= x*y
// x <= b, y <= d --> x*y <= b*d (uses the fact that x is pos (a is not negative) or y is pos (c is not negative))
TRACE("interval_bug", tout << "(P, P) #" << call_id << "\n";);
@ -1232,7 +1233,7 @@ void interval_manager<C>::power_jst(interval const & a, unsigned n, interval_dep
else if (upper_is_neg(a)) {
// [l, u]^n = [u^n, l^n] if u < 0
// l <= x <= u < 0 --> x^n <= l^n (use lower and upper bound -- need the fact that x is negative)
// x <= u < 0 --> u^n <= x^n
// x <= u < 0 --> u^n <= x^n
b_deps.m_lower_deps = DEP_IN_UPPER1;
if (lower_is_inf(a))
b_deps.m_upper_deps = 0;
@ -1252,7 +1253,7 @@ void interval_manager<C>::power_jst(interval const & a, unsigned n, interval_dep
b_deps.m_lower_deps = 0;
else
b_deps.m_lower_deps = DEP_IN_LOWER1;
if (upper_is_inf(a))
b_deps.m_upper_deps = 0;
else
@ -1298,15 +1299,15 @@ void interval_manager<C>::power(interval const & a, unsigned n, interval & b) {
else if (upper_is_neg(a)) {
// [l, u]^n = [u^n, l^n] if u < 0
// l <= x <= u < 0 --> x^n <= l^n (use lower and upper bound -- need the fact that x is negative)
// x <= u < 0 --> u^n <= x^n
// x <= u < 0 --> u^n <= x^n
SASSERT(!upper_is_inf(a));
bool lower_a_open = lower_is_open(a), upper_a_open = upper_is_open(a);
bool lower_a_inf = lower_is_inf(a);
m().set(lower(b), lower(a));
m().set(upper(b), upper(a));
m().swap(lower(b), upper(b)); // we use a swap because a and b can be aliased
round_to_minus_inf();
m().power(lower(b), n, lower(b));
@ -1337,7 +1338,7 @@ void interval_manager<C>::power(interval const & a, unsigned n, interval & b) {
round_to_plus_inf();
::power(m(), un1, un1_kind, n);
::power(m(), un2, un2_kind, n);
if (::gt(m(), un1, un1_kind, un2, un2_kind) || (::eq(m(), un1, un1_kind, un2, un2_kind) && !lower_is_open(a) && upper_is_open(a))) {
m().swap(upper(b), un1);
set_upper_is_inf(b, un1_kind == EN_PLUS_INFINITY);
@ -1348,7 +1349,7 @@ void interval_manager<C>::power(interval const & a, unsigned n, interval & b) {
set_upper_is_inf(b, un2_kind == EN_PLUS_INFINITY);
set_upper_is_open(b, upper_is_open(a));
}
m().reset(lower(b));
set_lower_is_inf(b, false);
set_lower_is_open(b, false);
@ -1359,8 +1360,8 @@ void interval_manager<C>::power(interval const & a, unsigned n, interval & b) {
if (lower_is_inf(a)) {
reset_lower(b);
}
else {
m().power(lower(a), n, lower(b));
else {
m().power(lower(a), n, lower(b));
set_lower_is_inf(b, false);
set_lower_is_open(b, lower_is_open(a));
}
@ -1409,7 +1410,7 @@ void interval_manager<C>::nth_root(interval const & a, unsigned n, numeral const
set_lower_is_open(b, lower_is_open(a) && m().eq(lo, hi));
m().set(lower(b), lo);
}
if (upper_is_inf(a)) {
m().reset(upper(b));
set_upper_is_inf(b, true);
@ -1423,7 +1424,7 @@ void interval_manager<C>::nth_root(interval const & a, unsigned n, numeral const
set_upper_is_open(b, upper_is_open(a) && m().eq(lo, hi));
m().set(upper(b), hi);
}
TRACE("interval_nth_root", display(tout, a); tout << " --> "; display(tout, b); tout << "\n";);
TRACE("interval_nth_root", display(tout, a); tout << " --> "; display(tout, b); tout << "\n";);
}
template<typename C>
@ -1523,17 +1524,17 @@ void interval_manager<C>::inv(interval const & a, interval & b) {
numeral & new_l_val = m_result_lower;
numeral & new_u_val = m_result_upper;
ext_numeral_kind new_l_kind, new_u_kind;
if (is_P1(a)) {
// 0 < l <= x --> 1/x <= 1/l
// 0 < l <= x <= u --> 1/u <= 1/x (use lower and upper bounds)
round_to_minus_inf();
m().set(new_l_val, upper(a)); new_l_kind = upper_kind(a);
::inv(m(), new_l_val, new_l_kind);
SASSERT(new_l_kind == EN_NUMERAL);
bool new_l_open = upper_is_open(a);
if (lower_is_zero(a)) {
SASSERT(lower_is_open(a));
m().reset(upper(b));
@ -1542,15 +1543,15 @@ void interval_manager<C>::inv(interval const & a, interval & b) {
}
else {
round_to_plus_inf();
m().set(new_u_val, lower(a));
m().set(new_u_val, lower(a));
m().inv(new_u_val);
m().swap(upper(b), new_u_val);
set_upper_is_inf(b, false);
set_upper_is_open(b, lower_is_open(a));
}
m().swap(lower(b), new_l_val);
set_lower_is_inf(b, false);
m().swap(lower(b), new_l_val);
set_lower_is_inf(b, false);
set_lower_is_open(b, new_l_open);
}
else if (is_N1(a)) {
@ -1577,8 +1578,8 @@ void interval_manager<C>::inv(interval const & a, interval & b) {
set_lower_is_inf(b, false);
set_lower_is_open(b, upper_is_open(a));
}
m().swap(upper(b), new_u_val);
m().swap(upper(b), new_u_val);
set_upper_is_inf(b, false);
set_upper_is_open(b, new_u_open);
}
@ -1610,12 +1611,12 @@ void interval_manager<C>::div_jst(interval const & i1, interval const & i2, inte
// x <= b <= 0, c <= y <= d < 0 --> b/c <= x/y
// a <= x <= b <= 0, y <= d < 0 --> x/y <= a/d
r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2;
r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2;
}
else {
// a <= x, a < 0, 0 < c <= y --> a/c <= x/y
// x <= b <= 0, 0 < c <= y <= d --> x/y <= b/d
r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2;
r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2;
r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
}
}
@ -1634,7 +1635,7 @@ void interval_manager<C>::div_jst(interval const & i1, interval const & i2, inte
}
}
else {
SASSERT(is_P(i1));
SASSERT(is_P(i1));
if (is_N1(i2)) {
// b > 0, x <= b, c <= y <= d < 0 --> b/d <= x/y
// 0 <= a <= x, c <= y <= d < 0 --> x/y <= a/c
@ -1665,7 +1666,7 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
#endif
SASSERT(!contains_zero(i2));
SASSERT(&i1 != &r);
if (is_zero(i1)) {
TRACE("interval_bug", tout << "div #" << call_id << "\n"; display(tout, i1); tout << "\n"; display(tout, i2); tout << "\n";);
@ -1682,12 +1683,12 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
numeral const & b = upper(i1); ext_numeral_kind b_k = upper_kind(i1);
numeral const & c = lower(i2); ext_numeral_kind c_k = lower_kind(i2);
numeral const & d = upper(i2); ext_numeral_kind d_k = upper_kind(i2);
bool a_o = lower_is_open(i1);
bool b_o = upper_is_open(i1);
bool c_o = lower_is_open(i2);
bool d_o = upper_is_open(i2);
numeral & new_l_val = m_result_lower;
numeral & new_u_val = m_result_upper;
ext_numeral_kind new_l_kind, new_u_kind;
@ -1698,7 +1699,7 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
tout << "c: "; m().display(tout, c); tout << "\n";
tout << "d: "; m().display(tout, d); tout << "\n";
);
if (is_N(i1)) {
if (is_N1(i2)) {
// x <= b <= 0, c <= y <= d < 0 --> b/c <= x/y
@ -1707,7 +1708,7 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
set_lower_is_open(r, is_N0(i1) ? false : b_o || c_o);
set_upper_is_open(r, a_o || d_o);
round_to_minus_inf();
::div(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
if (m().is_zero(d)) {
@ -1725,10 +1726,10 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
// x <= b <= 0, 0 < c <= y <= d --> x/y <= b/d
TRACE("interval_bug", tout << "(N, P) #" << call_id << "\n";);
SASSERT(is_P1(i2));
set_upper_is_open(r, is_N0(i1) ? false : (b_o || d_o));
set_lower_is_open(r, a_o || c_o);
if (m().is_zero(c)) {
SASSERT(c_o);
m().reset(new_l_val);
@ -1747,10 +1748,10 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
// 0 < a <= x <= b < 0, y <= d < 0 --> b/d <= x/y
// 0 < a <= x <= b < 0, y <= d < 0 --> x/y <= a/d
TRACE("interval_bug", tout << "(M, N) #" << call_id << "\n";);
set_lower_is_open(r, b_o || d_o);
set_upper_is_open(r, a_o || d_o);
set_upper_is_open(r, a_o || d_o);
if (m().is_zero(d)) {
SASSERT(d_o);
m().reset(new_l_val); m().reset(new_u_val);
@ -1771,10 +1772,10 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
TRACE("interval_bug", tout << "(M, P) #" << call_id << "\n";);
SASSERT(is_P1(i2));
set_lower_is_open(r, a_o || c_o);
set_upper_is_open(r, b_o || c_o);
set_upper_is_open(r, b_o || c_o);
if (m().is_zero(c)) {
SASSERT(c_o);
m().reset(new_l_val); m().reset(new_u_val);
@ -1790,15 +1791,15 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
}
}
else {
SASSERT(is_P(i1));
SASSERT(is_P(i1));
if (is_N1(i2)) {
// b > 0, x <= b, c <= y <= d < 0 --> b/d <= x/y
// 0 <= a <= x, c <= y <= d < 0 --> x/y <= a/c
TRACE("interval_bug", tout << "(P, N) #" << call_id << "\n";);
set_upper_is_open(r, is_P0(i1) ? false : a_o || c_o);
set_lower_is_open(r, b_o || d_o);
if (m().is_zero(d)) {
SASSERT(d_o);
m().reset(new_l_val);
@ -1816,10 +1817,10 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
// 0 <= a <= x, 0 < c <= y <= d --> a/d <= x/y
// b > 0 x <= b, 0 < c <= y --> x/y <= b/c
TRACE("interval_bug", tout << "(P, P) #" << call_id << "\n";);
set_lower_is_open(r, is_P0(i1) ? false : a_o || d_o);
set_upper_is_open(r, b_o || c_o);
round_to_minus_inf();
::div(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
if (m().is_zero(c)) {
@ -1833,7 +1834,7 @@ void interval_manager<C>::div(interval const & i1, interval const & i2, interval
}
}
}
m().swap(lower(r), new_l_val);
m().swap(upper(r), new_u_val);
set_lower_is_inf(r, new_l_kind == EN_MINUS_INFINITY);
@ -1851,16 +1852,16 @@ void interval_manager<C>::pi_series(int x, numeral & r, bool up) {
_scoped_numeral<numeral_manager> f(m());
set_rounding(up);
m().set(r, 4, 8*x + 1);
set_rounding(!up);
set_rounding(!up);
m().set(f, 2, 8*x + 4);
set_rounding(up);
m().sub(r, f, r);
set_rounding(!up);
set_rounding(!up);
m().set(f, 1, 8*x + 5);
set_rounding(up);
m().sub(r, f, r);
set_rounding(!up);
m().set(f, 1, 8*x + 6);
m().set(f, 1, 8*x + 6);
set_rounding(up);
m().sub(r, f, r);
m().set(f, 1, 16);
@ -1870,16 +1871,16 @@ void interval_manager<C>::pi_series(int x, numeral & r, bool up) {
template<typename C>
void interval_manager<C>::pi(unsigned n, interval & r) {
// Compute an interval that contains pi using the series
// Compute an interval that contains pi using the series
// P[0] + P[1] + ... + P[n]
// where
// P[n] := 1/16^x (4/(8x + 1) - 2/(8x + 4) - 1/(8x + 5) - 1/(8x + 6))
//
//
// The size of the interval is 1/15 * 1/(16^n)
//
// Lower is P[0] + P[1] + ... + P[n]
// Upper is Lower + 1/15 * 1/(16^n)
// compute size of the resulting interval
round_to_plus_inf(); // overestimate size of the interval
_scoped_numeral<numeral_manager> len(m());
@ -1888,7 +1889,7 @@ void interval_manager<C>::pi(unsigned n, interval & r) {
m().power(len, n, len);
m().set(p, 1, 15);
m().mul(p, len, len);
// compute lower bound
numeral & l_val = m_result_lower;
m().reset(l_val);
@ -1897,7 +1898,7 @@ void interval_manager<C>::pi(unsigned n, interval & r) {
round_to_minus_inf();
m().add(l_val, p, l_val);
}
// computer upper bound
numeral & u_val = m_result_upper;
if (m().precise()) {
@ -1959,7 +1960,7 @@ void interval_manager<C>::e_series(unsigned k, bool upper, numeral & o) {
template<typename C>
void interval_manager<C>::e(unsigned k, interval & r) {
// Store in r lower and upper bounds for Euler's constant.
//
//
// The procedure uses the series
//
// V = 1 + 1/1 + 1/2! + 1/3! + ... + 1/k!
@ -1976,9 +1977,9 @@ void interval_manager<C>::e(unsigned k, interval & r) {
fact(k+1, error);
round_to_plus_inf();
m().inv(error); // error == 1/(k+1)!
m().set(aux, 4);
m().set(aux, 4);
m().mul(aux, error, error); // error == 4/(k+1)!
if (m().precise()) {
m().set(hi, lo);
m().add(hi, error, hi);

383
src/math/polynomial/algebraic_numbers.cpp
View file

File diff suppressed because it is too large Load diff

719
src/math/polynomial/polynomial.cpp
View file

File diff suppressed because it is too large Load diff

315
src/math/polynomial/upolynomial.cpp
View file

@ -8,12 +8,12 @@ Module Name:
Abstract:
Goodies for creating and handling univariate polynomials.
A dense representation is much better for Root isolation algorithms,
A dense representation is much better for Root isolation algorithms,
encoding algebraic numbers, factorization, etc.
We also use integers as the coefficients of univariate polynomials.
Author:
Leonardo (leonardo) 2011-11-29
@ -26,21 +26,22 @@ Notes:
#include"polynomial_primes.h"
#include"buffer.h"
#include"cooperate.h"
#include"common_msgs.h"
namespace upolynomial {
core_manager::factors::factors(core_manager & upm):
m_upm(upm),
m_total_factors(0),
m_total_degree(0) {
nm().set(m_constant, 1);
core_manager::factors::factors(core_manager & upm):
m_upm(upm),
m_total_factors(0),
m_total_degree(0) {
nm().set(m_constant, 1);
}
core_manager::factors::~factors() {
clear();
nm().del(m_constant);
}
void core_manager::factors::clear() {
for (unsigned i = 0; i < m_factors.size(); ++ i) {
m_upm.reset(m_factors[i]);
@ -104,19 +105,19 @@ namespace upolynomial {
}
}
numeral_manager & core_manager::factors::nm() const {
return m_upm.m();
}
numeral_manager & core_manager::factors::nm() const {
return m_upm.m();
}
void core_manager::factors::set_constant(numeral const & constant) {
nm().set(m_constant, constant);
}
void core_manager::factors::set_constant(numeral const & constant) {
nm().set(m_constant, constant);
}
void core_manager::factors::set_degree(unsigned i, unsigned degree) {
void core_manager::factors::set_degree(unsigned i, unsigned degree) {
m_total_degree -= m_upm.degree(m_factors[i])*m_degrees[i];
m_total_factors -= m_degrees[i];
m_degrees[i] = degree;
m_total_factors += degree;
m_total_factors -= m_degrees[i];
m_degrees[i] = degree;
m_total_factors += degree;
m_total_degree += m_upm.degree(m_factors[i])*degree;
}
@ -125,7 +126,7 @@ namespace upolynomial {
m_total_degree += m_upm.degree(p)*m_degrees[i];
m_factors[i].swap(p);
}
void core_manager::factors::swap(factors & other) {
m_factors.swap(other.m_factors);
m_degrees.swap(other.m_degrees);
@ -154,8 +155,8 @@ namespace upolynomial {
}
void core_manager::checkpoint() {
if (!m_limit.inc())
throw upolynomial_exception("canceled");
if (!m_limit.inc())
throw upolynomial_exception(Z3_CANCELED_MSG);
cooperate("upolynomial");
}
@ -209,7 +210,7 @@ namespace upolynomial {
set_size(sz, buffer);
}
void core_manager::get_primitive_and_content(unsigned f_sz, numeral const * f, numeral_vector & pp, numeral & cont) {
void core_manager::get_primitive_and_content(unsigned f_sz, numeral const * f, numeral_vector & pp, numeral & cont) {
SASSERT(f_sz > 0);
m().gcd(f_sz, f, cont);
SASSERT(m().is_pos(cont));
@ -221,7 +222,7 @@ namespace upolynomial {
for (unsigned i = 0; i < f_sz; i++) {
if (!m().is_zero(f[i])) {
m().div(f[i], cont, pp[i]);
}
}
else {
m().set(pp[i], 0);
}
@ -252,7 +253,7 @@ namespace upolynomial {
neg_core(sz, p, m_basic_tmp);
buffer.swap(m_basic_tmp);
}
// buffer := p1 + p2
void core_manager::add_core(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & buffer) {
SASSERT(!is_alias(p1, buffer));
@ -334,7 +335,7 @@ namespace upolynomial {
numeral const & b_j = p2[j];
if (m().is_zero(b_j))
continue;
m().addmul(buffer[i+j], a_i, b_j, buffer[i+j]);
m().addmul(buffer[i+j], a_i, b_j, buffer[i+j]);
}
}
set_size(new_sz, buffer);
@ -378,7 +379,7 @@ namespace upolynomial {
return;
for (unsigned i = 0; i < sz; i++) {
#ifdef Z3DEBUG
scoped_numeral old_p_i(m());
scoped_numeral old_p_i(m());
old_p_i = p[i];
#endif
// Actual code
@ -387,7 +388,7 @@ namespace upolynomial {
#ifdef Z3DEBUG
scoped_numeral tmp(m());
m().mul(g, p[i], tmp);
CTRACE("div_bug", !m().eq(tmp, old_p_i), tout << "old(p[i]): " << m().to_string(old_p_i) << ", g: " << m().to_string(g) << ", p[i]: " <<
CTRACE("div_bug", !m().eq(tmp, old_p_i), tout << "old(p[i]): " << m().to_string(old_p_i) << ", g: " << m().to_string(g) << ", p[i]: " <<
m().to_string(p[i]) << ", tmp: " << m().to_string(tmp) << "\n";);
SASSERT(tmp == old_p_i);
#endif
@ -428,7 +429,7 @@ namespace upolynomial {
}
// Pseudo division
void core_manager::div_rem_core(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2,
void core_manager::div_rem_core(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2,
unsigned & d, numeral_vector & q, numeral_vector & r) {
SASSERT(!is_alias(p1, q)); SASSERT(!is_alias(p2, q));
SASSERT(!is_alias(p1, r)); SASSERT(!is_alias(p2, r));
@ -438,7 +439,7 @@ namespace upolynomial {
set(sz1, p1, q);
if (field()) {
div(q, *p2);
}
}
reset(r);
return;
}
@ -464,7 +465,7 @@ namespace upolynomial {
set_size(qsz, q);
return;
}
unsigned m_n = sz1 - sz2; // m-n
unsigned m_n = sz1 - sz2; // m-n
if (field()) {
numeral & ratio = a_m;
m().div(r[sz1 - 1], b_n, ratio);
@ -492,7 +493,7 @@ namespace upolynomial {
}
// Pseudo division
void core_manager::div_rem(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, unsigned & d,
void core_manager::div_rem(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, unsigned & d,
numeral_vector & q, numeral_vector & r) {
numeral_vector & _r = m_div_tmp1;
numeral_vector & _q = m_div_tmp2;
@ -677,15 +678,15 @@ namespace upolynomial {
scoped_numeral a1(m());
scoped_numeral a2(m());
m().mul(b2, inv2, a1); // a1 is the multiplicator for coefficients of C1
m().mul(b1, inv1, a2); // a2 is the multiplicator for coefficients of C2
m().mul(b1, inv1, a2); // a2 is the multiplicator for coefficients of C2
TRACE("CRA", tout << "a1: " << a1 << ", a2: " << a2 << "\n";);
// new bound
scoped_numeral new_bound(m());
m().mul(b1, b2, new_bound);
m().mul(b1, b2, new_bound);
scoped_numeral lower(m());
scoped_numeral upper(m());
scoped_numeral new_a(m()), tmp1(m()), tmp2(m()), tmp3(m());
m().div(new_bound, 2, upper);
m().div(new_bound, 2, upper);
m().set(lower, upper);
m().neg(lower);
TRACE("CRA", tout << "lower: " << lower << ", upper: " << upper << "\n";);
@ -700,7 +701,7 @@ namespace upolynomial {
R.push_back(numeral()); \
m().set(R.back(), new_a); \
}
numeral zero(0);
unsigned i = 0;
unsigned sz1 = C1.size();
@ -718,7 +719,7 @@ namespace upolynomial {
m().set(b2, new_bound);
R.swap(C2);
}
void core_manager::mod_gcd(unsigned sz_u, numeral const * u,
unsigned sz_v, numeral const * v,
numeral_vector & result) {
@ -726,8 +727,8 @@ namespace upolynomial {
SASSERT(sz_u > 0 && sz_v > 0);
SASSERT(!m().modular());
scoped_numeral c_u(m()), c_v(m());
numeral_vector & pp_u = m_mgcd_tmp[0];
numeral_vector & pp_v = m_mgcd_tmp[1];
numeral_vector & pp_u = m_mgcd_tmp[0];
numeral_vector & pp_v = m_mgcd_tmp[1];
get_primitive_and_content(sz_u, u, pp_u, c_u);
get_primitive_and_content(sz_v, v, pp_v, c_v);
scoped_numeral c_g(m());
@ -745,7 +746,7 @@ namespace upolynomial {
numeral_vector & v_Zp = m_mgcd_tmp[3];
numeral_vector & q = m_mgcd_tmp[4];
numeral_vector & C = m_mgcd_tmp[5];
for (unsigned i = 0; i < NUM_BIG_PRIMES; i++) {
m().set(p, polynomial::g_big_primes[i]);
TRACE("mgcd", tout << "trying prime: " << p << "\n";);
@ -797,8 +798,8 @@ namespace upolynomial {
TRACE("mgcd", tout << "candidate:\n"; display_star(tout, candidate); tout << "\n";);
SASSERT(candidate.size() > 0);
numeral const & lc_candidate = candidate[candidate.size() - 1];
if (m().divides(lc_candidate, lc_g) &&
divides(pp_u, candidate) &&
if (m().divides(lc_candidate, lc_g) &&
divides(pp_u, candidate) &&
divides(pp_v, candidate)) {
TRACE("mgcd", tout << "found GCD\n";);
mul(candidate, c_g);
@ -845,9 +846,9 @@ namespace upolynomial {
else {
flip_sign_if_lm_neg(buffer);
}
TRACE("upolynomial", tout << "GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
TRACE("upolynomial", tout << "GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
display(tout, buffer); tout << "\n";);
return;
return;
}
rem(A.size(), A.c_ptr(), B.size(), B.c_ptr(), R);
normalize(R);
@ -868,7 +869,7 @@ namespace upolynomial {
flip_sign_if_lm_neg(buffer);
return;
}
if (!modular())
if (!modular())
mod_gcd(sz1, p1, sz2, p2, buffer);
else
euclid_gcd(sz1, p1, sz2, p2, buffer);
@ -913,9 +914,9 @@ namespace upolynomial {
else {
flip_sign_if_lm_neg(buffer);
}
TRACE("upolynomial", tout << "subresultant GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
TRACE("upolynomial", tout << "subresultant GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
display(tout, buffer); tout << "\n";);
return;
return;
}
rem(A.size(), A.c_ptr(), B.size(), B.c_ptr(), d, R);
unsigned pseudo_div_d = A.size() - B.size();
@ -928,7 +929,7 @@ namespace upolynomial {
m().power(B[B.size() - 1], pseudo_div_d + 1 - d, aux);
mul(R, aux);
}
d = pseudo_div_d;
d = pseudo_div_d;
TRACE("upolynomial", tout << "R: "; display(tout, R); tout << "\nd: " << d << "\n";);
// aux <- g*h^d
m().power(h, d, aux);
@ -961,7 +962,7 @@ namespace upolynomial {
// buffer := square-free(p)
void core_manager::square_free(unsigned sz, numeral const * p, numeral_vector & buffer) {
SASSERT(!is_alias(p, buffer));
SASSERT(!is_alias(p, buffer));
if (sz <= 1) {
set(sz, p, buffer);
return;
@ -1000,11 +1001,11 @@ namespace upolynomial {
return;
}
if (k == 1 || sz == 0 || (sz == 1 && m().is_one(p[0]))) {
if (k == 1 || sz == 0 || (sz == 1 && m().is_one(p[0]))) {
set(sz, p, r);
return;
}
numeral_vector & result = m_pw_tmp;
set(sz, p, result);
for (unsigned i = 1; i < k; i++)
@ -1039,18 +1040,18 @@ namespace upolynomial {
m().mul(p[d], lc_inv, p[d]);
}
}
}
}
// Extended GCD
void core_manager::ext_gcd(unsigned szA, numeral const * A, unsigned szB, numeral const * B,
void core_manager::ext_gcd(unsigned szA, numeral const * A, unsigned szB, numeral const * B,
numeral_vector & U, numeral_vector & V, numeral_vector & D) {
SASSERT(!is_alias(A, U)); SASSERT(!is_alias(A, V)); SASSERT(!is_alias(A, D));
SASSERT(!is_alias(B, U)); SASSERT(!is_alias(B, V)); SASSERT(!is_alias(B, D));
SASSERT(field());
scoped_numeral_vector V1(m()), V3(m()), Q(m()), R(m()), T(m()), V1Q(m());
// since we are in a field define gcd(A, B) to be monic
// if AU + BV = D and D is not monic we make it monic, and then divide U and V by the same inverse
@ -1058,13 +1059,13 @@ namespace upolynomial {
// U <- 1
reset(U); U.push_back(numeral()); m().set(U.back(), 1);
// D <- A
set(szA, A, D);
set(szA, A, D);
mk_monic(szA, D.c_ptr());
// V1 <- 0
reset(V1);
reset(V1);
// V3 <- B
set(szB, B, V3);
set(szB, B, V3);
while (true) {
if (V3.empty()) {
// V3 is the zero polynomial
@ -1088,10 +1089,10 @@ namespace upolynomial {
mul(V, lc_inv);
return;
}
// D = QV3 + R
div_rem(D.size(), D.c_ptr(), V3.size(), V3.c_ptr(), Q, R);
// T <- U - V1Q
mul(V1.size(), V1.c_ptr(), Q.size(), Q.c_ptr(), V1Q);
sub(U.size(), U.c_ptr(), V1Q.size(), V1Q.c_ptr(), T);
@ -1104,8 +1105,8 @@ namespace upolynomial {
// V3 <- R
V3.swap(R);
}
}
}
// Display p
void core_manager::display(std::ostream & out, unsigned sz, numeral const * p, char const * var_name, bool use_star) const {
bool displayed = false;
@ -1144,7 +1145,7 @@ namespace upolynomial {
if (!displayed)
out << "0";
}
static void display_smt2_mumeral(std::ostream & out, numeral_manager & m, mpz const & n) {
if (m.is_neg(n)) {
out << "(- ";
@ -1172,7 +1173,7 @@ namespace upolynomial {
out << "(^ " << var_name << " " << k << ")";
}
else {
out << "(* ";
out << "(* ";
display_smt2_mumeral(out, m, n);
out << " ";
if (k == 1)
@ -1189,12 +1190,12 @@ namespace upolynomial {
out << "0";
return;
}
if (sz == 1) {
display_smt2_mumeral(out, m(), p[0]);
return;
}
unsigned non_zero_idx = UINT_MAX;
unsigned num_non_zeros = 0;
for (unsigned i = 0; i < sz; i++) {
@ -1208,7 +1209,7 @@ namespace upolynomial {
SASSERT(non_zero_idx != UINT_MAX && non_zero_idx >= 1);
display_smt2_monomial(out, m(), p[non_zero_idx], non_zero_idx, var_name);
}
out << "(+";
unsigned i = sz;
while (i > 0) {
@ -1230,7 +1231,7 @@ namespace upolynomial {
}
return true;
}
void upolynomial_sequence::push(unsigned sz, numeral * p) {
m_begins.push_back(m_seq_coeffs.size());
m_szs.push_back(sz);
@ -1253,8 +1254,8 @@ namespace upolynomial {
_scoped_numeral_vector<numeral_manager>(m.m()) {
}
scoped_upolynomial_sequence::~scoped_upolynomial_sequence() {
m_manager.reset(*this);
scoped_upolynomial_sequence::~scoped_upolynomial_sequence() {
m_manager.reset(*this);
}
manager::~manager() {
@ -1355,7 +1356,7 @@ namespace upolynomial {
return r;
}
// Return the descartes bound for (0, +oo)
unsigned manager::descartes_bound(unsigned sz, numeral const * p) {
return sign_changes(sz, p);
@ -1382,9 +1383,9 @@ namespace upolynomial {
unsigned num_vars = 0;
// a0 a1 a2 a3
// a0 a0+a1 a0+a1+a2 a0+a1+a2+a3
// a0 2a0+a1 3a0+2a1+a2
// a0 3a0+a1
// a0
// a0 2a0+a1 3a0+2a1+a2
// a0 3a0+a1
// a0
for (unsigned i = 0; i < sz; i++) {
checkpoint();
unsigned k;
@ -1496,12 +1497,12 @@ namespace upolynomial {
}
// Given b of the form c/(2^k)
// p(x) := (2^k)^n * p(x + c/(2^k)) where b = c/2^k
// p(x) := (2^k)^n * p(x + c/(2^k)) where b = c/2^k
//
// Given p: a_n * x^n + ... + a_0
//
//
// buffer := a_n * (2^k * x + c)^n + a_{n-1} * 2^k * (2^k * x + c)^{n-1} + ... + a_1 * (2^k)^{n-1} * (2^k * x + c) + a_0 * (2^k)^n
//
//
// We perform the transformation in two steps:
// 1) a_n * x^n + a_{n-1} * 2^k * x^{n-1} + ... + a_1 * (2^k)^{n-1} + a_0 * (2^k)^n
// Let d_n, ..., d_0 be the coefficients after this step
@ -1509,33 +1510,33 @@ namespace upolynomial {
// d_n * x^n + ... + d_0
// 2) d_n * (2^k * x + c)^n + ... + d_0
// This step is like the translation with integers, but the coefficients must be adjusted by 2^k in each iteration.
//
//
// That is, it is a special translation such as:
// a_n*(b*x + c)^n + a_{n-1}*(b*x + c)^{n-1} + ... + a_1*(b*x + c) + a_0
// This kind of translation can be computed by applying the following recurrence:
// To simplify the notation, we fix n = 4
// Moreover we use a_i[k] to represent the coefficient a_i at iteration k
// a_i[0] = a_i
//
//
// a_4*(b*x + c)^4 + a_3*(b*x + c)^3 + a_2*(b*x + c)^2 + a_1*(b*x + c) + a_0
// -->
// a_4*b*x*(b*x + c)^3 + (a_3 + a_4*c)*(b*x + c)^3 + a_2*(b*x + c)^2 + a_1*(b*x + c) + a_0
// Thus a_4[1] = a_4[0]*b, a_3[1] = a_3[0] + a_4[0]*c, a_2[1] = a_2[0], a_1[1] = a_1[0], a_0[1] = a_0[0]
//
//
// a_4[1]*x*(b*x + c)^3 + a_3[1]*(b*x + c)^3 + a_2[1]*(b*x + c)^2 + a_1[1]*(b*x + c) + a_0[1]
// -->
// a_4[1]*b*x^2*(b*x + c)^2 + (a_3[1]*b + a_4[1]*c)*x*(b*x + c)^2 + (a_2[1] + a_3[1]*c)*(b*x + c)^2 + a_1[1]*(b*x + c) + a_0[1]
// Thus a_4[2] = a_4[1]*b, a_3[2] = a_3[1]*b + a_4[1]*c, a_2[2] = a_2[1] + a_3[1]*c, a_1[2] = a_1[1], a_0[2] = a_0[1]
//
//
// a_4[2]*x^2*(b*x + c)^2 + a_3[2]*x*(b*x + c)^2 + a_2[2]*(b*x + c)^2 + a_1[2]*(b*x + c) + a_0[2]
// -->
// a_4[2]*b*x^3*(b*x + c) + (a_3[2]*b + a_4[2]*c)*x^2*(b*x + c) + (a_2[2]*b + a_3[2]*c)*x*(b*x + c) + (a_1[2] + a_2[2]*c)*(b*x + c) + a_0[2]
// Thus a_4[3] = a_4[2]*b, a_3[3] = a_3[2]*b + a_4[2]*c, a_2[3] = a_2[2]*b + a_3[2]*c, a_1[3] = a_1[2] + a_2[2]*c, a_0[3] = a_1[3]
//
//
// a_4[3]*x^3*(b*x + c) + a_3[3]*x^2*(b*x + c) + a_2[3]*x*(b*x + c) + a_1[3]*(b*x + c) + a_0[3]
// --->
// a_4[3]*b*x^4 + (a_3[3]*b + a_4[3]*c)*x^3 + (a_2[3]*b + a_3[3]*c)*x^2 + (a_1[3]*b + a_2[3]*c)*x + (a_0[3] + a_1[3]*c)
//
// a_4[3]*b*x^4 + (a_3[3]*b + a_4[3]*c)*x^3 + (a_2[3]*b + a_3[3]*c)*x^2 + (a_1[3]*b + a_2[3]*c)*x + (a_0[3] + a_1[3]*c)
//
void manager::translate_bq(unsigned sz, numeral * p, mpbq const & b) {
if (sz <= 1)
return;
@ -1617,7 +1618,7 @@ namespace upolynomial {
}
}
// p(x) := p(2^k * x)
// p(x) := p(2^k * x)
void manager::compose_p_2k_x(unsigned sz, numeral * p, unsigned k) {
// (2^k)^n*a_n*x^n + (2^k)^(n-1)*x^{n-1} + ... + a_0
if (sz <= 1)
@ -1629,11 +1630,11 @@ namespace upolynomial {
}
}
// p(x) := (2^k)^n * p(x/2^k)
// p(x) := (2^k)^n * p(x/2^k)
//
// Let old(p) be of the form:
// a_n * x^n + a_{n-1}*x^{n-1} + ... + a_1 * x + a_0
//
//
// Then p is of the form:
// a_n * x^n + a_{n-1} * 2^k * x^{n-1} + a_{n-2} * (2^k)^2 * x^{n-2} + ... + a_0 * (2^k)^n
void manager::compose_2kn_p_x_div_2k(unsigned sz, numeral * p, unsigned k) {
@ -1646,7 +1647,7 @@ namespace upolynomial {
m().mul2k(p[i], k_i);
}
}
// p(x) := a^n * p(x/a)
// See compose_2kn_p_x_div_2k
void manager::compose_an_p_x_div_a(unsigned sz, numeral * p, numeral const & a) {
@ -1675,13 +1676,13 @@ namespace upolynomial {
m().mul(b_i, b, b_i);
}
}
// Let b be of the form c/(2^k), then this operation is equivalent to:
// (2^k)^n*p(c*x/(2^k))
//
//
// Let old(p) be of the form:
// a_n * x^n + a_{n-1}*x^{n-1} + ... + a_1 * x + a_0
//
//
// Then p is of the form:
// a_n * c^n * x^n + a_{n-1} * c^{n-1} * 2^k * x^{n-1} + ... + a_1 * c * (2^k)^(n-1) * x + a_0 * (2^k)^n
void manager::compose_p_b_x(unsigned sz, numeral * p, mpbq const & b) {
@ -1705,12 +1706,12 @@ namespace upolynomial {
// Let q be of the form b/c, then this operation is equivalent to:
// p(x) := c^n*p(b*x/c)
//
//
// If u is a root of old(p), then u*(c/b) is a root of new p.
//
//
// Let old(p) be of the form:
// a_n * x^n + a_{n-1}*x^{n-1} + ... + a_1 * x + a_0
//
//
// Then p is of the form:
// a_n * b^n * x^n + a_{n-1} * b^{n-1} * c * x^{n-1} + ... + a_1 * b * c^(n-1) * x + a_0 * c^n
void manager::compose_p_q_x(unsigned sz, numeral * p, mpq const & q) {
@ -1730,14 +1731,14 @@ namespace upolynomial {
}
}
}
// Evaluate the sign of p(b)
int manager::eval_sign_at(unsigned sz, numeral const * p, mpbq const & b) {
// Actually, given b = c/2^k, we compute the sign of (2^k)^n*p(b)
// Original Horner Sequence
// ((a_n * b + a_{n-1})*b + a_{n-2})*b + a_{n-3} ...
// Variation of the Horner Sequence for (2^k)^n*p(b)
// ((a_n * c + a_{n-1}*2_k)*c + a_{n-2}*(2_k)^2)*c + a_{n-3}*(2_k)^3 ... + a_0*(2_k)^n
// ((a_n * c + a_{n-1}*2_k)*c + a_{n-2}*(2_k)^2)*c + a_{n-3}*(2_k)^3 ... + a_0*(2_k)^n
if (sz == 0)
return 0;
if (sz == 1)
@ -1770,7 +1771,7 @@ namespace upolynomial {
// Original Horner Sequence
// ((a_n * b + a_{n-1})*b + a_{n-2})*b + a_{n-3} ...
// Variation of the Horner Sequence for (d^n)*p(b)
// ((a_n * c + a_{n-1}*d)*c + a_{n-2}*d^2)*c + a_{n-3}*d^3 ... + a_0*d^n
// ((a_n * c + a_{n-1}*d)*c + a_{n-2}*d^2)*c + a_{n-3}*d^3 ... + a_0*d^n
if (sz == 0)
return 0;
if (sz == 1)
@ -1797,7 +1798,7 @@ namespace upolynomial {
}
return sign_of(r);
}
// Evaluate the sign of p(b)
int manager::eval_sign_at(unsigned sz, numeral const * p, mpz const & b) {
// Using Horner Sequence
@ -1906,7 +1907,7 @@ namespace upolynomial {
void manager::root_upper_bound(unsigned sz, numeral const * p, numeral & r) {
// Using Cauchy's Inequality
// We have that any root u of p must satisfy
// We have that any root u of p must satisfy
// |u| < (max(p) + min(p))/min(p)
// |u| < (max(p) + |a_n|)/|a_n|
// where: max(p) is the maximum coefficient in absolute value.
@ -1932,7 +1933,7 @@ namespace upolynomial {
init = true;
continue;
}
if (m().gt(c, max))
if (m().gt(c, max))
m().set(max, c);
if (m().lt(c, min))
m().set(min, c);
@ -1946,26 +1947,26 @@ namespace upolynomial {
m().div(r2, a_n, r2);
m().add(r2, numeral(1), r2);
// use the best bound
if (m().lt(r2, r))
if (m().lt(r2, r))
swap(r, r2);
SASSERT(m().le(r, r2));
}
/**
\brief Find positive root upper bound using Knuth's approach.
Given p(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_0
If a_n is positive,
Let B = max({ (-a_{n-k}/a_n)^{1/k} | 1 <= k <= n, a_{n-k} < 0 })
Then, 2*B is a bound for the positive roots
Similarly, if a_n is negative
Let B = max({ (-a_{n-k}/a_n)^{1/k} | 1 <= k <= n, a_{n-k} > 0 })
Then, 2*B is a bound for the positive roots
This procedure returns a k s.t. 2*B <= 2^k
The procedure actually computes
max of log2(abs(a_{n-k})) + 1 - log2(abs(a_n))/k
@ -1974,12 +1975,12 @@ namespace upolynomial {
Let u > 0 be a root of p(x).
If u <= B, then we are done. So, let us assume u > B
Assume, a_n > 0 (the proof is similar for a_n < 0)
Then: a_n * u^n + a_{n-1} * u^{n-1} + ... + a0 = 0
u^n = 1/a_n (-a_{n-1} * u^{n-1} - ... - a_0)
Note that, if we remove the monomials s.t. a_i is positive, then
Note that, if we remove the monomials s.t. a_i is positive, then
we are increasing the value of (-a_{n-1} * u^{n-1} - ... - a_0).
Thus,
u^n <= 1/a_n(SUM_{a_i < 0, 0 <= i < n} (-a_i * u^i))
@ -1990,7 +1991,7 @@ namespace upolynomial {
1 <= 1/a_n(SUM_{a_{n-k} < 0, n >= k > 0} (-a_{n-k} * u^{n-k})) = 1/a_n(SUM_{a_{n-k} < 0, 1 <= k <= n} (-a_i * u^{n-k})) < Sum_{1 <= k < +oo}(B/u)^k
Since u > B, we have that Sum_{1 <= k < +oo}(B/u)^k = B/(u - B)
Thus, we have
1 < B/(u - B)
1 < B/(u - B)
and u < 2B
*/
unsigned manager::knuth_positive_root_upper_bound(unsigned sz, numeral const * p) {
@ -2044,14 +2045,14 @@ namespace upolynomial {
We essentially compute the upper bound for the roots of x^n*p(1/x) where n is the degree of p.
Remark: alpha is a nonzero root of p(x) iff 1/alpha is a root of x^n*p(1/x).
Thus, if 2^k is upper bound for the root of x^n*p(1/x). Then we have that
-2^k < 1/alpha < 2^k
-2^k < 1/alpha < 2^k
and consequently
alpha < -1/2^k or 1/2^k < alpha
/pre p is not the zero polynomial.
*/
unsigned manager::nonzero_root_lower_bound(unsigned sz, numeral const * p) {
SASSERT(sz > 0);
SASSERT(sz > 0);
// consume zeros
unsigned i = 0;
while (true) {
@ -2082,7 +2083,7 @@ namespace upolynomial {
struct manager::drs_frame {
unsigned m_parent_idx; // position of the parent frame, UINT_MAX if it doesn't have a parent frame
unsigned m_size:30; // size of the polynomial associated with this frame
unsigned m_first:1; // first time visiting the frame?
unsigned m_first:1; // first time visiting the frame?
unsigned m_left:1; // true if the frame is the left child of the parent frame.
drs_frame(unsigned pidx, unsigned sz, bool left):
m_parent_idx(pidx),
@ -2115,14 +2116,14 @@ namespace upolynomial {
// I don't really need the following test, because 0 - 1 == UINT_MAX
unsigned parent_idx = frame_stack.empty() ? UINT_MAX : frame_stack.size() - 1;
numeral_vector & p_aux = m_push_tmp;
// Possible optimization: the coefficients of the parent frame are not needed anymore.
// So, we could reuse/steal them.
// left child
set(sz, p, p_aux);
compose_2n_p_x_div_2(p_aux.size(), p_aux.c_ptr());
normalize(p_aux);
normalize(p_aux);
for (unsigned i = 0; i < sz; i++) {
p_stack.push_back(numeral());
m().set(p_stack.back(), p_aux[i]);
@ -2149,7 +2150,7 @@ namespace upolynomial {
unsigned idx = frame_stack.size() - 1;
while (idx != UINT_MAX) {
drs_frame const & fr = frame_stack[idx];
TRACE("upolynomial",
TRACE("upolynomial",
tout << "normalizing...\n";
tout << "idx: " << idx << ", left: " << fr.m_left << ", l: " << bqm.to_string(l) << ", u: " << bqm.to_string(u) << "\n";);
if (fr.m_left) {
@ -2170,7 +2171,7 @@ namespace upolynomial {
swap(lowers.back(), l);
swap(uppers.back(), u);
}
// 1/2 is a root of the polynomial associated with the top frame.
// Apply transformations for obtaining root of the original polynomial:
// We use the following transformations:
@ -2259,7 +2260,7 @@ namespace upolynomial {
push_child_frames(q.size(), q.c_ptr(), p_stack, frame_stack);
}
else {
push_child_frames(sz, p, p_stack, frame_stack);
push_child_frames(sz, p, p_stack, frame_stack);
}
}
}
@ -2314,7 +2315,7 @@ namespace upolynomial {
TRACE("upolynomial", tout << "searching at (-1, 0) using:\n"; display(tout, sz, p); tout << "\n";);
old_roots_sz = roots.size();
old_lowers_sz = lowers.size();
drs_isolate_0_1_roots(sz, p, bqm, roots, lowers, uppers);
drs_isolate_0_1_roots(sz, p, bqm, roots, lowers, uppers);
SASSERT(lowers.size() == uppers.size());
adjust_neg(bqm, roots, old_roots_sz, neg_k);
adjust_neg(bqm, lowers, old_lowers_sz, neg_k);
@ -2390,7 +2391,7 @@ namespace upolynomial {
}
// Isolate roots in an interval (-2^neg_k, 2^pos_k) using an approach based on Descartes rule of signs.
void manager::sturm_isolate_roots_core(unsigned sz, numeral * p, unsigned neg_k, unsigned pos_k,
void manager::sturm_isolate_roots_core(unsigned sz, numeral * p, unsigned neg_k, unsigned pos_k,
mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers) {
SASSERT(!has_zero_roots(sz, p));
scoped_upolynomial_sequence seq(*this);
@ -2450,10 +2451,10 @@ namespace upolynomial {
bqm.swap(curr_lower, f.m_lower);
bqm.swap(curr_upper, f.m_upper);
pop_ss_frame(bqm, s);
SASSERT(lower_sv > upper_sv + 1);
SASSERT(lower_sv > upper_sv + 1);
bqm.add(curr_lower, curr_upper, mid);
bqm.div2(mid);
TRACE("upolynomial",
TRACE("upolynomial",
tout << "depth: " << s.size() << "\n";
tout << "lower_sv: " << lower_sv << "\n";
tout << "upper_sv: " << upper_sv << "\n";
@ -2492,8 +2493,8 @@ namespace upolynomial {
}
void manager::sqf_isolate_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers) {
bqm.reset(roots);
bqm.reset(lowers);
bqm.reset(roots);
bqm.reset(lowers);
bqm.reset(uppers);
if (has_zero_roots(sz, p)) {
roots.push_back(mpbq(0));
@ -2516,10 +2517,10 @@ namespace upolynomial {
TRACE("upolynomial", tout << "square free part:\n"; display(tout, sqf_p); tout << "\n";);
sqf_isolate_roots(sqf_p.size(), sqf_p.c_ptr(), bqm, roots, lowers, uppers);
}
// Keep expanding the Sturm sequence starting at seq
void manager::sturm_seq_core(upolynomial_sequence & seq) {
scoped_numeral_vector r(m());
scoped_numeral_vector r(m());
while (true) {
unsigned sz = seq.size();
srem(seq.size(sz-2), seq.coeffs(sz-2), seq.size(sz-1), seq.coeffs(sz-1), r);
@ -2539,7 +2540,7 @@ namespace upolynomial {
void manager::sturm_seq(unsigned sz, numeral const * p, upolynomial_sequence & seq) {
reset(seq);
scoped_numeral_vector p_prime(m());
scoped_numeral_vector p_prime(m());
seq.push(m(), sz, p);
derivative(sz, p, p_prime);
seq.push(p_prime.size(), p_prime.c_ptr());
@ -2548,7 +2549,7 @@ namespace upolynomial {
void manager::sturm_tarski_seq(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, upolynomial_sequence & seq) {
reset(seq);
scoped_numeral_vector p1p2(m());
scoped_numeral_vector p1p2(m());
seq.push(m(), sz1, p1);
derivative(sz1, p1, p1p2);
mul(p1p2.size(), p1p2.c_ptr(), sz2, p2, p1p2);
@ -2558,7 +2559,7 @@ namespace upolynomial {
void manager::fourier_seq(unsigned sz, numeral const * p, upolynomial_sequence & seq) {
reset(seq);
scoped_numeral_vector p_prime(m());
scoped_numeral_vector p_prime(m());
seq.push(m(), sz, p);
if (sz == 0)
return;
@ -2570,16 +2571,16 @@ namespace upolynomial {
seq.push(p_prime.size(), p_prime.c_ptr());
}
}
/**
We say an interval (a, b) of a polynomial p is ISOLATING if p has only one root in the
We say an interval (a, b) of a polynomial p is ISOLATING if p has only one root in the
interval (a, b).
We say an isolating interval (a, b) of a square free polynomial p is REFINEABLE if
sign(p(a)) = -sign(p(b))
Not every isolating interval (a, b) of a square free polynomial p is refineable, because
sign(p(a)) or sign(p(b)) may be zero.
sign(p(a)) or sign(p(b)) may be zero.
Refinable intervals of square free polynomials are useful, because we can increase precision
("squeeze" the interval) by just evaluating p at (a+b)/2
@ -2590,7 +2591,7 @@ namespace upolynomial {
The method returns TRUE if it produced a REFINABLE interval (a', b'). The new interval
is stored in input variables a and b.
The method returns FALSE if it found the root a' in the interval (a, b). The root is
The method returns FALSE if it found the root a' in the interval (a, b). The root is
stored in the input variable a.
\pre p MUST BE SQUARE FREE.
@ -2631,7 +2632,7 @@ namespace upolynomial {
}
}
}
if (sign_a != 0 && sign_b == 0 ) {
if (sign_a != 0 && sign_b == 0 ) {
// CASE 3
scoped_mpbq new_b(bqm);
// new_b <- (a+b)/2
@ -2671,7 +2672,7 @@ namespace upolynomial {
// b1 <- (a+b)/2
bqm.add(a, b, b1);
bqm.div2(b1);
// a2 <- (a+b)/2
// a2 <- (a+b)/2
bqm.set(a2, b1);
int sign_b1 = eval_sign_at(sz, p, b1);
int sign_a2 = sign_b1;
@ -2691,13 +2692,13 @@ namespace upolynomial {
bqm.div2(new_b2);
while (true) {
TRACE("upolynomial",
TRACE("upolynomial",
tout << "CASE 4\na1: " << bqm.to_string(a1) << ", b1: " << bqm.to_string(b1) << ", new_a1: " << bqm.to_string(new_a1) << "\n";
tout << "a2: " << bqm.to_string(a2) << ", b2: " << bqm.to_string(b2) << ", new_b2: " << bqm.to_string(new_b2) << "\n";);
int sign_new_a1 = eval_sign_at(sz, p, new_a1);
if (sign_new_a1 == 0) {
// found root
swap(new_a1, a);
swap(new_a1, a);
result = false;
goto end;
}
@ -2716,7 +2717,7 @@ namespace upolynomial {
result = false;
goto end;
}
if (sign_new_b2 == -sign_a2) {
// found interval
swap(a2, a);
@ -2739,7 +2740,7 @@ namespace upolynomial {
bqm.add(b2, a2, new_b2);
bqm.div2(new_b2);
}
end:
return result;
}
@ -2748,11 +2749,11 @@ namespace upolynomial {
Given a square-free polynomial p, and a refinable interval (a,b), then "squeeze" (a,b).
That is, return a new interval (a',b') s.t. b' - a' = (b - a)/2
See isolating2refinable for a definition of refinable interval.
Return TRUE, if interval was squeezed, and new interval is stored in (a,b).
Return FALSE, if the actual root was found, it is stored in a.
The arguments sign_a and sign_b must contain the values returned by
The arguments sign_a and sign_b must contain the values returned by
eval_sign_at(sz, p, a) and eval_sign_at(sz, p, b).
*/
bool manager::refine_core(unsigned sz, numeral const * p, int sign_a, mpbq_manager & bqm, mpbq & a, mpbq & b) {
@ -2786,7 +2787,7 @@ namespace upolynomial {
}
// Keeps reducing the interval until b - a < 1/2^k or a root is found.
// See comments in refine_core above.
// See comments in refine_core above.
//
// Return TRUE, if interval was squeezed, and new interval is stored in (a,b).
// Return FALSE, if the actual root was found, it is stored in a.
@ -2821,7 +2822,7 @@ namespace upolynomial {
SASSERT(sign_a != 0 && sign_b != 0);
SASSERT(sign_a == -sign_b);
bool found_d = false;
TRACE("convert_bug",
TRACE("convert_bug",
tout << "a: " << m().to_string(a.numerator()) << "/" << m().to_string(a.denominator()) << "\n";
tout << "b: " << m().to_string(b.numerator()) << "/" << m().to_string(b.denominator()) << "\n";
tout << "sign_a: " << sign_a << "\n";
@ -2839,7 +2840,7 @@ namespace upolynomial {
bqm.mul2(upper);
if (m_manager.is_neg(a.numerator()))
::swap(lower, upper);
TRACE("convert_bug",
TRACE("convert_bug",
tout << "a: "; m().display(tout, a.numerator()); tout << "/"; m().display(tout, a.denominator()); tout << "\n";
tout << "lower: "; bqm.display(tout, lower); tout << ", upper: "; bqm.display(tout, upper); tout << "\n";);
SASSERT(bqm.lt(lower, a));
@ -2848,7 +2849,7 @@ namespace upolynomial {
bqm.refine_upper(a, lower, upper);
}
SASSERT(bqm.lt(upper, b));
while (true) {
while (true) {
int sign_upper = eval_sign_at(sz, p, upper);
if (sign_upper == 0) {
// found root
@ -2872,7 +2873,7 @@ namespace upolynomial {
bqm.refine_upper(a, lower, upper);
}
}
if (!found_d) {
if (bqm.to_mpbq(b, lower)) {
// found d
@ -2893,7 +2894,7 @@ namespace upolynomial {
SASSERT(bqm.lt(c, lower));
SASSERT(bqm.lt(lower, upper));
SASSERT(bqm.lt(lower, b));
while (true) {
while (true) {
int sign_lower = eval_sign_at(sz, p, lower);
if (sign_lower == 0) {
// found root
@ -2946,8 +2947,8 @@ namespace upolynomial {
bool manager::normalize_interval(unsigned sz, numeral const * p, mpbq_manager & m, mpbq & a, mpbq & b) {
return normalize_interval_core(sz, p, INT_MIN, m, a, b);
}
}
unsigned manager::get_root_id(unsigned sz, numeral const * p, mpbq const & l) {
scoped_upolynomial_sequence seq(*this);
sturm_seq(sz, p, seq);
@ -2963,7 +2964,7 @@ namespace upolynomial {
m_manager.neg(new_c);
r.set_constant(new_c);
}
void manager::factor_2_sqf_pp(numeral_vector & p, factors & r, unsigned k) {
SASSERT(p.size() == 3); // p has degree 2
TRACE("factor", tout << "factor square free (degree == 2):\n"; display(tout, p); tout << "\n";);
@ -2980,7 +2981,7 @@ namespace upolynomial {
m().mul(a, c, ac);
m().addmul(b2, numeral(-4), ac, disc);
// discriminant must be different from 0, since p is square free
SASSERT(!m().is_zero(disc));
SASSERT(!m().is_zero(disc));
scoped_numeral disc_sqrt(m());
if (!m().is_perfect_square(disc, disc_sqrt)) {
// p is irreducible
@ -3053,13 +3054,13 @@ namespace upolynomial {
scoped_numeral_vector pp(m());
get_primitive_and_content(sz, p, pp, content);
r.set_constant(content);
//
//
scoped_numeral_vector & C = pp;
// Let C be a primitive polynomial of the form: P_1^1 * P_2^2 * ... * P_k^k, where each P_i is square free
scoped_numeral_vector C_prime(m());
derivative(C, C_prime);
scoped_numeral_vector A(m()), B(m()), D(m());
gcd(C, C_prime, B);
gcd(C, C_prime, B);
bool result = true;
if (is_const(B)) {
@ -3071,7 +3072,7 @@ namespace upolynomial {
}
else {
// B is of the form P_2 * P_3^2 * ... * P_k^{k-1}
VERIFY(exact_div(C, B, A));
VERIFY(exact_div(C, B, A));
TRACE("factor_bug", tout << "C: "; display(tout, C); tout << "\nB: "; display(tout, B); tout << "\nA: "; display(tout, A); tout << "\n";);
// A is of the form P_1 * P_2 * ... * P_k
unsigned j = 1;
@ -3084,7 +3085,7 @@ namespace upolynomial {
// D is of the form P_{j+1} * P_{j+2} * ... * P_k
VERIFY(exact_div(A, D, C));
// C is of the form P_j
if (!is_const(C)) {
if (!is_const(C)) {
flip_factor_sign_if_lm_neg(C, r, j);
if (!factor_sqf_pp(C, r, j, params))
result = false;
@ -3111,7 +3112,7 @@ namespace upolynomial {
bool manager::factor(unsigned sz, numeral const * p, factors & r, factor_params const & params) {
bool result = factor_core(sz, p, r, params);
#ifndef _EXTERNAL_RELEASE
#ifndef _EXTERNAL_RELEASE
IF_VERBOSE(FACTOR_VERBOSE_LVL, verbose_stream() << "(polynomial-factorization :distinct-factors " << r.distinct_factors() << ")" << std::endl;);
#endif
return result;

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

File diff suppressed because it is too large Load diff

126
src/math/subpaving/subpaving_t_def.h
View file

@ -31,8 +31,8 @@ namespace subpaving {
template<typename C>
class breadth_first_node_selector : public context_t<C>::node_selector {
typedef typename context_t<C>::node node;
public:
breadth_first_node_selector(context_t<C> * ctx):
public:
breadth_first_node_selector(context_t<C> * ctx):
context_t<C>::node_selector(ctx) {
}
@ -47,8 +47,8 @@ public:
template<typename C>
class depth_first_node_selector : public context_t<C>::node_selector {
typedef typename context_t<C>::node node;
public:
depth_first_node_selector(context_t<C> * ctx):
public:
depth_first_node_selector(context_t<C> * ctx):
context_t<C>::node_selector(ctx) {
}
@ -73,12 +73,12 @@ class round_robing_var_selector : public context_t<C>::var_selector {
x = 0;
}
public:
public:
round_robing_var_selector(context_t<C> * ctx, bool only_non_def = true):
context_t<C>::var_selector(ctx),
m_only_non_def(only_non_def) {
}
// Return the next variable to branch.
virtual var operator()(typename context_t<C>::node * n) {
typename context_t<C>::numeral_manager & nm = this->ctx()->nm();
@ -86,7 +86,7 @@ public:
var x = this->ctx()->splitting_var(n);
if (x == null_var)
x = 0;
else
else
next(x);
var start = x;
do {
@ -112,13 +112,13 @@ public:
- All variables x s.t. lower(x) == upper(x) are ignored.
- If node contains an unbounded variable (-oo, oo), then return the smallest unbounded variable.
- Otherwise, select the smallest variable x with the largest width, where
width is defined as:
If x has lower and upper bound, then width = upper(x) - lower(x).
If x has only lower, width = penalty/max(|lower(x)|, 1)
If x has only upper, width = penalty/max(|upper(x)|, 1)
penaly is a parameter of this class.
If x has only upper, width = penalty/max(|upper(x)|, 1)
penaly is a parameter of this class.
This strategy guarantees fairness.
*/
@ -127,13 +127,13 @@ class largest_interval_var_selector : public context_t<C>::var_selector {
unsigned m_penalty;
bool m_only_non_def;
public:
public:
largest_interval_var_selector(context_t<C> * ctx, unsigned unbounded_penalty = 10, bool only_non_def = true):
context_t<C>::var_selector(ctx),
m_penalty(unbounded_penalty),
m_only_non_def(only_non_def) {
}
// Return the next variable to branch.
virtual var operator()(typename context_t<C>::node * n) {
var y = null_var;
@ -149,7 +149,7 @@ public:
typename context_t<C>::bound * u = n->upper(x);
// variables without lower and upper bounds are selected immediately.
if (l == 0 && u == 0)
return x;
return x;
if (l != 0 && u != 0) {
// ignore variables s.t. lower(x) == upper(x)
if (nm.eq(l->value(), u->value()))
@ -165,7 +165,7 @@ public:
nm.set(width, l->value());
else
nm.set(width, u->value());
C::round_to_plus_inf(nm);
C::round_to_plus_inf(nm);
if (nm.is_neg(width))
nm.neg(width);
if (nm.lt(width, one))
@ -215,7 +215,7 @@ public:
nm.set(delta, static_cast<int>(m_delta));
nm.set(mid, upper->value());
C::round_to_minus_inf(nm);
nm.sub(mid, delta, mid);
nm.sub(mid, delta, mid);
// mid == upper - delta
}
else if (upper == 0) {
@ -224,7 +224,7 @@ public:
nm.set(delta, static_cast<int>(m_delta));
nm.set(mid, lower->value());
C::round_to_plus_inf(nm);
nm.add(mid, delta, mid);
nm.add(mid, delta, mid);
// mid == lower + delta
}
else {
@ -237,10 +237,10 @@ public:
throw subpaving::exception();
// mid == (lower + upper)/2
}
this->mk_decided_bound(x, mid, false, m_left_open, left);
this->mk_decided_bound(x, mid, true, !m_left_open, right);
TRACE("subpaving_int_split",
TRACE("subpaving_int_split",
tout << "LEFT:\n"; this->ctx()->display_bounds(tout, left);
tout << "\nRIGHT:\n"; this->ctx()->display_bounds(tout, right););
}
@ -342,7 +342,7 @@ void context_t<C>::node::push(bound * b) {
}
}
/**
/**
\brief Return the most recent variable that was used for splitting on node n.
*/
template<typename C>
@ -395,11 +395,11 @@ void context_t<C>::polynomial::display(std::ostream & out, numeral_manager & nm,
out << nm.to_rational_string(m_c);
first = false;
}
for (unsigned i = 0; i < m_size; i++) {
if (first)
first = false;
else
else
out << " + ";
if (!nm.is_one(a(i))) {
out << nm.to_rational_string(a(i));
@ -460,7 +460,7 @@ context_t<C>::~context_t() {
template<typename C>
void context_t<C>::checkpoint() {
if (!m_limit.inc())
throw default_exception("canceled");
throw default_exception(Z3_CANCELED_MSG);
if (memory::get_allocation_size() > m_max_memory)
throw default_exception(Z3_MAX_MEMORY_MSG);
cooperate("subpaving");
@ -497,7 +497,7 @@ void context_t<C>::updt_params(params_ref const & p) {
m_max_memory = megabytes_to_bytes(p.get_uint("max_memory", UINT_MAX));
unsigned prec = p.get_uint("nth_root_precision", 8192);
if (prec == 0)
if (prec == 0)
prec = 1;
nm().set(m_nth_root_prec, static_cast<int>(prec));
nm().inv(m_nth_root_prec);
@ -513,7 +513,7 @@ void context_t<C>::collect_param_descrs(param_descrs & d) {
}
template<typename C>
void context_t<C>::display_params(std::ostream & out) const {
void context_t<C>::display_params(std::ostream & out) const {
out << "max_nodes " << m_max_nodes << "\n";
out << "max_depth " << m_max_depth << "\n";
out << "epsilon " << nm().to_rational_string(m_epsilon) << "\n";
@ -626,7 +626,7 @@ void context_t<C>::display_definition(std::ostream & out, definition const * d,
template<typename C>
void context_t<C>::display(std::ostream & out, constraint * c, bool use_star) const {
if (c->get_kind() == constraint::CLAUSE)
if (c->get_kind() == constraint::CLAUSE)
static_cast<clause*>(c)->display(out, nm(), *m_display_proc);
else
display_definition(out, static_cast<definition*>(c), use_star);
@ -639,7 +639,7 @@ void context_t<C>::display_bounds(std::ostream & out, node * n) const {
bound * l = n->lower(x);
bound * u = n->upper(x);
if (l != 0) {
display(out, l);
display(out, l);
out << " ";
}
if (u != 0) {
@ -782,18 +782,18 @@ typename context_t<C>::ineq * context_t<C>::mk_ineq(var x, numeral const & k, bo
}
template<typename C>
void context_t<C>::inc_ref(ineq * a) {
void context_t<C>::inc_ref(ineq * a) {
TRACE("subpaving_ref_count", tout << "inc-ref: " << a << " " << a->m_ref_count << "\n";);
if (a)
a->m_ref_count++;
if (a)
a->m_ref_count++;
}
template<typename C>
void context_t<C>::dec_ref(ineq * a) {
if (a) {
TRACE("subpaving_ref_count",
tout << "dec-ref: " << a << " " << a->m_ref_count << "\n";
a->display(tout, nm());
TRACE("subpaving_ref_count",
tout << "dec-ref: " << a << " " << a->m_ref_count << "\n";
a->display(tout, nm());
tout << "\n";);
SASSERT(a->m_ref_count > 0);
a->m_ref_count--;
@ -804,7 +804,7 @@ void context_t<C>::dec_ref(ineq * a) {
}
}
}
template<typename C>
void context_t<C>::add_clause_core(unsigned sz, ineq * const * atoms, bool lemma, bool watch) {
SASSERT(lemma || watch);
@ -846,7 +846,7 @@ void context_t<C>::del_clause(clause * c) {
SASSERT(c->m_num_jst == 0); // We cannot delete a clause that is being used to justify some bound
bool watch = c->watched();
var prev_x = null_var;
unsigned sz = c->size();
unsigned sz = c->size();
for (unsigned i = 0; i < sz; i++) {
var x = c->m_atoms[i]->x();
if (watch) {
@ -883,10 +883,10 @@ typename context_t<C>::node * context_t<C>::mk_node(node * parent) {
else
r = new (mem) node(parent, m_node_id_gen.mk());
m_var_selector->new_node_eh(r);
// Add node in the leaf dlist
push_front(r);
m_num_nodes++;
return r;
}
@ -903,9 +903,9 @@ void context_t<C>::del_node(node * n) {
// recycle id
m_node_id_gen.recycle(n->id());
// disconnect n from list of leaves.
// disconnect n from list of leaves.
remove_from_leaf_dlist(n);
// disconnect n from parent
node * p = n->parent();
bound * b = n->trail_stack();
@ -1118,7 +1118,7 @@ void context_t<C>::del_definitions() {
if (d == 0)
continue;
switch (d->get_kind()) {
case constraint::MONOMIAL:
case constraint::MONOMIAL:
del_monomial(static_cast<monomial*>(d));
break;
case constraint::POLYNOMIAL:
@ -1130,7 +1130,7 @@ void context_t<C>::del_definitions() {
}
}
}
template<typename C>
void context_t<C>::display_constraints(std::ostream & out, bool use_star) const {
// display definitions
@ -1160,9 +1160,9 @@ void context_t<C>::display_constraints(std::ostream & out, bool use_star) const
// -----------------------------------
template<typename C>
void context_t<C>::set_conflict(var x, node * n) {
void context_t<C>::set_conflict(var x, node * n) {
m_num_conflicts++;
n->set_conflict(x);
n->set_conflict(x);
remove_from_leaf_dlist(n);
}
@ -1222,12 +1222,12 @@ bool context_t<C>::relevant_new_bound(var x, numeral const & k, bool lower, bool
if (m_zero_epsilon && curr_lower != 0 && (nm().lt(k, curr_lower->value()) || ((curr_lower->is_open() || !open) && nm().eq(k, curr_lower->value())))) {
// new lower bound does not improve existing bound
TRACE("subpaving_relevant_bound", tout << "irrelevant because does not improve existing bound.\n";);
return false;
return false;
}
if (curr_upper == 0 && nm().lt(m_max_bound, k)) {
// new lower bound exceeds the :max-bound threshold
TRACE("subpaving_relevant_bound", tout << "irrelevant because exceeds :max-bound threshold.\n";);
return false;
return false;
}
if (!m_zero_epsilon && curr_lower != 0) {
// check if:
@ -1250,10 +1250,10 @@ bool context_t<C>::relevant_new_bound(var x, numeral const & k, bool lower, bool
nm().set(delta, min);
nm().mul(delta, m_epsilon, delta);
nm().add(curr_lower->value(), delta, delta);
TRACE("subpaving_relevant_bound_bug",
tout << "k: "; nm().display(tout, k);
TRACE("subpaving_relevant_bound_bug",
tout << "k: "; nm().display(tout, k);
tout << ", delta: "; nm().display(tout, delta); tout << "\n";
tout << "curr_lower: "; nm().display(tout, curr_lower->value());
tout << "curr_lower: "; nm().display(tout, curr_lower->value());
tout << ", min: "; nm().display(tout, min); tout << "\n";);
if (nm().le(k, delta)) {
TRACE("subpaving_relevant_bound", tout << "irrelevant because does not improve existing bound to at least ";
@ -1272,7 +1272,7 @@ bool context_t<C>::relevant_new_bound(var x, numeral const & k, bool lower, bool
if (m_zero_epsilon && curr_upper != 0 && (nm().lt(curr_upper->value(), k) || ((curr_upper->is_open() || !open) && nm().eq(k, curr_upper->value())))) {
// new upper bound does not improve existing bound
TRACE("subpaving_relevant_bound", tout << "irrelevant because does not improve existing bound.\n";);
return false;
return false;
}
if (curr_lower == 0 && nm().lt(k, m_minus_max_bound)) {
// new upper bound exceeds the -:max-bound threshold
@ -1343,7 +1343,7 @@ bool context_t<C>::conflicting_bounds(var x, node * n) const {
/**
\brief Return the truth value of the inequality t in node n.
The result may be l_true (True), l_false (False), or l_undef(Unknown).
*/
template<typename C>
@ -1399,7 +1399,7 @@ void context_t<C>::propagate_clause(clause * c, node * n) {
ineq * a = (*c)[j];
TRACE("propagate_clause", tout << "propagating inequality: "; display(tout, a); tout << "\n";);
propagate_bound(a->x(), a->value(), a->is_lower(), a->is_open(), n, justification(c));
// A clause can propagate only once.
// A clause can propagate only once.
// So, we can safely set its timestamp again to avoid another useless visit.
c->set_visited(m_timestamp);
}
@ -1411,7 +1411,7 @@ void context_t<C>::propagate_polynomial(var x, node * n, var y) {
polynomial * p = get_polynomial(x);
unsigned sz = p->size();
interval & r = m_i_tmp1; r.set_mutable();
interval & v = m_i_tmp2;
interval & v = m_i_tmp2;
interval & av = m_i_tmp3; av.set_mutable();
if (x == y) {
for (unsigned i = 0; i < sz; i++) {
@ -1482,7 +1482,7 @@ void context_t<C>::propagate_polynomial(var x, node * n) {
}
}
TRACE("propagate_polynomial", tout << "unbounded_var: "; display(tout, unbounded_var); tout << "\n";);
if (unbounded_var != null_var) {
propagate_polynomial(x, n, unbounded_var);
}
@ -1528,10 +1528,10 @@ void context_t<C>::propagate_monomial(var x, node * n) {
// x must be zero
numeral & zero = m_tmp1;
nm().set(zero, 0);
propagate_bound(x, zero, true, false, n, justification(x));
propagate_bound(x, zero, true, false, n, justification(x));
if (inconsistent(n))
return;
propagate_bound(x, zero, false, false, n, justification(x));
propagate_bound(x, zero, false, false, n, justification(x));
}
// no need to downward propagation
return;
@ -1572,11 +1572,11 @@ void context_t<C>::propagate_monomial_upward(var x, node * n) {
monomial * m = get_monomial(x);
unsigned sz = m->size();
interval & r = m_i_tmp1; r.set_mutable();
interval & y = m_i_tmp2;
interval & y = m_i_tmp2;
interval & yk = m_i_tmp3; yk.set_mutable();
for (unsigned i = 0; i < sz; i++) {
y.set_constant(n, m->x(i));
im().power(y, m->degree(i), yk);
im().power(y, m->degree(i), yk);
if (i == 0)
im().set(r, yk);
else
@ -1606,37 +1606,37 @@ void context_t<C>::propagate_monomial_downward(var x, node * n, unsigned j) {
monomial * m = get_monomial(x);
SASSERT(j < m->size());
unsigned sz = m->size();
interval & r = m_i_tmp3;
if (sz > 1) {
interval & d = m_i_tmp1; d.set_mutable();
interval & y = m_i_tmp2;
interval & y = m_i_tmp2;
interval & yk = m_i_tmp3; yk.set_mutable();
bool first = true;
for (unsigned i = 0; i < sz; i++) {
if (i == j)
continue;
y.set_constant(n, m->x(i));
im().power(y, m->degree(i), yk);
im().power(y, m->degree(i), yk);
if (first)
im().set(d, yk);
else
im().mul(d, yk, r);
}
interval & aux = m_i_tmp2;
interval & aux = m_i_tmp2;
aux.set_constant(n, x);
im().div(aux, d, r);
}
else {
SASSERT(sz == 1);
SASSERT(j == 0);
interval & aux = m_i_tmp2;
interval & aux = m_i_tmp2;
aux.set_constant(n, x);
im().set(r, aux);
}
unsigned d = m->degree(j);
if (d > 1) {
if (d % 2 == 0 && im().lower_is_neg(r))
if (d % 2 == 0 && im().lower_is_neg(r))
return; // If d is even, we can't take the nth-root when lower(r) is negative.
im().xn_eq_y(r, d, m_nth_root_prec, r);
}
@ -1814,7 +1814,7 @@ void context_t<C>::operator()() {
if (m_num_nodes > m_max_nodes)
break;
node * n = (*m_node_selector)(m_leaf_head, m_leaf_tail);
if (n == 0)
if (n == 0)
break;
TRACE("subpaving_main", tout << "selected node: #" << n->id() << ", depth: " << n->depth() << "\n";);
remove_from_leaf_dlist(n);

47
src/math/subpaving/tactic/expr2subpaving.cpp
View file

@ -8,8 +8,8 @@ Module Name:
Abstract:
Translator from Z3 expressions into generic subpaving data-structure.
Author:
Leonardo (leonardo) 2012-08-08
@ -24,6 +24,7 @@ Notes:
#include"cooperate.h"
#include"arith_decl_plugin.h"
#include"scoped_numeral_buffer.h"
#include"common_msgs.h"
struct expr2subpaving::imp {
struct frame {
@ -43,9 +44,9 @@ struct expr2subpaving::imp {
expr_ref_vector m_var2expr;
typedef svector<subpaving::var> var_vector;
obj_map<expr, unsigned> m_cache;
var_vector m_cached_vars;
var_vector m_cached_vars;
scoped_mpz_vector m_cached_numerators;
scoped_mpz_vector m_cached_denominators;
@ -69,17 +70,17 @@ struct expr2subpaving::imp {
m_expr2var = e2v;
m_expr2var_owner = false;
}
}
~imp() {
reset_cache();
if (m_expr2var_owner)
dealloc(m_expr2var);
}
ast_manager & m() { return m_manager; }
subpaving::context & s() { return m_subpaving; }
unsynch_mpq_manager & qm() const { return m_qm; }
@ -94,7 +95,7 @@ struct expr2subpaving::imp {
void checkpoint() {
if (m().canceled())
throw default_exception("canceled");
throw default_exception(Z3_CANCELED_MSG);
cooperate("expr2subpaving");
}
@ -111,7 +112,7 @@ struct expr2subpaving::imp {
}
return x;
}
void found_non_simplified() {
throw default_exception("you must apply simplifier before internalizing expressions into the subpaving module.");
}
@ -128,7 +129,7 @@ struct expr2subpaving::imp {
SASSERT(!m_cache.contains(t));
SASSERT(m_cached_numerators.size() == m_cached_vars.size());
SASSERT(m_cached_denominators.size() == m_cached_vars.size());
if (t->get_ref_count() <= 1)
if (t->get_ref_count() <= 1)
return;
unsigned idx = m_cached_vars.size();
m_cache.insert(t, idx);
@ -196,7 +197,7 @@ struct expr2subpaving::imp {
sbuffer<subpaving::power> pws;
for (unsigned i = 0; i < sz; i++) {
expr * arg = margs[i];
unsigned k;
unsigned k;
as_power(arg, arg, k);
subpaving::var x_arg = process(arg, depth+1, n_arg, d_arg);
qm().power(n_arg, k, n_arg);
@ -216,7 +217,7 @@ struct expr2subpaving::imp {
cache_result(t, x, n, d);
return x;
}
typedef _scoped_numeral_buffer<unsynch_mpz_manager> mpz_buffer;
typedef sbuffer<subpaving::var> var_buffer;
@ -291,13 +292,13 @@ struct expr2subpaving::imp {
switch (t->get_decl_kind()) {
case OP_NUM:
return process_num(t, depth, n, d);
case OP_ADD:
case OP_ADD:
return process_add(t, depth, n, d);
case OP_MUL:
case OP_MUL:
return process_mul(t, depth, n, d);
case OP_POWER:
return process_power(t, depth, n, d);
case OP_TO_REAL:
case OP_TO_REAL:
return process(t->get_arg(0), depth+1, n, d);
case OP_SUB:
case OP_UMINUS:
@ -310,7 +311,7 @@ struct expr2subpaving::imp {
case OP_REM:
case OP_IRRATIONAL_ALGEBRAIC_NUM:
throw default_exception("you must apply arithmetic purifier before internalizing expressions into the subpaving module.");
case OP_SIN:
case OP_SIN:
case OP_COS:
case OP_TAN:
case OP_ASIN:
@ -350,12 +351,12 @@ struct expr2subpaving::imp {
return process_arith_app(to_app(t), depth, n, d);
}
bool is_var(expr * t) const {
return m_expr2var->is_var(t);
bool is_var(expr * t) const {
return m_expr2var->is_var(t);
}
subpaving::var internalize_term(expr * t, mpz & n, mpz & d) {
return process(t, 0, n, d);
}
@ -376,11 +377,11 @@ ast_manager & expr2subpaving::m() const {
subpaving::context & expr2subpaving::s() const {
return m_imp->s();
}
bool expr2subpaving::is_var(expr * t) const {
return m_imp->is_var(t);
}
subpaving::var expr2subpaving::internalize_term(expr * t, mpz & n, mpz & d) {
return m_imp->internalize_term(t, n, d);

324
src/muz/bmc/dl_bmc_engine.cpp
View file

@ -60,11 +60,11 @@ namespace datalog {
expr_ref fml(m.mk_app(q, T), m);
b.assert_expr(fml);
res = b.m_solver.check();
if (res == l_true) {
res = get_model();
}
b.m_solver.pop(1);
b.m_solver.pop(1);
++m_bit_width;
}
return res;
@ -74,11 +74,11 @@ namespace datalog {
sort_ref mk_index_sort() {
return sort_ref(m_bv.mk_sort(m_bit_width), m);
}
var_ref mk_index_var() {
return var_ref(m.mk_var(0, mk_index_sort()), m);
}
void compile() {
sort_ref index_sort = mk_index_sort();
var_ref var = mk_index_var();
@ -92,17 +92,17 @@ namespace datalog {
// Assert: forall level . p(T) => body of rule i + equalities for head of rule i
func_decl_ref pr = mk_q_func_decl(p);
expr_ref pred = expr_ref(m.mk_app(pr, var.get()), m);
expr_ref_vector rules(m), sub(m), conjs(m);
expr_ref_vector rules(m), sub(m), conjs(m);
expr_ref trm(m), rule_body(m), rule_i(m);
for (unsigned i = 0; i < rls.size(); ++i) {
sub.reset();
for (unsigned i = 0; i < rls.size(); ++i) {
sub.reset();
conjs.reset();
rule& r = *rls[i];
rule& r = *rls[i];
rule_i = m.mk_app(mk_q_rule(p, i), var.get());
rules.push_back(rule_i);
mk_qrule_vars(r, i, sub);
// apply substitution to body.
var_subst vs(m, false);
for (unsigned k = 0; k < p->get_arity(); ++k) {
@ -137,21 +137,21 @@ namespace datalog {
b.assert_expr(trm);
}
}
void setup() {
b.m_fparams.m_relevancy_lvl = 2;
b.m_fparams.m_model = true;
b.m_fparams.m_model_compact = true;
b.m_fparams.m_mbqi = true;
}
void mk_qrule_vars(datalog::rule const& r, unsigned rule_id, expr_ref_vector& sub) {
ptr_vector<sort> sorts;
r.get_vars(m, sorts);
// populate substitution of bound variables.
sub.reset();
sub.resize(sorts.size());
for (unsigned k = 0; k < r.get_decl()->get_arity(); ++k) {
expr* arg = r.get_head()->get_arg(k);
if (is_var(arg)) {
@ -177,52 +177,52 @@ namespace datalog {
if (sorts[j] && !sub[j].get()) {
sub[j] = mk_q_var(r.get_decl(), sorts[j], rule_id, idx++);
}
}
}
}
expr_ref mk_q_var(func_decl* pred, sort* s, unsigned rule_id, unsigned idx) {
std::stringstream _name;
_name << pred->get_name() << "#" << rule_id << "_" << idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
var_ref var = mk_index_var();
return expr_ref(m.mk_app(m.mk_func_decl(nm, mk_index_sort(), s), var), m);
}
}
expr_ref mk_q_arg(func_decl* pred, unsigned idx, bool is_current) {
SASSERT(idx < pred->get_arity());
std::stringstream _name;
_name << pred->get_name() << "#" << idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
expr_ref var(mk_index_var(), m);
if (!is_current) {
var = m_bv.mk_bv_sub(var, mk_q_one());
}
return expr_ref(m.mk_app(m.mk_func_decl(nm, mk_index_sort(), pred->get_domain(idx)), var), m);
}
expr_ref mk_q_one() {
return mk_q_num(1);
}
expr_ref mk_q_num(unsigned i) {
return expr_ref(m_bv.mk_numeral(i, m_bit_width), m);
}
func_decl_ref mk_q_func_decl(func_decl* f) {
std::stringstream _name;
_name << f->get_name() << "#";
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return func_decl_ref(m.mk_func_decl(nm, mk_index_sort(), f->get_range()), m);
}
func_decl_ref mk_q_rule(func_decl* f, unsigned rule_id) {
std::stringstream _name;
_name << f->get_name() << "#" << rule_id;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return func_decl_ref(m.mk_func_decl(nm, mk_index_sort(), m.mk_bool_sort()), m);
}
expr_ref eval_q(model_ref& model, func_decl* f, unsigned i) {
func_decl_ref fn = mk_q_func_decl(f);
expr_ref t(m), result(m);
@ -230,7 +230,7 @@ namespace datalog {
model->eval(t, result);
return result;
}
expr_ref eval_q(model_ref& model, expr* t, unsigned i) {
expr_ref tmp(m), result(m), num(m);
var_subst vs(m, false);
@ -239,8 +239,8 @@ namespace datalog {
vs(t, 1, nums, tmp);
model->eval(tmp, result);
return result;
}
}
lbool get_model() {
rule_manager& rm = b.m_ctx.get_rule_manager();
func_decl_ref q = mk_q_func_decl(b.m_query_pred);
@ -250,7 +250,7 @@ namespace datalog {
rule_unifier unifier(b.m_ctx);
rational num;
unsigned level, bv_size;
b.m_solver.get_model(md);
func_decl* pred = b.m_query_pred;
dl_decl_util util(m);
@ -261,7 +261,7 @@ namespace datalog {
level = num.get_unsigned();
SASSERT(m.is_true(eval_q(md, b.m_query_pred, level)));
TRACE("bmc", model_smt2_pp(tout, m, *md, 0););
rule_ref r0(rm), r1(rm), r2(rm);
while (true) {
TRACE("bmc", tout << "Predicate: " << pred->get_name() << "\n";);
@ -269,8 +269,8 @@ namespace datalog {
rule_vector const& rls = b.m_rules.get_predicate_rules(pred);
rule* r = 0;
unsigned i = 0;
for (; i < rls.size(); ++i) {
rule_i = m.mk_app(mk_q_rule(pred, i), mk_q_num(level).get());
for (; i < rls.size(); ++i) {
rule_i = m.mk_app(mk_q_rule(pred, i), mk_q_num(level).get());
TRACE("bmc", rls[i]->display(b.m_ctx, tout << "Checking rule " << mk_pp(rule_i, m) << " "););
if (m.is_true(eval_q(md, rule_i, level))) {
r = rls[i];
@ -280,13 +280,13 @@ namespace datalog {
SASSERT(r);
mk_qrule_vars(*r, i, sub);
// we have rule, we have variable names of rule.
// extract values for the variables in the rule.
for (unsigned j = 0; j < sub.size(); ++j) {
expr_ref vl = eval_q(md, sub[j].get(), i);
if (vl) {
// vl can be 0 if the interpretation does not assign a value to it.
sub[j] = vl;
sub[j] = vl;
}
else {
sub[j] = m.mk_var(j, m.get_sort(sub[j].get()));
@ -295,7 +295,7 @@ namespace datalog {
svector<std::pair<unsigned, unsigned> > positions;
vector<expr_ref_vector> substs;
expr_ref fml(m), concl(m);
rm.to_formula(*r, fml);
r2 = r;
rm.substitute(r2, sub.size(), sub.c_ptr());
@ -308,15 +308,15 @@ namespace datalog {
unifier.apply(*r0.get(), 0, *r2.get(), r1);
rm.to_formula(*r1.get(), concl);
scoped_proof _sp(m);
p = r->get_proof();
if (!p) {
p = m.mk_asserted(fml);
}
proof* premises[2] = { pr, p };
positions.push_back(std::make_pair(0, 1));
substs.push_back(sub1);
substs.push_back(sub);
pr = m.mk_hyper_resolve(2, premises, concl, positions, substs);
@ -339,7 +339,7 @@ namespace datalog {
}
r0 = r2;
}
if (level == 0) {
SASSERT(r->get_uninterpreted_tail_size() == 0);
break;
@ -370,7 +370,7 @@ namespace datalog {
lbool check() {
setup();
for (unsigned i = 0; ; ++i) {
IF_VERBOSE(1, verbose_stream() << "level: " << i << "\n";);
IF_VERBOSE(1, verbose_stream() << "level: " << i << "\n";);
b.checkpoint();
expr_ref_vector fmls(m);
compile(b.m_rules, fmls, i);
@ -392,12 +392,12 @@ namespace datalog {
for (unsigned i = 0; i < level_p->get_arity(); ++i) {
std::stringstream _name;
_name << query_pred->get_name() << "#" << level << "_" << i;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
vars.push_back(m.mk_const(nm, level_p->get_domain(i)));
}
return expr_ref(m.mk_app(level_p, vars.size(), vars.c_ptr()), m);
}
void compile(rule_set const& rules, expr_ref_vector& result, unsigned level) {
bool_rewriter br(m);
rule_set::decl2rules::iterator it = rules.begin_grouped_rules();
@ -405,16 +405,16 @@ namespace datalog {
for (; it != end; ++it) {
func_decl* p = it->m_key;
rule_vector const& rls = *it->m_value;
// Assert: p_level(vars) => r1_level(vars) \/ r2_level(vars) \/ r3_level(vars) \/ ...
// Assert: r_i_level(vars) => exists aux_vars . body of rule i for level
// Assert: r_i_level(vars) => exists aux_vars . body of rule i for level
func_decl_ref level_pred = mk_level_predicate(p, level);
expr_ref_vector rules(m);
expr_ref body(m), head(m);
for (unsigned i = 0; i < rls.size(); ++i) {
rule& r = *rls[i];
func_decl_ref rule_i = mk_level_rule(p, i, level);
for (unsigned i = 0; i < rls.size(); ++i) {
rule& r = *rls[i];
func_decl_ref rule_i = mk_level_rule(p, i, level);
rules.push_back(apply_vars(rule_i));
ptr_vector<sort> rule_vars;
@ -440,10 +440,10 @@ namespace datalog {
head = m.mk_app(rule_i, args.size(), args.c_ptr());
for (unsigned i = 0; i < r.get_tail_size(); ++i) {
conjs.push_back(r.get_tail(i));
}
}
br.mk_and(conjs.size(), conjs.c_ptr(), body);
replace_by_level_predicates(level, body);
replace_by_level_predicates(level, body);
body = skolemize_vars(r, args, rule_vars, body);
body = m.mk_implies(head, body);
body = bind_vars(body, head);
@ -496,7 +496,7 @@ namespace datalog {
// find the rule that was triggered by evaluating the arguments to prop.
rule_vector const& rls = b.m_rules.get_predicate_rules(pred);
rule* r = 0;
for (unsigned i = 0; i < rls.size(); ++i) {
for (unsigned i = 0; i < rls.size(); ++i) {
func_decl_ref rule_i = mk_level_rule(pred, i, level);
TRACE("bmc", rls[i]->display(b.m_ctx, tout << "Checking rule " << mk_pp(rule_i, m) << " "););
prop_r = m.mk_app(rule_i, prop->get_num_args(), prop->get_args());
@ -537,7 +537,7 @@ namespace datalog {
prs.push_back(get_proof(md, head_j, to_app(prop_body), level-1));
positions.push_back(std::make_pair(j+1,0));
substs.push_back(expr_ref_vector(m));
}
}
pr = m.mk_hyper_resolve(sz+1, prs.c_ptr(), prop, positions, substs);
return pr;
}
@ -557,16 +557,16 @@ namespace datalog {
func_decl_ref mk_level_predicate(func_decl* p, unsigned level) {
std::stringstream _name;
_name << p->get_name() << "#" << level;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return func_decl_ref(m.mk_func_decl(nm, p->get_arity(), p->get_domain(), m.mk_bool_sort()), m);
}
func_decl_ref mk_level_rule(func_decl* p, unsigned rule_idx, unsigned level) {
std::stringstream _name;
_name << "rule:" << p->get_name() << "#" << level << "_" << rule_idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return func_decl_ref(m.mk_func_decl(nm, p->get_arity(), p->get_domain(), m.mk_bool_sort()), m);
}
}
expr_ref apply_vars(func_decl* p) {
expr_ref_vector vars(m);
@ -617,7 +617,7 @@ namespace datalog {
func_decl_ref mk_body_func(rule& r, ptr_vector<sort> const& args, unsigned index, sort* s) {
std::stringstream _name;
_name << r.get_decl()->get_name() << "@" << index;
symbol name(_name.str().c_str());
symbol name(_name.str().c_str());
func_decl* f = m.mk_func_decl(name, args.size(), args.c_ptr(), s);
return func_decl_ref(f, m);
}
@ -689,16 +689,16 @@ namespace datalog {
struct level_replacer_cfg : public default_rewriter_cfg {
level_replacer m_r;
level_replacer_cfg(nonlinear& nl, unsigned level):
level_replacer_cfg(nonlinear& nl, unsigned level):
m_r(nl, level) {}
bool rewrite_patterns() const { return false; }
br_status reduce_app(func_decl * f, unsigned num, expr * const * args, expr_ref & result, proof_ref & result_pr) {
return m_r.mk_app_core(f, num, args, result);
}
bool reduce_quantifier(quantifier * old_q,
expr * new_body,
expr * const * new_patterns,
bool reduce_quantifier(quantifier * old_q,
expr * new_body,
expr * const * new_patterns,
expr * const * new_no_patterns,
expr_ref & result,
proof_ref & result_pr) {
@ -725,9 +725,9 @@ namespace datalog {
fml = tmp;
}
};
// --------------------------------------------------------------------------
// Basic non-linear BMC based on compiling into data-types (it is inefficient)
@ -747,7 +747,7 @@ namespace datalog {
setup();
declare_datatypes();
compile();
return check_query();
return check_query();
}
private:
@ -765,7 +765,7 @@ namespace datalog {
func_decl_ref mk_predicate(func_decl* pred) {
std::stringstream _name;
_name << pred->get_name() << "#";
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
sort* pred_trace_sort = m_pred2sort.find(pred);
return func_decl_ref(m.mk_func_decl(nm, pred_trace_sort, m_path_sort, m.mk_bool_sort()), m);
}
@ -773,30 +773,30 @@ namespace datalog {
func_decl_ref mk_rule(func_decl* p, unsigned rule_idx) {
std::stringstream _name;
_name << "rule:" << p->get_name() << "#" << rule_idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
sort* pred_trace_sort = m_pred2sort.find(p);
return func_decl_ref(m.mk_func_decl(nm, pred_trace_sort, m_path_sort, m.mk_bool_sort()), m);
return func_decl_ref(m.mk_func_decl(nm, pred_trace_sort, m_path_sort, m.mk_bool_sort()), m);
}
expr_ref mk_var(func_decl* pred, sort*s, unsigned idx, expr* path_arg, expr* trace_arg) {
std::stringstream _name;
_name << pred->get_name() << "#V_" << idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
func_decl_ref fn(m);
fn = m.mk_func_decl(nm, m_pred2sort.find(pred), m_path_sort, s);
return expr_ref(m.mk_app(fn, trace_arg, path_arg), m);
}
expr_ref mk_arg(func_decl* pred, unsigned idx, expr* path_arg, expr* trace_arg) {
SASSERT(idx < pred->get_arity());
std::stringstream _name;
_name << pred->get_name() << "#X_" << idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
func_decl_ref fn(m);
fn = m.mk_func_decl(nm, m_pred2sort.find(pred), m_path_sort, pred->get_domain(idx));
return expr_ref(m.mk_app(fn, trace_arg, path_arg), m);
}
void mk_subst(datalog::rule& r, expr* path, app* trace, expr_ref_vector& sub) {
datatype_util dtu(m);
ptr_vector<sort> sorts;
@ -840,52 +840,52 @@ namespace datalog {
}
}
}
/**
\brief compile Horn rule into co-Horn implication.
forall args . R(path_var, rule_i(trace_vars)) => Body[X(path_var, rule_i(trace_vars)), Y(S_j(path_var), trace_vars_j)]
*/
*/
void compile() {
datatype_util dtu(m);
rule_set::decl2rules::iterator it = b.m_rules.begin_grouped_rules();
rule_set::decl2rules::iterator end = b.m_rules.end_grouped_rules();
for (; it != end; ++it) {
func_decl* p = it->m_key;
rule_vector const& rls = *it->m_value;
// Assert: p_level => r1_level \/ r2_level \/ r3_level \/ ...
// where: r_i_level = body of rule i for level + equalities for head of rule i
expr_ref rule_body(m), tmp(m), pred(m), trace_arg(m), fml(m);
var_ref path_var(m), trace_var(m);
expr_ref_vector rules(m), sub(m), conjs(m), vars(m), patterns(m);
sort* pred_sort = m_pred2sort.find(p);
path_var = m.mk_var(0, m_path_sort);
trace_var = m.mk_var(1, pred_sort);
trace_var = m.mk_var(1, pred_sort);
// sort* sorts[2] = { pred_sort, m_path_sort };
ptr_vector<func_decl> const& cnstrs = *dtu.get_datatype_constructors(pred_sort);
ptr_vector<func_decl> const& succs = *dtu.get_datatype_constructors(m_path_sort);
SASSERT(cnstrs.size() == rls.size());
pred = m.mk_app(mk_predicate(p), trace_var.get(), path_var.get());
for (unsigned i = 0; i < rls.size(); ++i) {
sub.reset();
for (unsigned i = 0; i < rls.size(); ++i) {
sub.reset();
conjs.reset();
vars.reset();
rule& r = *rls[i];
rule& r = *rls[i];
func_decl_ref rule_pred_i = mk_rule(p, i);
// Create cnstr_rule_i(Vars)
func_decl* cnstr = cnstrs[i];
rules.push_back(m.mk_app(rule_pred_i, trace_var.get(), path_var.get()));
unsigned arity = cnstr->get_arity();
for (unsigned j = 0; j < arity; ++j) {
vars.push_back(m.mk_var(arity-j,cnstr->get_domain(j)));
vars.push_back(m.mk_var(arity-j,cnstr->get_domain(j)));
}
trace_arg = m.mk_app(cnstr, vars.size(), vars.c_ptr());
mk_subst(r, path_var, to_app(trace_arg), sub);
// apply substitution to body.
var_subst vs(m, false);
for (unsigned k = 0; k < p->get_arity(); ++k) {
@ -926,28 +926,28 @@ namespace datalog {
names.push_back(symbol("path"));
SASSERT(names.size() == q_sorts.size());
SASSERT(vars.size() == names.size());
symbol qid = r.name(), skid;
symbol qid = r.name(), skid;
tmp = m.mk_app(mk_predicate(p), trace_arg.get(), path_var.get());
patterns.reset();
patterns.push_back(m.mk_pattern(to_app(tmp)));
fml = m.mk_implies(tmp, rule_body);
fml = m.mk_forall(vars.size(), q_sorts.c_ptr(), names.c_ptr(), fml, 1, qid, skid, 1, patterns.c_ptr());
b.assert_expr(fml);
b.assert_expr(fml);
}
}
}
}
void declare_datatypes() {
rule_set::decl2rules::iterator it = b.m_rules.begin_grouped_rules();
rule_set::decl2rules::iterator end = b.m_rules.end_grouped_rules();
datatype_util dtu(m);
ptr_vector<datatype_decl> dts;
obj_map<func_decl, unsigned> pred_idx;
for (unsigned i = 0; it != end; ++it, ++i) {
pred_idx.insert(it->m_key, i);
}
it = b.m_rules.begin_grouped_rules();
for (; it != end; ++it) {
rule_vector const& rls = *it->m_value;
@ -971,27 +971,27 @@ namespace datalog {
_name << "?";
symbol is_name(_name.str().c_str());
cnstrs.push_back(mk_constructor_decl(name, is_name, accs.size(), accs.c_ptr()));
}
}
dts.push_back(mk_datatype_decl(pred->get_name(), cnstrs.size(), cnstrs.c_ptr()));
}
sort_ref_vector new_sorts(m);
family_id dfid = m.mk_family_id("datatype");
datatype_decl_plugin* dtp = static_cast<datatype_decl_plugin*>(m.get_plugin(dfid));
VERIFY (dtp->mk_datatypes(dts.size(), dts.c_ptr(), new_sorts));
it = b.m_rules.begin_grouped_rules();
for (unsigned i = 0; it != end; ++it, ++i) {
for (unsigned i = 0; it != end; ++it, ++i) {
m_pred2sort.insert(it->m_key, new_sorts[i].get());
m_sort2pred.insert(new_sorts[i].get(), it->m_key);
m_pinned.push_back(new_sorts[i].get());
m_pinned.push_back(new_sorts[i].get());
}
if (new_sorts.size() > 0) {
TRACE("bmc", dtu.display_datatype(new_sorts[0].get(), tout););
}
del_datatype_decls(dts.size(), dts.c_ptr());
// declare path data-type.
{
new_sorts.reset();
@ -1004,9 +1004,9 @@ namespace datalog {
rule* r = *it;
unsigned sz = r->get_uninterpreted_tail_size();
max_arity = std::max(sz, max_arity);
}
}
cnstrs.push_back(mk_constructor_decl(symbol("Z#"), symbol("Z#?"), 0, 0));
for (unsigned i = 0; i + 1 < max_arity; ++i) {
std::stringstream _name;
_name << "succ#" << i;
@ -1019,13 +1019,13 @@ namespace datalog {
type_ref tr(0);
accs.push_back(mk_accessor_decl(name, tr));
cnstrs.push_back(mk_constructor_decl(name, is_name, accs.size(), accs.c_ptr()));
}
}
dts.push_back(mk_datatype_decl(symbol("Path"), cnstrs.size(), cnstrs.c_ptr()));
VERIFY (dtp->mk_datatypes(dts.size(), dts.c_ptr(), new_sorts));
m_path_sort = new_sorts[0].get();
}
}
proof_ref get_proof(model_ref& md, app* trace, app* path) {
datatype_util dtu(m);
rule_manager& rm = b.m_ctx.get_rule_manager();
@ -1044,7 +1044,7 @@ namespace datalog {
proof_ref pr(m);
expr_ref fml(m), head(m), tmp(m);
app_ref path1(m);
var_subst vs(m, false);
mk_subst(*rules[i], path, trace, sub);
rm.to_formula(*rules[i], fml);
@ -1061,30 +1061,30 @@ namespace datalog {
rule_ref rl(b.m_ctx.get_rule_manager());
rl = rules[i];
b.m_ctx.get_rule_manager().substitute(rl, sub.size(), sub.c_ptr());
substs.push_back(sub);
for (unsigned j = 0; j < sz; ++j) {
if (j == 0) {
path1 = path;
path1 = path;
}
else {
path1 = m.mk_app(succs[j], path);
}
prs.push_back(get_proof(md, to_app(trace->get_arg(j)), path1));
positions.push_back(std::make_pair(j+1,0));
substs.push_back(expr_ref_vector(m));
substs.push_back(expr_ref_vector(m));
}
head = rl->get_head();
pr = m.mk_hyper_resolve(sz+1, prs.c_ptr(), head, positions, substs);
pr = m.mk_hyper_resolve(sz+1, prs.c_ptr(), head, positions, substs);
return pr;
}
}
UNREACHABLE();
return proof_ref(0, m);
}
// instantiation of algebraic data-types takes care of the rest.
lbool check_query() {
sort* trace_sort = m_pred2sort.find(b.m_query_pred);
@ -1102,10 +1102,10 @@ namespace datalog {
}
b.m_solver.get_model(md);
mk_answer(md, trace, path);
return l_true;
return l_true;
}
}
bool check_model(model_ref& md, expr* trace) {
expr_ref trace_val(m), eq(m);
md->eval(trace, trace_val);
@ -1122,10 +1122,10 @@ namespace datalog {
eq = m.mk_not(eq);
b.assert_expr(eq);
}
return is_sat != l_false;
return is_sat != l_false;
}
void mk_answer(model_ref& md, expr_ref& trace, expr_ref& path) {
void mk_answer(model_ref& md, expr_ref& trace, expr_ref& path) {
IF_VERBOSE(2, model_smt2_pp(verbose_stream(), m, *md, 0););
md->eval(trace, trace);
md->eval(path, path);
@ -1139,7 +1139,7 @@ namespace datalog {
apply(m, b.m_ctx.get_proof_converter().get(), pr);
b.m_answer = pr;
}
};
// --------------------------------------------------------------------------
@ -1155,7 +1155,7 @@ namespace datalog {
lbool check() {
setup();
for (unsigned i = 0; ; ++i) {
IF_VERBOSE(1, verbose_stream() << "level: " << i << "\n";);
IF_VERBOSE(1, verbose_stream() << "level: " << i << "\n";);
b.checkpoint();
compile(i);
lbool res = check(i);
@ -1183,9 +1183,9 @@ namespace datalog {
b.m_solver.get_model(md);
func_decl* pred = b.m_query_pred;
SASSERT(m.is_true(md->get_const_interp(to_app(level_query)->get_decl())));
TRACE("bmc", model_smt2_pp(tout, m, *md, 0););
rule_ref r0(rm), r1(rm), r2(rm);
while (true) {
TRACE("bmc", tout << "Predicate: " << pred->get_name() << "\n";);
@ -1193,7 +1193,7 @@ namespace datalog {
rule_vector const& rls = b.m_rules.get_predicate_rules(pred);
rule* r = 0;
unsigned i = 0;
for (; i < rls.size(); ++i) {
for (; i < rls.size(); ++i) {
expr_ref rule_i = mk_level_rule(pred, i, level);
TRACE("bmc", rls[i]->display(b.m_ctx, tout << "Checking rule " << mk_pp(rule_i, m) << " "););
if (m.is_true(md->get_const_interp(to_app(rule_i)->get_decl()))) {
@ -1204,13 +1204,13 @@ namespace datalog {
SASSERT(r);
mk_rule_vars(*r, level, i, sub);
// we have rule, we have variable names of rule.
// extract values for the variables in the rule.
for (unsigned j = 0; j < sub.size(); ++j) {
expr* vl = md->get_const_interp(to_app(sub[j].get())->get_decl());
if (vl) {
// vl can be 0 if the interpretation does not assign a value to it.
sub[j] = vl;
sub[j] = vl;
}
else {
sub[j] = m.mk_var(j, m.get_sort(sub[j].get()));
@ -1219,7 +1219,7 @@ namespace datalog {
svector<std::pair<unsigned, unsigned> > positions;
vector<expr_ref_vector> substs;
expr_ref fml(m), concl(m);
rm.to_formula(*r, fml);
r2 = r;
rm.substitute(r2, sub.size(), sub.c_ptr());
@ -1239,12 +1239,12 @@ namespace datalog {
apply_subst(sub, sub2);
unifier.apply(*r0.get(), 0, *r2.get(), r1);
rm.to_formula(*r1.get(), concl);
scoped_proof _sp(m);
proof* premises[2] = { pr, p };
positions.push_back(std::make_pair(0, 1));
substs.push_back(sub1);
substs.push_back(sub);
pr = m.mk_hyper_resolve(2, premises, concl, positions, substs);
@ -1263,7 +1263,7 @@ namespace datalog {
}
r0 = r2;
}
if (level == 0) {
SASSERT(r->get_uninterpreted_tail_size() == 0);
break;
@ -1276,8 +1276,8 @@ namespace datalog {
apply(m, b.m_ctx.get_proof_converter().get(), pr);
b.m_answer = pr;
}
void setup() {
b.m_fparams.m_relevancy_lvl = 0;
b.m_fparams.m_model = true;
@ -1285,54 +1285,54 @@ namespace datalog {
b.m_fparams.m_mbqi = false;
// m_fparams.m_auto_config = false;
}
lbool check(unsigned level) {
expr_ref level_query = mk_level_predicate(b.m_query_pred, level);
expr* q = level_query.get();
return b.m_solver.check(1, &q);
}
expr_ref mk_level_predicate(func_decl* p, unsigned level) {
return mk_level_predicate(p->get_name(), level);
}
expr_ref mk_level_predicate(symbol const& name, unsigned level) {
std::stringstream _name;
_name << name << "#" << level;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return expr_ref(m.mk_const(nm, m.mk_bool_sort()), m);
}
expr_ref mk_level_arg(func_decl* pred, unsigned idx, unsigned level) {
SASSERT(idx < pred->get_arity());
std::stringstream _name;
_name << pred->get_name() << "#" << level << "_" << idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return expr_ref(m.mk_const(nm, pred->get_domain(idx)), m);
}
expr_ref mk_level_var(func_decl* pred, sort* s, unsigned rule_id, unsigned idx, unsigned level) {
std::stringstream _name;
_name << pred->get_name() << "#" << level << "_" << rule_id << "_" << idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return expr_ref(m.mk_const(nm, s), m);
}
expr_ref mk_level_rule(func_decl* p, unsigned rule_idx, unsigned level) {
std::stringstream _name;
_name << "rule:" << p->get_name() << "#" << level << "_" << rule_idx;
symbol nm(_name.str().c_str());
symbol nm(_name.str().c_str());
return expr_ref(m.mk_const(nm, m.mk_bool_sort()), m);
}
void mk_rule_vars(rule& r, unsigned level, unsigned rule_id, expr_ref_vector& sub) {
ptr_vector<sort> sorts;
r.get_vars(m, sorts);
// populate substitution of bound variables.
sub.reset();
sub.resize(sorts.size());
for (unsigned k = 0; k < r.get_decl()->get_arity(); ++k) {
expr* arg = r.get_head()->get_arg(k);
if (is_var(arg)) {
@ -1361,34 +1361,34 @@ namespace datalog {
}
}
}
void compile(unsigned level) {
rule_set::decl2rules::iterator it = b.m_rules.begin_grouped_rules();
rule_set::decl2rules::iterator end = b.m_rules.end_grouped_rules();
for (; it != end; ++it) {
func_decl* p = it->m_key;
rule_vector const& rls = *it->m_value;
// Assert: p_level => r1_level \/ r2_level \/ r3_level \/ ...
// Assert: r_i_level => body of rule i for level + equalities for head of rule i
expr_ref level_pred = mk_level_predicate(p, level);
expr_ref_vector rules(m), sub(m), conjs(m);
expr_ref rule_body(m), tmp(m);
for (unsigned i = 0; i < rls.size(); ++i) {
sub.reset();
for (unsigned i = 0; i < rls.size(); ++i) {
sub.reset();
conjs.reset();
rule& r = *rls[i];
expr_ref rule_i = mk_level_rule(p, i, level);
rule& r = *rls[i];
expr_ref rule_i = mk_level_rule(p, i, level);
rules.push_back(rule_i);
if (level == 0 && r.get_uninterpreted_tail_size() > 0) {
tmp = m.mk_not(rule_i);
b.assert_expr(tmp);
continue;
}
mk_rule_vars(r, level, i, sub);
// apply substitution to body.
var_subst vs(m, false);
for (unsigned k = 0; k < p->get_arity(); ++k) {
@ -1418,11 +1418,11 @@ namespace datalog {
}
}
};
bmc::bmc(context& ctx):
bmc::bmc(context& ctx):
engine_base(ctx.get_manager(), "bmc"),
m_ctx(ctx),
m(ctx.get_manager()),
m_ctx(ctx),
m(ctx.get_manager()),
m_solver(m, m_fparams),
m_rules(ctx),
m_query_pred(m),
@ -1440,9 +1440,9 @@ namespace datalog {
rule_set& rules0 = m_ctx.get_rules();
datalog::rule_set old_rules(rules0);
rule_manager.mk_query(query, rules0);
expr_ref bg_assertion = m_ctx.get_background_assertion();
expr_ref bg_assertion = m_ctx.get_background_assertion();
apply_default_transformation(m_ctx);
if (m_ctx.xform_slice()) {
datalog::rule_transformer transformer(m_ctx);
datalog::mk_slice* slice = alloc(datalog::mk_slice, m_ctx);
@ -1490,7 +1490,7 @@ namespace datalog {
}
}
void bmc::assert_expr(expr* e) {
void bmc::assert_expr(expr* e) {
TRACE("bmc", tout << mk_pp(e, m) << "\n";);
m_solver.assert_expr(e);
}
@ -1510,7 +1510,7 @@ namespace datalog {
void bmc::checkpoint() {
if (m.canceled()) {
throw default_exception("bmc canceled");
throw default_exception(Z3_CANCELED_MSG);
}
}
@ -1526,7 +1526,7 @@ namespace datalog {
m_solver.reset_statistics();
}
expr_ref bmc::get_answer() {
expr_ref bmc::get_answer() {
return m_answer;
}

52
src/muz/duality/duality_dl_interface.cpp
View file

@ -62,7 +62,7 @@ namespace Duality {
context ctx;
RPFP::LogicSolver *ls;
RPFP *rpfp;
DualityStatus status;
std::vector<expr> clauses;
std::vector<std::vector<RPFP::label_struct> > clause_labels;
@ -104,14 +104,14 @@ namespace Duality {
//
// Check if the new rules are weaker so that we can
// Check if the new rules are weaker so that we can
// re-use existing context.
//
//
#if 0
void dl_interface::check_reset() {
// TODO
datalog::rule_ref_vector const& new_rules = m_ctx.get_rules().get_rules();
datalog::rule_ref_vector const& old_rules = m_old_rules.get_rules();
datalog::rule_ref_vector const& old_rules = m_old_rules.get_rules();
bool is_subsumed = !old_rules.empty();
for (unsigned i = 0; is_subsumed && i < new_rules.size(); ++i) {
is_subsumed = false;
@ -155,7 +155,7 @@ namespace Duality {
expr_ref_vector rules(m_ctx.get_manager());
svector< ::symbol> names;
svector< ::symbol> names;
unsigned_vector bounds;
// m_ctx.get_rules_as_formulas(rules, names);
@ -189,7 +189,7 @@ namespace Duality {
rules.push_back(f);
}
}
else
else
m_ctx.get_raw_rule_formulas(rules, names, bounds);
// get all the rules as clauses
@ -199,7 +199,7 @@ namespace Duality {
expr e(_d->ctx,rules[i].get());
clauses.push_back(e);
}
std::vector<sort> b_sorts;
std::vector<symbol> b_names;
used_vars uv;
@ -216,7 +216,7 @@ namespace Duality {
#if 0
// turn the query into a clause
expr q(_d->ctx,m_ctx.bind_variables(query,false));
std::vector<sort> b_sorts;
std::vector<symbol> b_names;
if (q.is_quantifier() && !q.is_quantifier_forall()) {
@ -235,14 +235,14 @@ namespace Duality {
qc = _d->ctx.make_quant(Forall,b_sorts,b_names,qc);
clauses.push_back(qc);
bounds.push_back(UINT_MAX);
// get the background axioms
unsigned num_asserts = m_ctx.get_num_assertions();
for (unsigned i = 0; i < num_asserts; ++i) {
expr e(_d->ctx,m_ctx.get_assertion(i));
_d->rpfp->AssertAxiom(e);
}
// make sure each predicate is the head of at least one clause
func_decl_set heads;
for(unsigned i = 0; i < clauses.size(); i++){
@ -297,14 +297,14 @@ namespace Duality {
// populate the edge-to-clause map
for(unsigned i = 0; i < _d->rpfp->edges.size(); ++i)
_d->map[_d->rpfp->edges[i]] = i;
// create a solver object
Solver *rs = Solver::Create("duality", _d->rpfp);
if(old_rs)
rs->LearnFrom(old_rs); // new solver gets hints from old solver
// set its options
IF_VERBOSE(1, rs->SetOption("report","1"););
rs->SetOption("full_expand",m_ctx.get_params().duality_full_expand() ? "1" : "0");
@ -328,12 +328,12 @@ namespace Duality {
ans = rs->Solve();
}
catch (Duality::solver::cancel_exception &exn){
throw default_exception("duality canceled");
throw default_exception(Z3_CANCELED_MSG);
}
catch (Duality::Solver::Incompleteness &exn){
throw default_exception("incompleteness");
}
// profile!
if(m_ctx.get_params().duality_profile())
@ -343,7 +343,7 @@ namespace Duality {
_d->status = ans ? StatusModel : StatusRefutation;
_d->cex.swap(rs->GetCounterexample()); // take ownership of cex
_d->old_rs = rs; // save this for later hints
if(old_data){
dealloc(old_data); // this deallocates the old solver if there is one
}
@ -392,7 +392,7 @@ namespace Duality {
context &ctx = d->dd()->ctx;
RPFP::Node &node = *root;
RPFP::Edge &edge = *node.Outgoing;
// first, prove the children (that are actually used)
for(unsigned i = 0; i < edge.Children.size(); i++){
@ -434,19 +434,19 @@ namespace Duality {
}
}
out << " )\n";
out << " (labels";
std::vector<symbol> labels;
tree->GetLabels(&edge,labels);
for(unsigned j = 0; j < labels.size(); j++){
out << " " << labels[j];
}
out << " )\n";
// reference the proofs of all the children, in syntactic order
// "true" means the child is not needed
out << " (ref ";
for(unsigned i = 0; i < edge.Children.size(); i++){
if(!tree->Empty(edge.Children[i]))
@ -468,7 +468,7 @@ namespace Duality {
ast_manager &m = m_ctx.get_manager();
model_ref md = get_model();
out << "(fixedpoint \n";
model_smt2_pp(out, m, *md.get(), 0);
model_smt2_pp(out, m, *md.get(), 0);
out << ")\n";
}
else if(_d->status == StatusRefutation){
@ -550,22 +550,22 @@ namespace Duality {
}
return md;
}
static proof_ref extract_proof(dl_interface *d, RPFP *tree, RPFP::Node *root) {
context &ctx = d->dd()->ctx;
ast_manager &mgr = ctx.m();
RPFP::Node &node = *root;
RPFP::Edge &edge = *node.Outgoing;
RPFP::Edge *orig_edge = edge.map;
// first, prove the children (that are actually used)
proof_ref_vector prems(mgr);
::vector<expr_ref_vector> substs;
int orig_clause = d->dd()->map[orig_edge];
expr &t = d->dd()->clauses[orig_clause];
prems.push_back(mgr.mk_asserted(ctx.uncook(t)));
substs.push_back(expr_ref_vector(mgr));
if (t.is_quantifier() && t.is_quantifier_forall()) {
int bound = t.get_quantifier_num_bound();
@ -599,12 +599,12 @@ namespace Duality {
for(unsigned i = 0; i < edge.F.IndParams.size(); i++)
args.push_back(tree->Eval(&edge,edge.F.IndParams[i]));
expr conc = f(args);
::vector< ::proof *> pprems;
for(unsigned i = 0; i < prems.size(); i++)
pprems.push_back(prems[i].get());
proof_ref res(mgr.mk_hyper_resolve(pprems.size(),&pprems[0], ctx.uncook(conc), pos, substs),mgr);
return res;

390
src/muz/pdr/pdr_context.cpp
View file

File diff suppressed because it is too large Load diff

40
src/smt/tactic/smt_tactic.cpp
View file

@ -59,7 +59,7 @@ void extract_clauses_and_dependencies(goal_ref const& g, expr_ref_vector& clause
clause.push_back(m.mk_not(d));
}
else {
// must normalize assumption
// must normalize assumption
expr * b = 0;
if (!dep2bool.find(d, b)) {
b = m.mk_fresh_const(0, m.mk_bool_sort());
@ -79,7 +79,7 @@ void extract_clauses_and_dependencies(goal_ref const& g, expr_ref_vector& clause
cls = mk_or(m, clause.size(), clause.c_ptr());
clauses.push_back(cls);
}
}
}
}
class smt_tactic : public tactic {
@ -96,7 +96,7 @@ class smt_tactic : public tactic {
public:
smt_tactic(params_ref const & p):
m_params_ref(p),
m_ctx(0),
m_ctx(0),
m_callback(0) {
updt_params_core(p);
TRACE("smt_tactic", tout << this << "\np: " << p << "\n";);
@ -118,20 +118,20 @@ public:
m_candidate_models = p.get_bool("candidate_models", false);
m_fail_if_inconclusive = p.get_bool("fail_if_inconclusive", true);
}
virtual void updt_params(params_ref const & p) {
TRACE("smt_tactic", tout << this << "\nupdt_params: " << p << "\n";);
updt_params_core(p);
fparams().updt_params(p);
SASSERT(p.get_bool("auto_config", fparams().m_auto_config) == fparams().m_auto_config);
}
virtual void collect_param_descrs(param_descrs & r) {
r.insert("candidate_models", CPK_BOOL, "(default: false) create candidate models even when quantifier or theory reasoning is incomplete.");
r.insert("fail_if_inconclusive", CPK_BOOL, "(default: true) fail if found unsat (sat) for under (over) approximated goal.");
smt_params_helper::collect_param_descrs(r);
}
virtual void collect_statistics(statistics & st) const {
if (m_ctx)
@ -146,7 +146,7 @@ public:
virtual void reset_statistics() {
m_stats.reset();
}
virtual void set_logic(symbol const & l) {
m_logic = l;
}
@ -154,7 +154,7 @@ public:
virtual void set_progress_callback(progress_callback * callback) {
m_callback = callback;
}
struct scoped_init_ctx {
smt_tactic & m_owner;
smt_params m_params; // smt-setup overwrites parameters depending on the current assertions.
@ -168,41 +168,41 @@ public:
new_ctx->set_progress_callback(o.m_callback);
}
o.m_ctx = new_ctx;
}
~scoped_init_ctx() {
smt::kernel * d = m_owner.m_ctx;
m_owner.m_ctx = 0;
if (d)
dealloc(d);
}
};
virtual void operator()(goal_ref const & in,
goal_ref_buffer & result,
model_converter_ref & mc,
virtual void operator()(goal_ref const & in,
goal_ref_buffer & result,
model_converter_ref & mc,
proof_converter_ref & pc,
expr_dependency_ref & core) {
try {
SASSERT(in->is_well_sorted());
ast_manager & m = in->m();
TRACE("smt_tactic", tout << this << "\nAUTO_CONFIG: " << fparams().m_auto_config << " HIDIV0: " << fparams().m_hi_div0 << " "
TRACE("smt_tactic", tout << this << "\nAUTO_CONFIG: " << fparams().m_auto_config << " HIDIV0: " << fparams().m_hi_div0 << " "
<< " PREPROCESS: " << fparams().m_preprocess << "\n";
tout << "RELEVANCY: " << fparams().m_relevancy_lvl << "\n";
tout << "fail-if-inconclusive: " << m_fail_if_inconclusive << "\n";
tout << "params_ref: " << m_params_ref << "\n";
tout << "nnf: " << fparams().m_nnf_cnf << "\n";);
TRACE("smt_tactic_detail", in->display(tout););
TRACE("smt_tactic_memory", tout << "wasted_size: " << m.get_allocator().get_wasted_size() << "\n";);
TRACE("smt_tactic_memory", tout << "wasted_size: " << m.get_allocator().get_wasted_size() << "\n";);
scoped_init_ctx init(*this, m);
SASSERT(m_ctx != 0);
expr_ref_vector clauses(m);
expr2expr_map bool2dep;
ptr_vector<expr> assumptions;
expr2expr_map bool2dep;
ptr_vector<expr> assumptions;
ref<filter_model_converter> fmc;
if (in->unsat_core_enabled()) {
extract_clauses_and_dependencies(in, clauses, assumptions, bool2dep, fmc);
@ -225,9 +225,9 @@ public:
}
}
if (m_ctx->canceled()) {
throw tactic_exception("smt_tactic canceled");
throw tactic_exception(Z3_CANCELED_MSG);
}
lbool r;
if (assumptions.empty())
r = m_ctx->setup_and_check();
@ -315,7 +315,7 @@ tactic * mk_smt_tactic(params_ref const & p) {
}
tactic * mk_smt_tactic_using(bool auto_config, params_ref const & _p) {
params_ref p = _p;
params_ref p = _p;
p.set_bool("auto_config", auto_config);
tactic * r = mk_smt_tactic(p);
TRACE("smt_tactic", tout << "auto_config: " << auto_config << "\nr: " << r << "\np: " << p << "\n";);

2
src/tactic/nlsat_smt/nl_purify_tactic.cpp
View file

@ -289,7 +289,7 @@ private:
void check_point() {
if (m.canceled()) {
throw tactic_exception("canceled");
throw tactic_exception(Z3_CANCELED_MSG);
}
}