diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt
index bde6ed93a..539c68712 100644
--- a/src/util/lp/CMakeLists.txt
+++ b/src/util/lp/CMakeLists.txt
@@ -6,6 +6,7 @@ z3_add_component(lp
     core_solver_pretty_printer.cpp
     dense_matrix.cpp
     eta_matrix.cpp
+    gomory.cpp
     indexed_vector.cpp
     int_solver.cpp
     lar_solver.cpp
diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp
new file mode 100644
index 000000000..55098dccf
--- /dev/null
+++ b/src/util/lp/gomory.cpp
@@ -0,0 +1,222 @@
+/*++
+  Copyright (c) 2017 Microsoft Corporation
+
+  Module Name:
+
+  <name>
+
+  Abstract:
+
+  <abstract>
+
+  Author:
+  Nikolaj Bjorner (nbjorner)
+  Lev Nachmanson (levnach)
+
+  Revision History:
+
+
+  --*/
+#include "util/lp/gomory.h"
+#include "util/lp/int_solver.h"
+#include "util/lp/lar_solver.h"
+namespace lp {
+
+class gomory::imp {
+    lar_term   &          m_t; // the term to return in the cut
+    mpq        &          m_k; // the right side of the cut
+    explanation&          m_ex; // the conflict explanation
+    unsigned              m_inf_col; // a basis column which has to be an integer but has a not integral value
+    const row_strip<mpq>& m_row;
+    const int_solver&           m_int_solver;
+
+
+    const impq & get_value(unsigned j) const { return m_int_solver.get_value(j); }
+    bool is_real(unsigned j) const { return m_int_solver.is_real(j); }
+    bool at_lower(unsigned j) const { return m_int_solver.at_lower(j); }
+    bool at_upper(unsigned j) const { return m_int_solver.at_upper(j); }
+    const impq & lower_bound(unsigned j) const { return m_int_solver.lower_bound(j); }
+    const impq & upper_bound(unsigned j) const  { return m_int_solver.upper_bound(j); }
+    constraint_index column_lower_bound_constraint(unsigned j) const { return m_int_solver.column_lower_bound_constraint(j); }
+    constraint_index column_upper_bound_constraint(unsigned j) const { return m_int_solver.column_upper_bound_constraint(j); }
+    bool column_is_fixed(unsigned j) const { return m_int_solver.m_lar_solver->column_is_fixed(j); }
+    void int_case_in_gomory_cut(const mpq & a, unsigned j,
+                                            mpq & lcm_den, const mpq& f0, const mpq& one_minus_f0) {
+        lp_assert(is_int(j) && !a.is_int());
+        mpq fj =  int_solver::fractional_part(a);
+        lp_assert(fj.is_pos());
+        TRACE("gomory_cut_detail", 
+              tout << a << " j=" << j << " k = " << m_k;
+              tout << ", fj: " << fj << ", ";
+              tout << "f0: " << f0 << ", ";
+              tout << "1 - f0: " << 1 - f0 << ", ";
+              tout << (at_lower(j)?"at_lower":"at_upper")<< std::endl;
+              );
+        mpq new_a;
+        mpq one_minus_fj = 1 - fj;
+        if (at_lower(j)) {
+            new_a = fj < one_minus_f0? fj / one_minus_f0 : one_minus_fj / f0;
+            m_k.addmul(new_a, lower_bound(j).x);
+            m_ex.push_justification(column_lower_bound_constraint(j), new_a);
+        }
+        else {
+            lp_assert(at_upper(j));
+            // the upper terms are inverted: therefore we have the minus
+            new_a = - (fj < f0? fj / f0 : one_minus_fj / one_minus_f0);
+            m_k.addmul(new_a, upper_bound(j).x);
+            m_ex.push_justification(column_upper_bound_constraint(j), new_a);
+        }
+        TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << m_k << "\n";);
+        m_t.add_monomial(new_a, j);
+        lcm_den = lcm(lcm_den, denominator(new_a));
+    }
+
+    void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f0, const mpq& one_minus_f0) {
+        TRACE("gomory_cut_detail_real", tout << "real\n";);
+        mpq new_a;
+        if (at_lower(x_j)) {
+            if (a.is_pos()) {
+                new_a  =  a / one_minus_f0;
+            }
+            else {
+                new_a  =  a / f0;
+                new_a.neg();
+            }
+            m_k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation than
+            // k += lower_bound(x_j).x * new_a;  
+            m_ex.push_justification(column_lower_bound_constraint(x_j), new_a);
+        }
+        else {
+            lp_assert(at_upper(x_j));
+            if (a.is_pos()) {
+                new_a =   a / f0; 
+                new_a.neg(); // the upper terms are inverted.
+            }
+            else {
+                new_a =   a / one_minus_f0; 
+            }
+            m_k.addmul(new_a, upper_bound(x_j).x); //  k += upper_bound(x_j).x * new_a; 
+            m_ex.push_justification(column_upper_bound_constraint(x_j), new_a);
+        }
+        TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << m_k << "\n";);
+        m_t.add_monomial(new_a, x_j);
+    }
+
+    lia_move report_conflict_from_gomory_cut() {
+        lp_assert(m_k.is_pos());
+        // conflict 0 >= k where k is positive
+        m_k.neg(); // returning 0 <= -k
+        return lia_move::conflict;
+    }
+
+    void adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) {
+        lp_assert(!m_t.is_empty());
+        auto pol = m_t.coeffs_as_vector();
+        m_t.clear();
+        if (pol.size() == 1) {
+            TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
+            unsigned v = pol[0].second;
+            lp_assert(is_int(v));
+            const mpq& a = pol[0].first;
+            m_k /= a;
+            if (a.is_pos()) { // we have av >= k
+                if (!m_k.is_int())
+                    m_k = ceil(m_k);
+                // switch size
+                m_t.add_monomial(- mpq(1), v);
+                m_k.neg();
+            } else {
+                if (!m_k.is_int())
+                    m_k = floor(m_k);
+                m_t.add_monomial(mpq(1), v);
+            }
+        } else {
+            TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;);
+            lcm_den = lcm(lcm_den, denominator(m_k));
+            lp_assert(lcm_den.is_pos());
+            if (!lcm_den.is_one()) {
+                // normalize coefficients of integer parameters to be integers.
+                for (auto & pi: pol) {
+                    pi.first *= lcm_den;
+                    SASSERT(!is_int(pi.second) || pi.first.is_int());
+                }
+                m_k *= lcm_den;
+            }
+            // negate everything to return -pol <= -m_k
+            for (const auto & pi: pol)
+                m_t.add_monomial(-pi.first, pi.second);
+            m_k.neg();
+        }
+        TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;);
+        lp_assert(m_k.is_int());
+    }
+public:
+    lia_move create_cut() {
+        TRACE("gomory_cut",
+              tout << "applying cut at:\n"; m_int_solver.m_lar_solver->print_row(m_row, tout); tout << std::endl;
+              for (auto & p : m_row) {
+                  m_int_solver.m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout);
+              }
+              tout << "inf_col = " << m_inf_col << std::endl;
+              );
+        
+        // gomory will be   t <= k and the current solution has a property t > k
+        m_k = 1;
+        mpq lcm_den(1);
+        bool some_int_columns = false;
+        mpq f0  = int_solver::fractional_part(get_value(m_inf_col));
+        mpq one_min_f0 = 1 - f0;
+        for (const auto & p : m_row) {
+            unsigned j = p.var();
+            if (column_is_fixed(j)) {
+                m_ex.push_justification(column_lower_bound_constraint(j));
+                m_ex.push_justification(column_upper_bound_constraint(j));
+                continue;
+            }
+            if (j == m_inf_col) {
+                lp_assert(p.coeff() == one_of_type<mpq>());
+                TRACE("gomory_cut_detail", tout << "seeing basic var";);
+                continue;
+            }
+            // make the format compatible with the format used in: Integrating Simplex with DPLL(T)
+            mpq a = - p.coeff();
+            if (is_real(j)) 
+                real_case_in_gomory_cut(a, j, f0, one_min_f0);
+            else if (!a.is_int()) { // fj will be zero and no monomial will be added
+                some_int_columns = true;
+                int_case_in_gomory_cut(a, j, lcm_den, f0, one_min_f0);
+            }
+        }
+
+        if (m_t.is_empty())
+            return report_conflict_from_gomory_cut();
+        if (some_int_columns)
+            adjust_term_and_k_for_some_ints_case_gomory(lcm_den);
+        lp_assert(m_int_solver.current_solution_is_inf_on_cut());
+        m_int_solver.m_lar_solver->subs_term_columns(m_t, m_k);
+        TRACE("gomory_cut", tout<<"gomory cut:"; m_int_solver.m_lar_solver->print_term(m_t, tout); tout << " <= " << m_k << std::endl;);
+        return lia_move::cut;
+    }
+    imp(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, const int_solver& int_slv ) :
+        m_t(t),
+        m_k(k),
+        m_ex(ex),
+        m_inf_col(basic_inf_int_j),
+        m_row(row),
+        m_int_solver(int_slv)
+    {
+    }
+    
+};
+
+lia_move gomory::create_cut() {
+    return m_imp->create_cut();
+}
+
+gomory::gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, const int_solver& s) {
+    m_imp = alloc(imp, t, k, ex, basic_inf_int_j, row, s);
+}
+
+gomory::~gomory() { dealloc(m_imp); }
+
+}
diff --git a/src/util/lp/gomory.h b/src/util/lp/gomory.h
new file mode 100644
index 000000000..b7946d6b0
--- /dev/null
+++ b/src/util/lp/gomory.h
@@ -0,0 +1,36 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+#include "util/lp/lar_term.h"
+#include "util/lp/lia_move.h"
+#include "util/lp/explanation.h"
+#include "util/lp/static_matrix.h"
+
+namespace lp {
+class int_solver;
+class gomory {
+    class imp;
+    imp                  *m_imp;
+public :
+    gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, const int_solver& s);
+    lia_move create_cut();
+    ~gomory();
+};
+}
diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp
index a77c202a0..cc01a0038 100644
--- a/src/util/lp/int_solver.cpp
+++ b/src/util/lp/int_solver.cpp
@@ -8,6 +8,7 @@
 #include "util/lp/lp_utils.h"
 #include <utility>
 #include "util/lp/monomial.h"
+#include "util/lp/gomory.h"
 namespace lp {
 
 
@@ -101,10 +102,7 @@ bool int_solver::is_gomory_cut_target(const row_strip<mpq>& row) {
     unsigned j;
     for (const auto & p : row) {
         j = p.var();
-        if (is_base(j)) continue;
-        if (!at_bound(j))
-            return false;
-        if (!is_zero(get_value(j).y)) {
+        if (!is_base(j) && (!at_bound(j) || !is_zero(get_value(j).y))) {
             TRACE("gomory_cut", tout << "row is not gomory cut target:\n";
                   display_column(tout, j);
                   tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";);
@@ -115,36 +113,6 @@ bool int_solver::is_gomory_cut_target(const row_strip<mpq>& row) {
 }
 
 
-void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0) {
-    TRACE("gomory_cut_detail_real", tout << "real\n";);
-    mpq new_a;
-    if (at_low(x_j)) {
-        if (a.is_pos()) {
-            new_a  =  a / one_minus_f_0;
-        }
-        else {
-            new_a  =  a / f_0;
-            new_a.neg();
-        }
-        m_k->addmul(new_a, lower_bound(x_j).x); // is it a faster operation than
-        // k += lower_bound(x_j).x * new_a;  
-        m_ex->push_justification(column_lower_bound_constraint(x_j), new_a);
-    }
-    else {
-        lp_assert(at_upper(x_j));
-        if (a.is_pos()) {
-            new_a =   a / f_0; 
-            new_a.neg(); // the upper terms are inverted.
-        }
-        else {
-            new_a =   a / one_minus_f_0; 
-        }
-        m_k->addmul(new_a, upper_bound(x_j).x); //  k += upper_bound(x_j).x * new_a; 
-        m_ex->push_justification(column_upper_bound_constraint(x_j), new_a);
-    }
-    TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << *m_k << "\n";);
-    m_t->add_monomial(new_a, x_j);
-}
 
 constraint_index int_solver::column_upper_bound_constraint(unsigned j) const {
     return m_lar_solver->get_column_upper_bound_witness(j);
@@ -155,99 +123,6 @@ constraint_index int_solver::column_lower_bound_constraint(unsigned j) const {
 }
 
 
-void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j,
-                                        mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) {
-    lp_assert(is_int(x_j));
-    lp_assert(!a.is_int());
-    mpq f_j =  fractional_part(a);
-    TRACE("gomory_cut_detail", 
-          tout << a << " x_j" << x_j << " k = " << *m_k << "\n";
-          tout << "f_j: " << f_j << "\n";
-          tout << "f_0: " << f_0 << "\n";
-          tout << "1 - f_0: " << 1 - f_0 << "\n";
-          tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl;
-          );
-    lp_assert (!f_j.is_zero());
-    mpq new_a;
-    if (at_low(x_j)) {
-        if (f_j <= one_minus_f_0) {
-            new_a = f_j / one_minus_f_0;
-        }
-        else {
-            new_a = (1 - f_j) / f_0;
-        }
-        m_k->addmul(new_a, lower_bound(x_j).x);
-        m_ex->push_justification(column_lower_bound_constraint(x_j), new_a);
-    }
-    else {
-        lp_assert(at_upper(x_j));
-        if (f_j <= f_0) {
-            new_a = f_j / f_0;
-        }
-        else {
-            new_a = (mpq(1) - f_j) / one_minus_f_0;
-        }
-        new_a.neg(); // the upper terms are inverted
-        m_k->addmul(new_a, upper_bound(x_j).x);
-        m_ex->push_justification(column_upper_bound_constraint(x_j), new_a);
-    }
-    TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << *m_k << "\n";);
-    m_t->add_monomial(new_a, x_j);
-    lcm_den = lcm(lcm_den, denominator(new_a));
-}
-
-lia_move int_solver::report_conflict_from_gomory_cut() {
-    TRACE("empty_pol",);
-    lp_assert(m_k->is_pos());
-    // conflict 0 >= k where k is positive
-    m_k->neg(); // returning 0 <= -k
-    return lia_move::conflict;
-}
-
-void int_solver::gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & pol,
-                                           lar_term & t,
-                                           mpq &k,
-                                           bool some_ints,
-                                           mpq & lcm_den) {
-    if (!some_ints)
-        return;
-        
-    t.clear();
-    if (pol.size() == 1) {
-        unsigned v = pol[0].second;
-        lp_assert(is_int(v));
-        bool k_is_int = k.is_int();
-        const mpq& a = pol[0].first;
-        k /= a;
-        if (a.is_pos()) { // we have av >= k
-            if (!k_is_int)
-                k = ceil(k);
-            // switch size
-            t.add_monomial(- mpq(1), v);
-            k.neg();
-        } else {
-            if (!k_is_int)
-                k = floor(k);
-            t.add_monomial(mpq(1), v);
-        }
-    } else if (some_ints) {
-        lcm_den = lcm(lcm_den, denominator(k));
-        lp_assert(lcm_den.is_pos());
-        if (!lcm_den.is_one()) {
-            // normalize coefficients of integer parameters to be integers.
-            for (auto & pi: pol) {
-                pi.first *= lcm_den;
-                SASSERT(!is_int(pi.second) || pi.first.is_int());
-            }
-            k *= lcm_den;
-        }
-        // negate everything to return -pol <= -k
-        for (const auto & pi: pol)
-            t.add_monomial(-pi.first, pi.second);
-        k.neg();
-    }
-}
-
 bool int_solver::current_solution_is_inf_on_cut() const {
     const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x;
     impq v = m_t->apply(x);
@@ -259,115 +134,18 @@ bool int_solver::current_solution_is_inf_on_cut() const {
     return v * sign > (*m_k) * sign;
 }
 
-void int_solver::adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) {
-    lp_assert(!m_t->is_empty());
-    auto pol = m_t->coeffs_as_vector();
-    m_t->clear();
-    if (pol.size() == 1) {
-        TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
-        unsigned v = pol[0].second;
-        lp_assert(is_int(v));
-        const mpq& a = pol[0].first;
-        (*m_k) /= a;
-        if (a.is_pos()) { // we have av >= k
-            if (!(*m_k).is_int())
-                (*m_k) = ceil((*m_k));
-            // switch size
-            m_t->add_monomial(- mpq(1), v);
-            (*m_k).neg();
-        } else {
-            if (!(*m_k).is_int())
-                (*m_k) = floor((*m_k));
-            m_t->add_monomial(mpq(1), v);
-        }
-    } else {
-        TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;);
-        lcm_den = lcm(lcm_den, denominator((*m_k)));
-        lp_assert(lcm_den.is_pos());
-        if (!lcm_den.is_one()) {
-            // normalize coefficients of integer parameters to be integers.
-            for (auto & pi: pol) {
-                pi.first *= lcm_den;
-                SASSERT(!is_int(pi.second) || pi.first.is_int());
-            }
-            (*m_k) *= lcm_den;
-        }
-        // negate everything to return -pol <= -(*m_k)
-        for (const auto & pi: pol)
-            m_t->add_monomial(-pi.first, pi.second);
-        (*m_k).neg();
-    }
-    TRACE("gomory_cut_detail", tout << "k = " << (*m_k) << std::endl;);
-    lp_assert((*m_k).is_int());
-}
-
-
-
-
 lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip<mpq> & row) {
-
     lp_assert(column_is_int_inf(inf_col));
 
-    TRACE("gomory_cut",
-          tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl;
-          for (auto & p : row) {
-              m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout);
-          }
-          tout << "inf_col = " << inf_col << std::endl;
-          );
-        
-    // gomory will be   t <= k and the current solution has a property t > k
-    *m_k = 1;
-    mpq lcm_den(1);
-    unsigned x_j;
-    mpq a;
-    bool some_int_columns = false;
-    mpq f_0  = int_solver::fractional_part(get_value(inf_col));
-    mpq one_min_f_0 = 1 - f_0;
-    for (const auto & p : row) {
-        x_j = p.var();
-        if (x_j == inf_col)
-            continue;
-        // make the format compatible with the format used in: Integrating Simplex with DPLL(T)
-        a = p.coeff();
-        a.neg();  
-        if (is_real(x_j)) 
-            real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0);
-        else if (!a.is_int()) { // f_j will be zero and no monomial will be added
-            some_int_columns = true;
-            int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0);
-        }
-    }
-
-    if (m_t->is_empty())
-        return report_conflict_from_gomory_cut();
-    if (some_int_columns)
-        adjust_term_and_k_for_some_ints_case_gomory(lcm_den);
-
-    lp_assert(current_solution_is_inf_on_cut());
-    m_lar_solver->subs_term_columns(*m_t);
-    TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;);
-    return lia_move::cut;
-}
-
-int int_solver::find_free_var_in_gomory_row(const row_strip<mpq>& row) {
-    unsigned j;
-    for (const auto & p : row) {
-        j = p.var();
-        if (!is_base(j) && is_free(j))
-            return static_cast<int>(j);
-    }
-    return -1;
+    gomory gc(*m_t, *m_k, *m_ex, inf_col, row, *this);
+    return gc.create_cut();
 }
 
 lia_move int_solver::proceed_with_gomory_cut(unsigned j) {
     const row_strip<mpq>& row = m_lar_solver->get_row(row_of_basic_column(j));
 
-    if (-1 != find_free_var_in_gomory_row(row))  
-        return lia_move::undef;
-
     if (!is_gomory_cut_target(row)) 
-        return lia_move::undef;
+        return create_branch_on_column(j);
 
     *m_upper = true;
     return mk_gomory_cut(j, row);
@@ -394,19 +172,25 @@ typedef monomial mono;
 
 
 // this will allow to enable and disable tracking of the pivot rows
-struct pivoted_rows_tracking_control {
-    lar_solver * m_lar_solver;
-    bool m_track_pivoted_rows;
-    pivoted_rows_tracking_control(lar_solver* ls) :
+struct check_return_helper {
+    lar_solver *     m_lar_solver;
+    const lia_move & m_r;
+    bool             m_track_pivoted_rows;
+    check_return_helper(lar_solver* ls, const lia_move& r) :
         m_lar_solver(ls),
+        m_r(r),
         m_track_pivoted_rows(ls->get_track_pivoted_rows())
     {
         TRACE("pivoted_rows", tout << "pivoted rows = " << ls->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;);
         m_lar_solver->set_track_pivoted_rows(false);
     }
-    ~pivoted_rows_tracking_control() {
+    ~check_return_helper() {
         TRACE("pivoted_rows", tout << "pivoted rows = " << m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;);
         m_lar_solver->set_track_pivoted_rows(m_track_pivoted_rows);
+        if (m_r == lia_move::cut || m_r == lia_move::branch) {
+            int_solver * s = m_lar_solver->get_int_solver();
+            m_lar_solver->adjust_cut_for_terms(*(s->m_t), *(s->m_k));
+        }
     }
 };
 
@@ -626,7 +410,7 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
     lia_move r = run_gcd_test();
     if (r != lia_move::undef) return r;
 
-    pivoted_rows_tracking_control pc(m_lar_solver);
+    check_return_helper pc(m_lar_solver, r);
 
     if(settings().m_int_pivot_fixed_vars_from_basis)
         m_lar_solver->pivot_fixed_vars_from_basis();
@@ -1126,7 +910,7 @@ bool int_solver::at_bound(unsigned j) const {
     }
 }
 
-bool int_solver::at_low(unsigned j) const {
+bool int_solver::at_lower(unsigned j) const {
     auto & mpq_solver = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
     switch (mpq_solver.m_column_types[j] ) {
     case column_type::fixed:
diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h
index ec708918d..dfe51711c 100644
--- a/src/util/lp/int_solver.h
+++ b/src/util/lp/int_solver.h
@@ -53,6 +53,13 @@ public:
     bool move_non_basic_column_to_bounds(unsigned j);
     lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex);    
     bool is_base(unsigned j) const;
+    bool is_real(unsigned j) const;
+    const impq & lower_bound(unsigned j) const;
+    const impq & upper_bound(unsigned j) const;
+    bool is_int(unsigned j) const;
+    const impq & get_value(unsigned j) const;
+    bool at_lower(unsigned j) const;
+    bool at_upper(unsigned j) const;
     
 private:
 
@@ -79,10 +86,7 @@ private:
     void add_to_explanation_from_fixed_or_boxed_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;
+private:
     bool is_boxed(unsigned j) const;
     bool is_fixed(unsigned j) const;
     bool is_free(unsigned j) const;
@@ -91,7 +95,6 @@ private:
     void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val);
     bool non_basic_columns_are_at_bounds() const;
     bool is_feasible() const;
-    const impq & get_value(unsigned j) const;
     bool column_is_int_inf(unsigned j) const;
     void trace_inf_rows() const;
     lia_move branch_or_sat();
@@ -104,14 +107,9 @@ private:
     bool move_non_basic_columns_to_bounds();
     void branch_infeasible_int_var(unsigned);
     lia_move mk_gomory_cut(unsigned inf_col, const row_strip<mpq>& row);
-    lia_move report_conflict_from_gomory_cut();
-    void adjust_term_and_k_for_some_ints_case_gomory(mpq& lcm_den);
     lia_move proceed_with_gomory_cut(unsigned j);
-    int find_free_var_in_gomory_row(const row_strip<mpq>& );
     bool is_gomory_cut_target(const row_strip<mpq>&);
     bool at_bound(unsigned j) const;
-    bool at_low(unsigned j) const;
-    bool at_upper(unsigned j) const;
     bool has_low(unsigned j) const;
     bool has_upper(unsigned j) const;
     unsigned row_of_basic_column(unsigned j) const;
@@ -126,17 +124,13 @@ public:
         lp_assert(is_rational(n));
         return n.x - floor(n.x);
     }
-private:
-    void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0);
-    void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0);
     constraint_index column_upper_bound_constraint(unsigned j) const;
     constraint_index column_lower_bound_constraint(unsigned j) const;
-    void display_row_info(std::ostream & out, unsigned row_index) const;
-    void gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k, bool num_ints, mpq &lcm_den);
     bool current_solution_is_inf_on_cut() const;
-public:
+
     bool shift_var(unsigned j, unsigned range);
 private:
+    void display_row_info(std::ostream & out, unsigned row_index) const;
     unsigned random();
     bool has_inf_int() const;
     lia_move create_branch_on_column(int j);
@@ -161,5 +155,5 @@ public:
     bool hnf_has_var_with_non_integral_value() const;
     bool hnf_cutter_is_full() const;
     void patch_nbasic_column(unsigned j, bool patch_only_int_vals);
-};
+  };
 }
diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp
index 1340d1826..09c3bd5b9 100644
--- a/src/util/lp/lar_solver.cpp
+++ b/src/util/lp/lar_solver.cpp
@@ -1656,6 +1656,7 @@ var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
         if (m_settings.bound_propagation())
             m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
     }
+    CTRACE("add_term_lar_solver", !m_v.is_zero(), print_term(*m_terms.back(), tout););
     lp_assert(m_var_register.size() == A_r().column_count());
     return ret;
 }
@@ -2265,6 +2266,16 @@ void lar_solver::set_cut_strategy(unsigned cut_frequency) {
     } 
 }
 
+void lar_solver::adjust_cut_for_terms(const lar_term& t, mpq & rs) {
+    for (const auto& p : t) {
+        if (!is_term(p.var())) continue;
+        const lar_term & p_term = get_term(p.var());
+        if (p_term.m_v.is_zero()) continue;
+        rs -= p.coeff() * p_term.m_v;
+    }
+}
+
+
 } // namespace lp
 
 
diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h
index 9afb70c72..3c0ed4fbf 100644
--- a/src/util/lp/lar_solver.h
+++ b/src/util/lp/lar_solver.h
@@ -535,7 +535,7 @@ public:
         return m_columns_to_ul_pairs()[j].lower_bound_witness();
     }
 
-    void subs_term_columns(lar_term& t) {
+    void subs_term_columns(lar_term& t, mpq & rs) {
         vector<std::pair<unsigned,unsigned>> columns_to_subs;
         for (const auto & m : t.m_coeffs) {
             unsigned tj = adjust_column_index_to_term_index(m.first);
@@ -545,9 +545,12 @@ public:
         for (const auto & p : columns_to_subs) {
             auto it = t.m_coeffs.find(p.first);
             lp_assert(it != t.m_coeffs.end());
+            const lar_term& lt = get_term(p.second);
             mpq v = it->second;
             t.m_coeffs.erase(it);
             t.m_coeffs[p.second] = v;
+            if (lt.m_v.is_zero()) continue;
+            rs -= v * lt.m_v;
         }
     }
 
@@ -584,5 +587,6 @@ public:
     lar_term get_term_to_maximize(unsigned ext_j) const;
     void set_cut_strategy(unsigned cut_frequency);
     bool sum_first_coords(const lar_term& t, mpq & val) const;
+    void adjust_cut_for_terms(const lar_term& t, mpq & rs);
 };
 }