Refactor transcendental solver (#5514)

The transcendental solver has grown over time, and a refactoring was due.
This PR splits the transcendental solver into five components:

a utility to compute taylor approximations
a common state for transcendental solvers
a solver for exponential function
a solver for sine function
a solver that wraps everything to a transcendental solver.

15 files changed, 2491 insertions, 1901 deletions

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 366e72af0..70cb68431 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -387,8 +387,16 @@ libcvc4_add_sources(
   theory/arith/nl/stats.h
   theory/arith/nl/strategy.cpp
   theory/arith/nl/strategy.h
-  theory/arith/nl/transcendental_solver.cpp
-  theory/arith/nl/transcendental_solver.h
+  theory/arith/nl/transcendental/exponential_solver.cpp
+  theory/arith/nl/transcendental/exponential_solver.h
+  theory/arith/nl/transcendental/sine_solver.cpp
+  theory/arith/nl/transcendental/sine_solver.h
+  theory/arith/nl/transcendental/taylor_generator.cpp
+  theory/arith/nl/transcendental/taylor_generator.h
+  theory/arith/nl/transcendental/transcendental_solver.cpp
+  theory/arith/nl/transcendental/transcendental_solver.h
+  theory/arith/nl/transcendental/transcendental_state.cpp
+  theory/arith/nl/transcendental/transcendental_state.h
   theory/arith/normal_form.cpp
   theory/arith/normal_form.h
   theory/arith/operator_elim.cpp
diff --git a/src/theory/arith/nl/nl_lemma_utils.h b/src/theory/arith/nl/nl_lemma_utils.h
index ac52e7d52..24302339c 100644
--- a/src/theory/arith/nl/nl_lemma_utils.h
+++ b/src/theory/arith/nl/nl_lemma_utils.h
@@ -98,6 +98,25 @@ struct SortNlModel
   bool operator()(Node i, Node j);
 };
 
+/**
+ * Wrapper for std::sort that uses SortNlModel to sort an iterator range.
+ */
+template <typename It>
+void sortByNlModel(It begin,
+                   It end,
+                   NlModel* model,
+                   bool concrete = true,
+                   bool absolute = false,
+                   bool reverse = false)
+{
+  SortNlModel smv;
+  smv.d_nlm = model;
+  smv.d_isConcrete = concrete;
+  smv.d_isAbsolute = absolute;
+  smv.d_reverse_order = reverse;
+  std::sort(begin, end, smv);
+}
+
 struct SortNonlinearDegree
 {
   SortNonlinearDegree(const std::map<Node, unsigned>& m) : d_mdegree(m) {}
diff --git a/src/theory/arith/nl/nonlinear_extension.h b/src/theory/arith/nl/nonlinear_extension.h
index 21b978a55..b563384a5 100644
--- a/src/theory/arith/nl/nonlinear_extension.h
+++ b/src/theory/arith/nl/nonlinear_extension.h
@@ -38,7 +38,7 @@
 #include "theory/arith/nl/nl_model.h"
 #include "theory/arith/nl/stats.h"
 #include "theory/arith/nl/strategy.h"
-#include "theory/arith/nl/transcendental_solver.h"
+#include "theory/arith/nl/transcendental/transcendental_solver.h"
 #include "theory/ext_theory.h"
 #include "theory/uf/equality_engine.h"
 
@@ -253,7 +253,7 @@ class NonlinearExtension
    * This is the subsolver responsible for running the procedure for
    * transcendental functions.
    */
-  TranscendentalSolver d_trSlv;
+  transcendental::TranscendentalSolver d_trSlv;
   /**
    * Holds common lookup data for the checks implemented in the "nl-ext"
    * solvers (from Cimatti et al., TACAS 2017).
diff --git a/src/theory/arith/nl/transcendental/exponential_solver.cpp b/src/theory/arith/nl/transcendental/exponential_solver.cpp
new file mode 100644
index 000000000..4275b2487
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/exponential_solver.cpp
@@ -0,0 +1,207 @@
+/*********************                                                        */
+/*! \file transcendental_solver.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King, 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 solver for handling transcendental functions.
+ **/
+
+#include "theory/arith/nl/transcendental/exponential_solver.h"
+
+#include <cmath>
+#include <set>
+
+#include "expr/node_algorithm.h"
+#include "expr/node_builder.h"
+#include "options/arith_options.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/arith_utilities.h"
+#include "theory/rewriter.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+ExponentialSolver::ExponentialSolver(TranscendentalState* tstate)
+    : d_data(tstate)
+{
+}
+
+ExponentialSolver::~ExponentialSolver() {}
+
+void ExponentialSolver::checkInitialRefine()
+{
+  NodeManager* nm = NodeManager::currentNM();
+  Trace("nl-ext")
+      << "Get initial refinement lemmas for transcendental functions..."
+      << std::endl;
+  for (std::pair<const Kind, std::vector<Node> >& tfl : d_data->d_funcMap)
+  {
+    if (tfl.first != Kind::EXPONENTIAL)
+    {
+      continue;
+    }
+    Assert(tfl.first == Kind::EXPONENTIAL);
+    for (const Node& t : tfl.second)
+    {
+      // initial refinements
+      if (d_tf_initial_refine.find(t) == d_tf_initial_refine.end())
+      {
+        d_tf_initial_refine[t] = true;
+        Node lem;
+        // ( exp(x) > 0 ) ^ ( x=0 <=> exp( x ) = 1 ) ^ ( x < 0 <=> exp( x ) <
+        // 1 ) ^ ( x <= 0 V exp( x ) > x + 1 )
+        lem = nm->mkNode(
+            Kind::AND,
+            nm->mkNode(Kind::GT, t, d_data->d_zero),
+            nm->mkNode(Kind::EQUAL,
+                       t[0].eqNode(d_data->d_zero),
+                       t.eqNode(d_data->d_one)),
+            nm->mkNode(Kind::EQUAL,
+                       nm->mkNode(Kind::LT, t[0], d_data->d_zero),
+                       nm->mkNode(Kind::LT, t, d_data->d_one)),
+            nm->mkNode(
+                Kind::OR,
+                nm->mkNode(Kind::LEQ, t[0], d_data->d_zero),
+                nm->mkNode(
+                    Kind::GT, t, nm->mkNode(Kind::PLUS, t[0], d_data->d_one))));
+        if (!lem.isNull())
+        {
+          d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_INIT_REFINE);
+        }
+      }
+    }
+  }
+}
+
+void ExponentialSolver::checkMonotonic()
+{
+  Trace("nl-ext") << "Get monotonicity lemmas for transcendental functions..."
+                  << std::endl;
+
+  auto it = d_data->d_funcMap.find(Kind::EXPONENTIAL);
+  if (it == d_data->d_funcMap.end())
+  {
+    Trace("nl-ext-exp") << "No exponential terms" << std::endl;
+    return;
+  }
+
+  // sort arguments of all transcendentals
+  std::vector<Node> tf_args;
+  std::map<Node, Node> tf_arg_to_term;
+
+  for (const Node& tf : it->second)
+  {
+    Node a = tf[0];
+    Node mvaa = d_data->d_model.computeAbstractModelValue(a);
+    if (mvaa.isConst())
+    {
+      Trace("nl-ext-exp") << "...tf term : " << a << std::endl;
+      tf_args.push_back(a);
+      tf_arg_to_term[a] = tf;
+    }
+  }
+
+  if (tf_args.empty())
+  {
+    return;
+  }
+
+  sortByNlModel(
+      tf_args.begin(), tf_args.end(), &d_data->d_model, true, false, true);
+
+  Node targ, targval, t, tval;
+  for (const auto& sarg : tf_args)
+  {
+    Node sargval = d_data->d_model.computeAbstractModelValue(sarg);
+    Assert(sargval.isConst());
+    Node s = tf_arg_to_term[sarg];
+    Node sval = d_data->d_model.computeAbstractModelValue(s);
+    Assert(sval.isConst());
+
+    // store the concavity region
+    d_data->d_tf_region[s] = 1;
+    Trace("nl-ext-concavity") << ", arg model value = " << sargval << std::endl;
+
+    if (!tval.isNull() && sval.getConst<Rational>() > tval.getConst<Rational>())
+    {
+      NodeManager* nm = NodeManager::currentNM();
+      Node mono_lem = nm->mkNode(Kind::IMPLIES,
+                                 nm->mkNode(Kind::GEQ, targ, sarg),
+                                 nm->mkNode(Kind::GEQ, t, s));
+      Trace("nl-ext-exp") << "Monotonicity lemma : " << mono_lem << std::endl;
+
+      d_data->d_im.addPendingArithLemma(mono_lem,
+                                        InferenceId::NL_T_MONOTONICITY);
+    }
+    // store the previous values
+    targ = sarg;
+    targval = sargval;
+    t = s;
+    tval = sval;
+  }
+}
+
+void ExponentialSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx)
+{
+  NodeManager* nm = NodeManager::currentNM();
+  // compute tangent plane
+  // Figure 3: T( x )
+  // We use zero slope tangent planes, since the concavity of the Taylor
+  // approximation cannot be easily established.
+  // Tangent plane is valid in the interval [c,u).
+  Node lem = nm->mkNode(Kind::IMPLIES,
+                        nm->mkNode(Kind::GEQ, e[0], c),
+                        nm->mkNode(Kind::GEQ, e, poly_approx));
+  Trace("nl-ext-exp") << "*** Tangent plane lemma (pre-rewrite): " << lem
+                      << std::endl;
+  lem = Rewriter::rewrite(lem);
+  Trace("nl-ext-exp") << "*** Tangent plane lemma : " << lem << std::endl;
+  Assert(d_data->d_model.computeAbstractModelValue(lem) == d_data->d_false);
+  // Figure 3 : line 9
+  d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_TANGENT, nullptr, true);
+}
+
+void ExponentialSolver::doSecantLemmas(TNode e,
+                                       TNode c,
+                                       TNode poly_approx,
+                                       unsigned d)
+{
+  d_data->doSecantLemmas(getSecantBounds(e, c, d), c, poly_approx, e, d, 1);
+}
+
+std::pair<Node, Node> ExponentialSolver::getSecantBounds(TNode e,
+                                                         TNode c,
+                                                         unsigned d)
+{
+  std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, c, d);
+
+  // Check if we already have neighboring secant points
+  if (bounds.first.isNull())
+  {
+    // pick c-1
+    bounds.first = Rewriter::rewrite(
+        NodeManager::currentNM()->mkNode(Kind::MINUS, c, d_data->d_one));
+  }
+  if (bounds.second.isNull())
+  {
+    // pick c+1
+    bounds.second = Rewriter::rewrite(
+        NodeManager::currentNM()->mkNode(Kind::PLUS, c, d_data->d_one));
+  }
+  return bounds;
+}
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl/transcendental/exponential_solver.h b/src/theory/arith/nl/transcendental/exponential_solver.h
new file mode 100644
index 000000000..f757e5892
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/exponential_solver.h
@@ -0,0 +1,110 @@
+/*********************                                                        */
+/*! \file exponential_solver.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** 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 Solving for handling exponential function.
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__EXPONENTIAL_SOLVER_H
+#define CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__EXPONENTIAL_SOLVER_H
+
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+#include "expr/node.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+#include "theory/arith/nl/transcendental/transcendental_state.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+/** Transcendental solver class
+ *
+ * This class implements model-based refinement schemes
+ * for transcendental functions, described in:
+ *
+ * - "Satisfiability Modulo Transcendental
+ * Functions via Incremental Linearization" by Cimatti
+ * et al., CADE 2017.
+ *
+ * It's main functionality are methods that implement lemma schemas below,
+ * which return a set of lemmas that should be sent on the output channel.
+ */
+class ExponentialSolver
+{
+ public:
+  ExponentialSolver(TranscendentalState* tstate);
+  ~ExponentialSolver();
+
+  void initLastCall(const std::vector<Node>& assertions,
+                    const std::vector<Node>& false_asserts,
+                    const std::vector<Node>& xts);
+
+  /**
+   * check initial refine
+   *
+   * Constructs a set of valid theory lemmas, based on
+   * simple facts about the exponential function.
+   * This mostly follows the initial axioms described in
+   * Section 4 of "Satisfiability
+   * Modulo Transcendental Functions via Incremental
+   * Linearization" by Cimatti et al., CADE 2017.
+   *
+   * Examples:
+   *
+   * exp( x )>0
+   * x<0 => exp( x )<1
+   */
+  void checkInitialRefine();
+
+  /**
+   * check monotonicity
+   *
+   * Constructs a set of valid theory lemmas, based on a
+   * lemma scheme that ensures that applications
+   * of the exponential function respect monotonicity.
+   *
+   * Examples:
+   *
+   * x > y => exp( x ) > exp( y )
+   */
+  void checkMonotonic();
+
+  /** Sent tangent lemma around c for e */
+  void doTangentLemma(TNode e, TNode c, TNode poly_approx);
+
+  /** Sent secant lemmas around c for e */
+  void doSecantLemmas(TNode e, TNode c, TNode poly_approx, unsigned d);
+
+ private:
+  /** Generate bounds for secant lemmas */
+  std::pair<Node, Node> getSecantBounds(TNode e, TNode c, unsigned d);
+
+  /** Holds common state for transcendental solvers */
+  TranscendentalState* d_data;
+
+  /** The transcendental functions we have done initial refinements on */
+  std::map<Node, bool> d_tf_initial_refine;
+
+}; /* class ExponentialSolver */
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__TRANSCENDENTAL_SOLVER_H */
diff --git a/src/theory/arith/nl/transcendental/sine_solver.cpp b/src/theory/arith/nl/transcendental/sine_solver.cpp
new file mode 100644
index 000000000..012b8e7ec
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/sine_solver.cpp
@@ -0,0 +1,328 @@
+/*********************                                                        */
+/*! \file transcendental_solver.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King, 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 solver for handling transcendental functions.
+ **/
+
+#include "theory/arith/nl/transcendental/sine_solver.h"
+
+#include <cmath>
+#include <set>
+
+#include "expr/node_algorithm.h"
+#include "expr/node_builder.h"
+#include "options/arith_options.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/arith_utilities.h"
+#include "theory/rewriter.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+SineSolver::SineSolver(TranscendentalState* tstate) : d_data(tstate) {}
+
+SineSolver::~SineSolver() {}
+
+void SineSolver::checkInitialRefine()
+{
+  NodeManager* nm = NodeManager::currentNM();
+  Trace("nl-ext")
+      << "Get initial refinement lemmas for transcendental functions..."
+      << std::endl;
+  for (std::pair<const Kind, std::vector<Node> >& tfl : d_data->d_funcMap)
+  {
+    if (tfl.first != Kind::SINE)
+    {
+      continue;
+    }
+    for (const Node& t : tfl.second)
+    {
+      // initial refinements
+      if (d_tf_initial_refine.find(t) == d_tf_initial_refine.end())
+      {
+        d_tf_initial_refine[t] = true;
+        Node lem;
+        Node symn = nm->mkNode(Kind::SINE,
+                               nm->mkNode(Kind::MULT, d_data->d_neg_one, t[0]));
+        symn = Rewriter::rewrite(symn);
+        // Can assume it is its own master since phase is split over 0,
+        // hence  -pi <= t[0] <= pi implies -pi <= -t[0] <= pi.
+        d_data->d_trMaster[symn] = symn;
+        d_data->d_trSlaves[symn].insert(symn);
+        Assert(d_data->d_trSlaves.find(t) != d_data->d_trSlaves.end());
+        std::vector<Node> children;
+
+        lem =
+            nm->mkNode(Kind::AND,
+                       // bounds
+                       nm->mkNode(Kind::AND,
+                                  nm->mkNode(Kind::LEQ, t, d_data->d_one),
+                                  nm->mkNode(Kind::GEQ, t, d_data->d_neg_one)),
+                       // symmetry
+                       nm->mkNode(Kind::PLUS, t, symn).eqNode(d_data->d_zero),
+                       // sign
+                       nm->mkNode(Kind::EQUAL,
+                                  nm->mkNode(Kind::LT, t[0], d_data->d_zero),
+                                  nm->mkNode(Kind::LT, t, d_data->d_zero)),
+                       // zero val
+                       nm->mkNode(Kind::EQUAL,
+                                  nm->mkNode(Kind::GT, t[0], d_data->d_zero),
+                                  nm->mkNode(Kind::GT, t, d_data->d_zero)));
+        lem = nm->mkNode(
+            Kind::AND,
+            lem,
+            // zero tangent
+            nm->mkNode(Kind::AND,
+                       nm->mkNode(Kind::IMPLIES,
+                                  nm->mkNode(Kind::GT, t[0], d_data->d_zero),
+                                  nm->mkNode(Kind::LT, t, t[0])),
+                       nm->mkNode(Kind::IMPLIES,
+                                  nm->mkNode(Kind::LT, t[0], d_data->d_zero),
+                                  nm->mkNode(Kind::GT, t, t[0]))),
+            // pi tangent
+            nm->mkNode(
+                Kind::AND,
+                nm->mkNode(
+                    Kind::IMPLIES,
+                    nm->mkNode(Kind::LT, t[0], d_data->d_pi),
+                    nm->mkNode(Kind::LT,
+                               t,
+                               nm->mkNode(Kind::MINUS, d_data->d_pi, t[0]))),
+                nm->mkNode(
+                    Kind::IMPLIES,
+                    nm->mkNode(Kind::GT, t[0], d_data->d_pi_neg),
+                    nm->mkNode(
+                        Kind::GT,
+                        t,
+                        nm->mkNode(Kind::MINUS, d_data->d_pi_neg, t[0])))));
+        if (!lem.isNull())
+        {
+          d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_INIT_REFINE);
+        }
+      }
+    }
+  }
+}
+
+void SineSolver::checkMonotonic()
+{
+  Trace("nl-ext") << "Get monotonicity lemmas for transcendental functions..."
+                  << std::endl;
+
+  auto it = d_data->d_funcMap.find(Kind::SINE);
+  if (it == d_data->d_funcMap.end())
+  {
+    Trace("nl-ext-exp") << "No sine terms" << std::endl;
+    return;
+  }
+
+  // sort arguments of all transcendentals
+  std::vector<Node> tf_args;
+  std::map<Node, Node> tf_arg_to_term;
+
+  for (const Node& tf : it->second)
+  {
+    Node a = tf[0];
+    Node mvaa = d_data->d_model.computeAbstractModelValue(a);
+    if (mvaa.isConst())
+    {
+      Trace("nl-ext-tf-mono-debug") << "...tf term : " << a << std::endl;
+      tf_args.push_back(a);
+      tf_arg_to_term[a] = tf;
+    }
+  }
+
+  if (tf_args.empty())
+  {
+    return;
+  }
+
+  sortByNlModel(
+      tf_args.begin(), tf_args.end(), &d_data->d_model, true, false, true);
+
+  std::vector<Node> mpoints = {d_data->d_pi,
+                               d_data->d_pi_2,
+                               d_data->d_zero,
+                               d_data->d_pi_neg_2,
+                               d_data->d_pi_neg};
+  std::vector<Node> mpoints_vals;
+
+  // get model values for points
+  for (const auto& point : mpoints)
+  {
+    mpoints_vals.emplace_back(d_data->d_model.computeAbstractModelValue(point));
+    Assert(mpoints_vals.back().isConst());
+  }
+
+  unsigned mdir_index = 0;
+  int monotonic_dir = -1;
+  Node mono_bounds[2];
+  Node targ, targval, t, tval;
+  for (const auto& sarg : tf_args)
+  {
+    Node sargval = d_data->d_model.computeAbstractModelValue(sarg);
+    Assert(sargval.isConst());
+    Node s = tf_arg_to_term[sarg];
+    Node sval = d_data->d_model.computeAbstractModelValue(s);
+    Assert(sval.isConst());
+
+    // increment to the proper monotonicity region
+    bool increment = true;
+    while (increment && mdir_index < mpoints.size())
+    {
+      increment = false;
+      Node pval = mpoints_vals[mdir_index];
+      Assert(pval.isConst());
+      if (sargval.getConst<Rational>() < pval.getConst<Rational>())
+      {
+        increment = true;
+        Trace("nl-ext-tf-mono")
+            << "...increment at " << sarg << " since model value is less than "
+            << mpoints[mdir_index] << std::endl;
+      }
+      if (increment)
+      {
+        tval = Node::null();
+        mono_bounds[1] = mpoints[mdir_index];
+        mdir_index++;
+        monotonic_dir = regionToMonotonicityDir(mdir_index);
+        if (mdir_index < mpoints.size())
+        {
+          mono_bounds[0] = mpoints[mdir_index];
+        }
+        else
+        {
+          mono_bounds[0] = Node::null();
+        }
+      }
+    }
+    // store the concavity region
+    d_data->d_tf_region[s] = mdir_index;
+    Trace("nl-ext-concavity")
+        << "Transcendental function " << s << " is in region #" << mdir_index;
+    Trace("nl-ext-concavity") << ", arg model value = " << sargval << std::endl;
+
+    if (!tval.isNull())
+    {
+      NodeManager* nm = NodeManager::currentNM();
+      Node mono_lem;
+      if (monotonic_dir == 1
+          && sval.getConst<Rational>() > tval.getConst<Rational>())
+      {
+        mono_lem = nm->mkNode(Kind::IMPLIES,
+                              nm->mkNode(Kind::GEQ, targ, sarg),
+                              nm->mkNode(Kind::GEQ, t, s));
+      }
+      else if (monotonic_dir == -1
+               && sval.getConst<Rational>() < tval.getConst<Rational>())
+      {
+        mono_lem = nm->mkNode(Kind::IMPLIES,
+                              nm->mkNode(Kind::LEQ, targ, sarg),
+                              nm->mkNode(Kind::LEQ, t, s));
+      }
+      if (!mono_lem.isNull())
+      {
+        if (!mono_bounds[0].isNull())
+        {
+          Assert(!mono_bounds[1].isNull());
+          mono_lem = nm->mkNode(
+              Kind::IMPLIES,
+              nm->mkNode(Kind::AND,
+                         mkBounded(mono_bounds[0], targ, mono_bounds[1]),
+                         mkBounded(mono_bounds[0], sarg, mono_bounds[1])),
+              mono_lem);
+        }
+        Trace("nl-ext-tf-mono")
+            << "Monotonicity lemma : " << mono_lem << std::endl;
+
+        d_data->d_im.addPendingArithLemma(mono_lem,
+                                          InferenceId::NL_T_MONOTONICITY);
+      }
+    }
+    // store the previous values
+    targ = sarg;
+    targval = sargval;
+    t = s;
+    tval = sval;
+  }
+}
+
+void SineSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx, int region)
+{
+  NodeManager* nm = NodeManager::currentNM();
+
+  // compute tangent plane
+  // Figure 3: T( x )
+  // We use zero slope tangent planes, since the concavity of the Taylor
+  // approximation cannot be easily established.
+  int concavity = regionToConcavity(region);
+  int mdir = regionToMonotonicityDir(region);
+  Node lem = nm->mkNode(
+      Kind::IMPLIES,
+      nm->mkNode(
+          Kind::AND,
+          nm->mkNode(Kind::GEQ,
+                     e[0],
+                     mdir == concavity ? Node(c) : regionToLowerBound(region)),
+          nm->mkNode(Kind::LEQ,
+                     e[0],
+                     mdir != concavity ? Node(c) : regionToUpperBound(region))),
+      nm->mkNode(concavity == 1 ? Kind::GEQ : Kind::LEQ, e, poly_approx));
+
+  Trace("nl-ext-sine") << "*** Tangent plane lemma (pre-rewrite): " << lem
+                       << std::endl;
+  lem = Rewriter::rewrite(lem);
+  Trace("nl-ext-sine") << "*** Tangent plane lemma : " << lem << std::endl;
+  Assert(d_data->d_model.computeAbstractModelValue(lem) == d_data->d_false);
+  // Figure 3 : line 9
+  d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_TANGENT, nullptr, true);
+}
+
+void SineSolver::doSecantLemmas(
+    TNode e, TNode c, TNode poly_approx, unsigned d, int region)
+{
+  d_data->doSecantLemmas(getSecantBounds(e, c, d, region),
+                         c,
+                         poly_approx,
+                         e,
+                         d,
+                         regionToConcavity(region));
+}
+
+std::pair<Node, Node> SineSolver::getSecantBounds(TNode e,
+                                                  TNode c,
+                                                  unsigned d,
+                                                  int region)
+{
+  std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, c, d);
+
+  // Check if we already have neighboring secant points
+  if (bounds.first.isNull())
+  {
+    // lower boundary point for this concavity region
+    bounds.first = regionToLowerBound(region);
+  }
+  if (bounds.second.isNull())
+  {
+    // upper boundary point for this concavity region
+    bounds.second = regionToUpperBound(region);
+  }
+  return bounds;
+}
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl/transcendental/sine_solver.h b/src/theory/arith/nl/transcendental/sine_solver.h
new file mode 100644
index 000000000..394c2d927
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/sine_solver.h
@@ -0,0 +1,177 @@
+/*********************                                                        */
+/*! \file sine_solver.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** 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 Solving for handling exponential function.
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__SINE_SOLVER_H
+#define CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__SINE_SOLVER_H
+
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+#include "expr/node.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+#include "theory/arith/nl/transcendental/transcendental_state.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+/** Transcendental solver class
+ *
+ * This class implements model-based refinement schemes
+ * for transcendental functions, described in:
+ *
+ * - "Satisfiability Modulo Transcendental
+ * Functions via Incremental Linearization" by Cimatti
+ * et al., CADE 2017.
+ *
+ * It's main functionality are methods that implement lemma schemas below,
+ * which return a set of lemmas that should be sent on the output channel.
+ */
+class SineSolver
+{
+ public:
+  SineSolver(TranscendentalState* tstate);
+  ~SineSolver();
+
+  void initLastCall(const std::vector<Node>& assertions,
+                    const std::vector<Node>& false_asserts,
+                    const std::vector<Node>& xts);
+
+  /**
+   * check initial refine
+   *
+   * Constructs a set of valid theory lemmas, based on
+   * simple facts about the sine function.
+   * This mostly follows the initial axioms described in
+   * Section 4 of "Satisfiability
+   * Modulo Transcendental Functions via Incremental
+   * Linearization" by Cimatti et al., CADE 2017.
+   *
+   * Examples:
+   *
+   * sin( x ) = -sin( -x )
+   * ( PI > x > 0 ) => 0 < sin( x ) < 1
+   */
+  void checkInitialRefine();
+
+  /**
+   * check monotonicity
+   *
+   * Constructs a set of valid theory lemmas, based on a
+   * lemma scheme that ensures that applications
+   * of the sine function respect monotonicity.
+   *
+   * Examples:
+   *
+   * PI/2 > x > y > 0 => sin( x ) > sin( y )
+   * PI > x > y > PI/2 => sin( x ) < sin( y )
+   */
+  void checkMonotonic();
+
+  /** Sent tangent lemma around c for e */
+  void doTangentLemma(TNode e, TNode c, TNode poly_approx, int region);
+
+  /** Sent secant lemmas around c for e */
+  void doSecantLemmas(
+      TNode e, TNode c, TNode poly_approx, unsigned d, int region);
+
+ private:
+  std::pair<Node, Node> getSecantBounds(TNode e,
+                                        TNode c,
+                                        unsigned d,
+                                        int region);
+
+  /** region to lower bound
+   *
+   * Returns the term corresponding to the lower
+   * bound of the region of transcendental function
+   * with kind k. Returns Node::null if the region
+   * is invalid, or there is no lower bound for the
+   * region.
+   */
+  Node regionToLowerBound(int region)
+  {
+    switch (region)
+    {
+      case 1: return d_data->d_pi_2;
+      case 2: return d_data->d_zero;
+      case 3: return d_data->d_pi_neg_2;
+      case 4: return d_data->d_pi_neg;
+      default: return Node();
+    }
+  }
+
+  /** region to upper bound
+   *
+   * Returns the term corresponding to the upper
+   * bound of the region of transcendental function
+   * with kind k. Returns Node::null if the region
+   * is invalid, or there is no upper bound for the
+   * region.
+   */
+  Node regionToUpperBound(int region)
+  {
+    switch (region)
+    {
+      case 1: return d_data->d_pi;
+      case 2: return d_data->d_pi_2;
+      case 3: return d_data->d_zero;
+      case 4: return d_data->d_pi_neg_2;
+      default: return Node();
+    }
+  }
+
+  int regionToMonotonicityDir(int region)
+  {
+    switch (region)
+    {
+      case 1:
+      case 4: return -1;
+      case 2:
+      case 3: return 1;
+      default: return 0;
+    }
+  }
+  int regionToConcavity(int region)
+  {
+    switch (region)
+    {
+      case 1:
+      case 2: return -1;
+      case 3:
+      case 4: return 1;
+      default: return 0;
+    }
+  }
+
+  /** Holds common state for transcendental solvers */
+  TranscendentalState* d_data;
+
+  /** The transcendental functions we have done initial refinements on */
+  std::map<Node, bool> d_tf_initial_refine;
+
+}; /* class SineSolver */
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__TRANSCENDENTAL_SOLVER_H */
diff --git a/src/theory/arith/nl/transcendental/taylor_generator.cpp b/src/theory/arith/nl/transcendental/taylor_generator.cpp
new file mode 100644
index 000000000..a373339e9
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/taylor_generator.cpp
@@ -0,0 +1,317 @@
+/*********************                                                        */
+/*! \file taylor_generator.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** 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 Generate taylor approximations transcendental lemmas.
+ **/
+
+#include "theory/arith/nl/transcendental/taylor_generator.h"
+
+#include "theory/arith/arith_utilities.h"
+#include "theory/rewriter.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+TaylorGenerator::TaylorGenerator()
+    : d_nm(NodeManager::currentNM()),
+      d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType())),
+      d_taylor_real_fv_base(d_nm->mkBoundVar("a", d_nm->realType())),
+      d_taylor_real_fv_base_rem(d_nm->mkBoundVar("b", d_nm->realType()))
+{
+}
+
+TNode TaylorGenerator::getTaylorVariable() { return d_taylor_real_fv; }
+
+std::pair<Node, Node> TaylorGenerator::getTaylor(TNode fa, std::uint64_t n)
+{
+  NodeManager* nm = NodeManager::currentNM();
+  Node d_zero = nm->mkConst(Rational(0));
+
+  Assert(n > 0);
+  Node fac;  // what term we cache for fa
+  if (fa[0] == d_zero)
+  {
+    // optimization : simpler to compute (x-fa[0])^n if we are centered around 0
+    fac = fa;
+  }
+  else
+  {
+    // otherwise we use a standard factor a in (x-a)^n
+    fac = nm->mkNode(fa.getKind(), d_taylor_real_fv_base);
+  }
+  Node taylor_rem;
+  Node taylor_sum;
+  // check if we have already computed this Taylor series
+  auto itt = s_taylor_sum[fac].find(n);
+  if (itt == s_taylor_sum[fac].end())
+  {
+    Node i_exp_base;
+    if (fa[0] == d_zero)
+    {
+      i_exp_base = d_taylor_real_fv;
+    }
+    else
+    {
+      i_exp_base = Rewriter::rewrite(
+          nm->mkNode(Kind::MINUS, d_taylor_real_fv, d_taylor_real_fv_base));
+    }
+    Node i_derv = fac;
+    Node i_fact = nm->mkConst(Rational(1));
+    Node i_exp = i_fact;
+    int i_derv_status = 0;
+    unsigned counter = 0;
+    std::vector<Node> sum;
+    do
+    {
+      counter++;
+      if (fa.getKind() == Kind::EXPONENTIAL)
+      {
+        // unchanged
+      }
+      else if (fa.getKind() == Kind::SINE)
+      {
+        if (i_derv_status % 2 == 1)
+        {
+          Node pi = nm->mkNullaryOperator(nm->realType(), Kind::PI);
+          Node pi_2 = Rewriter::rewrite(nm->mkNode(
+              Kind::MULT, pi, nm->mkConst(Rational(1) / Rational(2))));
+
+          Node arg = nm->mkNode(Kind::PLUS, pi_2, d_taylor_real_fv_base);
+          i_derv = nm->mkNode(Kind::SINE, arg);
+        }
+        else
+        {
+          i_derv = fa;
+        }
+        if (i_derv_status >= 2)
+        {
+          i_derv = nm->mkNode(Kind::MINUS, d_zero, i_derv);
+        }
+        i_derv = Rewriter::rewrite(i_derv);
+        i_derv_status = i_derv_status == 3 ? 0 : i_derv_status + 1;
+      }
+      if (counter == (n + 1))
+      {
+        TNode x = d_taylor_real_fv_base;
+        i_derv = i_derv.substitute(x, d_taylor_real_fv_base_rem);
+      }
+      Node curr = nm->mkNode(
+          Kind::MULT, nm->mkNode(Kind::DIVISION, i_derv, i_fact), i_exp);
+      if (counter == (n + 1))
+      {
+        taylor_rem = curr;
+      }
+      else
+      {
+        sum.push_back(curr);
+        i_fact = Rewriter::rewrite(
+            nm->mkNode(Kind::MULT, nm->mkConst(Rational(counter)), i_fact));
+        i_exp = Rewriter::rewrite(nm->mkNode(Kind::MULT, i_exp_base, i_exp));
+      }
+    } while (counter <= n);
+    taylor_sum = sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum);
+
+    if (fac[0] != d_taylor_real_fv_base)
+    {
+      TNode x = d_taylor_real_fv_base;
+      taylor_sum = taylor_sum.substitute(x, fac[0]);
+    }
+
+    // cache
+    s_taylor_sum[fac][n] = taylor_sum;
+    d_taylor_rem[fac][n] = taylor_rem;
+  }
+  else
+  {
+    taylor_sum = itt->second;
+    Assert(d_taylor_rem[fac].find(n) != d_taylor_rem[fac].end());
+    taylor_rem = d_taylor_rem[fac][n];
+  }
+
+  // must substitute for the argument if we were using a different lookup
+  if (fa[0] != fac[0])
+  {
+    TNode x = d_taylor_real_fv_base;
+    taylor_sum = taylor_sum.substitute(x, fa[0]);
+  }
+  return std::pair<Node, Node>(taylor_sum, taylor_rem);
+}
+
+void TaylorGenerator::getPolynomialApproximationBounds(
+    Kind k, unsigned d, std::vector<Node>& pbounds)
+{
+  if (d_poly_bounds[k][d].empty())
+  {
+    NodeManager* nm = NodeManager::currentNM();
+    Node tft = nm->mkNode(k, nm->mkConst(Rational(0)));
+    // n is the Taylor degree we are currently considering
+    unsigned n = 2 * d;
+    // n must be even
+    std::pair<Node, Node> taylor = getTaylor(tft, n);
+    Trace("nl-ext-tftp-debug2")
+        << "Taylor for " << k << " is : " << taylor.first << std::endl;
+    Node taylor_sum = Rewriter::rewrite(taylor.first);
+    Trace("nl-ext-tftp-debug2")
+        << "Taylor for " << k << " is (post-rewrite) : " << taylor_sum
+        << std::endl;
+    Assert(taylor.second.getKind() == Kind::MULT);
+    Assert(taylor.second.getNumChildren() == 2);
+    Assert(taylor.second[0].getKind() == Kind::DIVISION);
+    Trace("nl-ext-tftp-debug2")
+        << "Taylor remainder for " << k << " is " << taylor.second << std::endl;
+    // ru is x^{n+1}/(n+1)!
+    Node ru = nm->mkNode(Kind::DIVISION, taylor.second[1], taylor.second[0][1]);
+    ru = Rewriter::rewrite(ru);
+    Trace("nl-ext-tftp-debug2")
+        << "Taylor remainder factor is (post-rewrite) : " << ru << std::endl;
+    if (k == Kind::EXPONENTIAL)
+    {
+      pbounds.push_back(taylor_sum);
+      pbounds.push_back(taylor_sum);
+      pbounds.push_back(Rewriter::rewrite(
+          nm->mkNode(Kind::MULT,
+                     taylor_sum,
+                     nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru))));
+      pbounds.push_back(
+          Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru)));
+    }
+    else
+    {
+      Assert(k == Kind::SINE);
+      Node l = Rewriter::rewrite(nm->mkNode(Kind::MINUS, taylor_sum, ru));
+      Node u = Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
+      pbounds.push_back(l);
+      pbounds.push_back(l);
+      pbounds.push_back(u);
+      pbounds.push_back(u);
+    }
+    Trace("nl-ext-tf-tplanes")
+        << "Polynomial approximation for " << k << " is: " << std::endl;
+    Trace("nl-ext-tf-tplanes") << " Lower (pos): " << pbounds[0] << std::endl;
+    Trace("nl-ext-tf-tplanes") << " Upper (pos): " << pbounds[2] << std::endl;
+    Trace("nl-ext-tf-tplanes") << " Lower (neg): " << pbounds[1] << std::endl;
+    Trace("nl-ext-tf-tplanes") << " Upper (neg): " << pbounds[3] << std::endl;
+    d_poly_bounds[k][d].insert(
+        d_poly_bounds[k][d].end(), pbounds.begin(), pbounds.end());
+  }
+  else
+  {
+    pbounds.insert(
+        pbounds.end(), d_poly_bounds[k][d].begin(), d_poly_bounds[k][d].end());
+  }
+}
+
+void TaylorGenerator::getPolynomialApproximationBoundForArg(
+    Kind k, Node c, unsigned d, std::vector<Node>& pbounds)
+{
+  getPolynomialApproximationBounds(k, d, pbounds);
+  Assert(c.isConst());
+  if (k == Kind::EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
+  {
+    NodeManager* nm = NodeManager::currentNM();
+    Node tft = nm->mkNode(k, nm->mkConst(Rational(0)));
+    bool success = false;
+    unsigned ds = d;
+    TNode ttrf = getTaylorVariable();
+    TNode tc = c;
+    do
+    {
+      success = true;
+      unsigned n = 2 * ds;
+      std::pair<Node, Node> taylor = getTaylor(tft, n);
+      // check that 1-c^{n+1}/(n+1)! > 0
+      Node ru =
+          nm->mkNode(Kind::DIVISION, taylor.second[1], taylor.second[0][1]);
+      Node rus = ru.substitute(ttrf, tc);
+      rus = Rewriter::rewrite(rus);
+      Assert(rus.isConst());
+      if (rus.getConst<Rational>() > 1)
+      {
+        success = false;
+        ds = ds + 1;
+      }
+    } while (!success);
+    if (ds > d)
+    {
+      Trace("nl-ext-exp-taylor")
+          << "*** Increase Taylor bound to " << ds << " > " << d << " for ("
+          << k << " " << c << ")" << std::endl;
+      // must use sound upper bound
+      std::vector<Node> pboundss;
+      getPolynomialApproximationBounds(k, ds, pboundss);
+      pbounds[2] = pboundss[2];
+    }
+  }
+}
+
+std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlModel& model)
+{
+  // compute the model value of the argument
+  Node c = model.computeAbstractModelValue(tf[0]);
+  Assert(c.isConst());
+  int csign = c.getConst<Rational>().sgn();
+  Kind k = tf.getKind();
+  if (csign == 0)
+  {
+    NodeManager* nm = NodeManager::currentNM();
+    // at zero, its trivial
+    if (k == Kind::SINE)
+    {
+      Node zero = nm->mkConst(Rational(0));
+      return std::pair<Node, Node>(zero, zero);
+    }
+    Assert(k == Kind::EXPONENTIAL);
+    Node one = nm->mkConst(Rational(1));
+    return std::pair<Node, Node>(one, one);
+  }
+  bool isNeg = csign == -1;
+
+  std::vector<Node> pbounds;
+  getPolynomialApproximationBoundForArg(k, c, d, pbounds);
+
+  std::vector<Node> bounds;
+  TNode tfv = getTaylorVariable();
+  TNode tfs = tf[0];
+  for (unsigned d2 = 0; d2 < 2; d2++)
+  {
+    int index = d2 == 0 ? (isNeg ? 1 : 0) : (isNeg ? 3 : 2);
+    Node pab = pbounds[index];
+    if (!pab.isNull())
+    {
+      // { x -> M_A(tf[0]) }
+      // Notice that we compute the model value of tfs first, so that
+      // the call to rewrite below does not modify the term, where notice that
+      // rewrite( x*x { x -> M_A(t) } ) = M_A(t)*M_A(t)
+      // is not equal to
+      // M_A( x*x { x -> t } ) = M_A( t*t )
+      // where M_A denotes the abstract model.
+      Node mtfs = model.computeAbstractModelValue(tfs);
+      pab = pab.substitute(tfv, mtfs);
+      pab = Rewriter::rewrite(pab);
+      Assert(pab.isConst());
+      bounds.push_back(pab);
+    }
+    else
+    {
+      bounds.push_back(Node::null());
+    }
+  }
+  return std::pair<Node, Node>(bounds[0], bounds[1]);
+}
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl/transcendental/taylor_generator.h b/src/theory/arith/nl/transcendental/taylor_generator.h
new file mode 100644
index 000000000..c647f7fd2
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/taylor_generator.h
@@ -0,0 +1,118 @@
+/*********************                                                        */
+/*! \file taylor_generator.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** 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 Generate taylor approximations transcendental lemmas.
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TAYLOR_GENERATOR_H
+#define CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TAYLOR_GENERATOR_H
+
+#include "expr/node.h"
+#include "theory/arith/nl/nl_model.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+class TaylorGenerator
+{
+ public:
+  TaylorGenerator();
+
+  /**
+   * Return the variable used as x in getTaylor().
+   */
+  TNode getTaylorVariable();
+
+  /**
+   * Get Taylor series of degree n for function fa centered around point fa[0].
+   *
+   * Return value is ( P_{n,f(a)}( x ), R_{n+1,f(a)}( x ) ) where
+   * the first part of the pair is the Taylor series expansion :
+   *    P_{n,f(a)}( x ) = sum_{i=0}^n (f^i( a )/i!)*(x-a)^i
+   * and the second part of the pair is the Taylor series remainder :
+   *    R_{n+1,f(a),b}( x ) = (f^{n+1}( b )/(n+1)!)*(x-a)^{n+1}
+   *
+   * The above values are cached for each (f,n) for a fixed variable "a".
+   * To compute the Taylor series for fa, we compute the Taylor series
+   *   for ( fa.getKind(), n ) then substitute { a -> fa[0] } if fa[0]!=0.
+   * We compute P_{n,f(0)}( x )/R_{n+1,f(0),b}( x ) for ( fa.getKind(), n )
+   *   if fa[0]=0.
+   * In the latter case, note we compute the exponential x^{n+1}
+   * instead of (x-a)^{n+1}, which can be done faster.
+   */
+  std::pair<Node, Node> getTaylor(TNode fa, std::uint64_t n);
+
+  /**
+   * polynomial approximation bounds
+   *
+   * This adds P_l+[x], P_l-[x], P_u+[x], P_u-[x] to pbounds, where x is
+   * d_taylor_real_fv. These are polynomial approximations of the Taylor series
+   * of <k>( 0 ) for degree 2*d where k is SINE or EXPONENTIAL.
+   * These correspond to P_l and P_u from Figure 3 of Cimatti et al., CADE 2017,
+   * for positive/negative (+/-) values of the argument of <k>( 0 ).
+   *
+   * Notice that for certain bounds (e.g. upper bounds for exponential), the
+   * Taylor approximation for a fixed degree is only sound up to a given
+   * upper bound on the argument. To obtain sound lower/upper bounds for a
+   * given <k>( c ), use the function below.
+   */
+  void getPolynomialApproximationBounds(Kind k,
+                                        unsigned d,
+                                        std::vector<Node>& pbounds);
+
+  /**
+   * polynomial approximation bounds
+   *
+   * This computes polynomial approximations P_l+[x], P_l-[x], P_u+[x], P_u-[x]
+   * that are sound (lower, upper) bounds for <k>( c ). Notice that these
+   * polynomials may depend on c. In particular, for P_u+[x] for <k>( c ) where
+   * c>0, we return the P_u+[x] from the function above for the minimum degree
+   * d' >= d such that (1-c^{2*d'+1}/(2*d'+1)!) is positive.
+   */
+  void getPolynomialApproximationBoundForArg(Kind k,
+                                             Node c,
+                                             unsigned d,
+                                             std::vector<Node>& pbounds);
+
+  /** get transcendental function model bounds
+   *
+   * This returns the current lower and upper bounds of transcendental
+   * function application tf based on Taylor of degree 2*d, which is dependent
+   * on the model value of its argument.
+   */
+  std::pair<Node, Node> getTfModelBounds(Node tf, unsigned d, NlModel& model);
+
+ private:
+  NodeManager* d_nm;
+  const Node d_taylor_real_fv;
+  const Node d_taylor_real_fv_base;
+  const Node d_taylor_real_fv_base_rem;
+  std::unordered_map<Node,
+                     std::unordered_map<std::uint64_t, Node>,
+                     NodeHashFunction>
+      s_taylor_sum;
+  std::unordered_map<Node,
+                     std::unordered_map<std::uint64_t, Node>,
+                     NodeHashFunction>
+      d_taylor_rem;
+  std::map<Kind, std::map<unsigned, std::vector<Node>>> d_poly_bounds;
+};
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__TRANSCENDENTAL_SOLVER_H */
diff --git a/src/theory/arith/nl/transcendental/transcendental_solver.cpp b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
new file mode 100644
index 000000000..1f76cd833
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
@@ -0,0 +1,402 @@
+/*********************                                                        */
+/*! \file transcendental_solver.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King, 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 solver for handling transcendental functions.
+ **/
+
+#include "theory/arith/nl/transcendental/transcendental_solver.h"
+
+#include <cmath>
+#include <set>
+
+#include "expr/node_algorithm.h"
+#include "expr/node_builder.h"
+#include "options/arith_options.h"
+#include "theory/arith/arith_msum.h"
+#include "theory/arith/arith_utilities.h"
+#include "theory/arith/nl/transcendental/taylor_generator.h"
+#include "theory/rewriter.h"
+
+using namespace CVC4::kind;
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+TranscendentalSolver::TranscendentalSolver(InferenceManager& im, NlModel& m)
+    : d_tstate(im, m), d_expSlv(&d_tstate), d_sineSlv(&d_tstate)
+{
+  d_taylor_degree = options::nlExtTfTaylorDegree();
+}
+
+TranscendentalSolver::~TranscendentalSolver() {}
+
+void TranscendentalSolver::initLastCall(const std::vector<Node>& assertions,
+                                        const std::vector<Node>& false_asserts,
+                                        const std::vector<Node>& xts)
+{
+  d_tstate.init(assertions, false_asserts, xts);
+}
+
+bool TranscendentalSolver::preprocessAssertionsCheckModel(
+    std::vector<Node>& assertions)
+{
+  std::vector<Node> pvars;
+  std::vector<Node> psubs;
+  for (const std::pair<const Node, Node>& tb : d_tstate.d_trMaster)
+  {
+    pvars.push_back(tb.first);
+    psubs.push_back(tb.second);
+  }
+
+  // initialize representation of assertions
+  std::vector<Node> passertions;
+  for (const Node& a : assertions)
+
+  {
+    Node pa = a;
+    if (!pvars.empty())
+    {
+      pa = arithSubstitute(pa, pvars, psubs);
+      pa = Rewriter::rewrite(pa);
+    }
+    if (!pa.isConst() || !pa.getConst<bool>())
+    {
+      Trace("nl-ext-cm-assert") << "- assert : " << pa << std::endl;
+      passertions.push_back(pa);
+    }
+  }
+  // get model bounds for all transcendental functions
+  Trace("nl-ext-cm") << "----- Get bounds for transcendental functions..."
+                     << std::endl;
+  for (std::pair<const Kind, std::vector<Node> >& tfs : d_tstate.d_funcMap)
+  {
+    for (const Node& tf : tfs.second)
+    {
+      Trace("nl-ext-cm") << "- Term: " << tf << std::endl;
+      bool success = true;
+      // tf is Figure 3 : tf( x )
+      std::pair<Node, Node> bounds;
+      if (tfs.first == Kind::PI)
+      {
+        bounds = {d_tstate.d_pi_bound[0], d_tstate.d_pi_bound[1]};
+      }
+      else
+      {
+        bounds = d_tstate.d_taylor.getTfModelBounds(
+            tf, d_taylor_degree, d_tstate.d_model);
+        if (bounds.first != bounds.second)
+        {
+          d_tstate.d_model.setUsedApproximate();
+        }
+      }
+      if (!bounds.first.isNull() && !bounds.second.isNull())
+      {
+        // for each function in the congruence classe
+        for (const Node& ctf : d_tstate.d_funcCongClass[tf])
+        {
+          // each term in congruence classes should be master terms
+          Assert(d_tstate.d_trSlaves.find(ctf) != d_tstate.d_trSlaves.end());
+          // we set the bounds for each slave of tf
+          for (const Node& stf : d_tstate.d_trSlaves[ctf])
+          {
+            Trace("nl-ext-cm")
+                << "...bound for " << stf << " : [" << bounds.first << ", "
+                << bounds.second << "]" << std::endl;
+            success = d_tstate.d_model.addCheckModelBound(
+                stf, bounds.first, bounds.second);
+          }
+        }
+      }
+      else
+      {
+        Trace("nl-ext-cm") << "...no bound for " << tf << std::endl;
+      }
+      if (!success)
+      {
+        // a bound was conflicting
+        Trace("nl-ext-cm") << "...failed to set bound for " << tf << std::endl;
+        Trace("nl-ext-cm") << "-----" << std::endl;
+        return false;
+      }
+    }
+  }
+  // replace the assertions
+  assertions = std::move(passertions);
+  return true;
+}
+
+void TranscendentalSolver::incrementTaylorDegree() { d_taylor_degree++; }
+unsigned TranscendentalSolver::getTaylorDegree() const
+{
+  return d_taylor_degree;
+}
+
+void TranscendentalSolver::processSideEffect(const NlLemma& se)
+{
+  for (const std::tuple<Node, unsigned, Node>& sp : se.d_secantPoint)
+  {
+    Node tf = std::get<0>(sp);
+    unsigned d = std::get<1>(sp);
+    Node c = std::get<2>(sp);
+    d_tstate.d_secant_points[tf][d].push_back(c);
+  }
+}
+
+void TranscendentalSolver::checkTranscendentalInitialRefine()
+{
+  d_expSlv.checkInitialRefine();
+  d_sineSlv.checkInitialRefine();
+}
+
+void TranscendentalSolver::checkTranscendentalMonotonic()
+{
+  d_expSlv.checkMonotonic();
+  d_sineSlv.checkMonotonic();
+}
+
+void TranscendentalSolver::checkTranscendentalTangentPlanes()
+{
+  Trace("nl-ext") << "Get tangent plane lemmas for transcendental functions..."
+                  << std::endl;
+  // this implements Figure 3 of "Satisfiaility Modulo Transcendental Functions
+  // via Incremental Linearization" by Cimatti et al
+  for (const std::pair<const Kind, std::vector<Node> >& tfs :
+       d_tstate.d_funcMap)
+  {
+    Kind k = tfs.first;
+    if (k == PI)
+    {
+      // We do not use Taylor approximation for PI currently.
+      // This is because the convergence is extremely slow, and hence an
+      // initial approximation is superior.
+      continue;
+    }
+
+    // we substitute into the Taylor sum P_{n,f(0)}( x )
+
+    for (const Node& tf : tfs.second)
+    {
+      // tf is Figure 3 : tf( x )
+      Trace("nl-ext-tftp") << "Compute tangent planes " << tf << std::endl;
+      // go until max degree is reached, or we don't meet bound criteria
+      for (unsigned d = 1; d <= d_taylor_degree; d++)
+      {
+        Trace("nl-ext-tftp") << "- run at degree " << d << "..." << std::endl;
+        unsigned prev =
+            d_tstate.d_im.numPendingLemmas() + d_tstate.d_im.numWaitingLemmas();
+        if (checkTfTangentPlanesFun(tf, d))
+        {
+          Trace("nl-ext-tftp") << "...fail, #lemmas = "
+                               << (d_tstate.d_im.numPendingLemmas()
+                                   + d_tstate.d_im.numWaitingLemmas() - prev)
+                               << std::endl;
+          break;
+        }
+        else
+        {
+          Trace("nl-ext-tftp") << "...success" << std::endl;
+        }
+      }
+    }
+  }
+}
+
+bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d)
+{
+  NodeManager* nm = NodeManager::currentNM();
+  Kind k = tf.getKind();
+  // this should only be run on master applications
+  Assert(d_tstate.d_trSlaves.find(tf) != d_tstate.d_trSlaves.end());
+
+  // Figure 3 : c
+  Node c = d_tstate.d_model.computeAbstractModelValue(tf[0]);
+  int csign = c.getConst<Rational>().sgn();
+  if (csign == 0)
+  {
+    // no secant/tangent plane is necessary
+    return true;
+  }
+  Assert(csign == 1 || csign == -1);
+
+  // Figure 3: P_l, P_u
+  // mapped to for signs of c
+  std::map<int, Node> poly_approx_bounds[2];
+  std::vector<Node> pbounds;
+  d_tstate.d_taylor.getPolynomialApproximationBoundForArg(k, c, d, pbounds);
+  poly_approx_bounds[0][1] = pbounds[0];
+  poly_approx_bounds[0][-1] = pbounds[1];
+  poly_approx_bounds[1][1] = pbounds[2];
+  poly_approx_bounds[1][-1] = pbounds[3];
+
+  // Figure 3 : v
+  Node v = d_tstate.d_model.computeAbstractModelValue(tf);
+
+  // check value of tf
+  Trace("nl-ext-tftp-debug") << "Process tangent plane refinement for " << tf
+                             << ", degree " << d << "..." << std::endl;
+  Trace("nl-ext-tftp-debug") << "  value in model : " << v << std::endl;
+  Trace("nl-ext-tftp-debug") << "  arg value in model : " << c << std::endl;
+
+  // compute the concavity
+  int region = -1;
+  std::unordered_map<Node, int, NodeHashFunction>::iterator itr =
+      d_tstate.d_tf_region.find(tf);
+  if (itr != d_tstate.d_tf_region.end())
+  {
+    region = itr->second;
+    Trace("nl-ext-tftp-debug") << "  region is : " << region << std::endl;
+  }
+  // Figure 3 : conc
+  int concavity = regionToConcavity(k, itr->second);
+  Trace("nl-ext-tftp-debug") << "  concavity is : " << concavity << std::endl;
+  if (concavity == 0)
+  {
+    // no secant/tangent plane is necessary
+    return true;
+  }
+
+  // Figure 3: P
+  Node poly_approx;
+
+  // compute whether this is a tangent refinement or a secant refinement
+  bool is_tangent = false;
+  bool is_secant = false;
+  std::pair<Node, Node> mvb =
+      d_tstate.d_taylor.getTfModelBounds(tf, d, d_tstate.d_model);
+  // this is the approximated value of tf(c), which is a value such that:
+  //    M_A(tf(c)) >= poly_appox_c >= tf(c) or
+  //    M_A(tf(c)) <= poly_appox_c <= tf(c)
+  // In other words, it is a better approximation of the true value of tf(c)
+  // in the case that we add a refinement lemma. We use this value in the
+  // refinement schemas below.
+  Node poly_approx_c;
+  for (unsigned r = 0; r < 2; r++)
+  {
+    Node pab = poly_approx_bounds[r][csign];
+    Node v_pab = r == 0 ? mvb.first : mvb.second;
+    if (!v_pab.isNull())
+    {
+      Trace("nl-ext-tftp-debug2")
+          << "...model value of " << pab << " is " << v_pab << std::endl;
+
+      Assert(v_pab.isConst());
+      Node comp = nm->mkNode(r == 0 ? LT : GT, v, v_pab);
+      Trace("nl-ext-tftp-debug2") << "...compare : " << comp << std::endl;
+      Node compr = Rewriter::rewrite(comp);
+      Trace("nl-ext-tftp-debug2") << "...got : " << compr << std::endl;
+      if (compr == d_tstate.d_true)
+      {
+        poly_approx_c = Rewriter::rewrite(v_pab);
+        // beyond the bounds
+        if (r == 0)
+        {
+          poly_approx = poly_approx_bounds[r][csign];
+          is_tangent = concavity == 1;
+          is_secant = concavity == -1;
+        }
+        else
+        {
+          poly_approx = poly_approx_bounds[r][csign];
+          is_tangent = concavity == -1;
+          is_secant = concavity == 1;
+        }
+        if (Trace.isOn("nl-ext-tftp"))
+        {
+          Trace("nl-ext-tftp") << "*** Outside boundary point (";
+          Trace("nl-ext-tftp") << (r == 0 ? "low" : "high") << ") ";
+          printRationalApprox("nl-ext-tftp", v_pab);
+          Trace("nl-ext-tftp") << ", will refine..." << std::endl;
+          Trace("nl-ext-tftp")
+              << "    poly_approx = " << poly_approx << std::endl;
+          Trace("nl-ext-tftp")
+              << "    is_tangent = " << is_tangent << std::endl;
+          Trace("nl-ext-tftp") << "    is_secant = " << is_secant << std::endl;
+        }
+        break;
+      }
+      else
+      {
+        Trace("nl-ext-tftp")
+            << "  ...within " << (r == 0 ? "low" : "high") << " bound : ";
+        printRationalApprox("nl-ext-tftp", v_pab);
+        Trace("nl-ext-tftp") << std::endl;
+      }
+    }
+  }
+
+  // Figure 3: P( c )
+  if (is_tangent || is_secant)
+  {
+    Trace("nl-ext-tftp-debug2")
+        << "...poly approximation at c is " << poly_approx_c << std::endl;
+  }
+  else
+  {
+    // we may want to continue getting better bounds
+    return false;
+  }
+
+  if (is_tangent)
+  {
+    if (k == Kind::EXPONENTIAL)
+    {
+      d_expSlv.doTangentLemma(tf, c, poly_approx_c);
+    }
+    else if (k == Kind::SINE)
+    {
+      d_sineSlv.doTangentLemma(tf, c, poly_approx_c, region);
+    }
+  }
+  else if (is_secant)
+  {
+    if (k == EXPONENTIAL)
+    {
+      d_expSlv.doSecantLemmas(tf, c, poly_approx_c, d);
+    }
+    else if (k == Kind::SINE)
+    {
+      d_sineSlv.doSecantLemmas(tf, c, poly_approx_c, d, region);
+    }
+  }
+  return true;
+}
+
+int TranscendentalSolver::regionToConcavity(Kind k, int region)
+{
+  if (k == EXPONENTIAL)
+  {
+    if (region == 1)
+    {
+      return 1;
+    }
+  }
+  else if (k == SINE)
+  {
+    if (region == 1 || region == 2)
+    {
+      return -1;
+    }
+    else if (region == 3 || region == 4)
+    {
+      return 1;
+    }
+  }
+  return 0;
+}
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl/transcendental/transcendental_solver.h b/src/theory/arith/nl/transcendental/transcendental_solver.h
new file mode 100644
index 000000000..b0db79b6a
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/transcendental_solver.h
@@ -0,0 +1,204 @@
+
+/*********************                                                        */
+/*! \file transcendental_solver.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** 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 Solving for handling transcendental functions.
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TRANSCENDENTAL_SOLVER_H
+#define CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TRANSCENDENTAL_SOLVER_H
+
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+#include "expr/node.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+#include "theory/arith/nl/transcendental/exponential_solver.h"
+#include "theory/arith/nl/transcendental/sine_solver.h"
+#include "theory/arith/nl/transcendental/transcendental_state.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+/** Transcendental solver class
+ *
+ * This class implements model-based refinement schemes
+ * for transcendental functions, described in:
+ *
+ * - "Satisfiability Modulo Transcendental
+ * Functions via Incremental Linearization" by Cimatti
+ * et al., CADE 2017.
+ *
+ * It's main functionality are methods that implement lemma schemas below,
+ * which return a set of lemmas that should be sent on the output channel.
+ */
+class TranscendentalSolver
+{
+ public:
+  TranscendentalSolver(InferenceManager& im, NlModel& m);
+  ~TranscendentalSolver();
+
+  /** init last call
+   *
+   * This is called at the beginning of last call effort check, where
+   * assertions are the set of assertions belonging to arithmetic,
+   * false_asserts is the subset of assertions that are false in the current
+   * model, and xts is the set of extended function terms that are active in
+   * the current context.
+   *
+   * This call may add lemmas to lems based on registering term
+   * information (for example, purification of sine terms).
+   */
+  void initLastCall(const std::vector<Node>& assertions,
+                    const std::vector<Node>& false_asserts,
+                    const std::vector<Node>& xts);
+  /** increment taylor degree */
+  void incrementTaylorDegree();
+  /** get taylor degree */
+  unsigned getTaylorDegree() const;
+  /** preprocess assertions check model
+   *
+   * This modifies the given assertions in preparation for running a call
+   * to check model.
+   *
+   * This method returns false if a bound for a transcendental function
+   * was conflicting.
+   */
+  bool preprocessAssertionsCheckModel(std::vector<Node>& assertions);
+  /** Process side effects in lemma se */
+  void processSideEffect(const NlLemma& se);
+  //-------------------------------------------- lemma schemas
+  // Relays to exp and sine solvers.
+  void checkTranscendentalInitialRefine();
+
+  // Relays to exp and sine solvers.
+  void checkTranscendentalMonotonic();
+
+  /** check transcendental tangent planes
+   *
+   * Constructs a set of valid theory lemmas, based on
+   * computing an "incremental linearization" of
+   * transcendental functions based on the model values
+   * of transcendental functions and their arguments.
+   * It is based on Figure 3 of "Satisfiability
+   * Modulo Transcendental Functions via Incremental
+   * Linearization" by Cimatti et al., CADE 2017.
+   * This schema is not terminating in general.
+   * It is not enabled by default, and can
+   * be enabled by --nl-ext-tf-tplanes.
+   *
+   * Example:
+   *
+   * Assume we have a term sin(y) where M( y ) = 1 where M is the current model.
+   * Note that:
+   *   sin(1) ~= .841471
+   *
+   * The Taylor series and remainder of sin(y) of degree 7 is
+   *   P_{7,sin(0)}( x ) = x + (-1/6)*x^3 + (1/20)*x^5
+   *   R_{7,sin(0),b}( x ) = (-1/5040)*x^7
+   *
+   * This gives us lower and upper bounds :
+   *   P_u( x ) = P_{7,sin(0)}( x ) + R_{7,sin(0),b}( x )
+   *     ...where note P_u( 1 ) = 4243/5040 ~= .841865
+   *   P_l( x ) = P_{7,sin(0)}( x ) - R_{7,sin(0),b}( x )
+   *     ...where note P_l( 1 ) = 4241/5040 ~= .841468
+   *
+   * Assume that M( sin(y) ) > P_u( 1 ).
+   * Since the concavity of sine in the region 0 < x < PI/2 is -1,
+   * we add a tangent plane refinement.
+   * The tangent plane at the point 1 in P_u is
+   * given by the formula:
+   *   T( x ) = P_u( 1 ) + ((d/dx)(P_u(x)))( 1 )*( x - 1 )
+   * We add the lemma:
+   *   ( 0 < y < PI/2 ) => sin( y ) <= T( y )
+   * which is:
+   *   ( 0 < y < PI/2 ) => sin( y ) <= (391/720)*(y - 2737/1506)
+   *
+   * Assume that M( sin(y) ) < P_u( 1 ).
+   * Since the concavity of sine in the region 0 < x < PI/2 is -1,
+   * we add a secant plane refinement for some constants ( l, u )
+   * such that 0 <= l < M( y ) < u <= PI/2. Assume we choose
+   * l = 0 and u = M( PI/2 ) = 150517/47912.
+   * The secant planes at point 1 for P_l
+   * are given by the formulas:
+   *   S_l( x ) = (x-l)*(P_l( l )-P_l(c))/(l-1) + P_l( l )
+   *   S_u( x ) = (x-u)*(P_l( u )-P_l(c))/(u-1) + P_l( u )
+   * We add the lemmas:
+   *   ( 0 < y < 1 ) => sin( y ) >= S_l( y )
+   *   ( 1 < y < PI/2 ) => sin( y ) >= S_u( y )
+   * which are:
+   *   ( 0 < y < 1 ) => (sin y) >= 4251/5040*y
+   *   ( 1 < y < PI/2 ) => (sin y) >= c1*(y+c2)
+   *     where c1, c2 are rationals (for brevity, omitted here)
+   *     such that c1 ~= .277 and c2 ~= 2.032.
+   */
+  void checkTranscendentalTangentPlanes();
+
+ private:
+  /** check transcendental function refinement for tf
+   *
+   * This method is called by the above method for each "master"
+   * transcendental function application that occurs in an assertion in the
+   * current context. For example, an application like sin(t) is not a master
+   * if we have introduced the constraints:
+   *   t=y+2*pi*n ^ -pi <= y <= pi ^ sin(t) = sin(y).
+   * See d_trMaster/d_trSlaves for more detail.
+   *
+   * This runs Figure 3 of Cimatti et al., CADE 2017 for transcendental
+   * function application tf for Taylor degree d. It may add a secant or
+   * tangent plane lemma to lems, which includes the side effect of the lemma
+   * (if one exists).
+   *
+   * It returns false if the bounds are not precise enough to add a
+   * secant or tangent plane lemma.
+   */
+  bool checkTfTangentPlanesFun(Node tf, unsigned d);
+  //-------------------------------------------- end lemma schemas
+
+  /** get concavity
+   *
+   * Returns whether we are concave (+1) or convex (-1)
+   * in region of transcendental function with kind k,
+   * where region is defined above.
+   * Returns 0 if region is invalid.
+   */
+  int regionToConcavity(Kind k, int region);
+
+  /** taylor degree
+   *
+   * Indicates that the degree of the polynomials in the Taylor approximation of
+   * all transcendental functions is 2*d_taylor_degree. This value is set
+   * initially to options::nlExtTfTaylorDegree() and may be incremented
+   * if the option options::nlExtTfIncPrecision() is enabled.
+   */
+  unsigned d_taylor_degree;
+
+  /** Common state for transcendental solver */
+  transcendental::TranscendentalState d_tstate;
+  /** The solver responsible for the exponential function */
+  transcendental::ExponentialSolver d_expSlv;
+  /** The solver responsible for the sine function */
+  transcendental::SineSolver d_sineSlv;
+}; /* class TranscendentalSolver */
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__TRANSCENDENTAL_SOLVER_H */
diff --git a/src/theory/arith/nl/transcendental/transcendental_state.cpp b/src/theory/arith/nl/transcendental/transcendental_state.cpp
new file mode 100644
index 000000000..b6882d606
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/transcendental_state.cpp
@@ -0,0 +1,384 @@
+/*********************                                                        */
+/*! \file transcendental_state.cpp
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King, 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 solver for handling transcendental functions.
+ **/
+
+#include "theory/arith/nl/transcendental/transcendental_state.h"
+
+#include "theory/arith/arith_utilities.h"
+#include "theory/arith/nl/transcendental/taylor_generator.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+namespace {
+
+/**
+ * Ensure a is in the main phase:
+ *   -pi <= a <= pi
+ */
+inline Node mkValidPhase(TNode a, TNode pi)
+{
+  return mkBounded(
+      NodeManager::currentNM()->mkNode(Kind::MULT, mkRationalNode(-1), pi),
+      a,
+      pi);
+}
+}  // namespace
+
+TranscendentalState::TranscendentalState(InferenceManager& im, NlModel& model)
+    : d_im(im), d_model(model)
+{
+  d_true = NodeManager::currentNM()->mkConst(true);
+  d_false = NodeManager::currentNM()->mkConst(false);
+  d_zero = NodeManager::currentNM()->mkConst(Rational(0));
+  d_one = NodeManager::currentNM()->mkConst(Rational(1));
+  d_neg_one = NodeManager::currentNM()->mkConst(Rational(-1));
+}
+
+void TranscendentalState::init(const std::vector<Node>& assertions,
+                               const std::vector<Node>& false_asserts,
+                               const std::vector<Node>& xts)
+{
+  d_funcCongClass.clear();
+  d_funcMap.clear();
+  d_tf_region.clear();
+
+  NodeManager* nm = NodeManager::currentNM();
+
+  // register the extended function terms
+  std::vector<Node> trNeedsMaster;
+  bool needPi = false;
+  // for computing congruence
+  std::map<Kind, ArgTrie> argTrie;
+  for (unsigned i = 0, xsize = xts.size(); i < xsize; i++)
+  {
+    Node a = xts[i];
+    Kind ak = a.getKind();
+    bool consider = true;
+    // if is an unpurified application of SINE, or it is a transcendental
+    // applied to a trancendental, purify.
+    if (isTranscendentalKind(ak))
+    {
+      // if we've already computed master for a
+      if (d_trMaster.find(a) != d_trMaster.end())
+      {
+        // a master has at least one slave
+        consider = (d_trSlaves.find(a) != d_trSlaves.end());
+      }
+      else
+      {
+        if (ak == Kind::SINE)
+        {
+          // always not a master
+          consider = false;
+        }
+        else
+        {
+          for (const Node& ac : a)
+          {
+            if (isTranscendentalKind(ac.getKind()))
+            {
+              consider = false;
+              break;
+            }
+          }
+        }
+        if (!consider)
+        {
+          // wait to assign a master below
+          trNeedsMaster.push_back(a);
+        }
+        else
+        {
+          d_trMaster[a] = a;
+          d_trSlaves[a].insert(a);
+        }
+      }
+    }
+    if (ak == Kind::EXPONENTIAL || ak == Kind::SINE)
+    {
+      needPi = needPi || (ak == Kind::SINE);
+      // if we didn't indicate that it should be purified above
+      if (consider)
+      {
+        std::vector<Node> repList;
+        for (const Node& ac : a)
+        {
+          Node r = d_model.computeConcreteModelValue(ac);
+          repList.push_back(r);
+        }
+        Node aa = argTrie[ak].add(a, repList);
+        if (aa != a)
+        {
+          // apply congruence to pairs of terms that are disequal and congruent
+          Assert(aa.getNumChildren() == a.getNumChildren());
+          Node mvaa = d_model.computeAbstractModelValue(a);
+          Node mvaaa = d_model.computeAbstractModelValue(aa);
+          if (mvaa != mvaaa)
+          {
+            std::vector<Node> exp;
+            for (unsigned j = 0, size = a.getNumChildren(); j < size; j++)
+            {
+              exp.push_back(a[j].eqNode(aa[j]));
+            }
+            Node expn = exp.size() == 1 ? exp[0] : nm->mkNode(Kind::AND, exp);
+            Node cong_lemma = nm->mkNode(Kind::OR, expn.negate(), a.eqNode(aa));
+            d_im.addPendingArithLemma(cong_lemma, InferenceId::NL_CONGRUENCE);
+          }
+        }
+        else
+        {
+          // new representative of congruence class
+          d_funcMap[ak].push_back(a);
+        }
+        // add to congruence class
+        d_funcCongClass[aa].push_back(a);
+      }
+    }
+    else if (ak == Kind::PI)
+    {
+      Assert(consider);
+      needPi = true;
+      d_funcMap[ak].push_back(a);
+      d_funcCongClass[a].push_back(a);
+    }
+  }
+  // initialize pi if necessary
+  if (needPi && d_pi.isNull())
+  {
+    mkPi();
+    getCurrentPiBounds();
+  }
+
+  if (d_im.hasUsed())
+  {
+    return;
+  }
+
+  // process SINE phase shifting
+  for (const Node& a : trNeedsMaster)
+  {
+    // should not have processed this already
+    Assert(d_trMaster.find(a) == d_trMaster.end());
+    Kind k = a.getKind();
+    Assert(k == Kind::SINE || k == Kind::EXPONENTIAL);
+    Node y =
+        nm->mkSkolem("y", nm->realType(), "phase shifted trigonometric arg");
+    Node new_a = nm->mkNode(k, y);
+    d_trSlaves[new_a].insert(new_a);
+    d_trSlaves[new_a].insert(a);
+    d_trMaster[a] = new_a;
+    d_trMaster[new_a] = new_a;
+    Node lem;
+    if (k == Kind::SINE)
+    {
+      Trace("nl-ext-tf") << "Basis sine : " << new_a << " for " << a
+                         << std::endl;
+      Assert(!d_pi.isNull());
+      Node shift = nm->mkSkolem("s", nm->integerType(), "number of shifts");
+      // TODO : do not introduce shift here, instead needs model-based
+      // refinement for constant shifts (cvc4-projects #1284)
+      lem = nm->mkNode(
+          Kind::AND,
+          transcendental::mkValidPhase(y, d_pi),
+          nm->mkNode(
+              Kind::ITE,
+              transcendental::mkValidPhase(a[0], d_pi),
+              a[0].eqNode(y),
+              a[0].eqNode(nm->mkNode(
+                  Kind::PLUS,
+                  y,
+                  nm->mkNode(
+                      Kind::MULT, nm->mkConst(Rational(2)), shift, d_pi)))),
+          new_a.eqNode(a));
+    }
+    else
+    {
+      // do both equalities to ensure that new_a becomes a preregistered term
+      lem = nm->mkNode(Kind::AND, a.eqNode(new_a), a[0].eqNode(y));
+    }
+    // note we must do preprocess on this lemma
+    Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : purify : " << lem
+                          << std::endl;
+    NlLemma nlem(
+        lem, LemmaProperty::PREPROCESS, nullptr, InferenceId::NL_T_PURIFY_ARG);
+    d_im.addPendingArithLemma(nlem);
+  }
+
+  if (Trace.isOn("nl-ext-mv"))
+  {
+    Trace("nl-ext-mv") << "Arguments of trancendental functions : "
+                       << std::endl;
+    for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
+    {
+      Kind k = tfl.first;
+      if (k == Kind::SINE || k == Kind::EXPONENTIAL)
+      {
+        for (const Node& tf : tfl.second)
+        {
+          Node v = tf[0];
+          d_model.computeConcreteModelValue(v);
+          d_model.computeAbstractModelValue(v);
+          d_model.printModelValue("nl-ext-mv", v);
+        }
+      }
+    }
+  }
+}
+
+void TranscendentalState::mkPi()
+{
+  NodeManager* nm = NodeManager::currentNM();
+  if (d_pi.isNull())
+  {
+    d_pi = nm->mkNullaryOperator(nm->realType(), Kind::PI);
+    d_pi_2 = Rewriter::rewrite(
+        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(1) / Rational(2))));
+    d_pi_neg_2 = Rewriter::rewrite(
+        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1) / Rational(2))));
+    d_pi_neg = Rewriter::rewrite(
+        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1))));
+    // initialize bounds
+    d_pi_bound[0] = nm->mkConst(Rational(103993) / Rational(33102));
+    d_pi_bound[1] = nm->mkConst(Rational(104348) / Rational(33215));
+  }
+}
+
+void TranscendentalState::getCurrentPiBounds()
+{
+  NodeManager* nm = NodeManager::currentNM();
+  Node pi_lem = nm->mkNode(Kind::AND,
+                           nm->mkNode(Kind::GEQ, d_pi, d_pi_bound[0]),
+                           nm->mkNode(Kind::LEQ, d_pi, d_pi_bound[1]));
+  d_im.addPendingArithLemma(pi_lem, InferenceId::NL_T_PI_BOUND);
+}
+
+std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e,
+                                                                  TNode c,
+                                                                  unsigned d)
+{
+  // bounds are the minimum and maximum previous secant points
+  // should not repeat secant points: secant lemmas should suffice to
+  // rule out previous assignment
+  Assert(
+      std::find(d_secant_points[e][d].begin(), d_secant_points[e][d].end(), c)
+      == d_secant_points[e][d].end());
+  // Insert into the (temporary) vector. We do not update this vector
+  // until we are sure this secant plane lemma has been processed. We do
+  // this by mapping the lemma to a side effect below.
+  std::vector<Node> spoints = d_secant_points[e][d];
+  spoints.push_back(c);
+
+  // sort
+  sortByNlModel(spoints.begin(), spoints.end(), &d_model);
+  // get the resulting index of c
+  unsigned index =
+      std::find(spoints.begin(), spoints.end(), c) - spoints.begin();
+
+  // bounds are the next closest upper/lower bound values
+  return {index > 0 ? spoints[index - 1] : Node(),
+          index < spoints.size() - 1 ? spoints[index + 1] : Node()};
+}
+
+Node TranscendentalState::mkSecantPlane(
+    TNode arg, TNode b, TNode c, TNode approx, TNode approx_c)
+{
+  NodeManager* nm = NodeManager::currentNM();
+  // Figure 3 : P(l), P(u), for s = 0,1
+  Node approx_b =
+      Rewriter::rewrite(approx.substitute(d_taylor.getTaylorVariable(), b));
+  // Figure 3: S_l( x ), S_u( x ) for s = 0,1
+  Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, b, c));
+  Assert(rcoeff_n.isConst());
+  Rational rcoeff = rcoeff_n.getConst<Rational>();
+  Assert(rcoeff.sgn() != 0);
+  return nm->mkNode(Kind::PLUS,
+                    approx_b,
+                    nm->mkNode(Kind::MULT,
+                               nm->mkNode(Kind::MINUS, approx_b, approx_c),
+                               nm->mkConst(rcoeff.inverse()),
+                               nm->mkNode(Kind::MINUS, arg, b)));
+}
+
+NlLemma TranscendentalState::mkSecantLemma(
+    TNode lower, TNode upper, int concavity, TNode tf, TNode splane)
+{
+  NodeManager* nm = NodeManager::currentNM();
+  // With respect to Figure 3, this is slightly different.
+  // In particular, we chose b to be the model value of bounds[s],
+  // which is a constant although bounds[s] may not be (e.g. if it
+  // contains PI).
+  // To ensure that c...b does not cross an inflection point,
+  // we guard with the symbolic version of bounds[s].
+  // This leads to lemmas e.g. of this form:
+  //   ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x -
+  //   b ) + P( b ) )
+  // where b = (PI/2)^M, the current value of PI/2 in the model.
+  // This is sound since we are guarded by the symbolic
+  // representation of PI/2.
+  Node antec_n = nm->mkNode(Kind::AND,
+                            nm->mkNode(Kind::GEQ, tf[0], lower),
+                            nm->mkNode(Kind::LEQ, tf[0], upper));
+  Node lem = nm->mkNode(
+      Kind::IMPLIES,
+      antec_n,
+      nm->mkNode(concavity == 1 ? Kind::LEQ : Kind::GEQ, tf, splane));
+  Trace("nl-ext-tftp-debug2")
+      << "*** Secant plane lemma (pre-rewrite) : " << lem << std::endl;
+  lem = Rewriter::rewrite(lem);
+  Trace("nl-ext-tftp-lemma") << "*** Secant plane lemma : " << lem << std::endl;
+  Assert(d_model.computeAbstractModelValue(lem) == d_false);
+  return NlLemma(lem, LemmaProperty::NONE, nullptr, InferenceId::NL_T_SECANT);
+}
+
+void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds,
+                                         TNode c,
+                                         TNode approx_c,
+                                         TNode tf,
+                                         unsigned d,
+                                         int concavity)
+{
+  // take the model value of l or u (since may contain PI)
+  // Make secant from bounds.first to c
+  Node lval = d_model.computeAbstractModelValue(bounds.first);
+  if (lval != c)
+  {
+    Node splane = mkSecantPlane(tf[0], lval, c, bounds.first, approx_c);
+    NlLemma nlem = mkSecantLemma(lval, c, concavity, tf, splane);
+    // The side effect says that if lem is added, then we should add the
+    // secant point c for (tf,d).
+    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, c));
+    d_im.addPendingArithLemma(nlem, true);
+  }
+
+  // Make secant from c to bounds.second
+  Node uval = d_model.computeAbstractModelValue(bounds.second);
+  if (c != uval)
+  {
+    Node splane = mkSecantPlane(tf[0], c, uval, approx_c, bounds.second);
+    NlLemma nlem = mkSecantLemma(c, uval, concavity, tf, splane);
+    // The side effect says that if lem is added, then we should add the
+    // secant point c for (tf,d).
+    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, c));
+    d_im.addPendingArithLemma(nlem, true);
+  }
+}
+
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
diff --git a/src/theory/arith/nl/transcendental/transcendental_state.h b/src/theory/arith/nl/transcendental/transcendental_state.h
new file mode 100644
index 000000000..2cf6400a3
--- /dev/null
+++ b/src/theory/arith/nl/transcendental/transcendental_state.h
@@ -0,0 +1,213 @@
+/*********************                                                        */
+/*! \file transcendental_state.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds, Tim King
+ ** 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 Utilities for transcendental lemmas.
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TRANSCENDENTAL_STATE_H
+#define CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TRANSCENDENTAL_STATE_H
+
+#include "expr/node.h"
+#include "theory/arith/inference_manager.h"
+#include "theory/arith/nl/nl_model.h"
+#include "theory/arith/nl/transcendental/taylor_generator.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+namespace nl {
+namespace transcendental {
+
+/**
+ * Holds common state and utilities for transcendental solvers.
+ *
+ * This includes common lookups and caches as well as generic utilities for
+ * secant plane lemmas and taylor approximations.
+ */
+struct TranscendentalState
+{
+  TranscendentalState(InferenceManager& im, NlModel& model);
+
+  /** init last call
+   *
+   * This is called at the beginning of last call effort check, where
+   * assertions are the set of assertions belonging to arithmetic,
+   * false_asserts is the subset of assertions that are false in the current
+   * model, and xts is the set of extended function terms that are active in
+   * the current context.
+   *
+   * This call may add lemmas to lems based on registering term
+   * information (for example, purification of sine terms).
+   */
+  void init(const std::vector<Node>& assertions,
+            const std::vector<Node>& false_asserts,
+            const std::vector<Node>& xts);
+
+  /** Initialize members for pi-related values */
+  void mkPi();
+  /** Get current bounds for pi as a lemma */
+  void getCurrentPiBounds();
+
+  /**
+   * Get the two closest secant points from the once stored already.
+   * "closest" is determined according to the current model.
+   * @param e The transcendental term (like (exp t))
+   * @param c The point currently under consideration (probably the model of t)
+   * @param d The taylor degree.
+   */
+  std::pair<Node, Node> getClosestSecantPoints(TNode e, TNode c, unsigned d);
+
+  /**
+   * Construct a secant plane between b and c
+   * @param arg The argument of the transcendental term
+   * @param b Left secant point
+   * @param c Right secant point
+   * @param approx Approximation for b (not yet substituted)
+   * @param approx_c Approximation for c (already substituted)
+   */
+  Node mkSecantPlane(TNode arg, TNode b, TNode c, TNode approx, TNode approx_c);
+
+  /**
+   * Construct a secant lemma between lower and upper for tf.
+   * @param lower Left secant point
+   * @param upper Right secant point
+   * @param concavity Concavity of the current regions
+   * @param tf Current transcendental term
+   * @param splane Secant plane as computed by mkSecantPlane()
+   */
+  NlLemma mkSecantLemma(
+      TNode lower, TNode upper, int concavity, TNode tf, TNode splane);
+
+  /**
+   * Construct and send secant lemmas (if appropriate)
+   * @param bounds Secant bounds
+   * @param c Current point
+   * @param approx_c Approximation for c
+   * @param tf Current transcendental term
+   * @param d Current taylor degree
+   * @param concavity Concavity in region of c
+   */
+  void doSecantLemmas(const std::pair<Node, Node>& bounds,
+                      TNode c,
+                      TNode approx_c,
+                      TNode tf,
+                      unsigned d,
+                      int concavity)
+    ;
+
+    Node d_true;
+    Node d_false;
+    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;
+    /** Utility to compute taylor approximations */
+    TaylorGenerator d_taylor;
+
+    /**
+     * Some transcendental functions f(t) are "purified", e.g. we add
+     * t = y ^ f(t) = f(y) where y is a fresh variable. Those that are not
+     * purified we call "master terms".
+     *
+     * The maps below maintain a master/slave relationship over
+     * transcendental functions (SINE, EXPONENTIAL, PI), where above
+     * f(y) is the master of itself and of f(t).
+     *
+     * This is used for ensuring that the argument y of SINE we process is on
+     * the interval [-pi .. pi], and that exponentials are not applied to
+     * arguments that contain transcendental functions.
+     */
+    std::map<Node, Node> d_trMaster;
+    std::map<Node, std::unordered_set<Node, NodeHashFunction>> d_trSlaves;
+
+    /** concavity region for transcendental functions
+     *
+     * This stores an integer that identifies an interval in
+     * which the current model value for an argument of an
+     * application of a transcendental function resides.
+     *
+     * For exp( x ):
+     *   region #1 is -infty < x < infty
+     * For sin( x ):
+     *   region #0 is pi < x < infty (this is an invalid region)
+     *   region #1 is pi/2 < x <= pi
+     *   region #2 is 0 < x <= pi/2
+     *   region #3 is -pi/2 < x <= 0
+     *   region #4 is -pi < x <= -pi/2
+     *   region #5 is -infty < x <= -pi (this is an invalid region)
+     * All regions not listed above, as well as regions 0 and 5
+     * for SINE are "invalid". We only process applications
+     * of transcendental functions whose arguments have model
+     * values that reside in valid regions.
+     */
+    std::unordered_map<Node, int, NodeHashFunction> d_tf_region;
+    /**
+     * Maps representives of a congruence class to the members of that class.
+     *
+     * In detail, a congruence class is a set of terms of the form
+     *   { f(t1), ..., f(tn) }
+     * such that t1 = ... = tn in the current context. We choose an arbitrary
+     * term among these to be the repesentative of this congruence class.
+     *
+     * Moreover, notice we compute congruence classes only over terms that
+     * are transcendental function applications that are "master terms",
+     * see d_trMaster/d_trSlave.
+     */
+    std::map<Node, std::vector<Node>> d_funcCongClass;
+    /**
+     * A list of all functions for each kind in { EXPONENTIAL, SINE, POW, PI }
+     * that are representives of their congruence class.
+     */
+    std::map<Kind, std::vector<Node>> d_funcMap;
+
+    /** secant points (sorted list) for transcendental functions
+     *
+     * This is used for tangent plane refinements for
+     * transcendental functions. This is the set
+     * "get-previous-secant-points" in "Satisfiability
+     * Modulo Transcendental Functions via Incremental
+     * Linearization" by Cimatti et al., CADE 2017, for
+     * each transcendental function application. We store this set for each
+     * Taylor degree.
+     */
+    std::unordered_map<Node,
+                       std::map<unsigned, std::vector<Node>>,
+                       NodeHashFunction>
+        d_secant_points;
+
+    /** PI
+     *
+     * Note that PI is a (symbolic, non-constant) nullary operator. This is
+     * because its value cannot be computed exactly. We constraint PI to
+     * concrete lower and upper bounds stored in d_pi_bound below.
+     */
+    Node d_pi;
+    /** PI/2 */
+    Node d_pi_2;
+    /** -PI/2 */
+    Node d_pi_neg_2;
+    /** -PI */
+    Node d_pi_neg;
+    /** the concrete lower and upper bounds for PI */
+    Node d_pi_bound[2];
+  };
+
+}  // namespace transcendental
+}  // namespace transcendental
+}  // namespace nl
+}  // namespace arith
+}  // namespace theory
+
+#endif /* CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TRANSCENDENTAL_STATE_H */
diff --git a/src/theory/arith/nl/transcendental_solver.cpp b/src/theory/arith/nl/transcendental_solver.cpp
deleted file mode 100644
index 2e25f1642..000000000
--- a/src/theory/arith/nl/transcendental_solver.cpp
+++ /dev/null
@@ -1,1474 +0,0 @@
-/*********************                                                        */
-/*! \file transcendental_solver.cpp
- ** \verbatim
- ** Top contributors (to current version):
- **   Andrew Reynolds, Tim King, 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 solver for handling transcendental functions.
- **/
-
-#include "theory/arith/nl/transcendental_solver.h"
-
-#include <cmath>
-#include <set>
-
-#include "expr/node_algorithm.h"
-#include "expr/node_builder.h"
-#include "options/arith_options.h"
-#include "theory/arith/arith_msum.h"
-#include "theory/arith/arith_utilities.h"
-#include "theory/rewriter.h"
-
-using namespace CVC4::kind;
-
-namespace CVC4 {
-namespace theory {
-namespace arith {
-namespace nl {
-
-TranscendentalSolver::TranscendentalSolver(InferenceManager& im, NlModel& m) : d_im(im), d_model(m)
-{
-  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_taylor_real_fv = nm->mkBoundVar("x", nm->realType());
-  d_taylor_real_fv_base = nm->mkBoundVar("a", nm->realType());
-  d_taylor_real_fv_base_rem = nm->mkBoundVar("b", nm->realType());
-  d_taylor_degree = options::nlExtTfTaylorDegree();
-}
-
-TranscendentalSolver::~TranscendentalSolver() {}
-
-void TranscendentalSolver::initLastCall(const std::vector<Node>& assertions,
-                                        const std::vector<Node>& false_asserts,
-                                        const std::vector<Node>& xts)
-{
-  d_funcCongClass.clear();
-  d_funcMap.clear();
-  d_tf_region.clear();
-
-  NodeManager* nm = NodeManager::currentNM();
-
-  // register the extended function terms
-  std::vector<Node> trNeedsMaster;
-  bool needPi = false;
-  // for computing congruence
-  std::map<Kind, ArgTrie> argTrie;
-  for (unsigned i = 0, xsize = xts.size(); i < xsize; i++)
-  {
-    Node a = xts[i];
-    Kind ak = a.getKind();
-    bool consider = true;
-    // if is an unpurified application of SINE, or it is a transcendental
-    // applied to a trancendental, purify.
-    if (isTranscendentalKind(ak))
-    {
-      // if we've already computed master for a
-      if (d_trMaster.find(a) != d_trMaster.end())
-      {
-        // a master has at least one slave
-        consider = (d_trSlaves.find(a) != d_trSlaves.end());
-      }
-      else
-      {
-        if (ak == SINE)
-        {
-          // always not a master
-          consider = false;
-        }
-        else
-        {
-          for (const Node& ac : a)
-          {
-            if (isTranscendentalKind(ac.getKind()))
-            {
-              consider = false;
-              break;
-            }
-          }
-        }
-        if (!consider)
-        {
-          // wait to assign a master below
-          trNeedsMaster.push_back(a);
-        }
-        else
-        {
-          d_trMaster[a] = a;
-          d_trSlaves[a].insert(a);
-        }
-      }
-    }
-    if (ak == EXPONENTIAL || ak == SINE)
-    {
-      needPi = needPi || (ak == SINE);
-      // if we didn't indicate that it should be purified above
-      if (consider)
-      {
-        std::vector<Node> repList;
-        for (const Node& ac : a)
-        {
-          Node r = d_model.computeConcreteModelValue(ac);
-          repList.push_back(r);
-        }
-        Node aa = argTrie[ak].add(a, repList);
-        if (aa != a)
-        {
-          // apply congruence to pairs of terms that are disequal and congruent
-          Assert(aa.getNumChildren() == a.getNumChildren());
-          Node mvaa = d_model.computeAbstractModelValue(a);
-          Node mvaaa = d_model.computeAbstractModelValue(aa);
-          if (mvaa != mvaaa)
-          {
-            std::vector<Node> exp;
-            for (unsigned j = 0, size = a.getNumChildren(); j < size; j++)
-            {
-              exp.push_back(a[j].eqNode(aa[j]));
-            }
-            Node expn = exp.size() == 1 ? exp[0] : nm->mkNode(AND, exp);
-            Node cong_lemma = nm->mkNode(OR, expn.negate(), a.eqNode(aa));
-            d_im.addPendingArithLemma(cong_lemma, InferenceId::NL_CONGRUENCE);
-          }
-        }
-        else
-        {
-          // new representative of congruence class
-          d_funcMap[ak].push_back(a);
-        }
-        // add to congruence class
-        d_funcCongClass[aa].push_back(a);
-      }
-    }
-    else if (ak == PI)
-    {
-      Assert(consider);
-      needPi = true;
-      d_funcMap[ak].push_back(a);
-      d_funcCongClass[a].push_back(a);
-    }
-  }
-  // initialize pi if necessary
-  if (needPi && d_pi.isNull())
-  {
-    mkPi();
-    getCurrentPiBounds();
-  }
-
-  if (d_im.hasUsed())
-  {
-    return;
-  }
-
-  // process SINE phase shifting
-  for (const Node& a : trNeedsMaster)
-  {
-    // should not have processed this already
-    Assert(d_trMaster.find(a) == d_trMaster.end());
-    Kind k = a.getKind();
-    Assert(k == SINE || k == EXPONENTIAL);
-    Node y =
-        nm->mkSkolem("y", nm->realType(), "phase shifted trigonometric arg");
-    Node new_a = nm->mkNode(k, y);
-    d_trSlaves[new_a].insert(new_a);
-    d_trSlaves[new_a].insert(a);
-    d_trMaster[a] = new_a;
-    d_trMaster[new_a] = new_a;
-    Node lem;
-    if (k == SINE)
-    {
-      Trace("nl-ext-tf") << "Basis sine : " << new_a << " for " << a
-                         << std::endl;
-      Assert(!d_pi.isNull());
-      Node shift = nm->mkSkolem("s", nm->integerType(), "number of shifts");
-      // TODO : do not introduce shift here, instead needs model-based
-      // refinement for constant shifts (cvc4-projects #1284)
-      lem = nm->mkNode(
-          AND,
-          mkValidPhase(y, d_pi),
-          nm->mkNode(
-              ITE,
-              mkValidPhase(a[0], d_pi),
-              a[0].eqNode(y),
-              a[0].eqNode(nm->mkNode(
-                  PLUS,
-                  y,
-                  nm->mkNode(MULT, nm->mkConst(Rational(2)), shift, d_pi)))),
-          new_a.eqNode(a));
-    }
-    else
-    {
-      // do both equalities to ensure that new_a becomes a preregistered term
-      lem = nm->mkNode(AND, a.eqNode(new_a), a[0].eqNode(y));
-    }
-    // note we must do preprocess on this lemma
-    Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : purify : " << lem
-                          << std::endl;
-    NlLemma nlem(lem, LemmaProperty::PREPROCESS, nullptr, InferenceId::NL_T_PURIFY_ARG);
-    d_im.addPendingArithLemma(nlem);
-  }
-
-  if (Trace.isOn("nl-ext-mv"))
-  {
-    Trace("nl-ext-mv") << "Arguments of trancendental functions : "
-                       << std::endl;
-    for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
-    {
-      Kind k = tfl.first;
-      if (k == SINE || k == EXPONENTIAL)
-      {
-        for (const Node& tf : tfl.second)
-        {
-          Node v = tf[0];
-          d_model.computeConcreteModelValue(v);
-          d_model.computeAbstractModelValue(v);
-          d_model.printModelValue("nl-ext-mv", v);
-        }
-      }
-    }
-  }
-}
-
-bool TranscendentalSolver::preprocessAssertionsCheckModel(
-    std::vector<Node>& assertions)
-{
-  std::vector<Node> pvars;
-  std::vector<Node> psubs;
-  for (const std::pair<const Node, Node>& tb : d_trMaster)
-  {
-    pvars.push_back(tb.first);
-    psubs.push_back(tb.second);
-  }
-
-  // initialize representation of assertions
-  std::vector<Node> passertions;
-  for (const Node& a : assertions)
-
-  {
-    Node pa = a;
-    if (!pvars.empty())
-    {
-      pa = arithSubstitute(pa, pvars, psubs);
-      pa = Rewriter::rewrite(pa);
-    }
-    if (!pa.isConst() || !pa.getConst<bool>())
-    {
-      Trace("nl-ext-cm-assert") << "- assert : " << pa << std::endl;
-      passertions.push_back(pa);
-    }
-  }
-  // get model bounds for all transcendental functions
-  Trace("nl-ext-cm") << "----- Get bounds for transcendental functions..."
-                     << std::endl;
-  for (std::pair<const Kind, std::vector<Node> >& tfs : d_funcMap)
-  {
-    Kind k = tfs.first;
-    for (const Node& tf : tfs.second)
-    {
-      Trace("nl-ext-cm") << "- Term: " << tf << std::endl;
-      bool success = true;
-      // tf is Figure 3 : tf( x )
-      Node bl;
-      Node bu;
-      if (k == PI)
-      {
-        bl = d_pi_bound[0];
-        bu = d_pi_bound[1];
-      }
-      else
-      {
-        std::pair<Node, Node> bounds = getTfModelBounds(tf, d_taylor_degree);
-        bl = bounds.first;
-        bu = bounds.second;
-        if (bl != bu)
-        {
-          d_model.setUsedApproximate();
-        }
-      }
-      if (!bl.isNull() && !bu.isNull())
-      {
-        // for each function in the congruence classe
-        for (const Node& ctf : d_funcCongClass[tf])
-        {
-          // each term in congruence classes should be master terms
-          Assert(d_trSlaves.find(ctf) != d_trSlaves.end());
-          // we set the bounds for each slave of tf
-          for (const Node& stf : d_trSlaves[ctf])
-          {
-            Trace("nl-ext-cm") << "...bound for " << stf << " : [" << bl << ", "
-                               << bu << "]" << std::endl;
-            success = d_model.addCheckModelBound(stf, bl, bu);
-          }
-        }
-      }
-      else
-      {
-        Trace("nl-ext-cm") << "...no bound for " << tf << std::endl;
-      }
-      if (!success)
-      {
-        // a bound was conflicting
-        Trace("nl-ext-cm") << "...failed to set bound for " << tf << std::endl;
-        Trace("nl-ext-cm") << "-----" << std::endl;
-        return false;
-      }
-    }
-  }
-  // replace the assertions
-  assertions = passertions;
-  return true;
-}
-
-void TranscendentalSolver::incrementTaylorDegree() { d_taylor_degree++; }
-unsigned TranscendentalSolver::getTaylorDegree() const
-{
-  return d_taylor_degree;
-}
-
-void TranscendentalSolver::processSideEffect(const NlLemma& se)
-{
-  for (const std::tuple<Node, unsigned, Node>& sp : se.d_secantPoint)
-  {
-    Node tf = std::get<0>(sp);
-    unsigned d = std::get<1>(sp);
-    Node c = std::get<2>(sp);
-    d_secant_points[tf][d].push_back(c);
-  }
-}
-
-void TranscendentalSolver::mkPi()
-{
-  NodeManager* nm = NodeManager::currentNM();
-  if (d_pi.isNull())
-  {
-    d_pi = nm->mkNullaryOperator(nm->realType(), PI);
-    d_pi_2 = Rewriter::rewrite(
-        nm->mkNode(MULT, d_pi, nm->mkConst(Rational(1) / Rational(2))));
-    d_pi_neg_2 = Rewriter::rewrite(
-        nm->mkNode(MULT, d_pi, nm->mkConst(Rational(-1) / Rational(2))));
-    d_pi_neg =
-        Rewriter::rewrite(nm->mkNode(MULT, d_pi, nm->mkConst(Rational(-1))));
-    // initialize bounds
-    d_pi_bound[0] = nm->mkConst(Rational(103993) / Rational(33102));
-    d_pi_bound[1] = nm->mkConst(Rational(104348) / Rational(33215));
-  }
-}
-
-void TranscendentalSolver::getCurrentPiBounds()
-{
-  NodeManager* nm = NodeManager::currentNM();
-  Node pi_lem = nm->mkNode(AND,
-                           nm->mkNode(GEQ, d_pi, d_pi_bound[0]),
-                           nm->mkNode(LEQ, d_pi, d_pi_bound[1]));
-  d_im.addPendingArithLemma(pi_lem, InferenceId::NL_T_PI_BOUND);
-}
-
-void TranscendentalSolver::checkTranscendentalInitialRefine()
-{
-  NodeManager* nm = NodeManager::currentNM();
-  Trace("nl-ext")
-      << "Get initial refinement lemmas for transcendental functions..."
-      << std::endl;
-  for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
-  {
-    Kind k = tfl.first;
-    for (const Node& t : tfl.second)
-    {
-      // initial refinements
-      if (d_tf_initial_refine.find(t) == d_tf_initial_refine.end())
-      {
-        d_tf_initial_refine[t] = true;
-        Node lem;
-        if (k == SINE)
-        {
-          Node symn = nm->mkNode(SINE, nm->mkNode(MULT, d_neg_one, t[0]));
-          symn = Rewriter::rewrite(symn);
-          // Can assume it is its own master since phase is split over 0,
-          // hence  -pi <= t[0] <= pi implies -pi <= -t[0] <= pi.
-          d_trMaster[symn] = symn;
-          d_trSlaves[symn].insert(symn);
-          Assert(d_trSlaves.find(t) != d_trSlaves.end());
-          std::vector<Node> children;
-
-          lem = nm->mkNode(AND,
-                           // bounds
-                           nm->mkNode(AND,
-                                      nm->mkNode(LEQ, t, d_one),
-                                      nm->mkNode(GEQ, t, d_neg_one)),
-                           // symmetry
-                           nm->mkNode(PLUS, t, symn).eqNode(d_zero),
-                           // sign
-                           nm->mkNode(EQUAL,
-                                      nm->mkNode(LT, t[0], d_zero),
-                                      nm->mkNode(LT, t, d_zero)),
-                           // zero val
-                           nm->mkNode(EQUAL,
-                                      nm->mkNode(GT, t[0], d_zero),
-                                      nm->mkNode(GT, t, d_zero)));
-          lem = nm->mkNode(
-              AND,
-              lem,
-              // zero tangent
-              nm->mkNode(AND,
-                         nm->mkNode(IMPLIES,
-                                    nm->mkNode(GT, t[0], d_zero),
-                                    nm->mkNode(LT, t, t[0])),
-                         nm->mkNode(IMPLIES,
-                                    nm->mkNode(LT, t[0], d_zero),
-                                    nm->mkNode(GT, t, t[0]))),
-              // pi tangent
-              nm->mkNode(
-                  AND,
-                  nm->mkNode(IMPLIES,
-                             nm->mkNode(LT, t[0], d_pi),
-                             nm->mkNode(LT, t, nm->mkNode(MINUS, d_pi, t[0]))),
-                  nm->mkNode(
-                      IMPLIES,
-                      nm->mkNode(GT, t[0], d_pi_neg),
-                      nm->mkNode(GT, t, nm->mkNode(MINUS, d_pi_neg, t[0])))));
-        }
-        else if (k == EXPONENTIAL)
-        {
-          // ( exp(x) > 0 ) ^ ( x=0 <=> exp( x ) = 1 ) ^ ( x < 0 <=> exp( x ) <
-          // 1 ) ^ ( x <= 0 V exp( x ) > x + 1 )
-          lem = nm->mkNode(
-              AND,
-              nm->mkNode(GT, t, d_zero),
-              nm->mkNode(EQUAL, t[0].eqNode(d_zero), t.eqNode(d_one)),
-              nm->mkNode(EQUAL,
-                         nm->mkNode(LT, t[0], d_zero),
-                         nm->mkNode(LT, t, d_one)),
-              nm->mkNode(OR,
-                         nm->mkNode(LEQ, t[0], d_zero),
-                         nm->mkNode(GT, t, nm->mkNode(PLUS, t[0], d_one))));
-        }
-        if (!lem.isNull())
-        {
-          d_im.addPendingArithLemma(lem, InferenceId::NL_T_INIT_REFINE);
-        }
-      }
-    }
-  }
-}
-
-void TranscendentalSolver::checkTranscendentalMonotonic()
-{
-  Trace("nl-ext") << "Get monotonicity lemmas for transcendental functions..."
-                  << std::endl;
-
-  // sort arguments of all transcendentals
-  std::map<Kind, std::vector<Node> > sorted_tf_args;
-  std::map<Kind, std::map<Node, Node> > tf_arg_to_term;
-
-  for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
-  {
-    Kind k = tfl.first;
-    if (k == EXPONENTIAL || k == SINE)
-    {
-      for (const Node& tf : tfl.second)
-      {
-        Node a = tf[0];
-        Node mvaa = d_model.computeAbstractModelValue(a);
-        if (mvaa.isConst())
-        {
-          Trace("nl-ext-tf-mono-debug") << "...tf term : " << a << std::endl;
-          sorted_tf_args[k].push_back(a);
-          tf_arg_to_term[k][a] = tf;
-        }
-      }
-    }
-  }
-
-  SortNlModel smv;
-  smv.d_nlm = &d_model;
-  // sort by concrete values
-  smv.d_isConcrete = true;
-  smv.d_reverse_order = true;
-  for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
-  {
-    Kind k = tfl.first;
-    if (!sorted_tf_args[k].empty())
-    {
-      std::sort(sorted_tf_args[k].begin(), sorted_tf_args[k].end(), smv);
-      Trace("nl-ext-tf-mono") << "Sorted transcendental function list for " << k
-                              << " : " << std::endl;
-      for (unsigned i = 0; i < sorted_tf_args[k].size(); i++)
-      {
-        Node targ = sorted_tf_args[k][i];
-        Node mvatarg = d_model.computeAbstractModelValue(targ);
-        Trace("nl-ext-tf-mono")
-            << "  " << targ << " -> " << mvatarg << std::endl;
-        Node t = tf_arg_to_term[k][targ];
-        Node mvat = d_model.computeAbstractModelValue(t);
-        Trace("nl-ext-tf-mono") << "     f-val : " << mvat << std::endl;
-      }
-      std::vector<Node> mpoints;
-      std::vector<Node> mpoints_vals;
-      if (k == SINE)
-      {
-        mpoints.push_back(d_pi);
-        mpoints.push_back(d_pi_2);
-        mpoints.push_back(d_zero);
-        mpoints.push_back(d_pi_neg_2);
-        mpoints.push_back(d_pi_neg);
-      }
-      else if (k == EXPONENTIAL)
-      {
-        mpoints.push_back(Node::null());
-      }
-      if (!mpoints.empty())
-      {
-        // get model values for points
-        for (unsigned i = 0; i < mpoints.size(); i++)
-        {
-          Node mpv;
-          if (!mpoints[i].isNull())
-          {
-            mpv = d_model.computeAbstractModelValue(mpoints[i]);
-            Assert(mpv.isConst());
-          }
-          mpoints_vals.push_back(mpv);
-        }
-
-        unsigned mdir_index = 0;
-        int monotonic_dir = -1;
-        Node mono_bounds[2];
-        Node targ, targval, t, tval;
-        for (unsigned i = 0, size = sorted_tf_args[k].size(); i < size; i++)
-        {
-          Node sarg = sorted_tf_args[k][i];
-          Node sargval = d_model.computeAbstractModelValue(sarg);
-          Assert(sargval.isConst());
-          Node s = tf_arg_to_term[k][sarg];
-          Node sval = d_model.computeAbstractModelValue(s);
-          Assert(sval.isConst());
-
-          // increment to the proper monotonicity region
-          bool increment = true;
-          while (increment && mdir_index < mpoints.size())
-          {
-            increment = false;
-            if (mpoints[mdir_index].isNull())
-            {
-              increment = true;
-            }
-            else
-            {
-              Node pval = mpoints_vals[mdir_index];
-              Assert(pval.isConst());
-              if (sargval.getConst<Rational>() < pval.getConst<Rational>())
-              {
-                increment = true;
-                Trace("nl-ext-tf-mono") << "...increment at " << sarg
-                                        << " since model value is less than "
-                                        << mpoints[mdir_index] << std::endl;
-              }
-            }
-            if (increment)
-            {
-              tval = Node::null();
-              mono_bounds[1] = mpoints[mdir_index];
-              mdir_index++;
-              monotonic_dir = regionToMonotonicityDir(k, mdir_index);
-              if (mdir_index < mpoints.size())
-              {
-                mono_bounds[0] = mpoints[mdir_index];
-              }
-              else
-              {
-                mono_bounds[0] = Node::null();
-              }
-            }
-          }
-          // store the concavity region
-          d_tf_region[s] = mdir_index;
-          Trace("nl-ext-concavity") << "Transcendental function " << s
-                                    << " is in region #" << mdir_index;
-          Trace("nl-ext-concavity")
-              << ", arg model value = " << sargval << std::endl;
-
-          if (!tval.isNull())
-          {
-            NodeManager* nm = NodeManager::currentNM();
-            Node mono_lem;
-            if (monotonic_dir == 1
-                && sval.getConst<Rational>() > tval.getConst<Rational>())
-            {
-              mono_lem = nm->mkNode(
-                  IMPLIES, nm->mkNode(GEQ, targ, sarg), nm->mkNode(GEQ, t, s));
-            }
-            else if (monotonic_dir == -1
-                     && sval.getConst<Rational>() < tval.getConst<Rational>())
-            {
-              mono_lem = nm->mkNode(
-                  IMPLIES, nm->mkNode(LEQ, targ, sarg), nm->mkNode(LEQ, t, s));
-            }
-            if (!mono_lem.isNull())
-            {
-              if (!mono_bounds[0].isNull())
-              {
-                Assert(!mono_bounds[1].isNull());
-                mono_lem = nm->mkNode(
-                    IMPLIES,
-                    nm->mkNode(AND,
-                               mkBounded(mono_bounds[0], targ, mono_bounds[1]),
-                               mkBounded(mono_bounds[0], sarg, mono_bounds[1])),
-                    mono_lem);
-              }
-              Trace("nl-ext-tf-mono")
-                  << "Monotonicity lemma : " << mono_lem << std::endl;
-
-              d_im.addPendingArithLemma(mono_lem, InferenceId::NL_T_MONOTONICITY);
-            }
-          }
-          // store the previous values
-          targ = sarg;
-          targval = sargval;
-          t = s;
-          tval = sval;
-        }
-      }
-    }
-  }
-}
-
-void TranscendentalSolver::checkTranscendentalTangentPlanes()
-{
-  Trace("nl-ext") << "Get tangent plane lemmas for transcendental functions..."
-                  << std::endl;
-  // this implements Figure 3 of "Satisfiaility Modulo Transcendental Functions
-  // via Incremental Linearization" by Cimatti et al
-  for (std::pair<const Kind, std::vector<Node> >& tfs : d_funcMap)
-  {
-    Kind k = tfs.first;
-    if (k == PI)
-    {
-      // We do not use Taylor approximation for PI currently.
-      // This is because the convergence is extremely slow, and hence an
-      // initial approximation is superior.
-      continue;
-    }
-    Trace("nl-ext-tftp-debug2") << "Taylor variables: " << std::endl;
-    Trace("nl-ext-tftp-debug2")
-        << "          taylor_real_fv : " << d_taylor_real_fv << std::endl;
-    Trace("nl-ext-tftp-debug2")
-        << "     taylor_real_fv_base : " << d_taylor_real_fv_base << std::endl;
-    Trace("nl-ext-tftp-debug2")
-        << " taylor_real_fv_base_rem : " << d_taylor_real_fv_base_rem
-        << std::endl;
-    Trace("nl-ext-tftp-debug2") << std::endl;
-
-    // we substitute into the Taylor sum P_{n,f(0)}( x )
-
-    for (const Node& tf : tfs.second)
-    {
-      // tf is Figure 3 : tf( x )
-      Trace("nl-ext-tftp") << "Compute tangent planes " << tf << std::endl;
-      // go until max degree is reached, or we don't meet bound criteria
-      for (unsigned d = 1; d <= d_taylor_degree; d++)
-      {
-        Trace("nl-ext-tftp") << "- run at degree " << d << "..." << std::endl;
-        unsigned prev = d_im.numPendingLemmas() + d_im.numWaitingLemmas();
-        if (checkTfTangentPlanesFun(tf, d))
-        {
-          Trace("nl-ext-tftp")
-              << "...fail, #lemmas = "
-              << (d_im.numPendingLemmas() + d_im.numWaitingLemmas() - prev)
-              << std::endl;
-          break;
-        }
-        else
-        {
-          Trace("nl-ext-tftp") << "...success" << std::endl;
-        }
-      }
-    }
-  }
-}
-
-bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf,
-                                                   unsigned d)
-{
-  NodeManager* nm = NodeManager::currentNM();
-  Kind k = tf.getKind();
-  // this should only be run on master applications
-  Assert(d_trSlaves.find(tf) != d_trSlaves.end());
-
-  // Figure 3 : c
-  Node c = d_model.computeAbstractModelValue(tf[0]);
-  int csign = c.getConst<Rational>().sgn();
-  if (csign == 0)
-  {
-    // no secant/tangent plane is necessary
-    return true;
-  }
-  Assert(csign == 1 || csign == -1);
-
-  // Figure 3: P_l, P_u
-  // mapped to for signs of c
-  std::map<int, Node> poly_approx_bounds[2];
-  std::vector<Node> pbounds;
-  getPolynomialApproximationBoundForArg(k, c, d, pbounds);
-  poly_approx_bounds[0][1] = pbounds[0];
-  poly_approx_bounds[0][-1] = pbounds[1];
-  poly_approx_bounds[1][1] = pbounds[2];
-  poly_approx_bounds[1][-1] = pbounds[3];
-
-  // Figure 3 : v
-  Node v = d_model.computeAbstractModelValue(tf);
-
-  // check value of tf
-  Trace("nl-ext-tftp-debug") << "Process tangent plane refinement for " << tf
-                             << ", degree " << d << "..." << std::endl;
-  Trace("nl-ext-tftp-debug") << "  value in model : " << v << std::endl;
-  Trace("nl-ext-tftp-debug") << "  arg value in model : " << c << std::endl;
-
-  std::vector<Node> taylor_vars;
-  taylor_vars.push_back(d_taylor_real_fv);
-
-  // compute the concavity
-  int region = -1;
-  std::unordered_map<Node, int, NodeHashFunction>::iterator itr =
-      d_tf_region.find(tf);
-  if (itr != d_tf_region.end())
-  {
-    region = itr->second;
-    Trace("nl-ext-tftp-debug") << "  region is : " << region << std::endl;
-  }
-  // Figure 3 : conc
-  int concavity = regionToConcavity(k, itr->second);
-  Trace("nl-ext-tftp-debug") << "  concavity is : " << concavity << std::endl;
-  if (concavity == 0)
-  {
-    // no secant/tangent plane is necessary
-    return true;
-  }
-  // bounds for which we are this concavity
-  // Figure 3: < l, u >
-  Node bounds[2];
-  if (k == SINE)
-  {
-    bounds[0] = regionToLowerBound(k, region);
-    Assert(!bounds[0].isNull());
-    bounds[1] = regionToUpperBound(k, region);
-    Assert(!bounds[1].isNull());
-  }
-
-  // Figure 3: P
-  Node poly_approx;
-
-  // compute whether this is a tangent refinement or a secant refinement
-  bool is_tangent = false;
-  bool is_secant = false;
-  std::pair<Node, Node> mvb = getTfModelBounds(tf, d);
-  // this is the approximated value of tf(c), which is a value such that:
-  //    M_A(tf(c)) >= poly_appox_c >= tf(c) or
-  //    M_A(tf(c)) <= poly_appox_c <= tf(c)
-  // In other words, it is a better approximation of the true value of tf(c)
-  // in the case that we add a refinement lemma. We use this value in the
-  // refinement schemas below.
-  Node poly_approx_c;
-  for (unsigned r = 0; r < 2; r++)
-  {
-    Node pab = poly_approx_bounds[r][csign];
-    Node v_pab = r == 0 ? mvb.first : mvb.second;
-    if (!v_pab.isNull())
-    {
-      Trace("nl-ext-tftp-debug2")
-          << "...model value of " << pab << " is " << v_pab << std::endl;
-
-      Assert(v_pab.isConst());
-      Node comp = nm->mkNode(r == 0 ? LT : GT, v, v_pab);
-      Trace("nl-ext-tftp-debug2") << "...compare : " << comp << std::endl;
-      Node compr = Rewriter::rewrite(comp);
-      Trace("nl-ext-tftp-debug2") << "...got : " << compr << std::endl;
-      if (compr == d_true)
-      {
-        poly_approx_c = v_pab;
-        // beyond the bounds
-        if (r == 0)
-        {
-          poly_approx = poly_approx_bounds[r][csign];
-          is_tangent = concavity == 1;
-          is_secant = concavity == -1;
-        }
-        else
-        {
-          poly_approx = poly_approx_bounds[r][csign];
-          is_tangent = concavity == -1;
-          is_secant = concavity == 1;
-        }
-        if (Trace.isOn("nl-ext-tftp"))
-        {
-          Trace("nl-ext-tftp") << "*** Outside boundary point (";
-          Trace("nl-ext-tftp") << (r == 0 ? "low" : "high") << ") ";
-          printRationalApprox("nl-ext-tftp", v_pab);
-          Trace("nl-ext-tftp") << ", will refine..." << std::endl;
-          Trace("nl-ext-tftp")
-              << "    poly_approx = " << poly_approx << std::endl;
-          Trace("nl-ext-tftp")
-              << "    is_tangent = " << is_tangent << std::endl;
-          Trace("nl-ext-tftp") << "    is_secant = " << is_secant << std::endl;
-        }
-        break;
-      }
-      else
-      {
-        Trace("nl-ext-tftp")
-            << "  ...within " << (r == 0 ? "low" : "high") << " bound : ";
-        printRationalApprox("nl-ext-tftp", v_pab);
-        Trace("nl-ext-tftp") << std::endl;
-      }
-    }
-  }
-
-  // Figure 3: P( c )
-  if (is_tangent || is_secant)
-  {
-    Trace("nl-ext-tftp-debug2")
-        << "...poly approximation at c is " << poly_approx_c << std::endl;
-  }
-  else
-  {
-    // we may want to continue getting better bounds
-    return false;
-  }
-
-  if (is_tangent)
-  {
-    // compute tangent plane
-    // Figure 3: T( x )
-    // We use zero slope tangent planes, since the concavity of the Taylor
-    // approximation cannot be easily established.
-    Node tplane = poly_approx_c;
-
-    Node lem = nm->mkNode(concavity == 1 ? GEQ : LEQ, tf, tplane);
-    std::vector<Node> antec;
-    int mdir = regionToMonotonicityDir(k, region);
-    for (unsigned i = 0; i < 2; i++)
-    {
-      // Tangent plane is valid in the interval [c,u) if the slope of the
-      // function matches its concavity, and is valid in (l, c] otherwise.
-      Node use_bound = (mdir == concavity) == (i == 0) ? c : bounds[i];
-      if (!use_bound.isNull())
-      {
-        Node ant = nm->mkNode(i == 0 ? GEQ : LEQ, tf[0], use_bound);
-        antec.push_back(ant);
-      }
-    }
-    if (!antec.empty())
-    {
-      Node antec_n = antec.size() == 1 ? antec[0] : nm->mkNode(AND, antec);
-      lem = nm->mkNode(IMPLIES, antec_n, lem);
-    }
-    Trace("nl-ext-tftp-debug2")
-        << "*** Tangent plane lemma (pre-rewrite): " << lem << std::endl;
-    lem = Rewriter::rewrite(lem);
-    Trace("nl-ext-tftp-lemma")
-        << "*** Tangent plane lemma : " << lem << std::endl;
-    Assert(d_model.computeAbstractModelValue(lem) == d_false);
-    // Figure 3 : line 9
-    d_im.addPendingArithLemma(lem, InferenceId::NL_T_TANGENT, nullptr, true);
-  }
-  else if (is_secant)
-  {
-    // bounds are the minimum and maximum previous secant points
-    // should not repeat secant points: secant lemmas should suffice to
-    // rule out previous assignment
-    Assert(std::find(
-               d_secant_points[tf][d].begin(), d_secant_points[tf][d].end(), c)
-           == d_secant_points[tf][d].end());
-    // Insert into the (temporary) vector. We do not update this vector
-    // until we are sure this secant plane lemma has been processed. We do
-    // this by mapping the lemma to a side effect below.
-    std::vector<Node> spoints = d_secant_points[tf][d];
-    spoints.push_back(c);
-
-    // sort
-    SortNlModel smv;
-    smv.d_nlm = &d_model;
-    smv.d_isConcrete = true;
-    std::sort(spoints.begin(), spoints.end(), smv);
-    // get the resulting index of c
-    unsigned index =
-        std::find(spoints.begin(), spoints.end(), c) - spoints.begin();
-    // bounds are the next closest upper/lower bound values
-    if (index > 0)
-    {
-      bounds[0] = spoints[index - 1];
-    }
-    else
-    {
-      // otherwise, we use the lower boundary point for this concavity
-      // region
-      if (k == SINE)
-      {
-        Assert(!bounds[0].isNull());
-      }
-      else if (k == EXPONENTIAL)
-      {
-        // pick c-1
-        bounds[0] = Rewriter::rewrite(nm->mkNode(MINUS, c, d_one));
-      }
-    }
-    if (index < spoints.size() - 1)
-    {
-      bounds[1] = spoints[index + 1];
-    }
-    else
-    {
-      // otherwise, we use the upper boundary point for this concavity
-      // region
-      if (k == SINE)
-      {
-        Assert(!bounds[1].isNull());
-      }
-      else if (k == EXPONENTIAL)
-      {
-        // pick c+1
-        bounds[1] = Rewriter::rewrite(nm->mkNode(PLUS, c, d_one));
-      }
-    }
-    Trace("nl-ext-tftp-debug2") << "...secant bounds are : " << bounds[0]
-                                << " ... " << bounds[1] << std::endl;
-
-    // the secant plane may be conjunction of 1-2 guarded inequalities
-    std::vector<Node> lemmaConj;
-    for (unsigned s = 0; s < 2; s++)
-    {
-      // compute secant plane
-      Assert(!poly_approx.isNull());
-      Assert(!bounds[s].isNull());
-      // take the model value of l or u (since may contain PI)
-      Node b = d_model.computeAbstractModelValue(bounds[s]);
-      Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds[s]
-                                  << " is " << b << std::endl;
-      Assert(b.isConst());
-      if (c != b)
-      {
-        // Figure 3 : P(l), P(u), for s = 0,1
-        Node poly_approx_b;
-        std::vector<Node> taylor_subs;
-        taylor_subs.push_back(b);
-        Assert(taylor_vars.size() == taylor_subs.size());
-        poly_approx_b = poly_approx.substitute(taylor_vars.begin(),
-                                               taylor_vars.end(),
-                                               taylor_subs.begin(),
-                                               taylor_subs.end());
-        // Figure 3: S_l( x ), S_u( x ) for s = 0,1
-        Node splane;
-        Node rcoeff_n = Rewriter::rewrite(nm->mkNode(MINUS, b, c));
-        Assert(rcoeff_n.isConst());
-        Rational rcoeff = rcoeff_n.getConst<Rational>();
-        Assert(rcoeff.sgn() != 0);
-        poly_approx_b = Rewriter::rewrite(poly_approx_b);
-        poly_approx_c = Rewriter::rewrite(poly_approx_c);
-        splane = nm->mkNode(
-            PLUS,
-            poly_approx_b,
-            nm->mkNode(MULT,
-                       nm->mkNode(MINUS, poly_approx_b, poly_approx_c),
-                       nm->mkConst(Rational(1) / rcoeff),
-                       nm->mkNode(MINUS, tf[0], b)));
-
-        Node lem = nm->mkNode(concavity == 1 ? LEQ : GEQ, tf, splane);
-        // With respect to Figure 3, this is slightly different.
-        // In particular, we chose b to be the model value of bounds[s],
-        // which is a constant although bounds[s] may not be (e.g. if it
-        // contains PI).
-        // To ensure that c...b does not cross an inflection point,
-        // we guard with the symbolic version of bounds[s].
-        // This leads to lemmas e.g. of this form:
-        //   ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x -
-        //   b ) + P( b ) )
-        // where b = (PI/2)^M, the current value of PI/2 in the model.
-        // This is sound since we are guarded by the symbolic
-        // representation of PI/2.
-        Node antec_n =
-            nm->mkNode(AND,
-                       nm->mkNode(GEQ, tf[0], s == 0 ? bounds[s] : c),
-                       nm->mkNode(LEQ, tf[0], s == 0 ? c : bounds[s]));
-        lem = nm->mkNode(IMPLIES, antec_n, lem);
-        Trace("nl-ext-tftp-debug2")
-            << "*** Secant plane lemma (pre-rewrite) : " << lem << std::endl;
-        lem = Rewriter::rewrite(lem);
-        Trace("nl-ext-tftp-lemma")
-            << "*** Secant plane lemma : " << lem << std::endl;
-        lemmaConj.push_back(lem);
-        Assert(d_model.computeAbstractModelValue(lem) == d_false);
-      }
-    }
-    // Figure 3 : line 22
-    Assert(!lemmaConj.empty());
-    Node lem =
-        lemmaConj.size() == 1 ? lemmaConj[0] : nm->mkNode(AND, lemmaConj);
-    NlLemma nlem(lem, LemmaProperty::NONE, nullptr, InferenceId::NL_T_SECANT);
-    // The side effect says that if lem is added, then we should add the
-    // secant point c for (tf,d).
-    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, c));
-    d_im.addPendingArithLemma(nlem, true);
-  }
-  return true;
-}
-
-int TranscendentalSolver::regionToMonotonicityDir(Kind k, int region)
-{
-  if (k == EXPONENTIAL)
-  {
-    if (region == 1)
-    {
-      return 1;
-    }
-  }
-  else if (k == SINE)
-  {
-    if (region == 1 || region == 4)
-    {
-      return -1;
-    }
-    else if (region == 2 || region == 3)
-    {
-      return 1;
-    }
-  }
-  return 0;
-}
-
-int TranscendentalSolver::regionToConcavity(Kind k, int region)
-{
-  if (k == EXPONENTIAL)
-  {
-    if (region == 1)
-    {
-      return 1;
-    }
-  }
-  else if (k == SINE)
-  {
-    if (region == 1 || region == 2)
-    {
-      return -1;
-    }
-    else if (region == 3 || region == 4)
-    {
-      return 1;
-    }
-  }
-  return 0;
-}
-
-Node TranscendentalSolver::regionToLowerBound(Kind k, int region)
-{
-  if (k == SINE)
-  {
-    if (region == 1)
-    {
-      return d_pi_2;
-    }
-    else if (region == 2)
-    {
-      return d_zero;
-    }
-    else if (region == 3)
-    {
-      return d_pi_neg_2;
-    }
-    else if (region == 4)
-    {
-      return d_pi_neg;
-    }
-  }
-  return Node::null();
-}
-
-Node TranscendentalSolver::regionToUpperBound(Kind k, int region)
-{
-  if (k == SINE)
-  {
-    if (region == 1)
-    {
-      return d_pi;
-    }
-    else if (region == 2)
-    {
-      return d_pi_2;
-    }
-    else if (region == 3)
-    {
-      return d_zero;
-    }
-    else if (region == 4)
-    {
-      return d_pi_neg_2;
-    }
-  }
-  return Node::null();
-}
-
-Node TranscendentalSolver::getDerivative(Node n, Node x)
-{
-  NodeManager* nm = NodeManager::currentNM();
-  Assert(x.isVar());
-  // only handle the cases of the taylor expansion of d
-  if (n.getKind() == EXPONENTIAL)
-  {
-    if (n[0] == x)
-    {
-      return n;
-    }
-  }
-  else if (n.getKind() == SINE)
-  {
-    if (n[0] == x)
-    {
-      Node na = nm->mkNode(MINUS, d_pi_2, n[0]);
-      Node ret = nm->mkNode(SINE, na);
-      ret = Rewriter::rewrite(ret);
-      return ret;
-    }
-  }
-  else if (n.getKind() == PLUS)
-  {
-    std::vector<Node> dchildren;
-    for (unsigned i = 0; i < n.getNumChildren(); i++)
-    {
-      // PLUS is flattened in rewriter, recursion depth is bounded by 1
-      Node dc = getDerivative(n[i], x);
-      if (dc.isNull())
-      {
-        return dc;
-      }
-      else
-      {
-        dchildren.push_back(dc);
-      }
-    }
-    return nm->mkNode(PLUS, dchildren);
-  }
-  else if (n.getKind() == MULT)
-  {
-    Assert(n[0].isConst());
-    Node dc = getDerivative(n[1], x);
-    if (!dc.isNull())
-    {
-      return nm->mkNode(MULT, n[0], dc);
-    }
-  }
-  else if (n.getKind() == NONLINEAR_MULT)
-  {
-    unsigned xcount = 0;
-    std::vector<Node> children;
-    unsigned xindex = 0;
-    for (unsigned i = 0, size = n.getNumChildren(); i < size; i++)
-    {
-      if (n[i] == x)
-      {
-        xcount++;
-        xindex = i;
-      }
-      children.push_back(n[i]);
-    }
-    if (xcount == 0)
-    {
-      return d_zero;
-    }
-    else
-    {
-      children[xindex] = nm->mkConst(Rational(xcount));
-    }
-    return nm->mkNode(MULT, children);
-  }
-  else if (n.isVar())
-  {
-    return n == x ? d_one : d_zero;
-  }
-  else if (n.isConst())
-  {
-    return d_zero;
-  }
-  Trace("nl-ext-debug") << "No derivative computed for " << n;
-  Trace("nl-ext-debug") << " for d/d{" << x << "}" << std::endl;
-  return Node::null();
-}
-
-std::pair<Node, Node> TranscendentalSolver::getTaylor(Node fa, unsigned n)
-{
-  NodeManager* nm = NodeManager::currentNM();
-  Assert(n > 0);
-  Node fac;  // what term we cache for fa
-  if (fa[0] == d_zero)
-  {
-    // optimization : simpler to compute (x-fa[0])^n if we are centered around 0
-    fac = fa;
-  }
-  else
-  {
-    // otherwise we use a standard factor a in (x-a)^n
-    fac = nm->mkNode(fa.getKind(), d_taylor_real_fv_base);
-  }
-  Node taylor_rem;
-  Node taylor_sum;
-  // check if we have already computed this Taylor series
-  std::unordered_map<unsigned, Node>::iterator itt = d_taylor_sum[fac].find(n);
-  if (itt == d_taylor_sum[fac].end())
-  {
-    Node i_exp_base;
-    if (fa[0] == d_zero)
-    {
-      i_exp_base = d_taylor_real_fv;
-    }
-    else
-    {
-      i_exp_base = Rewriter::rewrite(
-          nm->mkNode(MINUS, d_taylor_real_fv, d_taylor_real_fv_base));
-    }
-    Node i_derv = fac;
-    Node i_fact = d_one;
-    Node i_exp = d_one;
-    int i_derv_status = 0;
-    unsigned counter = 0;
-    std::vector<Node> sum;
-    do
-    {
-      counter++;
-      if (fa.getKind() == EXPONENTIAL)
-      {
-        // unchanged
-      }
-      else if (fa.getKind() == SINE)
-      {
-        if (i_derv_status % 2 == 1)
-        {
-          Node arg = nm->mkNode(PLUS, d_pi_2, d_taylor_real_fv_base);
-          i_derv = nm->mkNode(SINE, arg);
-        }
-        else
-        {
-          i_derv = fa;
-        }
-        if (i_derv_status >= 2)
-        {
-          i_derv = nm->mkNode(MINUS, d_zero, i_derv);
-        }
-        i_derv = Rewriter::rewrite(i_derv);
-        i_derv_status = i_derv_status == 3 ? 0 : i_derv_status + 1;
-      }
-      if (counter == (n + 1))
-      {
-        TNode x = d_taylor_real_fv_base;
-        i_derv = i_derv.substitute(x, d_taylor_real_fv_base_rem);
-      }
-      Node curr = nm->mkNode(MULT, nm->mkNode(DIVISION, i_derv, i_fact), i_exp);
-      if (counter == (n + 1))
-      {
-        taylor_rem = curr;
-      }
-      else
-      {
-        sum.push_back(curr);
-        i_fact = Rewriter::rewrite(
-            nm->mkNode(MULT, nm->mkConst(Rational(counter)), i_fact));
-        i_exp = Rewriter::rewrite(nm->mkNode(MULT, i_exp_base, i_exp));
-      }
-    } while (counter <= n);
-    taylor_sum = sum.size() == 1 ? sum[0] : nm->mkNode(PLUS, sum);
-
-    if (fac[0] != d_taylor_real_fv_base)
-    {
-      TNode x = d_taylor_real_fv_base;
-      taylor_sum = taylor_sum.substitute(x, fac[0]);
-    }
-
-    // cache
-    d_taylor_sum[fac][n] = taylor_sum;
-    d_taylor_rem[fac][n] = taylor_rem;
-  }
-  else
-  {
-    taylor_sum = itt->second;
-    Assert(d_taylor_rem[fac].find(n) != d_taylor_rem[fac].end());
-    taylor_rem = d_taylor_rem[fac][n];
-  }
-
-  // must substitute for the argument if we were using a different lookup
-  if (fa[0] != fac[0])
-  {
-    TNode x = d_taylor_real_fv_base;
-    taylor_sum = taylor_sum.substitute(x, fa[0]);
-  }
-  return std::pair<Node, Node>(taylor_sum, taylor_rem);
-}
-
-void TranscendentalSolver::getPolynomialApproximationBounds(
-    Kind k, unsigned d, std::vector<Node>& pbounds)
-{
-  if (d_poly_bounds[k][d].empty())
-  {
-    NodeManager* nm = NodeManager::currentNM();
-    Node tft = nm->mkNode(k, d_zero);
-    // n is the Taylor degree we are currently considering
-    unsigned n = 2 * d;
-    // n must be even
-    std::pair<Node, Node> taylor = getTaylor(tft, n);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor for " << k << " is : " << taylor.first << std::endl;
-    Node taylor_sum = Rewriter::rewrite(taylor.first);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor for " << k << " is (post-rewrite) : " << taylor_sum
-        << std::endl;
-    Assert(taylor.second.getKind() == MULT);
-    Assert(taylor.second.getNumChildren() == 2);
-    Assert(taylor.second[0].getKind() == DIVISION);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor remainder for " << k << " is " << taylor.second << std::endl;
-    // ru is x^{n+1}/(n+1)!
-    Node ru = nm->mkNode(DIVISION, taylor.second[1], taylor.second[0][1]);
-    ru = Rewriter::rewrite(ru);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor remainder factor is (post-rewrite) : " << ru << std::endl;
-    if (k == EXPONENTIAL)
-    {
-      pbounds.push_back(taylor_sum);
-      pbounds.push_back(taylor_sum);
-      pbounds.push_back(Rewriter::rewrite(
-          nm->mkNode(MULT, taylor_sum, nm->mkNode(PLUS, d_one, ru))));
-      pbounds.push_back(Rewriter::rewrite(nm->mkNode(PLUS, taylor_sum, ru)));
-    }
-    else
-    {
-      Assert(k == SINE);
-      Node l = Rewriter::rewrite(nm->mkNode(MINUS, taylor_sum, ru));
-      Node u = Rewriter::rewrite(nm->mkNode(PLUS, taylor_sum, ru));
-      pbounds.push_back(l);
-      pbounds.push_back(l);
-      pbounds.push_back(u);
-      pbounds.push_back(u);
-    }
-    Trace("nl-ext-tf-tplanes")
-        << "Polynomial approximation for " << k << " is: " << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Lower (pos): " << pbounds[0] << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Upper (pos): " << pbounds[2] << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Lower (neg): " << pbounds[1] << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Upper (neg): " << pbounds[3] << std::endl;
-    d_poly_bounds[k][d].insert(
-        d_poly_bounds[k][d].end(), pbounds.begin(), pbounds.end());
-  }
-  else
-  {
-    pbounds.insert(
-        pbounds.end(), d_poly_bounds[k][d].begin(), d_poly_bounds[k][d].end());
-  }
-}
-
-void TranscendentalSolver::getPolynomialApproximationBoundForArg(
-    Kind k, Node c, unsigned d, std::vector<Node>& pbounds)
-{
-  getPolynomialApproximationBounds(k, d, pbounds);
-  Assert(c.isConst());
-  if (k == EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
-  {
-    NodeManager* nm = NodeManager::currentNM();
-    Node tft = nm->mkNode(k, d_zero);
-    bool success = false;
-    unsigned ds = d;
-    TNode ttrf = d_taylor_real_fv;
-    TNode tc = c;
-    do
-    {
-      success = true;
-      unsigned n = 2 * ds;
-      std::pair<Node, Node> taylor = getTaylor(tft, n);
-      // check that 1-c^{n+1}/(n+1)! > 0
-      Node ru = nm->mkNode(DIVISION, taylor.second[1], taylor.second[0][1]);
-      Node rus = ru.substitute(ttrf, tc);
-      rus = Rewriter::rewrite(rus);
-      Assert(rus.isConst());
-      if (rus.getConst<Rational>() > d_one.getConst<Rational>())
-      {
-        success = false;
-        ds = ds + 1;
-      }
-    } while (!success);
-    if (ds > d)
-    {
-      Trace("nl-ext-exp-taylor")
-          << "*** Increase Taylor bound to " << ds << " > " << d << " for ("
-          << k << " " << c << ")" << std::endl;
-      // must use sound upper bound
-      std::vector<Node> pboundss;
-      getPolynomialApproximationBounds(k, ds, pboundss);
-      pbounds[2] = pboundss[2];
-    }
-  }
-}
-
-std::pair<Node, Node> TranscendentalSolver::getTfModelBounds(Node tf,
-                                                             unsigned d)
-{
-  // compute the model value of the argument
-  Node c = d_model.computeAbstractModelValue(tf[0]);
-  Assert(c.isConst());
-  int csign = c.getConst<Rational>().sgn();
-  Kind k = tf.getKind();
-  if (csign == 0)
-  {
-    // at zero, its trivial
-    if (k == SINE)
-    {
-      return std::pair<Node, Node>(d_zero, d_zero);
-    }
-    Assert(k == EXPONENTIAL);
-    return std::pair<Node, Node>(d_one, d_one);
-  }
-  bool isNeg = csign == -1;
-
-  std::vector<Node> pbounds;
-  getPolynomialApproximationBoundForArg(k, c, d, pbounds);
-
-  std::vector<Node> bounds;
-  TNode tfv = d_taylor_real_fv;
-  TNode tfs = tf[0];
-  for (unsigned d2 = 0; d2 < 2; d2++)
-  {
-    int index = d2 == 0 ? (isNeg ? 1 : 0) : (isNeg ? 3 : 2);
-    Node pab = pbounds[index];
-    if (!pab.isNull())
-    {
-      // { x -> M_A(tf[0]) }
-      // Notice that we compute the model value of tfs first, so that
-      // the call to rewrite below does not modify the term, where notice that
-      // rewrite( x*x { x -> M_A(t) } ) = M_A(t)*M_A(t)
-      // is not equal to
-      // M_A( x*x { x -> t } ) = M_A( t*t )
-      // where M_A denotes the abstract model.
-      Node mtfs = d_model.computeAbstractModelValue(tfs);
-      pab = pab.substitute(tfv, mtfs);
-      pab = Rewriter::rewrite(pab);
-      Assert (pab.isConst());
-      bounds.push_back(pab);
-    }
-    else
-    {
-      bounds.push_back(Node::null());
-    }
-  }
-  return std::pair<Node, Node>(bounds[0], bounds[1]);
-}
-
-Node TranscendentalSolver::mkValidPhase(Node a, Node pi)
-{
-  return mkBounded(
-      NodeManager::currentNM()->mkNode(MULT, mkRationalNode(-1), pi), a, pi);
-}
-
-}  // namespace nl
-}  // namespace arith
-}  // namespace theory
-}  // namespace CVC4
diff --git a/src/theory/arith/nl/transcendental_solver.h b/src/theory/arith/nl/transcendental_solver.h
deleted file mode 100644
index 004221be2..000000000
--- a/src/theory/arith/nl/transcendental_solver.h
+++ /dev/null
@@ -1,423 +0,0 @@
-/*********************                                                        */
-/*! \file transcendental_solver.h
- ** \verbatim
- ** Top contributors (to current version):
- **   Andrew Reynolds, Tim King
- ** 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 Solving for handling transcendental functions.
- **/
-
-#ifndef CVC4__THEORY__ARITH__NL__TRANSCENDENTAL_SOLVER_H
-#define CVC4__THEORY__ARITH__NL__TRANSCENDENTAL_SOLVER_H
-
-#include <map>
-#include <unordered_map>
-#include <unordered_set>
-#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 {
-
-/** Transcendental solver class
- *
- * This class implements model-based refinement schemes
- * for transcendental functions, described in:
- *
- * - "Satisfiability Modulo Transcendental
- * Functions via Incremental Linearization" by Cimatti
- * et al., CADE 2017.
- *
- * It's main functionality are methods that implement lemma schemas below,
- * which return a set of lemmas that should be sent on the output channel.
- */
-class TranscendentalSolver
-{
- public:
-  TranscendentalSolver(InferenceManager& im, NlModel& m);
-  ~TranscendentalSolver();
-
-  /** init last call
-   *
-   * This is called at the beginning of last call effort check, where
-   * assertions are the set of assertions belonging to arithmetic,
-   * false_asserts is the subset of assertions that are false in the current
-   * model, and xts is the set of extended function terms that are active in
-   * the current context.
-   *
-   * This call may add lemmas to lems based on registering term
-   * information (for example, purification of sine terms).
-   */
-  void initLastCall(const std::vector<Node>& assertions,
-                    const std::vector<Node>& false_asserts,
-                    const std::vector<Node>& xts);
-  /** increment taylor degree */
-  void incrementTaylorDegree();
-  /** get taylor degree */
-  unsigned getTaylorDegree() const;
-  /** preprocess assertions check model
-   *
-   * This modifies the given assertions in preparation for running a call
-   * to check model.
-   *
-   * This method returns false if a bound for a transcendental function
-   * was conflicting.
-   */
-  bool preprocessAssertionsCheckModel(std::vector<Node>& assertions);
-  /** Process side effects in lemma se */
-  void processSideEffect(const NlLemma& se);
-  //-------------------------------------------- lemma schemas
-  /** check transcendental initial refine
-   *
-   * Constructs a set of valid theory lemmas, based on
-   * simple facts about transcendental functions.
-   * This mostly follows the initial axioms described in
-   * Section 4 of "Satisfiability
-   * Modulo Transcendental Functions via Incremental
-   * Linearization" by Cimatti et al., CADE 2017.
-   *
-   * Examples:
-   *
-   * sin( x ) = -sin( -x )
-   * ( PI > x > 0 ) => 0 < sin( x ) < 1
-   * exp( x )>0
-   * x<0 => exp( x )<1
-   */
-  void checkTranscendentalInitialRefine();
-
-  /** check transcendental monotonic
-   *
-   * Constructs a set of valid theory lemmas, based on a
-   * lemma scheme that ensures that applications
-   * of transcendental functions respect monotonicity.
-   *
-   * Examples:
-   *
-   * x > y => exp( x ) > exp( y )
-   * PI/2 > x > y > 0 => sin( x ) > sin( y )
-   * PI > x > y > PI/2 => sin( x ) < sin( y )
-   */
-  void checkTranscendentalMonotonic();
-
-  /** check transcendental tangent planes
-   *
-   * Constructs a set of valid theory lemmas, based on
-   * computing an "incremental linearization" of
-   * transcendental functions based on the model values
-   * of transcendental functions and their arguments.
-   * It is based on Figure 3 of "Satisfiability
-   * Modulo Transcendental Functions via Incremental
-   * Linearization" by Cimatti et al., CADE 2017.
-   * This schema is not terminating in general.
-   * It is not enabled by default, and can
-   * be enabled by --nl-ext-tf-tplanes.
-   *
-   * Example:
-   *
-   * Assume we have a term sin(y) where M( y ) = 1 where M is the current model.
-   * Note that:
-   *   sin(1) ~= .841471
-   *
-   * The Taylor series and remainder of sin(y) of degree 7 is
-   *   P_{7,sin(0)}( x ) = x + (-1/6)*x^3 + (1/20)*x^5
-   *   R_{7,sin(0),b}( x ) = (-1/5040)*x^7
-   *
-   * This gives us lower and upper bounds :
-   *   P_u( x ) = P_{7,sin(0)}( x ) + R_{7,sin(0),b}( x )
-   *     ...where note P_u( 1 ) = 4243/5040 ~= .841865
-   *   P_l( x ) = P_{7,sin(0)}( x ) - R_{7,sin(0),b}( x )
-   *     ...where note P_l( 1 ) = 4241/5040 ~= .841468
-   *
-   * Assume that M( sin(y) ) > P_u( 1 ).
-   * Since the concavity of sine in the region 0 < x < PI/2 is -1,
-   * we add a tangent plane refinement.
-   * The tangent plane at the point 1 in P_u is
-   * given by the formula:
-   *   T( x ) = P_u( 1 ) + ((d/dx)(P_u(x)))( 1 )*( x - 1 )
-   * We add the lemma:
-   *   ( 0 < y < PI/2 ) => sin( y ) <= T( y )
-   * which is:
-   *   ( 0 < y < PI/2 ) => sin( y ) <= (391/720)*(y - 2737/1506)
-   *
-   * Assume that M( sin(y) ) < P_u( 1 ).
-   * Since the concavity of sine in the region 0 < x < PI/2 is -1,
-   * we add a secant plane refinement for some constants ( l, u )
-   * such that 0 <= l < M( y ) < u <= PI/2. Assume we choose
-   * l = 0 and u = M( PI/2 ) = 150517/47912.
-   * The secant planes at point 1 for P_l
-   * are given by the formulas:
-   *   S_l( x ) = (x-l)*(P_l( l )-P_l(c))/(l-1) + P_l( l )
-   *   S_u( x ) = (x-u)*(P_l( u )-P_l(c))/(u-1) + P_l( u )
-   * We add the lemmas:
-   *   ( 0 < y < 1 ) => sin( y ) >= S_l( y )
-   *   ( 1 < y < PI/2 ) => sin( y ) >= S_u( y )
-   * which are:
-   *   ( 0 < y < 1 ) => (sin y) >= 4251/5040*y
-   *   ( 1 < y < PI/2 ) => (sin y) >= c1*(y+c2)
-   *     where c1, c2 are rationals (for brevity, omitted here)
-   *     such that c1 ~= .277 and c2 ~= 2.032.
-   */
-  void checkTranscendentalTangentPlanes();
- private:
-  /** check transcendental function refinement for tf
-   *
-   * This method is called by the above method for each "master"
-   * transcendental function application that occurs in an assertion in the
-   * current context. For example, an application like sin(t) is not a master
-   * if we have introduced the constraints:
-   *   t=y+2*pi*n ^ -pi <= y <= pi ^ sin(t) = sin(y).
-   * See d_trMaster/d_trSlaves for more detail.
-   *
-   * This runs Figure 3 of Cimatti et al., CADE 2017 for transcendental
-   * function application tf for Taylor degree d. It may add a secant or
-   * tangent plane lemma to lems, which includes the side effect of the lemma
-   * (if one exists).
-   *
-   * It returns false if the bounds are not precise enough to add a
-   * secant or tangent plane lemma.
-   */
-  bool checkTfTangentPlanesFun(Node tf, unsigned d);
-  //-------------------------------------------- end lemma schemas
-  /** polynomial approximation bounds
-   *
-   * This adds P_l+[x], P_l-[x], P_u+[x], P_u-[x] to pbounds, where x is
-   * d_taylor_real_fv. These are polynomial approximations of the Taylor series
-   * of <k>( 0 ) for degree 2*d where k is SINE or EXPONENTIAL.
-   * These correspond to P_l and P_u from Figure 3 of Cimatti et al., CADE 2017,
-   * for positive/negative (+/-) values of the argument of <k>( 0 ).
-   *
-   * Notice that for certain bounds (e.g. upper bounds for exponential), the
-   * Taylor approximation for a fixed degree is only sound up to a given
-   * upper bound on the argument. To obtain sound lower/upper bounds for a
-   * given <k>( c ), use the function below.
-   */
-  void getPolynomialApproximationBounds(Kind k,
-                                        unsigned d,
-                                        std::vector<Node>& pbounds);
-  /** polynomial approximation bounds
-   *
-   * This computes polynomial approximations P_l+[x], P_l-[x], P_u+[x], P_u-[x]
-   * that are sound (lower, upper) bounds for <k>( c ). Notice that these
-   * polynomials may depend on c. In particular, for P_u+[x] for <k>( c ) where
-   * c>0, we return the P_u+[x] from the function above for the minimum degree
-   * d' >= d such that (1-c^{2*d'+1}/(2*d'+1)!) is positive.
-   */
-  void getPolynomialApproximationBoundForArg(Kind k,
-                                             Node c,
-                                             unsigned d,
-                                             std::vector<Node>& pbounds);
-  /** get transcendental function model bounds
-   *
-   * This returns the current lower and upper bounds of transcendental
-   * function application tf based on Taylor of degree 2*d, which is dependent
-   * on the model value of its argument.
-   */
-  std::pair<Node, Node> getTfModelBounds(Node tf, unsigned d);
-  /** get monotonicity direction
-   *
-   * Returns whether the slope is positive (+1) or negative(-1)
-   * in region of transcendental function with kind k.
-   * Returns 0 if region is invalid.
-   */
-  int regionToMonotonicityDir(Kind k, int region);
-  /** get concavity
-   *
-   * Returns whether we are concave (+1) or convex (-1)
-   * in region of transcendental function with kind k,
-   * where region is defined above.
-   * Returns 0 if region is invalid.
-   */
-  int regionToConcavity(Kind k, int region);
-  /** region to lower bound
-   *
-   * Returns the term corresponding to the lower
-   * bound of the region of transcendental function
-   * with kind k. Returns Node::null if the region
-   * is invalid, or there is no lower bound for the
-   * region.
-   */
-  Node regionToLowerBound(Kind k, int region);
-  /** region to upper bound
-   *
-   * Returns the term corresponding to the upper
-   * bound of the region of transcendental function
-   * with kind k. Returns Node::null if the region
-   * is invalid, or there is no upper bound for the
-   * region.
-   */
-  Node regionToUpperBound(Kind k, int region);
-  /** get derivative
-   *
-   * Returns d/dx n. Supports cases of n
-   * for transcendental functions applied to x,
-   * multiplication, addition, constants and variables.
-   * Returns Node::null() if derivative is an
-   * unhandled case.
-   */
-  Node getDerivative(Node n, Node x);
-
-  void mkPi();
-  void getCurrentPiBounds();
-  /** Make the node -pi <= a <= pi */
-  static Node mkValidPhase(Node a, Node pi);
-
-  /** The inference manager that we push conflicts and lemmas to. */
-  InferenceManager& d_im;
-  /** 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_true;
-  Node d_false;
-  /**
-   * Some transcendental functions f(t) are "purified", e.g. we add
-   * t = y ^ f(t) = f(y) where y is a fresh variable. Those that are not
-   * purified we call "master terms".
-   *
-   * The maps below maintain a master/slave relationship over
-   * transcendental functions (SINE, EXPONENTIAL, PI), where above
-   * f(y) is the master of itself and of f(t).
-   *
-   * This is used for ensuring that the argument y of SINE we process is on the
-   * interval [-pi .. pi], and that exponentials are not applied to arguments
-   * that contain transcendental functions.
-   */
-  std::map<Node, Node> d_trMaster;
-  std::map<Node, std::unordered_set<Node, NodeHashFunction>> d_trSlaves;
-  /** The transcendental functions we have done initial refinements on */
-  std::map<Node, bool> d_tf_initial_refine;
-
-  /** concavity region for transcendental functions
-   *
-   * This stores an integer that identifies an interval in
-   * which the current model value for an argument of an
-   * application of a transcendental function resides.
-   *
-   * For exp( x ):
-   *   region #1 is -infty < x < infty
-   * For sin( x ):
-   *   region #0 is pi < x < infty (this is an invalid region)
-   *   region #1 is pi/2 < x <= pi
-   *   region #2 is 0 < x <= pi/2
-   *   region #3 is -pi/2 < x <= 0
-   *   region #4 is -pi < x <= -pi/2
-   *   region #5 is -infty < x <= -pi (this is an invalid region)
-   * All regions not listed above, as well as regions 0 and 5
-   * for SINE are "invalid". We only process applications
-   * of transcendental functions whose arguments have model
-   * values that reside in valid regions.
-   */
-  std::unordered_map<Node, int, NodeHashFunction> d_tf_region;
-  /** cache of the above function */
-  std::map<Kind, std::map<unsigned, std::vector<Node>>> d_poly_bounds;
-
-  /**
-   * Maps representives of a congruence class to the members of that class.
-   *
-   * In detail, a congruence class is a set of terms of the form
-   *   { f(t1), ..., f(tn) }
-   * such that t1 = ... = tn in the current context. We choose an arbitrary
-   * term among these to be the repesentative of this congruence class.
-   *
-   * Moreover, notice we compute congruence classes only over terms that
-   * are transcendental function applications that are "master terms",
-   * see d_trMaster/d_trSlave.
-   */
-  std::map<Node, std::vector<Node>> d_funcCongClass;
-  /**
-   * A list of all functions for each kind in { EXPONENTIAL, SINE, POW, PI }
-   * that are representives of their congruence class.
-   */
-  std::map<Kind, std::vector<Node>> d_funcMap;
-
-  // tangent plane bounds
-  std::map<Node, std::map<Node, Node>> d_tangent_val_bound[4];
-
-  /** secant points (sorted list) for transcendental functions
-   *
-   * This is used for tangent plane refinements for
-   * transcendental functions. This is the set
-   * "get-previous-secant-points" in "Satisfiability
-   * Modulo Transcendental Functions via Incremental
-   * Linearization" by Cimatti et al., CADE 2017, for
-   * each transcendental function application. We store this set for each
-   * Taylor degree.
-   */
-  std::unordered_map<Node,
-                     std::map<unsigned, std::vector<Node>>,
-                     NodeHashFunction>
-      d_secant_points;
-
-  /** get Taylor series of degree n for function fa centered around point fa[0].
-   *
-   * Return value is ( P_{n,f(a)}( x ), R_{n+1,f(a)}( x ) ) where
-   * the first part of the pair is the Taylor series expansion :
-   *    P_{n,f(a)}( x ) = sum_{i=0}^n (f^i( a )/i!)*(x-a)^i
-   * and the second part of the pair is the Taylor series remainder :
-   *    R_{n+1,f(a),b}( x ) = (f^{n+1}( b )/(n+1)!)*(x-a)^{n+1}
-   *
-   * The above values are cached for each (f,n) for a fixed variable "a".
-   * To compute the Taylor series for fa, we compute the Taylor series
-   *   for ( fa.getKind(), n ) then substitute { a -> fa[0] } if fa[0]!=0.
-   * We compute P_{n,f(0)}( x )/R_{n+1,f(0),b}( x ) for ( fa.getKind(), n )
-   *   if fa[0]=0.
-   * In the latter case, note we compute the exponential x^{n+1}
-   * instead of (x-a)^{n+1}, which can be done faster.
-   */
-  std::pair<Node, Node> getTaylor(Node fa, unsigned n);
-
-  /** internal variables used for constructing (cached) versions of the Taylor
-   * series above.
-   */
-  Node d_taylor_real_fv;           // x above
-  Node d_taylor_real_fv_base;      // a above
-  Node d_taylor_real_fv_base_rem;  // b above
-
-  /** cache of sum and remainder terms for getTaylor */
-  std::unordered_map<Node, std::unordered_map<unsigned, Node>, NodeHashFunction>
-      d_taylor_sum;
-  std::unordered_map<Node, std::unordered_map<unsigned, Node>, NodeHashFunction>
-      d_taylor_rem;
-  /** taylor degree
-   *
-   * Indicates that the degree of the polynomials in the Taylor approximation of
-   * all transcendental functions is 2*d_taylor_degree. This value is set
-   * initially to options::nlExtTfTaylorDegree() and may be incremented
-   * if the option options::nlExtTfIncPrecision() is enabled.
-   */
-  unsigned d_taylor_degree;
-  /** PI
-   *
-   * Note that PI is a (symbolic, non-constant) nullary operator. This is
-   * because its value cannot be computed exactly. We constraint PI to concrete
-   * lower and upper bounds stored in d_pi_bound below.
-   */
-  Node d_pi;
-  /** PI/2 */
-  Node d_pi_2;
-  /** -PI/2 */
-  Node d_pi_neg_2;
-  /** -PI */
-  Node d_pi_neg;
-  /** the concrete lower and upper bounds for PI */
-  Node d_pi_bound[2];
-}; /* class TranscendentalSolver */
-
-}  // namespace nl
-}  // namespace arith
-}  // namespace theory
-}  // namespace CVC4
-
-#endif /* CVC4__THEORY__ARITH__TRANSCENDENTAL_SOLVER_H */