From eca78aa9c6f70596ba0e9fb8fc509962409b5dcf Mon Sep 17 00:00:00 2001
From: Leonardo de Moura <leonardo@microsoft.com>
Date: Thu, 10 Jan 2013 08:52:25 -0800
Subject: [PATCH] Fix incorrect assertions and bug
Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>
---
src/math/realclosure/mpz_matrix.cpp | 11 ++++++++++-
src/math/realclosure/mpz_matrix.h | 4 +++-
src/math/realclosure/realclosure.cpp | 11 ++++++++---
3 files changed, 21 insertions(+), 5 deletions(-)
diff --git a/src/math/realclosure/mpz_matrix.cpp b/src/math/realclosure/mpz_matrix.cpp
index fcbad5695..836ac35c2 100644
--- a/src/math/realclosure/mpz_matrix.cpp
+++ b/src/math/realclosure/mpz_matrix.cpp
@@ -349,7 +349,7 @@ void mpz_matrix_manager::permute_rows(mpz_matrix const & A, unsigned const * p,
B.swap(C);
}
-unsigned mpz_matrix_manager::linear_independent_rows(mpz_matrix const & _A, unsigned * r) {
+unsigned mpz_matrix_manager::linear_independent_rows(mpz_matrix const & _A, unsigned * r, mpz_matrix & B) {
unsigned r_sz = 0;
scoped_mpz_matrix A(*this);
scoped_mpz g(nm());
@@ -389,6 +389,15 @@ unsigned mpz_matrix_manager::linear_independent_rows(mpz_matrix const & _A, unsi
k2++;
}
std::sort(r, r + r_sz);
+ // Copy linear independent rows to B
+ mpz_matrix & C = A;
+ mk(r_sz, _A.n, C);
+ for (unsigned i = 0; i < r_sz; i++ ) {
+ for (unsigned j = 0; j < _A.n; j++) {
+ nm().set(C(i, j), _A(r[i], j));
+ }
+ }
+ B.swap(C);
return r_sz;
}
diff --git a/src/math/realclosure/mpz_matrix.h b/src/math/realclosure/mpz_matrix.h
index 5986e3a21..b48585c60 100644
--- a/src/math/realclosure/mpz_matrix.h
+++ b/src/math/realclosure/mpz_matrix.h
@@ -108,8 +108,10 @@ public:
\remark The vector r must have at least A.n() capacity
The numer of linear independent rows is returned.
+
+ Store the new matrix in B.
*/
- unsigned linear_independent_rows(mpz_matrix const & A, unsigned * r);
+ unsigned linear_independent_rows(mpz_matrix const & A, unsigned * r, mpz_matrix & B);
// method for debugging purposes
void display(std::ostream & out, mpz_matrix const & A, unsigned cell_width=4) const;
diff --git a/src/math/realclosure/realclosure.cpp b/src/math/realclosure/realclosure.cpp
index 30fceede3..d970136c0 100644
--- a/src/math/realclosure/realclosure.cpp
+++ b/src/math/realclosure/realclosure.cpp
@@ -1761,8 +1761,9 @@ namespace realclosure {
// Solve
// new_M_s * sc_cardinalities = new_taqrs
VERIFY(mm().solve(new_M_s, sc_cardinalities.c_ptr(), new_taqrs.c_ptr()));
+ TRACE("rcf_sign_det", tout << "solution: "; for (unsigned i = 0; i < sc_cardinalities.size(); i++) { tout << sc_cardinalities[i] << " "; } tout << "\n";);
// The solution must contain only positive values <= num_roots
- DEBUG_CODE(for (unsigned j = 0; j < sc_cardinalities.size(); j++) { SASSERT(0 < sc_cardinalities[j] && sc_cardinalities[j] <= num_roots); });
+ DEBUG_CODE(for (unsigned j = 0; j < sc_cardinalities.size(); j++) { SASSERT(0 <= sc_cardinalities[j] && sc_cardinalities[j] <= num_roots); });
// We should keep q only if it discriminated something.
// That is,
// If !use_q2, then There is an i s.t. sc_cardinalities[2*i] > 0 && sc_cardinalities[2*i] > 0
@@ -1799,9 +1800,9 @@ namespace realclosure {
SASSERT(new_scs.empty());
// Update M_s
mm().filter_cols(new_M_s, cols_to_keep.size(), cols_to_keep.c_ptr(), M_s);
- SASSERT(new_M_s.n() == cols_to_keep.size());
+ SASSERT(M_s.n() == cols_to_keep.size());
new_row_idxs.resize(cols_to_keep.size(), 0);
- unsigned new_num_rows = mm().linear_independent_rows(M_s, new_row_idxs.c_ptr());
+ unsigned new_num_rows = mm().linear_independent_rows(M_s, new_row_idxs.c_ptr(), M_s);
SASSERT(new_num_rows == cols_to_keep.size());
// Update taqrs and prs
prs.reset();
@@ -1927,6 +1928,9 @@ namespace realclosure {
\brief Root isolation for polynomials where 0 is not a root.
*/
void nz_isolate_roots(unsigned n, value * const * p, numeral_vector & roots) {
+ TRACE("rcf_isolate",
+ tout << "nz_isolate_roots\n";
+ display_poly(tout, n, p); tout << "\n";);
SASSERT(n > 0);
SASSERT(!is_zero(p[0]));
SASSERT(!is_zero(p[n-1]));
@@ -1943,6 +1947,7 @@ namespace realclosure {
\brief Root isolation entry point.
*/
void isolate_roots(unsigned n, numeral const * p, numeral_vector & roots) {
+ TRACE("rcf_isolate_bug", tout << "isolate_roots: "; for (unsigned i = 0; i < n; i++) { display(tout, p[i]); tout << " "; } tout << "\n";);
SASSERT(n > 0);
SASSERT(!is_zero(p[n-1]));
if (n == 1) {