diff --git a/examples/python/data/horn1.smt2 b/examples/python/data/horn1.smt2
new file mode 100644
index 000000000..20d043534
--- /dev/null
+++ b/examples/python/data/horn1.smt2
@@ -0,0 +1,50 @@
+(declare-rel Goal (Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool))
+(declare-rel Invariant (Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool))
+(declare-var A Bool)
+(declare-var B Bool)
+(declare-var C Bool)
+(declare-var D Bool)
+(declare-var E Bool)
+(declare-var F Bool)
+(declare-var G Bool)
+(declare-var H Bool)
+(declare-var I Bool)
+(declare-var J Bool)
+(declare-var K Bool)
+(declare-var L Bool)
+(declare-var M Bool)
+(declare-var N Bool)
+(declare-var O Bool)
+(declare-var P Bool)
+(declare-var Q Bool)
+(declare-var R Bool)
+(declare-var S Bool)
+(declare-var T Bool)
+(declare-var U Bool)
+(declare-var V Bool)
+(declare-var W Bool)
+(declare-var X Bool)
+(rule (=> (not (or L K J I H G F E D C B A)) (Invariant L K J I H G F E D C B A)))
+(rule (let ((a!1 (and (Invariant X W V U T S R Q P O N M)
+                (=> (not (and true)) (not F))
+                (=> (not (and true)) (not E))
+                (=> (not (and W)) (not D))
+                (=> (not (and W)) (not C))
+                (=> (not (and U)) (not B))
+                (=> (not (and U)) (not A))
+                (= L (xor F X))
+                (= K (xor E W))
+                (= J (xor D V))
+                (= I (xor C U))
+                (= H (xor B T))
+                (= G (xor A S))
+                (=> D (not E))
+                (=> C (not E))
+                (=> B (not C))
+                (=> A (not C))
+                ((_ at-most 5) L K J I H G))))
+  (=> a!1 (Invariant L K J I H G F E D C B A))))
+(rule (=> (and (Invariant L K J I H G F E D C B A) L (not K) J (not I) H G)
+    (Goal L K J I H G F E D C B A)))
+
+(query Goal)
\ No newline at end of file
diff --git a/examples/python/data/horn2.smt2 b/examples/python/data/horn2.smt2
new file mode 100644
index 000000000..478c39d5f
--- /dev/null
+++ b/examples/python/data/horn2.smt2
@@ -0,0 +1,44 @@
+(declare-rel Invariant (Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool))
+(declare-rel Goal (Bool Bool Bool Bool Bool Bool Bool Bool Bool Bool))
+(declare-var A Bool)
+(declare-var B Bool)
+(declare-var C Bool)
+(declare-var D Bool)
+(declare-var E Bool)
+(declare-var F Bool)
+(declare-var G Bool)
+(declare-var H Bool)
+(declare-var I Bool)
+(declare-var J Bool)
+(declare-var K Bool)
+(declare-var L Bool)
+(declare-var M Bool)
+(declare-var N Bool)
+(declare-var O Bool)
+(declare-var P Bool)
+(declare-var Q Bool)
+(declare-var R Bool)
+(declare-var S Bool)
+(declare-var T Bool)
+(rule (=> (not (or J I H G F E D C B A)) (Invariant J I H G F E D C B A)))
+(rule (let ((a!1 (and (Invariant T S R Q P O N M L K)
+                (=> (not (and true)) (not E))
+                (=> (not (and T)) (not D))
+                (=> (not (and S)) (not C))
+                (=> (not (and R)) (not B))
+                (=> (not (and Q)) (not A))
+                (= J (xor E T))
+                (= I (xor D S))
+                (= H (xor C R))
+                (= G (xor B Q))
+                (= F (xor A P))
+                (=> D (not E))
+                (=> C (not D))
+                (=> B (not C))
+                (=> A (not B))
+                ((_ at-most 3) J I H G F))))
+  (=> a!1 (Invariant J I H G F E D C B A))))
+(rule (=> (and (Invariant J I H G F E D C B A) (not J) (not I) (not H) (not G) F)
+    (Goal J I H G F E D C B A)))
+
+(query Goal)
diff --git a/examples/python/mini_ic3.py b/examples/python/mini_ic3.py
new file mode 100644
index 000000000..5a8ea566b
--- /dev/null
+++ b/examples/python/mini_ic3.py
@@ -0,0 +1,438 @@
+from z3 import *
+import heapq
+
+
+# Simplistic (and fragile) converter from
+# a class of Horn clauses corresponding to 
+# a transition system into a transition system
+# representation as <init, trans, goal>
+# It assumes it is given three Horn clauses
+# of the form:
+#  init(x) => Invariant(x)
+#  Invariant(x) and trans(x,x') => Invariant(x')
+#  Invariant(x) and goal(x) => Goal(x)
+# where Invariant and Goal are uninterpreted predicates
+    
+class Horn2Transitions:
+    def __init__(self):
+        self.trans = True
+        self.init = True
+        self.goal = True
+        self.index = 0
+        
+    def parse(self, file):
+        fp = Fixedpoint()
+        goals = fp.parse_file(file)
+        for r in fp.get_rules():
+            if not is_quantifier(r):
+                continue
+            b = r.body()
+            if not is_implies(b):
+                continue
+            f = b.arg(0)
+            g = b.arg(1)
+            if self.is_goal(f, g):
+                continue
+            if self.is_transition(f, g):
+                continue
+            if self.is_init(f, g):
+                continue
+
+    def is_pred(self, p, name):
+        return is_app(p) and p.decl().name() == name
+    
+    def is_goal(self, body, head):
+        if not self.is_pred(head, "Goal"):
+            return False
+        pred, inv = self.is_body(body)
+        if pred is None:
+            return False
+        self.goal = self.subst_vars("x", inv, pred)
+        return True
+
+    def is_body(self, body):
+        if not is_and(body):
+            return None, None
+        fmls = [f for f in body.children() if self.is_inv(f) is None]
+        inv = None
+        for f in body.children():
+            if self.is_inv(f) is not None:
+                inv = f;
+                break
+        return And(fmls), inv
+
+    def is_inv(self, f):
+        if self.is_pred(f, "Invariant"):
+            return f
+        return None
+
+    def is_transition(self, body, head):
+        pred, inv0 = self.is_body(body)
+        if pred is None:
+            return False
+        inv1 = self.is_inv(head)
+        if inv1 is None:
+            return False
+        pred = self.subst_vars("x",  inv0, pred)
+        self.xs = self.vars
+        pred = self.subst_vars("xn", inv1, pred)
+        self.xns = self.vars
+        self.trans = pred
+        return True
+
+    def is_init(self, body, head):
+        for f in body.children():
+            if self.is_inv(f) is not None:
+               return False
+        inv = self.is_inv(head)
+        if inv is None:
+            return False
+        self.init = self.subst_vars("x", inv, body)
+        return True
+    
+    def subst_vars(self, prefix, inv, fml):
+        subst = self.mk_subst(prefix, inv)
+        self.vars = [ v for (k,v) in subst ]
+        return substitute(fml, subst)
+
+    def mk_subst(self, prefix, inv):
+        self.index = 0
+        return [(f, self.mk_bool(prefix)) for f in inv.children()]
+
+    def mk_bool(self, prefix):
+        self.index += 1
+        return Bool("%s%d" % (prefix, self.index))
+
+# Produce a finite domain solver.
+# The theory QF_FD covers bit-vector formulas
+# and pseudo-Boolean constraints.
+# By default cardinality and pseudo-Boolean 
+# constraints are converted to clauses. To override
+# this default for cardinality constraints
+# we set sat.cardinality.solver to True
+
+def fd_solver():
+    s = SolverFor("QF_FD")
+    s.set("sat.cardinality.solver", True)
+    return s
+
+
+# negate, avoid double negation
+def negate(f):
+    if is_not(f):
+        return f.arg(0)
+    else:
+        return Not(f)
+
+def cube2clause(cube):
+    return Or([negate(f) for f in cube])
+
+class State:
+    def __init__(self, s):
+        self.R = set([])
+        self.solver = s
+
+    def add(self, clause):
+        if clause not in self.R:
+           self.R |= { clause }
+           self.solver.add(clause)
+    
+class Goal:
+    def __init__(self, cube, parent, level):
+        self.level = level
+        self.cube = cube
+        self.parent = parent
+
+def is_seq(f):
+    return isinstance(f, list) or isinstance(f, tuple) or isinstance(f, AstVector)
+
+# Check if the initial state is bad
+def check_disjoint(a, b):
+    s = fd_solver()
+    s.add(a)
+    s.add(b)
+    return unsat == s.check()
+
+
+# Remove clauses that are subsumed
+def prune(R):
+    removed = set([])
+    s = fd_solver()
+    for f1 in R:
+        s.push()
+        for f2 in R:
+            if f2 not in removed:
+               s.add(Not(f2) if f1.eq(f2) else f2)
+        if s.check() == unsat:
+            removed |= { f1 }
+        s.pop()
+    return R - removed
+                
+class MiniIC3:
+    def __init__(self, init, trans, goal, x0, xn):
+        self.x0 = x0
+        self.xn = xn
+        self.init = init
+        self.bad = goal
+        self.trans = trans
+        self.min_cube_solver = fd_solver()
+        self.min_cube_solver.add(Not(trans))
+        self.goals = []
+        s = State(fd_solver())
+        s.add(init)
+        s.solver.add(trans)
+        self.states = [s]
+        self.s_bad = fd_solver()
+        self.s_good = fd_solver()
+        self.s_bad.add(self.bad)
+        self.s_good.add(Not(self.bad))        
+        
+    def next(self, f):
+        if is_seq(f):
+           return [self.next(f1) for f1 in f]
+        return substitute(f, zip(self.x0, self.xn))    
+    
+    def prev(self, f):
+        if is_seq(f):
+           return [self.prev(f1) for f1 in f]
+        return substitute(f, zip(self.xn, self.x0))    
+    
+    def add_solver(self):
+        s = fd_solver()
+        s.add(self.trans)
+        self.states += [State(s)]        
+
+    def R(self, i):
+        return And(self.states[i].R)
+
+    # Check if there are two states next to each other that have the same clauses.
+    def is_valid(self):
+        i = 1
+        while i + 1 < len(self.states):
+            if not (self.states[i].R - self.states[i+1].R):
+               return And(prune(self.states[i].R))
+            i += 1
+        return None
+
+    def value2literal(self, m, x):
+        value = m.eval(x)
+        if is_true(value):
+            return x
+        if is_false(value):
+            return Not(x)
+        return None
+
+    def values2literals(self, m, xs):
+        p = [self.value2literal(m, x) for x in xs]
+        return [x for x in p if x is not None]
+
+    def project0(self, m):
+        return self.values2literals(m, self.x0)
+
+    def projectN(self, m):
+        return self.values2literals(m, self.xn)
+
+    # Determine if there is a cube for the current state 
+    # that is potentially reachable.
+    def unfold(self):
+        core = []
+        self.s_bad.push()
+        R = self.R(len(self.states)-1)
+        self.s_bad.add(R)
+        is_sat = self.s_bad.check()
+        if is_sat == sat:
+           m = self.s_bad.model()
+           props = self.project0(m)
+           self.s_good.push()
+           self.s_good.add(R)
+           is_sat2 = self.s_good.check(props)
+           assert is_sat2 == unsat
+           core = self.s_good.unsat_core()
+           self.s_good.pop()
+        self.s_bad.pop()
+        return is_sat, core
+
+    # Block a cube by asserting the clause corresponding to its negation
+    def block_cube(self, i, cube):
+        self.assert_clause(i, cube2clause(cube))
+
+    # Add a clause to levels 0 until i
+    def assert_clause(self, i, clause):
+        for j in range(i + 1):
+            self.states[j].add(clause)
+
+    # minimize cube that is core of Dual solver.
+    # this assumes that props & cube => Trans    
+    def minimize_cube(self, cube, lits):
+        is_sat = self.min_cube_solver.check(lits + [c for c in cube])
+        assert is_sat == unsat
+        core = self.min_cube_solver.unsat_core()
+        assert core
+        return [c for c in core if c in set(cube)]
+
+    # push a goal on a heap
+    def push_heap(self, goal):
+        heapq.heappush(self.goals, (goal.level, goal))
+
+    # A state s0 and level f0 such that
+    # not(s0) is f0-1 inductive
+    def ic3_blocked(self, s0, f0):
+        self.push_heap(Goal(self.next(s0), None, f0))
+        while self.goals:
+            f, g = heapq.heappop(self.goals)
+            sys.stdout.write("%d." % f)
+            sys.stdout.flush()
+            # Not(g.cube) is f-1 invariant
+            if f == 0:
+               print("")
+               return g
+            cube, f, is_sat = self.is_inductive(f, g.cube)
+            if is_sat == unsat:
+               self.block_cube(f, self.prev(cube))
+               if f < f0:
+                  self.push_heap(Goal(g.cube, g.parent, f + 1))
+            elif is_sat == sat:
+               self.push_heap(Goal(cube, g, f - 1))
+               self.push_heap(g)
+            else:
+               return is_sat
+        print("")
+        return None
+
+    # Rudimentary generalization:
+    # If the cube is already unsat with respect to transition relation
+    # extract a core (not necessarily minimal)
+    # otherwise, just return the cube.
+    def generalize(self, cube, f):
+        s = self.states[f - 1].solver
+        if unsat == s.check(cube):
+            return s.unsat_core(), f
+        return cube, f
+
+    # Check if the negation of cube is inductive at level f
+    def is_inductive(self, f, cube):
+        s = self.states[f - 1].solver
+        s.push()
+        s.add(self.prev(Not(And(cube))))
+        is_sat = s.check(cube)
+        if is_sat == sat:
+           m = s.model()
+        s.pop()
+        if is_sat == sat:
+           cube = self.next(self.minimize_cube(self.project0(m), self.projectN(m)))
+        elif is_sat == unsat:
+           cube, f = self.generalize(cube, f)
+        return cube, f, is_sat
+                        
+    def run(self):
+        if not check_disjoint(self.init, self.bad):
+           return "goal is reached in initial state"
+        level = 0
+        while True:
+            inv = self.is_valid()
+            if inv is not None:
+                return inv
+            is_sat, cube = self.unfold()
+            if is_sat == unsat:
+               level += 1
+               print("Unfold %d" % level)
+               sys.stdout.flush()
+               self.add_solver()
+            elif is_sat == sat:
+               cex = self.ic3_blocked(cube, level)
+               if cex is not None:
+                  return cex
+            else:
+               return is_sat  
+
+def test(file):
+    h2t = Horn2Transitions()
+    h2t.parse(file)
+    mp = MiniIC3(h2t.init, h2t.trans, h2t.goal, h2t.xs, h2t.xns)
+    result = mp.run()    
+    if isinstance(result, Goal):
+       g = result
+       print("Trace")
+       while g:
+          print(g.level, g.cube)
+          g = g.parent
+       return
+    if isinstance(result, ExprRef):
+       print("Invariant:\n%s " % result)
+       return
+    print(result)
+
+test("data/horn1.smt2")
+test("data/horn2.smt2")
+
+
+
+"""
+# TBD: Quip variant of IC3
+
+must = True
+may = False
+
+class QGoal:
+    def __init__(self, cube, parent, level, must):
+        self.level = level
+        self.cube = cube
+        self.parent = parent
+        self.must = must
+
+class Quip(MiniIC3):
+
+    # prev & tras -> r', such that r' intersects with cube
+    def add_reachable(self, prev, cube):
+        s = fd_solver()
+        s.add(self.trans)
+        s.add(prev)
+        s.add(Or(cube))
+        is_sat = s.check()
+        assert is_sat == sat
+        m = s.model();
+        result = self.values2literals(m, cube)
+        assert result
+        self.reachable.add(result)
+
+    # A state s0 and level f0 such that
+    # not(s0) is f0-1 inductive
+    def quip_blocked(self, s0, f0):
+        self.push_heap(QGoal(self.next(s0), None, f0, must))
+        while self.goals:
+           f, g = heapq.heappop(self.goals)
+           sys.stdout.write("%d." % f)
+           sys.stdout.flush()
+           if f == 0:
+              if g.must:
+                 print("")
+                 return g
+              self.add_reachable(self.init, p.parent.cube)
+              continue
+
+        # TBD
+        return None
+
+                        
+    def run(self):
+        if not check_disjoint(self.init, self.bad):
+           return "goal is reached in initial state"
+        level = 0
+        while True:
+            inv = self.is_valid()
+            if inv is not None:
+                return inv
+            is_sat, cube = self.unfold()
+            if is_sat == unsat:
+               level += 1
+               print("Unfold %d" % level)
+               sys.stdout.flush()
+               self.add_solver()
+            elif is_sat == sat:
+               cex = self.quipie_blocked(cube, level)
+               if cex is not None:
+                  return cex
+            else:
+               return is_sat  
+
+"""
\ No newline at end of file
diff --git a/scripts/mk_project.py b/scripts/mk_project.py
index ca62f5c5f..1ec5f05b5 100644
--- a/scripts/mk_project.py
+++ b/scripts/mk_project.py
@@ -10,7 +10,7 @@ from mk_util import *
 # Z3 Project definition
 def init_project_def():
     set_version(4, 8, 0, 0)
-    add_lib('util', [])
+    add_lib('util', [], includes2install = ['z3_version.h'])
     add_lib('polynomial', ['util'], 'math/polynomial')
     add_lib('sat', ['util'])
     add_lib('nlsat', ['polynomial', 'sat'])
diff --git a/scripts/mk_util.py b/scripts/mk_util.py
index 770e118ee..ebe017739 100644
--- a/scripts/mk_util.py
+++ b/scripts/mk_util.py
@@ -2805,8 +2805,8 @@ def get_full_version_string(major, minor, build, revision):
 # Update files with the version number
 def mk_version_dot_h(major, minor, build, revision):
     c = get_component(UTIL_COMPONENT)
-    version_template = os.path.join(c.src_dir, 'version.h.in')
-    version_header_output = os.path.join(c.src_dir, 'version.h')
+    version_template = os.path.join(c.src_dir, 'z3_version.h.in')
+    version_header_output = os.path.join(c.src_dir, 'z3_version.h')
     # Note the substitution names are what is used by the CMake
     # builds system. If you change these you should change them
     # in the CMake build too
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 826f87e8c..c497c19ee 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -166,6 +166,8 @@ foreach (header ${libz3_public_headers})
   set_property(TARGET libz3 APPEND PROPERTY
     PUBLIC_HEADER "${CMAKE_SOURCE_DIR}/src/api/${header}")
 endforeach()
+set_property(TARGET libz3 APPEND PROPERTY
+    PUBLIC_HEADER "${CMAKE_CURRENT_BINARY_DIR}/util/z3_version.h")
 
 install(TARGETS libz3
   EXPORT Z3_EXPORTED_TARGETS
diff --git a/src/api/api_context.cpp b/src/api/api_context.cpp
index cc2a13aed..4b3b85399 100644
--- a/src/api/api_context.cpp
+++ b/src/api/api_context.cpp
@@ -19,7 +19,7 @@ Revision History:
 --*/
 #include<typeinfo>
 #include "api/api_context.h"
-#include "util/version.h"
+#include "util/z3_version.h"
 #include "ast/ast_pp.h"
 #include "ast/ast_ll_pp.h"
 #include "api/api_log_macros.h"
diff --git a/src/api/api_log.cpp b/src/api/api_log.cpp
index 1bdbb8735..0f531b98e 100644
--- a/src/api/api_log.cpp
+++ b/src/api/api_log.cpp
@@ -19,7 +19,7 @@ Revision History:
 #include "api/z3.h"
 #include "api/api_log_macros.h"
 #include "util/util.h"
-#include "util/version.h"
+#include "util/z3_version.h"
 
 std::ostream * g_z3_log = nullptr;
 bool g_z3_log_enabled   = false;
diff --git a/src/api/api_solver.cpp b/src/api/api_solver.cpp
index 9ad51aaf4..5a4537a4a 100644
--- a/src/api/api_solver.cpp
+++ b/src/api/api_solver.cpp
@@ -179,6 +179,7 @@ extern "C" {
         LOG_Z3_solver_from_file(c, s, file_name);
         char const* ext = get_extension(file_name);
         std::ifstream is(file_name);
+        init_solver(c, s);
         if (!is) {
             SET_ERROR_CODE(Z3_FILE_ACCESS_ERROR, nullptr);
         }
@@ -371,6 +372,21 @@ extern "C" {
         Z3_CATCH_RETURN(0);
     }
 
+    Z3_ast_vector Z3_API Z3_solver_get_non_units(Z3_context c, Z3_solver s) {
+        Z3_TRY;
+        LOG_Z3_solver_get_non_units(c, s);
+        RESET_ERROR_CODE();
+        init_solver(c, s);
+        Z3_ast_vector_ref * v = alloc(Z3_ast_vector_ref, *mk_c(c), mk_c(c)->m());
+        mk_c(c)->save_object(v);
+        expr_ref_vector fmls = to_solver_ref(s)->get_non_units(mk_c(c)->m());
+        for (expr* f : fmls) {
+            v->m_ast_vector.push_back(f);
+        }
+        RETURN_Z3(of_ast_vector(v));
+        Z3_CATCH_RETURN(0);
+    }
+
     static Z3_lbool _solver_check(Z3_context c, Z3_solver s, unsigned num_assumptions, Z3_ast const assumptions[]) {
         for (unsigned i = 0; i < num_assumptions; i++) {
             if (!is_expr(to_ast(assumptions[i]))) {
diff --git a/src/api/c++/z3++.h b/src/api/c++/z3++.h
index e1f263e17..a283571d0 100644
--- a/src/api/c++/z3++.h
+++ b/src/api/c++/z3++.h
@@ -2038,6 +2038,8 @@ namespace z3 {
         stats statistics() const { Z3_stats r = Z3_solver_get_statistics(ctx(), m_solver); check_error(); return stats(ctx(), r); }
         expr_vector unsat_core() const { Z3_ast_vector r = Z3_solver_get_unsat_core(ctx(), m_solver); check_error(); return expr_vector(ctx(), r); }
         expr_vector assertions() const { Z3_ast_vector r = Z3_solver_get_assertions(ctx(), m_solver); check_error(); return expr_vector(ctx(), r); }
+        expr_vector non_units() const { Z3_ast_vector r = Z3_solver_get_non_units(ctx(), m_solver); check_error(); return expr_vector(ctx(), r); }
+        expr_vector units() const { Z3_ast_vector r = Z3_solver_get_units(ctx(), m_solver); check_error(); return expr_vector(ctx(), r); }
         expr proof() const { Z3_ast r = Z3_solver_get_proof(ctx(), m_solver); check_error(); return expr(ctx(), r); }
         friend std::ostream & operator<<(std::ostream & out, solver const & s);
 
diff --git a/src/api/ml/z3.ml b/src/api/ml/z3.ml
index 7415fc507..231587729 100644
--- a/src/api/ml/z3.ml
+++ b/src/api/ml/z3.ml
@@ -1815,9 +1815,10 @@ struct
     | _ -> UNKNOWN
 
   let get_model x =
-  let q = Z3native.solver_get_model (gc x) x in
-    try if Z3native.is_null_model q then None else Some q with | _ -> None
-    
+    try 
+       let q = Z3native.solver_get_model (gc x) x in
+       if Z3native.is_null_model q then None else Some q 
+    with | _ -> None
 
   let get_proof x =
     let q = Z3native.solver_get_proof (gc x) x in
@@ -1953,8 +1954,10 @@ struct
     | _ -> Solver.UNKNOWN
 
   let get_model (x:optimize) =
-    let q = Z3native.optimize_get_model (gc x) x in
-    if Z3native.is_null_model q then None else Some q
+    try
+      let q = Z3native.optimize_get_model (gc x) x in
+      if Z3native.is_null_model q then None else Some q
+    with | _ -> None
 
   let get_lower (x:handle) = Z3native.optimize_get_lower (gc x.opt) x.opt x.h
   let get_upper (x:handle) = Z3native.optimize_get_upper (gc x.opt) x.opt x.h
diff --git a/src/api/python/z3/z3.py b/src/api/python/z3/z3.py
index 0eeeeaecc..8b37fb802 100644
--- a/src/api/python/z3/z3.py
+++ b/src/api/python/z3/z3.py
@@ -1258,6 +1258,11 @@ def Consts(names, sort):
         names = names.split(" ")
     return [Const(name, sort) for name in names]
 
+def FreshConst(sort, prefix='c'):
+    """Create a fresh constant of a specified sort"""
+    ctx = _get_ctx(sort.ctx)
+    return _to_expr_ref(Z3_mk_fresh_const(ctx.ref(), prefix, sort.ast), ctx)
+    
 def Var(idx, s):
     """Create a Z3 free variable. Free variables are used to create quantified formulas.
 
@@ -2176,6 +2181,8 @@ class ArithRef(ExprRef):
         >>> (x * y).sort()
         Real
         """
+        if isinstance(other, BoolRef):
+           return If(other, self, 0)
         a, b = _coerce_exprs(self, other)
         return ArithRef(_mk_bin(Z3_mk_mul, a, b), self.ctx)
 
@@ -4278,7 +4285,7 @@ def get_map_func(a):
         _z3_assert(is_map(a), "Z3 array map expression expected.")
     return FuncDeclRef(Z3_to_func_decl(a.ctx_ref(), Z3_get_decl_ast_parameter(a.ctx_ref(), a.decl().ast, 0)), a.ctx)
 
-def ArraySort(d, r):
+def ArraySort(*sig):
     """Return the Z3 array sort with the given domain and range sorts.
 
     >>> A = ArraySort(IntSort(), BoolSort())
@@ -4292,12 +4299,23 @@ def ArraySort(d, r):
     >>> AA
     Array(Int, Array(Int, Bool))
     """
+    sig = _get_args(sig)
     if __debug__:
-        _z3_assert(is_sort(d), "Z3 sort expected")
-        _z3_assert(is_sort(r), "Z3 sort expected")
-        _z3_assert(d.ctx == r.ctx, "Context mismatch")
+        _z3_assert(len(sig) > 1, "At least two arguments expected")
+    arity = len(sig) - 1
+    r = sig[arity]
+    d = sig[0]
+    if __debug__:
+        for s in sig:
+            _z3_assert(is_sort(s), "Z3 sort expected")
+            _z3_assert(s.ctx == r.ctx, "Context mismatch")
     ctx = d.ctx
-    return ArraySortRef(Z3_mk_array_sort(ctx.ref(), d.ast, r.ast), ctx)
+    if len(sig) == 2:
+        return ArraySortRef(Z3_mk_array_sort(ctx.ref(), d.ast, r.ast), ctx)
+    dom = (Sort * arity)()
+    for i in range(arity):
+        dom[i] = sig[i].ast    
+    return ArraySortRef(Z3_mk_array_sort_n(ctx.ref(), arity, dom, r.ast), ctx)
 
 def Array(name, dom, rng):
     """Return an array constant named `name` with the given domain and range sorts.
@@ -6603,7 +6621,12 @@ class Solver(Z3PPObject):
             _handle_parse_error(e, self.ctx)        
     
     def cube(self, vars = None):
-        """Get set of cubes"""
+        """Get set of cubes
+        The method takes an optional set of variables that restrict which
+        variables may be used as a starting point for cubing.
+        If vars is not None, then the first case split is based on a variable in
+        this set.
+        """
         self.cube_vs = AstVector(None, self.ctx)
         if vars is not None:
            for v in vars:
@@ -6619,19 +6642,15 @@ class Solver(Z3PPObject):
                 return
 
     def cube_vars(self):
+        """Access the set of variables that were touched by the most recently generated cube.
+        This set of variables can be used as a starting point for additional cubes.
+        The idea is that variables that appear in clauses that are reduced by the most recent
+        cube are likely more useful to cube on."""
         return self.cube_vs
 
     def proof(self):
         """Return a proof for the last `check()`. Proof construction must be enabled."""
         return _to_expr_ref(Z3_solver_get_proof(self.ctx.ref(), self.solver), self.ctx)
-
-    def from_file(self, filename):
-        """Parse assertions from a file"""
-        Z3_solver_from_file(self.ctx.ref(), self.solver, filename)
-
-    def from_string(self, s):
-        """Parse assertions from a string"""
-        Z3_solver_from_string(self.ctx.ref(), self.solver, s)
         
     def assertions(self):
         """Return an AST vector containing all added constraints.
@@ -6652,6 +6671,11 @@ class Solver(Z3PPObject):
         """
         return AstVector(Z3_solver_get_units(self.ctx.ref(), self.solver), self.ctx)
 
+    def non_units(self):
+        """Return an AST vector containing all atomic formulas in solver state that are not units.
+        """
+        return AstVector(Z3_solver_get_non_units(self.ctx.ref(), self.solver), self.ctx)
+
     def statistics(self):
         """Return statistics for the last `check()`.
 
@@ -8040,7 +8064,7 @@ def substitute(t, *m):
     """
     if isinstance(m, tuple):
         m1 = _get_args(m)
-        if isinstance(m1, list):
+        if isinstance(m1, list) and all(isinstance(p, tuple) for p in m1):
             m = m1
     if __debug__:
         _z3_assert(is_expr(t), "Z3 expression expected")
diff --git a/src/api/z3_api.h b/src/api/z3_api.h
index 2657e558d..03bce5d5e 100644
--- a/src/api/z3_api.h
+++ b/src/api/z3_api.h
@@ -6121,6 +6121,14 @@ extern "C" {
     */
     Z3_ast_vector Z3_API Z3_solver_get_units(Z3_context c, Z3_solver s);
 
+
+    /**
+       \brief Return the set of non units in the solver state.
+
+       def_API('Z3_solver_get_non_units', AST_VECTOR, (_in(CONTEXT), _in(SOLVER)))
+    */
+    Z3_ast_vector Z3_API Z3_solver_get_non_units(Z3_context c, Z3_solver s);
+
     /**
        \brief Check whether the assertions in a given solver are consistent or not.
 
diff --git a/src/ast/arith_decl_plugin.cpp b/src/ast/arith_decl_plugin.cpp
index db3604a99..cdf18875c 100644
--- a/src/ast/arith_decl_plugin.cpp
+++ b/src/ast/arith_decl_plugin.cpp
@@ -21,6 +21,7 @@ Revision History:
 #include "math/polynomial/algebraic_numbers.h"
 #include "util/id_gen.h"
 #include "ast/ast_smt2_pp.h"
+#include "util/gparams.h"
 
 struct arith_decl_plugin::algebraic_numbers_wrapper {
     unsynch_mpq_manager           m_qmanager;
@@ -487,7 +488,7 @@ func_decl * arith_decl_plugin::mk_func_decl(decl_kind k, unsigned num_parameters
         if (arity != 1 || domain[0] != m_int_decl || num_parameters != 1 || !parameters[0].is_int()) {
             m_manager->raise_exception("invalid divides application. Expects integer parameter and one argument of sort integer");
         }
-        return m_manager->mk_func_decl(symbol("divides"), 1, &m_int_decl, m_manager->mk_bool_sort(), 
+        return m_manager->mk_func_decl(symbol("divisible"), 1, &m_int_decl, m_manager->mk_bool_sort(), 
                                        func_decl_info(m_family_id, k, num_parameters, parameters));
     }
 
@@ -512,7 +513,7 @@ func_decl * arith_decl_plugin::mk_func_decl(decl_kind k, unsigned num_parameters
         if (num_args != 1 || m_manager->get_sort(args[0]) != m_int_decl || num_parameters != 1 || !parameters[0].is_int()) {
             m_manager->raise_exception("invalid divides application. Expects integer parameter and one argument of sort integer");
         }
-        return m_manager->mk_func_decl(symbol("divides"), 1, &m_int_decl, m_manager->mk_bool_sort(), 
+        return m_manager->mk_func_decl(symbol("divisible"), 1, &m_int_decl, m_manager->mk_bool_sort(), 
                                        func_decl_info(m_family_id, k, num_parameters, parameters));
     }
     if (m_manager->int_real_coercions() && use_coercion(k)) {
@@ -549,8 +550,9 @@ void arith_decl_plugin::get_op_names(svector<builtin_name>& op_names, symbol con
     op_names.push_back(builtin_name("*",OP_MUL));
     op_names.push_back(builtin_name("/",OP_DIV));
     op_names.push_back(builtin_name("div",OP_IDIV));
-    // clashes with user-defined functions
-    // op_names.push_back(builtin_name("divides",OP_IDIVIDES));
+    if (gparams::get_value("smtlib2_compliant") == "true") {
+        op_names.push_back(builtin_name("divisible",OP_IDIVIDES));
+    }
     op_names.push_back(builtin_name("rem",OP_REM));
     op_names.push_back(builtin_name("mod",OP_MOD));
     op_names.push_back(builtin_name("to_real",OP_TO_REAL));
diff --git a/src/ast/ast.cpp b/src/ast/ast.cpp
index 11a15492c..65d3c7821 100644
--- a/src/ast/ast.cpp
+++ b/src/ast/ast.cpp
@@ -1656,6 +1656,12 @@ bool ast_manager::are_distinct(expr* a, expr* b) const {
     return false;
 }
 
+func_decl* ast_manager::get_rec_fun_decl(quantifier* q) const {
+    SASSERT(is_rec_fun_def(q)); 
+    return to_app(to_app(q->get_pattern(0))->get_arg(0))->get_decl(); 
+}
+
+
 void ast_manager::register_plugin(family_id id, decl_plugin * plugin) {
     SASSERT(m_plugins.get(id, 0) == 0);
     m_plugins.setx(id, plugin, 0);
diff --git a/src/ast/ast.h b/src/ast/ast.h
index 89df04961..c1193dfbd 100644
--- a/src/ast/ast.h
+++ b/src/ast/ast.h
@@ -1632,6 +1632,7 @@ public:
 
     bool is_rec_fun_def(quantifier* q) const { return q->get_qid() == m_rec_fun; }
     bool is_lambda_def(quantifier* q) const { return q->get_qid() == m_lambda_def; }
+    func_decl* get_rec_fun_decl(quantifier* q) const;
     
     symbol const& rec_fun_qid() const { return m_rec_fun; }
 
diff --git a/src/ast/seq_decl_plugin.cpp b/src/ast/seq_decl_plugin.cpp
index 07c6fd06a..06f30cbdf 100644
--- a/src/ast/seq_decl_plugin.cpp
+++ b/src/ast/seq_decl_plugin.cpp
@@ -842,7 +842,9 @@ void seq_decl_plugin::get_sort_names(svector<builtin_name> & sort_names, symbol
 }
 
 app* seq_decl_plugin::mk_string(symbol const& s) {
-    parameter param(s);
+    zstring canonStr(s.bare_str());
+    symbol canonSym(canonStr.encode().c_str());
+    parameter param(canonSym);
     func_decl* f = m_manager->mk_const_decl(m_stringc_sym, m_string,
                                             func_decl_info(m_family_id, OP_STRING_CONST, 1, &param));
     return m_manager->mk_const(f);
diff --git a/src/cmd_context/basic_cmds.cpp b/src/cmd_context/basic_cmds.cpp
index e65eb1b32..ea5994ece 100644
--- a/src/cmd_context/basic_cmds.cpp
+++ b/src/cmd_context/basic_cmds.cpp
@@ -17,7 +17,7 @@ Notes:
 --*/
 #include "util/gparams.h"
 #include "util/env_params.h"
-#include "util/version.h"
+#include "util/z3_version.h"
 #include "ast/ast_smt_pp.h"
 #include "ast/ast_smt2_pp.h"
 #include "ast/ast_pp_dot.h"
diff --git a/src/muz/base/dl_rule_set.cpp b/src/muz/base/dl_rule_set.cpp
index 5114fdca4..80af80266 100644
--- a/src/muz/base/dl_rule_set.cpp
+++ b/src/muz/base/dl_rule_set.cpp
@@ -31,7 +31,7 @@ namespace datalog {
     rule_dependencies::rule_dependencies(const rule_dependencies & o, bool reversed):
         m_context(o.m_context) {
         if (reversed) {
-            for (auto & kv : o) { 
+            for (auto & kv : o) {
                 func_decl * pred = kv.m_key;
                 item_set & orig_items = *kv.get_value();
 
@@ -114,8 +114,8 @@ namespace datalog {
                 app* a = to_app(e);
                 d = a->get_decl();
                 if (m_context.is_predicate(d)) {
-                    // insert d and ensure the invariant 
-                    // that every predicate is present as 
+                    // insert d and ensure the invariant
+                    // that every predicate is present as
                     // a key in m_data
                     s.insert(d);
                     ensure_key(d);
@@ -148,7 +148,7 @@ namespace datalog {
             item_set& itms = *kv.get_value();
             set_intersection(itms, allowed);
         }
-        for (func_decl* f : to_remove) 
+        for (func_decl* f : to_remove)
             remove_m_data_entry(f);
     }
 
@@ -253,18 +253,18 @@ namespace datalog {
     //
     // -----------------------------------
 
-    rule_set::rule_set(context & ctx) 
-          : m_context(ctx), 
-            m_rule_manager(ctx.get_rule_manager()), 
-            m_rules(m_rule_manager), 
+    rule_set::rule_set(context & ctx)
+          : m_context(ctx),
+            m_rule_manager(ctx.get_rule_manager()),
+            m_rules(m_rule_manager),
             m_deps(ctx),
             m_stratifier(nullptr),
             m_refs(ctx.get_manager()) {
     }
 
-    rule_set::rule_set(const rule_set & other) 
-        : m_context(other.m_context), 
-          m_rule_manager(other.m_rule_manager), 
+    rule_set::rule_set(const rule_set & other)
+        : m_context(other.m_context),
+          m_rule_manager(other.m_rule_manager),
           m_rules(m_rule_manager),
           m_deps(other.m_context),
           m_stratifier(nullptr),
@@ -353,10 +353,27 @@ namespace datalog {
                 break;                                  \
             }                                           \
         }                                               \
-        
+
         DEL_VECTOR(*rules);
         DEL_VECTOR(m_rules);
-    }    
+    }
+
+    void rule_set::replace_rule(rule * r, rule * other) {
+        TRACE("dl", r->display(m_context, tout << "replace:"););
+        func_decl* d = r->get_decl();
+        rule_vector* rules = m_head2rules.find(d);
+#define REPLACE_VECTOR(_v)                              \
+        for (unsigned i = (_v).size(); i > 0; ) {       \
+            --i;                                        \
+            if ((_v)[i] == r) {                         \
+                (_v)[i] = other;                        \
+                break;                                  \
+            }                                           \
+        }                                               \
+
+        REPLACE_VECTOR(*rules);
+        REPLACE_VECTOR(m_rules);
+    }
 
     void rule_set::ensure_closed() {
         if (!is_closed()) {
@@ -365,7 +382,7 @@ namespace datalog {
     }
 
     bool rule_set::close() {
-        SASSERT(!is_closed()); //the rule_set is not already closed        
+        SASSERT(!is_closed()); //the rule_set is not already closed
         m_deps.populate(*this);
         m_stratifier = alloc(rule_stratifier, m_deps);
         if (!stratified_negation()) {
@@ -426,7 +443,7 @@ namespace datalog {
         inherit_predicates(src);
     }
 
-    const rule_vector & rule_set::get_predicate_rules(func_decl * pred) const { 
+    const rule_vector & rule_set::get_predicate_rules(func_decl * pred) const {
         decl2rules::obj_map_entry * e = m_head2rules.find_core(pred);
         if (!e) {
             return m_empty_rule_vector;
@@ -519,7 +536,7 @@ namespace datalog {
                 out << "\n";
                 non_empty = false;
             }
-            
+
             for (func_decl * first : *strat) {
                 const func_decl_set & deps = m_deps.get_deps(first);
                 for (func_decl * dep : deps) {
@@ -545,8 +562,8 @@ namespace datalog {
     unsigned rule_stratifier::get_predicate_strat(func_decl * pred) const {
         unsigned num;
         if (!m_pred_strat_nums.find(pred, num)) {
-            //the number of the predicate is not stored, therefore it did not appear 
-            //in the algorithm and therefore it does not depend on anything and nothing 
+            //the number of the predicate is not stored, therefore it did not appear
+            //in the algorithm and therefore it does not depend on anything and nothing
             //depends on it. So it is safe to assign zero strate to it, although it is
             //not strictly true.
             num = 0;
@@ -641,7 +658,7 @@ namespace datalog {
         }
 
 
-        // We put components whose indegree is zero to m_strats and assign its 
+        // We put components whose indegree is zero to m_strats and assign its
         // m_components entry to zero.
         unsigned comp_cnt = m_components.size();
         for (unsigned i = 0; i < comp_cnt; i++) {
@@ -680,7 +697,7 @@ namespace datalog {
             strats_index++;
         }
         //we have managed to topologicaly order all the components
-        SASSERT(std::find_if(m_components.begin(), m_components.end(), 
+        SASSERT(std::find_if(m_components.begin(), m_components.end(),
             std::bind1st(std::not_equal_to<item_set*>(), (item_set*)0)) == m_components.end());
 
         //reverse the strats array, so that the only the later components would depend on earlier ones
@@ -713,7 +730,7 @@ namespace datalog {
             }
             out << "\n";
         }
-        
+
     }
 
 };
diff --git a/src/muz/base/dl_rule_set.h b/src/muz/base/dl_rule_set.h
index 736dd8888..e870e369c 100644
--- a/src/muz/base/dl_rule_set.h
+++ b/src/muz/base/dl_rule_set.h
@@ -77,7 +77,7 @@ namespace datalog {
            \brief Number of predicates that depend on \c f.
          */
         unsigned out_degree(func_decl * f) const;
-        
+
         /**
            \brief If the rependency graph is acyclic, put all elements into \c res
              ordered so that elements can depend only on elements that are before them.
@@ -131,7 +131,7 @@ namespace datalog {
             it must exist for the whole lifetime of the \c stratifier object.
         */
         rule_stratifier(const rule_dependencies & deps)
-            : m_deps(deps), m_next_preorder(0) 
+            : m_deps(deps), m_next_preorder(0)
         {
             process();
         }
@@ -145,7 +145,7 @@ namespace datalog {
         const comp_vector & get_strats() const { return m_strats; }
 
         unsigned get_predicate_strat(func_decl * pred) const;
-        
+
         void display( std::ostream & out ) const;
     };
 
@@ -203,6 +203,10 @@ namespace datalog {
            \brief Remove rule \c r from the rule set.
         */
         void del_rule(rule * r);
+        /**
+           \brief Replace a rule \c r with the rule \c other
+        */
+        void replace_rule(rule * r, rule * other);
 
         /**
            \brief Add all rules from a different rule_set.
@@ -276,8 +280,7 @@ namespace datalog {
     inline std::ostream& operator<<(std::ostream& out, rule_set const& r) { r.display(out); return out; }
 
 
-    
+
 };
 
 #endif /* DL_RULE_SET_H_ */
-
diff --git a/src/muz/base/fp_params.pyg b/src/muz/base/fp_params.pyg
index 8cee99540..18eb85662 100644
--- a/src/muz/base/fp_params.pyg
+++ b/src/muz/base/fp_params.pyg
@@ -177,5 +177,6 @@ def_module_params('fp',
                           ('spacer.dump_threshold', DOUBLE, 5.0, 'Threshold in seconds on dumping benchmarks'),
                           ('spacer.gpdr', BOOL, False, 'Use GPDR solving strategy for non-linear CHC'),
                           ('spacer.gpdr.bfs', BOOL, True, 'Use BFS exploration strategy for expanding model search'),
+                          ('spacer.use_bg_invs', BOOL, False, 'Enable external background invariants'),
 
                           ))
diff --git a/src/muz/spacer/spacer_context.cpp b/src/muz/spacer/spacer_context.cpp
index f0e86f1a5..dd6656954 100644
--- a/src/muz/spacer/spacer_context.cpp
+++ b/src/muz/spacer/spacer_context.cpp
@@ -493,7 +493,8 @@ lemma::lemma (ast_manager &manager, expr * body, unsigned lvl) :
     m_pob(nullptr), m_ctp(nullptr),
     m_lvl(lvl), m_init_lvl(m_lvl),
     m_bumped(0), m_weakness(WEAKNESS_MAX),
-    m_external(false), m_blocked(false) {
+    m_external(false), m_blocked(false),
+    m_background(false) {
     SASSERT(m_body);
     normalize(m_body, m_body);
 }
@@ -505,7 +506,8 @@ lemma::lemma(pob_ref const &p) :
     m_pob(p), m_ctp(nullptr),
     m_lvl(p->level()), m_init_lvl(m_lvl),
     m_bumped(0), m_weakness(p->weakness()),
-    m_external(false), m_blocked(false) {
+    m_external(false), m_blocked(false),
+    m_background(false) {
     SASSERT(m_pob);
     m_pob->get_skolems(m_zks);
     add_binding(m_pob->get_binding());
@@ -519,8 +521,8 @@ lemma::lemma(pob_ref const &p, expr_ref_vector &cube, unsigned lvl) :
     m_pob(p), m_ctp(nullptr),
     m_lvl(p->level()), m_init_lvl(m_lvl),
     m_bumped(0), m_weakness(p->weakness()),
-    m_external(false), m_blocked(false)
-{
+    m_external(false), m_blocked(false),
+    m_background(false) {
     if (m_pob) {
         m_pob->get_skolems(m_zks);
         add_binding(m_pob->get_binding());
@@ -921,10 +923,10 @@ void pred_transformer::simplify_formulas()
 {m_frames.simplify_formulas ();}
 
 
-expr_ref pred_transformer::get_formulas(unsigned level) const
+expr_ref pred_transformer::get_formulas(unsigned level, bool bg) const
 {
     expr_ref_vector res(m);
-    m_frames.get_frame_geq_lemmas (level, res);
+    m_frames.get_frame_geq_lemmas (level, res, bg);
     return mk_and(res);
 }
 
@@ -935,6 +937,7 @@ bool pred_transformer::propagate_to_next_level (unsigned src_level)
 /// \brief adds a lemma to the solver and to child solvers
 void pred_transformer::add_lemma_core(lemma* lemma, bool ground_only)
 {
+    SASSERT(!lemma->is_background());
     unsigned lvl = lemma->level();
     expr* l = lemma->get_expr();
     SASSERT(!lemma->is_ground() || is_clause(m, l));
@@ -975,8 +978,9 @@ void pred_transformer::add_lemma_core(lemma* lemma, bool ground_only)
                                       next_level(lvl), ground_only); }
 }
 
-bool pred_transformer::add_lemma (expr *e, unsigned lvl) {
+bool pred_transformer::add_lemma (expr *e, unsigned lvl, bool bg) {
     lemma_ref lem = alloc(lemma, m, e, lvl);
+    lem->set_background(bg);
     return m_frames.add_lemma(lem.get());
 }
 
@@ -1217,15 +1221,18 @@ expr_ref pred_transformer::get_origin_summary (model &mdl,
 }
 
 
-void pred_transformer::add_cover(unsigned level, expr* property)
+void pred_transformer::add_cover(unsigned level, expr* property, bool bg)
 {
+    SASSERT(!bg || is_infty_level(level));
     // replace bound variables by local constants.
     expr_ref result(property, m), v(m), c(m);
     expr_substitution sub(m);
+    proof_ref pr(m);
+    pr = m.mk_asserted(m.mk_true());
     for (unsigned i = 0; i < sig_size(); ++i) {
         c = m.mk_const(pm.o2n(sig(i), 0));
         v = m.mk_var(i, sig(i)->get_range());
-        sub.insert(v, c);
+        sub.insert(v, c, pr);
     }
     scoped_ptr<expr_replacer> rep = mk_default_expr_replacer(m);
     rep->set_substitution(&sub);
@@ -1236,13 +1243,38 @@ void pred_transformer::add_cover(unsigned level, expr* property)
     expr_ref_vector lemmas(m);
     flatten_and(result, lemmas);
     for (unsigned i = 0, sz = lemmas.size(); i < sz; ++i) {
-        add_lemma(lemmas.get(i), level);
+        add_lemma(lemmas.get(i), level, bg);
     }
 }
 
 void pred_transformer::propagate_to_infinity (unsigned level)
 {m_frames.propagate_to_infinity (level);}
 
+// compute a conjunction of all background facts
+void pred_transformer::get_pred_bg_invs(expr_ref_vector& out) {
+    expr_ref inv(m), tmp1(m), tmp2(m);
+    ptr_vector<func_decl> preds;
+    for (auto kv : m_pt_rules) {
+        expr* tag = kv.m_value->tag();
+        datalog::rule const &r = kv.m_value->rule();
+        find_predecessors (r, preds);
+
+        for (unsigned i = 0, preds_sz = preds.size(); i < preds_sz; i++) {
+            func_decl* pre = preds[i];
+            pred_transformer &pt = ctx.get_pred_transformer(pre);
+            const lemma_ref_vector &invs = pt.get_bg_invs();
+            CTRACE("spacer", !invs.empty(),
+                   tout << "add-bg-invariant: " << mk_pp (pre, m) << "\n";);
+            for (auto inv : invs) {
+                // tag -> inv1 ...  tag -> invn
+                tmp1 = m.mk_implies(tag, inv->get_expr());
+                pm.formula_n2o(tmp1, tmp2, i);
+                out.push_back(tmp2);
+                TRACE("spacer", tout << tmp2 << "\n";);
+            }
+        }
+    }
+}
 
 
 /// \brief Returns true if the obligation is already blocked by current lemmas
@@ -1344,6 +1376,14 @@ lbool pred_transformer::is_reachable(pob& n, expr_ref_vector* core,
 
     expr_ref_vector post (m), reach_assumps (m);
     post.push_back (n.post ());
+    flatten_and(post);
+
+    // if equality propagation is disabled in arithmetic, expand
+    // equality literals into two inequalities to increase the space
+    // for interpolation
+    if (!ctx.use_eq_prop()) {
+        expand_literals(m, post);
+    }
 
     // populate reach_assumps
     if (n.level () > 0 && !m_all_init) {
@@ -1472,7 +1512,7 @@ bool pred_transformer::is_invariant(unsigned level, lemma* lem,
     expr_ref lemma_expr(m);
     lemma_expr = lem->get_expr();
 
-    expr_ref_vector conj(m), aux(m);
+    expr_ref_vector cand(m), aux(m), conj(m);
     expr_ref gnd_lemma(m);
 
 
@@ -1482,8 +1522,8 @@ bool pred_transformer::is_invariant(unsigned level, lemma* lem,
         lemma_expr = gnd_lemma.get();
     }
 
-    conj.push_back(mk_not(m, lemma_expr));
-    flatten_and (conj);
+    cand.push_back(mk_not(m, lemma_expr));
+    flatten_and (cand);
 
     prop_solver::scoped_level _sl(*m_solver, level);
     prop_solver::scoped_subset_core _sc (*m_solver, true);
@@ -1494,9 +1534,12 @@ bool pred_transformer::is_invariant(unsigned level, lemma* lem,
     if (ctx.use_ctp()) {mdl_ref_ptr = &mdl;}
     m_solver->set_core(core);
     m_solver->set_model(mdl_ref_ptr);
-    expr * bg = m_extend_lit.get ();
-    lbool r = m_solver->check_assumptions (conj, aux, m_transition_clause,
-                                          1, &bg, 1);
+
+    conj.push_back(m_extend_lit);
+    if (ctx.use_bg_invs()) get_pred_bg_invs(conj);
+
+    lbool r = m_solver->check_assumptions (cand, aux, m_transition_clause,
+                                           conj.size(), conj.c_ptr(), 1);
     if (r == l_false) {
         solver_level = m_solver->uses_level ();
         lem->reset_ctp();
@@ -1527,6 +1570,7 @@ bool pred_transformer::check_inductive(unsigned level, expr_ref_vector& state,
     m_solver->set_core(&core);
     m_solver->set_model (nullptr);
     expr_ref_vector aux (m);
+    if (ctx.use_bg_invs()) get_pred_bg_invs(conj);
     conj.push_back (m_extend_lit);
     lbool res = m_solver->check_assumptions (state, aux,
                                             m_transition_clause,
@@ -1941,14 +1985,27 @@ void pred_transformer::update_solver_with_rfs(prop_solver *solver,
 }
 
 /// pred_transformer::frames
-
-
 bool pred_transformer::frames::add_lemma(lemma *new_lemma)
 {
     TRACE("spacer", tout << "add-lemma: " << pp_level(new_lemma->level()) << " "
           << m_pt.head()->get_name() << " "
           << mk_pp(new_lemma->get_expr(), m_pt.get_ast_manager()) << "\n";);
 
+    if (new_lemma->is_background()) {
+        SASSERT (is_infty_level(new_lemma->level()));
+
+        for (auto &l : m_bg_invs) {
+            if (l->get_expr() == new_lemma->get_expr()) return false;
+        }
+        TRACE("spacer", tout << "add-external-lemma: "
+              << pp_level(new_lemma->level()) << " "
+              << m_pt.head()->get_name() << " "
+              << mk_pp(new_lemma->get_expr(), m_pt.get_ast_manager()) << "\n";);
+
+        m_bg_invs.push_back(new_lemma);
+        return true;
+    }
+
     unsigned i = 0;
     for (auto *old_lemma : m_lemmas) {
         if (old_lemma->get_expr() == new_lemma->get_expr()) {
@@ -2295,6 +2352,7 @@ void context::updt_params() {
     m_use_restarts = m_params.spacer_restarts();
     m_restart_initial_threshold = m_params.spacer_restart_initial_threshold();
     m_pdr_bfs = m_params.spacer_gpdr_bfs();
+    m_use_bg_invs = m_params.spacer_use_bg_invs();
 
     if (m_use_gpdr) {
         // set options to be compatible with GPDR
@@ -2423,36 +2481,38 @@ expr_ref context::get_cover_delta(int level, func_decl* p_orig, func_decl* p)
     if (m_rels.find(p, pt)) {
         return pt->get_cover_delta(p_orig, level);
     } else {
-        IF_VERBOSE(10, verbose_stream() << "did not find predicate " << p->get_name() << "\n";);
+        IF_VERBOSE(10, verbose_stream() << "did not find predicate "
+                   << p->get_name() << "\n";);
         return expr_ref(m.mk_true(), m);
     }
 }
 
-void context::add_cover(int level, func_decl* p, expr* property)
+void context::add_cover(int level, func_decl* p, expr* property, bool bg)
 {
+    scoped_proof _pf_(m);
+
     pred_transformer* pt = nullptr;
     if (!m_rels.find(p, pt)) {
         pt = alloc(pred_transformer, *this, get_manager(), p);
         m_rels.insert(p, pt);
-        IF_VERBOSE(10, verbose_stream() << "did not find predicate " << p->get_name() << "\n";);
+        IF_VERBOSE(10, verbose_stream() << "did not find predicate "
+                   << p->get_name() << "\n";);
     }
     unsigned lvl = (level == -1)?infty_level():((unsigned)level);
-    pt->add_cover(lvl, property);
+    pt->add_cover(lvl, property, bg);
 }
 
 void context::add_invariant (func_decl *p, expr *property)
-{add_cover (infty_level(), p, property);}
+{add_cover (infty_level(), p, property, true);}
 
-expr_ref context::get_reachable(func_decl *p)
-{
+expr_ref context::get_reachable(func_decl *p) {
     pred_transformer* pt = nullptr;
     if (!m_rels.find(p, pt))
     { return expr_ref(m.mk_false(), m); }
     return pt->get_reachable();
 }
 
-bool context::validate()
-{
+bool context::validate() {
     if (!m_validate_result) { return true; }
 
     std::stringstream msg;
@@ -2483,7 +2543,7 @@ bool context::validate()
         model_ref model;
         vector<relation_info> rs;
         model_converter_ref mc;
-        get_level_property(m_inductive_lvl, refs, rs);
+        get_level_property(m_inductive_lvl, refs, rs, use_bg_invs());
         inductive_property ex(m, mc, rs);
         ex.to_model(model);
         var_subst vs(m, false);
@@ -2624,13 +2684,13 @@ void context::init_lemma_generalizers()
 }
 
 void context::get_level_property(unsigned lvl, expr_ref_vector& res,
-                                 vector<relation_info>& rs) const {
+                                 vector<relation_info>& rs, bool with_bg) const {
     for (auto const& kv : m_rels) {
         pred_transformer* r = kv.m_value;
         if (r->head() == m_query_pred) {
             continue;
         }
-        expr_ref conj = r->get_formulas(lvl);
+        expr_ref conj = r->get_formulas(lvl, with_bg);
         m_pm.formula_n2o(0, false, conj);
         res.push_back(conj);
         ptr_vector<func_decl> sig(r->head()->get_arity(), r->sig());
@@ -2662,7 +2722,7 @@ lbool context::solve(unsigned from_lvl)
             IF_VERBOSE(1, {
                     expr_ref_vector refs(m);
                     vector<relation_info> rs;
-                    get_level_property(m_inductive_lvl, refs, rs);
+                    get_level_property(m_inductive_lvl, refs, rs, use_bg_invs());
                     model_converter_ref mc;
                     inductive_property ex(m, mc, rs);
                     verbose_stream() << ex.to_string();
@@ -2844,7 +2904,7 @@ model_ref context::get_model()
     model_ref model;
     expr_ref_vector refs(m);
     vector<relation_info> rs;
-    get_level_property(m_inductive_lvl, refs, rs);
+    get_level_property(m_inductive_lvl, refs, rs, use_bg_invs());
     inductive_property ex(m, const_cast<model_converter_ref&>(m_mc), rs);
     ex.to_model (model);
     return model;
@@ -2877,7 +2937,7 @@ expr_ref context::mk_unsat_answer() const
 {
     expr_ref_vector refs(m);
     vector<relation_info> rs;
-    get_level_property(m_inductive_lvl, refs, rs);
+    get_level_property(m_inductive_lvl, refs, rs, use_bg_invs());
     inductive_property ex(m, const_cast<model_converter_ref&>(m_mc), rs);
     return ex.to_expr();
 }
diff --git a/src/muz/spacer/spacer_context.h b/src/muz/spacer/spacer_context.h
index 48c27f96d..0d8b2daf6 100644
--- a/src/muz/spacer/spacer_context.h
+++ b/src/muz/spacer/spacer_context.h
@@ -128,8 +128,9 @@ class lemma {
     unsigned m_init_lvl;       // level at which lemma was created
     unsigned m_bumped:16;
     unsigned m_weakness:16;
-    unsigned m_external:1;
-    unsigned m_blocked:1;
+    unsigned m_external:1;    // external lemma from another solver
+    unsigned m_blocked:1;     // blocked by CTP
+    unsigned m_background:1;  // background assumed fact
 
     void mk_expr_core();
     void mk_cube_core();
@@ -163,6 +164,9 @@ public:
     void set_external(bool ext){m_external = ext;}
     bool external() { return m_external;}
 
+    void set_background(bool v) {m_background = v;}
+    bool is_background() {return m_background;}
+
     bool is_blocked() {return m_blocked;}
     void set_blocked(bool v) {m_blocked=v;}
 
@@ -222,6 +226,7 @@ class pred_transformer {
         pred_transformer &m_pt;            // parent pred_transformer
         lemma_ref_vector m_pinned_lemmas;  // all created lemmas
         lemma_ref_vector m_lemmas;         // active lemmas
+        lemma_ref_vector m_bg_invs;        // background (assumed) invariants
         unsigned m_size;                   // num of frames
 
         bool m_sorted;                     // true if m_lemmas is sorted by m_lt
@@ -230,7 +235,8 @@ class pred_transformer {
         void sort ();
 
     public:
-        frames (pred_transformer &pt) : m_pt (pt), m_size(0), m_sorted (true) {}
+        frames (pred_transformer &pt) : m_pt (pt),
+                                        m_size(0), m_sorted (true) {}
         ~frames() {}
         void simplify_formulas ();
 
@@ -245,17 +251,25 @@ class pred_transformer {
                 }
             }
         }
-        void get_frame_geq_lemmas (unsigned level, expr_ref_vector &out) const {
+        void get_frame_geq_lemmas (unsigned level, expr_ref_vector &out,
+                                   bool with_bg = false) const {
             for (auto &lemma : m_lemmas) {
                 if (lemma->level() >= level) {
                     out.push_back(lemma->get_expr());
                 }
             }
+            if (with_bg) {
+                for (auto &lemma : m_bg_invs)
+                    out.push_back(lemma->get_expr());
+            }
         }
 
-        unsigned size () const {return m_size;}
-        unsigned lemma_size () const {return m_lemmas.size ();}
-        void add_frame () {m_size++;}
+        const lemma_ref_vector& get_bg_invs() const {return m_bg_invs;}
+        unsigned size() const {return m_size;}
+        unsigned lemma_size() const {return m_lemmas.size ();}
+        unsigned bg_invs_size() const {return m_bg_invs.size();}
+
+        void add_frame() {m_size++;}
         void inherit_frames (frames &other) {
             for (auto &other_lemma : other.m_lemmas) {
                 lemma_ref new_lemma = alloc(lemma, m_pt.get_ast_manager(),
@@ -265,6 +279,7 @@ class pred_transformer {
                 add_lemma(new_lemma.get());
             }
             m_sorted = false;
+            m_bg_invs.append(other.m_bg_invs);
         }
 
         bool add_lemma (lemma *new_lemma);
@@ -418,6 +433,11 @@ class pred_transformer {
 
     app_ref mk_fresh_rf_tag ();
 
+    // get tagged formulae of all of the background invariants for all of the
+    // predecessors of the current transformer
+    void get_pred_bg_invs(expr_ref_vector &out);
+    const lemma_ref_vector &get_bg_invs() const {return m_frames.get_bg_invs();}
+
 public:
     pred_transformer(context& ctx, manager& pm, func_decl* head);
     ~pred_transformer() {}
@@ -448,7 +468,7 @@ public:
     }
     unsigned get_num_levels() const {return m_frames.size ();}
     expr_ref get_cover_delta(func_decl* p_orig, int level);
-    void     add_cover(unsigned level, expr* property);
+    void     add_cover(unsigned level, expr* property, bool bg = false);
     expr_ref get_reachable();
 
     std::ostream& display(std::ostream& strm) const;
@@ -484,7 +504,7 @@ public:
     bool propagate_to_next_level(unsigned level);
     void propagate_to_infinity(unsigned level);
     /// \brief  Add a lemma to the current context and all users
-    bool add_lemma(expr * lemma, unsigned lvl);
+    bool add_lemma(expr * e, unsigned lvl, bool bg);
     bool add_lemma(lemma* lem) {return m_frames.add_lemma(lem);}
     expr* get_reach_case_var (unsigned idx) const;
     bool has_rfs () const { return !m_reach_facts.empty () ;}
@@ -527,7 +547,7 @@ public:
     bool check_inductive(unsigned level, expr_ref_vector& state,
                          unsigned& assumes_level, unsigned weakness = UINT_MAX);
 
-    expr_ref get_formulas(unsigned level) const;
+    expr_ref get_formulas(unsigned level, bool bg = false) const;
 
     void simplify_formulas();
 
@@ -958,6 +978,7 @@ class context {
     bool                 m_simplify_formulas_pre;
     bool                 m_simplify_formulas_post;
     bool                 m_pdr_bfs;
+    bool                 m_use_bg_invs;
     unsigned             m_push_pob_max_depth;
     unsigned             m_max_level;
     unsigned             m_restart_initial_threshold;
@@ -992,7 +1013,8 @@ class context {
 
     // Generate inductive property
     void get_level_property(unsigned lvl, expr_ref_vector& res,
-                            vector<relation_info> & rs) const;
+                            vector<relation_info> & rs,
+                            bool with_bg = false) const;
 
 
     // Initialization
@@ -1027,18 +1049,20 @@ public:
 
 
     const fp_params &get_params() const { return m_params; }
-    bool use_native_mbp () {return m_use_native_mbp;}
-    bool use_ground_pob () {return m_ground_pob;}
-    bool use_instantiate () {return m_instantiate;}
-    bool weak_abs() {return m_weak_abs;}
-    bool use_qlemmas () {return m_use_qlemmas;}
-    bool use_euf_gen() {return m_use_euf_gen;}
-    bool simplify_pob() {return m_simplify_pob;}
-    bool use_ctp() {return m_use_ctp;}
-    bool use_inc_clause() {return m_use_inc_clause;}
-    unsigned blast_term_ite_inflation() {return m_blast_term_ite_inflation;}
-    bool elim_aux() {return m_elim_aux;}
-    bool reach_dnf() {return m_reach_dnf;}
+    bool use_eq_prop() const {return m_use_eq_prop;}
+    bool use_native_mbp() const {return m_use_native_mbp;}
+    bool use_ground_pob() const {return m_ground_pob;}
+    bool use_instantiate() const {return m_instantiate;}
+    bool weak_abs() const {return m_weak_abs;}
+    bool use_qlemmas() const {return m_use_qlemmas;}
+    bool use_euf_gen() const {return m_use_euf_gen;}
+    bool simplify_pob() const {return m_simplify_pob;}
+    bool use_ctp() const {return m_use_ctp;}
+    bool use_inc_clause() const {return m_use_inc_clause;}
+    unsigned blast_term_ite_inflation() const {return m_blast_term_ite_inflation;}
+    bool elim_aux() const {return m_elim_aux;}
+    bool reach_dnf() const {return m_reach_dnf;}
+    bool use_bg_invs() const {return m_use_bg_invs;}
 
     ast_manager&      get_ast_manager() const {return m;}
     manager&          get_manager() {return m_pm;}
@@ -1081,7 +1105,7 @@ public:
     unsigned get_num_levels(func_decl* p);
 
     expr_ref get_cover_delta(int level, func_decl* p_orig, func_decl* p);
-    void add_cover(int level, func_decl* pred, expr* property);
+    void add_cover(int level, func_decl* pred, expr* property, bool bg = false);
     expr_ref get_reachable (func_decl* p);
     void add_invariant (func_decl *pred, expr* property);
     model_ref get_model();
diff --git a/src/muz/spacer/spacer_prop_solver.cpp b/src/muz/spacer/spacer_prop_solver.cpp
index 004ea6754..4899ad692 100644
--- a/src/muz/spacer/spacer_prop_solver.cpp
+++ b/src/muz/spacer/spacer_prop_solver.cpp
@@ -94,6 +94,7 @@ void prop_solver::add_level()
 
 void prop_solver::ensure_level(unsigned lvl)
 {
+    if (is_infty_level(lvl)) return;
     while (lvl >= level_cnt()) {
         add_level();
     }
diff --git a/src/muz/spacer/spacer_util.h b/src/muz/spacer/spacer_util.h
index 8da574d0e..0328d8640 100644
--- a/src/muz/spacer/spacer_util.h
+++ b/src/muz/spacer/spacer_util.h
@@ -48,7 +48,8 @@ namespace spacer {
     }
 
     inline bool is_infty_level(unsigned lvl) {
-        return lvl == infty_level ();
+        // XXX: level is 16 bits in class pob
+        return lvl >= 65535;
     }
 
     inline unsigned next_level(unsigned lvl) {
diff --git a/src/muz/transforms/CMakeLists.txt b/src/muz/transforms/CMakeLists.txt
index e92aa1c2f..62272450c 100644
--- a/src/muz/transforms/CMakeLists.txt
+++ b/src/muz/transforms/CMakeLists.txt
@@ -25,6 +25,7 @@ z3_add_component(transforms
     dl_mk_array_eq_rewrite.cpp
     dl_mk_array_instantiation.cpp
     dl_mk_elim_term_ite.cpp
+    dl_mk_synchronize.cpp
   COMPONENT_DEPENDENCIES
     dataflow
     hilbert
diff --git a/src/muz/transforms/dl_mk_synchronize.cpp b/src/muz/transforms/dl_mk_synchronize.cpp
new file mode 100644
index 000000000..348c9b5de
--- /dev/null
+++ b/src/muz/transforms/dl_mk_synchronize.cpp
@@ -0,0 +1,376 @@
+/*++
+Copyright (c) 2017-2018 Saint-Petersburg State University
+
+Module Name:
+
+    dl_mk_synchronize.h
+
+Abstract:
+
+    Rule transformer that attempts to merge recursive iterations
+    relaxing the shape of the inductive invariant.
+
+Author:
+
+    Dmitry Mordvinov (dvvrd) 2017-05-24
+    Lidiia Chernigovskaia (LChernigovskaya) 2017-10-20
+
+Revision History:
+
+--*/
+#include "muz/transforms/dl_mk_synchronize.h"
+#include <algorithm>
+
+namespace datalog {
+
+    typedef mk_synchronize::item_set_vector item_set_vector;
+
+    mk_synchronize::mk_synchronize(context& ctx, unsigned priority):
+        rule_transformer::plugin(priority, false),
+        m_ctx(ctx),
+        m(ctx.get_manager()),
+        rm(ctx.get_rule_manager())
+    {}
+
+    bool mk_synchronize::is_recursive(rule &r, func_decl &decl) const {
+        func_decl *hdecl = r.get_head()->get_decl();
+        // AG: shouldn't decl appear in the body?
+        if (hdecl == &decl)  return true;
+        auto & strata = m_stratifier->get_strats();
+        unsigned num_of_stratum = m_stratifier->get_predicate_strat(hdecl);
+        return strata[num_of_stratum]->contains(&decl);
+    }
+
+    bool mk_synchronize::has_recursive_premise(app * app) const {
+        func_decl* app_decl = app->get_decl();
+        if (m_deps->get_deps(app_decl).contains(app_decl)) {
+            return true;
+        }
+        rule_stratifier::comp_vector const & strata = m_stratifier->get_strats();
+        unsigned num_of_stratum = m_stratifier->get_predicate_strat(app_decl);
+        return strata[num_of_stratum]->size() > 1;
+    }
+
+    item_set_vector mk_synchronize::add_merged_decls(ptr_vector<app> & apps) {
+        unsigned sz = apps.size();
+        item_set_vector merged_decls;
+        merged_decls.resize(sz);
+        auto & strata = m_stratifier->get_strats();
+        for (unsigned j = 0; j < sz; ++j) {
+            unsigned nos;
+            nos = m_stratifier->get_predicate_strat(apps[j]->get_decl());
+            merged_decls[j] = strata[nos];
+        }
+        return merged_decls;
+    }
+
+    void mk_synchronize::add_new_rel_symbols(unsigned idx,
+                                             item_set_vector const & decls,
+                                             ptr_vector<func_decl> & decls_buf,
+                                             bool & was_added) {
+        if (idx >= decls.size()) {
+            string_buffer<> buffer;
+            ptr_vector<sort> domain;
+            for (auto &d : decls_buf) {
+                buffer << d->get_name() << "!!";
+                domain.append(d->get_arity(), d->get_domain());
+            }
+
+            symbol new_name = symbol(buffer.c_str());
+
+            if (!m_cache.contains(new_name)) {
+                was_added = true;
+                func_decl* orig = decls_buf[0];
+                func_decl* product_pred = m_ctx.mk_fresh_head_predicate(new_name,
+                    symbol::null, domain.size(), domain.c_ptr(), orig);
+                m_cache.insert(new_name, product_pred);
+            }
+            return;
+        }
+
+        // -- compute Cartesian product of decls, and create a new
+        // -- predicate for each element of the product
+        for (auto &p : *decls[idx]) {
+            decls_buf[idx] = p;
+            add_new_rel_symbols(idx + 1, decls, decls_buf, was_added);
+        }
+    }
+
+    void mk_synchronize::replace_applications(rule & r, rule_set & rules,
+                                              ptr_vector<app> & apps) {
+        app_ref replacing = product_application(apps);
+
+        ptr_vector<app> new_tail;
+        svector<bool> new_tail_neg;
+        unsigned n = r.get_tail_size() - apps.size() + 1;
+        unsigned tail_idx = 0;
+        new_tail.resize(n);
+        new_tail_neg.resize(n);
+        new_tail[0] = replacing;
+        new_tail_neg[0] = false;
+
+        for (unsigned i = 0; i < r.get_positive_tail_size(); ++i) {
+            app* tail = r.get_tail(i);
+            if (!apps.contains(tail)) {
+                ++tail_idx;
+                new_tail[tail_idx] = tail;
+                new_tail_neg[tail_idx] = false;
+            }
+        }
+        for (unsigned i = r.get_positive_tail_size(); i < r.get_uninterpreted_tail_size(); ++i) {
+            ++tail_idx;
+            new_tail[tail_idx] = r.get_tail(i);
+            new_tail_neg[tail_idx] = true;
+        }
+        for (unsigned i = r.get_uninterpreted_tail_size(); i < r.get_tail_size(); ++i) {
+            ++tail_idx;
+            new_tail[tail_idx] = r.get_tail(i);
+            new_tail_neg[tail_idx] = false;
+        }
+
+        rule_ref new_rule(rm);
+        new_rule = rm.mk(r.get_head(), tail_idx + 1,
+            new_tail.c_ptr(), new_tail_neg.c_ptr(), symbol::null, false);
+        rules.replace_rule(&r, new_rule.get());
+    }
+
+    rule_ref mk_synchronize::rename_bound_vars_in_rule(rule * r,
+                                                       unsigned & var_idx) {
+        // AG: shift all variables in a rule so that lowest var index is var_idx?
+        // AG: update var_idx in the process?
+        ptr_vector<sort> sorts;
+        r->get_vars(m, sorts);
+        expr_ref_vector revsub(m);
+        revsub.resize(sorts.size());
+        for (unsigned i = 0; i < sorts.size(); ++i) {
+            if (sorts[i]) {
+                revsub[i] = m.mk_var(var_idx++, sorts[i]);
+            }
+        }
+
+        rule_ref new_rule(rm);
+        new_rule = rm.mk(r);
+        rm.substitute(new_rule, revsub.size(), revsub.c_ptr());
+        return new_rule;
+    }
+
+    vector<rule_ref_vector> mk_synchronize::rename_bound_vars(item_set_vector const & heads,
+                                                              rule_set & rules) {
+        // AG: is every item_set in heads corresponds to rules that are merged?
+        // AG: why are bound variables renamed in the first place?
+        // AG: the data structure seems too complex
+        vector<rule_ref_vector> result;
+        unsigned var_idx = 0;
+        for (auto item : heads) {
+            rule_ref_vector dst_vector(rm);
+            for (auto *head : *item) {
+                for (auto *r : rules.get_predicate_rules(head)) {
+                    rule_ref new_rule = rename_bound_vars_in_rule(r, var_idx);
+                    dst_vector.push_back(new_rule.get());
+                }
+            }
+            result.push_back(dst_vector);
+        }
+        return result;
+    }
+
+    void mk_synchronize::add_rec_tail(vector< ptr_vector<app> > & recursive_calls,
+                                      app_ref_vector & new_tail,
+                                      svector<bool> & new_tail_neg,
+                                      unsigned & tail_idx) {
+        unsigned max_sz = 0;
+        for (auto &rc : recursive_calls)
+            max_sz= std::max(rc.size(), max_sz);
+
+        unsigned n = recursive_calls.size();
+        ptr_vector<app> merged_recursive_calls;
+
+        for (unsigned j = 0; j < max_sz; ++j) {
+            merged_recursive_calls.reset();
+            merged_recursive_calls.resize(n);
+            for (unsigned i = 0; i < n; ++i) {
+                unsigned sz = recursive_calls[i].size();
+                merged_recursive_calls[i] =
+                    j < sz ? recursive_calls[i][j] : recursive_calls[i][sz - 1];
+            }
+            ++tail_idx;
+            new_tail[tail_idx] = product_application(merged_recursive_calls);
+            new_tail_neg[tail_idx] = false;
+        }
+    }
+
+    void mk_synchronize::add_non_rec_tail(rule & r, app_ref_vector & new_tail,
+                                          svector<bool> & new_tail_neg,
+                                          unsigned & tail_idx) {
+        for (unsigned i = 0, sz = r.get_positive_tail_size(); i < sz; ++i) {
+            app* tail = r.get_tail(i);
+            if (!is_recursive(r, *tail)) {
+                ++tail_idx;
+                new_tail[tail_idx] = tail;
+                new_tail_neg[tail_idx] = false;
+            }
+        }
+        for (unsigned i = r.get_positive_tail_size(),
+                 sz = r.get_uninterpreted_tail_size() ; i < sz; ++i) {
+            ++tail_idx;
+            new_tail[tail_idx] = r.get_tail(i);
+            new_tail_neg[tail_idx] = true;
+        }
+        for (unsigned i = r.get_uninterpreted_tail_size(),
+                 sz = r.get_tail_size(); i < sz; ++i) {
+            ++tail_idx;
+            new_tail[tail_idx] = r.get_tail(i);
+            new_tail_neg[tail_idx] = r.is_neg_tail(i);
+        }
+    }
+
+    app_ref mk_synchronize::product_application(ptr_vector<app> const &apps) {
+        unsigned args_num = 0;
+        string_buffer<> buffer;
+
+        // AG: factor out into mk_name
+        for (auto *app : apps) {
+            buffer << app->get_decl()->get_name() << "!!";
+            args_num += app->get_num_args();
+        }
+
+        symbol name = symbol(buffer.c_str());
+        SASSERT(m_cache.contains(name));
+        func_decl * pred = m_cache[name];
+
+        ptr_vector<expr> args;
+        args.resize(args_num);
+        unsigned idx = 0;
+        for (auto *a : apps) {
+            for (unsigned i = 0, sz = a->get_num_args(); i < sz; ++i, ++idx)
+                args[idx] = a->get_arg(i);
+        }
+
+        return app_ref(m.mk_app(pred, args_num, args.c_ptr()), m);
+    }
+
+    rule_ref mk_synchronize::product_rule(rule_ref_vector const & rules) {
+        unsigned n = rules.size();
+
+        string_buffer<> buffer;
+        bool first_rule = true;
+        for (auto it = rules.begin(); it != rules.end(); ++it, first_rule = false) {
+            if (!first_rule) {
+                buffer << "+";
+            }
+            buffer << (*it)->name();
+        }
+
+        ptr_vector<app> heads;
+        heads.resize(n);
+        for (unsigned i = 0; i < n; ++i) {
+            heads[i] = rules[i]->get_head();
+        }
+        app_ref product_head = product_application(heads);
+        unsigned product_tail_length = 0;
+        bool has_recursion = false;
+        vector< ptr_vector<app> > recursive_calls;
+        recursive_calls.resize(n);
+        for (unsigned i = 0; i < n; ++i) {
+            rule& rule = *rules[i];
+            product_tail_length += rule.get_tail_size();
+            for (unsigned j = 0; j < rule.get_positive_tail_size(); ++j) {
+                app* tail = rule.get_tail(j);
+                if (is_recursive(rule, *tail)) {
+                    has_recursion = true;
+                    recursive_calls[i].push_back(tail);
+                }
+            }
+            if (recursive_calls[i].empty()) {
+                recursive_calls[i].push_back(rule.get_head());
+            }
+        }
+
+        app_ref_vector new_tail(m);
+        svector<bool> new_tail_neg;
+        new_tail.resize(product_tail_length);
+        new_tail_neg.resize(product_tail_length);
+        unsigned tail_idx = -1;
+        if (has_recursion) {
+            add_rec_tail(recursive_calls, new_tail, new_tail_neg, tail_idx);
+        }
+
+        for (rule_vector::const_iterator it = rules.begin(); it != rules.end(); ++it) {
+            rule& rule = **it;
+            add_non_rec_tail(rule, new_tail, new_tail_neg, tail_idx);
+        }
+
+        rule_ref new_rule(rm);
+        new_rule = rm.mk(product_head, tail_idx + 1,
+            new_tail.c_ptr(), new_tail_neg.c_ptr(), symbol(buffer.c_str()), false);
+        rm.fix_unbound_vars(new_rule, false);
+        return new_rule;
+    }
+
+    void mk_synchronize::merge_rules(unsigned idx, rule_ref_vector & buf,
+                                     vector<rule_ref_vector> const & merged_rules,
+                                     rule_set & all_rules) {
+        if (idx >= merged_rules.size()) {
+            rule_ref product = product_rule(buf);
+            all_rules.add_rule(product.get());
+            return;
+        }
+
+        for (auto *r : merged_rules[idx]) {
+            buf[idx] = r;
+            merge_rules(idx + 1, buf, merged_rules, all_rules);
+        }
+    }
+
+    void mk_synchronize::merge_applications(rule & r, rule_set & rules) {
+        ptr_vector<app> non_recursive_preds;
+        obj_hashtable<app> apps;
+        for (unsigned i = 0; i < r.get_positive_tail_size(); ++i) {
+            app* t = r.get_tail(i);
+            if (!is_recursive(r, *t) && has_recursive_premise(t)) {
+                apps.insert(t);
+            }
+        }
+        if (apps.size() < 2) return;
+        for (auto *a : apps) non_recursive_preds.push_back(a);
+
+        item_set_vector merged_decls = add_merged_decls(non_recursive_preds);
+
+        unsigned n = non_recursive_preds.size();
+        ptr_vector<func_decl> decls_buf;
+        decls_buf.resize(n);
+        bool was_added = false;
+        add_new_rel_symbols(0, merged_decls, decls_buf, was_added);
+        if (was_added){
+            rule_ref_vector rules_buf(rm);
+            rules_buf.resize(n);
+            vector<rule_ref_vector> renamed_rules = rename_bound_vars(merged_decls, rules);
+            merge_rules(0, rules_buf, renamed_rules, rules);
+        }
+
+        replace_applications(r, rules, non_recursive_preds);
+        m_deps->populate(rules);
+        m_stratifier = alloc(rule_stratifier, *m_deps);
+    }
+
+    rule_set * mk_synchronize::operator()(rule_set const & source) {
+        rule_set* rules = alloc(rule_set, m_ctx);
+        rules->inherit_predicates(source);
+
+        for (auto *r : source) { rules->add_rule(r); }
+
+        m_deps = alloc(rule_dependencies, m_ctx);
+        m_deps->populate(*rules);
+        m_stratifier = alloc(rule_stratifier, *m_deps);
+
+        unsigned current_rule = 0;
+        while (current_rule < rules->get_num_rules()) {
+            rule *r = rules->get_rule(current_rule);
+            merge_applications(*r, *rules);
+            ++current_rule;
+        }
+
+        return rules;
+    }
+
+};
diff --git a/src/muz/transforms/dl_mk_synchronize.h b/src/muz/transforms/dl_mk_synchronize.h
new file mode 100644
index 000000000..6f4b3ca40
--- /dev/null
+++ b/src/muz/transforms/dl_mk_synchronize.h
@@ -0,0 +1,134 @@
+/*++
+Copyright (c) 2017-2018 Saint-Petersburg State University
+
+Module Name:
+
+    dl_mk_synchronize.h
+
+Abstract:
+
+    Rule transformer that attempts to merge recursive iterations
+    relaxing the shape of the inductive invariant.
+
+Example:
+
+    Q(z)  :- A(x), B(y), phi1(x,y,z).
+    A(x)  :- A(x'), phi2(x,x').
+    A(x)  :- phi3(x).
+    B(y)  :- C(y'), phi4(y,y').
+    C(y)  :- B(y'), phi5(y,y').
+    B(y)  :- phi6(y).
+
+    Transformed clauses:
+
+    Q(z)    :- AB(x,y), phi1(x,y,z).
+    AB(x,y) :- AC(x',y'), phi2(x,x'), phi4(y,y').
+    AC(x,y) :- AB(x',y'), phi2(x,x'), phi5(y,y').
+    AB(x,y) :- AC(x, y'), phi3(x), phi4(y,y').
+    AC(x,y) :- AB(x, y'), phi3(x), phi5(y,y').
+    AB(x,y) :- AB(x',y), phi2(x,x'), phi6(y).
+    AB(x,y) :- phi3(x), phi6(y).
+
+Author:
+
+    Dmitry Mordvinov (dvvrd) 2017-05-24
+    Lidiia Chernigovskaia (LChernigovskaya) 2017-10-20
+
+Revision History:
+
+--*/
+#ifndef DL_MK_SYNCHRONIZE_H_
+#define DL_MK_SYNCHRONIZE_H_
+
+#include"muz/base/dl_context.h"
+#include"muz/base/dl_rule_set.h"
+#include"util/uint_set.h"
+#include"muz/base/dl_rule_transformer.h"
+#include"muz/transforms/dl_mk_rule_inliner.h"
+
+namespace datalog {
+
+    /**
+       \brief Implements a synchronous product transformation.
+    */
+    class mk_synchronize : public rule_transformer::plugin {
+    public:
+        typedef ptr_vector<rule_stratifier::item_set> item_set_vector;
+    private:
+        context&        m_ctx;
+        ast_manager&    m;
+        rule_manager&   rm;
+
+        scoped_ptr<rule_dependencies> m_deps;
+        scoped_ptr<rule_stratifier> m_stratifier;
+
+        // symbol table that maps new predicate names to corresponding
+        // func_decl
+        map<symbol, func_decl*, symbol_hash_proc, symbol_eq_proc> m_cache;
+
+        /// returns true if decl is recursive in the given rule
+        /// requires that decl appears in the tail of r
+        bool is_recursive(rule &r, func_decl &decl) const;
+        bool is_recursive(rule &r, expr &e) const {
+            SASSERT(is_app(&e));
+            return is_app(&e) && is_recursive(r, *to_app(&e)->get_decl());
+        }
+
+        /// returns true if the top-level predicate of app is recursive
+        bool has_recursive_premise(app * app) const;
+
+        item_set_vector add_merged_decls(ptr_vector<app> & apps);
+
+        /**
+            Compute Cartesian product of decls and create a new
+            predicate for each element. For example, if decls is
+
+            ( (a, b), (c, d) ) )
+
+            create new predicates: a!!c, a!!d, b!!c, and b!!d
+        */
+        void add_new_rel_symbols(unsigned idx, item_set_vector const & decls,
+                                 ptr_vector<func_decl> & buf, bool & was_added);
+
+        /**
+           Convert vector of predicate apps into a single app. For example,
+           (Foo a) (Bar b) becomes (Foo!!Bar a b)
+         */
+        app_ref product_application(ptr_vector<app> const & apps);
+
+        /**
+            Replace a given rule by a rule in which conjunction of
+            predicates is replaced by a single product predicate
+         */
+        void replace_applications(rule & r, rule_set & rules,
+                                  ptr_vector<app> & apps);
+
+        rule_ref rename_bound_vars_in_rule(rule * r, unsigned & var_idx);
+        vector<rule_ref_vector> rename_bound_vars(item_set_vector const & heads,
+                                                  rule_set & rules);
+
+        void add_rec_tail(vector< ptr_vector<app> > & recursive_calls,
+                          app_ref_vector & new_tail,
+                          svector<bool> & new_tail_neg, unsigned & tail_idx);
+        void add_non_rec_tail(rule & r, app_ref_vector & new_tail,
+                              svector<bool> & new_tail_neg,
+                              unsigned & tail_idx);
+
+        rule_ref product_rule(rule_ref_vector const & rules);
+
+        void merge_rules(unsigned idx, rule_ref_vector & buf,
+            vector<rule_ref_vector> const & merged_rules, rule_set & all_rules);
+        void merge_applications(rule & r, rule_set & rules);
+
+    public:
+        /**
+           \brief Create synchronous product transformer.
+         */
+        mk_synchronize(context & ctx, unsigned priority = 22500);
+
+        rule_set * operator()(rule_set const & source);
+    };
+
+};
+
+#endif /* DL_MK_SYNCHRONIZE_H_ */
diff --git a/src/opt/opt_context.cpp b/src/opt/opt_context.cpp
index e5c0bcddb..b1296f71e 100644
--- a/src/opt/opt_context.cpp
+++ b/src/opt/opt_context.cpp
@@ -351,7 +351,7 @@ namespace opt {
     void context::get_model_core(model_ref& mdl) {
         mdl = m_model;
         fix_model(mdl);
-        mdl->set_model_completion(true);
+        if (mdl) mdl->set_model_completion(true);
         TRACE("opt", tout << *mdl;);
     }
 
@@ -1084,11 +1084,7 @@ namespace opt {
             }
             term = m_arith.mk_add(args.size(), args.c_ptr());
         }
-        else if (m_arith.is_arith_expr(term) && !is_mul_const(term)) {
-            TRACE("opt", tout << "Purifying " << term << "\n";);
-            term = purify(fm, term);
-        }
-        else if (m.is_ite(term)) {
+        else if (m.is_ite(term) || !is_mul_const(term)) {
             TRACE("opt", tout << "Purifying " << term << "\n";);
             term = purify(fm, term);
         }
diff --git a/src/qe/qsat.cpp b/src/qe/qsat.cpp
index c87b1c2eb..2ad5b9b96 100644
--- a/src/qe/qsat.cpp
+++ b/src/qe/qsat.cpp
@@ -1266,9 +1266,9 @@ namespace qe {
                 in->reset();
                 in->inc_depth();
                 result.push_back(in.get());
-                if (in->models_enabled()) {
+                if (in->models_enabled()) {                    
                     model_converter_ref mc;
-                    mc = model2model_converter(m_model.get());
+                    mc = model2model_converter(m_model_save.get());
                     mc = concat(m_pred_abs.fmc(), mc.get());
                     in->add(mc.get());
                 }
diff --git a/src/sat/ba_solver.cpp b/src/sat/ba_solver.cpp
index 4226b9791..de6635d09 100644
--- a/src/sat/ba_solver.cpp
+++ b/src/sat/ba_solver.cpp
@@ -1551,10 +1551,10 @@ namespace sat {
         if (k == 1 && lit == null_literal) {
             literal_vector _lits(lits);
             s().mk_clause(_lits.size(), _lits.c_ptr(), learned);
-            return 0;
+            return nullptr;
         }
         if (!learned && clausify(lit, lits.size(), lits.c_ptr(), k)) {
-            return 0;
+            return nullptr;
         }
         void * mem = m_allocator.allocate(card::get_obj_size(lits.size()));
         card* c = new (mem) card(next_id(), lit, lits, k);
@@ -1615,7 +1615,7 @@ namespace sat {
         bool units = true;
         for (wliteral wl : wlits) units &= wl.first == 1;
         if (k == 0 && lit == null_literal) {
-            return 0;
+            return nullptr;
         }
         if (units || k == 1) {
             literal_vector lits;
@@ -3405,7 +3405,7 @@ namespace sat {
             return;
         }
         for (wliteral l : p1) {
-            SASSERT(m_weights[l.second.index()] == 0);
+            SASSERT(m_weights.size() <= l.second.index() || m_weights[l.second.index()] == 0);
             m_weights.setx(l.second.index(), l.first, 0);
             mark_visited(l.second);  
         }
@@ -3837,8 +3837,8 @@ namespace sat {
         reset_active_var_set();
         m_wlits.reset();
         uint64_t sum = 0;
-        if (m_bound == 1) return 0;
-        if (m_overflow) return 0;
+        if (m_bound == 1) return nullptr;
+        if (m_overflow) return nullptr;
         
         for (bool_var v : m_active_vars) {
             int coeff = get_int_coeff(v);
@@ -3850,7 +3850,7 @@ namespace sat {
         }
 
         if (m_overflow || sum >= UINT_MAX/2) {
-            return 0;
+            return nullptr;
         }
         else {
             return add_pb_ge(null_literal, m_wlits, m_bound, true);
@@ -3905,7 +3905,7 @@ namespace sat {
             ++k;
         }
         if (k == 1) {
-            return 0;
+            return nullptr;
         }
         while (!m_wlits.empty()) {
             wliteral wl = m_wlits.back();
@@ -3928,7 +3928,7 @@ namespace sat {
                 ++num_max_level;
             }
         }
-        if (m_overflow) return 0;
+        if (m_overflow) return nullptr;
 
         if (slack >= k) {
 #if 0
@@ -3937,7 +3937,7 @@ namespace sat {
             std::cout << "not asserting\n";
             display(std::cout, m_A, true);
 #endif
-            return 0;
+            return nullptr;
         }
 
         // produce asserting cardinality constraint
diff --git a/src/sat/sat_lookahead.cpp b/src/sat/sat_lookahead.cpp
index 3833e2a52..4ca2c5f84 100644
--- a/src/sat/sat_lookahead.cpp
+++ b/src/sat/sat_lookahead.cpp
@@ -16,6 +16,8 @@ Author:
 
 Notes:
 
+
+
 --*/
 
 #include <cmath>
diff --git a/src/sat/sat_params.pyg b/src/sat/sat_params.pyg
index 89776c479..113a8133e 100644
--- a/src/sat/sat_params.pyg
+++ b/src/sat/sat_params.pyg
@@ -50,17 +50,39 @@ def_module_params('sat',
 	                  ('unit_walk', BOOL, False, 'use unit-walk search instead of CDCL'),
                           ('unit_walk_threads', UINT, 0, 'number of unit-walk search threads to find satisfiable solution'),
                           ('lookahead.cube.cutoff', SYMBOL, 'depth', 'cutoff type used to create lookahead cubes: depth, freevars, psat, adaptive_freevars, adaptive_psat'),
+                          # - depth: the maximal cutoff is fixed to the value of lookahead.cube.depth. 
+                          #          So if the value is 10, at most 1024 cubes will be generated of length 10.
+                          # - freevars: cutoff based on a variable fraction of lookahead.cube.freevars. 
+                          #          Cut if the number of current unassigned variables drops below a fraction of number of initial variables.
+                          # - psat:  Let psat_heur := (Sum_{clause C} (psat.clause_base ^ {-|C|+1})) / |freevars|^psat.var_exp
+                          #          Cut if the value of psat_heur exceeds psat.trigger
+                          # - adaptive_freevars: Cut if the number of current unassigned variables drops below a fraction of free variables
+                          #          at the time of the last conflict. The fraction is increased every time the a cutoff is created.
+                          # - adative_psat: Cut based on psat_heur in an adaptive way.
                           ('lookahead.cube.fraction', DOUBLE, 0.4, 'adaptive fraction to create lookahead cubes. Used when lookahead.cube.cutoff is adaptive_freevars or adaptive_psat'),
                           ('lookahead.cube.depth', UINT, 1, 'cut-off depth to create cubes. Used when lookahead.cube.cutoff is depth.'),
-                          ('lookahead.cube.freevars', DOUBLE, 0.8, 'cube free fariable fraction. Used when lookahead.cube.cutoff is freevars'),
+                          ('lookahead.cube.freevars', DOUBLE, 0.8, 'cube free variable fraction. Used when lookahead.cube.cutoff is freevars'),
                           ('lookahead.cube.psat.var_exp', DOUBLE, 1, 'free variable exponent for PSAT cutoff'),
                           ('lookahead.cube.psat.clause_base', DOUBLE, 2, 'clause base for PSAT cutoff'),
                           ('lookahead.cube.psat.trigger', DOUBLE, 5, 'trigger value to create lookahead cubes for PSAT cutoff. Used when lookahead.cube.cutoff is psat'),
-                          ('lookahead_search', BOOL, False, 'use lookahead solver'),
                           ('lookahead.preselect', BOOL, False, 'use pre-selection of subset of variables for branching'),
                           ('lookahead_simplify', BOOL, False, 'use lookahead solver during simplification'),
                           ('lookahead.use_learned', BOOL, False, 'use learned clauses when selecting lookahead literal'),
                           ('lookahead_simplify.bca', BOOL, True, 'add learned binary clauses as part of lookahead simplification'),
                           ('lookahead.global_autarky', BOOL, False, 'prefer to branch on variables that occur in clauses that are reduced'),
-                          ('lookahead.reward', SYMBOL, 'march_cu', 'select lookahead heuristic: ternary, heule_schur (Heule Schur), heuleu (Heule Unit), unit, or march_cu')))
+                          ('lookahead.reward', SYMBOL, 'march_cu', 'select lookahead heuristic: ternary, heule_schur (Heule Schur), heuleu (Heule Unit), unit, or march_cu'))
+                          # reward function used to determine which literal to cube on.
+                          # - ternary: reward function useful for random 3-SAT instances. Used by Heule and Knuth in March.
+                          # - heule_schur: reward function based on "Schur Number 5", Heule, AAAI 2018
+                          #   The score of a literal lit is: 
+                          #           Sum_{C in Clauses | lit in C} 2 ^ (- |C|+1)
+                          #           * Sum_{lit' in C | lit' != lit} lit_occs(~lit')
+                          #           / | C |
+                          #   where lit_occs(lit) is the number of clauses containing lit.
+                          # - heuleu: The score of a literal lit is: Sum_{C in Clauses | lit in C} 2 ^ (-|C| + 1)
+                          # - unit: heule_schur + also counts number of unit clauses.
+                          # - march_cu: default reward function used in a version of March
+                          # Each reward function also comes with its own variant of "mix_diff", which 
+                          # is the function for combining reward metrics for the positive and negative variant of a literal. 
+                          )
 
diff --git a/src/sat/sat_solver/inc_sat_solver.cpp b/src/sat/sat_solver/inc_sat_solver.cpp
index ff55598c2..5de64b496 100644
--- a/src/sat/sat_solver/inc_sat_solver.cpp
+++ b/src/sat/sat_solver/inc_sat_solver.cpp
@@ -259,7 +259,7 @@ public:
         return m_num_scopes;
     }
 
-    void assert_expr_core2(expr * t, expr * a) override {
+    void assert_expr_core2(expr * t, expr * a) override {        
         if (a) {
             m_asmsf.push_back(a);
             assert_expr_core(m.mk_implies(a, t));
@@ -311,7 +311,10 @@ public:
     expr_ref_vector cube(expr_ref_vector& vs, unsigned backtrack_level) override {
         if (!is_internalized()) {
             lbool r = internalize_formulas();
-            if (r != l_true) return expr_ref_vector(m);
+            if (r != l_true) {
+                IF_VERBOSE(0, verbose_stream() << "internalize produced " << r << "\n");
+                return expr_ref_vector(m);
+            }
         }
         convert_internalized();
         obj_hashtable<expr> _vs;
@@ -329,6 +332,7 @@ public:
             return result;
         }
         if (result == l_true) {
+            IF_VERBOSE(1, verbose_stream() << "formulas are SAT\n");
             return expr_ref_vector(m);
         }        
         expr_ref_vector fmls(m);
@@ -345,6 +349,7 @@ public:
                 vs.push_back(x);
             }
         }
+        if (fmls.empty()) { IF_VERBOSE(0, verbose_stream() << "no literals were produced in cube\n"); }
         return fmls;
     }
     
@@ -473,6 +478,7 @@ public:
     }
 
     void convert_internalized() {
+        m_solver.pop_to_base_level();
         if (!is_internalized() && m_fmls_head > 0) {
             internalize_formulas();
         }
diff --git a/src/shell/main.cpp b/src/shell/main.cpp
index bb1c19b47..1c8b6908d 100644
--- a/src/shell/main.cpp
+++ b/src/shell/main.cpp
@@ -26,7 +26,7 @@ Revision History:
 #include "shell/smtlib_frontend.h"
 #include "shell/z3_log_frontend.h"
 #include "util/warning.h"
-#include "util/version.h"
+#include "util/z3_version.h"
 #include "shell/dimacs_frontend.h"
 #include "shell/datalog_frontend.h"
 #include "shell/opt_frontend.h"
diff --git a/src/smt/CMakeLists.txt b/src/smt/CMakeLists.txt
index fb1997fb4..26b099240 100644
--- a/src/smt/CMakeLists.txt
+++ b/src/smt/CMakeLists.txt
@@ -13,6 +13,7 @@ z3_add_component(smt
     old_interval.cpp
     qi_queue.cpp
     smt_almost_cg_table.cpp
+    smt_arith_value.cpp
     smt_case_split_queue.cpp
     smt_cg_table.cpp
     smt_checker.cpp
diff --git a/src/smt/smt_arith_value.cpp b/src/smt/smt_arith_value.cpp
new file mode 100644
index 000000000..41679851a
--- /dev/null
+++ b/src/smt/smt_arith_value.cpp
@@ -0,0 +1,97 @@
+/*++
+Copyright (c) 2018 Microsoft Corporation
+
+Module Name:
+
+    smt_arith_value.cpp
+
+Abstract:
+
+    Utility to extract arithmetic values from context.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2018-12-08.
+
+Revision History:
+
+--*/
+
+#include "smt/smt_arith_value.h"
+#include "smt/theory_lra.h"
+#include "smt/theory_arith.h"
+
+namespace smt {
+
+    arith_value::arith_value(context& ctx): 
+        m_ctx(ctx), m(ctx.get_manager()), a(m) {}
+
+    bool arith_value::get_lo(expr* e, rational& lo, bool& is_strict) {
+        if (!m_ctx.e_internalized(e)) return false;
+        family_id afid = a.get_family_id();
+        is_strict = false;
+        enode* next = m_ctx.get_enode(e), *n = next;
+        bool found = false;
+        bool is_strict1; 
+        rational lo1;
+        theory* th = m_ctx.get_theory(afid);
+        theory_mi_arith* tha = dynamic_cast<theory_mi_arith*>(th);
+        theory_i_arith*  thi = dynamic_cast<theory_i_arith*>(th);
+        theory_lra*      thr = dynamic_cast<theory_lra*>(th);
+        do {
+            if ((tha && tha->get_lower(next, lo1, is_strict1)) ||
+                (thi && thi->get_lower(next, lo1, is_strict1)) ||
+                (thr && thr->get_lower(next, lo1, is_strict1))) {
+                if (!found || lo1 > lo || (lo == lo1 && is_strict1)) lo = lo1, is_strict = is_strict1;
+                found = true;
+            }
+            next = next->get_next();
+        }
+        while (n != next);
+        return found;
+    }
+
+    bool arith_value::get_up(expr* e, rational& up, bool& is_strict) {
+        if (!m_ctx.e_internalized(e)) return false;
+        family_id afid = a.get_family_id();
+        is_strict = false;
+        enode* next = m_ctx.get_enode(e), *n = next;
+        bool found = false, is_strict1;
+        rational up1;
+        theory* th = m_ctx.get_theory(afid);
+        theory_mi_arith* tha = dynamic_cast<theory_mi_arith*>(th);
+        theory_i_arith*  thi = dynamic_cast<theory_i_arith*>(th);
+        theory_lra*      thr = dynamic_cast<theory_lra*>(th);
+        do {
+            if ((tha && tha->get_upper(next, up1, is_strict1)) ||
+                (thi && thi->get_upper(next, up1, is_strict1)) ||
+                (thr && thr->get_upper(next, up1, is_strict1))) {
+                if (!found || up1 < up || (up1 == up && is_strict1)) up = up1, is_strict = is_strict1;
+                found = true;
+            }
+            next = next->get_next();
+        }
+        while (n != next);
+        return found;
+    }
+
+    bool arith_value::get_value(expr* e, rational& val) {
+        if (!m_ctx.e_internalized(e)) return false;
+        expr_ref _val(m);
+        enode* next = m_ctx.get_enode(e), *n = next;
+        family_id afid = a.get_family_id();
+        theory* th = m_ctx.get_theory(afid);
+        theory_mi_arith* tha = dynamic_cast<theory_mi_arith*>(th);
+        theory_i_arith* thi = dynamic_cast<theory_i_arith*>(th);
+        theory_lra* thr = dynamic_cast<theory_lra*>(th);           
+        do { 
+            e = next->get_owner();
+            if (tha && tha->get_value(next, _val) && a.is_numeral(_val, val)) return true;
+            if (thi && thi->get_value(next, _val) && a.is_numeral(_val, val)) return true;
+            if (thr && thr->get_value(next, val)) return true;
+            next = next->get_next();
+        }
+        while (next != n);
+        return false;
+    }
+};
diff --git a/src/smt/smt_arith_value.h b/src/smt/smt_arith_value.h
new file mode 100644
index 000000000..4e22b44f8
--- /dev/null
+++ b/src/smt/smt_arith_value.h
@@ -0,0 +1,37 @@
+
+/*++
+Copyright (c) 2018 Microsoft Corporation
+
+Module Name:
+
+    smt_arith_value.h
+
+Abstract:
+
+    Utility to extract arithmetic values from context.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2018-12-08.
+
+Revision History:
+
+--*/
+#pragma once
+
+#include "ast/arith_decl_plugin.h"
+#include "smt/smt_context.h"
+
+
+namespace smt {
+    class arith_value {
+        context& m_ctx;        
+        ast_manager& m;
+        arith_util a;
+    public:
+        arith_value(context& ctx);
+        bool get_lo(expr* e, rational& lo, bool& strict);
+        bool get_up(expr* e, rational& up, bool& strict);
+        bool get_value(expr* e, rational& value);
+    };
+};
diff --git a/src/smt/smt_context_pp.cpp b/src/smt/smt_context_pp.cpp
index 2a46fd07f..fb67d91d6 100644
--- a/src/smt/smt_context_pp.cpp
+++ b/src/smt/smt_context_pp.cpp
@@ -426,6 +426,7 @@ namespace smt {
         std::stringstream strm;
         strm << "lemma_" << (++m_lemma_id) << ".smt2";
         std::ofstream out(strm.str());
+        TRACE("lemma", tout << strm.str() << "\n";);
         display_lemma_as_smt_problem(out, num_antecedents, antecedents, consequent, logic);
         out.close();
         return m_lemma_id;
@@ -466,6 +467,7 @@ namespace smt {
         std::stringstream strm;
         strm << "lemma_" << (++m_lemma_id) << ".smt2";
         std::ofstream out(strm.str());
+        TRACE("lemma", tout << strm.str() << "\n";);
         display_lemma_as_smt_problem(out, num_antecedents, antecedents, num_eq_antecedents, eq_antecedents, consequent, logic);
         out.close();
         return m_lemma_id;
diff --git a/src/smt/smt_farkas_util.cpp b/src/smt/smt_farkas_util.cpp
index ff415ad0c..f7aaea61b 100644
--- a/src/smt/smt_farkas_util.cpp
+++ b/src/smt/smt_farkas_util.cpp
@@ -312,6 +312,13 @@ namespace smt {
     }
     
     expr_ref farkas_util::get() {
+        TRACE("arith", 
+              for (unsigned i = 0; i < m_coeffs.size(); ++i) {
+                  tout << m_coeffs[i] << " * (" << mk_pp(m_ineqs[i].get(), m) << ") ";
+              }
+              tout << "\n";
+              );
+
         m_normalize_factor = rational::one();
         expr_ref res(m);
         if (m_coeffs.empty()) {
@@ -330,13 +337,12 @@ namespace smt {
             partition_ineqs();
             expr_ref_vector lits(m);
             unsigned lo = 0;
-            for (unsigned i = 0; i < m_his.size(); ++i) {
-                unsigned hi = m_his[i];
+            for (unsigned hi : m_his) {
                 lits.push_back(extract_consequence(lo, hi));
                 lo = hi;
             }
             bool_rewriter(m).mk_or(lits.size(), lits.c_ptr(), res);
-            IF_VERBOSE(2, { if (lits.size() > 1) { verbose_stream() << "combined lemma: " << mk_pp(res, m) << "\n"; } });
+            IF_VERBOSE(2, { if (lits.size() > 1) { verbose_stream() << "combined lemma: " << res << "\n"; } });
         }
         else {
             res = extract_consequence(0, m_coeffs.size());
diff --git a/src/smt/smt_internalizer.cpp b/src/smt/smt_internalizer.cpp
index 1414cb522..618450f9d 100644
--- a/src/smt/smt_internalizer.cpp
+++ b/src/smt/smt_internalizer.cpp
@@ -216,13 +216,17 @@ namespace smt {
         SASSERT(m_manager.is_bool(n));
         if (is_gate(m_manager, n)) {
             switch(to_app(n)->get_decl_kind()) {
-            case OP_AND:
-                UNREACHABLE();
+            case OP_AND: {
+                for (expr * arg : *to_app(n)) {
+                    internalize(arg, true);
+                    literal lit = get_literal(arg);
+                    mk_root_clause(1, &lit, pr);
+                }
+                break;
+            }
             case OP_OR: {
                 literal_buffer lits;
-                unsigned num = to_app(n)->get_num_args();
-                for (unsigned i = 0; i < num; i++) {
-                    expr * arg = to_app(n)->get_arg(i); 
+                for (expr * arg : *to_app(n)) { 
                     internalize(arg, true);
                     lits.push_back(get_literal(arg));
                 }
@@ -294,8 +298,7 @@ namespace smt {
         sort * s = m_manager.get_sort(n->get_arg(0));
         sort_ref u(m_manager.mk_fresh_sort("distinct-elems"), m_manager);
         func_decl_ref f(m_manager.mk_fresh_func_decl("distinct-aux-f", "", 1, &s, u), m_manager);
-        for (unsigned i = 0; i < num_args; i++) {
-            expr * arg  = n->get_arg(i);
+        for (expr * arg : *n) {
             app_ref fapp(m_manager.mk_app(f, arg), m_manager);
             app_ref val(m_manager.mk_fresh_const("unique-value", u), m_manager);
             enode * e   = mk_enode(val, false, false, true);
@@ -826,9 +829,7 @@ namespace smt {
     void context::internalize_uninterpreted(app * n) {
         SASSERT(!e_internalized(n));
         // process args
-        unsigned num = n->get_num_args();
-        for (unsigned i = 0; i < num; i++) {
-            expr * arg = n->get_arg(i); 
+        for (expr * arg : *n) {
             internalize(arg, false);
             SASSERT(e_internalized(arg));
         }
@@ -1542,10 +1543,9 @@ namespace smt {
     void context::add_and_rel_watches(app * n) {
         if (relevancy()) {
             relevancy_eh * eh = m_relevancy_propagator->mk_and_relevancy_eh(n);
-            unsigned num = n->get_num_args();
-            for (unsigned i = 0; i < num; i++) {
+            for (expr * arg : *n) {
                 // if one child is assigned to false, the and-parent must be notified
-                literal l = get_literal(n->get_arg(i));
+                literal l = get_literal(arg);
                 add_rel_watch(~l, eh);
             }
         }
@@ -1554,10 +1554,9 @@ namespace smt {
     void context::add_or_rel_watches(app * n) {
         if (relevancy()) {
             relevancy_eh * eh = m_relevancy_propagator->mk_or_relevancy_eh(n);
-            unsigned num = n->get_num_args();
-            for (unsigned i = 0; i < num; i++) {
+            for (expr * arg : *n) {
                 // if one child is assigned to true, the or-parent must be notified
-                literal l = get_literal(n->get_arg(i));
+                literal l = get_literal(arg);
                 add_rel_watch(l, eh);
             }
         }
@@ -1588,9 +1587,8 @@ namespace smt {
         TRACE("mk_and_cnstr", tout << "l: "; display_literal(tout, l); tout << "\n";);
         literal_buffer buffer;
         buffer.push_back(l);
-        unsigned num_args = n->get_num_args();
-        for (unsigned i = 0; i < num_args; i++) {
-            literal l_arg = get_literal(n->get_arg(i));
+        for (expr * arg : *n) {
+            literal l_arg = get_literal(arg);
             TRACE("mk_and_cnstr", tout << "l_arg: "; display_literal(tout, l_arg); tout << "\n";);
             mk_gate_clause(~l, l_arg);
             buffer.push_back(~l_arg);
@@ -1602,9 +1600,8 @@ namespace smt {
         literal l = get_literal(n);
         literal_buffer buffer;
         buffer.push_back(~l);
-        unsigned num_args = n->get_num_args();
-        for (unsigned i = 0; i < num_args; i++) {
-            literal l_arg = get_literal(n->get_arg(i));
+        for (expr * arg : *n) {
+            literal l_arg = get_literal(arg);
             mk_gate_clause(l, ~l_arg);
             buffer.push_back(l_arg);
         }
diff --git a/src/smt/smt_model_checker.cpp b/src/smt/smt_model_checker.cpp
index 0fea4d13d..c3af41dcf 100644
--- a/src/smt/smt_model_checker.cpp
+++ b/src/smt/smt_model_checker.cpp
@@ -47,6 +47,7 @@ namespace smt {
         m_model_finder(mf),
         m_max_cexs(1),
         m_iteration_idx(0),
+        m_has_rec_fun(false),
         m_curr_model(nullptr),
         m_pinned_exprs(m) {
     }
@@ -351,9 +352,7 @@ namespace smt {
     bool model_checker::check_rec_fun(quantifier* q, bool strict_rec_fun) {
         TRACE("model_checker", tout << mk_pp(q, m) << "\n";);
         SASSERT(q->get_num_patterns() == 2); // first pattern is the function, second is the body.
-        expr* fn = to_app(q->get_pattern(0))->get_arg(0);
-        SASSERT(is_app(fn));
-        func_decl* f = to_app(fn)->get_decl();
+        func_decl* f = m.get_rec_fun_decl(q);
 
         expr_ref_vector args(m);
         unsigned num_decls = q->get_num_decls();
@@ -433,7 +432,7 @@ namespace smt {
         TRACE("model_checker", tout << "model checker result: " << (num_failures == 0) << "\n";);
         m_max_cexs += m_params.m_mbqi_max_cexs;
 
-        if (num_failures == 0 && !m_context->validate_model()) {
+        if (num_failures == 0 && (!m_context->validate_model() || has_rec_under_quantifiers())) {
             num_failures = 1;
             // this time force expanding recursive function definitions
             // that are not forced true in the current model.
@@ -450,6 +449,43 @@ namespace smt {
         return num_failures == 0;
     }
 
+    struct has_rec_fun_proc {
+        obj_hashtable<func_decl>& m_rec_funs;
+        bool m_has_rec_fun;
+
+        bool has_rec_fun() const { return m_has_rec_fun; }
+
+        has_rec_fun_proc(obj_hashtable<func_decl>& rec_funs):
+            m_rec_funs(rec_funs),
+            m_has_rec_fun(false) {}
+
+        void operator()(app* fn) {
+            m_has_rec_fun |= m_rec_funs.contains(fn->get_decl());
+        }
+        void operator()(expr*) {}
+    };
+
+    bool model_checker::has_rec_under_quantifiers() {
+        if (!m_has_rec_fun) {
+            return false;
+        }
+        obj_hashtable<func_decl> rec_funs;
+        for (quantifier * q : *m_qm) {
+            if (m.is_rec_fun_def(q)) {
+                rec_funs.insert(m.get_rec_fun_decl(q));
+            }            
+        }
+        expr_fast_mark1 visited;
+        has_rec_fun_proc proc(rec_funs);
+        for (quantifier * q : *m_qm) {
+            if (!m.is_rec_fun_def(q)) {
+                quick_for_each_expr(proc, visited, q);
+                if (proc.has_rec_fun()) return true;
+            }            
+        }
+        return false;
+    }
+
     // 
     // (repeated from defined_names.cpp)
     // NB. The pattern for lambdas is incomplete.
@@ -479,6 +515,7 @@ namespace smt {
             }
             found_relevant = true;
             if (m.is_rec_fun_def(q)) {
+                m_has_rec_fun = true;
                 if (!check_rec_fun(q, strict_rec_fun)) {
                     TRACE("model_checker", tout << "checking recursive function failed\n";);
                     num_failures++;
diff --git a/src/smt/smt_model_checker.h b/src/smt/smt_model_checker.h
index 40a58efea..57edf3034 100644
--- a/src/smt/smt_model_checker.h
+++ b/src/smt/smt_model_checker.h
@@ -51,8 +51,10 @@ namespace smt {
         scoped_ptr<context>                         m_aux_context; // Auxiliary context used for model checking quantifiers.
         unsigned                                    m_max_cexs;
         unsigned                                    m_iteration_idx;
+        bool                                        m_has_rec_fun;
         proto_model *                               m_curr_model;
         obj_map<expr, expr *>                       m_value2expr;
+
         friend class instantiation_set;
 
         void init_aux_context();
@@ -64,6 +66,7 @@ namespace smt {
         bool add_blocking_clause(model * cex, expr_ref_vector & sks);
         bool check(quantifier * q);
         bool check_rec_fun(quantifier* q, bool strict_rec_fun);
+        bool has_rec_under_quantifiers();
         void check_quantifiers(bool strict_rec_fun, bool& found_relevant, unsigned& num_failures);
 
         struct instance {
diff --git a/src/smt/theory_arith.h b/src/smt/theory_arith.h
index 4b3ee52f6..a2c6c1191 100644
--- a/src/smt/theory_arith.h
+++ b/src/smt/theory_arith.h
@@ -1071,6 +1071,8 @@ namespace smt {
 
         bool get_lower(enode* n, expr_ref& r);
         bool get_upper(enode* n, expr_ref& r);
+        bool get_lower(enode* n, rational& r, bool &is_strict);
+        bool get_upper(enode* n, rational& r, bool &is_strict);
         bool to_expr(inf_numeral const& val, bool is_int, expr_ref& r);
 
 
diff --git a/src/smt/theory_arith_core.h b/src/smt/theory_arith_core.h
index f96b6228b..bce029753 100644
--- a/src/smt/theory_arith_core.h
+++ b/src/smt/theory_arith_core.h
@@ -484,7 +484,6 @@ namespace smt {
     void theory_arith<Ext>::mk_idiv_mod_axioms(expr * dividend, expr * divisor) {
         if (!m_util.is_zero(divisor)) {
             ast_manager & m = get_manager();
-            bool is_numeral = m_util.is_numeral(divisor);
             // if divisor is zero, then idiv and mod are uninterpreted functions.
             expr_ref div(m), mod(m), zero(m), abs_divisor(m), one(m);
             expr_ref eqz(m), eq(m), lower(m), upper(m);
@@ -3303,6 +3302,21 @@ namespace smt {
         return b && to_expr(b->get_value(), is_int(v), r);
     }
 
+
+    template<typename Ext>
+    bool theory_arith<Ext>::get_lower(enode * n, rational& r, bool& is_strict) {
+        theory_var v = n->get_th_var(get_id());
+        bound* b = (v == null_theory_var) ? nullptr : lower(v);
+        return b && (r = b->get_value().get_rational().to_rational(), is_strict = b->get_value().get_infinitesimal().is_pos(), true);
+    }
+
+    template<typename Ext>
+    bool theory_arith<Ext>::get_upper(enode * n, rational& r, bool& is_strict) {
+        theory_var v = n->get_th_var(get_id());
+        bound* b = (v == null_theory_var) ? nullptr : upper(v);
+        return b && (r = b->get_value().get_rational().to_rational(), is_strict = b->get_value().get_infinitesimal().is_neg(), true);
+    }
+
     // -----------------------------------
     //
     // Backtracking
diff --git a/src/smt/theory_arith_int.h b/src/smt/theory_arith_int.h
index ee3bd5e2e..afe527a98 100644
--- a/src/smt/theory_arith_int.h
+++ b/src/smt/theory_arith_int.h
@@ -396,7 +396,9 @@ namespace smt {
         for (; it != end; ++it) {
             if (!it->is_dead() && it->m_var != b && is_free(it->m_var)) {
                 theory_var v  = it->m_var;
-                expr * bound  = m_util.mk_ge(get_enode(v)->get_owner(), m_util.mk_numeral(rational::zero(), is_int(v)));
+                expr* e = get_enode(v)->get_owner();
+                bool _is_int = m_util.is_int(e);
+                expr * bound  = m_util.mk_ge(e, m_util.mk_numeral(rational::zero(), _is_int));
                 context & ctx = get_context();
                 ctx.internalize(bound, true);
                 ctx.mark_as_relevant(bound);
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 267412da6..5b1de851e 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -55,7 +55,7 @@ std::ostream& operator<<(std::ostream& out, bound_kind const& k) {
 }
 
 class bound { 
-    smt::bool_var     m_bv;
+    smt::bool_var    m_bv;
     smt::theory_var  m_var;
     bool             m_is_int;
     rational         m_value;
@@ -129,6 +129,7 @@ class theory_lra::imp {
 
     struct scope {
         unsigned m_bounds_lim;
+        unsigned m_idiv_lim;
         unsigned m_asserted_qhead;            
         unsigned m_asserted_atoms_lim;
         unsigned m_underspecified_lim;
@@ -154,6 +155,7 @@ class theory_lra::imp {
     ast_manager&         m;
     theory_arith_params& m_arith_params;
     arith_util           a;
+    bool                 m_has_int;
     arith_eq_adapter     m_arith_eq_adapter;
     vector<rational>     m_columns;
       
@@ -163,13 +165,13 @@ class theory_lra::imp {
         expr_ref_vector     m_terms;                     
         vector<rational>    m_coeffs;
         svector<theory_var> m_vars;
-        rational            m_coeff;
+        rational            m_offset;
         ptr_vector<expr>    m_terms_to_internalize;
         internalize_state(ast_manager& m): m_terms(m) {}
         void reset() {
             m_terms.reset();
             m_coeffs.reset();
-            m_coeff.reset();
+            m_offset.reset();
             m_vars.reset();
             m_terms_to_internalize.reset();
         }
@@ -195,7 +197,7 @@ class theory_lra::imp {
         expr_ref_vector&     terms() { return m_st.m_terms; }                     
         vector<rational>&    coeffs() { return m_st.m_coeffs; }
         svector<theory_var>& vars() { return m_st.m_vars; }
-        rational&            coeff() { return m_st.m_coeff; }
+        rational&            offset() { return m_st.m_offset; }
         ptr_vector<expr>&    terms_to_internalize() { return m_st.m_terms_to_internalize; }            
         void push(expr* e, rational c) { m_st.m_terms.push_back(e); m_st.m_coeffs.push_back(c); }
         void set_back(unsigned i) { 
@@ -214,6 +216,10 @@ class theory_lra::imp {
     svector<theory_var>      m_term_index2theory_var;   // reverse map from lp_solver variables to theory variables  
     var_coeffs               m_left_side;              // constraint left side
     mutable std::unordered_map<lp::var_index, rational> m_variable_values; // current model
+    lp::var_index m_one_var;
+    lp::var_index m_zero_var;
+    lp::var_index m_rone_var;
+    lp::var_index m_rzero_var;
 
     enum constraint_source {
         inequality_source,
@@ -229,6 +235,7 @@ class theory_lra::imp {
     svector<delayed_atom>  m_asserted_atoms;        
     expr*                  m_not_handled;
     ptr_vector<app>        m_underspecified;
+    ptr_vector<expr>       m_idiv_terms;
     unsigned_vector        m_var_trail;
     vector<ptr_vector<lp_api::bound> > m_use_list;        // bounds where variables are used.
 
@@ -328,6 +335,31 @@ class theory_lra::imp {
         }
     }
 
+    lp::var_index add_const(int c, lp::var_index& var, bool is_int) {
+        if (var != UINT_MAX) {
+            return var;
+        }
+        app_ref cnst(a.mk_numeral(rational(c), is_int), m);
+        TRACE("arith", tout << "add " << cnst << "\n";);
+        enode* e = mk_enode(cnst);
+        theory_var v = mk_var(cnst);
+        var = m_solver->add_var(v, true);
+        m_theory_var2var_index.setx(v, var, UINT_MAX);
+        m_var_index2theory_var.setx(var, v, UINT_MAX);
+        m_var_trail.push_back(v);
+        add_def_constraint(m_solver->add_var_bound(var, lp::GE, rational(c)));
+        add_def_constraint(m_solver->add_var_bound(var, lp::LE, rational(c)));
+        return var;
+    }
+
+    lp::var_index get_one(bool is_int) {
+        return add_const(1, is_int ? m_one_var : m_rone_var, is_int);
+    }
+
+    lp::var_index get_zero(bool is_int) {
+        return add_const(0, is_int ? m_zero_var : m_rzero_var, is_int);
+    }
+
 
     void found_not_handled(expr* n) {
         m_not_handled = n;
@@ -372,7 +404,7 @@ class theory_lra::imp {
         expr_ref_vector & terms = st.terms();
         svector<theory_var>& vars = st.vars();
         vector<rational>& coeffs = st.coeffs();
-        rational& coeff = st.coeff();
+        rational& offset = st.offset();
         rational r;
         expr* n1, *n2;
         unsigned index = 0;
@@ -412,7 +444,7 @@ class theory_lra::imp {
                 ++index;
             }
             else if (a.is_numeral(n, r)) {
-                coeff += coeffs[index]*r;
+                offset += coeffs[index]*r;
                 ++index;
             }
             else if (a.is_uminus(n, n1)) {
@@ -424,7 +456,6 @@ class theory_lra::imp {
                 app* t = to_app(n);
                 internalize_args(t);
                 mk_enode(t);
-				
                 theory_var v = mk_var(n);
                 coeffs[vars.size()] = coeffs[index];
                 vars.push_back(v);
@@ -442,6 +473,7 @@ class theory_lra::imp {
                 }
                 else if (a.is_idiv(n, n1, n2)) {
                     if (!a.is_numeral(n2, r) || r.is_zero()) found_not_handled(n);
+                    m_idiv_terms.push_back(n);
                     app * mod = a.mk_mod(n1, n2);
                     ctx().internalize(mod, false);
                     if (ctx().relevancy()) ctx().add_relevancy_dependency(n, mod);
@@ -451,6 +483,7 @@ class theory_lra::imp {
                     if (!is_num) {
                         found_not_handled(n);
                     }
+#if 0
                     else {
                         app_ref div(a.mk_idiv(n1, n2), m);
                         mk_enode(div);
@@ -461,7 +494,8 @@ class theory_lra::imp {
                         // abs(r) > v >= 0
                         assert_idiv_mod_axioms(u, v, w, r);
                     }
-                    if (!ctx().relevancy() && !is_num) mk_idiv_mod_axioms(n1, n2);                    
+#endif
+                    if (!ctx().relevancy()) mk_idiv_mod_axioms(n1, n2);                    
                 }
                 else if (a.is_rem(n, n1, n2)) {
                     if (!a.is_numeral(n2, r) || r.is_zero()) found_not_handled(n);
@@ -544,6 +578,7 @@ class theory_lra::imp {
     }
 
     enode * mk_enode(app * n) {
+        TRACE("arith", tout << expr_ref(n, m) << "\n";);
         if (ctx().e_internalized(n)) {
             return get_enode(n);
         }
@@ -624,6 +659,7 @@ class theory_lra::imp {
         }
         if (result == UINT_MAX) {
             result = m_solver->add_var(v, is_int(v));
+            m_has_int |= is_int(v);
             m_theory_var2var_index.setx(v, result, UINT_MAX);
             m_var_index2theory_var.setx(result, v, UINT_MAX);
             m_var_trail.push_back(v);
@@ -677,7 +713,7 @@ class theory_lra::imp {
         m_constraint_sources.setx(index, inequality_source, null_source);
         m_inequalities.setx(index, lit, null_literal);
         ++m_stats.m_add_rows;
-        TRACE("arith", m_solver->print_constraint(index, tout); tout << "\n";);
+        TRACE("arith", m_solver->print_constraint(index, tout) << "\n";);
     }
 
     void add_def_constraint(lp::constraint_index index) {
@@ -692,9 +728,7 @@ class theory_lra::imp {
         ++m_stats.m_add_rows;
     }
         
-    void internalize_eq(theory_var v1, theory_var v2) {     
-        enode* n1 = get_enode(v1);
-        enode* n2 = get_enode(v2);
+    void internalize_eq(theory_var v1, theory_var v2) {  
         app_ref term(m.mk_fresh_const("eq", a.mk_real()), m);
         scoped_internalize_state st(*this);
         st.vars().push_back(v1);
@@ -707,8 +741,8 @@ class theory_lra::imp {
         add_def_constraint(m_solver->add_var_bound(vi, lp::GE, rational::zero()));
         TRACE("arith", 
               {
-                  expr*  o1 = n1->get_owner();
-                  expr*  o2 = n2->get_owner();                  
+                  expr*  o1 = get_enode(v1)->get_owner();
+                  expr*  o2 = get_enode(v2)->get_owner();                  
                   tout << "v" << v1 << " = " << "v" << v2 << ": "
                        << mk_pp(o1, m) << " = " << mk_pp(o2, m) << "\n";
               });
@@ -733,10 +767,19 @@ class theory_lra::imp {
     }
        
     bool is_unit_var(scoped_internalize_state& st) {
-        return st.coeff().is_zero() && st.vars().size() == 1 && st.coeffs()[0].is_one();
+        return st.offset().is_zero() && st.vars().size() == 1 && st.coeffs()[0].is_one();
+    }
+
+    bool is_one(scoped_internalize_state& st) {
+        return st.offset().is_one() && st.vars().empty();
+    }
+
+    bool is_zero(scoped_internalize_state& st) {
+        return st.offset().is_zero() && st.vars().empty();
     }
 
     theory_var internalize_def(app* term, scoped_internalize_state& st) {
+        TRACE("arith", tout << expr_ref(term, m) << "\n";);
         if (ctx().e_internalized(term)) {
             IF_VERBOSE(0, verbose_stream() << "repeated term\n";);
             return mk_var(term, false);
@@ -766,13 +809,24 @@ class theory_lra::imp {
         if (is_unit_var(st)) {
             return st.vars()[0];
         }
+        else if (is_one(st)) {
+            return get_one(a.is_int(term));
+        }
+        else if (is_zero(st)) {
+            return get_zero(a.is_int(term));
+        }
         else {
             init_left_side(st);
             theory_var v = mk_var(term);
             lp::var_index vi = m_theory_var2var_index.get(v, UINT_MAX);
-            TRACE("arith", tout << mk_pp(term, m) << " " << v << " " << vi << "\n";);
+            TRACE("arith", tout << mk_pp(term, m) << " v" << v << "\n";);
             if (vi == UINT_MAX) {
-                vi = m_solver->add_term(m_left_side, st.coeff());
+                rational const& offset = st.offset();
+                if (!offset.is_zero()) {
+                    m_left_side.push_back(std::make_pair(offset, get_one(a.is_int(term))));
+                }
+                SASSERT(!m_left_side.empty());
+                vi = m_solver->add_term(m_left_side);
                 m_theory_var2var_index.setx(v, vi, UINT_MAX);
                 if (m_solver->is_term(vi)) {
                     m_term_index2theory_var.setx(m_solver->adjust_term_index(vi), v, UINT_MAX);
@@ -782,7 +836,7 @@ class theory_lra::imp {
                 }
                 m_var_trail.push_back(v);
                 TRACE("arith_verbose", tout << "v" << v << " := " << mk_pp(term, m) << " slack: " << vi << " scopes: " << m_scopes.size() << "\n";
-                      m_solver->print_term(m_solver->get_term(vi), tout); tout << "\n";);
+                      m_solver->print_term(m_solver->get_term(vi), tout) << "\n";);
             }
             rational val;
             if (a.is_numeral(term, val)) {
@@ -798,9 +852,14 @@ public:
         th(th), m(m), 
         m_arith_params(ap), 
         a(m), 
+        m_has_int(false),
         m_arith_eq_adapter(th, ap, a),            
         m_internalize_head(0),
-        m_not_handled(0),
+        m_one_var(UINT_MAX),
+        m_zero_var(UINT_MAX),
+        m_rone_var(UINT_MAX),
+        m_rzero_var(UINT_MAX),
+        m_not_handled(nullptr),
         m_asserted_qhead(0), 
         m_assume_eq_head(0),
         m_num_conflicts(0),
@@ -871,16 +930,18 @@ public:
         }
         return true;
     }
+
+    bool is_arith(enode* n) {
+        return n && n->get_th_var(get_id()) != null_theory_var;
+    }
         
     void internalize_eq_eh(app * atom, bool_var) {
-        expr* lhs = 0, *rhs = 0;
+        TRACE("arith_verbose", tout << mk_pp(atom, m) << "\n";);
+        expr* lhs = nullptr, *rhs = nullptr;
         VERIFY(m.is_eq(atom, lhs, rhs));
         enode * n1 = get_enode(lhs);
         enode * n2 = get_enode(rhs);
-        if (n1->get_th_var(get_id()) != null_theory_var &&
-            n2->get_th_var(get_id()) != null_theory_var &&
-            n1 != n2) {
-            TRACE("arith_verbose", tout << mk_pp(atom, m) << "\n";);
+        if (is_arith(n1) && is_arith(n2) && n1 != n2) {
             m_arith_eq_adapter.mk_axioms(n1, n2);
         }
     }
@@ -910,6 +971,7 @@ public:
         scope& s = m_scopes.back();
         s.m_bounds_lim = m_bounds_trail.size();
         s.m_asserted_qhead = m_asserted_qhead;
+        s.m_idiv_lim = m_idiv_terms.size();
         s.m_asserted_atoms_lim = m_asserted_atoms.size();
         s.m_not_handled = m_not_handled;
         s.m_underspecified_lim = m_underspecified.size();
@@ -935,6 +997,7 @@ public:
             }
             m_theory_var2var_index[m_var_trail[i]] = UINT_MAX;
         }
+        m_idiv_terms.shrink(m_scopes[old_size].m_idiv_lim);
         m_asserted_atoms.shrink(m_scopes[old_size].m_asserted_atoms_lim);
         m_asserted_qhead = m_scopes[old_size].m_asserted_qhead;
         m_underspecified.shrink(m_scopes[old_size].m_underspecified_lim);
@@ -1030,37 +1093,74 @@ public:
         add_def_constraint(m_solver->add_var_bound(vi, lp::LE, rational::zero()));
         add_def_constraint(m_solver->add_var_bound(get_var_index(v), lp::GE, rational::zero()));
         add_def_constraint(m_solver->add_var_bound(get_var_index(v), lp::LT, abs(r)));
+        TRACE("arith", m_solver->print_constraints(tout << term << "\n"););
     }
 
     void mk_idiv_mod_axioms(expr * p, expr * q) {
         if (a.is_zero(q)) {
             return;
         }
+        TRACE("arith", tout << expr_ref(p, m) << " " << expr_ref(q, m) << "\n";);
         // if q is zero, then idiv and mod are uninterpreted functions.
         expr_ref div(a.mk_idiv(p, q), m);
         expr_ref mod(a.mk_mod(p, q), m);
         expr_ref zero(a.mk_int(0), m);
-        literal q_ge_0     = mk_literal(a.mk_ge(q, zero));
-        literal q_le_0     = mk_literal(a.mk_le(q, zero));
-        //            literal eqz        = th.mk_eq(q, zero, false);
         literal eq         = th.mk_eq(a.mk_add(a.mk_mul(q, div), mod), p, false);
         literal mod_ge_0   = mk_literal(a.mk_ge(mod, zero));
-        // q >= 0 or p = (p mod q) + q * (p div q)
-        // q <= 0 or p = (p mod q) + q * (p div q)
-        // q >= 0 or (p mod q) >= 0
-        // q <= 0 or (p mod q) >= 0
-        // q <= 0 or (p mod q) <  q
-        // q >= 0 or (p mod q) < -q
-        // enable_trace("mk_bool_var");
-        mk_axiom(q_ge_0, eq);
-        mk_axiom(q_le_0, eq);
-        mk_axiom(q_ge_0, mod_ge_0);
-        mk_axiom(q_le_0, mod_ge_0);
-        mk_axiom(q_le_0, ~mk_literal(a.mk_ge(a.mk_sub(mod, q), zero)));
-        mk_axiom(q_ge_0, ~mk_literal(a.mk_ge(a.mk_add(mod, q), zero)));
-        rational k;
-        if (m_arith_params.m_arith_enum_const_mod && a.is_numeral(q, k) && 
-            k.is_pos() && k < rational(8)) {
+        literal div_ge_0   = mk_literal(a.mk_ge(div, zero));
+        literal div_le_0   = mk_literal(a.mk_le(div, zero));
+        literal p_ge_0     = mk_literal(a.mk_ge(p, zero));
+        literal p_le_0     = mk_literal(a.mk_le(p, zero));
+
+        rational k(0);
+        expr_ref upper(m);
+
+        if (a.is_numeral(q, k)) {
+            if (k.is_pos()) { 
+                upper = a.mk_numeral(k - 1, true);
+            }
+            else if (k.is_neg()) {
+                upper = a.mk_numeral(-k - 1, true);
+            }
+        }
+        else {
+            k = rational::zero();
+        }
+
+        if (!k.is_zero()) {
+            mk_axiom(eq);
+            mk_axiom(mod_ge_0);
+            mk_axiom(mk_literal(a.mk_le(mod, upper)));
+            if (k.is_pos()) {
+                mk_axiom(~p_ge_0, div_ge_0);
+                mk_axiom(~p_le_0, div_le_0);
+            }
+            else {
+                mk_axiom(~p_ge_0, div_le_0);
+                mk_axiom(~p_le_0, div_ge_0);
+            }
+        }
+        else {
+            // q >= 0 or p = (p mod q) + q * (p div q)
+            // q <= 0 or p = (p mod q) + q * (p div q)
+            // q >= 0 or (p mod q) >= 0
+            // q <= 0 or (p mod q) >= 0
+            // q <= 0 or (p mod q) <  q
+            // q >= 0 or (p mod q) < -q
+            literal q_ge_0 = mk_literal(a.mk_ge(q, zero));
+            literal q_le_0 = mk_literal(a.mk_le(q, zero));
+            mk_axiom(q_ge_0, eq);
+            mk_axiom(q_le_0, eq);
+            mk_axiom(q_ge_0, mod_ge_0);
+            mk_axiom(q_le_0, mod_ge_0);
+            mk_axiom(q_le_0, ~mk_literal(a.mk_ge(a.mk_sub(mod, q), zero)));
+            mk_axiom(q_ge_0, ~mk_literal(a.mk_ge(a.mk_add(mod, q), zero)));        
+            mk_axiom(q_le_0, ~p_ge_0, div_ge_0); 
+            mk_axiom(q_le_0, ~p_le_0, div_le_0); 
+            mk_axiom(q_ge_0, ~p_ge_0, div_le_0); 
+            mk_axiom(q_ge_0, ~p_le_0, div_ge_0); 
+        }
+        if (m_arith_params.m_arith_enum_const_mod && k.is_pos() && k < rational(8)) {
             unsigned _k = k.get_unsigned();
             literal_buffer lits;
             for (unsigned j = 0; j < _k; ++j) {
@@ -1152,7 +1252,6 @@ public:
             m_todo_terms.pop_back();
             if (m_solver->is_term(vi)) {
                 const lp::lar_term& term = m_solver->get_term(vi);
-                result += term.m_v * coeff;
                 for (const auto & i:  term) {
                     m_todo_terms.push_back(std::make_pair(i.var(), coeff * i.coeff()));
                 }                    
@@ -1189,7 +1288,6 @@ public:
             m_todo_terms.pop_back();
             if (m_solver->is_term(wi)) {
                 const lp::lar_term& term = m_solver->get_term(wi);
-                result += term.m_v * coeff;
                 for (const auto & i : term) {
                     if (m_variable_values.count(i.var()) > 0) {
                         result += m_variable_values[i.var()] * coeff * i.coeff();
@@ -1208,10 +1306,10 @@ public:
     }
 
     void init_variable_values() {
+        reset_variable_values();
         if (!m.canceled() && m_solver.get() && th.get_num_vars() > 0) {
-            reset_variable_values();
+            TRACE("arith", tout << "update variable values\n";);
             m_solver->get_model(m_variable_values);
-            TRACE("arith", display(tout););
         }
     }
 
@@ -1314,6 +1412,7 @@ public:
     }
 
     final_check_status final_check_eh() {
+        IF_VERBOSE(2, verbose_stream() << "final-check\n");
         m_use_nra_model = false;
         lbool is_sat = l_true;
         if (m_solver->get_status() != lp::lp_status::OPTIMAL) {
@@ -1328,7 +1427,7 @@ public:
             }
             if (assume_eqs()) {
                 return FC_CONTINUE;
-            }
+            }            
 
             switch (check_lia()) {
             case l_true:
@@ -1340,7 +1439,7 @@ public:
                 st = FC_CONTINUE;
                 break;
             }
-                
+            
             switch (check_nra()) {
             case l_true:
                 break;
@@ -1378,6 +1477,18 @@ public:
         u_map<rational> coeffs;
         term2coeffs(term, coeffs, rational::one(), offset);
         offset.neg();
+        TRACE("arith", 
+              m_solver->print_term(term, tout << "term: ") << "\n";
+              for (auto const& kv : coeffs) {
+                  tout << "v" << kv.m_key << " * " << kv.m_value << "\n";
+              }
+              tout << offset << "\n";
+              rational g(0);
+              for (auto const& kv : coeffs) {
+                  g = gcd(g, kv.m_value);
+              }
+              tout << "gcd: " << g << "\n";
+              );
         if (is_int) {
             // 3x + 6y >= 5 -> x + 3y >= 5/3, then x + 3y >= 2
             // 3x + 6y <= 5 -> x + 3y <= 1
@@ -1385,10 +1496,12 @@ public:
             rational g = gcd_reduce(coeffs);
             if (!g.is_one()) {
                 if (lower_bound) {
-                    offset = div(offset + g - rational::one(), g);
+                    TRACE("arith", tout << "lower: " << offset << " / " << g << " = " << offset / g << " >= " << ceil(offset / g) << "\n";);
+                    offset = ceil(offset / g);
                 }
                 else {
-                    offset = div(offset, g);
+                    TRACE("arith", tout << "upper: " << offset << " / " << g << " = " << offset / g << " <= " << floor(offset / g) << "\n";);
+                    offset = floor(offset / g);
                 }
             }
         }
@@ -1408,27 +1521,247 @@ public:
         }
 
         TRACE("arith", tout << t << ": " << atom << "\n";
-              m_solver->print_term(term, tout << "bound atom: "); tout << (lower_bound?" >= ":" <= ") << k << "\n";);
+              m_solver->print_term(term, tout << "bound atom: ") << (lower_bound?" >= ":" <= ") << k << "\n";);
         ctx().internalize(atom, true);
         ctx().mark_as_relevant(atom.get());
         return atom;
     }
 
+    bool make_sure_all_vars_have_bounds() {
+        if (!m_has_int) {
+            return true;
+        }
+        unsigned nv = std::min(th.get_num_vars(), m_theory_var2var_index.size());
+        bool all_bounded = true;
+        for (unsigned v = 0; v < nv; ++v) {
+            lp::var_index vi = m_theory_var2var_index[v];
+            if (vi == UINT_MAX)
+                continue;
+            if (!m_solver->is_term(vi) && !var_has_bound(vi, true) && !var_has_bound(vi, false)) {
+                lp::lar_term term;
+                term.add_monomial(rational::one(), vi);
+                app_ref b = mk_bound(term, rational::zero(), true);
+                TRACE("arith", tout << "added bound " << b << "\n";);
+                IF_VERBOSE(2, verbose_stream() << "bound: " << b << "\n");
+                all_bounded = false;
+            }
+        }
+        return all_bounded;
+    }
+
+    /**
+     * n = (div p q)
+     *
+     * (div p q) * q + (mod p q) = p, (mod p q) >= 0
+     *
+     * 0 < q => (p/q <= v(p)/v(q) => n <= floor(v(p)/v(q)))
+     * 0 < q => (v(p)/v(q) <= p/q => v(p)/v(q) - 1 < n) 
+     * 
+     */
+    bool check_idiv_bounds() {
+        if (m_idiv_terms.empty()) {
+            return true;
+        }
+        bool all_divs_valid = true;        
+        init_variable_values();
+        for (expr* n : m_idiv_terms) {
+            expr* p = nullptr, *q = nullptr;
+            VERIFY(a.is_idiv(n, p, q));
+            theory_var v  = mk_var(n);
+            theory_var v1 = mk_var(p);
+            theory_var v2 = mk_var(q);
+            rational r1 = get_value(v1);
+            rational r2;
+
+            if (!r1.is_int() || r1.is_neg()) {
+                // TBD
+                // r1 = 223/4, r2 = 2, r = 219/8 
+                // take ceil(r1), floor(r1), ceil(r2), floor(r2), for floor(r2) > 0
+                // then 
+                //      p/q <= ceil(r1)/floor(r2) => n <= div(ceil(r1), floor(r2))
+                //      p/q >= floor(r1)/ceil(r2) => n >= div(floor(r1), ceil(r2))
+                continue;
+            }
+
+            if (a.is_numeral(q, r2) && r2.is_pos()) {
+                if (get_value(v) == div(r1, r2)) continue;
+            
+                rational div_r = div(r1, r2);
+                // p <= q * div(r1, q) + q - 1 => div(p, q) <= div(r1, r2)
+                // p >= q * div(r1, q) => div(r1, q) <= div(p, q)
+                rational mul(1);
+                rational hi = r2 * div_r + r2 - 1;
+                rational lo = r2 * div_r;
+
+                // used to normalize inequalities so they 
+                // don't appear as 8*x >= 15, but x >= 2
+                expr *n1 = nullptr, *n2 = nullptr;
+                if (a.is_mul(p, n1, n2) && is_numeral(n1, mul) && mul.is_pos()) {
+                    p = n2;
+                    hi = floor(hi/mul);
+                    lo = ceil(lo/mul);
+                }
+                literal p_le_r1  = mk_literal(a.mk_le(p, a.mk_numeral(hi, true)));
+                literal p_ge_r1  = mk_literal(a.mk_ge(p, a.mk_numeral(lo, true)));
+                literal n_le_div = mk_literal(a.mk_le(n, a.mk_numeral(div_r, true)));
+                literal n_ge_div = mk_literal(a.mk_ge(n, a.mk_numeral(div_r, true)));
+                mk_axiom(~p_le_r1, n_le_div); 
+                mk_axiom(~p_ge_r1, n_ge_div);
+
+                all_divs_valid = false;
+
+                TRACE("arith",
+                      tout << r1 << " div " << r2 << "\n";
+                      literal_vector lits;
+                      lits.push_back(~p_le_r1);
+                      lits.push_back(n_le_div);
+                      ctx().display_literals_verbose(tout, lits) << "\n\n";
+                      lits[0] = ~p_ge_r1;
+                      lits[1] = n_ge_div;
+                      ctx().display_literals_verbose(tout, lits) << "\n";);                      
+                continue;
+            }
+#if 0
+
+            // TBD similar for non-linear division.
+            // better to deal with in nla_solver:
+
+            all_divs_valid = false;
+
+        
+            //  
+            //    p/q <= r1/r2 => n <= div(r1, r2)
+            // <=>
+            //    p*r2 <= q*r1 => n <= div(r1, r2)
+            //
+            //    p/q >= r1/r2 => n >= div(r1, r2)
+            // <=>
+            //    p*r2 >= r1*q => n >= div(r1, r2)
+            // 
+            expr_ref zero(a.mk_int(0), m);
+            expr_ref divc(a.mk_numeral(div(r1, r2), true), m);
+            expr_ref pqr(a.mk_sub(a.mk_mul(a.mk_numeral(r2, true), p), a.mk_mul(a.mk_numeral(r1, true), q)), m);            
+            literal pq_lhs   = ~mk_literal(a.mk_le(pqr, zero));
+            literal pq_rhs   = ~mk_literal(a.mk_ge(pqr, zero));
+            literal n_le_div = mk_literal(a.mk_le(n, divc));
+            literal n_ge_div = mk_literal(a.mk_ge(n, divc)); 
+            mk_axiom(pq_lhs, n_le_div); 
+            mk_axiom(pq_rhs, n_ge_div);
+            TRACE("arith",
+                  literal_vector lits;
+                  lits.push_back(pq_lhs);
+                  lits.push_back(n_le_div);
+                  ctx().display_literals_verbose(tout, lits) << "\n";
+                  lits[0] = pq_rhs;
+                  lits[1] = n_ge_div;
+                  ctx().display_literals_verbose(tout, lits) << "\n";);
+#endif
+        }
+        
+        return all_divs_valid;
+    }
+
+    expr_ref var2expr(lp::var_index v) {
+        std::ostringstream name;
+        name << "v" << m_solver->local2external(v);
+        return expr_ref(m.mk_const(symbol(name.str().c_str()), a.mk_int()), m);
+    }
+
+    expr_ref multerm(rational const& r, expr* e) {
+        if (r.is_one()) return expr_ref(e, m);
+        return expr_ref(a.mk_mul(a.mk_numeral(r, true), e), m);
+    }
+
+    expr_ref term2expr(lp::lar_term const& term) {
+        expr_ref t(m);
+        expr_ref_vector ts(m);
+        for (auto const& p : term) {
+            lp::var_index wi = p.var();
+            if (m_solver->is_term(wi)) {
+                ts.push_back(multerm(p.coeff(), term2expr(m_solver->get_term(wi))));
+            }
+            else {
+                ts.push_back(multerm(p.coeff(), var2expr(wi)));
+            }
+        }
+        if (ts.size() == 1) {
+            t = ts.back();
+        }
+        else {
+            t = a.mk_add(ts.size(), ts.c_ptr());
+        }
+        return t;
+    }
+
+    expr_ref constraint2fml(lp::constraint_index ci) {
+        lp::lar_base_constraint const& c = *m_solver->constraints()[ci];
+        expr_ref fml(m);
+        expr_ref_vector ts(m);
+        rational rhs = c.m_right_side;
+        for (auto cv : c.get_left_side_coefficients()) {
+            ts.push_back(multerm(cv.first, var2expr(cv.second)));
+        }
+        switch (c.m_kind) {
+        case lp::LE: fml = a.mk_le(a.mk_add(ts.size(), ts.c_ptr()), a.mk_numeral(rhs, true)); break;
+        case lp::LT: fml = a.mk_lt(a.mk_add(ts.size(), ts.c_ptr()), a.mk_numeral(rhs, true)); break;
+        case lp::GE: fml = a.mk_ge(a.mk_add(ts.size(), ts.c_ptr()), a.mk_numeral(rhs, true)); break;
+        case lp::GT: fml = a.mk_gt(a.mk_add(ts.size(), ts.c_ptr()), a.mk_numeral(rhs, true)); break;
+        case lp::EQ: fml = m.mk_eq(a.mk_add(ts.size(), ts.c_ptr()), a.mk_numeral(rhs, true)); break;
+        }
+        return fml;
+    }
+
+    void dump_cut_lemma(std::ostream& out, lp::lar_term const& term, lp::mpq const& k, lp::explanation const& ex, bool upper) {
+        m_solver->print_term(term, out << "bound: "); 
+        out << (upper?" <= ":" >= ") << k << "\n";
+        for (auto const& p : term) {
+            lp::var_index wi = p.var();
+            out << p.coeff() << " * ";
+            if (m_solver->is_term(wi)) {
+                m_solver->print_term(m_solver->get_term(wi), out) << "\n";
+            }
+            else {
+                out << "v" << m_solver->local2external(wi) << "\n";
+            }
+        }
+        for (auto const& ev : ex.m_explanation) {
+            m_solver->print_constraint(ev.second, out << ev.first << ": ");
+        }
+        expr_ref_vector fmls(m);
+        for (auto const& ev : ex.m_explanation) {
+            fmls.push_back(constraint2fml(ev.second));
+        }        
+        expr_ref t(term2expr(term), m);
+        if (upper) {
+            fmls.push_back(m.mk_not(a.mk_ge(t, a.mk_numeral(k, true))));
+        }
+        else {
+            fmls.push_back(m.mk_not(a.mk_le(t, a.mk_numeral(k, true))));
+        }
+        ast_pp_util visitor(m);
+        visitor.collect(fmls);
+        visitor.display_decls(out);
+        visitor.display_asserts(out, fmls, true);
+        out << "(check-sat)\n";            
+    }
+
     lbool check_lia() {
         if (m.canceled()) {
             TRACE("arith", tout << "canceled\n";);
             return l_undef;
         }
-        lp::lar_term term;
-        lp::mpq k;
-        lp::explanation ex; // TBD, this should be streamlined accross different explanations
-        bool upper;
-        switch(m_lia->check(term, k, ex, upper)) {
+        if (!check_idiv_bounds()) {
+            TRACE("arith", tout << "idiv bounds check\n";);
+            return l_false;
+        }
+        m_explanation.reset();
+        switch (m_lia->check()) {
         case lp::lia_move::sat:
             return l_true;
         case lp::lia_move::branch: {
             TRACE("arith", tout << "branch\n";);
-            app_ref b = mk_bound(term, k, !upper);
+            app_ref b = mk_bound(m_lia->get_term(), m_lia->get_offset(), !m_lia->is_upper());
+            IF_VERBOSE(2, verbose_stream() << "branch " << b << "\n";);
             // branch on term >= k + 1
             // branch on term <= k
             // TBD: ctx().force_phase(ctx().get_literal(b));
@@ -1440,11 +1773,13 @@ public:
             TRACE("arith", tout << "cut\n";);
             ++m_stats.m_gomory_cuts;
             // m_explanation implies term <= k
-            app_ref b = mk_bound(term, k, !upper);
+            app_ref b = mk_bound(m_lia->get_term(), m_lia->get_offset(), !m_lia->is_upper());
+            IF_VERBOSE(2, verbose_stream() << "cut " << b << "\n");
+            TRACE("arith", dump_cut_lemma(tout, m_lia->get_term(), m_lia->get_offset(), m_lia->get_explanation(), m_lia->is_upper()););
             m_eqs.reset();
             m_core.reset();
             m_params.reset();
-            for (auto const& ev : ex.m_explanation) {
+            for (auto const& ev : m_lia->get_explanation().m_explanation) {
                 if (!ev.first.is_zero()) { 
                     set_evidence(ev.second);
                 }
@@ -1457,8 +1792,9 @@ public:
             return l_false;
         }
         case lp::lia_move::conflict:
+            TRACE("arith", tout << "conflict\n";);
             // ex contains unsat core
-            m_explanation = ex.m_explanation;
+            m_explanation = m_lia->get_explanation().m_explanation;
             set_conflict1();
             return l_false;
         case lp::lia_move::undef:
@@ -2062,18 +2398,18 @@ public:
         SASSERT(!bounds.empty());
         if (bounds.size() == 1) return;
         if (m_unassigned_bounds[v] == 0) return;
-
+        bool v_is_int = is_int(v);
         literal lit1(bv, !is_true);
         literal lit2 = null_literal;
         bool find_glb = (is_true == (k == lp_api::lower_t));
+        TRACE("arith", tout << "find_glb: " << find_glb << " is_true: " << is_true << " k: " << k << " is_lower: " << (k == lp_api::lower_t) << "\n";);
         if (find_glb) {
             rational glb;
-            lp_api::bound* lb = 0;
-            for (unsigned i = 0; i < bounds.size(); ++i) {
-                lp_api::bound* b2 = bounds[i];
+            lp_api::bound* lb = nullptr;
+            for (lp_api::bound* b2 : bounds) {
                 if (b2 == &b) continue;
                 rational const& val2 = b2->get_value();
-                if ((is_true ? val2 < val : val2 <= val) && (!lb || glb < val2)) {
+                if (((is_true || v_is_int) ? val2 < val : val2 <= val) && (!lb || glb < val2)) {
                     lb = b2;
                     glb = val2;
                 }
@@ -2084,12 +2420,11 @@ public:
         }
         else {
             rational lub;
-            lp_api::bound* ub = 0;
-            for (unsigned i = 0; i < bounds.size(); ++i) {
-                lp_api::bound* b2 = bounds[i];
+            lp_api::bound* ub = nullptr;
+            for (lp_api::bound* b2 : bounds) {
                 if (b2 == &b) continue;
                 rational const& val2 = b2->get_value();
-                if ((is_true ? val < val2 : val <= val2) && (!ub || val2 < lub)) {
+                if (((is_true || v_is_int) ? val < val2 : val <= val2) && (!ub || val2 < lub)) {
                     ub = b2;
                     lub = val2;
                 }
@@ -2107,7 +2442,7 @@ public:
         m_params.reset();
         m_core.reset();
         m_eqs.reset();
-        m_core.push_back(lit2);
+        m_core.push_back(lit1);
         m_params.push_back(parameter(symbol("farkas")));
         m_params.push_back(parameter(rational(1)));
         m_params.push_back(parameter(rational(1)));
@@ -2383,6 +2718,18 @@ public:
         }
     }
 
+    bool var_has_bound(lp::var_index vi, bool is_lower) {
+        bool is_strict = false;
+        rational b;
+        lp::constraint_index ci;
+        if (is_lower) {
+            return m_solver->has_lower_bound(vi, ci, b, is_strict);
+        }
+        else {
+            return m_solver->has_upper_bound(vi, ci, b, is_strict);
+        }        
+    }
+
     bool has_upper_bound(lp::var_index vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, false); }
 
     bool has_lower_bound(lp::var_index vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, true); }
@@ -2606,7 +2953,7 @@ public:
             m_todo_terms.push_back(std::make_pair(vi, rational::one()));
 
             TRACE("arith", tout << "v" << v << " := w" << vi << "\n";
-                  m_solver->print_term(m_solver->get_term(vi), tout); tout << "\n";);
+                  m_solver->print_term(m_solver->get_term(vi), tout) << "\n";);
 
             m_nra->am().set(r, 0);
             while (!m_todo_terms.empty()) {
@@ -2614,13 +2961,13 @@ public:
                 vi = m_todo_terms.back().first;
                 m_todo_terms.pop_back();
                 lp::lar_term const& term = m_solver->get_term(vi);
-                TRACE("arith", m_solver->print_term(term, tout); tout << "\n";);
+                TRACE("arith", m_solver->print_term(term, tout) << "\n";);
                 scoped_anum r1(m_nra->am());
-                rational c1 = term.m_v * wcoeff;
+                rational c1(0);
                 m_nra->am().set(r1, c1.to_mpq());
                 m_nra->am().add(r, r1, r);                
                 for (auto const & arg : term) {
-                    lp::var_index wi = m_solver->local2external(arg.var());
+                    lp::var_index wi = arg.var();
                     c1 = arg.coeff() * wcoeff;
                     if (m_solver->is_term(wi)) {
                         m_todo_terms.push_back(std::make_pair(wi, c1));
@@ -2658,13 +3005,23 @@ public:
         }
     }
 
-    bool get_value(enode* n, expr_ref& r) {
+    bool get_value(enode* n, rational& val) {
         theory_var v = n->get_th_var(get_id());            
         if (!can_get_bound(v)) return false;
         lp::var_index vi = m_theory_var2var_index[v];
-        rational val;
         if (m_solver->has_value(vi, val)) {
+            TRACE("arith", tout << expr_ref(n->get_owner(), m) << " := " << val << "\n";);
             if (is_int(n) && !val.is_int()) return false;
+            return true;
+        }
+        else {
+            return false;
+        }
+    }    
+
+    bool get_value(enode* n, expr_ref& r) {
+        rational val;
+        if (get_value(n, val)) {
             r = a.mk_numeral(val, is_int(n));
             return true;
         }
@@ -2673,7 +3030,7 @@ public:
         }
     }    
 
-    bool get_lower(enode* n, expr_ref& r) {
+    bool get_lower(enode* n, rational& val, bool& is_strict) {
         theory_var v = n->get_th_var(get_id());
         if (!can_get_bound(v)) {
             TRACE("arith", tout << "cannot get lower for " << v << "\n";);
@@ -2681,29 +3038,36 @@ public:
         }
         lp::var_index vi = m_theory_var2var_index[v];
         lp::constraint_index ci;
-        rational val;
+        return m_solver->has_lower_bound(vi, ci, val, is_strict);
+    }
+
+    bool get_lower(enode* n, expr_ref& r) {
         bool is_strict;
-        if (m_solver->has_lower_bound(vi, ci, val, is_strict)) {
+        rational val;
+        if (get_lower(n, val, is_strict) && !is_strict) {
             r = a.mk_numeral(val, is_int(n));
             return true;
         }
-        TRACE("arith", m_solver->print_constraints(tout << "does not have lower bound " << vi << "\n"););
         return false;
     }
 
-    bool get_upper(enode* n, expr_ref& r) {
+    bool get_upper(enode* n, rational& val, bool& is_strict) {
         theory_var v = n->get_th_var(get_id());
         if (!can_get_bound(v))
             return false;
         lp::var_index vi = m_theory_var2var_index[v];
         lp::constraint_index ci;
-        rational val;
+        return m_solver->has_upper_bound(vi, ci, val, is_strict);
+
+    }
+
+    bool get_upper(enode* n, expr_ref& r) {
         bool is_strict;
-        if (m_solver->has_upper_bound(vi, ci, val, is_strict)) {
+        rational val;
+        if (get_upper(n, val, is_strict) && !is_strict) {
             r = a.mk_numeral(val, is_int(n));
             return true;
         }
-        TRACE("arith", m_solver->print_constraints(tout << "does not have upper bound " << vi << "\n"););
         return false;
     }
 
@@ -2875,7 +3239,6 @@ public:
             coeffs.find(w, c0);
             coeffs.insert(w, c0 + ti.coeff() * coeff);
         }
-        offset += coeff * term.m_v;
     }
 
     app_ref coeffs2app(u_map<rational> const& coeffs, rational const& offset, bool is_int) {
@@ -2915,15 +3278,17 @@ public:
 
     rational gcd_reduce(u_map<rational>& coeffs) {
         rational g(0);
-         for (auto const& kv : coeffs) {
-             g = gcd(g, kv.m_value);
-         }
-         if (!g.is_one() && !g.is_zero()) {
-             for (auto& kv : coeffs) {
-                 kv.m_value /= g;
-             }             
-         }
-         return g;
+        for (auto const& kv : coeffs) {
+            g = gcd(g, kv.m_value);
+        }
+        if (g.is_zero())
+            return rational::one();
+        if (!g.is_one()) {
+            for (auto& kv : coeffs) {
+                kv.m_value /= g;
+            }             
+        }
+        return g;
     }
 
     app_ref mk_obj(theory_var v) {
@@ -2951,7 +3316,7 @@ public:
         }
         if (!ctx().b_internalized(b)) {
             fm.hide(b->get_decl());
-            bool_var bv = ctx().mk_bool_var(b);
+            bool_var bv =  ctx().mk_bool_var(b);
             ctx().set_var_theory(bv, get_id());
             // ctx().set_enode_flag(bv, true);
             lp_api::bound_kind bkind = lp_api::bound_kind::lower_t;
@@ -3132,6 +3497,9 @@ void theory_lra::init_model(model_generator & m) {
 model_value_proc * theory_lra::mk_value(enode * n, model_generator & mg) {
     return m_imp->mk_value(n, mg);
 }
+bool theory_lra::get_value(enode* n, rational& r) {
+    return m_imp->get_value(n, r);
+}
 bool theory_lra::get_value(enode* n, expr_ref& r) {
     return m_imp->get_value(n, r);
 }
@@ -3141,6 +3509,12 @@ bool theory_lra::get_lower(enode* n, expr_ref& r) {
 bool theory_lra::get_upper(enode* n, expr_ref& r) {
     return m_imp->get_upper(n, r);
 }
+bool theory_lra::get_lower(enode* n, rational& r, bool& is_strict) {
+    return m_imp->get_lower(n, r, is_strict);
+}
+bool theory_lra::get_upper(enode* n, rational& r, bool& is_strict) {
+    return m_imp->get_upper(n, r, is_strict);
+}
 
 bool theory_lra::validate_eq_in_model(theory_var v1, theory_var v2, bool is_true) const {
     return m_imp->validate_eq_in_model(v1, v2, is_true);
diff --git a/src/smt/theory_lra.h b/src/smt/theory_lra.h
index 074b11ba7..23adaa557 100644
--- a/src/smt/theory_lra.h
+++ b/src/smt/theory_lra.h
@@ -78,8 +78,11 @@ namespace smt {
         model_value_proc * mk_value(enode * n, model_generator & mg) override;
 
         bool get_value(enode* n, expr_ref& r) override;
+        bool get_value(enode* n, rational& r);
         bool get_lower(enode* n, expr_ref& r);
         bool get_upper(enode* n, expr_ref& r);
+        bool get_lower(enode* n, rational& r, bool& is_strict);
+        bool get_upper(enode* n, rational& r, bool& is_strict);
 
         bool validate_eq_in_model(theory_var v1, theory_var v2, bool is_true) const override;
                 
diff --git a/src/smt/theory_seq.cpp b/src/smt/theory_seq.cpp
index b00c1565c..857691b8b 100644
--- a/src/smt/theory_seq.cpp
+++ b/src/smt/theory_seq.cpp
@@ -4588,10 +4588,10 @@ bool theory_seq::lower_bound2(expr* _e, rational& lo) {
     theory_mi_arith* tha = get_th_arith<theory_mi_arith>(ctx, m_autil.get_family_id(), e);
     if (!tha) {
         theory_i_arith* thi = get_th_arith<theory_i_arith>(ctx, m_autil.get_family_id(), e);
-        if (!thi || !thi->get_lower(ctx.get_enode(e), _lo)) return false;
+        if (!thi || !thi->get_lower(ctx.get_enode(e), _lo) || !m_autil.is_numeral(_lo, lo)) return false;
     }
     enode *ee = ctx.get_enode(e);
-    if (!tha->get_lower(ee, _lo) || m_autil.is_numeral(_lo, lo)) {
+    if (tha && (!tha->get_lower(ee, _lo) || m_autil.is_numeral(_lo, lo))) {
         enode *next = ee->get_next();
         bool flag = false;
         while (next != ee) {
diff --git a/src/solver/mus.cpp b/src/solver/mus.cpp
index 094b27ed3..4ae93a52e 100644
--- a/src/solver/mus.cpp
+++ b/src/solver/mus.cpp
@@ -217,13 +217,12 @@ struct mus::imp {
         }
 
         expr_set mss_set;
-        for (unsigned i = 0; i < mss.size(); ++i) {
-            mss_set.insert(mss[i]);
+        for (expr* e : mss) {
+            mss_set.insert(e);
         }
-        expr_set::iterator it = min_core.begin(), end = min_core.end();
-        for (; it != end; ++it) {
-            if (mss_set.contains(*it) && min_lit != *it) {
-                unknown.push_back(*it);
+        for (expr * e : min_core) {
+            if (mss_set.contains(e) && min_lit != e) {
+                unknown.push_back(e);
             }
         }
         core_literal = min_lit;
diff --git a/src/solver/solver.cpp b/src/solver/solver.cpp
index 4044c4a85..0e2128990 100644
--- a/src/solver/solver.cpp
+++ b/src/solver/solver.cpp
@@ -178,10 +178,19 @@ lbool solver::preferred_sat(expr_ref_vector const& asms, vector<expr_ref_vector>
     return check_sat(0, nullptr);
 }
 
-bool solver::is_literal(ast_manager& m, expr* e) {
-    return is_uninterp_const(e) || (m.is_not(e, e) && is_uninterp_const(e));
+
+static bool is_m_atom(ast_manager& m, expr* f) {
+    if (!is_app(f)) return true;
+    app* _f = to_app(f);
+    family_id bfid = m.get_basic_family_id();
+    if (_f->get_family_id() != bfid) return true;
+    if (_f->get_num_args() > 0 && m.is_bool(_f->get_arg(0))) return false;    
+    return m.is_eq(f) || m.is_distinct(f);
 }
 
+bool solver::is_literal(ast_manager& m, expr* e) {
+    return is_m_atom(m, e) || (m.is_not(e, e) && is_m_atom(m, e));
+}
 
 void solver::assert_expr(expr* f) {
     expr_ref fml(f, get_manager());
@@ -256,3 +265,40 @@ expr_ref_vector solver::get_units(ast_manager& m) {
 
     return result;
 }
+
+
+expr_ref_vector solver::get_non_units(ast_manager& m) {
+    expr_ref_vector result(m), fmls(m);
+    get_assertions(fmls);
+    family_id bfid = m.get_basic_family_id();
+    expr_mark marked;
+    unsigned sz0 = fmls.size();
+    for (unsigned i = 0; i < fmls.size(); ++i) {
+        expr* f = fmls.get(i);
+        if (marked.is_marked(f)) continue;
+        marked.mark(f);
+        if (!is_app(f)) {
+            if (i >= sz0) result.push_back(f);
+            continue;
+        }
+        app* _f = to_app(f);
+        if (_f->get_family_id() == bfid) {
+            // basic objects are true/false/and/or/not/=/distinct 
+            // and proof objects (that are not Boolean).
+            if (i < sz0 && m.is_not(f) && is_m_atom(m, _f->get_arg(0))) {
+                marked.mark(_f->get_arg(0));
+            }
+            else if (_f->get_num_args() > 0 && m.is_bool(_f->get_arg(0))) {
+                fmls.append(_f->get_num_args(), _f->get_args());
+            }
+            else if (i >= sz0 && is_m_atom(m, f)) {
+                result.push_back(f);
+            }
+        }
+        
+        else {
+            if (i >= sz0) result.push_back(f);
+        }
+    }
+    return result;
+}
diff --git a/src/solver/solver.h b/src/solver/solver.h
index c371be284..5329161cd 100644
--- a/src/solver/solver.h
+++ b/src/solver/solver.h
@@ -236,6 +236,8 @@ public:
     */
     expr_ref_vector get_units(ast_manager& m);
 
+    expr_ref_vector get_non_units(ast_manager& m);
+
     class scoped_push {
         solver& s;
         bool    m_nopop;
diff --git a/src/test/lp/gomory_test.h b/src/test/lp/gomory_test.h
index 03ad5b187..972466dc3 100644
--- a/src/test/lp/gomory_test.h
+++ b/src/test/lp/gomory_test.h
@@ -1,4 +1,5 @@
 namespace lp {
+#include "util/lp/lp_utils.h"
 struct gomory_test {
     gomory_test(
         std::function<std::string (unsigned)> name_function_p,
@@ -88,7 +89,7 @@ struct gomory_test {
         lp_assert(is_int(x_j));
         lp_assert(!a.is_int());
              lp_assert(f_0 > zero_of_type<mpq>() && f_0 < one_of_type<mpq>());
-        mpq f_j =  int_solver::fractional_part(a);
+        mpq f_j =  fractional_part(a);
         TRACE("gomory_cut_detail", 
               tout << a << " x_j = " << x_j << ", k = " << k << "\n";
               tout << "f_j: " << f_j << "\n";
@@ -184,7 +185,6 @@ struct gomory_test {
     }
 
     void print_term(lar_term & t, std::ostream & out) {
-        lp_assert(is_zero(t.m_v));
         vector<std::pair<mpq, unsigned>>  row;
         for (auto p : t.m_coeffs)
             row.push_back(std::make_pair(p.second, p.first));
@@ -206,7 +206,7 @@ struct gomory_test {
         unsigned x_j;
         mpq a;
         bool some_int_columns = false;
-        mpq f_0  = int_solver::fractional_part(get_value(inf_col));
+        mpq f_0  = fractional_part(get_value(inf_col));
         mpq one_min_f_0 = 1 - f_0;
         for ( auto pp : row) {
             a = pp.first;
diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp
index 6e418fe68..ff9de0e58 100644
--- a/src/test/lp/lp.cpp
+++ b/src/test/lp/lp.cpp
@@ -2667,13 +2667,20 @@ void test_term() {
     lar_solver solver;
     unsigned _x = 0;
     unsigned _y = 1;
+    unsigned _one = 2;
     var_index x = solver.add_var(_x, false);
     var_index y = solver.add_var(_y, false);
+    var_index one = solver.add_var(_one, false);
+
+    vector<std::pair<mpq, var_index>> term_one;
+    term_one.push_back(std::make_pair((int)1, one));
+    solver.add_constraint(term_one, lconstraint_kind::EQ, mpq(0));
 
     vector<std::pair<mpq, var_index>> term_ls;
     term_ls.push_back(std::pair<mpq, var_index>((int)1, x));
     term_ls.push_back(std::pair<mpq, var_index>((int)1, y));
-    var_index z = solver.add_term(term_ls, mpq(3));
+    term_ls.push_back(std::make_pair((int)3, one));
+    var_index z = solver.add_term(term_ls);
 
     vector<std::pair<mpq, var_index>> ls;
     ls.push_back(std::pair<mpq, var_index>((int)1, x));
@@ -2743,10 +2750,10 @@ void test_bound_propagation_one_small_sample1() {
     vector<std::pair<mpq, var_index>> coeffs;
     coeffs.push_back(std::pair<mpq, var_index>(1, a));
     coeffs.push_back(std::pair<mpq, var_index>(-1, c));
-    ls.add_term(coeffs, zero_of_type<mpq>());
+    ls.add_term(coeffs);
     coeffs.pop_back();
     coeffs.push_back(std::pair<mpq, var_index>(-1, b));
-    ls.add_term(coeffs, zero_of_type<mpq>());
+    ls.add_term(coeffs);
     coeffs.clear();
     coeffs.push_back(std::pair<mpq, var_index>(1, a));
     coeffs.push_back(std::pair<mpq, var_index>(-1, b));
@@ -3485,12 +3492,12 @@ void test_maximize_term() {
     vector<std::pair<mpq, var_index>> term_ls;
     term_ls.push_back(std::pair<mpq, var_index>((int)1, x));
     term_ls.push_back(std::pair<mpq, var_index>((int)-1, y));
-    unsigned term_x_min_y = solver.add_term(term_ls, mpq(0));
+    unsigned term_x_min_y = solver.add_term(term_ls);
     term_ls.clear();
     term_ls.push_back(std::pair<mpq, var_index>((int)2, x));
     term_ls.push_back(std::pair<mpq, var_index>((int)2, y));
     
-    unsigned term_2x_pl_2y = solver.add_term(term_ls, mpq(0));
+    unsigned term_2x_pl_2y = solver.add_term(term_ls);
     solver.add_var_bound(term_x_min_y,  LE, zero_of_type<mpq>());
     solver.add_var_bound(term_2x_pl_2y, LE, mpq((int)5));
     solver.find_feasible_solution();
@@ -3502,8 +3509,7 @@ void test_maximize_term() {
         std::cout<< "v[" << p.first << "] = " << p.second << std::endl;
     }
     std::cout << "calling int_solver\n";
-    lar_term t; mpq k; explanation ex; bool upper;
-    lia_move lm = i_solver.check(t, k, ex, upper);
+    lia_move lm = i_solver.check();
     VERIFY(lm == lia_move::sat);
     impq term_max;
     lp_status st = solver.maximize_term(term_2x_pl_2y, term_max);
diff --git a/src/test/pb2bv.cpp b/src/test/pb2bv.cpp
index 493d81bb7..623f7c9f1 100644
--- a/src/test/pb2bv.cpp
+++ b/src/test/pb2bv.cpp
@@ -18,6 +18,7 @@ Copyright (c) 2015 Microsoft Corporation
 #include "ast/rewriter/th_rewriter.h"
 #include "tactic/fd_solver/fd_solver.h"
 #include "solver/solver.h"
+#include "ast/arith_decl_plugin.h"
 
 static void test1() {
     ast_manager m;
@@ -194,9 +195,20 @@ static void test3() {
     }
 }
 
+static void test4() {
+    ast_manager m;
+    reg_decl_plugins(m);
+    arith_util arith(m);
+    expr_ref a(m.mk_const(symbol("a"), arith.mk_int()), m);
+    expr_ref b(m.mk_const(symbol("b"), arith.mk_int()), m);
+    expr_ref eq(m.mk_eq(a,b), m);
+    std::cout << "is_atom: " << is_atom(m, eq) << "\n";
+}
+
 void tst_pb2bv() {
     test1();
     test2();
     test3();
+    test4();
 }
 
diff --git a/src/util/CMakeLists.txt b/src/util/CMakeLists.txt
index 85b6f955c..b6abb785f 100644
--- a/src/util/CMakeLists.txt
+++ b/src/util/CMakeLists.txt
@@ -1,11 +1,12 @@
-if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/version.h")
-  message(FATAL_ERROR "\"${CMAKE_CURRENT_SOURCE_DIR}/version.h\""
+if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/z3_version.h")
+  message(FATAL_ERROR "\"${CMAKE_CURRENT_SOURCE_DIR}/z3_version.h\""
           ${z3_polluted_tree_msg}
   )
 endif()
 
 set(Z3_FULL_VERSION "\"${Z3_FULL_VERSION_STR}\"")
-configure_file(version.h.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/version.h)
+configure_file(z3_version.h.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/z3_version.h)
+
 
 z3_add_component(util
   SOURCES
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/bound_propagator.cpp b/src/util/lp/bound_propagator.cpp
index c4fa2aefa..a5c7c976a 100644
--- a/src/util/lp/bound_propagator.cpp
+++ b/src/util/lp/bound_propagator.cpp
@@ -17,10 +17,6 @@ const impq & bound_propagator::get_upper_bound(unsigned j) const {
 }
 void bound_propagator::try_add_bound(mpq  v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
     j = m_lar_solver.adjust_column_index_to_term_index(j);    
-    if (m_lar_solver.is_term(j)) {
-        // lp treats terms as not having a free coefficient, restoring it below for the outside consumption
-        v += m_lar_solver.get_term(j).m_v;
-    }
 
     lconstraint_kind kind = is_low? GE : LE;
     if (strict)
diff --git a/src/util/lp/column_info.h b/src/util/lp/column_info.h
index 407f40dfc..2a38900c1 100644
--- a/src/util/lp/column_info.h
+++ b/src/util/lp/column_info.h
@@ -69,16 +69,6 @@ public:
         m_column_index(static_cast<unsigned>(-1))
     {}
     
-    column_info(unsigned column_index) :
-        m_lower_bound_is_set(false),
-        m_lower_bound_is_strict(false),
-        m_upper_bound_is_set (false),
-        m_upper_bound_is_strict (false),
-        m_is_fixed(false),
-        m_cost(numeric_traits<T>::zero()),
-        m_column_index(column_index) {
-    }
-
     column_info(const column_info & ci) {
         m_name = ci.m_name;
         m_lower_bound_is_set = ci.m_lower_bound_is_set;
diff --git a/src/util/lp/column_namer.h b/src/util/lp/column_namer.h
index e6e8e53a2..51baf445f 100644
--- a/src/util/lp/column_namer.h
+++ b/src/util/lp/column_namer.h
@@ -33,29 +33,6 @@ public:
         print_linear_combination_of_column_indices(coeff, out);
     }
     
-    template <typename T>
-    void print_linear_combination_of_column_indices_only(const vector<std::pair<T, unsigned>> & coeffs, std::ostream & out) const {
-        bool first = true;
-        for (const auto & it : coeffs) {
-            auto val = it.first;
-            if (first) {
-                first = false;
-            } else {
-                if (numeric_traits<T>::is_pos(val)) {
-                    out << " + ";
-                } else {
-                    out << " - ";
-                    val = -val;
-                }
-            }
-            if (val == -numeric_traits<T>::one())
-                out << " - ";
-            else if (val != numeric_traits<T>::one())
-                out << T_to_string(val);
-        
-            out << "v" << it.second;
-        }
-    }
 
     
     template <typename T>
diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp
new file mode 100644
index 000000000..ad1c02625
--- /dev/null
+++ b/src/util/lp/gomory.cpp
@@ -0,0 +1,341 @@
+/*++
+  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"
+#include "util/lp/lp_utils.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 non integral value
+    const row_strip<mpq>& m_row;
+    const int_solver&     m_int_solver;
+    mpq                   m_lcm_den;
+    mpq                   m_f;
+    mpq                   m_one_minus_f;
+    mpq                   m_fj;
+    mpq                   m_one_minus_fj;
+    
+    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(unsigned j) {
+        lp_assert(is_int(j) && m_fj.is_pos());
+        TRACE("gomory_cut_detail", 
+              tout << " k = " << m_k;
+              tout << ", fj: " << m_fj << ", ";
+              tout << (at_lower(j)?"at_lower":"at_upper")<< std::endl;
+              );
+        mpq new_a;
+        if (at_lower(j)) {
+            new_a = m_fj <= m_one_minus_f ? m_fj / m_one_minus_f : ((1 - m_fj) / m_f);
+            lp_assert(new_a.is_pos());
+            m_k.addmul(new_a, lower_bound(j).x);
+            m_ex.push_justification(column_lower_bound_constraint(j));            
+        }
+        else {
+            lp_assert(at_upper(j));
+            // the upper terms are inverted: therefore we have the minus
+            new_a = - (m_fj <= m_f ? m_fj / m_f  : ((1 - m_fj) / m_one_minus_f));
+            lp_assert(new_a.is_neg());
+            m_k.addmul(new_a, upper_bound(j).x);
+            m_ex.push_justification(column_upper_bound_constraint(j));
+        }
+        m_t.add_monomial(new_a, j);
+        m_lcm_den = lcm(m_lcm_den, denominator(new_a));
+        TRACE("gomory_cut_detail", tout << "v" << j << " new_a = " << new_a << ", k = " << m_k << ", m_lcm_den = " << m_lcm_den << "\n";);
+    }
+
+    void real_case_in_gomory_cut(const mpq & a, unsigned j) {
+        TRACE("gomory_cut_detail_real", tout << "real\n";);
+        mpq new_a;
+        if (at_lower(j)) {
+            if (a.is_pos()) {
+                new_a = a / m_one_minus_f;
+            }
+            else {
+                new_a = - a / m_f;
+            }
+            m_k.addmul(new_a, lower_bound(j).x); // is it a faster operation than
+            // k += lower_bound(j).x * new_a;  
+            m_ex.push_justification(column_lower_bound_constraint(j));
+        }
+        else {
+            lp_assert(at_upper(j));
+            if (a.is_pos()) {
+                new_a =  - a / m_f; 
+            }
+            else {
+                new_a =   a / m_one_minus_f; 
+            }
+            m_k.addmul(new_a, upper_bound(j).x); //  k += upper_bound(j).x * new_a; 
+            m_ex.push_justification(column_upper_bound_constraint(j));
+        }
+        TRACE("gomory_cut_detail_real", tout << a << "*v" << j << " k: " << m_k << "\n";);
+        m_t.add_monomial(new_a, 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() {
+        lp_assert(!m_t.is_empty());
+        // k = 1 + sum of m_t at bounds
+        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 {
+            m_lcm_den = lcm(m_lcm_den, denominator(m_k));
+            lp_assert(m_lcm_den.is_pos());
+            TRACE("gomory_cut_detail", tout << "pol.size() > 1 den: " << m_lcm_den << std::endl;);
+            if (!m_lcm_den.is_one()) {
+                // normalize coefficients of integer parameters to be integers.
+                for (auto & pi: pol) {
+                    pi.first *= m_lcm_den;
+                    SASSERT(!is_int(pi.second) || pi.first.is_int());
+                }
+                m_k *= m_lcm_den;
+            }
+            // negate everything to return -pol <= -m_k
+            for (const auto & pi: pol) {
+                TRACE("gomory_cut", tout << pi.first << "* " << "v" << pi.second << "\n";);
+                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());
+    }
+
+    std::string var_name(unsigned j) const {
+        return std::string("x") + std::to_string(j);
+    }
+
+    std::ostream& dump_coeff_val(std::ostream & out, const mpq & a) const {
+        if (a.is_int()) {
+            out << a;
+        } 
+        else if ( a >= zero_of_type<mpq>())
+            out << "(/ " << numerator(a) << " " << denominator(a) << ")";
+        else {
+            out << "(- ( / " <<   numerator(-a) << " " << denominator(-a) << "))";
+        }
+        return out;
+    }
+    
+    template <typename T>
+    void dump_coeff(std::ostream & out, const T& c) const {
+        out << "( * ";
+        dump_coeff_val(out, c.coeff());
+        out << " " << var_name(c.var()) << ")";
+    }
+    
+    std::ostream& dump_row_coefficients(std::ostream & out) const {
+        mpq lc(1);
+        for (const auto& p : m_row) {
+            lc = lcm(lc, denominator(p.coeff())); 
+        }        
+        for (const auto& p : m_row) {
+            dump_coeff_val(out << " (* ", p.coeff()*lc) << " " << var_name(p.var()) << ")";
+        }
+        return out;
+    }
+    
+    void dump_the_row(std::ostream& out) const {
+        out << "; the row, excluding fixed vars\n";
+        out << "(assert ( = ( +";
+        dump_row_coefficients(out) << ") 0))\n";
+    }
+
+    void dump_declaration(std::ostream& out, unsigned v) const {
+        out << "(declare-const " << var_name(v) << (is_int(v) ? " Int" : " Real") << ")\n";
+    }
+    
+    void dump_declarations(std::ostream& out) const {
+        // for a column j the var name is vj
+        for (const auto & p : m_row) {
+            dump_declaration(out, p.var());
+        }
+        for (const auto& p : m_t) {
+            unsigned v = p.var();
+            if (m_int_solver.m_lar_solver->is_term(v)) {
+                dump_declaration(out, v);
+            }
+        }
+    }
+
+    void dump_lower_bound_expl(std::ostream & out, unsigned j) const {
+        out << "(assert (>= " << var_name(j) << " " << lower_bound(j).x << "))\n";
+    }
+    void dump_upper_bound_expl(std::ostream & out, unsigned j) const {
+        out << "(assert (<= " << var_name(j) << " " << upper_bound(j).x << "))\n";
+    }
+    
+    void dump_explanations(std::ostream& out) const {
+        for (const auto & p : m_row) {            
+            unsigned j = p.var();
+            if (j == m_inf_col || (!is_real(j) && p.coeff().is_int())) {
+                continue;
+            }
+            else if (at_lower(j)) {
+                dump_lower_bound_expl(out, j);
+            } else {
+                lp_assert(at_upper(j));
+                dump_upper_bound_expl(out, j);
+            }
+        }
+    }
+
+    std::ostream& dump_term_coefficients(std::ostream & out) const {
+        for (const auto& p : m_t) {
+            dump_coeff(out, p);
+        }
+        return out;
+    }
+    
+    std::ostream& dump_term_sum(std::ostream & out) const {
+        return dump_term_coefficients(out << "(+ ") << ")";
+    }
+    
+    std::ostream& dump_term_le_k(std::ostream & out) const {
+        return dump_term_sum(out << "(<= ") << " " << m_k << ")";
+    }
+
+    void dump_the_cut_assert(std::ostream & out) const {
+        dump_term_le_k(out << "(assert (not ") << "))\n";
+    }
+
+
+    void dump_cut_and_constraints_as_smt_lemma(std::ostream& out) const {
+        dump_declarations(out);
+        dump_the_row(out);
+        dump_explanations(out);
+        dump_the_cut_assert(out);
+        out << "(check-sat)\n";
+    }
+public:
+    lia_move create_cut() {
+        TRACE("gomory_cut",
+              tout << "applying cut at:\n"; print_linear_combination_of_column_indices_only(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;
+        m_t.clear();
+        mpq m_lcm_den(1);
+        bool some_int_columns = false;
+        mpq m_f  = fractional_part(get_value(m_inf_col));
+        TRACE("gomory_cut_detail", tout << "m_f: " << m_f << ", ";
+              tout << "1 - m_f: " << 1 - m_f << ", get_value(m_inf_col).x - m_f = " << get_value(m_inf_col).x - m_f;);
+        lp_assert(m_f.is_pos() && (get_value(m_inf_col).x - m_f).is_int());  
+
+        mpq one_min_m_f = 1 - m_f;
+        for (const auto & p : m_row) {
+            unsigned j = p.var();
+            if (j == m_inf_col) {
+                lp_assert(p.coeff() == one_of_type<mpq>());
+                TRACE("gomory_cut_detail", tout << "seeing basic var";);
+                continue;
+            }
+
+             // use -p.coeff() to make the format compatible with the format used in: Integrating Simplex with DPLL(T)
+            if (is_real(j)) {  
+                real_case_in_gomory_cut(- p.coeff(), j);
+            } else {
+                if (p.coeff().is_int()) {
+                    // m_fj will be zero and no monomial will be added
+                    continue;
+                }
+                some_int_columns = true;
+                m_fj = fractional_part(-p.coeff());
+				m_one_minus_fj = 1 - m_fj;
+                int_case_in_gomory_cut(j);
+            }
+        }
+
+        if (m_t.is_empty())
+            return report_conflict_from_gomory_cut();
+        if (some_int_columns)
+            adjust_term_and_k_for_some_ints_case_gomory();
+        lp_assert(m_int_solver.current_solution_is_inf_on_cut());
+        TRACE("gomory_cut_detail", dump_cut_and_constraints_as_smt_lemma(tout););
+        m_int_solver.m_lar_solver->subs_term_columns(m_t);
+        TRACE("gomory_cut", print_linear_combination_of_column_indices_only(m_t, tout << "gomory cut:"); 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),
+        m_lcm_den(1),
+        m_f(fractional_part(get_value(basic_inf_int_j).x)),
+        m_one_minus_f(1 - m_f) {}
+    
+};
+
+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..0967c6cc6 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,221 +123,31 @@ 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);
-    mpq sign = *m_upper ? one_of_type<mpq>()  : -one_of_type<mpq>();
-    CTRACE("current_solution_is_inf_on_cut", v * sign <= (*m_k) * sign,
-           tout << "m_upper = " << *m_upper << std::endl;
-           tout << "v = " << v << ", k = " << (*m_k) << std::endl;
+    impq v = m_t.apply(x);
+    mpq sign = m_upper ? one_of_type<mpq>()  : -one_of_type<mpq>();
+    CTRACE("current_solution_is_inf_on_cut", v * sign <= m_k * sign,
+           tout << "m_upper = " << m_upper << std::endl;
+           tout << "v = " << v << ", k = " << m_k << std::endl;
           );
-    return v * sign > (*m_k) * sign;
+    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;
+    m_upper = true;
     return mk_gomory_cut(j, row);
 }
 
@@ -394,17 +172,19 @@ 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);
     }
@@ -589,21 +369,21 @@ lia_move int_solver::make_hnf_cut() {
 #else
     vector<mpq> x0;
 #endif
-    lia_move r =  m_hnf_cutter.create_cut(*m_t, *m_k, *m_ex, *m_upper, x0);
+    lia_move r =  m_hnf_cutter.create_cut(m_t, m_k, m_ex, m_upper, x0);
 
     if (r == lia_move::cut) {      
         TRACE("hnf_cut",
-              m_lar_solver->print_term(*m_t, tout << "cut:"); 
-              tout << " <= " << *m_k << std::endl;
+              m_lar_solver->print_term(m_t, tout << "cut:"); 
+              tout << " <= " << m_k << std::endl;
               for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
                   m_lar_solver->print_constraint(i, tout);
               }              
               );
         lp_assert(current_solution_is_inf_on_cut());
         settings().st().m_hnf_cuts++;
-        m_ex->clear();        
+        m_ex.clear();        
         for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
-             m_ex->push_justification(i);
+             m_ex.push_justification(i);
         }
     } 
     return r;
@@ -619,14 +399,17 @@ lia_move int_solver::hnf_cut() {
     return lia_move::undef;
 }
 
-lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
+lia_move int_solver::check() {
     if (!has_inf_int()) return lia_move::sat;
 
-    m_t = &t;  m_k = &k;  m_ex = &ex; m_upper = &upper;
+    m_t.clear();
+    m_k.reset();
+    m_ex.clear();
+    m_upper = false;
     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();
@@ -862,8 +645,8 @@ bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, uns
 void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j) {
     constraint_index lc, uc;
     m_lar_solver->get_bound_constraint_witnesses_for_column(j, lc, uc);
-    m_ex->m_explanation.push_back(std::make_pair(mpq(1), lc));
-    m_ex->m_explanation.push_back(std::make_pair(mpq(1), uc));
+    m_ex.m_explanation.push_back(std::make_pair(mpq(1), lc));
+    m_ex.m_explanation.push_back(std::make_pair(mpq(1), uc));
 }
 void int_solver::fill_explanation_from_fixed_columns(const row_strip<mpq> & row) {
     for (const auto & c : row) {
@@ -1126,7 +909,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:
@@ -1258,20 +1041,20 @@ const impq& int_solver::lower_bound(unsigned j) const {
 
 lia_move int_solver::create_branch_on_column(int j) {
     TRACE("check_main_int", tout << "branching" << std::endl;);
-    lp_assert(m_t->is_empty());
+    lp_assert(m_t.is_empty());
     lp_assert(j != -1);
-    m_t->add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
+    m_t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
     if (is_free(j)) {
-        *m_upper = true;
-        *m_k = mpq(0);
+        m_upper = true;
+        m_k = mpq(0);
     } else {
-        *m_upper = left_branch_is_more_narrow_than_right(j);
-        *m_k = *m_upper? floor(get_value(j)) : ceil(get_value(j));        
+        m_upper = left_branch_is_more_narrow_than_right(j);
+        m_k = m_upper? floor(get_value(j)) : ceil(get_value(j));        
     }
 
     TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n";
           display_column(tout, j);
-          tout << "k = " << *m_k << std::endl;
+          tout << "k = " << m_k << std::endl;
           );
     return lia_move::branch;
 
diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h
index ec708918d..17ce20481 100644
--- a/src/util/lp/int_solver.h
+++ b/src/util/lp/int_solver.h
@@ -39,20 +39,31 @@ public:
     // fields
     lar_solver          *m_lar_solver;
     unsigned            m_number_of_calls;
-    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
-    bool                *m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise
+    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
+    bool                m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise
     hnf_cutter          m_hnf_cutter;
     // methods
     int_solver(lar_solver* lp);
 
     // main function to check that the solution provided by lar_solver is valid for integral values,
     // or provide a way of how it can be adjusted.
-    lia_move check(lar_term& t, mpq& k, explanation& ex, bool & upper);
+    lia_move check();
+    lar_term const& get_term() const { return m_t; }
+    mpq const& get_offset() const { return m_k; }
+    explanation const& get_explanation() const { return m_ex; }
+    bool is_upper() const { return m_upper; }
+
     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 +90,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 +99,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,39 +111,22 @@ 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;
-    inline static bool is_rational(const impq & n) {
-        return is_zero(n.y);  
-    }
 
 public:
     void display_column(std::ostream & out, unsigned j) const;
-    inline static
-    mpq fractional_part(const impq & n) {
-        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 +151,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_constraints.h b/src/util/lp/lar_constraints.h
index ac15028bb..6305089b4 100644
--- a/src/util/lp/lar_constraints.h
+++ b/src/util/lp/lar_constraints.h
@@ -75,7 +75,7 @@ struct lar_term_constraint: public lar_base_constraint {
     }
     unsigned size() const override { return m_term->size();}
     lar_term_constraint(const lar_term *t, lconstraint_kind kind, const mpq& right_side) : lar_base_constraint(kind, right_side), m_term(t) { }
-    mpq get_free_coeff_of_left_side() const override { return m_term->m_v;}
+    // mpq get_free_coeff_of_left_side() const override { return m_term->m_v;}
 
 };
 
diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp
index 654eb7017..5b3028a98 100644
--- a/src/util/lp/lar_solver.cpp
+++ b/src/util/lp/lar_solver.cpp
@@ -27,7 +27,7 @@ void clear() {lp_assert(false); // not implemented
 }
 
 
-lar_solver::lar_solver() : m_status(lp_status::OPTIMAL),
+lar_solver::lar_solver() : m_status(lp_status::UNKNOWN),
                            m_infeasible_column_index(-1),
                            m_terms_start_index(1000000),
                            m_mpq_lar_core_solver(m_settings, *this),
@@ -71,7 +71,7 @@ bool lar_solver::sizes_are_correct() const {
 }
     
  
-void lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out) const {
+std::ostream& lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out) const {
     out << "implied bound\n";
     unsigned v = be.m_j;
     if (is_term(v)) {
@@ -83,6 +83,7 @@ void lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out
     }
     out << " " << lconstraint_kind_string(be.kind()) << " "  << be.m_bound << std::endl;
     out << "end of implied bound" << std::endl;
+    return out;
 }
     
 bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be, const vector<std::pair<mpq, unsigned>> & explanation) const {
@@ -136,7 +137,7 @@ bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be,
             kind = static_cast<lconstraint_kind>(-kind);
         }
         rs_of_evidence /= ratio;
-        rs_of_evidence += t->m_v * ratio;
+        // rs_of_evidence += t->m_v * ratio;
     }
 
     return kind == be.kind() && rs_of_evidence == be.m_bound;
@@ -601,7 +602,7 @@ void lar_solver::register_monoid_in_map(std::unordered_map<var_index, mpq> & coe
 
 
 void lar_solver::substitute_terms_in_linear_expression(const vector<std::pair<mpq, var_index>>& left_side_with_terms,
-                                                       vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const {
+                                                       vector<std::pair<mpq, var_index>> &left_side) const {
     std::unordered_map<var_index, mpq> coeffs;
     for (auto & t : left_side_with_terms) {
         unsigned j = t.second;
@@ -612,7 +613,6 @@ void lar_solver::substitute_terms_in_linear_expression(const vector<std::pair<mp
             for (auto & p : term.coeffs()){
                 register_monoid_in_map(coeffs, t.first * p.second , p.first);
             }
-            free_coeff += t.first * term.m_v;
         }
     }
 
@@ -781,7 +781,7 @@ void lar_solver::solve_with_core_solver() {
         update_x_and_inf_costs_for_columns_with_changed_bounds();
     m_mpq_lar_core_solver.solve();
     set_status(m_mpq_lar_core_solver.m_r_solver.get_status());
-    lp_assert(m_status != lp_status::OPTIMAL || all_constraints_hold());
+    lp_assert((m_settings.random_next() % 100) != 0 || m_status != lp_status::OPTIMAL || all_constraints_hold());
 }
 
     
@@ -909,13 +909,8 @@ bool lar_solver::try_to_set_fixed(column_info<mpq> & ci) {
     return false;
 }
 
-column_type lar_solver::get_column_type(const column_info<mpq> & ci) {
-    auto ret = ci.get_column_type_no_flipping();
-    if (ret == column_type::boxed) { // changing boxed to fixed because of the no span
-        if (ci.get_lower_bound() == ci.get_upper_bound())
-            ret = column_type::fixed;
-    }
-    return ret;
+column_type lar_solver::get_column_type(unsigned j) const{
+    return m_mpq_lar_core_solver.m_column_types[j];
 }
 
 std::string lar_solver::get_column_name(unsigned j) const {
@@ -954,9 +949,9 @@ bool lar_solver::constraint_holds(const lar_base_constraint & constr, std::unord
     mpq left_side_val = get_left_side_val(constr, var_map);
     switch (constr.m_kind) {
     case LE: return left_side_val <= constr.m_right_side;
-    case LT: return left_side_val < constr.m_right_side;
+    case LT: return left_side_val <  constr.m_right_side;
     case GE: return left_side_val >= constr.m_right_side;
-    case GT: return left_side_val > constr.m_right_side;
+    case GT: return left_side_val >  constr.m_right_side;
     case EQ: return left_side_val == constr.m_right_side;
     default:
         lp_unreachable();
@@ -975,8 +970,10 @@ bool lar_solver::the_relations_are_of_same_type(const vector<std::pair<mpq, unsi
             flip_kind(m_constraints[con_ind]->m_kind);
         if (kind == GT || kind == LT)
             strict = true;
-        if (kind == GE || kind == GT) n_of_G++;
-        else if (kind == LE || kind == LT) n_of_L++;
+        if (kind == GE || kind == GT) 
+            n_of_G++;
+        else if (kind == LE || kind == LT) 
+            n_of_L++;
     }
     the_kind_of_sum = n_of_G ? GE : (n_of_L ? LE : EQ);
     if (strict)
@@ -1116,7 +1113,7 @@ bool lar_solver::has_upper_bound(var_index var, constraint_index& ci, mpq& value
 bool lar_solver::has_value(var_index var, mpq& value) const {
     if (is_term(var)) {
         lar_term const& t = get_term(var);
-        value = t.m_v;
+        value = 0;
         for (auto const& cv : t) {
             impq const& r = get_column_value(cv.var());
             if (!numeric_traits<mpq>::is_zero(r.y)) return false;
@@ -1177,6 +1174,7 @@ void lar_solver::get_model(std::unordered_map<var_index, mpq> & variable_values)
         std::unordered_set<impq> set_of_different_pairs; 
         std::unordered_set<mpq> set_of_different_singles;
         delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(delta);
+        TRACE("get_model", tout << "delta=" << delta << "size = " << m_mpq_lar_core_solver.m_r_x.size() << std::endl;);
         for (i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++ ) {
             const numeric_pair<mpq> & rp = m_mpq_lar_core_solver.m_r_x[i];
             set_of_different_pairs.insert(rp);
@@ -1208,42 +1206,40 @@ std::string lar_solver::get_variable_name(var_index vi) const {
 }
 
 // ********** print region start
-void lar_solver::print_constraint(constraint_index ci, std::ostream & out) const {
+std::ostream& lar_solver::print_constraint(constraint_index ci, std::ostream & out) const {
     if (ci >= m_constraints.size()) {
         out << "constraint " << T_to_string(ci) << " is not found";
         out << std::endl;
-        return;
+        return out;
     }
 
-    print_constraint(m_constraints[ci], out);
+    return print_constraint(m_constraints[ci], out);
 }
 
-void lar_solver::print_constraints(std::ostream& out) const  {
+std::ostream& lar_solver::print_constraints(std::ostream& out) const  {
     out << "number of constraints = " << m_constraints.size() << std::endl;
     for (auto c : m_constraints) {
         print_constraint(c, out);
     }
+    return out;
 }
 
-void lar_solver::print_terms(std::ostream& out) const  {
+std::ostream& lar_solver::print_terms(std::ostream& out) const  {
     for (auto it : m_terms) {
-        print_term(*it, out);
-        out << "\n";
+        print_term(*it, out) << "\n";
     }
+    return out;
 }
 
-void lar_solver::print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const {
+std::ostream& lar_solver::print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const {
     print_linear_combination_of_column_indices(c->get_left_side_coefficients(), out);
     mpq free_coeff = c->get_free_coeff_of_left_side();
     if (!is_zero(free_coeff))
         out << " + " << free_coeff;
-        
+    return out;
 }
 
-void lar_solver::print_term(lar_term const& term, std::ostream & out) const {
-    if (!numeric_traits<mpq>::is_zero(term.m_v)) {
-        out << term.m_v << " + ";
-    }
+std::ostream& lar_solver::print_term(lar_term const& term, std::ostream & out) const {
     bool first = true;
     for (const auto p : term) {
         mpq val = p.coeff();
@@ -1263,14 +1259,12 @@ void lar_solver::print_term(lar_term const& term, std::ostream & out) const {
             out << T_to_string(val);
         out << this->get_column_name(p.var());
     }
-
+    return out;
 }
 
-void lar_solver::print_term_as_indices(lar_term const& term, std::ostream & out) const {
-    if (!numeric_traits<mpq>::is_zero(term.m_v)) {
-        out << term.m_v << " + ";
-    }
-    print_linear_combination_of_column_indices_only(term.coeffs_as_vector(), out);
+std::ostream& lar_solver::print_term_as_indices(lar_term const& term, std::ostream & out) const {
+    print_linear_combination_of_column_indices_only(term, out);
+    return out;
 }
 
 mpq lar_solver::get_left_side_val(const lar_base_constraint &  cns, const std::unordered_map<var_index, mpq> & var_map) const {
@@ -1284,9 +1278,10 @@ mpq lar_solver::get_left_side_val(const lar_base_constraint &  cns, const std::u
     return ret;
 }
 
-void lar_solver::print_constraint(const lar_base_constraint * c, std::ostream & out) const {
+std::ostream& lar_solver::print_constraint(const lar_base_constraint * c, std::ostream & out) const {
     print_left_side_of_constraint(c, out);
     out << " " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side << std::endl;
+    return out;
 }
 
 void lar_solver::fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list) {
@@ -1492,7 +1487,7 @@ bool lar_solver::term_is_int(const lar_term * t) const {
     for (auto const & p :  t->m_coeffs)
         if (! (column_is_int(p.first)  && p.second.is_int()))
             return false;
-    return t->m_v.is_int();
+    return true;
 }
 
 bool lar_solver::var_is_int(var_index v) const {
@@ -1593,17 +1588,13 @@ void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
 }
 
 
-var_index lar_solver::add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
-                                         const mpq &m_v) {
-    push_and_register_term(new lar_term(coeffs, m_v));
+var_index lar_solver::add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs) {
+    push_and_register_term(new lar_term(coeffs));
     return m_terms_start_index + m_terms.size() - 1;
 }
 
 #if Z3DEBUG_CHECK_UNIQUE_TERMS
-bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs, const mpq& v) {
-    if (coeffs.empty()) {
-        return is_zero(v);
-    }
+bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs) {
 
     for (const auto & p : coeffs) {
         if (column_is_real(p.second))
@@ -1638,12 +1629,11 @@ void lar_solver::push_and_register_term(lar_term* t) {
 }
 
 // terms
-var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
-                               const mpq &m_v) {
+var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs) {
     if (strategy_is_undecided())
-        return add_term_undecided(coeffs, m_v);
+        return add_term_undecided(coeffs);
 
-    push_and_register_term(new lar_term(coeffs, m_v));
+    push_and_register_term(new lar_term(coeffs));
     unsigned adjusted_term_index = m_terms.size() - 1;
     var_index ret = m_terms_start_index + adjusted_term_index;
     if (use_tableau() && !coeffs.empty()) {
@@ -1656,7 +1646,7 @@ var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
 }
 
 void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
-    TRACE("dump_terms", print_term(*term, tout); tout << std::endl;);
+    TRACE("dump_terms", print_term(*term, tout) << std::endl;);
     register_new_ext_var_index(term_ext_index, term_is_int(term));
     // j will be a new variable
     unsigned j = A_r().column_count();
@@ -1738,9 +1728,8 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k
     //    lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int());
     unsigned term_j;
     if (m_var_register.external_is_used(j, term_j)) {
-        mpq rs = right_side - m_terms[adjusted_term_index]->m_v;
         m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side));
-        update_column_type_and_bound(term_j, kind, rs, ci);
+        update_column_type_and_bound(term_j, kind, right_side, ci);
     }
     else {
         add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side);
@@ -1749,11 +1738,10 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k
 
 constraint_index lar_solver::add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm) {
     vector<std::pair<mpq, var_index>> left_side;
-    mpq rs = -right_side_parm;
-    substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs);
-    unsigned term_index = add_term(left_side, zero_of_type<mpq>());
+    substitute_terms_in_linear_expression(left_side_with_terms, left_side);
+    unsigned term_index = add_term(left_side);
     constraint_index ci = m_constraints.size();
-    add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci);
+    add_var_bound_on_constraint_for_term(term_index, kind_par, right_side_parm, ci);
     return ci;
 }
 
@@ -1762,7 +1750,7 @@ void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned ter
 
     add_row_from_term_no_constraint(term, term_j);
     unsigned j = A_r().column_count() - 1;
-    update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
+    update_column_type_and_bound(j, kind, right_side, m_constraints.size());
     m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
     lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
 }
@@ -1964,7 +1952,6 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
 
             if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
                 m_status = lp_status::INFEASIBLE;
-                lp_assert(false);
                 m_infeasible_column_index = j;
             }
             else {
@@ -2261,6 +2248,7 @@ void lar_solver::set_cut_strategy(unsigned cut_frequency) {
     } 
 }
 
+
 } // namespace lp
 
 
diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h
index 283c13c38..4189bad4e 100644
--- a/src/util/lp/lar_solver.h
+++ b/src/util/lp/lar_solver.h
@@ -164,13 +164,11 @@ public:
 
 
     // terms
-    var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
-                       const mpq &m_v);
+    var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs);
 
-    var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
-                                 const mpq &m_v);
+    var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs);
 
-    bool term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs, const mpq& v);
+    bool term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs);
     void push_and_register_term(lar_term* t);
 
     void add_row_for_term(const lar_term * term, unsigned term_ext_index);
@@ -208,7 +206,10 @@ public:
     void update_lower_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
 
     void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
+
     //end of init region
+
+
     lp_settings & settings();
 
     lp_settings const & settings() const;
@@ -227,9 +228,7 @@ public:
     bool use_lu() const;
     
     bool sizes_are_correct() const;
- 
-    void print_implied_bound(const implied_bound& be, std::ostream & out) const;
-    
+     
     bool implied_bound_is_correctly_explained(implied_bound const & be, const vector<std::pair<mpq, unsigned>> & explanation) const;
     
     void analyze_new_bounds_on_row(
@@ -238,8 +237,7 @@ public:
 
     void analyze_new_bounds_on_row_tableau(
         unsigned row_index,
-        bound_propagator & bp
-                                           );
+        bound_propagator & bp);
 
     
     void substitute_basis_var_in_terms_for_row(unsigned i);
@@ -331,7 +329,7 @@ public:
 
 
     void substitute_terms_in_linear_expression( const vector<std::pair<mpq, var_index>>& left_side_with_terms,
-                                                vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const;
+                                                vector<std::pair<mpq, var_index>> &left_side) const;
 
 
     void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j);
@@ -397,7 +395,7 @@ public:
 
     bool try_to_set_fixed(column_info<mpq> & ci);
 
-    column_type get_column_type(const column_info<mpq> & ci);
+    column_type get_column_type(unsigned j) const;
 
     std::string get_column_name(unsigned j) const;
 
@@ -436,30 +434,33 @@ public:
         int inf_sign) const;
 
 
-
     void get_model(std::unordered_map<var_index, mpq> & variable_values) const;
 
     void get_model_do_not_care_about_diff_vars(std::unordered_map<var_index, mpq> & variable_values) const;
 
     std::string get_variable_name(var_index vi) const;
 
-    // ********** print region start
-    void print_constraint(constraint_index ci, std::ostream & out) const;
+    // print utilities
 
-    void print_constraints(std::ostream& out) const ;
+    std::ostream& print_constraint(constraint_index ci, std::ostream & out) const;
 
-    void print_terms(std::ostream& out) const;
+    std::ostream& print_constraints(std::ostream& out) const ;
 
-    void print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const;
+    std::ostream& print_terms(std::ostream& out) const;
 
-    void print_term(lar_term const& term, std::ostream & out) const;
+    std::ostream& print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const;
 
-    void print_term_as_indices(lar_term const& term, std::ostream & out) const;
+    std::ostream& print_term(lar_term const& term, std::ostream & out) const;
 
+    std::ostream& print_term_as_indices(lar_term const& term, std::ostream & out) const;
+
+    std::ostream& print_constraint(const lar_base_constraint * c, std::ostream & out) const;
+
+    std::ostream& print_implied_bound(const implied_bound& be, std::ostream & out) const;
+
+    
     mpq get_left_side_val(const lar_base_constraint &  cns, const std::unordered_map<var_index, mpq> & var_map) const;
 
-    void print_constraint(const lar_base_constraint * c, std::ostream & out) const;
-
     void fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list);
 
     void random_update(unsigned sz, var_index const * vars);
diff --git a/src/util/lp/lar_term.h b/src/util/lp/lar_term.h
index 519847848..e9259b8c0 100644
--- a/src/util/lp/lar_term.h
+++ b/src/util/lp/lar_term.h
@@ -21,9 +21,9 @@
 #include "util/lp/indexed_vector.h"
 namespace lp {
 struct lar_term {
-    // the term evaluates to sum of m_coeffs + m_v
+    // the term evaluates to sum of m_coeffs 
     std::unordered_map<unsigned, mpq> m_coeffs;
-    mpq m_v;
+    // mpq m_v;
     lar_term() {}
     void add_monomial(const mpq& c, unsigned j) {
         auto it = m_coeffs.find(j);
@@ -37,7 +37,7 @@ struct lar_term {
     }
 
     bool is_empty() const {
-        return m_coeffs.size() == 0 && is_zero(m_v);
+        return m_coeffs.size() == 0; // && is_zero(m_v);
     }
     
     unsigned size() const { return static_cast<unsigned>(m_coeffs.size()); }
@@ -46,8 +46,7 @@ struct lar_term {
         return m_coeffs;
     }
     
-    lar_term(const vector<std::pair<mpq, unsigned>>& coeffs,
-             const mpq & v) : m_v(v) {
+    lar_term(const vector<std::pair<mpq, unsigned>>& coeffs) {
         for (const auto & p : coeffs) {
             add_monomial(p.first, p.second);
         }
@@ -87,7 +86,7 @@ struct lar_term {
 
     template <typename T>
     T apply(const vector<T>& x) const {
-        T ret = T(m_v);
+        T ret(0);
         for (const auto & t : m_coeffs) {
             ret += t.second * x[t.first];
         }
@@ -96,7 +95,6 @@ struct lar_term {
    
     void clear() {
         m_coeffs.clear();
-        m_v = zero_of_type<mpq>();
     }
 
     struct ival {
diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h
index 41b6fe31d..9a6549917 100644
--- a/src/util/lp/lp_core_solver_base.h
+++ b/src/util/lp/lp_core_solver_base.h
@@ -577,7 +577,7 @@ public:
     }
 
     void print_column_info(unsigned j, std::ostream & out) const {
-        out << "j = " << j << ", name = "<< column_name(j);
+        out << "j = " << j << ",\tname = "<< column_name(j) << "\t";
         switch (m_column_types[j]) {
         case column_type::fixed:
         case column_type::boxed:
@@ -596,11 +596,11 @@ public:
             lp_assert(false);
         }
         //        out << "basis heading = " << m_basis_heading[j] << std::endl;
-        out << " x =                " << m_x[j];
+        out << "\tx = " << m_x[j];
         if (m_basis_heading[j] >= 0)
             out << " base\n";
         else
-            out << " nbas\n";
+           out << " \n";
     }
 
     bool column_is_free(unsigned j) const { return this->m_column_type[j] == free; }
diff --git a/src/util/lp/lp_primal_core_solver_def.h b/src/util/lp/lp_primal_core_solver_def.h
index 1e9edbd31..872922f60 100644
--- a/src/util/lp/lp_primal_core_solver_def.h
+++ b/src/util/lp/lp_primal_core_solver_def.h
@@ -1238,6 +1238,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_column
         break;
     case column_type::free_column:
         out << "( _" << this->m_x[j] << "_)" << std::endl;
+        break;
     default:
         lp_unreachable();
     }
diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h
index dd19df23a..1bbefd154 100644
--- a/src/util/lp/lp_settings.h
+++ b/src/util/lp/lp_settings.h
@@ -357,7 +357,7 @@ public:
     }
     
 #ifdef Z3DEBUG
-    static unsigned ddd; // used for debugging    
+static unsigned ddd; // used for debugging    
 #endif
 }; // end of lp_settings class
 
@@ -433,7 +433,15 @@ inline void ensure_increasing(vector<unsigned> & v) {
         }
 }
 
+inline static bool is_rational(const impq & n) { return is_zero(n.y); }
 
+inline static mpq fractional_part(const impq & n) {
+    lp_assert(is_rational(n));
+    return n.x - floor(n.x);
+}
+inline static mpq fractional_part(const mpq & n) {
+    return n - floor(n);
+}
 
 #if Z3DEBUG
 bool D();
diff --git a/src/util/lp/lp_solver_def.h b/src/util/lp/lp_solver_def.h
index 10c7a6feb..9b385dee6 100644
--- a/src/util/lp/lp_solver_def.h
+++ b/src/util/lp/lp_solver_def.h
@@ -24,7 +24,7 @@ Revision History:
 namespace lp {
 template <typename T, typename X> column_info<T> * lp_solver<T, X>::get_or_create_column_info(unsigned column) {
     auto it = m_map_from_var_index_to_column_info.find(column);
-    return (it == m_map_from_var_index_to_column_info.end())? (m_map_from_var_index_to_column_info[column] = new column_info<T>(static_cast<unsigned>(-1))) : it->second;
+    return (it == m_map_from_var_index_to_column_info.end())? (m_map_from_var_index_to_column_info[column] = new column_info<T>()) : it->second;
 }
 
 template <typename T, typename X>
diff --git a/src/util/lp/lp_utils.h b/src/util/lp/lp_utils.h
index e776092a2..573fc319c 100644
--- a/src/util/lp/lp_utils.h
+++ b/src/util/lp/lp_utils.h
@@ -50,11 +50,34 @@ bool contains(const std::unordered_map<A, B> & map, const A& key) {
 
 namespace lp {
 
-
-   inline void throw_exception(std::string && str) {
-        throw default_exception(std::move(str));
+template <typename T>
+void print_linear_combination_of_column_indices_only(const T & coeffs, std::ostream & out) {
+    bool first = true;
+    for (const auto & it : coeffs) {
+        auto val = it.coeff();
+        if (first) {
+            first = false;
+        } else {
+            if (val.is_pos()) {
+                out << " + ";
+            } else {
+                out << " - ";
+                val = -val;
+            }
+        }
+        if (val == 1)
+            out << " ";
+        else 
+            out << T_to_string(val);
+        
+        out << "x" << it.var();
     }
-    typedef z3_exception exception;
+}
+
+inline void throw_exception(std::string && str) {
+    throw default_exception(std::move(str));
+}
+typedef z3_exception exception;
 
 #define lp_assert(_x_) { SASSERT(_x_); }
 inline void lp_unreachable() { lp_assert(false); }
diff --git a/src/util/mpz.h b/src/util/mpz.h
index 65f89f078..2df03666f 100644
--- a/src/util/mpz.h
+++ b/src/util/mpz.h
@@ -152,7 +152,7 @@ class mpz_manager {
     
     // make sure that n is a big number and has capacity equal to at least c.
     void allocate_if_needed(mpz & n, unsigned c) {
-        c = std::max(c, m_init_cell_capacity);
+        if (m_init_cell_capacity > c) c = m_init_cell_capacity;
         if (n.m_ptr == nullptr || capacity(n) < c) {
             deallocate(n);
             n.m_val             = 1;
diff --git a/src/util/version.h.cmake.in b/src/util/z3_version.h.cmake.in
similarity index 100%
rename from src/util/version.h.cmake.in
rename to src/util/z3_version.h.cmake.in
diff --git a/src/util/version.h.in b/src/util/z3_version.h.in
similarity index 100%
rename from src/util/version.h.in
rename to src/util/z3_version.h.in