diff --git a/src/math/lp/core_solver_pretty_printer_def.h b/src/math/lp/core_solver_pretty_printer_def.h index 27aa6c75d..b8048b241 100644 --- a/src/math/lp/core_solver_pretty_printer_def.h +++ b/src/math/lp/core_solver_pretty_printer_def.h @@ -279,9 +279,9 @@ template void core_solver_pretty_printer::print() print_row(i); } m_out << std::endl; - if (m_core_solver.inf_set().size()) { + if (m_core_solver.inf_heap().size()) { m_out << "inf columns: "; - print_vector(m_core_solver.inf_set(), m_out); + print_vector(m_core_solver.inf_heap(), m_out); m_out << std::endl; } } diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h index 550b6fe36..942c87b3b 100644 --- a/src/math/lp/lar_core_solver_def.h +++ b/src/math/lp/lar_core_solver_def.h @@ -85,8 +85,8 @@ unsigned lar_core_solver::get_number_of_non_ints() const { void lar_core_solver::solve() { TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";); lp_assert(m_r_solver.non_basic_columns_are_set_correctly()); - lp_assert(m_r_solver.inf_set_is_correct()); - TRACE("find_feas_stats", tout << "infeasibles = " << m_r_solver.inf_set_size() << ", int_infs = " << get_number_of_non_ints() << std::endl;); + lp_assert(m_r_solver.inf_heap_is_correct()); + TRACE("find_feas_stats", tout << "infeasibles = " << m_r_solver.inf_heap_size() << ", int_infs = " << get_number_of_non_ints() << std::endl;); if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) { m_r_solver.set_status(lp_status::OPTIMAL); TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";); @@ -117,11 +117,9 @@ void lar_core_solver::solve() { } lp_assert(r_basis_is_OK()); lp_assert(m_r_solver.non_basic_columns_are_set_correctly()); - lp_assert(m_r_solver.inf_set_is_correct()); - - TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";); -} - + lp_assert(m_r_solver.inf_heap_is_correct()); + TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";); } +} // namespace lp diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 62712914c..5f544dff1 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -244,6 +244,14 @@ namespace lp { set.erase(j); } + void lar_solver::clean_popped_elements_for_heap(unsigned n, lpvar_heap& heap) { + vector to_remove; + for (unsigned j : heap) + if (j >= n) + to_remove.push_back(j); + for (unsigned j : to_remove) + heap.erase(j); + } void lar_solver::pop(unsigned k) { @@ -271,7 +279,7 @@ namespace lp { unsigned m = A_r().row_count(); clean_popped_elements(m, m_rows_with_changed_bounds); - clean_inf_set_of_r_solver_after_pop(); + clean_inf_heap_of_r_solver_after_pop(); lp_assert( m_settings.simplex_strategy() == simplex_strategy_enum::undecided || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); @@ -328,7 +336,7 @@ namespace lp { void lar_solver::set_costs_to_zero(const lar_term& term) { auto& rslv = m_mpq_lar_core_solver.m_r_solver; - auto& jset = m_mpq_lar_core_solver.m_r_solver.inf_set(); // hijack this set that should be empty right now + auto& jset = m_mpq_lar_core_solver.m_r_solver.inf_heap(); // hijack this set that should be empty right now lp_assert(jset.empty()); for (lar_term::ival p : term) { @@ -674,9 +682,9 @@ namespace lp { void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) { if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { if (costs_are_used()) { - bool was_infeas = m_mpq_lar_core_solver.m_r_solver.inf_set_contains(j); + bool was_infeas = m_mpq_lar_core_solver.m_r_solver.inf_heap_contains(j); m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j); - if (was_infeas != m_mpq_lar_core_solver.m_r_solver.inf_set_contains(j)) + if (was_infeas != m_mpq_lar_core_solver.m_r_solver.inf_heap_contains(j)) m_basic_columns_with_changed_cost.insert(j); } else { @@ -1293,12 +1301,12 @@ namespace lp { lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); } - void lar_solver::clean_inf_set_of_r_solver_after_pop() { + void lar_solver::clean_inf_heap_of_r_solver_after_pop() { vector became_feas; - clean_popped_elements(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.inf_set()); + clean_popped_elements_for_heap(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.inf_heap()); std::unordered_set basic_columns_with_changed_cost; m_inf_index_copy.reset(); - for (auto j : m_mpq_lar_core_solver.m_r_solver.inf_set()) + for (auto j : m_mpq_lar_core_solver.m_r_solver.inf_heap()) m_inf_index_copy.push_back(j); for (auto j : m_inf_index_copy) { if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { @@ -1316,16 +1324,16 @@ namespace lp { lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis_heading[j] < 0); m_mpq_lar_core_solver.m_r_solver.m_d[j] -= m_mpq_lar_core_solver.m_r_solver.m_costs[j]; m_mpq_lar_core_solver.m_r_solver.m_costs[j] = zero_of_type(); - m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j); + m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_heap(j); } became_feas.clear(); - for (unsigned j : m_mpq_lar_core_solver.m_r_solver.inf_set()) { + for (unsigned j : m_mpq_lar_core_solver.m_r_solver.inf_heap()) { lp_assert(m_mpq_lar_core_solver.m_r_heading[j] >= 0); if (m_mpq_lar_core_solver.m_r_solver.column_is_feasible(j)) became_feas.push_back(j); } for (unsigned j : became_feas) - m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j); + m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_heap(j); } @@ -1504,7 +1512,7 @@ namespace lp { m_mpq_lar_core_solver.m_r_x.resize(j + 1); m_mpq_lar_core_solver.m_r_lower_bounds.increase_size_by_one(); m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one(); - m_mpq_lar_core_solver.m_r_solver.inf_set_increase_size_by_one(); + m_mpq_lar_core_solver.m_r_solver.inf_heap_increase_size_by_one(); m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1); m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1); lp_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index 5b421d70d..f2017acd7 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -94,7 +94,7 @@ class lar_solver : public column_namer { // these are basic columns with the value changed, so the corresponding row in the tableau // does not sum to zero anymore u_set m_incorrect_columns; - // copy of m_r_solver.inf_set() + // copy of m_r_solver.inf_heap() unsigned_vector m_inf_index_copy; stacked_value m_term_count; vector m_terms; @@ -178,7 +178,9 @@ class lar_solver : public column_namer { bp); } + static void clean_popped_elements_for_heap(unsigned n, lpvar_heap& set); static void clean_popped_elements(unsigned n, u_set& set); + bool maximize_term_on_tableau(const lar_term & term, impq &term_max); bool costs_are_zeros_for_r_solver() const; @@ -230,7 +232,7 @@ class lar_solver : public column_namer { void remove_last_column_from_basis_tableau(unsigned j); void remove_last_column_from_tableau(); void pop_tableau(); - void clean_inf_set_of_r_solver_after_pop(); + void clean_inf_heap_of_r_solver_after_pop(); inline bool column_value_is_integer(unsigned j) const { return get_column_value(j).is_int(); } bool model_is_int_feasible() const; @@ -386,14 +388,15 @@ public: return tv::is_term(idx)? m_term_register.local_to_external(idx) : m_var_register.local_to_external(idx); } + bool column_corresponds_to_term(unsigned) const; inline unsigned row_count() const { return A_r().row_count(); } bool var_is_registered(var_index vj) const; - void clear_inf_set() { - m_mpq_lar_core_solver.m_r_solver.inf_set().clear(); + void clear_inf_heap() { + m_mpq_lar_core_solver.m_r_solver.inf_heap().clear(); } - inline void remove_column_from_inf_set(unsigned j) { - m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j); + inline void remove_column_from_inf_heap(unsigned j) { + m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_heap(j); } template void change_basic_columns_dependend_on_a_given_nb_column_report(unsigned j, diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp index f1ae95ea0..28f92d8e2 100644 --- a/src/math/lp/lp_core_solver_base.cpp +++ b/src/math/lp/lp_core_solver_base.cpp @@ -63,8 +63,8 @@ template bool lp::lp_core_solver_base >::colu template bool lp::lp_core_solver_base>::pivot_column_tableau(unsigned int, unsigned int); template bool lp::lp_core_solver_base::pivot_column_tableau(unsigned int, unsigned int); template void lp::lp_core_solver_base >::transpose_rows_tableau(unsigned int, unsigned int); -template bool lp::lp_core_solver_base >::inf_set_is_correct() const; -template bool lp::lp_core_solver_base::inf_set_is_correct() const; +template bool lp::lp_core_solver_base >::inf_heap_is_correct() const; +template bool lp::lp_core_solver_base::inf_heap_is_correct() const; template bool lp::lp_core_solver_base >::remove_from_basis(unsigned int); diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index fb0c28507..f38587573 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -28,9 +28,13 @@ Revision History: #include "math/lp/permutation_matrix.h" #include "math/lp/column_namer.h" #include "math/lp/u_set.h" - +#include "util/heap.h" namespace lp { +struct lpvar_lt { + bool operator()(lpvar v1, lpvar v2) const { return v1 < v2; } +}; +typedef heap lpvar_heap; template X dot_product(const vector & a, const vector & b) { lp_assert(a.size() == b.size()); @@ -51,24 +55,24 @@ private: public: bool current_x_is_feasible() const { TRACE("feas", - if (m_inf_set.size()) { - tout << "column " << *m_inf_set.begin() << " is infeasible" << std::endl; - print_column_info(*m_inf_set.begin(), tout); + if (m_inf_heap.size()) { + tout << "column " << *m_inf_heap.begin() << " is infeasible" << std::endl; + print_column_info(*m_inf_heap.begin(), tout); } else { tout << "x is feasible\n"; } ); - return m_inf_set.size() == 0; + return m_inf_heap.empty(); } - bool current_x_is_infeasible() const { return m_inf_set.size() != 0; } + bool current_x_is_infeasible() const { return m_inf_heap.size() != 0; } private: - u_set m_inf_set; + lpvar_heap m_inf_heap; public: - const u_set& inf_set() const { return m_inf_set; } - u_set& inf_set() { return m_inf_set; } - void inf_set_increase_size_by_one() { m_inf_set.increase_size_by_one(); } - bool inf_set_contains(unsigned j) const { return m_inf_set.contains(j); } - unsigned inf_set_size() const { return m_inf_set.size(); } + const lpvar_heap& inf_heap() const { return m_inf_heap; } + lpvar_heap& inf_heap() { return m_inf_heap; } + void inf_heap_increase_size_by_one() { m_inf_heap.reserve(m_inf_heap.size() + 1); } + bool inf_heap_contains(unsigned j) const { return m_inf_heap.contains(j); } + unsigned inf_heap_size() const { return m_inf_heap.size(); } indexed_vector m_pivot_row; // this is the real pivot row of the simplex tableu static_matrix & m_A; // the matrix A // vector const & m_b; // the right side @@ -255,7 +259,7 @@ public: bool calc_current_x_is_feasible_include_non_basis() const; - bool inf_set_is_correct() const; + bool inf_heap_is_correct() const; bool column_is_dual_feasible(unsigned j) const; @@ -526,8 +530,8 @@ public: swap(m_basis_heading, m_basis[i], m_basis[ii]); } - bool column_is_in_inf_set(unsigned j) const { - return m_inf_set.contains(j); + bool column_is_in_inf_heap(unsigned j) const { + return m_inf_heap.contains(j); } bool column_is_base(unsigned j) const { @@ -560,29 +564,26 @@ public: void track_column_feasibility(unsigned j) { if (column_is_feasible(j)) - remove_column_from_inf_set(j); + remove_column_from_inf_heap(j); else - insert_column_into_inf_set(j); + insert_column_into_inf_heap(j); } - void insert_column_into_inf_set(unsigned j) { + void insert_column_into_inf_heap(unsigned j) { TRACE("lar_solver", tout << "j = " << j << "\n";); - m_inf_set.insert(j); + if (!m_inf_heap.contains(j)) + m_inf_heap.insert(j); lp_assert(!column_is_feasible(j)); } - void remove_column_from_inf_set(unsigned j) { + void remove_column_from_inf_heap(unsigned j) { TRACE("lar_solver", tout << "j = " << j << "\n";); - m_inf_set.erase(j); + if (m_inf_heap.contains(j)) + m_inf_heap.erase(j); lp_assert(column_is_feasible(j)); } - void resize_inf_set(unsigned size) { + void clear_inf_heap() { TRACE("lar_solver",); - m_inf_set.resize(size); - } - - void clear_inf_set() { - TRACE("lar_solver",); - m_inf_set.clear(); + m_inf_heap.clear(); } bool costs_on_nbasis_are_zeros() const { diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h index 8619c926e..319966091 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -42,7 +42,7 @@ lp_core_solver_base(static_matrix & A, m_total_iterations(0), m_iters_with_no_cost_growing(0), m_status(lp_status::FEASIBLE), - m_inf_set(A.column_count()), + m_inf_heap(std::max(static_cast(1024), A.column_count())), m_pivot_row(A.column_count()), m_A(A), m_basis(basis), @@ -250,9 +250,9 @@ template bool lp_core_solver_base::calc_current_x return true; } -template bool lp_core_solver_base::inf_set_is_correct() const { +template bool lp_core_solver_base::inf_heap_is_correct() const { for (unsigned j = 0; j < this->m_n(); j++) { - bool belongs_to_set = m_inf_set.contains(j); + bool belongs_to_set = m_inf_heap.contains(j); bool is_feas = column_is_feasible(j); if (is_feas == belongs_to_set) { TRACE("lp_core", tout << "incorrectly set column in inf set "; print_column_info(j, tout) << "\n";); diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index 207428985..b2e402c08 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -33,6 +33,8 @@ Revision History: #include #include #include +#include "util/heap.h" + namespace lp { // This core solver solves (Ax=b, lower_bound_values \leq x \leq @@ -152,19 +154,28 @@ public: UNREACHABLE(); // unreachable return false; } - - unsigned get_number_of_basic_vars_that_might_become_inf( - unsigned j) const { // consider looking at the signs here: todo +/** + * Return the number of base non-free variables depending on the column j, + * different from bj, + * but take the min with the (bound+1). + * This function is used to select the pivot variable. +*/ + unsigned get_num_of_not_free_basic_dependent_vars( + unsigned j, unsigned bound, unsigned bj) const { // consider looking at the signs here: todo unsigned r = 0; for (const auto &cc : this->m_A.m_columns[j]) { - unsigned k = this->m_basis[cc.var()]; - if (this->m_column_types[k] != column_type::free_column) - r++; + unsigned basic_for_row = this->m_basis[cc.var()]; + if (basic_for_row == bj) + continue; + + // std::cout << this->m_A.m_rows[cc.var()] << std::endl; + if (this->m_column_types[basic_for_row] != column_type::free_column) + if (r++ > bound) return r; } return r; } - int find_beneficial_column_in_row_tableau_rows_bland_mode(int i, T &a_ent) { + int find_beneficial_entering_in_row_tableau_rows_bland_mode(int i, T &a_ent) { int j = -1; unsigned bj = this->m_basis[i]; bool bj_needs_to_grow = needs_to_grow(bj); @@ -190,15 +201,15 @@ public: return j; } - int find_beneficial_column_in_row_tableau_rows(int i, T &a_ent) { + int find_beneficial_entering_tableau_rows(int i, T &a_ent) { if (m_bland_mode_tableau) - return find_beneficial_column_in_row_tableau_rows_bland_mode(i, a_ent); + return find_beneficial_entering_in_row_tableau_rows_bland_mode(i, a_ent); // a short row produces short infeasibility explanation and benefits at // least one pivot operation int choice = -1; int nchoices = 0; - unsigned num_of_non_free_basics = 1000000; - unsigned len = 100000000; + unsigned min_non_free_so_far = -1; + unsigned best_col_sz = -1; unsigned bj = this->m_basis[i]; bool bj_needs_to_grow = needs_to_grow(bj); for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) { @@ -213,17 +224,18 @@ public: if (!monoid_can_increase(rc)) continue; } - unsigned damage = get_number_of_basic_vars_that_might_become_inf(j); - if (damage < num_of_non_free_basics) { - num_of_non_free_basics = damage; - len = this->m_A.m_columns[j].size(); + unsigned not_free = get_num_of_not_free_basic_dependent_vars(j, min_non_free_so_far, bj); + unsigned col_sz = this->m_A.m_columns[j].size(); + if (not_free < min_non_free_so_far || (not_free == min_non_free_so_far && col_sz < best_col_sz)) { + min_non_free_so_far = not_free; + best_col_sz = this->m_A.m_columns[j].size(); choice = k; nchoices = 1; - } else if (damage == num_of_non_free_basics && - this->m_A.m_columns[j].size() <= len && - (this->m_settings.random_next() % (++nchoices))) { - choice = k; - len = this->m_A.m_columns[j].size(); + } else if (not_free == min_non_free_so_far && + col_sz == best_col_sz) { + if (this->m_settings.random_next(++nchoices) == 0){ + choice = k; + } } } @@ -282,7 +294,7 @@ public: unsigned j, vector &leavings, T m, X &t, T &abs_of_d_of_leaving); - X get_max_bound(vector &b); + #ifdef Z3DEBUG void check_Ax_equal_b(); @@ -291,14 +303,6 @@ public: void check_correctness(); #endif - // from page 183 of Istvan Maros's book - // the basis structures have not changed yet - void update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving); - - // return 0 if the reduced cost at entering is close enough to the refreshed - // 1 if it is way off, and 2 if it is unprofitable - int refresh_reduced_cost_at_entering_and_check_that_it_is_off( - unsigned entering); void backup_and_normalize_costs(); @@ -356,15 +360,13 @@ public: } int find_smallest_inf_column() { - int j = -1; - for (unsigned k : this->inf_set()) { - if (k < static_cast(j)) { - j = k; - } - } - return j; + if (this->inf_heap().empty()) + return -1; + return this->inf_heap().min_value(); } + + const X &get_val_for_leaving(unsigned j) const { lp_assert(!this->column_is_feasible(j)); switch (this->m_column_types[j]) { @@ -405,7 +407,7 @@ public: } } T a_ent; - int entering = find_beneficial_column_in_row_tableau_rows( + int entering = find_beneficial_entering_tableau_rows( this->m_basis_heading[leaving], a_ent); if (entering == -1) { this->set_status(lp_status::INFEASIBLE); @@ -414,7 +416,7 @@ public: const X &new_val_for_leaving = get_val_for_leaving(leaving); X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent; this->m_x[leaving] = new_val_for_leaving; - this->remove_column_from_inf_set(leaving); + this->remove_column_from_inf_heap(leaving); advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta); if (this->current_x_is_feasible()) this->set_status(lp_status::OPTIMAL); diff --git a/src/math/lp/lp_primal_core_solver_tableau_def.h b/src/math/lp/lp_primal_core_solver_tableau_def.h index 898abb152..0c4f6216e 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -30,7 +30,7 @@ template void lp_primal_core_solver::one_iteratio else { advance_on_entering_tableau(entering); } - lp_assert(this->inf_set_is_correct()); + lp_assert(this->inf_heap_is_correct()); } template void lp_primal_core_solver::advance_on_entering_tableau(int entering) { @@ -256,7 +256,7 @@ template void lp_primal_core_solver::init_run_tab this->m_basis_sort_counter = 0; // to initiate the sort of the basis // this->set_total_iterations(0); this->iters_with_no_cost_growing() = 0; - lp_assert(this->inf_set_is_correct()); + lp_assert(this->inf_heap_is_correct()); if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) return; if (this->m_settings.backup_costs) diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h index c213333e0..727bc3531 100644 --- a/src/math/lp/lp_settings.h +++ b/src/math/lp/lp_settings.h @@ -218,6 +218,8 @@ public: unsigned hnf_cut_period() const { return m_hnf_cut_period; } void set_hnf_cut_period(unsigned period) { m_hnf_cut_period = period; } unsigned random_next() { return m_rand(); } + unsigned random_next(unsigned u ) { return m_rand(u); } + void set_random_seed(unsigned s) { m_rand.set_seed(s); } bool bound_progation() const { diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index 4d1cc6edb..4380b940b 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -1436,7 +1436,7 @@ void core::patch_monomials() { } else { m_lar_solver.pop(); restore_tableau(); - m_lar_solver.clear_inf_set(); + m_lar_solver.clear_inf_heap(); } SASSERT(m_lar_solver.ax_is_correct()); } diff --git a/src/math/lp/nla_tangent_lemmas.cpp b/src/math/lp/nla_tangent_lemmas.cpp index 299d8031f..4ba9eeccc 100644 --- a/src/math/lp/nla_tangent_lemmas.cpp +++ b/src/math/lp/nla_tangent_lemmas.cpp @@ -75,8 +75,8 @@ private: c().negate_relation(lemma, m_jy, m_y.rat_sign()*pl.y); #if Z3DEBUG SASSERT(c().val(m_x) == m_xy.x && c().val(m_y) == m_xy.y); - int mult_sign = nla::rat_sign(pl.x - m_xy.x)*nla::rat_sign(pl.y - m_xy.y); - SASSERT((mult_sign == 1) == m_below); + // int mult_sign = nla::rat_sign(pl.x - m_xy.x)*nla::rat_sign(pl.y - m_xy.y); + SASSERT((nla::rat_sign(pl.x - m_xy.x)*nla::rat_sign(pl.y - m_xy.y) == 1) == m_below); // If "mult_sign is 1" then (a - x)(b-y) > 0 and ab - bx - ay + xy > 0 // or -ab + bx + ay < xy or -ay - bx + xy > -ab // val(j) stands for xy. So, finally we have -ay - bx + j > - ab