Skip to content

Commit 7c10fb8

Browse files
NikolajBjornerNikolajBjorner
committed
fix #2615
Signed-off-by: Nikolaj Bjorner <[email protected]>
1 parent f9b6e4e commit 7c10fb8

File tree

4 files changed

+92
-78
lines changed

4 files changed

+92
-78
lines changed

‎src/math/polynomial/algebraic_numbers.cpp

Lines changed: 64 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -16,16 +16,16 @@ Module Name:
1616
Notes:
1717
1818
--*/
19-
#include "math/polynomial/algebraic_numbers.h"
20-
#include "math/polynomial/upolynomial.h"
2119
#include "util/mpbq.h"
2220
#include "util/basic_interval.h"
23-
#include "math/polynomial/sexpr2upolynomial.h"
2421
#include "util/scoped_ptr_vector.h"
2522
#include "util/mpbqi.h"
2623
#include "util/timeit.h"
27-
#include "math/polynomial/algebraic_params.hpp"
2824
#include "util/common_msgs.h"
25+
#include "math/polynomial/algebraic_numbers.h"
26+
#include "math/polynomial/upolynomial.h"
27+
#include "math/polynomial/sexpr2upolynomial.h"
28+
#include "math/polynomial/algebraic_params.hpp"
2929

3030
namespace algebraic_numbers {
3131

@@ -124,6 +124,10 @@ namespace algebraic_numbers {
124124
~imp() {
125125
}
126126

127+
bool acell_inv(algebraic_cell const& c) {
128+
return c.m_sign_lower == (upm().eval_sign_at(c.m_p_sz, c.m_p, lower(&c)) == polynomial::sign_neg);
129+
}
130+
127131
void checkpoint() {
128132
if (!m_limit.inc())
129133
throw algebraic_exception(Z3_CANCELED_MSG);
@@ -366,24 +370,25 @@ namespace algebraic_numbers {
366370
return c;
367371
}
368372

369-
int sign_lower(algebraic_cell * c) {
370-
return c->m_sign_lower == 0 ? 1 : -1;
373+
polynomial::sign sign_lower(algebraic_cell * c) const {
374+
return c->m_sign_lower == 0 ? polynomial::sign_pos : polynomial::sign_neg;
371375
}
372376

373-
mpbq & lower(algebraic_cell * c) {
374-
return c->m_interval.lower();
375-
}
377+
mpbq const & lower(algebraic_cell const * c) const { return c->m_interval.lower(); }
376378

377-
mpbq & upper(algebraic_cell * c) {
378-
return c->m_interval.upper();
379-
}
379+
mpbq const & upper(algebraic_cell const * c) const { return c->m_interval.upper(); }
380+
381+
mpbq & lower(algebraic_cell * c) { return c->m_interval.lower(); }
382+
383+
mpbq & upper(algebraic_cell * c) { return c->m_interval.upper(); }
380384

381385
void update_sign_lower(algebraic_cell * c) {
382386
polynomial::sign sl = upm().eval_sign_at(c->m_p_sz, c->m_p, lower(c));
383387
// The isolating intervals are refinable. Thus, the polynomial has opposite signs at lower and upper.
384388
SASSERT(sl != polynomial::sign_zero);
385389
SASSERT(upm().eval_sign_at(c->m_p_sz, c->m_p, upper(c)) == -sl);
386390
c->m_sign_lower = sl == polynomial::sign_neg;
391+
SASSERT(acell_inv(*c));
387392
}
388393

389394
// Make sure the GCD of the coefficients is one and the leading coefficient is positive
@@ -393,6 +398,7 @@ namespace algebraic_numbers {
393398
if (upm().m().is_neg(c->m_p[c->m_p_sz-1])) {
394399
upm().neg(c->m_p_sz, c->m_p);
395400
c->m_sign_lower = !(c->m_sign_lower);
401+
SASSERT(acell_inv(*c));
396402
}
397403
}
398404

@@ -469,6 +475,8 @@ namespace algebraic_numbers {
469475
target->m_sign_lower = source->m_sign_lower;
470476
target->m_not_rational = source->m_not_rational;
471477
target->m_i = source->m_i;
478+
SASSERT(acell_inv(*source));
479+
SASSERT(acell_inv(*target));
472480
}
473481

474482
void set(numeral & a, unsigned sz, mpz const * p, mpbq const & lower, mpbq const & upper, bool minimal) {
@@ -499,7 +507,7 @@ namespace algebraic_numbers {
499507
update_sign_lower(c);
500508
normalize_coeffs(c);
501509
}
502-
SASSERT(sign_lower(a.to_algebraic()) == upm().eval_sign_at(a.to_algebraic()->m_p_sz, a.to_algebraic()->m_p, a.to_algebraic()->m_interval.lower()));
510+
SASSERT(acell_inv(*a.to_algebraic()));
503511
}
504512
TRACE("algebraic", tout << "a: "; display_root(tout, a); tout << "\n";);
505513
}
@@ -519,6 +527,7 @@ namespace algebraic_numbers {
519527
algebraic_cell * c = new (mem) algebraic_cell();
520528
a.m_cell = TAG(void *, c, ROOT);
521529
copy(c, b.to_algebraic());
530+
SASSERT(acell_inv(*c));
522531
}
523532
}
524533
else {
@@ -532,6 +541,7 @@ namespace algebraic_numbers {
532541
del_poly(a.to_algebraic());
533542
del_interval(a.to_algebraic());
534543
copy(a.to_algebraic(), b.to_algebraic());
544+
SASSERT(acell_inv(*a.to_algebraic()));
535545
}
536546
}
537547
}
@@ -693,6 +703,7 @@ namespace algebraic_numbers {
693703
algebraic_cell * c = a.to_algebraic();
694704
if (!upm().normalize_interval_core(c->m_p_sz, c->m_p, sign_lower(c), bqm(), lower(c), upper(c)))
695705
reset(a);
706+
SASSERT(acell_inv(*c));
696707
}
697708
}
698709

@@ -727,7 +738,9 @@ namespace algebraic_numbers {
727738
Return FALSE, if actual root was found.
728739
*/
729740
bool refine_core(algebraic_cell * c) {
730-
return upm().refine_core(c->m_p_sz, c->m_p, sign_lower(c), bqm(), lower(c), upper(c));
741+
bool r = upm().refine_core(c->m_p_sz, c->m_p, sign_lower(c), bqm(), lower(c), upper(c));
742+
SASSERT(acell_inv(*c));
743+
return r;
731744
}
732745

733746
/**
@@ -746,15 +759,17 @@ namespace algebraic_numbers {
746759
if (a.is_basic())
747760
return false;
748761
algebraic_cell * c = a.to_algebraic();
749-
if (!refine_core(c)) {
762+
if (refine_core(c)) {
763+
return true;
764+
}
765+
else {
750766
// root was found
751767
scoped_mpq r(qm());
752768
to_mpq(qm(), lower(c), r);
753769
del(c);
754770
a.m_cell = mk_basic_cell(r);
755771
return false;
756772
}
757-
return true;
758773
}
759774

760775
bool refine(numeral & a, unsigned k) {
@@ -776,6 +791,7 @@ namespace algebraic_numbers {
776791
a.m_cell = mk_basic_cell(r);
777792
return false;
778793
}
794+
SASSERT(acell_inv(*c));
779795
return true;
780796
}
781797

@@ -1573,6 +1589,7 @@ namespace algebraic_numbers {
15731589
upm().p_minus_x(c->m_p_sz, c->m_p);
15741590
bqim().neg(c->m_interval);
15751591
update_sign_lower(c);
1592+
SASSERT(acell_inv(*c));
15761593
}
15771594
}
15781595

@@ -1583,38 +1600,41 @@ namespace algebraic_numbers {
15831600
if (a.is_basic())
15841601
return;
15851602
algebraic_cell * cell_a = a.to_algebraic();
1586-
mpbq & lower = cell_a->m_interval.lower();
1587-
mpbq & upper = cell_a->m_interval.upper();
1588-
if (!bqm().is_zero(lower) && !bqm().is_zero(upper))
1603+
SASSERT(acell_inv(*cell_a));
1604+
mpbq & _lower = cell_a->m_interval.lower();
1605+
mpbq & _upper = cell_a->m_interval.upper();
1606+
if (!bqm().is_zero(_lower) && !bqm().is_zero(_upper))
15891607
return;
1590-
int sign_l = sign_lower(cell_a);
1591-
SASSERT(sign_l != 0);
1592-
int sign_u = -sign_l;
1608+
auto sign_l = sign_lower(cell_a);
1609+
SASSERT(!polynomial::is_zero(sign_l));
1610+
auto sign_u = -sign_l;
15931611

1594-
#define REFINE_LOOP(BOUND, TARGET_SIGN) \
1595-
while (true) { \
1596-
bqm().div2(BOUND); \
1612+
#define REFINE_LOOP(BOUND, TARGET_SIGN) \
1613+
while (true) { \
1614+
bqm().div2(BOUND); \
15971615
polynomial::sign new_sign = upm().eval_sign_at(cell_a->m_p_sz, cell_a->m_p, BOUND); \
1598-
if (new_sign == polynomial::sign_zero) { \
1599-
/* found actual root */ \
1600-
scoped_mpq r(qm()); \
1601-
to_mpq(qm(), BOUND, r); \
1602-
set(a, r); \
1603-
return; \
1604-
} \
1605-
if (new_sign == TARGET_SIGN) \
1606-
return; \
1607-
}
1608-
1609-
if (bqm().is_zero(lower)) {
1610-
bqm().set(lower, upper);
1611-
REFINE_LOOP(lower, sign_l);
1616+
if (new_sign == polynomial::sign_zero) { \
1617+
/* found actual root */ \
1618+
scoped_mpq r(qm()); \
1619+
to_mpq(qm(), BOUND, r); \
1620+
set(a, r); \
1621+
break; \
1622+
} \
1623+
if (new_sign == TARGET_SIGN) { \
1624+
break; \
1625+
} \
1626+
}
1627+
1628+
if (bqm().is_zero(_lower)) {
1629+
bqm().set(_lower, _upper);
1630+
REFINE_LOOP(_lower, sign_l);
16121631
}
16131632
else {
1614-
SASSERT(bqm().is_zero(upper));
1615-
bqm().set(upper, lower);
1616-
REFINE_LOOP(upper, sign_u);
1633+
SASSERT(bqm().is_zero(_upper));
1634+
bqm().set(_upper, _lower);
1635+
REFINE_LOOP(_upper, sign_u);
16171636
}
1637+
SASSERT(acell_inv(*cell_a));
16181638
}
16191639

16201640
void inv(numeral & a) {
@@ -1642,6 +1662,8 @@ namespace algebraic_numbers {
16421662
// convert isolating interval back as a binary rational bound
16431663
upm().convert_q2bq_interval(cell_a->m_p_sz, cell_a->m_p, inv_lower, inv_upper, bqm(), lower(cell_a), upper(cell_a));
16441664
TRACE("algebraic_bug", tout << "after inv: "; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
1665+
update_sign_lower(cell_a);
1666+
SASSERT(acell_inv(*cell_a));
16451667
}
16461668
}
16471669

@@ -2118,7 +2140,6 @@ namespace algebraic_numbers {
21182140
// compute the resultants
21192141
polynomial_ref q_i(pm());
21202142
std::stable_sort(xs.begin(), xs.end(), var_degree_lt(*this, x2v));
2121-
// std::cout << "R: " << R << "\n";
21222143
for (unsigned i = 0; i < xs.size(); i++) {
21232144
checkpoint();
21242145
polynomial::var x_i = xs[i];
@@ -2127,7 +2148,6 @@ namespace algebraic_numbers {
21272148
SASSERT(!v_i.is_basic());
21282149
algebraic_cell * c = v_i.to_algebraic();
21292150
q_i = pm().to_polynomial(c->m_p_sz, c->m_p, x_i);
2130-
// std::cout << "q_i: " << q_i << std::endl;
21312151
pm().resultant(R, q_i, x_i, R);
21322152
SASSERT(!pm().is_zero(R));
21332153
}
@@ -2136,7 +2156,6 @@ namespace algebraic_numbers {
21362156
upm().to_numeral_vector(R, _R);
21372157
unsigned k = upm().nonzero_root_lower_bound(_R.size(), _R.c_ptr());
21382158
TRACE("anum_eval_sign", tout << "R: " << R << "\nk: " << k << "\nri: "<< ri << "\n";);
2139-
// std::cout << "R: " << R << "\n";
21402159
scoped_mpbq mL(bqm()), L(bqm());
21412160
bqm().set(mL, -1);
21422161
bqm().set(L, 1);

‎src/math/polynomial/polynomial.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,7 @@ namespace polynomial {
4848
inline sign operator-(sign s) { switch (s) { case sign_neg: return sign_pos; case sign_pos: return sign_neg; default: return sign_zero; } };
4949
inline sign to_sign(int s) { return s == 0 ? sign_zero : (s > 0 ? sign_pos : sign_neg); }
5050
inline sign operator*(sign a, sign b) { return to_sign((int)a * (int)b); }
51+
inline bool is_zero(sign s) { return s == sign_zero; }
5152

5253
int lex_compare(monomial const * m1, monomial const * m2);
5354
int lex_compare2(monomial const * m1, monomial const * m2, var min_var);

0 commit comments

Comments
 (0)