Refactor operator elimination in arithmetic (#4519)

This is a major refactor of how operators are eliminated in arithmetic. Currently there were (at least) two things wrong:
(1) ppRewriteTerm sent lemmas on the output channel. This behavior is incompatible with how preprocessing works. In particular, this caused unconstrained simplification to be unaware of terms from such lemmas, leading to incorrect "sat" answers.
(2) Lemmas used to eliminate certain "div-like" terms were processed in a context-independent way. However, lemmas should be cached in a user-context-dependent way. This was leading to incorrect "sat" answers in incremental.

The solution to these issues is to eliminate operators via the construction of witness terms. No lemmas are sent out, and instead these lemmas are the consequence of term formula removal in the standard way.

As a result of the refactor, 2 quantifiers regressions time out due to infinite branch and bound issues (one only during --check-unsat-cores). These appear to be random and I've changed the options to avoid these issues. 3 others now have check-model warnings, which I've added --quiet to. Improving check-model will be addressed on a future PR.

This PR is not required for SMT COMP since we have workarounds that avoid both the incorrect behaviors in our scripts.

Also notice that --rewrite-divk is effectively now enabled by default always.

Fixes #4484, fixes #4486, fixes #4481.

4 files changed, 464 insertions, 421 deletions

diff --git a/src/options/arith_options.toml b/src/options/arith_options.toml
index 4865ccead..d16513cfb 100644
--- a/src/options/arith_options.toml
+++ b/src/options/arith_options.toml
@@ -276,14 +276,6 @@ header = "options/arith_options.h"
   help       = "use quick explain to minimize the sum of infeasibility conflicts"
 
 [[option]]
-  name       = "rewriteDivk"
-  category   = "regular"
-  long       = "rewrite-divk"
-  type       = "bool"
-  default    = "false"
-  help       = "rewrite division and mod when by a constant into linear terms"
-
-[[option]]
   name       = "trySolveIntStandardEffort"
   category   = "regular"
   long       = "se-solve-int"
diff --git a/src/smt/set_defaults.cpp b/src/smt/set_defaults.cpp
index 899da4d4a..e06363883 100644
--- a/src/smt/set_defaults.cpp
+++ b/src/smt/set_defaults.cpp
@@ -1097,11 +1097,6 @@ void setDefaults(SmtEngine& smte, LogicInfo& logic)
     {
       options::quantSplit.set(false);
     }
-    // rewrite divk
-    if (!options::rewriteDivk.wasSetByUser())
-    {
-      options::rewriteDivk.set(true);
-    }
     // do not do macros
     if (!options::macrosQuant.wasSetByUser())
     {
@@ -1136,11 +1131,6 @@ void setDefaults(SmtEngine& smte, LogicInfo& logic)
   }
   if (options::cegqi())
   {
-    // must rewrite divk
-    if (!options::rewriteDivk.wasSetByUser())
-    {
-      options::rewriteDivk.set(true);
-    }
     if (options::incrementalSolving())
     {
       // cannot do nested quantifier elimination in incremental mode
diff --git a/src/theory/arith/theory_arith_private.cpp b/src/theory/arith/theory_arith_private.cpp
index 7fdab2034..09b6d742a 100644
--- a/src/theory/arith/theory_arith_private.cpp
+++ b/src/theory/arith/theory_arith_private.cpp
@@ -34,6 +34,7 @@
 #include "expr/node.h"
 #include "expr/node_algorithm.h"
 #include "expr/node_builder.h"
+#include "expr/proof_skolem_cache.h"
 #include "options/arith_options.h"
 #include "options/smt_options.h"  // for incrementalSolving()
 #include "preprocessing/util/ite_utilities.h"
@@ -226,100 +227,6 @@ static void resolve(ConstraintCPVec& buf, ConstraintP c, const ConstraintCPVec&
   // return safeConstructNary(nb);
 }
 
-// Integer division axioms:
-// These concenrate the high level constructs needed to constrain the functions:
-// div, mod, div_total and mod_total.
-//
-// The high level constraint.
-// (for all ((m Int) (n Int))
-//   (=> (distinct n 0)
-//       (let ((q (div m n)) (r (mod m n)))
-//         (and (= m (+ (* n q) r))
-//              (<= 0 r (- (abs n) 1))))))
-//
-// We now add division by 0 functions.
-// (for all ((m Int) (n Int))
-//   (let ((q (div m n)) (r (mod m n)))
-//     (ite (= n 0)
-//          (and (= q (div_0_func m)) (= r (mod_0_func m)))
-//          (and (= m (+ (* n q) r))
-//               (<= 0 r (- (abs n) 1)))))))
-//
-// We now express div and mod in terms of div_total and mod_total.
-// (for all ((m Int) (n Int))
-//   (let ((q (div m n)) (r (mod m n)))
-//     (ite (= n 0)
-//          (and (= q (div_0_func m)) (= r (mod_0_func m)))
-//          (and (= q (div_total m n)) (= r (mod_total m n))))))
-//
-// Alternative div_total and mod_total without the abs function:
-// (for all ((m Int) (n Int))
-//   (let ((q (div_total m n)) (r (mod_total m n)))
-//     (ite (= n 0)
-//          (and (= q 0) (= r 0))
-//          (and (r = m - n * q)
-//               (ite (> n 0)
-//                    (n*q <= m < n*q + n)
-//                    (n*q <= m < n*q - n))))))
-
-
-// Returns a formula that entails that q equals the total integer division of m
-// by n.
-// (for all ((m Int) (n Int))
-//   (let ((q (div_total m n)))
-//     (ite (= n 0)
-//          (= q 0)
-//          (ite (> n 0)
-//               (n*q <= m < n*q + n)
-//               (n*q <= m < n*q - n)))))
-Node mkAxiomForTotalIntDivision(Node m, Node n, Node q) {
-  NodeManager* nm = NodeManager::currentNM();
-  Node zero = mkRationalNode(0);
-  // (n*q <= m < n*q + n) is equivalent to (0 <= m - n*q < n).
-  Node diff = nm->mkNode(kind::MINUS, m, mkMult(n, q));
-  Node when_n_is_positive = mkInRange(diff, zero, n);
-  Node when_n_is_negative = mkInRange(diff, zero, nm->mkNode(kind::UMINUS, n));
-  return n.eqNode(zero).iteNode(
-      q.eqNode(zero), nm->mkNode(kind::GT, n, zero)
-                          .iteNode(when_n_is_positive, when_n_is_negative));
-}
-
-// Returns a formula that entails that r equals the integer division total of m
-// by n assuming q is equal to (div_total m n).
-// (for all ((m Int) (n Int))
-//   (let ((q (div_total m n)) (r (mod_total m n)))
-//     (ite (= n 0)
-//          (= r 0)
-//          (= r (- m (n * q))))))
-Node mkAxiomForTotalIntMod(Node m, Node n, Node q, Node r) {
-  NodeManager* nm = NodeManager::currentNM();
-  Node diff = nm->mkNode(kind::MINUS, m, mkMult(n, q));
-  return mkOnZeroIte(n, r, mkRationalNode(0), diff);
-}
-
-// Returns an expression that constrains a term to be the result of the
-// [total] integer division of num by den. Equivalently:
-// (and (=> (den > 0) (den*term <= num < den*term + den))
-//      (=> (den < 0) (den*term <= num < den*term - den))
-//      (=> (den = 0) (term = 0)))
-// static Node mkIntegerDivisionConstraint(Node term, Node num, Node den) {
-//   // We actually encode:
-//   // (and (=> (den > 0) (0 <= num - den*term < den))
-//   //      (=> (den < 0) (0 <= num - den*term < -den))
-//   //      (=> (den = 0) (term = 0)))
-//   NodeManager* nm = NodeManager::currentNM();
-//   Node zero = nm->mkConst(Rational(0));
-//   Node den_is_positive = nm->mkNode(kind::GT, den, zero);
-//   Node den_is_negative = nm->mkNode(kind::LT, den, zero);
-//   Node diff = nm->mkNode(kind::MINUS, num, mkMult(den, term));
-//   Node when_den_positive = den_positive.impNode(mkInRange(diff, zero, den));
-//   Node when_den_negative = den_negative.impNode(
-//       mkInRange(diff, zero, nm->mkNode(kind::UMINUS, den)));
-//   Node when_den_is_zero = (den.eqNode(zero)).impNode(term.eqNode(zero));
-//   return mk->mkNode(kind::AND, when_den_positive, when_den_negative,
-//                     when_den_is_zero);
-// }
-
 void TheoryArithPrivate::setMasterEqualityEngine(eq::EqualityEngine* eq) {
   d_congruenceManager.setMasterEqualityEngine(eq);
 }
@@ -1176,126 +1083,443 @@ Node TheoryArithPrivate::ppRewriteTerms(TNode n) {
   if(Theory::theoryOf(n) != THEORY_ARITH) {
     return n;
   }
+  // Eliminate operators recursively. Notice we must do this here since other
+  // theories may generate lemmas that involve non-standard operators. For
+  // example, quantifier instantiation may use TO_INTEGER terms; SyGuS may
+  // introduce non-standard arithmetic terms appearing in grammars.
+  // call eliminate operators
+  Node nn = eliminateOperators(n);
+  if (nn != n)
+  {
+    // since elimination may introduce new operators to eliminate, we must
+    // recursively eliminate result
+    return eliminateOperatorsRec(nn);
+  }
+  return n;
+}
+
+void TheoryArithPrivate::checkNonLinearLogic(Node term)
+{
+  if (getLogicInfo().isLinear())
+  {
+    Trace("arith-logic") << "ERROR: Non-linear term in linear logic: " << term
+                         << std::endl;
+    std::stringstream serr;
+    serr << "A non-linear fact was asserted to arithmetic in a linear logic."
+         << std::endl;
+    serr << "The fact in question: " << term << endl;
+    throw LogicException(serr.str());
+  }
+}
+
+struct ArithElimOpAttributeId
+{
+};
+typedef expr::Attribute<ArithElimOpAttributeId, Node> ArithElimOpAttribute;
+
+Node TheoryArithPrivate::eliminateOperatorsRec(Node n)
+{
+  ArithElimOpAttribute aeoa;
+  Trace("arith-elim") << "Begin elim: " << n << std::endl;
+  NodeManager* nm = NodeManager::currentNM();
+  std::unordered_map<Node, Node, TNodeHashFunction> visited;
+  std::unordered_map<Node, Node, TNodeHashFunction>::iterator it;
+  std::vector<Node> visit;
+  Node cur;
+  visit.push_back(n);
+  do
+  {
+    cur = visit.back();
+    visit.pop_back();
+    it = visited.find(cur);
+    if (Theory::theoryOf(cur) != THEORY_ARITH)
+    {
+      visited[cur] = cur;
+    }
+    else if (it == visited.end())
+    {
+      if (cur.hasAttribute(aeoa))
+      {
+        visited[cur] = cur.getAttribute(aeoa);
+      }
+      else
+      {
+        visited[cur] = Node::null();
+        visit.push_back(cur);
+        for (const Node& cn : cur)
+        {
+          visit.push_back(cn);
+        }
+      }
+    }
+    else if (it->second.isNull())
+    {
+      Node ret = cur;
+      bool childChanged = false;
+      std::vector<Node> children;
+      if (cur.getMetaKind() == metakind::PARAMETERIZED)
+      {
+        children.push_back(cur.getOperator());
+      }
+      for (const Node& cn : cur)
+      {
+        it = visited.find(cn);
+        Assert(it != visited.end());
+        Assert(!it->second.isNull());
+        childChanged = childChanged || cn != it->second;
+        children.push_back(it->second);
+      }
+      if (childChanged)
+      {
+        ret = nm->mkNode(cur.getKind(), children);
+      }
+      Node retElim = eliminateOperators(ret);
+      if (retElim != ret)
+      {
+        // recursively eliminate operators in result, since some eliminations
+        // are defined in terms of other non-standard operators.
+        ret = eliminateOperatorsRec(retElim);
+      }
+      cur.setAttribute(aeoa, ret);
+      visited[cur] = ret;
+    }
+  } while (!visit.empty());
+  Assert(visited.find(n) != visited.end());
+  Assert(!visited.find(n)->second.isNull());
+  return visited[n];
+}
 
+Node TheoryArithPrivate::eliminateOperators(Node node)
+{
   NodeManager* nm = NodeManager::currentNM();
 
-  switch(Kind k = n.getKind()) {
-
-  case kind::TO_INTEGER:
-  case kind::IS_INTEGER: {
-    Node toIntSkolem;
-    NodeMap::const_iterator it = d_to_int_skolem.find( n[0] );
-    if( it==d_to_int_skolem.end() ){
-      toIntSkolem = nm->mkSkolem("toInt", nm->integerType(),
-                            "a conversion of a Real term to its Integer part");
-      d_to_int_skolem[n[0]] = toIntSkolem;
-      // n[0] - 1 < toIntSkolem <= n[0]
-      // -1 < toIntSkolem - n[0] <= 0
-      // 0 <= n[0] - toIntSkolem < 1
-      Node one = mkRationalNode(1);
-      Node lem = mkAxiomForTotalIntDivision(n[0], one, toIntSkolem);
-      outputLemma(lem);
-    }else{
-      toIntSkolem = (*it).second;
+  Kind k = node.getKind();
+  switch (k)
+  {
+    case kind::TANGENT:
+    case kind::COSECANT:
+    case kind::SECANT:
+    case kind::COTANGENT:
+    {
+      // these are eliminated by rewriting
+      return Rewriter::rewrite(node);
+      break;
     }
-    if(k == kind::IS_INTEGER) {
-      return nm->mkNode(kind::EQUAL, n[0], toIntSkolem);
+    case kind::TO_INTEGER:
+    case kind::IS_INTEGER:
+    {
+      Node toIntSkolem;
+      NodeMap::const_iterator it = d_to_int_skolem.find(node[0]);
+      if (it == d_to_int_skolem.end())
+      {
+        // node[0] - 1 < toIntSkolem <= node[0]
+        // -1 < toIntSkolem - node[0] <= 0
+        // 0 <= node[0] - toIntSkolem < 1
+        Node v = nm->mkBoundVar(nm->integerType());
+        Node one = mkRationalNode(1);
+        Node zero = mkRationalNode(0);
+        Node diff = nm->mkNode(kind::MINUS, node[0], v);
+        Node lem = mkInRange(diff, zero, one);
+        toIntSkolem = ProofSkolemCache::mkSkolem(
+            v, lem, "toInt", "a conversion of a Real term to its Integer part");
+        toIntSkolem = ProofSkolemCache::getWitnessForm(toIntSkolem);
+        d_to_int_skolem[node[0]] = toIntSkolem;
+      }
+      else
+      {
+        toIntSkolem = (*it).second;
+      }
+      if (k == kind::IS_INTEGER)
+      {
+        return nm->mkNode(kind::EQUAL, node[0], toIntSkolem);
+      }
+      Assert(k == kind::TO_INTEGER);
+      return toIntSkolem;
     }
-    Assert(k == kind::TO_INTEGER);
-    return toIntSkolem;
-  }
-  case kind::INTS_DIVISION:
-  case kind::INTS_MODULUS:
-  case kind::DIVISION:
-    // these should be removed during expand definitions
-    Assert(false);
-    break;
-  
-  case kind::INTS_DIVISION_TOTAL: 
-  case kind::INTS_MODULUS_TOTAL: {
-    Node den = Rewriter::rewrite(n[1]);
-    if(!options::rewriteDivk() && den.isConst()) {
-      return n;
-    }
-    Node num = Rewriter::rewrite(n[0]);
-    Node intVar;
-    Node rw = nm->mkNode(k, num, den);
-    NodeMap::const_iterator it = d_int_div_skolem.find( rw );
-    if( it==d_int_div_skolem.end() ){
-      intVar = nm->mkSkolem("linearIntDiv", nm->integerType(), "the result of an intdiv-by-k term");
-      d_int_div_skolem[rw] = intVar;
-      Node lem;
-      if (den.isConst()) {
-        const Rational& rat = den.getConst<Rational>();
-        Assert(!num.isConst());
-        if(rat != 0) {
-          if(rat > 0) {
-            lem = nm->mkNode(kind::AND, nm->mkNode(kind::LEQ, nm->mkNode(kind::MULT, den, intVar), num), 
-                                        nm->mkNode(kind::LT, num, nm->mkNode(kind::MULT, den, nm->mkNode(kind::PLUS, intVar, nm->mkConst(Rational(1))))));
-          } else {
-            lem = nm->mkNode(kind::AND, nm->mkNode(kind::LEQ, nm->mkNode(kind::MULT, den, intVar), num), 
-                                        nm->mkNode(kind::LT, num, nm->mkNode(kind::MULT, den, nm->mkNode(kind::PLUS, intVar, nm->mkConst(Rational(-1))))));
+
+    case kind::INTS_DIVISION_TOTAL:
+    case kind::INTS_MODULUS_TOTAL:
+    {
+      Node den = Rewriter::rewrite(node[1]);
+      Node num = Rewriter::rewrite(node[0]);
+      Node intVar;
+      Node rw = nm->mkNode(k, num, den);
+      NodeMap::const_iterator it = d_int_div_skolem.find(rw);
+      if (it == d_int_div_skolem.end())
+      {
+        Node v = nm->mkBoundVar(nm->integerType());
+        Node lem;
+        Node leqNum = nm->mkNode(LEQ, nm->mkNode(MULT, den, v), num);
+        if (den.isConst())
+        {
+          const Rational& rat = den.getConst<Rational>();
+          if (num.isConst() || rat == 0)
+          {
+            // just rewrite
+            return Rewriter::rewrite(node);
+          }
+          if (rat > 0)
+          {
+            lem = nm->mkNode(
+                AND,
+                leqNum,
+                nm->mkNode(
+                    LT,
+                    num,
+                    nm->mkNode(MULT,
+                               den,
+                               nm->mkNode(PLUS, v, nm->mkConst(Rational(1))))));
+          }
+          else
+          {
+            lem = nm->mkNode(
+                AND,
+                leqNum,
+                nm->mkNode(
+                    LT,
+                    num,
+                    nm->mkNode(
+                        MULT,
+                        den,
+                        nm->mkNode(PLUS, v, nm->mkConst(Rational(-1))))));
           }
         }
-      }else{
-        lem = nm->mkNode(kind::AND,
-                nm->mkNode(kind::IMPLIES, NodeManager::currentNM()->mkNode( kind::GT, den, nm->mkConst(Rational(0)) ),
-                  nm->mkNode(kind::AND, nm->mkNode(kind::LEQ, nm->mkNode(kind::MULT, den, intVar), num), 
-                                        nm->mkNode(kind::LT, num, nm->mkNode(kind::MULT, den, nm->mkNode(kind::PLUS, intVar, nm->mkConst(Rational(1))))))),
-                nm->mkNode(kind::IMPLIES, NodeManager::currentNM()->mkNode( kind::LT, den, nm->mkConst(Rational(0)) ),
-                  nm->mkNode(kind::AND, nm->mkNode(kind::LEQ, nm->mkNode(kind::MULT, den, intVar), num), 
-                                        nm->mkNode(kind::LT, num, nm->mkNode(kind::MULT, den, nm->mkNode(kind::PLUS, intVar, nm->mkConst(Rational(-1))))))));                
-      }    
-      if( !lem.isNull() ){
-        outputLemma(lem);
+        else
+        {
+          checkNonLinearLogic(node);
+          lem = nm->mkNode(
+              AND,
+              nm->mkNode(
+                  IMPLIES,
+                  nm->mkNode(GT, den, nm->mkConst(Rational(0))),
+                  nm->mkNode(
+                      AND,
+                      leqNum,
+                      nm->mkNode(
+                          LT,
+                          num,
+                          nm->mkNode(
+                              MULT,
+                              den,
+                              nm->mkNode(PLUS, v, nm->mkConst(Rational(1))))))),
+              nm->mkNode(
+                  IMPLIES,
+                  nm->mkNode(LT, den, nm->mkConst(Rational(0))),
+                  nm->mkNode(
+                      AND,
+                      leqNum,
+                      nm->mkNode(
+                          LT,
+                          num,
+                          nm->mkNode(
+                              MULT,
+                              den,
+                              nm->mkNode(
+                                  PLUS, v, nm->mkConst(Rational(-1))))))));
+        }
+        intVar = ProofSkolemCache::mkSkolem(
+            v, lem, "linearIntDiv", "the result of an intdiv-by-k term");
+        intVar = ProofSkolemCache::getWitnessForm(intVar);
+        d_int_div_skolem[rw] = intVar;
       }
-    }else{
-      intVar = (*it).second;
+      else
+      {
+        intVar = (*it).second;
+      }
+      if (k == kind::INTS_MODULUS_TOTAL)
+      {
+        Node nn =
+            nm->mkNode(kind::MINUS, num, nm->mkNode(kind::MULT, den, intVar));
+        return nn;
+      }
+      else
+      {
+        return intVar;
+      }
+      break;
     }
-    if( k==kind::INTS_MODULUS_TOTAL ){
-      Node node = nm->mkNode(kind::MINUS, num, nm->mkNode(kind::MULT, den, intVar));
-      return node;
-    }else{
-      return intVar;
+    case kind::DIVISION_TOTAL:
+    {
+      Node num = Rewriter::rewrite(node[0]);
+      Node den = Rewriter::rewrite(node[1]);
+      if (den.isConst())
+      {
+        // No need to eliminate here, can eliminate via rewriting later.
+        // Moreover, rewriting may change the type of this node from real to
+        // int, which impacts certain issues with subtyping.
+        return node;
+      }
+      checkNonLinearLogic(node);
+      Node var;
+      Node rw = nm->mkNode(k, num, den);
+      NodeMap::const_iterator it = d_div_skolem.find(rw);
+      if (it == d_div_skolem.end())
+      {
+        Node v = nm->mkBoundVar(nm->realType());
+        Node lem = nm->mkNode(IMPLIES,
+                              den.eqNode(nm->mkConst(Rational(0))).negate(),
+                              nm->mkNode(MULT, den, v).eqNode(num));
+        var = ProofSkolemCache::mkSkolem(
+            v, lem, "nonlinearDiv", "the result of a non-linear div term");
+        var = ProofSkolemCache::getWitnessForm(var);
+        d_div_skolem[rw] = var;
+      }
+      else
+      {
+        var = (*it).second;
+      }
+      return var;
+      break;
     }
-    break;
-  }
-  case kind::DIVISION_TOTAL: {
-    Node num = Rewriter::rewrite(n[0]);
-    Node den = Rewriter::rewrite(n[1]);
-    Assert(!den.isConst());
-    Node var;
-    Node rw = nm->mkNode(k, num, den);
-    NodeMap::const_iterator it = d_div_skolem.find( rw );
-    if( it==d_div_skolem.end() ){
-      var = nm->mkSkolem("nonlinearDiv", nm->realType(), "the result of a non-linear div term");
-      d_div_skolem[rw] = var;
-      outputLemma(nm->mkNode(kind::IMPLIES, den.eqNode(nm->mkConst(Rational(0))).negate(), nm->mkNode(kind::MULT, den, var).eqNode(num)));
-    }else{
-      var = (*it).second;
+    case kind::DIVISION:
+    {
+      Node num = Rewriter::rewrite(node[0]);
+      Node den = Rewriter::rewrite(node[1]);
+      Node ret = nm->mkNode(kind::DIVISION_TOTAL, num, den);
+      if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
+      {
+        checkNonLinearLogic(node);
+        Node divByZeroNum = getArithSkolemApp(num, ArithSkolemId::DIV_BY_ZERO);
+        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
+        ret = nm->mkNode(kind::ITE, denEq0, divByZeroNum, ret);
+      }
+      return ret;
+      break;
     }
-    return var;
-    break;
-  }
-  default:
-    break;
-  }
 
-  for(TNode::const_iterator i = n.begin(); i != n.end(); ++i) {
-    Node rewritten = ppRewriteTerms(*i);
-    if(rewritten != *i) {
-      NodeBuilder<> b(n.getKind());
-      b.append(n.begin(), i);
-      b << rewritten;
-      for(++i; i != n.end(); ++i) {
-        b << ppRewriteTerms(*i);
+    case kind::INTS_DIVISION:
+    {
+      // partial function: integer div
+      Node num = Rewriter::rewrite(node[0]);
+      Node den = Rewriter::rewrite(node[1]);
+      Node ret = nm->mkNode(kind::INTS_DIVISION_TOTAL, num, den);
+      if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
+      {
+        checkNonLinearLogic(node);
+        Node intDivByZeroNum =
+            getArithSkolemApp(num, ArithSkolemId::INT_DIV_BY_ZERO);
+        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
+        ret = nm->mkNode(kind::ITE, denEq0, intDivByZeroNum, ret);
       }
-      rewritten = b;
-      return rewritten;
+      return ret;
+      break;
     }
-  }
 
-  return n;
+    case kind::INTS_MODULUS:
+    {
+      // partial function: mod
+      Node num = Rewriter::rewrite(node[0]);
+      Node den = Rewriter::rewrite(node[1]);
+      Node ret = nm->mkNode(kind::INTS_MODULUS_TOTAL, num, den);
+      if (!den.isConst() || den.getConst<Rational>().sgn() == 0)
+      {
+        checkNonLinearLogic(node);
+        Node modZeroNum = getArithSkolemApp(num, ArithSkolemId::MOD_BY_ZERO);
+        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
+        ret = nm->mkNode(kind::ITE, denEq0, modZeroNum, ret);
+      }
+      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::SQRT:
+    case kind::ARCSINE:
+    case kind::ARCCOSINE:
+    case kind::ARCTANGENT:
+    case kind::ARCCOSECANT:
+    case kind::ARCSECANT:
+    case kind::ARCCOTANGENT:
+    {
+      checkNonLinearLogic(node);
+      // eliminate inverse functions here
+      NodeMap::const_iterator it = d_nlin_inverse_skolem.find(node);
+      if (it == d_nlin_inverse_skolem.end())
+      {
+        Node var = nm->mkBoundVar(nm->realType());
+        Node lem;
+        if (k == kind::SQRT)
+        {
+          Node skolemApp = getArithSkolemApp(node[0], ArithSkolemId::SQRT);
+          Node uf = skolemApp.eqNode(var);
+          Node nonNeg = nm->mkNode(
+              kind::AND, nm->mkNode(kind::MULT, var, var).eqNode(node[0]), uf);
+
+          // sqrt(x) reduces to:
+          // witness y. ite(x >= 0.0, y * y = x ^ y = Uf(x), y = Uf(x))
+          //
+          // Uf(x) makes sure that the reduction still behaves like a function,
+          // otherwise the reduction of (x = 1) ^ (sqrt(x) != sqrt(1)) would be
+          // satisfiable. On the original formula, this would require that we
+          // simultaneously interpret sqrt(1) as 1 and -1, which is not a valid
+          // model.
+          lem = nm->mkNode(
+              kind::ITE,
+              nm->mkNode(kind::GEQ, node[0], nm->mkConst(Rational(0))),
+              nonNeg,
+              uf);
+        }
+        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));
+          }
+
+          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);
+          lem = nm->mkNode(AND, rlem, invTerm.eqNode(node[0]));
+        }
+        Assert(!lem.isNull());
+        Node ret = ProofSkolemCache::mkSkolem(
+            var,
+            lem,
+            "tfk",
+            "Skolem to eliminate a non-standard transcendental function");
+        ret = ProofSkolemCache::getWitnessForm(ret);
+        d_nlin_inverse_skolem[node] = ret;
+        return ret;
+      }
+      return (*it).second;
+      break;
+    }
+
+    default: break;
+  }
+  return node;
 }
 
 Node TheoryArithPrivate::ppRewrite(TNode atom) {
@@ -1312,6 +1536,7 @@ Node TheoryArithPrivate::ppRewrite(TNode atom) {
       Node rewritten = Rewriter::rewrite(leq.andNode(geq));
       Debug("arith::preprocess")
           << "arith::preprocess() : returning " << rewritten << endl;
+      // don't need to rewrite terms since rewritten is not a non-standard op
       return rewritten;
     }
   }
@@ -1441,11 +1666,6 @@ void TheoryArithPrivate::setupVariable(const Variable& x){
   //setupInitialValue(varN);
 
   markSetup(n);
-
-  if(x.isDivLike()){
-    setupDivLike(x);
-  }
-
 }
 
 void TheoryArithPrivate::setupVariableList(const VarList& vl){
@@ -1508,52 +1728,6 @@ void TheoryArithPrivate::cautiousSetupPolynomial(const Polynomial& p){
   }
 }
 
-void TheoryArithPrivate::setupDivLike(const Variable& v){
-  Assert(v.isDivLike());
-
-  if(getLogicInfo().isLinear()){
-    throw LogicException(
-        "A non-linear fact (involving div/mod/divisibility) was asserted to "
-        "arithmetic in a linear logic;\nif you only use division (or modulus) "
-        "by a constant value, or if you only use the divisibility-by-k "
-        "predicate, try using the --rewrite-divk option.");
-  }
-
-  Node vnode = v.getNode();
-  Assert(
-      isSetup(vnode));  // Otherwise there is some invariant breaking recursion
-  Polynomial m = Polynomial::parsePolynomial(vnode[0]);
-  Polynomial n = Polynomial::parsePolynomial(vnode[1]);
-
-  cautiousSetupPolynomial(m);
-  cautiousSetupPolynomial(n);
-
-  Node lem;
-  switch(vnode.getKind()){
-  case DIVISION:
-  case INTS_DIVISION:
-  case INTS_MODULUS:
-    // these should be removed during expand definitions
-    Assert(false);
-    break;
-  case DIVISION_TOTAL:
-    lem = axiomIteForTotalDivision(vnode);
-    break;
-  case INTS_DIVISION_TOTAL:
-  case INTS_MODULUS_TOTAL:
-    lem = axiomIteForTotalIntDivision(vnode);
-    break;
-  default:
-    /* intentionally blank */
-    break;
-  }
-
-  if(!lem.isNull()){
-    Debug("arith::div") << lem << endl;
-    outputLemma(lem);
-  }
-}
-
 Node TheoryArithPrivate::axiomIteForTotalDivision(Node div_tot){
   Assert(div_tot.getKind() == DIVISION_TOTAL);
 
@@ -5056,157 +5230,15 @@ const BoundsInfo& TheoryArithPrivate::boundsInfo(ArithVar basic) const{
 
 Node TheoryArithPrivate::expandDefinition(Node node)
 {
-  NodeManager* nm = NodeManager::currentNM();
-
-  // 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)
+  // call eliminate operators
+  Node nn = eliminateOperators(node);
+  if (nn != node)
   {
-    node = Rewriter::rewrite(node);
-    k = node.getKind();
+    // since elimination may introduce new operators to eliminate, we must
+    // recursively eliminate result
+    return eliminateOperatorsRec(nn);
   }
-
-  switch (k)
-  {
-    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)
-      {
-        Node divByZeroNum = getArithSkolemApp(num, ArithSkolemId::DIV_BY_ZERO);
-        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
-        ret = nm->mkNode(kind::ITE, denEq0, divByZeroNum, ret);
-      }
-      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)
-      {
-        Node intDivByZeroNum =
-            getArithSkolemApp(num, ArithSkolemId::INT_DIV_BY_ZERO);
-        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
-        ret = nm->mkNode(kind::ITE, denEq0, intDivByZeroNum, ret);
-      }
-      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)
-      {
-        Node modZeroNum = getArithSkolemApp(num, ArithSkolemId::MOD_BY_ZERO);
-        Node denEq0 = nm->mkNode(kind::EQUAL, den, nm->mkConst(Rational(0)));
-        ret = nm->mkNode(kind::ITE, denEq0, modZeroNum, ret);
-      }
-      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::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->mkBoundVar(nm->realType());
-        Node lem;
-        if (k == kind::SQRT)
-        {
-          Node skolemApp = getArithSkolemApp(node[0], ArithSkolemId::SQRT);
-          Node uf = skolemApp.eqNode(var);
-          Node nonNeg = nm->mkNode(
-              kind::AND, nm->mkNode(kind::MULT, var, var).eqNode(node[0]), uf);
-
-          // sqrt(x) reduces to:
-          // witness y. ite(x >= 0.0, y * y = x ^ Uf(x), Uf(x))
-          //
-          // Uf(x) makes sure that the reduction still behaves like a function,
-          // otherwise the reduction of (x = 1) ^ (sqrt(x) != sqrt(1)) would be
-          // satisfiable. On the original formula, this would require that we
-          // simultaneously interpret sqrt(1) as 1 and -1, which is not a valid
-          // model.
-          lem = nm->mkNode(
-              kind::ITE,
-              nm->mkNode(kind::GEQ, node[0], nm->mkConst(Rational(0))),
-              nonNeg,
-              uf);
-        }
-        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));
-          }
-
-          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);
-          lem = nm->mkNode(AND, rlem, invTerm.eqNode(node[0]));
-        }
-        Assert(!lem.isNull());
-        Node ret = nm->mkNode(WITNESS, nm->mkNode(BOUND_VAR_LIST, var), lem);
-        d_nlin_inverse_skolem[node] = ret;
-        return ret;
-      }
-      return (*it).second;
-      break;
-    }
-
-    default: return node; break;
-  }
-
-  Unreachable();
+  return node;
 }
 
 Node TheoryArithPrivate::getArithSkolem(ArithSkolemId asi)
diff --git a/src/theory/arith/theory_arith_private.h b/src/theory/arith/theory_arith_private.h
index b5debe76d..8198dbcf1 100644
--- a/src/theory/arith/theory_arith_private.h
+++ b/src/theory/arith/theory_arith_private.h
@@ -189,9 +189,6 @@ public:
     d_setupNodes.insert(n);
   }
 private:
-
-  void setupDivLike(const Variable& x);
-
   void setupVariable(const Variable& x);
   void setupVariableList(const VarList& vl);
   void setupPolynomial(const Polynomial& poly);
@@ -422,7 +419,39 @@ private:
   // handle linear /, div, mod, and also is_int, to_int
   Node ppRewriteTerms(TNode atom);
 
-public:
+  /**
+   * Called when a non-linear term n is given to this class. Throw an exception
+   * if the logic is linear.
+   */
+  void checkNonLinearLogic(Node term);
+  /**
+   * Eliminate operators in term n. If n has top symbol that is not a core
+   * one (including division, int division, mod, to_int, is_int, syntactic sugar
+   * transcendental functions), then we replace it by a form that eliminates
+   * that operator. This may involve the introduction of witness terms.
+   *
+   * One exception to the above rule is that we may leave certain applications
+   * like (/ 4 1) unchanged, since replacing this by 4 changes its type from
+   * real to int. This is important for some subtyping issues during
+   * expandDefinition. Moreover, applications like this can be eliminated
+   * trivially later by rewriting.
+   *
+   * This method is called both during expandDefinition and during ppRewrite.
+   *
+   * @param n The node to eliminate operators from.
+   * @return The (single step) eliminated form of n.
+   */
+  Node eliminateOperators(Node n);
+  /**
+   * Recursively ensure that n has no non-standard operators. This applies
+   * the above method on all subterms of n.
+   *
+   * @param n The node to eliminate operators from.
+   * @return The eliminated form of n.
+   */
+  Node eliminateOperatorsRec(Node n);
+
+ public:
   TheoryArithPrivate(TheoryArith& containing, context::Context* c, context::UserContext* u, OutputChannel& out, Valuation valuation, const LogicInfo& logicInfo);
   ~TheoryArithPrivate();