diff --git a/src/math/polysat/fixplex_def.h b/src/math/polysat/fixplex_def.h
index 252b042d3..6a5db1ef7 100644
--- a/src/math/polysat/fixplex_def.h
+++ b/src/math/polysat/fixplex_def.h
@@ -608,7 +608,7 @@ namespace polysat {
             return;
         if (m_var_is_touched.get(v, false))
             return;
-        m_var_is_touched.set(v, true);
+        m_var_is_touched.setx(v, true, false);
         for (auto idx : m_var2ineqs[v]) {
             if (!m_ineqs[idx].is_active) {
                 m_ineqs[idx].is_active = true;
diff --git a/src/test/fixplex.cpp b/src/test/fixplex.cpp
index 3ef9774b4..5240c605c 100644
--- a/src/test/fixplex.cpp
+++ b/src/test/fixplex.cpp
@@ -1,3 +1,4 @@
+
 #include "math/polysat/fixplex.h"
 #include "math/polysat/fixplex_def.h"
 #include "ast/bv_decl_plugin.h"
@@ -367,10 +368,168 @@ namespace polysat {
 
         std::cout << "number of failures: " << num_bad << "\n";
     }
+
+    static void test_lp(
+        vector<svector<std::pair<unsigned, uint64_t>>> const& rows,
+        svector<ineq> const& ineqs,
+        vector<mod_interval<uint64_t>> const& bounds) {
+
+        unsigned num_vars = 0;
+        for (auto const& row : rows)
+            for (auto [v, c] : row)
+                num_vars = std::max(num_vars, v);
+        for (auto const& i : ineqs)
+            num_vars = std::max(num_vars, std::max(i.v, i.w));
+        num_vars = std::max(bounds.size(), num_vars + 1);
+
+        ast_manager m;
+        reg_decl_plugins(m);
+        bv_util bv(m);
+        smt_params pa;
+        smt::kernel solver(m, pa);
+        expr_ref_vector variables(m);
+        for (unsigned i = 0; i < num_vars; ++i)
+            variables.push_back(m.mk_fresh_const("v", bv.mk_sort(32)));
+        auto mk_ult = [&](expr* a, expr* b) {
+            return m.mk_not(bv.mk_ule(b, a));
+        };
+
+        scoped_fp fp;
+        for (auto& row : rows) {
+            vector<rational> coeffs;
+            svector<var_t> vars;
+            expr_ref term(m);
+            term = bv.mk_numeral(0, 32);
+            for (auto [v, c] : row) {
+                coeffs.push_back(rational(c, rational::ui64()));
+                vars.push_back(v);
+                term = bv.mk_bv_add(term, bv.mk_bv_mul(bv.mk_numeral(c, 32), variables.get(v)));
+            }
+            fp.add_row(row[0].first, row.size(), vars.data(), coeffs.data());
+            solver.assert_expr(m.mk_eq(bv.mk_numeral(0, 32), term));
+        }
+        unsigned index = 0;
+        for (auto const& i : ineqs) {
+            ++index;
+            if (i.strict) {
+                fp.add_lt(i.v, i.w, index);
+                solver.assert_expr(mk_ult(variables.get(i.v), variables.get(i.w)));
+            }
+            else {
+                fp.add_le(i.v, i.w, index);
+                solver.assert_expr(bv.mk_ule(variables.get(i.v), variables.get(i.w)));
+            }
+        }
+        unsigned v = 0;
+        for (auto const& b : bounds) {
+            ++index;
+            if (b.is_free())
+                continue;
+
+            fp.set_bounds(v, rational(b.lo, rational::ui64()), rational(b.hi, rational::ui64()), index);
+
+            expr_ref A(mk_ult(variables.get(v), bv.mk_numeral(b.hi, 32)), m);
+            expr_ref B(bv.mk_ule(bv.mk_numeral(b.lo, 32), variables.get(v)), m);
+
+            if (b.lo < b.hi) 
+                solver.assert_expr(m.mk_and(A, B));            
+            else if (b.hi == 0)
+                solver.assert_expr(B);
+            else 
+                solver.assert_expr(m.mk_or(A, B));            
+        }
+
+        solver.display(std::cout);
+        std::cout << fp << "\n";
+
+        lbool r1 = solver.check();
+        lbool r2 = fp.make_feasible();
+
+        std::cout << r1 << " " << r2 << "\n" << fp << "\n";
+
+        if (r2 == l_true) {
+            for (auto const& row : rows) {
+                uint64_t sum = 0;
+                for (auto col : row) 
+                    sum += col.second * fp.value(col.first);
+                SASSERT(sum == 0);
+            }
+            for (unsigned i = 0; i < bounds.size(); ++i) {
+                uint64_t val = fp.value(i);
+                uint64_t lo = bounds[i].lo;
+                uint64_t hi = bounds[i].hi;
+                SASSERT(!(lo < hi || hi == 0) || lo <= val && val < hi);
+                SASSERT(!(hi > lo) || val < hi || lo <= val);
+            }
+            for (auto const& i : ineqs) {
+                SASSERT(fp.value(i.v) <= fp.value(i.w));
+                SASSERT(!i.strict || fp.value(i.v) < fp.value(i.w));
+            }
+        }
+        if (r1 == r2) {
+            std::cout << "agree\n";
+        }
+        else if (r1 == l_false && r2 == l_true) {
+            std::cout << "BUG should be unsat\n";
+            solver.display(std::cout);
+            std::cout << fp << "\n";
+        }
+        else if (r1 == l_true && r2 == l_false) {
+            std::cout << "BUG should be sat\n";
+            solver.display(std::cout);
+            std::cout << fp << "\n";
+        }
+        else {
+            std::cout << r2 << " missed with verdict " << r1 << "\n";
+        }
+    }
+
+
+    static void test_lps(random_gen& r, unsigned num_vars, unsigned num_rows, unsigned num_vars_per_row, unsigned num_ineqs) {
+        vector<svector<std::pair<unsigned, uint64_t>>> rows;
+        svector<ineq> ineqs;
+        vector<mod_interval<uint64_t>> bounds;
+        for (unsigned i = 0; i < num_vars; ++i) {
+            uint64_t l = (r() % 2 == 0) ? r() % 100 : (0 - r() % 10);
+            uint64_t h = (r() % 2 == 0) ? r() % 100 : (0 - r() % 10);
+            bounds.push_back(mod_interval<uint64_t>(l, h));
+        }
+        for (unsigned i = 0; i < num_ineqs; ++i) {
+            var_t v = r() % num_vars;
+            var_t w = r() % num_vars;
+            ineqs.push_back(ineq(v, w, 0, 0 == r() % 2));
+        }
+        for (unsigned i = 0; i < num_rows; ++i) {
+            svector<std::pair<unsigned, uint64_t>> row;
+            uint64_t coeff = (r() % 2 == 0) ? r() % 100 : (0 - r() % 10);
+            row.push_back(std::make_pair(i, coeff));
+            for (unsigned j = 0; j + 1 < num_vars_per_row; ++j) {
+                var_t v = (r() % (num_vars - num_vars_per_row)) + num_vars_per_row;
+                coeff = (r() % 2 == 0) ? r() % 100 : (0 - r() % 10);
+                row.push_back(std::make_pair(v, coeff));
+            }
+            rows.push_back(row);
+        }
+
+        test_lp(rows, ineqs, bounds);
+    }
+
+    static void test_lps() {
+        random_gen r;
+        for (unsigned i = 0; i < 10000; ++i)
+            test_lps(r, 6, 0, 0, 5);
+        for (unsigned i = 0; i < 10000; ++i)
+            test_lps(r, 6, 3, 3, 0);
+        for (unsigned i = 0; i < 10000; ++i)
+            test_lps(r, 6, 3, 3, 3);
+
+    }
 }
 
 void tst_fixplex() {
 
+    polysat::test_lps();
+
     polysat::test_ineq_propagation1();
     polysat::test_ineq_propagation2();
     polysat::test_ineqs();