mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-23 17:15:31 +00:00
Code Activity

"canceled" -> Z3_CANCELED_MSG

Browse source 
Relates to #431
This commit is contained in:
Christoph M. Wintersteiger 2016-02-04 13:52:43 +00:00
parent 9b979b6e1e
commit 4e37821dde
13 changed files with 1632 additions and 1625 deletions
Download patch file Download diff file Expand all files Collapse all files

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