diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index b5d64d345..9a1284b37 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -422,8 +422,9 @@ class theory_lra::imp {
         m_theory_var2var_index.setx(v, var, UINT_MAX);
         m_var_index2theory_var.setx(var, v, UINT_MAX);
         m_var_trail.push_back(v);
-        add_def_constraint(m_solver->add_var_bound(var, lp::GE, rational(c)));
-        add_def_constraint(m_solver->add_var_bound(var, lp::LE, rational(c)));
+        lp::explanation e;
+        add_def_constraint(m_solver->add_var_bound(var, lp::GE, rational(c), e));
+        add_def_constraint(m_solver->add_var_bound(var, lp::LE, rational(c), e));
         return var;
     }
 
@@ -799,7 +800,15 @@ class theory_lra::imp {
         m_definitions.setx(index, v, null_theory_var);
         ++m_stats.m_add_rows;
     }
-        
+
+    void process_conflict() {
+        NOT_IMPLEMENTED_YET();
+    }
+
+    bool is_infeasible() const {
+        return m_solver->get_status() == lp::lp_status::INFEASIBLE;
+    }
+    
     void internalize_eq(theory_var v1, theory_var v2) {  
         app_ref term(m.mk_fresh_const("eq", a.mk_real()), m);
         scoped_internalize_state st(*this);
@@ -809,8 +818,14 @@ class theory_lra::imp {
         st.coeffs().push_back(rational::minus_one());
         theory_var z = internalize_linearized_def(term, st);      
         lp::var_index vi = register_theory_var_in_lar_solver(z);
-        add_def_constraint(m_solver->add_var_bound(vi, lp::LE, rational::zero()));
-        add_def_constraint(m_solver->add_var_bound(vi, lp::GE, rational::zero()));
+        add_def_constraint(m_solver->add_var_bound(vi, lp::LE, rational::zero(), m_explanation));
+        if (is_infeasible()) {
+            process_conflict(); // exit here?
+        }
+        add_def_constraint(m_solver->add_var_bound(vi, lp::GE, rational::zero(), m_explanation));
+        if (is_infeasible()) {
+            process_conflict();
+        }
         TRACE("arith", 
               {
                   expr*  o1 = get_enode(v1)->get_owner();
@@ -1226,10 +1241,12 @@ public:
                                           get_enode(w)->get_owner())));
         theory_var z = internalize_def(term);
         lp::var_index vi = register_theory_var_in_lar_solver(z);
-        add_def_constraint(m_solver->add_var_bound(vi, lp::GE, rational::zero()));
-        add_def_constraint(m_solver->add_var_bound(vi, lp::LE, rational::zero()));
-        add_def_constraint(m_solver->add_var_bound(register_theory_var_in_lar_solver(v), lp::GE, rational::zero()));
-        add_def_constraint(m_solver->add_var_bound(register_theory_var_in_lar_solver(v), lp::LT, abs(r)));
+        lp::explanation e;
+        add_def_constraint(m_solver->add_var_bound(vi, lp::GE, rational::zero(), e));
+        add_def_constraint(m_solver->add_var_bound(vi, lp::LE, rational::zero(), e));
+        add_def_constraint(m_solver->add_var_bound(register_theory_var_in_lar_solver(v), lp::GE, rational::zero(), e));
+        add_def_constraint(m_solver->add_var_bound(register_theory_var_in_lar_solver(v), lp::LT, abs(r), e));
+        SASSERT(m_solver->get_status() != lp::lp_status::INFEASIBLE);
         TRACE("arith", m_solver->print_constraints(tout << term << "\n"););
     }
 
@@ -2971,12 +2988,13 @@ public:
         lp::constraint_index ci;
         if (is_int && !is_true) {
             rational bound = b.get_value(false).get_rational();
-            ci = m_solver->add_var_bound(vi, k, bound);
+            ci = m_solver->add_var_bound(vi, k, bound, m_explanation);
         }
         else {
-            ci = m_solver->add_var_bound(vi, k, b.get_value());
+            ci = m_solver->add_var_bound(vi, k, b.get_value(), m_explanation);
         }
         if (m_solver->get_status() == lp::lp_status::INFEASIBLE) {
+            NOT_IMPLEMENTED_YET();
             return;
         }
         TRACE("arith", tout << "v" << b.get_var() << "\n";);
diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp
index 31b767f1a..4ae037a96 100644
--- a/src/test/lp/lp.cpp
+++ b/src/test/lp/lp.cpp
@@ -2682,7 +2682,8 @@ void test_term() {
 
     vector<std::pair<mpq, var_index>> term_one;
     term_one.push_back(std::make_pair(mpq(1), one));
-    solver.add_constraint(term_one, lconstraint_kind::EQ, mpq(0));
+    explanation e;
+    solver.add_constraint(term_one, lconstraint_kind::EQ, mpq(0), e);
 
     vector<std::pair<mpq, var_index>> term_ls;
     term_ls.push_back(std::pair<mpq, var_index>(mpq(1), x));
@@ -2694,13 +2695,13 @@ void test_term() {
     ls.push_back(std::pair<mpq, var_index>(mpq(1), x));
     ls.push_back(std::pair<mpq, var_index>(mpq(1), y));
     ls.push_back(std::pair<mpq, var_index>(mpq(1), z));
-    
-    solver.add_constraint(ls, lconstraint_kind::EQ, mpq(0));
+
+    solver.add_constraint(ls, lconstraint_kind::EQ, mpq(0), e);
     ls.clear();
     ls.push_back(std::pair<mpq, var_index>(mpq(1), x));
-    solver.add_constraint(ls, lconstraint_kind::LT, mpq(0));
+    solver.add_constraint(ls, lconstraint_kind::LT, mpq(0), e);
     ls.push_back(std::pair<mpq, var_index>(mpq(2), y));
-    solver.add_constraint(ls, lconstraint_kind::GT, mpq(0));
+    solver.add_constraint(ls, lconstraint_kind::GT, mpq(0), e);
     auto status = solver.solve();
     std::cout << lp_status_to_string(status) << std::endl;
     std::unordered_map<var_index, mpq> model;
@@ -2717,6 +2718,7 @@ void test_term() {
 
 void test_evidence_for_total_inf_simple(argument_parser & args_parser) {
     lar_solver solver;
+    explanation e;
     var_index x = solver.add_var(0, false);
     var_index y = solver.add_var(1, false);
     solver.add_var_bound(x, LE, mpq(-1));
@@ -2725,10 +2727,11 @@ void test_evidence_for_total_inf_simple(argument_parser & args_parser) {
     
     ls.push_back(std::pair<mpq, var_index>(mpq(1), x));
     ls.push_back(std::pair<mpq, var_index>(mpq(1), y));
-    solver.add_constraint(ls, GE, mpq(1));
+
+    solver.add_constraint(ls, GE, mpq(1), e);
     ls.pop_back();
     ls.push_back(std::pair<mpq, var_index>(- mpq(1), y));
-    solver.add_constraint(ls, lconstraint_kind::GE, mpq(0));
+    solver.add_constraint(ls, lconstraint_kind::GE, mpq(0), e);
     auto status = solver.solve();
     std::cout << lp_status_to_string(status) << std::endl;
     std::unordered_map<var_index, mpq> model;
@@ -2772,7 +2775,7 @@ void test_bound_propagation_one_small_sample1() {
     coeffs.push_back(std::pair<mpq, var_index>(mpq(-1), c));
     ls.add_constraint(coeffs, LE, zero_of_type<mpq>());
     vector<implied_bound> ev;
-    ls.add_var_bound(a, LE, mpq(1));
+    ls.add_var_bound(a, LE, mpq(1), e);
     ls.solve();
     my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
@@ -2823,9 +2826,10 @@ void test_bound_propagation_one_row() {
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(mpq(1), x0));
     c.push_back(std::pair<mpq, var_index>(mpq(-1), x1));
-    ls.add_constraint(c, EQ, one_of_type<mpq>());
+    explanation e;
+    ls.add_constraint(c, EQ, one_of_type<mpq>(), e);
     vector<implied_bound> ev;
-    ls.add_var_bound(x0, LE, mpq(1));
+    ls.add_var_bound(x0, LE, mpq(1), e);
     ls.solve();
     my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
@@ -2837,11 +2841,12 @@ void test_bound_propagation_one_row_with_bounded_vars() {
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(mpq(1), x0));
     c.push_back(std::pair<mpq, var_index>(mpq(-1), x1));
-    ls.add_constraint(c, EQ, one_of_type<mpq>());
+    explanation e;
+    ls.add_constraint(c, EQ, one_of_type<mpq>(), e);
     vector<implied_bound> ev;
-    ls.add_var_bound(x0, GE, mpq(-3));
-    ls.add_var_bound(x0, LE, mpq(3));
-    ls.add_var_bound(x0, LE, mpq(1));
+    ls.add_var_bound(x0, GE, mpq(-3), e);
+    ls.add_var_bound(x0, LE, mpq(3), e);
+    ls.add_var_bound(x0, LE, mpq(1), e);
     ls.solve();
     my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
@@ -2853,9 +2858,10 @@ void test_bound_propagation_one_row_mixed() {
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(mpq(1), x0));
     c.push_back(std::pair<mpq, var_index>(mpq(-1), x1));
-    ls.add_constraint(c, EQ, one_of_type<mpq>());
+    explanation e;
+    ls.add_constraint(c, EQ, one_of_type<mpq>(), e);
     vector<implied_bound> ev;
-    ls.add_var_bound(x1, LE, mpq(1));
+    ls.add_var_bound(x1, LE, mpq(1), e);
     ls.solve();
     my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
@@ -2870,15 +2876,16 @@ void test_bound_propagation_two_rows() {
     c.push_back(std::pair<mpq, var_index>(mpq(1), x));
     c.push_back(std::pair<mpq, var_index>(mpq(2), y));
     c.push_back(std::pair<mpq, var_index>(mpq(3), z));
-    ls.add_constraint(c, GE, one_of_type<mpq>());
+    explanation e;
+    ls.add_constraint(c, GE, one_of_type<mpq>(), e);
     c.clear();
     c.push_back(std::pair<mpq, var_index>(mpq(3), x));
     c.push_back(std::pair<mpq, var_index>(mpq(2), y));
     c.push_back(std::pair<mpq, var_index>(mpq(y), z));
-    ls.add_constraint(c, GE, one_of_type<mpq>());
-    ls.add_var_bound(x, LE, mpq(2));
+    ls.add_constraint(c, GE, one_of_type<mpq>(), e);
+    ls.add_var_bound(x, LE, mpq(2), e);
     vector<implied_bound> ev;
-    ls.add_var_bound(y, LE, mpq(1));
+    ls.add_var_bound(y, LE, mpq(1), e);
     ls.solve();
     my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
@@ -2891,14 +2898,15 @@ void test_total_case_u() {
     unsigned y = ls.add_var(1, false);
     unsigned z = ls.add_var(2, false);
     vector<std::pair<mpq, var_index>> c;
+    explanation e;
     c.push_back(std::pair<mpq, var_index>(mpq(1), x));
     c.push_back(std::pair<mpq, var_index>(mpq(2), y));
     c.push_back(std::pair<mpq, var_index>(mpq(3), z));
-    ls.add_constraint(c, LE, one_of_type<mpq>());
-    ls.add_var_bound(x, GE, zero_of_type<mpq>());
-    ls.add_var_bound(y, GE, zero_of_type<mpq>());
+    ls.add_constraint(c, LE, one_of_type<mpq>(), e);
+    ls.add_var_bound(x, GE, zero_of_type<mpq>(), e);
+    ls.add_var_bound(y, GE, zero_of_type<mpq>(), e);
     vector<implied_bound> ev;
-    ls.add_var_bound(z, GE, zero_of_type<mpq>());
+    ls.add_var_bound(z, GE, zero_of_type<mpq>(), e);
     ls.solve();
     my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
@@ -2917,15 +2925,16 @@ void test_total_case_l(){
     unsigned y = ls.add_var(1, false);
     unsigned z = ls.add_var(2, false);
     vector<std::pair<mpq, var_index>> c;
+    explanation e;
     c.push_back(std::pair<mpq, var_index>(mpq(1), x));
     c.push_back(std::pair<mpq, var_index>(mpq(2), y));
     c.push_back(std::pair<mpq, var_index>(mpq(3), z));
-    ls.add_constraint(c, GE, one_of_type<mpq>());
-    ls.add_var_bound(x, LE, one_of_type<mpq>());
-    ls.add_var_bound(y, LE, one_of_type<mpq>());
+    ls.add_constraint(c, GE, one_of_type<mpq>(), e);
+    ls.add_var_bound(x, LE, one_of_type<mpq>(), e);
+    ls.add_var_bound(y, LE, one_of_type<mpq>(), e);
     ls.settings().presolve_with_double_solver_for_lar = true;
     vector<implied_bound> ev;
-    ls.add_var_bound(z, LE, zero_of_type<mpq>());
+    ls.add_var_bound(z, LE, zero_of_type<mpq>(), e);
     ls.solve();
     my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
@@ -3507,8 +3516,9 @@ void test_maximize_term() {
     term_ls.push_back(std::pair<mpq, var_index>(mpq(2), y));
     
     unsigned term_2x_pl_2y = solver.add_term(term_ls);
-    solver.add_var_bound(term_x_min_y,  LE, zero_of_type<mpq>());
-    solver.add_var_bound(term_2x_pl_2y, LE, mpq(5));
+    explanation e;
+    solver.add_var_bound(term_x_min_y,  LE, zero_of_type<mpq>(), e);
+    solver.add_var_bound(term_2x_pl_2y, LE, mpq(5), e);
     solver.find_feasible_solution();
     lp_assert(solver.get_status() == lp_status::OPTIMAL);
     solver.print_constraints(std::cout);
diff --git a/src/test/lp/smt_reader.h b/src/test/lp/smt_reader.h
index 4bce99765..4f456c3f6 100644
--- a/src/test/lp/smt_reader.h
+++ b/src/test/lp/smt_reader.h
@@ -386,17 +386,17 @@ namespace lp {
             return ret;
         }
         
-        void add_constraint_to_solver(lar_solver * solver, formula_constraint & fc) {
+        void add_constraint_to_solver(lar_solver * solver, formula_constraint & fc, explanation& e) {
             vector<std::pair<mpq, var_index>> ls;
             for (auto & it : fc.m_coeffs) {
                 ls.push_back(std::make_pair(it.first, solver->add_var(register_name(it.second), false)));
             }
-            solver->add_constraint(ls, fc.m_kind, fc.m_right_side);
+            solver->add_constraint(ls, fc.m_kind, fc.m_right_side, e);
         }
-
+        explanation e;
         void fill_lar_solver(lar_solver * solver) {
             for (formula_constraint & fc : m_constraints)
-                add_constraint_to_solver(solver, fc);
+                add_constraint_to_solver(solver, fc, e);
         }
 
 
diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp
index 401c17450..97e582504 100644
--- a/src/util/lp/lar_solver.cpp
+++ b/src/util/lp/lar_solver.cpp
@@ -28,7 +28,6 @@ void clear() {lp_assert(false); // not implemented
 
 
 lar_solver::lar_solver() : m_status(lp_status::UNKNOWN),
-                           m_infeasible_column_index(-1),
                            m_terms_start_index(1000000),
                            m_mpq_lar_core_solver(m_settings, *this),
                            m_int_solver(nullptr)
@@ -303,16 +302,8 @@ lp_status lar_solver::solve() {
     return m_status;
 }
 
-void lar_solver::fill_explanation_from_infeasible_column(explanation & evidence) const{
-
-    // this is the case when the lower bound is in conflict with the upper one
-    const ul_pair & ul =  m_columns_to_ul_pairs[m_infeasible_column_index];
-    evidence.push_justification(ul.upper_bound_witness(),  numeric_traits<mpq>::one());
-    evidence.push_justification(ul.lower_bound_witness(), -numeric_traits<mpq>::one());
-}
-
-    
 unsigned lar_solver::get_total_iterations() const { return m_mpq_lar_core_solver.m_r_solver.total_iterations(); }
+
 vector<unsigned> lar_solver::get_list_of_all_var_indices() const {
     vector<unsigned> ret;
     for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_heading.size(); j++)
@@ -323,7 +314,6 @@ void lar_solver::push() {
     m_simplex_strategy = m_settings.simplex_strategy();
     m_simplex_strategy.push();
     m_columns_to_ul_pairs.push();
-    m_infeasible_column_index.push();
     m_mpq_lar_core_solver.push();
     m_term_count = m_terms.size();
     m_term_count.push();
@@ -1207,10 +1197,6 @@ bool lar_solver::has_value(var_index var, mpq& value) const {
 
 void lar_solver::get_infeasibility_explanation(explanation& exp) const {
     exp.clear();
-    if (m_infeasible_column_index != -1) {
-        fill_explanation_from_infeasible_column(exp);
-        return;
-    }
     if (m_mpq_lar_core_solver.get_infeasible_sum_sign() == 0) {
         return;
     }
@@ -1781,17 +1767,17 @@ bool lar_solver::bound_is_integer_for_integer_column(unsigned j, const mpq & rig
     return right_side.is_int();
 }
 
-constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side) {
+constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side, explanation& e) {
     TRACE("lar_solver", tout << "j = " << j << std::endl;);
     constraint_index ci = m_constraints.size();
     if (!is_term(j)) { // j is a var
         lp_assert(bound_is_integer_for_integer_column(j, right_side));
         auto vc = new lar_var_constraint(j, kind, right_side);
         m_constraints.push_back(vc);
-        update_column_type_and_bound(j, kind, right_side, ci);
+        update_column_type_and_bound(j, kind, right_side, ci, e);
     }
     else {
-        add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
+        add_var_bound_on_constraint_for_term(j, kind, right_side, ci, e);
     }
     lp_assert(sizes_are_correct());
     return ci;
@@ -1816,57 +1802,44 @@ bool lar_solver::compare_values(impq const& lhs, lconstraint_kind k, const mpq &
     }
 }
 
-void lar_solver::update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
-    switch (m_mpq_lar_core_solver.m_column_types[j]) {
-    case column_type::free_column:
-        update_free_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::boxed:
-        update_boxed_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::lower_bound:
-        update_lower_bound_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::upper_bound:
-        update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::fixed:
-        update_fixed_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    default:
-        lp_assert(false); // cannot be here
-    }
+void lar_solver::update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side,
+                                              constraint_index constr_index,
+                                              explanation& e) {
+    if (column_has_upper_bound(j))
+        update_column_type_and_bound_with_ub(j, kind, right_side, constr_index, e);
+    else 
+        update_column_type_and_bound_with_no_ub(j, kind, right_side, constr_index, e);
 }
 
-void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci, explanation & e) {
     lp_assert(is_term(j));
     unsigned adjusted_term_index = adjust_term_index(j);
     //    lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int());
     unsigned term_j;
     if (m_var_register.external_is_used(j, term_j)) {
         m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side));
-        update_column_type_and_bound(term_j, kind, right_side, ci);
+        update_column_type_and_bound(term_j, kind, right_side, ci, e);
     }
     else {
-        add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side);
+        add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side, e);
     }
 }
 
-constraint_index lar_solver::add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm) {
+constraint_index lar_solver::add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm, explanation & e) {
     vector<std::pair<mpq, var_index>> left_side;
     substitute_terms_in_linear_expression(left_side_with_terms, left_side);
     unsigned term_index = add_term(left_side);
     constraint_index ci = m_constraints.size();
-    add_var_bound_on_constraint_for_term(term_index, kind_par, right_side_parm, ci);
+    add_var_bound_on_constraint_for_term(term_index, kind_par, right_side_parm, ci, e);
     return ci;
 }
 
 void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
-                                                                    lconstraint_kind kind, const mpq & right_side) {
+                                                                    lconstraint_kind kind, const mpq & right_side, explanation & e) {
 
     add_row_from_term_no_constraint(term, term_j);
     unsigned j = A_r().column_count() - 1;
-    update_column_type_and_bound(j, kind, right_side, m_constraints.size());
+    update_column_type_and_bound(j, kind, right_side, m_constraints.size(), e);
     m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
     lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
 }
@@ -1992,7 +1965,8 @@ void lar_solver::update_free_column_type_and_bound(var_index j, lconstraint_kind
     m_columns_with_changed_bound.insert(j);
 }
 
-void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci,
+                                                          explanation& e) {
     lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
     mpq y_of_bound(0);
     switch (kind) {
@@ -2019,7 +1993,8 @@ void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstrai
             m_columns_with_changed_bound.insert(j);
             if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                e.push_justification(ci);
+                e.push_justification(get_column_lower_bound_witness(j));
             }
             else {
                 m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_lower_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed;
@@ -2032,7 +2007,7 @@ void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstrai
             if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
                 set_lower_bound_witness(j, ci);
-                m_infeasible_column_index = j;
+                //                m_infeasible_column_index = j;
             }
             else {
                 m_mpq_lar_core_solver.m_r_lower_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
@@ -2068,7 +2043,7 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
 
             if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                //                m_infeasible_column_index = j;
             }
             else {
                 if (m_mpq_lar_core_solver.m_r_lower_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])
@@ -2088,7 +2063,7 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
             }
             if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
             }
             else if (low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
@@ -2100,12 +2075,12 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
             auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
             if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_upper_bound_witness(j, ci);
             }
             else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_lower_bound_witness(j, ci);
             }
             else {
@@ -2124,6 +2099,240 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
 
     }
 }
+void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const mpq & right_side, unsigned constraint_index, lp::explanation& e) {
+    SASSERT(column_has_upper_bound(j));
+    if (column_has_lower_bound(j)) {
+        update_bound_with_ub_lb(j, kind, right_side, constraint_index, e);
+    } else {
+        update_bound_with_ub_no_lb(j, kind, right_side, constraint_index, e);
+    }
+}
+
+void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const mpq & right_side, unsigned constraint_index, lp::explanation& e) {
+    SASSERT(!column_has_upper_bound(j));
+    if (column_has_lower_bound(j)) {
+        update_bound_with_no_ub_lb(j, kind, right_side, constraint_index, e);
+    } else {
+        update_bound_with_no_ub_no_lb(j, kind, right_side, constraint_index, e);
+    }
+}
+
+
+void lar_solver::update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci, explanation& e) {
+    lp_assert(column_has_lower_bound(j) && column_has_upper_bound(j));
+
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+	{
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_lower_bound_witness(j));
+            }
+            if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return;
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+            set_upper_bound_witness(j, ci);
+            m_columns_with_changed_bound.insert(j);
+	}
+	break;
+    case GT:
+        y_of_bound = 1;
+    case GE:
+	{
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_upper_bound_witness(j));
+            }
+            if (low < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                return;
+            }
+            m_mpq_lar_core_solver.m_r_lower_bounds[j] = low;
+            m_columns_with_changed_bound.insert(j);
+            set_lower_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_column_types[j] = (low == m_mpq_lar_core_solver.m_r_upper_bounds[j]? column_type::fixed : column_type::boxed);
+	}
+	break;
+    case EQ:
+	{
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_upper_bound_witness(j));
+            }
+            if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_lower_bound_witness(j));
+            }    
+            set_upper_bound_witness(j, ci);
+            set_lower_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
+            m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            break;
+	}
+
+    default:
+        lp_unreachable();
+    }
+}
+void lar_solver::update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci, explanation& e) {
+    lp_assert(column_has_lower_bound(j) && !column_has_upper_bound(j));
+
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+	{
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_lower_bound_witness(j));
+            }
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+            set_upper_bound_witness(j, ci);
+            m_columns_with_changed_bound.insert(j);
+            m_mpq_lar_core_solver.m_column_types[j] = (up == m_mpq_lar_core_solver.m_r_lower_bounds[j]? column_type::fixed : column_type::boxed);
+	}
+	break;
+    case GT:
+        y_of_bound = 1;
+    case GE:
+	{
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            if (low < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                return;
+            }
+            m_mpq_lar_core_solver.m_r_lower_bounds[j] = low;
+            m_columns_with_changed_bound.insert(j);
+            set_lower_bound_witness(j, ci);
+	}
+	break;
+    case EQ:
+	{
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_lower_bound_witness(j));
+            }
+                
+            set_upper_bound_witness(j, ci);
+            set_lower_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
+            m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            break;
+	}
+
+    default:
+        lp_unreachable();
+    }
+
+}
+
+void lar_solver::update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci, explanation& e) {
+    lp_assert(!column_has_lower_bound(j) && column_has_upper_bound(j));
+
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+	{
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return;
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+            set_upper_bound_witness(j, ci);
+            m_columns_with_changed_bound.insert(j);
+	}
+	break;
+    case GT:
+        y_of_bound = 1;
+    case GE:
+	{
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_upper_bound_witness(j));
+            }
+            m_mpq_lar_core_solver.m_r_lower_bounds[j] = low;
+            m_columns_with_changed_bound.insert(j);
+            set_lower_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_column_types[j] = (low == m_mpq_lar_core_solver.m_r_upper_bounds[j]? column_type::fixed : column_type::boxed);
+	}
+	break;
+    case EQ:
+	{
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                set_status(lp_status::INFEASIBLE);
+                e.push_justification(ci);
+                e.push_justification(get_column_upper_bound_witness(j));
+            }
+                
+            set_upper_bound_witness(j, ci);
+            set_lower_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
+            m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            break;
+	}
+
+    default:
+        lp_unreachable();
+    }
+}
+void lar_solver::update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci, explanation&) {
+    lp_assert(!column_has_lower_bound(j) && !column_has_upper_bound(j));
+    m_columns_with_changed_bound.insert(j);
+
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+	{
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+            set_upper_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound;
+	}
+	break;
+    case GT:
+        y_of_bound = 1;
+    case GE:
+	{
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_lower_bounds[j] = low;
+            m_columns_with_changed_bound.insert(j);
+            set_lower_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_column_types[j] = column_type::lower_bound;
+	}
+	break;
+    case EQ:
+	{
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            set_upper_bound_witness(j, ci);
+            set_lower_bound_witness(j, ci);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
+            m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            break;
+	}
+
+    default:
+        lp_unreachable();
+    }
+}
+
+
 void lar_solver::update_lower_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
     lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::lower_bound);
     mpq y_of_bound(0);
@@ -2139,7 +2348,7 @@ void lar_solver::update_lower_bound_column_type_and_bound(var_index j, lconstrai
 
             if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
             }
             else {
                 m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_lower_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed;
@@ -2163,7 +2372,7 @@ void lar_solver::update_lower_bound_column_type_and_bound(var_index j, lconstrai
             auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
             if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_upper_bound_witness(j, ci);
             }
             else {
@@ -2192,7 +2401,7 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
     case LT:
         if (v <= m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
             m_status = lp_status::INFEASIBLE;
-            m_infeasible_column_index = j;
+            // m_infeasible_column_index = j;
             set_upper_bound_witness(j, ci);
         }
         break;
@@ -2200,7 +2409,7 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
 	{
             if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_upper_bound_witness(j, ci);
             }
 	}
@@ -2209,7 +2418,7 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
 	{
             if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_lower_bound_witness(j, ci);
             }
 	}
@@ -2218,7 +2427,7 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
 	{
             if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_lower_bound_witness(j, ci);
             }
 	}
@@ -2227,12 +2436,12 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
 	{
             if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_upper_bound_witness(j, ci);
             }
             else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                m_infeasible_column_index = j;
+                // m_infeasible_column_index = j;
                 set_lower_bound_witness(j, ci);
             }
             break;
@@ -2266,18 +2475,19 @@ bool lar_solver::tighten_term_bounds_by_delta(unsigned term_index, const impq& d
             return false;
         }
     }
+    explanation e;
     TRACE("cube", tout << "can tighten";);
     if (slv.column_has_upper_bound(j)) {
         if (!is_zero(delta.y))
-            add_var_bound(tj, lconstraint_kind::LT, slv.m_upper_bounds[j].x - delta.x);
+            add_var_bound(tj, lconstraint_kind::LT, slv.m_upper_bounds[j].x - delta.x, e);
         else 
-            add_var_bound(tj, lconstraint_kind::LE, slv.m_upper_bounds[j].x - delta.x);
+            add_var_bound(tj, lconstraint_kind::LE, slv.m_upper_bounds[j].x - delta.x, e);
     }
     if (slv.column_has_lower_bound(j)) {
         if (!is_zero(delta.y))
-            add_var_bound(tj, lconstraint_kind::GT, slv.m_lower_bounds[j].x + delta.x);
+            add_var_bound(tj, lconstraint_kind::GT, slv.m_lower_bounds[j].x + delta.x, e);
         else 
-            add_var_bound(tj, lconstraint_kind::GE, slv.m_lower_bounds[j].x + delta.x);
+            add_var_bound(tj, lconstraint_kind::GE, slv.m_lower_bounds[j].x + delta.x, e);
     }
     return true;
 }
diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h
index 2ac1ddf5d..c50c28ff6 100644
--- a/src/util/lp/lar_solver.h
+++ b/src/util/lp/lar_solver.h
@@ -193,19 +193,25 @@ public:
 
     void add_basic_var_to_core_fields();
 
-    constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side) ;
+    constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side, explanation&) ;
 
     bool compare_values(var_index j, lconstraint_kind kind, const mpq & right_side);
 
     bool compare_values(impq const& lhs, lconstraint_kind k, const mpq & rhs);
 
-    void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
+    void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, explanation&);
+    void update_column_type_and_bound_with_ub(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, explanation&);
+    void update_column_type_and_bound_with_no_ub(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, explanation&);
+    void update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, explanation&);
+    void update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, explanation&);
+    void update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, explanation&);
+    void update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, explanation&);
 
-    void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
+    void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci, explanation&);
 
 
     void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
-                                                            lconstraint_kind kind, const mpq & right_side);
+                                                            lconstraint_kind kind, const mpq & right_side, explanation&);
 
     unsigned row_of_basic_column(unsigned) const;
     
@@ -222,7 +228,7 @@ public:
 
     void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind);
 
-    void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
+    void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci, explanation &);
     
     void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
     void update_lower_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
@@ -316,9 +322,6 @@ public:
     
     lp_status solve();
 
-    void fill_explanation_from_infeasible_column(explanation & evidence) const;
-
-    
     unsigned get_total_iterations() const;
     // see http://research.microsoft.com/projects/z3/smt07.pdf
     // This method searches for a feasible solution with as many different values of variables, reverenced in vars, as it can find
@@ -448,7 +451,7 @@ public:
 
     bool all_constrained_variables_are_registered(const vector<std::pair<mpq, var_index>>& left_side);
 
-    constraint_index add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm);
+    constraint_index add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm, explanation &);
     bool all_constraints_hold() const;
     bool constraint_holds(const lar_base_constraint & constr, std::unordered_map<var_index, mpq> & var_map) const;
     bool the_relations_are_of_same_type(const vector<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum) const;
diff --git a/src/util/lp/mps_reader.h b/src/util/lp/mps_reader.h
index bc2a8432e..72616509b 100644
--- a/src/util/lp/mps_reader.h
+++ b/src/util/lp/mps_reader.h
@@ -817,6 +817,7 @@ public:
     }
     
     void fill_lar_solver_on_row(row * row, lar_solver *solver)  {
+        explanation e;
         if (row->m_name != m_cost_row_name) {
             auto kind = get_lar_relation_from_row(row->m_type);
             vector<std::pair<mpq, var_index>> ls;
@@ -824,7 +825,7 @@ public:
                 var_index i = solver->add_var(get_var_index(s.first), false);
                 ls.push_back(std::make_pair(s.second, i));
             }
-            solver->add_constraint(ls, kind, row->m_right_side);
+            solver->add_constraint(ls, kind, row->m_right_side, e);
         } else {
             // ignore the cost row
         }
@@ -841,21 +842,24 @@ public:
         vector<std::pair<mpq, var_index>> ls;
         var_index i = solver->add_var(col->m_index, false);
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
-        solver->add_constraint(ls, GE, b->m_low);
+        explanation e;
+        solver->add_constraint(ls, GE, b->m_low, e);
     }
 
     void create_upper_constraint_for_var(column* col, bound * b, lar_solver *solver) {
         var_index i = solver->add_var(col->m_index, false);
         vector<std::pair<mpq, var_index>> ls;
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
-        solver->add_constraint(ls, LE, b->m_upper);
+        explanation e;
+        solver->add_constraint(ls, LE, b->m_upper, e);
     }
 
     void create_equality_contraint_for_var(column* col, bound * b, lar_solver *solver) {
         var_index i = solver->add_var(col->m_index, false);
         vector<std::pair<mpq, var_index>> ls;
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
-        solver->add_constraint(ls, EQ, b->m_fixed_value);
+        explanation e;
+        solver->add_constraint(ls, EQ, b->m_fixed_value, e);
     }
 
     void fill_lar_solver_on_columns(lar_solver * solver) {