From 26a3d2ca31072116a8f5c2cd7c0447f43acdc87e Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Tue, 21 Jan 2014 08:40:28 -0800 Subject: [PATCH] add stand-alone simplex Signed-off-by: Nikolaj Bjorner --- scripts/mk_project.py | 3 +- src/math/simplex/simplex.cpp | 26 ++ src/math/simplex/simplex.h | 157 ++++++++ src/math/simplex/simplex_def.h | 480 ++++++++++++++++++++++ src/math/simplex/sparse_matrix.h | 254 ++++++++++++ src/math/simplex/sparse_matrix_def.h | 580 +++++++++++++++++++++++++++ src/test/simplex.cpp | 24 ++ src/util/mpq_inf.cpp | 2 +- src/util/mpq_inf.h | 54 ++- 9 files changed, 1560 insertions(+), 20 deletions(-) create mode 100644 src/math/simplex/simplex.cpp create mode 100644 src/math/simplex/simplex.h create mode 100644 src/math/simplex/simplex_def.h create mode 100644 src/math/simplex/sparse_matrix.h create mode 100644 src/math/simplex/sparse_matrix_def.h create mode 100644 src/test/simplex.cpp diff --git a/scripts/mk_project.py b/scripts/mk_project.py index 6215e0fa0..952c22c20 100644 --- a/scripts/mk_project.py +++ b/scripts/mk_project.py @@ -15,6 +15,7 @@ def init_project_def(): add_lib('sat', ['util']) add_lib('nlsat', ['polynomial', 'sat']) add_lib('hilbert', ['util'], 'math/hilbert') + add_lib('simplex', ['util'], 'math/simplex') add_lib('interval', ['util'], 'math/interval') add_lib('realclosure', ['interval'], 'math/realclosure') add_lib('subpaving', ['interval'], 'math/subpaving') @@ -75,7 +76,7 @@ def init_project_def(): add_lib('api', ['portfolio', 'user_plugin', 'smtparser', 'realclosure','interp','opt'], includes2install=['z3.h', 'z3_v1.h', 'z3_macros.h'] + API_files) add_exe('shell', ['api', 'sat', 'extra_cmds','opt'], exe_name='z3') - add_exe('test', ['api', 'fuzzing'], exe_name='test-z3', install=False) + add_exe('test', ['api', 'fuzzing', 'simplex'], exe_name='test-z3', install=False) add_dll('api_dll', ['api', 'sat', 'extra_cmds'], 'api/dll', reexports=['api'], dll_name='libz3', diff --git a/src/math/simplex/simplex.cpp b/src/math/simplex/simplex.cpp new file mode 100644 index 000000000..494c0b6bb --- /dev/null +++ b/src/math/simplex/simplex.cpp @@ -0,0 +1,26 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + simplex.h + +Abstract: + + Multi-precision simplex tableau. + +Author: + + Nikolaj Bjorner (nbjorner) 2014-01-15 + +Notes: + +--*/ + +#include"simplex.h" +#include"sparse_matrix_def.h" +#include"simplex_def.h" +namespace simplex { + template class simplex; + template class simplex; +}; diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h new file mode 100644 index 000000000..0cb4852f2 --- /dev/null +++ b/src/math/simplex/simplex.h @@ -0,0 +1,157 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + simplex.h + +Abstract: + + Multi-precision simplex tableau. + + - It uses code from theory_arith where applicable. + + - It is detached from the theory class and ASTs. + + - It uses non-shared mpz/mpq's avoiding global locks and operations on rationals. + + - It follows the same sparse tableau layout (no LU yet). + + - It does not include features for non-linear arithmetic. + + - Branch/bound/cuts is external. + +Author: + + Nikolaj Bjorner (nbjorner) 2014-01-15 + +Notes: + +--*/ + +#ifndef _SIMPLEX_H_ +#define _SIMPLEX_H_ + +#include "sparse_matrix.h" +#include "mpq_inf.h" +#include "heap.h" +#include "lbool.h" + +namespace simplex { + + template + class simplex { + + typedef unsigned var_t; + typedef typename Ext::eps_numeral eps_numeral; + typedef typename Ext::numeral numeral; + typedef typename Ext::manager manager; + typedef typename Ext::eps_manager eps_manager; + typedef typename Ext::scoped_numeral scoped_numeral; + typedef _scoped_numeral scoped_eps_numeral; + + typedef typename _scoped_numeral_vector scoped_eps_numeral_vector; + + typedef sparse_matrix matrix; + + struct var_lt { + bool operator()(var_t v1, var_t v2) const { return v1 < v2; } + }; + typedef heap var_heap; + + enum pivot_strategy_t { + S_BLAND, + S_GREATEST_ERROR, + S_LEAST_ERROR, + S_DEFAULT + }; + + struct var_info { + unsigned m_base2row:29; + unsigned m_is_base:1; + unsigned m_lower_valid:1; + unsigned m_upper_valid:1; + eps_numeral m_value; + eps_numeral m_lower; + eps_numeral m_upper; + numeral m_base_coeff; + var_info(): + m_base2row(0), + m_is_base(false), + m_lower_valid(false), + m_upper_valid(false) + {} + }; + + static const var_t null_var = UINT_MAX; + matrix M; + unsigned m_max_iterations; + volatile bool m_cancel; + var_heap m_to_patch; + mutable manager m; + mutable eps_manager em; + vector m_vars; + svector m_row2base; + bool m_bland; + random_gen m_random; + + public: + simplex(): + m_max_iterations(UINT_MAX), + m_cancel(false), + m_to_patch(1024), + m_bland(false) {} + + typedef typename matrix::row row; + typedef typename matrix::row_iterator row_iterator; + typedef typename matrix::col_iterator col_iterator; + + void ensure_var(var_t v); + row add_row(unsigned num_vars, var_t base, var_t const* vars, numeral const* coeffs); + void del_row(row const& r); + void set_lower(var_t var, eps_numeral const& b); + void set_upper(var_t var, eps_numeral const& b); + void unset_lower(var_t var); + void unset_upper(var_t var); + void set_value(var_t var, eps_numeral const& b); + void set_cancel(bool f) { m_cancel = f; } + void set_max_iterations(unsigned m) { m_max_iterations = m; } + lbool make_feasible(); + lbool optimize(var_t var); + eps_numeral const& get_value(var_t v); + void display(std::ostream& out) const; + + private: + + var_t select_var_to_fix(); + pivot_strategy_t pivot_strategy(); + var_t select_smallest_var() { return m_to_patch.empty()?null_var:m_to_patch.erase_min(); } + var_t select_error_var(bool least); + // row get_infeasible_row() { } + void check_blands_rule(var_t v) { } + bool make_var_feasible(var_t x_i); + void update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, eps_numeral const& new_value); + void update_value(var_t v, eps_numeral const& delta); + void update_value_core(var_t v, eps_numeral const& delta); + var_t select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); + var_t select_blands_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); + template + var_t select_pivot_core(var_t x_i, scoped_numeral& out_a_ij); + int get_num_non_free_dep_vars(var_t x_j, int best_so_far); + + bool below_lower(var_t v) const; + bool above_upper(var_t v) const; + bool above_lower(var_t v) const; + bool below_upper(var_t v) const; + bool outside_bounds(var_t v) const { return below_lower(v) || above_upper(v); } + bool is_free(var_t v) const { return !m_vars[v].m_lower_valid && !m_vars[v].m_upper_valid; } + bool is_non_free(var_t v) const { return !is_free(v); } + unsigned get_num_vars() const { return m_vars.size(); } + bool is_base(var_t x) const { return m_vars[x].m_is_base; } + + bool well_formed() const; + }; + +}; + +#endif diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h new file mode 100644 index 000000000..b92ec7551 --- /dev/null +++ b/src/math/simplex/simplex_def.h @@ -0,0 +1,480 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + simplex_def.h + +Abstract: + +Author: + + Nikolaj Bjorner (nbjorner) 2014-01-15 + +Notes: + +--*/ + +#ifndef _SIMPLEX_DEF_H_ +#define _SIMPLEX_DEF_H_ + + +namespace simplex { + + template + typename simplex::row + simplex::add_row(var_t base, unsigned num_vars, var_t const* vars, numeral const* coeffs) { + DEBUG_CODE( + bool found = false; + for (unsigned i = 0; !found && i < num_vars; ++i) found = vars[i] == base; + SASSERT(found); + ); + row r = M.mk_row(); + for (unsigned i = 0; i < num_vars; ++i) { + M.add(r, coeffs[i], vars[i]); + } + while (m_row2base.size() <= r.id()) { + m_row2base.push_back(null_var); + } + m_row2base[r.id()] = base; + m_vars[base].m_base2row = r.id(); + m_vars[base].m_is_base = true; + return r; + } + + template + void simplex::del_row(row const& r) { + m_is_base[m_row2base[r.id()]] = false; + M.del(r); + } + + template + void simplex::set_lower(var_t var, eps_numeral const& b) { + var_info& vi = m_vars[var]; + em.set(vi.m_lower, b); + vi.m_lower_valid = true; + if (em.lt(vi.m_value, b)) { + scoped_eps_numeral delta(em); + em.sub(b, vi.m_value, delta); + update_value(var, delta); + } + } + + template + void simplex::set_upper(var_t var, eps_numeral const& b) { + var_info& vi = m_vars[var]; + em.set(vi.m_upper, b); + vi.m_upper_valid = true; + if (em.gt(vi.m_value, b)) { + scoped_eps_numeral delta(em); + em.sub(b, vi.m_value, delta); + update_value(var, delta); + } + } + + template + void simplex::unset_lower(var_t var) { + m_vars[var].m_lower_valid = false; + } + + template + void simplex::unset_upper(var_t var) { + m_vars[var].m_upper_valid = false; + } + + template + void simplex::set_value(var_t var, eps_numeral const& b) { + scoped_eps_numeral delta(em); + em.sub(b, m_vars[var].m_value, delta); + update_value(var, delta); + } + + template + typename simplex::eps_numeral const& + simplex::get_value(var_t v) { + return m_vars[v].m_value; + } + + template + void simplex::display(std::ostream& out) const { + M.display(out); + for (unsigned i = 0; i < m_vars.size(); ++i) { + var_info const& vi = m_vars[i]; + out << "v" << i << " "; + if (vi.m_is_base) out << "b:" << vi.m_base2row << " "; + em.display(out, vi.m_value); + out << " ["; + if (vi.m_lower_valid) em.display(out, vi.m_lower) else out << "-oo"; + out << ":"; + if (vi.m_upper_valid) em.display(out, vi.m_upper) else out << "oo"; + out << "]"; + out << "\n"; + } + } + + template + void simplex::ensure_var(var_t v) { + while (v >= m_vars.size()) { + M.ensure_var(m_vars.size()); + m_vars.push_back(var_info()); + } + } + + template + lbool simplex::make_feasible() { + unsigned num_iterations = 0; + var_t v = null_var; + while ((v = select_var_to_fix()) != null_var) { + if (m_cancel || num_iterations > m_max_iterations) { + return l_undef; + } + check_blands_rule(v); + if (!make_var_feasible(v)) { + return l_false; + } + ++num_iterations; + SASSERT(well_formed()); + } + return l_true; + } + + /** + \brief Make x_j the new base variable for row of x_i. + x_i is assumed base variable of row r_i. + x_j is assumed to have coefficient a_ij in r_i. + + a_ii*x_i + a_ij*x_j + r_i = 0 + + current value of x_i is v_i + new value of x_i is new_value + a_ii*(x_i + new_value - x_i) + a_ij*((x_i - new_value)*a_ii/a_ij + x_j) + r_i = 0 + + Let r_k be a row where x_j has coefficient x_kj != 0. + r_k <- r_k * a_ij - r_i * a_kj + */ + template + void simplex::update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, eps_numeral const& new_value) { + SASSERT(is_base(x_i)); + SASSERT(!is_base(x_j)); + var_info& x_iI = m_vars[x_i]; + var_info& x_jI = m_vars[x_j]; + scoped_eps_numeral theta(em); + theta = x_iI.m_value; + theta -= new_value; + numeral const& a_ii = x_iI.m_base_coeff; + em.mul(theta, a_ii, theta); + em.div(theta, a_ij, theta); + update_value(x_j, theta); + SASSERT(em.eq(x_iI.m_value, new_value)); + unsigned r_i = x_iI.m_base2row; + x_jI.m_base2row = r_i; + m_row2base[r_i] = x_j; + x_jI.m_is_base = true; + x_iI.m_is_base = false; + if (outside_bounds(x_j)) { + m_to_patch.insert(x_j); + } + col_iterator it = M.col_begin(x_j), end = M.col_end(x_j); + scoped_numeral a_kj(m), g(m); + for (; it != end; ++it) { + row r_k = it.get_row(); + if (r_k != r_i) { + a_kj = it.get_row_entry().m_coeff; + a_kj.neg(); + M.mul(r_k, a_ij); + M.add(r_k, a_kj, r_i); + var_t s = m_row2base[r_k.id()]; + numeral& coeff = m_vars[s].m_base_coeff; + m.mul(coeff, a_ij, coeff); + M.gcd_normalize(r_k, g); + if (!m.is_one(g)) { + m.div(coeff, g, coeff); + } + } + } + } + + template + void simplex::update_value(var_t v, eps_numeral const& delta) { + update_value_core(v, delta); + col_iterator it = M.col_begin(v), end = M.col_end(v); + + // v <- v + delta + // s*s_coeff + v*v_coeff + R = 0 + // -> + // (v + delta)*v_coeff + (s - delta*v_coeff/s_coeff)*v + R = 0 + for (; it != end; ++it) { + row r = it.get_row(); + var_t s = m_row2base[r.id()]; + var_info& si = m_vars[s]; + scoped_eps_numeral delta2(em); + numeral const& coeff = it.get_row_entry().m_coeff; + em.mul(delta, coeff, delta2); + em.div(delta2, si.m_base_coeff, delta2); + delta2.neg(); + update_value_core(s, delta2); + // TBD m.add(si.m_base_coeff, delta2, si.m_base_coeff); + } + } + + template + void simplex::update_value_core(var_t v, eps_numeral const& delta) { + eps_numeral& val = m_vars[v].m_value; + em.add(val, delta, val); + if (is_base(v) && outside_bounds(v)) { + m_to_patch.insert(v); + } + } + + template + bool simplex::below_lower(var_t v) const { + var_info const& vi = m_vars[v]; + return vi.m_lower_valid && em.lt(vi.m_value, vi.m_lower); + } + + template + bool simplex::above_upper(var_t v) const { + var_info const& vi = m_vars[v]; + return vi.m_upper_valid && em.gt(vi.m_value, vi.m_upper); + } + + template + bool simplex::above_lower(var_t v) const { + var_info const& vi = m_vars[v]; + return !vi.m_lower_valid || em.gt(vi.m_value, vi.m_lower); + } + + template + bool simplex::below_upper(var_t v) const { + var_info const& vi = m_vars[v]; + return !vi.m_upper_valid || em.lt(vi.m_value, vi.m_upper); + } + + template + bool simplex::make_var_feasible(var_t x_i) { + scoped_numeral a_ij(m); + scoped_eps_numeral value(em); + bool is_below; + if (below_lower(x_i)) { + is_below = true; + value = m_vars[x_i].m_lower; + } + else if (above_upper(x_i)) { + is_below = false; + value = m_vars[x_i].m_upper; + } + else { + // x_i is already feasible + return true; + } + var_t x_j = select_pivot(x_i, is_below, a_ij); + if (x_j != null_var) { + update_and_pivot(x_i, x_j, a_ij, value); + } + return x_j != null_var; + } + + /** + \brief Wrapper for select_blands_pivot_core and select_pivot_core + */ + template + typename simplex::var_t + simplex::select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij) { + if (m_bland) { + return select_blands_pivot(x_i, is_below, out_a_ij); + } + if (is_below) { + return select_pivot_core(x_i, out_a_ij); + } + else { + return select_pivot_core(x_i, out_a_ij); + } + } + + /** + \brief Select a variable x_j in the row r defining the base var x_i, + s.t. x_j can be used to patch the error in x_i. Return null_theory_var + if there is no x_j. Otherwise, return x_j and store its coefficient + in out_a_ij. + + The argument is_below is true (false) if x_i is below its lower + bound (above its upper bound). + */ + template + template + typename simplex::var_t + simplex::select_pivot_core(var_t x_i, scoped_numeral & out_a_ij) { + SASSERT(is_base(x_i)); + var_t max = get_num_vars(); + var_t result = max; + row r = row(m_vars[x_i].m_base2row); + int n; + unsigned best_col_sz = UINT_MAX; + int best_so_far = INT_MAX; + + row_iterator it = M.row_begin(r), end = M.row_end(r); + + for (; it != end; ++it) { + var_t x_j = it->m_var; + numeral const & a_ij = it->m_coeff; + + bool is_neg = is_below ? m.is_neg(a_ij) : m.is_pos(a_ij); + bool is_pos = !is_neg; + if (x_i != x_j && ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { + int num = get_num_non_free_dep_vars(x_j, best_so_far); + unsigned col_sz = M.column_size(x_j); + if (num < best_so_far || (num == best_so_far && col_sz < best_col_sz)) { + result = x_j; + out_a_ij = a_ij; + best_so_far = num; + best_col_sz = col_sz; + n = 1; + } + else if (num == best_so_far && col_sz == best_col_sz) { + n++; + if (m_random()%n == 0) { + result = x_j; + out_a_ij = a_ij; + } + } + } + } + return result < max ? result : null_var; + } + + /** + \brief Return the number of base variables that are non free and are v dependent. + The function adds 1 to the result if v is non free. + The function returns with a partial result r if r > best_so_far. + This function is used to select the pivot variable. + */ + template + int simplex::get_num_non_free_dep_vars(var_t x_j, int best_so_far) { + int result = is_non_free(x_j); + col_iterator it = M.col_begin(x_j), end = M.col_end(x_j); + for (; it != end; ++it) { + var_t s = m_row2base[it.get_row().id()]; + result += is_non_free(s); + if (result > best_so_far) + return result; + } + return result; + } + + + /** + \brief Using Bland's rule, select a variable x_j in the row r defining the base var x_i, + s.t. x_j can be used to patch the error in x_i. Return null_theory_var + if there is no x_j. Otherwise, return x_j and store its coefficient + in out_a_ij. + */ + template + typename simplex::var_t + simplex::select_blands_pivot(var_t x_i, bool is_below, scoped_numeral & out_a_ij) { + SASSERT(is_base(x_i)); + unsigned max = get_num_vars(); + var_t result = max; + row r(m_vars[x_i].m_base2row); + row_iterator it = M.row_begin(r), end = M.row_end(r); + for (; it != end; ++it) { + var_t x_j = it->m_var; + numeral const & a_ij = it->m_coeff; + bool is_neg = is_below ? m.is_neg(a_ij) : m.is_pos(a_ij); + bool is_pos = !is_neg; + if (x_i != x_j && ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { + SASSERT(!is_base(x_j)); + if (x_j < result) { + result = x_j; + out_a_ij = a_ij; + } + } + } + return result < max ? result : null_var; + } + + template + lbool simplex::optimize(var_t v) { + // TBD SASSERT(is_feasible()); + // pick row for v and check if primal + // bounds are slack. + // return l_false for unbounded. + // return l_undef for canceled. + // return l_true for optimal. + return l_true; + } + + template + typename simplex::pivot_strategy_t + simplex::pivot_strategy() { + if (m_bland) { + return S_BLAND; + } + return S_DEFAULT; + } + + template + typename simplex::var_t + simplex::select_var_to_fix() { + switch (pivot_strategy()) { + case S_BLAND: + return select_smallest_var(); + case S_GREATEST_ERROR: + return select_error_var(false); + case S_LEAST_ERROR: + return select_error_var(true); + default: + return select_smallest_var(); + } + } + + template + typename simplex::var_t + simplex::select_error_var(bool least) { + var_t best = null_var; + scoped_eps_numeral best_error(em); + scoped_eps_numeral curr_error(em); + typename var_heap::iterator it = m_to_patch.begin(); + typename var_heap::iterator end = m_to_patch.end(); + for (; it != end; ++it) { + var_t v = *it; + var_info const& vi = m_vars[v]; + if (below_lower(v)) + em.sub(vi.m_lower, vi.m_value, curr_error); + else if (above_upper(v)) + em.sub(vi.m_value, vi.m_lower, curr_error); + else + continue; + SASSERT(is_pos(curr_error)); + if ((best == null_var) || + (!least && curr_error > best_error) || + (least && curr_error < best_error)) { + best = v; + best_error = curr_error; + } + } + if (best == null_var) + m_to_patch.clear(); // all variables are satisfied + else + m_to_patch.erase(best); + return best; + } + + template + bool simplex::well_formed() const { + SASSERT(M.well_formed()); + for (unsigned i = 0; i < m_row2base.size(); ++i) { + var_t s = m_row2base[i]; + SASSERT(i == m_vars[s].m_base2row); // unless row is deleted. + // + // TBD: extract coefficient of base variable and compare + // with m_vars[s].m_base_coeff; + // + } + return true; + } + + +}; + +#endif + diff --git a/src/math/simplex/sparse_matrix.h b/src/math/simplex/sparse_matrix.h new file mode 100644 index 000000000..cecd1c44a --- /dev/null +++ b/src/math/simplex/sparse_matrix.h @@ -0,0 +1,254 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + sparse_matrix.h + +Abstract: + + +Author: + + Nikolaj Bjorner (nbjorner) 2014-01-15 + +Notes: + +--*/ + +#ifndef _SPARSE_MATRIX_H_ +#define _SPARSE_MATRIX_H_ + +#include "mpq_inf.h" + +namespace simplex { + + template + class sparse_matrix { + public: + typedef typename Ext::numeral numeral; + typedef typename Ext::scoped_numeral scoped_numeral; + typedef typename Ext::manager manager; + typedef unsigned var_t; + + struct row_entry { + numeral m_coeff; + var_t m_var; + row_entry(numeral& c, var_t v): m_coeff(c), m_var(v) {} + }; + + private: + + static const int dead_id = -1; + + /** + \brief A row_entry is: m_var*m_coeff + + m_col_idx points to the place in the + column where the variable occurs. + */ + struct _row_entry : public row_entry { + union { + int m_col_idx; + int m_next_free_row_entry_idx; + }; + _row_entry(numeral const & c, var_t v): row_entry(c, v), m_col_idx(0) {} + _row_entry() : row_entry(numeral(), dead_id), m_col_idx(0) {} + bool is_dead() const { return m_var == dead_id; } + }; + + /** + \brief A column entry points to the row and the row_entry within the row + that has a non-zero coefficient on the variable associated + with the column entry. + */ + struct col_entry { + int m_row_id; + union { + int m_row_idx; + int m_next_free_col_entry_idx; + }; + col_entry(int r, int i): m_row_id(r), m_row_idx(i) {} + col_entry(): m_row_id(0), m_row_idx(0) {} + bool is_dead() const { return m_row_id == dead_id; } + }; + + struct column; + + /** + \brief A row contains a base variable and set of + row_entries. The base variable must occur in the set of + row_entries with coefficient 1. + */ + struct _row { + vector<_row_entry> m_entries; + unsigned m_size; // the real size, m_entries contains dead row_entries. + int m_first_free_idx; // first available position. + _row(); + unsigned size() const { return m_size; } + unsigned num_entries() const { return m_entries.size(); } + void reset(manager& m); + _row_entry & add_row_entry(unsigned & pos_idx); + void del_row_entry(unsigned idx); + void compress(manager& m, vector & cols); + void compress_if_needed(manager& m, vector & cols); + void save_var_pos(svector & result_map) const; + void reset_var_pos(svector & result_map) const; + bool is_coeff_of(var_t v, numeral const & expected) const; + int get_idx_of(var_t v) const; + }; + + /** + \brief A column stores in which rows a variable occurs. + The column may have free/dead entries. The field m_first_free_idx + is a reference to the first free/dead entry. + */ + struct column { + svector m_entries; + unsigned m_size; + int m_first_free_idx; + + column():m_size(0), m_first_free_idx(-1) {} + unsigned size() const { return m_size; } + unsigned num_entries() const { return m_entries.size(); } + void reset(); + void compress(vector<_row> & rows); + void compress_if_needed(vector<_row> & rows); + void compress_singleton(vector<_row> & rows, unsigned singleton_pos); + col_entry const * get_first_col_entry() const; + col_entry & add_col_entry(int & pos_idx); + void del_col_entry(unsigned idx); + }; + + manager m; + vector<_row> m_rows; + svector m_dead_rows; // rows to recycle + vector m_columns; // per var + svector m_var_pos; // temporary map from variables to positions in row + + bool well_formed_row(unsigned row_id) const; + bool well_formed_column(unsigned column_id) const; + void del_row_entry(_row& r, unsigned pos); + + public: + + sparse_matrix() {} + ~sparse_matrix(); + + class row { + unsigned m_id; + public: + row(unsigned r):m_id(r) {} + row():m_id(UINT_MAX) {} + bool operator!=(row const& other) const { + return m_id != other.m_id; + } + unsigned id() const { return m_id; } + }; + + void ensure_var(var_t v); + + row mk_row(); + void add(row r, numeral const& n, var_t var); + void add(row r, numeral const& n, row src); + void mul(row r, numeral const& n); + void neg(row r); + void del(row r); + + void gcd_normalize(row const& r, scoped_numeral& g); + + class row_iterator { + friend sparse_matrix; + unsigned m_curr; + _row & m_row; + void move_to_used() { + while (m_curr < m_row.m_entries.size() && + m_row.m_entries[m_curr].is_dead()) { + ++m_curr; + } + } + row_iterator(_row & r, bool begin): + m_curr(0), m_row(r) { + if (begin) { + move_to_used(); + } + else { + m_curr = m_row.m_entries.size(); + } + } + public: + row_entry & operator*() const { return m_row.m_entries[m_curr]; } + row_entry * operator->() const { return &(operator*()); } + row_iterator & operator++() { ++m_curr; move_to_used(); return *this; } + row_iterator operator++(int) { row_iterator tmp = *this; ++*this; return tmp; } + bool operator==(row_iterator const & it) const { return m_curr == it.m_curr; } + bool operator!=(row_iterator const & it) const { return m_curr != it.m_curr; } + }; + + row_iterator row_begin(row const& r) { return row_iterator(m_rows[r.id()], true); } + row_iterator row_end(row const& r) { return row_iterator(m_rows[r.id()], false); } + + unsigned column_size(var_t v) const { return m_columns[v].size(); } + + class col_iterator { + friend sparse_matrix; + unsigned m_curr; + column const& m_col; + vector<_row> const& m_rows; + void move_to_used() { + while (m_curr < m_col.m_entries.size() && m_col.m_entries[m_curr].is_dead()) { + ++m_curr; + } + } + col_iterator(column const& c, vector<_row> const& r, bool begin): + m_curr(0), m_col(c), m_rows(r) { + if (begin) { + move_to_used(); + } + else { + m_curr = m_col.m_entries.size(); + } + } + public: + row get_row() { return row(m_col.m_entries[m_curr].m_row_id); } + row_entry const& get_row_entry() { + col_entry const& c = m_col.m_entries[m_curr]; + int row_id = c.m_row_id; + return m_rows[row_id].m_entries[c.m_row_idx]; + } + + col_iterator & operator++() { ++m_curr; move_to_used(); return *this; } + col_iterator operator++(int) { col_iterator tmp = *this; ++*this; return tmp; } + bool operator==(col_iterator const & it) const { return m_curr == it.m_curr; } + bool operator!=(col_iterator const & it) const { return m_curr != it.m_curr; } + }; + + col_iterator col_begin(int v) const { return col_iterator(m_columns[v], m_rows, true); } + col_iterator col_end(int v) const { return col_iterator(m_columns[v], m_rows, false); } + + void display(std::ostream& out) const; + bool well_formed() const; + + + }; + + struct mpz_ext { + typedef mpz numeral; + typedef scoped_mpz scoped_numeral; + typedef unsynch_mpz_manager manager; + typedef mpq_inf eps_numeral; + typedef unsynch_mpq_inf_manager eps_manager; + }; + + struct mpq_ext { + typedef mpq numeral; + typedef scoped_mpq scoped_numeral; + typedef unsynch_mpq_manager manager; + typedef mpq_inf eps_numeral; + typedef unsynch_mpq_inf_manager eps_manager; + }; + +}; + + +#endif diff --git a/src/math/simplex/sparse_matrix_def.h b/src/math/simplex/sparse_matrix_def.h new file mode 100644 index 000000000..4e46c823d --- /dev/null +++ b/src/math/simplex/sparse_matrix_def.h @@ -0,0 +1,580 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + sparse_matrix_def.h + +Abstract: + + +Author: + + Nikolaj Bjorner (nbjorner) 2014-01-15 + +Notes: + + mainly hoisted from theory_arith.h and theory_arith_aux.h + +--*/ + +#ifndef _SPARSE_MATRIX_DEF_H_ +#define _SPARSE_MATRIX_DEF_H_ + +#include "sparse_matrix.h" +#include "uint_set.h" + +namespace simplex { + + // ----------------------------------- + // + // Rows + // + // ----------------------------------- + + template + sparse_matrix::_row::_row(): + m_size(0), + m_first_free_idx(-1) { + } + + template + void sparse_matrix::_row::reset(manager& m) { + for (unsigned i = 0; i < m_entries.size(); ++i) { + m.reset(m_entries[i].m_coeff); + } + m_entries.reset(); + m_size = 0; + m_first_free_idx = -1; + } + + /** + \brief Add a new row_entry. The result is a reference to the new row_entry. + The position of the new row_entry in the + row is stored in pos_idx. + */ + template + typename sparse_matrix::_row_entry & + sparse_matrix::_row::add_row_entry(unsigned & pos_idx) { + m_size++; + if (m_first_free_idx == -1) { + pos_idx = m_entries.size(); + m_entries.push_back(_row_entry()); + return m_entries.back(); + } + else { + SASSERT(m_first_free_idx >= 0); + pos_idx = static_cast(m_first_free_idx); + _row_entry & result = m_entries[pos_idx]; + SASSERT(result.is_dead()); + m_first_free_idx = result.m_next_free_row_entry_idx; + return result; + } + } + + /** + \brief Delete row_entry at position idx. + */ + template + void sparse_matrix::_row::del_row_entry(unsigned idx) { + _row_entry & t = m_entries[idx]; + SASSERT(!t.is_dead()); + t.m_next_free_row_entry_idx = m_first_free_idx; + t.m_var = dead_id; + m_size--; + m_first_free_idx = idx; + SASSERT(t.is_dead()); + } + + /** + \brief Remove holes (i.e., dead entries) from the row. + */ + template + void sparse_matrix::_row::compress(manager& m, vector & cols) { + unsigned i = 0; + unsigned j = 0; + unsigned sz = m_entries.size(); + for (; i < sz; i++) { + _row_entry & t1 = m_entries[i]; + if (!t1.is_dead()) { + if (i != j) { + _row_entry & t2 = m_entries[j]; + t2.m_coeff.swap(t1.m_coeff); + t2.m_var = t1.m_var; + t2.m_col_idx = t1.m_col_idx; + SASSERT(!t2.is_dead()); + column & col = cols[t2.m_var]; + col.m_entries[t2.m_col_idx].m_row_idx = j; + } + j++; + } + } + SASSERT(j == m_size); + // + // alternative: update the free-list to point to the + // tail and avoid shrinking. + // if m.does not allocate memory (for wrapper around + // double), also bypass this step. + // + for (unsigned i = m_size; i < m_entries.size(); ++i) { + m.reset(m_entries[i].m_coeff); + } + m_entries.shrink(m_size); + m_first_free_idx = -1; + } + + template + void sparse_matrix::_row::compress_if_needed(manager& m, vector & cols) { + if (size() *2 < num_entries()) { + compress(m, cols); + } + } + + /** + \brief Fill the map var -> pos/idx + */ + template + inline void sparse_matrix::_row::save_var_pos(svector & result_map) const { + typename vector<_row_entry>::const_iterator it = m_entries.begin(); + typename vector<_row_entry>::const_iterator end = m_entries.end(); + unsigned idx = 0; + for (; it != end; ++it, ++idx) { + if (!it->is_dead()) { + result_map[it->m_var] = idx; + } + } + } + + /** + \brief Reset the map var -> pos/idx. That is for all variables v in the row, set result[v] = -1 + This method can be viewed as the "inverse" of save_var_pos. + */ + template + inline void sparse_matrix::_row::reset_var_pos(svector & result_map) const { + typename vector<_row_entry>::const_iterator it = m_entries.begin(); + typename vector<_row_entry>::const_iterator end = m_entries.end(); + for (; it != end; ++it) { + if (!it->is_dead()) { + result_map[it->m_var] = -1; + } + } + } + + template + int sparse_matrix::_row::get_idx_of(var_t v) const { + typename vector<_row_entry>::const_iterator it = m_entries.begin(); + typename vector<_row_entry>::const_iterator end = m_entries.end(); + for (unsigned idx = 0; it != end; ++it, ++idx) { + if (!it->is_dead() && it->m_var == v) + return idx; + } + return -1; + } + + // ----------------------------------- + // + // Columns + // + // ----------------------------------- + + template + void sparse_matrix::column::reset() { + m_entries.reset(); + m_size = 0; + m_first_free_idx = -1; + } + + /** + \brief Remove holes (i.e., dead entries) from the column. + */ + template + void sparse_matrix::column::compress(vector<_row> & rows) { + unsigned i = 0; + unsigned j = 0; + unsigned sz = m_entries.size(); + for (; i < sz; i++) { + col_entry & e1 = m_entries[i]; + if (!e1.is_dead()) { + if (i != j) { + m_entries[j] = e1; + _row & r = rows[e1.m_row_id]; + r.m_entries[e1.m_row_idx].m_col_idx = j; + } + j++; + } + } + SASSERT(j == m_size); + m_entries.shrink(m_size); + m_first_free_idx = -1; + } + + /** + \brief Invoke compress if the column contains too many holes (i.e., dead entries). + */ + template + inline void sparse_matrix::column::compress_if_needed(vector<_row> & rows) { + if (size() * 2 < num_entries()) { + compress(rows); + } + } + + /** + \brief Special version of compress, that is used when the column contain + only one entry located at position singleton_pos. + */ + template + void sparse_matrix::column::compress_singleton(vector<_row> & rows, unsigned singleton_pos) { + SASSERT(m_size == 1); + if (singleton_pos != 0) { + col_entry & s = m_entries[singleton_pos]; + m_entries[0] = s; + row & r = rows[s.m_row_id]; + r[s.m_row_idx].m_col_idx = 0; + } + m_first_free_idx = -1; + m_entries.shrink(1); + } + + template + const typename sparse_matrix::col_entry * + sparse_matrix::column::get_first_col_entry() const { + typename svector::const_iterator it = m_entries.begin(); + typename svector::const_iterator end = m_entries.end(); + for (; it != end; ++it) { + if (!it->is_dead()) { + return it; + } + } + return 0; + } + + template + typename sparse_matrix::col_entry & + sparse_matrix::column::add_col_entry(int & pos_idx) { + m_size++; + if (m_first_free_idx == -1) { + pos_idx = m_entries.size(); + m_entries.push_back(col_entry()); + return m_entries.back(); + } + else { + pos_idx = m_first_free_idx; + col_entry & result = m_entries[pos_idx]; + SASSERT(result.is_dead()); + m_first_free_idx = result.m_next_free_col_entry_idx; + return result; + } + } + + template + void sparse_matrix::column::del_col_entry(unsigned idx) { + col_entry & c = m_entries[idx]; + SASSERT(!c.is_dead()); + c.m_row_id = dead_id; + c.m_next_free_col_entry_idx = m_first_free_idx; + m_first_free_idx = idx; + m_size--; + } + + // ----------------------------------- + // + // Matrix + // + // ----------------------------------- + + template + sparse_matrix::~sparse_matrix() { + for (unsigned i = 0; i < m_rows.size(); ++i) { + _row& r = m_rows[i]; + for (unsigned j = 0; j < r.m_entries.size(); ++j) { + m.reset(r.m_entries[j].m_coeff); + } + } + } + + template + void sparse_matrix::ensure_var(var_t v) { + while (m_columns.size() <= v) { + m_columns.push_back(column()); + m_var_pos.push_back(-1); + } + } + + template + typename sparse_matrix::row + sparse_matrix::mk_row() { + if (m_dead_rows.empty()) { + row r(m_rows.size()); + m_rows.push_back(_row()); + return r; + } + else { + row r(m_dead_rows.back()); + m_dead_rows.pop_back(); + return r; + } + } + + template + void sparse_matrix::add(row dst, numeral const& n, var_t v) { + _row& r = m_rows[dst.id()]; + column& c = m_columns[v]; + unsigned r_idx; + int c_idx; + _row_entry & r_entry = r.add_row_entry(r_idx); + col_entry& c_entry = c.add_col_entry(c_idx); + r_entry.m_var = v; + m.set(r_entry.m_coeff, n); + r_entry.m_col_idx = c_idx; + c_entry.m_row_id = dst.id(); + c_entry.m_row_idx = r_idx; + } + + /** + \brief Set row1 <- row1 + row2 * n + */ + template + void sparse_matrix::add(row row1, numeral const& n, row row2) { + // m_stats.m_add_rows++; + _row & r1 = m_rows[row1.id()]; + _row & r2 = m_rows[row2.id()]; + + r1.save_var_pos(m_var_pos); + + // + // loop over variables in row2, + // add terms in row2 to row1. + // + +#define ADD_ROW(_SET_COEFF_, _ADD_COEFF_) \ + row_iterator it = row_begin(row2); \ + row_iterator end = row_end(row2); \ + for (; it != end; ++it) { \ + var_t v = it->m_var; \ + int pos = m_var_pos[v]; \ + if (pos == -1) { \ + /* variable v is not in row1 */ \ + unsigned row_idx; \ + _row_entry & r_entry = r1.add_row_entry(row_idx); \ + r_entry.m_var = v; \ + m.set(r_entry.m_coeff, it->m_coeff); \ + _SET_COEFF_; \ + column & c = m_columns[v]; \ + int col_idx; \ + col_entry & c_entry = c.add_col_entry(col_idx); \ + r_entry.m_col_idx = col_idx; \ + c_entry.m_row_id = row1.id(); \ + c_entry.m_row_idx = row_idx; \ + } \ + else { \ + /* variable v is in row1 */ \ + _row_entry & r_entry = r1.m_entries[pos]; \ + SASSERT(r_entry.m_var == v); \ + _ADD_COEFF_; \ + if (m.is_zero(r_entry.m_coeff)) { \ + del_row_entry(r1, pos); \ + } \ + } \ + } \ + + ((void) 0) + + if (m.is_one(n)) { + ADD_ROW({}, + m.add(r_entry.m_coeff, it->m_coeff, r_entry.m_coeff)); + } + else if (m.is_minus_one(n)) { + ADD_ROW(m.neg(r_entry.m_coeff), + m.sub(r_entry.m_coeff, it->m_coeff, r_entry.m_coeff)); + } + else { + scoped_numeral tmp(m); + ADD_ROW(m.mul(r_entry.m_coeff, n, r_entry.m_coeff), + m.mul(it->m_coeff, n, tmp); + m.add(r_entry.m_coeff, tmp, r_entry.m_coeff)); + } + + r1.reset_var_pos(m_var_pos); + r1.compress_if_needed(m, m_columns); + } + + + template + void sparse_matrix::del_row_entry(_row& r, unsigned pos) { + _row_entry & r_entry = r.m_entries[pos]; + var_t v = r_entry.m_var; + int col_idx = r_entry.m_col_idx; + r.del_row_entry(pos); + column & c = m_columns[v]; + c.del_col_entry(col_idx); + c.compress_if_needed(m_rows); + } + + /** + \brief Set row <- -row + */ + template + void sparse_matrix::neg(row r) { + row_iterator it = row_begin(r); + row_iterator end = row_end(r); + for (; it != end; ++it) { + m.neg(it->m_coeff); + } + } + + /** + \brief Set row <- n*row + */ + template + void sparse_matrix::mul(row r, numeral const& n) { + SASSERT(!m.is_zero(n)); + if (m.is_one(n)) { + // no op + } + else if (m.is_minus_one(n)) { + neg(r); + } + else { + row_iterator it = row_begin(r); + row_iterator end = row_end(r); + for (; it != end; ++it) { + m.mul(it->m_coeff, n, it->m_coeff); + } + } + } + + /** + \brief Delete row. + */ + template + void sparse_matrix::del(row r) { + _row& rw = r.id(); + for (unsigned i = 0; i < rw.m_entries.size(); ++i) { + _row_entry& e = rw.m_entries[i]; + if (!e.is_dead()) { + del_row_entry(rw, i); + } + } + SASSERT(rw.m_size == 0); + m_dead_rows.push_back(r.id()); + } + + /** + \brief normalize coefficients by dividing with they coefficients. + return the gcd. + */ + template + void sparse_matrix::gcd_normalize(row const& r, scoped_numeral& g) { + g.reset(); + row_iterator it = row_begin(r), end = row_end(r); + for (; it != end && !m.is_one(g); ++it) { + if (m.is_zero(g)) g = it->m_coeff; + else m.gcd(g, it->m_coeff, g); + } + if (m.is_zero(g)) { + g = numeral(1); + } + if (!m.is_one(g)) { + row_iterator it2 = row_begin(r); + for (; it2 != end; ++it2) { + m.div(it2->m_coeff, g, it2->m_coeff); + } + } + } + + /** + \brief well_formed check + */ + template + bool sparse_matrix::well_formed() const { + for (unsigned i = 0; i < m_rows.size(); ++i) { + well_formed_row(i); + } + for (unsigned i = 0; i < m_columns.size(); ++i) { + well_formed_column(i); + } + return true; + } + + /** + \brief well_formed row check + */ + template + bool sparse_matrix::well_formed_row(unsigned row_id) const { + uint_set vars, dead; + _row const& r = m_rows[row_id]; + for (unsigned i = 0; i < r.num_entries(); ++i) { + _row_entry const& e = r.m_entries[i]; + if (e.is_dead()) { + dead.insert(i); + continue; + } + SASSERT(!vars.contains(e.m_var)); + SASSERT(!m.is_zero(e.m_coeff)); + SASSERT(e.m_var != dead_id); + col_entry const& c = m_columns[e.m_var].m_entries[e.m_col_idx]; + SASSERT(c.m_row_id == row_id); + SASSERT(c.m_row_idx == i); + vars.insert(e.m_var); + } + int idx = r.m_first_free_idx; + while (idx != -1) { + SASSERT(dead.contains(idx)); + dead.remove(idx); + idx = r.m_entries[idx].m_next_free_row_entry_idx; + } + SASSERT(dead.empty()); + return true; + } + + /** + \brief well_formed column check + */ + template + bool sparse_matrix::well_formed_column(var_t v) const { + uint_set rows, dead; + column const& col = m_columns[v]; + for (unsigned i = 0; i < col.num_entries(); ++i) { + col_entry const& c = col.m_entries[i]; + if (c.is_dead()) { + dead.insert(i); + continue; + } + SASSERT(!rows.contains(c.m_row_id)); + _row const& row = m_rows[c.m_row_id]; + _row_entry const& r = row.m_entries[c.m_row_idx]; + SASSERT(r.m_var == v); + SASSERT(r.m_col_idx == i); + rows.insert(c.m_row_id); + } + int idx = col.m_first_free_idx; + while (idx != -1) { + SASSERT(dead.contains(idx)); + dead.remove(idx); + idx = col.m_entries[idx].m_next_free_col_entry_idx; + } + SASSERT(dead.empty()); + return true; + } + + /** + \brief display method + */ + template + void sparse_matrix::display(std::ostream& out) const { + for (unsigned i = 0; i < m_rows.size(); ++i) { + if (m_rows[i].size() == 0) continue; + row_iterator it = row_begin(row(i)), end = row_end(row(i)); + for (; it != end; ++it) { + m.display(out, it->m_coeff) + out << "*v" << it->m_var << " "; + } + out << "\n"; + } + } + + + +}; + +#endif diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp new file mode 100644 index 000000000..ca69cec1f --- /dev/null +++ b/src/test/simplex.cpp @@ -0,0 +1,24 @@ +#include "sparse_matrix.h" +#include "sparse_matrix_def.h" +#include "simplex.h" +#include "simplex_def.h" + +typedef simplex::simplex Simplex; + +void tst_simplex() { + simplex::sparse_matrix M; + Simplex S; + + S.make_feasible(); + + unsynch_mpz_manager m; + scoped_mpz_vector coeffs(m); + svector vars; + for (unsigned i = 0; i < 5; ++i) { + S.ensure_var(i); + vars.push_back(i); + coeffs.push_back(mpz(i+1)); + } + + Simplex::row r = S.add_row(1, coeffs.size(), vars.c_ptr(), coeffs.c_ptr()); +} diff --git a/src/util/mpq_inf.cpp b/src/util/mpq_inf.cpp index 6f6d805f2..cb3232d6c 100644 --- a/src/util/mpq_inf.cpp +++ b/src/util/mpq_inf.cpp @@ -19,7 +19,7 @@ Revision History: #include"mpq_inf.h" template -std::string mpq_inf_manager::to_string(mpq_inf const & a) const { +std::string mpq_inf_manager::to_string(mpq_inf const & a) { if (m.is_zero(a.second)) return m.to_string(a.first); diff --git a/src/util/mpq_inf.h b/src/util/mpq_inf.h index 217b58567..eef066a6b 100644 --- a/src/util/mpq_inf.h +++ b/src/util/mpq_inf.h @@ -26,12 +26,12 @@ typedef std::pair mpq_inf; template class mpq_inf_manager { - mpq_manager & m; + mpq_manager m; double m_inf; public: typedef mpq_inf numeral; - mpq_inf_manager(mpq_manager & _m, double inf = 0.0001):m(_m) { + mpq_inf_manager(double inf = 0.0001) { set_inf(inf); } @@ -83,6 +83,14 @@ public: } bool is_int(mpq_inf const & a) const { return m.is_int(a.first) && m.is_zero(a.second); } + + bool is_pos(mpq_inf const & a) const { + return m.is_pos(a.first) || (m.is_zero(a.first) && m.is_pos(a.second)); + } + + bool is_neg(mpq_inf const & a) const { + return m.is_neg(a.first) || (m.is_zero(a.first) && m.is_neg(a.second)); + } bool is_rational(mpq_inf const & a) const { return m.is_zero(a.second); } @@ -104,15 +112,15 @@ public: return m.is_zero(a.first) && m.is_zero(a.second); } - bool eq(mpq_inf const & a, mpq_inf const & b) const { + bool eq(mpq_inf const & a, mpq_inf const & b) { return m.eq(a.first, b.first) && m.eq(a.second, b.second); } - bool eq(mpq_inf const & a, mpq const & b) const { + bool eq(mpq_inf const & a, mpq const & b) { return m.eq(a.first, b) && m.is_zero(a.second); } - bool eq(mpq_inf const & a, mpq const & b, inf_kind k) const { + bool eq(mpq_inf const & a, mpq const & b, inf_kind k) { if (!m.eq(a.first, b)) return false; switch (k) { @@ -124,15 +132,15 @@ public: return false; } - bool lt(mpq_inf const & a, mpq_inf const & b) const { + bool lt(mpq_inf const & a, mpq_inf const & b) { return m.lt(a.first, b.first) || (m.lt(a.second, b.second) && m.eq(a.first, b.first)); } - bool lt(mpq_inf const & a, mpq const & b) const { + bool lt(mpq_inf const & a, mpq const & b) { return m.lt(a.first, b) || (m.is_neg(a.second) && m.eq(a.first, b)); } - bool lt(mpq_inf const & a, mpq const & b, inf_kind k) const { + bool lt(mpq_inf const & a, mpq const & b, inf_kind k) { if (m.lt(a.first, b)) return true; if (m.eq(a.first, b)) { @@ -146,13 +154,13 @@ public: return false; } - bool gt(mpq_inf const & a, mpq_inf const & b) const { return lt(b, a); } + bool gt(mpq_inf const & a, mpq_inf const & b) { return lt(b, a); } - bool gt(mpq_inf const & a, mpq const & b) const { + bool gt(mpq_inf const & a, mpq const & b) { return m.gt(a.first, b) || (m.is_pos(a.second) && m.eq(a.first, b)); } - bool gt(mpq_inf const & a, mpq const & b, inf_kind k) const { + bool gt(mpq_inf const & a, mpq const & b, inf_kind k) { if (m.gt(a.first, b)) return true; if (m.eq(a.first, b)) { @@ -166,17 +174,17 @@ public: return false; } - bool le(mpq_inf const & a, mpq_inf const & b) const { return !gt(a, b); } + bool le(mpq_inf const & a, mpq_inf const & b) { return !gt(a, b); } - bool le(mpq_inf const & a, mpq const & b) const { return !gt(a, b); } + bool le(mpq_inf const & a, mpq const & b) { return !gt(a, b); } - bool le(mpq_inf const & a, mpq const & b, inf_kind k) const { return !gt(a, b, k); } + bool le(mpq_inf const & a, mpq const & b, inf_kind k) { return !gt(a, b, k); } - bool ge(mpq_inf const & a, mpq_inf const & b) const { return !lt(a, b); } + bool ge(mpq_inf const & a, mpq_inf const & b) { return !lt(a, b); } - bool ge(mpq_inf const & a, mpq const & b) const { return !lt(a, b); } + bool ge(mpq_inf const & a, mpq const & b) { return !lt(a, b); } - bool ge(mpq_inf const & a, mpq const & b, inf_kind k) const { return !lt(a, b, k); } + bool ge(mpq_inf const & a, mpq const & b, inf_kind k) { return !lt(a, b, k); } void add(mpq_inf const & a, mpq_inf const & b, mpq_inf & c) { m.add(a.first, b.first, c.first); @@ -208,6 +216,16 @@ public: m.mul(b, a.second, c.second); } + void div(mpq_inf const & a, mpq const & b, mpq_inf & c) { + m.div(a.first, b, c.first); + m.div(a.second, b, c.second); + } + + void div(mpq_inf const & a, mpz const & b, mpq_inf & c) { + m.div(a.first, b, c.first); + m.div(a.second, b, c.second); + } + void inc(mpq_inf & a) { m.inc(a.first); } @@ -246,7 +264,7 @@ public: } } - std::string to_string(mpq_inf const & a) const; + std::string to_string(mpq_inf const & a); }; typedef mpq_inf_manager synch_mpq_inf_manager;