From 9c8b428ffbffb984538e3417d26a5a3cbf2dd6a7 Mon Sep 17 00:00:00 2001
From: Leonardo de Moura <leonardo@microsoft.com>
Date: Tue, 8 Jan 2013 17:58:34 -0800
Subject: [PATCH] Add matrix operations needed for implementing non-naive sign
determination
Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>
---
src/math/realclosure/mpz_matrix.cpp | 409 ++++++++++++++++++++++++++++
src/math/realclosure/mpz_matrix.h | 151 ++++++++++
src/test/rcf.cpp | 86 +++++-
3 files changed, 644 insertions(+), 2 deletions(-)
create mode 100644 src/math/realclosure/mpz_matrix.cpp
create mode 100644 src/math/realclosure/mpz_matrix.h
diff --git a/src/math/realclosure/mpz_matrix.cpp b/src/math/realclosure/mpz_matrix.cpp
new file mode 100644
index 000000000..fb9327a9e
--- /dev/null
+++ b/src/math/realclosure/mpz_matrix.cpp
@@ -0,0 +1,409 @@
+/*++
+Copyright (c) 2013 Microsoft Corporation
+
+Module Name:
+
+ mpz_matrix.h
+
+Abstract:
+
+ Matrix with integer coefficients. This is not a general purpose
+ module for handling matrices with integer coefficients. Instead,
+ it is a custom package that only contains operations needed to
+ implement Sign Determination (Algorithm 10.11) in the Book:
+ "Algorithms in real algebraic geometry", Basu, Pollack, Roy
+
+ Design choices:
+ - Dense representation. The matrices in Alg 10.11 are small and dense.
+ - Integer coefficients instead of rational coefficients (it only complicates the solver a little bit).
+ Remark: in Algorithm 10.11, the coefficients of the input matrices are always in {-1, 0, 1}.
+ During solving, bigger coefficients are produced, but they are usually very small. It may be
+ an overkill to use mpz instead of int. We use mpz just to be safe.
+ Remark: We do not use rational arithmetic. The solver is slightly more complicated with integers, but is saves space.
+
+Author:
+
+ Leonardo (leonardo) 2013-01-07
+
+Notes:
+
+--*/
+#include"mpz_matrix.h"
+#include"buffer.h"
+
+mpz_matrix_manager::mpz_matrix_manager(unsynch_mpz_manager & nm, small_object_allocator & a):
+ m_nm(nm),
+ m_allocator(a) {
+}
+
+mpz_matrix_manager::~mpz_matrix_manager() {
+}
+
+void mpz_matrix_manager::mk(unsigned m, unsigned n, mpz_matrix & A) {
+ SASSERT(m > 0 && n > 0);
+ del(A);
+ A.m = m;
+ A.n = n;
+ A.a_ij = new (m_allocator) mpz[m*n];
+}
+
+void mpz_matrix_manager::del(mpz_matrix & A) {
+ if (A.a_ij != 0) {
+ for (unsigned i = 0; i < A.m; i++)
+ for (unsigned j = 0; j < A.n; j++)
+ nm().del(A(i,j));
+ unsigned sz = sizeof(mpz) * A.m * A.n;
+ m_allocator.deallocate(sz, A.a_ij);
+ A.m = 0;
+ A.n = 0;
+ A.a_ij = 0;
+ }
+}
+
+void mpz_matrix_manager::set(mpz_matrix & A, mpz_matrix const & B) {
+ if (&A == &B)
+ return;
+ if (A.m != B.m || A.n != B.n) {
+ del(A);
+ mk(B.m, B.n, A);
+ }
+ SASSERT(A.m == B.m && A.n == B.n);
+ for (unsigned i = 0; i < B.m; i++)
+ for (unsigned j = 0; j < B.n; j++)
+ nm().set(A(i, j), B(i, j));
+}
+
+void mpz_matrix_manager::tensor_product(mpz_matrix const & A, mpz_matrix const & B, mpz_matrix & C) {
+ scoped_mpz_matrix _C(*this);
+ mk(A.m * B.m, A.n * B.n, _C);
+ for (unsigned i = 0; i < _C.m(); i++)
+ for (unsigned j = 0; j < _C.n(); j++)
+ nm().mul(A(i / B.m, j / B.n),
+ B(i % B.m, j % B.n),
+ _C(i, j));
+ C.swap(_C);
+}
+
+void mpz_matrix_manager::swap_rows(mpz_matrix & A, unsigned i, unsigned j) {
+ if (i != j) {
+ for (unsigned k = 0; k < A.n; k++)
+ ::swap(A(i, k), A(j, k));
+ }
+}
+
+// If b_i == 0, then method just divides the given row by its GCD
+// If b_i != 0
+// If the GCD of the row divides *b_i
+// divide the row and *b_i by the GCD
+// Else
+// If int_solver == true ==> return false (the system is unsolvable)
+bool mpz_matrix_manager::normalize_row(mpz * A_i, unsigned n, mpz * b_i, bool int_solver) {
+ scoped_mpz g(nm());
+ bool first = true;
+ for (unsigned j = 0; j < n; j++) {
+ if (nm().is_zero(A_i[j]))
+ continue;
+ if (first) {
+ nm().set(g, A_i[j]);
+ nm().abs(g);
+ first = false;
+ }
+ else {
+ nm().gcd(g, A_i[j], g);
+ }
+ if (nm().is_one(g))
+ return true;
+ }
+ if (first)
+ return true; // zero row
+ if (!nm().is_one(g)) {
+ if (b_i) {
+ if (nm().divides(g, *b_i)) {
+ for (unsigned j = 0; j < n; j++) {
+ nm().div(A_i[j], g, A_i[j]);
+ }
+ nm().div(*b_i, g, *b_i);
+ }
+ else {
+ if (int_solver)
+ return false; // system does not have an integer solution
+ }
+ }
+ else {
+ for (unsigned j = 0; j < n; j++) {
+ nm().div(A_i[j], g, A_i[j]);
+ }
+ }
+ }
+ return true;
+}
+
+/*
+ Given a matrix of the form
+
+ k2
+ |
+ V
+ X X ... X X ... X
+ 0 X ... X X ... X
+ ... ... X X ... X
+k1=> 0 0 ... 0 X ... X
+ 0 0 ... 0 X ... X
+ ... ... 0 X ... X
+ 0 0 ... 0 X ... X
+
+ It will "zero" the elements a_{k1+1, k2} ... a_{m, k2} by addining multiples of the row k1 to multiples of the
+ rows k1+1, ..., m
+
+ The resultant matrix will look like
+
+ k2
+ |
+ V
+ X X ... X X ... X
+ 0 X ... X X ... X
+ ... ... X X ... X
+k1=> 0 0 ... 0 X ... X
+ 0 0 ... 0 0 ... X
+ ... ... 0 0 ... X
+ 0 0 ... 0 0 ... X
+
+
+ If b != 0, then the transformations are also applied to b.
+ If int_solver == true and b != 0, then the method returns false if when
+ performing the transformations it detected that it is impossible to
+ solve the integer system of equations A x = b.
+*/
+bool mpz_matrix_manager::eliminate(mpz_matrix & A, mpz * b, unsigned k1, unsigned k2, bool int_solver) {
+ // check if first k2-1 positions of row k1 are 0
+ DEBUG_CODE(for (unsigned j = 0; j < k2; j++) { SASSERT(nm().is_zero(A(k1, j))); });
+ mpz & a_kk = A(k1, k2);
+ SASSERT(!nm().is_zero(a_kk));
+ scoped_mpz t1(nm()), t2(nm());
+ scoped_mpz a_ik_prime(nm()), a_kk_prime(nm()), lcm_a_kk_a_ik(nm());
+ // for all rows below pivot
+ for (unsigned i = k1+1; i < A.m; i++) {
+ // check if first k-1 positions of row k are 0
+ DEBUG_CODE(for (unsigned j = 0; j < k2; j++) { SASSERT(nm().is_zero(A(i, j))); });
+ mpz & a_ik = A(i, k2);
+ if (!nm().is_zero(a_ik)) {
+ // a_ik' = lcm(a_kk, a_ik)/a_kk
+ // a_kk' = lcm(a_kk, a_ik)/a_ik
+ nm().lcm(a_kk, a_ik, lcm_a_kk_a_ik);
+ nm().div(lcm_a_kk_a_ik, a_kk, a_ik_prime);
+ nm().div(lcm_a_kk_a_ik, a_ik, a_kk_prime);
+ for (unsigned j = k2+1; j < A.n; j++) {
+ // a_ij <- a_kk' * a_ij - a_ik' * a_kj
+ nm().mul(a_ik_prime, A(k1, j), t1);
+ nm().mul(a_kk_prime, A(i, j), t2);
+ nm().sub(t2, t1, A(i, j));
+ }
+ if (b) {
+ // b_i <- a_kk' * b_i - a_ik' * b_k
+ nm().mul(a_ik_prime, b[k1], t1);
+ nm().mul(a_kk_prime, b[i], t2);
+ nm().sub(t2, t1, b[i]);
+ }
+ // a_ik <- 0
+ nm().set(A(i, k2), 0);
+ // normalize row i
+ if (!normalize_row(A.row(i), A.n, b ? &(b[i]) : 0, int_solver))
+ return false;
+ }
+ SASSERT(nm().is_zero(A(i, k2)));
+ }
+ return true;
+}
+
+bool mpz_matrix_manager::solve_core(mpz_matrix const & _A, mpz * b, bool int_solver) {
+ SASSERT(_A.n == _A.m);
+ scoped_mpz_matrix A(*this);
+ set(A, _A);
+ for (unsigned k = 0; k < A.m(); k++) {
+ TRACE("mpz_matrix",
+ tout << "k: " << k << "\n" << A;
+ tout << "b:";
+ for (unsigned i = 0; i < A.m(); i++) {
+ tout << " ";
+ nm().display(tout, b[i]);
+ }
+ tout << "\n";);
+ // find pivot
+ unsigned i = k;
+ for (; i < A.m(); i++) {
+ if (!nm().is_zero(A(i, k)))
+ break;
+ }
+ if (i == A.m())
+ return false; // matrix is singular
+ // swap rows k and i
+ swap_rows(A, k, i);
+ swap(b[k], b[i]);
+ //
+ if (!eliminate(A, b, k, k, int_solver))
+ return false;
+ }
+ // Back substitution
+ unsigned k = A.m();
+ while (k > 0) {
+ --k;
+ DEBUG_CODE(for (unsigned j = 0; j < A.n(); j++) { SASSERT(j == k || nm().is_zero(A(k, j))); });
+ SASSERT(!nm().is_zero(A(k, k)));
+ if (nm().divides(A(k, k), b[k])) {
+ nm().div(b[k], A(k, k), b[k]);
+ nm().set(A(k, k), 1);
+ }
+ else {
+ if (int_solver)
+ return false; // no integer solution
+ if (nm().is_neg(A(k, k))) {
+ nm().neg(A(k, k));
+ nm().neg(b[k]);
+ }
+ }
+ if (!int_solver) {
+ // REMARK:
+ // For the sign determination algorithm, we only use int_solver == true.
+ //
+ // TODO: implement backward substitution when int_solver == false
+ // In this case, A(k, k) may not be 1.
+ NOT_IMPLEMENTED_YET();
+ }
+ SASSERT(!int_solver || nm().is_one(A(k, k)));
+ // back substitute
+ unsigned i = k;
+ while (i > 0) {
+ --i;
+ // Assuming int_solver == true
+ SASSERT(int_solver); // See comment above
+ // b_i <- b_i - a_ik * b_k
+ nm().submul(b[i], A(i, k), b[k], b[i]);
+ nm().set(A(i, k), 0);
+ }
+ }
+ return true;
+}
+
+bool mpz_matrix_manager::solve(mpz_matrix const & A, mpz * b, mpz const * c) {
+ for (unsigned i = 0; i < A.n; i++)
+ nm().set(b[i], c[i]);
+ return solve_core(A, b, true);
+}
+
+bool mpz_matrix_manager::solve(mpz_matrix const & A, int * b, int const * c) {
+ scoped_mpz_matrix _b(*this);
+ mk(A.n, 1, _b);
+ for (unsigned i = 0; i < A.n; i++)
+ nm().set(_b(i,0), c[i]);
+ bool r = solve_core(A, _b.A.a_ij, true);
+ if (r) {
+ for (unsigned i = 0; i < A.n; i++)
+ b[i] = _b.get_int(i, 0);
+ }
+ return r;
+}
+
+void mpz_matrix_manager::filter_cols(mpz_matrix const & A, unsigned num_cols, unsigned const * cols, mpz_matrix & B) {
+ SASSERT(num_cols <= A.n);
+ // Check pre-condition:
+ // - All elements in cols are smaller than A.n
+ // - cols is sorted
+ // - cols does not contain repeated elements
+ DEBUG_CODE({
+ for (unsigned i = 0; i < num_cols; i ++) {
+ SASSERT(cols[i] < A.n);
+ SASSERT(i == 0 || cols[i-1] < cols[i]);
+ }
+ });
+ if (num_cols == A.n) {
+ // keep everything
+ set(B, A);
+ }
+ else {
+ SASSERT(num_cols < A.n);
+ scoped_mpz_matrix C(*this);
+ mk(A.m, num_cols, C);
+ for (unsigned i = 0; i < A.m; i++)
+ for (unsigned j = 0; j < num_cols; j++)
+ nm().set(C(i, j), A(i, cols[j]));
+ B.swap(C);
+ }
+}
+
+void mpz_matrix_manager::permute_rows(mpz_matrix const & A, unsigned const * p, mpz_matrix & B) {
+ // Check if p is really a permutation
+ DEBUG_CODE({
+ buffer<bool> seen;
+ seen.resize(A.m, false);
+ for (unsigned i = 0; i < A.m; i++) {
+ SASSERT(p[i] < A.m);
+ SASSERT(!seen[p[i]]);
+ seen[p[i]] = true;
+ }
+ });
+ scoped_mpz_matrix C(*this);
+ mk(A.m, A.n, C);
+ for (unsigned i = 0; i < A.m; i++)
+ for (unsigned j = 0; j < A.n; j++)
+ nm().set(C(i, j), A(p[i], j));
+ B.swap(C);
+}
+
+void mpz_matrix_manager::linear_independent_rows(mpz_matrix const & _A, unsigned_vector & r) {
+ r.reset();
+ scoped_mpz_matrix A(*this);
+ scoped_mpz g(nm());
+ scoped_mpz t1(nm()), t2(nm());
+ scoped_mpz a_ik_prime(nm()), a_kk_prime(nm()), lcm_a_kk_a_ik(nm());
+ set(A, _A);
+ sbuffer<unsigned, 128> rows;
+ rows.resize(A.m(), 0);
+ for (unsigned i = 0; i < A.m(); i++)
+ rows[i] = i;
+ for (unsigned k1 = 0, k2 = 0; k1 < A.m(); k1++) {
+ TRACE("mpz_matrix", tout << "k1: " << k1 << ", k2: " << k2 << "\n" << A;);
+ // find pivot
+ unsigned pivot = UINT_MAX;
+ for (unsigned i = k1; i < A.m(); i++) {
+ if (!nm().is_zero(A(i, k2))) {
+ if (pivot == UINT_MAX) {
+ pivot = i;
+ }
+ else {
+ if (rows[i] < rows[pivot])
+ pivot = i;
+ }
+ }
+ }
+ if (pivot == UINT_MAX)
+ continue;
+ // swap rows k and pivot
+ swap_rows(A, k1, pivot);
+ std::swap(rows[k1], rows[pivot]);
+ //
+ r.push_back(rows[k1]);
+ if (r.size() >= A.n())
+ break;
+ eliminate(A, 0, k1, k2, false);
+ k2++;
+ }
+ std::sort(r.begin(), r.end());
+}
+
+void mpz_matrix_manager::display(std::ostream & out, mpz_matrix const & A, unsigned cell_width) const {
+ out << A.m << " x " << A.n << " Matrix\n";
+ for (unsigned i = 0; i < A.m; i++) {
+ for (unsigned j = 0; j < A.n; j++) {
+ if (j > 0)
+ out << " ";
+ std::string s = nm().to_string(A(i, j));
+ if (s.size() < cell_width) {
+ unsigned space = cell_width - s.size();
+ for (unsigned k = 0; k < space; k++)
+ out << " ";
+ }
+ out << s;
+ }
+ out << "\n";
+ }
+}
diff --git a/src/math/realclosure/mpz_matrix.h b/src/math/realclosure/mpz_matrix.h
new file mode 100644
index 000000000..9ce9aaac7
--- /dev/null
+++ b/src/math/realclosure/mpz_matrix.h
@@ -0,0 +1,151 @@
+/*++
+Copyright (c) 2013 Microsoft Corporation
+
+Module Name:
+
+ mpz_matrix.h
+
+Abstract:
+
+ Matrix with integer coefficients. This is not a general purpose
+ module for handling matrices with integer coefficients. Instead,
+ it is a custom package that only contains operations needed to
+ implement Sign Determination (Algorithm 10.11) in the Book:
+ "Algorithms in real algebraic geometry", Basu, Pollack, Roy
+
+ Design choices:
+ - Dense representation. The matrices in Alg 10.11 are small and dense.
+ - Integer coefficients instead of rational coefficients (it only complicates the solver a little bit).
+ Remark: in Algorithm 10.11, the coefficients of the input matrices are always in {-1, 0, 1}.
+ During solving, bigger coefficients are produced, but they are usually very small. It may be
+ an overkill to use mpz instead of int. We use mpz just to be safe.
+ Remark: We do not use rational arithmetic. The solver is slightly more complicated with integers, but is saves space.
+
+Author:
+
+ Leonardo (leonardo) 2013-01-07
+
+Notes:
+
+--*/
+#ifndef _MPZ_MATRIX_H_
+#define _MPZ_MATRIX_H_
+
+#include"mpz.h"
+
+/**
+ \brief A mxn matrix.
+ Remark: Algorithm 10.11 only uses square matrices, but supporting
+ arbitrary matrices does not increase the complexity of this module.
+*/
+class mpz_matrix {
+ friend class mpz_matrix_manager;
+ friend class scoped_mpz_matrix;
+ unsigned m;
+ unsigned n;
+ mpz * a_ij;
+public:
+ mpz_matrix():m(0), n(0), a_ij(0) {}
+ mpz const & operator()(unsigned i, unsigned j) const {
+ SASSERT(i < m);
+ SASSERT(j < n);
+ return a_ij[i*n + j]; }
+ mpz & operator()(unsigned i, unsigned j) { SASSERT(i < m); SASSERT(j < n); return a_ij[i*n + j]; }
+ void swap(mpz_matrix & B) { std::swap(m, B.m); std::swap(n, B.n); std::swap(a_ij, B.a_ij); }
+ mpz * row(unsigned i) const { SASSERT(i < m); return a_ij + i*n; }
+};
+
+class mpz_matrix_manager {
+ unsynch_mpz_manager & m_nm;
+ small_object_allocator & m_allocator;
+ static void swap_rows(mpz_matrix & A, unsigned i, unsigned j);
+ bool normalize_row(mpz * A_i, unsigned n, mpz * b_i, bool int_solver);
+ bool eliminate(mpz_matrix & A, mpz * b, unsigned k1, unsigned k2, bool int_solver);
+ bool solve_core(mpz_matrix const & A, mpz * b, bool int_solver);
+public:
+ mpz_matrix_manager(unsynch_mpz_manager & nm, small_object_allocator & a);
+ ~mpz_matrix_manager();
+ unsynch_mpz_manager & nm() const { return m_nm; }
+ void mk(unsigned m, unsigned n, mpz_matrix & A);
+ void del(mpz_matrix & r);
+ void set(mpz_matrix & A, mpz_matrix const & B);
+ void tensor_product(mpz_matrix const & A, mpz_matrix const & B, mpz_matrix & C);
+ /**
+ \brief Solve A*b = c
+
+ Return false if the system does not have an integer solution.
+
+ \pre A is a square matrix
+ \pre b and c are vectors of size A.n (== A.m)
+ */
+ bool solve(mpz_matrix const & A, mpz * b, mpz const * c);
+ /**
+ \brief Solve A*b = c
+
+ Return false if the system does not have an integer solution.
+
+ \pre A is a square matrix
+ \pre b and c are vectors of size A.n (== A.m)
+ */
+ bool solve(mpz_matrix const & A, int * b, int const * c);
+ /**
+ \brief Store in B that contains the subset cols of columns of A.
+
+ \pre num_cols <= A.n
+ \pre Forall i < num_cols, cols[i] < A.n
+ \pre Forall 0 < i < num_cols, cols[i-1] < cols[i]
+ */
+ void filter_cols(mpz_matrix const & A, unsigned num_cols, unsigned const * cols, mpz_matrix & B);
+ /**
+ \brief Store in B the matrix obtained after applying the given permutation to the rows of A.
+ */
+ void permute_rows(mpz_matrix const & A, unsigned const * p, mpz_matrix & B);
+ /**
+ \brief Store in r the row (ids) of A that are linear independent.
+
+ \remark If there is an option between rows i and j,
+ this method will give preference to the row that occurs first.
+ */
+ void linear_independent_rows(mpz_matrix const & A, unsigned_vector & r);
+
+ // method for debugging purposes
+ void display(std::ostream & out, mpz_matrix const & A, unsigned cell_width=4) const;
+};
+
+class scoped_mpz_matrix {
+ friend class mpz_matrix_manager;
+ mpz_matrix_manager & m_manager;
+ mpz_matrix A;
+public:
+ scoped_mpz_matrix(mpz_matrix_manager & m):m_manager(m) {}
+ scoped_mpz_matrix(mpz_matrix const & B, mpz_matrix_manager & m):m_manager(m) { m_manager.set(A, B); }
+ ~scoped_mpz_matrix() { m_manager.del(A); }
+ mpz_matrix_manager & mm() const { return m_manager; }
+ unsynch_mpz_manager & nm() const { return mm().nm(); }
+
+ unsigned m() const { return A.m; }
+ unsigned n() const { return A.n; }
+ mpz * row(unsigned i) const { return A.row(i); }
+
+ operator mpz_matrix const &() const { return A; }
+ operator mpz_matrix &() { return A; }
+ mpz_matrix const & get() const { return A; }
+ mpz_matrix & get() { return A; }
+
+ void swap(mpz_matrix & B) { A.swap(B); }
+
+ void set(unsigned i, unsigned j, mpz const & v) { nm().set(A(i, j), v); }
+ void set(unsigned i, unsigned j, int v) { nm().set(A(i, j), v); }
+
+ mpz const & operator()(unsigned i, unsigned j) const { return A(i, j); }
+ mpz & operator()(unsigned i, unsigned j) { return A(i, j); }
+
+ int get_int(unsigned i, unsigned j) const { SASSERT(nm().is_int(A(i, j))); return nm().get_int(A(i, j)); }
+};
+
+inline std::ostream & operator<<(std::ostream & out, scoped_mpz_matrix const & m) {
+ m.mm().display(out, m);
+ return out;
+}
+
+#endif
diff --git a/src/test/rcf.cpp b/src/test/rcf.cpp
index b50ebc56c..0583c54a0 100644
--- a/src/test/rcf.cpp
+++ b/src/test/rcf.cpp
@@ -17,6 +17,7 @@ Notes:
--*/
#include"realclosure.h"
+#include"mpz_matrix.h"
static void tst1() {
unsynch_mpq_manager qm;
@@ -66,6 +67,87 @@ static void tst1() {
std::cout << interval_pp((a + eps)/(a - eps)) << std::endl;
}
-void tst_rcf() {
- tst1();
+static void tst2() {
+ enable_trace("mpz_matrix");
+ unsynch_mpq_manager nm;
+ small_object_allocator allocator;
+ mpz_matrix_manager mm(nm, allocator);
+ scoped_mpz_matrix A(mm);
+ mm.mk(3, 3, A);
+ // Matrix
+ // 1 1 1
+ // 0 1 -1
+ // 0 1 1
+ A.set(0, 0, 1); A.set(0, 1, 1); A.set(0, 2, 1);
+ A.set(1, 0, 0); A.set(1, 1, 1); A.set(1, 2, -1);
+ A.set(2, 0, 0); A.set(2, 1, 1); A.set(2, 2, 1);
+ std::cout << A;
+ {
+ int b[3];
+ int c[3] = { 10, -2, 8 };
+ std::cout << "solve: " << mm.solve(A, b, c) << "\n";
+ for (unsigned i = 0; i < 3; i++) std::cout << b[i] << " ";
+ std::cout << "\n";
+ }
+ scoped_mpz_matrix A2(mm);
+ mm.tensor_product(A, A, A2);
+ std::cout << A2;
+ scoped_mpz_matrix B(mm);
+ unsigned cols[] = { 1, 3, 7, 8 };
+ mm.filter_cols(A2, 4, cols, B);
+ std::cout << B;
+ scoped_mpz_matrix C(mm);
+ unsigned perm[] = { 8, 7, 6, 5, 4, 3, 2, 1, 0 };
+ mm.permute_rows(B, perm, C);
+ std::cout << C;
+}
+
+static void tst_solve(unsigned n, int _A[], int _b[], int _c[], bool solved) {
+ unsynch_mpq_manager nm;
+ small_object_allocator allocator;
+ mpz_matrix_manager mm(nm, allocator);
+ scoped_mpz_matrix A(mm);
+ mm.mk(n, n, A);
+ for (unsigned i = 0; i < n; i++)
+ for (unsigned j = 0; j < n; j++)
+ A.set(i, j, _A[i*n + j]);
+ svector<int> b;
+ b.resize(n, 0);
+ if (mm.solve(A, b.c_ptr(), _c)) {
+ SASSERT(solved);
+ for (unsigned i = 0; i < n; i++) {
+ SASSERT(b[i] == _b[i]);
+ }
+ }
+ else {
+ SASSERT(!solved);
+ }
+}
+
+static void tst_lin_indep(unsigned m, unsigned n, int _A[], unsigned ex_sz, unsigned ex_r[]) {
+ unsynch_mpq_manager nm;
+ small_object_allocator allocator;
+ mpz_matrix_manager mm(nm, allocator);
+ scoped_mpz_matrix A(mm);
+ mm.mk(m, n, A);
+ for (unsigned i = 0; i < m; i++)
+ for (unsigned j = 0; j < n; j++)
+ A.set(i, j, _A[i*n + j]);
+ unsigned_vector r;
+ mm.linear_independent_rows(A, r);
+ SASSERT(r.size() == ex_sz);
+ for (unsigned i = 0; i < ex_sz; i++) {
+ SASSERT(r[i] == ex_r[i]);
+ }
+}
+
+
+void tst_rcf() {
+ // tst1();
+ tst2();
+ { int A[] = {0, 1, 1, 1, 0, 1, 1, 1, -1}; int c[] = {10, 4, -4}; int b[] = {-2, 4, 6}; tst_solve(3, A, b, c, true); }
+ { int A[] = {1, 1, 1, 0, 1, 1, 0, 1, 1}; int c[] = {3, 2, 2}; int b[] = {1, 1, 1}; tst_solve(3, A, b, c, false); }
+ { int A[] = {1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, -1}; unsigned r[] = {0, 1, 4}; tst_lin_indep(5, 3, A, 3, r); }
+ { int A[] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, -1}; unsigned r[] = {0, 4}; tst_lin_indep(5, 3, A, 2, r); }
+ { int A[] = {1, 1, 1, 1, 1, 1, 1, 0, 1, 2, 1, 2, 3, 1, 3}; unsigned r[] = {0, 2}; tst_lin_indep(5, 3, A, 2, r); }
}