From 9ba4026bc6b78ff3ba9b56444bfe1b411f4a3863 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Thu, 31 May 2018 10:37:22 -0700 Subject: [PATCH] avoid going creating hnf_cuts if all involved vars have integral values Signed-off-by: Lev Nachmanson add explanations to hnf cuts Signed-off-by: Lev Nachmanson nits and virtual methods (#68) * local Signed-off-by: Nikolaj Bjorner * virtual method in bound propagator Signed-off-by: Nikolaj Bjorner cleanup from std::cout Signed-off-by: Lev Nachmanson handle the case when the number of terms is greater than the number of variables in hnf Signed-off-by: Lev Nachmanson method name's fix Signed-off-by: Lev Nachmanson restore hnf_cutter to work with m_row_count <= m_column_count Signed-off-by: Lev Nachmanson tune addition of rational numbers (#70) * log quantifiers only if present Signed-off-by: Nikolaj Bjorner * merge and fix some warnings Signed-off-by: Nikolaj Bjorner * set new arith as default for LIA Signed-off-by: Nikolaj Bjorner * local Signed-off-by: Nikolaj Bjorner * virtual method in bound propagator Signed-off-by: Nikolaj Bjorner * prepare for mixed integer-real Signed-off-by: Nikolaj Bjorner * fix default tactic usage Signed-off-by: Nikolaj Bjorner give shorter explanations, call hnf only when have a not integral var Signed-off-by: Lev Nachmanson overhaul of mpq (#71) * log quantifiers only if present Signed-off-by: Nikolaj Bjorner * merge and fix some warnings Signed-off-by: Nikolaj Bjorner * set new arith as default for LIA Signed-off-by: Nikolaj Bjorner * local Signed-off-by: Nikolaj Bjorner * virtual method in bound propagator Signed-off-by: Nikolaj Bjorner * prepare for mixed integer-real Signed-off-by: Nikolaj Bjorner * fix default tactic usage Signed-off-by: Nikolaj Bjorner * overhaul of mpz/mpq Signed-off-by: Nikolaj Bjorner * disabled temporary setting Signed-off-by: Nikolaj Bjorner * remove prints Signed-off-by: Nikolaj Bjorner fix for 32 bit build (#72) * log quantifiers only if present Signed-off-by: Nikolaj Bjorner * merge and fix some warnings Signed-off-by: Nikolaj Bjorner * set new arith as default for LIA Signed-off-by: Nikolaj Bjorner * local Signed-off-by: Nikolaj Bjorner * virtual method in bound propagator Signed-off-by: Nikolaj Bjorner * prepare for mixed integer-real Signed-off-by: Nikolaj Bjorner * fix default tactic usage Signed-off-by: Nikolaj Bjorner * overhaul of mpz/mpq Signed-off-by: Nikolaj Bjorner * disabled temporary setting Signed-off-by: Nikolaj Bjorner * remove prints Signed-off-by: Nikolaj Bjorner * customize for 64 bit Signed-off-by: Nikolaj Bjorner yes (#74) * log quantifiers only if present Signed-off-by: Nikolaj Bjorner * merge and fix some warnings Signed-off-by: Nikolaj Bjorner * set new arith as default for LIA Signed-off-by: Nikolaj Bjorner * local Signed-off-by: Nikolaj Bjorner * virtual method in bound propagator Signed-off-by: Nikolaj Bjorner * prepare for mixed integer-real Signed-off-by: Nikolaj Bjorner * fix default tactic usage Signed-off-by: Nikolaj Bjorner * overhaul of mpz/mpq Signed-off-by: Nikolaj Bjorner * disabled temporary setting Signed-off-by: Nikolaj Bjorner * remove prints Signed-off-by: Nikolaj Bjorner * customize for 64 bit Signed-off-by: Nikolaj Bjorner * customize for 64 bit Signed-off-by: Nikolaj Bjorner * more refactor Signed-off-by: Nikolaj Bjorner fix the merge Signed-off-by: Lev Nachmanson fixes in maximize_term untested Signed-off-by: Lev Nachmanson fix compilation (#75) * log quantifiers only if present Signed-off-by: Nikolaj Bjorner * merge and fix some warnings Signed-off-by: Nikolaj Bjorner * set new arith as default for LIA Signed-off-by: Nikolaj Bjorner * local Signed-off-by: Nikolaj Bjorner * virtual method in bound propagator Signed-off-by: Nikolaj Bjorner * prepare for mixed integer-real Signed-off-by: Nikolaj Bjorner * fix default tactic usage Signed-off-by: Nikolaj Bjorner * overhaul of mpz/mpq Signed-off-by: Nikolaj Bjorner * disabled temporary setting Signed-off-by: Nikolaj Bjorner * remove prints Signed-off-by: Nikolaj Bjorner * customize for 64 bit Signed-off-by: Nikolaj Bjorner * customize for 64 bit Signed-off-by: Nikolaj Bjorner * more refactor Signed-off-by: Nikolaj Bjorner * merge Signed-off-by: Nikolaj Bjorner * relax check Signed-off-by: Nikolaj Bjorner * change for gcc Signed-off-by: Nikolaj Bjorner * merge Signed-off-by: Nikolaj Bjorner --- src/smt/params/smt_params_helper.pyg | 2 +- src/smt/smt_setup.cpp | 5 +- src/smt/theory_lra.cpp | 23 +- src/tactic/smtlogics/qflia_tactic.cpp | 12 - src/test/lp/lp.cpp | 87 +-- src/test/simplex.cpp | 1 - src/util/lp/bound_propagator.h | 2 +- src/util/lp/explanation.h | 1 + src/util/lp/hnf.h | 87 ++- src/util/lp/hnf_cutter.h | 50 +- src/util/lp/int_solver.cpp | 32 +- src/util/lp/int_solver.h | 10 +- src/util/lp/lar_core_solver.h | 2 - src/util/lp/lar_solver.cpp | 98 +-- src/util/lp/lar_solver.h | 12 +- src/util/lp/lp_core_solver_base.cpp | 1 + src/util/lp/lp_core_solver_base.h | 12 +- src/util/lp/lp_core_solver_base_def.h | 39 +- src/util/lp/lp_dual_core_solver_def.h | 6 - src/util/lp/lp_primal_core_solver.h | 3 +- .../lp/lp_primal_core_solver_tableau_def.h | 1 - src/util/lp/lp_primal_simplex_def.h | 1 - src/util/lp/lp_settings.h | 5 +- src/util/lp/lp_settings_def.h | 4 - src/util/lp/lu_def.h | 5 - src/util/lp/mps_reader.h | 1 - src/util/lp/permutation_matrix_def.h | 1 - src/util/lp/scaler_def.h | 2 - src/util/lp/square_sparse_matrix.h | 4 - src/util/lp/square_sparse_matrix_def.h | 9 - src/util/lp/stacked_map.h | 255 ------- src/util/lp/stacked_unordered_set.h | 107 --- src/util/lp/static_matrix_def.h | 7 - src/util/lp/test_bound_analyzer.h | 2 - src/util/lp/var_register.h | 2 + src/util/mpff.cpp | 6 +- src/util/mpfx.cpp | 2 +- src/util/mpq.cpp | 106 +++ src/util/mpq.h | 71 +- src/util/mpz.cpp | 689 ++++++++++++------ src/util/mpz.h | 439 ++++------- 41 files changed, 933 insertions(+), 1271 deletions(-) delete mode 100644 src/util/lp/stacked_map.h delete mode 100644 src/util/lp/stacked_unordered_set.h diff --git a/src/smt/params/smt_params_helper.pyg b/src/smt/params/smt_params_helper.pyg index 816764896a..9cb485668f 100644 --- a/src/smt/params/smt_params_helper.pyg +++ b/src/smt/params/smt_params_helper.pyg @@ -40,7 +40,7 @@ def_module_params(module_name='smt', ('bv.reflect', BOOL, True, 'create enode for every bit-vector term'), ('bv.enable_int2bv', BOOL, True, 'enable support for int2bv and bv2int operators'), ('arith.random_initial_value', BOOL, False, 'use random initial values in the simplex-based procedure for linear arithmetic'), - ('arith.solver', UINT, 2, 'arithmetic solver: 0 - no solver, 1 - bellman-ford based solver (diff. logic only), 2 - simplex based solver, 3 - floyd-warshall based solver (diff. logic only) and no theory combination 4 - utvpi, 5 - infinitary lra, 6 - lra solver'), + ('arith.solver', UINT, 6, 'arithmetic solver: 0 - no solver, 1 - bellman-ford based solver (diff. logic only), 2 - simplex based solver, 3 - floyd-warshall based solver (diff. logic only) and no theory combination 4 - utvpi, 5 - infinitary lra, 6 - lra solver'), ('arith.nl', BOOL, True, '(incomplete) nonlinear arithmetic support based on Groebner basis and interval propagation'), ('arith.nl.gb', BOOL, True, 'groebner Basis computation, this option is ignored when arith.nl=false'), ('arith.nl.branching', BOOL, True, 'branching on integer variables in non linear clusters'), diff --git a/src/smt/smt_setup.cpp b/src/smt/smt_setup.cpp index 9390d91ba5..6a689c8d98 100644 --- a/src/smt/smt_setup.cpp +++ b/src/smt/smt_setup.cpp @@ -818,10 +818,7 @@ namespace smt { m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params)); break; default: - if (m_params.m_arith_int_only && int_only) - setup_i_arith(); - else - m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params)); + setup_i_arith(); break; } } diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 07a8aecaed..7e2ee2a7e9 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -751,12 +751,14 @@ class theory_lra::imp { init_solver(); } - bool internalize_atom(app * atom, bool gate_ctx) { - return internalize_atom_strict(atom, gate_ctx); - + void internalize_is_int(app * n) { + SASSERT(a.is_is_int(n)); + (void) mk_enode(n); + if (!ctx().relevancy()) + mk_is_int_axiom(n); } - bool internalize_atom_strict(app * atom, bool gate_ctx) { + bool internalize_atom(app * atom, bool gate_ctx) { SASSERT(!ctx().b_internalized(atom)); bool_var bv = ctx().mk_bool_var(atom); ctx().set_var_theory(bv, get_id()); @@ -772,6 +774,10 @@ class theory_lra::imp { v = internalize_def(to_app(n1)); k = lp_api::lower_t; } + else if (a.is_is_int(atom)) { + internalize_is_int(atom); + return true; + } else { TRACE("arith", tout << "Could not internalize " << mk_pp(atom, m) << "\n";); found_not_handled(atom); @@ -1527,7 +1533,7 @@ class theory_lra::imp { return m_imp.bound_is_interesting(j, kind, v); } - virtual void consume(rational const& v, unsigned j) { + void consume(rational const& v, lp::constraint_index j) override { m_imp.set_evidence(j); m_imp.m_explanation.push_back(std::make_pair(v, j)); } @@ -2619,12 +2625,17 @@ class theory_lra::imp { coeff = rational::zero(); } lp::impq term_max; - if (m_solver->maximize_term(coeffs, term_max)) { + lp::lp_status st = m_solver->maximize_term(coeffs, term_max); + if (st == lp::lp_status::OPTIMAL) { blocker = mk_gt(v); inf_rational val(term_max.x + coeff, term_max.y); return inf_eps(rational::zero(), val); + } if (st == lp::lp_status::FEASIBLE) { + lp_assert(false); // not implemented + // todo: what to return? } else { + SASSERT(st == lp::lp_status::UNBOUNDED); TRACE("arith", tout << "Unbounded " << v << "\n";); has_shared = false; blocker = m.mk_false(); diff --git a/src/tactic/smtlogics/qflia_tactic.cpp b/src/tactic/smtlogics/qflia_tactic.cpp index b6b0cd6bc1..46b766bd40 100644 --- a/src/tactic/smtlogics/qflia_tactic.cpp +++ b/src/tactic/smtlogics/qflia_tactic.cpp @@ -228,18 +228,6 @@ tactic * mk_qflia_tactic(ast_manager & m, params_ref const & p) { #endif main_p); - // - - - // tactic * st = using_params(and_then(preamble_st, - // or_else(mk_ilp_model_finder_tactic(m), - // mk_pb_tactic(m), - // and_then(fail_if_not(mk_is_quasi_pb_probe()), - // using_params(mk_lia2sat_tactic(m), quasi_pb_p), - // mk_fail_if_undecided_tactic()), - // mk_bounded_tactic(m), - // mk_smt_tactic())), - // main_p); st->updt_params(p); return st; diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 140aebf3d2..e762398a86 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -44,10 +44,8 @@ #include "util/lp/numeric_pair.h" #include "util/lp/binary_heap_upair_queue.h" #include "util/lp/stacked_value.h" -#include "util/lp/stacked_unordered_set.h" #include "util/lp/int_set.h" #include "util/stopwatch.h" -#include "util/lp/stacked_map.h" #include #include "test/lp/gomory_test.h" #include "util/lp/matrix.h" @@ -55,9 +53,16 @@ #include "util/lp/square_sparse_matrix_def.h" #include "util/lp/lu_def.h" #include "util/lp/general_matrix.h" +#include "util/lp/bound_propagator.h" namespace lp { unsigned seed = 1; + class my_bound_propagator : public bound_propagator { + public: + my_bound_propagator(lar_solver & ls): bound_propagator(ls) {} + void consume(mpq const& v, lp::constraint_index j) override {} + }; + random_gen g_rand; static unsigned my_random() { return g_rand(); @@ -1942,44 +1947,6 @@ void setup_args_parser(argument_parser & parser) { struct fff { int a; int b;}; -void test_stacked_map_itself() { - vector v(3,0); - for(auto u : v) - std::cout << u << std::endl; - - std::unordered_map foo; - fff l; - l.a = 0; - l.b =1; - foo[1] = l; - int r = 1; - int k = foo[r].a; - std::cout << k << std::endl; - - stacked_map m; - m[0] = 3; - m[1] = 4; - m.push(); - m[1] = 5; - m[2] = 2; - m.pop(); - m.erase(2); - m[2] = 3; - m.erase(1); - m.push(); - m[3] = 100; - m[4] = 200; - m.erase(1); - m.push(); - m[5] = 300; - m[6] = 400; - m[5] = 301; - m.erase(5); - m[3] = 122; - - m.pop(2); - m.pop(); -} void test_stacked_unsigned() { std::cout << "test stacked unsigned" << std::endl; @@ -2038,36 +2005,10 @@ void test_stacked_vector() { } -void test_stacked_set() { -#ifdef Z3DEBUG - std::cout << "test_stacked_set" << std::endl; - stacked_unordered_set s; - s.insert(1); - s.insert(2); - s.insert(3); - std::unordered_set scopy = s(); - s.push(); - s.insert(4); - s.pop(); - lp_assert(s() == scopy); - s.push(); - s.push(); - s.insert(4); - s.insert(5); - s.push(); - s.insert(4); - s.pop(3); - lp_assert(s() == scopy); -#endif -} void test_stacked() { - std::cout << "test_stacked_map()" << std::endl; - test_stacked_map_itself(); test_stacked_value(); test_stacked_vector(); - test_stacked_set(); - } char * find_home_dir() { @@ -2815,7 +2756,7 @@ void test_bound_propagation_one_small_sample1() { vector ev; ls.add_var_bound(a, LE, mpq(1)); ls.solve(); - bound_propagator bp(ls); + my_bound_propagator bp(ls); ls.propagate_bounds_for_touched_rows(bp); std::cout << " bound ev from test_bound_propagation_one_small_sample1" << std::endl; for (auto & be : bp.m_ibounds) { @@ -2868,7 +2809,7 @@ void test_bound_propagation_one_row() { vector ev; ls.add_var_bound(x0, LE, mpq(1)); ls.solve(); - bound_propagator bp(ls); + my_bound_propagator bp(ls); ls.propagate_bounds_for_touched_rows(bp); } void test_bound_propagation_one_row_with_bounded_vars() { @@ -2884,7 +2825,7 @@ void test_bound_propagation_one_row_with_bounded_vars() { ls.add_var_bound(x0, LE, mpq(3)); ls.add_var_bound(x0, LE, mpq(1)); ls.solve(); - bound_propagator bp(ls); + my_bound_propagator bp(ls); ls.propagate_bounds_for_touched_rows(bp); } void test_bound_propagation_one_row_mixed() { @@ -2898,7 +2839,7 @@ void test_bound_propagation_one_row_mixed() { vector ev; ls.add_var_bound(x1, LE, mpq(1)); ls.solve(); - bound_propagator bp(ls); + my_bound_propagator bp(ls); ls.propagate_bounds_for_touched_rows(bp); } @@ -2921,7 +2862,7 @@ void test_bound_propagation_two_rows() { vector ev; ls.add_var_bound(y, LE, mpq(1)); ls.solve(); - bound_propagator bp(ls); + my_bound_propagator bp(ls); ls.propagate_bounds_for_touched_rows(bp); } @@ -2941,7 +2882,7 @@ void test_total_case_u() { vector ev; ls.add_var_bound(z, GE, zero_of_type()); ls.solve(); - bound_propagator bp(ls); + my_bound_propagator bp(ls); ls.propagate_bounds_for_touched_rows(bp); } bool contains_j_kind(unsigned j, lconstraint_kind kind, const mpq & rs, const vector & ev) { @@ -2968,7 +2909,7 @@ void test_total_case_l(){ vector ev; ls.add_var_bound(z, LE, zero_of_type()); ls.solve(); - bound_propagator bp(ls); + my_bound_propagator bp(ls); ls.propagate_bounds_for_touched_rows(bp); lp_assert(ev.size() == 4); lp_assert(contains_j_kind(x, GE, - one_of_type(), ev)); diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp index 6620dc83bc..782be7ea5e 100644 --- a/src/test/simplex.cpp +++ b/src/test/simplex.cpp @@ -4,7 +4,6 @@ Copyright (c) 2015 Microsoft Corporation --*/ -#include "util/lp/sparse_matrix.h" #include "math/simplex/sparse_matrix_def.h" #include "math/simplex/simplex.h" #include "math/simplex/simplex_def.h" diff --git a/src/util/lp/bound_propagator.h b/src/util/lp/bound_propagator.h index 61e3045f3b..a1f0301aad 100644 --- a/src/util/lp/bound_propagator.h +++ b/src/util/lp/bound_propagator.h @@ -22,6 +22,6 @@ class bound_propagator { lp::lconstraint_kind kind, const rational & bval) {return true;} unsigned number_of_found_bounds() const { return m_ibounds.size(); } - virtual void consume(mpq const& v, unsigned j) { std::cout << "doh\n"; } + virtual void consume(mpq const& v, lp::constraint_index j) = 0; }; } diff --git a/src/util/lp/explanation.h b/src/util/lp/explanation.h index c3d876afe7..f811e1d0ec 100644 --- a/src/util/lp/explanation.h +++ b/src/util/lp/explanation.h @@ -20,6 +20,7 @@ Revision History: #pragma once namespace lp { struct explanation { + void clear() { m_explanation.clear(); } vector> m_explanation; void push_justification(constraint_index j, const mpq& v) { m_explanation.push_back(std::make_pair(v, j)); diff --git a/src/util/lp/hnf.h b/src/util/lp/hnf.h index 7a1b448c7f..e0d7b0f20d 100644 --- a/src/util/lp/hnf.h +++ b/src/util/lp/hnf.h @@ -32,7 +32,7 @@ Revision History: namespace lp { namespace hnf_calc { - // d = u * a + b * b and the sum of abs(u) + abs(v) is minimal, d is positive + // d = u * a + v * b and the sum of abs(u) + abs(v) is minimal, d is positive inline void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq & v) { if (is_zero(a)) { @@ -47,7 +47,11 @@ void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq d = a; return; } +#if 1 + d = gcd(a, b, u, v); +#else extended_gcd(a, b, d, u, v); +#endif if (is_neg(d)) { d = -d; u = -u; @@ -103,7 +107,8 @@ void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq -template bool prepare_pivot_for_lower_triangle(M &m, unsigned r) { +template +bool prepare_pivot_for_lower_triangle(M &m, unsigned r) { for (unsigned i = r; i < m.row_count(); i++) { for (unsigned j = r; j < m.column_count(); j++) { if (!is_zero(m[i][j])) { @@ -120,13 +125,13 @@ template bool prepare_pivot_for_lower_triangle(M &m, unsigned r) { return false; } -template void pivot_column_non_fractional(M &m, unsigned r) { +template +void pivot_column_non_fractional(M &m, unsigned r) { lp_assert(!is_zero(m[r][r])); - lp_assert(m.row_count() <= m.column_count()); for (unsigned j = r + 1; j < m.column_count(); j++) { for (unsigned i = r + 1; i < m.row_count(); i++) { m[i][j] = - (r > 0 )? + (r > 0) ? (m[r][r]*m[i][j] - m[i][r]*m[r][j]) / m[r-1][r-1] : (m[r][r]*m[i][j] - m[i][r]*m[r][j]); lp_assert(is_int(m[i][j])); @@ -140,8 +145,8 @@ template void pivot_column_non_fractional(M &m, unsigned r) { } // returns the rank of the matrix -template unsigned to_lower_triangle_non_fractional(M &m) { - lp_assert(m.row_count() <= m.column_count()); +template +unsigned to_lower_triangle_non_fractional(M &m) { unsigned i = 0; for (; i < m.row_count(); i++) { if (!prepare_pivot_for_lower_triangle(m, i)) { @@ -152,6 +157,8 @@ template unsigned to_lower_triangle_non_fractional(M &m) { lp_assert(i == m.row_count()); return i; } + +// returns gcd of values below diagonal i,i template mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) { mpq g = zero_of_type(); @@ -170,28 +177,28 @@ mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) { return g; } - -// it fills "r" - the basic rows of m -template mpq determinant_of_rectangular_matrix(const M& m, svector & basis_rows) { - if (m.column_count() < m.row_count()) - throw "not implemented"; // consider working with the transposed m, or create a "transposed" version if needed - // The plan is to transform m to the lower triangular form by using non-fractional Gaussian Elimination by columns. - // Then the elements of the following elements of the last non-zero row of the matrix - // m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r. - // The gcd of these minors is the return value - auto mc = m; - unsigned rank = to_lower_triangle_non_fractional(mc); +// It fills "r" - the basic rows of m. +// The plan is to transform m to the lower triangular form by using non-fractional Gaussian Elimination by columns. +// Then the trailing after the diagonal elements of the following elements of the last non-zero row of the matrix, +// namely, m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r. +// The gcd of these minors is the return value. + +template +mpq determinant_of_rectangular_matrix(const M& m, svector & basis_rows) { + auto m_copy = m; + unsigned rank = to_lower_triangle_non_fractional(m_copy); if (rank == 0) return one_of_type(); for (unsigned i = 0; i < rank; i++) { - basis_rows.push_back(mc.adjust_row(i)); + basis_rows.push_back(m_copy.adjust_row(i)); } - TRACE("hnf_calc", tout << "basis_rows = "; print_vector(basis_rows, tout); mc.print(tout, "mc = ");); - return gcd_of_row_starting_from_diagonal(mc, rank - 1); + TRACE("hnf_calc", tout << "basis_rows = "; print_vector(basis_rows, tout); m_copy.print(tout, "m_copy = ");); + return gcd_of_row_starting_from_diagonal(m_copy, rank - 1); } -template mpq determinant(const M& m) { +template +mpq determinant(const M& m) { lp_assert(m.row_count() == m.column_count()); auto mc = m; svector basis_rows; @@ -229,10 +236,8 @@ class hnf { mpq mod_R(const mpq & a) const { mpq t = a % m_R; - t = is_neg(t) ? t + m_R : t; - if(is_neg(t)){ - std::cout << "a=" << a << ", m_R= " << m_R << std::endl; - } + t = is_neg(t) ? t + m_R : t; + CTRACE("hnf", is_neg(t), tout << "a=" << a << ", m_R= " << m_R << std::endl;); return t; } @@ -319,8 +324,8 @@ class hnf { } } - void handle_column_ij_in_row_i(unsigned i, unsigned j) { - lp_assert(is_correct_modulo()); + void handle_column_ij_in_row_i(unsigned i, unsigned j) { + lp_assert(is_correct_modulo()); const mpq& aii = m_H[i][i]; const mpq& aij = m_H[i][j]; mpq p,q,r; @@ -470,15 +475,14 @@ class hnf { bool is_correct_final() const { if (!is_correct()) { - std::cout << "m_H = "; m_H.print(std::cout, 17); - - std::cout << "\nm_A_orig * m_U = "; (m_A_orig * m_U).print(std::cout, 17); - - std::cout << "is_correct() does not hold" << std::endl; + TRACE("hnf_calc", + tout << "m_H = "; m_H.print(tout, 17); + tout << "\nm_A_orig * m_U = "; (m_A_orig * m_U).print(tout, 17); + tout << "is_correct() does not hold" << std::endl;); return false; } if (!is_correct_form()) { - std::cout << "is_correct_form() does not hold" << std::endl; + TRACE("hnf_calc", tout << "is_correct_form() does not hold" << std::endl;); return false; } return true; @@ -538,10 +542,6 @@ class hnf { copy_buffer_to_col_i_W_modulo(); } - void endl() const { - std::cout << std::endl; - } - void fix_row_under_diagonal_W_modulo() { mpq d, u, v; if (is_zero(m_W[m_i][m_i])) { @@ -616,12 +616,11 @@ class hnf { #endif calculate_by_modulo(); #ifdef Z3DEBUG - if (m_H != m_W) { - std::cout << "A = "; m_A_orig.print(std::cout, 4); endl(); - std::cout << "H = "; m_H.print(std::cout, 4); endl(); - std::cout << "W = "; m_W.print(std::cout, 4); endl(); - lp_assert(false); - } + CTRACE("hnf_calc", m_H != m_W, + tout << "A = "; m_A_orig.print(tout, 4); tout << std::endl; + tout << "H = "; m_H.print(tout, 4); tout << std::endl; + tout << "W = "; m_W.print(tout, 4); tout << std::endl;); + lp_assert (m_H == m_W); #endif } diff --git a/src/util/lp/hnf_cutter.h b/src/util/lp/hnf_cutter.h index c2a8672f71..1035183d7d 100644 --- a/src/util/lp/hnf_cutter.h +++ b/src/util/lp/hnf_cutter.h @@ -29,10 +29,12 @@ class hnf_cutter { var_register m_var_register; general_matrix m_A; vector m_terms; + svector m_constraints_for_explanation; vector m_right_sides; unsigned m_row_count; unsigned m_column_count; lp_settings & m_settings; + public: hnf_cutter(lp_settings & settings) : m_row_count(0), m_column_count(0), m_settings(settings) {} @@ -41,25 +43,29 @@ class hnf_cutter { } const vector& terms() const { return m_terms; } + const svector& constraints_for_explanation() const { + return m_constraints_for_explanation; + } const vector & right_sides() const { return m_right_sides; } void clear() { // m_A will be filled from scratch in init_matrix_A m_var_register.clear(); m_terms.clear(); - m_right_sides.clear(); + m_constraints_for_explanation.clear(); + m_right_sides.clear(); m_row_count = m_column_count = 0; } - void add_term(const lar_term* t, const mpq &rs) { + void add_term(const lar_term* t, const mpq &rs, constraint_index ci) { m_terms.push_back(t); + m_right_sides.push_back(rs); + m_constraints_for_explanation.push_back(ci); for (const auto &p : *t) { m_var_register.add_var(p.var()); } - m_right_sides.push_back(rs); if (m_terms.size() <= m_var_register.size()) { // capture the maximal acceptable sizes m_row_count = m_terms.size(); m_column_count = m_var_register.size(); } - } void print(std::ostream & out) { out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl; @@ -93,8 +99,9 @@ class hnf_cutter { if (basis_rows.size() == m_right_sides.size()) return m_right_sides; vector b; - for (unsigned i : basis_rows) + for (unsigned i : basis_rows) { b.push_back(m_right_sides[i]); + } return b; } @@ -161,19 +168,30 @@ class hnf_cutter { return ret; } #endif + void shrink_explanation(const svector& basis_rows) { + svector new_expl; + for (unsigned i : basis_rows) { + new_expl.push_back(m_constraints_for_explanation[i]); + } + m_constraints_for_explanation = new_expl; + } + lia_move create_cut(lar_term& t, mpq& k, explanation& ex, bool & upper #ifdef Z3DEBUG , const vector & x0 #endif ) { + // we suppose that x0 has at least one non integer element init_matrix_A(); svector basis_rows; mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows); if (m_settings.get_cancel_flag()) return lia_move::undef; - if (basis_rows.size() < m_A.row_count()) + if (basis_rows.size() < m_A.row_count()) { m_A.shrink_to_rank(basis_rows); + shrink_explanation(basis_rows); + } hnf h(m_A, d); // general_matrix A_orig = m_A; @@ -184,22 +202,10 @@ class hnf_cutter { find_h_minus_1_b(h.W(), b); // lp_assert(bcopy == h.W().take_first_n_columns(b.size()) * b); int cut_row = find_cut_row_index(b); - if (cut_row == -1) { + if (cut_row == -1) return lia_move::undef; - } - // test region - /* - general_matrix U(m_A.column_count()); - vector rt(m_A.column_count()); - for (unsigned i = 0; i < U.row_count(); i++) { - get_ei_H_minus_1(i, h.W(), rt); - vector ft = rt * A_orig; - for (unsigned j = 0; j < ft.size(); j++) - U[i][j] = ft[j]; - } - std::cout << "U reverse = "; U.print(std::cout, 12); std::cout << std::endl; - */ - // end test region + // the matrix is not square - we can get + // all integers in b's projection vector row(m_A.column_count()); get_ei_H_minus_1(cut_row, h.W(), row); @@ -209,5 +215,7 @@ class hnf_cutter { upper = true; return lia_move::cut; } + + svector vars() const { return m_var_register.vars(); } }; } diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index 73fdacbaf7..020b241717 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -409,6 +409,8 @@ struct pivoted_rows_tracking_control { } }; + + impq int_solver::get_cube_delta_for_term(const lar_term& t) const { if (t.size() == 2) { bool seen_minus = false; @@ -528,34 +530,37 @@ lia_move int_solver::gomory_cut() { } -bool int_solver::try_add_term_to_A_for_hnf(unsigned i) { +void int_solver::try_add_term_to_A_for_hnf(unsigned i) { mpq rs; const lar_term* t = m_lar_solver->terms()[i]; +#if Z3DEBUG for (const auto & p : *t) { if (!is_int(p.var())) { lp_assert(false); - return false; // todo : the mix case! } } - bool has_bounds; - if (m_lar_solver->get_equality_and_right_side_for_term_on_corrent_x(i, rs, has_bounds)) { - m_hnf_cutter.add_term(t, rs); - return true; +#endif + constraint_index ci; + if (m_lar_solver->get_equality_and_right_side_for_term_on_current_x(i, rs, ci)) { + m_hnf_cutter.add_term(t, rs, ci); } - return !has_bounds; } -bool int_solver::hnf_matrix_is_empty() const { return true; } +bool int_solver::hnf_has_non_integral_var() const { + for (unsigned j : m_hnf_cutter.vars()) + if (get_value(j).is_int() == false) + return true; + return false; +} bool int_solver::init_terms_for_hnf_cut() { m_hnf_cutter.clear(); - for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) { + for (unsigned i = 0; i < m_lar_solver->terms().size() && m_hnf_cutter.row_count() < settings().limit_on_rows_for_hnf_cutter; i++) { try_add_term_to_A_for_hnf(i); } - return m_hnf_cutter.row_count() < settings().limit_on_rows_for_hnf_cutter; + return hnf_has_non_integral_var(); } - lia_move int_solver::make_hnf_cut() { if (!init_terms_for_hnf_cut()) { return lia_move::undef; @@ -574,6 +579,10 @@ lia_move int_solver::make_hnf_cut() { if (r == lia_move::cut) { lp_assert(current_solution_is_inf_on_cut()); settings().st().m_hnf_cuts++; + m_ex->clear(); + for (unsigned i : m_hnf_cutter.constraints_for_explanation()) { + m_ex->push_justification(i); + } } return r; } @@ -599,7 +608,6 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { r = patch_nbasic_columns(); if (r != lia_move::undef) return r; - ++m_branch_cut_counter; r = find_cube(); diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index c2aa323606..0bc7d5ef05 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -51,7 +51,9 @@ class int_solver { // or provide a way of how it can be adjusted. lia_move check(lar_term& t, mpq& k, explanation& ex, bool & upper); bool move_non_basic_column_to_bounds(unsigned j); - lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex); + lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex); + bool is_base(unsigned j) const; + private: // how to tighten bounds for integer variables. @@ -75,14 +77,12 @@ class int_solver { mpq const & consts); void fill_explanation_from_fixed_columns(const row_strip & row); void add_to_explanation_from_fixed_or_boxed_column(unsigned j); - void patch_nbasic_column(unsigned j); lia_move patch_nbasic_columns(); bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m); const impq & lower_bound(unsigned j) const; const impq & upper_bound(unsigned j) const; bool is_int(unsigned j) const; bool is_real(unsigned j) const; - bool is_base(unsigned j) const; bool is_boxed(unsigned j) const; bool is_fixed(unsigned j) const; bool is_free(unsigned j) const; @@ -156,6 +156,8 @@ class int_solver { lia_move make_hnf_cut(); bool init_terms_for_hnf_cut(); bool hnf_matrix_is_empty() const; - bool try_add_term_to_A_for_hnf(unsigned term_index); + void try_add_term_to_A_for_hnf(unsigned term_index); + bool hnf_has_non_integral_var() const; + void patch_nbasic_column(unsigned j); }; } diff --git a/src/util/lp/lar_core_solver.h b/src/util/lp/lar_core_solver.h index 24f1ddad3a..5f55185287 100644 --- a/src/util/lp/lar_core_solver.h +++ b/src/util/lp/lar_core_solver.h @@ -26,7 +26,6 @@ Revision History: #include "util/lp/indexed_vector.h" #include "util/lp/binary_heap_priority_queue.h" #include "util/lp/breakpoint.h" -#include "util/lp/stacked_unordered_set.h" #include "util/lp/lp_primal_core_solver.h" #include "util/lp/stacked_vector.h" #include "util/lp/lar_solution_signature.h" @@ -611,7 +610,6 @@ class lar_core_solver { } if (no_r_lu()) { // it is the case where m_d_solver gives a degenerated basis, we need to roll back - // std::cout << "no_r_lu" << std::endl; catch_up_in_lu_in_reverse(changes_of_basis, m_r_solver); m_r_solver.find_feasible_solution(); m_d_basis = m_r_basis; diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index d6d5563033..def7edb154 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -370,7 +370,6 @@ void lar_solver::pop(unsigned k) { unsigned m = A_r().row_count(); clean_popped_elements(m, m_rows_with_changed_bounds); clean_inf_set_of_r_solver_after_pop(); - m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN; lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided || (!use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); @@ -508,14 +507,45 @@ bool lar_solver::maximize_term_on_corrected_r_solver(const vector> & term, +lp_status lar_solver::maximize_term(const vector> & term, impq &term_max) { lp_assert(m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()); m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false; - return maximize_term_on_corrected_r_solver(term, term_max); + if (!maximize_term_on_corrected_r_solver(term, term_max)) + return lp_status::UNBOUNDED; + + bool change = false; + for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_x.size(); j++) { + if (!column_is_int(j)) + continue; + if (column_value_is_integer(j)) + continue; + if (m_int_solver->is_base(j)) { + if (!remove_from_basis(j)) // consider a special version of remove_from_basis that would not remove inf_int columns + return lp_status::FEASIBLE; // it should not happen + } + m_int_solver->patch_nbasic_column(j); + if (!column_value_is_integer(j)) + return lp_status::FEASIBLE; + change = true; + } + if (change) { + term_max = zero_of_type(); + for (const auto& p : term) { + term_max += p.first * m_mpq_lar_core_solver.m_r_x[p.second]; + } + } + return lp_status::OPTIMAL; } @@ -792,18 +822,11 @@ void lar_solver::add_last_rows_to_lu(lp_primal_core_solver & s) { bool lar_solver::x_is_correct() const { if (m_mpq_lar_core_solver.m_r_x.size() != A_r().column_count()) { - // std::cout << "the size is off " << m_r_solver.m_x.size() << ", " << A().column_count() << std::endl; return false; } for (unsigned i = 0; i < A_r().row_count(); i++) { numeric_pair delta = A_r().dot_product_with_row(i, m_mpq_lar_core_solver.m_r_x); if (!delta.is_zero()) { - // std::cout << "x is off ("; - // std::cout << "m_b[" << i << "] = " << m_b[i] << " "; - // std::cout << "left side = " << A().dot_product_with_row(i, m_r_solver.m_x) << ' '; - // std::cout << "delta = " << delta << ' '; - // std::cout << "iters = " << total_iterations() << ")" << std::endl; - // std::cout << "row " << i << " is off" << std::endl; return false; } } @@ -906,7 +929,6 @@ bool lar_solver::all_constraints_hold() const { for (unsigned i = 0; i < m_constraints.size(); i++) { if (!constraint_holds(*m_constraints[i], var_map)) { - print_constraint(i, std::cout); return false; } } @@ -973,13 +995,6 @@ bool lar_solver::the_left_sides_sum_to_zero(const vector= m_terms_start_index) - throw 0; // todo : what is the right way to exit? auto it = m_ext_vars_to_columns.find(ext_j); if (it != m_ext_vars_to_columns.end()) { return it->second.internal_j(); @@ -1580,7 +1594,6 @@ bool lar_solver::term_coeffs_are_ok(const vector> & co bool g_is_set = false; for (const auto & p : coeffs) { if (!p.first.is_int()) { - std::cout << "p.first = " << p.first << " is not an int\n"; return false; } if (!g_is_set) { @@ -1593,9 +1606,6 @@ bool lar_solver::term_coeffs_are_ok(const vector> & co if (g == one_of_type()) return true; - std::cout << "term is not ok: g = " << g << std::endl; - this->print_linear_combination_of_column_indices_only(coeffs, std::cout); - std::cout << " rs = " << v << std::endl; return false; } #endif @@ -2207,36 +2217,34 @@ void lar_solver::round_to_integer_solution() { } } -bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & rs, bool & has_bounds) const { +bool lar_solver::get_equality_and_right_side_for_term_on_current_x(unsigned term_index, mpq & rs, constraint_index& ci) const { unsigned tj = term_index + m_terms_start_index; auto it = m_ext_vars_to_columns.find(tj); - has_bounds = false; if (it == m_ext_vars_to_columns.end()) return false; unsigned j = it->second.internal_j(); - auto & slv = m_mpq_lar_core_solver.m_r_solver; impq term_val; bool term_val_ready = false; - if (slv.column_has_upper_bound(j)) { - has_bounds = true; - const impq & b = slv.m_upper_bounds[j]; - lp_assert(is_zero(b.y) && is_int(b.x)); + mpq b; + bool is_strict; + if (has_upper_bound(j, ci, b, is_strict)) { + lp_assert(!is_strict); + lp_assert(b.is_int()); term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); term_val_ready = true; - if (term_val == b) { - rs = b.x; + if (term_val.x == b) { + rs = b; return true; } } - if (slv.column_has_lower_bound(j)) { - has_bounds = true; + if (has_lower_bound(j, ci, b, is_strict)) { + lp_assert(!is_strict); if (!term_val_ready) term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); - const impq & b = slv.m_lower_bounds[j]; - lp_assert(is_zero(b.y) && is_int(b.x)); + lp_assert(b.is_int()); - if (term_val == b) { - rs = b.x; + if (term_val.x == b) { + rs = b; return true; } } diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index c7b62ee1de..ce975738f7 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -36,7 +36,6 @@ Revision History: #include #include "util/lp/stacked_value.h" #include "util/lp/stacked_vector.h" -#include "util/lp/stacked_unordered_set.h" #include "util/lp/implied_bound.h" #include "util/lp/bound_analyzer_on_row.h" #include "util/lp/conversion_helper.h" @@ -326,8 +325,8 @@ public : impq &term_max); // starting from a given feasible state look for the maximum of the term // return true if found and false if unbounded - bool maximize_term(const vector> & term, - impq &term_max); + lp_status maximize_term(const vector> & term, + impq &term_max); @@ -435,9 +434,9 @@ public : mpq sum_of_right_sides_of_explanation(const vector> & explanation) const; - bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict); + bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const; - bool has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict); + bool has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const; void get_infeasibility_explanation(vector> & explanation) const; @@ -587,6 +586,7 @@ public : unsigned column_count() const { return A_r().column_count(); } const vector & r_basis() const { return m_mpq_lar_core_solver.r_basis(); } const vector & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); } - bool get_equality_and_right_side_for_term_on_corrent_x(unsigned i, mpq &rs, bool & has_bounds) const; + bool get_equality_and_right_side_for_term_on_current_x(unsigned i, mpq &rs, constraint_index& ci) const; + bool remove_from_basis(unsigned); }; } diff --git a/src/util/lp/lp_core_solver_base.cpp b/src/util/lp/lp_core_solver_base.cpp index 613404e68f..00c1322c2f 100644 --- a/src/util/lp/lp_core_solver_base.cpp +++ b/src/util/lp/lp_core_solver_base.cpp @@ -145,3 +145,4 @@ template bool lp::lp_core_solver_base >::infe template bool lp::lp_core_solver_base::infeasibility_costs_are_correct() const; template bool lp::lp_core_solver_base::infeasibility_costs_are_correct() const; template void lp::lp_core_solver_base >::calculate_pivot_row(unsigned int); +template bool lp::lp_core_solver_base >::remove_from_basis(unsigned int); diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index c43ca7142f..b8b7b7c0dc 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -225,7 +225,6 @@ class lp_core_solver_base { lp_assert(m_A.is_correct()); if (m_using_infeas_costs) { if (infeasibility_costs_are_correct() == false) { - std::cout << "infeasibility_costs_are_correct() does not hold" << std::endl; return false; } } @@ -234,9 +233,6 @@ class lp_core_solver_base { for (unsigned j = 0; j < n; j++) { if (m_basis_heading[j] >= 0) { if (!is_zero(m_d[j])) { - - std::cout << "case a\n"; - print_column_info(j, std::cout); return false; } } else { @@ -245,8 +241,6 @@ class lp_core_solver_base { d -= this->m_costs[this->m_basis[cc.m_i]] * this->m_A.get_val(cc); } if (m_d[j] != d) { - std::cout << "case b\n"; - print_column_info(j, std::cout); return false; } } @@ -453,6 +447,7 @@ class lp_core_solver_base { void init_lu(); int pivots_in_column_and_row_are_different(int entering, int leaving) const; void pivot_fixed_vars_from_basis(); + bool remove_from_basis(unsigned j); bool pivot_column_general(unsigned j, unsigned j_basic, indexed_vector & w); bool pivot_for_tableau_on_basis(); bool pivot_row_for_tableau_on_basis(unsigned row); @@ -552,7 +547,6 @@ class lp_core_solver_base { bool non_basic_columns_are_set_correctly() const { for (unsigned j : this->m_nbasis) if (!column_is_feasible(j)) { - print_column_info(j, std::cout); return false; } return true; @@ -601,10 +595,6 @@ class lp_core_solver_base { out << " base\n"; else out << " nbas\n"; - - /* - std::cout << "cost = " << m_costs[j] << std::endl; - std:: cout << "m_d = " << m_d[j] << std::endl;*/ } bool column_is_free(unsigned j) const { return this->m_column_type[j] == free; } diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index 4577b26ecd..dd3de910f9 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -228,11 +228,6 @@ A_mult_x_is_off() const { for (unsigned i = 0; i < m_m(); i++) { X delta = m_b[i] - m_A.dot_product_with_row(i, m_x); if (delta != numeric_traits::zero()) { - std::cout << "precise x is off ("; - std::cout << "m_b[" << i << "] = " << m_b[i] << " "; - std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' '; - std::cout << "delta = " << delta << ' '; - std::cout << "iters = " << total_iterations() << ")" << std::endl; return true; } } @@ -265,11 +260,6 @@ A_mult_x_is_off_on_index(const vector & index) const { for (unsigned i : index) { X delta = m_b[i] - m_A.dot_product_with_row(i, m_x); if (delta != numeric_traits::zero()) { - // std::cout << "x is off ("; - // std::cout << "m_b[" << i << "] = " << m_b[i] << " "; - // std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' '; - // std::cout << "delta = " << delta << ' '; - // std::cout << "iters = " << total_iterations() << ")" << std::endl; return true; } } @@ -416,13 +406,10 @@ column_is_dual_feasible(unsigned j) const { case column_type::lower_bound: return x_is_at_lower_bound(j) && d_is_not_negative(j); case column_type::upper_bound: - LP_OUT(m_settings, "upper_bound type should be switched to lower_bound" << std::endl); lp_assert(false); // impossible case case column_type::free_column: return numeric_traits::is_zero(m_d[j]); default: - LP_OUT(m_settings, "column = " << j << std::endl); - LP_OUT(m_settings, "unexpected column type = " << column_type_to_string(m_column_types[j]) << std::endl); lp_unreachable(); } lp_unreachable(); @@ -530,9 +517,6 @@ template bool lp_core_solver_base::inf_set_is_cor bool is_feas = column_is_feasible(j); if (is_feas == belongs_to_set) { - print_column_info(j, std::cout); - std::cout << "belongs_to_set = " << belongs_to_set << std::endl; - std::cout <<( is_feas? "feas":"inf") << std::endl; return false; } } @@ -714,22 +698,18 @@ template bool lp_core_solver_base:: lp_assert(m_basis.size() == m_A.row_count()); lp_assert(m_nbasis.size() <= m_A.column_count() - m_A.row_count()); // for the dual the size of non basis can be smaller if (!basis_has_no_doubles()) { - // std::cout << "basis_has_no_doubles" << std::endl; return false; } if (!non_basis_has_no_doubles()) { - // std::cout << "non_basis_has_no_doubles" << std::endl; return false; } if (!basis_is_correctly_represented_in_heading()) { - // std::cout << "basis_is_correctly_represented_in_heading" << std::endl; return false; } if (!non_basis_is_correctly_represented_in_heading()) { - // std::cout << "non_basis_is_correctly_represented_in_heading" << std::endl; return false; } @@ -989,6 +969,17 @@ template void lp_core_solver_base::pivot_fixed_v } } +template bool lp_core_solver_base::remove_from_basis(unsigned basic_j) { + indexed_vector w(m_basis.size()); // the buffer + unsigned i = m_basis_heading[basic_j]; + for (auto &c : m_A.m_rows[i]) { + if (pivot_column_general(c.var(), basic_j, w)) + return true; + } + return false; +} + + template bool lp_core_solver_base::infeasibility_costs_are_correct() const { if (! this->m_using_infeas_costs) @@ -996,13 +987,9 @@ lp_core_solver_base::infeasibility_costs_are_correct() const { lp_assert(costs_on_nbasis_are_zeros()); for (unsigned j :this->m_basis) { if (!infeasibility_cost_is_correct_for_column(j)) { - std::cout << "infeasibility_cost_is_correct_for_column does not hold\n"; - print_column_info(j, std::cout); return false; } if (!is_zero(m_d[j])) { - std::cout << "m_d is not zero\n"; - print_column_info(j, std::cout); return false; } } @@ -1052,9 +1039,9 @@ void lp_core_solver_base::calculate_pivot_row(unsigned i) { m_pivot_row.clear(); m_pivot_row.resize(m_n()); if (m_settings.use_tableau()) { - unsigned basis_j = m_basis[i]; + unsigned basic_j = m_basis[i]; for (auto & c : m_A.m_rows[i]) { - if (c.m_j != basis_j) + if (c.m_j != basic_j) m_pivot_row.set_value(c.get_val(), c.m_j); } return; diff --git a/src/util/lp/lp_dual_core_solver_def.h b/src/util/lp/lp_dual_core_solver_def.h index 0b47e71788..86e6231fc4 100644 --- a/src/util/lp/lp_dual_core_solver_def.h +++ b/src/util/lp/lp_dual_core_solver_def.h @@ -535,12 +535,6 @@ template bool lp_dual_core_solver::snap_runaway_n template bool lp_dual_core_solver::problem_is_dual_feasible() const { for (unsigned j : this->non_basis()){ if (!this->column_is_dual_feasible(j)) { - // std::cout << "column " << j << " is not dual feasible" << std::endl; - // std::cout << "m_d[" << j << "] = " << this->m_d[j] << std::endl; - // std::cout << "x[" << j << "] = " << this->m_x[j] << std::endl; - // std::cout << "type = " << column_type_to_string(this->m_column_type[j]) << std::endl; - // std::cout << "bounds = " << this->m_lower_bounds[j] << "," << this->m_upper_bounds[j] << std::endl; - // std::cout << "total_iterations = " << this->total_iterations() << std::endl; return false; } } diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 26fc550b3d..424823a7c6 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -504,7 +504,8 @@ class lp_primal_core_solver:public lp_core_solver_base { } X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent; advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta ); - lp_assert(this->m_x[leaving] == new_val_for_leaving); + X xleaving = this->m_x[leaving]; + lp_assert(xleaving == new_val_for_leaving); if (this->current_x_is_feasible()) this->set_status(lp_status::OPTIMAL); } diff --git a/src/util/lp/lp_primal_core_solver_tableau_def.h b/src/util/lp/lp_primal_core_solver_tableau_def.h index b1d52c6743..f85c131a32 100644 --- a/src/util/lp/lp_primal_core_solver_tableau_def.h +++ b/src/util/lp/lp_primal_core_solver_tableau_def.h @@ -322,7 +322,6 @@ template int lp_primal_core_solver::find_leaving_ return m_leaving_candidates[k]; } template void lp_primal_core_solver::init_run_tableau() { - // print_matrix(&(this->m_A), std::cout); CASSERT("A_off", this->A_mult_x_is_off() == false); lp_assert(basis_columns_are_set_correctly()); this->m_basis_sort_counter = 0; // to initiate the sort of the basis diff --git a/src/util/lp/lp_primal_simplex_def.h b/src/util/lp/lp_primal_simplex_def.h index 5366144d83..722e4bd799 100644 --- a/src/util/lp/lp_primal_simplex_def.h +++ b/src/util/lp/lp_primal_simplex_def.h @@ -278,7 +278,6 @@ template bool lp_primal_simplex::bounds_hold(std: } if (!it.second->bounds_hold(sol_it->second)) { - // std::cout << "bounds do not hold for " << it.second->get_name() << std::endl; it.second->bounds_hold(sol_it->second); return false; } diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index d8c877e8ee..553e58c627 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -214,7 +214,7 @@ struct lp_settings { lp_settings() : m_default_resource_limit(*this), m_resource_limit(&m_default_resource_limit), - m_debug_out( &std::cout), + m_debug_out(&std::cout), m_message_out(&std::cout), reps_in_scaler(20), pivot_epsilon(0.00000001), @@ -231,7 +231,6 @@ struct lp_settings { drop_tolerance ( 1e-14), tolerance_for_artificials ( 1e-4), can_be_taken_to_basis_tolerance ( 0.00001), - percent_of_entering_to_check ( 5),// we try to find a profitable column in a percentage of the columns use_scaling ( true), scaling_maximum ( 1), @@ -265,7 +264,7 @@ struct lp_settings { m_chase_cut_solver_cycle_on_var(10), m_int_pivot_fixed_vars_from_basis(false), m_int_patch_only_integer_values(true), - limit_on_rows_for_hnf_cutter(100) + limit_on_rows_for_hnf_cutter(75) {} void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; } diff --git a/src/util/lp/lp_settings_def.h b/src/util/lp/lp_settings_def.h index a35d081140..02963718e0 100644 --- a/src/util/lp/lp_settings_def.h +++ b/src/util/lp/lp_settings_def.h @@ -74,14 +74,12 @@ bool vectors_are_equal(T * a, vector &b, unsigned n) { if (numeric_traits::precise()) { for (unsigned i = 0; i < n; i ++){ if (!numeric_traits::is_zero(a[i] - b[i])) { - // std::cout << "a[" << i <<"]" << a[i] << ", " << "b[" << i <<"]" << b[i] << std::endl; return false; } } } else { for (unsigned i = 0; i < n; i ++){ if (std::abs(numeric_traits::get_double(a[i] - b[i])) > 0.000001) { - // std::cout << "a[" << i <<"]" << a[i] << ", " << "b[" << i <<"]" << b[i] << std::endl; return false; } } @@ -97,7 +95,6 @@ bool vectors_are_equal(const vector & a, const vector &b) { if (numeric_traits::precise()) { for (unsigned i = 0; i < n; i ++){ if (!numeric_traits::is_zero(a[i] - b[i])) { - // std::cout << "a[" << i <<"]" << a[i] << ", " << "b[" << i <<"]" << b[i] << std::endl; return false; } } @@ -112,7 +109,6 @@ bool vectors_are_equal(const vector & a, const vector &b) { } if (fabs(da - db) > 0.000001) { - // std::cout << "a[" << i <<"] = " << a[i] << ", but " << "b[" << i <<"] = " << b[i] << std::endl; return false; } } diff --git a/src/util/lp/lu_def.h b/src/util/lp/lu_def.h index 89c959e4d4..be4cd724d3 100644 --- a/src/util/lp/lu_def.h +++ b/src/util/lp/lu_def.h @@ -144,8 +144,6 @@ lu::lu(const M& A, m_row_eta_work_vector(A.row_count()), m_refactor_counter(0) { lp_assert(A.row_count() == A.column_count()); - std::cout << "m_U = \n"; - print_matrix(&m_U, std::cout); create_initial_factorization(); #ifdef Z3DEBUG lp_assert(is_correct()); @@ -659,13 +657,10 @@ template bool lu::is_correct() { #ifdef Z3DEBUG if (get_status() != LU_status::OK) { - std::cout << " status" << std::endl; return false; } dense_matrix left_side = get_left_side(); - std::cout << "ls = "; print_matrix(left_side, std::cout); dense_matrix right_side = get_right_side(); - std::cout << "rs = "; print_matrix(right_side, std::cout); return left_side == right_side; #else return true; diff --git a/src/util/lp/mps_reader.h b/src/util/lp/mps_reader.h index ca6188dd79..09762cd5ef 100644 --- a/src/util/lp/mps_reader.h +++ b/src/util/lp/mps_reader.h @@ -303,7 +303,6 @@ class mps_reader { do { read_line(); if (m_line.find("RHS") == 0) { - // cout << "found RHS" << std::endl; break; } if (m_line.size() < 22) { diff --git a/src/util/lp/permutation_matrix_def.h b/src/util/lp/permutation_matrix_def.h index bb2ca6c2a1..76af12ed87 100644 --- a/src/util/lp/permutation_matrix_def.h +++ b/src/util/lp/permutation_matrix_def.h @@ -64,7 +64,6 @@ void permutation_matrix::apply_from_left(vector & w, lp_settings & ) { // L * deb_w = clone_vector(w, row_count()); // deb.apply_from_left(deb_w); #endif - // std::cout << " apply_from_left " << std::endl; lp_assert(m_X_buffer.size() == w.size()); unsigned i = size(); while (i-- > 0) { diff --git a/src/util/lp/scaler_def.h b/src/util/lp/scaler_def.h index 3f21f19539..2710f89bf3 100644 --- a/src/util/lp/scaler_def.h +++ b/src/util/lp/scaler_def.h @@ -214,7 +214,6 @@ template void scaler::scale_rows() { } template void scaler::scale_row(unsigned i) { - std::cout << "t" << "\n"; T row_max = std::max(m_A.get_max_abs_in_row(i), abs(convert_struct::convert(m_b[i]))); T alpha = numeric_traits::one(); if (numeric_traits::is_zero(row_max)) { @@ -244,7 +243,6 @@ template void scaler::scale_column(unsigned i) if (numeric_traits::is_zero(column_max)){ return; // the column has zeros only } - std::cout << "f"; if (column_max < m_scaling_minimum) { do { alpha *= 2; diff --git a/src/util/lp/square_sparse_matrix.h b/src/util/lp/square_sparse_matrix.h index 7c32064006..d84a8a289e 100644 --- a/src/util/lp/square_sparse_matrix.h +++ b/src/util/lp/square_sparse_matrix.h @@ -254,11 +254,7 @@ class square_sparse_matrix void delete_column(int i); void swap_columns(unsigned a, unsigned b) { - // cout << "swaapoiiin" << std::endl; - // dense_matrix d(*this); m_column_permutation.transpose_from_left(a, b); - // d.swap_columns(a, b); - // lp_assert(*this == d); } void swap_rows(unsigned a, unsigned b) { diff --git a/src/util/lp/square_sparse_matrix_def.h b/src/util/lp/square_sparse_matrix_def.h index 4e9600ffc2..791bdb6ae5 100644 --- a/src/util/lp/square_sparse_matrix_def.h +++ b/src/util/lp/square_sparse_matrix_def.h @@ -1214,7 +1214,6 @@ bool square_sparse_matrix::is_upper_triangular_and_maximums_are_set_correc return false; if (aj == ai) { if (iv.m_value != 1) { - // std::cout << "value at diagonal = " << iv.m_value << std::endl; return false; } } @@ -1245,18 +1244,14 @@ void square_sparse_matrix::check_column_vs_rows(unsigned col) { for (indexed_value & column_iv : mc) { indexed_value & row_iv = column_iv_other(column_iv); if (row_iv.m_index != col) { - // std::cout << "m_other in row does not belong to column " << col << ", but to column " << row_iv.m_index << std::endl; lp_assert(false); } if (& row_iv_other(row_iv) != &column_iv) { - // std::cout << "row and col do not point to each other" << std::endl; lp_assert(false); } if (row_iv.m_value != column_iv.m_value) { - // std::cout << "the data from col " << col << " for row " << column_iv.m_index << " is different in the column " << std::endl; - // std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl; lp_assert(false); } } @@ -1269,18 +1264,14 @@ void square_sparse_matrix::check_row_vs_columns(unsigned row) { indexed_value & column_iv = row_iv_other(row_iv); if (column_iv.m_index != row) { - // std::cout << "col_iv does not point to correct row " << row << " but to " << column_iv.m_index << std::endl; lp_assert(false); } if (& row_iv != & column_iv_other(column_iv)) { - // std::cout << "row and col do not point to each other" << std::endl; lp_assert(false); } if (row_iv.m_value != column_iv.m_value) { - // std::cout << "the data from col " << column_iv.m_index << " for row " << row << " is different in the column " << std::endl; - // std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl; lp_assert(false); } } diff --git a/src/util/lp/stacked_map.h b/src/util/lp/stacked_map.h deleted file mode 100644 index e1a8963fbf..0000000000 --- a/src/util/lp/stacked_map.h +++ /dev/null @@ -1,255 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ - -#pragma once -// this class implements a map with some stack functionality -#include -#include -#include -#include -#include -namespace lp { -template, - class _Alloc = std::allocator > > -class stacked_map { - struct delta { - std::unordered_set m_new; - std::unordered_map m_original_changed; - // std::unordered_map m_deb_copy; - }; - std::map m_map; - std::stack m_stack; -public: - class ref { - stacked_map & m_map; - const A & m_key; - public: - ref(stacked_map & m, const A & key) :m_map(m), m_key(key) {} - ref & operator=(const B & b) { - m_map.emplace_replace(m_key, b); - return *this; - } - ref & operator=(const ref & b) { lp_assert(false); return *this; } - operator const B&() const { - auto it = m_map.m_map.find(m_key); - lp_assert(it != m_map.m_map.end()); - return it->second; - } - }; - typename std::map < A, B, _Pr, _Alloc>::iterator upper_bound(const A& k) { - return m_map.upper_bound(k); - } - typename std::map < A, B, _Pr, _Alloc>::const_iterator upper_bound(const A& k) const { - return m_map.upper_bound(k); - } - typename std::map < A, B, _Pr, _Alloc>::iterator lower_bound(const A& k) { - return m_map.lower_bound(k); - } - typename std::map < A, B, _Pr, _Alloc>::const_iterator lower_bound(const A& k) const { - return m_map.lower_bound(k); - } - typename std::map < A, B, _Pr, _Alloc>::iterator end() { - return m_map.end(); - } - typename std::map < A, B, _Pr, _Alloc>::const_iterator end() const { - return m_map.end(); - } - typename std::map < A, B, _Pr, _Alloc>::reverse_iterator rend() { - return m_map.rend(); - } - typename std::map < A, B, _Pr, _Alloc>::const_reverse_iterator rend() const { - return m_map.rend(); - } - typename std::map < A, B, _Pr, _Alloc>::iterator begin() { - return m_map.begin(); - } - typename std::map < A, B, _Pr, _Alloc>::const_iterator begin() const { - return m_map.begin(); - } - typename std::map < A, B, _Pr, _Alloc>::reverse_iterator rbegin() { - return m_map.rbegin(); - } - typename std::map < A, B, _Pr, _Alloc>::const_reverse_iterator rbegin() const { - return m_map.rbegin(); - } -private: - void emplace_replace(const A & a, const B & b) { - if (!m_stack.empty()) { - delta & d = m_stack.top(); - auto it = m_map.find(a); - if (it == m_map.end()) { - d.m_new.insert(a); - m_map.emplace(a, b); - } else if (it->second != b) { - auto nit = d.m_new.find(a); - if (nit == d.m_new.end()) { // we do not have the old key - auto & orig_changed= d.m_original_changed; - auto itt = orig_changed.find(a); - if (itt == orig_changed.end()) { - orig_changed.emplace(a, it->second); - } else if (itt->second == b) { - orig_changed.erase(itt); - } - } - it->second = b; - } - } else { // there is no stack: just emplace or replace - m_map[a] = b; - } - } -public: - ref operator[] (const A & a) { - return ref(*this, a); - } - - const B & operator[]( const A & a) const { - auto it = m_map.find(a); - if (it == m_map.end()) { - lp_assert(false); - } - - return it->second; - } - - bool try_get_value(const A& key, B& val) const { - auto it = m_map.find(key); - if (it == m_map.end()) - return false; - - val = it->second; - return true; - } - bool try_get_value(const A&& key, B& val) const { - auto it = m_map.find(std::move(key)); - if (it == m_map.end()) - return false; - - val = it->second; - return true; - } - - unsigned size() const { - return static_cast(m_map.size()); - } - - bool contains(const A & key) const { - return m_map.find(key) != m_map.end(); - } - - bool contains(const A && key) const { - return m_map.find(std::move(key)) != m_map.end(); - } - - void push() { - delta d; - // d.m_deb_copy = m_map; - m_stack.push(d); - } - - void revert() { - if (m_stack.empty()) return; - - delta & d = m_stack.top(); - for (auto & t : d.m_new) { - m_map.erase(t); - } - for (auto & t: d.m_original_changed) { - m_map[t.first] = t.second; - } - d.clear(); - } - - void pop() { - pop(1); - } - void pop(unsigned k) { - while (k-- > 0) { - if (m_stack.empty()) - return; - delta & d = m_stack.top(); - for (auto & t : d.m_new) { - m_map.erase(t); - } - for (auto & t: d.m_original_changed) { - m_map[t.first] = t.second; - } - // lp_assert(d.m_deb_copy == m_map); - m_stack.pop(); - } - } - void erase(const typename std::map < A, B, _Pr, _Alloc>::iterator & it) { - erase(it->first); - } - void erase(const typename std::map < A, B, _Pr, _Alloc>::const_iterator & it) { - erase(it->first); - } - void erase(const typename std::map < A, B, _Pr, _Alloc>::reverse_iterator & it) { - erase(it->first); - } - void erase(const A & key) { - if (m_stack.empty()) { - m_map.erase(key); - return; - } - - delta & d = m_stack.top(); - auto it = m_map.find(key); - if (it == m_map.end()) { - lp_assert(d.m_new.find(key) == d.m_new.end()); - return; - } - auto &orig_changed = d.m_original_changed; - auto nit = d.m_new.find(key); - if (nit == d.m_new.end()) { // key is old - if (orig_changed.find(key) == orig_changed.end()) - orig_changed.emplace(it->first, it->second); // need to restore - } else { // k is new - lp_assert(orig_changed.find(key) == orig_changed.end()); - d.m_new.erase(nit); - } - - m_map.erase(it); - } - - void clear() { - if (m_stack.empty()) { - m_map.clear(); - return; - } - - delta & d = m_stack.top(); - auto & oc = d.m_original_changed; - for (auto & p : m_map) { - const auto & it = oc.find(p.first); - if (it == oc.end() && d.m_new.find(p.first) == d.m_new.end()) - oc.emplace(p.first, p.second); - } - m_map.clear(); - } - - unsigned stack_size() const { return m_stack.size(); } - - bool empty() const { return m_map.empty(); } - - const std::map& operator()() const { return m_map;} -}; -} diff --git a/src/util/lp/stacked_unordered_set.h b/src/util/lp/stacked_unordered_set.h deleted file mode 100644 index 5824aeadeb..0000000000 --- a/src/util/lp/stacked_unordered_set.h +++ /dev/null @@ -1,107 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ - -#pragma once -// this class implements an unordered_set with some stack functionality -#include -#include -#include -namespace lp { - -template , - typename KeyEqual = std::equal_to, - typename Allocator = std::allocator - > class stacked_unordered_set { - struct delta { - std::unordered_set m_inserted; - std::unordered_set m_erased; - std::unordered_set m_deb_copy; - }; - std::unordered_set m_set; - std::stack m_stack; -public: - void insert(const A & a) { - if (m_stack.empty()) { - m_set.insert(a); - } else if (m_set.find(a) == m_set.end()) { - m_set.insert(a); - size_t in_erased = m_stack.top().m_erased.erase(a); - if (in_erased == 0) { - m_stack.top().m_inserted.insert(a); - } - } - } - - void erase(const A &a) { - if (m_stack.empty()) { - m_set.erase(a); - return; - } - auto erased = m_set.erase(a); - if (erased == 1) { - auto was_new = m_stack.top().m_inserted.erase(a); - if (was_new == 0) { - m_stack.top().m_erased.insert(a); - } - } - } - - unsigned size() const { - return m_set.size(); - } - - bool contains(A & key) const { - return m_set.find(key) != m_set.end(); - } - - bool contains(A && key) const { - return m_set.find(std::move(key)) != m_set.end(); - } - - void push() { - delta d; - d.m_deb_copy = m_set; - m_stack.push(d); - } - - void pop() { - pop(1); - } - void pop(unsigned k) { - while (k-- > 0) { - if (m_stack.empty()) - return; - delta & d = m_stack.top(); - for (auto & t : d.m_inserted) { - m_set.erase(t); - } - for (auto & t : d.m_erased) { - m_set.insert(t); - } - lp_assert(d.m_deb_copy == m_set); - m_stack.pop(); - } - } - - const std::unordered_set& operator()() const { return m_set;} -}; -} - diff --git a/src/util/lp/static_matrix_def.h b/src/util/lp/static_matrix_def.h index a246a5aead..42249b4d8f 100644 --- a/src/util/lp/static_matrix_def.h +++ b/src/util/lp/static_matrix_def.h @@ -260,10 +260,6 @@ template void static_matrix::check_consistency for (auto & t : by_rows) { auto ic = by_cols.find(t.first); - if (ic == by_cols.end()){ - //std::cout << "rows have pair (" << t.first.first <<"," << t.first.second - // << "), but columns don't " << std::endl; - } lp_assert(ic != by_cols.end()); lp_assert(t.second == ic->second); } @@ -347,7 +343,6 @@ template bool static_matrix::is_correct() const { std::unordered_set s; for (auto & rc : r) { if (s.find(rc.m_j) != s.end()) { - std::cout << "found column " << rc.m_j << " twice in a row " << i << std::endl; return false; } s.insert(rc.m_j); @@ -358,7 +353,6 @@ template bool static_matrix::is_correct() const { if (rc.get_val() != get_val(m_columns[rc.m_j][rc.m_offset])) return false; if (is_zero(rc.get_val())) { - std::cout << "found zero column " << rc.m_j << " in row " << i << std::endl; return false; } @@ -370,7 +364,6 @@ template bool static_matrix::is_correct() const { std::unordered_set s; for (auto & cc : c) { if (s.find(cc.m_i) != s.end()) { - std::cout << "found row " << cc.m_i << " twice in a column " << j << std::endl; return false; } s.insert(cc.m_i); diff --git a/src/util/lp/test_bound_analyzer.h b/src/util/lp/test_bound_analyzer.h index 3765f480ea..7cee0cf56f 100644 --- a/src/util/lp/test_bound_analyzer.h +++ b/src/util/lp/test_bound_analyzer.h @@ -115,7 +115,6 @@ public : } } default: - // std::cout << " got an upper bound with " << T_to_string(l) << "\n"; m_implied_bounds.push_back(implied_bound(l, j, false, is_pos(m_coeffs[i]), m_row_or_term_index, strict)); } } @@ -207,7 +206,6 @@ public : } } default: - // std::cout << " got a lower bound with " << T_to_string(l) << "\n"; m_implied_bounds.push_back(implied_bound(l, j, true, is_pos(m_coeffs[i]), m_row_or_term_index, strict)); } } diff --git a/src/util/lp/var_register.h b/src/util/lp/var_register.h index e24667e946..7528e8d497 100644 --- a/src/util/lp/var_register.h +++ b/src/util/lp/var_register.h @@ -34,6 +34,8 @@ class var_register { return ret; } + const svector & vars() const { return m_local_vars_to_external; } + unsigned local_var_to_user_var(unsigned local_var) const { return m_local_vars_to_external[local_var]; } diff --git a/src/util/mpff.cpp b/src/util/mpff.cpp index 64e07f9f07..e00a25b1bd 100644 --- a/src/util/mpff.cpp +++ b/src/util/mpff.cpp @@ -1070,7 +1070,7 @@ bool mpff_manager::is_power_of_two(mpff const & a) const { template void mpff_manager::significand_core(mpff const & n, mpz_manager & m, mpz & t) { - m.set(t, m_precision, sig(n)); + m.set_digits(t, m_precision, sig(n)); } void mpff_manager::significand(mpff const & n, unsynch_mpz_manager & m, mpz & t) { @@ -1090,10 +1090,10 @@ void mpff_manager::to_mpz_core(mpff const & n, mpz_manager & m, mpz & t) to_buffer(0, n); unsigned * b = m_buffers[0].c_ptr(); shr(m_precision, b, -exp, m_precision, b); - m.set(t, m_precision, b); + m.set_digits(t, m_precision, b); } else { - m.set(t, m_precision, sig(n)); + m.set_digits(t, m_precision, sig(n)); if (exp > 0) { _scoped_numeral > p(m); m.set(p, 2); diff --git a/src/util/mpfx.cpp b/src/util/mpfx.cpp index e17a5e766f..2ebc548403 100644 --- a/src/util/mpfx.cpp +++ b/src/util/mpfx.cpp @@ -705,7 +705,7 @@ template void mpfx_manager::to_mpz_core(mpfx const & n, mpz_manager & m, mpz & t) { SASSERT(is_int(n)); unsigned * w = words(n); - m.set(t, m_int_part_sz, w+m_frac_part_sz); + m.set_digits(t, m_int_part_sz, w+m_frac_part_sz); if (is_neg(n)) m.neg(t); } diff --git a/src/util/mpq.cpp b/src/util/mpq.cpp index bc6f992135..4c636f1afb 100644 --- a/src/util/mpq.cpp +++ b/src/util/mpq.cpp @@ -315,6 +315,112 @@ unsigned mpq_manager::prev_power_of_two(mpq const & a) { return prev_power_of_two(_tmp); } + +template +template +void mpq_manager::lin_arith_op(mpq const& a, mpq const& b, mpq& c, mpz& g, mpz& tmp1, mpz& tmp2, mpz& tmp3) { + gcd(a.m_den, b.m_den, g); + if (is_one(g)) { + mul(a.m_num, b.m_den, tmp1); + mul(b.m_num, a.m_den, tmp2); + if (SUB) sub(tmp1, tmp2, c.m_num); else add(tmp1, tmp2, c.m_num); + mul(a.m_den, b.m_den, c.m_den); + } + else { + div(a.m_den, g, tmp3); + mul(tmp3, b.m_den, c.m_den); + mul(tmp3, b.m_num, tmp2); + div(b.m_den, g, tmp3); + mul(tmp3, a.m_num, tmp1); + if (SUB) sub(tmp1, tmp2, tmp3); else add(tmp1, tmp2, tmp3); + gcd(tmp3, g, tmp1); + if (is_one(tmp1)) { + set(c.m_num, tmp3); + } + else { + div(tmp3, tmp1, c.m_num); + div(c.m_den, tmp1, c.m_den); + } + } +} + +template +void mpq_manager::rat_mul(mpq const & a, mpq const & b, mpq & c, mpz& g1, mpz& g2, mpz& tmp1, mpz& tmp2) { + gcd(a.m_den, b.m_num, g1); + gcd(a.m_num, b.m_den, g2); + div(a.m_num, g2, tmp1); + div(b.m_num, g1, tmp2); + mul(tmp1, tmp2, c.m_num); + div(b.m_den, g2, tmp1); + div(a.m_den, g1, tmp2); + mul(tmp1, tmp2, c.m_den); +} + +template +void mpq_manager::rat_mul(mpz const & a, mpq const & b, mpq & c) { + STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); + mul(a, b.m_num, c.m_num); + set(c.m_den, b.m_den); + normalize(c); + STRACE("rat_mpq", tout << to_string(c) << "\n";); +} + + +template +void mpq_manager::rat_mul(mpq const & a, mpq const & b, mpq & c) { + STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); + if (SYNCH) { + mpz g1, g2, tmp1, tmp2; + rat_mul(a, b, c, g1, g2, tmp1, tmp2); + del(g1); + del(g2); + del(tmp1); + del(tmp2); + } + else { + mpz& g1 = m_n_tmp, &g2 = m_addmul_tmp.m_num, &tmp1 = m_add_tmp1, &tmp2 = m_add_tmp2; + rat_mul(a, b, c, g1, g2, tmp1, tmp2); + } + STRACE("rat_mpq", tout << to_string(c) << "\n";); +} + +template +void mpq_manager::rat_add(mpq const & a, mpq const & b, mpq & c) { + STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " + " << to_string(b) << " == ";); + if (SYNCH) { + mpz_stack tmp1, tmp2, tmp3, g; + lin_arith_op(a, b, c, g, tmp1, tmp2, tmp3); + del(tmp1); + del(tmp2); + del(tmp3); + del(g); + } + else { + mpz& g = m_n_tmp, &tmp1 = m_add_tmp1, &tmp2 = m_add_tmp2, &tmp3 = m_addmul_tmp.m_num; + lin_arith_op(a, b, c, g, tmp1, tmp2, tmp3); + } + STRACE("rat_mpq", tout << to_string(c) << "\n";); +} + +template +void mpq_manager::rat_sub(mpq const & a, mpq const & b, mpq & c) { + STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " - " << to_string(b) << " == ";); + if (SYNCH) { + mpz tmp1, tmp2, tmp3, g; + lin_arith_op(a, b, c, g, tmp1, tmp2, tmp3); + del(tmp1); + del(tmp2); + del(tmp3); + del(g); + } + else { + mpz& g = m_n_tmp, &tmp1 = m_add_tmp1, &tmp2 = m_add_tmp2, &tmp3 = m_addmul_tmp.m_num; + lin_arith_op(a, b, c, g, tmp1, tmp2, tmp3); + } + STRACE("rat_mpq", tout << to_string(c) << "\n";); +} + + template class mpq_manager; template class mpq_manager; diff --git a/src/util/mpq.h b/src/util/mpq.h index 01ad06db42..1bccabc745 100644 --- a/src/util/mpq.h +++ b/src/util/mpq.h @@ -74,27 +74,7 @@ class mpq_manager : public mpz_manager { } } - void rat_add(mpq const & a, mpq const & b, mpq & c) { - STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " + " << to_string(b) << " == ";); - if (SYNCH) { - mpz tmp1, tmp2; - mul(a.m_num, b.m_den, tmp1); - mul(b.m_num, a.m_den, tmp2); - mul(a.m_den, b.m_den, c.m_den); - add(tmp1, tmp2, c.m_num); - normalize(c); - del(tmp1); - del(tmp2); - } - else { - mul(a.m_num, b.m_den, m_add_tmp1); - mul(b.m_num, a.m_den, m_add_tmp2); - mul(a.m_den, b.m_den, c.m_den); - add(m_add_tmp1, m_add_tmp2, c.m_num); - normalize(c); - } - STRACE("rat_mpq", tout << to_string(c) << "\n";); - } + void rat_add(mpq const & a, mpq const & b, mpq & c); void rat_add(mpq const & a, mpz const & b, mpq & c) { STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " + " << to_string(b) << " == ";); @@ -115,46 +95,19 @@ class mpq_manager : public mpz_manager { STRACE("rat_mpq", tout << to_string(c) << "\n";); } - void rat_sub(mpq const & a, mpq const & b, mpq & c) { - STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " - " << to_string(b) << " == ";); - if (SYNCH) { - mpz tmp1, tmp2; - mul(a.m_num, b.m_den, tmp1); - mul(b.m_num, a.m_den, tmp2); - mul(a.m_den, b.m_den, c.m_den); - sub(tmp1, tmp2, c.m_num); - normalize(c); - del(tmp1); - del(tmp2); - } - else { - mul(a.m_num, b.m_den, m_add_tmp1); - mul(b.m_num, a.m_den, m_add_tmp2); - mul(a.m_den, b.m_den, c.m_den); - sub(m_add_tmp1, m_add_tmp2, c.m_num); - normalize(c); - } - STRACE("rat_mpq", tout << to_string(c) << "\n";); - } + void rat_sub(mpq const & a, mpq const & b, mpq & c); - void rat_mul(mpq const & a, mpq const & b, mpq & c) { - STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); - mul(a.m_num, b.m_num, c.m_num); - mul(a.m_den, b.m_den, c.m_den); - normalize(c); - STRACE("rat_mpq", tout << to_string(c) << "\n";); - } + void rat_mul(mpq const & a, mpq const & b, mpq & c); - void rat_mul(mpz const & a, mpq const & b, mpq & c) { - STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); - mul(a, b.m_num, c.m_num); - set(c.m_den, b.m_den); - normalize(c); - STRACE("rat_mpq", tout << to_string(c) << "\n";); - } + void rat_mul(mpz const & a, mpq const & b, mpq & c); bool rat_lt(mpq const & a, mpq const & b); + template + void lin_arith_op(mpq const& a, mpq const& b, mpq& c, mpz& g, mpz& tmp1, mpz& tmp2, mpz& tmp3); + + void rat_mul(mpq const & a, mpq const & b, mpq & c, mpz& g1, mpz& g2, mpz& tmp1, mpz& tmp2); + public: typedef mpq numeral; typedef mpq rational; @@ -746,10 +699,12 @@ class mpq_manager : public mpz_manager { reset_denominator(a); } - void set(mpz & a, unsigned sz, digit_t const * digits) { mpz_manager::set(a, sz, digits); } + void set(mpz & a, unsigned sz, digit_t const * digits) { + mpz_manager::set_digits(a, sz, digits); + } void set(mpq & a, unsigned sz, digit_t const * digits) { - mpz_manager::set(a.m_num, sz, digits); + mpz_manager::set_digits(a.m_num, sz, digits); reset_denominator(a); } diff --git a/src/util/mpz.cpp b/src/util/mpz.cpp index 861a31cfb2..d306fd8495 100644 --- a/src/util/mpz.cpp +++ b/src/util/mpz.cpp @@ -30,6 +30,7 @@ Revision History: #else #error No multi-precision library selected. #endif +#include // Available GCD algorithms // #define EUCLID_GCD @@ -45,9 +46,6 @@ Revision History: #define LEHMER_GCD #endif - -#include - #if defined(__GNUC__) #define _trailing_zeros32(X) __builtin_ctz(X) #else @@ -146,12 +144,7 @@ mpz_manager::mpz_manager(): else { m_init_cell_capacity = 6; } - for (unsigned i = 0; i < 2; i++) { - m_tmp[i] = allocate(m_init_cell_capacity); - m_arg[i] = allocate(m_init_cell_capacity); - m_arg[i]->m_size = 1; - } - set(m_int_min, -static_cast(INT_MIN)); + set(m_int_min, -static_cast(INT_MIN)); #else // GMP mpz_init(m_tmp); @@ -160,8 +153,6 @@ mpz_manager::mpz_manager(): mpz_set_ui(m_two32, UINT_MAX); mpz_set_ui(m_tmp, 1); mpz_add(m_two32, m_two32, m_tmp); - m_arg[0] = allocate(); - m_arg[1] = allocate(); mpz_init(m_uint64_max); unsigned max_l = static_cast(UINT64_MAX); unsigned max_h = static_cast(UINT64_MAX >> 32); @@ -193,16 +184,10 @@ mpz_manager::~mpz_manager() { del(m_two64); #ifndef _MP_GMP del(m_int_min); - for (unsigned i = 0; i < 2; i++) { - deallocate(m_tmp[i]); - deallocate(m_arg[i]); - } #else mpz_clear(m_tmp); mpz_clear(m_tmp2); mpz_clear(m_two32); - deallocate(m_arg[0]); - deallocate(m_arg[1]); mpz_clear(m_uint64_max); mpz_clear(m_int64_max); mpz_clear(m_int64_min); @@ -212,11 +197,47 @@ mpz_manager::~mpz_manager() { } template -void mpz_manager::set_big_i64(mpz & c, int64_t v) { +void mpz_manager::del(mpz & a) { + if (a.m_ptr) { + deallocate(a.m_owner == mpz_self, a.m_ptr); + a.m_ptr = nullptr; + a.m_kind = mpz_small; + a.m_owner = mpz_self; + } +} + +template +void mpz_manager::add(mpz const & a, mpz const & b, mpz & c) { + STRACE("mpz", tout << "[mpz] " << to_string(a) << " + " << to_string(b) << " == ";); + if (is_small(a) && is_small(b)) { + set_i64(c, i64(a) + i64(b)); + } + else { + big_add(a, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + +template +void mpz_manager::sub(mpz const & a, mpz const & b, mpz & c) { + STRACE("mpz", tout << "[mpz] " << to_string(a) << " - " << to_string(b) << " == ";); + if (is_small(a) && is_small(b)) { + set_i64(c, i64(a) - i64(b)); + } + else { + big_sub(a, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + +template +void mpz_manager::set_big_i64(mpz & c, int64 v) { #ifndef _MP_GMP - if (is_small(c)) { + if (c.m_ptr == nullptr) { c.m_ptr = allocate(m_init_cell_capacity); + c.m_owner = mpz_self; } + c.m_kind = mpz_ptr; SASSERT(capacity(c) >= m_init_cell_capacity); uint64_t _v; if (v < 0) { @@ -239,9 +260,11 @@ void mpz_manager::set_big_i64(mpz & c, int64_t v) { c.m_ptr->m_size = digits(c)[1] == 0 ? 1 : 2; } #else - if (is_small(c)) { + if (c.m_ptr == nullptr) { c.m_ptr = allocate(); + c.m_owner = mpz_self; } + c.m_kind = mpz_ptr; uint64_t _v; bool sign; if (v < 0) { @@ -253,9 +276,11 @@ void mpz_manager::set_big_i64(mpz & c, int64_t v) { sign = false; } mpz_set_ui(*c.m_ptr, static_cast(_v)); + MPZ_BEGIN_CRITICAL(); mpz_set_ui(m_tmp, static_cast(_v >> 32)); mpz_mul(m_tmp, m_tmp, m_two32); mpz_add(*c.m_ptr, *c.m_ptr, m_tmp); + MPZ_END_CRITICAL(); if (sign) mpz_neg(*c.m_ptr, *c.m_ptr); #endif @@ -264,9 +289,11 @@ void mpz_manager::set_big_i64(mpz & c, int64_t v) { template void mpz_manager::set_big_ui64(mpz & c, uint64_t v) { #ifndef _MP_GMP - if (is_small(c)) { + if (c.m_ptr == nullptr) { c.m_ptr = allocate(m_init_cell_capacity); + c.m_owner = mpz_self; } + c.m_kind = mpz_ptr; SASSERT(capacity(c) >= m_init_cell_capacity); c.m_val = 1; if (sizeof(digit_t) == sizeof(uint64_t)) { @@ -281,54 +308,71 @@ void mpz_manager::set_big_ui64(mpz & c, uint64_t v) { c.m_ptr->m_size = digits(c)[1] == 0 ? 1 : 2; } #else - if (is_small(c)) { + if (c.m_ptr == nullptr) { c.m_ptr = allocate(); + c.m_owner = mpz_self; } + c.m_kind = mpz_ptr; mpz_set_ui(*c.m_ptr, static_cast(v)); + MPZ_BEGIN_CRITICAL(); mpz_set_ui(m_tmp, static_cast(v >> 32)); mpz_mul(m_tmp, m_tmp, m_two32); mpz_add(*c.m_ptr, *c.m_ptr, m_tmp); + MPZ_END_CRITICAL(); #endif } -#ifndef _MP_GMP +#ifdef _MP_GMP + template -template -void mpz_manager::set(mpz & a, int sign, unsigned sz) { -#if 0 - static unsigned max_sz = 0; - if (sz > max_sz) { - max_sz = sz; - verbose_stream() << "max_sz: " << max_sz << "\n"; +mpz_manager::ensure_mpz_t(mpz_const& a) { + if (is_small(a)) { + m_result = &m_local; + mpz_init(m_local); + mpz_set_si(m_local, a.m_val); + } + else { + m_result = a.m_ptr; + } +} + +template +mpz_manager::ensure_mpz_t::~ensure_mpz_t() { + if (m_result == &m_local) { + mpz_clear(m_local); } +} + #endif + +#ifndef _MP_GMP +template +void mpz_manager::set(mpz_cell& src, mpz & a, int sign, unsigned sz) { unsigned i = sz; - for (; i > 0; --i) { - if (m_tmp[IDX]->m_digits[i-1] != 0) - break; - } + for (; i > 0 && src.m_digits[i-1] == 0; --i) ; if (i == 0) { - // m_tmp[IDX] is zero + // src is zero reset(a); return; } - if (i == 1 && m_tmp[IDX]->m_digits[0] <= INT_MAX) { - // m_tmp[IDX] fits is a fixnum - del(a); - a.m_val = sign < 0 ? -static_cast(m_tmp[IDX]->m_digits[0]) : static_cast(m_tmp[IDX]->m_digits[0]); + unsigned d = src.m_digits[0]; + if (i == 1 && d <= INT_MAX) { + // src fits is a fixnum + a.m_val = sign < 0 ? -static_cast(d) : static_cast(d); + a.m_kind = mpz_small; return; } + set_digits(a, i, src.m_digits); a.m_val = sign; - std::swap(a.m_ptr, m_tmp[IDX]); - a.m_ptr->m_size = i; - if (!m_tmp[IDX]) // 'a' was a small number - m_tmp[IDX] = allocate(m_init_cell_capacity); + + SASSERT(a.m_kind == mpz_ptr); } #endif + template void mpz_manager::set(mpz & a, char const * val) { reset(a); @@ -353,7 +397,7 @@ void mpz_manager::set(mpz & a, char const * val) { } template -void mpz_manager::set(mpz & target, unsigned sz, digit_t const * digits) { +void mpz_manager::set_digits(mpz & target, unsigned sz, digit_t const * digits) { // remove zero digits while (sz > 0 && digits[sz - 1] == 0) sz--; @@ -364,114 +408,299 @@ void mpz_manager::set(mpz & target, unsigned sz, digit_t const * digits) else { #ifndef _MP_GMP target.m_val = 1; // number is positive. - if (is_small(target)) { + if (target.m_ptr == nullptr) { unsigned c = sz < m_init_cell_capacity ? m_init_cell_capacity : sz; target.m_ptr = allocate(c); target.m_ptr->m_size = sz; target.m_ptr->m_capacity = c; + target.m_kind = mpz_ptr; + target.m_owner = mpz_self; memcpy(target.m_ptr->m_digits, digits, sizeof(digit_t) * sz); } + else if (capacity(target) < sz) { + SASSERT(sz > m_init_cell_capacity); + mpz_cell* ptr = allocate(sz); + memcpy(ptr->m_digits, digits, sizeof(digit_t) * sz); + ptr->m_size = sz; + ptr->m_capacity = sz; + deallocate(target); + target.m_val = 1; + target.m_ptr = ptr; + target.m_kind = mpz_ptr; + target.m_owner = mpz_self; + } else { - if (capacity(target) < sz) { - SASSERT(sz > m_init_cell_capacity); - deallocate(target.m_ptr); - target.m_ptr = allocate(sz); - target.m_ptr->m_size = sz; - target.m_ptr->m_capacity = sz; - memcpy(target.m_ptr->m_digits, digits, sizeof(digit_t) * sz); - } - else { - target.m_ptr->m_size = sz; + target.m_ptr->m_size = sz; + if (target.m_ptr->m_digits != digits) memcpy(target.m_ptr->m_digits, digits, sizeof(digit_t) * sz); - } + target.m_kind = mpz_ptr; } #else - mk_big(target); - // reset - mpz_set_ui(*target.m_ptr, digits[sz - 1]); - SASSERT(sz > 0); - unsigned i = sz - 1; - while (i > 0) { - --i; - mpz_mul_2exp(*target.m_ptr, *target.m_ptr, 32); - mpz_set_ui(m_tmp, digits[i]); - mpz_add(*target.m_ptr, *target.m_ptr, m_tmp); - } + mk_big(target); + // reset + mpz_set_ui(*target.m_ptr, digits[sz - 1]); + SASSERT(sz > 0); + unsigned i = sz - 1; + MPZ_BEGIN_CRITICAL(); + while (i > 0) { + --i; + mpz_mul_2exp(*target.m_ptr, *target.m_ptr, 32); + mpz_set_ui(m_tmp, digits[i]); + mpz_add(*target.m_ptr, *target.m_ptr, m_tmp); + } + MPZ_END_CRITICAL(); #endif } } +template +void mpz_manager::mul(mpz const & a, mpz const & b, mpz & c) { + STRACE("mpz", tout << "[mpz] " << to_string(a) << " * " << to_string(b) << " == ";); + if (is_small(a) && is_small(b)) { + set_i64(c, i64(a) * i64(b)); + } + else { + big_mul(a, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + +// d <- a + b*c +template +void mpz_manager::addmul(mpz const & a, mpz const & b, mpz const & c, mpz & d) { + if (is_one(b)) { + add(a, c, d); + } + else if (is_minus_one(b)) { + sub(a, c, d); + } + else { + mpz tmp; + mul(b,c,tmp); + add(a,tmp,d); + del(tmp); + } +} + + +// d <- a - b*c +template +void mpz_manager::submul(mpz const & a, mpz const & b, mpz const & c, mpz & d) { + if (is_one(b)) { + sub(a, c, d); + } + else if (is_minus_one(b)) { + add(a, c, d); + } + else { + mpz tmp; + mul(b,c,tmp); + sub(a,tmp,d); + del(tmp); + } +} + + +template +void mpz_manager::machine_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r) { + STRACE("mpz", tout << "[mpz-ext] divrem(" << to_string(a) << ", " << to_string(b) << ") == ";); + if (is_small(a) && is_small(b)) { + int64 _a = i64(a); + int64 _b = i64(b); + set_i64(q, _a / _b); + set_i64(r, _a % _b); + } + else { + big_div_rem(a, b, q, r); + } + STRACE("mpz", tout << "(" << to_string(q) << ", " << to_string(r) << ")\n";); +} + +template +void mpz_manager::machine_div(mpz const & a, mpz const & b, mpz & c) { + STRACE("mpz", tout << "[mpz-ext] machine-div(" << to_string(a) << ", " << to_string(b) << ") == ";); + if (is_small(a) && is_small(b)) { + set_i64(c, i64(a) / i64(b)); + } + else { + big_div(a, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + +template +void mpz_manager::rem(mpz const & a, mpz const & b, mpz & c) { + STRACE("mpz", tout << "[mpz-ext] rem(" << to_string(a) << ", " << to_string(b) << ") == ";); + if (is_small(a) && is_small(b)) { + set_i64(c, i64(a) % i64(b)); + } + else { + big_rem(a, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + + +template +void mpz_manager::div_gcd(mpz const& a, mpz const& b, mpz & c) { + STRACE("mpz", tout << "[mpz-ext] div(" << to_string(a) << ", " << to_string(b) << ") == ";); + if (is_one(b)) { + set(c, a); + } + else { + machine_div(a, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + +template +void mpz_manager::div(mpz const & a, mpz const & b, mpz & c) { + STRACE("mpz", tout << "[mpz-ext] div(" << to_string(a) << ", " << to_string(b) << ") == ";); + if (is_one(b)) { + set(c, a); + } + else if (is_neg(a)) { + mpz tmp; + machine_div_rem(a, b, c, tmp); + if (!is_zero(tmp)) { + if (is_neg(b)) + add(c, mk_z(1), c); + else + sub(c, mk_z(1), c); + } + del(tmp); + } + else { + machine_div(a, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + +template +void mpz_manager::mod(mpz const & a, mpz const & b, mpz & c) { + STRACE("mpz", tout << "[mpz-ext] mod(" << to_string(a) << ", " << to_string(b) << ") == ";); + rem(a, b, c); + if (is_neg(c)) { + if (is_pos(b)) + add(c, b, c); + else + sub(c, b, c); + } + STRACE("mpz", tout << to_string(c) << "\n";); +} + +template +void mpz_manager::neg(mpz & a) { + STRACE("mpz", tout << "[mpz] 0 - " << to_string(a) << " == ";); + if (is_small(a) && a.m_val == INT_MIN) { + // neg(INT_MIN) is not a small int + set_big_i64(a, - static_cast(INT_MIN)); + return; + } +#ifndef _MP_GMP + a.m_val = -a.m_val; +#else + if (is_small(a)) { + a.m_val = -a.m_val; + } + else { + mpz_neg(*a.m_ptr, *a.m_ptr); + } +#endif + STRACE("mpz", tout << to_string(a) << "\n";); +} + +template +void mpz_manager::abs(mpz & a) { + if (is_small(a)) { + if (a.m_val < 0) { + if (a.m_val == INT_MIN) { + // abs(INT_MIN) is not a small int + set_big_i64(a, - static_cast(INT_MIN)); + } + else + a.m_val = -a.m_val; + } + } + else { +#ifndef _MP_GMP + a.m_val = 1; +#else + mpz_abs(*a.m_ptr, *a.m_ptr); +#endif + } +} + + +// TBD: replace use of 'tmp' by 'c'. #ifndef _MP_GMP template template void mpz_manager::big_add_sub(mpz const & a, mpz const & b, mpz & c) { - int sign_a; - int sign_b; - mpz_cell * cell_a; - mpz_cell * cell_b; - get_sign_cell<0>(a, sign_a, cell_a); - get_sign_cell<1>(b, sign_b, cell_b); + sign_cell ca(*this, a), cb(*this, b); + int sign_b = cb.sign(); + mpz_stack tmp; if (SUB) sign_b = -sign_b; size_t real_sz; - if (sign_a == sign_b) { - unsigned sz = std::max(cell_a->m_size, cell_b->m_size)+1; - ensure_tmp_capacity<0>(sz); - m_mpn_manager.add(cell_a->m_digits, cell_a->m_size, - cell_b->m_digits, cell_b->m_size, - m_tmp[0]->m_digits, sz, &real_sz); + if (ca.sign() == sign_b) { + unsigned sz = std::max(ca.cell()->m_size, cb.cell()->m_size)+1; + allocate_if_needed(tmp, sz); + m_mpn_manager.add(ca.cell()->m_digits, ca.cell()->m_size, + cb.cell()->m_digits, cb.cell()->m_size, + tmp.m_ptr->m_digits, sz, &real_sz); SASSERT(real_sz <= sz); - set<0>(c, sign_a, static_cast(real_sz)); + set(*tmp.m_ptr, c, ca.sign(), static_cast(real_sz)); } else { digit_t borrow; - int r = m_mpn_manager.compare(cell_a->m_digits, cell_a->m_size, - cell_b->m_digits, cell_b->m_size); + int r = m_mpn_manager.compare(ca.cell()->m_digits, ca.cell()->m_size, + cb.cell()->m_digits, cb.cell()->m_size); if (r == 0) { reset(c); } else if (r < 0) { // a < b - unsigned sz = cell_b->m_size; - ensure_tmp_capacity<0>(sz); - m_mpn_manager.sub(cell_b->m_digits, - cell_b->m_size, - cell_a->m_digits, - cell_a->m_size, - m_tmp[0]->m_digits, + unsigned sz = cb.cell()->m_size; + allocate_if_needed(tmp, sz); + m_mpn_manager.sub(cb.cell()->m_digits, + cb.cell()->m_size, + ca.cell()->m_digits, + ca.cell()->m_size, + tmp.m_ptr->m_digits, &borrow); SASSERT(borrow == 0); - set<0>(c, sign_b, sz); + set(*tmp.m_ptr, c, sign_b, sz); } else { // a > b - unsigned sz = cell_a->m_size; - ensure_tmp_capacity<0>(sz); - m_mpn_manager.sub(cell_a->m_digits, - cell_a->m_size, - cell_b->m_digits, - cell_b->m_size, - m_tmp[0]->m_digits, + unsigned sz = ca.cell()->m_size; + allocate_if_needed(tmp, sz); + m_mpn_manager.sub(ca.cell()->m_digits, + ca.cell()->m_size, + cb.cell()->m_digits, + cb.cell()->m_size, + tmp.m_ptr->m_digits, &borrow); SASSERT(borrow == 0); - set<0>(c, sign_a, sz); + set(*tmp.m_ptr, c, ca.sign(), sz); } } + del(tmp); } #endif + + template void mpz_manager::big_add(mpz const & a, mpz const & b, mpz & c) { #ifndef _MP_GMP big_add_sub(a, b, c); #else // GMP version - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_add(*c.m_ptr, *arg0, *arg1); + mpz_add(*c.m_ptr, a1(), b1()); #endif } @@ -481,97 +710,87 @@ void mpz_manager::big_sub(mpz const & a, mpz const & b, mpz & c) { big_add_sub(a, b, c); #else // GMP version - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_sub(*c.m_ptr, *arg0, *arg1); + mpz_sub(*c.m_ptr, a1(), b1()); #endif } template void mpz_manager::big_mul(mpz const & a, mpz const & b, mpz & c) { #ifndef _MP_GMP - int sign_a; - int sign_b; - mpz_cell * cell_a; - mpz_cell * cell_b; - get_sign_cell<0>(a, sign_a, cell_a); - get_sign_cell<1>(b, sign_b, cell_b); - unsigned sz = cell_a->m_size + cell_b->m_size; - ensure_tmp_capacity<0>(sz); - m_mpn_manager.mul(cell_a->m_digits, - cell_a->m_size, - cell_b->m_digits, - cell_b->m_size, - m_tmp[0]->m_digits); - set<0>(c, sign_a == sign_b ? 1 : -1, sz); + // TBD replace tmp by c. + mpz_stack tmp; + sign_cell ca(*this, a), cb(*this, b); + unsigned sz = ca.cell()->m_size + cb.cell()->m_size; + allocate_if_needed(tmp, sz); + m_mpn_manager.mul(ca.cell()->m_digits, + ca.cell()->m_size, + cb.cell()->m_digits, + cb.cell()->m_size, + tmp.m_ptr->m_digits); + set(*tmp.m_ptr, c, ca.sign() == cb.sign() ? 1 : -1, sz); + del(tmp); #else // GMP version - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_mul(*c.m_ptr, *arg0, *arg1); + mpz_mul(*c.m_ptr, a1(), b1()); +#endif +} + + +template +void mpz_manager::big_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r) { +#ifndef _MP_GMP + quot_rem_core(a, b, q, r); +#else + // GMP version + ensure_mpz a1(a), b1(b); + mk_big(q); + mk_big(r); + mpz_tdiv_qr(*q.m_ptr, *r.m_ptr, a1(), b1()); #endif } #ifndef _MP_GMP template template -void mpz_manager::quot_rem_core(mpz const & a, mpz const & b, mpz & q, mpz & r) { +void mpz_manager::quot_rem_core(mpz const & a, mpz const & b, mpz & q, mpz & r) +{ /* +26 / +7 = +3, remainder is +5 -26 / +7 = -3, remainder is -5 +26 / -7 = -3, remainder is +5 -26 / -7 = +3, remainder is -5 */ - int sign_a; - int sign_b; - mpz_cell * cell_a; - mpz_cell * cell_b; - get_sign_cell<0>(a, sign_a, cell_a); - get_sign_cell<1>(b, sign_b, cell_b); - if (cell_b->m_size > cell_a->m_size) { + + mpz_stack q1, r1; + sign_cell ca(*this, a), cb(*this, b); + if (cb.cell()->m_size > ca.cell()->m_size) { if (MODE == REM_ONLY || MODE == QUOT_AND_REM) set(r, a); if (MODE == QUOT_ONLY || MODE == QUOT_AND_REM) reset(q); return; } - unsigned q_sz = cell_a->m_size - cell_b->m_size + 1; - unsigned r_sz = cell_b->m_size; - ensure_tmp_capacity<0>(q_sz); - ensure_tmp_capacity<1>(r_sz); - m_mpn_manager.div(cell_a->m_digits, cell_a->m_size, - cell_b->m_digits, cell_b->m_size, - m_tmp[0]->m_digits, - m_tmp[1]->m_digits); + unsigned q_sz = ca.cell()->m_size - cb.cell()->m_size + 1; + unsigned r_sz = cb.cell()->m_size; + allocate_if_needed(q1, q_sz); + allocate_if_needed(r1, r_sz); + m_mpn_manager.div(ca.cell()->m_digits, ca.cell()->m_size, + cb.cell()->m_digits, cb.cell()->m_size, + q1.m_ptr->m_digits, + r1.m_ptr->m_digits); if (MODE == QUOT_ONLY || MODE == QUOT_AND_REM) - set<0>(q, sign_a == sign_b ? 1 : -1, q_sz); + set(*q1.m_ptr, q, ca.sign() == cb.sign() ? 1 : -1, q_sz); if (MODE == REM_ONLY || MODE == QUOT_AND_REM) - set<1>(r, sign_a, r_sz); + set(*r1.m_ptr, r, ca.sign(), r_sz); + del(q1); + del(r1); } #endif -template -void mpz_manager::big_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r) { -#ifndef _MP_GMP - quot_rem_core(a, b, q, r); -#else - // GMP version - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); - mk_big(q); - mk_big(r); - mpz_tdiv_qr(*q.m_ptr, *r.m_ptr, *arg0, *arg1); -#endif -} - template void mpz_manager::big_div(mpz const & a, mpz const & b, mpz & c) { #ifndef _MP_GMP @@ -580,12 +799,9 @@ void mpz_manager::big_div(mpz const & a, mpz const & b, mpz & c) { SASSERT(is_zero(dummy)); #else // GMP version - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_tdiv_q(*c.m_ptr, *arg0, *arg1); + mpz_tdiv_q(*c.m_ptr, a1(), b1()); #endif } @@ -597,18 +813,16 @@ void mpz_manager::big_rem(mpz const & a, mpz const & b, mpz & c) { SASSERT(is_zero(dummy)); #else // GMP version - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_tdiv_r(*c.m_ptr, *arg0, *arg1); + mpz_tdiv_r(*c.m_ptr, a1(), b1()); #endif } template void mpz_manager::gcd(mpz const & a, mpz const & b, mpz & c) { static_assert(sizeof(a.m_val) == sizeof(int), "size mismatch"); + static_assert(sizeof(mpz) <= 16, "mpz size overflow"); if (is_small(a) && is_small(b) && a.m_val != INT_MIN && b.m_val != INT_MIN) { int _a = a.m_val; int _b = b.m_val; @@ -619,12 +833,9 @@ void mpz_manager::gcd(mpz const & a, mpz const & b, mpz & c) { } else { #ifdef _MP_GMP - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_gcd(*c.m_ptr, *arg0, *arg1); + mpz_gcd(*c.m_ptr, a1(), b1()); return; #endif if (is_zero(a)) { @@ -879,6 +1090,14 @@ unsigned mpz_manager::size_info(mpz const & a) { #endif } +template +mpz_manager::sign_cell::sign_cell(mpz_manager& m, mpz const& a): + m_local(reinterpret_cast(m_bytes)), m_a(a) { + m_local.m_ptr->m_capacity = capacity; + m.get_sign_cell(a, m_sign, m_cell, m_local.m_ptr); +} + + template struct mpz_manager::sz_lt { mpz_manager & m; @@ -891,6 +1110,7 @@ struct mpz_manager::sz_lt { } }; + template void mpz_manager::gcd(unsigned sz, mpz const * as, mpz & g) { #if 0 @@ -1070,8 +1290,8 @@ void mpz_manager::bitwise_or(mpz const & a, mpz const & b, mpz & c) { SASSERT(is_nonneg(b)); TRACE("mpz", tout << "is_small(a): " << is_small(a) << ", is_small(b): " << is_small(b) << "\n";); if (is_small(a) && is_small(b)) { - del(c); c.m_val = a.m_val | b.m_val; + c.m_kind = mpz_small; } else { #ifndef _MP_GMP @@ -1104,12 +1324,9 @@ void mpz_manager::bitwise_or(mpz const & a, mpz const & b, mpz & c) { } del(a1); del(b1); del(a2); del(b2); del(m); del(tmp); #else - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_ior(*c.m_ptr, *arg0, *arg1); + mpz_ior(*c.m_ptr, a1(), b1()); #endif } } @@ -1119,8 +1336,8 @@ void mpz_manager::bitwise_and(mpz const & a, mpz const & b, mpz & c) { SASSERT(is_nonneg(a)); SASSERT(is_nonneg(b)); if (is_small(a) && is_small(b)) { - del(c); c.m_val = a.m_val & b.m_val; + c.m_kind = mpz_small; } else { #ifndef _MP_GMP @@ -1142,12 +1359,9 @@ void mpz_manager::bitwise_and(mpz const & a, mpz const & b, mpz & c) { } del(a1); del(b1); del(a2); del(b2); del(m); del(tmp); #else - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_and(*c.m_ptr, *arg0, *arg1); + mpz_and(*c.m_ptr, a1(), b1()); #endif } } @@ -1187,12 +1401,9 @@ void mpz_manager::bitwise_xor(mpz const & a, mpz const & b, mpz & c) { } del(a1); del(b1); del(a2); del(b2); del(m); del(tmp); #else - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); + ensure_mpz a1(a), b1(b); mk_big(c); - mpz_xor(*c.m_ptr, *arg0, *arg1); + mpz_xor(*c.m_ptr, a1(), b1()); #endif } } @@ -1238,24 +1449,27 @@ void mpz_manager::big_set(mpz & target, mpz const & source) { if (&target == &source) return; target.m_val = source.m_val; - if (is_small(target)) { + if (target.m_ptr == nullptr) { + target.m_ptr = allocate(capacity(source)); + target.m_ptr->m_size = size(source); + target.m_ptr->m_capacity = capacity(source); + target.m_kind = mpz_ptr; + target.m_owner = mpz_self; + memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source)); + } + else if (capacity(target) < size(source)) { + deallocate(target); target.m_ptr = allocate(capacity(source)); target.m_ptr->m_size = size(source); target.m_ptr->m_capacity = capacity(source); + target.m_kind = mpz_ptr; + target.m_owner = mpz_self; memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source)); } else { - if (capacity(target) < size(source)) { - deallocate(target.m_ptr); - target.m_ptr = allocate(capacity(source)); - target.m_ptr->m_size = size(source); - target.m_ptr->m_capacity = capacity(source); - memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source)); - } - else { - target.m_ptr->m_size = size(source); - memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source)); - } + target.m_ptr->m_size = size(source); + memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source)); + target.m_kind = mpz_ptr; } #else // GMP version @@ -1267,19 +1481,14 @@ void mpz_manager::big_set(mpz & target, mpz const & source) { template int mpz_manager::big_compare(mpz const & a, mpz const & b) { #ifndef _MP_GMP - int sign_a; - int sign_b; - mpz_cell * cell_a; - mpz_cell * cell_b; - get_sign_cell<0>(a, sign_a, cell_a); - get_sign_cell<1>(b, sign_b, cell_b); - - if (sign_a > 0) { + sign_cell ca(*this, a), cb(*this, b); + + if (ca.sign() > 0) { // a is positive - if (sign_b > 0) { + if (cb.sign() > 0) { // a & b are positive - return m_mpn_manager.compare(cell_a->m_digits, cell_a->m_size, - cell_b->m_digits, cell_b->m_size); + return m_mpn_manager.compare(ca.cell()->m_digits, ca.cell()->m_size, + cb.cell()->m_digits, cb.cell()->m_size); } else { // b is negative @@ -1288,23 +1497,20 @@ int mpz_manager::big_compare(mpz const & a, mpz const & b) { } else { // a is negative - if (sign_b > 0) { + if (cb.sign() > 0) { // b is positive return -1; // a < b } else { // a & b are negative - return m_mpn_manager.compare(cell_b->m_digits, cell_b->m_size, - cell_a->m_digits, cell_a->m_size); + return m_mpn_manager.compare(cb.cell()->m_digits, cb.cell()->m_size, + ca.cell()->m_digits, ca.cell()->m_size); } } #else // GMP version - mpz_t * arg0; - mpz_t * arg1; - get_arg<0>(a, arg0); - get_arg<1>(b, arg1); - return mpz_cmp(*arg0, *arg1); + ensure_mpz a1(a), b1(b); + return mpz_cmp(a1(), b1()); #endif } @@ -1369,6 +1575,7 @@ uint64_t mpz_manager::get_uint64(mpz const & a) const { return mpz_get_ui(*a.m_ptr); } else { + MPZ_BEGIN_CRITICAL(); mpz_manager * _this = const_cast(this); mpz_set(_this->m_tmp, *a.m_ptr); mpz_mod(_this->m_tmp, m_tmp, m_two32); @@ -1376,6 +1583,7 @@ uint64_t mpz_manager::get_uint64(mpz const & a) const { mpz_set(_this->m_tmp, *a.m_ptr); mpz_div(_this->m_tmp, m_tmp, m_two32); r += static_cast(mpz_get_ui(m_tmp)) << static_cast(32); + MPZ_END_CRITICAL(); return r; } #endif @@ -1400,11 +1608,13 @@ int64_t mpz_manager::get_int64(mpz const & a) const { return mpz_get_si(*a.m_ptr); } else { + MPZ_BEGIN_CRITICAL(); mpz_manager * _this = const_cast(this); mpz_mod(_this->m_tmp, *a.m_ptr, m_two32); int64_t r = static_cast(mpz_get_ui(m_tmp)); mpz_div(_this->m_tmp, *a.m_ptr, m_two32); r += static_cast(mpz_get_si(m_tmp)) << static_cast(32); + MPZ_END_CRITICAL(); return r; } #endif @@ -1512,8 +1722,8 @@ void mpz_manager::power(mpz const & a, unsigned p, mpz & b) { if (is_small(a)) { if (a.m_val == 2) { if (p < 8 * sizeof(int) - 1) { - del(b); b.m_val = 1 << p; + b.m_kind = mpz_small; } else { unsigned sz = p/(8 * sizeof(digit_t)) + 1; @@ -1526,6 +1736,7 @@ void mpz_manager::power(mpz const & a, unsigned p, mpz & b) { b.m_ptr->m_digits[i] = 0; b.m_ptr->m_digits[sz-1] = 1 << shift; b.m_val = 1; + b.m_kind = mpz_ptr; } return; } @@ -1608,8 +1819,10 @@ void mpz_manager::ensure_capacity(mpz & a, unsigned capacity) { return; if (capacity < m_init_cell_capacity) capacity = m_init_cell_capacity; - if (is_small(a)) { + if (a.m_ptr == nullptr) { a.m_ptr = allocate(capacity); + a.m_owner = mpz_self; + a.m_kind = mpz_ptr; SASSERT(a.m_ptr->m_capacity == capacity); if (a.m_val == INT_MIN) { unsigned intmin_sz = m_int_min.m_ptr->m_size; @@ -1629,17 +1842,17 @@ void mpz_manager::ensure_capacity(mpz & a, unsigned capacity) { a.m_ptr->m_size = 1; } } - else { - if (a.m_ptr->m_capacity >= capacity) - return; + else if (a.m_ptr->m_capacity < capacity) { mpz_cell * new_cell = allocate(capacity); SASSERT(new_cell->m_capacity == capacity); unsigned old_sz = a.m_ptr->m_size; new_cell->m_size = old_sz; for (unsigned i = 0; i < old_sz; i++) new_cell->m_digits[i] = a.m_ptr->m_digits[i]; - deallocate(a.m_ptr); + deallocate(a); a.m_ptr = new_cell; + a.m_owner = mpz_self; + a.m_kind = mpz_ptr; } } @@ -1662,8 +1875,8 @@ void mpz_manager::normalize(mpz & a) { if (i == 1 && ds[0] <= INT_MAX) { // a is small int val = a.m_val < 0 ? -static_cast(ds[0]) : static_cast(ds[0]); - del(a); a.m_val = val; + a.m_kind = mpz_small; return; } // adjust size @@ -1731,10 +1944,9 @@ void mpz_manager::machine_div2k(mpz & a, unsigned k) { c->m_size = new_sz; normalize(a); #else + ensure_mpz a1(a); MPZ_BEGIN_CRITICAL(); - mpz_t * arg0; - get_arg<0>(a, arg0); - mpz_tdiv_q_2exp(m_tmp, *arg0, k); + mpz_tdiv_q_2exp(m_tmp, a1(), k); mk_big(a); mpz_swap(*a.m_ptr, m_tmp); MPZ_END_CRITICAL(); @@ -1796,10 +2008,9 @@ void mpz_manager::mul2k(mpz & a, unsigned k) { normalize(a); TRACE("mpz_mul2k", tout << "mul2k result:\n" << to_string(a) << "\n";); #else - mpz_t * arg0; - get_arg<0>(a, arg0); + ensure_mpz a1(a); mk_big(a); - mpz_mul_2exp(*a.m_ptr, *arg0, k); + mpz_mul_2exp(*a.m_ptr, a1(), k); #endif } @@ -1902,8 +2113,10 @@ unsigned mpz_manager::mlog2(mpz const & a) { else return (sz - 1)*32 + ::log2(static_cast(ds[sz-1])); #else + MPZ_BEGIN_CRITICAL(); mpz_neg(m_tmp, *a.m_ptr); unsigned r = mpz_sizeinbase(m_tmp, 2); + MPZ_END_CRITICAL(); SASSERT(r > 0); return r - 1; #endif @@ -2087,6 +2300,7 @@ bool mpz_manager::decompose(mpz const & a, svector & digits) { return a.m_val < 0; #else bool r = is_neg(a); + MPZ_BEGIN_CRITICAL(); mpz_set(m_tmp, *a.m_ptr); mpz_abs(m_tmp, m_tmp); while (mpz_sgn(m_tmp) != 0) { @@ -2095,6 +2309,7 @@ bool mpz_manager::decompose(mpz const & a, svector & digits) { digits.push_back(v); mpz_tdiv_q_2exp(m_tmp, m_tmp, 32); } + MPZ_END_CRITICAL(); return r; #endif } diff --git a/src/util/mpz.h b/src/util/mpz.h index 92c6d0d108..967ed60263 100644 --- a/src/util/mpz.h +++ b/src/util/mpz.h @@ -64,6 +64,7 @@ class mpz_cell { digit_t m_digits[0]; friend class mpz_manager; friend class mpz_manager; + friend class mpz_stack; }; #else #include @@ -72,17 +73,25 @@ class mpz_cell { /** \brief Multi-precision integer. - If m_ptr == 0, the it is a small number and the value is stored at m_val. - Otherwise, m_val contains the sign (-1 negative, 1 positive), and m_ptr points to a mpz_cell that - store the value. <<< This last statement is true only in Windows. + If m_kind == mpz_small, it is a small number and the value is stored in m_val. + If m_kind == mpz_ptr, the value is stored in m_ptr and m_ptr != nullptr. + m_val contains the sign (-1 negative, 1 positive) + under winodws, m_ptr points to a mpz_cell that store the value. */ + +enum mpz_kind { mpz_small = 0, mpz_ptr = 1}; +enum mpz_owner { mpz_self = 0, mpz_ext = 1}; + class mpz { - int m_val; #ifndef _MP_GMP - mpz_cell * m_ptr; + typedef mpz_cell mpz_type; #else - mpz_t * m_ptr; + typedef mpz_t mpz_type; #endif + int m_val; + unsigned m_kind:1; + unsigned m_owner:1; + mpz_type * m_ptr; friend class mpz_manager; friend class mpz_manager; friend class mpq_manager; @@ -90,18 +99,37 @@ class mpz { friend class mpq; friend class mpbq; friend class mpbq_manager; + friend class mpz_stack; mpz & operator=(mpz const & other) { UNREACHABLE(); return *this; } public: - mpz(int v):m_val(v), m_ptr(nullptr) {} - mpz():m_val(0), m_ptr(nullptr) {} - mpz(mpz && other) : m_val(other.m_val), m_ptr(nullptr) { + mpz(int v):m_val(v), m_kind(mpz_small), m_owner(mpz_self), m_ptr(nullptr) {} + mpz():m_val(0), m_kind(mpz_small), m_owner(mpz_self), m_ptr(nullptr) {} + mpz(mpz_type* ptr): m_val(0), m_kind(mpz_small), m_owner(mpz_ext), m_ptr(ptr) { SASSERT(ptr);} + mpz(mpz && other) : m_val(other.m_val), m_kind(mpz_small), m_owner(mpz_self), m_ptr(nullptr) { std::swap(m_ptr, other.m_ptr); + unsigned o = m_owner; m_owner = other.m_owner; other.m_owner = o; + unsigned k = m_kind; m_kind = other.m_kind; other.m_kind = k; } void swap(mpz & other) { std::swap(m_val, other.m_val); std::swap(m_ptr, other.m_ptr); + unsigned o = m_owner; m_owner = other.m_owner; other.m_owner = o; + unsigned k = m_kind; m_kind = other.m_kind; other.m_kind = k; + } +}; + +#ifndef _MP_GMP +class mpz_stack : public mpz { + static const unsigned capacity = 8; + unsigned char m_bytes[sizeof(mpz_cell) + sizeof(digit_t) * capacity]; +public: + mpz_stack():mpz(reinterpret_cast(m_bytes)) { + m_ptr->m_capacity = capacity; } }; +#else +class mpz_stack : public mpz {}; +#endif inline void swap(mpz & m1, mpz & m2) { m1.swap(m2); } @@ -115,57 +143,56 @@ class mpz_manager { #ifndef _MP_GMP unsigned m_init_cell_capacity; - mpz_cell * m_tmp[2]; - mpz_cell * m_arg[2]; mpz m_int_min; - static unsigned cell_size(unsigned capacity) { return sizeof(mpz_cell) + sizeof(digit_t) * capacity; } + static unsigned cell_size(unsigned capacity) { + return sizeof(mpz_cell) + sizeof(digit_t) * capacity; + } mpz_cell * allocate(unsigned capacity) { SASSERT(capacity >= m_init_cell_capacity); + MPZ_BEGIN_CRITICAL(); mpz_cell * cell = reinterpret_cast(m_allocator.allocate(cell_size(capacity))); + MPZ_END_CRITICAL(); cell->m_capacity = capacity; return cell; } + + void deallocate(mpz& n) { + if (n.m_ptr) { + deallocate(n.m_owner == mpz_self, n.m_ptr); + n.m_ptr = nullptr; + n.m_kind = mpz_small; + } + } // make sure that n is a big number and has capacity equal to at least c. void allocate_if_needed(mpz & n, unsigned c) { - if (c < m_init_cell_capacity) - c = m_init_cell_capacity; - if (is_small(n)) { + c = std::max(c, m_init_cell_capacity); + if (n.m_ptr == nullptr || capacity(n) < c) { + deallocate(n); n.m_val = 1; + n.m_kind = mpz_ptr; + n.m_owner = mpz_self; n.m_ptr = allocate(c); - n.m_ptr->m_capacity = c; } - else if (capacity(n) < c) { - deallocate(n.m_ptr); - n.m_val = 1; - n.m_ptr = allocate(c); - n.m_ptr->m_capacity = c; - } - } - - void deallocate(mpz_cell * ptr) { - m_allocator.deallocate(cell_size(ptr->m_capacity), ptr); } - /** - \brief Make sure that m_tmp[IDX] can hold the given number of digits - */ - template - void ensure_tmp_capacity(unsigned capacity) { - if (m_tmp[IDX]->m_capacity >= capacity) - return; - deallocate(m_tmp[IDX]); - unsigned new_capacity = (3 * capacity + 1) >> 1; - m_tmp[IDX] = allocate(new_capacity); - SASSERT(m_tmp[IDX]->m_capacity >= capacity); + void deallocate(bool is_heap, mpz_cell * ptr) { + if (is_heap) { + MPZ_BEGIN_CRITICAL(); + m_allocator.deallocate(cell_size(ptr->m_capacity), ptr); + MPZ_END_CRITICAL(); + } } // Expand capacity of a while preserving its content. void ensure_capacity(mpz & a, unsigned sz); void normalize(mpz & a); + + void clear(mpz& n) { } + #else // GMP code mpz_t m_tmp, m_tmp2; @@ -175,22 +202,34 @@ class mpz_manager { mpz_t m_int64_max; mpz_t m_int64_min; - mpz_t * allocate() { + mpz_t * allocate() { + MPZ_BEGIN_CRITICAL(); mpz_t * cell = reinterpret_cast(m_allocator.allocate(sizeof(mpz_t))); + MPZ_END_CRITICAL(); mpz_init(*cell); return cell; } - void deallocate(mpz_t * ptr) { mpz_clear(*ptr); m_allocator.deallocate(sizeof(mpz_t), ptr); } + void deallocate(bool is_heap, mpz_t * ptr) { + mpz_clear(*ptr); + if (is_heap) { + MPZ_BEGIN_CRITICAL(); + m_allocator.deallocate(sizeof(mpz_t), ptr); + MPZ_END_CRITICAL(); + } + } + + void clear(mpz& n) { if (n.m_ptr) mpz_clear(n.m_ptr); } #endif + + mpz m_two64; /** - \brief Set \c a with the value stored at m_tmp[IDX], and the given sign. - \c sz is an overapproximation of the size of the number stored at \c tmp. + \brief Set \c a with the value stored at src, and the given sign. + \c sz is an overapproximation of the size of the number stored at \c src. */ - template - void set(mpz & a, int sign, unsigned sz); + void set(mpz_cell& src, mpz & a, int sign, unsigned sz); static int64_t i64(mpz const & a) { return static_cast(a.m_val); } @@ -198,19 +237,19 @@ class mpz_manager { void set_i64(mpz & c, int64_t v) { if (v >= INT_MIN && v <= INT_MAX) { - del(c); c.m_val = static_cast(v); + c.m_kind = mpz_small; } else { - MPZ_BEGIN_CRITICAL(); set_big_i64(c, v); - MPZ_END_CRITICAL(); } } void set_big_ui64(mpz & c, uint64_t v); + #ifndef _MP_GMP + static unsigned capacity(mpz const & c) { return c.m_ptr->m_capacity; } static unsigned size(mpz const & c) { return c.m_ptr->m_size; } @@ -241,16 +280,28 @@ class mpz_manager { return ((static_cast(digits(a)[1]) << 32) | (static_cast(digits(a)[0]))); } - template - void get_sign_cell(mpz const & a, int & sign, mpz_cell * & cell) { + class sign_cell { + static const unsigned capacity = 2; + unsigned char m_bytes[sizeof(mpz_cell) + sizeof(digit_t) * capacity]; + mpz m_local; + mpz const& m_a; + int m_sign; + mpz_cell* m_cell; + public: + sign_cell(mpz_manager& m, mpz const& a); + int sign() { return m_sign; } + mpz_cell const* cell() { return m_cell; } + }; + + void get_sign_cell(mpz const & a, int & sign, mpz_cell * & cell, mpz_cell* reserve) { if (is_small(a)) { if (a.m_val == INT_MIN) { sign = -1; cell = m_int_min.m_ptr; } else { - cell = m_arg[IDX]; - SASSERT(cell->m_size == 1); + cell = reserve; + cell->m_size = 1; if (a.m_val < 0) { sign = -1; cell->m_digits[0] = -a.m_val; @@ -266,26 +317,29 @@ class mpz_manager { cell = a.m_ptr; } } + #else // GMP code - template - void get_arg(mpz const & a, mpz_t * & result) { - if (is_small(a)) { - result = m_arg[IDX]; - mpz_set_si(*result, a.m_val); - } - else { - result = a.m_ptr; - } - } + class ensure_mpz_t { + mpz_t m_local; + mpz_t* m_result; + public: + ensure_mpz_t(mpz const& a); + ~ensure_mpz_t(); + mpz_t& operator()() { return *m_result; } + }; void mk_big(mpz & a) { - if (a.m_ptr == 0) { + if (a.m_ptr == null) { a.m_val = 0; a.m_ptr = allocate(); + a.m_owner = mpz_self; } + a.m_kind = mpz_ptr; } + + #endif #ifndef _MP_GMP @@ -333,245 +387,49 @@ class mpz_manager { ~mpz_manager(); - static bool is_small(mpz const & a) { return a.m_ptr == nullptr; } + static bool is_small(mpz const & a) { return a.m_kind == mpz_small; } static mpz mk_z(int val) { return mpz(val); } - void del(mpz & a) { - if (a.m_ptr != nullptr) { - MPZ_BEGIN_CRITICAL(); - deallocate(a.m_ptr); - MPZ_END_CRITICAL(); - a.m_ptr = nullptr; - } - } + void del(mpz & a); - void add(mpz const & a, mpz const & b, mpz & c) { - STRACE("mpz", tout << "[mpz] " << to_string(a) << " + " << to_string(b) << " == ";); - if (is_small(a) && is_small(b)) { - set_i64(c, i64(a) + i64(b)); - } - else { - MPZ_BEGIN_CRITICAL(); - big_add(a, b, c); - MPZ_END_CRITICAL(); - } - STRACE("mpz", tout << to_string(c) << "\n";); - } + void add(mpz const & a, mpz const & b, mpz & c); - void sub(mpz const & a, mpz const & b, mpz & c) { - STRACE("mpz", tout << "[mpz] " << to_string(a) << " - " << to_string(b) << " == ";); - if (is_small(a) && is_small(b)) { - set_i64(c, i64(a) - i64(b)); - } - else { - MPZ_BEGIN_CRITICAL(); - big_sub(a, b, c); - MPZ_END_CRITICAL(); - } - STRACE("mpz", tout << to_string(c) << "\n";); - } + void sub(mpz const & a, mpz const & b, mpz & c); void inc(mpz & a) { add(a, mpz(1), a); } void dec(mpz & a) { add(a, mpz(-1), a); } - void mul(mpz const & a, mpz const & b, mpz & c) { - STRACE("mpz", tout << "[mpz] " << to_string(a) << " * " << to_string(b) << " == ";); - if (is_small(a) && is_small(b)) { - set_i64(c, i64(a) * i64(b)); - } - else { - MPZ_BEGIN_CRITICAL(); - big_mul(a, b, c); - MPZ_END_CRITICAL(); - } - STRACE("mpz", tout << to_string(c) << "\n";); - } + void mul(mpz const & a, mpz const & b, mpz & c); // d <- a + b*c - void addmul(mpz const & a, mpz const & b, mpz const & c, mpz & d) { - if (is_one(b)) { - add(a, c, d); - } - else if (is_minus_one(b)) { - sub(a, c, d); - } - else { - mpz tmp; - mul(b,c,tmp); - add(a,tmp,d); - del(tmp); - } - } - + void addmul(mpz const & a, mpz const & b, mpz const & c, mpz & d); // d <- a - b*c - void submul(mpz const & a, mpz const & b, mpz const & c, mpz & d) { - if (is_one(b)) { - sub(a, c, d); - } - else if (is_minus_one(b)) { - add(a, c, d); - } - else { - mpz tmp; - mul(b,c,tmp); - sub(a,tmp,d); - del(tmp); - } - } + void submul(mpz const & a, mpz const & b, mpz const & c, mpz & d); + void machine_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r); - void machine_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r) { - STRACE("mpz", tout << "[mpz-ext] divrem(" << to_string(a) << ", " << to_string(b) << ") == ";); - if (is_small(a) && is_small(b)) { - int64_t _a = i64(a); - int64_t _b = i64(b); - set_i64(q, _a / _b); - set_i64(r, _a % _b); - } - else { - MPZ_BEGIN_CRITICAL(); - big_div_rem(a, b, q, r); - MPZ_END_CRITICAL(); - } - STRACE("mpz", tout << "(" << to_string(q) << ", " << to_string(r) << ")\n";); - } + void machine_div(mpz const & a, mpz const & b, mpz & c); - void machine_div(mpz const & a, mpz const & b, mpz & c) { - STRACE("mpz", tout << "[mpz-ext] machine-div(" << to_string(a) << ", " << to_string(b) << ") == ";); - if (is_small(a) && is_small(b)) { - set_i64(c, i64(a) / i64(b)); - } - else { - MPZ_BEGIN_CRITICAL(); - big_div(a, b, c); - MPZ_END_CRITICAL(); - } - STRACE("mpz", tout << to_string(c) << "\n";); - } + void rem(mpz const & a, mpz const & b, mpz & c); - void rem(mpz const & a, mpz const & b, mpz & c) { - STRACE("mpz", tout << "[mpz-ext] rem(" << to_string(a) << ", " << to_string(b) << ") == ";); - if (is_small(a) && is_small(b)) { - set_i64(c, i64(a) % i64(b)); - } - else { - MPZ_BEGIN_CRITICAL(); - big_rem(a, b, c); - MPZ_END_CRITICAL(); - } - STRACE("mpz", tout << to_string(c) << "\n";); - } - - void div(mpz const & a, mpz const & b, mpz & c) { - STRACE("mpz", tout << "[mpz-ext] div(" << to_string(a) << ", " << to_string(b) << ") == ";); - if (is_neg(a)) { - mpz tmp; - machine_div_rem(a, b, c, tmp); - if (!is_zero(tmp)) { - if (is_neg(b)) - add(c, mk_z(1), c); - else - sub(c, mk_z(1), c); - } - del(tmp); - } - else { - machine_div(a, b, c); - } - STRACE("mpz", tout << to_string(c) << "\n";); - } + void div_gcd(mpz const & a, mpz const & b, mpz & c); - void mod(mpz const & a, mpz const & b, mpz & c) { - STRACE("mpz", tout << "[mpz-ext] mod(" << to_string(a) << ", " << to_string(b) << ") == ";); - rem(a, b, c); - if (is_neg(c)) { - if (is_pos(b)) - add(c, b, c); - else - sub(c, b, c); - } - STRACE("mpz", tout << to_string(c) << "\n";); - } + void div(mpz const & a, mpz const & b, mpz & c); - void neg(mpz & a) { - STRACE("mpz", tout << "[mpz] 0 - " << to_string(a) << " == ";); - if (is_small(a) && a.m_val == INT_MIN) { - // neg(INT_MIN) is not a small int - MPZ_BEGIN_CRITICAL(); - set_big_i64(a, - static_cast(INT_MIN)); - MPZ_END_CRITICAL(); - return; - } -#ifndef _MP_GMP - a.m_val = -a.m_val; -#else - if (is_small(a)) { - a.m_val = -a.m_val; - } - else { - mpz_neg(*a.m_ptr, *a.m_ptr); - } -#endif - STRACE("mpz", tout << to_string(a) << "\n";); - } + void mod(mpz const & a, mpz const & b, mpz & c); - void abs(mpz & a) { - if (is_small(a)) { - if (a.m_val < 0) { - if (a.m_val == INT_MIN) { - // abs(INT_MIN) is not a small int - MPZ_BEGIN_CRITICAL(); - set_big_i64(a, - static_cast(INT_MIN)); - MPZ_END_CRITICAL(); - } - else - a.m_val = -a.m_val; - } - } - else { -#ifndef _MP_GMP - a.m_val = 1; -#else - mpz_abs(*a.m_ptr, *a.m_ptr); -#endif - } - } + void neg(mpz & a); - static bool is_pos(mpz const & a) { -#ifndef _MP_GMP - return a.m_val > 0; -#else - if (is_small(a)) - return a.m_val > 0; - else - return mpz_sgn(*a.m_ptr) > 0; -#endif - } + void abs(mpz & a); - static bool is_neg(mpz const & a) { -#ifndef _MP_GMP - return a.m_val < 0; -#else - if (is_small(a)) - return a.m_val < 0; - else - return mpz_sgn(*a.m_ptr) < 0; -#endif - } + static bool is_pos(mpz const & a) { return sign(a) > 0; } + + static bool is_neg(mpz const & a) { return sign(a) < 0; } - static bool is_zero(mpz const & a) { -#ifndef _MP_GMP - return a.m_val == 0; -#else - if (is_small(a)) - return a.m_val == 0; - else - return mpz_sgn(*a.m_ptr) == 0; -#endif - } + static bool is_zero(mpz const & a) { return sign(a) == 0; } static int sign(mpz const & a) { #ifndef _MP_GMP @@ -593,10 +451,7 @@ class mpz_manager { return a.m_val == b.m_val; } else { - MPZ_BEGIN_CRITICAL(); - bool res = big_compare(a, b) == 0; - MPZ_END_CRITICAL(); - return res; + return big_compare(a, b) == 0; } } @@ -614,10 +469,7 @@ class mpz_manager { return a.m_val < b.m_val; } else { - MPZ_BEGIN_CRITICAL(); - bool res = big_compare(a, b) < 0; - MPZ_END_CRITICAL(); - return res; + return big_compare(a, b) < 0; } } @@ -661,25 +513,24 @@ class mpz_manager { void set(mpz & target, mpz const & source) { if (is_small(source)) { - del(target); target.m_val = source.m_val; + target.m_kind = mpz_small; } else { - MPZ_BEGIN_CRITICAL(); big_set(target, source); - MPZ_END_CRITICAL(); } } void set(mpz & target, mpz && source) { - del(target); target.m_val = source.m_val; std::swap(target.m_ptr, source.m_ptr); + auto o = target.m_owner; target.m_owner = source.m_owner; source.m_owner = o; + auto k = target.m_kind; target.m_kind = source.m_kind; source.m_kind = k; } void set(mpz & a, int val) { - del(a); a.m_val = val; + a.m_kind = mpz_small; } void set(mpz & a, unsigned val) { @@ -697,17 +548,15 @@ class mpz_manager { void set(mpz & a, uint64_t val) { if (val < INT_MAX) { - del(a); a.m_val = static_cast(val); + a.m_kind = mpz_small; } else { - MPZ_BEGIN_CRITICAL(); set_big_ui64(a, val); - MPZ_END_CRITICAL(); } } - void set(mpz & target, unsigned sz, digit_t const * digits); + void set_digits(mpz & target, unsigned sz, digit_t const * digits); mpz dup(const mpz & source) { mpz temp; @@ -716,13 +565,15 @@ class mpz_manager { } void reset(mpz & a) { - del(a); a.m_val = 0; + a.m_kind = mpz_small; } void swap(mpz & a, mpz & b) { std::swap(a.m_val, b.m_val); std::swap(a.m_ptr, b.m_ptr); + auto o = a.m_owner; a.m_owner = b.m_owner; b.m_owner = o; + auto k = a.m_kind; a.m_kind = b.m_kind; b.m_kind = k; } bool is_uint64(mpz const & a) const;