This commit merges into trunk the branch branches/arithmetic/integers2 from r2650 to r2779.

- This excludes revision 2777.  This revision had some strange performance implications and was delaying the merge.
- This includes the new DioSolver. The DioSolver can discover conflicts, produce substitutions, and produce cuts.
- The DioSolver can be disabled at command line using --disable-dio-solver.
- This includes a number of changes to the arithmetic normal form.
- The Integer class features a number of new number theoretic function.
- This commit includes a few rather loud warning. I will do my best to take care of them today.

16 files changed, 2578 insertions, 326 deletions

diff --git a/src/theory/arith/arith_rewriter.cpp b/src/theory/arith/arith_rewriter.cpp
index 66223b479..a309b9403 100644
--- a/src/theory/arith/arith_rewriter.cpp
+++ b/src/theory/arith/arith_rewriter.cpp
@@ -182,37 +182,43 @@ RewriteResponse ArithRewriter::postRewriteAtomConstantRHS(TNode t){
   TNode left  = t[0];
   TNode right = t[1];
 
-  Comparison cmp = Comparison::mkComparison(t.getKind(), Polynomial::parsePolynomial(left), Constant(right));
+  Comparison cmp = Comparison::mkNormalComparison(t.getKind(), Polynomial::parsePolynomial(left), Constant(right));
 
-  if(cmp.isBoolean()){
-    return RewriteResponse(REWRITE_DONE, cmp.getNode());
-  }
+  Assert(cmp.isNormalForm());
+  return RewriteResponse(REWRITE_DONE, cmp.getNode());
 
-  if(cmp.getLeft().containsConstant()){
-    Monomial constantHead = cmp.getLeft().getHead();
-    Assert(constantHead.isConstant());
 
-    Constant constant = constantHead.getConstant();
+  // Comparison cmp = Comparison::mkComparison(t.getKind(), Polynomial::parsePolynomial(left), Constant(right));
 
-    Constant negativeConstantHead = -constant;
+  // if(cmp.isBoolean()){
+  //   return RewriteResponse(REWRITE_DONE, cmp.getNode());
+  // }
 
-    cmp = cmp.addConstant(negativeConstantHead);
-  }
-  Assert(!cmp.getLeft().containsConstant());
+  // if(cmp.getLeft().containsConstant()){
+  //   Monomial constantHead = cmp.getLeft().getHead();
+  //   Assert(constantHead.isConstant());
 
-  if(!cmp.getLeft().getHead().coefficientIsOne()){
-    Monomial constantHead = cmp.getLeft().getHead();
-    Assert(!constantHead.isConstant());
-    Constant constant = constantHead.getConstant();
+  //   Constant constant = constantHead.getConstant();
 
-    Constant inverse = Constant::mkConstant(constant.getValue().inverse());
+  //   Constant negativeConstantHead = -constant;
 
-    cmp = cmp.multiplyConstant(inverse);
-  }
-  Assert(cmp.getLeft().getHead().coefficientIsOne());
+  //   cmp = cmp.addConstant(negativeConstantHead);
+  // }
+  // Assert(!cmp.getLeft().containsConstant());
 
-  Assert(cmp.isBoolean() || cmp.isNormalForm());
-  return RewriteResponse(REWRITE_DONE, cmp.getNode());
+  // if(!cmp.getLeft().getHead().coefficientIsOne()){
+  //   Monomial constantHead = cmp.getLeft().getHead();
+  //   Assert(!constantHead.isConstant());
+  //   Constant constant = constantHead.getConstant();
+
+  //   Constant inverse = Constant::mkConstant(constant.getValue().inverse());
+
+  //   cmp = cmp.multiplyConstant(inverse);
+  // }
+  // Assert(cmp.getLeft().getHead().coefficientIsOne());
+
+  // Assert(cmp.isBoolean() || cmp.isNormalForm());
+  // return RewriteResponse(REWRITE_DONE, cmp.getNode());
 }
 
 RewriteResponse ArithRewriter::postRewriteAtom(TNode atom){
diff --git a/src/theory/arith/arith_utilities.h b/src/theory/arith/arith_utilities.h
index 42f4704df..0168ee043 100644
--- a/src/theory/arith/arith_utilities.h
+++ b/src/theory/arith/arith_utilities.h
@@ -49,11 +49,17 @@ typedef __gnu_cxx::hash_map<ArithVar, Node> ArithVarToNodeMap;
 typedef __gnu_cxx::hash_set<Node, NodeHashFunction> NodeSet;
 typedef context::CDSet<Node, NodeHashFunction> CDNodeSet;
 
+typedef context::CDSet<ArithVar> CDArithVarSet;
+
 
 inline Node mkRationalNode(const Rational& q){
   return NodeManager::currentNM()->mkConst<Rational>(q);
 }
 
+inline Node mkIntegerNode(const Integer& z){
+  return NodeManager::currentNM()->mkConst<Integer>(z);
+}
+
 inline Node mkBoolNode(bool b){
   return NodeManager::currentNM()->mkConst<bool>(b);
 }
diff --git a/src/theory/arith/delta_rational.h b/src/theory/arith/delta_rational.h
index fdc3bee3b..682d13720 100644
--- a/src/theory/arith/delta_rational.h
+++ b/src/theory/arith/delta_rational.h
@@ -118,6 +118,38 @@ public:
     return *(this);
   }
 
+  bool isIntegral() const {
+    if(getInfinitesimalPart().sgn() == 0){
+      return getNoninfinitesimalPart().isIntegral();
+    }else{
+      return false;
+    }
+  }
+
+  Integer floor() const {
+    if(getNoninfinitesimalPart().isIntegral()){
+      if(getInfinitesimalPart().sgn() >= 0){
+        return getNoninfinitesimalPart().getNumerator();
+      }else{
+        return getNoninfinitesimalPart().getNumerator() - Integer(1);
+      }
+    }else{
+      return getNoninfinitesimalPart().floor();
+    }
+  }
+
+  Integer ceiling() const {
+    if(getNoninfinitesimalPart().isIntegral()){
+      if(getInfinitesimalPart().sgn() <= 0){
+        return getNoninfinitesimalPart().getNumerator();
+      }else{
+        return getNoninfinitesimalPart().getNumerator() + Integer(1);
+      }
+    }else{
+      return getNoninfinitesimalPart().ceiling();
+    }
+  }
+
   std::string toString() const;
 
 };
diff --git a/src/theory/arith/dio_solver.cpp b/src/theory/arith/dio_solver.cpp
index 151859f21..dead34807 100644
--- a/src/theory/arith/dio_solver.cpp
+++ b/src/theory/arith/dio_solver.cpp
@@ -26,91 +26,766 @@ namespace CVC4 {
 namespace theory {
 namespace arith {
 
-/*
-static void printEquation(vector<Rational>& coeffs,
-                          vector<ArithVar>& variables,
-                          ostream& out) {
-  Assert(coeffs.size() == variables.size());
-  vector<Rational>::const_iterator i = coeffs.begin();
-  vector<ArithVar>::const_iterator j = variables.begin();
-  while(i != coeffs.end()) {
-    out << *i << "*" << *j;
-    ++i;
-    ++j;
-    if(i != coeffs.end()) {
-      out << " + ";
+inline Node makeIntegerVariable(){
+  NodeManager* curr = NodeManager::currentNM();
+  return curr->mkVar(curr->integerType());
+}
+
+DioSolver::DioSolver(context::Context* ctxt) :
+  d_lastUsedVariable(ctxt,0),
+  d_inputConstraints(ctxt),
+  d_nextInputConstraintToEnqueue(ctxt, 0),
+  d_trail(ctxt),
+  d_subs(ctxt),
+  d_currentF(),
+  d_savedQueue(ctxt),
+  d_savedQueueIndex(ctxt, 0),
+  d_conflictHasBeenRaised(ctxt, false),
+  d_maxInputCoefficientLength(ctxt, 0),
+  d_usedDecomposeIndex(ctxt, false),
+  d_lastPureSubstitution(ctxt, 0),
+  d_pureSubstitionIter(ctxt, 0)
+{}
+
+DioSolver::Statistics::Statistics() :
+  d_conflictCalls("theory::arith::dio::conflictCalls",0),
+  d_cutCalls("theory::arith::dio::cutCalls",0),
+  d_cuts("theory::arith::dio::cuts",0),
+  d_conflicts("theory::arith::dio::conflicts",0),
+  d_conflictTimer("theory::arith::dio::conflictTimer"),
+  d_cutTimer("theory::arith::dio::cutTimer")
+{
+  StatisticsRegistry::registerStat(&d_conflictCalls);
+  StatisticsRegistry::registerStat(&d_cutCalls);
+
+  StatisticsRegistry::registerStat(&d_cuts);
+  StatisticsRegistry::registerStat(&d_conflicts);
+
+  StatisticsRegistry::registerStat(&d_conflictTimer);
+  StatisticsRegistry::registerStat(&d_cutTimer);
+}
+
+DioSolver::Statistics::~Statistics(){
+  StatisticsRegistry::unregisterStat(&d_conflictCalls);
+  StatisticsRegistry::unregisterStat(&d_cutCalls);
+
+  StatisticsRegistry::unregisterStat(&d_cuts);
+  StatisticsRegistry::unregisterStat(&d_conflicts);
+
+  StatisticsRegistry::unregisterStat(&d_conflictTimer);
+  StatisticsRegistry::unregisterStat(&d_cutTimer);
+}
+
+size_t DioSolver::allocateVariableInPool() {
+  Assert(d_lastUsedVariable <= d_variablePool.size());
+  if(d_lastUsedVariable == d_variablePool.size()){
+    Assert(d_lastUsedVariable == d_variablePool.size());
+    Node intVar = makeIntegerVariable();
+    d_variablePool.push_back(Variable(intVar));
+  }
+  size_t res = d_lastUsedVariable;
+  d_lastUsedVariable = d_lastUsedVariable + 1;
+  return res;
+}
+
+
+  Node DioSolver::nextPureSubstitution(){
+    Assert(hasMorePureSubstitutions());
+    SubIndex curr = d_pureSubstitionIter;
+    d_pureSubstitionIter = d_pureSubstitionIter + 1;
+
+    Assert(d_subs[curr].d_fresh.isNull());
+    Variable v = d_subs[curr].d_eliminated;
+
+    SumPair sp = d_trail[d_subs[curr].d_constraint].d_eq;
+    Polynomial p = sp.getPolynomial();
+    Constant c = -sp.getConstant();
+    Polynomial cancelV = p + Polynomial::mkPolynomial(v);
+    Node eq = NodeManager::currentNM()->mkNode(kind::EQUAL, v.getNode(), cancelV.getNode());
+    return eq;
+  }
+
+
+bool DioSolver::debugEqualityInInputEquations(Node eq){
+  typedef context::CDList<InputConstraint>::const_iterator const_iterator;
+  const_iterator i=d_inputConstraints.begin(), end = d_inputConstraints.end();
+  for(; i != end; ++i){
+    Node reason_i = (*i).d_reason;
+    if(eq == reason_i){
+      return true;
     }
   }
-  out << " = 0";
+  return false;
 }
-*/
 
-bool DioSolver::solve() {
-  Trace("integers") << "DioSolver::solve()" << endl;
-  if(Debug.isOn("integers")) {
-    ScopedDebug dtab("tableau");
-    d_tableau.printTableau();
+bool DioSolver::acceptableOriginalNodes(Node n){
+  Kind k = n.getKind();
+  if(k == kind::EQUAL){
+    return true;
+  }else if(k == kind::AND){
+    Node ub = n[0];
+    Node lb = n[1];
+    Kind kub = simplifiedKind(ub);
+    Kind klb = simplifiedKind(lb);
+    return (kub == kind::LEQ || kub==kind::LT) && (klb == kind::GEQ || klb == kind::GT);
+  }else{
+    return false;
   }
-  for(ArithVarSet::const_iterator i = d_tableau.beginBasic();
-      i != d_tableau.endBasic();
-      ++i) {
-    Debug("integers") << "going through row " << *i << endl;
+}
+
+void DioSolver::pushInputConstraint(const Comparison& eq, Node reason){
+  Assert(!debugEqualityInInputEquations(reason));
+  Assert(eq.isIntegral());
+  Assert(eq.getNode().getKind() == kind::EQUAL);
+  Assert(acceptableOriginalNodes(reason));
+
+  SumPair sp = SumPair::comparisonToSumPair(eq);
+  uint32_t length = sp.maxLength();
+  if(length > d_maxInputCoefficientLength){
+    d_maxInputCoefficientLength = length;
+  }
+
+  size_t varIndex = allocateVariableInPool();
+  Variable proofVariable(d_variablePool[varIndex]);
+
+  TrailIndex posInTrail = d_trail.size();
+  d_trail.push_back(Constraint(sp,Polynomial(Monomial(VarList(proofVariable)))));
+
+  size_t posInConstraintList = d_inputConstraints.size();
+  d_inputConstraints.push_back(InputConstraint(reason, posInTrail));
 
-    Integer m = 1;
-    for(Tableau::RowIterator j = d_tableau.rowIterator(*i); !j.atEnd(); ++j) {
-      Debug("integers") << "  column " << (*j).getCoefficient() << " * " << (*j).getColVar() << endl;
-      m *= (*j).getCoefficient().getDenominator();
+  d_varToInputConstraintMap[proofVariable.getNode()] = posInConstraintList;
+}
+
+
+DioSolver::TrailIndex DioSolver::scaleEqAtIndex(DioSolver::TrailIndex i, const Integer& g){
+  Assert(g != 0);
+  Constant invg = Constant::mkConstant(Rational(Integer(1),g));
+  const SumPair& sp = d_trail[i].d_eq;
+  const Polynomial& proof = d_trail[i].d_proof;
+
+  SumPair newSP = sp * invg;
+  Polynomial newProof = proof * invg;
+
+  Assert(newSP.isIntegral());
+  Assert(newSP.gcd() == 1);
+
+  TrailIndex j = d_trail.size();
+
+  d_trail.push_back(Constraint(newSP, newProof));
+
+  Debug("arith::dio") << "scaleEqAtIndex(" << i <<","<<g<<")"<<endl;
+  Debug("arith::dio") << "derived "<< newSP.getNode()
+                      <<" with proof " << newProof.getNode() << endl;
+  return j;
+}
+
+Node DioSolver::proveIndex(TrailIndex i){
+  Assert(inRange(i));
+  const Polynomial& proof = d_trail[i].d_proof;
+  Assert(!proof.isConstant());
+
+  NodeBuilder<> nb(kind::AND);
+  for(Polynomial::iterator iter = proof.begin(), end = proof.end(); iter!= end; ++iter){
+    Monomial m = (*iter);
+    Assert(!m.isConstant());
+    VarList vl = m.getVarList();
+    Assert(vl.singleton());
+    Variable v = vl.getHead();
+
+    Node input = proofVariableToReason(v);
+    Assert(acceptableOriginalNodes(input));
+    if(input.getKind() == kind::AND){
+      nb << input[0] << input[1];
+    }else{
+      nb << input;
     }
-    Assert(m >= 1);
-    Debug("integers") << "final multiplier is " << m << endl;
-
-    Integer gcd = 0;
-    Rational c = 0;
-    Debug("integers") << "going throw row " << *i << " to find gcd" << endl;
-    for(Tableau::RowIterator j = d_tableau.rowIterator(*i); !j.atEnd(); ++j) {
-      Rational r = (*j).getCoefficient();
-      ArithVar v = (*j).getColVar();
-      r *= m;
-      Debug("integers") << "  column " << r << " * " << v << endl;
-      Assert(r.getDenominator() == 1);
-      bool foundConstant = false;
-      // The tableau doesn't store constants.  Constants int he input
-      // are mapped to slack variables that are constrained with
-      // bounds in the partial model.  So we detect and accumulate all
-      // variables whose upper bound equals their lower bound, which
-      // catches these input constants as well as any (contextually)
-      // eliminated variables.
-      if(d_partialModel.hasUpperBound(v) && d_partialModel.hasLowerBound(v)) {
-        DeltaRational b = d_partialModel.getLowerBound(v);
-        if(b.getInfinitesimalPart() == 0 && b == d_partialModel.getUpperBound(v)) {
-          r *= b.getNoninfinitesimalPart();
-          Debug("integers") << "    var " << v << " == [" << b << "], so found a constant part " << r << endl;
-          c += r;
-          foundConstant = true;
+  }
+
+  Node result = (nb.getNumChildren() == 1) ? nb[0] : (Node)nb;
+  Debug("arith::dio") << "Proof at " << i << " is "
+                      << d_trail[i].d_eq.getNode() << endl
+                      << d_trail[i].d_proof.getNode() << endl
+                      << " which becomes " << result << endl;
+  return result;
+}
+
+bool DioSolver::anyCoefficientExceedsMaximum(TrailIndex j) const{
+  uint32_t length = d_trail[j].d_eq.maxLength();
+  uint32_t nmonos = d_trail[j].d_eq.getPolynomial().numMonomials();
+
+  bool result =
+    nmonos >= 2 &&
+    length > d_maxInputCoefficientLength + MAX_GROWTH_RATE;
+  if(Debug.isOn("arith::dio::max") && result){
+    Debug("arith::dio::max") << "about to drop:";
+    debugPrintTrail(j);
+  }
+  return result;
+}
+
+void DioSolver::enqueueInputConstraints(){
+  Assert(d_currentF.empty());
+  while(d_savedQueueIndex < d_savedQueue.size()){
+    d_currentF.push_back(d_savedQueue[d_savedQueueIndex]);
+    d_savedQueueIndex = d_savedQueueIndex + 1;
+  }
+
+  while(d_nextInputConstraintToEnqueue < d_inputConstraints.size()  && !inConflict()){
+    size_t curr = d_nextInputConstraintToEnqueue;
+    d_nextInputConstraintToEnqueue = d_nextInputConstraintToEnqueue + 1;
+
+    TrailIndex i = d_inputConstraints[curr].d_trailPos;
+    TrailIndex j = applyAllSubstitutionsToIndex(i);
+
+    if(!(triviallySat(j)  || anyCoefficientExceedsMaximum(j))){
+      if(triviallyUnsat(j)){
+        raiseConflict(j);
+      }else{
+        TrailIndex k = reduceByGCD(j);
+
+        if(!inConflict()){
+          if(triviallyUnsat(k)){
+            raiseConflict(k);
+          }else if(!triviallySat(k)){
+            pushToQueueBack(k);
+          }
         }
       }
-      if(!foundConstant) {
-        // calculate gcd of all (NON-constant) coefficients
-        gcd = (gcd == 0) ? r.getNumerator().abs() : gcd.gcd(r.getNumerator());
-        Debug("integers") << "    gcd now " << gcd << endl;
+    }
+  }
+}
+
+
+/*TODO Currently linear in the size of the Queue
+ *It is not clear if am O(log n) strategy would be better.
+ *Right before this in the algorithm is a substition which could potentially
+ *effect the key values of everything in the queue.
+ */
+void DioSolver::moveMinimumByAbsToQueueFront(){
+  Assert(!queueEmpty());
+
+
+  //Select the minimum element.
+  size_t indexInQueue = 0;
+  Monomial minMonomial = d_trail[d_currentF[indexInQueue]].d_minimalMonomial;
+
+  size_t N = d_currentF.size();
+  for(size_t i=1; i < N; ++i){
+    Monomial curr = d_trail[d_currentF[i]].d_minimalMonomial;
+    if(curr.absLessThan(minMonomial)){
+      indexInQueue = i;
+      minMonomial = curr;
+    }
+  }
+
+  TrailIndex tmp = d_currentF[indexInQueue];
+  d_currentF[indexInQueue] = d_currentF.front();
+  d_currentF.front() = tmp;
+}
+
+bool DioSolver::queueEmpty() const{
+  return d_currentF.empty();
+}
+
+Node DioSolver::columnGcdIsOne() const{
+  std::hash_map<Node, Integer, NodeHashFunction> gcdMap;
+
+  std::deque<TrailIndex>::const_iterator iter, end;
+  for(iter = d_currentF.begin(), end = d_currentF.end(); iter != end; ++iter){
+    TrailIndex curr = *iter;
+    Polynomial p = d_trail[curr].d_eq.getPolynomial();
+    Polynomial::iterator monoIter = p.begin(), monoEnd = p.end();
+    for(; monoIter != monoEnd; ++monoIter){
+      Monomial m = *monoIter;
+      VarList vl = m.getVarList();
+      Node vlNode = vl.getNode();
+
+      Constant c = m.getConstant();
+      Integer zc = c.getValue().getNumerator();
+      if(gcdMap.find(vlNode) == gcdMap.end()){
+        // have not seen vl yet.
+        // gcd is c
+        Assert(c.isIntegral());
+        Assert(!m.absCoefficientIsOne());
+        gcdMap.insert(make_pair(vlNode, zc.abs()));
+      }else{
+        const Integer& currentGcd = gcdMap[vlNode];
+        Integer newGcd = currentGcd.gcd(zc);
+        if(newGcd == 1){
+          return vlNode;
+        }else{
+          gcdMap[vlNode] = newGcd;
+        }
       }
     }
-    Debug("integers") << "addEquation: gcd is " << gcd << ", c is " << c << endl;
-    // The gcd must divide c for this equation to be satisfiable.
-    // If c is not an integer, there's no way it can.
-    if(c.getDenominator() == 1 && gcd == gcd.gcd(c.getNumerator())) {
-      Trace("integers") << "addEquation: this eqn is fine" << endl;
-    } else {
-      Trace("integers") << "addEquation: eqn not satisfiable, returning false" << endl;
-      return false;
+  }
+  return Node::null();
+}
+
+void DioSolver::saveQueue(){
+  std::deque<TrailIndex>::const_iterator iter, end;
+  for(iter = d_currentF.begin(), end = d_currentF.end(); iter != end; ++iter){
+    d_savedQueue.push_back(*iter);
+  }
+}
+
+DioSolver::TrailIndex DioSolver::impliedGcdOfOne(){
+  Node canReduce = columnGcdIsOne();
+  if(canReduce.isNull()){
+    return 0;
+  }else{
+    VarList vl = VarList::parseVarList(canReduce);
+
+    TrailIndex current;
+    Integer currentCoeff, currentGcd;
+
+    //step 1 find the first equation containing vl
+    //Set current and currentCoefficient
+    std::deque<TrailIndex>::const_iterator iter, end;
+    for(iter = d_currentF.begin(), end = d_currentF.end(); true; ++iter){
+      Assert(iter != end);
+      current = *iter;
+      Constant coeff = d_trail[current].d_eq.getPolynomial().getCoefficient(vl);
+      if(!coeff.isZero()){
+        currentCoeff = coeff.getValue().getNumerator();
+        currentGcd = currentCoeff.abs();
+
+        ++iter;
+        break;
+      }
+    }
+
+    //For the rest of the equations keep reducing until the coefficient is one
+    for(; iter != end; ++iter){
+      TrailIndex inQueue = *iter;
+      Constant iqc = d_trail[inQueue].d_eq.getPolynomial().getCoefficient(vl);
+      if(!iqc.isZero()){
+        Integer inQueueCoeff = iqc.getValue().getNumerator();
+
+        //mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, mpz_t a, mpz_t b);
+        Integer g, s, t;
+        // g = a*s + b*t
+        Integer::extendedGcd(g, s, t, currentCoeff, inQueueCoeff);
+
+        Debug("arith::dio") << "extendedReduction" << endl;
+        Debug("arith::dio") << g << " = " << s <<"*"<< currentCoeff << " + " << t <<"*"<< inQueueCoeff << endl;
+        
+        Assert(g <= currentGcd);
+        if(g < currentGcd){
+          if(s.sgn() == 0){
+            Debug("arith::dio") << "extendedReduction drop" << endl;
+            Assert(inQueueCoeff.divides(currentGcd));
+            current = *iter;
+            currentCoeff = inQueueCoeff;
+            currentGcd = inQueueCoeff;
+          }else{
+            Debug("arith::dio") << "extendedReduction combine" << endl;
+          
+            TrailIndex next = combineEqAtIndexes(current, s, inQueue, t);
+
+            Assert(d_trail[next].d_eq.getPolynomial().getCoefficient(vl).getValue().getNumerator() == g);
+
+            current = next;
+            currentCoeff = g;
+            currentGcd = g;
+            if(currentGcd == 1){
+              return current;
+            }
+          }
+        }
+      }
+    }
+    //This is not reachble as it is assured that the gcd of the column is 1
+    Unreachable();
+  }
+}
+
+bool DioSolver::processEquations(bool allowDecomposition){
+  Assert(!inConflict());
+
+  enqueueInputConstraints();
+  while(! queueEmpty() && !inConflict()){
+    moveMinimumByAbsToQueueFront();
+
+    TrailIndex minimum = d_currentF.front();
+    TrailIndex reduceIndex;
+
+    Assert(inRange(minimum));
+    Assert(!inConflict());
+
+    Debug("arith::dio") << "processEquations " << minimum << " : " << d_trail[minimum].d_eq.getNode() << endl;
+
+    Assert(queueConditions(minimum));
+
+    bool canDirectlySolve = d_trail[minimum].d_minimalMonomial.absCoefficientIsOne();
+
+    std::pair<SubIndex, TrailIndex> p;
+    if(canDirectlySolve){
+      d_currentF.pop_front();
+      p = solveIndex(minimum);
+      reduceIndex = minimum;
+    }else{
+      TrailIndex implied = impliedGcdOfOne();
+      
+      if(implied != 0){
+        p = solveIndex(implied);
+        reduceIndex = implied;
+      }else if(allowDecomposition){
+        d_currentF.pop_front();
+        p = decomposeIndex(minimum);
+        reduceIndex = minimum;
+      }else {
+        // Cannot make progress without decomposeIndex
+        saveQueue();
+        break;
+      }
+    }
+
+    SubIndex subIndex = p.first;
+    TrailIndex next = p.second;
+    subAndReduceCurrentFByIndex(subIndex);
+
+    if(next != reduceIndex){
+      if(triviallyUnsat(next)){
+        raiseConflict(next);
+      }else if(! triviallySat(next) ){
+        pushToQueueBack(next);
+      }
     }
   }
 
-  return true;
+  d_currentF.clear();
+  return inConflict();
 }
 
-void DioSolver::getSolution() {
-  Unimplemented();
+Node DioSolver::processEquationsForConflict(){
+  TimerStat::CodeTimer codeTimer(d_statistics.d_conflictTimer);
+  ++(d_statistics.d_conflictCalls);
+
+  Assert(!inConflict());
+  if(processEquations(false)){
+    ++(d_statistics.d_conflicts);
+    return proveIndex(getConflictIndex());
+  }else{
+    return Node::null();
+  }
+}
+
+SumPair DioSolver::processEquationsForCut(){
+  TimerStat::CodeTimer codeTimer(d_statistics.d_cutTimer);
+  ++(d_statistics.d_cutCalls);
+
+  Assert(!inConflict());
+  if(processEquations(true)){
+    ++(d_statistics.d_cuts);
+    return purifyIndex(getConflictIndex());
+  }else{
+    return SumPair::mkZero();
+  }
+}
+
+
+SumPair DioSolver::purifyIndex(TrailIndex i){
+#warning "This uses the substition trail to reverse the substitions from the sum term. Using the proof term should be more efficient."
+
+  SumPair curr = d_trail[i].d_eq;
+
+  Constant negOne = Constant::mkConstant(-1);
+
+  for(uint32_t revIter = d_subs.size(); revIter > 0; --revIter){
+    uint32_t i = revIter - 1;
+    Node freshNode = d_subs[i].d_fresh;
+    if(freshNode.isNull()){
+      continue;
+    }else{
+      Variable var(freshNode);
+      Polynomial vsum = curr.getPolynomial();
+
+      Constant a = vsum.getCoefficient(VarList(var));
+      if(!a.isZero()){
+        const SumPair& sj = d_trail[d_subs[i].d_constraint].d_eq;
+        Assert(sj.getPolynomial().getCoefficient(VarList(var)).isOne());
+        SumPair newSi = (curr * negOne) + (sj * a);
+        Assert(newSi.getPolynomial().getCoefficient(VarList(var)).isZero());
+        curr = newSi;
+      }
+    }
+  }
+  return curr;
+}
+
+DioSolver::TrailIndex DioSolver::combineEqAtIndexes(DioSolver::TrailIndex i, const Integer& q, DioSolver::TrailIndex j, const Integer& r){
+  Constant cq = Constant::mkConstant(q);
+  Constant cr = Constant::mkConstant(r);
+
+  const SumPair& si = d_trail[i].d_eq;
+  const SumPair& sj = d_trail[j].d_eq;
+
+  Debug("arith::dio") << "combineEqAtIndexes(" << i <<","<<q<<","<<j<<","<<r<<")"<<endl;
+  Debug("arith::dio") << "d_facts[i] = " << si.getNode() << endl
+                      << "d_facts[j] = " << sj.getNode() << endl;
+
+
+  SumPair newSi = (si * cq) + (sj * cr);
+
+
+  const Polynomial& pi = d_trail[i].d_proof;
+  const Polynomial& pj = d_trail[j].d_proof;
+  Polynomial newPi = (pi * cq) + (pj * cr);
+
+  TrailIndex k = d_trail.size();
+  d_trail.push_back(Constraint(newSi, newPi));
+
+
+  Debug("arith::dio") << "derived "<< newSi.getNode()
+                      <<" with proof " << newPi.getNode() << endl;
+
+  return k;
+
+}
+
+void DioSolver::printQueue(){
+  Debug("arith::dio") << "DioSolver::printQueue()" << endl;
+  for(TrailIndex i = 0, last = d_trail.size(); i < last; ++i){
+    Debug("arith::dio") << "d_trail[i].d_eq = " << d_trail[i].d_eq.getNode() << endl;
+    Debug("arith::dio") << "d_trail[i].d_proof = " << d_trail[i].d_proof.getNode() << endl;
+  }
+
+  Debug("arith::dio") << "DioSolver::printSubs()" << endl;
+  for(SubIndex si=0, sN=d_subs.size(); si < sN; ++si){
+    Debug("arith::dio") << "d_subs[i] = {"
+                        << "d_fresh="<< d_subs[si].d_fresh <<","
+                        << "d_eliminated="<< d_subs[si].d_eliminated.getNode() <<","
+                        << "d_constraint="<< d_subs[si].d_constraint <<"}" << endl;
+    Debug("arith::dio") << "d_trail[d_subs[i].d_constraint].d_eq="
+                        << d_trail[d_subs[si].d_constraint].d_eq.getNode() << endl;
+  }
+}
+
+DioSolver::TrailIndex DioSolver::applyAllSubstitutionsToIndex(DioSolver::TrailIndex trailIndex){
+  TrailIndex currentIndex = trailIndex;
+  for(SubIndex subIter = 0, siEnd = d_subs.size(); subIter < siEnd; ++subIter){
+    currentIndex = applySubstitution(subIter, currentIndex);
+  }
+  return currentIndex;
+}
+
+bool DioSolver::debugSubstitutionApplies(DioSolver::SubIndex si, DioSolver::TrailIndex ti){
+  Variable var = d_subs[si].d_eliminated;
+  TrailIndex subIndex = d_subs[si].d_constraint;
+
+  const SumPair& curr = d_trail[ti].d_eq;
+  Polynomial vsum = curr.getPolynomial();
+
+  Constant a = vsum.getCoefficient(VarList(var));
+  return !a.isZero();
+}
+
+bool DioSolver::debugAnySubstitionApplies(DioSolver::TrailIndex i){
+  for(SubIndex subIter = 0, siEnd = d_subs.size(); subIter < siEnd; ++subIter){
+    if(debugSubstitutionApplies(subIter, i)){
+      return true;
+    }
+  }
+  return false;
+}
+
+std::pair<DioSolver::SubIndex, DioSolver::TrailIndex> DioSolver::solveIndex(DioSolver::TrailIndex i){
+  const SumPair& si = d_trail[i].d_eq;
+
+  Debug("arith::dio") << "before solveIndex("<<i<<":"<<si.getNode()<< ")" << endl;
+
+  const Polynomial& p = si.getPolynomial();
+  Assert(p.isIntegral());
+
+  Assert(p.selectAbsMinimum() == d_trail[i].d_minimalMonomial);
+  const Monomial& av = d_trail[i].d_minimalMonomial;
+
+  VarList vl = av.getVarList();
+  Assert(vl.singleton());
+  Variable var = vl.getHead();
+  Constant a = av.getConstant();
+  Integer a_abs = a.getValue().getNumerator().abs();
+
+  Assert(a_abs == 1);
+
+  TrailIndex ci = !a.isNegative() ? scaleEqAtIndex(i, Integer(-1)) : i;
+
+  SubIndex subBy = d_subs.size();
+  d_subs.push_back(Substitution(Node::null(), var, ci));
+
+  Debug("arith::dio") << "after solveIndex " <<  d_trail[ci].d_eq.getNode() << " for " << av.getNode() << endl;
+  Assert(d_trail[ci].d_eq.getPolynomial().getCoefficient(vl) == Constant::mkConstant(-1));
+
+  return make_pair(subBy, i);
+}
+
+std::pair<DioSolver::SubIndex, DioSolver::TrailIndex> DioSolver::decomposeIndex(DioSolver::TrailIndex i){
+  const SumPair& si = d_trail[i].d_eq;
+
+  d_usedDecomposeIndex = true;
+
+  Debug("arith::dio") << "before decomposeIndex("<<i<<":"<<si.getNode()<< ")" << endl;
+
+  const Polynomial& p = si.getPolynomial();
+  Assert(p.isIntegral());
+
+  Assert(p.selectAbsMinimum() == d_trail[i].d_minimalMonomial);
+  const Monomial& av = d_trail[i].d_minimalMonomial;
+
+  VarList vl = av.getVarList();
+  Assert(vl.singleton());
+  Variable var = vl.getHead();
+  Constant a = av.getConstant();
+  Integer a_abs = a.getValue().getNumerator().abs();
+
+  Assert(a_abs > 1);
+  
+  //It is not sufficient to reduce the case where abs(a) == 1 to abs(a) > 1.
+  //We need to handle both cases seperately to ensure termination.
+  Node qr = SumPair::computeQR(si, a.getValue().getNumerator());
+
+  Assert(qr.getKind() == kind::PLUS);
+  Assert(qr.getNumChildren() == 2);
+  SumPair q = SumPair::parseSumPair(qr[0]);
+  SumPair r = SumPair::parseSumPair(qr[1]);
+
+  Assert(q.getPolynomial().getCoefficient(vl) == Constant::mkConstant(1));
+
+  Assert(!r.isZero());
+  Node freshNode = makeIntegerVariable();
+  Variable fresh(freshNode);
+  SumPair fresh_one=SumPair::mkSumPair(fresh);
+  SumPair fresh_a = fresh_one * a;
+
+  SumPair newSI = SumPair(fresh_one) - q;
+  // this normalizes the coefficient of var to -1
+
+
+  TrailIndex ci = d_trail.size();
+  d_trail.push_back(Constraint(newSI, Polynomial::mkZero()));
+
+  Debug("arith::dio") << "Decompose ci(" << ci <<":" <<  d_trail[ci].d_eq.getNode()
+                      << ") for " << av.getNode() << endl;
+  Assert(d_trail[ci].d_eq.getPolynomial().getCoefficient(vl) == Constant::mkConstant(-1));
+
+  SumPair newFact = r + fresh_a;
+
+  TrailIndex nextIndex = d_trail.size();
+  d_trail.push_back(Constraint(newFact, d_trail[i].d_proof));
+
+  SubIndex subBy = d_subs.size();
+  d_subs.push_back(Substitution(freshNode, var, ci));
+
+  Debug("arith::dio") << "Decompose nextIndex " <<  d_trail[nextIndex].d_eq.getNode() << endl;
+  return make_pair(subBy, nextIndex);
+}
+
+
+DioSolver::TrailIndex DioSolver::applySubstitution(DioSolver::SubIndex si, DioSolver::TrailIndex ti){
+  Variable var = d_subs[si].d_eliminated;
+  TrailIndex subIndex = d_subs[si].d_constraint;
+
+  const SumPair& curr = d_trail[ti].d_eq;
+  Polynomial vsum = curr.getPolynomial();
+
+  Constant a = vsum.getCoefficient(VarList(var));
+  Assert(a.isIntegral());
+  if(!a.isZero()){
+    Integer one(1);
+    TrailIndex afterSub = combineEqAtIndexes(ti, one, subIndex, a.getValue().getNumerator());
+    Assert(d_trail[afterSub].d_eq.getPolynomial().getCoefficient(VarList(var)).isZero());
+    return afterSub;
+  }else{
+    return ti;
+  }
+}
+
+
+DioSolver::TrailIndex DioSolver::reduceByGCD(DioSolver::TrailIndex ti){
+  const SumPair& sp = d_trail[ti].d_eq;
+  Polynomial vsum = sp.getPolynomial();
+  Constant c = sp.getConstant();
+
+  Debug("arith::dio") << "reduceByGCD " << vsum.getNode() << endl;
+  Assert(!vsum.isConstant());
+  Integer g = vsum.gcd();
+  Assert(g >= 1);
+  Debug("arith::dio") << "gcd("<< vsum.getNode() <<")=" << g << " " << c.getValue() << endl;
+  if(g.divides(c.getValue().getNumerator())){
+    if(g > 1){
+      return scaleEqAtIndex(ti, g);
+    }else{
+      return ti;
+    }
+  }else{
+    raiseConflict(ti);
+    return ti;
+  }
+}
+
+bool DioSolver::triviallySat(TrailIndex i){
+  const SumPair& eq = d_trail[i].d_eq;
+  if(eq.isConstant()){
+    return eq.getConstant().isZero();
+  }else{
+    return false;
+  }
+}
+
+bool DioSolver::triviallyUnsat(DioSolver::TrailIndex i){
+  const SumPair& eq = d_trail[i].d_eq;
+  if(eq.isConstant()){
+    return !eq.getConstant().isZero();
+  }else{
+    return false;
+  }
+}
+
+
+bool DioSolver::gcdIsOne(DioSolver::TrailIndex i){
+  const SumPair& eq = d_trail[i].d_eq;
+  return eq.gcd() == Integer(1);
+}
+
+void DioSolver::debugPrintTrail(DioSolver::TrailIndex i) const{
+  const SumPair& eq = d_trail[i].d_eq;
+  const Polynomial& proof = d_trail[i].d_proof;
+
+  cout << "d_trail["<<i<<"].d_eq = " << eq.getNode() << endl;
+  cout << "d_trail["<<i<<"].d_proof = " << proof.getNode() << endl;
+}
+
+void DioSolver::subAndReduceCurrentFByIndex(DioSolver::SubIndex subIndex){
+  size_t N = d_currentF.size();
+
+  size_t readIter = 0, writeIter = 0;
+  for(; readIter < N && !inConflict(); ++readIter){
+    TrailIndex curr = d_currentF[readIter];
+    TrailIndex nextTI = applySubstitution(subIndex, curr);
+    if(nextTI == curr){
+      d_currentF[writeIter] = curr;
+      ++writeIter;
+    }else{
+      Assert(nextTI != curr);
+
+      if(triviallyUnsat(nextTI)){
+        raiseConflict(nextTI);
+      }else if(!(triviallySat(nextTI) || anyCoefficientExceedsMaximum(nextTI))){
+        TrailIndex nextNextTI = reduceByGCD(nextTI);
+
+        if(!inConflict()){
+          Assert(queueConditions(nextNextTI));
+          d_currentF[writeIter] = nextNextTI;
+          ++writeIter;
+        }
+      }
+    }
+  }
+  if(!inConflict() && writeIter < N){
+    d_currentF.resize(writeIter);
+  }
 }
 
 }/* CVC4::theory::arith namespace */
diff --git a/src/theory/arith/dio_solver.h b/src/theory/arith/dio_solver.h
index c91ddd994..da41d94cd 100644
--- a/src/theory/arith/dio_solver.h
+++ b/src/theory/arith/dio_solver.h
@@ -28,6 +28,8 @@
 #include "theory/arith/arith_utilities.h"
 #include "util/rational.h"
 
+#include "util/stats.h"
+
 #include <vector>
 #include <utility>
 
@@ -36,33 +38,384 @@ namespace theory {
 namespace arith {
 
 class DioSolver {
-  context::Context* d_context;
-  const Tableau& d_tableau;
-  ArithPartialModel& d_partialModel;
+private:
+  typedef size_t TrailIndex;
+  typedef size_t InputConstraintIndex;
+  typedef size_t SubIndex;
+
+  std::vector<Variable> d_variablePool;
+  context::CDO<size_t> d_lastUsedVariable;
+
+  /**
+   * The set of input constraints is stored in a CDList.
+   * Each constraint point to an element of the trail.
+   */
+  struct InputConstraint {
+    Node d_reason;
+    TrailIndex d_trailPos;
+    InputConstraint(Node reason, TrailIndex pos) : d_reason(reason), d_trailPos(pos) {}
+  };
+  context::CDList<InputConstraint> d_inputConstraints;
+
+  /**
+   * This is the next input constraint to handle.
+   */
+  context::CDO<size_t> d_nextInputConstraintToEnqueue;
+
+
+  /**
+   * We maintain a map from the variables associated with proofs to an input constraint.
+   * These variables can then be used in polynomial manipulations.
+   */
+  typedef std::hash_map<Node, InputConstraintIndex, NodeHashFunction> NodeToInputConstraintIndexMap;
+  NodeToInputConstraintIndexMap d_varToInputConstraintMap;
+
+  Node proofVariableToReason(const Variable& v) const{
+    Assert(d_varToInputConstraintMap.find(v.getNode()) != d_varToInputConstraintMap.end());
+    InputConstraintIndex pos = (*(d_varToInputConstraintMap.find(v.getNode()))).second;
+    Assert(pos < d_inputConstraints.size());
+    return d_inputConstraints[pos].d_reason;
+  }
+
+  /**
+   * The main work horse of the algorithm, the trail of constraints.
+   * Each constraint is a SumPair that implicitly represents an equality against 0.
+   *   d_trail[i].d_eq = (+ c (+ [(* coeff var)])) representing (+ [(* coeff var)]) = -c
+   * Each constraint has a proof in terms of a linear combination of the input constraints.
+   *   d_trail[i].d_proof
+   *
+   * Each Constraint also a monomial in d_eq.getPolynomial()
+   * of minimal absolute value by the coefficients.
+   * d_trail[i].d_minimalMonomial
+   *
+   * See Alberto's paper for how linear proofs are maintained for the abstract
+   * state machine in rules (7), (8) and (9).
+   */
+  struct Constraint {
+    SumPair d_eq;
+    Polynomial d_proof;
+    Monomial d_minimalMonomial;
+    Constraint(const SumPair& eq, const Polynomial& p) :
+      d_eq(eq), d_proof(p), d_minimalMonomial(d_eq.getPolynomial().selectAbsMinimum())
+    {}
+  };
+  context::CDList<Constraint> d_trail;
+
+  /** Compare by d_minimal. */
+  struct TrailMinimalCoefficientOrder {
+    const context::CDList<Constraint>& d_trail;
+    TrailMinimalCoefficientOrder(const context::CDList<Constraint>& trail):
+      d_trail(trail)
+    {}
+
+    bool operator()(TrailIndex i, TrailIndex j){
+      return d_trail[i].d_minimalMonomial.absLessThan(d_trail[j].d_minimalMonomial);
+    }
+  };
+
+  /**
+   * A substitution is stored as a constraint in the trail together with
+   * the variable to be eliminated, and a fresh variable if one was introduced.
+   * The variable d_subs[i].d_eliminated is substituted using the implicit equality in
+   * d_trail[d_subs[i].d_constraint]
+   *  - d_subs[i].d_eliminated is normalized to have coefficient -1 in
+   *    d_trail[d_subs[i].d_constraint].
+   *  - d_subs[i].d_fresh is either Node::null() or it is variable it is normalized
+   *    to have coefficient 1 in d_trail[d_subs[i].d_constraint].
+   */
+  struct Substitution {
+    Node d_fresh;
+    Variable d_eliminated;
+    TrailIndex d_constraint;
+    Substitution(Node f, const Variable& e, TrailIndex c) :
+      d_fresh(f), d_eliminated(e), d_constraint(c)
+    {}
+  };
+  context::CDList<Substitution> d_subs;
+
+  /**
+   * This is the queue of constraints to be processed in the current context level.
+   * This is to be empty upon entering solver and cleared upon leaving the solver.
+   *
+   * All elements in currentF:
+   * - are fully substituted according to d_subs.
+   * - !isConstant().
+   * - If the element is (+ constant (+ [(* coeff var)] )), then the gcd(coeff) = 1
+   */
+  std::deque<TrailIndex> d_currentF;
+  context::CDList<TrailIndex> d_savedQueue;
+  context::CDO<size_t> d_savedQueueIndex;
+
+  context::CDO<bool> d_conflictHasBeenRaised;
+  TrailIndex d_conflictIndex;
+
+  /**
+   * Drop derived constraints with a coefficient length larger than
+   * the maximum input constraints length than 2**MAX_GROWTH_RATE.
+   */
+  context::CDO<uint32_t> d_maxInputCoefficientLength;
+  static const uint32_t MAX_GROWTH_RATE = 3;
+
+  /** Returns true if the element on the trail should be dropped.*/
+  bool anyCoefficientExceedsMaximum(TrailIndex j) const;
+
+  /**
+   * Is true if decomposeIndex has been used in this context.
+   */
+  context::CDO<bool> d_usedDecomposeIndex;
+
+  context::CDO<SubIndex> d_lastPureSubstitution;
+  context::CDO<SubIndex> d_pureSubstitionIter;
 
 public:
 
   /** Construct a Diophantine equation solver with the given context. */
-  DioSolver(context::Context* ctxt, const Tableau& tab, ArithPartialModel& pmod) :
-    d_context(ctxt),
-    d_tableau(tab),
-    d_partialModel(pmod) {
+  DioSolver(context::Context* ctxt);
+
+  /** Returns true if the substitutions use no new variables. */
+  bool hasMorePureSubstitutions() const{
+    return d_pureSubstitionIter < d_lastPureSubstitution;
+  }
+
+  Node nextPureSubstitution();
+
+  /**
+   * Adds an equality to the queue of the DioSolver.
+   * orig is blamed in a conflict.
+   * orig can either be of the form (= p c) or (and ub lb).
+   * where ub is either (leq p c) or (not (> p [- c 1])), and
+   * where lb is either (geq p c) or (not (< p [+ c 1]))
+   */
+  void pushInputConstraint(const Comparison& eq, Node reason);
+
+  /**
+   * Processes the queue looking for any conflict.
+   * If a conflict is found, this returns conflict.
+   * Otherwise, it returns null.
+   * The conflict is guarenteed to be over literals given in addEquality.
+   */
+  Node processEquationsForConflict();
+
+  /**
+   * Processes the queue looking for an integer unsatisfiable cutting plane.
+   * If such a plane is found this returns an entailed plane using no
+   * fresh variables.
+   */
+  SumPair processEquationsForCut();
+
+private:
+  /** Returns true if the TrailIndex refers to a element in the trail. */
+  bool inRange(TrailIndex i) const{
+    return i < d_trail.size();
+  }
+
+  Node columnGcdIsOne() const;
+
+
+  /**
+   * Returns true if the context dependent flag for conflicts
+   * has been raised.
+   */
+  bool inConflict() const{
+    return d_conflictHasBeenRaised;
+  }
+
+  /** Raises a conflict at the index ti. */
+  void raiseConflict(TrailIndex ti){
+    Assert(!inConflict());
+    d_conflictHasBeenRaised = true;
+    d_conflictIndex = ti;
+  }
+
+  /** Returns the conflict index. */
+  TrailIndex getConflictIndex() const{
+    Assert(inConflict());
+    return d_conflictIndex;
+  }
+
+  /**
+   * Allocates a "unique" variables from the pool of integer variables.
+   * This variable is fresh with respect to the context.
+   * Returns index of the variable in d_variablePool;
+   */
+  size_t allocateVariableInPool();
+
+  /**
+   * Returns true if the node can be accepted as a reason according to the
+   * kinds.
+   */
+  bool acceptableOriginalNodes(Node n);
+
+  /** Empties the unproccessed input constraints into the queue. */
+  void enqueueInputConstraints();
+
+  /**
+   * Returns true if an input equality is in the map.
+   * This is expensive and is only for debug assertions.
+   */
+  bool debugEqualityInInputEquations(Node eq);
+
+  /** Applies the substitution at subIndex to currentF. */
+  void subAndReduceCurrentFByIndex(SubIndex d_subIndex);
+
+  /**
+   * Takes as input a TrailIndex i and an integer that divides d_trail[i].d_eq, and
+   * returns a TrailIndex j s.t.
+   *   d_trail[j].d_eq = (1/g) d_trail[i].d_eq
+   * and
+   *   d_trail[j].d_proof = (1/g) d_trail[i].d_proof.
+   *
+   * g must be non-zero.
+   *
+   * This corresponds to an application of Alberto's rule (7).
+   */
+  TrailIndex scaleEqAtIndex(TrailIndex i, const Integer& g);
+
+
+  /**
+   * Takes as input TrailIndex's i and j and Integer's q and r and a TrailIndex k s.t.
+   *   d_trail[k].d_eq == d_trail[i].d_eq * q + d_trail[j].d_eq * r
+   * and
+   *   d_trail[k].d_proof == d_trail[i].d_proof * q + d_trail[j].d_proof * r
+   *
+   * This corresponds to an application of Alberto's rule (8).
+   */
+  TrailIndex combineEqAtIndexes(TrailIndex i, const Integer& q, TrailIndex j, const Integer& r);
+
+  /**
+   * Decomposes the equation at index ti of trail by the variable
+   * with the lowest coefficient.
+   * This corresponds to an application of Alberto's rule (9).
+   *
+   * Returns a pair of a SubIndex and a TrailIndex.
+   * The SubIndex is the index of a newly introduced substition.
+   */
+  std::pair<SubIndex, TrailIndex> decomposeIndex(TrailIndex ti);
+
+  /** Solves the index at ti for the value in minimumMonomial. */
+  std::pair<SubIndex, TrailIndex> solveIndex(TrailIndex ti);
+
+  /** Prints the queue for debugging purposes to Debug("arith::dio"). */
+  void printQueue();
+
+  /**
+   * Exhaustively applies all substitutions discovered to an element of the trail.
+   * Returns a TrailIndex corresponding to the substitutions being applied.
+   */
+  TrailIndex applyAllSubstitutionsToIndex(TrailIndex i);
+
+  /**
+   * Applies a substitution to an element in the trail.
+   */
+  TrailIndex applySubstitution(SubIndex s, TrailIndex i);
+
+  /**
+   * Reduces the trail node at i by the gcd of the variables.
+   * Returns the new trail element.
+   *
+   * This raises the conflict flag if unsat is detected.
+   */
+  TrailIndex reduceByGCD(TrailIndex i);
+
+  /**
+   * Returns true if i'th element in the trail is trivially true.
+   * (0 = 0)
+   */
+  bool triviallySat(TrailIndex t);
+
+  /**
+   * Returns true if i'th element in the trail is trivially unsatisfiable.
+   * (1 = 0)
+   */
+  bool triviallyUnsat(TrailIndex t);
+
+  /** Returns true if the gcd of the i'th element of the trail is 1.*/
+  bool gcdIsOne(TrailIndex t);
+
+  bool debugAnySubstitionApplies(TrailIndex t);
+  bool debugSubstitutionApplies(SubIndex si, TrailIndex ti);
+
+
+  /** Returns true if the queue of nodes to process is empty. */
+  bool queueEmpty() const;
+
+  bool queueConditions(TrailIndex t){
+    /* debugPrintTrail(t); */
+    
+    /* std::cout << !inConflict() << std::endl; */
+    /* std::cout << gcdIsOne(t) << std::endl; */
+    /* std::cout << !debugAnySubstitionApplies(t) << std::endl; */
+    /* std::cout << !triviallySat(t) << std::endl; */
+    /* std::cout << !triviallyUnsat(t) << std::endl; */
+
+    return
+      !inConflict() &&
+      gcdIsOne(t) &&
+      !debugAnySubstitionApplies(t) &&
+      !triviallySat(t) &&
+      !triviallyUnsat(t);
+  }
+
+  void pushToQueueBack(TrailIndex t){
+    Assert(queueConditions(t));
+    d_currentF.push_back(t);
+  }
+
+  void pushToQueueFront(TrailIndex t){
+    Assert(queueConditions(t));
+    d_currentF.push_front(t);
   }
 
   /**
-   * Solve the set of Diophantine equations in the tableau.
+   * Moves the minimum Constraint by absolute value of the minimum coefficient to
+   * the front of the queue.
+   */
+  void moveMinimumByAbsToQueueFront();
+
+  void saveQueue();
+
+  TrailIndex impliedGcdOfOne();
+
+
+  /**
+   * Processing the current set of equations.
    *
-   * @return true if the set of equations was solved; false if no
-   * solution
+   * decomposeIndex() rule is only applied if allowDecomposition is true.
    */
-  bool solve();
+  bool processEquations(bool allowDecomposition);
 
   /**
-   * Get the parameterized solution to this set of Diophantine
-   * equations.  Must be preceded by a call to solve() that returned
-   * true. */
-  void getSolution();
+   * Constructs a proof from any d_trail[i] in terms of input literals.
+   */
+  Node proveIndex(TrailIndex i);
+
+  /**
+   * Returns the SumPair in d_trail[i].d_eq with all of the fresh variables purified out.
+   */
+  SumPair purifyIndex(TrailIndex i);
+
+
+  void debugPrintTrail(TrailIndex i) const;
+
+public:
+  /** These fields are designed to be accessable to TheoryArith methods. */
+  class Statistics {
+  public:
+
+    IntStat d_conflictCalls;
+    IntStat d_cutCalls;
+
+    IntStat d_cuts;
+    IntStat d_conflicts;
+
+    TimerStat d_conflictTimer;
+    TimerStat d_cutTimer;
+
+    Statistics();
+    ~Statistics();
+  };
 
+  Statistics d_statistics;
 };/* class DioSolver */
 
 }/* CVC4::theory::arith namespace */
diff --git a/src/theory/arith/normal_form.cpp b/src/theory/arith/normal_form.cpp
index e7df14df7..a4dc78c9f 100644
--- a/src/theory/arith/normal_form.cpp
+++ b/src/theory/arith/normal_form.cpp
@@ -21,9 +21,10 @@
 #include <list>
 
 using namespace std;
-using namespace CVC4;
-using namespace CVC4::theory;
-using namespace CVC4::theory::arith;
+
+namespace CVC4 {
+namespace theory{
+namespace arith {
 
 bool VarList::isSorted(iterator start, iterator end) {
   return __gnu_cxx::is_sorted(start, end);
@@ -131,9 +132,6 @@ Monomial Monomial::operator*(const Monomial& mono) const {
 
 vector<Monomial> Monomial::sumLikeTerms(const std::vector<Monomial> & monos) {
   Assert(isSorted(monos));
-
-  Debug("blah") << "start sumLikeTerms" << std::endl;
-  printList(monos);
   vector<Monomial> outMonomials;
   typedef vector<Monomial>::const_iterator iterator;
   for(iterator rangeIter = monos.begin(), end=monos.end(); rangeIter != end;) {
@@ -150,9 +148,6 @@ vector<Monomial> Monomial::sumLikeTerms(const std::vector<Monomial> & monos) {
       outMonomials.push_back(nonZero);
     }
   }
-  Debug("blah") << "outmonomials" << std::endl;
-  printList(monos);
-  Debug("blah") << "end sumLikeTerms" << std::endl;
 
   Assert(isStrictlySorted(outMonomials));
   return outMonomials;
@@ -169,8 +164,6 @@ void Monomial::printList(const std::vector<Monomial>& list) {
   }
 }
 Polynomial Polynomial::operator+(const Polynomial& vl) const {
-  this->printList();
-  vl.printList();
 
   std::vector<Monomial> sortedMonos;
   merge_ranges(begin(), end(), vl.begin(), vl.end(), sortedMonos);
@@ -178,7 +171,6 @@ Polynomial Polynomial::operator+(const Polynomial& vl) const {
   std::vector<Monomial> combined = Monomial::sumLikeTerms(sortedMonos);
 
   Polynomial result = mkPolynomial(combined);
-  result.printList();
   return result;
 }
 
@@ -211,9 +203,47 @@ Polynomial Polynomial::operator*(const Polynomial& poly) const {
   return res;
 }
 
+Monomial Polynomial::selectAbsMinimum() const {
+  iterator iter = begin(), myend = end();
+  Assert(iter != myend);
+
+  Monomial min = *iter;
+  ++iter;
+  for(; iter != end(); ++iter){
+    Monomial curr = *iter;
+    if(curr.absLessThan(min)){
+      min = curr;
+    }
+  }
+  return min;
+}
+
+Integer Polynomial::gcd() const {
+  //We'll use the standardization that gcd(0, 0) = 0
+  //So that the gcd of the zero polynomial is gcd{0} = 0
+  Assert(isIntegral());
+  iterator i=begin(), e=end();
+  Assert(i!=e);
+
+  Integer d = (*i).getConstant().getValue().getNumerator().abs();
+  ++i;
+  for(; i!=e; ++i){
+    Integer c = (*i).getConstant().getValue().getNumerator();
+    d = d.gcd(c);
+  }
+  return d;
+}
+
+Integer Polynomial::denominatorLCM() const {
+  Integer tmp(1);
+  for(iterator i=begin(), e=end(); i!=e; ++i){
+    const Constant& c = (*i).getConstant();
+    tmp = tmp.lcm(c.getValue().getDenominator());
+  }
+  return tmp;
+}
 
 Node Comparison::toNode(Kind k, const Polynomial& l, const Constant& r) {
-  Assert(!l.isConstant());
   Assert(isRelationOperator(k));
   switch(k) {
   case kind::GEQ:
@@ -306,22 +336,27 @@ bool Comparison::pbComparison(Kind k, TNode left, const Rational& right, bool& r
   return false;
 }
 
+// Comparison Comparison::mkComparison(Kind k, const Polynomial& left, const Constant& right) {
+//   Assert(isRelationOperator(k));
+//   if(left.isConstant()) {
+//     const Rational& lConst =  left.getNode().getConst<Rational>();
+//     const Rational& rConst = right.getNode().getConst<Rational>();
+//     bool res = evaluateConstantPredicate(k, lConst, rConst);
+//     return Comparison(res);
+//   }
+
+//   if(left.getNode().getType().isPseudoboolean()) {
+//     bool result;
+//     if(pbComparison(k, left.getNode(), right.getNode().getConst<Rational>(), result)) {
+//       return Comparison(result);
+//     }
+//   }
+
+//   return Comparison(toNode(k, left, right), k, left, right);
+// }
+
 Comparison Comparison::mkComparison(Kind k, const Polynomial& left, const Constant& right) {
   Assert(isRelationOperator(k));
-  if(left.isConstant()) {
-    const Rational& lConst =  left.getNode().getConst<Rational>();
-    const Rational& rConst = right.getNode().getConst<Rational>();
-    bool res = evaluateConstantPredicate(k, lConst, rConst);
-    return Comparison(res);
-  }
-
-  if(left.getNode().getType().isPseudoboolean()) {
-    bool result;
-    if(pbComparison(k, left.getNode(), right.getNode().getConst<Rational>(), result)) {
-      return Comparison(result);
-    }
-  }
-
   return Comparison(toNode(k, left, right), k, left, right);
 }
 
@@ -334,9 +369,231 @@ Comparison Comparison::addConstant(const Constant& constant) const {
   return mkComparison(oper, newLeft, newRight);
 }
 
+bool Comparison::constantInLefthand() const{
+  return getLeft().containsConstant();
+}
+
+Comparison Comparison::cancelLefthandConstant() const {
+  if(constantInLefthand()){
+    Monomial constantHead = getLeft().getHead();
+    Assert(constantHead.isConstant());
+
+    Constant constant = constantHead.getConstant();
+    Constant negativeConstantHead = -constant;
+    return addConstant(negativeConstantHead);
+  }else{
+    return *this;
+  }
+}
+
 Comparison Comparison::multiplyConstant(const Constant& constant) const {
   Assert(!isBoolean());
   Kind newOper = (constant.getValue() < 0) ? reverseRelationKind(oper) : oper;
 
   return mkComparison(newOper, left*Monomial(constant), right*constant);
 }
+
+Integer Comparison::denominatorLCM() const {
+  // Get the coefficients to be all integers.
+  Integer leftDenominatorLCM = left.denominatorLCM();
+  Integer rightDenominator = right.getValue().getDenominator();
+  Integer denominatorLCM = leftDenominatorLCM.lcm(rightDenominator);
+  Assert(denominatorLCM.sgn() > 0);
+  return denominatorLCM;
+}
+
+Comparison Comparison::multiplyByDenominatorLCM() const{
+  return multiplyConstant(Constant::mkConstant(denominatorLCM()));
+}
+
+Comparison Comparison::normalizeLeadingCoefficientPositive() const {
+  if(getLeft().getHead().getConstant().isNegative()){
+    return multiplyConstant(Constant::mkConstant(-1));
+  }else{
+    return *this;
+  }
+}
+
+bool Comparison::isIntegral() const {
+  return getRight().isIntegral() && getLeft().isIntegral();
+}
+
+bool Comparison::isConstant() const {
+  return getLeft().isConstant();
+}
+
+Comparison Comparison::evaluateConstant() const {
+  Assert(left.isConstant());
+  const Rational& rConst = getRight().getValue();
+  const Rational& lConst = getLeft().getHead().getConstant().getValue();
+  bool res = evaluateConstantPredicate(getKind(), lConst, rConst);
+  return Comparison(res);
+}
+
+Comparison Comparison::divideByLefthandGCD() const {
+  Assert(isIntegral());
+
+  Integer zr = getRight().getValue().getNumerator();
+  Integer gcd = getLeft().gcd();
+  Polynomial newLeft = getLeft().exactDivide(gcd);
+
+  Constant newRight = Constant::mkConstant(Rational(zr,gcd));
+  return mkComparison(getKind(), newLeft, newRight);
+}
+
+Comparison Comparison::divideByLeadingCoefficient() const {
+  //Handle the rational/mixed case
+  Monomial head = getLeft().getHead();
+  Constant leadingCoefficient = head.getConstant();
+  Assert(!leadingCoefficient.isZero());
+
+  Constant inverse = leadingCoefficient.inverse();
+
+  return multiplyConstant(inverse);
+}
+
+Comparison Comparison::tightenIntegralConstraint() const {
+  Assert(getLeft().isIntegral());
+
+  if(getRight().isIntegral()){
+    return *this;
+  }else{
+    switch(getKind()){
+    case kind::EQUAL:
+      //If the gcd of the left hand side does not cleanly divide the right hand side,
+      //this is unsatisfiable in the theory of Integers.
+      return Comparison(false);
+    case kind::GEQ:
+    case kind::GT:
+      {
+        //(>= l (/ n d))
+        //(>= l (ceil (/ n d)))
+        //This also hold for GT as (ceil (/ n d)) > (/ n d)
+        Integer ceilr = getRight().getValue().ceiling();
+        Constant newRight = Constant::mkConstant(ceilr);
+        return Comparison(toNode(kind::GEQ, getLeft(), newRight),kind::GEQ, getLeft(),newRight);
+      }
+    case kind::LEQ:
+    case kind::LT:
+      {
+        //(<= l (/ n d))
+        //(<= l (floor (/ n d)))
+        //This also hold for LT as (floor (/ n d)) < (/ n d)
+        Integer floor = getRight().getValue().floor();
+        Constant newRight = Constant::mkConstant(floor);
+        return Comparison(toNode(kind::LEQ, getLeft(), newRight),kind::LEQ, getLeft(),newRight);
+      }
+    default:
+      Unreachable();
+    }
+  }
+}
+
+bool Comparison::isIntegerNormalForm() const{
+  if(constantInLefthand()){ return false; }
+  else if(getLeft().getHead().getConstant().isNegative()){ return false; }
+  else if(!isIntegral()){ return false; }
+  else {
+    return getLeft().gcd() == 1;
+  }
+}
+bool Comparison::isMixedNormalForm() const {
+  if(constantInLefthand()){ return false; }
+  else if(allIntegralVariables()) { return false; }
+  else{
+    return getLeft().getHead().getConstant().getValue() == 1;
+  }
+}
+
+Comparison Comparison::normalize(Comparison c) {
+  if(c.isConstant()){
+    return c.evaluateConstant();
+  }else{
+    Comparison c0 = c.constantInLefthand() ? c.cancelLefthandConstant() : c;
+    Comparison c1 = c0.normalizeLeadingCoefficientPositive();
+    if(c1.allIntegralVariables()){
+      //All Integer Variable Case
+      Comparison integer0 = c1.multiplyByDenominatorLCM();
+      Comparison integer1 = integer0.divideByLefthandGCD();
+      Comparison integer2 = integer1.tightenIntegralConstraint();
+      Assert(integer2.isBoolean() || integer2.isIntegerNormalForm());
+      return integer2;
+    }else{
+      //Mixed case
+      Comparison mixed = c1.divideByLeadingCoefficient();
+      Assert(mixed.isMixedNormalForm());
+      return mixed;
+    }
+  }
+}
+
+Comparison Comparison::mkNormalComparison(Kind k, const Polynomial& left, const Constant& right) {
+  Comparison cmp = mkComparison(k,left,right);
+  Comparison normalized = cmp.normalize(cmp);
+  Assert(normalized.isNormalForm());
+  return normalized;
+}
+
+Node Polynomial::computeQR(const Polynomial& p, const Integer& div){
+  Assert(p.isIntegral());
+  std::vector<Monomial> q_vec, r_vec;
+  Integer tmp_q, tmp_r;
+  for(iterator iter = p.begin(), pend = p.end(); iter != pend; ++iter){
+    Monomial curr = *iter;
+    VarList vl = curr.getVarList();
+    Constant c = curr.getConstant();
+
+    const Integer& a = c.getValue().getNumerator();
+    Integer::floorQR(tmp_q, tmp_r, a, div);
+    Constant q=Constant::mkConstant(tmp_q);
+    Constant r=Constant::mkConstant(tmp_r);
+    if(!q.isZero()){
+      q_vec.push_back(Monomial::mkMonomial(q, vl));
+    }
+    if(!r.isZero()){
+      r_vec.push_back(Monomial::mkMonomial(r, vl));
+    }
+  }
+
+  Polynomial p_q = Polynomial::mkPolynomial(q_vec);
+  Polynomial p_r = Polynomial::mkPolynomial(r_vec);
+
+  return NodeManager::currentNM()->mkNode(kind::PLUS, p_q.getNode(), p_r.getNode());
+}
+
+Node SumPair::computeQR(const SumPair& sp, const Integer& div){
+  Assert(sp.isIntegral());
+
+  const Integer& constant = sp.getConstant().getValue().getNumerator();
+
+  Integer constant_q, constant_r;
+  Integer::floorQR(constant_q, constant_r, constant, div);
+
+  Node p_qr = Polynomial::computeQR(sp.getPolynomial(), div);
+  Assert(p_qr.getKind() == kind::PLUS);
+  Assert(p_qr.getNumChildren() == 2);
+
+  Polynomial p_q = Polynomial::parsePolynomial(p_qr[0]);
+  Polynomial p_r = Polynomial::parsePolynomial(p_qr[1]);
+
+  SumPair sp_q(p_q, Constant::mkConstant(constant_q));
+  SumPair sp_r(p_r, Constant::mkConstant(constant_r));
+
+  return NodeManager::currentNM()->mkNode(kind::PLUS, sp_q.getNode(), sp_r.getNode());
+}
+
+Constant Polynomial::getCoefficient(const VarList& vl) const{
+  //TODO improve to binary search...
+  for(iterator iter=begin(), myend=end(); iter != myend; ++iter){
+    Monomial m = *iter;
+    VarList curr = m.getVarList();
+    if(curr == vl){
+      return m.getConstant();
+    }
+  }
+  return Constant::mkConstant(0);
+}
+
+} //namespace arith
+} //namespace theory
+} //namespace CVC4
diff --git a/src/theory/arith/normal_form.h b/src/theory/arith/normal_form.h
index a1a33fed0..71d7c96f4 100644
--- a/src/theory/arith/normal_form.h
+++ b/src/theory/arith/normal_form.h
@@ -50,6 +50,7 @@ namespace arith {
  * variable := n
  *   where
  *     n.getMetaKind() == metakind::VARIABLE
+ *     n.getType() \in {Integer, Real}
  *
  * constant := n
  *   where
@@ -70,13 +71,24 @@ namespace arith {
  *     isStrictlySorted monoOrder [monomial]
  *     forall (\x -> x != 0) [monomial]
  *
- * restricted_cmp := (|><| polynomial constant)
+ * rational_restricted_cmp := (|><| qpolynomial constant)
  *   where
  *     |><| is GEQ, EQ, or EQ
- *     not (exists constantMonomial (monomialList polynomial))
- *     monomialCoefficient (head (monomialList polynomial)) == 1
+ *     not (exists constantMonomial (monomialList qpolynomial))
+ *     (exists realMonomial (monomialList qpolynomial))
+ *     monomialCoefficient (head (monomialList qpolynomial)) == 1
  *
- * comparison := TRUE | FALSE | restricted_cmp | (not restricted_cmp)
+ * integer_restricted_cmp := (|><| zpolynomial constant)
+ *   where
+ *     |><| is GEQ, EQ, or EQ
+ *     not (exists constantMonomial (monomialList zpolynomial))
+ *     (forall integerMonomial (monomialList qpolynomial))
+ *     the denominator of all coefficients and the constant is 1
+ *     the gcd of all numerators of coefficients is 1
+ *
+ * comparison := TRUE | FALSE
+ *   | rational_restricted_cmp | (not rational_restricted_cmp)
+ *   | integer_restricted_cmp | (not integer_restricted_cmp)
  *
  * Normal Form for terms := polynomial
  * Normal Form for atoms := comparison
@@ -146,6 +158,11 @@ namespace arith {
  *
  *  monoOrder m0 m1 = var_listOrder (monomialVarList m0) (monomialVarList m1)
  *
+ *  integerMonomial mono =
+ *    forall varHasTypeInteger (monomialVarList mono)
+ *
+ *  realMonomial mono = not (integerMonomial mono)
+ *
  *  constantMonomial monomial =
  *    match monomial with
  *        constant -> true
@@ -193,6 +210,10 @@ public:
 
   bool isNormalForm() { return isMember(getNode()); }
 
+  bool isIntegral() const {
+    return getNode().getType().isInteger();
+  }
+
   bool operator<(const Variable& v) const { return getNode() < v.getNode();}
   bool operator==(const Variable& v) const { return getNode() == v.getNode();}
 
@@ -223,7 +244,14 @@ public:
     return getNode().getConst<Rational>();
   }
 
-  bool isZero() const { return getValue() == 0; }
+  bool isIntegral() const { return getValue().isIntegral(); }
+
+  int sgn() const { return getValue().sgn(); }
+
+  bool isZero() const { return sgn() == 0; }
+  bool isNegative() const { return sgn() < 0; }
+  bool isPositive() const { return sgn() > 0; }
+
   bool isOne() const { return getValue() == 1; }
 
   Constant operator*(const Constant& other) const {
@@ -236,6 +264,33 @@ public:
     return mkConstant(-getValue());
   }
 
+  Constant inverse() const{
+    Assert(!isZero());
+    return mkConstant(getValue().inverse());
+  }
+
+  bool operator<(const Constant& other) const {
+    return getValue() < other.getValue();
+  }
+
+  bool operator==(const Constant& other) const {
+    //Rely on node uniqueness.
+    return getNode() == other.getNode();
+  }
+
+  Constant abs() const {
+    if(isNegative()){
+      return -(*this);
+    }else{
+      return (*this);
+    }
+  }
+
+  uint32_t length() const{
+    Assert(isIntegral());
+    return getValue().getNumerator().length();
+  }
+
 };/* class Constant */
 
 
@@ -359,6 +414,11 @@ public:
     return iterator(internalEnd());
   }
 
+  Variable getHead() const {
+    Assert(!empty());
+    return *(begin());
+  }
+
   VarList(Variable v) : NodeWrapper(v.getNode()) {
     Assert(isSorted(begin(), end()));
   }
@@ -412,6 +472,16 @@ public:
 
   bool operator==(const VarList& vl) const { return cmp(vl) == 0; }
 
+  bool isIntegral() const {
+    for(iterator i = begin(), e=end(); i != e; ++i ){
+      Variable var = *i;
+      if(!var.isIntegral()){
+        return false;
+      }
+    }
+    return true;
+  }
+
 private:
   bool isSorted(iterator start, iterator end);
 
@@ -471,6 +541,9 @@ public:
   /** Makes a monomial with no restrictions on c and vl. */
   static Monomial mkMonomial(const Constant& c, const VarList& vl);
 
+  static Monomial mkMonomial(const Variable& v){
+    return Monomial(VarList(v));
+  }
 
   static Monomial parseMonomial(Node n);
 
@@ -495,6 +568,10 @@ public:
     return constant.isOne();
   }
 
+  bool absCoefficientIsOne() const {
+    return coefficientIsOne() || constant.getValue() == -1;
+  }
+
   Monomial operator*(const Monomial& mono) const;
 
 
@@ -519,11 +596,40 @@ public:
   }
 
   /**
+   * The variable product
+   */
+  bool integralVariables() const {
+    return getVarList().isIntegral();
+  }
+
+  /**
+   * The coefficient of the monomial is integral.
+   */
+  bool integralCoefficient() const {
+    return getConstant().isIntegral();
+  }
+
+  /**
+   * A Monomial is an "integral" monomial if the constant is integral.
+   */
+  bool isIntegral() const {
+    return integralCoefficient() && integralVariables();
+  }
+
+  /**
    * Given a sorted list of monomials, this function transforms this
    * into a strictly sorted list of monomials that does not contain zero.
    */
   static std::vector<Monomial> sumLikeTerms(const std::vector<Monomial>& monos);
 
+  bool absLessThan(const Monomial& other) const{
+    return getConstant().abs() < other.getConstant().abs();
+  }
+
+  uint32_t coefficientLength() const{
+    return getConstant().length();
+  }
+
   void print() const;
   static void printList(const std::vector<Monomial>& list);
 
@@ -638,6 +744,9 @@ public:
     Assert( Monomial::isStrictlySorted(m) );
   }
 
+  static Polynomial mkPolynomial(const Variable& v){
+    return Monomial::mkMonomial(v);
+  }
 
   static Polynomial mkPolynomial(const std::vector<Monomial>& m) {
     if(m.size() == 0) {
@@ -696,12 +805,103 @@ public:
     }
   }
 
+  /** A Polynomial is an "integral" polynomial if all of the monomials are integral. */
+  bool allIntegralVariables() const {
+    for(iterator i = begin(), e=end(); i!=e; ++i){
+      if(!(*i).integralVariables()){
+        return false;
+      }
+    }
+    return true;
+  }
+
+  /**
+   * A Polynomial is an "integral" polynomial if all of the monomials are integral
+   * and all of the coefficients are Integral. */
+  bool isIntegral() const {
+    for(iterator i = begin(), e=end(); i!=e; ++i){
+      if(!(*i).isIntegral()){
+        return false;
+      }
+    }
+    return true;
+  }
+
+  /**
+   * Selects a minimal monomial in the polynomial by the absolute value of
+   * the coefficient.
+   */
+  Monomial selectAbsMinimum() const;
+
+  /**
+   * Returns the Least Common Multiple of the denominators of the coefficients
+   * of the monomials.
+   */
+  Integer denominatorLCM() const;
+
+  /**
+   * Returns the GCD of the coefficients of the monomials.
+   * Requires this to be an isIntegral() polynomial.
+   */
+  Integer gcd() const;
+
+  Polynomial exactDivide(const Integer& z) const {
+    Assert(isIntegral());
+    Constant invz = Constant::mkConstant(Rational(1,z));
+    Polynomial prod = (*this) * Monomial(invz);
+    Assert(prod.isIntegral());
+    return prod;
+  }
+
   Polynomial operator+(const Polynomial& vl) const;
 
   Polynomial operator*(const Monomial& mono) const;
 
   Polynomial operator*(const Polynomial& poly) const;
 
+  /**
+   * Viewing the integer polynomial as a list [(* coeff_i mono_i)]
+   * The quotient and remainder of p divided by the non-zero integer z is:
+   *   q := [(* floor(coeff_i/z) mono_i )]
+   *   r := [(* rem(coeff_i/z) mono_i)]
+   * computeQR(p,z) returns the node (+ q r).
+   *
+   * q and r are members of the Polynomial class.
+   * For example:
+   * computeQR( p = (+ 5 (* 3 x) (* 8 y)) , z = 2) returns
+   *   (+ (+ 2 x (* 4 y)) (+ 1 x))
+   */
+  static Node computeQR(const Polynomial& p, const Integer& z);
+
+  /** Returns the coefficient assiociated with the VarList in the polynomial. */
+  Constant getCoefficient(const VarList& vl) const;
+
+  uint32_t maxLength() const{
+    iterator i = begin(), e=end();
+    if( i == e){
+      return 1;
+    }else{
+      uint32_t max = (*i).coefficientLength();
+      ++i;
+      for(; i!=e; ++i){      
+        uint32_t curr = (*i).coefficientLength();
+        if(curr > max){
+          max = curr;
+        }
+      }
+      return max;
+    }
+  }
+
+  uint32_t numMonomials() const {
+    if( getNode().getKind() == kind::PLUS ){
+      return getNode().getNumChildren();
+    }else if(isZero()){
+      return 0;
+    }else{
+      return 1;
+    }
+  }
 };/* class Polynomial */
 
 
@@ -732,6 +932,8 @@ private:
    */
   static bool pbComparison(Kind k, TNode left, const Rational& right, bool& result);
 
+  Kind getKind() const { return oper; }
+
 public:
   Comparison(bool val) :
     NodeWrapper(NodeManager::currentNM()->mkConst(val)),
@@ -750,30 +952,112 @@ public:
 
   static Comparison mkComparison(Kind k, const Polynomial& left, const Constant& right);
 
+  /**
+   * Returns the left hand side of the comparison.
+   * This is a polynomial that always contains all of the variables.
+   */
+  const Polynomial& getLeft() const { return left; }
+
+  /**
+   * Returns the right hand constatn of the comparison.
+   * If in normal form, this is the only constant term.
+   */
+  const Constant& getRight() const { return right; }
+
+  /** Returns true if the comparison is a boolean constant. */
   bool isBoolean() const {
     return (oper == kind::CONST_BOOLEAN);
   }
 
+  /** Returns true if all of the variables are integers. */
+  bool allIntegralVariables() const {
+    return getLeft().allIntegralVariables();
+  }
+
+  /**
+   * Returns true if the comparison is either a boolean term,
+   * in integer normal form or mixed normal form.
+   */
   bool isNormalForm() const {
     if(isBoolean()) {
       return true;
-    } else if(left.containsConstant()) {
-      return false;
-    } else if(left.getHead().getConstant().isOne()) {
-      return true;
-    } else {
-      return false;
+    } else if(allIntegralVariables()) {
+      return isIntegerNormalForm();
+    } else{
+      return isMixedNormalForm();
     }
   }
 
-  const Polynomial& getLeft() const { return left; }
-  const Constant& getRight() const { return right; }
+private:
+  /** Normal form check if at least one variable is real. */
+  bool isMixedNormalForm() const;
+
+  /** Normal form check is all variables are real.*/
+  bool isIntegerNormalForm() const;
+
+public:
+  /**
+   * Returns true if the left hand side is the sum of a non-zero constant
+   * and a polynomial.*/
+  bool constantInLefthand() const;
+
+  /**
+   * Returns a polynomial that is equivalent to the original but does not contain
+   * a constant in the top level sum on the left hand side.
+   */
+  Comparison cancelLefthandConstant() const;
+
+
+  /** Returns true if the polynomial is a constant. */
+  bool isConstant() const;
+
+  /** Reduces a constant comparison to a boolean comaprison.*/
+  Comparison evaluateConstant() const;
+
+
+  /**
+   * Returns true if all of the variables are integers, the coefficients are integers,
+   * and the right hand coefficient is an integer.
+   */
+  bool isIntegral() const;
+
+  /**
+   * Returns the Least Common Multiple of the monomials
+   * on the lefthand side and the constant on the right.
+   */
+  Integer denominatorLCM() const;
+
+  /** Multiplies the comparison by the denominatorLCM(). */
+  Comparison multiplyByDenominatorLCM() const;
+
+  /** If the leading coefficient is negative, multiply by -1. */
+  Comparison normalizeLeadingCoefficientPositive() const;
+
+  /** Divides the Comaprison by the gcd of the lefthand.*/
+  Comparison divideByLefthandGCD() const;
+
+  /** Divides the Comparison by the leading coefficient. */
+  Comparison divideByLeadingCoefficient() const;
+
+  /**
+   * If the left hand is integral and the right hand is not,
+   * the comparison is tightened to the nearest integer.
+   */
+  Comparison tightenIntegralConstraint() const;
 
   Comparison addConstant(const Constant& constant) const;
+
+  /* Multiply a non-boolean comparison by a constant term. */
   Comparison multiplyConstant(const Constant& constant) const;
 
   static Comparison parseNormalForm(TNode n);
 
+  /* Returns a logically equivalent comparison in normal form. */
+  static Comparison normalize(Comparison c);
+
+  /** Makes a comparison that is in normal form. */
+  static Comparison mkNormalComparison(Kind k, const Polynomial& left, const Constant& right);
+
   inline static bool isNormalAtom(TNode n){
     Comparison parse = Comparison::parseNormalForm(n);
     return parse.isNormalForm();
@@ -781,6 +1065,132 @@ public:
 
 };/* class Comparison */
 
+/**
+ * SumPair is a utility class that extends polynomials for use in computations.
+ * A SumPair is always a combination of (+ p c) where
+ *  c is a constant and p is a polynomial such that p = 0 or !p.containsConstant().
+ *
+ * These are a useful utility for representing the equation p = c as (+ p -c) where the pair
+ * is known to implicitly be equal to 0.
+ *
+ * SumPairs do not have unique representations due to the potential for p = 0.
+ * This makes them inappropraite for normal forms.
+ */
+class SumPair : public NodeWrapper {
+private:
+  static Node toNode(const Polynomial& p, const Constant& c){
+    return NodeManager::currentNM()->mkNode(kind::PLUS, p.getNode(), c.getNode());
+  }
+
+  SumPair(TNode n) :
+    NodeWrapper(n)
+  {
+    Assert(isNormalForm());
+  }
+
+public:
+
+  SumPair(const Polynomial& p):
+    NodeWrapper(toNode(p, Constant::mkConstant(0)))
+  {
+    Assert(isNormalForm());
+  }
+
+  SumPair(const Polynomial& p, const Constant& c):
+    NodeWrapper(toNode(p, c))
+  {
+    Assert(isNormalForm());
+  }
+
+  static bool isMember(TNode n) {
+    if(n.getKind() == kind::PLUS && n.getNumChildren() == 2){
+      if(Constant::isMember(n[1])){
+        if(Polynomial::isMember(n[0])){
+          Polynomial p = Polynomial::parsePolynomial(n[0]);
+          return p.isZero() || (!p.containsConstant());
+        }else{
+          return false;
+        }
+      }else{
+        return false;
+      }
+    }else{
+      return false;
+    }
+  }
+
+  bool isNormalForm() const {
+    return isMember(getNode());
+  }
+
+  Polynomial getPolynomial() const {
+    return Polynomial::parsePolynomial(getNode()[0]);
+  }
+
+  Constant getConstant() const {
+    return Constant::mkConstant((getNode())[1]);
+  }
+
+  SumPair operator+(const SumPair& other) const {
+    return SumPair(getPolynomial() + other.getPolynomial(),
+                   getConstant() + other.getConstant());
+  }
+
+  SumPair operator*(const Constant& c) const {
+    return SumPair(getPolynomial() * c, getConstant() * c);
+  }
+
+  SumPair operator-(const SumPair& other) const {
+    return (*this) + (other * Constant::mkConstant(-1));
+  }
+
+  static SumPair mkSumPair(const Variable& var){
+    return SumPair(Polynomial::mkPolynomial(var));
+  }
+
+  static SumPair parseSumPair(TNode n){
+    return SumPair(n);
+  }
+
+  bool isIntegral() const{
+    return getConstant().isIntegral() && getPolynomial().isIntegral();
+  }
+
+  bool isConstant() const {
+    return getPolynomial().isZero();
+  }
+
+  bool isZero() const {
+    return getConstant().isZero() && isConstant();
+  }
+
+  /**
+   * Returns the greatest common divisor of gcd(getPolynomial()) and getConstant().
+   * The SumPair must be integral.
+   */
+  Integer gcd() const {
+    Assert(isIntegral());
+    return (getPolynomial().gcd()).gcd(getConstant().getValue().getNumerator());
+  }
+
+  uint32_t maxLength() const {
+    Assert(isIntegral());
+    return std::max(getPolynomial().maxLength(), getConstant().length());
+  }
+
+  static SumPair mkZero() {
+    return SumPair(Polynomial::mkZero(), Constant::mkConstant(0));
+  }
+
+  static Node computeQR(const SumPair& sp, const Integer& div);
+
+
+  static SumPair comparisonToSumPair(const Comparison& cmp){
+    return SumPair(cmp.getLeft(), - cmp.getRight());
+  }
+
+};/* class SumPair */
+
 }/* CVC4::theory::arith namespace */
 }/* CVC4::theory namespace */
 }/* CVC4 namespace */
diff --git a/src/theory/arith/partial_model.cpp b/src/theory/arith/partial_model.cpp
index 73f80a70d..4a126ec55 100644
--- a/src/theory/arith/partial_model.cpp
+++ b/src/theory/arith/partial_model.cpp
@@ -28,6 +28,13 @@ using namespace CVC4::theory;
 using namespace CVC4::theory::arith;
 
 
+bool ArithPartialModel::boundsAreEqual(ArithVar x){
+  if(hasLowerBound(x) && hasUpperBound(x)){
+    return d_upperBound[x] == d_lowerBound[x];
+  }else{
+    return false;
+  }
+}
 
 void ArithPartialModel::zeroDifferenceDetected(ArithVar x){
   Assert(d_dm.isDifferenceSlack(x));
diff --git a/src/theory/arith/partial_model.h b/src/theory/arith/partial_model.h
index 6e85f7e80..cf0fc7d4e 100644
--- a/src/theory/arith/partial_model.h
+++ b/src/theory/arith/partial_model.h
@@ -116,7 +116,7 @@ private:
   void zeroDifferenceDetected(ArithVar x);
 
 public:
-
+  bool boundsAreEqual(ArithVar x);
 
   void setUpperBound(ArithVar x, const DeltaRational& r);
   void setLowerBound(ArithVar x, const DeltaRational& r);
diff --git a/src/theory/arith/theory_arith.cpp b/src/theory/arith/theory_arith.cpp
index ca0763a0c..824bb59ed 100644
--- a/src/theory/arith/theory_arith.cpp
+++ b/src/theory/arith/theory_arith.cpp
@@ -58,14 +58,17 @@ static const uint32_t RESET_START = 2;
 
 TheoryArith::TheoryArith(context::Context* c, context::UserContext* u, OutputChannel& out, Valuation valuation) :
   Theory(THEORY_ARITH, c, u, out, valuation),
+  d_hasDoneWorkSinceCut(false),
   d_atomsInContext(c),
   d_learner(d_pbSubstitutions),
   d_nextIntegerCheckVar(0),
+  d_constantIntegerVariables(c),
+  d_CivIterator(c,0),
+  d_varsInDioSolver(c),
   d_partialModel(c, d_differenceManager),
-  d_userVariables(),
   d_diseq(c),
   d_tableau(),
-  d_diosolver(c, d_tableau, d_partialModel),
+  d_diosolver(c),
   d_pbSubstitutions(u),
   d_restartsCounter(0),
   d_rowHasBeenAdded(false),
@@ -405,7 +408,7 @@ void TheoryArith::preRegisterTerm(TNode n) {
 }
 
 
-ArithVar TheoryArith::requestArithVar(TNode x, bool basic){
+ArithVar TheoryArith::requestArithVar(TNode x, bool slack){
   Assert(isLeaf(x) || x.getKind() == PLUS);
   Assert(!d_arithvarNodeMap.hasArithVar(x));
   Assert(x.getType().isReal());// real or integer
@@ -413,13 +416,22 @@ ArithVar TheoryArith::requestArithVar(TNode x, bool basic){
   ArithVar varX = d_variables.size();
   d_variables.push_back(Node(x));
   Debug("integers") << "isInteger[[" << x << "]]: " << x.getType().isInteger() << endl;
-  d_integerVars.push_back(x.getType().isPseudoboolean() ? 2 : (x.getType().isInteger() ? 1 : 0));
+
+  if(slack){
+    //The type computation is not quite accurate for Rationals that are integral.
+    //We'll use the isIntegral check from the polynomial package instead.
+    Polynomial p = Polynomial::parsePolynomial(x);
+    d_variableTypes.push_back(p.isIntegral() ? ATInteger : ATReal);
+  }else{
+    d_variableTypes.push_back(nodeToArithType(x));
+  }
+
+  d_slackVars.push_back(slack);
 
   d_simplex.increaseMax();
 
   d_arithvarNodeMap.setArithVar(x,varX);
 
-  d_userVariables.init(varX, !basic);
   d_tableau.increaseSize();
 
   Debug("arith::arithvar") << x << " |-> " << varX << endl;
@@ -577,11 +589,114 @@ bool TheoryArith::canSafelyAvoidEqualitySetup(TNode equality){
   return d_arithvarNodeMap.hasArithVar(equality[0]);
 }
 
+Comparison TheoryArith::mkIntegerEqualityFromAssignment(ArithVar v){
+  const DeltaRational& beta = d_partialModel.getAssignment(v);
+
+  Assert(beta.isIntegral());
+  Constant betaAsConstant = Constant::mkConstant(beta.floor());
+
+  TNode var = d_arithvarNodeMap.asNode(v);
+  Polynomial varAsPolynomial = Polynomial::parsePolynomial(var);
+  return Comparison::mkComparison(EQUAL, varAsPolynomial, betaAsConstant);
+}
+
+Node TheoryArith::dioCutting(){
+  context::Context::ScopedPush speculativePush(getContext());
+  //DO NOT TOUCH THE OUTPUTSTREAM
+
+  //TODO: Improve this
+  for(ArithVar v = 0, end = d_variables.size(); v != end; ++v){
+    if(isInteger(v)){
+      const DeltaRational& dr = d_partialModel.getAssignment(v);
+      if(d_partialModel.equalsUpperBound(v, dr) || d_partialModel.equalsLowerBound(v, dr)){
+        if(!d_partialModel.boundsAreEqual(v)){
+          // If the bounds are equal this is already in the dioSolver
+          //Add v = dr as a speculation.
+          Comparison eq = mkIntegerEqualityFromAssignment(v);
+          Assert(!eq.isBoolean());
+          d_diosolver.pushInputConstraint(eq, eq.getNode());
+          // It does not matter what the explanation of eq is.
+          // It cannot be used in a conflict
+        }
+      }
+    }
+  }
+
+  SumPair plane = d_diosolver.processEquationsForCut();
+  if(plane.isZero()){
+    return Node::null();
+  }else{
+    Polynomial p = plane.getPolynomial();
+    Constant c = plane.getConstant() * Constant::mkConstant(-1);
+    Integer gcd = p.gcd();
+    Assert(p.isIntegral());
+    Assert(c.isIntegral());
+    Assert(gcd > 1);
+    Assert(!gcd.divides(c.getValue().getNumerator()));
+    Comparison leq = Comparison::mkComparison(LEQ, p, c);
+    Comparison geq = Comparison::mkComparison(GEQ, p, c);
+    Node lemma = NodeManager::currentNM()->mkNode(OR, leq.getNode(), geq.getNode());
+    Node rewrittenLemma = Rewriter::rewrite(lemma);
+    Debug("arith::dio") << "dioCutting found the plane: " << plane.getNode() << endl;
+    Debug("arith::dio") << "resulting in the cut: " << lemma << endl;
+    Debug("arith::dio") << "rewritten " << rewrittenLemma << endl;
+    return rewrittenLemma;
+  }
+}
+
+Node TheoryArith::callDioSolver(){
+  while(d_CivIterator < d_constantIntegerVariables.size()){
+    ArithVar v = d_constantIntegerVariables[d_CivIterator];
+    d_CivIterator = d_CivIterator + 1;
+
+    Debug("arith::dio")  << v << endl;
+
+    Assert(isInteger(v));
+    Assert(d_partialModel.boundsAreEqual(v));
+
+    if(d_varsInDioSolver.find(v) != d_varsInDioSolver.end()){
+      continue;
+    }else{
+      d_varsInDioSolver.insert(v);
+    }
+
+    TNode lb = d_partialModel.getLowerConstraint(v);
+    TNode ub = d_partialModel.getUpperConstraint(v);
+
+    Node orig = Node::null();
+    if(lb == ub){
+      Assert(lb.getKind() == EQUAL);
+      orig = lb;
+    }else{
+      NodeBuilder<> nb(AND);
+      nb << ub << lb;
+      orig = nb;
+    }
+
+    Assert(d_partialModel.assignmentIsConsistent(v));
+
+    Comparison eq = mkIntegerEqualityFromAssignment(v);
+
+    if(eq.isBoolean()){
+      //This can only be a conflict
+      Assert(!eq.getNode().getConst<bool>());
+
+      //This should be handled by the normal form earlier in the case of equality
+      Assert(orig.getKind() != EQUAL);
+      return orig;
+    }else{
+      d_diosolver.pushInputConstraint(eq, orig);
+    }
+  }
+
+  return d_diosolver.processEquationsForConflict();
+}
+
 Node TheoryArith::assertionCases(TNode assertion){
-  Kind simpKind = simplifiedKind(assertion);
-  Assert(simpKind != UNDEFINED_KIND);
-  if(simpKind == EQUAL || simpKind == DISTINCT){
-    Node eq = (simpKind == DISTINCT) ? assertion[0] : assertion;
+  Kind simpleKind = simplifiedKind(assertion);
+  Assert(simpleKind != UNDEFINED_KIND);
+  if(simpleKind == EQUAL || simpleKind == DISTINCT){
+    Node eq = (simpleKind == DISTINCT) ? assertion[0] : assertion;
 
     if(!isSetup(eq)){
       //The previous code was equivalent to:
@@ -592,35 +707,71 @@ Node TheoryArith::assertionCases(TNode assertion){
     }
   }
 
-  ArithVar x_i = determineLeftVariable(assertion, simpKind);
-  DeltaRational c_i = determineRightConstant(assertion, simpKind);
+  ArithVar x_i = determineLeftVariable(assertion, simpleKind);
+  DeltaRational c_i = determineRightConstant(assertion, simpleKind);
+
+  bool tightened = false;
+
+  //If the variable is an integer tighen the constraint.
+  if(isInteger(x_i)){
+    if(simpleKind == LT){
+      tightened = true;
+      c_i = DeltaRational(c_i.floor());
+    }else if(simpleKind == GT){
+      tightened = true;
+      c_i = DeltaRational(c_i.ceiling());
+    }
+  }
 
   Debug("arith::assertions")  << "arith assertion @" << getContext()->getLevel()
                               <<"(" << assertion
 			     << " \\-> "
-			     << x_i<<" "<< simpKind <<" "<< c_i << ")" << std::endl;
-  switch(simpKind){
+			     << x_i<<" "<< simpleKind <<" "<< c_i << ")" << std::endl;
+
+  switch(simpleKind){
   case LEQ:
-    if (d_partialModel.hasLowerBound(x_i) && d_partialModel.getLowerBound(x_i) == c_i) {
-      Node diseq = assertion[0].eqNode(assertion[1]).notNode();
-      if (d_diseq.find(diseq) != d_diseq.end()) {
-        Node lb = d_partialModel.getLowerConstraint(x_i);
-        return disequalityConflict(diseq, lb , assertion);
+  case LT:
+    if(simpleKind == LEQ || (simpleKind == LT && tightened)){
+      if (d_partialModel.hasLowerBound(x_i) && d_partialModel.getLowerBound(x_i) == c_i) {
+        //If equal
+        TNode left  = getSide<false>(assertion, simpleKind);
+        TNode right = getSide<true>(assertion, simpleKind);
+
+        Node diseq = left.eqNode(right).notNode();
+        if (d_diseq.find(diseq) != d_diseq.end()) {
+          Node lb = d_partialModel.getLowerConstraint(x_i);
+          return disequalityConflict(diseq, lb , assertion);
+        }
+
+        if(isInteger(x_i)){
+          d_constantIntegerVariables.push_back(x_i);
+        }
       }
     }
-  case LT:
     return  d_simplex.AssertUpper(x_i, c_i, assertion);
   case GEQ:
-    if (d_partialModel.hasUpperBound(x_i) && d_partialModel.getUpperBound(x_i) == c_i) {
-      Node diseq = assertion[0].eqNode(assertion[1]).notNode();
-      if (d_diseq.find(diseq) != d_diseq.end()) {
-        Node ub = d_partialModel.getUpperConstraint(x_i);
-        return disequalityConflict(diseq, assertion, ub);
+  case GT:
+    if(simpleKind == GEQ || (simpleKind == GT && tightened)){
+      if (d_partialModel.hasUpperBound(x_i) && d_partialModel.getUpperBound(x_i) == c_i) {
+        //If equal
+        TNode left  = getSide<false>(assertion, simpleKind);
+        TNode right = getSide<true>(assertion, simpleKind);
+
+        Node diseq = left.eqNode(right).notNode();
+        if (d_diseq.find(diseq) != d_diseq.end()) {
+          Node ub = d_partialModel.getUpperConstraint(x_i);
+          return disequalityConflict(diseq, assertion, ub);
+        }
+        if(isInteger(x_i)){
+          d_constantIntegerVariables.push_back(x_i);
+        }
       }
     }
-  case GT:
     return d_simplex.AssertLower(x_i, c_i, assertion);
   case EQUAL:
+    if(isInteger(x_i)){
+      d_constantIntegerVariables.push_back(x_i);
+    }
     return d_simplex.AssertEquality(x_i, c_i, assertion);
   case DISTINCT:
     {
@@ -649,11 +800,34 @@ Node TheoryArith::assertionCases(TNode assertion){
   }
 }
 
+/**
+ * Looks for the next integer variable without an integer assignment in a round robin fashion.
+ * Changes the value of d_nextIntegerCheckVar.
+ *
+ * If this returns false, d_nextIntegerCheckVar does not have an integer assignment.
+ * If this returns true, all integer variables have an integer assignment.
+ */
+bool TheoryArith::hasIntegerModel(){
+  if(d_variables.size() > 0){
+    const ArithVar rrEnd = d_nextIntegerCheckVar;
+    do {
+      //Do not include slack variables
+      if(isInteger(d_nextIntegerCheckVar) && !isSlackVariable(d_nextIntegerCheckVar)) { // integer
+        const DeltaRational& d = d_partialModel.getAssignment(d_nextIntegerCheckVar);
+        if(!d.isIntegral()){
+          return false;
+        }
+      }
+    } while((d_nextIntegerCheckVar = (1 + d_nextIntegerCheckVar == d_variables.size() ? 0 : 1 + d_nextIntegerCheckVar)) != rrEnd);
+  }
+  return true;
+}
 
 
 void TheoryArith::check(Effort effortLevel){
   Debug("arith") << "TheoryArith::check begun" << std::endl;
 
+  d_hasDoneWorkSinceCut = d_hasDoneWorkSinceCut || !done();
   while(!done()){
 
     Node assertion = get();
@@ -672,6 +846,7 @@ void TheoryArith::check(Effort effortLevel){
     debugPrintAssertions();
   }
 
+  bool emmittedConflictOrSplit = false;
   Node possibleConflict = d_simplex.updateInconsistentVars();
   if(possibleConflict != Node::null()){
     d_partialModel.revertAssignmentChanges();
@@ -679,146 +854,184 @@ void TheoryArith::check(Effort effortLevel){
     Debug("arith::conflict") << "conflict   " << possibleConflict << endl;
 
     d_out->conflict(possibleConflict);
+    emmittedConflictOrSplit = true;
   }else{
     d_partialModel.commitAssignmentChanges();
-
-    if (fullEffort(effortLevel)) {
-      splitDisequalities();
-    }
   }
 
-  if(fullEffort(effortLevel) && d_integerVars.size() > 0) {
-    const ArithVar rrEnd = d_nextIntegerCheckVar;
-    do {
-      ArithVar v = d_nextIntegerCheckVar;
-      short type = d_integerVars[v];
-      if(type > 0) { // integer
-        const DeltaRational& d = d_partialModel.getAssignment(v);
-        const Rational& r = d.getNoninfinitesimalPart();
-        const Rational& i = d.getInfinitesimalPart();
-        Trace("integers") << "integers: assignment to [[" << d_arithvarNodeMap.asNode(v) << "]] is " << r << "[" << i << "]" << endl;
-        if(type == 2) {
-          // pseudoboolean
-          if(r.getDenominator() == 1 && i.getNumerator() == 0 &&
-             (r.getNumerator() == 0 || r.getNumerator() == 1)) {
-            // already pseudoboolean; skip
-            continue;
-          }
+  if(!emmittedConflictOrSplit && fullEffort(effortLevel)){
+    emmittedConflictOrSplit = splitDisequalities();
+  }
 
-          TNode var = d_arithvarNodeMap.asNode(v);
-          Node zero = NodeManager::currentNM()->mkConst(Integer(0));
-          Node one = NodeManager::currentNM()->mkConst(Integer(1));
-          Node eq0 = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::EQUAL, var, zero));
-          Node eq1 = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::EQUAL, var, one));
-          Node lem = NodeManager::currentNM()->mkNode(kind::OR, eq0, eq1);
-          Trace("pb") << "pseudobooleans: branch & bound: " << lem << endl;
-          Trace("integers") << "pseudobooleans: branch & bound: " << lem << endl;
-          //d_out->lemma(lem);
-        }
-        if(r.getDenominator() == 1 && i.getNumerator() == 0) {
-          // already an integer assignment; skip
-          continue;
-        }
+  if(!emmittedConflictOrSplit && fullEffort(effortLevel) && !hasIntegerModel()){
 
-        // otherwise, try the Diophantine equation solver
-        //bool result = d_diosolver.solve();
-        //Debug("integers") << "the dio solver returned " << (result ? "true" : "false") << endl;
-
-        // branch and bound
-        if(r.getDenominator() == 1) {
-          // r is an integer, but the infinitesimal might not be
-          if(i.getNumerator() < 0) {
-            // lemma: v <= r - 1 || v >= r
-
-            TNode var = d_arithvarNodeMap.asNode(v);
-            Node nrMinus1 = NodeManager::currentNM()->mkConst(r - 1);
-            Node nr = NodeManager::currentNM()->mkConst(r);
-            Node leq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::LEQ, var, nrMinus1));
-            Node geq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::GEQ, var, nr));
-
-            Node lem = NodeManager::currentNM()->mkNode(kind::OR, leq, geq);
-            Trace("integers") << "integers: branch & bound: " << lem << endl;
-            if(d_valuation.isSatLiteral(lem[0])) {
-              Debug("integers") << "    " << lem[0] << " == " << d_valuation.getSatValue(lem[0]) << endl;
-            } else {
-              Debug("integers") << "    " << lem[0] << " is not assigned a SAT literal" << endl;
-            }
-            if(d_valuation.isSatLiteral(lem[1])) {
-              Debug("integers") << "    " << lem[1] << " == " << d_valuation.getSatValue(lem[1]) << endl;
-            } else {
-              Debug("integers") << "    " << lem[1] << " is not assigned a SAT literal" << endl;
-            }
-            ++(d_statistics.d_externalBranchAndBounds);
-            d_out->lemma(lem);
-
-            // split only on one var
-            break;
-          } else if(i.getNumerator() > 0) {
-            // lemma: v <= r || v >= r + 1
-
-            TNode var = d_arithvarNodeMap.asNode(v);
-            Node nr = NodeManager::currentNM()->mkConst(r);
-            Node nrPlus1 = NodeManager::currentNM()->mkConst(r + 1);
-            Node leq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::LEQ, var, nr));
-            Node geq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::GEQ, var, nrPlus1));
-
-            Node lem = NodeManager::currentNM()->mkNode(kind::OR, leq, geq);
-            Trace("integers") << "integers: branch & bound: " << lem << endl;
-            if(d_valuation.isSatLiteral(lem[0])) {
-              Debug("integers") << "    " << lem[0] << " == " << d_valuation.getSatValue(lem[0]) << endl;
-            } else {
-              Debug("integers") << "    " << lem[0] << " is not assigned a SAT literal" << endl;
-            }
-            if(d_valuation.isSatLiteral(lem[1])) {
-              Debug("integers") << "    " << lem[1] << " == " << d_valuation.getSatValue(lem[1]) << endl;
-            } else {
-              Debug("integers") << "    " << lem[1] << " is not assigned a SAT literal" << endl;
-            }
-            ++(d_statistics.d_externalBranchAndBounds);
-            d_out->lemma(lem);
+    if(!emmittedConflictOrSplit){
+      possibleConflict = callDioSolver();
+      if(possibleConflict != Node::null()){
+        Debug("arith::conflict") << "dio conflict   " << possibleConflict << endl;
+        d_out->conflict(possibleConflict);
+        emmittedConflictOrSplit = true;
+      }
+    }
 
-            // split only on one var
-            break;
-          } else {
-            Unreachable();
-          }
-        } else {
-          // lemma: v <= floor(r) || v >= ceil(r)
-
-          TNode var = d_arithvarNodeMap.asNode(v);
-          Node floor = NodeManager::currentNM()->mkConst(r.floor());
-          Node ceiling = NodeManager::currentNM()->mkConst(r.ceiling());
-          Node leq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::LEQ, var, floor));
-          Node geq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::GEQ, var, ceiling));
-
-          Node lem = NodeManager::currentNM()->mkNode(kind::OR, leq, geq);
-          Trace("integers") << "integers: branch & bound: " << lem << endl;
-          if(d_valuation.isSatLiteral(lem[0])) {
-            Debug("integers") << "    " << lem[0] << " == " << d_valuation.getSatValue(lem[0]) << endl;
-          } else {
-            Debug("integers") << "    " << lem[0] << " is not assigned a SAT literal" << endl;
-          }
-          if(d_valuation.isSatLiteral(lem[1])) {
-            Debug("integers") << "    " << lem[1] << " == " << d_valuation.getSatValue(lem[1]) << endl;
-          } else {
-            Debug("integers") << "    " << lem[1] << " is not assigned a SAT literal" << endl;
-          }
-          ++(d_statistics.d_externalBranchAndBounds);
-          d_out->lemma(lem);
+    if(!emmittedConflictOrSplit && d_hasDoneWorkSinceCut){
+      Node possibleLemma = dioCutting();
+      if(!possibleLemma.isNull()){
+        Debug("arith") << "dio cut   " << possibleLemma << endl;
+        emmittedConflictOrSplit = true;
+        d_hasDoneWorkSinceCut = false;
+        d_out->lemma(possibleLemma);
+      }
+    }
 
-          // split only on one var
-          break;
-        }
-      }// if(arithvar is integer-typed)
-    } while((d_nextIntegerCheckVar = (1 + d_nextIntegerCheckVar == d_variables.size() ? 0 : 1 + d_nextIntegerCheckVar)) != rrEnd);
-  }// if(full effort)
+    if(!emmittedConflictOrSplit) {
+      Node possibleLemma = roundRobinBranch();
+      if(!possibleLemma.isNull()){
+        ++(d_statistics.d_externalBranchAndBounds);
+        emmittedConflictOrSplit = true;
+        d_out->lemma(possibleLemma);
+      }
+    }
+  }//if !emmittedConflictOrSplit && fullEffort(effortLevel) && !hasIntegerModel()
 
   if(Debug.isOn("paranoid:check_tableau")){ d_simplex.debugCheckTableau(); }
   if(Debug.isOn("arith::print_model")) { debugPrintModel(); }
   Debug("arith") << "TheoryArith::check end" << std::endl;
 }
 
-void TheoryArith::splitDisequalities(){
+/** Returns true if the roundRobinBranching() issues a lemma. */
+Node TheoryArith::roundRobinBranch(){
+  if(hasIntegerModel()){
+    return Node::null();
+  }else{
+    ArithVar v = d_nextIntegerCheckVar;
+
+    Assert(isInteger(v));
+    Assert(!isSlackVariable(v));
+
+    const DeltaRational& d = d_partialModel.getAssignment(v);
+    const Rational& r = d.getNoninfinitesimalPart();
+    const Rational& i = d.getInfinitesimalPart();
+    Trace("integers") << "integers: assignment to [[" << d_arithvarNodeMap.asNode(v) << "]] is " << r << "[" << i << "]" << endl;
+
+    Assert(! (r.getDenominator() == 1 && i.getNumerator() == 0));
+    Assert(!d.isIntegral());
+
+    TNode var = d_arithvarNodeMap.asNode(v);
+    Integer floor_d = d.floor();
+    Integer ceil_d = d.ceiling();
+
+    Node leq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::LEQ, var, mkIntegerNode(floor_d)));
+    Node geq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::GEQ, var, mkIntegerNode(ceil_d)));
+
+
+    Node lem = NodeManager::currentNM()->mkNode(kind::OR, leq, geq);
+    Trace("integers") << "integers: branch & bound: " << lem << endl;
+    if(d_valuation.isSatLiteral(lem[0])) {
+      Debug("integers") << "    " << lem[0] << " == " << d_valuation.getSatValue(lem[0]) << endl;
+    } else {
+      Debug("integers") << "    " << lem[0] << " is not assigned a SAT literal" << endl;
+    }
+    if(d_valuation.isSatLiteral(lem[1])) {
+      Debug("integers") << "    " << lem[1] << " == " << d_valuation.getSatValue(lem[1]) << endl;
+    } else {
+      Debug("integers") << "    " << lem[1] << " is not assigned a SAT literal" << endl;
+    }
+    return lem;
+
+  //   // branch and bound
+  //   if(r.getDenominator() == 1) {
+  //     // r is an integer, but the infinitesimal might not be
+  //     if(i.getNumerator() < 0) {
+  //       // lemma: v <= r - 1 || v >= r
+
+  //       TNode var = d_arithvarNodeMap.asNode(v);
+  //       Node nrMinus1 = NodeManager::currentNM()->mkConst(r - 1);
+  //       Node nr = NodeManager::currentNM()->mkConst(r);
+  //       Node leq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::LEQ, var, nrMinus1));
+  //       Node geq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::GEQ, var, nr));
+
+  //       Node lem = NodeManager::currentNM()->mkNode(kind::OR, leq, geq);
+  //       Trace("integers") << "integers: branch & bound: " << lem << endl;
+  //       if(d_valuation.isSatLiteral(lem[0])) {
+  //         Debug("integers") << "    " << lem[0] << " == " << d_valuation.getSatValue(lem[0]) << endl;
+  //       } else {
+  //         Debug("integers") << "    " << lem[0] << " is not assigned a SAT literal" << endl;
+  //       }
+  //       if(d_valuation.isSatLiteral(lem[1])) {
+  //         Debug("integers") << "    " << lem[1] << " == " << d_valuation.getSatValue(lem[1]) << endl;
+  //       } else {
+  //         Debug("integers") << "    " << lem[1] << " is not assigned a SAT literal" << endl;
+  //       }
+  //       return lem;
+  //     } else if(i.getNumerator() > 0) {
+  //         // lemma: v <= r || v >= r + 1
+
+  //         TNode var = d_arithvarNodeMap.asNode(v);
+  //         Node nr = NodeManager::currentNM()->mkConst(r);
+  //         Node nrPlus1 = NodeManager::currentNM()->mkConst(r + 1);
+  //         Node leq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::LEQ, var, nr));
+  //         Node geq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::GEQ, var, nrPlus1));
+
+  //         Node lem = NodeManager::currentNM()->mkNode(kind::OR, leq, geq);
+  //         Trace("integers") << "integers: branch & bound: " << lem << endl;
+  //         if(d_valuation.isSatLiteral(lem[0])) {
+  //           Debug("integers") << "    " << lem[0] << " == " << d_valuation.getSatValue(lem[0]) << endl;
+  //         } else {
+  //           Debug("integers") << "    " << lem[0] << " is not assigned a SAT literal" << endl;
+  //         }
+  //         if(d_valuation.isSatLiteral(lem[1])) {
+  //           Debug("integers") << "    " << lem[1] << " == " << d_valuation.getSatValue(lem[1]) << endl;
+  //         } else {
+  //           Debug("integers") << "    " << lem[1] << " is not assigned a SAT literal" << endl;
+  //         }
+  //         ++(d_statistics.d_externalBranchAndBounds);
+  //         d_out->lemma(lem);
+  //         result = true;
+
+  //         // split only on one var
+  //         break;
+  //       } else {
+  //         Unreachable();
+  //       }
+  //     } else {
+  //       // lemma: v <= floor(r) || v >= ceil(r)
+
+  //       TNode var = d_arithvarNodeMap.asNode(v);
+  //       Node floor = NodeManager::currentNM()->mkConst(r.floor());
+  //       Node ceiling = NodeManager::currentNM()->mkConst(r.ceiling());
+  //       Node leq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::LEQ, var, floor));
+  //       Node geq = Rewriter::rewrite(NodeManager::currentNM()->mkNode(kind::GEQ, var, ceiling));
+
+  //       Node lem = NodeManager::currentNM()->mkNode(kind::OR, leq, geq);
+  //       Trace("integers") << "integers: branch & bound: " << lem << endl;
+  //       if(d_valuation.isSatLiteral(lem[0])) {
+  //         Debug("integers") << "    " << lem[0] << " == " << d_valuation.getSatValue(lem[0]) << endl;
+  //       } else {
+  //         Debug("integers") << "    " << lem[0] << " is not assigned a SAT literal" << endl;
+  //       }
+  //       if(d_valuation.isSatLiteral(lem[1])) {
+  //         Debug("integers") << "    " << lem[1] << " == " << d_valuation.getSatValue(lem[1]) << endl;
+  //       } else {
+  //         Debug("integers") << "    " << lem[1] << " is not assigned a SAT literal" << endl;
+  //       }
+  //       ++(d_statistics.d_externalBranchAndBounds);
+  //       d_out->lemma(lem);
+  //       result = true;
+
+  //       // split only on one var
+  //       break;
+  //     }
+  //   }// if(arithvar is integer-typed)
+  // } while((d_nextIntegerCheckVar = (1 + d_nextIntegerCheckVar == d_variables.size() ? 0 : 1 + d_nextIntegerCheckVar)) != rrEnd);
+
+  // return result;
+  }
+}
+
+bool TheoryArith::splitDisequalities(){
+  bool splitSomething = false;
+
   context::CDSet<Node, NodeHashFunction>::iterator it = d_diseq.begin();
   context::CDSet<Node, NodeHashFunction>::iterator it_end = d_diseq.end();
   for(; it != it_end; ++ it) {
@@ -839,8 +1052,10 @@ void TheoryArith::splitDisequalities(){
       Node lemma = NodeBuilder<3>(OR) << eq << ltNode << gtNode;
       ++(d_statistics.d_statDisequalitySplits);
       d_out->lemma(lemma);
+      splitSomething = true;
     }
   }
+  return splitSomething;
 }
 
 /**
@@ -852,12 +1067,12 @@ void TheoryArith::debugPrintAssertions() {
   for (ArithVar i = 0; i < d_variables.size(); ++ i) {
     if (d_partialModel.hasLowerBound(i)) {
       Node lConstr = d_partialModel.getLowerConstraint(i);
-      Debug("arith::print_assertions") << lConstr.toString() << endl;
+      Debug("arith::print_assertions") << lConstr << endl;
     }
 
     if (d_partialModel.hasUpperBound(i)) {
       Node uConstr = d_partialModel.getUpperConstraint(i);
-      Debug("arith::print_assertions") << uConstr.toString() << endl;
+      Debug("arith::print_assertions") << uConstr << endl;
     }
   }
   context::CDSet<Node, NodeHashFunction>::iterator it = d_diseq.begin();
@@ -1162,7 +1377,7 @@ void TheoryArith::presolve(){
     for(VarIter i = d_variables.begin(), end = d_variables.end(); i != end; ++i){
       Node variableNode = *i;
       ArithVar var = d_arithvarNodeMap.asArithVar(variableNode);
-      if(d_userVariables.isMember(var) &&
+      if(!isSlackVariable(var) &&
          !d_atomDatabase.hasAnyAtoms(variableNode) &&
          !variableNode.getType().isInteger()){
 	//The user variable is unconstrained.
diff --git a/src/theory/arith/theory_arith.h b/src/theory/arith/theory_arith.h
index 6255efbbc..6b14ae6ff 100644
--- a/src/theory/arith/theory_arith.h
+++ b/src/theory/arith/theory_arith.h
@@ -12,8 +12,7 @@
  ** information.\endverbatim
  **
  ** \brief Arithmetic theory.
- **
- ** Arithmetic theory.
+ ** ** Arithmetic theory.
  **/
 
 #include "cvc4_private.h"
@@ -59,6 +58,8 @@ namespace arith {
 class TheoryArith : public Theory {
 private:
 
+  bool d_hasDoneWorkSinceCut;
+
   /**
    * The set of atoms that are currently in the context.
    * This is exactly the union of preregistered atoms and
@@ -108,10 +109,40 @@ private:
 
 
   /**
-   * List of the types of variables in the system.
-   * "True" means integer, "false" means (non-integer) real.
+   * (For the moment) the type hierarchy goes as:
+   * PsuedoBoolean <: Integer <: Real
+   * The type number of a variable is an integer representing the most specific
+   * type of the variable. The possible values of type number are:
    */
-  std::vector<short> d_integerVars;
+  enum ArithType
+    {
+      ATReal = 0,
+      ATInteger = 1,
+      ATPsuedoBoolean = 2
+   };
+
+  std::vector<ArithType> d_variableTypes;
+  inline ArithType nodeToArithType(TNode x) const {
+    return x.getType().isPseudoboolean() ? ATPsuedoBoolean :
+      (x.getType().isInteger() ? ATInteger : ATReal);
+  }
+
+  /** Returns true if x is of type Integer or PsuedoBoolean. */
+  inline bool isInteger(ArithVar x) const {
+    return d_variableTypes[x] >= ATInteger;
+  }
+  /** Returns true if x is of type PsuedoBoolean. */
+  inline bool isPsuedoBoolean(ArithVar x) const {
+    return d_variableTypes[x] == ATPsuedoBoolean;
+  }
+
+  /** This is the set of variables initially introduced as slack variables. */
+  std::vector<bool> d_slackVars;
+
+  /** Returns true if the variable was initially introduced as a slack variable. */
+  inline bool isSlackVariable(ArithVar x) const{
+    return d_slackVars[x];
+  }
 
   /**
    * On full effort checks (after determining LA(Q) satisfiability), we
@@ -123,6 +154,21 @@ private:
   ArithVar d_nextIntegerCheckVar;
 
   /**
+   * Queue of Integer variables that are known to be equal to a constant.
+   */
+  context::CDList<ArithVar> d_constantIntegerVariables;
+  /** Iterator over d_constantIntegerVariables. */
+  context::CDO<unsigned int> d_CivIterator;
+
+  Node callDioSolver();
+  Node dioCutting();
+
+  Comparison mkIntegerEqualityFromAssignment(ArithVar v);
+
+  #warning "DO NOT COMMIT TO TRUNK, USE MORE EFFICIENT CHECK INSTEAD"
+  CDArithVarSet d_varsInDioSolver;
+
+  /**
    * If ArithVar v maps to the node n in d_removednode,
    * then n = (= asNode(v) rhs) where rhs is a term that
    * can be used to determine the value of n using getValue().
@@ -136,11 +182,6 @@ private:
   ArithPartialModel d_partialModel;
 
   /**
-   * Set of ArithVars that were introduced via preregisteration.
-   */
-  ArithVarSet d_userVariables;
-
-  /**
    * List of all of the inequalities asserted in the current context.
    */
   context::CDSet<Node, NodeHashFunction> d_diseq;
@@ -246,15 +287,37 @@ private:
    */
   ArithVar determineLeftVariable(TNode assertion, Kind simpleKind);
 
-  /** Splits the disequalities in d_diseq that are violated using lemmas on demand. */
-  void splitDisequalities();
+  /**
+   * Splits the disequalities in d_diseq that are violated using lemmas on demand.
+   * returns true if any lemmas were issued.
+   * returns false if all disequalities are satisified in the current model.
+   */
+  bool splitDisequalities();
+
+
+
+  /**
+   * Looks for the next integer variable without an integer assignment in a round robin fashion.
+   * Changes the value of d_nextIntegerCheckVar.
+   *
+   * If this returns false, d_nextIntegerCheckVar does not have an integer assignment.
+   * If this returns true, all integer variables have an integer assignment.
+   */
+  bool hasIntegerModel();
+
+  /**
+   * Issues branches for non-slack integer variables with non-integer assignments.
+   * Returns a cut for a lemma.
+   * If there is an integer model, this returns Node::null().
+   */
+  Node roundRobinBranch();
 
   /**
    * This requests a new unique ArithVar value for x.
    * This also does initial (not context dependent) set up for a variable,
    * except for setting up the initial.
    */
-  ArithVar requestArithVar(TNode x, bool basic);
+  ArithVar requestArithVar(TNode x, bool slack);
 
   /** Initial (not context dependent) sets up for a variable.*/
   void setupInitialValue(ArithVar x);
diff --git a/src/util/bitvector.h b/src/util/bitvector.h
index f05ebaf17..d7f0e13a5 100644
--- a/src/util/bitvector.h
+++ b/src/util/bitvector.h
@@ -91,15 +91,18 @@ public:
   }
 
   BitVector operator ~() const {
+    //is this right? it looks like a no-op?
     return BitVector(d_size, d_value);
   }
 
   BitVector concat (const BitVector& other) const {
-    return BitVector(d_size + other.d_size, (d_value * Integer(2).pow(other.d_size)) + other.d_value);
+    return BitVector(d_size + other.d_size, (d_value.multiplyByPow2(other.d_size)) + other.d_value);
+    //return BitVector(d_size + other.d_size, (d_value * Integer(2).pow(other.d_size)) + other.d_value);
   }
 
   BitVector extract(unsigned high, unsigned low) const {
-    return BitVector(high - low + 1, (d_value % (Integer(2).pow(high + 1))) / Integer(2).pow(low));
+    return BitVector(high - low + 1, d_value.extractBitRange(high - low + 1, low));
+    //return BitVector(high - low + 1, (d_value % (Integer(2).pow(high + 1))) / Integer(2).pow(low));
   }
 
   size_t hash() const {
diff --git a/src/util/integer_cln_imp.h b/src/util/integer_cln_imp.h
index f8ffc0d65..06459e3e1 100644
--- a/src/util/integer_cln_imp.h
+++ b/src/util/integer_cln_imp.h
@@ -167,6 +167,7 @@ public:
     return *this;
   }
 
+  /*
   Integer operator/(const Integer& y) const {
     return Integer( cln::floor1(d_value, y.d_value) );
   }
@@ -182,6 +183,65 @@ public:
     d_value = cln::floor2(d_value, y.d_value).remainder;
     return *this;
   }
+  */
+
+  /**
+   * Return this*(2^pow).
+   */
+  Integer multiplyByPow2(uint32_t pow) const {
+    cln::cl_I ipow(pow);
+    return Integer( d_value << ipow);
+  }
+
+  /** See CLN Documentation. */
+  Integer extractBitRange(uint32_t bitCount, uint32_t low) const {
+    cln::cl_byte range(bitCount, low);
+    return Integer(cln::ldb(d_value, range));
+  }
+
+  /**
+   * Returns the floor(this / y)
+   */
+  Integer floorDivideQuotient(const Integer& y) const {
+    return Integer( cln::floor1(d_value, y.d_value) );
+  }
+
+  /**
+   * Returns r == this - floor(this/y)*y
+   */
+  Integer floorDivideRemainder(const Integer& y) const {
+    return Integer( cln::floor2(d_value, y.d_value).remainder );
+  }
+   /**
+   * Computes a floor quoient and remainder for x divided by y.
+   */
+  static void floorQR(Integer& q, Integer& r, const Integer& x, const Integer& y) {
+    cln::cl_I_div_t res = cln::floor2(x.d_value, y.d_value);
+    q.d_value = res.quotient;
+    r.d_value = res.remainder;
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideQuotient(const Integer& y) const {
+    return Integer( cln::ceiling1(d_value, y.d_value) );
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideRemainder(const Integer& y) const {
+    return Integer( cln::ceiling2(d_value, y.d_value).remainder );
+  }
+
+  /**
+   * If y divides *this, then exactQuotient returns (this/y)
+   */
+  Integer exactQuotient(const Integer& y) const {
+    Assert(y.divides(*this));
+    return Integer( cln::exquo(d_value, y.d_value) );
+  }
 
   /**
    * Raise this Integer to the power <code>exp</code>.
@@ -208,6 +268,22 @@ public:
   }
 
   /**
+   * Return the least common multiple of this integer with another.
+   */
+  Integer lcm(const Integer& y) const {
+    cln::cl_I result = cln::lcm(d_value, y.d_value);
+    return Integer(result);
+  }
+
+  /**
+   * Return true if *this exactly divides y.
+   */
+  bool divides(const Integer& y) const {
+    cln::cl_I result = cln::rem(y.d_value, d_value);
+    return cln::zerop(result);
+  }
+
+  /**
    * Return the absolute value of this integer.
    */
   Integer abs() const {
@@ -243,6 +319,12 @@ public:
     return output;
   }
 
+  int sgn() const {
+    cln::cl_I sgn = cln::signum(d_value);
+    Assert(sgn == 0 || sgn == -1 || sgn == 1);
+    return cln::cl_I_to_int(sgn);
+  }
+
   //friend std::ostream& operator<<(std::ostream& os, const Integer& n);
 
   long getLong() const {
@@ -281,6 +363,27 @@ public:
     return cln::logbitp(n, d_value);
   }
 
+  /**
+   * If x != 0, returns the unique n s.t. 2^{n-1} <= abs(x) < 2^{n}.
+   * If x == 0, returns 1.
+   */
+  size_t length() const {
+    int s = sgn();
+    if(s == 0){
+      return 1;
+    }else if(s < 0){
+      return cln::integer_length(-d_value);
+    }else{
+      return cln::integer_length(d_value);
+    }
+  }
+
+/*   cl_I xgcd (const cl_I& a, const cl_I& b, cl_I* u, cl_I* v) */
+/* This function ("extended gcd") returns the greatest common divisor g of a and b and at the same time the representation of g as an integral linear combination of a and b: u and v with u*a+v*b = g, g >= 0. u and v will be normalized to be of smallest possible absolute value, in the following sense: If a and b are non-zero, and abs(a) != abs(b), u and v will satisfy the inequalities abs(u) <= abs(b)/(2*g), abs(v) <= abs(a)/(2*g). */
+  static void extendedGcd(Integer& g, Integer& s, Integer& t, const Integer& a, const Integer& b){    
+    g.d_value = cln::xgcd(a.d_value, b.d_value, &s.d_value, &t.d_value);
+  }
+
   friend class CVC4::Rational;
 };/* class Integer */
 
diff --git a/src/util/integer_gmp_imp.h b/src/util/integer_gmp_imp.h
index 16ca8313b..161666df5 100644
--- a/src/util/integer_gmp_imp.h
+++ b/src/util/integer_gmp_imp.h
@@ -135,20 +135,82 @@ public:
     return *this;
   }
 
-  Integer operator/(const Integer& y) const {
-    return Integer( d_value / y.d_value );
+  /**
+   * Return this*(2^pow).
+   */
+  Integer multiplyByPow2(uint32_t pow) const{
+    mpz_class result;
+    mpz_mul_2exp(result.get_mpz_t(), d_value.get_mpz_t(), pow);
+    return Integer( result );
   }
-  Integer& operator/=(const Integer& y) {
-    d_value /= y.d_value;
-    return *this;
+
+  /** See GMP Documentation. */
+  Integer extractBitRange(uint32_t bitCount, uint32_t low) const {
+    // bitCount = high-low+1
+    uint32_t high = low + bitCount-1;
+    //— Function: void mpz_fdiv_r_2exp (mpz_t r, mpz_t n, mp_bitcnt_t b)
+    mpz_class rem, div;
+    mpz_fdiv_r_2exp(rem.get_mpz_t(), d_value.get_mpz_t(), high+1);
+    mpz_fdiv_q_2exp(div.get_mpz_t(), rem.get_mpz_t(), low);
+
+    return Integer(div);
   }
 
-  Integer operator%(const Integer& y) const {
-    return Integer( d_value % y.d_value );
+  /**
+   * Returns the floor(this / y)
+   */
+  Integer floorDivideQuotient(const Integer& y) const {
+    mpz_class q;
+    mpz_fdiv_q(q.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( q );
   }
-  Integer& operator%=(const Integer& y) {
-    d_value %= y.d_value;
-    return *this;
+
+  /**
+   * Returns r == this - floor(this/y)*y
+   */
+  Integer floorDivideRemainder(const Integer& y) const {
+    mpz_class r;
+    mpz_fdiv_r(r.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( r );
+  }
+
+  /**
+   * Computes a floor quoient and remainder for x divided by y.
+   */
+  static void floorQR(Integer& q, Integer& r, const Integer& x, const Integer& y) {
+    mpz_fdiv_qr(q.d_value.get_mpz_t(), r.d_value.get_mpz_t(), x.d_value.get_mpz_t(), y.d_value.get_mpz_t());
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideQuotient(const Integer& y) const {
+    mpz_class q;
+    mpz_cdiv_q(q.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( q );
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideRemainder(const Integer& y) const {
+    mpz_class r;
+    mpz_cdiv_r(r.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( r );
+  }
+
+  /**
+   * If y divides *this, then exactQuotient returns (this/y)
+   */
+  Integer exactQuotient(const Integer& y) const {
+    Assert(y.divides(*this));
+    mpz_class q;
+    mpz_divexact(q.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( q );
+  }
+
+  int sgn() const {
+    return mpz_sgn(d_value.get_mpz_t());
   }
 
   /**
@@ -172,6 +234,24 @@ public:
   }
 
   /**
+   * Return the least common multiple of this integer with another.
+   */
+  Integer lcm(const Integer& y) const {
+    mpz_class result;
+    mpz_lcm(result.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer(result);
+  }
+
+  /**
+   * All non-zero integers z, z.divide(0)
+   * ! zero.divides(zero)
+   */
+  bool divides(const Integer& y) const {
+    int res = mpz_divisible_p(y.d_value.get_mpz_t(), d_value.get_mpz_t());
+    return res != 0;
+  }
+
+  /**
    * Return the absolute value of this integer.
    */
   Integer abs() const {
@@ -217,6 +297,24 @@ public:
     return mpz_tstbit(d_value.get_mpz_t(), n);
   }
 
+  /**
+   * If x != 0, returns the smallest n s.t. 2^{n-1} <= abs(x) < 2^{n}.
+   * If x == 0, returns 1.
+   */
+  size_t length() const {
+    if(sgn() == 0){
+      return 1;
+    }else{
+      return mpz_sizeinbase(d_value.get_mpz_t(),2);
+    }
+  }
+
+  static void extendedGcd(Integer& g, Integer& s, Integer& t, const Integer& a, const Integer& b){
+    //mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, mpz_t a, mpz_t b);
+    mpz_gcdext (g.d_value.get_mpz_t(), s.d_value.get_mpz_t(), t.d_value.get_mpz_t(), a.d_value.get_mpz_t(), b.d_value.get_mpz_t());
+  }
+
+
   friend class CVC4::Rational;
 };/* class Integer */
 
diff --git a/src/util/rational_cln_imp.h b/src/util/rational_cln_imp.h
index 2f2c14ed8..885e6b628 100644
--- a/src/util/rational_cln_imp.h
+++ b/src/util/rational_cln_imp.h
@@ -192,6 +192,18 @@ public:
     }
   }
 
+  Rational abs() const {
+    if(sgn() < 0){
+      return -(*this);
+    }else{
+      return *this;
+    }
+  }
+
+  bool isIntegral() const{
+    return getDenominator() == 1;
+  }
+
   Integer floor() const {
     return Integer(cln::floor1(d_value));
   }
diff --git a/src/util/rational_gmp_imp.h b/src/util/rational_gmp_imp.h
index 37c3c8364..4635ce881 100644
--- a/src/util/rational_gmp_imp.h
+++ b/src/util/rational_gmp_imp.h
@@ -169,6 +169,14 @@ public:
     return mpq_sgn(d_value.get_mpq_t());
   }
 
+  Rational abs() const {
+    if(sgn() < 0){
+      return -(*this);
+    }else{
+      return *this;
+    }
+  }
+
   Integer floor() const {
     mpz_class q;
     mpz_fdiv_q(q.get_mpz_t(), d_value.get_num_mpz_t(), d_value.get_den_mpz_t());
@@ -244,6 +252,10 @@ public:
     return (*this);
   }
 
+  bool isIntegral() const{
+    return getDenominator() == 1;
+  }
+
   /** Returns a string representing the rational in the given base. */
   std::string toString(int base = 10) const {
     return d_value.get_str(base);