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

fix #2613 fix #2612

Browse source 
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2019-10-05 16:57:51 -07:00
parent 016732aa59
commit 39edf73e78
16 changed files with 222 additions and 162 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -257,7 +257,7 @@ namespace algebraic_numbers {
SASSERT(bqm().ge(upper(c), candidate));
if (bqm().lt(lower(c), candidate) && upm().eval_sign_at(c->m_p_sz, c->m_p, candidate) == 0) {
if (bqm().lt(lower(c), candidate) && upm().eval_sign_at(c->m_p_sz, c->m_p, candidate) == polynomial::sign_zero) {
m_wrapper.set(a, candidate);
return true;
}
@ -320,7 +320,7 @@ namespace algebraic_numbers {
SASSERT(bqm().ge(upper(c), candidate));
// Find if candidate is an actual root
if (bqm().lt(lower(c), candidate) && upm().eval_sign_at(c->m_p_sz, c->m_p, candidate) == 0) {
if (bqm().lt(lower(c), candidate) && upm().eval_sign_at(c->m_p_sz, c->m_p, candidate) == polynomial::sign_zero) {
saved_a.restore_if_too_small();
set(a, candidate);
return true;
@ -379,11 +379,11 @@ namespace algebraic_numbers {
}
void update_sign_lower(algebraic_cell * c) {
int sl = upm().eval_sign_at(c->m_p_sz, c->m_p, lower(c));
polynomial::sign sl = upm().eval_sign_at(c->m_p_sz, c->m_p, lower(c));
// The isolating intervals are refinable. Thus, the polynomial has opposite signs at lower and upper.
SASSERT(sl != 0);
SASSERT(sl != polynomial::sign_zero);
SASSERT(upm().eval_sign_at(c->m_p_sz, c->m_p, upper(c)) == -sl);
c->m_sign_lower = sl < 0;
c->m_sign_lower = sl == polynomial::sign_neg;
}
// Make sure the GCD of the coefficients is one and the leading coefficient is positive
@ -1594,8 +1594,8 @@ namespace algebraic_numbers {
#define REFINE_LOOP(BOUND, TARGET_SIGN) \
while (true) { \
bqm().div2(BOUND); \
int new_sign = upm().eval_sign_at(cell_a->m_p_sz, cell_a->m_p, BOUND); \
if (new_sign == 0) { \
polynomial::sign new_sign = upm().eval_sign_at(cell_a->m_p_sz, cell_a->m_p, BOUND); \
if (new_sign == polynomial::sign_zero) { \
/* found actual root */ \
scoped_mpq r(qm()); \
to_mpq(qm(), BOUND, r); \
@ -1695,8 +1695,8 @@ namespace algebraic_numbers {
if (bqm().ge(l, b))
return 1;
// b is in the isolating interval (l, u)
int sign_b = upm().eval_sign_at(c->m_p_sz, c->m_p, b);
if (sign_b == 0)
auto sign_b = upm().eval_sign_at(c->m_p_sz, c->m_p, b);
if (sign_b == polynomial::sign_zero)
return 0;
return sign_b == sign_lower(c) ? 1 : -1;
}
@ -1979,7 +1979,7 @@ namespace algebraic_numbers {
};
polynomial::var_vector m_eval_sign_vars;
int eval_sign_at(polynomial_ref const & p, polynomial::var2anum const & x2v) {
polynomial::sign eval_sign_at(polynomial_ref const & p, polynomial::var2anum const & x2v) {
polynomial::manager & ext_pm = p.m();
TRACE("anum_eval_sign", tout << "evaluating sign of: " << p << "\n";);
while (true) {
@ -1990,7 +1990,7 @@ namespace algebraic_numbers {
scoped_mpq r(qm());
ext_pm.eval(p, x2v_basic, r);
TRACE("anum_eval_sign", tout << "all variables are assigned to rationals, value of p: " << r << "\n";);
return qm().sign(r);
return polynomial::to_sign(qm().sign(r));
}
catch (const opt_var2basic::failed &) {
// continue
@ -2004,13 +2004,13 @@ namespace algebraic_numbers {
if (ext_pm.is_zero(p_prime)) {
// polynomial vanished after substituting rational values.
return 0;
return polynomial::sign_zero;
}
if (is_const(p_prime)) {
// polynomial became the constant polynomial after substitution.
SASSERT(size(p_prime) == 1);
return ext_pm.m().sign(ext_pm.coeff(p_prime, 0));
return polynomial::to_sign(ext_pm.m().sign(ext_pm.coeff(p_prime, 0)));
}
// Try to find sign using intervals
@ -2026,7 +2026,7 @@ namespace algebraic_numbers {
ext_pm.eval(p_prime, x2v_interval, ri);
TRACE("anum_eval_sign", tout << "evaluating using intervals: " << ri << "\n";);
if (!bqim().contains_zero(ri)) {
return bqim().is_pos(ri) ? 1 : -1;
return bqim().is_pos(ri) ? polynomial::sign_pos : polynomial::sign_neg;
}
// refine intervals if magnitude > m_min_magnitude
bool refined = false;
@ -2067,7 +2067,7 @@ namespace algebraic_numbers {
// Remark: m_zero_accuracy == 0 means use precise computation.
if (m_zero_accuracy > 0) {
// assuming the value is 0, since the result is in (-1/2^k, 1/2^k), where m_zero_accuracy = k
return 0;
return polynomial::sign_zero;
}
#if 0
// Evaluating sign using algebraic arithmetic
@ -2143,7 +2143,7 @@ namespace algebraic_numbers {
bqm().div2k(mL, k);
bqm().div2k(L, k);
if (bqm().lt(mL, ri.lower()) && bqm().lt(ri.upper(), L))
return 0;
return polynomial::sign_zero;
// keep refining ri until ri is inside (-L, L) or
// ri does not contain zero.
@ -2166,14 +2166,13 @@ namespace algebraic_numbers {
TRACE("anum_eval_sign", tout << "evaluating using intervals: " << ri << "\n";
tout << "zero lower bound is: " << L << "\n";);
if (!bqim().contains_zero(ri)) {
return bqim().is_pos(ri) ? 1 : -1;
return bqim().is_pos(ri) ? polynomial::sign_pos : polynomial::sign_neg;
}
if (bqm().lt(mL, ri.lower()) && bqm().lt(ri.upper(), L))
return 0;
return polynomial::sign_zero;
for (unsigned i = 0; i < xs.size(); i++) {
polynomial::var x = xs[i];
for (auto x : xs) {
SASSERT(x2v.contains(x));
anum const & v = x2v(x);
SASSERT(!v.is_basic());
@ -2242,18 +2241,14 @@ namespace algebraic_numbers {
unsigned sz = roots.size();
unsigned j = 0;
// std::cout << "p: " << p << "\n";
// std::cout << "sz: " << sz << "\n";
for (unsigned i = 0; i < sz; i++) {
checkpoint();
// display_root(std::cout, roots[i]); std::cout << std::endl;
ext_var2num ext_x2v(m_wrapper, x2v, x, roots[i]);
TRACE("isolate_roots", tout << "filter_roots i: " << i << ", ext_x2v: x" << x << " -> "; display_root(tout, roots[i]); tout << "\n";);
int sign = eval_sign_at(p, ext_x2v);
polynomial::sign sign = eval_sign_at(p, ext_x2v);
TRACE("isolate_roots", tout << "filter_roots i: " << i << ", result sign: " << sign << "\n";);
if (sign != 0)
continue;
// display_decimal(std::cout, roots[i], 10); std::cout << " is root" << std::endl;
if (i != j)
set(roots[j], roots[i]);
j++;
@ -2453,7 +2448,7 @@ namespace algebraic_numbers {
}
}
int eval_at_mpbq(polynomial_ref const & p, polynomial::var2anum const & x2v, polynomial::var x, mpbq const & v) {
polynomial::sign eval_at_mpbq(polynomial_ref const & p, polynomial::var2anum const & x2v, polynomial::var x, mpbq const & v) {
scoped_mpq qv(qm());
to_mpq(qm(), v, qv);
scoped_anum av(m_wrapper);
@ -2568,13 +2563,13 @@ namespace algebraic_numbers {
#define DEFAULT_PRECISION 2
void isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, svector<int> & signs) {
void isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, svector<polynomial::sign> & signs) {
isolate_roots(p, x2v, roots);
unsigned num_roots = roots.size();
if (num_roots == 0) {
anum zero;
ext2_var2num ext_x2v(m_wrapper, x2v, zero);
int s = eval_sign_at(p, ext_x2v);
polynomial::sign s = eval_sign_at(p, ext_x2v);
signs.push_back(s);
}
else {
@ -2601,8 +2596,8 @@ namespace algebraic_numbers {
TRACE("isolate_roots_bug", tout << "w: "; display_root(tout, w); tout << "\n";);
{
ext2_var2num ext_x2v(m_wrapper, x2v, w);
int s = eval_sign_at(p, ext_x2v);
SASSERT(s != 0);
auto s = eval_sign_at(p, ext_x2v);
SASSERT(s != polynomial::sign_zero);
signs.push_back(s);
}
@ -2611,16 +2606,16 @@ namespace algebraic_numbers {
numeral & curr = roots[i];
select(prev, curr, w);
ext2_var2num ext_x2v(m_wrapper, x2v, w);
int s = eval_sign_at(p, ext_x2v);
SASSERT(s != 0);
auto s = eval_sign_at(p, ext_x2v);
SASSERT(s != polynomial::sign_zero);
signs.push_back(s);
}
int_gt(roots[num_roots - 1], w);
{
ext2_var2num ext_x2v(m_wrapper, x2v, w);
int s = eval_sign_at(p, ext_x2v);
SASSERT(s != 0);
auto s = eval_sign_at(p, ext_x2v);
SASSERT(s != polynomial::sign_zero);
signs.push_back(s);
}
}
@ -2879,7 +2874,7 @@ namespace algebraic_numbers {
m_imp->isolate_roots(p, x2v, roots);
}
void manager::isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, svector<int> & signs) {
void manager::isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, svector<polynomial::sign> & signs) {
m_imp->isolate_roots(p, x2v, roots, signs);
}
@ -3037,7 +3032,7 @@ namespace algebraic_numbers {
l = rational(_l);
}
int manager::eval_sign_at(polynomial_ref const & p, polynomial::var2anum const & x2v) {
polynomial::sign manager::eval_sign_at(polynomial_ref const & p, polynomial::var2anum const & x2v) {
SASSERT(&(x2v.m()) == this);
return m_imp->eval_sign_at(p, x2v);
}

4
src/math/polynomial/algebraic_numbers.h
View file

@ -173,7 +173,7 @@ namespace algebraic_numbers {
/**
\brief Isolate the roots of the given polynomial, and compute its sign between them.
*/
void isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, svector<int> & signs);
void isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, svector<polynomial::sign> & signs);
/**
\brief Store in r the i-th root of p.
@ -304,7 +304,7 @@ namespace algebraic_numbers {
Return 0 if p(alpha_1, ..., alpha_n) == 0
Return positive number if p(alpha_1, ..., alpha_n) > 0
*/
int eval_sign_at(polynomial_ref const & p, polynomial::var2anum const & x2v);
polynomial::sign eval_sign_at(polynomial_ref const & p, polynomial::var2anum const & x2v);
void get_polynomial(numeral const & a, svector<mpz> & r);

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

@ -1068,12 +1068,12 @@ namespace polynomial {
g.reserve(std::min(sz1, sz2));
r1.reserve(sz2); // r1 has at most num_args2 arguments
r2.reserve(sz1); // r2 has at most num_args1 arguments
bool found = false;
unsigned i1 = 0;
unsigned i2 = 0;
unsigned j1 = 0;
unsigned j2 = 0;
unsigned j3 = 0;
bool found = false;
unsigned i1 = 0;
unsigned i2 = 0;
unsigned j1 = 0;
unsigned j2 = 0;
unsigned j3 = 0;
while (true) {
if (i1 == sz1) {
if (found) {
@ -2501,6 +2501,32 @@ namespace polynomial {
return p;
}
void gcd_simplify(polynomial * p) {
if (m_manager.finite()) return;
auto& m = m_manager.m();
unsigned sz = p->size();
if (sz == 0)
return;
unsigned g = 0;
for (unsigned i = 0; i < sz; i++) {
if (!m.is_int(p->a(i))) {
return;
}
int j = m.get_int(p->a(i));
if (j == INT_MIN || j == 1 || j == -1)
return;
g = u_gcd(abs(j), g);
if (g == 1)
return;
}
scoped_mpz r(m), gg(m);
m.set(gg, g);
for (unsigned i = 0; i < sz; ++i) {
m.div_gcd(p->a(i), gg, r);
m.set(p->a(i), r);
}
}
polynomial * mk_zero() {
return m_zero;
}
@ -7041,6 +7067,10 @@ namespace polynomial {
return m_imp->hash(p);
}
void manager::gcd_simplify(polynomial * p) {
m_imp->gcd_simplify(p);
}
polynomial * manager::coeff(polynomial const * p, var x, unsigned k) {
return m_imp->coeff(p, x, k);
}

11
src/math/polynomial/polynomial.h
View file

@ -44,6 +44,11 @@ namespace polynomial {
typedef svector<var> var_vector;
class monomial;
typedef enum { sign_neg = -1, sign_zero = 0, sign_pos = 1} sign;
inline sign operator-(sign s) { switch (s) { case sign_neg: return sign_pos; case sign_pos: return sign_neg; default: return sign_zero; } };
inline sign to_sign(int s) { return s == 0 ? sign_zero : (s > 0 ? sign_pos : sign_neg); }
inline sign operator*(sign a, sign b) { return to_sign((int)a * (int)b); }
int lex_compare(monomial const * m1, monomial const * m2);
int lex_compare2(monomial const * m1, monomial const * m2, var min_var);
int graded_lex_compare(monomial const * m1, monomial const * m2);
@ -278,6 +283,12 @@ namespace polynomial {
*/
static unsigned id(polynomial const * p);
/**
\brief Normalize coefficients by dividing by their gcd
*/
void gcd_simplify(polynomial* p);
/**
\brief Return true if \c m is the unit monomial.
*/

1
src/math/polynomial/polynomial_cache.cpp
View file

@ -139,6 +139,7 @@ namespace polynomial {
polynomial * mk_unique(polynomial * p) {
if (m_in_cache.get(pid(p), false))
return p;
// m.gcd_simplify(p);
polynomial * p_prime = m_poly_table.insert_if_not_there(p);
if (p == p_prime) {
m_cached_polys.push_back(p_prime);

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

@ -1327,25 +1327,25 @@ namespace upolynomial {
div(sz, p, 2, two_x_1, buffer);
}
int manager::sign_of(numeral const & c) {
polynomial::sign manager::sign_of(numeral const & c) {
if (m().is_zero(c))
return 0;
return polynomial::sign_zero;
if (m().is_pos(c))
return 1;
return polynomial::sign_pos;
else
return -1;
return polynomial::sign_neg;
}
// Return the number of sign changes in the coefficients of p
unsigned manager::sign_changes(unsigned sz, numeral const * p) {
unsigned r = 0;
int prev_sign = 0;
auto prev_sign = polynomial::sign_zero;
unsigned i = 0;
for (; i < sz; i++) {
int sign = sign_of(p[i]);
if (sign == 0)
auto sign = sign_of(p[i]);
if (sign == polynomial::sign_zero)
continue;
if (sign != prev_sign && prev_sign != 0)
if (sign != prev_sign && prev_sign != polynomial::sign_zero)
r++;
prev_sign = sign;
}
@ -1729,14 +1729,14 @@ namespace upolynomial {
}
// Evaluate the sign of p(b)
int manager::eval_sign_at(unsigned sz, numeral const * p, mpbq const & b) {
polynomial::sign 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
if (sz == 0)
return 0;
return polynomial::sign_zero;
if (sz == 1)
return sign_of(p[0]);
numeral const & c = b.numerator();
@ -1762,14 +1762,14 @@ namespace upolynomial {
}
// Evaluate the sign of p(b)
int manager::eval_sign_at(unsigned sz, numeral const * p, mpq const & b) {
polynomial::sign manager::eval_sign_at(unsigned sz, numeral const * p, mpq const & b) {
// Actually, given b = c/d, we compute the sign of (d^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 (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
if (sz == 0)
return 0;
return polynomial::sign_zero;
if (sz == 1)
return sign_of(p[0]);
numeral const & c = b.numerator();
@ -1796,11 +1796,11 @@ namespace upolynomial {
}
// Evaluate the sign of p(b)
int manager::eval_sign_at(unsigned sz, numeral const * p, mpz const & b) {
polynomial::sign manager::eval_sign_at(unsigned sz, numeral const * p, mpz const & b) {
// Using Horner Sequence
// ((a_n * b + a_{n-1})*b + a_{n-2})*b + a_{n-3} ...
if (sz == 0)
return 0;
return polynomial::sign_zero;
if (sz == 1)
return sign_of(p[0]);
scoped_numeral r(m());
@ -1817,21 +1817,21 @@ namespace upolynomial {
return sign_of(r);
}
int manager::eval_sign_at_zero(unsigned sz, numeral const * p) {
polynomial::sign manager::eval_sign_at_zero(unsigned sz, numeral const * p) {
if (sz == 0)
return 0;
return polynomial::sign_zero;
return sign_of(p[0]);
}
int manager::eval_sign_at_plus_inf(unsigned sz, numeral const * p) {
polynomial::sign manager::eval_sign_at_plus_inf(unsigned sz, numeral const * p) {
if (sz == 0)
return 0;
return polynomial::sign_zero;
return sign_of(p[sz-1]);
}
int manager::eval_sign_at_minus_inf(unsigned sz, numeral const * p) {
polynomial::sign manager::eval_sign_at_minus_inf(unsigned sz, numeral const * p) {
if (sz == 0)
return 0;
return polynomial::sign_zero;
unsigned degree = sz - 1;
if (degree % 2 == 0)
return sign_of(p[sz-1]);

22
src/math/polynomial/upolynomial.h
View file

@ -554,7 +554,7 @@ namespace upolynomial {
numeral_vector m_tr_tmp;
numeral_vector m_push_tmp;
int sign_of(numeral const & c);
polynomial::sign sign_of(numeral const & c);
struct drs_frame;
void pop_top_frame(numeral_vector & p_stack, svector<drs_frame> & frame_stack);
void push_child_frames(unsigned sz, numeral const * p, numeral_vector & p_stack, svector<drs_frame> & frame_stack);
@ -735,32 +735,32 @@ namespace upolynomial {
/**
\brief Evaluate the sign of p(b)
*/
int eval_sign_at(unsigned sz, numeral const * p, mpbq const & b);
polynomial::sign eval_sign_at(unsigned sz, numeral const * p, mpbq const & b);
/**
\brief Evaluate the sign of p(b)
*/
polynomial::sign eval_sign_at(unsigned sz, numeral const * p, mpq const & b);
/**
\brief Evaluate the sign of p(b)
*/
int eval_sign_at(unsigned sz, numeral const * p, mpq const & b);
/**
\brief Evaluate the sign of p(b)
*/
int eval_sign_at(unsigned sz, numeral const * p, mpz const & b);
polynomial::sign eval_sign_at(unsigned sz, numeral const * p, mpz const & b);
/**
\brief Evaluate the sign of p(0)
*/
int eval_sign_at_zero(unsigned sz, numeral const * p);
polynomial::sign eval_sign_at_zero(unsigned sz, numeral const * p);
/**
\brief Evaluate the sign of p(+oo)
*/
int eval_sign_at_plus_inf(unsigned sz, numeral const * p);
polynomial::sign eval_sign_at_plus_inf(unsigned sz, numeral const * p);
/**
\brief Evaluate the sign of p(-oo)
*/
int eval_sign_at_minus_inf(unsigned sz, numeral const * p);
polynomial::sign eval_sign_at_minus_inf(unsigned sz, numeral const * p);
/**
\brief Evaluate the sign variations in the polynomial sequence at -oo