Split NlSolver in multiple subsolvers (#5315)

The NlSolver started as one place for nonlinear reasoning, but has grown significantly since. This PR splits the NlSolver class into multiple smaller classes.

20 files changed, 2360 insertions, 1965 deletions

diff --git a/src/theory/arith/nl/nl_constraint.cpp b/src/theory/arith/nl/ext/constraint.cpp
index d4aed15f1..f395ee2d4 100644
--- a/src/theory/arith/nl/nl_constraint.cpp
+++ b/src/theory/arith/nl/ext/constraint.cpp
@@ -1,5 +1,5 @@
 /*********************                                                        */
-/*! \file nl_constraint.cpp
+/*! \file constraint.cpp
  ** \verbatim
  ** Top contributors (to current version):
  **   Andrew Reynolds, Tim King
@@ -12,7 +12,7 @@
  ** \brief Implementation of utilities for non-linear constraints
  **/
 
-#include "theory/arith/nl/nl_constraint.h"
+#include "theory/arith/nl/ext/constraint.h"
 
 #include "theory/arith/arith_msum.h"
 #include "theory/arith/arith_utilities.h"
diff --git a/src/theory/arith/nl/nl_constraint.h b/src/theory/arith/nl/ext/constraint.h
index b126eeb64..699d5b350 100644
--- a/src/theory/arith/nl/nl_constraint.h
+++ b/src/theory/arith/nl/ext/constraint.h
@@ -1,5 +1,5 @@
 /*********************                                                        */
-/*! \file nl_constraint.h
+/*! \file constraint.h
  ** \verbatim
  ** Top contributors (to current version):
  **   Andrew Reynolds
@@ -12,15 +12,15 @@
  ** \brief Utilities for non-linear constraints
  **/
 
-#ifndef CVC4__THEORY__ARITH__NL__NL_CONSTRAINT_H
-#define CVC4__THEORY__ARITH__NL__NL_CONSTRAINT_H
+#ifndef CVC4__THEORY__ARITH__NL__EXT__CONSTRAINT_H
+#define CVC4__THEORY__ARITH__NL__EXT__CONSTRAINT_H
 
 #include <map>
 #include <vector>
 
 #include "expr/kind.h"
 #include "expr/node.h"
-#include "theory/arith/nl/nl_monomial.h"
+#include "theory/arith/nl/ext/monomial.h"
 
 namespace CVC4 {
 namespace theory {
diff --git a/src/theory/arith/nl/ext/ext_state.cpp b/src/theory/arith/nl/ext/ext_state.cpp
new file mode 100644
index 000000000..3ac4699a7
--- /dev/null
+++ b/src/theory/arith/nl/ext/ext_state.cpp
@@ -0,0 +1,95 @@
+/*********************                                                        */
+/*! \file shared_check_data.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Common data shared by multiple checks
+ **/
+
+#include "theory/arith/nl/ext/ext_state.h"
+
+#include <vector>
+
+#include "expr/node.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/ext/monomial.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+ExtState::ExtState(InferenceManager& im, NlModel& model, context::Context* c)
+    : d_im(im), d_model(model)
+{
+  d_false = NodeManager::currentNM()->mkConst(false);
+  d_true = NodeManager::currentNM()->mkConst(true);
+  d_zero = NodeManager::currentNM()->mkConst(Rational(0));
+  d_one = NodeManager::currentNM()->mkConst(Rational(1));
+  d_neg_one = NodeManager::currentNM()->mkConst(Rational(-1));
+}
+
+void ExtState::init(const std::vector<Node>& xts)
+{
+  d_ms_vars.clear();
+  d_ms.clear();
+  d_mterms.clear();
+  d_tplane_refine.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 == Kind::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);
+        }
+      }
+      // mark processed if has a "one" factor (will look at reduced monomial)?
+    }
+  }
+
+  // register constants
+  d_mdb.registerMonomial(d_one);
+
+  // 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;
+}
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl/ext/ext_state.h b/src/theory/arith/nl/ext/ext_state.h
new file mode 100644
index 000000000..64c9c6b94
--- /dev/null
+++ b/src/theory/arith/nl/ext/ext_state.h
@@ -0,0 +1,67 @@
+/*********************                                                        */
+/*! \file shared_check_data.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Common data shared by multiple checks
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__EXT__SHARED_CHECK_DATA_H
+#define CVC4__THEORY__ARITH__NL__EXT__SHARED_CHECK_DATA_H
+
+#include <vector>
+
+#include "expr/node.h"
+#include "expr/proof.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/ext/monomial.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+struct ExtState
+{
+  ExtState(InferenceManager& im, NlModel& model, context::Context* c);
+
+  void init(const std::vector<Node>& xts);
+
+  Node d_false;
+  Node d_true;
+  Node d_zero;
+  Node d_one;
+  Node d_neg_one;
+
+  /** The inference manager that we push conflicts and lemmas to. */
+  InferenceManager& d_im;
+  /** Reference to the non-linear model object */
+  NlModel& d_model;
+
+  // information about monomials
+  std::vector<Node> d_ms;
+  std::vector<Node> d_ms_vars;
+  std::vector<Node> d_mterms;
+
+  /** Context-independent database of monomial information */
+  MonomialDb d_mdb;
+
+  // ( 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;
+  /** the set of monomials we should apply tangent planes to */
+  std::unordered_set<Node, NodeHashFunction> d_tplane_refine;
+};
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif
diff --git a/src/theory/arith/nl/ext/factoring_check.cpp b/src/theory/arith/nl/ext/factoring_check.cpp
new file mode 100644
index 000000000..e329db121
--- /dev/null
+++ b/src/theory/arith/nl/ext/factoring_check.cpp
@@ -0,0 +1,181 @@
+/*********************                                                        */
+/*! \file factoring_check.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 factoring check
+ **/
+
+#include "theory/arith/nl/ext/factoring_check.h"
+
+#include "expr/node.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+FactoringCheck::FactoringCheck(InferenceManager& im, NlModel& model)
+    : d_im(im), d_model(model)
+{
+  d_zero = NodeManager::currentNM()->mkConst(Rational(0));
+  d_one = NodeManager::currentNM()->mkConst(Rational(1));
+}
+
+void FactoringCheck::check(const std::vector<Node>& asserts,
+                           const std::vector<Node>& false_asserts)
+{
+  NodeManager* nm = NodeManager::currentNM();
+  Trace("nl-ext") << "Get factoring lemmas..." << std::endl;
+  for (const Node& lit : asserts)
+  {
+    bool polarity = lit.getKind() != Kind::NOT;
+    Node atom = lit.getKind() == Kind::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() == Kind::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(Kind::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(Kind::PLUS, itf->second);
+          sum = Rewriter::rewrite(sum);
+          Trace("nl-ext-factor")
+              << "* Factored sum for " << x << " : " << sum << std::endl;
+          Node kf = getFactorSkolem(sum);
+          std::vector<Node> poly;
+          poly.push_back(nm->mkNode(Kind::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(Kind::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(Kind::OR, lemma_disj);
+          Trace("nl-ext-factor") << "...lemma is " << flem << std::endl;
+          d_im.addPendingArithLemma(flem, InferenceId::NL_FACTOR);
+        }
+      }
+    }
+  }
+}
+
+Node FactoringCheck::getFactorSkolem(Node n)
+{
+  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));
+    d_im.addPendingArithLemma(k_eq, InferenceId::NL_FACTOR);
+    d_factor_skolem[n] = k;
+    return k;
+  }
+  return itf->second;
+}
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
\ No newline at end of file
diff --git a/src/theory/arith/nl/ext/factoring_check.h b/src/theory/arith/nl/ext/factoring_check.h
new file mode 100644
index 000000000..8474ac610
--- /dev/null
+++ b/src/theory/arith/nl/ext/factoring_check.h
@@ -0,0 +1,68 @@
+/*********************                                                        */
+/*! \file factoring_check.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Check for factoring lemma
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__EXT__FACTORING_CHECK_H
+#define CVC4__THEORY__ARITH__NL__EXT__FACTORING_CHECK_H
+
+#include <vector>
+
+#include "expr/node.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+class FactoringCheck
+{
+ public:
+  FactoringCheck(InferenceManager& im, NlModel& model);
+
+  /** 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.
+   */
+  void check(const std::vector<Node>& asserts,
+             const std::vector<Node>& false_asserts);
+
+ private:
+  /** The inference manager that we push conflicts and lemmas to. */
+  InferenceManager& d_im;
+  /** Reference to the non-linear model object */
+  NlModel& d_model;
+  /** maps nodes to their factor skolems */
+  std::map<Node, Node> d_factor_skolem;
+
+  Node d_zero;
+  Node d_one;
+
+  Node getFactorSkolem(Node n);
+};
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif
diff --git a/src/theory/arith/nl/nl_monomial.cpp b/src/theory/arith/nl/ext/monomial.cpp
index 613d204f6..0b46cc88e 100644
--- a/src/theory/arith/nl/nl_monomial.cpp
+++ b/src/theory/arith/nl/ext/monomial.cpp
@@ -1,5 +1,5 @@
 /*********************                                                        */
-/*! \file nl_monomial.cpp
+/*! \file monomial.cpp
  ** \verbatim
  ** Top contributors (to current version):
  **   Andrew Reynolds, Tim King, Andres Noetzli
@@ -12,7 +12,7 @@
  ** \brief Implementation of utilities for monomials
  **/
 
-#include "theory/arith/nl/nl_monomial.h"
+#include "theory/arith/nl/ext/monomial.h"
 
 #include "theory/arith/arith_utilities.h"
 #include "theory/arith/nl/nl_lemma_utils.h"
diff --git a/src/theory/arith/nl/nl_monomial.h b/src/theory/arith/nl/ext/monomial.h
index 0875919bc..93a291ca0 100644
--- a/src/theory/arith/nl/nl_monomial.h
+++ b/src/theory/arith/nl/ext/monomial.h
@@ -1,5 +1,5 @@
 /*********************                                                        */
-/*! \file nl_monomial.h
+/*! \file monomial.h
  ** \verbatim
  ** Top contributors (to current version):
  **   Andrew Reynolds, Tim King
@@ -12,8 +12,8 @@
  ** \brief Utilities for monomials
  **/
 
-#ifndef CVC4__THEORY__ARITH__NL__NL_MONOMIAL_H
-#define CVC4__THEORY__ARITH__NL__NL_MONOMIAL_H
+#ifndef CVC4__THEORY__ARITH__NL__EXT__MONOMIAL_H
+#define CVC4__THEORY__ARITH__NL__EXT__MONOMIAL_H
 
 #include <map>
 #include <vector>
diff --git a/src/theory/arith/nl/ext/monomial_bounds_check.cpp b/src/theory/arith/nl/ext/monomial_bounds_check.cpp
new file mode 100644
index 000000000..fe4b1b83a
--- /dev/null
+++ b/src/theory/arith/nl/ext/monomial_bounds_check.cpp
@@ -0,0 +1,500 @@
+/*********************                                                        */
+/*! \file monomial_bounds_check.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Check for monomial bound inference lemmas
+ **/
+
+#include "theory/arith/nl/ext/monomial_bounds_check.h"
+
+#include "expr/node.h"
+#include "options/arith_options.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/arith_utilities.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+namespace {
+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() == Kind::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;
+}
+}  // namespace
+
+MonomialBoundsCheck::MonomialBoundsCheck(ExtState* data)
+    : d_data(data), d_cdb(d_data->d_mdb)
+{
+}
+
+void MonomialBoundsCheck::init()
+{
+  d_ci.clear();
+  d_ci_exp.clear();
+  d_ci_max.clear();
+}
+
+void MonomialBoundsCheck::checkBounds(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_data->d_mdb.sortByDegree(d_data->d_ms);
+  // all monomials
+  d_data->d_mterms.insert(d_data->d_mterms.end(),
+                          d_data->d_ms_vars.begin(),
+                          d_data->d_ms_vars.end());
+  d_data->d_mterms.insert(
+      d_data->d_mterms.end(), d_data->d_ms.begin(), d_data->d_ms.end());
+
+  const std::map<Node, std::map<Node, ConstraintInfo> >& cim =
+      d_cdb.getConstraints();
+
+  NodeManager* nm = NodeManager::currentNM();
+  // register constraints
+  Trace("nl-ext-debug") << "Register bound constraints..." << std::endl;
+  for (const Node& lit : asserts)
+  {
+    bool polarity = lit.getKind() != Kind::NOT;
+    Node atom = lit.getKind() == Kind::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 == Kind::EQUAL)
+        {
+          // we will take the strict inequality in the direction of the
+          // model
+          Node lhs = ArithMSum::mkCoeffTerm(coeff, x);
+          Node query = nm->mkNode(Kind::GT, lhs, rhs);
+          Node query_mv = d_data->d_model.computeAbstractModelValue(query);
+          if (query_mv == d_data->d_true)
+          {
+            exp = query;
+            type = Kind::GT;
+          }
+          else
+          {
+            Assert(query_mv == d_data->d_false);
+            exp = nm->mkNode(Kind::LT, lhs, rhs);
+            type = Kind::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 == Kind::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(Kind::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_data->d_tplane_refine.insert(x);
+        }
+      }
+    }
+  }
+  // reflexive constraints
+  Node null_coeff;
+  for (unsigned j = 0; j < d_data->d_mterms.size(); j++)
+  {
+    Node n = d_data->d_mterms[j];
+    d_ci[n][null_coeff][n] = Kind::EQUAL;
+    d_ci_exp[n][null_coeff][n] = d_data->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_data->d_mdb.getContainsParentMap();
+  for (unsigned k = 0; k < d_data->d_mterms.size(); k++)
+  {
+    Node x = d_data->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_data->d_mdb.getContainsDiff(x, y);
+          // x <k> t => m*x <k'> t  where y = m*x
+          // get the sign of mult
+          Node mmv = d_data->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(Kind::MULT, mult, t);
+          Node infer_rhs = nm->mkNode(Kind::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_data->d_model.computeAbstractModelValue(infer);
+          Trace("nl-ext-bound-debug")
+              << "       ...infer model value is " << infer_mv << std::endl;
+          if (infer_mv == d_data->d_false)
+          {
+            Node exp = nm->mkNode(
+                Kind::AND,
+                nm->mkNode(
+                    mmv_sign == 1 ? Kind::GT : Kind::LT, mult, d_data->d_zero),
+                d_ci_exp[x][coeff][rhs]);
+            Node iblem = nm->mkNode(Kind::IMPLIES, exp, infer);
+            Node pr_iblem = iblem;
+            iblem = Rewriter::rewrite(iblem);
+            bool introNewTerms = hasNewMonomials(iblem, d_data->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;
+            d_data->d_im.addPendingArithLemma(
+                iblem, InferenceId::NL_INFER_BOUNDS_NT, introNewTerms);
+          }
+        }
+      }
+    }
+  }
+}
+
+void MonomialBoundsCheck::checkResBounds()
+{
+  NodeManager* nm = NodeManager::currentNM();
+  Trace("nl-ext") << "Get monomial resolution inferred bound lemmas..."
+                  << std::endl;
+  size_t nmterms = d_data->d_mterms.size();
+  for (unsigned j = 0; j < nmterms; j++)
+  {
+    Node a = d_data->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_data->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_data->d_mono_diff[a].find(b);
+      if (ita == d_data->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_data->d_mono_diff[b].find(a);
+      Assert(itb != d_data->d_mono_diff[b].end());
+      Node mv_a = d_data->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_data->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(Kind::MULT, itb->second, rhs_a);
+          rhs_a_res_base = Rewriter::rewrite(rhs_a_res_base);
+          if (hasNewMonomials(rhs_a_res_base, d_data->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(Kind::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_data->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(Kind::GT, pivot_factor, d_data->d_zero));
+                }
+                else
+                {
+                  exp.push_back(
+                      nm->mkNode(Kind::LT, pivot_factor, d_data->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 != Kind::UNDEFINED_KIND)
+              {
+                Node conc = nm->mkNode(jk, rhs_a_res, rhs_b_res);
+                Node conc_mv = d_data->d_model.computeAbstractModelValue(conc);
+                if (conc_mv == d_data->d_false)
+                {
+                  Node rblem = nm->mkNode(Kind::IMPLIES, nm->mkAnd(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;
+                  d_data->d_im.addPendingArithLemma(
+                      rblem, InferenceId::NL_RES_INFER_BOUNDS);
+                }
+              }
+              exp.pop_back();
+              exp.pop_back();
+              exp.pop_back();
+            }
+          }
+          exp.pop_back();
+        }
+      }
+    }
+  }
+}
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
\ No newline at end of file
diff --git a/src/theory/arith/nl/ext/monomial_bounds_check.h b/src/theory/arith/nl/ext/monomial_bounds_check.h
new file mode 100644
index 000000000..d919b1272
--- /dev/null
+++ b/src/theory/arith/nl/ext/monomial_bounds_check.h
@@ -0,0 +1,90 @@
+/*********************                                                        */
+/*! \file monomial_bounds_check.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Check for monomial bound inference lemmas
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__EXT__MONOMIAL_BOUNDS_CHECK_H
+#define CVC4__THEORY__ARITH__NL__EXT__MONOMIAL_BOUNDS_CHECK_H
+
+#include "expr/node.h"
+#include "theory/arith/nl/ext/constraint.h"
+#include "theory/arith/nl/ext/ext_state.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+class MonomialBoundsCheck
+{
+ public:
+  MonomialBoundsCheck(ExtState* data);
+
+  void init();
+
+  /** 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.
+   */
+  void checkBounds(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.
+   */
+  void checkResBounds();
+
+ private:
+  /** Basic data that is shared with other checks */
+  ExtState* d_data;
+
+  /** Context-independent database of constraint information */
+  ConstraintDb d_cdb;
+
+  // 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;
+};
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif
diff --git a/src/theory/arith/nl/ext/monomial_check.cpp b/src/theory/arith/nl/ext/monomial_check.cpp
new file mode 100644
index 000000000..3b8d52c13
--- /dev/null
+++ b/src/theory/arith/nl/ext/monomial_check.cpp
@@ -0,0 +1,740 @@
+/*********************                                                        */
+/*! \file monomial_check.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Check for some monomial lemmas
+ **/
+
+#include "theory/arith/nl/ext/monomial_check.h"
+
+#include "expr/node.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+MonomialCheck::MonomialCheck(ExtState* data) : d_data(data)
+{
+  d_order_points.push_back(d_data->d_neg_one);
+  d_order_points.push_back(d_data->d_zero);
+  d_order_points.push_back(d_data->d_one);
+}
+
+void MonomialCheck::init(const std::vector<Node>& xts)
+{
+  d_ms_proc.clear();
+  d_m_nconst_factor.clear();
+
+  for (unsigned i = 0, xsize = xts.size(); i < xsize; i++)
+  {
+    Node a = xts[i];
+    if (a.getKind() == Kind::NONLINEAR_MULT)
+    {
+      const std::vector<Node>& varList = d_data->d_mdb.getVariableList(a);
+      for (const Node& v : varList)
+      {
+        Node mvk = d_data->d_model.computeAbstractModelValue(v);
+        if (!mvk.isConst())
+        {
+          d_m_nconst_factor[a] = true;
+        }
+      }
+    }
+  }
+
+  for (unsigned j = 0; j < d_order_points.size(); j++)
+  {
+    Node c = d_order_points[j];
+    d_data->d_model.computeConcreteModelValue(c);
+    d_data->d_model.computeAbstractModelValue(c);
+  }
+}
+
+void MonomialCheck::checkSign()
+{
+  std::map<Node, int> signs;
+  Trace("nl-ext") << "Get monomial sign lemmas..." << std::endl;
+  for (unsigned j = 0; j < d_data->d_ms.size(); j++)
+  {
+    Node a = d_data->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_data->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);
+        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;
+      }
+    }
+  }
+}
+
+void MonomialCheck::checkMagnitude(unsigned c)
+{
+  // ensure information is setup
+  if (c == 0)
+  {
+    // sort by absolute values of abstract model values
+    assignOrderIds(d_data->d_ms_vars, d_order_vars, false, true);
+
+    // sort individual variable lists
+    Trace("nl-ext-proc") << "Assign order var lists..." << std::endl;
+    d_data->d_mdb.sortVariablesByModel(d_data->d_ms, d_data->d_model);
+  }
+
+  unsigned r = 1;
+  std::vector<ArithLemma> 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_data->d_ms.size(); j++)
+  {
+    Node a = d_data->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_data->d_one,
+                        d_data->d_one,
+                        b_exp_proc,
+                        exp,
+                        lemmas,
+                        cmp_infers);
+      }
+      else
+      {
+        const NodeMultiset& mea = d_data->d_mdb.getMonomialExponentMap(a);
+        if (c == 1)
+        {
+          // could compare not just against containing variables?
+          // compare magnitude against variables
+          for (unsigned k = 0; k < d_data->d_ms_vars.size(); k++)
+          {
+            Node v = d_data->d_ms_vars[k];
+            Node mvcv = d_data->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_data->d_ms.size(); k++)
+          {
+            Node b = d_data->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_data->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::unordered_set<Node, NodeHashFunction> 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.insert(itc2->second);
+              Trace("nl-ext-comp")
+                  << "...inference of " << itc->first << " > " << itc2->first
+                  << " was redundant." << std::endl;
+              break;
+            }
+          }
+        }
+      }
+    }
+  }
+  for (unsigned i = 0; i < lemmas.size(); i++)
+  {
+    if (r_lemmas.find(lemmas[i].d_node) == r_lemmas.end())
+    {
+      d_data->d_im.addPendingArithLemma(lemmas[i]);
+    }
+  }
+  // could only take maximal lower/minimial lower bounds?
+}
+
+// show a <> 0 by inequalities between variables in monomial a w.r.t 0
+int MonomialCheck::compareSign(
+    Node oa, Node a, unsigned a_index, int status, std::vector<Node>& exp)
+{
+  Trace("nl-ext-debug") << "Process " << a << " at index " << a_index
+                        << ", status is " << status << std::endl;
+  NodeManager* nm = NodeManager::currentNM();
+  Node mvaoa = d_data->d_model.computeAbstractModelValue(oa);
+  const std::vector<Node>& vla = d_data->d_mdb.getVariableList(a);
+  if (a_index == vla.size())
+  {
+    if (mvaoa.getConst<Rational>().sgn() != status)
+    {
+      Node lemma =
+          nm->mkAnd(exp).impNode(mkLit(oa, d_data->d_zero, status * 2));
+      d_data->d_im.addPendingArithLemma(lemma, InferenceId::NL_SIGN);
+    }
+    return status;
+  }
+  Assert(a_index < vla.size());
+  Node av = vla[a_index];
+  unsigned aexp = d_data->d_mdb.getExponent(a, av);
+  // take current sign in model
+  Node mvaav = d_data->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_data->d_zero).impNode(oa.eqNode(d_data->d_zero));
+      d_data->d_im.addPendingArithLemma(lemma, InferenceId::NL_SIGN);
+    }
+    return 0;
+  }
+  if (aexp % 2 == 0)
+  {
+    exp.push_back(av.eqNode(d_data->d_zero).negate());
+    return compareSign(oa, a, a_index + 1, status, exp);
+  }
+  exp.push_back(nm->mkNode(sgn == 1 ? Kind::GT : Kind::LT, av, d_data->d_zero));
+  return compareSign(oa, a, a_index + 1, status * sgn, exp);
+}
+
+bool MonomialCheck::compareMonomial(
+    Node oa,
+    Node a,
+    NodeMultiset& a_exp_proc,
+    Node ob,
+    Node b,
+    NodeMultiset& b_exp_proc,
+    std::vector<Node>& exp,
+    std::vector<ArithLemma>& 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;
+}
+
+// trying to show a ( >, = ) b   by inequalities between variables in
+// monomials a,b
+bool MonomialCheck::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<ArithLemma>& 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_data->d_mdb.getVariableList(a);
+  const std::vector<Node>& vlb = d_data->d_mdb.getVariableList(b);
+  if (a_index == vla.size() && b_index == vlb.size())
+  {
+    // finished, compare absolute value of abstract model values
+    int modelStatus = d_data->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_data->d_zero).negate());
+        }
+      }
+      NodeManager* nm = NodeManager::currentNM();
+      Node clem = nm->mkNode(
+          Kind::IMPLIES, nm->mkAnd(exp), mkLit(oa, ob, status, true));
+      Trace("nl-ext-comp-lemma") << "comparison lemma : " << clem << std::endl;
+      lem.emplace_back(
+          clem, LemmaProperty::NONE, nullptr, InferenceId::NL_COMPARISON);
+      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_data->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_data->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_data->d_one) != d_order_vars.end());
+  unsigned ovo = d_order_vars[d_data->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_data->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_data->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_data->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_data->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;
+}
+
+bool MonomialCheck::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;
+}
+
+void MonomialCheck::assignOrderIds(std::vector<Node>& vars,
+                                   NodeMultiset& order,
+                                   bool isConcrete,
+                                   bool isAbsolute)
+{
+  SortNlModel smv;
+  smv.d_nlm = &d_data->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_data->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_data->d_model.computeModelValue(
+              d_order_points[order_index], isConcrete);
+          if (d_data->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_data->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++;
+  }
+}
+Node MonomialCheck::mkLit(Node a, Node b, int status, bool isAbsolute) const
+{
+  if (status == 0)
+  {
+    Node a_eq_b = a.eqNode(b);
+    if (!isAbsolute)
+    {
+      return a_eq_b;
+    }
+    Node negate_b = NodeManager::currentNM()->mkNode(Kind::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 ? Kind::GEQ : Kind::GT;
+  if (!isAbsolute)
+  {
+    return nm->mkNode(greater_op, a, b);
+  }
+  // return nm->mkNode( greater_op, mkAbs( a ), mkAbs( b ) );
+  Node a_is_nonnegative = nm->mkNode(Kind::GEQ, a, d_data->d_zero);
+  Node b_is_nonnegative = nm->mkNode(Kind::GEQ, b, d_data->d_zero);
+  Node negate_a = nm->mkNode(Kind::UMINUS, a);
+  Node negate_b = nm->mkNode(Kind::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)));
+}
+
+void MonomialCheck::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_data->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_data->d_mdb.mkMonomialRemFactor(a, common);
+  }
+}
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
\ No newline at end of file
diff --git a/src/theory/arith/nl/ext/monomial_check.h b/src/theory/arith/nl/ext/monomial_check.h
new file mode 100644
index 000000000..89b5847a2
--- /dev/null
+++ b/src/theory/arith/nl/ext/monomial_check.h
@@ -0,0 +1,196 @@
+/*********************                                                        */
+/*! \file monomial_check.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Check for some monomial lemmas
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__EXT__MONOMIAL_CHECK_H
+#define CVC4__THEORY__ARITH__NL__EXT__MONOMIAL_CHECK_H
+
+#include "expr/node.h"
+#include "theory/arith/nl/ext/ext_state.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+class MonomialCheck
+{
+ public:
+  MonomialCheck(ExtState* data);
+
+  void init(const std::vector<Node>& xts);
+
+  /** 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
+   */
+  void checkSign();
+
+  /** 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.
+   */
+  void checkMagnitude(unsigned c);
+
+ private:
+  /** 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);
+  /** 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<ArithLemma>& 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<ArithLemma>& 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);
+  /** assign order ids */
+  void assignOrderIds(std::vector<Node>& vars,
+                      NodeMultiset& d_order,
+                      bool isConcrete,
+                      bool isAbsolute);
+  /** Make literal */
+  Node mkLit(Node a, Node b, int status, bool isAbsolute = false) const;
+  /** register monomial */
+  void setMonomialFactor(Node a, Node b, const NodeMultiset& common);
+
+  /** Basic data that is shared with other checks */
+  ExtState* d_data;
+
+  std::map<Node, bool> d_ms_proc;
+  // ordering, stores variables and 0,1,-1
+  std::map<Node, unsigned> d_order_vars;
+  std::vector<Node> d_order_points;
+
+  // 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;
+};
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif
diff --git a/src/theory/arith/nl/ext/split_zero_check.cpp b/src/theory/arith/nl/ext/split_zero_check.cpp
new file mode 100644
index 000000000..fcbf84be4
--- /dev/null
+++ b/src/theory/arith/nl/ext/split_zero_check.cpp
@@ -0,0 +1,54 @@
+/*********************                                                        */
+/*! \file split_zero_check.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 split zero check
+ **/
+
+#include "theory/arith/nl/ext/split_zero_check.h"
+
+#include "expr/node.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+SplitZeroCheck::SplitZeroCheck(ExtState* data, context::Context* ctx)
+    : d_data(data), d_zero_split(ctx)
+{
+}
+
+void SplitZeroCheck::check()
+{
+  for (unsigned i = 0; i < d_data->d_ms_vars.size(); i++)
+  {
+    Node v = d_data->d_ms_vars[i];
+    if (d_zero_split.insert(v))
+    {
+      Node eq = v.eqNode(d_data->d_zero);
+      eq = Rewriter::rewrite(eq);
+      d_data->d_im.addPendingPhaseRequirement(eq, true);
+      ArithLemma lem(eq.orNode(eq.negate()),
+                     LemmaProperty::NONE,
+                     nullptr,
+                     InferenceId::NL_SPLIT_ZERO);
+      d_data->d_im.addPendingArithLemma(lem);
+    }
+  }
+}
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
\ No newline at end of file
diff --git a/src/theory/arith/nl/ext/split_zero_check.h b/src/theory/arith/nl/ext/split_zero_check.h
new file mode 100644
index 000000000..1a9f0a44d
--- /dev/null
+++ b/src/theory/arith/nl/ext/split_zero_check.h
@@ -0,0 +1,53 @@
+/*********************                                                        */
+/*! \file split_zero_check.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Check for split zero lemma
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__EXT__SPLIT_ZERO_CHECK_H
+#define CVC4__THEORY__ARITH__NL__EXT__SPLIT_ZERO_CHECK_H
+
+#include "expr/node.h"
+#include "theory/arith/nl/ext/ext_state.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+class SplitZeroCheck
+{
+ public:
+  SplitZeroCheck(ExtState* data, context::Context* ctx);
+
+  /** 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.
+   */
+  void check();
+
+ private:
+  using NodeSet = context::CDHashSet<Node, NodeHashFunction>;
+
+  /** Basic data that is shared with other checks */
+  ExtState* d_data;
+  /** cache of terms t for which we have added the lemma ( t = 0 V t != 0 ). */
+  NodeSet d_zero_split;
+};
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif
diff --git a/src/theory/arith/nl/ext/tangent_plane_check.cpp b/src/theory/arith/nl/ext/tangent_plane_check.cpp
new file mode 100644
index 000000000..ff66be0e4
--- /dev/null
+++ b/src/theory/arith/nl/ext/tangent_plane_check.cpp
@@ -0,0 +1,198 @@
+/*********************                                                        */
+/*! \file tangent_plane_check.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 tangent_plane check
+ **/
+
+#include "theory/arith/nl/ext/tangent_plane_check.h"
+
+#include "expr/node.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+TangentPlaneCheck::TangentPlaneCheck(ExtState* data) : d_data(data) {}
+
+void TangentPlaneCheck::check(bool asWaitingLemmas)
+{
+  Trace("nl-ext") << "Get monomial tangent plane lemmas..." << std::endl;
+  NodeManager* nm = NodeManager::currentNM();
+  const std::map<Node, std::vector<Node> >& ccMap =
+      d_data->d_mdb.getContainsChildrenMap();
+  unsigned kstart = d_data->d_ms_vars.size();
+  for (unsigned k = kstart; k < d_data->d_mterms.size(); k++)
+  {
+    Node t = d_data->d_mterms[k];
+    // if this term requires a refinement
+    if (d_data->d_tplane_refine.find(t) == d_data->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_data->d_one)
+      {
+        Node tc_diff = d_data->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_data->d_model.computeAbstractModelValue(a);
+          Node b_v_c = d_data->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) ? Kind::GT : Kind::LT, curr_v, pt_v);
+                do_extend = Rewriter::rewrite(do_extend);
+                if (do_extend == d_data->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(Kind::MINUS,
+                                     nm->mkNode(Kind::PLUS,
+                                                nm->mkNode(Kind::MULT, b_v, a),
+                                                nm->mkNode(Kind::MULT, a_v, b)),
+                                     nm->mkNode(Kind::MULT, a_v, b_v));
+            // conjuncts of the tangent plane lemma
+            std::vector<Node> tplaneConj;
+            for (unsigned d = 0; d < 4; d++)
+            {
+              Node aa =
+                  nm->mkNode(d == 0 || d == 3 ? Kind::GEQ : Kind::LEQ, a, a_v);
+              Node ab =
+                  nm->mkNode(d == 1 || d == 3 ? Kind::GEQ : Kind::LEQ, b, b_v);
+              Node conc = nm->mkNode(d <= 1 ? Kind::LEQ : Kind::GEQ, t, tplane);
+              Node tlem = nm->mkNode(Kind::OR, aa.negate(), ab.negate(), conc);
+              Trace("nl-ext-tplanes")
+                  << "Tangent plane lemma : " << tlem << std::endl;
+              tplaneConj.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(Kind::LEQ, a, a_v);
+            Node b_leq_bv = nm->mkNode(Kind::LEQ, b, b_v);
+            Node a_geq_av = nm->mkNode(Kind::GEQ, a, a_v);
+            Node b_geq_bv = nm->mkNode(Kind::GEQ, b, b_v);
+
+            Node t_leq_tplane = nm->mkNode(Kind::LEQ, t, tplane);
+            Node a_leq_av_or_b_leq_bv =
+                nm->mkNode(Kind::OR, a_leq_av, b_leq_bv);
+            Node b_geq_bv_or_a_geq_av =
+                nm->mkNode(Kind::OR, b_geq_bv, a_geq_av);
+            Node ub_reverse1 = nm->mkNode(
+                Kind::OR, t_leq_tplane.negate(), a_leq_av_or_b_leq_bv);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << ub_reverse1
+                << std::endl;
+            tplaneConj.push_back(ub_reverse1);
+            Node ub_reverse2 = nm->mkNode(
+                Kind::OR, t_leq_tplane.negate(), b_geq_bv_or_a_geq_av);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << ub_reverse2
+                << std::endl;
+            tplaneConj.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(Kind::GEQ, t, tplane);
+            Node a_leq_av_or_b_geq_bv =
+                nm->mkNode(Kind::OR, a_leq_av, b_geq_bv);
+            Node a_geq_av_or_b_leq_bv =
+                nm->mkNode(Kind::OR, a_geq_av, b_leq_bv);
+            Node lb_reverse1 = nm->mkNode(
+                Kind::OR, t_geq_tplane.negate(), a_leq_av_or_b_geq_bv);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << lb_reverse1
+                << std::endl;
+            tplaneConj.push_back(lb_reverse1);
+            Node lb_reverse2 = nm->mkNode(
+                Kind::OR, t_geq_tplane.negate(), a_geq_av_or_b_leq_bv);
+            Trace("nl-ext-tplanes")
+                << "Tangent plane lemma (reverse) : " << lb_reverse2
+                << std::endl;
+            tplaneConj.push_back(lb_reverse2);
+
+            Node tlem = nm->mkAnd(tplaneConj);
+            d_data->d_im.addPendingArithLemma(
+                tlem, InferenceId::NL_TANGENT_PLANE, asWaitingLemmas);
+          }
+        }
+      }
+    }
+  }
+}
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
\ No newline at end of file
diff --git a/src/theory/arith/nl/ext/tangent_plane_check.h b/src/theory/arith/nl/ext/tangent_plane_check.h
new file mode 100644
index 000000000..748ab6372
--- /dev/null
+++ b/src/theory/arith/nl/ext/tangent_plane_check.h
@@ -0,0 +1,69 @@
+/*********************                                                        */
+/*! \file tangent_plane_check.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Gereon Kremer
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2020 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 Check for tangent_plane lemma
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__EXT__TANGENT_PLANE_CHECK_H
+#define CVC4__THEORY__ARITH__NL__EXT__TANGENT_PLANE_CHECK_H
+
+#include "expr/node.h"
+#include "theory/arith/nl/ext/ext_state.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+
+class TangentPlaneCheck
+{
+ public:
+  TangentPlaneCheck(ExtState* data);
+
+  /** 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
+   */
+  void check(bool asWaitingLemmas);
+
+ private:
+  /** Basic data that is shared with other checks */
+  ExtState* d_data;
+  /** tangent plane bounds */
+  std::map<Node, std::map<Node, Node> > d_tangent_val_bound[4];
+};
+
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif
diff --git a/src/theory/arith/nl/nl_solver.cpp b/src/theory/arith/nl/nl_solver.cpp
deleted file mode 100644
index 5ffba7229..000000000
--- a/src/theory/arith/nl/nl_solver.cpp
+++ /dev/null
@@ -1,1564 +0,0 @@
-/*********************                                                        */
-/*! \file nl_solver.cpp
- ** \verbatim
- ** Top contributors (to current version):
- **   Andrew Reynolds, Tim King, Ahmed Irfan
- ** This file is part of the CVC4 project.
- ** Copyright (c) 2009-2020 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/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 {
-namespace nl {
-
-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(InferenceManager& im, ArithState& astate, NlModel& model)
-    : d_im(im),
-      d_astate(astate),
-      d_model(model),
-      d_cdb(d_mdb),
-      d_zero_split(d_astate.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);
-  }
-}
-
-void NlSolver::checkSplitZero()
-{
-  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);
-      d_im.addPendingPhaseRequirement(eq, true);
-      Node lem = eq.orNode(eq.negate());
-      d_im.addPendingArithLemma(lem, InferenceId::NL_SPLIT_ZERO);
-    }
-  }
-}
-
-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)
-{
-  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));
-      d_im.addPendingArithLemma(lemma, InferenceId::NL_SIGN);
-    }
-    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));
-      d_im.addPendingArithLemma(lemma, InferenceId::NL_SIGN);
-    }
-    return 0;
-  }
-  if (aexp % 2 == 0)
-  {
-    exp.push_back(av.eqNode(d_zero).negate());
-    return compareSign(oa, a, a_index + 1, status, exp);
-  }
-  exp.push_back(nm->mkNode(sgn == 1 ? GT : LT, av, d_zero));
-  return compareSign(oa, a, a_index + 1, status * sgn, exp);
-}
-
-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<ArithLemma>& 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<ArithLemma>& 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.emplace_back(clem, LemmaProperty::NONE, nullptr, InferenceId::NL_COMPARISON);
-      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;
-}
-
-void NlSolver::checkMonomialSign()
-{
-  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);
-        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;
-      }
-    }
-  }
-}
-
-void 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<ArithLemma> 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::unordered_set<Node, NodeHashFunction> 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.insert(itc2->second);
-              Trace("nl-ext-comp")
-                  << "...inference of " << itc->first << " > " << itc2->first
-                  << " was redundant." << std::endl;
-              break;
-            }
-          }
-        }
-      }
-    }
-  }
-  for (unsigned i = 0; i < lemmas.size(); i++)
-  {
-    if (r_lemmas.find(lemmas[i].d_node) == r_lemmas.end())
-    {
-      d_im.addPendingArithLemma(lemmas[i]);
-    }
-  }
-  // could only take maximal lower/minimial lower bounds?
-}
-
-void NlSolver::checkTangentPlanes(bool asWaitingLemmas)
-{
-  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));
-            // conjuncts of the tangent plane lemma
-            std::vector<Node> tplaneConj;
-            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;
-              tplaneConj.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;
-            tplaneConj.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;
-            tplaneConj.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;
-            tplaneConj.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;
-            tplaneConj.push_back(lb_reverse2);
-
-            Node tlem = nm->mkNode(AND, tplaneConj);
-            d_im.addPendingArithLemma(
-                tlem, InferenceId::NL_TANGENT_PLANE, asWaitingLemmas);
-          }
-        }
-      }
-    }
-  }
-}
-
-void NlSolver::checkMonomialInferBounds(
-    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();
-
-  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;
-            d_im.addPendingArithLemma(iblem, InferenceId::NL_INFER_BOUNDS_NT, introNewTerms);
-          }
-        }
-      }
-    }
-  }
-}
-
-void NlSolver::checkFactoring(const std::vector<Node>& asserts, const std::vector<Node>& false_asserts)
-{
-  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);
-          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;
-          d_im.addPendingArithLemma(flem, InferenceId::NL_FACTOR);
-        }
-      }
-    }
-  }
-}
-
-Node NlSolver::getFactorSkolem(Node n)
-{
-  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));
-    d_im.addPendingArithLemma(k_eq, InferenceId::NL_FACTOR);
-    d_factor_skolem[n] = k;
-    return k;
-  }
-  return itf->second;
-}
-
-void NlSolver::checkMonomialInferResBounds()
-{
-  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;
-                  d_im.addPendingArithLemma(rblem, InferenceId::NL_RES_INFER_BOUNDS);
-                }
-              }
-              exp.pop_back();
-              exp.pop_back();
-              exp.pop_back();
-            }
-          }
-          exp.pop_back();
-        }
-      }
-    }
-  }
-}
-
-}  // namespace nl
-}  // namespace arith
-}  // namespace theory
-}  // namespace CVC4
diff --git a/src/theory/arith/nl/nl_solver.h b/src/theory/arith/nl/nl_solver.h
deleted file mode 100644
index 903220fb2..000000000
--- a/src/theory/arith/nl/nl_solver.h
+++ /dev/null
@@ -1,370 +0,0 @@
-/*********************                                                        */
-/*! \file nl_solver.h
- ** \verbatim
- ** Top contributors (to current version):
- **   Andrew Reynolds, Tim King, Tianyi Liang
- ** This file is part of the CVC4 project.
- ** Copyright (c) 2009-2020 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/nl_constraint.h"
-#include "theory/arith/nl/nl_lemma_utils.h"
-#include "theory/arith/nl/nl_model.h"
-#include "theory/arith/nl/nl_monomial.h"
-#include "theory/arith/inference_manager.h"
-
-namespace CVC4 {
-namespace theory {
-namespace arith {
-namespace nl {
-
-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(InferenceManager& im, ArithState& astate, 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.
-   */
-  void 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
-   */
-  void 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.
-   */
-  void 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.
-   */
-  void checkMonomialInferBounds(
-      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.
-   */
-  void 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.
-   */
-  void 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
-   */
-  void checkTangentPlanes(bool asWaitingLemmas);
-
-  //-------------------------------------------- end lemma schemas
- private:
-  /** The inference manager that we push conflicts and lemmas to. */
-  InferenceManager& d_im;
-  /** Reference to the state. */
-  ArithState& d_astate;
-  /** 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);
-  /** 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<ArithLemma>& 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<ArithLemma>& 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);
-}; /* class NlSolver */
-
-}  // namespace nl
-}  // namespace arith
-}  // namespace theory
-}  // namespace CVC4
-
-#endif /* CVC4__THEORY__ARITH__NL_SOLVER_H */
diff --git a/src/theory/arith/nl/nonlinear_extension.cpp b/src/theory/arith/nl/nonlinear_extension.cpp
index 8b7d8f4f5..a5ac8e3fe 100644
--- a/src/theory/arith/nl/nonlinear_extension.cpp
+++ b/src/theory/arith/nl/nonlinear_extension.cpp
@@ -47,7 +47,12 @@ NonlinearExtension::NonlinearExtension(TheoryArith& containing,
                   containing.getOutputChannel()),
       d_model(containing.getSatContext()),
       d_trSlv(d_im, d_model),
-      d_nlSlv(d_im, state, d_model),
+      d_extState(d_im, d_model, containing.getSatContext()),
+      d_factoringSlv(d_im, d_model),
+      d_monomialBoundsSlv(&d_extState),
+      d_monomialSlv(&d_extState),
+      d_splitZeroSlv(&d_extState, state.getUserContext()),
+      d_tangentPlaneSlv(&d_extState),
       d_cadSlv(d_im, d_model),
       d_icpSlv(d_im),
       d_iandSlv(d_im, state, d_model),
@@ -686,7 +691,7 @@ void NonlinearExtension::runStrategy(Theory::Effort effort,
       case InferStep::CAD_FULL: d_cadSlv.checkFull(); break;
       case InferStep::CAD_INIT: d_cadSlv.initLastCall(assertions); break;
       case InferStep::NL_FACTORING:
-        d_nlSlv.checkFactoring(assertions, false_asserts);
+        d_factoringSlv.check(assertions, false_asserts);
         break;
       case InferStep::IAND_INIT:
         d_iandSlv.initLastCall(assertions, false_asserts, xts);
@@ -698,30 +703,30 @@ void NonlinearExtension::runStrategy(Theory::Effort effort,
         d_icpSlv.check();
         break;
       case InferStep::NL_INIT:
-        d_nlSlv.initLastCall(assertions, false_asserts, xts);
+        d_extState.init(xts);
+        d_monomialBoundsSlv.init();
+        d_monomialSlv.init(xts);
         break;
       case InferStep::NL_MONOMIAL_INFER_BOUNDS:
-        d_nlSlv.checkMonomialInferBounds(assertions, false_asserts);
+        d_monomialBoundsSlv.checkBounds(assertions, false_asserts);
         break;
       case InferStep::NL_MONOMIAL_MAGNITUDE0:
-        d_nlSlv.checkMonomialMagnitude(0);
+        d_monomialSlv.checkMagnitude(0);
         break;
       case InferStep::NL_MONOMIAL_MAGNITUDE1:
-        d_nlSlv.checkMonomialMagnitude(1);
+        d_monomialSlv.checkMagnitude(1);
         break;
       case InferStep::NL_MONOMIAL_MAGNITUDE2:
-        d_nlSlv.checkMonomialMagnitude(2);
+        d_monomialSlv.checkMagnitude(2);
         break;
-      case InferStep::NL_MONOMIAL_SIGN: d_nlSlv.checkMonomialSign(); break;
+      case InferStep::NL_MONOMIAL_SIGN: d_monomialSlv.checkSign(); break;
       case InferStep::NL_RESOLUTION_BOUNDS:
-        d_nlSlv.checkMonomialInferResBounds();
-        break;
-      case InferStep::NL_SPLIT_ZERO: d_nlSlv.checkSplitZero(); break;
-      case InferStep::NL_TANGENT_PLANES:
-        d_nlSlv.checkTangentPlanes(false);
+        d_monomialBoundsSlv.checkResBounds();
         break;
+      case InferStep::NL_SPLIT_ZERO: d_splitZeroSlv.check(); break;
+      case InferStep::NL_TANGENT_PLANES: d_tangentPlaneSlv.check(false); break;
       case InferStep::NL_TANGENT_PLANES_WAITING:
-        d_nlSlv.checkTangentPlanes(true);
+        d_tangentPlaneSlv.check(true);
         break;
       case InferStep::TRANS_INIT:
         d_trSlv.initLastCall(assertions, false_asserts, xts);
diff --git a/src/theory/arith/nl/nonlinear_extension.h b/src/theory/arith/nl/nonlinear_extension.h
index cf45942c8..cd26a027f 100644
--- a/src/theory/arith/nl/nonlinear_extension.h
+++ b/src/theory/arith/nl/nonlinear_extension.h
@@ -28,12 +28,16 @@
 #include "expr/node.h"
 #include "theory/arith/inference_manager.h"
 #include "theory/arith/nl/cad_solver.h"
+#include "theory/arith/nl/ext/factoring_check.h"
+#include "theory/arith/nl/ext/monomial_bounds_check.h"
+#include "theory/arith/nl/ext/monomial_check.h"
+#include "theory/arith/nl/ext/split_zero_check.h"
+#include "theory/arith/nl/ext/tangent_plane_check.h"
 #include "theory/arith/nl/ext_theory_callback.h"
 #include "theory/arith/nl/iand_solver.h"
 #include "theory/arith/nl/icp/icp_solver.h"
 #include "theory/arith/nl/nl_lemma_utils.h"
 #include "theory/arith/nl/nl_model.h"
-#include "theory/arith/nl/nl_solver.h"
 #include "theory/arith/nl/stats.h"
 #include "theory/arith/nl/strategy.h"
 #include "theory/arith/nl/transcendental_solver.h"
@@ -256,12 +260,21 @@ 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, Cimatti et al., TACAS 2017.
+  /**
+   * Holds common lookup data for the checks implemented in the "nl-ext"
+   * solvers (from Cimatti et al., TACAS 2017).
    */
-  NlSolver d_nlSlv;
+  ExtState d_extState;
+  /** Solver for factoring lemmas. */
+  FactoringCheck d_factoringSlv;
+  /** Solver for lemmas about monomial bounds. */
+  MonomialBoundsCheck d_monomialBoundsSlv;
+  /** Solver for lemmas about monomials. */
+  MonomialCheck d_monomialSlv;
+  /** Solver for lemmas that split multiplication at zero. */
+  SplitZeroCheck d_splitZeroSlv;
+  /** Solver for tangent plane lemmas. */
+  TangentPlaneCheck d_tangentPlaneSlv;
   /** The CAD-based solver */
   CadSolver d_cadSlv;
   /** The ICP-based solver */