diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h
index 974aa7cfb..bcc76c64f 100644
--- a/src/math/simplex/simplex.h
+++ b/src/math/simplex/simplex.h
@@ -111,7 +111,7 @@ namespace simplex {
         typedef typename matrix::col_iterator col_iterator;
 
         void  ensure_var(var_t v);
-        row   add_row(unsigned num_vars, var_t base, var_t const* vars, numeral const* coeffs);
+        row   add_row(var_t base, unsigned num_vars, var_t const* vars, numeral const* coeffs);
         void  del_row(row const& r);
         void  set_lower(var_t var, eps_numeral const& b);
         void  set_upper(var_t var, eps_numeral const& b);
@@ -125,6 +125,9 @@ namespace simplex {
         eps_numeral const& get_value(var_t v);
         void display(std::ostream& out) const;
 
+        unsigned get_num_vars() const { return m_vars.size(); }
+
+
     private:
 
         var_t select_var_to_fix();
@@ -159,10 +162,11 @@ namespace simplex {
         bool outside_bounds(var_t v) const { return below_lower(v) || above_upper(v); }
         bool is_free(var_t v) const { return !m_vars[v].m_lower_valid && !m_vars[v].m_upper_valid; }
         bool is_non_free(var_t v) const { return !is_free(v); }
-        unsigned get_num_vars() const { return m_vars.size(); }
         bool is_base(var_t x) const { return m_vars[x].m_is_base; }
+        void add_patch(var_t v);
 
         bool well_formed() const;
+        bool well_formed_row(row const& r) const;
         bool is_feasible() const;
     };
 
diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h
index 46915695f..631b89ce7 100644
--- a/src/math/simplex/simplex_def.h
+++ b/src/math/simplex/simplex_def.h
@@ -30,13 +30,26 @@ namespace simplex {
             SASSERT(found);
             );
         scoped_numeral base_coeff(m);
+        scoped_eps_numeral value(em), tmp(em);
         row r = M.mk_row();
         for (unsigned i = 0; i < num_vars; ++i) {
-            M.add(r, coeffs[i], vars[i]);
-            if (vars[i] == base) {
-                m.set(base_coeff, coeffs[i]);
+            if (!m.is_zero(coeffs[i])) {
+                var_t v = vars[i];
+                M.add_var(r, coeffs[i], v);
+                if (v == base) {
+                    m.set(base_coeff, coeffs[i]);
+                }
+                else {
+                    SASSERT(!is_base(v));
+                    em.mul(m_vars[v].m_value, coeffs[i], tmp);
+                    em.add(value, tmp, value);
+                }
             }
         }
+        em.neg(value);
+        em.div(value, base_coeff, value);
+        SASSERT(!m.is_zero(base_coeff));
+        SASSERT(!m_vars[base].m_is_base);
         while (m_row2base.size() <= r.id()) {
             m_row2base.push_back(null_var);
         }
@@ -44,9 +57,20 @@ namespace simplex {
         m_vars[base].m_base2row = r.id();
         m_vars[base].m_is_base = true;
         m.set(m_vars[base].m_base_coeff, base_coeff); 
+        em.set(m_vars[base].m_value, value);
+        add_patch(base);
+        SASSERT(well_formed_row(r));
         return r;
     }
 
+    template<typename Ext>
+    void simplex<Ext>::add_patch(var_t v) {
+        SASSERT(is_base(v));
+        if (outside_bounds(v)) {
+            m_to_patch.insert(v);
+        }
+    }
+
     template<typename Ext>
     void simplex<Ext>::del_row(row const& r) {
         m_vars[m_row2base[r.id()]].m_is_base = false;
@@ -114,8 +138,12 @@ namespace simplex {
             if (vi.m_lower_valid) out << em.to_string(vi.m_lower); else out << "-oo";
             out << ":";
             if (vi.m_upper_valid) out << em.to_string(vi.m_upper); else out << "oo";
-            out << "]";
-            if (vi.m_is_base) out << " b:" << vi.m_base2row;
+            out << "] ";
+            if (vi.m_is_base) out << "b:" << vi.m_base2row << " ";
+            col_iterator it = M.col_begin(i), end = M.col_end(i);
+            for (; it != end; ++it) {
+                out << "r" << it.get_row().id() << " ";
+            }
             out << "\n";
         }
     }
@@ -134,6 +162,7 @@ namespace simplex {
         unsigned num_iterations = 0;
         unsigned num_repeated = 0;
         var_t v = null_var;
+        SASSERT(well_formed());
         while ((v = select_var_to_fix()) != null_var) {
             if (m_cancel || num_iterations > m_max_iterations) {
                 return l_undef;
@@ -143,8 +172,8 @@ namespace simplex {
                 return l_false;
             }
             ++num_iterations;
-            SASSERT(well_formed());
         }
+        SASSERT(well_formed());
         return l_true;
     }
 
@@ -188,18 +217,18 @@ namespace simplex {
         m.set(x_jI.m_base_coeff, a_ij);
         x_jI.m_is_base = true;
         x_iI.m_is_base = false;
-        if (outside_bounds(x_j)) {
-            m_to_patch.insert(x_j);
-        }
+        add_patch(x_j);
+        SASSERT(well_formed_row(row(r_i)));
+
         col_iterator it = M.col_begin(x_j), end = M.col_end(x_j);
         scoped_numeral a_kj(m), g(m);
         for (; it != end; ++it) {
             row r_k = it.get_row();
-            if (r_k != r_i) {
+            if (r_k.id() != r_i) {
                 a_kj = it.get_row_entry().m_coeff;
                 a_kj.neg();
                 M.mul(r_k, a_ij);
-                M.add(r_k, a_kj, r_i);
+                M.add(r_k, a_kj, row(r_i));
                 var_t s = m_row2base[r_k.id()];
                 numeral& coeff = m_vars[s].m_base_coeff;
                 m.mul(coeff, a_ij, coeff);
@@ -207,6 +236,7 @@ namespace simplex {
                 if (!m.is_one(g)) {
                     m.div(coeff, g, coeff);
                 }
+                SASSERT(well_formed_row(row(r_k)));
             }
         }
     }
@@ -240,8 +270,8 @@ namespace simplex {
     void simplex<Ext>::update_value_core(var_t v, eps_numeral const& delta) {
         eps_numeral& val = m_vars[v].m_value;
         em.add(val, delta, val);
-        if (is_base(v) && outside_bounds(v)) {
-            m_to_patch.insert(v);
+        if (is_base(v)) {
+            add_patch(v);
         }
     }
 
@@ -705,18 +735,7 @@ namespace simplex {
             var_t s = m_row2base[i];
             if (s == null_var) continue;
             SASSERT(i == m_vars[s].m_base2row); 
-            //
-            // TBD: extract coefficient of base variable and compare
-            // with m_vars[s].m_base_coeff;
-            //
-            // check that sum of assignments add up to 0 for every row.
-            row_iterator it = M.row_begin(row(i)), end = M.row_end(row(i));
-            scoped_eps_numeral sum(em), tmp(em);
-            for (; it != end; ++it) {
-                em.mul(m_vars[it->m_var].m_value, it->m_coeff, tmp);
-                sum += tmp;
-            }
-            SASSERT(em.is_zero(sum));
+            SASSERT(well_formed_row(row(i)));            
         }
         return true;
     }
@@ -729,6 +748,25 @@ namespace simplex {
         return true;
     }
 
+    template<typename Ext>
+    bool simplex<Ext>::well_formed_row(row const& r) const {
+
+        //
+        // TBD: extract coefficient of base variable and compare
+        // with m_vars[s].m_base_coeff;
+        //
+        // check that sum of assignments add up to 0 for every row.
+        row_iterator it = M.row_begin(r), end = M.row_end(r);
+        scoped_eps_numeral sum(em), tmp(em);
+        for (; it != end; ++it) {
+            em.mul(m_vars[it->m_var].m_value, it->m_coeff, tmp);
+            sum += tmp;
+        }
+        SASSERT(em.is_zero(sum));
+
+        return true;
+    }
+
 };
 
 #endif
diff --git a/src/math/simplex/sparse_matrix.h b/src/math/simplex/sparse_matrix.h
index 6a85d69f0..d16d1d1cf 100644
--- a/src/math/simplex/sparse_matrix.h
+++ b/src/math/simplex/sparse_matrix.h
@@ -92,8 +92,7 @@ namespace simplex {
             void del_row_entry(unsigned idx);
             void compress(manager& m, vector<column> & cols); 
             void compress_if_needed(manager& m, vector<column> & cols);
-            void save_var_pos(svector<int> & result_map) const;
-            void reset_var_pos(svector<int> & result_map) const;
+            void save_var_pos(svector<int> & result_map, unsigned_vector& idxs) const;
             bool is_coeff_of(var_t v, numeral const & expected) const;
             int get_idx_of(var_t v) const;
         };
@@ -125,6 +124,7 @@ namespace simplex {
         svector<unsigned>       m_dead_rows;        // rows to recycle
         vector<column>          m_columns;          // per var
         svector<int>            m_var_pos;          // temporary map from variables to positions in row
+        unsigned_vector         m_var_pos_idx;      // indices in m_var_pos
 
         bool well_formed_row(unsigned row_id) const;
         bool well_formed_column(unsigned column_id) const;
@@ -138,7 +138,7 @@ namespace simplex {
         class row {
             unsigned m_id;            
         public:
-            row(unsigned r):m_id(r) {}
+            explicit row(unsigned r):m_id(r) {}
             row():m_id(UINT_MAX) {}
             bool operator!=(row const& other) const {
                 return m_id != other.m_id;
@@ -149,7 +149,7 @@ namespace simplex {
         void ensure_var(var_t v);
         
         row mk_row();
-        void add(row r, numeral const& n, var_t var);
+        void add_var(row r, numeral const& n, var_t var);
         void add(row r, numeral const& n, row src);
         void mul(row r, numeral const& n);
         void neg(row r);
diff --git a/src/math/simplex/sparse_matrix_def.h b/src/math/simplex/sparse_matrix_def.h
index 20e9d676f..a170ee47f 100644
--- a/src/math/simplex/sparse_matrix_def.h
+++ b/src/math/simplex/sparse_matrix_def.h
@@ -134,31 +134,18 @@ namespace simplex {
        \brief Fill the map var -> pos/idx
     */
     template<typename Ext>
-    inline void sparse_matrix<Ext>::_row::save_var_pos(svector<int> & result_map) const {
+    inline void sparse_matrix<Ext>::_row::save_var_pos(svector<int> & result_map, unsigned_vector& idxs) const {
         typename vector<_row_entry>::const_iterator it  = m_entries.begin();
         typename vector<_row_entry>::const_iterator end = m_entries.end();
         unsigned idx = 0;
         for (; it != end; ++it, ++idx) {
             if (!it->is_dead()) {
                 result_map[it->m_var] = idx;
+                idxs.push_back(it->m_var);
             }
         }
     }
 
-    /**
-       \brief Reset the map var -> pos/idx. That is for all variables v in the row, set result[v] = -1
-       This method can be viewed as the "inverse" of save_var_pos.
-    */
-    template<typename Ext>
-    inline void sparse_matrix<Ext>::_row::reset_var_pos(svector<int> & result_map) const {
-        typename vector<_row_entry>::const_iterator it  = m_entries.begin();
-        typename vector<_row_entry>::const_iterator end = m_entries.end();
-        for (; it != end; ++it) {
-            if (!it->is_dead()) {
-                result_map[it->m_var] = -1;
-            }
-        }
-    }
 
     template<typename Ext>
     int sparse_matrix<Ext>::_row::get_idx_of(var_t v) const {
@@ -316,7 +303,7 @@ namespace simplex {
     }
 
     template<typename Ext>
-    void sparse_matrix<Ext>::add(row dst, numeral const& n, var_t v) {
+    void sparse_matrix<Ext>::add_var(row dst, numeral const& n, var_t v) {
         _row& r   = m_rows[dst.id()];
         column& c = m_columns[v];
         unsigned r_idx;
@@ -339,7 +326,7 @@ namespace simplex {
         _row & r1 = m_rows[row1.id()];
         _row & r2 = m_rows[row2.id()];
         
-        r1.save_var_pos(m_var_pos);
+        r1.save_var_pos(m_var_pos, m_var_pos_idx);
 
         // 
         // loop over variables in row2,
@@ -376,7 +363,6 @@ namespace simplex {
                 }                                                       \
             }                                                           \
         }                                                               \
-        
         ((void) 0)
 
         if (m.is_one(n)) {
@@ -394,7 +380,11 @@ namespace simplex {
                     m.add(r_entry.m_coeff, tmp, r_entry.m_coeff));
         }
         
-        r1.reset_var_pos(m_var_pos);
+        // reset m_var_pos:
+        for (unsigned i = 0; i < m_var_pos_idx.size(); ++i) {
+            m_var_pos[m_var_pos_idx[i]] = -1;
+        }
+        m_var_pos_idx.reset();
         r1.compress_if_needed(m, m_columns);
     }
     
diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp
index def6f8be3..986caa30b 100644
--- a/src/test/simplex.cpp
+++ b/src/test/simplex.cpp
@@ -3,11 +3,123 @@
 #include "simplex.h"
 #include "simplex_def.h"
 #include "mpq_inf.h"
+#include "vector.h"
+#include "rational.h"
 
+#define R rational
 typedef simplex::simplex<simplex::mpz_ext> Simplex;
+typedef simplex::sparse_matrix<simplex::mpz_ext> sparse_matrix;
+
+static vector<R> vec(int i, int j) {
+    vector<R> nv;
+    nv.resize(2);
+    nv[0] = R(i);
+    nv[1] = R(j);
+    return nv;
+}
+
+static vector<R> vec(int i, int j, int k) {
+    vector<R> nv = vec(i, j);
+    nv.push_back(R(k));
+    return nv;
+}
+
+static vector<R> vec(int i, int j, int k, int l) {
+    vector<R> nv = vec(i, j, k);
+    nv.push_back(R(l));
+    return nv;
+}
+
+static vector<R> vec(int i, int j, int k, int l, int x) {
+    vector<R> nv = vec(i, j, k, l);
+    nv.push_back(R(x));
+    return nv;
+}
+
+static vector<R> vec(int i, int j, int k, int l, int x, int y) {
+    vector<R> nv = vec(i, j, k, l, x);
+    nv.push_back(R(y));
+    return nv;
+}
+
+static vector<R> vec(int i, int j, int k, int l, int x, int y, int z) {
+    vector<R> nv = vec(i, j, k, l, x, y);
+    nv.push_back(R(z));
+    return nv;
+}
+
+
+
+void add_row(Simplex& S, vector<R> const& _v, R const& _b, bool is_eq = false) {
+    unsynch_mpz_manager m;
+    unsigned_vector vars;
+    vector<R> v(_v);
+    R b(_b);
+    R l(denominator(b));
+    scoped_mpz_vector coeffs(m);
+    for (unsigned i = 0; i < v.size(); ++i) {
+        l = lcm(l, denominator(v[i]));
+        vars.push_back(i);
+        S.ensure_var(i);
+    }
+    b *= l;
+    b.neg();
+    for (unsigned i = 0; i < v.size(); ++i) {
+        v[i] *= l;
+        coeffs.push_back(v[i].to_mpq().numerator());
+    }
+    unsigned nv = S.get_num_vars();
+    vars.push_back(nv);
+    vars.push_back(nv+1);
+    S.ensure_var(nv);
+    S.ensure_var(nv+1);
+    coeffs.push_back(mpz(-1));
+    coeffs.push_back(b.to_mpq().numerator());
+    mpq_inf one(mpq(1),mpq(0));
+    mpq_inf zero(mpq(0),mpq(0));
+    SASSERT(vars.size() == coeffs.size());
+    std::cout << coeffs.size() << " " << nv << "\n";
+    for (unsigned i = 0; i < vars.size(); ++i) std::cout << vars[i] << " ";
+    std::cout << "\n";
+    S.set_lower(nv, zero);
+    if (is_eq) S.set_upper(nv, zero);
+    S.set_lower(nv+1, one);
+    S.set_upper(nv+1, one);
+    S.add_row(nv, coeffs.size(), vars.c_ptr(), coeffs.c_ptr());
+}
+
+static void feas(Simplex& S) {
+    S.display(std::cout);
+    lbool is_sat = S.make_feasible();
+    std::cout << "feasible: " << is_sat << "\n";
+    S.display(std::cout);
+}
+
+static void test1() {
+    Simplex S;
+    add_row(S, vec(1,0), R(1));
+    add_row(S, vec(0,1), R(1));
+    add_row(S, vec(1,1), R(1));
+    feas(S);
+}
+
+static void test2() {
+    Simplex S;
+    add_row(S, vec(1, 0), R(1));
+    add_row(S, vec(0, 1), R(1));
+    add_row(S, vec(1, 1), R(1), true);
+    feas(S);
+}
+
+static void test3() {
+    Simplex S;
+    add_row(S, vec(-1, 0), R(-1));
+    add_row(S, vec(0, -1), R(-1));
+    add_row(S, vec(1, 1), R(1), true);
+    feas(S);
+}
 
 void tst_simplex() {
-    simplex::sparse_matrix<simplex::mpz_ext> M;
     Simplex S;
 
     std::cout << "simplex\n";
@@ -37,4 +149,8 @@ void tst_simplex() {
     is_sat = S.make_feasible();
     std::cout << "feasible: " << is_sat << "\n";
     S.display(std::cout);
+
+    test1();
+    test2();
+    test3();
 }