diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h
index 52c5b6f1a..966e2f0d7 100644
--- a/src/math/lp/lar_core_solver.h
+++ b/src/math/lp/lar_core_solver.h
@@ -20,9 +20,6 @@ Author:
 namespace lp {
 
 class lar_core_solver  {
-    // m_sign_of_entering is set to 1 if the entering variable needs
-    // to grow and is set to -1  otherwise
-    int m_sign_of_entering_delta;
     vector<std::pair<mpq, unsigned>> m_infeasible_linear_combination;
     int m_infeasible_sum_sign; // todo: get rid of this field
     vector<numeric_pair<mpq>> m_right_sides_dummy;
@@ -40,8 +37,6 @@ public:
     vector<unsigned>         m_r_basis;
     vector<unsigned>         m_r_nbasis;
     vector<int>              m_r_heading;
-    stacked_vector<unsigned> m_r_columns_nz;
-    stacked_vector<unsigned> m_r_rows_nz;
     
 
     lp_primal_core_solver<mpq, numeric_pair<mpq>> m_r_solver; // solver in rational numbers
@@ -82,16 +77,9 @@ public:
         m_r_solver.print_column_bound_info(m_r_solver.m_basis[row_index], out);        
     }
     
-    bool row_is_infeasible(unsigned row);
-
-    bool row_is_evidence(unsigned row);
-
-    bool find_evidence_row();
-
+    
     void prefix_r();
 
-    void prefix_d();
-
     unsigned m_m() const { return m_r_A.row_count();  }
 
     unsigned m_n() const { return m_r_A.column_count(); }
@@ -103,8 +91,6 @@ public:
     template <typename L>
     int get_sign(const L & v) { return v > zero_of_type<L>() ? 1 : (v < zero_of_type<L>() ? -1 : 0); }
 
-    void fill_evidence(unsigned row);
-
     unsigned get_number_of_non_ints() const;
 
     void solve();
@@ -131,31 +117,6 @@ public:
         
     }
 
-    template <typename K> 
-    void push_vector(stacked_vector<K> & pushed_vector, const vector<K> & vector) {
-        lp_assert(pushed_vector.size() <= vector.size());
-        for (unsigned i = 0; i < vector.size();i++) {
-            if (i == pushed_vector.size()) {
-                pushed_vector.push_back(vector[i]);
-            } else {
-                pushed_vector[i] = vector[i];
-            }
-        }
-        pushed_vector.push();
-    }
-
-    void pop_markowitz_counts(unsigned k) {
-        m_r_columns_nz.pop(k);
-        m_r_rows_nz.pop(k);
-        m_r_solver.m_columns_nz.resize(m_r_columns_nz.size());
-        m_r_solver.m_rows_nz.resize(m_r_rows_nz.size());
-        for (unsigned i = 0; i < m_r_columns_nz.size(); i++)
-            m_r_solver.m_columns_nz[i] = m_r_columns_nz[i];
-        for (unsigned i = 0; i < m_r_rows_nz.size(); i++)
-            m_r_solver.m_rows_nz[i] = m_r_rows_nz[i];
-    }
-
-    
     void pop(unsigned k) {
         // rationals
         m_r_lower_bounds.pop(k);
@@ -172,72 +133,7 @@ public:
     }
 
     
-    template <typename L>
-    bool is_zero_vector(const vector<L> & b) {
-        for (const L & m: b)
-            if (!is_zero(m)) return false;
-        return true;
-    }
-
-
-    bool update_xj_and_get_delta(unsigned j, non_basic_column_value_position pos_type, numeric_pair<mpq> & delta) {
-        auto & x = m_r_x[j];
-        switch (pos_type) {
-        case at_lower_bound:
-            if (x == m_r_solver.m_lower_bounds[j])
-                return false;
-            delta = m_r_solver.m_lower_bounds[j] - x;
-            m_r_solver.m_x[j] = m_r_solver.m_lower_bounds[j];
-            break;
-        case at_fixed:
-        case at_upper_bound:
-            if (x == m_r_solver.m_upper_bounds[j])
-                return false;
-            delta = m_r_solver.m_upper_bounds[j] - x;
-            x = m_r_solver.m_upper_bounds[j];
-            break;
-        case free_of_bounds: {
-            return false;
-        }
-        case not_at_bound:
-            switch (m_column_types[j]) {
-            case column_type::free_column:
-                return false;
-            case column_type::upper_bound:
-                delta = m_r_solver.m_upper_bounds[j] - x;
-                x = m_r_solver.m_upper_bounds[j];
-                break;
-            case column_type::lower_bound:
-                delta = m_r_solver.m_lower_bounds[j] - x;
-                x = m_r_solver.m_lower_bounds[j];
-                break;
-            case column_type::boxed:
-                if (x > m_r_solver.m_upper_bounds[j]) {
-                    delta = m_r_solver.m_upper_bounds[j] - x;
-                    x += m_r_solver.m_upper_bounds[j];
-                } else {
-                    delta = m_r_solver.m_lower_bounds[j] - x;
-                    x = m_r_solver.m_lower_bounds[j];
-                }
-                break;
-            case column_type::fixed:
-                delta = m_r_solver.m_lower_bounds[j] - x;
-                x = m_r_solver.m_lower_bounds[j];
-                break;
-
-            default:
-                lp_assert(false);
-            }
-            break;
-        default:
-            lp_unreachable();
-        }
-        m_r_solver.remove_column_from_inf_set(j);
-        return true;
-    }
-
     
-
     bool r_basis_is_OK() const {
 #ifdef Z3DEBUG
         
@@ -255,28 +151,6 @@ public:
     }
     
    
-    
-    void fill_basis_d(
-                      vector<unsigned>& basis_d,
-                      vector<int>& heading_d,
-                      vector<unsigned>& nbasis_d){
-        basis_d = m_r_basis;
-        heading_d = m_r_heading;
-        nbasis_d = m_r_nbasis;
-    }
-
-    template <typename L, typename K>
-    void extract_signature_from_lp_core_solver(const lp_primal_core_solver<L, K> & solver, lar_solution_signature & signature) {
-        signature.clear();
-        lp_assert(signature.size() == 0);
-        for (unsigned j = 0; j < solver.m_basis_heading.size(); j++) {
-            if (solver.m_basis_heading[j] < 0) {
-                signature[j] = solver.get_non_basic_column_value_position(j);
-            }
-        }
-    }
-
-    
     bool lower_bound_is_set(unsigned j) const {
         switch (m_column_types[j]) {
         case column_type::free_column:
diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h
index 419345f0c..2e205291e 100644
--- a/src/math/lp/lar_core_solver_def.h
+++ b/src/math/lp/lar_core_solver_def.h
@@ -46,7 +46,6 @@ void lar_core_solver::prefix_r() {
 }
 
 
-
 void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() {
     unsigned bj = m_r_basis[m_r_solver.m_inf_row_index_for_tableau];
     m_infeasible_sum_sign =  m_r_solver.inf_sign_of_column(bj);
diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h
index 96f8729a0..9fb4d4333 100644
--- a/src/math/lp/lp_core_solver_base.h
+++ b/src/math/lp/lp_core_solver_base.h
@@ -72,8 +72,6 @@ public:
     unsigned inf_set_size() const { return m_inf_set.size(); }
     bool using_infeas_costs() const { return m_using_infeas_costs; }
     void set_using_infeas_costs(bool val)  { m_using_infeas_costs = val; }
-    vector<unsigned>      m_columns_nz; // m_columns_nz[i] keeps an approximate value of non zeroes the i-th column
-    vector<unsigned>      m_rows_nz; // m_rows_nz[i] keeps an approximate value of non zeroes in the i-th row
     indexed_vector<T>     m_pivot_row; // this is the real pivot row of the simplex tableu
     static_matrix<T, X> & m_A; // the matrix A
     // vector<X> const &           m_b; // the right side
diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h
index fecf66e68..17b1ea494 100644
--- a/src/math/lp/lp_primal_core_solver.h
+++ b/src/math/lp/lp_primal_core_solver.h
@@ -41,54 +41,21 @@ namespace lp {
 template <typename T, typename X>
 class lp_primal_core_solver:public lp_core_solver_base<T, X> {
 public:
-    // m_sign_of_entering is set to 1 if the entering variable needs
-    // to grow and is set to -1  otherwise
-    unsigned       m_column_norm_update_counter;
-    T              m_enter_price_eps;
     int            m_sign_of_entering_delta;
-    indexed_vector<T> m_beta; // see Swietanowski working vector beta for column norms
-    T                 m_epsilon_of_reduced_cost;
     vector<T>         m_costs_backup;
     unsigned          m_inf_row_index_for_tableau;
     bool              m_bland_mode_tableau;
-    u_set           m_left_basis_tableau;
+    u_set             m_left_basis_tableau;
     unsigned          m_bland_mode_threshold;
     unsigned          m_left_basis_repeated;
     vector<unsigned>  m_leaving_candidates;
     
     std::list<unsigned> m_non_basis_list;
     void sort_non_basis();
-    void sort_non_basis_rational();
     int choose_entering_column(unsigned number_of_benefitial_columns_to_go_over);
     int choose_entering_column_tableau();
     int choose_entering_column_presize(unsigned number_of_benefitial_columns_to_go_over);
     
-    bool column_is_benefitial_for_entering_basis_on_sign_row_strategy(unsigned j, int sign) const {
-        // sign = 1 means the x of the basis column of the row has to grow to become feasible, when the coeff before j is neg, or x - has to diminish when the coeff is pos
-        // we have xbj = -aj * xj
-        lp_assert(this->m_basis_heading[j] < 0);
-        lp_assert(this->column_is_feasible(j));
-        switch (this->m_column_types[j]) {
-        case column_type::free_column: return true;
-        case column_type::fixed: return false;
-        case column_type::lower_bound:
-            if (sign < 0)
-                return true;
-            return !this->x_is_at_lower_bound(j);
-        case column_type::upper_bound:
-            if (sign > 0)
-                return true;
-            return !this->x_is_at_upper_bound(j);
-        case column_type::boxed:
-            if (sign < 0)
-                return !this->x_is_at_lower_bound(j);
-            return !this->x_is_at_upper_bound(j);
-        }
-
-        lp_assert(false); // cannot be here
-        return false;
-    }
-    
 
     bool needs_to_grow(unsigned bj) const {
         lp_assert(!this->column_is_feasible(bj));
@@ -272,10 +239,7 @@ public:
         a_ent = rc.coeff();
         return rc.var();
     }
-    static X positive_infinity() {
-        return convert_struct<X, unsigned>::convert(std::numeric_limits<unsigned>::max());
-    }
-
+    
     bool try_jump_to_another_bound_on_entering(unsigned entering, const X & theta, X & t, bool & unlimited);
     bool try_jump_to_another_bound_on_entering_unlimited(unsigned entering, X & t);
     int find_leaving_and_t_tableau(unsigned entering, X & t);
@@ -313,8 +277,6 @@ public:
 
     void get_bound_on_variable_and_update_leaving_precisely(unsigned j, vector<unsigned> & leavings, T m, X & t, T & abs_of_d_of_leaving);
 
-    vector<T> m_lower_bounds_dummy; // needed for the base class only
-
     X get_max_bound(vector<X> & b);
 
 #ifdef Z3DEBUG
@@ -334,10 +296,6 @@ public:
 
     void backup_and_normalize_costs();
 
-    void init_run();
-
-    void calc_working_vector_beta_for_column_norms();
-
     void advance_on_entering_and_leaving(int entering, int leaving, X & t);
     void advance_on_entering_and_leaving_tableau(int entering, int leaving, X & t);
     void advance_on_entering_equal_leaving(int entering, X & t);
@@ -368,7 +326,6 @@ public:
 
     // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;}
 
-    void one_iteration();
     void one_iteration_tableau();
 
     // this version assumes that the leaving already has the right value, and does not update it
@@ -658,12 +615,6 @@ public:
     void init_infeasibility_costs_for_changed_basis_only();
 
     void print_column(unsigned j, std::ostream & out);
-    // j is the basic column, x is the value at x[j]
-    // d is the coefficient before m_entering in the row with j as the basis column
-    template <typename L>
-    bool same_sign_with_entering_delta(const L & a) {
-        return (a > zero_of_type<L>() && m_sign_of_entering_delta > 0) || (a < zero_of_type<L>() && m_sign_of_entering_delta < 0);
-    }
     
     void print_bound_info_and_x(unsigned j, std::ostream & out);
     
@@ -730,8 +681,6 @@ public:
                                   column_type_array,
                                   lower_bound_values,
                                   upper_bound_values),
-        m_beta(A.row_count()),
-        m_epsilon_of_reduced_cost(T(1)/T(10000000)),
         m_bland_mode_threshold(1000) {
         this->set_status(lp_status::UNKNOWN);
     }
diff --git a/src/math/lp/lp_primal_core_solver_def.h b/src/math/lp/lp_primal_core_solver_def.h
index 04c32d41c..1c48f1163 100644
--- a/src/math/lp/lp_primal_core_solver_def.h
+++ b/src/math/lp/lp_primal_core_solver_def.h
@@ -32,7 +32,7 @@ namespace lp {
 // The right side b is given implicitly by x and the basis
 
 template <typename T, typename X>
-void lp_primal_core_solver<T, X>::sort_non_basis_rational() {
+void lp_primal_core_solver<T, X>::sort_non_basis() {
      std::sort(this->m_nbasis.begin(), this->m_nbasis.end(), [this](unsigned a, unsigned b) {
                 unsigned ca = this->m_A.number_of_non_zeroes_in_column(a);
                 unsigned cb = this->m_A.number_of_non_zeroes_in_column(b);
@@ -50,10 +50,6 @@ void lp_primal_core_solver<T, X>::sort_non_basis_rational() {
 }
 
 
-template <typename T, typename X>
-void lp_primal_core_solver<T, X>::sort_non_basis() {
-     sort_non_basis_rational();   
-}
 
 template <typename T, typename X>
 bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_basis(unsigned j) const {
@@ -249,15 +245,6 @@ lp_primal_core_solver<T, X>::get_bound_on_variable_and_update_leaving_precisely(
     }
 }
 
-template <typename T, typename X>    X lp_primal_core_solver<T, X>::get_max_bound(vector<X> & b) {
-    X ret = zero_of_type<X>();
-    for (auto & v : b) {
-        X a = abs(v);
-        if (a > ret) ret = a;
-    }
-    return ret;
-}
-
 #ifdef Z3DEBUG
 template <typename T, typename X>   void lp_primal_core_solver<T, X>::check_Ax_equal_b() {
     dense_matrix<T, X> d(this->m_A);
@@ -282,38 +269,6 @@ template <typename T, typename X>    void lp_primal_core_solver<T, X>::check_cor
 }
 #endif
 
-// from page 183 of Istvan Maros's book
-// the basis structures have not changed yet
-template <typename T, typename X>
-void lp_primal_core_solver<T, X>::update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving) {
-    // the basis heading has changed already
-#ifdef Z3DEBUG
-    auto & basis_heading = this->m_basis_heading;
-    lp_assert(basis_heading[entering] >= 0 && static_cast<unsigned>(basis_heading[entering]) < this->m_m());
-    lp_assert(basis_heading[leaving] < 0);
-#endif
-    T pivot = this->m_pivot_row[entering];
-    T dq = this->m_d[entering]/pivot;
-    for (auto j : this->m_pivot_row.m_index) {
-        //            for (auto j : this->m_nbasis)
-        if (this->m_basis_heading[j] >= 0) continue;
-        if (j != leaving)
-            this->m_d[j] -= dq * this->m_pivot_row[j];
-    }
-    this->m_d[leaving] = -dq;
-    if (this->current_x_is_infeasible()) {
-        this->m_d[leaving] -= this->m_costs[leaving];
-        this->m_costs[leaving] = zero_of_type<T>();
-    }
-    this->m_d[entering] = numeric_traits<T>::zero();
-}
-
-// 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
-template <typename T, typename X>    int lp_primal_core_solver<T, X>::refresh_reduced_cost_at_entering_and_check_that_it_is_off(unsigned entering) {
-    return 0;
-    
-}
 
 template <typename T, typename X>    void lp_primal_core_solver<T, X>::backup_and_normalize_costs() {
     if (this->m_look_for_feasible_solution_only)
@@ -321,9 +276,6 @@ template <typename T, typename X>    void lp_primal_core_solver<T, X>::backup_an
     m_costs_backup = this->m_costs;    
 }
 
-template <typename T, typename X>    void lp_primal_core_solver<T, X>::init_run() {
-    
-}
 
 
 template <typename T, typename X>
@@ -377,20 +329,6 @@ template <typename T, typename X>    void lp_primal_core_solver<T, X>::find_feas
     solve();
 }
 
-template <typename T, typename X> void lp_primal_core_solver<T, X>::one_iteration() {
-    unsigned number_of_benefitial_columns_to_go_over = get_number_of_non_basic_column_to_try_for_enter();
-    int entering = choose_entering_column(number_of_benefitial_columns_to_go_over);
-    if (entering == -1) {
-        decide_on_status_when_cannot_find_entering();
-    }
-    else {
-        advance_on_entering(entering);
-    }
-}
-
-
-
-
 template <typename T, typename X>
 void lp_primal_core_solver<T, X>::init_infeasibility_costs_for_changed_basis_only() {
     for (unsigned i :  this->m_ed.m_index)
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 d4355fa80..de2f1c208 100644
--- a/src/math/lp/lp_primal_core_solver_tableau_def.h
+++ b/src/math/lp/lp_primal_core_solver_tableau_def.h
@@ -283,7 +283,6 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tab
             return;
         if (this->m_settings.backup_costs)
             backup_and_normalize_costs();
-        m_epsilon_of_reduced_cost = zero_of_type<T>();
         
         if (this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows)
             init_tableau_rows();
diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h
index 197b7c04f..314930fc1 100644
--- a/src/math/lp/lp_settings.h
+++ b/src/math/lp/lp_settings.h
@@ -179,15 +179,9 @@ public:
     bool int_run_gcd_test() const { return m_int_run_gcd_test; }
     bool& int_run_gcd_test() { return m_int_run_gcd_test; }
     unsigned      reps_in_scaler { 20 };
-    // when the absolute value of an element is less than pivot_epsilon
-    // in pivoting, we treat it as a zero
-    // a quotation "if some choice of the entering variable leads to an eta matrix
-    // whose diagonal element in the eta column is less than e2 (entering_diag_epsilon) in magnitude, the this choice is rejected ...
     int           c_partial_pivoting { 10 }; // this is the constant c from page 410
     unsigned      depth_of_rook_search { 4 };
     bool          using_partial_pivoting { true };
-    // dissertation of Achim Koberstein
-    // if Bx - b is different at any component more that refactor_epsilon then we refactor
     
     unsigned     percent_of_entering_to_check { 5 }; // we try to find a profitable column in a percentage of the columns
     bool         use_scaling { true };
diff --git a/src/math/lp/numeric_pair.h b/src/math/lp/numeric_pair.h
index f59aa84ba..4a60c82ca 100644
--- a/src/math/lp/numeric_pair.h
+++ b/src/math/lp/numeric_pair.h
@@ -107,7 +107,6 @@ public:
 template <typename X, typename Y>
 struct convert_struct {
     static X convert(const Y & y){ return X(y);}
-    static bool is_epsilon_small(const X & x,  const double & y) { return std::abs(numeric_traits<X>::get_double(x)) < y; }
     static bool below_bound_numeric(const X &, const X &, const Y &) { /*lp_unreachable();*/ return false;}
     static bool above_bound_numeric(const X &, const X &, const Y &) { /*lp_unreachable();*/ return false; }
 };
@@ -316,16 +315,12 @@ struct convert_struct<double, numeric_pair<double>> {
 
 typedef numeric_pair<mpq> impq;
 
-template <typename X> bool is_epsilon_small(const X & v, const double& eps);   // forward definition { return convert_struct<X, double>::is_epsilon_small(v, eps);}
 
 template <typename T>
 struct convert_struct<numeric_pair<T>, double> {
     static numeric_pair<T> convert(const double & q) {
         return numeric_pair<T>(convert_struct<T, double>::convert(q), numeric_traits<T>::zero());
     }
-    static bool is_epsilon_small(const numeric_pair<T> & p, const double & eps) {
-        return convert_struct<T, double>::is_epsilon_small(p.x, eps) && convert_struct<T, double>::is_epsilon_small(p.y, eps);
-    }
     static bool below_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
         // lp_unreachable();
         return false;
@@ -340,10 +335,7 @@ struct convert_struct<numeric_pair<double>, double> {
     static numeric_pair<double> convert(const double & q) {
         return numeric_pair<double>(q, 0.0);
     }
-    static bool is_epsilon_small(const numeric_pair<double> & p, const double & eps) {
-        return std::abs(p.x) < eps && std::abs(p.y) < eps;
-    }
-
+    
     static int compare_on_coord(const double & x, const double & bound, const double eps) {
         if (bound == 0) return (x < - eps)? -1: (x > eps? 1 : 0); // it is an important special case
         double relative = (bound > 0)? - eps: eps;
@@ -369,9 +361,6 @@ struct convert_struct<numeric_pair<double>, double> {
 
 template <>
 struct convert_struct<double, double> {
-    static bool is_epsilon_small(const double& x, const double & eps) {
-        return x < eps && x > -eps;
-    }
     static double convert(const double & y){ return y;}
     static bool below_bound_numeric(const double & x, const double & bound, const double & eps) {
         if (bound == 0) return x < - eps;
@@ -384,8 +373,6 @@ struct convert_struct<double, double> {
         return x > bound * (1.0 + relative) + eps;
     }
 };
-
-template <typename X> bool is_epsilon_small(const X & v, const double &eps) { return convert_struct<X, double>::is_epsilon_small(v, eps);}
 template <typename X> bool below_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::below_bound_numeric(x, bound, eps);}
 template <typename X> bool above_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::above_bound_numeric(x, bound, eps);}
 template  <typename T>  T floor(const numeric_pair<T> & r) {