diff --git a/src/ast/sls/sls_arith_base.cpp b/src/ast/sls/sls_arith_base.cpp
index 2f1045ccc..5543ad58d 100644
--- a/src/ast/sls/sls_arith_base.cpp
+++ b/src/ast/sls/sls_arith_base.cpp
@@ -609,8 +609,8 @@ namespace sls {
         for (auto idx : vi.m_muls) {
             auto const& [w, coeff, monomial] = m_muls[idx];
             num_t prod(coeff);
-            for (auto w : monomial)
-                prod *= value(w);
+            for (auto [w, p]:monomial)
+                prod *= power_of(value(w), p);
             if (value(w) != prod) 
                 update(w, prod);            
         }
@@ -715,10 +715,19 @@ namespace sls {
             default: {
                 v = mk_var(e);
                 unsigned idx = m_muls.size();
-                m_muls.push_back({ v, c, m });
+                std::stable_sort(m.begin(), m.end(), [&](unsigned a, unsigned b) { return a < b; });
+                svector<std::pair<unsigned, unsigned>> mp;
+                for (unsigned i = 0; i < m.size(); ++i) {
+                    auto w = m[i];
+                    auto p = 1;
+                    while (i + 1 < m.size() && m[i + 1] == w)
+                        ++p, ++i;
+                    mp.push_back({ w, p });
+                }
+                m_muls.push_back({ v, c, mp });
                 num_t prod(c);
-                for (auto w : m)
-                    m_vars[w].m_muls.push_back(idx), prod *= value(w);
+                for (auto [w, p] : mp)
+                    m_vars[w].m_muls.push_back(idx), prod *= power_of(value(w), p);
                 m_vars[v].m_def_idx = idx;
                 m_vars[v].m_op = arith_op_kind::OP_MUL;
                 m_vars[v].m_value = prod;
@@ -940,8 +949,8 @@ namespace sls {
         case OP_MUL: {
             auto const& [w, coeff, monomial] = m_muls[vi.m_def_idx];
             num_t prod(coeff);
-            for (auto w : monomial)
-                prod *= value(w);
+            for (auto [w, p] : monomial)
+                prod *= power_of(value(w), p);
             update(v, prod);
             break;
         }
@@ -1142,50 +1151,130 @@ namespace sls {
     }
 
     template<typename num_t>
-    bool arith_base<num_t>::repair_square(mul_def const& md) {
+    bool arith_base<num_t>::repair_power(mul_def const& md) {
         auto const& [v, coeff, monomial] = md;
-        if (!is_int(v) || monomial.size() != 2 || monomial[0] != monomial[1])
+        if (!is_int(v))
             return false;
-        
-        num_t val = value(v);
-        val = div(val, coeff);
-        var_t w = monomial[0];
-        if (val < 0) 
-            update(w, num_t(ctx.rand(10)));        
-        else {
-            num_t root = sqrt(val);
-            if (ctx.rand(3) == 0)
+        for (auto [c, v, p] : monomial_iterator(md)) {           
+            num_t val1 = div(value(v), c);
+            if (val1 < 0 && p % 2 == 0)
+                continue;
+            num_t root = root_of(p, val1);
+            if (in_bounds(v, -root) && (((p % 2 == 0) && ctx.rand(3) == 0) || val1 < 0))
                 root = -root;
-            if (root * root == val)
-                update(w, root);
-            else
-                update(w, root + num_t(ctx.rand(3)) - 1);
+            if (update(v, root))
+                return true;
         }
-        verbose_stream() << "ROOT " << val << " v" << w << " := " << value(w) << "\n";
-        return true;
+        return false;
     }
 
+    // increment/decrement the value of one of the variables as long
+    // as it improves the repair distance.
     template<typename num_t>
-    bool arith_base<num_t>::repair_mul1(mul_def const& md) {
+    bool arith_base<num_t>::repair_mul_inc(mul_def const& md) {
+        num_t inc;
+        unsigned i = 0;
+        double sum_probs = 0;
+        unsigned sz = md.m_monomial.size();
+        m_probs.reserve(sz);
+        for (auto [c, v, p] : monomial_iterator(md)) {
+            int r = -1; // repair_delta(v, inc);
+            auto prob = r <= 0 ? 0.1 : r;
+            sum_probs += prob;
+            m_probs[i++] = prob;
+
+        }
+        i = sz;
+        double lim = sum_probs * ((double)ctx.rand() / random_gen().max_value());
+        do {
+            lim -= m_probs[--i];
+        } while (lim >= 0 && i > 0);
+
+        return false;
+    }
+
+    // update entries to +/- 1, except for one, which is set to be +/- value(v)
+    template<typename num_t>
+    bool arith_base<num_t>::repair_mul_ones(mul_def const& md) {
+        auto const& [v, coeff, monomial] = md;
+        auto val = value(v);
+        vector<num_t> coeffs(monomial.size(), num_t(1));
+        for (auto [c, v, p] : monomial_iterator(md)) {
+            if (c == 0)
+                continue;
+            auto new_value = root_of(p, abs(val));
+            int sign = ctx.rand(2) == 0 ? 1 : -1;
+            if (sign == -1)
+                new_value = -new_value;
+            if (!in_bounds(v, new_value)) {
+                new_value = -new_value;
+                if (p % 2 != 0)
+                    sign *= -1;
+            }
+            if (!in_bounds(v, new_value))
+                continue;
+            
+            unsigned i = 0;
+            bool failed = false;
+            for (auto [w, q] : monomial) {
+                if (w == v) 
+                    coeffs[i] = new_value;                                   
+                else if (q % 2 == 0) 
+                    coeffs[i] = num_t(ctx.rand(2) == 0 ? 1 : -1);
+                else {
+                    auto sign1 = ctx.rand(2) == 0 ? 1 : -1;
+                    coeffs[i] = num_t(sign1);
+                    if (!in_bounds(w, coeffs[i])) {
+                        sign1 = -sign1;
+                        coeffs[i] = num_t(sign1);
+                    }
+                    if (!in_bounds(w, coeffs[i]))
+                        failed = true;
+                    sign *= sign1;
+                }
+                ++i;
+            }
+            if (failed)
+                continue;
+
+            if ((sign == 1) != val > 0) {
+                bool fixed = false;
+                for (auto [w, q] : monomial) {
+                    if (q % 2 == 1 && in_bounds(w, -coeffs[i])) {
+                        coeffs[i] = -coeffs[i];
+                        fixed = true;
+                        break;
+                    }
+                }
+                if (!fixed)
+                    continue;
+            }           
+            i = 0;
+            for (auto [w, q] : monomial) {
+                if (!update(w, coeffs[i]))
+                    return false;
+                ++i;
+            }
+            return true;
+        }
+        return false;
+    }
+
+    // update one of the variables such that the product aligns with value(v)
+    template<typename num_t>
+    bool arith_base<num_t>::repair_mul_one(mul_def const& md) {
         auto const& [v, coeff, monomial] = md;
         if (!is_int(v))
             return false;
         num_t val = value(v);
-        val = div(val, coeff);
-        if (val == 0)
+        if (!divides(coeff, val))
             return false;
-        unsigned sz = monomial.size();
-        unsigned start = ctx.rand(sz);
-        for (unsigned i = 0; i < sz; ++i) {
-            unsigned j = (start + i) % sz;
-            auto w = monomial[j];
-            num_t product(1);
-            for (auto v : monomial)
-                if (v != w)
-                    product *= value(v);
-            if (product == 0 || !divides(product, val))
+        for (auto [c, v, p] : monomial_iterator(md)) {
+            if (c == 0 || !divides(c, val))
                 continue;
-            if (update(w, div(val, product)))
+            auto val1 = div(val, c);
+            val1 = root_of(p, val);
+            if (update(v, val1))
                 return true;
         }
         return false;
@@ -1197,8 +1286,8 @@ namespace sls {
         num_t product(coeff);
         num_t val = value(v);
         num_t new_value;
-        for (auto v : monomial)
-            product *= value(v);
+        for (auto [v, p]: monomial)
+            product *= power_of(value(v), p);
         if (product == val)
             return true;
         verbose_stream() << "repair mul " << mk_bounded_pp(m_vars[v].m_expr, m) << " := " << val << "(product: " << product << ")\n";
@@ -1206,72 +1295,72 @@ namespace sls {
         if (ctx.rand(20) == 0)
             return update(v, product);
         else if (val == 0) {
-            auto v = monomial[ctx.rand(sz)];
-            num_t zero(0);
-            return update(v, zero);
+            for (auto [c, v, p] : monomial_iterator(md))
+                if (update(v, num_t(0)))
+                    return true;
+            return false;
         }
-        else if (repair_square(md))
+        else if (abs(val) == 1 && repair_mul_one(md))
             return true;
-        else if (ctx.rand(4) != 0 && repair_mul1(md)) {
-#if 0
-            verbose_stream() << "mul1 " << val << " " << coeff << " ";
-            for (auto v : monomial)
-                verbose_stream() << "v" << v << " = " << value(v) << " ";
-            verbose_stream() << "\n";
-#endif
+        else if (repair_power(md))
             return true;
-        }
-        else if (is_int(v)) {
-#if 0
-            verbose_stream() << "repair mul2 - ";
-            for (auto v : monomial)
-                verbose_stream() << "v" << v << " = " << value(v) << " ";
-#endif
-            num_t n = div(val, coeff);
-            if (!divides(coeff, val) && ctx.rand(2) == 0)
-                n = div(val + coeff - 1, coeff);
-            auto const& fs = factor(abs(n));
-            vector<num_t> coeffs(sz, num_t(1));
-            vector<num_t> gcds(sz, num_t(0));
-            num_t sign(1);
-            for (auto c : coeffs)
-                sign *= c;
-            unsigned i = 0;            
-            for (auto w : monomial) {
-                for (auto idx : m_vars[w].m_muls) {
-                    auto const& [w1, coeff1, monomial1] = m_muls[idx];
-                    gcds[i] = gcd(gcds[i], abs(value(w1)));
-                }
-                auto const& vi = m_vars[w];
-                if (vi.m_lo && vi.m_lo->value >= 0)
-                    coeffs[i] = 1;
-                else if (vi.m_hi && vi.m_hi->value < 0)
-                    coeffs[i] = -1;
-                else 
-                    coeffs[i] = num_t(ctx.rand(2) == 0 ? 1 : -1);
-                ++i;
-            }
-            for (auto f : fs) 
-                coeffs[ctx.rand(sz)] *= f;
-            if ((sign == 0) != (n == 0))
-                coeffs[ctx.rand(sz)] *= -1;
-            verbose_stream() << "value " << val << " coeff: " << coeff << " coeffs: " << coeffs << " factors: " << fs << "\n";
-            i = 0;
-            for (auto w : monomial) {
-                if (!update(w, coeffs[i++])) {
-                    verbose_stream() << "failed to update v" << w << " to " << coeffs[i - 1] << "\n";
-                    return false;
-                }
-            }
-            verbose_stream() << "all updated for v" << v << " := " << value(v) << "\n";
+        else if (ctx.rand(4) != 0 && repair_mul_one(md)) 
+            return true;        
+        else if (repair_mul_factors(md))
             return true;
-        }
         else {
             NOT_IMPLEMENTED_YET();
         }
         return false;
     }
 
+    template<typename num_t>
+    bool arith_base<num_t>::repair_mul_factors(mul_def const& md) {
+        auto const& [v, coeff, monomial] = md;
+        auto val = value(v);
+        if (!is_int(v))
+            return false;
+        num_t n = div(val, coeff);
+        if (!divides(coeff, val) && ctx.rand(2) == 0)
+            n = div(val + coeff - 1, coeff);
+        auto const& fs = factor(abs(n));
+        unsigned sz = monomial.size();
+        vector<num_t> coeffs(sz, num_t(1));
+        vector<num_t> gcds(sz, num_t(0));
+        num_t sign(1);
+        for (auto c : coeffs)
+            sign *= c;
+        unsigned i = 0;
+        for (auto [w, p] : monomial) {
+            for (auto idx : m_vars[w].m_muls) {
+                auto const& [w1, coeff1, monomial1] = m_muls[idx];
+                gcds[i] = gcd(gcds[i], abs(value(w1)));
+            }
+            auto const& vi = m_vars[w];
+            if (vi.m_lo && vi.m_lo->value >= 0)
+                coeffs[i] = 1;
+            else if (vi.m_hi && vi.m_hi->value < 0)
+                coeffs[i] = -1;
+            else
+                coeffs[i] = num_t(ctx.rand(2) == 0 ? 1 : -1);
+            ++i;
+        }
+        for (auto f : fs)
+            coeffs[ctx.rand(sz)] *= f;
+        if ((sign == 0) != (n == 0))
+            coeffs[ctx.rand(sz)] *= -1;
+        verbose_stream() << "value " << val << " coeff: " << coeff << " coeffs: " << coeffs << " factors: " << fs << "\n";
+        i = 0;
+        for (auto [w, p] : monomial) {
+            if (!update(w, coeffs[i++])) {
+                verbose_stream() << "failed to update v" << w << " to " << coeffs[i - 1] << "\n";
+                return false;
+            }
+        }
+        verbose_stream() << "all updated for v" << v << " := " << value(v) << "\n";
+        return true;
+    }
+
     template<typename num_t>
     bool arith_base<num_t>::repair_rem(op_def const& od) {
         auto v1 = value(od.m_arg1);
@@ -1438,34 +1527,37 @@ namespace sls {
         return max_result;
     }
 
-#if 0
-    double sum_prob = 0;
-    unsigned i = 0;
-    clause const& c = get_clause(cls_idx);
-    for (literal lit : c) {
-        double prob = m_prob_break[m_breaks[lit.var()]];
-        m_probs[i++] = prob;
-        sum_prob += prob;
-    }
-    double lim = sum_prob * ((double)m_rand() / m_rand.max_value());
-    do {
-        lim -= m_probs[--i];
-    } while (lim >= 0 && i > 0);
-#endif
-
-    // Newton function for integer square root.
     template<typename num_t>
-    num_t arith_base<num_t>::sqrt(num_t n) {
-        if (n <= 1)
-            return n;
+    num_t arith_base<num_t>::power_of(num_t x, unsigned k) {
+        num_t r = x;
+        while (k > 1) {
+            if (k % 2 == 1) {
+                r = x * r;
+                --k;
+            }
+            x = x * x;
+            k /= 2;
+        }
+        return x * r;
+    }
 
-        auto x0 = div(n, num_t(2));
+    // Newton function for integer n'th root of a
+    // x_{k+1} = 1/k ((k-1)*x_k + a / x_k^{n-1})
+    template<typename num_t>
+    num_t arith_base<num_t>::root_of(unsigned k, num_t a) {
+        if (a <= 1)
+            return a;
+        if (k == 1)
+            return a;
+        SASSERT(k > 1);
+       
+        auto x0 = div(a, num_t(k));
 
-        auto x1 = div(x0 + div(n, x0), num_t(2));
+        auto x1 = div((x0 * num_t(k - 1)) + div(a, power_of(x0, k - 1)), num_t(k));
 
         while (x1 < x0)	{
             x0 = x1;
-            x1 = div(x0 + div(n, x0), num_t(2));
+            x1 = div((x0 * num_t(k - 1)) + div(a, power_of(x0, k - 1)), num_t(k));
         }
         return x0;
     }
@@ -1663,8 +1755,13 @@ namespace sls {
         
         for (auto md : m_muls) {
             out << "v" << md.m_var << " := ";
-            for (auto w : md.m_monomial)
-                out << "v" << w << " ";
+            for (auto [w, p] : md.m_monomial) {
+                out << "v" << w;
+                if (p > 1)
+                    out << "^" << p;
+                out << " ";
+            }
+
             out << "\n";
         }
         for (auto ad : m_adds) {
@@ -1683,6 +1780,7 @@ namespace sls {
         return out;
     }
 
+
     template<typename num_t>
     void arith_base<num_t>::invariant() {
         for (unsigned v = 0; v < ctx.num_bool_vars(); ++v) {
@@ -1694,14 +1792,18 @@ namespace sls {
         for (auto md : m_muls) {
             auto const& [w, coeff, monomial] = md;
             num_t prod(coeff);
-            for (auto v : monomial)
-                prod *= value(v);
+            for (auto [v,p] : monomial)
+                prod *= power_of(value(v), p);
             //verbose_stream() << "check " << w << " " << monomial << "\n";
             if (prod != value(w)) {
                 out << prod << " " << value(w) << "\n";
                 out << "v" << w << " := ";
-                for (auto w : monomial)
-                    out << "v" << w << " ";
+                for (auto [w, p] : monomial) {
+                    out << "v" << w;
+                    if (p > 1)
+                        out << "^" << p;
+                    out << " ";
+                }
                 out << "\n";
             }
             SASSERT(prod == value(w));
diff --git a/src/ast/sls/sls_arith_base.h b/src/ast/sls/sls_arith_base.h
index 488101f51..a3dba2791 100644
--- a/src/ast/sls/sls_arith_base.h
+++ b/src/ast/sls/sls_arith_base.h
@@ -81,7 +81,7 @@ namespace sls {
         struct mul_def {
             unsigned        m_var;
             num_t           m_coeff;
-            unsigned_vector m_monomial;
+            svector<std::pair<unsigned, unsigned>> m_monomial;
         };
 
         struct add_def : public linear_term {
@@ -111,8 +111,11 @@ namespace sls {
 
         unsigned get_num_vars() const { return m_vars.size(); }
 
-        bool repair_mul1(mul_def const& md);
-        bool repair_square(mul_def const& md);
+        bool repair_mul_one(mul_def const& md);
+        bool repair_power(mul_def const& md);
+        bool repair_mul_factors(mul_def const& md);
+        bool repair_mul_ones(mul_def const& md);
+        bool repair_mul_inc(mul_def const& md);        
         bool repair_mul(mul_def const& md);
         bool repair_add(add_def const& ad);
         bool repair_mod(op_def const& od);
@@ -130,7 +133,46 @@ namespace sls {
 
         vector<num_t> m_factors;
         vector<num_t> const& factor(num_t n);
-        num_t sqrt(num_t n);
+        num_t root_of(unsigned n, num_t a);
+        num_t power_of(num_t a, unsigned k);
+
+        struct monomial_elem {
+            num_t other_product;
+            var_t v;
+            unsigned p; // power
+        };
+
+        struct monomial_iterator_base {
+            arith_base& a;
+            mul_def const& md;
+            unsigned m_seed = 0;
+            unsigned m_idx;
+            monomial_iterator_base(arith_base& a, mul_def const& md, unsigned idx, unsigned seed) : a(a), md(md), m_seed(seed), m_idx(idx) {}
+            monomial_elem operator*() {
+                auto [v, p] = md.m_monomial[(m_idx + m_seed) % md.m_monomial.size()];
+                num_t prod(md.m_coeff);
+                for (auto [w, q] : md.m_monomial)
+                    if (w != v)
+                        prod *= a.power_of(a.value(w), q);
+                return { prod, v, p };
+            }
+            monomial_iterator_base& operator++() {
+                ++m_idx;
+                return *this;
+            }
+            bool operator==(monomial_iterator_base const& other) const { return m_idx == other.m_idx; }
+            bool operator!=(monomial_iterator_base const& other) const { return m_idx != other.m_idx; }
+        };
+
+        struct monomials {
+            arith_base& a;
+            mul_def const& md;
+            monomials(arith_base & a, mul_def const& md) : a(a), md(md) {}
+            monomial_iterator_base begin() { return monomial_iterator_base(a, md, 0, a.ctx.rand()); }
+            monomial_iterator_base end() { return monomial_iterator_base(a, md, md.m_monomial.size(), 0); }
+        };
+
+        monomials monomial_iterator(mul_def const& md) { return monomials(*this, md); }
 
         double reward(sat::literal lit);