Add remaining transcendental functions (#1551)

11 files changed, 501 insertions, 204 deletions

diff --git a/src/parser/smt2/smt2.cpp b/src/parser/smt2/smt2.cpp
index e4f6569b8..77b50af4c 100644
--- a/src/parser/smt2/smt2.cpp
+++ b/src/parser/smt2/smt2.cpp
@@ -61,6 +61,17 @@ void Smt2::addArithmeticOperators() {
   addOperator(kind::SINE, "sin");
   addOperator(kind::COSINE, "cos");
   addOperator(kind::TANGENT, "tan");
+  addOperator(kind::COSECANT, "csc");
+  addOperator(kind::SECANT, "sec");
+  addOperator(kind::COTANGENT, "cot");
+  addOperator(kind::ARCSINE, "arcsin");
+  addOperator(kind::ARCCOSINE, "arccos");
+  addOperator(kind::ARCTANGENT, "arctan");
+  addOperator(kind::ARCCOSECANT, "arccsc");
+  addOperator(kind::ARCSECANT, "arcsec");
+  addOperator(kind::ARCCOTANGENT, "arccot");
+
+  addOperator(kind::SQRT, "sqrt");
 }
 
 void Smt2::addBitvectorOperators() {
diff --git a/src/printer/smt2/smt2_printer.cpp b/src/printer/smt2/smt2_printer.cpp
index 24fd0924f..e06f8c062 100644
--- a/src/printer/smt2/smt2_printer.cpp
+++ b/src/printer/smt2/smt2_printer.cpp
@@ -425,6 +425,16 @@ void Smt2Printer::toStream(std::ostream& out,
   case kind::SINE:
   case kind::COSINE:
   case kind::TANGENT:
+  case kind::COSECANT:
+  case kind::SECANT:
+  case kind::COTANGENT:
+  case kind::ARCSINE:
+  case kind::ARCCOSINE:
+  case kind::ARCTANGENT:
+  case kind::ARCCOSECANT:
+  case kind::ARCSECANT:
+  case kind::ARCCOTANGENT:
+  case kind::SQRT:
   case kind::MINUS:
   case kind::UMINUS:
   case kind::LT:
@@ -891,6 +901,16 @@ static string smtKindString(Kind k)
   case kind::SINE: return "sin";
   case kind::COSINE: return "cos";
   case kind::TANGENT: return "tan";
+  case kind::COSECANT: return "csc";
+  case kind::SECANT: return "sec";
+  case kind::COTANGENT: return "cot";
+  case kind::ARCSINE: return "arcsin";
+  case kind::ARCCOSINE: return "arccos";
+  case kind::ARCTANGENT: return "arctan";
+  case kind::ARCCOSECANT: return "arccsc";
+  case kind::ARCSECANT: return "arcsec";
+  case kind::ARCCOTANGENT: return "arccot";
+  case kind::SQRT: return "sqrt";
   case kind::MINUS: return "-";
   case kind::UMINUS: return "-";
   case kind::LT: return "<";
diff --git a/src/theory/arith/arith_rewriter.cpp b/src/theory/arith/arith_rewriter.cpp
index b862e7604..a9761ade4 100644
--- a/src/theory/arith/arith_rewriter.cpp
+++ b/src/theory/arith/arith_rewriter.cpp
@@ -107,7 +107,15 @@ RewriteResponse ArithRewriter::preRewriteTerm(TNode t){
     case kind::SINE:
     case kind::COSINE:
     case kind::TANGENT:
-      return preRewriteTranscendental(t);
+    case kind::COSECANT:
+    case kind::SECANT:
+    case kind::COTANGENT:
+    case kind::ARCSINE:
+    case kind::ARCCOSINE:
+    case kind::ARCTANGENT:
+    case kind::ARCCOSECANT:
+    case kind::ARCSECANT:
+    case kind::ARCCOTANGENT: return preRewriteTranscendental(t);
     case kind::INTS_DIVISION:
     case kind::INTS_MODULUS:
       return RewriteResponse(REWRITE_DONE, t);
@@ -163,7 +171,15 @@ RewriteResponse ArithRewriter::postRewriteTerm(TNode t){
     case kind::SINE:
     case kind::COSINE:
     case kind::TANGENT:
-      return postRewriteTranscendental(t);
+    case kind::COSECANT:
+    case kind::SECANT:
+    case kind::COTANGENT:
+    case kind::ARCSINE:
+    case kind::ARCCOSINE:
+    case kind::ARCTANGENT:
+    case kind::ARCCOSECANT:
+    case kind::ARCSECANT:
+    case kind::ARCCOTANGENT: return postRewriteTranscendental(t);
     case kind::INTS_DIVISION:
     case kind::INTS_MODULUS:
       return RewriteResponse(REWRITE_DONE, t);
@@ -360,28 +376,30 @@ RewriteResponse ArithRewriter::preRewriteTranscendental(TNode t) {
 
 RewriteResponse ArithRewriter::postRewriteTranscendental(TNode t) { 
   Trace("arith-tf-rewrite") << "Rewrite transcendental function : " << t << std::endl;
+  NodeManager* nm = NodeManager::currentNM();
   switch( t.getKind() ){
   case kind::EXPONENTIAL: {
     if(t[0].getKind() == kind::CONST_RATIONAL){
-      Node one = NodeManager::currentNM()->mkConst(Rational(1));
+      Node one = nm->mkConst(Rational(1));
       if(t[0].getConst<Rational>().sgn()>=0 && t[0].getType().isInteger() && t[0]!=one){
-        return RewriteResponse(REWRITE_AGAIN, NodeManager::currentNM()->mkNode(kind::POW, NodeManager::currentNM()->mkNode( kind::EXPONENTIAL, one ), t[0]));
+        return RewriteResponse(
+            REWRITE_AGAIN,
+            nm->mkNode(kind::POW, nm->mkNode(kind::EXPONENTIAL, one), t[0]));
       }else{          
         return RewriteResponse(REWRITE_DONE, t);
       }
     }else if(t[0].getKind() == kind::PLUS ){
       std::vector<Node> product;
       for( unsigned i=0; i<t[0].getNumChildren(); i++ ){
-        product.push_back( NodeManager::currentNM()->mkNode( kind::EXPONENTIAL, t[0][i] ) );
+        product.push_back(nm->mkNode(kind::EXPONENTIAL, t[0][i]));
       }
-      return RewriteResponse(REWRITE_AGAIN, NodeManager::currentNM()->mkNode(kind::MULT, product));
+      return RewriteResponse(REWRITE_AGAIN, nm->mkNode(kind::MULT, product));
     }
   }
     break;
   case kind::SINE:
     if(t[0].getKind() == kind::CONST_RATIONAL){
       const Rational& rat = t[0].getConst<Rational>();
-      NodeManager* nm = NodeManager::currentNM();
       if(rat.sgn() == 0){
         return RewriteResponse(REWRITE_DONE, nm->mkConst(Rational(0)));
       }
@@ -433,26 +451,29 @@ RewriteResponse ArithRewriter::postRewriteTranscendental(TNode t) {
         if (r_abs > rone)
         {
           //add/substract 2*pi beyond scope
-          Node ra_div_two = NodeManager::currentNM()->mkNode(
+          Node ra_div_two = nm->mkNode(
               kind::INTS_DIVISION, mkRationalNode(r_abs + rone), ntwo);
           Node new_pi_factor;
           if( r.sgn()==1 ){
-            new_pi_factor = NodeManager::currentNM()->mkNode( kind::MINUS, pi_factor, NodeManager::currentNM()->mkNode( kind::MULT, ntwo, ra_div_two ) );
+            new_pi_factor =
+                nm->mkNode(kind::MINUS,
+                           pi_factor,
+                           nm->mkNode(kind::MULT, ntwo, ra_div_two));
           }else{
             Assert( r.sgn()==-1 );
-            new_pi_factor = NodeManager::currentNM()->mkNode( kind::PLUS, pi_factor, NodeManager::currentNM()->mkNode( kind::MULT, ntwo, ra_div_two ) );
+            new_pi_factor =
+                nm->mkNode(kind::PLUS,
+                           pi_factor,
+                           nm->mkNode(kind::MULT, ntwo, ra_div_two));
           }
-          Node new_arg =
-              NodeManager::currentNM()->mkNode(kind::MULT, new_pi_factor, pi);
+          Node new_arg = nm->mkNode(kind::MULT, new_pi_factor, pi);
           if (!rem.isNull())
           {
-            new_arg =
-                NodeManager::currentNM()->mkNode(kind::PLUS, new_arg, rem);
+            new_arg = nm->mkNode(kind::PLUS, new_arg, rem);
           }
           // sin( 2*n*PI + x ) = sin( x )
-          return RewriteResponse(
-              REWRITE_AGAIN_FULL,
-              NodeManager::currentNM()->mkNode(kind::SINE, new_arg));
+          return RewriteResponse(REWRITE_AGAIN_FULL,
+                                 nm->mkNode(kind::SINE, new_arg));
         }
         else if (r_abs == rone)
         {
@@ -465,9 +486,7 @@ RewriteResponse ArithRewriter::postRewriteTranscendental(TNode t) {
           {
             return RewriteResponse(
                 REWRITE_AGAIN_FULL,
-                NodeManager::currentNM()->mkNode(
-                    kind::UMINUS,
-                    NodeManager::currentNM()->mkNode(kind::SINE, rem)));
+                nm->mkNode(kind::UMINUS, nm->mkNode(kind::SINE, rem)));
           }
         }
         else if (rem.isNull())
@@ -498,16 +517,48 @@ RewriteResponse ArithRewriter::postRewriteTranscendental(TNode t) {
     }
     break;
   case kind::COSINE: {
-    return RewriteResponse(REWRITE_AGAIN_FULL, NodeManager::currentNM()->mkNode( kind::SINE, 
-                                                 NodeManager::currentNM()->mkNode( kind::MINUS, 
-                                                   NodeManager::currentNM()->mkNode( kind::MULT, 
-                                                     NodeManager::currentNM()->mkConst( Rational(1)/Rational(2) ),
-                                                     NodeManager::currentNM()->mkNullaryOperator( NodeManager::currentNM()->realType(), kind::PI ) ),
-                                                     t[0] ) ) );
-  } break;
+    return RewriteResponse(
+        REWRITE_AGAIN_FULL,
+        nm->mkNode(kind::SINE,
+                   nm->mkNode(kind::MINUS,
+                              nm->mkNode(kind::MULT,
+                                         nm->mkConst(Rational(1) / Rational(2)),
+                                         mkPi()),
+                              t[0])));
+  }
+  break;
   case kind::TANGENT:
-    return RewriteResponse(REWRITE_AGAIN_FULL, NodeManager::currentNM()->mkNode(kind::DIVISION, NodeManager::currentNM()->mkNode( kind::SINE, t[0] ), 
-                                                                                           NodeManager::currentNM()->mkNode( kind::COSINE, t[0] ) ));
+  {
+    return RewriteResponse(REWRITE_AGAIN_FULL,
+                           nm->mkNode(kind::DIVISION,
+                                      nm->mkNode(kind::SINE, t[0]),
+                                      nm->mkNode(kind::COSINE, t[0])));
+  }
+  break;
+  case kind::COSECANT:
+  {
+    return RewriteResponse(REWRITE_AGAIN_FULL,
+                           nm->mkNode(kind::DIVISION,
+                                      mkRationalNode(Rational(1)),
+                                      nm->mkNode(kind::SINE, t[0])));
+  }
+  break;
+  case kind::SECANT:
+  {
+    return RewriteResponse(REWRITE_AGAIN_FULL,
+                           nm->mkNode(kind::DIVISION,
+                                      mkRationalNode(Rational(1)),
+                                      nm->mkNode(kind::COSINE, t[0])));
+  }
+  break;
+  case kind::COTANGENT:
+  {
+    return RewriteResponse(REWRITE_AGAIN_FULL,
+                           nm->mkNode(kind::DIVISION,
+                                      nm->mkNode(kind::COSINE, t[0]),
+                                      nm->mkNode(kind::SINE, t[0])));
+  }
+  break;
   default:
     break;
   }
diff --git a/src/theory/arith/kinds b/src/theory/arith/kinds
index 34ae30f4c..3073d0078 100644
--- a/src/theory/arith/kinds
+++ b/src/theory/arith/kinds
@@ -31,6 +31,17 @@ operator EXPONENTIAL 1 "exponential"
 operator SINE 1 "sine"
 operator COSINE 1 "consine"
 operator TANGENT 1 "tangent"
+operator COSECANT 1 "cosecant"
+operator SECANT 1 "secant"
+operator COTANGENT 1 "cotangent"
+operator ARCSINE 1 "arc sine"
+operator ARCCOSINE 1 "arc consine"
+operator ARCTANGENT 1 "arc tangent"
+operator ARCCOSECANT 1 "arc cosecant"
+operator ARCSECANT 1 "arc secant"
+operator ARCCOTANGENT 1 "arc cotangent"
+
+operator SQRT 1 "square root"
 
 constant DIVISIBLE_OP \
         ::CVC4::Divisible \
@@ -105,6 +116,17 @@ typerule EXPONENTIAL ::CVC4::theory::arith::RealOperatorTypeRule
 typerule SINE ::CVC4::theory::arith::RealOperatorTypeRule
 typerule COSINE ::CVC4::theory::arith::RealOperatorTypeRule
 typerule TANGENT ::CVC4::theory::arith::RealOperatorTypeRule
+typerule COSECANT ::CVC4::theory::arith::RealOperatorTypeRule
+typerule SECANT ::CVC4::theory::arith::RealOperatorTypeRule
+typerule COTANGENT ::CVC4::theory::arith::RealOperatorTypeRule
+typerule ARCSINE ::CVC4::theory::arith::RealOperatorTypeRule
+typerule ARCCOSINE ::CVC4::theory::arith::RealOperatorTypeRule
+typerule ARCTANGENT ::CVC4::theory::arith::RealOperatorTypeRule
+typerule ARCCOSECANT ::CVC4::theory::arith::RealOperatorTypeRule
+typerule ARCSECANT ::CVC4::theory::arith::RealOperatorTypeRule
+typerule ARCCOTANGENT ::CVC4::theory::arith::RealOperatorTypeRule
+
+typerule SQRT ::CVC4::theory::arith::RealOperatorTypeRule
 
 nullaryoperator PI "pi"
 
diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp
index 65a7597f1..1522708d9 100644
--- a/src/theory/arith/nonlinear_extension.cpp
+++ b/src/theory/arith/nonlinear_extension.cpp
@@ -216,12 +216,14 @@ bool hasNewMonomials(Node n, const std::vector<Node>& existing) {
 
 NonlinearExtension::NonlinearExtension(TheoryArith& containing,
                                        eq::EqualityEngine* ee)
-    : d_lemmas(containing.getUserContext()),
+    : d_def_lemmas(containing.getUserContext()),
+      d_lemmas(containing.getUserContext()),
       d_zero_split(containing.getUserContext()),
       d_skolem_atoms(containing.getUserContext()),
       d_containing(containing),
       d_ee(ee),
-      d_needsLastCall(false) {
+      d_needsLastCall(false)
+{
   d_true = NodeManager::currentNM()->mkConst(true);
   d_false = NodeManager::currentNM()->mkConst(false);
   d_zero = NodeManager::currentNM()->mkConst(Rational(0));
@@ -1032,7 +1034,9 @@ Kind NonlinearExtension::transKinds(Kind k1, Kind k2) {
 }
 
 bool NonlinearExtension::isTranscendentalKind(Kind k) {
-  Assert(k != TANGENT && k != COSINE);  // eliminated
+  // many operators are eliminated during rewriting
+  Assert(k != TANGENT && k != COSINE && k != COSECANT && k != SECANT
+         && k != COTANGENT);
   return k == EXPONENTIAL || k == SINE || k == PI;
 }
  
@@ -1161,10 +1165,12 @@ bool NonlinearExtension::checkModelTf(const std::vector<Node>& assertions)
 
   if (check_assertions.empty())
   {
+    Trace("nl-ext-tf-check-model") << "...simple check succeeded." << std::endl;
     return true;
   }
   else
   {
+    Trace("nl-ext-tf-check-model") << "...simple check failed." << std::endl;
     // TODO (#1450) check model for general case
     return false;
   }
@@ -1250,7 +1256,6 @@ bool NonlinearExtension::simpleCheckModelTfLit(Node lit)
       return comp == d_true;
     }
   }
-
   Trace("nl-ext-tf-check-model-simple") << "  failed due to unknown literal."
                                         << std::endl;
   return false;
@@ -1288,11 +1293,13 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
   d_tf_check_model_bounds.clear();
 
   int lemmas_proc = 0;
-  std::vector<Node> lemmas;  
-  
+  std::vector<Node> lemmas;
+  NodeManager* nm = NodeManager::currentNM();
+
   Trace("nl-ext-mv") << "Extended terms : " << std::endl;
   // register the extended function terms
   std::map< Node, Node > mvarg_to_term;
+  std::vector<Node> trig_no_base;
   for( unsigned i=0; i<xts.size(); i++ ){
     Node a = xts[i];
     computeModelValue(a, 0);
@@ -1341,38 +1348,11 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
       {
         if( d_trig_is_base.find( a )==d_trig_is_base.end() ){
           consider = false;
-          if( d_trig_base.find( a )==d_trig_base.end() ){
-            Node y = NodeManager::currentNM()->mkSkolem("y",NodeManager::currentNM()->realType(),"phase shifted trigonometric arg");
-            Node new_a = NodeManager::currentNM()->mkNode( a.getKind(), y );
-            d_trig_is_base[new_a] = true;
-            d_trig_base[a] = new_a;
-            Trace("nl-ext-tf") << "Basis sine : " << new_a << " for " << a << std::endl;
-            if( d_pi.isNull() ){
-              mkPi();
-              getCurrentPiBounds( lemmas );
-            }
-            Node shift = NodeManager::currentNM()->mkSkolem( "s", NodeManager::currentNM()->integerType(), "number of shifts" );
-            // FIXME : do not introduce shift here, instead needs model-based
-            // refinement for constant shifts (#1284)
-            Node shift_lem = NodeManager::currentNM()->mkNode(
-                AND,
-                mkValidPhase(y, d_pi),
-                a[0].eqNode(NodeManager::currentNM()->mkNode(
-                    PLUS,
-                    y,
-                    NodeManager::currentNM()->mkNode(
-                        MULT,
-                        NodeManager::currentNM()->mkConst(Rational(2)),
-                        shift,
-                        d_pi))),
-                // particular case of above for shift=0
-                NodeManager::currentNM()->mkNode(
-                    IMPLIES, mkValidPhase(a[0], d_pi), a[0].eqNode(y)),
-                new_a.eqNode(a));
-            //must do preprocess on this one
-            Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : shift : " << shift_lem << std::endl;
-            d_containing.getOutputChannel().lemma(shift_lem, false, true);
-            lemmas_proc++;
+          trig_no_base.push_back(a);
+          if (d_pi.isNull())
+          {
+            mkPi();
+            getCurrentPiBounds(lemmas);
           }
         }
       }
@@ -1383,7 +1363,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
           //verify they have the same model value
           if( d_mv[1][a]!=d_mv[1][itrm->second] ){
             // if not, add congruence lemma
-            Node cong_lemma = NodeManager::currentNM()->mkNode(
+            Node cong_lemma = nm->mkNode(
                 IMPLIES, a[0].eqNode(itrm->second[0]), a.eqNode(itrm->second));
             lemmas.push_back( cong_lemma );
             //Assert( false );
@@ -1407,6 +1387,45 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
     return lemmas_proc;
   }
 
+  // process SINE phase shifting
+  for (const Node& a : trig_no_base)
+  {
+    if (d_trig_base.find(a) == d_trig_base.end())
+    {
+      Node y =
+          nm->mkSkolem("y", nm->realType(), "phase shifted trigonometric arg");
+      Node new_a = nm->mkNode(a.getKind(), y);
+      d_trig_is_base[new_a] = true;
+      d_trig_base[a] = new_a;
+      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");
+      // FIXME : do not introduce shift here, instead needs model-based
+      // refinement for constant shifts (#1284)
+      Node shift_lem = nm->mkNode(
+          AND,
+          mkValidPhase(y, d_pi),
+          a[0].eqNode(nm->mkNode(
+              PLUS,
+              y,
+              nm->mkNode(MULT, nm->mkConst(Rational(2)), shift, d_pi))),
+          // particular case of above for shift=0
+          nm->mkNode(IMPLIES, mkValidPhase(a[0], d_pi), a[0].eqNode(y)),
+          new_a.eqNode(a));
+      // must do preprocess on this one
+      Trace("nl-ext-lemma")
+          << "NonlinearExtension::Lemma : shift : " << shift_lem << std::endl;
+      d_containing.getOutputChannel().lemma(shift_lem, false, true);
+      lemmas_proc++;
+    }
+  }
+  if (lemmas_proc > 0)
+  {
+    Trace("nl-ext") << "  ...finished with " << lemmas_proc
+                    << " new lemmas SINE phase shifting." << std::endl;
+    return lemmas_proc;
+  }
 
   // register constants
   registerMonomial(d_one);
@@ -1742,6 +1761,24 @@ void NonlinearExtension::check(Theory::Effort e) {
   }
 }
 
+void NonlinearExtension::addDefinition(Node lem)
+{
+  Trace("nl-ext") << "NonlinearExtension::addDefinition : " << lem << std::endl;
+  d_def_lemmas.insert(lem);
+}
+
+void NonlinearExtension::presolve()
+{
+  Trace("nl-ext") << "NonlinearExtension::presolve, #defs = "
+                  << d_def_lemmas.size() << std::endl;
+  for (NodeSet::const_iterator it = d_def_lemmas.begin();
+       it != d_def_lemmas.end();
+       ++it)
+  {
+    flushLemma(*it);
+  }
+}
+
 void NonlinearExtension::assignOrderIds(std::vector<Node>& vars,
                                         NodeMultiset& order,
                                         unsigned orderType) {
@@ -3487,6 +3524,10 @@ std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes()
                   antec.size() == 1 ? antec[0] : nm->mkNode(AND, antec);
               lem = nm->mkNode(IMPLIES, antec_n, lem);
             }
+            Trace("nl-ext-tf-tplanes-debug")
+                << "*** Tangent plane lemma (pre-rewrite): " << lem
+                << std::endl;
+            lem = Rewriter::rewrite(lem);
             Trace("nl-ext-tf-tplanes") << "*** Tangent plane lemma : " << lem
                                        << std::endl;
             // Figure 3 : line 9
@@ -3607,6 +3648,10 @@ std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes()
                                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-tf-tplanes-debug")
+                    << "*** Secant plane lemma (pre-rewrite) : " << lem
+                    << std::endl;
+                lem = Rewriter::rewrite(lem);
                 Trace("nl-ext-tf-tplanes") << "*** Secant plane lemma : " << lem
                                            << std::endl;
                 // Figure 3 : line 22
diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h
index 34da28f6c..84acc0269 100644
--- a/src/theory/arith/nonlinear_extension.h
+++ b/src/theory/arith/nonlinear_extension.h
@@ -121,6 +121,23 @@ class NonlinearExtension {
   void check(Theory::Effort e);
   /** Does this class need a call to check(...) at last call effort? */
   bool needsCheckLastEffort() const { return d_needsLastCall; }
+  /** add definition
+   *
+   * This function notifies this class that lem is a formula that defines or
+   * constrains an auxiliary variable. For example, during
+   * TheoryArith::expandDefinitions, we replace a term like arcsin( x ) with an
+   * auxiliary variable k. The lemmas 0 <= k < pi and sin( x ) = k are added as
+   * definitions to this class.
+   */
+  void addDefinition(Node lem);
+  /** presolve
+   *
+   * This function is called during TheoryArith's presolve command.
+   * In this function, we send lemmas we accumulated during preprocessing,
+   * for instance, definitional lemmas from expandDefinitions are sent out
+   * on the output channel of TheoryArith in this function.
+   */
+  void presolve();
   /** Compare arithmetic terms i and j based an ordering.
    *
    * orderType = 0 : compare concrete model values
@@ -387,8 +404,11 @@ class NonlinearExtension {
   // ( x*y, x*z, y ) for each pair of monomials ( x*y, x*z ) with common factors
   std::map<Node, std::map<Node, Node> > d_mono_diff;
 
-  // cache of all lemmas sent
+  /** cache of definition lemmas (user-context-dependent) */
+  NodeSet d_def_lemmas;
+  /** cache of all lemmas sent on the output channel (user-context-dependent) */
   NodeSet d_lemmas;
+  /** cache of terms t for which we have added the lemma ( t = 0 V t != 0 ). */
   NodeSet d_zero_split;
   
   // literals with Skolems (need not be satisfied by model)
diff --git a/src/theory/arith/normal_form.cpp b/src/theory/arith/normal_form.cpp
index 30b9ca0b5..76782d8a5 100644
--- a/src/theory/arith/normal_form.cpp
+++ b/src/theory/arith/normal_form.cpp
@@ -95,7 +95,16 @@ bool Variable::isTranscendentalMember(Node n) {
   case kind::SINE:
   case kind::COSINE:
   case kind::TANGENT:
-    return Polynomial::isMember(n[0]);
+  case kind::COSECANT:
+  case kind::SECANT:
+  case kind::COTANGENT:
+  case kind::ARCSINE:
+  case kind::ARCCOSINE:
+  case kind::ARCTANGENT:
+  case kind::ARCCOSECANT:
+  case kind::ARCSECANT:
+  case kind::ARCCOTANGENT:
+  case kind::SQRT: return Polynomial::isMember(n[0]);
   case kind::PI:
     return true;
   default:
diff --git a/src/theory/arith/normal_form.h b/src/theory/arith/normal_form.h
index 21301da91..ba740146b 100644
--- a/src/theory/arith/normal_form.h
+++ b/src/theory/arith/normal_form.h
@@ -245,6 +245,16 @@ public:
     case kind::SINE:
     case kind::COSINE:
     case kind::TANGENT:
+    case kind::COSECANT:
+    case kind::SECANT:
+    case kind::COTANGENT:
+    case kind::ARCSINE:
+    case kind::ARCCOSINE:
+    case kind::ARCTANGENT:
+    case kind::ARCCOSECANT:
+    case kind::ARCSECANT:
+    case kind::ARCCOTANGENT:
+    case kind::SQRT:
     case kind::PI:
       return isTranscendentalMember(n);      
     case kind::ABS:
diff --git a/src/theory/arith/theory_arith.cpp b/src/theory/arith/theory_arith.cpp
index e354305d7..1390cbee6 100644
--- a/src/theory/arith/theory_arith.cpp
+++ b/src/theory/arith/theory_arith.cpp
@@ -42,8 +42,6 @@ TheoryArith::TheoryArith(context::Context* c, context::UserContext* u,
     getExtTheory()->addFunctionKind(kind::NONLINEAR_MULT);
     getExtTheory()->addFunctionKind(kind::EXPONENTIAL);
     getExtTheory()->addFunctionKind(kind::SINE);
-    getExtTheory()->addFunctionKind(kind::COSINE);
-    getExtTheory()->addFunctionKind(kind::TANGENT);
     getExtTheory()->addFunctionKind(kind::PI);
   }
 }
diff --git a/src/theory/arith/theory_arith_private.cpp b/src/theory/arith/theory_arith_private.cpp
index 6a4456844..fc0673d21 100644
--- a/src/theory/arith/theory_arith_private.cpp
+++ b/src/theory/arith/theory_arith_private.cpp
@@ -85,67 +85,84 @@ namespace arith {
 static Node toSumNode(const ArithVariables& vars, const DenseMap<Rational>& sum);
 static bool complexityBelow(const DenseMap<Rational>& row, uint32_t cap);
 
-
-TheoryArithPrivate::TheoryArithPrivate(TheoryArith& containing, context::Context* c, context::UserContext* u, OutputChannel& out, Valuation valuation, const LogicInfo& logicInfo) :
-  d_containing(containing),
-  d_nlIncomplete( false),
-  d_rowTracking(),
-  d_constraintDatabase(c, u, d_partialModel, d_congruenceManager, RaiseConflict(*this)),
-  d_qflraStatus(Result::SAT_UNKNOWN),
-  d_unknownsInARow(0),
-  d_hasDoneWorkSinceCut(false),
-  d_learner(u),
-  d_assertionsThatDoNotMatchTheirLiterals(c),
-  d_nextIntegerCheckVar(0),
-  d_constantIntegerVariables(c),
-  d_diseqQueue(c, false),
-  d_currentPropagationList(),
-  d_learnedBounds(c),
-  d_partialModel(c, DeltaComputeCallback(*this)),
-  d_errorSet(d_partialModel, TableauSizes(&d_tableau), BoundCountingLookup(*this)),
-  d_tableau(),
-  d_linEq(d_partialModel, d_tableau, d_rowTracking, BasicVarModelUpdateCallBack(*this)),
-  d_diosolver(c),
-  d_restartsCounter(0),
-  d_tableauSizeHasBeenModified(false),
-  d_tableauResetDensity(1.6),
-  d_tableauResetPeriod(10),
-  d_conflicts(c),
-  d_blackBoxConflict(c, Node::null()),
-  d_congruenceManager(c, d_constraintDatabase, SetupLiteralCallBack(*this), d_partialModel, RaiseEqualityEngineConflict(*this)),
-  d_cmEnabled(c, true),
-
-  d_dualSimplex(d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
-  d_fcSimplex(d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
-  d_soiSimplex(d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
-  d_attemptSolSimplex(d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
-  d_nonlinearExtension( NULL ),
-  d_pass1SDP(NULL),
-  d_otherSDP(NULL),
-  d_lastContextIntegerAttempted(c,-1),
-
-
-  d_DELTA_ZERO(0),
-  d_approxCuts(c),
-  d_fullCheckCounter(0),
-  d_cutCount(c, 0),
-  d_cutInContext(c),
-  d_likelyIntegerInfeasible(c, false),
-  d_guessedCoeffSet(c, false),
-  d_guessedCoeffs(),
-  d_treeLog(NULL),
-  d_replayVariables(),
-  d_replayConstraints(),
-  d_lhsTmp(),
-  d_approxStats(NULL),
-  d_attemptSolveIntTurnedOff(u, 0),
-  d_dioSolveResources(0),
-  d_solveIntMaybeHelp(0u),
-  d_solveIntAttempts(0u),
-  d_statistics(),
-  d_to_int_skolem(u),
-  d_div_skolem(u),
-  d_int_div_skolem(u)
+TheoryArithPrivate::TheoryArithPrivate(TheoryArith& containing,
+                                       context::Context* c,
+                                       context::UserContext* u,
+                                       OutputChannel& out,
+                                       Valuation valuation,
+                                       const LogicInfo& logicInfo)
+    : d_containing(containing),
+      d_nlIncomplete(false),
+      d_rowTracking(),
+      d_constraintDatabase(
+          c, u, d_partialModel, d_congruenceManager, RaiseConflict(*this)),
+      d_qflraStatus(Result::SAT_UNKNOWN),
+      d_unknownsInARow(0),
+      d_hasDoneWorkSinceCut(false),
+      d_learner(u),
+      d_assertionsThatDoNotMatchTheirLiterals(c),
+      d_nextIntegerCheckVar(0),
+      d_constantIntegerVariables(c),
+      d_diseqQueue(c, false),
+      d_currentPropagationList(),
+      d_learnedBounds(c),
+      d_partialModel(c, DeltaComputeCallback(*this)),
+      d_errorSet(
+          d_partialModel, TableauSizes(&d_tableau), BoundCountingLookup(*this)),
+      d_tableau(),
+      d_linEq(d_partialModel,
+              d_tableau,
+              d_rowTracking,
+              BasicVarModelUpdateCallBack(*this)),
+      d_diosolver(c),
+      d_restartsCounter(0),
+      d_tableauSizeHasBeenModified(false),
+      d_tableauResetDensity(1.6),
+      d_tableauResetPeriod(10),
+      d_conflicts(c),
+      d_blackBoxConflict(c, Node::null()),
+      d_congruenceManager(c,
+                          d_constraintDatabase,
+                          SetupLiteralCallBack(*this),
+                          d_partialModel,
+                          RaiseEqualityEngineConflict(*this)),
+      d_cmEnabled(c, true),
+
+      d_dualSimplex(
+          d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
+      d_fcSimplex(
+          d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
+      d_soiSimplex(
+          d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
+      d_attemptSolSimplex(
+          d_linEq, d_errorSet, RaiseConflict(*this), TempVarMalloc(*this)),
+      d_nonlinearExtension(NULL),
+      d_pass1SDP(NULL),
+      d_otherSDP(NULL),
+      d_lastContextIntegerAttempted(c, -1),
+
+      d_DELTA_ZERO(0),
+      d_approxCuts(c),
+      d_fullCheckCounter(0),
+      d_cutCount(c, 0),
+      d_cutInContext(c),
+      d_likelyIntegerInfeasible(c, false),
+      d_guessedCoeffSet(c, false),
+      d_guessedCoeffs(),
+      d_treeLog(NULL),
+      d_replayVariables(),
+      d_replayConstraints(),
+      d_lhsTmp(),
+      d_approxStats(NULL),
+      d_attemptSolveIntTurnedOff(u, 0),
+      d_dioSolveResources(0),
+      d_solveIntMaybeHelp(0u),
+      d_solveIntAttempts(0u),
+      d_statistics(),
+      d_to_int_skolem(u),
+      d_div_skolem(u),
+      d_int_div_skolem(u),
+      d_nlin_inverse_skolem(u)
 {
   if( options::nlExt() ){
     d_nonlinearExtension = new NonlinearExtension(
@@ -4419,6 +4436,11 @@ void TheoryArithPrivate::presolve(){
     Debug("arith::oldprop") << " lemma lemma duck " <<lem << endl;
     outputLemma(lem);
   }
+
+  if (options::nlExt())
+  {
+    d_nonlinearExtension->presolve();
+  }
 }
 
 EqualityStatus TheoryArithPrivate::getEqualityStatus(TNode a, TNode b) {
@@ -4863,90 +4885,178 @@ const BoundsInfo& TheoryArithPrivate::boundsInfo(ArithVar basic) const{
 Node TheoryArithPrivate::expandDefinition(LogicRequest &logicRequest, Node node) {
   NodeManager* nm = NodeManager::currentNM();
 
-  // eliminate here since involves division
-  if( node.getKind()==kind::TANGENT ){
-    node = nm->mkNode(kind::DIVISION, nm->mkNode( kind::SINE, node[0] ), 
-                                      nm->mkNode( kind::COSINE, node[0] ) );
+  // eliminate here since the rewritten form of these may introduce division
+  Kind k = node.getKind();
+  if (k == kind::TANGENT || k == kind::COSECANT || k == kind::SECANT
+      || k == kind::COTANGENT)
+  {
+    node = Rewriter::rewrite(node);
+    k = node.getKind();
   }
 
-  switch(node.getKind()) {
-  case kind::DIVISION: {
-    TNode num = node[0], den = node[1];
-    Node ret = nm->mkNode(kind::DIVISION_TOTAL, num, den);
-    if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
+  switch (k)
+  {
+    case kind::DIVISION:
     {
-      // partial function: division
-      if (d_divByZero.isNull())
+      TNode num = node[0], den = node[1];
+      Node ret = nm->mkNode(kind::DIVISION_TOTAL, num, den);
+      if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
       {
-        d_divByZero =
-            nm->mkSkolem("divByZero",
-                         nm->mkFunctionType(nm->realType(), nm->realType()),
-                         "partial real division",
-                         NodeManager::SKOLEM_EXACT_NAME);
-        logicRequest.widenLogic(THEORY_UF);
+        // partial function: division
+        if (d_divByZero.isNull())
+        {
+          d_divByZero =
+              nm->mkSkolem("divByZero",
+                           nm->mkFunctionType(nm->realType(), nm->realType()),
+                           "partial real division",
+                           NodeManager::SKOLEM_EXACT_NAME);
+          logicRequest.widenLogic(THEORY_UF);
+        }
+        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
+        Node divByZeroNum = nm->mkNode(kind::APPLY_UF, d_divByZero, num);
+        ret = nm->mkNode(kind::ITE, denEq0, divByZeroNum, ret);
       }
-      Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
-      Node divByZeroNum = nm->mkNode(kind::APPLY_UF, d_divByZero, num);
-      ret = nm->mkNode(kind::ITE, denEq0, divByZeroNum, ret);
+      return ret;
+      break;
     }
-    return ret;
-    break;
-  }
 
-  case kind::INTS_DIVISION: {
-    // partial function: integer div
-    TNode num = node[0], den = node[1];
-    Node ret = nm->mkNode(kind::INTS_DIVISION_TOTAL, num, den);
-    if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
+    case kind::INTS_DIVISION:
     {
-      if (d_intDivByZero.isNull())
+      // partial function: integer div
+      TNode num = node[0], den = node[1];
+      Node ret = nm->mkNode(kind::INTS_DIVISION_TOTAL, num, den);
+      if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
       {
-        d_intDivByZero = nm->mkSkolem(
-            "intDivByZero",
-            nm->mkFunctionType(nm->integerType(), nm->integerType()),
-            "partial integer division",
-            NodeManager::SKOLEM_EXACT_NAME);
-        logicRequest.widenLogic(THEORY_UF);
+        if (d_intDivByZero.isNull())
+        {
+          d_intDivByZero = nm->mkSkolem(
+              "intDivByZero",
+              nm->mkFunctionType(nm->integerType(), nm->integerType()),
+              "partial integer division",
+              NodeManager::SKOLEM_EXACT_NAME);
+          logicRequest.widenLogic(THEORY_UF);
+        }
+        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
+        Node intDivByZeroNum = nm->mkNode(kind::APPLY_UF, d_intDivByZero, num);
+        ret = nm->mkNode(kind::ITE, denEq0, intDivByZeroNum, ret);
       }
-      Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
-      Node intDivByZeroNum = nm->mkNode(kind::APPLY_UF, d_intDivByZero, num);
-      ret = nm->mkNode(kind::ITE, denEq0, intDivByZeroNum, ret);
+      return ret;
+      break;
     }
-    return ret;
-    break;
-  }
 
-  case kind::INTS_MODULUS: {
-    // partial function: mod
-    TNode num = node[0], den = node[1];
-    Node ret = nm->mkNode(kind::INTS_MODULUS_TOTAL, num, den);
-    if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
+    case kind::INTS_MODULUS:
     {
-      if (d_modZero.isNull())
+      // partial function: mod
+      TNode num = node[0], den = node[1];
+      Node ret = nm->mkNode(kind::INTS_MODULUS_TOTAL, num, den);
+      if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
       {
-        d_modZero = nm->mkSkolem(
-            "modZero",
-            nm->mkFunctionType(nm->integerType(), nm->integerType()),
-            "partial modulus",
-            NodeManager::SKOLEM_EXACT_NAME);
-        logicRequest.widenLogic(THEORY_UF);
+        if (d_modZero.isNull())
+        {
+          d_modZero = nm->mkSkolem(
+              "modZero",
+              nm->mkFunctionType(nm->integerType(), nm->integerType()),
+              "partial modulus",
+              NodeManager::SKOLEM_EXACT_NAME);
+          logicRequest.widenLogic(THEORY_UF);
+        }
+        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
+        Node modZeroNum = nm->mkNode(kind::APPLY_UF, d_modZero, num);
+        ret = nm->mkNode(kind::ITE, denEq0, modZeroNum, ret);
       }
-      Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
-      Node modZeroNum = nm->mkNode(kind::APPLY_UF, d_modZero, num);
-      ret = nm->mkNode(kind::ITE, denEq0, modZeroNum, ret);
+      return ret;
+      break;
     }
-    return ret;
-    break;
-  }
 
-  case kind::ABS: {
-    return nm->mkNode(kind::ITE, nm->mkNode(kind::LT, node[0], nm->mkConst(Rational(0))), nm->mkNode(kind::UMINUS, node[0]), node[0]);
-    break;
-  }
+    case kind::ABS:
+    {
+      return nm->mkNode(kind::ITE,
+                        nm->mkNode(kind::LT, node[0], nm->mkConst(Rational(0))),
+                        nm->mkNode(kind::UMINUS, node[0]),
+                        node[0]);
+      break;
+    }
+    case kind::SQRT:
+    case kind::ARCSINE:
+    case kind::ARCCOSINE:
+    case kind::ARCTANGENT:
+    case kind::ARCCOSECANT:
+    case kind::ARCSECANT:
+    case kind::ARCCOTANGENT:
+    {
+      // eliminate inverse functions here
+      NodeMap::const_iterator it = d_nlin_inverse_skolem.find(node);
+      if (it == d_nlin_inverse_skolem.end())
+      {
+        Node var = nm->mkSkolem("nonlinearInv",
+                                nm->realType(),
+                                "the result of a non-linear inverse function");
+        d_nlin_inverse_skolem[node] = var;
+        Node lem;
+        if (k == kind::SQRT)
+        {
+          lem = nm->mkNode(kind::MULT, node[0], node[0]).eqNode(var);
+        }
+        else
+        {
+          Node pi = mkPi();
+
+          // range of the skolem
+          Node rlem;
+          if (k == kind::ARCSINE || k == ARCTANGENT || k == ARCCOSECANT)
+          {
+            Node half = nm->mkConst(Rational(1) / Rational(2));
+            Node pi2 = nm->mkNode(kind::MULT, half, pi);
+            Node npi2 = nm->mkNode(kind::MULT, nm->mkConst(Rational(-1)), pi2);
+            // -pi/2 < var <= pi/2
+            rlem = nm->mkNode(
+                AND, nm->mkNode(LT, npi2, var), nm->mkNode(LEQ, var, pi2));
+          }
+          else
+          {
+            // 0 <= var < pi
+            rlem = nm->mkNode(AND,
+                              nm->mkNode(LEQ, nm->mkConst(Rational(0)), var),
+                              nm->mkNode(LT, var, pi));
+          }
+          if (options::nlExt())
+          {
+            d_nonlinearExtension->addDefinition(rlem);
+          }
 
-  default:
-    return node;
-    break;
+          Kind rk = k == kind::ARCSINE
+                        ? kind::SINE
+                        : (k == kind::ARCCOSINE
+                               ? kind::COSINE
+                               : (k == kind::ARCTANGENT
+                                      ? kind::TANGENT
+                                      : (k == kind::ARCCOSECANT
+                                             ? kind::COSECANT
+                                             : (k == kind::ARCSECANT
+                                                    ? kind::SECANT
+                                                    : kind::COTANGENT))));
+          Node invTerm = nm->mkNode(rk, var);
+          // since invTerm may introduce division,
+          // we must also call expandDefinition on the result
+          invTerm = expandDefinition(logicRequest, invTerm);
+          lem = invTerm.eqNode(node[0]);
+        }
+        Assert(!lem.isNull());
+        if (options::nlExt())
+        {
+          d_nonlinearExtension->addDefinition(lem);
+        }
+        else
+        {
+          d_nlIncomplete = true;
+        }
+        return var;
+      }
+      return (*it).second;
+      break;
+    }
+
+    default: return node; break;
   }
 
   Unreachable();
diff --git a/src/theory/arith/theory_arith_private.h b/src/theory/arith/theory_arith_private.h
index 459bb41d9..23712016d 100644
--- a/src/theory/arith/theory_arith_private.h
+++ b/src/theory/arith/theory_arith_private.h
@@ -848,16 +848,17 @@ private:
    * semantics.  Needed to deal with partial function "mod".
    */
   Node d_modZero;
-  
-  /** 
-   *  Maps for Skolems for to-integer, real/integer div-by-k.
-   *  Introduced during ppRewriteTerms.
+
+  /**
+   *  Maps for Skolems for to-integer, real/integer div-by-k, and inverse
+   *  non-linear operators that are introduced during ppRewriteTerms.
    */
   typedef context::CDHashMap< Node, Node, NodeHashFunction > NodeMap;
   NodeMap d_to_int_skolem;
   NodeMap d_div_skolem;
   NodeMap d_int_div_skolem;
-  
+  NodeMap d_nlin_inverse_skolem;
+
 };/* class TheoryArithPrivate */
 
 }/* CVC4::theory::arith namespace */