Split the non-linear solver (#4219)

This splits the "non-linear solver" from "NonlinearExtension". The non-linear solver is the module that implements the inference schemas whereas NonlinearExtension is the glue code that manages the solver(s) for non-linear.

This also involves moving utilities from the non-linear solver to their own file.

10 files changed, 2675 insertions, 1998 deletions

diff --git a/src/theory/arith/nl_constraint.cpp b/src/theory/arith/nl_constraint.cpp
new file mode 100644
index 000000000..8fb4535ea
--- /dev/null
+++ b/src/theory/arith/nl_constraint.cpp
@@ -0,0 +1,123 @@
+/*********************                                                        */
+/*! \file nl_constraint.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2019 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved.  See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** \brief Implementation of utilities for non-linear constraints
+ **/
+
+#include "theory/arith/nl_constraint.h"
+
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/arith_utilities.h"
+
+using namespace CVC4::kind;
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+
+ConstraintDb::ConstraintDb(MonomialDb& mdb) : d_mdb(mdb) {}
+
+void ConstraintDb::registerConstraint(Node atom)
+{
+  if (std::find(d_constraints.begin(), d_constraints.end(), atom)
+      != d_constraints.end())
+  {
+    return;
+  }
+  d_constraints.push_back(atom);
+  Trace("nl-ext-debug") << "Register constraint : " << atom << std::endl;
+  std::map<Node, Node> msum;
+  if (ArithMSum::getMonomialSumLit(atom, msum))
+  {
+    Trace("nl-ext-debug") << "got monomial sum: " << std::endl;
+    if (Trace.isOn("nl-ext-debug"))
+    {
+      ArithMSum::debugPrintMonomialSum(msum, "nl-ext-debug");
+    }
+    unsigned max_degree = 0;
+    std::vector<Node> all_m;
+    std::vector<Node> max_deg_m;
+    for (std::map<Node, Node>::iterator itm = msum.begin(); itm != msum.end();
+         ++itm)
+    {
+      if (!itm->first.isNull())
+      {
+        all_m.push_back(itm->first);
+        d_mdb.registerMonomial(itm->first);
+        Trace("nl-ext-debug2")
+            << "...process monomial " << itm->first << std::endl;
+        unsigned d = d_mdb.getDegree(itm->first);
+        if (d > max_degree)
+        {
+          max_degree = d;
+          max_deg_m.clear();
+        }
+        if (d >= max_degree)
+        {
+          max_deg_m.push_back(itm->first);
+        }
+      }
+    }
+    // isolate for each maximal degree monomial
+    for (unsigned i = 0; i < all_m.size(); i++)
+    {
+      Node m = all_m[i];
+      Node rhs, coeff;
+      int res = ArithMSum::isolate(m, msum, coeff, rhs, atom.getKind());
+      if (res != 0)
+      {
+        Kind type = atom.getKind();
+        if (res == -1)
+        {
+          type = reverseRelationKind(type);
+        }
+        Trace("nl-ext-constraint") << "Constraint : " << atom << " <=> ";
+        if (!coeff.isNull())
+        {
+          Trace("nl-ext-constraint") << coeff << " * ";
+        }
+        Trace("nl-ext-constraint")
+            << m << " " << type << " " << rhs << std::endl;
+        ConstraintInfo& ci = d_c_info[atom][m];
+        ci.d_rhs = rhs;
+        ci.d_coeff = coeff;
+        ci.d_type = type;
+      }
+    }
+    for (unsigned i = 0; i < max_deg_m.size(); i++)
+    {
+      Node m = max_deg_m[i];
+      d_c_info_maxm[atom][m] = true;
+    }
+  }
+  else
+  {
+    Trace("nl-ext-debug") << "...failed to get monomial sum." << std::endl;
+  }
+}
+
+const std::map<Node, std::map<Node, ConstraintInfo> >&
+ConstraintDb::getConstraints()
+{
+  return d_c_info;
+}
+
+bool ConstraintDb::isMaximal(Node atom, Node x) const
+{
+  std::map<Node, std::map<Node, bool> >::const_iterator itcm =
+      d_c_info_maxm.find(atom);
+  Assert(itcm != d_c_info_maxm.end());
+  return itcm->second.find(x) != itcm->second.end();
+}
+
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl_constraint.h b/src/theory/arith/nl_constraint.h
new file mode 100644
index 000000000..faa3cc812
--- /dev/null
+++ b/src/theory/arith/nl_constraint.h
@@ -0,0 +1,86 @@
+/*********************                                                        */
+/*! \file nl_constraint.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2019 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved.  See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** \brief Utilities for non-linear constraints
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL_CONSTRAINT_H
+#define CVC4__THEORY__ARITH__NL_CONSTRAINT_H
+
+#include <map>
+#include <vector>
+
+#include "expr/kind.h"
+#include "expr/node.h"
+#include "theory/arith/nl_monomial.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+
+/** constraint information
+ *
+ * The struct ( d_rhs, d_coeff, d_type ) represents that a literal is of the
+ * form (d_coeff * x) <d_type> d_rhs.
+ */
+struct ConstraintInfo
+{
+ public:
+  /** The term on the right hand side of the constraint */
+  Node d_rhs;
+  /** The coefficent */
+  Node d_coeff;
+  /** The type (relation) of the constraint */
+  Kind d_type;
+}; /* struct ConstraintInfo */
+
+/** A database for constraints */
+class ConstraintDb
+{
+ public:
+  ConstraintDb(MonomialDb& mdb);
+  ~ConstraintDb() {}
+  /** register constraint
+   *
+   * This ensures that atom is in the domain of the constraints maintained by
+   * this database.
+   */
+  void registerConstraint(Node atom);
+  /** get constraints
+   *
+   * Returns a map m such that whenever
+   * m[lit][x] = ( r, coeff, k ), then
+   * ( lit <=>  (coeff * x) <k> r )
+   */
+  const std::map<Node, std::map<Node, ConstraintInfo> >& getConstraints();
+  /** Returns true if m is of maximal degree in atom
+   *
+   * For example, for atom x^2 + x*y + y >=0, the monomials x^2 and x*y
+   * are of maximal degree (2).
+   */
+  bool isMaximal(Node atom, Node m) const;
+
+ private:
+  /** Reference to a monomial database */
+  MonomialDb& d_mdb;
+  /** List of all constraints */
+  std::vector<Node> d_constraints;
+  /** Is maximal degree */
+  std::map<Node, std::map<Node, bool> > d_c_info_maxm;
+  /** Constraint information */
+  std::map<Node, std::map<Node, ConstraintInfo> > d_c_info;
+};
+
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__NL_SOLVER_H */
diff --git a/src/theory/arith/nl_lemma_utils.h b/src/theory/arith/nl_lemma_utils.h
index 9ad9f2ca5..bd729dad9 100644
--- a/src/theory/arith/nl_lemma_utils.h
+++ b/src/theory/arith/nl_lemma_utils.h
@@ -71,9 +71,9 @@ struct SortNlModel
 
 struct SortNonlinearDegree
 {
-  SortNonlinearDegree(std::map<Node, unsigned>& m) : d_mdegree(m) {}
+  SortNonlinearDegree(const std::map<Node, unsigned>& m) : d_mdegree(m) {}
   /** pointer to the non-linear extension */
-  std::map<Node, unsigned>& d_mdegree;
+  const std::map<Node, unsigned>& d_mdegree;
   /** Get the degree of n in d_mdegree */
   unsigned getDegree(Node n) const;
   /**
diff --git a/src/theory/arith/nl_monomial.cpp b/src/theory/arith/nl_monomial.cpp
new file mode 100644
index 000000000..be472217d
--- /dev/null
+++ b/src/theory/arith/nl_monomial.cpp
@@ -0,0 +1,334 @@
+/*********************                                                        */
+/*! \file nl_monomial.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2019 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved.  See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** \brief Implementation of utilities for monomials
+ **/
+
+#include "theory/arith/nl_monomial.h"
+
+#include "theory/arith/arith_utilities.h"
+#include "theory/arith/nl_lemma_utils.h"
+#include "theory/rewriter.h"
+
+using namespace CVC4::kind;
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+
+// Returns a[key] if key is in a or value otherwise.
+unsigned getCountWithDefault(const NodeMultiset& a, Node key, unsigned value)
+{
+  NodeMultiset::const_iterator it = a.find(key);
+  return (it == a.end()) ? value : it->second;
+}
+// Given two multisets return the multiset difference a \ b.
+NodeMultiset diffMultiset(const NodeMultiset& a, const NodeMultiset& b)
+{
+  NodeMultiset difference;
+  for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a)
+  {
+    Node key = it_a->first;
+    const unsigned a_value = it_a->second;
+    const unsigned b_value = getCountWithDefault(b, key, 0);
+    if (a_value > b_value)
+    {
+      difference[key] = a_value - b_value;
+    }
+  }
+  return difference;
+}
+
+// Return a vector containing a[key] repetitions of key in a multiset a.
+std::vector<Node> ExpandMultiset(const NodeMultiset& a)
+{
+  std::vector<Node> expansion;
+  for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a)
+  {
+    expansion.insert(expansion.end(), it_a->second, it_a->first);
+  }
+  return expansion;
+}
+
+// status 0 : n equal, -1 : n superset, 1 : n subset
+void MonomialIndex::addTerm(Node n,
+                            const std::vector<Node>& reps,
+                            MonomialDb* nla,
+                            int status,
+                            unsigned argIndex)
+{
+  if (status == 0)
+  {
+    if (argIndex == reps.size())
+    {
+      d_monos.push_back(n);
+    }
+    else
+    {
+      d_data[reps[argIndex]].addTerm(n, reps, nla, status, argIndex + 1);
+    }
+  }
+  for (std::map<Node, MonomialIndex>::iterator it = d_data.begin();
+       it != d_data.end();
+       ++it)
+  {
+    if (status != 0 || argIndex == reps.size() || it->first != reps[argIndex])
+    {
+      // if we do not contain this variable, then if we were a superset,
+      // fail (-2), otherwise we are subset.  if we do contain this
+      // variable, then if we were equal, we are superset since variables
+      // are ordered, otherwise we remain the same.
+      int new_status =
+          std::find(reps.begin(), reps.end(), it->first) == reps.end()
+              ? (status >= 0 ? 1 : -2)
+              : (status == 0 ? -1 : status);
+      if (new_status != -2)
+      {
+        it->second.addTerm(n, reps, nla, new_status, argIndex);
+      }
+    }
+  }
+  // compare for subsets
+  for (unsigned i = 0; i < d_monos.size(); i++)
+  {
+    Node m = d_monos[i];
+    if (m != n)
+    {
+      // we are superset if we are equal and haven't traversed all variables
+      int cstatus = status == 0 ? (argIndex == reps.size() ? 0 : -1) : status;
+      Trace("nl-ext-mindex-debug") << "  compare " << n << " and " << m
+                                   << ", status = " << cstatus << std::endl;
+      if (cstatus <= 0 && nla->isMonomialSubset(m, n))
+      {
+        nla->registerMonomialSubset(m, n);
+        Trace("nl-ext-mindex-debug") << "...success" << std::endl;
+      }
+      else if (cstatus >= 0 && nla->isMonomialSubset(n, m))
+      {
+        nla->registerMonomialSubset(n, m);
+        Trace("nl-ext-mindex-debug") << "...success (rev)" << std::endl;
+      }
+    }
+  }
+}
+
+MonomialDb::MonomialDb()
+{
+  d_one = NodeManager::currentNM()->mkConst(Rational(1));
+}
+
+void MonomialDb::registerMonomial(Node n)
+{
+  if (std::find(d_monomials.begin(), d_monomials.end(), n) != d_monomials.end())
+  {
+    return;
+  }
+  d_monomials.push_back(n);
+  Trace("nl-ext-debug") << "Register monomial : " << n << std::endl;
+  Kind k = n.getKind();
+  if (k == NONLINEAR_MULT)
+  {
+    // get exponent count
+    unsigned nchild = n.getNumChildren();
+    for (unsigned i = 0; i < nchild; i++)
+    {
+      d_m_exp[n][n[i]]++;
+      if (i == 0 || n[i] != n[i - 1])
+      {
+        d_m_vlist[n].push_back(n[i]);
+      }
+    }
+    d_m_degree[n] = nchild;
+  }
+  else if (n == d_one)
+  {
+    d_m_exp[n].clear();
+    d_m_vlist[n].clear();
+    d_m_degree[n] = 0;
+  }
+  else
+  {
+    Assert(k != PLUS && k != MULT);
+    d_m_exp[n][n] = 1;
+    d_m_vlist[n].push_back(n);
+    d_m_degree[n] = 1;
+  }
+  std::sort(d_m_vlist[n].begin(), d_m_vlist[n].end());
+  Trace("nl-ext-mindex") << "Add monomial to index : " << n << std::endl;
+  d_m_index.addTerm(n, d_m_vlist[n], this);
+}
+
+void MonomialDb::registerMonomialSubset(Node a, Node b)
+{
+  Assert(isMonomialSubset(a, b));
+
+  const NodeMultiset& a_exponent_map = getMonomialExponentMap(a);
+  const NodeMultiset& b_exponent_map = getMonomialExponentMap(b);
+
+  std::vector<Node> diff_children =
+      ExpandMultiset(diffMultiset(b_exponent_map, a_exponent_map));
+  Assert(!diff_children.empty());
+
+  d_m_contain_parent[a].push_back(b);
+  d_m_contain_children[b].push_back(a);
+
+  Node mult_term = safeConstructNary(MULT, diff_children);
+  Node nlmult_term = safeConstructNary(NONLINEAR_MULT, diff_children);
+  d_m_contain_mult[a][b] = mult_term;
+  d_m_contain_umult[a][b] = nlmult_term;
+  Trace("nl-ext-mindex") << "..." << a << " is a subset of " << b
+                         << ", difference is " << mult_term << std::endl;
+}
+
+bool MonomialDb::isMonomialSubset(Node am, Node bm) const
+{
+  const NodeMultiset& a = getMonomialExponentMap(am);
+  const NodeMultiset& b = getMonomialExponentMap(bm);
+  for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a)
+  {
+    Node key = it_a->first;
+    const unsigned a_value = it_a->second;
+    const unsigned b_value = getCountWithDefault(b, key, 0);
+    if (a_value > b_value)
+    {
+      return false;
+    }
+  }
+  return true;
+}
+
+const NodeMultiset& MonomialDb::getMonomialExponentMap(Node monomial) const
+{
+  MonomialExponentMap::const_iterator it = d_m_exp.find(monomial);
+  Assert(it != d_m_exp.end());
+  return it->second;
+}
+
+unsigned MonomialDb::getExponent(Node monomial, Node v) const
+{
+  MonomialExponentMap::const_iterator it = d_m_exp.find(monomial);
+  if (it == d_m_exp.end())
+  {
+    return 0;
+  }
+  std::map<Node, unsigned>::const_iterator itv = it->second.find(v);
+  if (itv == it->second.end())
+  {
+    return 0;
+  }
+  return itv->second;
+}
+
+const std::vector<Node>& MonomialDb::getVariableList(Node monomial) const
+{
+  std::map<Node, std::vector<Node> >::const_iterator itvl =
+      d_m_vlist.find(monomial);
+  Assert(itvl != d_m_vlist.end());
+  return itvl->second;
+}
+
+unsigned MonomialDb::getDegree(Node monomial) const
+{
+  std::map<Node, unsigned>::const_iterator it = d_m_degree.find(monomial);
+  Assert(it != d_m_degree.end());
+  return it->second;
+}
+
+void MonomialDb::sortByDegree(std::vector<Node>& ms) const
+{
+  SortNonlinearDegree snlad(d_m_degree);
+  std::sort(ms.begin(), ms.end(), snlad);
+}
+
+void MonomialDb::sortVariablesByModel(std::vector<Node>& ms, NlModel& m)
+{
+  SortNlModel smv;
+  smv.d_nlm = &m;
+  smv.d_isConcrete = false;
+  smv.d_isAbsolute = true;
+  smv.d_reverse_order = true;
+  for (const Node& msc : ms)
+  {
+    std::sort(d_m_vlist[msc].begin(), d_m_vlist[msc].end(), smv);
+  }
+}
+
+const std::map<Node, std::vector<Node> >& MonomialDb::getContainsChildrenMap()
+{
+  return d_m_contain_children;
+}
+
+const std::map<Node, std::vector<Node> >& MonomialDb::getContainsParentMap()
+{
+  return d_m_contain_parent;
+}
+
+Node MonomialDb::getContainsDiff(Node a, Node b) const
+{
+  std::map<Node, std::map<Node, Node> >::const_iterator it =
+      d_m_contain_mult.find(a);
+  if (it == d_m_contain_umult.end())
+  {
+    return Node::null();
+  }
+  std::map<Node, Node>::const_iterator it2 = it->second.find(b);
+  if (it2 == it->second.end())
+  {
+    return Node::null();
+  }
+  return it2->second;
+}
+
+Node MonomialDb::getContainsDiffNl(Node a, Node b) const
+{
+  std::map<Node, std::map<Node, Node> >::const_iterator it =
+      d_m_contain_umult.find(a);
+  if (it == d_m_contain_umult.end())
+  {
+    return Node::null();
+  }
+  std::map<Node, Node>::const_iterator it2 = it->second.find(b);
+  if (it2 == it->second.end())
+  {
+    return Node::null();
+  }
+  return it2->second;
+}
+
+Node MonomialDb::mkMonomialRemFactor(Node n,
+                                     const NodeMultiset& n_exp_rem) const
+{
+  std::vector<Node> children;
+  const NodeMultiset& exponent_map = getMonomialExponentMap(n);
+  for (NodeMultiset::const_iterator itme2 = exponent_map.begin();
+       itme2 != exponent_map.end();
+       ++itme2)
+  {
+    Node v = itme2->first;
+    unsigned inc = itme2->second;
+    Trace("nl-ext-mono-factor")
+        << "..." << inc << " factors of " << v << std::endl;
+    unsigned count_in_n_exp_rem = getCountWithDefault(n_exp_rem, v, 0);
+    Assert(count_in_n_exp_rem <= inc);
+    inc -= count_in_n_exp_rem;
+    Trace("nl-ext-mono-factor")
+        << "......rem, now " << inc << " factors of " << v << std::endl;
+    children.insert(children.end(), inc, v);
+  }
+  Node ret = safeConstructNary(MULT, children);
+  ret = Rewriter::rewrite(ret);
+  Trace("nl-ext-mono-factor") << "...return : " << ret << std::endl;
+  return ret;
+}
+
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl_monomial.h b/src/theory/arith/nl_monomial.h
new file mode 100644
index 000000000..81665c4d9
--- /dev/null
+++ b/src/theory/arith/nl_monomial.h
@@ -0,0 +1,147 @@
+/*********************                                                        */
+/*! \file nl_monomial.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2019 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved.  See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** \brief Utilities for monomials
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL_MONOMIAL_H
+#define CVC4__THEORY__ARITH__NL_MONOMIAL_H
+
+#include <map>
+#include <vector>
+
+#include "expr/node.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+
+class MonomialDb;
+class NlModel;
+
+typedef std::map<Node, unsigned> NodeMultiset;
+typedef std::map<Node, NodeMultiset> MonomialExponentMap;
+
+/** An index data structure for node multisets (monomials) */
+class MonomialIndex
+{
+ public:
+  /**
+   * Add term to this trie. The argument status indicates what the status
+   * of n is with respect to the current node in the trie, where:
+   *   0 : n is equal, -1 : n is superset, 1 : n is subset
+   * of the node described by the current path in the trie.
+   */
+  void addTerm(Node n,
+               const std::vector<Node>& reps,
+               MonomialDb* nla,
+               int status = 0,
+               unsigned argIndex = 0);
+
+ private:
+  /** The children of this node */
+  std::map<Node, MonomialIndex> d_data;
+  /** The monomials at this node */
+  std::vector<Node> d_monos;
+}; /* class MonomialIndex */
+
+/** Context-independent database for monomial information */
+class MonomialDb
+{
+ public:
+  MonomialDb();
+  ~MonomialDb() {}
+  /** register monomial */
+  void registerMonomial(Node n);
+  /**
+   * Register monomial subset. This method is called when we infer that b is
+   * a subset of monomial a, e.g. x*y^2 is a subset of x^3*y^2*z.
+   */
+  void registerMonomialSubset(Node a, Node b);
+  /**
+   * returns true if the multiset containing the
+   * factors of monomial a is a subset of the multiset
+   * containing the factors of monomial b.
+   */
+  bool isMonomialSubset(Node a, Node b) const;
+  /** Returns the NodeMultiset for a registered monomial. */
+  const NodeMultiset& getMonomialExponentMap(Node monomial) const;
+  /** Returns the exponent of variable v in the given monomial */
+  unsigned getExponent(Node monomial, Node v) const;
+  /** Get the list of unique variables is the monomial */
+  const std::vector<Node>& getVariableList(Node monomial) const;
+  /** Get degree of monomial, e.g. the degree of x^2*y^2 = 4 */
+  unsigned getDegree(Node monomial) const;
+  /** Sort monomials in ms by their degree
+   *
+   * Updates ms so that degree(ms[i]) <= degree(ms[j]) for i <= j.
+   */
+  void sortByDegree(std::vector<Node>& ms) const;
+  /** Sort the variable lists based on model values
+   *
+   * This updates the variable lists of monomials in ms based on the absolute
+   * value of their current model values in m.
+   *
+   * In other words, for each i, getVariableList(ms[i]) returns
+   * v1, ..., vn where |m(v1)| <= ... <= |m(vn)| after this method is invoked.
+   */
+  void sortVariablesByModel(std::vector<Node>& ms, NlModel& m);
+  /** Get monomial contains children map
+   *
+   * This maps monomials to other monomials that are contained in them, e.g.
+   * x^2 * y may map to { x, x^2, y } if these three terms exist have been
+   * registered to this class.
+   */
+  const std::map<Node, std::vector<Node> >& getContainsChildrenMap();
+  /** Get monomial contains parent map, reverse of above */
+  const std::map<Node, std::vector<Node> >& getContainsParentMap();
+  /**
+   * Get contains difference. Return the difference of a and b or null if it
+   * does not exist. In other words, this returns a term equivalent to a/b
+   * that does not contain division.
+   */
+  Node getContainsDiff(Node a, Node b) const;
+  /**
+   * Get contains difference non-linear. Same as above, but stores terms of kind
+   * NONLINEAR_MULT instead of MULT.
+   */
+  Node getContainsDiffNl(Node a, Node b) const;
+  /** Make monomial remainder factor */
+  Node mkMonomialRemFactor(Node n, const NodeMultiset& n_exp_rem) const;
+
+ private:
+  /** commonly used terms */
+  Node d_one;
+  /** list of all monomials */
+  std::vector<Node> d_monomials;
+  /** Map from monomials to var^index. */
+  MonomialExponentMap d_m_exp;
+  /**
+   * Mapping from monomials to the list of variables that occur in it. For
+   * example, x*x*y*z -> { x, y, z }.
+   */
+  std::map<Node, std::vector<Node> > d_m_vlist;
+  /** Degree information */
+  std::map<Node, unsigned> d_m_degree;
+  /** monomial index, by sorted variables */
+  MonomialIndex d_m_index;
+  /** containment ordering */
+  std::map<Node, std::vector<Node> > d_m_contain_children;
+  std::map<Node, std::vector<Node> > d_m_contain_parent;
+  std::map<Node, std::map<Node, Node> > d_m_contain_mult;
+  std::map<Node, std::map<Node, Node> > d_m_contain_umult;
+};
+
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__NL_MONOMIAL_H */
diff --git a/src/theory/arith/nl_solver.cpp b/src/theory/arith/nl_solver.cpp
new file mode 100644
index 000000000..893d3dbd7
--- /dev/null
+++ b/src/theory/arith/nl_solver.cpp
@@ -0,0 +1,1585 @@
+/*********************                                                        */
+/*! \file nl_solver.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2019 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved.  See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** \brief Implementation of non-linear solver
+ **/
+
+#include "theory/arith/nl_solver.h"
+
+#include "options/arith_options.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/arith_utilities.h"
+#include "theory/arith/theory_arith.h"
+#include "theory/theory_model.h"
+
+using namespace CVC4::kind;
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+
+void debugPrintBound(const char* c, Node coeff, Node x, Kind type, Node rhs)
+{
+  Node t = ArithMSum::mkCoeffTerm(coeff, x);
+  Trace(c) << t << " " << type << " " << rhs;
+}
+
+bool hasNewMonomials(Node n, const std::vector<Node>& existing)
+{
+  std::set<Node> visited;
+
+  std::vector<Node> worklist;
+  worklist.push_back(n);
+  while (!worklist.empty())
+  {
+    Node current = worklist.back();
+    worklist.pop_back();
+    if (visited.find(current) == visited.end())
+    {
+      visited.insert(current);
+      if (current.getKind() == NONLINEAR_MULT)
+      {
+        if (std::find(existing.begin(), existing.end(), current)
+            == existing.end())
+        {
+          return true;
+        }
+      }
+      else
+      {
+        worklist.insert(worklist.end(), current.begin(), current.end());
+      }
+    }
+  }
+  return false;
+}
+
+NlSolver::NlSolver(TheoryArith& containing, NlModel& model)
+    : d_containing(containing),
+      d_model(model),
+      d_cdb(d_mdb),
+      d_zero_split(containing.getUserContext())
+{
+  NodeManager* nm = NodeManager::currentNM();
+  d_true = nm->mkConst(true);
+  d_false = nm->mkConst(false);
+  d_zero = nm->mkConst(Rational(0));
+  d_one = nm->mkConst(Rational(1));
+  d_neg_one = nm->mkConst(Rational(-1));
+  d_order_points.push_back(d_neg_one);
+  d_order_points.push_back(d_zero);
+  d_order_points.push_back(d_one);
+}
+
+NlSolver::~NlSolver() {}
+
+void NlSolver::initLastCall(const std::vector<Node>& assertions,
+                            const std::vector<Node>& false_asserts,
+                            const std::vector<Node>& xts)
+{
+  d_ms_vars.clear();
+  d_ms_proc.clear();
+  d_ms.clear();
+  d_mterms.clear();
+  d_m_nconst_factor.clear();
+  d_tplane_refine.clear();
+  d_ci.clear();
+  d_ci_exp.clear();
+  d_ci_max.clear();
+
+  Trace("nl-ext-mv") << "Extended terms : " << std::endl;
+  // for computing congruence
+  std::map<Kind, ArgTrie> argTrie;
+  for (unsigned i = 0, xsize = xts.size(); i < xsize; i++)
+  {
+    Node a = xts[i];
+    d_model.computeConcreteModelValue(a);
+    d_model.computeAbstractModelValue(a);
+    d_model.printModelValue("nl-ext-mv", a);
+    Kind ak = a.getKind();
+    if (ak == NONLINEAR_MULT)
+    {
+      d_ms.push_back(a);
+
+      // context-independent registration
+      d_mdb.registerMonomial(a);
+
+      const std::vector<Node>& varList = d_mdb.getVariableList(a);
+      for (const Node& v : varList)
+      {
+        if (std::find(d_ms_vars.begin(), d_ms_vars.end(), v) == d_ms_vars.end())
+        {
+          d_ms_vars.push_back(v);
+        }
+        Node mvk = d_model.computeAbstractModelValue(v);
+        if (!mvk.isConst())
+        {
+          d_m_nconst_factor[a] = true;
+        }
+      }
+      // mark processed if has a "one" factor (will look at reduced monomial)?
+    }
+  }
+
+  // register constants
+  d_mdb.registerMonomial(d_one);
+  for (unsigned j = 0; j < d_order_points.size(); j++)
+  {
+    Node c = d_order_points[j];
+    d_model.computeConcreteModelValue(c);
+    d_model.computeAbstractModelValue(c);
+  }
+
+  // register variables
+  Trace("nl-ext-mv") << "Variables in monomials : " << std::endl;
+  for (unsigned i = 0; i < d_ms_vars.size(); i++)
+  {
+    Node v = d_ms_vars[i];
+    d_mdb.registerMonomial(v);
+    d_model.computeConcreteModelValue(v);
+    d_model.computeAbstractModelValue(v);
+    d_model.printModelValue("nl-ext-mv", v);
+  }
+
+  Trace("nl-ext") << "We have " << d_ms.size() << " monomials." << std::endl;
+}
+
+void NlSolver::setMonomialFactor(Node a, Node b, const NodeMultiset& common)
+{
+  // Could not tell if this was being inserted intentionally or not.
+  std::map<Node, Node>& mono_diff_a = d_mono_diff[a];
+  if (mono_diff_a.find(b) == mono_diff_a.end())
+  {
+    Trace("nl-ext-mono-factor")
+        << "Set monomial factor for " << a << "/" << b << std::endl;
+    mono_diff_a[b] = d_mdb.mkMonomialRemFactor(a, common);
+  }
+}
+
+std::vector<Node> NlSolver::checkSplitZero()
+{
+  std::vector<Node> lemmas;
+  for (unsigned i = 0; i < d_ms_vars.size(); i++)
+  {
+    Node v = d_ms_vars[i];
+    if (d_zero_split.insert(v))
+    {
+      Node eq = v.eqNode(d_zero);
+      eq = Rewriter::rewrite(eq);
+      Node literal = d_containing.getValuation().ensureLiteral(eq);
+      d_containing.getOutputChannel().requirePhase(literal, true);
+      lemmas.push_back(literal.orNode(literal.negate()));
+    }
+  }
+  return lemmas;
+}
+
+void NlSolver::assignOrderIds(std::vector<Node>& vars,
+                              NodeMultiset& order,
+                              bool isConcrete,
+                              bool isAbsolute)
+{
+  SortNlModel smv;
+  smv.d_nlm = &d_model;
+  smv.d_isConcrete = isConcrete;
+  smv.d_isAbsolute = isAbsolute;
+  smv.d_reverse_order = false;
+  std::sort(vars.begin(), vars.end(), smv);
+
+  order.clear();
+  // assign ordering id's
+  unsigned counter = 0;
+  unsigned order_index = isConcrete ? 0 : 1;
+  Node prev;
+  for (unsigned j = 0; j < vars.size(); j++)
+  {
+    Node x = vars[j];
+    Node v = d_model.computeModelValue(x, isConcrete);
+    if (!v.isConst())
+    {
+      Trace("nl-ext-mvo") << "..do not assign order to " << x << " : " << v
+                          << std::endl;
+      // don't assign for non-constant values (transcendental function apps)
+      break;
+    }
+    Trace("nl-ext-mvo") << "  order " << x << " : " << v << std::endl;
+    if (v != prev)
+    {
+      // builtin points
+      bool success;
+      do
+      {
+        success = false;
+        if (order_index < d_order_points.size())
+        {
+          Node vv = d_model.computeModelValue(d_order_points[order_index],
+                                              isConcrete);
+          if (d_model.compareValue(v, vv, isAbsolute) <= 0)
+          {
+            counter++;
+            Trace("nl-ext-mvo") << "O[" << d_order_points[order_index]
+                                << "] = " << counter << std::endl;
+            order[d_order_points[order_index]] = counter;
+            prev = vv;
+            order_index++;
+            success = true;
+          }
+        }
+      } while (success);
+    }
+    if (prev.isNull() || d_model.compareValue(v, prev, isAbsolute) != 0)
+    {
+      counter++;
+    }
+    Trace("nl-ext-mvo") << "O[" << x << "] = " << counter << std::endl;
+    order[x] = counter;
+    prev = v;
+  }
+  while (order_index < d_order_points.size())
+  {
+    counter++;
+    Trace("nl-ext-mvo") << "O[" << d_order_points[order_index]
+                        << "] = " << counter << std::endl;
+    order[d_order_points[order_index]] = counter;
+    order_index++;
+  }
+}
+
+// show a <> 0 by inequalities between variables in monomial a w.r.t 0
+int NlSolver::compareSign(Node oa,
+                          Node a,
+                          unsigned a_index,
+                          int status,
+                          std::vector<Node>& exp,
+                          std::vector<Node>& lem)
+{
+  Trace("nl-ext-debug") << "Process " << a << " at index " << a_index
+                        << ", status is " << status << std::endl;
+  NodeManager* nm = NodeManager::currentNM();
+  Node mvaoa = d_model.computeAbstractModelValue(oa);
+  const std::vector<Node>& vla = d_mdb.getVariableList(a);
+  if (a_index == vla.size())
+  {
+    if (mvaoa.getConst<Rational>().sgn() != status)
+    {
+      Node lemma =
+          safeConstructNary(AND, exp).impNode(mkLit(oa, d_zero, status * 2));
+      lem.push_back(lemma);
+    }
+    return status;
+  }
+  Assert(a_index < vla.size());
+  Node av = vla[a_index];
+  unsigned aexp = d_mdb.getExponent(a, av);
+  // take current sign in model
+  Node mvaav = d_model.computeAbstractModelValue(av);
+  int sgn = mvaav.getConst<Rational>().sgn();
+  Trace("nl-ext-debug") << "Process var " << av << "^" << aexp
+                        << ", model sign = " << sgn << std::endl;
+  if (sgn == 0)
+  {
+    if (mvaoa.getConst<Rational>().sgn() != 0)
+    {
+      Node lemma = av.eqNode(d_zero).impNode(oa.eqNode(d_zero));
+      lem.push_back(lemma);
+    }
+    return 0;
+  }
+  if (aexp % 2 == 0)
+  {
+    exp.push_back(av.eqNode(d_zero).negate());
+    return compareSign(oa, a, a_index + 1, status, exp, lem);
+  }
+  exp.push_back(nm->mkNode(sgn == 1 ? GT : LT, av, d_zero));
+  return compareSign(oa, a, a_index + 1, status * sgn, exp, lem);
+}
+
+bool NlSolver::compareMonomial(
+    Node oa,
+    Node a,
+    NodeMultiset& a_exp_proc,
+    Node ob,
+    Node b,
+    NodeMultiset& b_exp_proc,
+    std::vector<Node>& exp,
+    std::vector<Node>& lem,
+    std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers)
+{
+  Trace("nl-ext-comp-debug")
+      << "Check |" << a << "| >= |" << b << "|" << std::endl;
+  unsigned pexp_size = exp.size();
+  if (compareMonomial(
+          oa, a, 0, a_exp_proc, ob, b, 0, b_exp_proc, 0, exp, lem, cmp_infers))
+  {
+    return true;
+  }
+  exp.resize(pexp_size);
+  Trace("nl-ext-comp-debug")
+      << "Check |" << b << "| >= |" << a << "|" << std::endl;
+  if (compareMonomial(
+          ob, b, 0, b_exp_proc, oa, a, 0, a_exp_proc, 0, exp, lem, cmp_infers))
+  {
+    return true;
+  }
+  return false;
+}
+
+Node NlSolver::mkLit(Node a, Node b, int status, bool isAbsolute)
+{
+  if (status == 0)
+  {
+    Node a_eq_b = a.eqNode(b);
+    if (!isAbsolute)
+    {
+      return a_eq_b;
+    }
+    Node negate_b = NodeManager::currentNM()->mkNode(UMINUS, b);
+    return a_eq_b.orNode(a.eqNode(negate_b));
+  }
+  else if (status < 0)
+  {
+    return mkLit(b, a, -status);
+  }
+  Assert(status == 1 || status == 2);
+  NodeManager* nm = NodeManager::currentNM();
+  Kind greater_op = status == 1 ? GEQ : GT;
+  if (!isAbsolute)
+  {
+    return nm->mkNode(greater_op, a, b);
+  }
+  // return nm->mkNode( greater_op, mkAbs( a ), mkAbs( b ) );
+  Node zero = mkRationalNode(0);
+  Node a_is_nonnegative = nm->mkNode(GEQ, a, zero);
+  Node b_is_nonnegative = nm->mkNode(GEQ, b, zero);
+  Node negate_a = nm->mkNode(UMINUS, a);
+  Node negate_b = nm->mkNode(UMINUS, b);
+  return a_is_nonnegative.iteNode(
+      b_is_nonnegative.iteNode(nm->mkNode(greater_op, a, b),
+                               nm->mkNode(greater_op, a, negate_b)),
+      b_is_nonnegative.iteNode(nm->mkNode(greater_op, negate_a, b),
+                               nm->mkNode(greater_op, negate_a, negate_b)));
+}
+
+bool NlSolver::cmp_holds(Node x,
+                         Node y,
+                         std::map<Node, std::map<Node, Node> >& cmp_infers,
+                         std::vector<Node>& exp,
+                         std::map<Node, bool>& visited)
+{
+  if (x == y)
+  {
+    return true;
+  }
+  else if (visited.find(x) != visited.end())
+  {
+    return false;
+  }
+  visited[x] = true;
+  std::map<Node, std::map<Node, Node> >::iterator it = cmp_infers.find(x);
+  if (it != cmp_infers.end())
+  {
+    for (std::map<Node, Node>::iterator itc = it->second.begin();
+         itc != it->second.end();
+         ++itc)
+    {
+      exp.push_back(itc->second);
+      if (cmp_holds(itc->first, y, cmp_infers, exp, visited))
+      {
+        return true;
+      }
+      exp.pop_back();
+    }
+  }
+  return false;
+}
+
+// trying to show a ( >, = ) b   by inequalities between variables in
+// monomials a,b
+bool NlSolver::compareMonomial(
+    Node oa,
+    Node a,
+    unsigned a_index,
+    NodeMultiset& a_exp_proc,
+    Node ob,
+    Node b,
+    unsigned b_index,
+    NodeMultiset& b_exp_proc,
+    int status,
+    std::vector<Node>& exp,
+    std::vector<Node>& lem,
+    std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers)
+{
+  Trace("nl-ext-comp-debug")
+      << "compareMonomial " << oa << " and " << ob << ", indices = " << a_index
+      << " " << b_index << std::endl;
+  Assert(status == 0 || status == 2);
+  const std::vector<Node>& vla = d_mdb.getVariableList(a);
+  const std::vector<Node>& vlb = d_mdb.getVariableList(b);
+  if (a_index == vla.size() && b_index == vlb.size())
+  {
+    // finished, compare absolute value of abstract model values
+    int modelStatus = d_model.compare(oa, ob, false, true) * -2;
+    Trace("nl-ext-comp") << "...finished comparison with " << oa << " <"
+                         << status << "> " << ob
+                         << ", model status = " << modelStatus << std::endl;
+    if (status != modelStatus)
+    {
+      Trace("nl-ext-comp-infer")
+          << "infer : " << oa << " <" << status << "> " << ob << std::endl;
+      if (status == 2)
+      {
+        // must state that all variables are non-zero
+        for (unsigned j = 0; j < vla.size(); j++)
+        {
+          exp.push_back(vla[j].eqNode(d_zero).negate());
+        }
+      }
+      NodeManager* nm = NodeManager::currentNM();
+      Node clem = nm->mkNode(
+          IMPLIES, safeConstructNary(AND, exp), mkLit(oa, ob, status, true));
+      Trace("nl-ext-comp-lemma") << "comparison lemma : " << clem << std::endl;
+      lem.push_back(clem);
+      cmp_infers[status][oa][ob] = clem;
+    }
+    return true;
+  }
+  // get a/b variable information
+  Node av;
+  unsigned aexp = 0;
+  unsigned avo = 0;
+  if (a_index < vla.size())
+  {
+    av = vla[a_index];
+    unsigned aexpTotal = d_mdb.getExponent(a, av);
+    Assert(a_exp_proc[av] <= aexpTotal);
+    aexp = aexpTotal - a_exp_proc[av];
+    if (aexp == 0)
+    {
+      return compareMonomial(oa,
+                             a,
+                             a_index + 1,
+                             a_exp_proc,
+                             ob,
+                             b,
+                             b_index,
+                             b_exp_proc,
+                             status,
+                             exp,
+                             lem,
+                             cmp_infers);
+    }
+    Assert(d_order_vars.find(av) != d_order_vars.end());
+    avo = d_order_vars[av];
+  }
+  Node bv;
+  unsigned bexp = 0;
+  unsigned bvo = 0;
+  if (b_index < vlb.size())
+  {
+    bv = vlb[b_index];
+    unsigned bexpTotal = d_mdb.getExponent(b, bv);
+    Assert(b_exp_proc[bv] <= bexpTotal);
+    bexp = bexpTotal - b_exp_proc[bv];
+    if (bexp == 0)
+    {
+      return compareMonomial(oa,
+                             a,
+                             a_index,
+                             a_exp_proc,
+                             ob,
+                             b,
+                             b_index + 1,
+                             b_exp_proc,
+                             status,
+                             exp,
+                             lem,
+                             cmp_infers);
+    }
+    Assert(d_order_vars.find(bv) != d_order_vars.end());
+    bvo = d_order_vars[bv];
+  }
+  // get "one" information
+  Assert(d_order_vars.find(d_one) != d_order_vars.end());
+  unsigned ovo = d_order_vars[d_one];
+  Trace("nl-ext-comp-debug") << "....vars : " << av << "^" << aexp << " " << bv
+                             << "^" << bexp << std::endl;
+
+  //--- cases
+  if (av.isNull())
+  {
+    if (bvo <= ovo)
+    {
+      Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl;
+      // can multiply b by <=1
+      exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true));
+      return compareMonomial(oa,
+                             a,
+                             a_index,
+                             a_exp_proc,
+                             ob,
+                             b,
+                             b_index + 1,
+                             b_exp_proc,
+                             bvo == ovo ? status : 2,
+                             exp,
+                             lem,
+                             cmp_infers);
+    }
+    Trace("nl-ext-comp-debug")
+        << "...failure, unmatched |b|>1 component." << std::endl;
+    return false;
+  }
+  else if (bv.isNull())
+  {
+    if (avo >= ovo)
+    {
+      Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl;
+      // can multiply a by >=1
+      exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true));
+      return compareMonomial(oa,
+                             a,
+                             a_index + 1,
+                             a_exp_proc,
+                             ob,
+                             b,
+                             b_index,
+                             b_exp_proc,
+                             avo == ovo ? status : 2,
+                             exp,
+                             lem,
+                             cmp_infers);
+    }
+    Trace("nl-ext-comp-debug")
+        << "...failure, unmatched |a|<1 component." << std::endl;
+    return false;
+  }
+  Assert(!av.isNull() && !bv.isNull());
+  if (avo >= bvo)
+  {
+    if (bvo < ovo && avo >= ovo)
+    {
+      Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl;
+      // do avo>=1 instead
+      exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true));
+      return compareMonomial(oa,
+                             a,
+                             a_index + 1,
+                             a_exp_proc,
+                             ob,
+                             b,
+                             b_index,
+                             b_exp_proc,
+                             avo == ovo ? status : 2,
+                             exp,
+                             lem,
+                             cmp_infers);
+    }
+    unsigned min_exp = aexp > bexp ? bexp : aexp;
+    a_exp_proc[av] += min_exp;
+    b_exp_proc[bv] += min_exp;
+    Trace("nl-ext-comp-debug") << "...take leading " << min_exp << " from "
+                               << av << " and " << bv << std::endl;
+    exp.push_back(mkLit(av, bv, avo == bvo ? 0 : 2, true));
+    bool ret = compareMonomial(oa,
+                               a,
+                               a_index,
+                               a_exp_proc,
+                               ob,
+                               b,
+                               b_index,
+                               b_exp_proc,
+                               avo == bvo ? status : 2,
+                               exp,
+                               lem,
+                               cmp_infers);
+    a_exp_proc[av] -= min_exp;
+    b_exp_proc[bv] -= min_exp;
+    return ret;
+  }
+  if (bvo <= ovo)
+  {
+    Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl;
+    // try multiply b <= 1
+    exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true));
+    return compareMonomial(oa,
+                           a,
+                           a_index,
+                           a_exp_proc,
+                           ob,
+                           b,
+                           b_index + 1,
+                           b_exp_proc,
+                           bvo == ovo ? status : 2,
+                           exp,
+                           lem,
+                           cmp_infers);
+  }
+  Trace("nl-ext-comp-debug")
+      << "...failure, leading |b|>|a|>1 component." << std::endl;
+  return false;
+}
+
+std::vector<Node> NlSolver::checkMonomialSign()
+{
+  std::vector<Node> lemmas;
+  std::map<Node, int> signs;
+  Trace("nl-ext") << "Get monomial sign lemmas..." << std::endl;
+  for (unsigned j = 0; j < d_ms.size(); j++)
+  {
+    Node a = d_ms[j];
+    if (d_ms_proc.find(a) == d_ms_proc.end())
+    {
+      std::vector<Node> exp;
+      if (Trace.isOn("nl-ext-debug"))
+      {
+        Node cmva = d_model.computeConcreteModelValue(a);
+        Trace("nl-ext-debug")
+            << "  process " << a << ", mv=" << cmva << "..." << std::endl;
+      }
+      if (d_m_nconst_factor.find(a) == d_m_nconst_factor.end())
+      {
+        signs[a] = compareSign(a, a, 0, 1, exp, lemmas);
+        if (signs[a] == 0)
+        {
+          d_ms_proc[a] = true;
+          Trace("nl-ext-debug")
+              << "...mark " << a << " reduced since its value is 0."
+              << std::endl;
+        }
+      }
+      else
+      {
+        Trace("nl-ext-debug")
+            << "...can't conclude sign lemma for " << a
+            << " since model value of a factor is non-constant." << std::endl;
+      }
+    }
+  }
+  return lemmas;
+}
+
+std::vector<Node> NlSolver::checkMonomialMagnitude(unsigned c)
+{
+  // ensure information is setup
+  if (c == 0)
+  {
+    // sort by absolute values of abstract model values
+    assignOrderIds(d_ms_vars, d_order_vars, false, true);
+
+    // sort individual variable lists
+    Trace("nl-ext-proc") << "Assign order var lists..." << std::endl;
+    d_mdb.sortVariablesByModel(d_ms, d_model);
+  }
+
+  unsigned r = 1;
+  std::vector<Node> lemmas;
+  // if (x,y,L) in cmp_infers, then x > y inferred as conclusion of L
+  // in lemmas
+  std::map<int, std::map<Node, std::map<Node, Node> > > cmp_infers;
+  Trace("nl-ext") << "Get monomial comparison lemmas (order=" << r
+                  << ", compare=" << c << ")..." << std::endl;
+  for (unsigned j = 0; j < d_ms.size(); j++)
+  {
+    Node a = d_ms[j];
+    if (d_ms_proc.find(a) == d_ms_proc.end()
+        && d_m_nconst_factor.find(a) == d_m_nconst_factor.end())
+    {
+      if (c == 0)
+      {
+        // compare magnitude against 1
+        std::vector<Node> exp;
+        NodeMultiset a_exp_proc;
+        NodeMultiset b_exp_proc;
+        compareMonomial(a,
+                        a,
+                        a_exp_proc,
+                        d_one,
+                        d_one,
+                        b_exp_proc,
+                        exp,
+                        lemmas,
+                        cmp_infers);
+      }
+      else
+      {
+        const NodeMultiset& mea = d_mdb.getMonomialExponentMap(a);
+        if (c == 1)
+        {
+          // could compare not just against containing variables?
+          // compare magnitude against variables
+          for (unsigned k = 0; k < d_ms_vars.size(); k++)
+          {
+            Node v = d_ms_vars[k];
+            Node mvcv = d_model.computeConcreteModelValue(v);
+            if (mvcv.isConst())
+            {
+              std::vector<Node> exp;
+              NodeMultiset a_exp_proc;
+              NodeMultiset b_exp_proc;
+              if (mea.find(v) != mea.end())
+              {
+                a_exp_proc[v] = 1;
+                b_exp_proc[v] = 1;
+                setMonomialFactor(a, v, a_exp_proc);
+                setMonomialFactor(v, a, b_exp_proc);
+                compareMonomial(a,
+                                a,
+                                a_exp_proc,
+                                v,
+                                v,
+                                b_exp_proc,
+                                exp,
+                                lemmas,
+                                cmp_infers);
+              }
+            }
+          }
+        }
+        else
+        {
+          // compare magnitude against other non-linear monomials
+          for (unsigned k = (j + 1); k < d_ms.size(); k++)
+          {
+            Node b = d_ms[k];
+            //(signs[a]==signs[b])==(r==0)
+            if (d_ms_proc.find(b) == d_ms_proc.end()
+                && d_m_nconst_factor.find(b) == d_m_nconst_factor.end())
+            {
+              const NodeMultiset& meb = d_mdb.getMonomialExponentMap(b);
+
+              std::vector<Node> exp;
+              // take common factors of monomials, set minimum of
+              // common exponents as processed
+              NodeMultiset a_exp_proc;
+              NodeMultiset b_exp_proc;
+              for (NodeMultiset::const_iterator itmea2 = mea.begin();
+                   itmea2 != mea.end();
+                   ++itmea2)
+              {
+                NodeMultiset::const_iterator itmeb2 = meb.find(itmea2->first);
+                if (itmeb2 != meb.end())
+                {
+                  unsigned min_exp = itmea2->second > itmeb2->second
+                                         ? itmeb2->second
+                                         : itmea2->second;
+                  a_exp_proc[itmea2->first] = min_exp;
+                  b_exp_proc[itmea2->first] = min_exp;
+                  Trace("nl-ext-comp") << "Common exponent : " << itmea2->first
+                                       << " : " << min_exp << std::endl;
+                }
+              }
+              if (!a_exp_proc.empty())
+              {
+                setMonomialFactor(a, b, a_exp_proc);
+                setMonomialFactor(b, a, b_exp_proc);
+              }
+              /*
+              if( !a_exp_proc.empty() ){
+                //reduction based on common exponents a > 0 => ( a * b
+              <> a * c <=> b <> c ), a < 0 => ( a * b <> a * c <=> b
+              !<> c )  ? }else{ compareMonomial( a, a, a_exp_proc, b,
+              b, b_exp_proc, exp, lemmas );
+              }
+              */
+              compareMonomial(
+                  a, a, a_exp_proc, b, b, b_exp_proc, exp, lemmas, cmp_infers);
+            }
+          }
+        }
+      }
+    }
+  }
+  // remove redundant lemmas, e.g. if a > b, b > c, a > c were
+  // inferred, discard lemma with conclusion a > c
+  Trace("nl-ext-comp") << "Compute redundancies for " << lemmas.size()
+                       << " lemmas." << std::endl;
+  // naive
+  std::vector<Node> r_lemmas;
+  for (std::map<int, std::map<Node, std::map<Node, Node> > >::iterator itb =
+           cmp_infers.begin();
+       itb != cmp_infers.end();
+       ++itb)
+  {
+    for (std::map<Node, std::map<Node, Node> >::iterator itc =
+             itb->second.begin();
+         itc != itb->second.end();
+         ++itc)
+    {
+      for (std::map<Node, Node>::iterator itc2 = itc->second.begin();
+           itc2 != itc->second.end();
+           ++itc2)
+      {
+        std::map<Node, bool> visited;
+        for (std::map<Node, Node>::iterator itc3 = itc->second.begin();
+             itc3 != itc->second.end();
+             ++itc3)
+        {
+          if (itc3->first != itc2->first)
+          {
+            std::vector<Node> exp;
+            if (cmp_holds(itc3->first, itc2->first, itb->second, exp, visited))
+            {
+              r_lemmas.push_back(itc2->second);
+              Trace("nl-ext-comp")
+                  << "...inference of " << itc->first << " > " << itc2->first
+                  << " was redundant." << std::endl;
+              break;
+            }
+          }
+        }
+      }
+    }
+  }
+  std::vector<Node> nr_lemmas;
+  for (unsigned i = 0; i < lemmas.size(); i++)
+  {
+    if (std::find(r_lemmas.begin(), r_lemmas.end(), lemmas[i])
+        == r_lemmas.end())
+    {
+      nr_lemmas.push_back(lemmas[i]);
+    }
+  }
+  // could only take maximal lower/minimial lower bounds?
+
+  Trace("nl-ext-comp") << nr_lemmas.size() << " / " << lemmas.size()
+                       << " were non-redundant." << std::endl;
+  return nr_lemmas;
+}
+
+std::vector<Node> NlSolver::checkTangentPlanes()
+{
+  std::vector<Node> lemmas;
+  Trace("nl-ext") << "Get monomial tangent plane lemmas..." << std::endl;
+  NodeManager* nm = NodeManager::currentNM();
+  const std::map<Node, std::vector<Node> >& ccMap =
+      d_mdb.getContainsChildrenMap();
+  unsigned kstart = d_ms_vars.size();
+  for (unsigned k = kstart; k < d_mterms.size(); k++)
+  {
+    Node t = d_mterms[k];
+    // if this term requires a refinement
+    if (d_tplane_refine.find(t) == d_tplane_refine.end())
+    {
+      continue;
+    }
+    Trace("nl-ext-tplanes")
+        << "Look at monomial requiring refinement : " << t << std::endl;
+    // get a decomposition
+    std::map<Node, std::vector<Node> >::const_iterator it = ccMap.find(t);
+    if (it == ccMap.end())
+    {
+      continue;
+    }
+    std::map<Node, std::map<Node, bool> > dproc;
+    for (unsigned j = 0; j < it->second.size(); j++)
+    {
+      Node tc = it->second[j];
+      if (tc != d_one)
+      {
+        Node tc_diff = d_mdb.getContainsDiffNl(tc, t);
+        Assert(!tc_diff.isNull());
+        Node a = tc < tc_diff ? tc : tc_diff;
+        Node b = tc < tc_diff ? tc_diff : tc;
+        if (dproc[a].find(b) == dproc[a].end())
+        {
+          dproc[a][b] = true;
+          Trace("nl-ext-tplanes")
+              << "  decomposable into : " << a << " * " << b << std::endl;
+          Node a_v_c = d_model.computeAbstractModelValue(a);
+          Node b_v_c = d_model.computeAbstractModelValue(b);
+          // points we will add tangent planes for
+          std::vector<Node> pts[2];
+          pts[0].push_back(a_v_c);
+          pts[1].push_back(b_v_c);
+          // if previously refined
+          bool prevRefine = d_tangent_val_bound[0][a].find(b)
+                            != d_tangent_val_bound[0][a].end();
+          // a_min, a_max, b_min, b_max
+          for (unsigned p = 0; p < 4; p++)
+          {
+            Node curr_v = p <= 1 ? a_v_c : b_v_c;
+            if (prevRefine)
+            {
+              Node pt_v = d_tangent_val_bound[p][a][b];
+              Assert(!pt_v.isNull());
+              if (curr_v != pt_v)
+              {
+                Node do_extend =
+                    nm->mkNode((p == 1 || p == 3) ? GT : LT, curr_v, pt_v);
+                do_extend = Rewriter::rewrite(do_extend);
+                if (do_extend == d_true)
+                {
+                  for (unsigned q = 0; q < 2; q++)
+                  {
+                    pts[p <= 1 ? 0 : 1].push_back(curr_v);
+                    pts[p <= 1 ? 1 : 0].push_back(
+                        d_tangent_val_bound[p <= 1 ? 2 + q : q][a][b]);
+                  }
+                }
+              }
+            }
+            else
+            {
+              d_tangent_val_bound[p][a][b] = curr_v;
+            }
+          }
+
+          for (unsigned p = 0; p < pts[0].size(); p++)
+          {
+            Node a_v = pts[0][p];
+            Node b_v = pts[1][p];
+
+            // tangent plane
+            Node tplane = nm->mkNode(
+                MINUS,
+                nm->mkNode(
+                    PLUS, nm->mkNode(MULT, b_v, a), nm->mkNode(MULT, a_v, b)),
+                nm->mkNode(MULT, a_v, b_v));
+            for (unsigned d = 0; d < 4; d++)
+            {
+              Node aa = nm->mkNode(d == 0 || d == 3 ? GEQ : LEQ, a, a_v);
+              Node ab = nm->mkNode(d == 1 || d == 3 ? GEQ : LEQ, b, b_v);
+              Node conc = nm->mkNode(d <= 1 ? LEQ : GEQ, t, tplane);
+              Node tlem = nm->mkNode(OR, aa.negate(), ab.negate(), conc);
+              Trace("nl-ext-tplanes")
+                  << "Tangent plane lemma : " << tlem << std::endl;
+              lemmas.push_back(tlem);
+            }
+
+            // tangent plane reverse implication
+
+            // t <= tplane -> ( (a <= a_v ^ b >= b_v) v
+            // (a >= a_v ^ b <= b_v) ).
+            // in clause form, the above becomes
+            // t <= tplane -> a <= a_v v b <= b_v.
+            // t <= tplane -> b >= b_v v a >= a_v.
+            Node a_leq_av = nm->mkNode(LEQ, a, a_v);
+            Node b_leq_bv = nm->mkNode(LEQ, b, b_v);
+            Node a_geq_av = nm->mkNode(GEQ, a, a_v);
+            Node b_geq_bv = nm->mkNode(GEQ, b, b_v);
+
+            Node t_leq_tplane = nm->mkNode(LEQ, t, tplane);
+            Node a_leq_av_or_b_leq_bv = nm->mkNode(OR, a_leq_av, b_leq_bv);
+            Node b_geq_bv_or_a_geq_av = nm->mkNode(OR, b_geq_bv, a_geq_av);
+            Node ub_reverse1 =
+                nm->mkNode(OR, t_leq_tplane.negate(), a_leq_av_or_b_leq_bv);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << ub_reverse1
+                << std::endl;
+            lemmas.push_back(ub_reverse1);
+            Node ub_reverse2 =
+                nm->mkNode(OR, t_leq_tplane.negate(), b_geq_bv_or_a_geq_av);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << ub_reverse2
+                << std::endl;
+            lemmas.push_back(ub_reverse2);
+
+            // t >= tplane -> ( (a <= a_v ^ b <= b_v) v
+            // (a >= a_v ^ b >= b_v) ).
+            // in clause form, the above becomes
+            // t >= tplane -> a <= a_v v b >= b_v.
+            // t >= tplane -> b >= b_v v a <= a_v
+            Node t_geq_tplane = nm->mkNode(GEQ, t, tplane);
+            Node a_leq_av_or_b_geq_bv = nm->mkNode(OR, a_leq_av, b_geq_bv);
+            Node a_geq_av_or_b_leq_bv = nm->mkNode(OR, a_geq_av, b_leq_bv);
+            Node lb_reverse1 =
+                nm->mkNode(OR, t_geq_tplane.negate(), a_leq_av_or_b_geq_bv);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << lb_reverse1
+                << std::endl;
+            lemmas.push_back(lb_reverse1);
+            Node lb_reverse2 =
+                nm->mkNode(OR, t_geq_tplane.negate(), a_geq_av_or_b_leq_bv);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << lb_reverse2
+                << std::endl;
+            lemmas.push_back(lb_reverse2);
+          }
+        }
+      }
+    }
+  }
+  Trace("nl-ext") << "...trying " << lemmas.size() << " tangent plane lemmas..."
+                  << std::endl;
+  return lemmas;
+}
+
+std::vector<Node> NlSolver::checkMonomialInferBounds(
+    std::vector<Node>& nt_lemmas,
+    const std::vector<Node>& asserts,
+    const std::vector<Node>& false_asserts)
+{
+  // sort monomials by degree
+  Trace("nl-ext-proc") << "Sort monomials by degree..." << std::endl;
+  d_mdb.sortByDegree(d_ms);
+  // all monomials
+  d_mterms.insert(d_mterms.end(), d_ms_vars.begin(), d_ms_vars.end());
+  d_mterms.insert(d_mterms.end(), d_ms.begin(), d_ms.end());
+
+  const std::map<Node, std::map<Node, ConstraintInfo> >& cim =
+      d_cdb.getConstraints();
+
+  std::vector<Node> lemmas;
+  NodeManager* nm = NodeManager::currentNM();
+  // register constraints
+  Trace("nl-ext-debug") << "Register bound constraints..." << std::endl;
+  for (const Node& lit : asserts)
+  {
+    bool polarity = lit.getKind() != NOT;
+    Node atom = lit.getKind() == NOT ? lit[0] : lit;
+    d_cdb.registerConstraint(atom);
+    bool is_false_lit =
+        std::find(false_asserts.begin(), false_asserts.end(), lit)
+        != false_asserts.end();
+    // add information about bounds to variables
+    std::map<Node, std::map<Node, ConstraintInfo> >::const_iterator itc =
+        cim.find(atom);
+    if (itc == cim.end())
+    {
+      continue;
+    }
+    for (const std::pair<const Node, ConstraintInfo>& itcc : itc->second)
+    {
+      Node x = itcc.first;
+      Node coeff = itcc.second.d_coeff;
+      Node rhs = itcc.second.d_rhs;
+      Kind type = itcc.second.d_type;
+      Node exp = lit;
+      if (!polarity)
+      {
+        // reverse
+        if (type == EQUAL)
+        {
+          // we will take the strict inequality in the direction of the
+          // model
+          Node lhs = ArithMSum::mkCoeffTerm(coeff, x);
+          Node query = nm->mkNode(GT, lhs, rhs);
+          Node query_mv = d_model.computeAbstractModelValue(query);
+          if (query_mv == d_true)
+          {
+            exp = query;
+            type = GT;
+          }
+          else
+          {
+            Assert(query_mv == d_false);
+            exp = nm->mkNode(LT, lhs, rhs);
+            type = LT;
+          }
+        }
+        else
+        {
+          type = negateKind(type);
+        }
+      }
+      // add to status if maximal degree
+      d_ci_max[x][coeff][rhs] = d_cdb.isMaximal(atom, x);
+      if (Trace.isOn("nl-ext-bound-debug2"))
+      {
+        Node t = ArithMSum::mkCoeffTerm(coeff, x);
+        Trace("nl-ext-bound-debug2") << "Add Bound: " << t << " " << type << " "
+                                     << rhs << " by " << exp << std::endl;
+      }
+      bool updated = true;
+      std::map<Node, Kind>::iterator its = d_ci[x][coeff].find(rhs);
+      if (its == d_ci[x][coeff].end())
+      {
+        d_ci[x][coeff][rhs] = type;
+        d_ci_exp[x][coeff][rhs] = exp;
+      }
+      else if (type != its->second)
+      {
+        Trace("nl-ext-bound-debug2")
+            << "Joining kinds : " << type << " " << its->second << std::endl;
+        Kind jk = joinKinds(type, its->second);
+        if (jk == UNDEFINED_KIND)
+        {
+          updated = false;
+        }
+        else if (jk != its->second)
+        {
+          if (jk == type)
+          {
+            d_ci[x][coeff][rhs] = type;
+            d_ci_exp[x][coeff][rhs] = exp;
+          }
+          else
+          {
+            d_ci[x][coeff][rhs] = jk;
+            d_ci_exp[x][coeff][rhs] =
+                nm->mkNode(AND, d_ci_exp[x][coeff][rhs], exp);
+          }
+        }
+        else
+        {
+          updated = false;
+        }
+      }
+      if (Trace.isOn("nl-ext-bound"))
+      {
+        if (updated)
+        {
+          Trace("nl-ext-bound") << "Bound: ";
+          debugPrintBound("nl-ext-bound", coeff, x, d_ci[x][coeff][rhs], rhs);
+          Trace("nl-ext-bound") << " by " << d_ci_exp[x][coeff][rhs];
+          if (d_ci_max[x][coeff][rhs])
+          {
+            Trace("nl-ext-bound") << ", is max degree";
+          }
+          Trace("nl-ext-bound") << std::endl;
+        }
+      }
+      // compute if bound is not satisfied, and store what is required
+      // for a possible refinement
+      if (options::nlExtTangentPlanes())
+      {
+        if (is_false_lit)
+        {
+          d_tplane_refine.insert(x);
+        }
+      }
+    }
+  }
+  // reflexive constraints
+  Node null_coeff;
+  for (unsigned j = 0; j < d_mterms.size(); j++)
+  {
+    Node n = d_mterms[j];
+    d_ci[n][null_coeff][n] = EQUAL;
+    d_ci_exp[n][null_coeff][n] = d_true;
+    d_ci_max[n][null_coeff][n] = false;
+  }
+
+  Trace("nl-ext") << "Get inferred bound lemmas..." << std::endl;
+  const std::map<Node, std::vector<Node> >& cpMap =
+      d_mdb.getContainsParentMap();
+  for (unsigned k = 0; k < d_mterms.size(); k++)
+  {
+    Node x = d_mterms[k];
+    Trace("nl-ext-bound-debug")
+        << "Process bounds for " << x << " : " << std::endl;
+    std::map<Node, std::vector<Node> >::const_iterator itm = cpMap.find(x);
+    if (itm == cpMap.end())
+    {
+      Trace("nl-ext-bound-debug") << "...has no parent monomials." << std::endl;
+      continue;
+    }
+    Trace("nl-ext-bound-debug")
+        << "...has " << itm->second.size() << " parent monomials." << std::endl;
+    // check derived bounds
+    std::map<Node, std::map<Node, std::map<Node, Kind> > >::iterator itc =
+        d_ci.find(x);
+    if (itc == d_ci.end())
+    {
+      continue;
+    }
+    for (std::map<Node, std::map<Node, Kind> >::iterator itcc =
+             itc->second.begin();
+         itcc != itc->second.end();
+         ++itcc)
+    {
+      Node coeff = itcc->first;
+      Node t = ArithMSum::mkCoeffTerm(coeff, x);
+      for (std::map<Node, Kind>::iterator itcr = itcc->second.begin();
+           itcr != itcc->second.end();
+           ++itcr)
+      {
+        Node rhs = itcr->first;
+        // only consider this bound if maximal degree
+        if (!d_ci_max[x][coeff][rhs])
+        {
+          continue;
+        }
+        Kind type = itcr->second;
+        for (unsigned j = 0; j < itm->second.size(); j++)
+        {
+          Node y = itm->second[j];
+          Node mult = d_mdb.getContainsDiff(x, y);
+          // x <k> t => m*x <k'> t  where y = m*x
+          // get the sign of mult
+          Node mmv = d_model.computeConcreteModelValue(mult);
+          Trace("nl-ext-bound-debug2")
+              << "Model value of " << mult << " is " << mmv << std::endl;
+          if (!mmv.isConst())
+          {
+            Trace("nl-ext-bound-debug")
+                << "     ...coefficient " << mult
+                << " is non-constant (probably transcendental)." << std::endl;
+            continue;
+          }
+          int mmv_sign = mmv.getConst<Rational>().sgn();
+          Trace("nl-ext-bound-debug2")
+              << "  sign of " << mmv << " is " << mmv_sign << std::endl;
+          if (mmv_sign == 0)
+          {
+            Trace("nl-ext-bound-debug")
+                << "     ...coefficient " << mult << " is zero." << std::endl;
+            continue;
+          }
+          Trace("nl-ext-bound-debug")
+              << "  from " << x << " * " << mult << " = " << y << " and " << t
+              << " " << type << " " << rhs << ", infer : " << std::endl;
+          Kind infer_type = mmv_sign == -1 ? reverseRelationKind(type) : type;
+          Node infer_lhs = nm->mkNode(MULT, mult, t);
+          Node infer_rhs = nm->mkNode(MULT, mult, rhs);
+          Node infer = nm->mkNode(infer_type, infer_lhs, infer_rhs);
+          Trace("nl-ext-bound-debug") << "     " << infer << std::endl;
+          infer = Rewriter::rewrite(infer);
+          Trace("nl-ext-bound-debug2")
+              << "     ...rewritten : " << infer << std::endl;
+          // check whether it is false in model for abstraction
+          Node infer_mv = d_model.computeAbstractModelValue(infer);
+          Trace("nl-ext-bound-debug")
+              << "       ...infer model value is " << infer_mv << std::endl;
+          if (infer_mv == d_false)
+          {
+            Node exp =
+                nm->mkNode(AND,
+                           nm->mkNode(mmv_sign == 1 ? GT : LT, mult, d_zero),
+                           d_ci_exp[x][coeff][rhs]);
+            Node iblem = nm->mkNode(IMPLIES, exp, infer);
+            Node pr_iblem = iblem;
+            iblem = Rewriter::rewrite(iblem);
+            bool introNewTerms = hasNewMonomials(iblem, d_ms);
+            Trace("nl-ext-bound-lemma")
+                << "*** Bound inference lemma : " << iblem
+                << " (pre-rewrite : " << pr_iblem << ")" << std::endl;
+            // Trace("nl-ext-bound-lemma") << "       intro new
+            // monomials = " << introNewTerms << std::endl;
+            if (!introNewTerms)
+            {
+              lemmas.push_back(iblem);
+            }
+            else
+            {
+              nt_lemmas.push_back(iblem);
+            }
+          }
+        }
+      }
+    }
+  }
+  return lemmas;
+}
+
+std::vector<Node> NlSolver::checkFactoring(
+    const std::vector<Node>& asserts, const std::vector<Node>& false_asserts)
+{
+  std::vector<Node> lemmas;
+  NodeManager* nm = NodeManager::currentNM();
+  Trace("nl-ext") << "Get factoring lemmas..." << std::endl;
+  for (const Node& lit : asserts)
+  {
+    bool polarity = lit.getKind() != NOT;
+    Node atom = lit.getKind() == NOT ? lit[0] : lit;
+    Node litv = d_model.computeConcreteModelValue(lit);
+    bool considerLit = false;
+    // Only consider literals that are in false_asserts.
+    considerLit = std::find(false_asserts.begin(), false_asserts.end(), lit)
+                  != false_asserts.end();
+
+    if (considerLit)
+    {
+      std::map<Node, Node> msum;
+      if (ArithMSum::getMonomialSumLit(atom, msum))
+      {
+        Trace("nl-ext-factor") << "Factoring for literal " << lit
+                               << ", monomial sum is : " << std::endl;
+        if (Trace.isOn("nl-ext-factor"))
+        {
+          ArithMSum::debugPrintMonomialSum(msum, "nl-ext-factor");
+        }
+        std::map<Node, std::vector<Node> > factor_to_mono;
+        std::map<Node, std::vector<Node> > factor_to_mono_orig;
+        for (std::map<Node, Node>::iterator itm = msum.begin();
+             itm != msum.end();
+             ++itm)
+        {
+          if (!itm->first.isNull())
+          {
+            if (itm->first.getKind() == NONLINEAR_MULT)
+            {
+              std::vector<Node> children;
+              for (unsigned i = 0; i < itm->first.getNumChildren(); i++)
+              {
+                children.push_back(itm->first[i]);
+              }
+              std::map<Node, bool> processed;
+              for (unsigned i = 0; i < itm->first.getNumChildren(); i++)
+              {
+                if (processed.find(itm->first[i]) == processed.end())
+                {
+                  processed[itm->first[i]] = true;
+                  children[i] = d_one;
+                  if (!itm->second.isNull())
+                  {
+                    children.push_back(itm->second);
+                  }
+                  Node val = nm->mkNode(MULT, children);
+                  if (!itm->second.isNull())
+                  {
+                    children.pop_back();
+                  }
+                  children[i] = itm->first[i];
+                  val = Rewriter::rewrite(val);
+                  factor_to_mono[itm->first[i]].push_back(val);
+                  factor_to_mono_orig[itm->first[i]].push_back(itm->first);
+                }
+              }
+            }
+          }
+        }
+        for (std::map<Node, std::vector<Node> >::iterator itf =
+                 factor_to_mono.begin();
+             itf != factor_to_mono.end();
+             ++itf)
+        {
+          Node x = itf->first;
+          if (itf->second.size() == 1)
+          {
+            std::map<Node, Node>::iterator itm = msum.find(x);
+            if (itm != msum.end())
+            {
+              itf->second.push_back(itm->second.isNull() ? d_one : itm->second);
+              factor_to_mono_orig[x].push_back(x);
+            }
+          }
+          if (itf->second.size() <= 1)
+          {
+            continue;
+          }
+          Node sum = nm->mkNode(PLUS, itf->second);
+          sum = Rewriter::rewrite(sum);
+          Trace("nl-ext-factor")
+              << "* Factored sum for " << x << " : " << sum << std::endl;
+          Node kf = getFactorSkolem(sum, lemmas);
+          std::vector<Node> poly;
+          poly.push_back(nm->mkNode(MULT, x, kf));
+          std::map<Node, std::vector<Node> >::iterator itfo =
+              factor_to_mono_orig.find(x);
+          Assert(itfo != factor_to_mono_orig.end());
+          for (std::map<Node, Node>::iterator itm = msum.begin();
+               itm != msum.end();
+               ++itm)
+          {
+            if (std::find(itfo->second.begin(), itfo->second.end(), itm->first)
+                == itfo->second.end())
+            {
+              poly.push_back(ArithMSum::mkCoeffTerm(
+                  itm->second, itm->first.isNull() ? d_one : itm->first));
+            }
+          }
+          Node polyn = poly.size() == 1 ? poly[0] : nm->mkNode(PLUS, poly);
+          Trace("nl-ext-factor")
+              << "...factored polynomial : " << polyn << std::endl;
+          Node conc_lit = nm->mkNode(atom.getKind(), polyn, d_zero);
+          conc_lit = Rewriter::rewrite(conc_lit);
+          if (!polarity)
+          {
+            conc_lit = conc_lit.negate();
+          }
+
+          std::vector<Node> lemma_disj;
+          lemma_disj.push_back(lit.negate());
+          lemma_disj.push_back(conc_lit);
+          Node flem = nm->mkNode(OR, lemma_disj);
+          Trace("nl-ext-factor") << "...lemma is " << flem << std::endl;
+          lemmas.push_back(flem);
+        }
+      }
+    }
+  }
+  return lemmas;
+}
+
+Node NlSolver::getFactorSkolem(Node n, std::vector<Node>& lemmas)
+{
+  std::map<Node, Node>::iterator itf = d_factor_skolem.find(n);
+  if (itf == d_factor_skolem.end())
+  {
+    NodeManager* nm = NodeManager::currentNM();
+    Node k = nm->mkSkolem("kf", n.getType());
+    Node k_eq = Rewriter::rewrite(k.eqNode(n));
+    lemmas.push_back(k_eq);
+    d_factor_skolem[n] = k;
+    return k;
+  }
+  return itf->second;
+}
+
+std::vector<Node> NlSolver::checkMonomialInferResBounds()
+{
+  std::vector<Node> lemmas;
+  NodeManager* nm = NodeManager::currentNM();
+  Trace("nl-ext") << "Get monomial resolution inferred bound lemmas..."
+                  << std::endl;
+  size_t nmterms = d_mterms.size();
+  for (unsigned j = 0; j < nmterms; j++)
+  {
+    Node a = d_mterms[j];
+    std::map<Node, std::map<Node, std::map<Node, Kind> > >::iterator itca =
+        d_ci.find(a);
+    if (itca == d_ci.end())
+    {
+      continue;
+    }
+    for (unsigned k = (j + 1); k < nmterms; k++)
+    {
+      Node b = d_mterms[k];
+      std::map<Node, std::map<Node, std::map<Node, Kind> > >::iterator itcb =
+          d_ci.find(b);
+      if (itcb == d_ci.end())
+      {
+        continue;
+      }
+      Trace("nl-ext-rbound-debug") << "resolution inferences : compare " << a
+                                   << " and " << b << std::endl;
+      // if they have common factors
+      std::map<Node, Node>::iterator ita = d_mono_diff[a].find(b);
+      if (ita == d_mono_diff[a].end())
+      {
+        continue;
+      }
+      Trace("nl-ext-rbound") << "Get resolution inferences for [a] " << a
+                             << " vs [b] " << b << std::endl;
+      std::map<Node, Node>::iterator itb = d_mono_diff[b].find(a);
+      Assert(itb != d_mono_diff[b].end());
+      Node mv_a = d_model.computeAbstractModelValue(ita->second);
+      Assert(mv_a.isConst());
+      int mv_a_sgn = mv_a.getConst<Rational>().sgn();
+      if (mv_a_sgn == 0)
+      {
+        // we don't compare monomials whose current model value is zero
+        continue;
+      }
+      Node mv_b = d_model.computeAbstractModelValue(itb->second);
+      Assert(mv_b.isConst());
+      int mv_b_sgn = mv_b.getConst<Rational>().sgn();
+      if (mv_b_sgn == 0)
+      {
+        // we don't compare monomials whose current model value is zero
+        continue;
+      }
+      Trace("nl-ext-rbound") << "  [a] factor is " << ita->second
+                             << ", sign in model = " << mv_a_sgn << std::endl;
+      Trace("nl-ext-rbound") << "  [b] factor is " << itb->second
+                             << ", sign in model = " << mv_b_sgn << std::endl;
+
+      std::vector<Node> exp;
+      // bounds of a
+      for (std::map<Node, std::map<Node, Kind> >::iterator itcac =
+               itca->second.begin();
+           itcac != itca->second.end();
+           ++itcac)
+      {
+        Node coeff_a = itcac->first;
+        for (std::map<Node, Kind>::iterator itcar = itcac->second.begin();
+             itcar != itcac->second.end();
+             ++itcar)
+        {
+          Node rhs_a = itcar->first;
+          Node rhs_a_res_base = nm->mkNode(MULT, itb->second, rhs_a);
+          rhs_a_res_base = Rewriter::rewrite(rhs_a_res_base);
+          if (hasNewMonomials(rhs_a_res_base, d_ms))
+          {
+            continue;
+          }
+          Kind type_a = itcar->second;
+          exp.push_back(d_ci_exp[a][coeff_a][rhs_a]);
+
+          // bounds of b
+          for (std::map<Node, std::map<Node, Kind> >::iterator itcbc =
+                   itcb->second.begin();
+               itcbc != itcb->second.end();
+               ++itcbc)
+          {
+            Node coeff_b = itcbc->first;
+            Node rhs_a_res = ArithMSum::mkCoeffTerm(coeff_b, rhs_a_res_base);
+            for (std::map<Node, Kind>::iterator itcbr = itcbc->second.begin();
+                 itcbr != itcbc->second.end();
+                 ++itcbr)
+            {
+              Node rhs_b = itcbr->first;
+              Node rhs_b_res = nm->mkNode(MULT, ita->second, rhs_b);
+              rhs_b_res = ArithMSum::mkCoeffTerm(coeff_a, rhs_b_res);
+              rhs_b_res = Rewriter::rewrite(rhs_b_res);
+              if (hasNewMonomials(rhs_b_res, d_ms))
+              {
+                continue;
+              }
+              Kind type_b = itcbr->second;
+              exp.push_back(d_ci_exp[b][coeff_b][rhs_b]);
+              if (Trace.isOn("nl-ext-rbound"))
+              {
+                Trace("nl-ext-rbound") << "* try bounds : ";
+                debugPrintBound("nl-ext-rbound", coeff_a, a, type_a, rhs_a);
+                Trace("nl-ext-rbound") << std::endl;
+                Trace("nl-ext-rbound") << "               ";
+                debugPrintBound("nl-ext-rbound", coeff_b, b, type_b, rhs_b);
+                Trace("nl-ext-rbound") << std::endl;
+              }
+              Kind types[2];
+              for (unsigned r = 0; r < 2; r++)
+              {
+                Node pivot_factor = r == 0 ? itb->second : ita->second;
+                int pivot_factor_sign = r == 0 ? mv_b_sgn : mv_a_sgn;
+                types[r] = r == 0 ? type_a : type_b;
+                if (pivot_factor_sign == (r == 0 ? 1 : -1))
+                {
+                  types[r] = reverseRelationKind(types[r]);
+                }
+                if (pivot_factor_sign == 1)
+                {
+                  exp.push_back(nm->mkNode(GT, pivot_factor, d_zero));
+                }
+                else
+                {
+                  exp.push_back(nm->mkNode(LT, pivot_factor, d_zero));
+                }
+              }
+              Kind jk = transKinds(types[0], types[1]);
+              Trace("nl-ext-rbound-debug")
+                  << "trans kind : " << types[0] << " + " << types[1] << " = "
+                  << jk << std::endl;
+              if (jk != UNDEFINED_KIND)
+              {
+                Node conc = nm->mkNode(jk, rhs_a_res, rhs_b_res);
+                Node conc_mv = d_model.computeAbstractModelValue(conc);
+                if (conc_mv == d_false)
+                {
+                  Node rblem = nm->mkNode(IMPLIES, nm->mkNode(AND, exp), conc);
+                  Trace("nl-ext-rbound-lemma-debug")
+                      << "Resolution bound lemma "
+                         "(pre-rewrite) "
+                         ": "
+                      << rblem << std::endl;
+                  rblem = Rewriter::rewrite(rblem);
+                  Trace("nl-ext-rbound-lemma")
+                      << "Resolution bound lemma : " << rblem << std::endl;
+                  lemmas.push_back(rblem);
+                }
+              }
+              exp.pop_back();
+              exp.pop_back();
+              exp.pop_back();
+            }
+          }
+          exp.pop_back();
+        }
+      }
+    }
+  }
+  return lemmas;
+}
+
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl_solver.h b/src/theory/arith/nl_solver.h
new file mode 100644
index 000000000..73ca97f00
--- /dev/null
+++ b/src/theory/arith/nl_solver.h
@@ -0,0 +1,368 @@
+/*********************                                                        */
+/*! \file nl_solver.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2019 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved.  See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** \brief Solver for standard non-linear constraints
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL_SOLVER_H
+#define CVC4__THEORY__ARITH__NL_SOLVER_H
+
+#include <map>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include "context/cdhashset.h"
+#include "context/cdinsert_hashmap.h"
+#include "context/cdlist.h"
+#include "context/cdqueue.h"
+#include "context/context.h"
+#include "expr/kind.h"
+#include "expr/node.h"
+#include "theory/arith/nl_constraint.h"
+#include "theory/arith/nl_lemma_utils.h"
+#include "theory/arith/nl_model.h"
+#include "theory/arith/nl_monomial.h"
+#include "theory/arith/theory_arith.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+
+typedef std::map<Node, unsigned> NodeMultiset;
+
+/** Non-linear solver class
+ *
+ * This class implements model-based refinement schemes
+ * for non-linear arithmetic, described in:
+ *
+ * - "Invariant Checking of NRA Transition Systems
+ * via Incremental Reduction to LRA with EUF" by
+ * Cimatti et al., TACAS 2017.
+ *
+ * - Section 5 of "Desiging Theory Solvers with
+ * Extensions" by Reynolds et al., FroCoS 2017.
+ */
+class NlSolver
+{
+  typedef std::map<Node, NodeMultiset> MonomialExponentMap;
+  typedef context::CDHashSet<Node, NodeHashFunction> NodeSet;
+
+ public:
+  NlSolver(TheoryArith& containing, NlModel& model);
+  ~NlSolver();
+
+  /** init last call
+   *
+   * This is called at the beginning of last call effort check, where
+   * assertions are the set of assertions belonging to arithmetic,
+   * false_asserts is the subset of assertions that are false in the current
+   * model, and xts is the set of extended function terms that are active in
+   * the current context.
+   */
+  void initLastCall(const std::vector<Node>& assertions,
+                    const std::vector<Node>& false_asserts,
+                    const std::vector<Node>& xts);
+  //-------------------------------------------- lemma schemas
+  /** check split zero
+   *
+   * Returns a set of theory lemmas of the form
+   *   t = 0 V t != 0
+   * where t is a term that exists in the current context.
+   */
+  std::vector<Node> checkSplitZero();
+
+  /** check monomial sign
+   *
+   * Returns a set of valid theory lemmas, based on a
+   * lemma schema which ensures that non-linear monomials
+   * respect sign information based on their facts.
+   * For more details, see Section 5 of "Design Theory
+   * Solvers with Extensions" by Reynolds et al., FroCoS 2017,
+   * Figure 5, this is the schema "Sign".
+   *
+   * Examples:
+   *
+   * x > 0 ^ y > 0 => x*y > 0
+   * x < 0 => x*y*y < 0
+   * x = 0 => x*y*z = 0
+   */
+  std::vector<Node> checkMonomialSign();
+
+  /** check monomial magnitude
+   *
+   * Returns a set of valid theory lemmas, based on a
+   * lemma schema which ensures that comparisons between
+   * non-linear monomials respect the magnitude of their
+   * factors.
+   * For more details, see Section 5 of "Design Theory
+   * Solvers with Extensions" by Reynolds et al., FroCoS 2017,
+   * Figure 5, this is the schema "Magnitude".
+   *
+   * Examples:
+   *
+   * |x|>|y| => |x*z|>|y*z|
+   * |x|>|y| ^ |z|>|w| ^ |x|>=1 => |x*x*z*u|>|y*w|
+   *
+   * Argument c indicates the class of inferences to perform for the
+   * (non-linear) monomials in the vector d_ms. 0 : compare non-linear monomials
+   * against 1, 1 : compare non-linear monomials against variables, 2 : compare
+   * non-linear monomials against other non-linear monomials.
+   */
+  std::vector<Node> checkMonomialMagnitude(unsigned c);
+
+  /** check monomial inferred bounds
+   *
+   * Returns a set of valid theory lemmas, based on a
+   * lemma schema that infers new constraints about existing
+   * terms based on mulitplying both sides of an existing
+   * constraint by a term.
+   * For more details, see Section 5 of "Design Theory
+   * Solvers with Extensions" by Reynolds et al., FroCoS 2017,
+   * Figure 5, this is the schema "Multiply".
+   *
+   * Examples:
+   *
+   * x > 0 ^ (y > z + w) => x*y > x*(z+w)
+   * x < 0 ^ (y > z + w) => x*y < x*(z+w)
+   *   ...where (y > z + w) and x*y are a constraint and term
+   *      that occur in the current context.
+   */
+  std::vector<Node> checkMonomialInferBounds(
+      std::vector<Node>& nt_lemmas,
+      const std::vector<Node>& asserts,
+      const std::vector<Node>& false_asserts);
+
+  /** check factoring
+   *
+   * Returns a set of valid theory lemmas, based on a
+   * lemma schema that states a relationship betwen monomials
+   * with common factors that occur in the same constraint.
+   *
+   * Examples:
+   *
+   * x*z+y*z > t => ( k = x + y ^ k*z > t )
+   *   ...where k is fresh and x*z + y*z > t is a
+   *      constraint that occurs in the current context.
+   */
+  std::vector<Node> checkFactoring(const std::vector<Node>& asserts,
+                                   const std::vector<Node>& false_asserts);
+
+  /** check monomial infer resolution bounds
+   *
+   * Returns a set of valid theory lemmas, based on a
+   * lemma schema which "resolves" upper bounds
+   * of one inequality with lower bounds for another.
+   * This schema is not enabled by default, and can
+   * be enabled by --nl-ext-rbound.
+   *
+   * Examples:
+   *
+   *  ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t
+   *  ...where s <= x*z and x*y <= t are constraints
+   *     that occur in the current context.
+   */
+  std::vector<Node> checkMonomialInferResBounds();
+
+  /** check tangent planes
+   *
+   * Returns a set of valid theory lemmas, based on an
+   * "incremental linearization" of non-linear monomials.
+   * This linearization is accomplished by adding constraints
+   * corresponding to "tangent planes" at the current
+   * model value of each non-linear monomial. In particular
+   * consider the definition for constants a,b :
+   *   T_{a,b}( x*y ) = b*x + a*y - a*b.
+   * The lemmas added by this function are of the form :
+   *  ( ( x>a ^ y<b) ^ (x<a ^ y>b) ) => x*y < T_{a,b}( x*y )
+   *  ( ( x>a ^ y>b) ^ (x<a ^ y<b) ) => x*y > T_{a,b}( x*y )
+   * It is inspired by "Invariant Checking of NRA Transition
+   * Systems via Incremental Reduction to LRA with EUF" by
+   * Cimatti et al., TACAS 2017.
+   * This schema is not terminating in general.
+   * It is not enabled by default, and can
+   * be enabled by --nl-ext-tplanes.
+   *
+   * Examples:
+   *
+   * ( ( x>2 ^ y>5) ^ (x<2 ^ y<5) ) => x*y > 5*x + 2*y - 10
+   * ( ( x>2 ^ y<5) ^ (x<2 ^ y>5) ) => x*y < 5*x + 2*y - 10
+   */
+  std::vector<Node> checkTangentPlanes();
+
+  //-------------------------------------------- end lemma schemas
+ private:
+  // The theory of arithmetic containing this extension.
+  TheoryArith& d_containing;
+  /** Reference to the non-linear model object */
+  NlModel& d_model;
+  /** commonly used terms */
+  Node d_zero;
+  Node d_one;
+  Node d_neg_one;
+  Node d_two;
+  Node d_true;
+  Node d_false;
+  /** Context-independent database of monomial information */
+  MonomialDb d_mdb;
+  /** Context-independent database of constraint information */
+  ConstraintDb d_cdb;
+
+  // ( x*y, x*z, y ) for each pair of monomials ( x*y, x*z ) with common factors
+  std::map<Node, std::map<Node, Node> > d_mono_diff;
+
+  /** cache of terms t for which we have added the lemma ( t = 0 V t != 0 ). */
+  NodeSet d_zero_split;
+
+  // ordering, stores variables and 0,1,-1
+  std::map<Node, unsigned> d_order_vars;
+  std::vector<Node> d_order_points;
+
+  // information about monomials
+  std::vector<Node> d_ms;
+  std::vector<Node> d_ms_vars;
+  std::map<Node, bool> d_ms_proc;
+  std::vector<Node> d_mterms;
+
+  // list of monomials with factors whose model value is non-constant in model
+  //  e.g. y*cos( x )
+  std::map<Node, bool> d_m_nconst_factor;
+  /** the set of monomials we should apply tangent planes to */
+  std::unordered_set<Node, NodeHashFunction> d_tplane_refine;
+  /** maps nodes to their factor skolems */
+  std::map<Node, Node> d_factor_skolem;
+  /** tangent plane bounds */
+  std::map<Node, std::map<Node, Node> > d_tangent_val_bound[4];
+  // term -> coeff -> rhs -> ( status, exp, b ),
+  //   where we have that : exp =>  ( coeff * term <status> rhs )
+  //   b is true if degree( term ) >= degree( rhs )
+  std::map<Node, std::map<Node, std::map<Node, Kind> > > d_ci;
+  std::map<Node, std::map<Node, std::map<Node, Node> > > d_ci_exp;
+  std::map<Node, std::map<Node, std::map<Node, bool> > > d_ci_max;
+
+  /** Make literal */
+  static Node mkLit(Node a, Node b, int status, bool isAbsolute = false);
+  /** register monomial */
+  void setMonomialFactor(Node a, Node b, const NodeMultiset& common);
+  /** assign order ids */
+  void assignOrderIds(std::vector<Node>& vars,
+                      NodeMultiset& d_order,
+                      bool isConcrete,
+                      bool isAbsolute);
+
+  /** Check whether we have already inferred a relationship between monomials
+   * x and y based on the information in cmp_infers. This computes the
+   * transitive closure of the relation stored in cmp_infers.
+   */
+  bool cmp_holds(Node x,
+                 Node y,
+                 std::map<Node, std::map<Node, Node> >& cmp_infers,
+                 std::vector<Node>& exp,
+                 std::map<Node, bool>& visited);
+  /** In the following functions, status states a relationship
+   * between two arithmetic terms, where:
+   * 0 : equal
+   * 1 : greater than or equal
+   * 2 : greater than
+   * -X : (greater -> less)
+   * TODO (#1287) make this an enum?
+   */
+  /** compute the sign of a.
+   *
+   * Calls to this function are such that :
+   *    exp => ( oa = a ^ a <status> 0 )
+   *
+   * This function iterates over the factors of a,
+   * where a_index is the index of the factor in a
+   * we are currently looking at.
+   *
+   * This function returns a status, which indicates
+   * a's relationship to 0.
+   * We add lemmas to lem of the form given by the
+   * lemma schema checkSign(...).
+   */
+  int compareSign(Node oa,
+                  Node a,
+                  unsigned a_index,
+                  int status,
+                  std::vector<Node>& exp,
+                  std::vector<Node>& lem);
+  /** compare monomials a and b
+   *
+   * Initially, a call to this function is such that :
+   *    exp => ( oa = a ^ ob = b )
+   *
+   * This function returns true if we can infer a valid
+   * arithmetic lemma of the form :
+   *    P => abs( a ) >= abs( b )
+   * where P is true and abs( a ) >= abs( b ) is false in the
+   * current model.
+   *
+   * This function is implemented by "processing" factors
+   * of monomials a and b until an inference of the above
+   * form can be made. For example, if :
+   *   a = x*x*y and b = z*w
+   * Assuming we are trying to show abs( a ) >= abs( c ),
+   * then if abs( M( x ) ) >= abs( M( z ) ) where M is the current model,
+   * then we can add abs( x ) >= abs( z ) to our explanation, and
+   * mark one factor of x as processed in a, and
+   * one factor of z as processed in b. The number of processed factors of a
+   * and b are stored in a_exp_proc and b_exp_proc respectively.
+   *
+   * cmp_infers stores information that is helpful
+   * in discarding redundant inferences.  For example,
+   * we do not want to infer abs( x ) >= abs( z ) if
+   * we have already inferred abs( x ) >= abs( y ) and
+   * abs( y ) >= abs( z ).
+   * It stores entries of the form (status,t1,t2)->F,
+   * which indicates that we constructed a lemma F that
+   * showed t1 <status> t2.
+   *
+   * We add lemmas to lem of the form given by the
+   * lemma schema checkMagnitude(...).
+   */
+  bool compareMonomial(
+      Node oa,
+      Node a,
+      NodeMultiset& a_exp_proc,
+      Node ob,
+      Node b,
+      NodeMultiset& b_exp_proc,
+      std::vector<Node>& exp,
+      std::vector<Node>& lem,
+      std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers);
+  /** helper function for above
+   *
+   * The difference is the inputs a_index and b_index, which are the indices of
+   * children (factors) in monomials a and b which we are currently looking at.
+   */
+  bool compareMonomial(
+      Node oa,
+      Node a,
+      unsigned a_index,
+      NodeMultiset& a_exp_proc,
+      Node ob,
+      Node b,
+      unsigned b_index,
+      NodeMultiset& b_exp_proc,
+      int status,
+      std::vector<Node>& exp,
+      std::vector<Node>& lem,
+      std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers);
+  /** Get factor skolem for n, add resulting lemmas to lemmas */
+  Node getFactorSkolem(Node n, std::vector<Node>& lemmas);
+}; /* class NlSolver */
+
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__NL_SOLVER_H */
diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp
index 61a8b18b0..b638d8a59 100644
--- a/src/theory/arith/nonlinear_extension.cpp
+++ b/src/theory/arith/nonlinear_extension.cpp
@@ -17,17 +17,10 @@
 
 #include "theory/arith/nonlinear_extension.h"
 
-#include <cmath>
-#include <set>
-
-#include "expr/node_algorithm.h"
-#include "expr/node_builder.h"
 #include "options/arith_options.h"
-#include "theory/arith/arith_msum.h"
 #include "theory/arith/arith_utilities.h"
 #include "theory/arith/theory_arith.h"
 #include "theory/ext_theory.h"
-#include "theory/quantifiers/quant_util.h"
 #include "theory/theory_model.h"
 
 using namespace CVC4::kind;
@@ -36,181 +29,25 @@ namespace CVC4 {
 namespace theory {
 namespace arith {
 
-namespace {
-
-// Return true if a collection c contains an elem k. Compatible with map and set
-// containers.
-template <class Container, class Key>
-bool Contains(const Container& c, const Key& k) {
-  return c.find(k) != c.end();
-}
-
-// Inserts value into the set/map c if the value was not present there. Returns
-// true if the value was inserted.
-template <class Container, class Value>
-bool InsertIfNotPresent(Container* c, const Value& value) {
-  return (c->insert(value)).second;
-}
-
-// Returns true if a vector c contains an elem t.
-template <class T>
-bool IsInVector(const std::vector<T>& c, const T& t) {
-  return std::find(c.begin(), c.end(), t) != c.end();
-}
-
-// Returns the a[key] and assertion fails in debug mode.
-inline unsigned getCount(const NodeMultiset& a, Node key) {
-  NodeMultiset::const_iterator it = a.find(key);
-  Assert(it != a.end());
-  return it->second;
-}
-
-// Returns a[key] if key is in a or value otherwise.
-unsigned getCountWithDefault(const NodeMultiset& a, Node key, unsigned value) {
-  NodeMultiset::const_iterator it = a.find(key);
-  return (it == a.end()) ? value : it->second;
-}
-
-// Returns true if for any key then a[key] <= b[key] where the value for any key
-// not present is interpreted as 0.
-bool isSubset(const NodeMultiset& a, const NodeMultiset& b) {
-  for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) {
-    Node key = it_a->first;
-    const unsigned a_value = it_a->second;
-    const unsigned b_value = getCountWithDefault(b, key, 0);
-    if (a_value > b_value) {
-      return false;
-    }
-  }
-  return true;
-}
-
-// Given two multisets return the multiset difference a \ b.
-NodeMultiset diffMultiset(const NodeMultiset& a, const NodeMultiset& b) {
-  NodeMultiset difference;
-  for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) {
-    Node key = it_a->first;
-    const unsigned a_value = it_a->second;
-    const unsigned b_value = getCountWithDefault(b, key, 0);
-    if (a_value > b_value) {
-      difference[key] = a_value - b_value;
-    }
-  }
-  return difference;
-}
-
-// Return a vector containing a[key] repetitions of key in a multiset a.
-std::vector<Node> ExpandMultiset(const NodeMultiset& a) {
-  std::vector<Node> expansion;
-  for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) {
-    expansion.insert(expansion.end(), it_a->second, it_a->first);
-  }
-  return expansion;
-}
-
-void debugPrintBound(const char* c, Node coeff, Node x, Kind type, Node rhs) {
-  Node t = ArithMSum::mkCoeffTerm(coeff, x);
-  Trace(c) << t << " " << type << " " << rhs;
-}
-
-bool hasNewMonomials(Node n, const std::vector<Node>& existing) {
-  std::set<Node> visited;
-
-  std::vector<Node> worklist;
-  worklist.push_back(n);
-  while (!worklist.empty()) {
-    Node current = worklist.back();
-    worklist.pop_back();
-    if (!Contains(visited, current)) {
-      visited.insert(current);
-      if (current.getKind() == NONLINEAR_MULT)
-      {
-        if (!IsInVector(existing, current)) {
-          return true;
-        }
-      } else {
-        worklist.insert(worklist.end(), current.begin(), current.end());
-      }
-    }
-  }
-  return false;
-}
-
-}  // namespace
-
 NonlinearExtension::NonlinearExtension(TheoryArith& containing,
                                        eq::EqualityEngine* ee)
     : d_lemmas(containing.getUserContext()),
-      d_zero_split(containing.getUserContext()),
       d_containing(containing),
       d_ee(ee),
       d_needsLastCall(false),
       d_model(containing.getSatContext()),
       d_trSlv(d_model),
+      d_nlSlv(containing, d_model),
       d_builtModel(containing.getSatContext(), false)
 {
   d_true = NodeManager::currentNM()->mkConst(true);
-  d_false = NodeManager::currentNM()->mkConst(false);
   d_zero = NodeManager::currentNM()->mkConst(Rational(0));
   d_one = NodeManager::currentNM()->mkConst(Rational(1));
   d_neg_one = NodeManager::currentNM()->mkConst(Rational(-1));
-  d_two = NodeManager::currentNM()->mkConst(Rational(2));
-  d_order_points.push_back(d_neg_one);
-  d_order_points.push_back(d_zero);
-  d_order_points.push_back(d_one);
 }
 
 NonlinearExtension::~NonlinearExtension() {}
 
-// Returns a reference to either map[key] if it exists in the map
-// or to a default value otherwise.
-//
-// Warning: sped_cial care must be taken if value is a temporary object.
-template <class MapType, class Key, class Value>
-const Value& FindWithDefault(const MapType& map, const Key& key,
-                             const Value& value) {
-  typename MapType::const_iterator it = map.find(key);
-  if (it == map.end()) {
-    return value;
-  }
-  return it->second;
-}
-
-const NodeMultiset& NonlinearExtension::getMonomialExponentMap(
-    Node monomial) const {
-  MonomialExponentMap::const_iterator it = d_m_exp.find(monomial);
-  Assert(it != d_m_exp.end());
-  return it->second;
-}
-
-bool NonlinearExtension::isMonomialSubset(Node a, Node b) const {
-  const NodeMultiset& a_exponent_map = getMonomialExponentMap(a);
-  const NodeMultiset& b_exponent_map = getMonomialExponentMap(b);
-
-  return isSubset(a_exponent_map, b_exponent_map);
-}
-
-void NonlinearExtension::registerMonomialSubset(Node a, Node b) {
-  Assert(isMonomialSubset(a, b));
-
-  const NodeMultiset& a_exponent_map = getMonomialExponentMap(a);
-  const NodeMultiset& b_exponent_map = getMonomialExponentMap(b);
-
-  std::vector<Node> diff_children =
-      ExpandMultiset(diffMultiset(b_exponent_map, a_exponent_map));
-  Assert(!diff_children.empty());
-
-  d_m_contain_parent[a].push_back(b);
-  d_m_contain_children[b].push_back(a);
-
-  Node mult_term = safeConstructNary(MULT, diff_children);
-  Node nlmult_term = safeConstructNary(NONLINEAR_MULT, diff_children);
-  d_m_contain_mult[a][b] = mult_term;
-  d_m_contain_umult[a][b] = nlmult_term;
-  Trace("nl-ext-mindex") << "..." << a << " is a subset of " << b
-                         << ", difference is " << mult_term << std::endl;
-}
-
 bool NonlinearExtension::getCurrentSubstitution(
     int effort, const std::vector<Node>& vars, std::vector<Node>& subs,
     std::map<Node, std::vector<Node> >& exp) {
@@ -239,8 +76,6 @@ bool NonlinearExtension::getCurrentSubstitution(
 
   // return true if the substitution is non-trivial
   return retVal;
-
-  // d_containing.getValuation().getModel()->getRepresentative( n );
 }
 
 std::pair<bool, Node> NonlinearExtension::isExtfReduced(
@@ -295,194 +130,6 @@ std::pair<bool, Node> NonlinearExtension::isExtfReduced(
   return std::make_pair(true, Node::null());
 }
 
-void NonlinearExtension::registerMonomial(Node n) {
-  if (!IsInVector(d_monomials, n)) {
-    d_monomials.push_back(n);
-    Trace("nl-ext-debug") << "Register monomial : " << n << std::endl;
-    if (n.getKind() == NONLINEAR_MULT)
-    {
-      // get exponent count
-      for (unsigned k = 0; k < n.getNumChildren(); k++) {
-        d_m_exp[n][n[k]]++;
-        if (k == 0 || n[k] != n[k - 1]) {
-          d_m_vlist[n].push_back(n[k]);
-        }
-      }
-      d_m_degree[n] = n.getNumChildren();
-    } else if (n == d_one) {
-      d_m_exp[n].clear();
-      d_m_vlist[n].clear();
-      d_m_degree[n] = 0;
-    } else {
-      Assert(!isArithKind(n.getKind()));
-      d_m_exp[n][n] = 1;
-      d_m_vlist[n].push_back(n);
-      d_m_degree[n] = 1;
-    }
-    // TODO: sort necessary here?
-    std::sort(d_m_vlist[n].begin(), d_m_vlist[n].end());
-    Trace("nl-ext-mindex") << "Add monomial to index : " << n << std::endl;
-    d_m_index.addTerm(n, d_m_vlist[n], this);
-  }
-}
-
-void NonlinearExtension::setMonomialFactor(Node a, Node b,
-                                           const NodeMultiset& common) {
-  // Could not tell if this was being inserted intentionally or not.
-  std::map<Node, Node>& mono_diff_a = d_mono_diff[a];
-  if (!Contains(mono_diff_a, b)) {
-    Trace("nl-ext-mono-factor")
-        << "Set monomial factor for " << a << "/" << b << std::endl;
-    mono_diff_a[b] = mkMonomialRemFactor(a, common);
-  }
-}
-
-void NonlinearExtension::registerConstraint(Node atom) {
-  if (!IsInVector(d_constraints, atom)) {
-    d_constraints.push_back(atom);
-    Trace("nl-ext-debug") << "Register constraint : " << atom << std::endl;
-    std::map<Node, Node> msum;
-    if (ArithMSum::getMonomialSumLit(atom, msum))
-    {
-      Trace("nl-ext-debug") << "got monomial sum: " << std::endl;
-      if (Trace.isOn("nl-ext-debug")) {
-        ArithMSum::debugPrintMonomialSum(msum, "nl-ext-debug");
-      }
-      unsigned max_degree = 0;
-      std::vector<Node> all_m;
-      std::vector<Node> max_deg_m;
-      for (std::map<Node, Node>::iterator itm = msum.begin(); itm != msum.end();
-           ++itm) {
-        if (!itm->first.isNull()) {
-          all_m.push_back(itm->first);
-          registerMonomial(itm->first);
-          Trace("nl-ext-debug2")
-              << "...process monomial " << itm->first << std::endl;
-          Assert(d_m_degree.find(itm->first) != d_m_degree.end());
-          unsigned d = d_m_degree[itm->first];
-          if (d > max_degree) {
-            max_degree = d;
-            max_deg_m.clear();
-          }
-          if (d >= max_degree) {
-            max_deg_m.push_back(itm->first);
-          }
-        }
-      }
-      // isolate for each maximal degree monomial
-      for (unsigned i = 0; i < all_m.size(); i++) {
-        Node m = all_m[i];
-        Node rhs, coeff;
-        int res = ArithMSum::isolate(m, msum, coeff, rhs, atom.getKind());
-        if (res != 0) {
-          Kind type = atom.getKind();
-          if (res == -1) {
-            type = reverseRelationKind(type);
-          }
-          Trace("nl-ext-constraint") << "Constraint : " << atom << " <=> ";
-          if (!coeff.isNull()) {
-            Trace("nl-ext-constraint") << coeff << " * ";
-          }
-          Trace("nl-ext-constraint")
-              << m << " " << type << " " << rhs << std::endl;
-          d_c_info[atom][m].d_rhs = rhs;
-          d_c_info[atom][m].d_coeff = coeff;
-          d_c_info[atom][m].d_type = type;
-        }
-      }
-      for (unsigned i = 0; i < max_deg_m.size(); i++) {
-        Node m = max_deg_m[i];
-        d_c_info_maxm[atom][m] = true;
-      }
-    } else {
-      Trace("nl-ext-debug") << "...failed to get monomial sum." << std::endl;
-    }
-  }
-}
-
-bool NonlinearExtension::isArithKind(Kind k) {
-  return k == PLUS || k == MULT || k == NONLINEAR_MULT;
-}
-
-Node NonlinearExtension::mkLit(Node a, Node b, int status, bool isAbsolute)
-{
-  if (status == 0) {
-    Node a_eq_b = a.eqNode(b);
-    if (!isAbsolute)
-    {
-      return a_eq_b;
-    }
-    else
-    {
-      // return mkAbs( a ).eqNode( mkAbs( b ) );
-      Node negate_b = NodeManager::currentNM()->mkNode(UMINUS, b);
-      return a_eq_b.orNode(a.eqNode(negate_b));
-    }
-  } else if (status < 0) {
-    return mkLit(b, a, -status);
-  } else {
-    Assert(status == 1 || status == 2);
-    NodeManager* nm = NodeManager::currentNM();
-    Kind greater_op = status == 1 ? GEQ : GT;
-    if (!isAbsolute)
-    {
-      return nm->mkNode(greater_op, a, b);
-    }
-    else
-    {
-      // return nm->mkNode( greater_op, mkAbs( a ), mkAbs( b ) );
-      Node zero = mkRationalNode(0);
-      Node a_is_nonnegative = nm->mkNode(GEQ, a, zero);
-      Node b_is_nonnegative = nm->mkNode(GEQ, b, zero);
-      Node negate_a = nm->mkNode(UMINUS, a);
-      Node negate_b = nm->mkNode(UMINUS, b);
-      return a_is_nonnegative.iteNode(
-          b_is_nonnegative.iteNode(nm->mkNode(greater_op, a, b),
-                                   nm->mkNode(greater_op, a, negate_b)),
-          b_is_nonnegative.iteNode(nm->mkNode(greater_op, negate_a, b),
-                                   nm->mkNode(greater_op, negate_a, negate_b)));
-    }
-  }
-}
-
-Node NonlinearExtension::mkAbs(Node a) {
-  if (a.isConst()) {
-    return mkRationalNode(a.getConst<Rational>().abs());
-  } else {
-    NodeManager* nm = NodeManager::currentNM();
-    Node a_is_nonnegative = nm->mkNode(GEQ, a, mkRationalNode(0));
-    return a_is_nonnegative.iteNode(a, nm->mkNode(UMINUS, a));
-  }
-}
-
-Node NonlinearExtension::mkValidPhase(Node a, Node pi) {
-  return mkBounded(
-      NodeManager::currentNM()->mkNode(MULT, mkRationalNode(-1), pi), a, pi);
-}
-
-Node NonlinearExtension::mkMonomialRemFactor(
-    Node n, const NodeMultiset& n_exp_rem) const {
-  std::vector<Node> children;
-  const NodeMultiset& exponent_map = getMonomialExponentMap(n);
-  for (NodeMultiset::const_iterator itme2 = exponent_map.begin();
-       itme2 != exponent_map.end(); ++itme2) {
-    Node v = itme2->first;
-    unsigned inc = itme2->second;
-    Trace("nl-ext-mono-factor")
-        << "..." << inc << " factors of " << v << std::endl;
-    unsigned count_in_n_exp_rem = getCountWithDefault(n_exp_rem, v, 0);
-    Assert(count_in_n_exp_rem <= inc);
-    inc -= count_in_n_exp_rem;
-    Trace("nl-ext-mono-factor")
-        << "......rem, now " << inc << " factors of " << v << std::endl;
-    children.insert(children.end(), inc, v);
-  }
-  Node ret = safeConstructNary(MULT, children);
-  ret = Rewriter::rewrite(ret);
-  Trace("nl-ext-mono-factor") << "...return : " << ret << std::endl;
-  return ret;
-}
-
 void NonlinearExtension::sendLemmas(const std::vector<Node>& out,
                                     bool preprocess,
                                     std::map<Node, NlLemmaSideEffect>& lemSE)
@@ -732,22 +379,6 @@ bool NonlinearExtension::checkModel(const std::vector<Node>& assertions,
   return ret;
 }
 
-
-std::vector<Node> NonlinearExtension::checkSplitZero() {
-  std::vector<Node> lemmas;
-  for (unsigned i = 0; i < d_ms_vars.size(); i++) {
-    Node v = d_ms_vars[i];
-    if (d_zero_split.insert(v)) {
-      Node eq = v.eqNode(d_zero);
-      eq = Rewriter::rewrite(eq);
-      Node literal = d_containing.getValuation().ensureLiteral(eq);
-      d_containing.getOutputChannel().requirePhase(literal, true);
-      lemmas.push_back(literal.orNode(literal.negate()));
-    }
-  }
-  return lemmas;
-}
-
 int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
                                       const std::vector<Node>& false_asserts,
                                       const std::vector<Node>& xts,
@@ -756,48 +387,8 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
                                       std::vector<Node>& wlems,
                                       std::map<Node, NlLemmaSideEffect>& lemSE)
 {
-  d_ms_vars.clear();
-  d_ms_proc.clear();
-  d_ms.clear();
-  d_mterms.clear();
-  d_m_nconst_factor.clear();
-  d_tplane_refine.clear();
-  d_ci.clear();
-  d_ci_exp.clear();
-  d_ci_max.clear();
-
-  Trace("nl-ext-mv") << "Extended terms : " << std::endl;
-  // for computing congruence
-  std::map<Kind, ArgTrie> argTrie;
-  for (unsigned i = 0, xsize = xts.size(); i < xsize; i++)
-  {
-    Node a = xts[i];
-    d_model.computeConcreteModelValue(a);
-    d_model.computeAbstractModelValue(a);
-    d_model.printModelValue("nl-ext-mv", a);
-    Kind ak = a.getKind();
-    if (ak == NONLINEAR_MULT)
-    {
-      d_ms.push_back( a );
-      
-      //context-independent registration
-      registerMonomial(a);
-
-      std::map<Node, std::vector<Node> >::iterator itvl = d_m_vlist.find(a);
-      Assert(itvl != d_m_vlist.end());
-      for (unsigned k = 0; k < itvl->second.size(); k++) {
-        if (!IsInVector(d_ms_vars, itvl->second[k])) {
-          d_ms_vars.push_back(itvl->second[k]);
-        }
-        Node mvk = d_model.computeAbstractModelValue(itvl->second[k]);
-        if( !mvk.isConst() ){
-          d_m_nconst_factor[a] = true;
-        }
-      }
-      // mark processed if has a "one" factor (will look at reduced monomial)?
-    }
-  }
-
+  // initialize the non-linear solver
+  d_nlSlv.initLastCall(assertions, false_asserts, xts);
   // initialize the trancendental function solver
   std::vector<Node> lemmas;
   d_trSlv.initLastCall(assertions, false_asserts, xts, lemmas, lemsPp);
@@ -810,30 +401,11 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
                     << " new lemmas during registration." << std::endl;
     return lems.size() + lemsPp.size();
   }
-  Trace("nl-ext") << "We have " << d_ms.size() << " monomials." << std::endl;
-
-  // register constants
-  registerMonomial(d_one);
-  for (unsigned j = 0; j < d_order_points.size(); j++) {
-    Node c = d_order_points[j];
-    d_model.computeConcreteModelValue(c);
-    d_model.computeAbstractModelValue(c);
-  }
-
-  // register variables
-  Trace("nl-ext-mv") << "Variables in monomials : " << std::endl;
-  for (unsigned i = 0; i < d_ms_vars.size(); i++) {
-    Node v = d_ms_vars[i];
-    registerMonomial(v);
-    d_model.computeConcreteModelValue(v);
-    d_model.computeAbstractModelValue(v);
-    d_model.printModelValue("nl-ext-mv", v);
-  }
 
   //----------------------------------- possibly split on zero
   if (options::nlExtSplitZero()) {
     Trace("nl-ext") << "Get zero split lemmas..." << std::endl;
-    lemmas = checkSplitZero();
+    lemmas = d_nlSlv.checkSplitZero();
     filterLemmas(lemmas, lems);
     if (!lems.empty())
     {
@@ -854,7 +426,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
   }
 
   //-----------------------------------lemmas based on sign (comparison to zero)
-  lemmas = checkMonomialSign();
+  lemmas = d_nlSlv.checkMonomialSign();
   filterLemmas(lemmas, lems);
   if (!lems.empty())
   {
@@ -874,23 +446,9 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
   }
 
   //-----------------------------------lemmas based on magnitude of non-zero monomials
-  Trace("nl-ext-proc") << "Assign order ids..." << std::endl;
-  // sort by absolute values of abstract model values
-  assignOrderIds(d_ms_vars, d_order_vars, false, true);
-
-  // sort individual variable lists
-  Trace("nl-ext-proc") << "Assign order var lists..." << std::endl;
-  SortNlModel smv;
-  smv.d_nlm = &d_model;
-  smv.d_isConcrete = false;
-  smv.d_isAbsolute = true;
-  smv.d_reverse_order = true;
-  for (unsigned j = 0; j < d_ms.size(); j++) {
-    std::sort(d_m_vlist[d_ms[j]].begin(), d_m_vlist[d_ms[j]].end(), smv);
-  }
   for (unsigned c = 0; c < 3; c++) {
     // c is effort level
-    lemmas = checkMonomialMagnitude( c );
+    lemmas = d_nlSlv.checkMonomialMagnitude(c);
     unsigned nlem = lemmas.size();
     filterLemmas(lemmas, lems);
     if (!lems.empty())
@@ -902,18 +460,11 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
     }
   }
 
-  // sort monomials by degree
-  Trace("nl-ext-proc") << "Sort monomials by degree..." << std::endl;
-  SortNonlinearDegree snlad(d_m_degree);
-  std::sort(d_ms.begin(), d_ms.end(), snlad);
-  // all monomials
-  d_mterms.insert(d_mterms.end(), d_ms_vars.begin(), d_ms_vars.end());
-  d_mterms.insert(d_mterms.end(), d_ms.begin(), d_ms.end());
-
   //-----------------------------------inferred bounds lemmas
   //  e.g. x >= t => y*x >= y*t
   std::vector< Node > nt_lemmas;
-  lemmas = checkMonomialInferBounds(nt_lemmas, assertions, false_asserts);
+  lemmas =
+      d_nlSlv.checkMonomialInferBounds(nt_lemmas, assertions, false_asserts);
   // Trace("nl-ext") << "Bound lemmas : " << lemmas.size() << ", " <<
   // nt_lemmas.size() << std::endl;  prioritize lemmas that do not
   // introduce new monomials
@@ -921,7 +472,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
 
   if (options::nlExtTangentPlanes() && options::nlExtTangentPlanesInterleave())
   {
-    lemmas = checkTangentPlanes();
+    lemmas = d_nlSlv.checkTangentPlanes();
     filterLemmas(lemmas, lems);
   }
 
@@ -944,7 +495,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
   //------------------------------------factoring lemmas
   //   x*y + x*z >= t => exists k. k = y + z ^ x*k >= t
   if( options::nlExtFactor() ){
-    lemmas = checkFactoring(assertions, false_asserts);
+    lemmas = d_nlSlv.checkFactoring(assertions, false_asserts);
     filterLemmas(lemmas, lems);
     if (!lems.empty())
     {
@@ -957,7 +508,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
   //------------------------------------resolution bound inferences
   //  e.g. ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t
   if (options::nlExtResBound()) {
-    lemmas = checkMonomialInferResBounds();
+    lemmas = d_nlSlv.checkMonomialInferResBounds();
     filterLemmas(lemmas, lems);
     if (!lems.empty())
     {
@@ -970,7 +521,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
   //------------------------------------tangent planes
   if (options::nlExtTangentPlanes() && !options::nlExtTangentPlanesInterleave())
   {
-    lemmas = checkTangentPlanes();
+    lemmas = d_nlSlv.checkTangentPlanes();
     filterLemmas(lemmas, wlems);
   }
   if (options::nlExtTfTangentPlanes())
@@ -1258,1187 +809,6 @@ void NonlinearExtension::presolve()
   Trace("nl-ext") << "NonlinearExtension::presolve" << std::endl;
 }
 
-void NonlinearExtension::assignOrderIds(std::vector<Node>& vars,
-                                        NodeMultiset& order,
-                                        bool isConcrete,
-                                        bool isAbsolute)
-{
-  SortNlModel smv;
-  smv.d_nlm = &d_model;
-  smv.d_isConcrete = isConcrete;
-  smv.d_isAbsolute = isAbsolute;
-  smv.d_reverse_order = false;
-  std::sort(vars.begin(), vars.end(), smv);
-
-  order.clear();
-  // assign ordering id's
-  unsigned counter = 0;
-  unsigned order_index = isConcrete ? 0 : 1;
-  Node prev;
-  for (unsigned j = 0; j < vars.size(); j++) {
-    Node x = vars[j];
-    Node v = d_model.computeModelValue(x, isConcrete);
-    if( !v.isConst() ){
-      Trace("nl-ext-mvo") << "..do not assign order to " << x << " : " << v << std::endl;
-      //don't assign for non-constant values (transcendental function apps)
-      break;
-    }
-    Trace("nl-ext-mvo") << "  order " << x << " : " << v << std::endl;
-    if (v != prev) {
-      // builtin points
-      bool success;
-      do {
-        success = false;
-        if (order_index < d_order_points.size()) {
-          Node vv = d_model.computeModelValue(d_order_points[order_index],
-                                              isConcrete);
-          if (d_model.compareValue(v, vv, isAbsolute) <= 0)
-          {
-            counter++;
-            Trace("nl-ext-mvo") << "O[" << d_order_points[order_index]
-                                << "] = " << counter << std::endl;
-            order[d_order_points[order_index]] = counter;
-            prev = vv;
-            order_index++;
-            success = true;
-          }
-        }
-      } while (success);
-    }
-    if (prev.isNull() || d_model.compareValue(v, prev, isAbsolute) != 0)
-    {
-      counter++;
-    }
-    Trace("nl-ext-mvo") << "O[" << x << "] = " << counter << std::endl;
-    order[x] = counter;
-    prev = v;
-  }
-  while (order_index < d_order_points.size()) {
-    counter++;
-    Trace("nl-ext-mvo") << "O[" << d_order_points[order_index]
-                        << "] = " << counter << std::endl;
-    order[d_order_points[order_index]] = counter;
-    order_index++;
-  }
-}
-
-bool NonlinearExtension::getApproximateSqrt(Node c,
-                                            Node& l,
-                                            Node& u,
-                                            unsigned iter) const
-{
-  Assert(c.isConst());
-  if (c == d_one || c == d_zero)
-  {
-    l = c;
-    u = c;
-    return true;
-  }
-  Rational rc = c.getConst<Rational>();
-
-  Rational rl = rc < Rational(1) ? rc : Rational(1);
-  Rational ru = rc < Rational(1) ? Rational(1) : rc;
-  unsigned count = 0;
-  Rational half = Rational(1) / Rational(2);
-  while (count < iter)
-  {
-    Rational curr = half * (rl + ru);
-    Rational curr_sq = curr * curr;
-    if (curr_sq == rc)
-    {
-      rl = curr_sq;
-      ru = curr_sq;
-      break;
-    }
-    else if (curr_sq < rc)
-    {
-      rl = curr;
-    }
-    else
-    {
-      ru = curr;
-    }
-    count++;
-  }
-
-  NodeManager* nm = NodeManager::currentNM();
-  l = nm->mkConst(rl);
-  u = nm->mkConst(ru);
-  return true;
-}
-
-// show a <> 0 by inequalities between variables in monomial a w.r.t 0
-int NonlinearExtension::compareSign(Node oa, Node a, unsigned a_index,
-                                    int status, std::vector<Node>& exp,
-                                    std::vector<Node>& lem) {
-  Trace("nl-ext-debug") << "Process " << a << " at index " << a_index
-                        << ", status is " << status << std::endl;
-  Node mvaoa = d_model.computeAbstractModelValue(oa);
-  if (a_index == d_m_vlist[a].size()) {
-    if (mvaoa.getConst<Rational>().sgn() != status)
-    {
-      Node lemma =
-          safeConstructNary(AND, exp).impNode(mkLit(oa, d_zero, status * 2));
-      lem.push_back(lemma);
-    }
-    return status;
-  } else {
-    Assert(a_index < d_m_vlist[a].size());
-    Node av = d_m_vlist[a][a_index];
-    unsigned aexp = d_m_exp[a][av];
-    // take current sign in model
-    Node mvaav = d_model.computeAbstractModelValue(av);
-    int sgn = mvaav.getConst<Rational>().sgn();
-    Trace("nl-ext-debug") << "Process var " << av << "^" << aexp
-                          << ", model sign = " << sgn << std::endl;
-    if (sgn == 0) {
-      if (mvaoa.getConst<Rational>().sgn() != 0)
-      {
-        Node lemma = av.eqNode(d_zero).impNode(oa.eqNode(d_zero));
-        lem.push_back(lemma);
-      }
-      return 0;
-    } else {
-      if (aexp % 2 == 0) {
-        exp.push_back(av.eqNode(d_zero).negate());
-        return compareSign(oa, a, a_index + 1, status, exp, lem);
-      } else {
-        exp.push_back(
-            NodeManager::currentNM()->mkNode(sgn == 1 ? GT : LT, av, d_zero));
-        return compareSign(oa, a, a_index + 1, status * sgn, exp, lem);
-      }
-    }
-  }
-}
-
-bool NonlinearExtension::compareMonomial(
-    Node oa, Node a, NodeMultiset& a_exp_proc, Node ob, Node b,
-    NodeMultiset& b_exp_proc, std::vector<Node>& exp, std::vector<Node>& lem,
-    std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers) {
-  Trace("nl-ext-comp-debug")
-      << "Check |" << a << "| >= |" << b << "|" << std::endl;
-  unsigned pexp_size = exp.size();
-  if (compareMonomial(oa, a, 0, a_exp_proc, ob, b, 0, b_exp_proc, 0, exp, lem,
-                      cmp_infers)) {
-    return true;
-  } else {
-    exp.resize(pexp_size);
-    Trace("nl-ext-comp-debug")
-        << "Check |" << b << "| >= |" << a << "|" << std::endl;
-    if (compareMonomial(ob, b, 0, b_exp_proc, oa, a, 0, a_exp_proc, 0, exp, lem,
-                        cmp_infers)) {
-      return true;
-    } else {
-      return false;
-    }
-  }
-}
-
-bool NonlinearExtension::cmp_holds(
-    Node x, Node y, std::map<Node, std::map<Node, Node> >& cmp_infers,
-    std::vector<Node>& exp, std::map<Node, bool>& visited) {
-  if (x == y) {
-    return true;
-  } else if (visited.find(x) == visited.end()) {
-    visited[x] = true;
-    std::map<Node, std::map<Node, Node> >::iterator it = cmp_infers.find(x);
-    if (it != cmp_infers.end()) {
-      for (std::map<Node, Node>::iterator itc = it->second.begin();
-           itc != it->second.end(); ++itc) {
-        exp.push_back(itc->second);
-        if (cmp_holds(itc->first, y, cmp_infers, exp, visited)) {
-          return true;
-        }
-        exp.pop_back();
-      }
-    }
-  }
-  return false;
-}
-
-// trying to show a ( >, = ) b   by inequalities between variables in
-// monomials a,b
-bool NonlinearExtension::compareMonomial(
-    Node oa, Node a, unsigned a_index, NodeMultiset& a_exp_proc, Node ob,
-    Node b, unsigned b_index, NodeMultiset& b_exp_proc, int status,
-    std::vector<Node>& exp, std::vector<Node>& lem,
-    std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers) {
-  Trace("nl-ext-comp-debug")
-      << "compareMonomial " << oa << " and " << ob << ", indices = " << a_index
-      << " " << b_index << std::endl;
-  Assert(status == 0 || status == 2);
-  if (a_index == d_m_vlist[a].size() && b_index == d_m_vlist[b].size()) {
-    // finished, compare absolute value of abstract model values
-    int modelStatus = d_model.compare(oa, ob, false, true) * -2;
-    Trace("nl-ext-comp") << "...finished comparison with " << oa << " <"
-                         << status << "> " << ob
-                         << ", model status = " << modelStatus << std::endl;
-    if (status != modelStatus) {
-      Trace("nl-ext-comp-infer")
-          << "infer : " << oa << " <" << status << "> " << ob << std::endl;
-      if (status == 2) {
-        // must state that all variables are non-zero
-        for (unsigned j = 0; j < d_m_vlist[a].size(); j++) {
-          exp.push_back(d_m_vlist[a][j].eqNode(d_zero).negate());
-        }
-      }
-      Node clem = NodeManager::currentNM()->mkNode(
-          IMPLIES, safeConstructNary(AND, exp), mkLit(oa, ob, status, true));
-      Trace("nl-ext-comp-lemma") << "comparison lemma : " << clem << std::endl;
-      lem.push_back(clem);
-      cmp_infers[status][oa][ob] = clem;
-    }
-    return true;
-  } else {
-    // get a/b variable information
-    Node av;
-    unsigned aexp = 0;
-    unsigned avo = 0;
-    if (a_index < d_m_vlist[a].size()) {
-      av = d_m_vlist[a][a_index];
-      Assert(a_exp_proc[av] <= d_m_exp[a][av]);
-      aexp = d_m_exp[a][av] - a_exp_proc[av];
-      if (aexp == 0) {
-        return compareMonomial(oa, a, a_index + 1, a_exp_proc, ob, b, b_index,
-                               b_exp_proc, status, exp, lem, cmp_infers);
-      }
-      Assert(d_order_vars.find(av) != d_order_vars.end());
-      avo = d_order_vars[av];
-    }
-    Node bv;
-    unsigned bexp = 0;
-    unsigned bvo = 0;
-    if (b_index < d_m_vlist[b].size()) {
-      bv = d_m_vlist[b][b_index];
-      Assert(b_exp_proc[bv] <= d_m_exp[b][bv]);
-      bexp = d_m_exp[b][bv] - b_exp_proc[bv];
-      if (bexp == 0) {
-        return compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index + 1,
-                               b_exp_proc, status, exp, lem, cmp_infers);
-      }
-      Assert(d_order_vars.find(bv) != d_order_vars.end());
-      bvo = d_order_vars[bv];
-    }
-    // get "one" information
-    Assert(d_order_vars.find(d_one) != d_order_vars.end());
-    unsigned ovo = d_order_vars[d_one];
-    Trace("nl-ext-comp-debug") << "....vars : " << av << "^" << aexp << " "
-                               << bv << "^" << bexp << std::endl;
-
-    //--- cases
-    if (av.isNull()) {
-      if (bvo <= ovo) {
-        Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl;
-        // can multiply b by <=1
-        exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true));
-        return compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index + 1,
-                               b_exp_proc, bvo == ovo ? status : 2, exp, lem,
-                               cmp_infers);
-      } else {
-        Trace("nl-ext-comp-debug")
-            << "...failure, unmatched |b|>1 component." << std::endl;
-        return false;
-      }
-    } else if (bv.isNull()) {
-      if (avo >= ovo) {
-        Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl;
-        // can multiply a by >=1
-        exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true));
-        return compareMonomial(oa, a, a_index + 1, a_exp_proc, ob, b, b_index,
-                               b_exp_proc, avo == ovo ? status : 2, exp, lem,
-                               cmp_infers);
-      } else {
-        Trace("nl-ext-comp-debug")
-            << "...failure, unmatched |a|<1 component." << std::endl;
-        return false;
-      }
-    } else {
-      Assert(!av.isNull() && !bv.isNull());
-      if (avo >= bvo) {
-        if (bvo < ovo && avo >= ovo) {
-          Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl;
-          // do avo>=1 instead
-          exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true));
-          return compareMonomial(oa, a, a_index + 1, a_exp_proc, ob, b, b_index,
-                                 b_exp_proc, avo == ovo ? status : 2, exp, lem,
-                                 cmp_infers);
-        } else {
-          unsigned min_exp = aexp > bexp ? bexp : aexp;
-          a_exp_proc[av] += min_exp;
-          b_exp_proc[bv] += min_exp;
-          Trace("nl-ext-comp-debug")
-              << "...take leading " << min_exp << " from " << av << " and "
-              << bv << std::endl;
-          exp.push_back(mkLit(av, bv, avo == bvo ? 0 : 2, true));
-          bool ret = compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index,
-                                     b_exp_proc, avo == bvo ? status : 2, exp,
-                                     lem, cmp_infers);
-          a_exp_proc[av] -= min_exp;
-          b_exp_proc[bv] -= min_exp;
-          return ret;
-        }
-      } else {
-        if (bvo <= ovo) {
-          Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl;
-          // try multiply b <= 1
-          exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true));
-          return compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index + 1,
-                                 b_exp_proc, bvo == ovo ? status : 2, exp, lem,
-                                 cmp_infers);
-        } else {
-          Trace("nl-ext-comp-debug")
-              << "...failure, leading |b|>|a|>1 component." << std::endl;
-          return false;
-        }
-      }
-    }
-  }
-}
-
-// status 0 : n equal, -1 : n superset, 1 : n subset
-void NonlinearExtension::MonomialIndex::addTerm(Node n,
-                                                const std::vector<Node>& reps,
-                                                NonlinearExtension* nla,
-                                                int status, unsigned argIndex) {
-  if (status == 0) {
-    if (argIndex == reps.size()) {
-      d_monos.push_back(n);
-    } else {
-      d_data[reps[argIndex]].addTerm(n, reps, nla, status, argIndex + 1);
-    }
-  }
-  for (std::map<Node, MonomialIndex>::iterator it = d_data.begin();
-       it != d_data.end(); ++it) {
-    if (status != 0 || argIndex == reps.size() || it->first != reps[argIndex]) {
-      // if we do not contain this variable, then if we were a superset,
-      // fail (-2), otherwise we are subset.  if we do contain this
-      // variable, then if we were equal, we are superset since variables
-      // are ordered, otherwise we remain the same.
-      int new_status =
-          std::find(reps.begin(), reps.end(), it->first) == reps.end()
-              ? (status >= 0 ? 1 : -2)
-              : (status == 0 ? -1 : status);
-      if (new_status != -2) {
-        it->second.addTerm(n, reps, nla, new_status, argIndex);
-      }
-    }
-  }
-  // compare for subsets
-  for (unsigned i = 0; i < d_monos.size(); i++) {
-    Node m = d_monos[i];
-    if (m != n) {
-      // we are superset if we are equal and haven't traversed all variables
-      int cstatus = status == 0 ? (argIndex == reps.size() ? 0 : -1) : status;
-      Trace("nl-ext-mindex-debug") << "  compare " << n << " and " << m
-                                   << ", status = " << cstatus << std::endl;
-      if (cstatus <= 0 && nla->isMonomialSubset(m, n)) {
-        nla->registerMonomialSubset(m, n);
-        Trace("nl-ext-mindex-debug") << "...success" << std::endl;
-      } else if (cstatus >= 0 && nla->isMonomialSubset(n, m)) {
-        nla->registerMonomialSubset(n, m);
-        Trace("nl-ext-mindex-debug") << "...success (rev)" << std::endl;
-      }
-    }
-  }
-}
-
-std::vector<Node> NonlinearExtension::checkMonomialSign() {
-  std::vector<Node> lemmas;
-  std::map<Node, int> signs;
-  Trace("nl-ext") << "Get monomial sign lemmas..." << std::endl;
-  for (unsigned j = 0; j < d_ms.size(); j++) {
-    Node a = d_ms[j];
-    if (d_ms_proc.find(a) == d_ms_proc.end()) {
-      std::vector<Node> exp;
-      if (Trace.isOn("nl-ext-debug"))
-      {
-        Node cmva = d_model.computeConcreteModelValue(a);
-        Trace("nl-ext-debug")
-            << "  process " << a << ", mv=" << cmva << "..." << std::endl;
-      }
-      if( d_m_nconst_factor.find( a )==d_m_nconst_factor.end() ){
-        signs[a] = compareSign(a, a, 0, 1, exp, lemmas);
-        if (signs[a] == 0) {
-          d_ms_proc[a] = true;
-          Trace("nl-ext-debug") << "...mark " << a
-                             << " reduced since its value is 0." << std::endl;
-        }
-      }else{
-        Trace("nl-ext-debug") << "...can't conclude sign lemma for " << a << " since model value of a factor is non-constant." << std::endl;
-        //TODO
-      }
-    }
-  }
-  return lemmas;
-}
-
-std::vector<Node> NonlinearExtension::checkMonomialMagnitude( unsigned c ) {
-  unsigned r = 1;
-  std::vector<Node> lemmas;
-// if (x,y,L) in cmp_infers, then x > y inferred as conclusion of L
-  // in lemmas
-  std::map<int, std::map<Node, std::map<Node, Node> > > cmp_infers;
-  Trace("nl-ext") << "Get monomial comparison lemmas (order=" << r
-                  << ", compare=" << c << ")..." << std::endl;
-  for (unsigned j = 0; j < d_ms.size(); j++) {
-    Node a = d_ms[j];
-    if (d_ms_proc.find(a) == d_ms_proc.end() && 
-        d_m_nconst_factor.find( a )==d_m_nconst_factor.end()) {
-      if (c == 0) {
-        // compare magnitude against 1
-        std::vector<Node> exp;
-        NodeMultiset a_exp_proc;
-        NodeMultiset b_exp_proc;
-        compareMonomial(a, a, a_exp_proc, d_one, d_one, b_exp_proc, exp,
-                        lemmas, cmp_infers);
-      } else {
-        std::map<Node, NodeMultiset>::iterator itmea = d_m_exp.find(a);
-        Assert(itmea != d_m_exp.end());
-        if (c == 1) {
-          // TODO : not just against containing variables?
-          // compare magnitude against variables
-          for (unsigned k = 0; k < d_ms_vars.size(); k++) {
-            Node v = d_ms_vars[k];
-            Node mvcv = d_model.computeConcreteModelValue(v);
-            if (mvcv.isConst())
-            {
-              std::vector<Node> exp;
-              NodeMultiset a_exp_proc;
-              NodeMultiset b_exp_proc;
-              if (itmea->second.find(v) != itmea->second.end()) {
-                a_exp_proc[v] = 1;
-                b_exp_proc[v] = 1;
-                setMonomialFactor(a, v, a_exp_proc);
-                setMonomialFactor(v, a, b_exp_proc);
-                compareMonomial(a, a, a_exp_proc, v, v, b_exp_proc, exp,
-                                lemmas, cmp_infers);
-              }
-            }
-          }
-        } else {
-          // compare magnitude against other non-linear monomials
-          for (unsigned k = (j + 1); k < d_ms.size(); k++) {
-            Node b = d_ms[k];
-            //(signs[a]==signs[b])==(r==0)
-            if (d_ms_proc.find(b) == d_ms_proc.end() && 
-                d_m_nconst_factor.find( b )==d_m_nconst_factor.end()) {
-              std::map<Node, NodeMultiset>::iterator itmeb =
-                  d_m_exp.find(b);
-              Assert(itmeb != d_m_exp.end());
-
-              std::vector<Node> exp;
-              // take common factors of monomials, set minimum of
-              // common exponents as processed
-              NodeMultiset a_exp_proc;
-              NodeMultiset b_exp_proc;
-              for (NodeMultiset::iterator itmea2 = itmea->second.begin();
-                   itmea2 != itmea->second.end(); ++itmea2) {
-                NodeMultiset::iterator itmeb2 =
-                    itmeb->second.find(itmea2->first);
-                if (itmeb2 != itmeb->second.end()) {
-                  unsigned min_exp = itmea2->second > itmeb2->second
-                                         ? itmeb2->second
-                                         : itmea2->second;
-                  a_exp_proc[itmea2->first] = min_exp;
-                  b_exp_proc[itmea2->first] = min_exp;
-                  Trace("nl-ext-comp")
-                      << "Common exponent : " << itmea2->first << " : "
-                      << min_exp << std::endl;
-                }
-              }
-              if (!a_exp_proc.empty()) {
-                setMonomialFactor(a, b, a_exp_proc);
-                setMonomialFactor(b, a, b_exp_proc);
-              }
-              /*
-              if( !a_exp_proc.empty() ){
-                //reduction based on common exponents a > 0 => ( a * b
-              <> a * c <=> b <> c ), a < 0 => ( a * b <> a * c <=> b
-              !<> c )  ? }else{ compareMonomial( a, a, a_exp_proc, b,
-              b, b_exp_proc, exp, lemmas );
-              }
-              */
-              compareMonomial(a, a, a_exp_proc, b, b, b_exp_proc, exp,
-                              lemmas, cmp_infers);
-            }
-          }
-        }
-      }
-    }
-  }
-  // remove redundant lemmas, e.g. if a > b, b > c, a > c were
-  // inferred, discard lemma with conclusion a > c
-  Trace("nl-ext-comp") << "Compute redundancies for " << lemmas.size()
-                       << " lemmas." << std::endl;
-  // naive
-  std::vector<Node> r_lemmas;
-  for (std::map<int, std::map<Node, std::map<Node, Node> > >::iterator itb =
-           cmp_infers.begin();
-       itb != cmp_infers.end(); ++itb) {
-    for (std::map<Node, std::map<Node, Node> >::iterator itc =
-             itb->second.begin();
-         itc != itb->second.end(); ++itc) {
-      for (std::map<Node, Node>::iterator itc2 = itc->second.begin();
-           itc2 != itc->second.end(); ++itc2) {
-        std::map<Node, bool> visited;
-        for (std::map<Node, Node>::iterator itc3 = itc->second.begin();
-             itc3 != itc->second.end(); ++itc3) {
-          if (itc3->first != itc2->first) {
-            std::vector<Node> exp;
-            if (cmp_holds(itc3->first, itc2->first, itb->second, exp,
-                          visited)) {
-              r_lemmas.push_back(itc2->second);
-              Trace("nl-ext-comp")
-                  << "...inference of " << itc->first << " > "
-                  << itc2->first << " was redundant." << std::endl;
-              break;
-            }
-          }
-        }
-      }
-    }
-  }
-  std::vector<Node> nr_lemmas;
-  for (unsigned i = 0; i < lemmas.size(); i++) {
-    if (std::find(r_lemmas.begin(), r_lemmas.end(), lemmas[i]) ==
-        r_lemmas.end()) {
-      nr_lemmas.push_back(lemmas[i]);
-    }
-  }
-  // TODO: only take maximal lower/minimial lower bounds?
-
-  Trace("nl-ext-comp") << nr_lemmas.size() << " / " << lemmas.size()
-                       << " were non-redundant." << std::endl;
-  return nr_lemmas;
-}
-
-std::vector<Node> NonlinearExtension::checkTangentPlanes() {
-  std::vector< Node > lemmas;
-  Trace("nl-ext") << "Get monomial tangent plane lemmas..." << std::endl;
-  NodeManager* nm = NodeManager::currentNM();
-  unsigned kstart = d_ms_vars.size();
-  for (unsigned k = kstart; k < d_mterms.size(); k++) {
-    Node t = d_mterms[k];
-    // if this term requires a refinement
-    if (d_tplane_refine.find(t) != d_tplane_refine.end())
-    {
-      Trace("nl-ext-tplanes")
-          << "Look at monomial requiring refinement : " << t << std::endl;
-      // get a decomposition
-      std::map<Node, std::vector<Node> >::iterator it =
-          d_m_contain_children.find(t);
-      if (it != d_m_contain_children.end()) {
-        std::map<Node, std::map<Node, bool> > dproc;
-        for (unsigned j = 0; j < it->second.size(); j++) {
-          Node tc = it->second[j];
-          if (tc != d_one) {
-            Node tc_diff = d_m_contain_umult[tc][t];
-            Assert(!tc_diff.isNull());
-            Node a = tc < tc_diff ? tc : tc_diff;
-            Node b = tc < tc_diff ? tc_diff : tc;
-            if (dproc[a].find(b) == dproc[a].end()) {
-              dproc[a][b] = true;
-              Trace("nl-ext-tplanes")
-                  << "  decomposable into : " << a << " * " << b << std::endl;
-              Node a_v_c = d_model.computeAbstractModelValue(a);
-              Node b_v_c = d_model.computeAbstractModelValue(b);
-              // points we will add tangent planes for
-              std::vector< Node > pts[2];
-              pts[0].push_back( a_v_c );
-              pts[1].push_back( b_v_c );
-              // if previously refined
-              bool prevRefine = d_tangent_val_bound[0][a].find( b )!=d_tangent_val_bound[0][a].end();
-              // a_min, a_max, b_min, b_max
-              for( unsigned p=0; p<4; p++ ){
-                Node curr_v = p<=1 ? a_v_c : b_v_c;
-                if( prevRefine ){
-                  Node pt_v = d_tangent_val_bound[p][a][b];
-                  Assert(!pt_v.isNull());
-                  if( curr_v!=pt_v ){
-                    Node do_extend =
-                        nm->mkNode((p == 1 || p == 3) ? GT : LT, curr_v, pt_v);
-                    do_extend = Rewriter::rewrite( do_extend );
-                    if( do_extend==d_true ){
-                      for( unsigned q=0; q<2; q++ ){
-                        pts[ p<=1 ? 0 : 1 ].push_back( curr_v );
-                        pts[ p<=1 ? 1 : 0 ].push_back( d_tangent_val_bound[ p<=1 ? 2+q : q ][a][b] );
-                      }
-                    }
-                  }
-                }else{
-                  d_tangent_val_bound[p][a][b] = curr_v;
-                }
-              }
-              
-              for( unsigned p=0; p<pts[0].size(); p++ ){
-                Node a_v = pts[0][p];
-                Node b_v = pts[1][p];
-              
-                // tangent plane
-                Node tplane = nm->mkNode(MINUS,
-                                         nm->mkNode(PLUS,
-                                                    nm->mkNode(MULT, b_v, a),
-                                                    nm->mkNode(MULT, a_v, b)),
-                                         nm->mkNode(MULT, a_v, b_v));
-                for (unsigned d = 0; d < 4; d++) {
-                  Node aa = nm->mkNode(d == 0 || d == 3 ? GEQ : LEQ, a, a_v);
-                  Node ab = nm->mkNode(d == 1 || d == 3 ? GEQ : LEQ, b, b_v);
-                  Node conc = nm->mkNode(d <= 1 ? LEQ : GEQ, t, tplane);
-                  Node tlem = nm->mkNode(OR, aa.negate(), ab.negate(), conc);
-                  Trace("nl-ext-tplanes")
-                      << "Tangent plane lemma : " << tlem << std::endl;
-                  lemmas.push_back(tlem);
-                }
-
-                // tangent plane reverse implication
-
-                // t <= tplane -> ( (a <= a_v ^ b >= b_v) v
-                // (a >= a_v ^ b <= b_v) ).
-                // in clause form, the above becomes
-                // t <= tplane -> a <= a_v v b <= b_v.
-                // t <= tplane -> b >= b_v v a >= a_v.
-                Node a_leq_av = nm->mkNode(LEQ, a, a_v);
-                Node b_leq_bv = nm->mkNode(LEQ, b, b_v);
-                Node a_geq_av = nm->mkNode(GEQ, a, a_v);
-                Node b_geq_bv = nm->mkNode(GEQ, b, b_v);
-
-                Node t_leq_tplane = nm->mkNode(LEQ, t, tplane);
-                Node a_leq_av_or_b_leq_bv = nm->mkNode(OR, a_leq_av, b_leq_bv);
-                Node b_geq_bv_or_a_geq_av = nm->mkNode(OR, b_geq_bv, a_geq_av);
-                Node ub_reverse1 =
-                    nm->mkNode(OR, t_leq_tplane.negate(), a_leq_av_or_b_leq_bv);
-                Trace("nl-ext-tplanes")
-                    << "Tangent plane lemma (reverse) : " << ub_reverse1
-                    << std::endl;
-                lemmas.push_back(ub_reverse1);
-                Node ub_reverse2 =
-                    nm->mkNode(OR, t_leq_tplane.negate(), b_geq_bv_or_a_geq_av);
-                Trace("nl-ext-tplanes")
-                    << "Tangent plane lemma (reverse) : " << ub_reverse2
-                    << std::endl;
-                lemmas.push_back(ub_reverse2);
-
-                // t >= tplane -> ( (a <= a_v ^ b <= b_v) v
-                // (a >= a_v ^ b >= b_v) ).
-                // in clause form, the above becomes
-                // t >= tplane -> a <= a_v v b >= b_v.
-                // t >= tplane -> b >= b_v v a <= a_v
-                Node t_geq_tplane = nm->mkNode(GEQ, t, tplane);
-                Node a_leq_av_or_b_geq_bv = nm->mkNode(OR, a_leq_av, b_geq_bv);
-                Node a_geq_av_or_b_leq_bv = nm->mkNode(OR, a_geq_av, b_leq_bv);
-                Node lb_reverse1 =
-                    nm->mkNode(OR, t_geq_tplane.negate(), a_leq_av_or_b_geq_bv);
-                Trace("nl-ext-tplanes")
-                    << "Tangent plane lemma (reverse) : " << lb_reverse1
-                    << std::endl;
-                lemmas.push_back(lb_reverse1);
-                Node lb_reverse2 =
-                    nm->mkNode(OR, t_geq_tplane.negate(), a_geq_av_or_b_leq_bv);
-                Trace("nl-ext-tplanes")
-                    << "Tangent plane lemma (reverse) : " << lb_reverse2
-                    << std::endl;
-                lemmas.push_back(lb_reverse2);
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-  Trace("nl-ext") << "...trying " << lemmas.size()
-                  << " tangent plane lemmas..." << std::endl;
-  return lemmas;
-}
-
-std::vector<Node> NonlinearExtension::checkMonomialInferBounds(
-    std::vector<Node>& nt_lemmas,
-    const std::vector<Node>& asserts,
-    const std::vector<Node>& false_asserts)
-{
-  std::vector< Node > lemmas; 
-  // register constraints
-  Trace("nl-ext-debug") << "Register bound constraints..." << std::endl;
-  for (const Node& lit : asserts)
-  {
-    bool polarity = lit.getKind() != NOT;
-    Node atom = lit.getKind() == NOT ? lit[0] : lit;
-    registerConstraint(atom);
-    bool is_false_lit =
-        std::find(false_asserts.begin(), false_asserts.end(), lit)
-        != false_asserts.end();
-    // add information about bounds to variables
-    std::map<Node, std::map<Node, ConstraintInfo> >::iterator itc =
-        d_c_info.find(atom);
-    std::map<Node, std::map<Node, bool> >::iterator itcm =
-        d_c_info_maxm.find(atom);
-    if (itc != d_c_info.end()) {
-      Assert(itcm != d_c_info_maxm.end());
-      for (std::map<Node, ConstraintInfo>::iterator itcc = itc->second.begin();
-           itcc != itc->second.end(); ++itcc) {
-        Node x = itcc->first;
-        Node coeff = itcc->second.d_coeff;
-        Node rhs = itcc->second.d_rhs;
-        Kind type = itcc->second.d_type;
-        Node exp = lit;
-        if (!polarity) {
-          // reverse
-          if (type == EQUAL)
-          {
-            // we will take the strict inequality in the direction of the
-            // model
-            Node lhs = ArithMSum::mkCoeffTerm(coeff, x);
-            Node query = NodeManager::currentNM()->mkNode(GT, lhs, rhs);
-            Node query_mv = d_model.computeAbstractModelValue(query);
-            if (query_mv == d_true) {
-              exp = query;
-              type = GT;
-            } else {
-              Assert(query_mv == d_false);
-              exp = NodeManager::currentNM()->mkNode(LT, lhs, rhs);
-              type = LT;
-            }
-          } else {
-            type = negateKind(type);
-          }
-        }
-        // add to status if maximal degree
-        d_ci_max[x][coeff][rhs] = itcm->second.find(x) != itcm->second.end();
-        if (Trace.isOn("nl-ext-bound-debug2")) {
-          Node t = ArithMSum::mkCoeffTerm(coeff, x);
-          Trace("nl-ext-bound-debug2")
-              << "Add Bound: " << t << " " << type << " " << rhs << " by "
-              << exp << std::endl;
-        }
-        bool updated = true;
-        std::map<Node, Kind>::iterator its = d_ci[x][coeff].find(rhs);
-        if (its == d_ci[x][coeff].end()) {
-          d_ci[x][coeff][rhs] = type;
-          d_ci_exp[x][coeff][rhs] = exp;
-        } else if (type != its->second) {
-          Trace("nl-ext-bound-debug2")
-              << "Joining kinds : " << type << " " << its->second << std::endl;
-          Kind jk = joinKinds(type, its->second);
-          if (jk == UNDEFINED_KIND)
-          {
-            updated = false;
-          } else if (jk != its->second) {
-            if (jk == type) {
-              d_ci[x][coeff][rhs] = type;
-              d_ci_exp[x][coeff][rhs] = exp;
-            } else {
-              d_ci[x][coeff][rhs] = jk;
-              d_ci_exp[x][coeff][rhs] = NodeManager::currentNM()->mkNode(
-                  AND, d_ci_exp[x][coeff][rhs], exp);
-            }
-          } else {
-            updated = false;
-          }
-        }
-        if (Trace.isOn("nl-ext-bound")) {
-          if (updated) {
-            Trace("nl-ext-bound") << "Bound: ";
-            debugPrintBound("nl-ext-bound", coeff, x, d_ci[x][coeff][rhs], rhs);
-            Trace("nl-ext-bound") << " by " << d_ci_exp[x][coeff][rhs];
-            if (d_ci_max[x][coeff][rhs]) {
-              Trace("nl-ext-bound") << ", is max degree";
-            }
-            Trace("nl-ext-bound") << std::endl;
-          }
-        }
-        // compute if bound is not satisfied, and store what is required
-        // for a possible refinement
-        if (options::nlExtTangentPlanes()) {
-          if (is_false_lit) {
-            d_tplane_refine.insert(x);
-          }
-        }
-      }
-    }
-  }
-  // reflexive constraints
-  Node null_coeff;
-  for (unsigned j = 0; j < d_mterms.size(); j++) {
-    Node n = d_mterms[j];
-    d_ci[n][null_coeff][n] = EQUAL;
-    d_ci_exp[n][null_coeff][n] = d_true;
-    d_ci_max[n][null_coeff][n] = false;
-  }
-
-  Trace("nl-ext") << "Get inferred bound lemmas..." << std::endl;
-  
-  for (unsigned k = 0; k < d_mterms.size(); k++) {
-    Node x = d_mterms[k];
-    Trace("nl-ext-bound-debug")
-        << "Process bounds for " << x << " : " << std::endl;
-    std::map<Node, std::vector<Node> >::iterator itm =
-        d_m_contain_parent.find(x);
-    if (itm != d_m_contain_parent.end()) {
-      Trace("nl-ext-bound-debug") << "...has " << itm->second.size()
-                                  << " parent monomials." << std::endl;
-      // check derived bounds
-      std::map<Node, std::map<Node, std::map<Node, Kind> > >::iterator itc =
-          d_ci.find(x);
-      if (itc != d_ci.end()) {
-        for (std::map<Node, std::map<Node, Kind> >::iterator itcc =
-                 itc->second.begin();
-             itcc != itc->second.end(); ++itcc) {
-          Node coeff = itcc->first;
-          Node t = ArithMSum::mkCoeffTerm(coeff, x);
-          for (std::map<Node, Kind>::iterator itcr = itcc->second.begin();
-               itcr != itcc->second.end(); ++itcr) {
-            Node rhs = itcr->first;
-            // only consider this bound if maximal degree
-            if (d_ci_max[x][coeff][rhs]) {
-              Kind type = itcr->second;
-              for (unsigned j = 0; j < itm->second.size(); j++) {
-                Node y = itm->second[j];
-                Assert(d_m_contain_mult[x].find(y)
-                       != d_m_contain_mult[x].end());
-                Node mult = d_m_contain_mult[x][y];
-                // x <k> t => m*x <k'> t  where y = m*x
-                // get the sign of mult
-                Node mmv = d_model.computeConcreteModelValue(mult);
-                Trace("nl-ext-bound-debug2")
-                    << "Model value of " << mult << " is " << mmv << std::endl;
-                if(mmv.isConst()){
-                  int mmv_sign = mmv.getConst<Rational>().sgn();
-                  Trace("nl-ext-bound-debug2")
-                      << "  sign of " << mmv << " is " << mmv_sign << std::endl;
-                  if (mmv_sign != 0) {
-                    Trace("nl-ext-bound-debug")
-                        << "  from " << x << " * " << mult << " = " << y
-                        << " and " << t << " " << type << " " << rhs
-                        << ", infer : " << std::endl;
-                    Kind infer_type =
-                        mmv_sign == -1 ? reverseRelationKind(type) : type;
-                    Node infer_lhs =
-                        NodeManager::currentNM()->mkNode(MULT, mult, t);
-                    Node infer_rhs =
-                        NodeManager::currentNM()->mkNode(MULT, mult, rhs);
-                    Node infer = NodeManager::currentNM()->mkNode(
-                        infer_type, infer_lhs, infer_rhs);
-                    Trace("nl-ext-bound-debug") << "     " << infer << std::endl;
-                    infer = Rewriter::rewrite(infer);
-                    Trace("nl-ext-bound-debug2")
-                        << "     ...rewritten : " << infer << std::endl;
-                    // check whether it is false in model for abstraction
-                    Node infer_mv = d_model.computeAbstractModelValue(infer);
-                    Trace("nl-ext-bound-debug")
-                        << "       ...infer model value is " << infer_mv
-                        << std::endl;
-                    if (infer_mv == d_false) {
-                      Node exp = NodeManager::currentNM()->mkNode(
-                          AND,
-                          NodeManager::currentNM()->mkNode(
-                              mmv_sign == 1 ? GT : LT, mult, d_zero),
-                          d_ci_exp[x][coeff][rhs]);
-                      Node iblem =
-                          NodeManager::currentNM()->mkNode(IMPLIES, exp, infer);
-                      Node pr_iblem = iblem;
-                      iblem = Rewriter::rewrite(iblem);
-                      bool introNewTerms = hasNewMonomials(iblem, d_ms);
-                      Trace("nl-ext-bound-lemma")
-                          << "*** Bound inference lemma : " << iblem
-                          << " (pre-rewrite : " << pr_iblem << ")" << std::endl;
-                      // Trace("nl-ext-bound-lemma") << "       intro new
-                      // monomials = " << introNewTerms << std::endl;
-                      if (!introNewTerms) {
-                        lemmas.push_back(iblem);
-                      } else {
-                        nt_lemmas.push_back(iblem);
-                      }
-                    }
-                  } else {
-                    Trace("nl-ext-bound-debug") << "     ...coefficient " << mult
-                                                << " is zero." << std::endl;
-                  }
-                }else{
-                  Trace("nl-ext-bound-debug") << "     ...coefficient " << mult
-                                              << " is non-constant (probably transcendental)." << std::endl;
-                }
-              }
-            }
-          }
-        }
-      }
-    } else {
-      Trace("nl-ext-bound-debug") << "...has no parent monomials." << std::endl;
-    }
-  }
-  return lemmas;
-}
-
-std::vector<Node> NonlinearExtension::checkFactoring(
-    const std::vector<Node>& asserts, const std::vector<Node>& false_asserts)
-{
-  std::vector< Node > lemmas; 
-  Trace("nl-ext") << "Get factoring lemmas..." << std::endl;
-  for (const Node& lit : asserts)
-  {
-    bool polarity = lit.getKind() != NOT;
-    Node atom = lit.getKind() == NOT ? lit[0] : lit;
-    Node litv = d_model.computeConcreteModelValue(lit);
-    bool considerLit = false;
-    // Only consider literals that are in false_asserts.
-    considerLit = std::find(false_asserts.begin(), false_asserts.end(), lit)
-                  != false_asserts.end();
-
-    if (considerLit)
-    {
-      std::map<Node, Node> msum;
-      if (ArithMSum::getMonomialSumLit(atom, msum))
-      {
-        Trace("nl-ext-factor") << "Factoring for literal " << lit << ", monomial sum is : " << std::endl;
-        if (Trace.isOn("nl-ext-factor")) {
-          ArithMSum::debugPrintMonomialSum(msum, "nl-ext-factor");
-        }
-        std::map< Node, std::vector< Node > > factor_to_mono;
-        std::map< Node, std::vector< Node > > factor_to_mono_orig;
-        for( std::map<Node, Node>::iterator itm = msum.begin(); itm != msum.end(); ++itm ){
-          if( !itm->first.isNull() ){
-            if( itm->first.getKind()==NONLINEAR_MULT ){
-              std::vector< Node > children;
-              for( unsigned i=0; i<itm->first.getNumChildren(); i++ ){
-                children.push_back( itm->first[i] );
-              }
-              std::map< Node, bool > processed;
-              for( unsigned i=0; i<itm->first.getNumChildren(); i++ ){
-                if( processed.find( itm->first[i] )==processed.end() ){
-                  processed[itm->first[i]] = true;
-                  children[i] = d_one;
-                  if( !itm->second.isNull() ){
-                    children.push_back( itm->second );
-                  }
-                  Node val = NodeManager::currentNM()->mkNode(MULT, children);
-                  if( !itm->second.isNull() ){
-                    children.pop_back();
-                  }
-                  children[i] = itm->first[i];
-                  val = Rewriter::rewrite( val );
-                  factor_to_mono[itm->first[i]].push_back( val );
-                  factor_to_mono_orig[itm->first[i]].push_back( itm->first );
-                }
-              }
-            }
-          }
-        }
-        for( std::map< Node, std::vector< Node > >::iterator itf = factor_to_mono.begin(); itf != factor_to_mono.end(); ++itf ){
-          Node x = itf->first;
-          if (itf->second.size() == 1)
-          {
-            std::map<Node, Node>::iterator itm = msum.find(x);
-            if (itm != msum.end())
-            {
-              itf->second.push_back(itm->second.isNull() ? d_one : itm->second);
-              factor_to_mono_orig[x].push_back(x);
-            }
-          }
-          if( itf->second.size()>1 ){
-            Node sum = NodeManager::currentNM()->mkNode(PLUS, itf->second);
-            sum = Rewriter::rewrite( sum );
-            Trace("nl-ext-factor")
-                << "* Factored sum for " << x << " : " << sum << std::endl;
-            Node kf = getFactorSkolem(sum, lemmas);
-            std::vector< Node > poly;
-            poly.push_back(NodeManager::currentNM()->mkNode(MULT, x, kf));
-            std::map<Node, std::vector<Node> >::iterator itfo =
-                factor_to_mono_orig.find(x);
-            Assert(itfo != factor_to_mono_orig.end());
-            for( std::map<Node, Node>::iterator itm = msum.begin(); itm != msum.end(); ++itm ){
-              if( std::find( itfo->second.begin(), itfo->second.end(), itm->first )==itfo->second.end() ){
-                poly.push_back(ArithMSum::mkCoeffTerm(
-                    itm->second, itm->first.isNull() ? d_one : itm->first));
-              }
-            }
-            Node polyn = poly.size() == 1
-                             ? poly[0]
-                             : NodeManager::currentNM()->mkNode(PLUS, poly);
-            Trace("nl-ext-factor") << "...factored polynomial : " << polyn << std::endl;
-            Node conc_lit = NodeManager::currentNM()->mkNode( atom.getKind(), polyn, d_zero );
-            conc_lit = Rewriter::rewrite( conc_lit );
-            if( !polarity ){
-              conc_lit = conc_lit.negate();
-            }
-            
-            std::vector< Node > lemma_disj;
-            lemma_disj.push_back( lit.negate() );
-            lemma_disj.push_back( conc_lit );
-            Node flem = NodeManager::currentNM()->mkNode(OR, lemma_disj);
-            Trace("nl-ext-factor") << "...lemma is " << flem << std::endl;
-            lemmas.push_back( flem ); 
-          }
-        }
-      }
-    }
-  }
-  return lemmas;
-}
-
-Node NonlinearExtension::getFactorSkolem( Node n, std::vector< Node >& lemmas ) {
-  std::map< Node, Node >::iterator itf = d_factor_skolem.find( n );
-  if( itf==d_factor_skolem.end() ){
-    Node k = NodeManager::currentNM()->mkSkolem( "kf", n.getType() );
-    Node k_eq = Rewriter::rewrite( k.eqNode( n ) );
-    lemmas.push_back( k_eq );
-    d_factor_skolem[n] = k;
-    return k;
-  }else{
-    return itf->second;
-  }  
-}
-
-std::vector<Node> NonlinearExtension::checkMonomialInferResBounds() {            
-  std::vector< Node > lemmas; 
-  Trace("nl-ext") << "Get monomial resolution inferred bound lemmas..." << std::endl;
-  for (unsigned j = 0; j < d_mterms.size(); j++) {
-    Node a = d_mterms[j];
-    std::map<Node, std::map<Node, std::map<Node, Kind> > >::iterator itca =
-        d_ci.find(a);
-    if (itca != d_ci.end()) {
-      for (unsigned k = (j + 1); k < d_mterms.size(); k++) {
-        Node b = d_mterms[k];
-        std::map<Node, std::map<Node, std::map<Node, Kind> > >::iterator
-            itcb = d_ci.find(b);
-        if (itcb != d_ci.end()) {
-          Trace("nl-ext-rbound-debug") << "resolution inferences : compare "
-                                       << a << " and " << b << std::endl;
-          // if they have common factors
-          std::map<Node, Node>::iterator ita = d_mono_diff[a].find(b);
-          if (ita != d_mono_diff[a].end()) {
-            Trace("nl-ext-rbound") << "Get resolution inferences for [a] " << a
-                                   << " vs [b] " << b << std::endl;
-            std::map<Node, Node>::iterator itb = d_mono_diff[b].find(a);
-            Assert(itb != d_mono_diff[b].end());
-            Node mv_a = d_model.computeAbstractModelValue(ita->second);
-            Assert(mv_a.isConst());
-            int mv_a_sgn = mv_a.getConst<Rational>().sgn();
-            if (mv_a_sgn == 0)
-            {
-              // we don't compare monomials whose current model value is zero
-              continue;
-            }
-            Node mv_b = d_model.computeAbstractModelValue(itb->second);
-            Assert(mv_b.isConst());
-            int mv_b_sgn = mv_b.getConst<Rational>().sgn();
-            if (mv_b_sgn == 0)
-            {
-              // we don't compare monomials whose current model value is zero
-              continue;
-            }
-            Trace("nl-ext-rbound")
-                << "  [a] factor is " << ita->second
-                << ", sign in model = " << mv_a_sgn << std::endl;
-            Trace("nl-ext-rbound")
-                << "  [b] factor is " << itb->second
-                << ", sign in model = " << mv_b_sgn << std::endl;
-
-            std::vector<Node> exp;
-            // bounds of a
-            for (std::map<Node, std::map<Node, Kind> >::iterator itcac =
-                     itca->second.begin();
-                 itcac != itca->second.end(); ++itcac) {
-              Node coeff_a = itcac->first;
-              for (std::map<Node, Kind>::iterator itcar =
-                       itcac->second.begin();
-                   itcar != itcac->second.end(); ++itcar) {
-                Node rhs_a = itcar->first;
-                Node rhs_a_res_base =
-                    NodeManager::currentNM()->mkNode(MULT, itb->second, rhs_a);
-                rhs_a_res_base = Rewriter::rewrite(rhs_a_res_base);
-                if (!hasNewMonomials(rhs_a_res_base, d_ms)) {
-                  Kind type_a = itcar->second;
-                  exp.push_back(d_ci_exp[a][coeff_a][rhs_a]);
-
-                  // bounds of b
-                  for (std::map<Node, std::map<Node, Kind> >::iterator itcbc =
-                           itcb->second.begin();
-                       itcbc != itcb->second.end(); ++itcbc) {
-                    Node coeff_b = itcbc->first;
-                    Node rhs_a_res =
-                        ArithMSum::mkCoeffTerm(coeff_b, rhs_a_res_base);
-                    for (std::map<Node, Kind>::iterator itcbr =
-                             itcbc->second.begin();
-                         itcbr != itcbc->second.end(); ++itcbr) {
-                      Node rhs_b = itcbr->first;
-                      Node rhs_b_res = NodeManager::currentNM()->mkNode(
-                          MULT, ita->second, rhs_b);
-                      rhs_b_res = ArithMSum::mkCoeffTerm(coeff_a, rhs_b_res);
-                      rhs_b_res = Rewriter::rewrite(rhs_b_res);
-                      if (!hasNewMonomials(rhs_b_res, d_ms)) {
-                        Kind type_b = itcbr->second;
-                        exp.push_back(d_ci_exp[b][coeff_b][rhs_b]);
-                        if (Trace.isOn("nl-ext-rbound")) {
-                          Trace("nl-ext-rbound") << "* try bounds : ";
-                          debugPrintBound("nl-ext-rbound", coeff_a, a, type_a,
-                                          rhs_a);
-                          Trace("nl-ext-rbound") << std::endl;
-                          Trace("nl-ext-rbound") << "               ";
-                          debugPrintBound("nl-ext-rbound", coeff_b, b, type_b,
-                                          rhs_b);
-                          Trace("nl-ext-rbound") << std::endl;
-                        }
-                        Kind types[2];
-                        for (unsigned r = 0; r < 2; r++) {
-                          Node pivot_factor =
-                              r == 0 ? itb->second : ita->second;
-                          int pivot_factor_sign =
-                              r == 0 ? mv_b_sgn : mv_a_sgn;
-                          types[r] = r == 0 ? type_a : type_b;
-                          if (pivot_factor_sign == (r == 0 ? 1 : -1)) {
-                            types[r] = reverseRelationKind(types[r]);
-                          }
-                          if (pivot_factor_sign == 1) {
-                            exp.push_back(NodeManager::currentNM()->mkNode(
-                                GT, pivot_factor, d_zero));
-                          } else {
-                            exp.push_back(NodeManager::currentNM()->mkNode(
-                                LT, pivot_factor, d_zero));
-                          }
-                        }
-                        Kind jk = transKinds(types[0], types[1]);
-                        Trace("nl-ext-rbound-debug")
-                            << "trans kind : " << types[0] << " + "
-                            << types[1] << " = " << jk << std::endl;
-                        if (jk != UNDEFINED_KIND)
-                        {
-                          Node conc = NodeManager::currentNM()->mkNode(
-                              jk, rhs_a_res, rhs_b_res);
-                          Node conc_mv =
-                              d_model.computeAbstractModelValue(conc);
-                          if (conc_mv == d_false) {
-                            Node rblem = NodeManager::currentNM()->mkNode(
-                                IMPLIES,
-                                NodeManager::currentNM()->mkNode(AND, exp),
-                                conc);
-                            Trace("nl-ext-rbound-lemma-debug")
-                                << "Resolution bound lemma "
-                                   "(pre-rewrite) "
-                                   ": "
-                                << rblem << std::endl;
-                            rblem = Rewriter::rewrite(rblem);
-                            Trace("nl-ext-rbound-lemma")
-                                << "Resolution bound lemma : " << rblem
-                                << std::endl;
-                            lemmas.push_back(rblem);
-                          }
-                        }
-                        exp.pop_back();
-                        exp.pop_back();
-                        exp.pop_back();
-                      }
-                    }
-                  }
-                  exp.pop_back();
-                }
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-  return lemmas;
-}
 
 }  // namespace arith
 }  // namespace theory
diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h
index bfc713b12..5aa2070a6 100644
--- a/src/theory/arith/nonlinear_extension.h
+++ b/src/theory/arith/nonlinear_extension.h
@@ -19,23 +19,15 @@
 #define CVC4__THEORY__ARITH__NONLINEAR_EXTENSION_H
 
 #include <stdint.h>
-
 #include <map>
-#include <queue>
-#include <set>
-#include <unordered_map>
-#include <utility>
 #include <vector>
 
-#include "context/cdhashset.h"
-#include "context/cdinsert_hashmap.h"
 #include "context/cdlist.h"
-#include "context/cdqueue.h"
-#include "context/context.h"
 #include "expr/kind.h"
 #include "expr/node.h"
 #include "theory/arith/nl_lemma_utils.h"
 #include "theory/arith/nl_model.h"
+#include "theory/arith/nl_solver.h"
 #include "theory/arith/theory_arith.h"
 #include "theory/arith/transcendental_solver.h"
 #include "theory/uf/equality_engine.h"
@@ -44,9 +36,6 @@ namespace CVC4 {
 namespace theory {
 namespace arith {
 
-typedef std::map<Node, unsigned> NodeMultiset;
-
-// TODO : refactor/document this class (#1287)
 /** Non-linear extension class
  *
  * This class implements model-based refinement schemes
@@ -72,6 +61,8 @@ typedef std::map<Node, unsigned> NodeMultiset;
  * where d_out is the output channel of TheoryArith.
  */
 class NonlinearExtension {
+  typedef context::CDHashSet<Node, NodeHashFunction> NodeSet;
+
  public:
   NonlinearExtension(TheoryArith& containing, eq::EqualityEngine* ee);
   ~NonlinearExtension();
@@ -187,34 +178,7 @@ class NonlinearExtension {
   bool modelBasedRefinement(std::vector<Node>& mlems,
                             std::vector<Node>& mlemsPp,
                             std::map<Node, NlLemmaSideEffect>& lemSE);
-  /** returns true if the multiset containing the
-   * factors of monomial a is a subset of the multiset
-   * containing the factors of monomial b.
-   */
-  bool isMonomialSubset(Node a, Node b) const;
-  void registerMonomialSubset(Node a, Node b);
-
-  typedef context::CDHashSet<Node, NodeHashFunction> NodeSet;
-
-  // monomial information (context-independent)
-  class MonomialIndex {
-   public:
-    // status 0 : n equal, -1 : n superset, 1 : n subset
-    void addTerm(Node n, const std::vector<Node>& reps, NonlinearExtension* nla,
-                 int status = 0, unsigned argIndex = 0);
-
-   private:
-    std::map<Node, MonomialIndex> d_data;
-    std::vector<Node> d_monos;
-  }; /* class MonomialIndex */
 
-  // constraint information (context-independent)
-  struct ConstraintInfo {
-   public:
-    Node d_rhs;
-    Node d_coeff;
-    Kind d_type;
-  }; /* struct ConstraintInfo */
 
   /** check last call
    *
@@ -249,24 +213,6 @@ class NonlinearExtension {
                     std::vector<Node>& lemsPp,
                     std::vector<Node>& wlems,
                     std::map<Node, NlLemmaSideEffect>& lemSE);
-  //---------------------------------------term utilities
-  static bool isArithKind(Kind k);
-  static Node mkLit(Node a, Node b, int status, bool isAbsolute = false);
-  static Node mkAbs(Node a);
-  static Node mkValidPhase(Node a, Node pi);
-  Node mkMonomialRemFactor(Node n, const NodeMultiset& n_exp_rem) const;
-  //---------------------------------------end term utilities
-
-  /** register monomial */
-  void registerMonomial(Node n);
-  void setMonomialFactor(Node a, Node b, const NodeMultiset& common);
-
-  void registerConstraint(Node atom);
-  /** assign order ids */
-  void assignOrderIds(std::vector<Node>& vars,
-                      NodeMultiset& d_order,
-                      bool isConcrete,
-                      bool isAbsolute);
 
   /** get assertions
    *
@@ -312,85 +258,6 @@ class NonlinearExtension {
                   std::vector<Node>& gs);
   //---------------------------end check model
 
-  /** In the following functions, status states a relationship
-  * between two arithmetic terms, where:
-  * 0 : equal
-  * 1 : greater than or equal
-  * 2 : greater than
-  * -X : (greater -> less)
-  * TODO (#1287) make this an enum?
-  */
-  /** compute the sign of a.
-  *
-  * Calls to this function are such that :
-  *    exp => ( oa = a ^ a <status> 0 )
-  *
-  * This function iterates over the factors of a,
-  * where a_index is the index of the factor in a
-  * we are currently looking at.
-  *
-  * This function returns a status, which indicates
-  * a's relationship to 0.
-  * We add lemmas to lem of the form given by the
-  * lemma schema checkSign(...).
-  */
-  int compareSign(Node oa, Node a, unsigned a_index, int status,
-                  std::vector<Node>& exp, std::vector<Node>& lem);
-  /** compare monomials a and b
-  *
-  * Initially, a call to this function is such that :
-  *    exp => ( oa = a ^ ob = b )
-  *
-  * This function returns true if we can infer a valid
-  * arithmetic lemma of the form :
-  *    P => abs( a ) >= abs( b )
-  * where P is true and abs( a ) >= abs( b ) is false in the
-  * current model.
-  *
-  * This function is implemented by "processing" factors
-  * of monomials a and b until an inference of the above
-  * form can be made. For example, if :
-  *   a = x*x*y and b = z*w
-  * Assuming we are trying to show abs( a ) >= abs( c ),
-  * then if abs( M( x ) ) >= abs( M( z ) ) where M is the current model,
-  * then we can add abs( x ) >= abs( z ) to our explanation, and
-  * mark one factor of x as processed in a, and
-  * one factor of z as processed in b. The number of processed factors of a
-  * and b are stored in a_exp_proc and b_exp_proc respectively.
-  *
-  * cmp_infers stores information that is helpful
-  * in discarding redundant inferences.  For example,
-  * we do not want to infer abs( x ) >= abs( z ) if
-  * we have already inferred abs( x ) >= abs( y ) and
-  * abs( y ) >= abs( z ).
-  * It stores entries of the form (status,t1,t2)->F,
-  * which indicates that we constructed a lemma F that
-  * showed t1 <status> t2.
-  *
-  * We add lemmas to lem of the form given by the
-  * lemma schema checkMagnitude(...).
-  */
-  bool compareMonomial(
-      Node oa, Node a, NodeMultiset& a_exp_proc, Node ob, Node b,
-      NodeMultiset& b_exp_proc, std::vector<Node>& exp, std::vector<Node>& lem,
-      std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers);
-  /** helper function for above
-   *
-   * The difference is the inputs a_index and b_index, which are the indices of
-   * children (factors) in monomials a and b which we are currently looking at.
-   */
-  bool compareMonomial(
-      Node oa, Node a, unsigned a_index, NodeMultiset& a_exp_proc, Node ob,
-      Node b, unsigned b_index, NodeMultiset& b_exp_proc, int status,
-      std::vector<Node>& exp, std::vector<Node>& lem,
-      std::map<int, std::map<Node, std::map<Node, Node> > >& cmp_infers);
-  /** Check whether we have already inferred a relationship between monomials
-   * x and y based on the information in cmp_infers. This computes the
-   * transitive closure of the relation stored in cmp_infers.
-   */
-  bool cmp_holds(Node x, Node y,
-                 std::map<Node, std::map<Node, Node> >& cmp_infers,
-                 std::vector<Node>& exp, std::map<Node, bool>& visited);
   /** Is n entailed with polarity pol in the current context? */
   bool isEntailed(Node n, bool pol);
 
@@ -412,61 +279,19 @@ class NonlinearExtension {
   /** Process side effect se */
   void processSideEffect(const NlLemmaSideEffect& se);
 
-  // Returns the NodeMultiset for an existing monomial.
-  const NodeMultiset& getMonomialExponentMap(Node monomial) const;
-
-  // Map from monomials to var^index.
-  typedef std::map<Node, NodeMultiset> MonomialExponentMap;
-  MonomialExponentMap d_m_exp;
-
-  /**
-   * Mapping from monomials to the list of variables that occur in it. For
-   * example, x*x*y*z -> { x, y, z }.
-   */
-  std::map<Node, std::vector<Node> > d_m_vlist;
-  NodeMultiset d_m_degree;
-  // monomial index, by sorted variables
-  MonomialIndex d_m_index;
-  // list of all monomials
-  std::vector<Node> d_monomials;
-  // containment ordering
-  std::map<Node, std::vector<Node> > d_m_contain_children;
-  std::map<Node, std::vector<Node> > d_m_contain_parent;
-  std::map<Node, std::map<Node, Node> > d_m_contain_mult;
-  std::map<Node, std::map<Node, Node> > d_m_contain_umult;
-  // ( x*y, x*z, y ) for each pair of monomials ( x*y, x*z ) with common factors
-  std::map<Node, std::map<Node, Node> > d_mono_diff;
-
   /** cache of all lemmas sent on the output channel (user-context-dependent) */
   NodeSet d_lemmas;
-  /** cache of terms t for which we have added the lemma ( t = 0 V t != 0 ). */
-  NodeSet d_zero_split;
-
   /** commonly used terms */
   Node d_zero;
   Node d_one;
   Node d_neg_one;
-  Node d_two;
   Node d_true;
-  Node d_false;
-
   // The theory of arithmetic containing this extension.
   TheoryArith& d_containing;
-
   // pointer to used equality engine
   eq::EqualityEngine* d_ee;
   // needs last call effort
   bool d_needsLastCall;
-
-  // if d_c_info[lit][x] = ( r, coeff, k ), then ( lit <=>  (coeff * x) <k> r )
-  std::map<Node, std::map<Node, ConstraintInfo> > d_c_info;
-  std::map<Node, std::map<Node, bool> > d_c_info_maxm;
-  std::vector<Node> d_constraints;
-
-  // per last-call effort
-
-  // model values/orderings
-
   /** The non-linear model object
    *
    * This class is responsible for computing model values for arithmetic terms
@@ -479,6 +304,12 @@ class NonlinearExtension {
    * transcendental functions.
    */
   TranscendentalSolver d_trSlv;
+  /** The nonlinear extension object
+   *
+   * This is the subsolver responsible for running the procedure for
+   * constraints involving nonlinear mulitplication.
+   */
+  NlSolver d_nlSlv;
   /**
    * The lemmas we computed during collectModelInfo. We store two vectors of
    * lemmas to be sent out on the output channel of TheoryArith. The first
@@ -495,182 +326,6 @@ class NonlinearExtension {
   std::map<Node, std::pair<Node, Node>> d_approximations;
   /** have we successfully built the model in this SAT context? */
   context::CDO<bool> d_builtModel;
-
-  // ordering, stores variables and 0,1,-1
-  std::map<Node, unsigned> d_order_vars;
-  std::vector<Node> d_order_points;
-  
-
- private:
-  //per last-call effort check
-  
-  //information about monomials
-  std::vector< Node > d_ms;
-  std::vector< Node > d_ms_vars;
-  std::map<Node, bool> d_ms_proc;
-  std::vector<Node> d_mterms;
-
-  //list of monomials with factors whose model value is non-constant in model 
-  //  e.g. y*cos( x )
-  std::map<Node, bool> d_m_nconst_factor;
-  /** the set of monomials we should apply tangent planes to */
-  std::unordered_set<Node, NodeHashFunction> d_tplane_refine;
-  // term -> coeff -> rhs -> ( status, exp, b ),
-  //   where we have that : exp =>  ( coeff * term <status> rhs )
-  //   b is true if degree( term ) >= degree( rhs )
-  std::map<Node, std::map<Node, std::map<Node, Kind> > > d_ci;
-  std::map<Node, std::map<Node, std::map<Node, Node> > > d_ci_exp;
-  std::map<Node, std::map<Node, std::map<Node, bool> > > d_ci_max;
-
-  // factor skolems
-  std::map< Node, Node > d_factor_skolem;
-  Node getFactorSkolem( Node n, std::vector< Node >& lemmas );
-
-  // tangent plane bounds
-  std::map< Node, std::map< Node, Node > > d_tangent_val_bound[4];
-
-  /** get approximate sqrt
-   *
-   * This approximates the square root of positive constant c. If this method
-   * returns true, then l and u are updated to constants such that
-   *   l <= sqrt( c ) <= u
-   * The argument iter is the number of iterations in the binary search to
-   * perform. By default, this is set to 15, which is usually enough to be
-   * precise in the majority of simple cases, whereas not prohibitively
-   * expensive to compute.
-   */
-  bool getApproximateSqrt(Node c, Node& l, Node& u, unsigned iter = 15) const;
-
- private:
-  //-------------------------------------------- lemma schemas
-  /** check split zero
-  *
-  * Returns a set of theory lemmas of the form
-  *   t = 0 V t != 0
-  * where t is a term that exists in the current context.
-  */
-  std::vector<Node> checkSplitZero();
-
-  /** check monomial sign
-  *
-  * Returns a set of valid theory lemmas, based on a
-  * lemma schema which ensures that non-linear monomials
-  * respect sign information based on their facts.
-  * For more details, see Section 5 of "Design Theory
-  * Solvers with Extensions" by Reynolds et al., FroCoS 2017,
-  * Figure 5, this is the schema "Sign".
-  *
-  * Examples:
-  *
-  * x > 0 ^ y > 0 => x*y > 0
-  * x < 0 => x*y*y < 0
-  * x = 0 => x*y*z = 0
-  */
-  std::vector<Node> checkMonomialSign();
-
-  /** check monomial magnitude
-  *
-  * Returns a set of valid theory lemmas, based on a
-  * lemma schema which ensures that comparisons between
-  * non-linear monomials respect the magnitude of their
-  * factors.
-  * For more details, see Section 5 of "Design Theory
-  * Solvers with Extensions" by Reynolds et al., FroCoS 2017,
-  * Figure 5, this is the schema "Magnitude".
-  *
-  * Examples:
-  *
-  * |x|>|y| => |x*z|>|y*z|
-  * |x|>|y| ^ |z|>|w| ^ |x|>=1 => |x*x*z*u|>|y*w|
-  *
-  * Argument c indicates the class of inferences to perform for the (non-linear)
-  * monomials in the vector d_ms.
-  *   0 : compare non-linear monomials against 1,
-  *   1 : compare non-linear monomials against variables,
-  *   2 : compare non-linear monomials against other non-linear monomials.
-  */
-  std::vector<Node> checkMonomialMagnitude( unsigned c );
-
-  /** check monomial inferred bounds
-  *
-  * Returns a set of valid theory lemmas, based on a
-  * lemma schema that infers new constraints about existing
-  * terms based on mulitplying both sides of an existing
-  * constraint by a term.
-  * For more details, see Section 5 of "Design Theory
-  * Solvers with Extensions" by Reynolds et al., FroCoS 2017,
-  * Figure 5, this is the schema "Multiply".
-  *
-  * Examples:
-  *
-  * x > 0 ^ (y > z + w) => x*y > x*(z+w)
-  * x < 0 ^ (y > z + w) => x*y < x*(z+w)
-  *   ...where (y > z + w) and x*y are a constraint and term
-  *      that occur in the current context.
-  */
-  std::vector<Node> checkMonomialInferBounds(
-      std::vector<Node>& nt_lemmas,
-      const std::vector<Node>& asserts,
-      const std::vector<Node>& false_asserts);
-
-  /** check factoring
-  *
-  * Returns a set of valid theory lemmas, based on a
-  * lemma schema that states a relationship betwen monomials
-  * with common factors that occur in the same constraint.
-  *
-  * Examples:
-  *
-  * x*z+y*z > t => ( k = x + y ^ k*z > t )
-  *   ...where k is fresh and x*z + y*z > t is a
-  *      constraint that occurs in the current context.
-  */
-  std::vector<Node> checkFactoring(const std::vector<Node>& asserts,
-                                   const std::vector<Node>& false_asserts);
-
-  /** check monomial infer resolution bounds
-  *
-  * Returns a set of valid theory lemmas, based on a
-  * lemma schema which "resolves" upper bounds
-  * of one inequality with lower bounds for another.
-  * This schema is not enabled by default, and can
-  * be enabled by --nl-ext-rbound.
-  *
-  * Examples:
-  *
-  *  ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t
-  *  ...where s <= x*z and x*y <= t are constraints
-  *     that occur in the current context.
-  */
-  std::vector<Node> checkMonomialInferResBounds();
-
-  /** check tangent planes
-  *
-  * Returns a set of valid theory lemmas, based on an
-  * "incremental linearization" of non-linear monomials.
-  * This linearization is accomplished by adding constraints
-  * corresponding to "tangent planes" at the current
-  * model value of each non-linear monomial. In particular
-  * consider the definition for constants a,b :
-  *   T_{a,b}( x*y ) = b*x + a*y - a*b.
-  * The lemmas added by this function are of the form :
-  *  ( ( x>a ^ y<b) ^ (x<a ^ y>b) ) => x*y < T_{a,b}( x*y )
-  *  ( ( x>a ^ y>b) ^ (x<a ^ y<b) ) => x*y > T_{a,b}( x*y )
-  * It is inspired by "Invariant Checking of NRA Transition
-  * Systems via Incremental Reduction to LRA with EUF" by
-  * Cimatti et al., TACAS 2017.
-  * This schema is not terminating in general.
-  * It is not enabled by default, and can
-  * be enabled by --nl-ext-tplanes.
-  *
-  * Examples:
-  *
-  * ( ( x>2 ^ y>5) ^ (x<2 ^ y<5) ) => x*y > 5*x + 2*y - 10
-  * ( ( x>2 ^ y<5) ^ (x<2 ^ y>5) ) => x*y < 5*x + 2*y - 10
-  */
-  std::vector<Node> checkTangentPlanes();
-
-  //-------------------------------------------- end lemma schemas
 }; /* class NonlinearExtension */
 
 }  // namespace arith
diff --git a/src/theory/arith/transcendental_solver.h b/src/theory/arith/transcendental_solver.h
index 8e539bbf0..88de49af3 100644
--- a/src/theory/arith/transcendental_solver.h
+++ b/src/theory/arith/transcendental_solver.h
@@ -47,6 +47,15 @@ class TranscendentalSolver
   ~TranscendentalSolver();
 
   /** init last call
+   *
+   * This is called at the beginning of last call effort check, where
+   * assertions are the set of assertions belonging to arithmetic,
+   * false_asserts is the subset of assertions that are false in the current
+   * model, and xts is the set of extended function terms that are active in
+   * the current context.
+   *
+   * This call may add lemmas to lems/lemsPp based on registering term
+   * information (for example, purification of sine terms).
    */
   void initLastCall(const std::vector<Node>& assertions,
                     const std::vector<Node>& false_asserts,