diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 4953cdb0a..de6606964 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -2822,6 +2822,8 @@ public:
         st.update("cube-success", m_solver->settings().st().m_cube_success);
         st.update("arith-patches", m_solver->settings().st().m_patches);
         st.update("arith-patches-success", m_solver->settings().st().m_patches_success);
+        st.update("arith-hnf-calls", m_solver->settings().st().m_hnf_cutter_calls);
+        st.update("arith-hnf-cuts", m_solver->settings().st().m_hnf_cuts);
     }        
 };
     
diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp
index ab4611dab..b40bd545a 100644
--- a/src/test/lp/lp.cpp
+++ b/src/test/lp/lp.cpp
@@ -3410,6 +3410,8 @@ void test_gomory_cut_1() {
     g.mk_gomory_cut(t,  k, expl, inf_col, row);
 }
 
+void call_hnf(general_matrix & A);
+
 void test_hnf_m_less_than_n() {
 #ifdef Z3DEBUG
     general_matrix A;
@@ -3432,9 +3434,7 @@ void test_hnf_m_less_than_n() {
     v.push_back(mpq(1));
     v.push_back(mpq(5));
     A.push_row(v);
-    unsigned r;
-    mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); 
-    hnf<general_matrix> h(A, d, r);
+    call_hnf(A);
 #endif
 }
 void test_hnf_m_greater_than_n() {
@@ -3456,9 +3456,7 @@ void test_hnf_m_greater_than_n() {
     v.push_back(mpq(12));
     v.push_back(mpq(55));
     A.push_row(v);
-    unsigned r;
-    mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); 
-    hnf<general_matrix> h(A, d, r);
+    call_hnf(A);
 #endif
 }
 
@@ -3508,11 +3506,11 @@ void test_determinant() {
         M[0][0] = 3; M[0][1] = -1; M[0][2] = 1;
         M[1][0] = 1; M[1][1] =  0; M[1][2] = 0;
         M[2][0] = 0; M[2][1] =  1; M[2][2] = 4;
-        unsigned r;
+        svector<unsigned> r;
         std::cout << "M = "; M.print(std::cout, 4); endl(std::cout);
         mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r);
         std::cout << "det M = " << d  << std::endl;
-        std::cout << "rank = " << r << std::endl;
+        std::cout << "rank = " << r.size() << std::endl;
     }
     {
         auto M = general_matrix(4, 6);
@@ -3520,11 +3518,11 @@ void test_determinant() {
         M[1][0] = 1; M[1][1] =  0; M[1][2] = 0; M[1][3] = 0; M[1][4] = 2; M[1][5] =   7;
         M[2][0] = 0; M[2][1] =  1; M[2][2] = 4; M[2][3] = 0; M[2][4] = 2; M[2][5] =   8;
         M[3][0] = 6; M[3][1] = -2; M[3][2] = 2; M[3][3] = 2; M[3][4] = 6; M[3][5] =  -2;
-        unsigned r;
+        svector<unsigned> r;
         std::cout << "M = "; M.print(std::cout, 4); endl(std::cout);
         mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r);
         std::cout << "det M = " << d  << std::endl;
-        std::cout << "rank = " << r << std::endl;
+        std::cout << "rank = " << r.size() << std::endl;
     }
 #endif
 }
@@ -3539,11 +3537,11 @@ void fill_general_matrix(general_matrix & M) {
             M[i][j] = mpq(static_cast<int>(my_random() % 13) - 6);
 }
 
-void call_hnf(general_matrix A) {
-    unsigned r;
+void call_hnf(general_matrix& A) {
+    svector<unsigned> r;
     mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r);
-    if (r == A.row_count())
-        hnf<general_matrix> h(A, d, r);
+    A.shrink_to_rank(r);
+    hnf<general_matrix> h(A, d);
 }
 
 
@@ -3559,9 +3557,7 @@ void test_hnf_1_2() {
     v.push_back(mpq(5));
     v.push_back(mpq(26));
     A.push_row(v);
-    unsigned r;
-    mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); 
-    hnf<general_matrix> h(A, d, r);
+    call_hnf(A);
     std::cout << "test_hnf_1_2 passed" << std::endl;
 }
 void test_hnf_2_2() {
@@ -3575,9 +3571,7 @@ void test_hnf_2_2() {
     v.push_back(mpq(2));
     v.push_back(mpq(11));
     A.push_row(v);
-    unsigned r;
-    mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); 
-    hnf<general_matrix> h(A, d, r);
+    call_hnf(A);
 
     std::cout << "test_hnf_2_2 passed" << std::endl;
 }
@@ -3676,8 +3670,41 @@ void test_hnf_5_5() {
     std::cout << "test_hnf_5_5 passed" << std::endl;
 }
 
+void test_small_generated_hnf() {
+    std::cout << "test_small_rank_hnf" << std::endl;
+    general_matrix A;
+    vector<mpq> v;
+    v.push_back(mpq(5));
+    v.push_back(mpq(26));
+    A.push_row(v);
+    v.clear();
+    v.push_back(zero_of_type<mpq>());
+    v.push_back(zero_of_type<mpq>());
+    A.push_row(v);
+    call_hnf(A);
+    std::cout << "test_small_rank_hnf passed" << std::endl;
+}
+void test_larger_generated_hnf() {
+    std::cout << "test_larger_generated_rank_hnf" << std::endl;
+    general_matrix A;
+    vector<mpq> v;
+    v.push_back(zero_of_type<mpq>());
+    v.push_back(zero_of_type<mpq>());
+    v.push_back(zero_of_type<mpq>());
+    A.push_row(v);
+    A.push_row(v);
+    v.clear();
+    v.push_back(mpq(5));
+    v.push_back(mpq(26));
+    v.push_back(mpq(3));
+    A.push_row(v);
+    call_hnf(A);
+    std::cout << "test_larger_generated_rank_hnf passed" << std::endl;
+}
+
 void test_hnf() {
-    test_determinant();
+    test_larger_generated_hnf();
+    test_small_generated_hnf();
     test_hnf_1_2();
     test_hnf_3_3();
     test_hnf_4_4();
diff --git a/src/util/lp/dense_matrix_def.h b/src/util/lp/dense_matrix_def.h
index 2ea79ae41..f615c412d 100644
--- a/src/util/lp/dense_matrix_def.h
+++ b/src/util/lp/dense_matrix_def.h
@@ -23,7 +23,6 @@ Revision History:
 #include "util/lp/numeric_pair.h"
 #include "util/lp/dense_matrix.h"
 namespace lp {
-template <typename T> void print_vector(const vector<T> & t, std::ostream & out);
 template <typename T, typename X> dense_matrix<T, X>::dense_matrix(unsigned m, unsigned n) : m_m(m), m_n(n), m_values(m * n, numeric_traits<T>::zero()) {
 }
 
diff --git a/src/util/lp/explanation.h b/src/util/lp/explanation.h
new file mode 100644
index 000000000..c3d876afe
--- /dev/null
+++ b/src/util/lp/explanation.h
@@ -0,0 +1,31 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+namespace lp {
+struct explanation {
+    vector<std::pair<mpq, constraint_index>> m_explanation;
+    void push_justification(constraint_index j, const mpq& v) {
+        m_explanation.push_back(std::make_pair(v, j));
+    }
+    void push_justification(constraint_index j) {
+        m_explanation.push_back(std::make_pair(one_of_type<mpq>(), j));
+    }
+};
+}
diff --git a/src/util/lp/general_matrix.h b/src/util/lp/general_matrix.h
index 2994d6101..73521424c 100644
--- a/src/util/lp/general_matrix.h
+++ b/src/util/lp/general_matrix.h
@@ -50,7 +50,7 @@ public:
 
     
     unsigned row_count() const { return m_data.size(); }
-    unsigned column_count() const { return m_data[0].size(); }
+    unsigned column_count() const { return m_data.size() > 0? m_data[0].size() : 0; }
 
     class ref_row {
         general_matrix& m_matrix;
@@ -70,11 +70,15 @@ public:
     ref_row operator[](unsigned i) { return ref_row(*this, m_data[adjust_row(i)]); }
     ref_row_const operator[](unsigned i) const { return ref_row_const(*this, m_data[adjust_row(i)]); }
 
+#ifdef Z3DEBUG
     void print(std::ostream & out, unsigned blanks = 0) const {
         print_matrix<mpq>(m_data, out, blanks);
     }
-
-    void clear() { m_data.clear(); }
+    void print(std::ostream & out, const char * ss) const {
+        std::string s(ss);
+        out << s;
+        print(out, s.size());
+    }
 
     void print_submatrix(std::ostream & out, unsigned k, unsigned blanks = 0) const {
         general_matrix m(row_count() - k, column_count() - k);
@@ -85,9 +89,29 @@ public:
         print_matrix<mpq>(m.m_data, out, blanks);
     }
 
+#endif
+
+    void clear() { m_data.clear(); }
+
+    bool row_is_initialized_correctly(const vector<mpq>& row) {
+        lp_assert(row.size() == column_count());
+        for (unsigned j = 0; j < row.size(); j ++)
+            lp_assert(is_zero(row[j]));
+        return true;
+    }
+    
+    template <typename T>
+    void init_row_from_container(int i, const T & c, std::function<unsigned (unsigned)> column_fix) {
+        auto & row = m_data[adjust_row(i)];
+        lp_assert(row_is_initialized_correctly(row));
+        for (const auto & p : c) {
+            unsigned j = adjust_column(column_fix(p.var()));
+            row[j] = p.coeff();
+        }
+    }
     
     void copy_column_to_indexed_vector(unsigned entering, indexed_vector<mpq> &w ) const {
-        lp_assert(false); //
+        lp_assert(false); // not implemented
     }
     general_matrix operator*(const general_matrix & m) const {
         lp_assert(m.row_count() == column_count());
@@ -191,5 +215,38 @@ public:
             v.resize(n);
         }
     }
+
+    void shrink_to_rank(const svector<unsigned>& basis_rows) {
+        if (basis_rows.size() == row_count()) return;
+        vector<vector<mpq>> data; // todo : not efficient code
+        for (unsigned i : basis_rows)
+            data.push_back(m_data[i]);
+        m_data = data;
+    }
+
+    // used for debug only
+    general_matrix take_first_n_columns(unsigned n) const {
+        lp_assert(n <= column_count());
+        if (n == column_count())
+            return *this;
+        general_matrix ret(row_count(), n);
+        for (unsigned i = 0; i < row_count(); i++)
+            for (unsigned j = 0; j < n; j++)
+                ret[i][j] = (*this)[i][j];
+        return ret;
+    }
+    inline
+    friend vector<mpq> operator*(const vector<mpq> & f, const general_matrix& a) {
+        vector<mpq> r(a.column_count());
+        for (unsigned j = 0; j < a.column_count(); j ++) {
+            mpq t = zero_of_type<mpq>();
+            for (unsigned i = 0; i < a.row_count(); i++) {
+                t += f[i] * a[i][j];
+            }
+            r[j] = t;
+        }
+        return r;
+    }
 };
+
 }
diff --git a/src/util/lp/hnf.h b/src/util/lp/hnf.h
index 8227deb46..aa31d5d01 100644
--- a/src/util/lp/hnf.h
+++ b/src/util/lp/hnf.h
@@ -104,7 +104,7 @@ void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq
 
 
 
-template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
+template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r, svector<unsigned> & basis_rows) {
     lp_assert(m.row_count() <= m.column_count());
     for (unsigned i = r; i < m.row_count(); i++) {
         for (unsigned j = r; j < m.column_count(); j++) {
@@ -112,6 +112,7 @@ template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
                 if (i != r) {
                     m.transpose_rows(i, r);
                 }
+                basis_rows.push_back(i);
                 if (j != r) {
                     m.transpose_columns(j, r);
                 }
@@ -137,12 +138,12 @@ template <typename M> void pivot_column_non_fractional(M &m, unsigned & r) {
 }
 
 // returns the rank of the matrix
-template <typename M> unsigned to_lower_triangle_non_fractional(M &m) {
+template <typename M> void to_lower_triangle_non_fractional(M &m, svector<unsigned> & basis_rows ) {
     lp_assert(m.row_count() <= m.column_count());
     unsigned i = 0;
     for (; i < m.row_count() - 1; i++) {
-        if (!prepare_pivot_for_lower_triangle(m, i)) {
-            return i;
+        if (!prepare_pivot_for_lower_triangle(m, i, basis_rows)) {
+            return;
         }
         pivot_column_non_fractional(m, i);
     }
@@ -150,10 +151,10 @@ template <typename M> unsigned to_lower_triangle_non_fractional(M &m) {
     // go over the last row and try to find a non-zero in the row to the right of diagonal
     for (unsigned j = i; j < m.column_count(); j++) {
         if (!is_zero(m[i][j])) {
-            return m.row_count();
+            basis_rows.push_back(i);
+            break;
         }
     }
-    return i;
 }
 template <typename M>
 mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) {
@@ -173,8 +174,8 @@ mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) {
 }
 
 
-// it returns "r" - the rank of the matrix and the gcd of r-minors
-template <typename M> mpq determinant_of_rectangular_matrix(const M& m, unsigned & r) {
+// it fills "r" - the basic rows of m
+template <typename M> mpq determinant_of_rectangular_matrix(const M& m, svector<unsigned> & basis_rows) {
     if (m.column_count() < m.row_count())
         throw "not implemented"; // consider working with the transposed m, or create a "transposed" version if needed
     // The plan is to transform m to the lower triangular form by using non-fractional Gaussian Elimination by columns.
@@ -182,19 +183,18 @@ template <typename M> mpq determinant_of_rectangular_matrix(const M& m, unsigned
     // m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r.
     // The gcd of these minors is the return value
     auto mc = m;
-    r = to_lower_triangle_non_fractional(mc);
-    if (r == 0)
+    to_lower_triangle_non_fractional(mc, basis_rows);
+    if (basis_rows.size() == 0)
         return one_of_type<mpq>();
-    lp_assert(!is_zero(gcd_of_row_starting_from_diagonal(mc, r - 1)));
-    return gcd_of_row_starting_from_diagonal(mc, r - 1);
+    return gcd_of_row_starting_from_diagonal(mc, basis_rows.size() - 1);
 }
 
 template <typename M> mpq determinant(const M& m) {
     lp_assert(m.row_count() == m.column_count());
     auto mc = m;
-    unsigned r;
-    mpq d = determinant_of_rectangular_matrix(mc, r);
-    return r < m.row_count() ? zero_of_type<mpq>() : d;
+    svector<unsigned> basis_rows;
+    mpq d = determinant_of_rectangular_matrix(mc, basis_rows);
+    return basis_rows.size() < m.row_count() ? zero_of_type<mpq>() : d;
 }
 
 } // end of namespace hnf_calc
@@ -204,7 +204,7 @@ class hnf {
     // fields
 
 #ifdef Z3DEBUG
-    M &          m_H;
+    M            m_H;
     M            m_U;
     M            m_U_reverse;
 #endif
@@ -215,7 +215,7 @@ class hnf {
     unsigned     m_m;
     unsigned     m_n;
     mpq          m_d; // it is a positive number and a multiple of gcd of r-minors of m_A_orig, where r is the rank of m_A_orig
-    mpq          m_r; // the rank of m_A
+    //  we suppose that the rank of m_A is equal to row_count(), and that row_count() <= column_count(), that is m_A has the full rank
     unsigned     m_i;
     unsigned     m_j;
     mpq          m_R;
@@ -391,7 +391,7 @@ class hnf {
 
     void work_on_columns_less_than_i_in_the_triangle(unsigned i) {
         const mpq & mii = m_H[i][i];
-        lp_assert(is_pos(mii));
+        if (is_zero(mii)) return;
         for (unsigned j = 0; j < i; j++) {
             const mpq & mij = m_H[i][j];
             if (!is_pos(mij) && - mij < mii)
@@ -444,7 +444,7 @@ class hnf {
             const mpq & hij = m_H[i][j];
             if (is_pos(hij))
                 return false;
-            if (- hij >= hii)
+            if (!is_zero(hii) && - hij >= hii)
                 return false;
         }
         
@@ -594,7 +594,7 @@ private:
     }
     
 public:
-    hnf(M & A, const mpq & d, unsigned r) :
+    hnf(M & A, const mpq & d) :
 #ifdef Z3DEBUG
         m_H(A),
 #endif
@@ -604,12 +604,10 @@ public:
         m_m(A.row_count()),
         m_n(A.column_count()),
         m_d(d),
-        m_r(r),
         m_R(m_d),
         m_half_R(floor(m_R / 2))
     {
-        lp_assert(m_m > 0 && m_n > 0);
-        if (is_zero(m_d))
+        if (m_m == 0 || m_n == 0 || is_zero(m_d))
             return;
 #ifdef Z3DEBUG
         prepare_U_and_U_reverse();
@@ -627,7 +625,7 @@ public:
 #endif
     }
 
-
+    const M & W() const { return m_W; }
     
 };
 
diff --git a/src/util/lp/hnf_cutter.h b/src/util/lp/hnf_cutter.h
index c0071d8a0..8677356f0 100644
--- a/src/util/lp/hnf_cutter.h
+++ b/src/util/lp/hnf_cutter.h
@@ -22,25 +22,170 @@ Revision History:
 #include "util/lp/hnf.h"
 #include "util/lp/general_matrix.h"
 #include "util/lp/var_register.h"
+#include "util/lp/lia_move.h"
+#include "util/lp/explanation.h"
+
 namespace lp  {
 class hnf_cutter {
-    var_register m_var_register;
-    general_matrix m_A;
-    vector<const lar_term*> m_terms;
+    var_register               m_var_register;
+    general_matrix             m_A;
+    vector<const lar_term*>    m_terms;
+    vector<mpq>                m_right_sides;
+    unsigned                   m_row_count;
+    unsigned                   m_column_count;
+    std::function<unsigned ()> m_random_next;
 public:
+    hnf_cutter(std::function<unsigned()> random) : m_random_next(random) {}
+
     void clear() {
-        m_A.clear();
         m_var_register.clear();
         m_terms.clear();
+        m_row_count = m_column_count = 0;
     }
-    void add_term_to_A_for_hnf(const lar_term* t, const mpq &) {
+    void add_term(const lar_term* t, const mpq &rs) {
         m_terms.push_back(t);
         for (const auto &p : *t) {
-            m_var_register.register_user_var(p.var());
+            m_var_register.add_var(p.var());
         }
+        m_right_sides.push_back(rs);
+        if (m_terms.size() <= m_var_register.size()) { // capture the maximal acceptable sizes
+            m_row_count = m_terms.size();
+            m_column_count = m_var_register.size();
+        }
+        
     }
     void print(std::ostream & out) {
         out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl;
     }
+
+    void initialize_row(unsigned i) {
+        m_A.init_row_from_container(i, * m_terms[i], [this](unsigned j) { return m_var_register.add_var(j);});
+    }
+
+    void init_matrix_A() {
+        m_A = general_matrix(m_row_count, m_column_count);
+        // use the last suitable counts to make the number
+        // of rows less than or equal to the number of columns
+        for (unsigned i = 0; i < m_row_count; i++)
+            initialize_row(i);
+    }
+
+    // todo: as we need only one row i with non integral b[i] need to optimize later
+    void find_h_minus_1_b(const general_matrix& H, vector<mpq> & b) {
+        // the solution will be put into b
+        for (unsigned i = 0; i < H.row_count() ;i++) {
+            for (unsigned j = 0; j < i; j++) {
+                b[i] -= H[i][j]*b[j];
+            }
+            b[i] /= H[i][i];
+            // consider return from here if b[i] is not an integer and return i
+        }
+    }
+
+    vector<mpq> create_b(const svector<unsigned> & basis_rows) {
+        if (basis_rows.size() == m_right_sides.size())
+            return m_right_sides;
+        vector<mpq> b;
+        for (unsigned i : basis_rows)
+            b.push_back(m_right_sides[i]);
+        return b;
+    }
+
+    int find_cut_row_index(const vector<mpq> & b) {
+        int ret = -1;
+        int n = 0;
+        for (int i = 0; i < static_cast<int>(b.size()); i++) {
+            if (!is_int(b[i])) {
+                if (n == 0 ) {
+                    lp_assert(ret == -1);
+                    n = 1;
+                    ret = i;
+                } else {
+                    if (m_random_next() % (++n) == 0) {
+                       ret = i;
+                    }
+                }
+            }
+        }
+        return ret;
+    }
+
+
+    // fills e_i*H_minus_1
+    void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector<mpq> & row) {
+        // we solve x = ei * H_min_1
+        // or x * H = ei
+        unsigned m = H.row_count();
+        for (unsigned k = i + 1; k < m; k++) {
+            row[k] = zero_of_type<mpq>();
+        }
+        row[i] = one_of_type<mpq>() / H[i][i];
+        for(int k = i - 1; k >= 0; k--) {
+            mpq t = zero_of_type<mpq>();
+            for (unsigned l = k + 1; l <= i; l++) {
+                t += H[l][k]*row[l];
+            }
+            row[k] = -t / H[k][k]; 
+        }
+
+        // test region
+        vector<mpq> ei(H.row_count(), zero_of_type<mpq>());
+        ei[i] = one_of_type<mpq>();
+        vector<mpq> pr = row * H;
+        pr.shrink(ei.size());
+        lp_assert(ei == pr);
+        // end test region
+        
+    }
+
+    void fill_term(const vector<mpq> & row, lar_term& t) {
+        for (unsigned j = 0; j < row.size(); j++) {
+            if (!is_zero(row[j]))
+                t.add_monomial(row[j], m_var_register.local_var_to_user_var(j));
+        }
+    }
+    
+    lia_move create_cut(lar_term& t, mpq& k, explanation& ex, bool & upper) {
+        init_matrix_A();
+        svector<unsigned> basis_rows;
+        mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows);
+        if (basis_rows.size() < m_A.row_count())
+            m_A.shrink_to_rank(basis_rows);
+        
+        hnf<general_matrix> h(m_A, d);
+        //  general_matrix A_orig = m_A;
+        
+        vector<mpq> b = create_b(basis_rows);
+        vector<mpq> bcopy = b;
+        find_h_minus_1_b(h.W(), b);
+        
+        lp_assert(bcopy == h.W().take_first_n_columns(b.size()) * b);
+        int cut_row = find_cut_row_index(b);
+        if (cut_row == -1) {
+            return lia_move::undef;
+        }
+        // test region
+        /*
+        general_matrix U(m_A.column_count());
+        vector<mpq> rt(m_A.column_count());
+        for (unsigned i = 0; i < U.row_count(); i++) {
+            get_ei_H_minus_1(i, h.W(), rt);
+            vector<mpq> ft = rt * A_orig;
+            for (unsigned j = 0; j < ft.size(); j++)
+                U[i][j] = ft[j];
+        }
+        std::cout << "U reverse = "; U.print(std::cout, 12); std::cout << std::endl;
+        */
+        // end test region
+        
+        vector<mpq> row(m_A.column_count());
+        get_ei_H_minus_1(cut_row, h.W(), row);
+        vector<mpq> f = row * m_A;
+
+        fill_term(f, t);
+        k = floor(b[cut_row]);
+        upper = true;
+        return lia_move::cut;
+    }
 };
 }
diff --git a/src/util/lp/indexed_vector.cpp b/src/util/lp/indexed_vector.cpp
index 348196195..180291705 100644
--- a/src/util/lp/indexed_vector.cpp
+++ b/src/util/lp/indexed_vector.cpp
@@ -42,12 +42,12 @@ template void lp::indexed_vector<double>::print(std::basic_ostream<char,struct s
 template void lp::indexed_vector<lp::numeric_pair<lp::mpq> >::print(std::ostream&);
 #endif
 }
-template void lp::print_vector<double>(vector<double> const&, std::ostream&);
-template void lp::print_vector<unsigned int>(vector<unsigned int> const&, std::ostream&);
-template void lp::print_vector<std::string>(vector<std::string> const&, std::ostream&);
-template void lp::print_vector<lp::numeric_pair<lp::mpq> >(vector<lp::numeric_pair<lp::mpq>> const&, std::ostream&);
+// template void lp::print_vector<double, vectro>(vector<double> const&, std::ostream&);
+// template void lp::print_vector<unsigned int>(vector<unsigned int> const&, std::ostream&);
+// template void lp::print_vector<std::string>(vector<std::string> const&, std::ostream&);
+// template void lp::print_vector<lp::numeric_pair<lp::mpq> >(vector<lp::numeric_pair<lp::mpq>> const&, std::ostream&);
 template void lp::indexed_vector<double>::resize(unsigned int);
-template void lp::print_vector< lp::mpq>(vector< lp::mpq> const &, std::basic_ostream<char, std::char_traits<char> > &);
-template void lp::print_vector<std::pair<lp::mpq, unsigned int> >(vector<std::pair<lp::mpq, unsigned int>> const&, std::ostream&);
+// template void lp::print_vector< lp::mpq>(vector< lp::mpq> const &, std::basic_ostream<char, std::char_traits<char> > &);
+// template void lp::print_vector<std::pair<lp::mpq, unsigned int> >(vector<std::pair<lp::mpq, unsigned int>> const&, std::ostream&);
 template void lp::indexed_vector<lp::numeric_pair<lp::mpq> >::erase_from_index(unsigned int);
 
diff --git a/src/util/lp/indexed_vector.h b/src/util/lp/indexed_vector.h
index cd313f446..d6ff4e76a 100644
--- a/src/util/lp/indexed_vector.h
+++ b/src/util/lp/indexed_vector.h
@@ -28,8 +28,6 @@ Revision History:
 #include <unordered_set>
 namespace lp {
 
-template <typename T> void print_vector(const vector<T> & t, std::ostream & out);
-template <typename T> void print_vector(const buffer<T> & t, std::ostream & out);
 template <typename T> void print_sparse_vector(const vector<T> & t, std::ostream & out);
 
 void print_vector_as_doubles(const vector<mpq> & t, std::ostream & out);
diff --git a/src/util/lp/indexed_vector_def.h b/src/util/lp/indexed_vector_def.h
index e4a517b55..2f7706089 100644
--- a/src/util/lp/indexed_vector_def.h
+++ b/src/util/lp/indexed_vector_def.h
@@ -22,21 +22,6 @@ Revision History:
 #include "util/lp/lp_settings.h"
 namespace lp {
 
-template <typename T>
-void print_vector(const vector<T> & t, std::ostream & out) {
-    for (unsigned i = 0; i < t.size(); i++)
-        out << t[i] << " ";
-    out << std::endl;
-}
-
-
-template <typename T>
-void print_vector(const buffer<T> & t, std::ostream & out) {
-    for (unsigned i = 0; i < t.size(); i++)
-        out << t[i] << " ";
-    out << std::endl;
-}
-
 template <typename T>
 void print_sparse_vector(const vector<T> & t, std::ostream & out) {
     for (unsigned i = 0; i < t.size(); i++) {
diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp
index f5d06c816..94121d27f 100644
--- a/src/util/lp/int_solver.cpp
+++ b/src/util/lp/int_solver.cpp
@@ -272,14 +272,20 @@ void int_solver::gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & po
 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);
-    TRACE(
-        "current_solution_is_inf_on_cut", tout << "v = " << v << " k = " << (*m_k) << std::endl;
-        if (v <=(*m_k)) {
-            tout << "v <= k - it should not happen!\n";
-        }
+    mpq sign = !(*m_upper) ? one_of_type<mpq>()  : -one_of_type<mpq>();
+    TRACE("current_solution_is_inf_on_cut",
+        if (is_pos(sign)) {
+                 tout << "v = " << v << " k = " << (*m_k) << std::endl;
+                if (v <=(*m_k)) {
+                    tout << "v <= k - it should not happen!\n";
+                }
+            } else {
+                if (v >= (*m_k)) {
+                    tout << "v > k - it should not happen!\n";
+                }
+            }
           );
-    
-    return v > (*m_k);
+    return v * sign > (*m_k) * sign;
 }
 
 void int_solver::adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) {
@@ -635,33 +641,44 @@ lia_move int_solver::gomory_cut() {
 }
 
 
-void int_solver::try_add_term_to_A_for_hnf(unsigned i) {
+bool int_solver::try_add_term_to_A_for_hnf(unsigned i) {
     mpq rs;
     const lar_term* t = m_lar_solver->terms()[i];
     for (const auto & p : *t) {
         if (!is_int(p.var()))
-            return; // todo : the mix case!
+            return false; // todo : the mix case!
+    }
+    if (m_lar_solver->get_equality_for_term_on_corrent_x(i, rs)) {
+        m_hnf_cutter.add_term(t, rs);
+        return true;
+    } else {
+        return false;
     }
-    if (!m_lar_solver->get_equality_for_term_on_corrent_x(i, rs))
-        return;
-    m_hnf_cutter.add_term_to_A_for_hnf(t, rs);
 }
 
 bool int_solver::hnf_matrix_is_empty() const { return true; }
 
 bool int_solver::prepare_matrix_A_for_hnf_cut() {
     m_hnf_cutter.clear();
-    for (unsigned i = 0; i < m_lar_solver->terms().size(); i++)
-        try_add_term_to_A_for_hnf(i);
-    m_hnf_cutter.print(std::cout);
-    return ! hnf_matrix_is_empty();
+    for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) {
+        bool r = try_add_term_to_A_for_hnf(i);
+        if (!r && settings().hnf_cutter_exit_if_x_is_not_on_bound_or_mixed )
+            return false;
+    }
+    return true;
 }
 
 
 lia_move int_solver::make_hnf_cut() {
-    if( !prepare_matrix_A_for_hnf_cut())
+    if (!prepare_matrix_A_for_hnf_cut()) {
         return lia_move::undef;
-    return lia_move::undef;
+    }
+    settings().st().m_hnf_cutter_calls++;
+    lia_move r =  m_hnf_cutter.create_cut(*m_t, *m_k, *m_ex, *m_upper);
+    CTRACE("hnf_cut", r == lia_move::cut, tout<< "cut:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;);
+    if (r == lia_move::cut)
+        settings().st().m_hnf_cuts++;
+    return r;
 }
 
 lia_move int_solver::hnf_cut() {
@@ -1009,10 +1026,12 @@ int_solver::int_solver(lar_solver* lar_slv) :
     m_lar_solver(lar_slv),
     m_branch_cut_counter(0),
     m_chase_cut_solver([this](unsigned j) {return m_lar_solver->get_column_name(j);},
-                 [this](unsigned j, std::ostream &o) {m_lar_solver->print_constraint(j, o);},
-                 [this]() {return m_lar_solver->A_r().column_count();},
-                 [this](unsigned j) {return get_value(j);},
-                 settings()) {
+                       [this](unsigned j, std::ostream &o) {m_lar_solver->print_constraint(j, o);},
+                       [this]() {return m_lar_solver->A_r().column_count();},
+                       [this](unsigned j) {return get_value(j);},
+                       settings()),
+    m_hnf_cutter([this](){ return settings().random_next(); })
+{
     m_lar_solver->set_int_solver(this);
 }
 
diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h
index 7975025d8..210f350e1 100644
--- a/src/util/lp/int_solver.h
+++ b/src/util/lp/int_solver.h
@@ -25,30 +25,15 @@ Revision History:
 #include "util/lp/chase_cut_solver.h"
 #include "util/lp/lar_constraints.h"
 #include "util/lp/hnf_cutter.h"
+#include "util/lp/lia_move.h"
+#include "util/lp/explanation.h"
 
 namespace lp {
 class lar_solver;
+
 template <typename T, typename X>
 struct lp_constraint;
-enum class lia_move {
-    sat,
-    branch,
-    cut,
-    conflict,
-    continue_with_check,
-    undef,
-    unsat
-};
 
-struct explanation {
-    vector<std::pair<mpq, constraint_index>> m_explanation;
-    void push_justification(constraint_index j, const mpq& v) {
-        m_explanation.push_back(std::make_pair(v, j));
-    }
-    void push_justification(constraint_index j) {
-        m_explanation.push_back(std::make_pair(one_of_type<mpq>(), j));
-    }
-};
 
 class int_solver {
 public:
@@ -182,6 +167,6 @@ public:
     lia_move make_hnf_cut();
     bool prepare_matrix_A_for_hnf_cut();
     bool hnf_matrix_is_empty() const;
-    void try_add_term_to_A_for_hnf(unsigned term_index);
+    bool try_add_term_to_A_for_hnf(unsigned term_index);
 };
 }
diff --git a/src/util/lp/lia_move.h b/src/util/lp/lia_move.h
new file mode 100644
index 000000000..ec4643e20
--- /dev/null
+++ b/src/util/lp/lia_move.h
@@ -0,0 +1,31 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+namespace lp {
+enum class lia_move {
+    sat,
+    branch,
+    cut,
+    conflict,
+    continue_with_check,
+    undef,
+    unsat
+};
+}
diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h
index 33df5a180..eb460d1ad 100644
--- a/src/util/lp/lp_settings.h
+++ b/src/util/lp/lp_settings.h
@@ -112,6 +112,8 @@ struct stats {
     unsigned m_cube_success;
     unsigned m_patches;
     unsigned m_patches_success;
+    unsigned m_hnf_cutter_calls;
+    unsigned m_hnf_cuts;
     stats() { reset(); }
     void reset() { memset(this, 0, sizeof(*this)); }
 };
@@ -141,22 +143,22 @@ private:
     random_gen m_rand;
 
 public:
-    unsigned reps_in_scaler;
+    unsigned       reps_in_scaler;
     // when the absolute value of an element is less than pivot_epsilon
     // in pivoting, we treat it as a zero
-    double pivot_epsilon;
+    double         pivot_epsilon;
     // see Chatal, page 115
-    double positive_price_epsilon;
+    double         positive_price_epsilon;
     // a quatation "if some choice of the entering vairable 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 ...
-    double entering_diag_epsilon;
-    int c_partial_pivoting; // this is the constant c from page 410
-    unsigned depth_of_rook_search;
-    bool using_partial_pivoting;
+    double         entering_diag_epsilon;
+    int            c_partial_pivoting; // this is the constant c from page 410
+    unsigned       depth_of_rook_search;
+    bool           using_partial_pivoting;
     // dissertation of Achim Koberstein
     // if Bx - b is different at any component more that refactor_epsilon then we refactor
-    double refactor_tolerance;
-    double pivot_tolerance;
+    double         refactor_tolerance;
+    double         pivot_tolerance;
     double zero_tolerance;
     double drop_tolerance;
     double tolerance_for_artificials;
@@ -175,6 +177,7 @@ public:
     double dual_feasibility_tolerance; // // page 71 of the PhD thesis of Achim Koberstein
     double primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
     double relative_primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
+    bool         hnf_cutter_exit_if_x_is_not_on_bound_or_mixed = true;
 
     bool m_bound_propagation;
     
diff --git a/src/util/lp/lp_utils.h b/src/util/lp/lp_utils.h
index fca4f782f..aa804919f 100644
--- a/src/util/lp/lp_utils.h
+++ b/src/util/lp/lp_utils.h
@@ -22,6 +22,13 @@ Revision History:
 #include "util/lp/numeric_pair.h"
 #include "util/debug.h"
 #include <unordered_map>
+template <typename C>
+void print_vector(const C & t, std::ostream & out) {
+    for (const auto & p : t)
+        out << p << " ";
+    out << std::endl;
+}
+
 template <typename A, typename B>
 bool try_get_value(const std::unordered_map<A,B> & map, const A& key, B & val) {
     const auto it = map.find(key);
diff --git a/src/util/lp/var_register.h b/src/util/lp/var_register.h
new file mode 100644
index 000000000..e24667e94
--- /dev/null
+++ b/src/util/lp/var_register.h
@@ -0,0 +1,48 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+namespace lp  {
+class var_register {
+    svector<unsigned> m_local_vars_to_external;
+    std::unordered_map<unsigned, unsigned> m_external_vars_to_local;
+public:
+    unsigned add_var(unsigned user_var) {
+        auto t = m_external_vars_to_local.find(user_var);
+        if (t != m_external_vars_to_local.end()) {
+            return t->second;
+        }
+
+        unsigned ret = size();
+        m_external_vars_to_local[user_var] = ret;
+        m_local_vars_to_external.push_back(user_var);
+        return ret;
+    }
+
+    unsigned local_var_to_user_var(unsigned local_var) const {
+        return m_local_vars_to_external[local_var];
+    }
+    unsigned size() const {
+        return m_local_vars_to_external.size();
+    }
+    void clear() {
+        m_local_vars_to_external.clear();
+        m_external_vars_to_local.clear();
+    }
+};
+}
diff --git a/src/util/vector.h b/src/util/vector.h
index c4443255a..edcee7449 100644
--- a/src/util/vector.h
+++ b/src/util/vector.h
@@ -188,6 +188,24 @@ public:
         m_data = nullptr;
     }
 
+    bool operator==(vector const & other) const {
+        if (this == &other) {
+            return true;
+        }
+        if (size() != other.size())
+            return false;
+        for (unsigned i = 0; i < size(); i++) {
+            if ((*this)[i] != other[i])
+                return false;
+        }
+        return true;
+    }
+
+    bool operator!=(vector const & other) const {
+        return !(*this == other);
+    }
+
+    
     vector & operator=(vector const & source) {
         if (this == &source) {
             return *this;