theory_arith_new.cpp

Go to the documentation of this file.
00001 /*****************************************************************************/
00002 /*!
00003  * \file theory_arith_new.cpp
00004  * 
00005  * Author: Dejan Jovanovic
00006  * 
00007  * Created: Thu Jun 14 13:42:00 2007
00008  *
00009  * <hr>
00010  *
00011  * License to use, copy, modify, sell and/or distribute this software
00012  * and its documentation for any purpose is hereby granted without
00013  * royalty, subject to the terms and conditions defined in the \ref
00014  * LICENSE file provided with this distribution.
00015  * 
00016  * <hr>
00017  * 
00018  */
00019 /*****************************************************************************/
00020 
00021 
00022 #include "theory_arith_new.h" 
00023 #include "arith_proof_rules.h"
00024 //#include "arith_expr.h"
00025 #include "arith_exception.h"
00026 #include "typecheck_exception.h"
00027 #include "eval_exception.h"
00028 #include "parser_exception.h"
00029 #include "smtlib_exception.h"
00030 #include "theory_core.h"
00031 #include "command_line_flags.h"
00032 
00033 
00034 using namespace std;
00035 using namespace CVC3;
00036 
00037 
00038 ///////////////////////////////////////////////////////////////////////////////
00039 // TheoryArithNew::FreeConst Methods                                            //
00040 ///////////////////////////////////////////////////////////////////////////////
00041 
00042 namespace CVC3 {
00043 
00044 ostream& operator<<(ostream& os, const TheoryArithNew::FreeConst& fc) {
00045   os << "FreeConst(r=" << fc.getConst() << ", " 
00046      << (fc.strict()? "strict" : "non-strict") << ")";
00047   return os;
00048 }
00049 
00050 ///////////////////////////////////////////////////////////////////////////////
00051 // TheoryArithNew::Ineq Methods                                                 //
00052 ///////////////////////////////////////////////////////////////////////////////
00053 
00054 ostream& operator<<(ostream& os, const TheoryArithNew::Ineq& ineq) {
00055   os << "Ineq(" << ineq.ineq().getExpr() << ", isolated on "
00056      << (ineq.varOnRHS()? "RHS" : "LHS") << ", const = "
00057      << ineq.getConst() << ")";
00058   return os;
00059 }
00060 } // End of namespace CVC3
00061 
00062 
00063 ///////////////////////////////////////////////////////////////////////////////
00064 // TheoryArithNew Private Methods                                               //
00065 ///////////////////////////////////////////////////////////////////////////////
00066 
00067 
00068 Theorem TheoryArithNew::isIntegerThm(const Expr& e) {
00069   // Quick check
00070   if(isReal(e.getType())) return Theorem();
00071   // Try harder
00072   return isIntegerDerive(Expr(IS_INTEGER, e), typePred(e));
00073 }
00074 
00075 
00076 Theorem TheoryArithNew::isIntegerDerive(const Expr& isIntE, const Theorem& thm) {
00077   const Expr& e = thm.getExpr();
00078   // We found it!
00079   if(e == isIntE) return thm;
00080 
00081   Theorem res;
00082   // If the theorem is an AND, look inside each child
00083   if(e.isAnd()) {
00084     int i, iend=e.arity();
00085     for(i=0; i<iend; ++i) {
00086       res = isIntegerDerive(isIntE, getCommonRules()->andElim(thm, i));
00087       if(!res.isNull()) return res;
00088     }
00089   }
00090   return res;
00091 }
00092 
00093 
00094 bool TheoryArithNew::kidsCanonical(const Expr& e) {
00095   if(isLeaf(e)) return true;
00096   bool res(true);
00097   for(int i=0; res && i<e.arity(); ++i) {
00098     Expr simp(canon(e[i]).getRHS());
00099     res = (e[i] == simp);
00100     IF_DEBUG(if(!res) debugger.getOS() << "\ne[" << i << "] = " << e[i]
00101        << "\nsimplified = " << simp << endl);
00102   }
00103   return res;
00104 }
00105 
00106 
00107 
00108 
00109 /**
00110  * Compute a canonical form for expression e and return a theorem that e is equal to its canonical form.
00111  * Note that canonical form for arith expressions is of the following form:
00112  * 
00113  * b + a_1 x_1 + a_2 x_2 + ... + a_n x_n (ONE SUM EXPRESION)
00114  * 
00115  * or just a rational b (RATIONAL EXPRESSION)
00116  * 
00117  * where a_i and b are rationals and x_i are arithmetical leaves (i.e. variables or non arithmetic terms)
00118  * 
00119  * @author Clark Barrett
00120  * @author Vijay Ganesh
00121  * @since Sat Feb  8 14:46:32 2003
00122  */  
00123 Theorem TheoryArithNew::canon(const Expr& e)
00124 {
00125   // Trace the arith canon of the expression in the debug version
00126   TRACE("arith", "canon(", e , ") {");
00127   
00128   // The invariant is that the kids of the expression should already be canonised 
00129   //DebugAssert(kidsCanonical(e), "TheoryArithNew::canon("+e.toString()+")");
00130   
00131   // The resulting theorem object 
00132     Theorem result;
00133     
00134     switch (e.getKind()) {
00135       
00136       // -E
00137       case UMINUS: {
00138         
00139         // Convert the unary minus to multiplication with -1
00140         result = d_rules->uMinusToMult(e[0]);
00141         
00142         // Get the resulting epression
00143           Expr e2 = result.getRHS();
00144           
00145           // Canonise the resulting expression and combine the theorems
00146           result = transitivityRule(result, canon(e2));
00147       
00148         // UMINUS case break
00149         break;
00150       }
00151       
00152         // e[0] + e[1]
00153       case PLUS: 
00154         
00155         // Call the canonPlus procedure to canonise +
00156         result = d_rules->canonPlus(e);
00157           
00158           // PLUS case break
00159           break;
00160       
00161       // e[0] - e[1]
00162       case MINUS: {
00163           
00164           // Arity of the minus should be exactly 2 (TODO: looks useless, maybe remove)
00165           DebugAssert(e.arity() == 2," ");
00166       
00167           // This produces e[0] + (-1)*e[1], we have to canonize it in 2 steps 
00168           Theorem minus_eq_sum = d_rules->minusToPlus(e[0], e[1]);
00169       
00170           // The resulting sum expression
00171           Expr sum(minus_eq_sum.getRHS());
00172           
00173           // Canonise the sum[1]
00174           Theorem thm(canon(sum[1]));
00175           
00176           // If the sum[1] sayed the same, canonise the sum expression
00177           if(thm.getLHS() == thm.getRHS()) 
00178             result = canonThm(minus_eq_sum);
00179           // The sum changed; do the work
00180           else {
00181             
00182             // Substitute and try to canonise again
00183             Theorem sum_eq_canon = canonThm(substitutivityRule(sum, 1, thm));
00184             
00185             // Combine the results, and thats it
00186             result = transitivityRule(minus_eq_sum, sum_eq_canon);
00187           }
00188           
00189           // MINUS case break 
00190           break;
00191       
00192       }
00193   
00194       // e[0] * e[1]
00195       case MULT:
00196     
00197         // Canonise using the canonMult() proof rule
00198           result = d_rules->canonMult(e);
00199     
00200           // MULT case break
00201           break;
00202   
00203       // e[0] DIV e[1]
00204       case DIVIDE: {
00205   
00206           // Division by 0 is OK (total extension, protected by TCCs)
00207         // if (e[1].isRational() && e[1].getRational() == 0)
00208       // throw ArithException("Divide by 0 error in "+e.toString());
00209           if (e[1].getKind() == PLUS)
00210             throw ArithException("Divide by a PLUS expression not handled in" + e.toString());
00211           
00212           // Canonise the divite using canonDivide() method 
00213           result = d_rules->canonDivide(e);
00214           
00215           // DIVIDE case break
00216           break;  
00217       }
00218     
00219       // e[0]^e[1]
00220       case POW :
00221     
00222         // If we are raising to a constant rational call canonPowConst() proof rule
00223         if(e[1].isRational()) result = d_rules->canonPowConst(e);
00224         // Othervise just return the same expression
00225         else result = reflexivityRule(e);
00226     
00227         // POW case break
00228         break;
00229     
00230       default:
00231         
00232           // How did we get here (TODO: check if need to throw and exception), anyhow, just return the reflexivity rule
00233           result = reflexivityRule(e);
00234      
00235         // Default break 
00236           break;
00237     }
00238   
00239   // Finished with the canonisation, trace the result in the debug version
00240   TRACE("arith","canon => ",result," }");
00241   
00242   // Return the resulting theorem 
00243   return result;
00244 }
00245 
00246 Theorem TheoryArithNew::canonSimplify(const Expr& e) {
00247     // Trace the function on the arith trace flag
00248     TRACE("arith", "canonSimplify(", e, ") {");
00249   
00250     // Assert that all the kids must be already canonical
00251     //DebugAssert(kidsCanonical(e), "TheoryArithNew::canonSimplify("+e.toString()+")");
00252     // Assert that all the leaves must be simplified
00253     //DebugAssert(leavesAreSimp(e), "Expected leaves to be simplified");
00254   
00255     // Make a copy of the expression
00256     Expr tmp(e);
00257   
00258     // Canonise the expresion
00259     Theorem thm = canon(e);
00260     
00261     // If the new simplified expression has a find, use it (TODO: Check whats going on here)
00262     if(thm.getRHS().hasFind())
00263       thm = transitivityRule(thm, find(thm.getRHS()));
00264   
00265     // We shouldn't rely on simplification in this function anymore
00266     DebugAssert(thm.getRHS() == simplifyExpr(thm.getRHS()), "canonSimplify("+e.toString()+")\n"+"canon(e) = "+thm.getRHS().toString()+"\nsimplify(canon(e)) = "+simplifyExpr(thm.getRHS()).toString());
00267 
00268   // Trace the resulting theorem
00269   TRACE("arith", "canonSimplify =>", thm, " }");
00270   
00271   // Return the theorem
00272   return thm;
00273 }
00274 
00275 /*! accepts a theorem, canonizes it, applies iffMP and substituvity to
00276  *  derive the canonized thm
00277  */
00278 Theorem 
00279 TheoryArithNew::canonPred(const Theorem& thm) {
00280   vector<Theorem> thms;
00281   DebugAssert(thm.getExpr().arity() == 2,
00282               "TheoryArithNew::canonPred: bad theorem: "+thm.toString());
00283   Expr e(thm.getExpr());
00284   thms.push_back(canonSimplify(e[0]));
00285   thms.push_back(canonSimplify(e[1]));
00286   return iffMP(thm, substitutivityRule(e.getOp(), thms));
00287 }
00288 
00289 /** 
00290  * Accepts an equivalence theorem, canonizes the subexpressions, applies transitivity and substituvity to derive the canonized theorem.
00291  * 
00292  * @param thm equivalence theorem
00293  * @return the canonised expression theorem
00294  */
00295 Theorem TheoryArithNew::canonPredEquiv(const Theorem& thm) {
00296   
00297   vector<Theorem> thms; // Vector to hold the simplified versions of e[0] and e[1]
00298   
00299   // Arity of the expression should be 2
00300   DebugAssert(thm.getRHS().arity() == 2, "TheoryArithNew::canonPredEquiv: bad theorem: " + thm.toString());
00301   
00302   // Get the expression of the theorem
00303   Expr e(thm.getRHS());
00304   
00305   // Simplify both sides of the equivalence
00306   thms.push_back(canonSimplify(e[0]));
00307   thms.push_back(canonSimplify(e[1]));
00308   
00309   // Return the theorem substituting e[0] and e[1] with their simplified versions
00310   return transitivityRule(thm, substitutivityRule(e.getOp(), thms));
00311 }
00312 
00313 /*! accepts an equivalence theorem whose RHS is a conjunction,
00314  *  canonizes it, applies iffMP and substituvity to derive the
00315  *  canonized thm
00316  */
00317 Theorem
00318 TheoryArithNew::canonConjunctionEquiv(const Theorem& thm) {
00319   vector<Theorem> thms;
00320   return thm;
00321 }
00322 
00323 /*! Psuedo-code for doSolve. (Input is an equation) (output is a Theorem)
00324  *  -# translate e to the form e' = 0
00325  *  -# if (e'.isRational()) then {if e' != 0 return false else true}
00326  *  -# a for loop checks if all the variables are integers. 
00327  *      - if not isolate a suitable real variable and call processRealEq().
00328  *      - if all variables are integers then isolate suitable variable 
00329  *         and call processIntEq(). 
00330  */
00331 Theorem TheoryArithNew::doSolve(const Theorem& thm)
00332 { 
00333   const Expr& e = thm.getExpr();
00334   TRACE("arith eq","doSolve(",e,") {");
00335   DebugAssert(thm.isRewrite(), "thm = "+thm.toString());
00336   Theorem eqnThm;
00337   vector<Theorem> thms;  
00338   // Move LHS to the RHS, if necessary
00339   if(e[0].isRational() && e[0].getRational() == 0)
00340     eqnThm = thm;
00341   else {
00342     eqnThm = iffMP(thm, d_rules->rightMinusLeft(e));
00343     eqnThm = canonPred(eqnThm);
00344   }
00345   // eqnThm is of the form 0 = e'
00346   // 'right' is of the form e'
00347   Expr right = eqnThm.getRHS();
00348   // Check for trivial equation
00349   if (right.isRational()) {
00350     Theorem result = iffMP(eqnThm, d_rules->constPredicate(eqnThm.getExpr()));
00351     TRACE("arith eq","doSolve => ",result," }");
00352     return result;
00353   }
00354 
00355   //normalize
00356   eqnThm = iffMP(eqnThm, normalize(eqnThm.getExpr()));
00357   right = eqnThm.getRHS();
00358   
00359   //eqn is of the form 0 = e' and is normalized where 'right' denotes e'
00360   //FIXME: change processRealEq to accept equations as well instead of theorems
00361   if(!isInteger(right)) {
00362     Theorem res;
00363     try {
00364       res = processRealEq(eqnThm);
00365     } catch(ArithException& e) {
00366       res = symmetryRule(eqnThm); // Flip to e' = 0
00367       TRACE("arith eq", "doSolve: failed to solve an equation: ", e, "");
00368       IF_DEBUG(debugger.counter("FAILED to solve real equalities")++);
00369       setIncomplete("Non-linear arithmetic inequalities");
00370     }
00371     IF_DEBUG(debugger.counter("solved real equalities")++);
00372     TRACE("arith eq", "doSolve [real] => ", res, " }");
00373     return res;
00374   }
00375   else {
00376     Theorem res = processIntEq(eqnThm);
00377     IF_DEBUG(debugger.counter("solved int equalities")++);
00378     TRACE("arith eq", "doSolve [int] => ", res, " }");
00379     return res;
00380   }
00381 }
00382 
00383 /*! pick a monomial for the input equation. This function is used only
00384  *  if the equation is an integer equation. Choose the monomial with
00385  *  the smallest absolute value of coefficient.
00386  */
00387 Expr
00388 TheoryArithNew::pickIntEqMonomial(const Expr& right)
00389 {
00390   DebugAssert(isPlus(right) && right.arity() > 2,
00391               "TheoryArithNew::pickIntEqMonomial right is wrong :-): " +
00392               right.toString());
00393   DebugAssert(right[0].isRational(),
00394               "TheoryArithNew::pickIntEqMonomial. right[0] must be const" +
00395               right.toString());
00396   DebugAssert(isInteger(right),
00397               "TheoryArithNew::pickIntEqMonomial: right is of type int: " +
00398               right.toString());
00399   //right is of the form "C + a1x1 + ... + anxn". min is initialized
00400   //to a1
00401   Expr::iterator i = right.begin();
00402   i++; //skip 'C'
00403   Rational min = isMult(*i) ? abs((*i)[0].getRational()) : 1;
00404   Expr pickedMon = *i;
00405   for(Expr::iterator iend = right.end(); i != iend; ++i) {
00406     Rational coeff = isMult(*i) ? abs((*i)[0].getRational()) : 1;
00407     if(min > coeff) {
00408       min = coeff;
00409       pickedMon = *i;
00410     }
00411   }
00412   return pickedMon;
00413 }
00414 
00415 /*! input is e1=e2<==>0=e' Theorem and some of the vars in e' are of
00416  * type REAL. isolate one of them and send back to framework. output
00417  * is "e1=e2 <==> var = e''" Theorem.
00418  */
00419 Theorem 
00420 TheoryArithNew::processRealEq(const Theorem& eqn)
00421 {
00422   Expr right = eqn.getRHS();
00423   // Find variable to isolate and store it in left.  Pick the largest
00424   // (according to the total ordering) variable.  FIXME: change from
00425   // total ordering to the ordering we devised for inequalities.
00426 
00427   // TODO: I have to pick a variable that appears as a variable in the
00428   // term but does not appear as a variable anywhere else.  The variable
00429   // must appear as a single leaf and not in a MULT expression with some
00430   // other variables and nor in a POW expression.
00431 
00432   bool found = false;
00433   
00434   Expr left;
00435   
00436   if (isPlus(right))  {
00437     for(int i = right.arity()-1; i>=0; --i) {
00438       Expr c = right[i];
00439       if(isRational(c))
00440         continue;
00441       if(!isInteger(c))  {
00442         if (isLeaf(c) || 
00443             ((isMult(c) && c.arity() == 2 && isLeaf(c[1])))) {
00444           int numoccurs = 0;
00445           Expr leaf = isLeaf(c) ? c : c[1];
00446           for (int j = 0; j < right.arity(); ++j) {
00447             if (j!= i
00448     && isLeafIn(leaf, right[j])
00449     // && leaf.subExprOf(right[j])
00450     ) {
00451               numoccurs++;
00452               break;
00453             }
00454           }
00455           if (!numoccurs) {
00456             left = c;
00457             found = true;
00458             break;
00459           }
00460         }
00461       }
00462     }
00463   }
00464   else if ((isMult(right) && right.arity() == 2 && isLeaf(right[1])) ||
00465            isLeaf(right)) {
00466     left = right;
00467     found = true;
00468   }
00469   
00470   if (!found) {
00471     // TODO:
00472     // throw an arithmetic exception that this cannot be done.
00473     throw 
00474       ArithException("Can't find a leaf for solve in "+eqn.toString());
00475   }
00476 
00477   Rational r = -1;
00478   if (isMult(left))  {
00479     DebugAssert(left.arity() == 2, "only leaf should be chosen as lhs");
00480     DebugAssert(left[0].getRational() != 0, "left = "+left.toString());
00481     r = -1/left[0].getRational();
00482     left = left[1];
00483   }
00484 
00485   DebugAssert(isReal(getBaseType(left)) && !isInteger(left),
00486               "TheoryArithNew::ProcessRealEq: left is integer:\n left = "
00487         +left.toString());
00488   // Normalize equation so that coefficient of the monomial
00489   // corresponding to "left" in eqn[1] is -1
00490   Theorem result(iffMP(eqn,
00491            d_rules->multEqn(eqn.getLHS(), eqn.getRHS(), rat(r))));
00492   result = canonPred(result);
00493 
00494   // Isolate left
00495   result = iffMP(result, d_rules->plusPredicate(result.getLHS(),
00496             result.getRHS(), left, EQ));
00497   result = canonPred(result);
00498   TRACE("arith","processRealEq => ",result," }");
00499   return result;
00500 }
00501 
00502 /*!
00503  * \param eqn is a single equation 0 = e
00504  * \return an equivalent Theorem (x = t AND 0 = e'), or just x = t
00505  */
00506 Theorem
00507 TheoryArithNew::processSimpleIntEq(const Theorem& eqn)
00508 {
00509   TRACE("arith eq", "processSimpleIntEq(", eqn.getExpr(), ") {");
00510   DebugAssert(eqn.isRewrite(),
00511               "TheoryArithNew::processSimpleIntEq: eqn must be equality" +
00512               eqn.getExpr().toString());
00513 
00514   Expr right = eqn.getRHS();
00515 
00516   DebugAssert(eqn.getLHS().isRational() && 0 == eqn.getLHS().getRational(),
00517               "TheoryArithNew::processSimpleIntEq: LHS must be 0:\n" +
00518               eqn.getExpr().toString());
00519 
00520   //recall that 0 = c case is already handled in doSolve() function.
00521   if(isMult(right)) {
00522     //here we take care of special case 0=c.x
00523     Expr c,x;
00524     separateMonomial(right, c, x);
00525     Theorem isIntx(isIntegerThm(x));
00526     DebugAssert(!isIntx.isNull(), "right = "+right.toString());
00527     Theorem res(iffMP(eqn, d_rules->intVarEqnConst(eqn.getExpr(), isIntx)));
00528     TRACE("arith eq", "processSimpleIntEq[0 = a*x] => ", res, " }");
00529     return res;
00530   } else if(isPlus(right)) {
00531     if(2 == right.arity()) {
00532       //we take care of special cases like 0 = c + a.x, 0 = c + x,
00533       Expr c,x;
00534       separateMonomial(right[1], c, x);
00535       Theorem isIntx(isIntegerThm(x));
00536       DebugAssert(!isIntx.isNull(), "right = "+right.toString()
00537       +"\n x = "+x.toString());
00538       Theorem res(iffMP(eqn, d_rules->intVarEqnConst(eqn.getExpr(), isIntx)));
00539       TRACE("arith eq", "processSimpleIntEq[0 = c + a*x] => ", res, " }");
00540       return res;
00541     }
00542     DebugAssert(right.arity() > 2,
00543                 "TheoryArithNew::processSimpleIntEq: RHS is not in correct form:"
00544                 +eqn.getExpr().toString());
00545     // Pick a suitable monomial. isolated can be of the form x, a.x, -a.x
00546     Expr isolated = pickIntEqMonomial(right);
00547     TRACE("arith eq", "processSimpleIntEq: isolated = ", isolated, "");
00548 
00549     // First, we compute the 'sign factor' with which to multiply the
00550     // eqn.  if the coeff of isolated is positive (i.e. 'isolated' is
00551     // of the form x or a.x where a>0 ) then r must be -1 and if coeff
00552     // of 'isolated' is negative, r=1.
00553     Rational r = isMult(isolated) ?
00554       ((isolated[0].getRational() > 0) ? -1 : 1) : -1;
00555     Theorem result;
00556     if (-1 == r) {
00557       // r=-1 and hence 'isolated' is 'x' or 'a.x' where a is
00558       // positive.  modify eqn (0=e') to the equation (0=canon(-1*e'))
00559       result = iffMP(eqn, d_rules->multEqn(eqn.getLHS(), right, rat(r)));
00560       result = canonPred(result);
00561       Expr e = result.getRHS();
00562 
00563       // Isolate the 'isolated'
00564       result = iffMP(result,
00565          d_rules->plusPredicate(result.getLHS(),result.getRHS(),
00566               isolated, EQ));
00567     } else {
00568       //r is 1 and hence isolated is -a.x. Make 'isolated' positive.
00569       const Rational& minusa = isolated[0].getRational();
00570       Rational a = -1*minusa;
00571       isolated = (a == 1)? isolated[1] : rat(a) * isolated[1];
00572       
00573       // Isolate the 'isolated'
00574       result = iffMP(eqn, d_rules->plusPredicate(eqn.getLHS(), 
00575              right,isolated,EQ));
00576     }
00577     // Canonize the result
00578     result = canonPred(result);
00579         
00580     //if isolated is 'x' or 1*x, then return result else continue processing.
00581     if(!isMult(isolated) || isolated[0].getRational() == 1) {   
00582       TRACE("arith eq", "processSimpleIntEq[x = rhs] => ", result, " }");
00583       return result;
00584     } else {
00585       DebugAssert(isMult(isolated) && isolated[0].getRational() >= 2,
00586                   "TheoryArithNew::processSimpleIntEq: isolated must be mult "
00587       "with coeff >= 2.\n isolated = " + isolated.toString());
00588 
00589       // Compute IS_INTEGER() for lhs and rhs
00590       Expr lhs = result.getLHS();
00591       Expr rhs = result.getRHS();
00592       Expr a, x;
00593       separateMonomial(lhs, a, x);
00594       Theorem isIntLHS = isIntegerThm(x);
00595       vector<Theorem> isIntRHS;
00596       if(!isPlus(rhs)) { // rhs is a MULT
00597   Expr c, v;
00598   separateMonomial(rhs, c, v);
00599   isIntRHS.push_back(isIntegerThm(v));
00600   DebugAssert(!isIntRHS.back().isNull(), "");
00601       } else { // rhs is a PLUS
00602   DebugAssert(isPlus(rhs), "rhs = "+rhs.toString());
00603   DebugAssert(rhs.arity() >= 2, "rhs = "+rhs.toString());
00604   Expr::iterator i=rhs.begin(), iend=rhs.end();
00605   ++i; // Skip the free constant
00606   for(; i!=iend; ++i) {
00607     Expr c, v;
00608     separateMonomial(*i, c, v);
00609     isIntRHS.push_back(isIntegerThm(v));
00610     DebugAssert(!isIntRHS.back().isNull(), "");
00611   }
00612       }
00613       // Derive (EXISTS (x:INT): x = t2 AND 0 = t3)
00614       result = d_rules->eqElimIntRule(result, isIntLHS, isIntRHS);
00615       // Skolemize the quantifier
00616       result = getCommonRules()->skolemize(result);
00617       // Canonize t2 and t3 generated by this rule
00618       DebugAssert(result.getExpr().isAnd() && result.getExpr().arity() == 2,
00619       "processSimpleIntEq: result = "+result.getExpr().toString());
00620       Theorem thm1 = canonPred(getCommonRules()->andElim(result, 0));
00621       Theorem thm2 = canonPred(getCommonRules()->andElim(result, 1));
00622       Theorem newRes = getCommonRules()->andIntro(thm1, thm2);
00623       if(newRes.getExpr() != result.getExpr()) result = newRes;
00624       
00625       TRACE("arith eq", "processSimpleIntEq => ", result, " }");
00626       return result;
00627     }
00628   } else {
00629     // eqn is 0 = x.  Flip it and return
00630     Theorem result = symmetryRule(eqn);
00631     TRACE("arith eq", "processSimpleIntEq[x = 0] => ", result, " }");
00632     return result;
00633   }
00634 }
00635 
00636 /*! input is e1=e2<==>0=e' Theorem and all of the vars in e' are of
00637  * type INT. isolate one of them and send back to framework. output
00638  * is "e1=e2 <==> var = e''" Theorem and some associated equations in
00639  * solved form.
00640  */
00641 Theorem 
00642 TheoryArithNew::processIntEq(const Theorem& eqn)
00643 {
00644   TRACE("arith eq", "processIntEq(", eqn.getExpr(), ") {");
00645   // Equations in the solved form.  
00646   std::vector<Theorem> solvedAndNewEqs;
00647   Theorem newEq(eqn), result;
00648   bool done(false);
00649   do {
00650     result = processSimpleIntEq(newEq);
00651     // 'result' is of the from (x1=t1)  AND 0=t'
00652     if(result.isRewrite()) {
00653       solvedAndNewEqs.push_back(result);
00654       done = true;
00655     }
00656     else if(!result.getExpr().isFalse()) {
00657       DebugAssert(result.getExpr().isAnd() && result.getExpr().arity() == 2,
00658       "TheoryArithNew::processIntEq("+eqn.getExpr().toString()
00659       +")\n result = "+result.getExpr().toString());
00660       solvedAndNewEqs.push_back(getCommonRules()->andElim(result, 0));
00661       newEq = getCommonRules()->andElim(result, 1);
00662     } else
00663       done = true;
00664   } while(!done);
00665   Theorem res;
00666   if(result.getExpr().isFalse()) res = result;
00667   else res = solvedForm(solvedAndNewEqs);
00668   TRACE("arith eq", "processIntEq => ", res.getExpr(), " }");
00669   return res;
00670 }
00671 
00672 /*!
00673  * Takes a vector of equations and for every equation x_i=t_i
00674  * substitutes t_j for x_j in t_i for all j>i.  This turns the system
00675  * of equations into a solved form.
00676  *
00677  * Assumption: variables x_j may appear in the RHS terms t_i ONLY for
00678  * i<j, but not for i>=j.
00679  */
00680 Theorem
00681 TheoryArithNew::solvedForm(const vector<Theorem>& solvedEqs) 
00682 {
00683   DebugAssert(solvedEqs.size() > 0, "TheoryArithNew::solvedForm()");
00684 
00685   // Trace code
00686   TRACE_MSG("arith eq", "TheoryArithNew::solvedForm:solvedEqs(\n  [");
00687   IF_DEBUG(if(debugger.trace("arith eq")) {
00688     for(vector<Theorem>::const_iterator j = solvedEqs.begin(),
00689     jend=solvedEqs.end(); j!=jend;++j)
00690       TRACE("arith eq", "", j->getExpr(), ",\n   ");
00691   });
00692   TRACE_MSG("arith eq", "  ]) {");
00693   // End of Trace code
00694   
00695   vector<Theorem>::const_reverse_iterator
00696     i = solvedEqs.rbegin(),
00697     iend = solvedEqs.rend();
00698   // Substitution map: a variable 'x' is mapped to a Theorem x=t.
00699   // This map accumulates the resulting solved form.
00700   ExprMap<Theorem> subst;
00701   for(; i!=iend; ++i) {
00702     if(i->isRewrite()) {
00703       Theorem thm = substAndCanonize(*i, subst);
00704       TRACE("arith eq", "solvedForm: subst["+i->getLHS().toString()+"] = ",
00705       thm.getExpr(), "");
00706       subst[i->getLHS()] = thm;
00707     }
00708     else {
00709       // This is the FALSE case: just return the contradiction
00710       DebugAssert(i->getExpr().isFalse(),
00711       "TheoryArithNew::solvedForm: an element of solvedEqs must "
00712       "be either EQ or FALSE: "+i->toString());
00713       return *i;
00714     }
00715   }
00716   // Now we've collected the solved form in 'subst'.  Wrap it up into
00717   // a conjunction and return.
00718   vector<Theorem> thms;
00719   for(ExprMap<Theorem>::iterator i=subst.begin(), iend=subst.end();
00720       i!=iend; ++i)
00721     thms.push_back(i->second);
00722   return getCommonRules()->andIntro(thms);
00723 }
00724 
00725 
00726 /*!
00727  * ASSUMPTION: 't' is a fully canonized arithmetic term, and every
00728  * element of subst is a fully canonized equation of the form x=e,
00729  * indexed by the LHS variable.
00730  */
00731 
00732 Theorem
00733 TheoryArithNew::substAndCanonize(const Expr& t, ExprMap<Theorem>& subst)
00734 {
00735   TRACE("arith eq", "substAndCanonize(", t, ") {");
00736   // Quick and dirty check: return immediately if subst is empty
00737   if(subst.empty()) {
00738     Theorem res(reflexivityRule(t));
00739     TRACE("arith eq", "substAndCanonize[subst empty] => ", res, " }");
00740     return res;
00741   }
00742   // Check if we can substitute 't' directly
00743   ExprMap<Theorem>::iterator i = subst.find(t), iend=subst.end();
00744   if(i!=iend) {
00745     TRACE("arith eq", "substAndCanonize[subst] => ", i->second, " }");
00746     return i->second;
00747   }
00748   // The base case: t is an i-leaf
00749   if(isLeaf(t)) {
00750     Theorem res(reflexivityRule(t));
00751     TRACE("arith eq", "substAndCanonize[i-leaf] => ", res, " }");
00752     return res;
00753   }
00754   // 't' is an arithmetic term; recurse into the children
00755   vector<Theorem> thms;
00756   vector<unsigned> changed;
00757   for(unsigned j=0, jend=t.arity(); j!=jend; ++j) {
00758     Theorem thm = substAndCanonize(t[j], subst);
00759     if(thm.getRHS() != t[j]) {
00760       thm = canonThm(thm);
00761       thms.push_back(thm);
00762       changed.push_back(j);
00763     }
00764   }
00765   // Do the actual substitution and canonize the result
00766   Theorem res;
00767   if(thms.size() > 0) {
00768     res = substitutivityRule(t, changed, thms);
00769     res = canonThm(res);
00770   }
00771   else
00772     res = reflexivityRule(t);
00773   TRACE("arith eq", "substAndCanonize => ", res, " }");
00774   return res;
00775 }
00776 
00777 
00778 /*!
00779  *  ASSUMPTION: 't' is a fully canonized equation of the form x = t,
00780  *  and so is every element of subst, indexed by the LHS variable.
00781  */
00782 
00783 Theorem
00784 TheoryArithNew::substAndCanonize(const Theorem& eq, ExprMap<Theorem>& subst)
00785 {
00786   // Quick and dirty check: return immediately if subst is empty
00787   if(subst.empty()) return eq;
00788 
00789   DebugAssert(eq.isRewrite(), "TheoryArithNew::substAndCanonize: t = "
00790         +eq.getExpr().toString());
00791   const Expr& t = eq.getRHS();
00792   // Do the actual substitution in the term t
00793   Theorem thm = substAndCanonize(t, subst);
00794   // Substitution had no result: return the original equation
00795   if(thm.getRHS() == t) return eq;
00796   // Otherwise substitute the result into the equation
00797   vector<Theorem> thms;
00798   vector<unsigned> changed;
00799   thms.push_back(thm);
00800   changed.push_back(1);
00801   return iffMP(eq, substitutivityRule(eq.getExpr(), changed, thms));
00802 }
00803 
00804 
00805 void TheoryArithNew::updateStats(const Rational& c, const Expr& v)
00806 {
00807   TRACE("arith ineq", "updateStats("+c.toString()+", ", v, ")");
00808   if(c > 0) {
00809     if(d_countRight.count(v) > 0) d_countRight[v] = d_countRight[v] + 1;
00810     else d_countRight[v] = 1;
00811   }
00812   else
00813     if(d_countLeft.count(v) > 0) d_countLeft[v] = d_countLeft[v] + 1;
00814     else d_countLeft[v] = 1;
00815 }
00816 
00817 
00818 void TheoryArithNew::updateStats(const Expr& monomial)
00819 {
00820   Expr c, m;
00821   separateMonomial(monomial, c, m);
00822   updateStats(c.getRational(), m);
00823 }
00824 
00825 
00826 void TheoryArithNew::addToBuffer(const Theorem& thm) {
00827   // First, turn the inequality into 0 < rhs
00828   // FIXME: check if this can be optimized away
00829   Theorem result(thm);
00830   const Expr& e = thm.getExpr();
00831   // int kind = e.getKind();
00832   if(!(e[0].isRational() && e[0].getRational() == 0)) {
00833     result = iffMP(result, d_rules->rightMinusLeft(e));
00834     result = canonPred(result);
00835   }
00836   TRACE("arith ineq", "addToBuffer(", result.getExpr(), ")");
00837   // Push it into the buffer
00838   d_buffer.push_back(thm);
00839 
00840   // Collect the statistics about variables
00841   const Expr& rhs = thm.getExpr()[1];
00842   if(isPlus(rhs))
00843     for(Expr::iterator i=rhs.begin(), iend=rhs.end(); i!=iend; ++i)
00844       updateStats(*i);
00845   else // It's a monomial
00846     updateStats(rhs);
00847 }
00848 
00849 
00850 Expr TheoryArithNew::computeNormalFactor(const Expr& right, NormalizationType type) {
00851     // Strategy: compute f = lcm(d1...dn)/gcd(c1...cn), where the RHS is of the form c1/d1*x1 + ... + cn/dn*xn or a rational
00852   
00853     // The expression must be canonised, i.e. it is a sum or a rational
00854     DebugAssert(isRational(right) || isPlus(right),"TheoryArithNew::computeNormalFactor: " + right.toString() + "wrong kind");
00855   
00856   // The factor we will compute
00857     Rational factor;
00858     
00859     // Sign of the factor 
00860     int sign = 0;
00861     
00862     // In case the expression is a PLUS expression find the gcd
00863     if(isPlus(right)) {
00864       
00865       // Vector of numerators and denominators
00866       vector<Rational> nums, denoms;
00867       
00868       // Pass throught the sum and pick up the integers
00869       for(int i = 0, iend = right.arity(); i < iend; i ++) {
00870       
00871           Rational c(1); // The rational constant of the current summand
00872                 
00873           // If rational or multiplicatino set c to be the apropriate rational
00874           switch(right[i].getKind()) {
00875             
00876             case RATIONAL_EXPR: 
00877             
00878               // Sum term is just a rational, so just add it 
00879               c = abs(right[i].getRational()); 
00880               break;
00881               
00882               case MULT:          
00883                 
00884                 // Sum term is involves multiplication
00885                 c = right[i][0].getRational();
00886                 
00887                 // If this is the first one (sign is still not set) set the sign depending on the sign of the term
00888                 if (!sign) {
00889                   
00890                   // If the normalization type is NORMALIZE_UNIT just return the invese of the first
00891                   if (type == NORMALIZE_UNIT) return rat(1/c);
00892                   
00893                   // Set the sign otherwise
00894                   sign = (c > 0 ? 1 : -1);                  
00895                 }
00896                 
00897                 // Constant should be positive for lcd and gcd computation
00898                 c = abs(c);
00899                 
00900                 break;
00901                 
00902             default: // Stays 1 or everything else
00903               
00904               break;
00905           }
00906           
00907           // Add the demoninator and the numerator to the apropriate vectors
00908           nums.push_back(c.getNumerator());
00909           denoms.push_back(c.getDenominator());   
00910       }
00911     
00912     // Compute the gcd of all the numerators
00913       Rational gcd_nums = gcd(nums);
00914     
00915       // x/0 == 0 in our model, as a total extension of arithmetic.  The
00916       // particular value of x/0 is irrelevant, since the DP is guarded
00917       // by the top-level TCCs, and it is safe to return any value in
00918       // cases when terms are undefined.
00919       factor = (gcd_nums == 0)? 0 : (sign * lcm(denoms) / gcd_nums);
00920     } else 
00921       // The expression is a rational, just return 1
00922       factor = 1;
00923   
00924   // Return the reational expression representing the factor
00925   return rat(factor);
00926 }
00927 
00928 
00929 bool TheoryArithNew::lessThanVar(const Expr& isolatedMonomial, const Expr& var2) 
00930 {
00931   DebugAssert(!isRational(var2) && !isRational(isolatedMonomial),
00932               "TheoryArithNew::findMaxVar: isolatedMonomial cannot be rational" + 
00933               isolatedMonomial.toString());
00934   Expr c, var0, var1;
00935   separateMonomial(isolatedMonomial, c, var0);
00936   separateMonomial(var2, c, var1);
00937   return var0 < var1;
00938 }
00939 
00940 
00941 void TheoryArithNew::separateMonomial(const Expr& e, Expr& c, Expr& var) {
00942   TRACE("separateMonomial", "separateMonomial(", e, ")");
00943   DebugAssert(!isPlus(e), 
00944         "TheoryArithNew::separateMonomial(e = "+e.toString()+")");
00945   if(isMult(e)) {
00946     DebugAssert(e.arity() >= 2,
00947     "TheoryArithNew::separateMonomial(e = "+e.toString()+")");
00948     c = e[0];
00949     if(e.arity() == 2) var = e[1];
00950     else {
00951       vector<Expr> kids = e.getKids();
00952       kids[0] = rat(1);
00953       var = multExpr(kids);
00954     }
00955   } else {
00956     c = rat(1);
00957     var = e;
00958   }
00959   DebugAssert(c.isRational(), "TheoryArithNew::separateMonomial(e = "
00960         +e.toString()+", c = "+c.toString()+")");
00961   DebugAssert(!isMult(var) || (var[0].isRational() && var[0].getRational()==1),
00962         "TheoryArithNew::separateMonomial(e = "
00963         +e.toString()+", var = "+var.toString()+")");
00964 }
00965 
00966 
00967 //returns true if e1 < e2, else false(i.e e2 < e1 or e1,e2 are not
00968 //comparable)
00969 bool TheoryArithNew::VarOrderGraph::lessThan(const Expr& e1, const Expr& e2) 
00970 {
00971   d_cache.clear();
00972   //returns true if e1 is in the subtree rooted at e2 implying e1 < e2
00973   return dfs(e1,e2);
00974 }
00975 
00976 //returns true if e1 is in the subtree rooted at e2 implying e1 < e2
00977 bool TheoryArithNew::VarOrderGraph::dfs(const Expr& e1, const Expr& e2)
00978 {
00979   if(e1 == e2)
00980     return true;
00981   if(d_cache.count(e2) > 0)
00982     return false;
00983   if(d_edges.count(e2) == 0)
00984     return false;
00985   d_cache[e2] = true;
00986   vector<Expr>& e2Edges = d_edges[e2];
00987   vector<Expr>::iterator i = e2Edges.begin();
00988   vector<Expr>::iterator iend = e2Edges.end();
00989   //if dfs finds e1 then i != iend else i is equal to iend
00990   for(; i != iend && !dfs(e1,*i); ++i);
00991   return (i != iend);
00992 }
00993 
00994 
00995 void TheoryArithNew::VarOrderGraph::selectSmallest(vector<Expr>& v1,
00996                                                vector<Expr>& v2) 
00997 {
00998   int v1Size = v1.size();
00999   vector<bool> v3(v1Size);
01000   for(int j=0; j < v1Size; ++j)
01001     v3[j] = false;
01002 
01003   for(int j=0; j < v1Size; ++j) {
01004     if(v3[j]) continue;
01005     for(int i =0; i < v1Size; ++i) {
01006       if((i == j) || v3[i]) 
01007         continue;
01008       if(lessThan(v1[i],v1[j])) {
01009         v3[j] = true;
01010         break;
01011       }
01012     }
01013   }
01014   vector<Expr> new_v1;
01015 
01016   for(int j = 0; j < v1Size; ++j) 
01017     if(!v3[j]) v2.push_back(v1[j]);
01018     else new_v1.push_back(v1[j]);
01019   v1 = new_v1;
01020 }
01021 
01022 
01023 void TheoryArithNew::VarOrderGraph::selectLargest(const vector<Expr>& v1,
01024                                                vector<Expr>& v2) 
01025 {
01026   int v1Size = v1.size();
01027   vector<bool> v3(v1Size);
01028   for(int j=0; j < v1Size; ++j)
01029     v3[j] = false;
01030 
01031   for(int j=0; j < v1Size; ++j) {
01032     if(v3[j]) continue;
01033     for(int i =0; i < v1Size; ++i) {
01034       if((i == j) || v3[i]) 
01035         continue;
01036       if(lessThan(v1[j],v1[i])) {
01037         v3[j] = true;
01038         break;
01039       }
01040     }
01041   }
01042   
01043   for(int j = 0; j < v1Size; ++j) 
01044     if(!v3[j]) v2.push_back(v1[j]);
01045 }
01046 
01047 ///////////////////////////////////////////////////////////////////////////////
01048 // TheoryArithNew Public Methods                                                //
01049 ///////////////////////////////////////////////////////////////////////////////
01050 
01051 
01052 TheoryArithNew::TheoryArithNew(TheoryCore* core)
01053   : TheoryArith(core, "ArithmeticNew"),
01054     d_diseq(core->getCM()->getCurrentContext()),
01055     d_diseqIdx(core->getCM()->getCurrentContext(), 0, 0),
01056     d_inModelCreation(core->getCM()->getCurrentContext(), false, 0),
01057     d_freeConstDB(core->getCM()->getCurrentContext()),
01058     d_buffer(core->getCM()->getCurrentContext()),
01059     // Initialize index to 0 at scope 0
01060     d_bufferIdx(core->getCM()->getCurrentContext(), 0, 0),
01061     d_bufferThres(&(core->getFlags()["ineq-delay"].getInt())),
01062     d_countRight(core->getCM()->getCurrentContext()),
01063     d_countLeft(core->getCM()->getCurrentContext()),
01064     d_sharedTerms(core->getCM()->getCurrentContext()),
01065     d_sharedVars(core->getCM()->getCurrentContext()),
01066     consistent(core->getCM()->getCurrentContext()),
01067     lowerBound(core->getCM()->getCurrentContext()),
01068     upperBound(core->getCM()->getCurrentContext()),
01069     beta(core->getCM()->getCurrentContext()),
01070     assertedExprCount(core->getCM()->getCurrentContext()),
01071     integer_lemma(core->getCM()->getCurrentContext())
01072 {
01073   IF_DEBUG(d_diseq.setName("CDList[TheoryArithNew::d_diseq]"));
01074   IF_DEBUG(d_buffer.setName("CDList[TheoryArithNew::d_buffer]"));
01075   IF_DEBUG(d_bufferIdx.setName("CDList[TheoryArithNew::d_bufferIdx]"));
01076 
01077   getEM()->newKind(REAL, "_REAL", true);
01078   getEM()->newKind(INT, "_INT", true);
01079   getEM()->newKind(SUBRANGE, "_SUBRANGE", true);
01080 
01081   getEM()->newKind(UMINUS, "_UMINUS");
01082   getEM()->newKind(PLUS, "_PLUS");
01083   getEM()->newKind(MINUS, "_MINUS");
01084   getEM()->newKind(MULT, "_MULT");
01085   getEM()->newKind(DIVIDE, "_DIVIDE");
01086   getEM()->newKind(POW, "_POW");
01087   getEM()->newKind(INTDIV, "_INTDIV");
01088   getEM()->newKind(MOD, "_MOD");
01089   getEM()->newKind(LT, "_LT");
01090   getEM()->newKind(LE, "_LE");
01091   getEM()->newKind(GT, "_GT");
01092   getEM()->newKind(GE, "_GE");
01093   getEM()->newKind(IS_INTEGER, "_IS_INTEGER");
01094   getEM()->newKind(NEGINF, "_NEGINF");
01095   getEM()->newKind(POSINF, "_POSINF");
01096   getEM()->newKind(DARK_SHADOW, "_DARK_SHADOW");
01097   getEM()->newKind(GRAY_SHADOW, "_GRAY_SHADOW");
01098 
01099   getEM()->newKind(REAL_CONST, "_REAL_CONST");
01100 
01101   vector<int> kinds;
01102   kinds.push_back(REAL);
01103   kinds.push_back(INT);
01104   kinds.push_back(SUBRANGE);
01105   kinds.push_back(IS_INTEGER);
01106   kinds.push_back(UMINUS);
01107   kinds.push_back(PLUS);
01108   kinds.push_back(MINUS);
01109   kinds.push_back(MULT);
01110   kinds.push_back(DIVIDE);
01111   kinds.push_back(POW);
01112   kinds.push_back(INTDIV);
01113   kinds.push_back(MOD);
01114   kinds.push_back(LT);
01115   kinds.push_back(LE);
01116   kinds.push_back(GT);
01117   kinds.push_back(GE);
01118   kinds.push_back(RATIONAL_EXPR);
01119   kinds.push_back(NEGINF);
01120   kinds.push_back(POSINF);
01121   kinds.push_back(DARK_SHADOW);
01122   kinds.push_back(GRAY_SHADOW);
01123   kinds.push_back(REAL_CONST);
01124 
01125   registerTheory(this, kinds, true);
01126 
01127   d_rules = createProofRules();
01128 
01129   d_realType = Type(getEM()->newLeafExpr(REAL));
01130   d_intType  = Type(getEM()->newLeafExpr(INT));
01131   consistent = SATISFIABLE;
01132   assertedExprCount = 0;
01133 }
01134 
01135 
01136 // Destructor: delete the proof rules class if it's present
01137 TheoryArithNew::~TheoryArithNew() {
01138   if(d_rules != NULL) delete d_rules;
01139   // Clear the inequality databases
01140   for(ExprMap<CDList<Ineq> *>::iterator i=d_inequalitiesRightDB.begin(),
01141   iend=d_inequalitiesRightDB.end(); i!=iend; ++i) {
01142     delete (i->second);
01143     free(i->second);
01144   }
01145   for(ExprMap<CDList<Ineq> *>::iterator i=d_inequalitiesLeftDB.begin(),
01146   iend=d_inequalitiesLeftDB.end(); i!=iend; ++i) {
01147     delete (i->second);
01148     free (i->second);
01149   }
01150 }
01151 
01152 void TheoryArithNew::collectVars(const Expr& e, vector<Expr>& vars,
01153             set<Expr>& cache) {
01154   // Check the cache first
01155   if(cache.count(e) > 0) return;
01156   // Not processed yet.  Cache the expression and proceed
01157   cache.insert(e);
01158   if(isLeaf(e)) vars.push_back(e);
01159   else
01160     for(Expr::iterator i=e.begin(), iend=e.end(); i!=iend; ++i)
01161       collectVars(*i, vars, cache);
01162 }
01163 
01164 void
01165 TheoryArithNew::processFiniteInterval(const Theorem& alphaLEax,
01166            const Theorem& bxLEbeta) {
01167   const Expr& ineq1(alphaLEax.getExpr());
01168   const Expr& ineq2(bxLEbeta.getExpr());
01169   DebugAssert(isLE(ineq1), "TheoryArithNew::processFiniteInterval: ineq1 = "
01170         +ineq1.toString());
01171   DebugAssert(isLE(ineq2), "TheoryArithNew::processFiniteInterval: ineq2 = "
01172         +ineq2.toString());
01173   // If the inequalities are not integer, just return (nothing to do)
01174   if(!isInteger(ineq1[0])
01175      || !isInteger(ineq1[1])
01176      || !isInteger(ineq2[0])
01177      || !isInteger(ineq2[1]))
01178     return;
01179 
01180   const Expr& ax = ineq1[1];
01181   const Expr& bx = ineq2[0];
01182   DebugAssert(!isPlus(ax) && !isRational(ax),
01183         "TheoryArithNew::processFiniteInterval:\n ax = "+ax.toString());
01184   DebugAssert(!isPlus(bx) && !isRational(bx),
01185         "TheoryArithNew::processFiniteInterval:\n bx = "+bx.toString());
01186   Expr a = isMult(ax)? ax[0] : rat(1);
01187   Expr b = isMult(bx)? bx[0] : rat(1);
01188 
01189   Theorem thm1(alphaLEax), thm2(bxLEbeta);
01190   // Multiply the inequalities by 'b' and 'a', and canonize them, if necessary
01191   if(a != b) {
01192     thm1 = canonPred(iffMP(alphaLEax, d_rules->multIneqn(ineq1, b)));
01193     thm2 = canonPred(iffMP(bxLEbeta, d_rules->multIneqn(ineq2, a)));
01194   }
01195   // Check that a*beta - b*alpha == c > 0
01196   const Expr& alphaLEt = thm1.getExpr();
01197   const Expr& alpha = alphaLEt[0];
01198   const Expr& t = alphaLEt[1];
01199   const Expr& beta = thm2.getExpr()[1];
01200   Expr c = canon(beta - alpha).getRHS();
01201 
01202   if(c.isRational() && c.getRational() >= 1) {
01203     // This is a finite interval.  First, derive t <= alpha + c:
01204     // canon(alpha+c) => (alpha+c == beta) ==> [symmetry] beta == alpha+c
01205     // Then substitute that in thm2
01206     Theorem bEQac = symmetryRule(canon(alpha + c));
01207     // Substitute beta == alpha+c for the second child of thm2
01208     vector<unsigned> changed;
01209     vector<Theorem> thms;
01210     changed.push_back(1);
01211     thms.push_back(bEQac);
01212     Theorem tLEac = substitutivityRule(thm2.getExpr(), changed, thms);
01213     tLEac = iffMP(thm2, tLEac);
01214     // Derive and enqueue the finite interval constraint
01215     Theorem isInta(isIntegerThm(alpha));
01216     Theorem isIntt(isIntegerThm(t));
01217     enqueueFact(d_rules->finiteInterval(thm1, tLEac, isInta, isIntt));
01218   }
01219 }
01220 
01221 
01222 void
01223 TheoryArithNew::processFiniteIntervals(const Expr& x) {
01224   // If x is not integer, do not bother
01225   if(!isInteger(x)) return;
01226   // Process every pair of the opposing inequalities for 'x'
01227   ExprMap<CDList<Ineq> *>::iterator iLeft, iRight;
01228   iLeft = d_inequalitiesLeftDB.find(x);
01229   if(iLeft == d_inequalitiesLeftDB.end()) return;
01230   iRight = d_inequalitiesRightDB.find(x);
01231   if(iRight == d_inequalitiesRightDB.end()) return;
01232   // There are some opposing inequalities; get the lists
01233   CDList<Ineq>& ineqsLeft = *(iLeft->second);
01234   CDList<Ineq>& ineqsRight = *(iRight->second);
01235   // Get the sizes of the lists
01236   size_t sizeLeft = ineqsLeft.size();
01237   size_t sizeRight = ineqsRight.size();
01238   // Process all the pairs of the opposing inequalities
01239   for(size_t l=0; l<sizeLeft; ++l)
01240     for(size_t r=0; r<sizeRight; ++r)
01241       processFiniteInterval(ineqsRight[r], ineqsLeft[l]);
01242 }
01243 
01244 /*! This function recursively decends expression tree <strong>without
01245  *  caching</strong> until it hits a node that is already setup.  Be
01246  *  careful on what expressions you are calling it.
01247  */
01248 void
01249 TheoryArithNew::setupRec(const Expr& e) {
01250   if(e.hasFind()) return;
01251   // First, set up the kids recursively
01252   for (int k = 0; k < e.arity(); ++k) {
01253     setupRec(e[k]);
01254   }
01255   // Create a find pointer for e
01256   e.setFind(reflexivityRule(e));
01257   // And call our own setup()   
01258   setup(e);
01259 }
01260 
01261 
01262 /*!
01263  * Keep track of all finitely bounded integer variables in shared
01264  * terms.
01265  *
01266  * When a new shared term t is reported, all of its variables x1...xn
01267  * are added to the set d_sharedVars.  
01268  *
01269  * For each newly added variable x, check all of its opposing
01270  * inequalities, find out all the finite bounds and assert the
01271  * corresponding GRAY_SHADOW constraints.
01272  *
01273  * When projecting integer inequalities, the database d_sharedVars is
01274  * consulted, and if the variable is in it, check the shadows for
01275  * being a finite interval, and produce the additional GRAY_SHADOW
01276  * constraints.
01277  */
01278 void TheoryArithNew::addSharedTerm(const Expr& e) {
01279   d_sharedTerms[e] = true;
01280   /*
01281   set<Expr> cache; // Aux. cache for collecting i-leaves
01282   vector<Expr> vars; // Vector of vars in e
01283   collectVars(e, vars, cache);
01284   for(vector<Expr>::iterator i=vars.begin(), iend=vars.end(); i!=iend; ++i) {
01285     if(d_sharedVars.count(*i) == 0) {
01286       TRACE("arith shared", "TheoryArithNew::addSharedTerm: new var = ", *i, "");
01287       processFiniteIntervals(*i);
01288       d_sharedVars[*i]=true;
01289     }
01290   }
01291   */
01292 }
01293 
01294 void TheoryArithNew::refineCounterExample()
01295 {
01296   d_inModelCreation = true;
01297   TRACE("model", "refineCounterExample[TheoryArithNew] ", "", "{");
01298   CDMap<Expr, bool>::iterator it = d_sharedTerms.begin(), it2,
01299     iend = d_sharedTerms.end();
01300   // Add equalities over all pairs of shared terms as suggested
01301   // splitters.  Notice, that we want to split on equality
01302   // (positively) first, to reduce the size of the model.
01303   for(; it!=iend; ++it) {
01304     // Copy by value: the elements in the pair from *it are NOT refs in CDMap
01305     Expr e1 = (*it).first;
01306     for(it2 = it, ++it2; it2!=iend; ++it2) {
01307       Expr e2 = (*it2).first;
01308       DebugAssert(isReal(getBaseType(e1)),
01309       "TheoryArithNew::refineCounterExample: e1 = "+e1.toString()
01310       +"\n type(e1) = "+e1.getType().toString());
01311       if(findExpr(e1) != findExpr(e2)) {
01312   DebugAssert(isReal(getBaseType(e2)),
01313         "TheoryArithNew::refineCounterExample: e2 = "+e2.toString()
01314         +"\n type(e2) = "+e2.getType().toString());
01315   Expr eq = simplifyExpr(e1.eqExpr(e2));
01316         if (!eq.isBoolConst())
01317           addSplitter(eq);
01318       }
01319     }
01320   }
01321   TRACE("model", "refineCounterExample[Theory::Arith] ", "", "}");
01322 }
01323 
01324 
01325 void
01326 TheoryArithNew::findRationalBound(const Expr& varSide, const Expr& ratSide, 
01327              const Expr& var,
01328              Rational &r)
01329 {
01330   Expr c, x;
01331   separateMonomial(varSide, c, x);
01332   
01333   DebugAssert(findExpr(c).isRational(), 
01334         "seperateMonomial failed"); 
01335   DebugAssert(findExpr(ratSide).isRational(), 
01336         "smallest variable in graph, should not have variables"
01337         " in inequalities: ");
01338   DebugAssert(x == var, "separateMonomial found different variable: " 
01339         + var.toString());
01340   r = findExpr(ratSide).getRational() / findExpr(c).getRational();
01341 } 
01342 
01343 bool
01344 TheoryArithNew::findBounds(const Expr& e, Rational& lub, Rational&  glb)
01345 {
01346   bool strictLB=false, strictUB=false;
01347   bool right = (d_inequalitiesRightDB.count(e) > 0
01348            && d_inequalitiesRightDB[e]->size() > 0);
01349   bool left = (d_inequalitiesLeftDB.count(e) > 0
01350          && d_inequalitiesLeftDB[e]->size() > 0);
01351   int numRight = 0, numLeft = 0;
01352   if(right) { //rationals less than e
01353     CDList<Ineq> * ratsLTe = d_inequalitiesRightDB[e];
01354     for(unsigned int i=0; i<ratsLTe->size(); i++) {
01355       DebugAssert((*ratsLTe)[i].varOnRHS(), "variable on wrong side!");
01356       Expr ineq = (*ratsLTe)[i].ineq().getExpr();
01357       Expr leftSide = ineq[0], rightSide = ineq[1];
01358       Rational r;
01359       findRationalBound(rightSide, leftSide, e , r);
01360       if(numRight==0 || r>glb){
01361   glb = r;
01362   strictLB = isLT(ineq);
01363       }
01364       numRight++;
01365     }
01366     TRACE("model", "   =>Lower bound ", glb.toString(), "");
01367   }
01368   if(left) {   //rationals greater than e
01369     CDList<Ineq> * ratsGTe = d_inequalitiesLeftDB[e];
01370     for(unsigned int i=0; i<ratsGTe->size(); i++) { 
01371       DebugAssert((*ratsGTe)[i].varOnLHS(), "variable on wrong side!");
01372       Expr ineq = (*ratsGTe)[i].ineq().getExpr();
01373       Expr leftSide = ineq[0], rightSide = ineq[1];
01374       Rational r;
01375       findRationalBound(leftSide, rightSide, e, r); 
01376       if(numLeft==0 || r<lub) {
01377   lub = r;
01378   strictUB = isLT(ineq);
01379       }
01380       numLeft++;
01381     }
01382     TRACE("model", "   =>Upper bound ", lub.toString(), "");
01383   }
01384   if(!left && !right) {
01385       lub = 0; 
01386       glb = 0;
01387   }
01388   if(!left && right) {lub = glb +2;}
01389   if(!right && left)  {glb =  lub-2;}
01390   DebugAssert(glb <= lub, "Greatest lower bound needs to be smaller "
01391         "than least upper bound"); 
01392   return strictLB;
01393 }
01394 
01395 void TheoryArithNew::assignVariables(std::vector<Expr>&v)
01396 {
01397   int count = 0;
01398   while (v.size() > 0) {
01399     std::vector<Expr> bottom;
01400     d_graph.selectSmallest(v, bottom);
01401     TRACE("model", "Finding variables to assign. Iteration # ", count, "");
01402     for(unsigned int i = 0; i<bottom.size(); i++) {
01403       Expr e = bottom[i];
01404       TRACE("model", "Found: ", e, "");
01405       // Check if it is already a concrete constant
01406       if(e.isRational()) continue;
01407       
01408       Rational lub, glb;
01409       bool strictLB;
01410       strictLB = findBounds(e, lub, glb);
01411       Rational mid;
01412       if(isInteger(e)) {
01413   if(strictLB && glb.isInteger())
01414     mid = glb + 1;
01415   else
01416     mid = ceil(glb);
01417       }
01418       else
01419   mid = (lub + glb)/2;
01420       TRACE("model", "Assigning mid = ", mid, " {");
01421       assignValue(e, rat(mid));
01422       TRACE("model", "Assigned find(e) = ", findExpr(e), " }");
01423       if(inconsistent()) return; // Punt immediately if failed
01424     }
01425     count++;
01426   }
01427 }
01428 
01429 void TheoryArithNew::computeModelBasic(const std::vector<Expr>& v)
01430 {
01431   d_inModelCreation = true;
01432   vector<Expr> reps;
01433   TRACE("model", "Arith=>computeModel ", "", "{");
01434   for(unsigned int i=0; i <v.size(); ++i) {
01435     const Expr& e = v[i];
01436     if(findExpr(e) == e) {
01437       TRACE("model", "arith variable:", e , "");
01438       reps.push_back(e);
01439     }
01440     else {
01441       TRACE("model", "arith variable:", e , "");
01442       TRACE("model", " ==> is defined by: ", findExpr(e) , "");
01443     }
01444   }
01445   assignVariables(reps);
01446   TRACE("model", "Arith=>computeModel", "", "}");
01447   d_inModelCreation = false;
01448 }
01449 
01450 // For any arith expression 'e', if the subexpressions are assigned
01451 // concrete values, then find(e) must already be a concrete value.
01452 void TheoryArithNew::computeModel(const Expr& e, vector<Expr>& vars) {
01453   DebugAssert(findExpr(e).isRational(), "TheoryArithNew::computeModel("
01454         +e.toString()+")\n e is not assigned concrete value.\n"
01455         +" find(e) = "+findExpr(e).toString());
01456   assignValue(simplify(e));
01457   vars.push_back(e);
01458 }
01459 
01460 /*! accepts a rewrite theorem over eqn|ineqn and normalizes it
01461  *  and returns a theorem to that effect. assumes e is non-trivial
01462  *  i.e. e is not '0=const' or 'const=0' or '0 <= const' etc.
01463  */
01464 Theorem TheoryArithNew::normalize(const Expr& e, NormalizationType type) {
01465   
01466   //e is an eqn or ineqn. e is not a trivial eqn or ineqn
01467     //trivial means 0 = const or 0 <= const.
01468   
01469     // Trace the normalization on the arithm trace flag
01470     TRACE("arith", "normalize(", e, ") {"); 
01471   
01472     // The constraint must be an equality or inequality
01473     DebugAssert(isIneq(e), "normalize: input must be an inequality: " + e.toString());
01474   
01475     // The right side must be a sum or a multiplication 
01476     DebugAssert(isPlus(e[1]) || (isMult(e[1]) && e[1][0].isRational()), "normalize: right side must be a sum or a mult: " + e.toString());
01477   
01478     // The left side must be 0
01479     DebugAssert(e[0].isRational() && 0 == e[0].getRational(), "normalize: e[0] must be 0" + e.toString());
01480   
01481     // Compute the factor to use for normalizing the inequality  
01482     Expr factor;
01483     // If a mult, just take the coefficient
01484     if (isMult(e[1])) factor = rat(1/e[1][0].getRational());
01485     // Otherwise compute it
01486     else factor = computeNormalFactor(e[1], type);
01487   
01488     // Trace the rezult on the arith flag
01489     TRACE("arith", "normalize: factor = ", factor, "");
01490     
01491     // The theorem we will be creating (reflexivity, just in case)
01492     Theorem thm;
01493   
01494     // Now multiply the equality by the factor, unless it is 1
01495     if(factor.getRational() != 1)
01496       switch(e.getKind()) {
01497         
01498         case LE:
01499         case LT:
01500         case GE:
01501         case GT:
01502                           
01503             // Multiply inequality by the factor  
01504             thm = d_rules->multIneqn(e, factor);
01505             
01506             // And canonize the result
01507             thm = canonPredEquiv(thm);
01508       
01509             // Inequalities case break
01510             break;
01511         
01512         default:
01513             
01514             // How did we get here
01515             FatalAssert(false, "normalize: control should not reach here" + e.toString());
01516             
01517             // Default case break
01518             break;
01519       }
01520     else 
01521       // If rational is 1 then nothing to be done
01522       thm = reflexivityRule(e);
01523 
01524   
01525   // Trace the resulting theorem on the arith trace flag
01526   TRACE("arith", "normalize => ", thm, " }");
01527   
01528   // Return the explaining theorem
01529   return(thm);
01530 }
01531 
01532 
01533 Theorem TheoryArithNew::normalize(const Theorem& eIffEqn, NormalizationType type) {
01534   return transitivityRule(eIffEqn, normalize(eIffEqn.getRHS(), type));
01535 }
01536 
01537 Theorem TheoryArithNew::rafineIntegerConstraints(const Theorem& thm) {
01538   
01539   // The resulting theorem
01540   Theorem result = thm;
01541   
01542   // Get the constraint
01543   const Expr& constr = result.getRHS();
01544   
01545   // Get the proof that this constraint is an integer constraint 
01546   Theorem isIntConstraintThm = isIntegerThm(constr[1]);
01547   
01548   // If not an integer, just return the same theorem
01549   if (isIntConstraintThm.isNull()) return result;
01550   
01551   // Trace the integer
01552   TRACE("integer", "TheoryArithNew::rafineIntegerConstraints(", constr, ")");
01553   TRACE("integer", "TheoryArithNew::rafineIntegerConstraints: int proof:", isIntConstraintThm, "");
01554   
01555   // Get the left-hand of the constraint (the rational)
01556   Rational r = constr[0].getRational();
01557   
01558   // If it is a strict inequality, make it non-strict
01559   if (constr.getKind() == LT || constr.getKind() == GT || !r.isInteger())
01560     result = transitivityRule(result, d_rules->rafineStrictInteger(isIntConstraintThm, constr));
01561 
01562   // Return the result    
01563   return result;
01564 }
01565 
01566 Theorem TheoryArithNew::rewrite(const Expr& e) {
01567   
01568     // Leaves are assumeed to be already simplified
01569     //DebugAssert(leavesAreSimp(e), "Expected leaves to be simplified");
01570   
01571     // Trace the rewrites on the arith flag
01572     TRACE("arith", "TheoryArithNew::rewrite(", e, ") {");
01573   
01574     // The resulting theorem object
01575     Theorem thm;
01576     
01577     // Are we in the NOT expression
01578     bool isNot = false;
01579   
01580     // If the expression is not a term, i.e a literal  
01581     if (!e.isTerm()) {
01582        
01583       // If the expression is not a literal just return the reflexivity theorem
01584       if (!e.isAbsLiteral()) { 
01585         
01586         // Set the expression REWRITE NORMAL FLAG
01587         e.setRewriteNormal();
01588         
01589         // Create the reflaxivity rule
01590           thm = reflexivityRule(e);
01591         
01592           // Trace the non-literal rewrite 
01593           TRACE("arith", "TheoryArithNew::rewrite[non-literal] => ", thm, " }");
01594         
01595           // Return the theorem
01596           return thm;
01597       }
01598     
01599       // Its a literal, switch the arithmetic relations for rewrite
01600       switch(e.getKind()) {
01601     
01602         // Equality
01603         case EQ: {
01604         
01605           // Rewrite equality as two inequalities
01606           thm = d_rules->eqToIneq(e);
01607           Expr andExpr = thm.getRHS();
01608     
01609           // Rewrite the left inequality
01610           Expr leExpr = andExpr[0];
01611           const Theorem& thmLeft = simplify(leExpr);
01612           
01613           // Rewrite the right inequality 
01614           Expr geExpr = andExpr[1]; 
01615           const Theorem& thmRight = simplify(geExpr);         
01616         
01617           // Do the substitution
01618           thm = transitivityRule(thm, substitutivityRule(andExpr, thmLeft, thmRight));
01619               
01620           // EQ case break
01621           break;
01622         }
01623     
01624         // Negation of an atomic formula
01625         case NOT:    
01626 
01627         // If here, the equality should have already been rewritten as two disequalitites
01628         DebugAssert(e[0].getKind() != EQ, "TheoryArithNew::rewrite, not expecting equality under negation");
01629           
01630         // Negate the inequality and continue with the normal case
01631         thm = d_rules->negatedInequality(e);
01632       
01633         // IMPORTANT, LE, LT, GE, GT MUST BE UNDER THIS
01634         isNot = true;
01635       
01636         case LT:
01637         case GT:
01638         case LE:
01639         case GE:
01640         
01641         // Move everything to the right side
01642         if (isNot)
01643           thm = transitivityRule(thm, d_rules->rightMinusLeft(thm.getRHS())); 
01644         else
01645           thm = d_rules->rightMinusLeft(e);
01646 
01647             // Canonise children again
01648             thm = canonPredEquiv(thm);
01649    
01650             // If a trivial inequation just return true or false
01651             if ((thm.getRHS())[1].isRational()) 
01652           thm = transitivityRule(thm, d_rules->constPredicate(thm.getRHS()));
01653             else { // It's a non-trivial inequation
01654               
01655               // Normalize the inequality 
01656           thm = normalize(thm, NORMALIZE_UNIT);
01657         
01658           // Get the left and right side
01659           const Expr& normalised = thm.getRHS();
01660           const Expr& rightSide = normalised[1]; 
01661           const Expr& leftSide  = normalised[0];
01662         
01663           // If the right side is a sum, move the rational to the right side
01664           if (isPlus(rightSide)) {
01665             // Check if the fist of the sum is a rational
01666             if (rightSide[0].isRational()) {
01667               // Add the negative constant to both sides
01668               thm = transitivityRule(thm, d_rules->plusPredicate(leftSide, rightSide, rat(-rightSide[0].getRational()), thm.getRHS().getKind()));    
01669               // Canonize the left side
01670               const Theorem& thmLeft  = d_rules->canonPlus(thm.getRHS()[0]);
01671               // Canonize the right side
01672               const Theorem& thmRight = d_rules->canonPlus(thm.getRHS()[1]); 
01673               // Sunstitute the canonised into the expression
01674               thm = transitivityRule(thm, substitutivityRule(thm.getRHS(), thmLeft, thmRight));
01675             }
01676           }
01677             }
01678           
01679             // If a strict inequality on integers, or bounded by a non-integer, rafine it to a non-strict one
01680 //            thm = rafineIntegerConstraints(thm);
01681             
01682             // EQ, LE, LT, GE, GT case break        
01683             break;
01684     
01685         case IS_INTEGER:
01686         
01687           // TRACE
01688           TRACE("integer", "Got ", e.toString(), "");
01689           TRACE("integer", "Type ", e[0].getType().toString(), ""); 
01690         
01691           // TODO: What with the integers
01692         thm = d_rules->dummyTheorem(e);
01693           
01694           // IS_INTEGER case break
01695           break;
01696     
01697         default:
01698             
01699             // How did we get here
01700             FatalAssert(false, "Theory_Arith::rewrite: control should not reach here");
01701       
01702             // default break
01703             break;
01704             
01705       } // End relation case
01706       
01707     } else { // Expression is a term 
01708       
01709       // Terms that don't contain formulas are canonised via the canon() function 
01710       if (e.isAtomic()) thm = canon(e);
01711       // Otherwise just return the same expression
01712       else thm = reflexivityRule(e);
01713     }
01714   
01715     // Compact the theorem into one rule
01716     //thm = d_rules->trustedRewrite(e, thm.getRHS());
01717   
01718     // Arith canonization is idempotent, the theory should stay the same
01719     if (theoryOf(thm.getRHS()) == this)
01720       thm.getRHS().setRewriteNormal();
01721     
01722     // Finished, trace the end of rewrite on the arith flag 
01723     TRACE("arith", "TheoryArithNew::rewrite => ", thm, " }");
01724   
01725     // Return the final theorem
01726     return thm;
01727 }
01728 
01729 
01730 void TheoryArithNew::setup(const Expr& e)
01731 {
01732   //If the expression is not a term but rather a predicate
01733     if (!e.isTerm()) {
01734       
01735       // If not or equality, just return
01736       if (e.isNot() || e.isEq()) return;
01737     
01738       // TODO: Integer::eqToIneq
01739       if (isIntPred(e)) return;
01740     
01741       // Check if the expression is canonised as (SUM t1 ... tn) @ b
01742       DebugAssert(isIneq(e) && e[0].isRational(), "TheoryArithNew::setup: b @ (sum t1 ... tn) :" + e.toString());
01743     
01744       // Add the sum to the notify list of e
01745       e[1].addToNotify(this, e);
01746       
01747     } else {
01748       
01749       // The arity of the expression 
01750       int arity(e.arity());
01751       
01752       // Go through all the children and add this expression to their notify lists 
01753       for (int k = 0 ; k < arity; k ++) {
01754         
01755         // Add the to the notify list of the children 
01756         e[k].addToNotify(this, e);
01757         
01758         // Trace this notification add
01759         TRACE("arith setup", "e["+int2string(k)+"]: ", *(e[k].getNotify()), "");
01760       }
01761     }
01762 }
01763 
01764 
01765 void TheoryArithNew::update(const Theorem& e, const Expr& d)
01766 {
01767   if (inconsistent()) return;
01768   IF_DEBUG(debugger.counter("arith update total")++);
01769   if (!d.hasFind()) return;
01770   if (isIneq(d)) {
01771     // Substitute e[1] for e[0] in d and enqueue new inequality
01772     DebugAssert(e.getLHS() == d[1], "Mismatch");
01773     Theorem thm = find(d);
01774     //    DebugAssert(thm.getRHS() == trueExpr(), "Expected find = true");
01775     vector<unsigned> changed;
01776     vector<Theorem> children;
01777     changed.push_back(1);
01778     children.push_back(e);
01779     Theorem thm2 = substitutivityRule(d, changed, children);
01780     if (thm.getRHS() == trueExpr()) {
01781       enqueueFact(iffMP(getCommonRules()->iffTrueElim(thm), thm2));
01782     }
01783     else {
01784       enqueueFact(getCommonRules()->iffFalseElim(
01785         transitivityRule(symmetryRule(thm2), thm)));
01786     }
01787     IF_DEBUG(debugger.counter("arith update ineq")++);
01788   }
01789   else if (find(d).getRHS() == d) {
01790     // Substitute e[1] for e[0] in d and assert the result equal to d
01791     Theorem thm = canonSimp(d);
01792     TRACE("arith", "TheoryArithNew::update(): thm = ", thm, "");
01793     DebugAssert(leavesAreSimp(thm.getRHS()), "updateHelper error: "
01794     +thm.getExpr().toString());
01795     assertEqualities(transitivityRule(thm, rewrite(thm.getRHS())));
01796     IF_DEBUG(debugger.counter("arith update find(d)=d")++);
01797   }
01798 }
01799 
01800 
01801 Theorem TheoryArithNew::solve(const Theorem& thm)
01802 {
01803   DebugAssert(thm.isRewrite() && thm.getLHS().isTerm(), "");
01804   const Expr& lhs = thm.getLHS();
01805   const Expr& rhs = thm.getRHS();
01806 
01807   // Check for already solved equalities.
01808 
01809   // Have to be careful about the types: integer variable cannot be
01810   // assigned a real term.  Also, watch for e[0] being a subexpression
01811   // of e[1]: this would create an unsimplifiable expression.
01812   if (isLeaf(lhs) && !isLeafIn(lhs, rhs) 
01813       && (lhs.getType() != intType() || isInteger(rhs))
01814       // && !e[0].subExprOf(e[1])
01815       )
01816     return thm;
01817 
01818   // Symmetric version is already solved
01819   if (isLeaf(rhs) && !isLeafIn(rhs, lhs)
01820       && (rhs.getType() != intType() || isInteger(lhs))
01821       // && !e[1].subExprOf(e[0])
01822       )
01823     return symmetryRule(thm);
01824 
01825   return doSolve(thm);
01826 }
01827 
01828  
01829 void TheoryArithNew::computeModelTerm(const Expr& e, std::vector<Expr>& v) {
01830   switch(e.getKind()) {
01831   case RATIONAL_EXPR: // Skip the constants
01832     break;
01833   case PLUS:
01834   case MULT:
01835   case DIVIDE:
01836   case POW: // This is not a variable; extract the variables from children
01837     for(Expr::iterator i=e.begin(), iend=e.end(); i!=iend; ++i)
01838       // getModelTerm(*i, v);
01839       v.push_back(*i);
01840     break;
01841   default: { // Otherwise it's a variable.  Check if it has a find pointer
01842     Expr e2(findExpr(e));
01843     if(e==e2) {
01844       TRACE("model", "TheoryArithNew::computeModelTerm(", e, "): a variable");
01845       // Leave it alone (it has no descendants)
01846       // v.push_back(e);
01847     } else {
01848       TRACE("model", "TheoryArithNew::computeModelTerm("+e.toString()
01849       +"): has find pointer to ", e2, "");
01850       v.push_back(e2);
01851     }
01852   }
01853   }
01854 }
01855 
01856 Expr TheoryArithNew::computeTypePred(const Type& t, const Expr& e) {
01857   Expr tExpr = t.getExpr();
01858   switch(tExpr.getKind()) {
01859   case INT:  
01860     return Expr(IS_INTEGER, e);
01861   case SUBRANGE: {
01862     std::vector<Expr> kids;
01863     kids.push_back(Expr(IS_INTEGER, e));
01864     kids.push_back(leExpr(tExpr[0], e));
01865     kids.push_back(leExpr(e, tExpr[1]));
01866     return andExpr(kids);
01867   }
01868   default:
01869     return e.getEM()->trueExpr();
01870   }
01871 }
01872 
01873 
01874 void TheoryArithNew::checkAssertEqInvariant(const Theorem& e)
01875 {
01876   return;
01877 }
01878 
01879 
01880 void TheoryArithNew::checkType(const Expr& e)
01881 {
01882   switch (e.getKind()) {
01883     case INT:
01884     case REAL:
01885       if (e.arity() > 0) {
01886         throw Exception("Ill-formed arithmetic type: "+e.toString());
01887       }
01888       break;
01889     case SUBRANGE:
01890       if (e.arity() != 2 ||
01891           !isIntegerConst(e[0]) ||
01892           !isIntegerConst(e[1]) ||
01893           e[0].getRational() > e[1].getRational()) {
01894         throw Exception("bad SUBRANGE type expression"+e.toString());
01895       }
01896       break;
01897     default:
01898       DebugAssert(false, "Unexpected kind in TheoryArithNew::checkType"
01899                   +getEM()->getKindName(e.getKind()));
01900   }
01901 }
01902 
01903 void TheoryArithNew::computeType(const Expr& e)
01904 {
01905   switch (e.getKind()) {
01906     case REAL_CONST:
01907       e.setType(d_realType);
01908       break;
01909     case RATIONAL_EXPR:
01910       if(e.getRational().isInteger())
01911         e.setType(d_intType);
01912       else
01913         e.setType(d_realType);
01914       break;
01915     case UMINUS:
01916     case PLUS:
01917     case MINUS:
01918     case MULT:
01919     case POW: {
01920       bool isInt = true;
01921       for(int k = 0; k < e.arity(); ++k) {
01922         if(d_realType != getBaseType(e[k]))
01923           throw TypecheckException("Expecting type REAL with `" +
01924                                    getEM()->getKindName(e.getKind()) + "',\n"+
01925                                    "but got a " + getBaseType(e[k]).toString()+ 
01926                                    " for:\n" + e.toString());
01927         if(isInt && !isInteger(e[k]))
01928           isInt = false;
01929       }
01930       if(isInt)
01931         e.setType(d_intType);
01932       else
01933         e.setType(d_realType);
01934       break;
01935     }
01936     case DIVIDE: {
01937       Expr numerator = e[0];
01938       Expr denominator = e[1];
01939       if (getBaseType(numerator) != d_realType ||
01940           getBaseType(denominator) != d_realType) {
01941         throw TypecheckException("Expecting only REAL types with `DIVIDE',\n"
01942                                  "but got " + getBaseType(numerator).toString()+ 
01943                                  " and " + getBaseType(denominator).toString() +
01944                                  " for:\n" + e.toString());
01945       }
01946       if(denominator.isRational() && 1 == denominator.getRational())
01947         e.setType(numerator.getType());
01948       else
01949         e.setType(d_realType);
01950       break;
01951     }
01952     case LT:
01953     case LE:
01954     case GT:
01955     case GE:
01956     case GRAY_SHADOW:
01957       // Need to know types for all exprs -Clark
01958       //    e.setType(boolType());
01959       //    break;
01960     case DARK_SHADOW:
01961       for (int k = 0; k < e.arity(); ++k) {
01962         if (d_realType != getBaseType(e[k]))
01963           throw TypecheckException("Expecting type REAL with `" +
01964                                    getEM()->getKindName(e.getKind()) + "',\n"+
01965                                    "but got a " + getBaseType(e[k]).toString()+ 
01966                                    " for:\n" + e.toString());
01967       }
01968       
01969       e.setType(boolType());
01970       break;
01971     case IS_INTEGER:
01972       if(d_realType != getBaseType(e[0]))
01973         throw TypecheckException("Expected type REAL, but got "
01974                                  +getBaseType(e[0]).toString()
01975                                  +"\n\nExpr = "+e.toString());
01976       e.setType(boolType());
01977       break;
01978     default:
01979       DebugAssert(false,"TheoryArithNew::computeType: unexpected expression:\n "
01980                   +e.toString());
01981       break;
01982   }
01983 }
01984 
01985 Type TheoryArithNew::computeBaseType(const Type& t) {
01986   IF_DEBUG(int kind = t.getExpr().getKind());
01987   DebugAssert(kind==INT || kind==REAL || kind==SUBRANGE,
01988         "TheoryArithNew::computeBaseType("+t.toString()+")");
01989   return realType();
01990 }
01991 
01992 Expr TheoryArithNew::computeTCC(const Expr& e) {
01993   Expr tcc(Theory::computeTCC(e));
01994   switch(e.getKind()) {
01995   case DIVIDE:
01996     DebugAssert(e.arity() == 2, "");
01997     return tcc.andExpr(!(e[1].eqExpr(rat(0))));
01998   default:
01999     return tcc;
02000   }
02001 }
02002 
02003 ///////////////////////////////////////////////////////////////////////////////
02004 //parseExprOp:
02005 //translating special Exprs to regular EXPR??
02006 ///////////////////////////////////////////////////////////////////////////////
02007 Expr TheoryArithNew::parseExprOp(const Expr& e) {
02008   TRACE("parser", "TheoryArithNew::parseExprOp(", e, ")");
02009   //std::cout << "Were here";
02010   // If the expression is not a list, it must have been already
02011   // parsed, so just return it as is.
02012   switch(e.getKind()) {
02013   case ID: {
02014     int kind = getEM()->getKind(e[0].getString());
02015     switch(kind) {
02016     case NULL_KIND: return e; // nothing to do
02017     case REAL:
02018     case INT:
02019     case NEGINF: 
02020     case POSINF: return getEM()->newLeafExpr(kind);
02021     default:
02022       DebugAssert(false, "Bad use of bare keyword: "+e.toString());
02023       return e;
02024     }
02025   }
02026   case RAW_LIST: break; // break out of switch, do the hard work
02027   default:
02028     return e;
02029   }
02030 
02031   DebugAssert(e.getKind() == RAW_LIST && e.arity() > 0,
02032         "TheoryArithNew::parseExprOp:\n e = "+e.toString());
02033   
02034   const Expr& c1 = e[0][0];
02035   int kind = getEM()->getKind(c1.getString());
02036   switch(kind) {
02037     case UMINUS: {
02038       if(e.arity() != 2)
02039   throw ParserException("UMINUS requires exactly one argument: "
02040       +e.toString());
02041       return uminusExpr(parseExpr(e[1])); 
02042     }
02043     case PLUS: {
02044       vector<Expr> k;
02045       Expr::iterator i = e.begin(), iend=e.end();
02046       // Skip first element of the vector of kids in 'e'.
02047       // The first element is the operator.
02048       ++i; 
02049       // Parse the kids of e and push them into the vector 'k'
02050       for(; i!=iend; ++i) k.push_back(parseExpr(*i));
02051       return plusExpr(k); 
02052     }
02053     case MINUS: {
02054       if(e.arity() == 2)
02055   return uminusExpr(parseExpr(e[1]));
02056       else if(e.arity() == 3)
02057   return minusExpr(parseExpr(e[1]), parseExpr(e[2]));
02058       else
02059   throw ParserException("MINUS requires one or two arguments:"
02060       +e.toString());
02061     }
02062     case MULT: {
02063       vector<Expr> k;
02064       Expr::iterator i = e.begin(), iend=e.end();
02065       // Skip first element of the vector of kids in 'e'.
02066       // The first element is the operator.
02067       ++i; 
02068       // Parse the kids of e and push them into the vector 'k'
02069       for(; i!=iend; ++i) k.push_back(parseExpr(*i));
02070       return multExpr(k);
02071     }
02072     case POW: { 
02073       return powExpr(parseExpr(e[1]), parseExpr(e[2])); 
02074     }
02075     case DIVIDE:
02076       { return divideExpr(parseExpr(e[1]), parseExpr(e[2]));  }
02077     case LT:  
02078       { return ltExpr(parseExpr(e[1]), parseExpr(e[2]));  }
02079     case LE:  
02080       { return leExpr(parseExpr(e[1]), parseExpr(e[2]));  }
02081     case GT:  
02082       { return gtExpr(parseExpr(e[1]), parseExpr(e[2]));  }
02083     case GE:  
02084       { return geExpr(parseExpr(e[1]), parseExpr(e[2]));  }
02085     case INTDIV:
02086     case MOD:
02087     case SUBRANGE: {
02088       vector<Expr> k;
02089       Expr::iterator i = e.begin(), iend=e.end();
02090       // Skip first element of the vector of kids in 'e'.
02091       // The first element is the operator.
02092       ++i; 
02093       // Parse the kids of e and push them into the vector 'k'
02094       for(; i!=iend; ++i) 
02095         k.push_back(parseExpr(*i));
02096       return Expr(kind, k, e.getEM());
02097     }
02098     case IS_INTEGER: {
02099       if(e.arity() != 2)
02100   throw ParserException("IS_INTEGER requires exactly one argument: "
02101       +e.toString());
02102       return Expr(IS_INTEGER, parseExpr(e[1]));
02103     }
02104     default:
02105       DebugAssert(false,
02106         "TheoryArithNew::parseExprOp: invalid input " + e.toString());
02107       break;
02108   }
02109   return e;
02110 }
02111 
02112 ///////////////////////////////////////////////////////////////////////////////
02113 // Pretty-printing                                                           //
02114 ///////////////////////////////////////////////////////////////////////////////
02115 
02116 
02117 ExprStream& TheoryArithNew::print(ExprStream& os, const Expr& e) {
02118   switch(os.lang()) {
02119     case SIMPLIFY_LANG:
02120       switch(e.getKind()) {
02121       case RATIONAL_EXPR:
02122   e.print(os);
02123   break;
02124       case SUBRANGE:
02125   os <<"ERROR:SUBRANGE:not supported in Simplify\n";
02126   break;
02127       case IS_INTEGER:
02128   os <<"ERROR:IS_INTEGER:not supported in Simplify\n";
02129   break;
02130       case PLUS:  {
02131   int i=0, iend=e.arity();
02132   os << "(+ "; 
02133   if(i!=iend) os << e[i];
02134   ++i;
02135   for(; i!=iend; ++i) os << " " << e[i];
02136   os <<  ")";
02137   break;
02138       }
02139       case MINUS:
02140   os << "(- " << e[0] << " " << e[1]<< ")";
02141   break;
02142       case UMINUS:
02143   os << "-" << e[0] ;
02144   break;
02145       case MULT:  {
02146   int i=0, iend=e.arity();
02147   os << "(* " ;
02148   if(i!=iend) os << e[i];
02149   ++i;
02150   for(; i!=iend; ++i) os << " " << e[i];
02151   os <<  ")";
02152   break;
02153       }
02154       case POW:
02155           os << "(" << push << e[1] << space << "^ " << e[0] << push << ")";
02156           break;
02157       case DIVIDE:
02158   os << "(" << push << e[0] << space << "/ " << e[1] << push << ")";
02159   break;
02160       case LT:
02161   if (isInt(e[0].getType()) || isInt(e[1].getType())) {
02162   }
02163   os << "(< " << e[0] << " " << e[1] <<")";
02164   break;
02165       case LE:
02166           os << "(<= " << e[0]  << " " << e[1] << ")";
02167           break;
02168       case GT:
02169   os << "(> " << e[0] << " " << e[1] << ")";
02170   break;
02171       case GE:
02172   os << "(>= " << e[0] << " " << e[1]  << ")";
02173   break;
02174       case DARK_SHADOW:
02175       case GRAY_SHADOW:
02176   os <<"ERROR:SHADOW:not supported in Simplify\n";
02177   break;
02178       default:
02179   // Print the top node in the default LISP format, continue with
02180   // pretty-printing for children.
02181           e.print(os);
02182     
02183           break;
02184       } 
02185       break; // end of case SIMPLIFY_LANG
02186       
02187     case PRESENTATION_LANG:
02188       switch(e.getKind()) {
02189         case REAL:
02190         case INT:
02191         case RATIONAL_EXPR:
02192         case NEGINF:
02193         case POSINF:
02194           e.print(os);
02195           break;
02196         case SUBRANGE:
02197           if(e.arity() != 2) e.printAST(os);
02198           else 
02199             os << "[" << push << e[0] << ".." << e[1] << push << "]";
02200           break;
02201         case IS_INTEGER:
02202     if(e.arity() == 1)
02203       os << "IS_INTEGER(" << push << e[0] << push << ")";
02204     else
02205       e.printAST(os);
02206     break;
02207         case PLUS:  {
02208           int i=0, iend=e.arity();
02209           os << "(" << push;
02210           if(i!=iend) os << e[i];
02211           ++i;
02212           for(; i!=iend; ++i) os << space << "+ " << e[i];
02213           os << push << ")";
02214           break;
02215         }
02216         case MINUS:
02217           os << "(" << push << e[0] << space << "- " << e[1] << push << ")";
02218           break;
02219         case UMINUS:
02220           os << "-(" << push << e[0] << push << ")";
02221           break;
02222         case MULT:  {
02223           int i=0, iend=e.arity();
02224           os << "(" << push;
02225           if(i!=iend) os << e[i];
02226           ++i;
02227           for(; i!=iend; ++i) os << space << "* " << e[i];
02228           os << push << ")";
02229           break;
02230         }
02231         case POW:
02232           os << "(" << push << e[1] << space << "^ " << e[0] << push << ")";
02233           break;
02234         case DIVIDE:
02235           os << "(" << push << e[0] << space << "/ " << e[1] << push << ")";
02236           break;
02237         case LT:
02238           if (isInt(e[0].getType()) || isInt(e[1].getType())) {
02239           }
02240           os << "(" << push << e[0] << space << "< " << e[1] << push << ")";
02241           break;
02242         case LE:
02243           os << "(" << push << e[0] << space << "<= " << e[1] << push << ")";
02244           break;
02245         case GT:
02246           os << "(" << push << e[0] << space << "> " << e[1] << push << ")";
02247           break;
02248         case GE:
02249           os << "(" << push << e[0] << space << ">= " << e[1] << push << ")";
02250           break;
02251         case DARK_SHADOW:
02252     os << "DARK_SHADOW(" << push << e[0] << ", " << space << e[1] << push << ")";
02253     break;
02254         case GRAY_SHADOW:
02255     os << "GRAY_SHADOW(" << push << e[0] << ","  << space << e[1]
02256        << "," << space << e[2] << "," << space << e[3] << push << ")";
02257     break;
02258         default:
02259           // Print the top node in the default LISP format, continue with
02260           // pretty-printing for children.
02261           e.printAST(os);
02262 
02263           break;
02264       } 
02265       break; // end of case PRESENTATION_LANG
02266     case SMTLIB_LANG: {
02267       switch(e.getKind()) {
02268         case REAL_CONST:
02269           printRational(os, e[0].getRational(), true);
02270           break;
02271         case RATIONAL_EXPR:
02272           printRational(os, e.getRational());
02273           break;
02274         case REAL:
02275           os << "Real";
02276           break;
02277         case INT:
02278           os << "Int";
02279           break;
02280         case SUBRANGE:
02281           throw SmtlibException("TheoryArithNew::print: SMTLIB: SUBRANGE not implemented");
02282 //           if(e.arity() != 2) e.print(os);
02283 //           else 
02284 //             os << "(" << push << "SUBRANGE" << space << e[0]
02285 //         << space << e[1] << push << ")";
02286           break;
02287         case IS_INTEGER:
02288     if(e.arity() == 1)
02289       os << "(" << push << "IsInt" << space << e[0] << push << ")";
02290     else
02291             throw SmtlibException("TheoryArithNew::print: SMTLIB: IS_INTEGER used unexpectedly");
02292     break;
02293         case PLUS:  {
02294           os << "(" << push << "+";
02295           Expr::iterator i = e.begin(), iend = e.end();
02296           for(; i!=iend; ++i) {
02297             os << space << (*i);
02298           }
02299           os << push << ")";
02300           break;
02301         }
02302         case MINUS: {
02303           os << "(" << push << "- " << e[0] << space << e[1] << push << ")";
02304           break;
02305         }
02306         case UMINUS: {
02307           os << "(" << push << "-" << space << e[0] << push << ")";
02308           break;
02309         }
02310         case MULT:  {
02311           int i=0, iend=e.arity();
02312           for(; i!=iend; ++i) {
02313             if (i < iend-1) {
02314               os << "(" << push << "*";
02315             }
02316             os << space << e[i];
02317           }
02318           for (i=0; i < iend-1; ++i) os << push << ")";
02319           break;
02320         }
02321         case POW:
02322           throw SmtlibException("TheoryArithNew::print: SMTLIB: POW not supported");
02323           //          os << "(" << push << "^ " << e[1] << space << e[0] << push << ")";
02324           break;
02325         case DIVIDE: {
02326           throw SmtlibException("TheoryArithNew::print: SMTLIB: unexpected use of DIVIDE");
02327           break;
02328         }
02329         case LT: {
02330           Rational r;
02331           os << "(" << push << "<" << space;
02332           os << e[0] << space << e[1] << push << ")";
02333           break;
02334         }
02335         case LE: {
02336           Rational r;
02337           os << "(" << push << "<=" << space;
02338           os << e[0] << space << e[1] << push << ")";
02339           break;
02340         }
02341         case GT: {
02342           Rational r;
02343           os << "(" << push << ">" << space;
02344           os << e[0] << space << e[1] << push << ")";
02345           break;
02346         }
02347         case GE: {
02348           Rational r;
02349           os << "(" << push << ">=" << space;
02350           os << e[0] << space << e[1] << push << ")";
02351           break;
02352         }
02353         case DARK_SHADOW:
02354           throw SmtlibException("TheoryArithNew::print: SMTLIB: DARK_SHADOW not supported");
02355     os << "(" << push << "DARK_SHADOW" << space << e[0]
02356        << space << e[1] << push << ")";
02357     break;
02358         case GRAY_SHADOW:
02359           throw SmtlibException("TheoryArithNew::print: SMTLIB: GRAY_SHADOW not supported");
02360     os << "GRAY_SHADOW(" << push << e[0] << ","  << space << e[1]
02361        << "," << space << e[2] << "," << space << e[3] << push << ")";
02362     break;
02363         default:
02364           throw SmtlibException("TheoryArithNew::print: SMTLIB: default not supported");
02365           // Print the top node in the default LISP format, continue with
02366           // pretty-printing for children.
02367           e.printAST(os);
02368 
02369           break;
02370       } 
02371       break; // end of case SMTLIB_LANG
02372     }
02373     case LISP_LANG:
02374       switch(e.getKind()) {
02375         case REAL:
02376         case INT:
02377         case RATIONAL_EXPR:
02378         case NEGINF:
02379         case POSINF:
02380           e.print(os);
02381           break;
02382         case SUBRANGE:
02383           if(e.arity() != 2) e.printAST(os);
02384           else 
02385             os << "(" << push << "SUBRANGE" << space << e[0]
02386          << space << e[1] << push << ")";
02387           break;
02388         case IS_INTEGER:
02389     if(e.arity() == 1)
02390       os << "(" << push << "IS_INTEGER" << space << e[0] << push << ")";
02391     else
02392       e.printAST(os);
02393     break;
02394         case PLUS:  {
02395           int i=0, iend=e.arity();
02396           os << "(" << push << "+";
02397           for(; i!=iend; ++i) os << space << e[i];
02398           os << push << ")";
02399           break;
02400         }
02401         case MINUS:
02402         //os << "(" << push << e[0] << space << "- " << e[1] << push << ")";
02403     os << "(" << push << "- " << e[0] << space << e[1] << push << ")";
02404           break;
02405         case UMINUS:
02406           os << "(" << push << "-" << space << e[0] << push << ")";
02407           break;
02408         case MULT:  {
02409           int i=0, iend=e.arity();
02410           os << "(" << push << "*";
02411           for(; i!=iend; ++i) os << space << e[i];
02412           os << push << ")";
02413           break;
02414         }
02415         case POW:
02416           os << "(" << push << "^ " << e[1] << space << e[0] << push << ")";
02417           break;
02418         case DIVIDE:
02419           os << "(" << push << "/ " << e[0] << space << e[1] << push << ")";
02420           break;
02421         case LT:
02422           os << "(" << push << "< " << e[0] << space << e[1] << push << ")";
02423           break;
02424         case LE:
02425           os << "(" << push << "<= " << e[0] << space << e[1] << push << ")";
02426           break;
02427         case GT:
02428           os << "(" << push << "> " << e[1]  << space << e[0] << push << ")";
02429           break;
02430         case GE:
02431           os << "(" << push << ">= " << e[0] << space << e[1] << push << ")";
02432           break;
02433         case DARK_SHADOW:
02434     os << "(" << push << "DARK_SHADOW" << space << e[0]
02435        << space << e[1] << push << ")";
02436     break;
02437         case GRAY_SHADOW:
02438     os << "(" << push << "GRAY_SHADOW" << space << e[0] << space
02439        << e[1] << space << e[2] << space << e[3] << push << ")";
02440     break;
02441         default:
02442           // Print the top node in the default LISP format, continue with
02443           // pretty-printing for children.
02444           e.printAST(os);
02445 
02446           break;
02447       } 
02448       break; // end of case LISP_LANG
02449     default:
02450      // Print the top node in the default LISP format, continue with
02451      // pretty-printing for children.
02452        e.printAST(os);
02453   }
02454   return os;
02455 }
02456 
02457 ///////////////////////////////////////////////////////////////////////////////
02458 // SIMPLEX SOLVER                                                            //
02459 ///////////////////////////////////////////////////////////////////////////////
02460 
02461 bool TheoryArithNew::isBasic(const Expr& x) const {
02462   // If it is on the right side of the tableaux the it is basic, non-basic otherwise
02463   return (tableaux.find(x) != tableaux.end());
02464 }
02465 
02466 void TheoryArithNew::pivot(const Expr& x_r, const Expr& x_s) {
02467   
02468   // Check that the variable is non-basic
02469   DebugAssert(isBasic(x_r), "TheoryArithNew::pivot, given variable should be basic" + x_r.toString());
02470   DebugAssert(!isBasic(x_s), "TheoryArithNew::pivot, given variable should be non-basic" + x_s.toString());
02471 
02472   // If trace is on for the simplex, print it out
02473   TRACE("simplex", "TheoryArithNew::pivot(", x_r.toString() + ", " + x_s.toString(), ")");
02474 
02475   // Get the iterator that points to the expression of x_r
02476   TebleauxMap::iterator it = tableaux.find(x_r); 
02477   
02478   // Get the expresiion of x_r
02479   Theorem x_r_Theorem = (*it).second;  
02480   
02481   // Erase the x_r expression from the tableaux
02482   tableaux.erase(x_r); // TODO: Add erase by iterator to ExprHashMap
02483   
02484   // Update the dependencies
02485   updateDependenciesRemove(x_r, x_r_Theorem.getExpr()[1]);
02486   
02487   // Construct the expresison (theorem) of x_s
02488   const Theorem& x_s_Theorem = pivotRule(x_r_Theorem, x_s); 
02489     
02490   // Replace all occurances of x_s with x_s_Expression
02491   substAndCanonizeTableaux(x_s_Theorem);
02492   
02493   // Update the dependencies
02494   updateDependenciesAdd(x_s, x_s_Theorem.getExpr()[1]);
02495   
02496   // Add the expression of x_s to the map
02497   tableaux[x_s] = x_s_Theorem;
02498 }
02499 
02500 void TheoryArithNew::update(const Expr& x_i, const EpsRational& v) {
02501   
02502   // Check that the variable is non-basic
02503   DebugAssert(tableaux.find(x_i) == tableaux.end(), "TheoryArithNew::update, given variable should be non-basic" + x_i.toString());
02504 
02505   // If trace is on for the simplex, print it out
02506   TRACE("simplex", "TheoryArithNew::update(", x_i.toString() + ", " + v.toString(), ")");
02507 
02508   // Remember the value of the x_i variable
02509   EpsRational old_value = getBeta(x_i);
02510 
02511   // If there are dependent variables then do the work 
02512   DependenciesMap::iterator find = dependenciesMap.find(x_i);
02513   if (find != dependenciesMap.end()) {
02514   
02515     // Go through all the variables that depend on x_i
02516     const set<Expr>& dependent = (*find).second;
02517     set<Expr>::const_iterator it     = dependent.begin();
02518     set<Expr>::const_iterator it_end = dependent.end();     
02519     // Fix the values of all the variables
02520     while (it != it_end) {
02521       
02522       // Get the current variable
02523       const Expr& x_j = *it;
02524       
02525       // Get the tableaux coefficient of the of x_i in the row of x_j (the right side ofg the tableaux expression)
02526       const Rational& a_ji = getTableauxEntry(x_j, x_i);
02527       
02528       // Update the beta valuation
02529       EpsRational b(getBeta(x_j)); // TODO: not necessary, beta's are all setup now
02530       EpsRational x_j_new = beta[x_j] = b + (v - old_value) * a_ji;
02531       
02532       // Check if the variable is sat or unsat under the new assignment 
02533       if (x_j_new < getLowerBound(x_j) || getUpperBound(x_j) < x_j_new)
02534         unsatBasicVariables.insert(x_j);
02535       else
02536         unsatBasicVariables.erase(x_j);
02537       
02538       // Go to the next one
02539       it ++;
02540     }
02541   }
02542   
02543   // Set the new value to x_i (x_i) is non_basic, no need to add to unsat set)
02544   beta[x_i] = v;
02545 }
02546       
02547 void TheoryArithNew::pivotAndUpdate(const Expr& x_i, const Expr& x_j, const EpsRational& v) {
02548   
02549   // Variable x_i should be a basic variable (left side of the tableaux) and x_j should be non_basic
02550   DebugAssert(tableaux.find(x_i) != tableaux.end(), "TheoryArithNew::pivotAndUpdate, " + x_i.toString() + " should be a basic variable");
02551   DebugAssert(tableaux.find(x_j) == tableaux.end(), "TheoryArithNew::pivotAndUpdate, " + x_j.toString() + " should be a non-basic variable");
02552 
02553   // If trace is on for the simplex, print it out
02554   TRACE("simplex", "TheoryArithNew::pivotAndUpdate(", x_i.toString() + ", " + x_j.toString() + ", " + v.toString(), ")");
02555 
02556   // The a_ij coefficient of the tableaux
02557   const Rational& a_ij = getTableauxEntry(x_i, x_j);
02558   
02559   // Let theta be the adjustment coefficient
02560   EpsRational theta((v - getBeta(x_i))/a_ij);
02561   
02562   // Set the new value to x_i (x_i becomes non-basic, hance sat)
02563   beta[x_i] = v;
02564   // x_i becomes non-basic, no need ==> it is sat
02565   unsatBasicVariables.erase(x_i);
02566   
02567   // Set the new value to x_j 
02568   EpsRational b = getBeta(x_j);
02569   EpsRational new_x_j = beta[x_j] = b + theta;
02570   
02571   // Check if the variable is sat or unsat under the new assignment 
02572   if (new_x_j < getLowerBound(x_j) || getUpperBound(x_j) < new_x_j)
02573     unsatBasicVariables.insert(x_j);
02574   else
02575     unsatBasicVariables.erase(x_j);
02576   
02577   // Go through all the variables that depend on x_j, and update their value (there will be at least one, i.e. x_i) // TODO: maybe optimise
02578   const set<Expr>& dependent = (*dependenciesMap.find(x_j)).second;
02579   set<Expr>::const_iterator it     = dependent.begin();
02580   set<Expr>::const_iterator it_end = dependent.end();     
02581   // Go throught all the basic variables
02582   while (it != it_end) {
02583     
02584     // Get the current variable
02585     const Expr& x_k = *it; 
02586     
02587     // If the basic variable is different from x_i update its value
02588     if (x_k != x_i) {
02589     
02590       // Get the a_kj coefficient from the tableaux
02591       const Rational& a_kj = getTableauxEntry(x_k, x_j);
02592       
02593       // Adjust the value (check fort sat/unsat, x_k is basic)
02594       b = getBeta(x_k);
02595       EpsRational x_k_new = beta[x_k] = b + theta * a_kj; 
02596       
02597       // Check if the variable is sat or unsat under the new assignment 
02598       if (x_k_new < getLowerBound(x_k) || getUpperBound(x_k) < x_k_new)
02599         unsatBasicVariables.insert(x_k);
02600       else
02601         unsatBasicVariables.erase(x_k);
02602     }
02603   
02604     // Go to the next one
02605     it ++;
02606   }
02607   
02608   // Last thing to do, pivot x_i and x_j
02609   pivot(x_i, x_j);
02610 }
02611       
02612 QueryResult TheoryArithNew::assertUpper(const Expr& x_i, const EpsRational& c, const Theorem& thm) {
02613   
02614   // If trace is on for the simplex, print it out
02615   TRACE("simplex", "TheoryArithNew::assertUpper(", x_i.toString() + ", " + c.toString(), ")");
02616 
02617   // Check if the given expression is actually a variable
02618   DebugAssert(isLeaf(x_i), "TheoryArithNew::assertUpper wronk kind, expected an arithmetic leaf (variable) got " + x_i.toString());
02619   
02620   // Check if the upper bound is worse then the last one
02621   if (getUpperBound(x_i) <= c) return (consistent == UNKNOWN? UNKNOWN : SATISFIABLE); 
02622   
02623   // If the new bound is lower then the lower bound */
02624   if (c < getLowerBound(x_i)) {
02625     // Create the explaining theorem
02626     explanation = d_rules->clashingBounds(getLowerBoundThm(x_i), thm);
02627     // Return unsatisfiable
02628     return UNSATISFIABLE;
02629   }
02630   
02631   // Since it is a tighter bound, propagate what can be 
02632   propagate = true;
02633 
02634   // Set the new upper bound */
02635   upperBound[x_i] = BoundInfo(c, thm);
02636 
02637   // If within the new bounds, return the previous state 
02638   if (getBeta(x_i) <= c) return (consistent == UNKNOWN? UNKNOWN : SATISFIABLE);
02639   
02640   // Otherwise, if the variable is non-basic then update it's value
02641   if (!isBasic(x_i)) update(x_i, c);
02642   // Else its basic, and non-sat, add to the unsat set
02643   else unsatBasicVariables.insert(x_i);
02644   
02645   // Everything went fine, return OK (not complete, means not UNSAT)
02646   return UNKNOWN;
02647 }
02648 
02649 QueryResult TheoryArithNew::assertLower(const Expr& x_i, const EpsRational& c, const Theorem& thm) {
02650 
02651   // Check if the given expression is actually a variable
02652   DebugAssert(isLeaf(x_i), "TheoryArithNew::assertLower wronk kind, expected an arithmetic leaf (variable) got " + x_i.toString());
02653 
02654   // If trace is on for the simplex, print it out
02655   TRACE("simplex", "TheoryArithNew::assertLower(", x_i.toString() + ", " + c.toString(), ")");
02656   
02657   // Check if the upper bound is worse then the last one
02658   if (c <= getLowerBound(x_i)) return (consistent == UNKNOWN? UNKNOWN : SATISFIABLE); 
02659   
02660   // Since it is a tighter bound, propagate what can be 
02661   propagate = true;
02662   
02663   // If the new bound is lower then the  bound */
02664   if (c > getUpperBound(x_i)) { 
02665     // Create the explaining theorem
02666     explanation = d_rules->clashingBounds(thm, getUpperBoundThm(x_i));
02667     // Return unsatisfiable
02668     return UNSATISFIABLE;
02669   }
02670   
02671   // Set the new upper bound */
02672   lowerBound[x_i] = BoundInfo(c, thm);
02673   
02674   // If within the new bounds, return the previous state 
02675   if (c <= getBeta(x_i)) return (consistent == UNKNOWN? UNKNOWN : SATISFIABLE);
02676   
02677   // Otherwise, if the variable is non-basic then update it's value
02678   if (!isBasic(x_i)) update(x_i, c);
02679   // Else its basic, and non-sat, add to the unsat set
02680   else unsatBasicVariables.insert(x_i);
02681   
02682   // Everything went fine, return OK (not complete, means not UNSAT)
02683   return UNKNOWN;
02684 }
02685       
02686 QueryResult TheoryArithNew::assertEqual(const Expr& x_i, const EpsRational& c, const Theorem& thm) {
02687   
02688   // Assert the upper bound
02689   consistent = assertUpper(x_i, c, thm); // TODO: thm: = |- <=
02690   
02691   // If the return value is UNSAT return UNSAT
02692   if (consistent == UNSATISFIABLE) return UNSATISFIABLE;
02693   
02694   // Assert the lower bound
02695   consistent = assertLower(x_i, c, thm); // TODO: thm: = |- >=
02696   
02697   // Return the consistency
02698   return consistent;
02699 }
02700 
02701 
02702 void TheoryArithNew::checkSat(bool fullEffort)
02703 { 
02704   // We should not be here if inconsistent
02705   DebugAssert(consistent != UNSATISFIABLE, "TheoryArithNew::checkSat : we are UNSAT before entering?!?!");
02706 
02707   // Trace the check sat if simplex flag is on
02708     TRACE("simplex", "TheoryArithNew::checkSat(fullEffort=",fullEffort? "true" : "false", ")");
02709   
02710     // Check if by backtracking we have more fresh variables than we expect
02711   if (assertedExprCount < assertedExpr.size()) 
02712     updateFreshVariables();
02713   
02714     // Do the simplex search if the database is not satisfiable (if inconsistent, we will not be here);
02715     if (consistent != SATISFIABLE || fullEffort)
02716       consistent = checkSatSimplex();
02717             
02718     // If the result is inconsistent set the inconsistent flag
02719     if (consistent == UNSATISFIABLE)
02720       setInconsistent(explanation);
02721     
02722     // If we are not - consistent, check the integer satisfiability
02723 //    if (consistent == SATISFIABLE) {
02724 //    // If we asserted a lemma and it hasn't been processed yet, we just don't do anything
02725 //    Theorem integer_lemma_thm = integer_lemma;
02726 //    if (!integer_lemma_thm.isNull()) {
02727 //      if (simplify(integer_lemma_thm.getExpr()).getRHS() == getEM()->trueExpr())
02728 //        consistent = checkSatInteger();
02729 //      else 
02730 //        consistent = UNKNOWN;
02731 //    } else
02732 //      // Lemma was not asserted, check integer sat
02733 //      consistent = checkSatInteger();   
02734 //    }
02735     
02736     // Trace the check sat if simplex flag is on
02737     TRACE("simplex", "TheoryArithNew::checkSat ==> ", consistent == UNSATISFIABLE? "UNSATISFIABLE" : "SATISFIABLE", "");
02738 }
02739 
02740 QueryResult TheoryArithNew::checkSatInteger() {
02741   
02742   // Trace the check sat if simplex flag is on
02743     TRACE("simplex", "TheoryArithNew::checkSatInteger()", "", "");
02744     
02745     // At this point we are REAL satisfiable, so we have to check the integers
02746   set<Expr>::iterator x_i_it     = intVariables.begin();
02747   set<Expr>::iterator x_i_it_end = intVariables.end();
02748 //  cerr << "Integer vars: ";
02749 //  while (x_i_it != x_i_it_end) {
02750 //    cerr << *x_i_it << " ";
02751 //    ++ x_i_it;
02752 //  }
02753 //  cerr << endl;
02754 //  
02755 //  x_i_it     = intVariables.begin();
02756 //  x_i_it_end = intVariables.end();
02757   while (x_i_it != x_i_it_end) {
02758     
02759     // Get the left-hand variable of the tableaux
02760     const Expr& x_i = (*x_i_it);
02761     
02762     // Only do work for unsat integer variables
02763     if (isInteger(x_i)) {
02764       
02765       // Get the current assignment of x_i
02766       const EpsRational& beta = getBeta(x_i);
02767       
02768       // Check if it is an integer
02769       if (beta.isInteger()) { ++ x_i_it; continue; }
02770       
02771       // Since the variable is not an integer, split on it being <= [beta] and >= [beta] + 1
02772       integer_lemma = d_rules->integerSplit(x_i, beta.getFloor()); 
02773       TRACE("integer", "TheoryArithNew::checkSatInteger ==> lemma : ", integer_lemma, "");
02774       enqueueFact(integer_lemma);
02775       
02776       // One split should be enough, return unknown
02777       return UNKNOWN;      
02778     }
02779     
02780     // Go to the next row
02781     ++ x_i_it;
02782   } 
02783   
02784   // We came through, huh, we are sat
02785   return SATISFIABLE; 
02786 }
02787 
02788 QueryResult TheoryArithNew::checkSatSimplex() {
02789   
02790   Expr x_i;                         // The basic variable we will use
02791   EpsRational x_i_Value;            // The current value of the variable x_i
02792   Expr x_j;                         // The non-basic variable we will use
02793   EpsRational x_j_Value;            // The current value of the variable x_j
02794   Rational a_ij;                    // The coefficient wwhen expanding x_i using x_j in the tableaux
02795   
02796   bool x_j_Found;                   // Did we find a suitable one
02797   EpsRational l_i;                  // Lower bound of x_i
02798   EpsRational u_i;                  // Upper bound of x_i
02799 
02800   // Trace the size of the tableaux and the unsat 
02801     TRACE("arith_atoms", "Tableaux size: ", tableaux.size(), "");
02802   TRACE("arith_atoms", "UNSAT vars: ", unsatBasicVariables.size(), "");
02803   
02804   // The main loop
02805   while (unsatBasicVariables.size()) {
02806    
02807     // The first unsat variable information
02808     x_i        = *unsatBasicVariables.begin();
02809     TebleauxMap::iterator it = tableaux.find(x_i); 
02810     x_i_Value  = getBeta(x_i);
02811     l_i        = getLowerBound(x_i);
02812     u_i        = getUpperBound(x_i);
02813     
02814     // If the variable doesn't obey the lower bound 
02815     if (l_i > x_i_Value) {
02816       
02817       // Did we find a suitable x_j, NOT YET
02818       x_j_Found = false; 
02819           
02820       // Find the smalles non-basic x_j that can improve x_i
02821       const Expr& x_i_RightSide = (*it).second.getExpr()[1];
02822       int non_basics_i_end     = x_i_RightSide.arity();
02823       for(int non_basics_i = 0; non_basics_i < non_basics_i_end; ++ non_basics_i) {
02824         
02825         // Get the info on the current variable
02826         x_j        = x_i_RightSide[non_basics_i][1]; 
02827         a_ij       = x_i_RightSide[non_basics_i][0].getRational(); 
02828         x_j_Value  = getBeta(x_j);           
02829         
02830         // If there is slack for improving x_i using x_j then do it               
02831         if (a_ij > 0) {
02832           if (x_j_Value < getUpperBound(x_j)) {
02833             x_j_Found = true;
02834             break;
02835           }
02836         } else {
02837           if (x_j_Value > getLowerBound(x_j)) {
02838             x_j_Found = true;
02839             break;
02840           }
02841         }     
02842       }
02843       
02844       // If none of the variables allows for slack, the tableaux is unsatisfiable
02845       if (!x_j_Found) {
02846         // Generate the explanation
02847         explanation = getLowerBoundExplanation(it);
02848         // Return unsat
02849         return UNSATISFIABLE;
02850       }
02851       
02852       // Otherwise do pivotAndUpdate and continue with the loop
02853       pivotAndUpdate(x_i, x_j, l_i);       
02854     } else
02855     // If the variable doesn't obey the upper bound 
02856     if (x_i_Value > u_i) {
02857       
02858       // Did we find a suitable x_j, NOT YET
02859       x_j_Found = false; 
02860           
02861       // Find the smalles non-basic x_j that can improve x_i
02862       const Expr& x_i_RightSide = (*it).second.getExpr()[1];
02863       int non_basics_i_end     = x_i_RightSide.arity();
02864       for(int non_basics_i = 0; non_basics_i < non_basics_i_end; ++ non_basics_i) {
02865         
02866         // Get the info on the current variable
02867         x_j        = x_i_RightSide[non_basics_i][1]; 
02868         a_ij       = x_i_RightSide[non_basics_i][0].getRational(); 
02869         x_j_Value  = getBeta(x_j);           
02870         
02871         // If there is slack for improving x_i using x_j then do it
02872         if (a_ij < 0) {
02873           if (x_j_Value < getUpperBound(x_j)) {
02874             x_j_Found = true;
02875             break;
02876           }
02877         } else 
02878           if (x_j_Value > getLowerBound(x_j)) {
02879             x_j_Found = true;
02880             break;
02881           }     
02882       }
02883       
02884       // If none of the variables allows for slack, the tableaux is unsatisfiable
02885       if (!x_j_Found) {
02886         // Generate the explanation
02887         explanation = getUpperBoundExplanation(it);
02888         // Return unsat
02889         return UNSATISFIABLE; 
02890       }
02891       
02892       // Otherwise do pivotAndUpdate and continue with the loop
02893       pivotAndUpdate(x_i, x_j, u_i);       
02894     } else
02895       // Due to backtracking a unsat might become sat, just erase it
02896       unsatBasicVariables.erase(x_i);
02897       
02898     
02899     // If trace is on, printout the current valuation and the tableaux
02900     TRACE("simplex", "TheoryArithNew::checkSatSimplex ==> new tableaux : \n", tableauxAsString(), "");
02901     TRACE("simplex", "TheoryArithNew::checkSatSimplex ==> new bounds   : \n", boundsAsString(), "");
02902     TRACE("simplex", "TheoryArithNew::checkSatSimplex ==> unsat   : \n", unsatAsString(), "");
02903     
02904   };
02905 
02906   // Since there is no more unsat constraints we are satisfiable
02907   return SATISFIABLE;
02908 }     
02909 
02910 Rational TheoryArithNew::getTableauxEntry(const Expr& x_i, const Expr& x_j) {
02911   return findCoefficient(x_j, tableaux[x_i].getExpr()[1]);
02912 }
02913 
02914 
02915 const Rational& TheoryArithNew::findCoefficient(const Expr& var, const Expr& expr) {
02916 
02917   // The zero coefficient 
02918   static Rational zeroCoefficient(0);
02919 
02920     // The given expression must be a sum
02921     DebugAssert(isPlus(expr), "TheoryArithNew::findCoefficient : expression must be a sum : " + expr.toString());
02922   
02923   // Go through the sum and find the coefficient
02924   int left = 0;
02925   int right = expr.arity() - 1;
02926   int i;
02927   while (left <= right) {
02928     
02929     // Take the middle one
02930     i = (left + right ) / 2;
02931     
02932     // Current expression
02933     const Expr& mult = expr[i];
02934     
02935     // Every summand must be a multiplication with a rational
02936     DebugAssert(isMult(mult) && isRational(mult[0]), "TheoryArithNew::findCoefficient : expression must be a sum of multiplications: " + expr.toString()); 
02937     
02938     // If the variable is the same return the coefficient
02939     int cmp = compare(mult[1], var);
02940     if (cmp == 0) 
02941       return mult[0].getRational();
02942     else 
02943       if (cmp > 0) 
02944         left = i + 1;
02945       else 
02946         right = i - 1;    
02947   }
02948     
02949   // Not found, return 0
02950   return zeroCoefficient; 
02951 }
02952 
02953 
02954 TheoryArithNew::EpsRational TheoryArithNew::getLowerBound(const Expr& x) const {
02955   // Try and find the lower bound in the map
02956   CDMap<Expr, BoundInfo>::iterator find = lowerBound.find(x);
02957   // If not found return +infinity as the default value, othervise return the value
02958   if (find == lowerBound.end()) return EpsRational::MinusInfinity;
02959   else return (*find).second.bound;
02960 }
02961 
02962 TheoryArithNew::EpsRational TheoryArithNew::getUpperBound(const Expr& x) const {
02963   // Try and find the upper bound in the map
02964   CDMap<Expr, BoundInfo>::iterator find = upperBound.find(x);
02965   // If not found return -infinity as the default value, othervise return the value
02966   if (find == upperBound.end()) return EpsRational::PlusInfinity;
02967   else return (*find).second.bound;
02968 }
02969 
02970 Theorem TheoryArithNew::getLowerBoundThm(const Expr& x) const {
02971   // Try and find the upper bound in the map
02972   CDMap<Expr, BoundInfo>::iterator find = lowerBound.find(x);
02973   // It has to be found
02974   DebugAssert(find != lowerBound.end(), "TheoryArithNew::getLowerBoundThm, theorem not found for " + x.toString()); 
02975   // Return the bound theorem 
02976   return (*find).second.theorem;
02977 }
02978 
02979 Theorem TheoryArithNew::getUpperBoundThm(const Expr& x) const {
02980   // Try and find the upper bound in the map
02981   CDMap<Expr, BoundInfo>::iterator find = upperBound.find(x);
02982   // It has to be found
02983   DebugAssert(find != upperBound.end(), "TheoryArithNew::getUpperBoundThm, theorem not found for " + x.toString()); 
02984   // Return the bound theorem 
02985   return (*find).second.theorem;
02986 }
02987 
02988       
02989 TheoryArithNew::EpsRational TheoryArithNew::getBeta(const Expr& x) {
02990   
02991   // Try to find the beta value in the map
02992   CDMap<Expr, EpsRational>::iterator find = beta.find(x);
02993   
02994   if (find == beta.end()) {
02995     // If not found return 0 (no need for sat/unsat, it's allways sat) 
02996     return beta[x] = EpsRational::Zero;
02997   
02998     // Check if the variable is sat or unsat under the new assignment 
02999     if (EpsRational::Zero < getLowerBound(x) || getUpperBound(x) < EpsRational::Zero)
03000       unsatBasicVariables.insert(x);
03001     else
03002       unsatBasicVariables.erase(x);
03003   }
03004   else 
03005     // If found, just return it
03006     return (*find).second;
03007 }
03008 
03009 
03010 void TheoryArithNew::assertFact(const Theorem& assertThm)
03011 {
03012   // Get the expression to be asserted
03013   const Expr& expr = assertThm.getExpr();
03014 
03015   // If tracing arithmetic, print out the expression to be asserted
03016   TRACE("simplex", "TheoryArithNew::assertFact(", expr, ")");
03017   TRACE("simplex", "asserted: ", assertedExpr.size(), "");
03018   TRACE("simplex", "counted: ", assertedExprCount, "");
03019   TRACE("arith_atoms", "Assert: ", expr.toString(), "");
03020             
03021   // Check if by backtracking we have more fresh variables than we expect
03022   if (assertedExprCount < assertedExpr.size()) 
03023     updateFreshVariables();
03024   
03025   // Get the constraint partsreturn 
03026   const Expr& leftSide   = expr[0];          // The left side of the constraint  
03027   const Expr& rightSide  = expr[1];          // The right side of the constraint
03028   int kind        = expr.getKind();   // The relational symbol, should be one of LT, LE, GT, GE. EQ was rewritten as LE and GE so its not here
03029   
03030   // The expression must be an inequality (sum a_1*x_1 .. a_n*x_n) @ b
03031   DebugAssert(isIneq(expr)         , "TheoryArithNew::assertFact wronk kind, expected inequality"                 + expr.toString());
03032   DebugAssert(isPlus(rightSide) || (isMult(rightSide) && isRational(rightSide[0]) && isLeaf(rightSide[1])), "TheoryArithNew::assertFact wronk kind, expected sum on the right side opr a simple 1*x"      + expr.toString());
03033   DebugAssert(isRational(leftSide), "TheoryArithNew::assertFact wronk kind, expected rational on the right side" + expr.toString());
03034 
03035   // Get the rational value on the right side 
03036   Rational leftSideValue = leftSide.getRational(); 
03037   
03038   // The rational bound on the constraint
03039   Rational c = leftSideValue;
03040   
03041   // The variable to be constrained
03042   Expr var;       
03043   
03044   // Theorem of the assertion
03045   Theorem assert = assertThm;                                 
03046   
03047   // For now we have that leftSide @ (c1x1 + c2x2 + ... cnxn) = rightSide
03048   // If rightSide is not a sum constraint is atomic (leftSide @ a*x)
03049   if (!isPlus(rightSide)) {
03050     
03051     // Tee left side in this case should be 1*x
03052     DebugAssert(rightSide[0].isRational() && rightSide[0].getRational() == 1, "TheoryArithNew::assertFact, left side should be multiplication by one");
03053     
03054     // Change the assert theorem to remove the 1*x to x
03055     assert = iffMP(assert, substitutivityRule(expr, reflexivityRule(leftSide), d_rules->oneElimination(rightSide)));
03056     
03057     // The variable will just be x1 (if var is not already present in the tableaux, it will be non-basic)
03058     var = rightSide[1]; // IMPORTANT!!!, when giving explanations, it should be of the same form
03059     
03060     // If an integer, add it to the integer set
03061     if (isInteger(var)) intVariables.insert(var);
03062           
03063   } else {
03064     // Get the theorem that corresponds to the introduction of the new variable var = leftSide 
03065     const Theorem& introductionThm = getVariableIntroThm(rightSide);
03066     
03067     // Take the new variable out
03068     var = introductionThm.getExpr()[0];
03069 
03070     // Change the theorem for the assertion so that it involves the new variable, i.e.  c < rightSide, var = rightSide |- c < var
03071     assert = iffMP(assertThm, substitutivityRule(expr, 1, symmetryRule(introductionThm)));
03072     
03073     // Insert all the integer variables into the integer set
03074     if (isInteger(var)) {
03075       intVariables.insert(var);
03076       int i_end = rightSide.arity();
03077       for(int i = 0; i < i_end; ++ i)
03078         intVariables.insert(rightSide[i][1]);
03079     } else {
03080       int i_end = rightSide.arity();
03081       for(int i = 0; i < i_end; ++ i)
03082         if (isInteger(rightSide[i][1])) intVariables.insert(rightSide[i][1]);
03083     }
03084   }
03085   
03086   // Make the BoundInfo object for theory propagation
03087   EpsRational bound;
03088   
03089   // By default we don't propagate
03090   propagate = false;
03091   
03092   // Remember the old bound
03093   EpsRational oldBound;
03094   
03095   // Finnaly assert the right constraint
03096   switch (kind) {
03097     case LT: 
03098       oldBound    = getLowerBound(var);
03099       // c < var, convert to c + epsilon <= var and assert it
03100       consistent = assertLower(var, bound = EpsRational(c, +1), assert);
03101       break;          
03102     case LE:
03103       oldBound    = getLowerBound(var);
03104       // c <= var, assert the lower bound
03105       consistent = assertLower(var, bound = c, assert);
03106       break;
03107     case GT: 
03108       oldBound    = getUpperBound(var);
03109       // c > var, convert to c - epsilon >= var and assert it
03110       consistent = assertUpper(var, bound = EpsRational(c, -1), assert);
03111       break;
03112     case GE:
03113       oldBound    = getUpperBound(var);
03114       // c >= var, assert the upper bound
03115       consistent = assertUpper(var, bound = c, assert);
03116       break;
03117     case EQ:
03118       // c == var, assert the equality
03119       consistent = assertEqual(var, bound = c, assert);
03120       // For now, don't propagate anything
03121       propagate = false; // TODO: some propagation is in place here (negations)     
03122       break;
03123     default:
03124       //How did we get here
03125           FatalAssert(false, "Theory_Arith::assertFact, control should not reach here");
03126           break;
03127   } 
03128   
03129   // If tracing arithmetic, print out the new tableaux and the current bounds
03130   TRACE("simplex", "TheoryArithNew::assertFact ==> ", consistent == UNSATISFIABLE? "UNSATISFIABLE" : consistent == UNKNOWN? "UNKNOWN" : "SATISFIABLE", "");
03131   TRACE("simplex", "TheoryArithNew::assertFact ==> new tableaux : \n", tableauxAsString(), "");
03132   TRACE("simplex", "TheoryArithNew::assertFact ==> new bounds   : \n", boundsAsString(), "");
03133   TRACE("simplex", "TheoryArithNew::assertFact ==> unsat   : \n", unsatAsString(), "");
03134   
03135   // If the result is inconsistent set the inconsistent flag
03136     if (consistent == UNSATISFIABLE)
03137       setInconsistent(explanation);
03138       
03139     // Not inconsistent, propagate all from this assertion
03140     if (propagate) 
03141       propagateTheory(assertThm.getExpr(), bound, oldBound);    
03142 }
03143 
03144 
03145 Theorem TheoryArithNew::getVariableIntroThm(const Expr& rightSide) {
03146   
03147   // Try to find the expression in the map
03148   TebleauxMap::iterator find = freshVariables.find(rightSide);
03149   
03150   // If not in tableaux, add it and assign it a new variable
03151   if (find == freshVariables.end()) {
03152   
03153     // Get the common rules
03154     CommonProofRules* c_rules = getCommonRules();
03155   
03156     // Create a new variable (\exists x . rightSide = x)
03157     Theorem thm = c_rules->varIntroRule(rightSide);
03158       
03159       // Skolemise it, to get an equality (\exists x . rightSide = x) <==> rightSide = c, then infer |- rightSide = c
03160       thm = c_rules->iffMP(thm, c_rules->skolemizeRewrite(thm.getExpr()));    
03161       
03162       // Reverse the equality into standard form
03163       thm = symmetryRule(thm);
03164       
03165       // Add the theorem to the variable introduction map (this is the theorem we return)
03166       Theorem returnThm = freshVariables[rightSide] = thm;
03167       
03168       // Now, flatten the equality modulo tableaux
03169       thm = substAndCanonizeModTableaux(thm);
03170       
03171       // Get the skolem constant that was introduced (|- c = leftSide)
03172       const Expr& var = thm.getExpr()[0];
03173             
03174       // Also add it to the tableaux
03175       tableaux[var] = thm;
03176 
03177     // Update the dependencies
03178     updateDependenciesAdd(var, thm.getExpr()[1]);
03179 
03180     // Add it to the expressions map
03181     assertedExpr.push_back(Expr(EQ, var, rightSide));
03182     assertedExprCount = assertedExprCount + 1;
03183       
03184       // Compute the value of the new variable using the tableaux
03185     updateValue(var, rightSide);
03186             
03187       // Return the variable
03188       return returnThm;
03189   }
03190   
03191   // Otherwise we have it, so return it
03192   return (*find).second;
03193 }
03194 
03195 void TheoryArithNew::updateValue(const Expr& var, const Expr& e) {
03196   
03197     // The initial value
03198     EpsRational varValue(0);
03199     
03200     // Go through the sum and compute the value
03201     int i_end = e.arity();
03202     for (int i = 0; i < i_end; ++ i)
03203       varValue += getBeta(e[i][1]) * e[i][0].getRational();
03204     
03205     // Update the beta  
03206     beta[var] = varValue;
03207     
03208     // Check if the variable is sat or unsat under the new assignment 
03209   if (varValue < getLowerBound(var) || getUpperBound(var) < varValue)
03210     unsatBasicVariables.insert(var);
03211   else
03212     unsatBasicVariables.erase(var);
03213 }
03214 
03215 string TheoryArithNew::tableauxAsString() const {
03216 
03217   // The string we are building 
03218   string str; 
03219 
03220   // Start from the begining
03221   TebleauxMap::const_iterator row     = tableaux.begin();
03222   TebleauxMap::const_iterator row_end = tableaux.end(); 
03223   
03224   // Print all the rows
03225   while (row != tableaux.end()) {
03226     // Print the equality
03227     str = str + ((*row).second).getExpr().toString() + "\n";
03228     
03229     // Continue to the next row
03230     row ++;
03231   }
03232   
03233   // Return the string representation 
03234   return str;
03235 }
03236 
03237 string TheoryArithNew::unsatAsString() const {
03238 
03239   // The string we are building 
03240   string str; 
03241 
03242   // Start from the begining
03243   set<Expr>::const_iterator it     = unsatBasicVariables.begin();
03244   set<Expr>::const_iterator it_end = unsatBasicVariables.end(); 
03245   
03246   // Print all the rows
03247   while (it != it_end) {
03248     // Print the equality
03249     str = str + (*it).toString() + " ";
03250     
03251     // Continue to the next row
03252     it ++;
03253   }
03254   
03255   // Return the string representation 
03256   return str;
03257 }
03258 
03259 
03260 string TheoryArithNew::boundsAsString() {
03261 
03262   // The string we are building 
03263   string str; 
03264 
03265   // Set containing all the variables
03266   set<Expr> all_variables;
03267   
03268   // Go throught the tableaux and pick up all the variables
03269   TebleauxMap::const_iterator it     = tableaux.begin();
03270   TebleauxMap::const_iterator it_end = tableaux.end();
03271   for(; it != it_end; it ++) {
03272   
03273     // Add the basic variable to the set
03274     all_variables.insert((*it).first);
03275     
03276     // Go through all the expression variables and add them to the set
03277     const Expr& rightSide = (*it).second.getExpr()[1];
03278     int term_i_end = rightSide.arity();
03279     for(int term_i = 0; term_i < term_i_end; ++ term_i)
03280       all_variables.insert(rightSide[term_i][1]); 
03281   } 
03282   
03283   // Go throught the bounds map and pickup all the variables
03284   CDMap<Expr, BoundInfo>::iterator bounds_it;
03285   for (bounds_it = lowerBound.begin(); bounds_it != lowerBound.end(); bounds_it ++)
03286     all_variables.insert((*bounds_it).first);
03287   for (bounds_it = upperBound.begin(); bounds_it != upperBound.end(); bounds_it ++)
03288     all_variables.insert((*bounds_it).first);
03289     
03290   // Start from the begining
03291   set<Expr>::iterator var_it     = all_variables.begin();
03292   set<Expr>::iterator var_it_end = all_variables.end(); 
03293   
03294   // Print all the rows
03295   while (var_it != var_it_end) {
03296     
03297     // The current variable
03298     const Expr& var = *var_it;
03299     
03300     // Print the equality
03301     str += getLowerBound(var).toString() + " <= " + var.toString() + "(" + getBeta(var).toString() + ") <= " + getUpperBound(var).toString() + "\n";
03302       
03303     // Continue to the next variable
03304     var_it ++;
03305   }
03306   
03307   // Return the string representation 
03308   return str;
03309 }
03310 
03311 // The infinity constant 
03312 const TheoryArithNew::EpsRational TheoryArithNew::EpsRational::PlusInfinity(PLUS_INFINITY);
03313 // The negative infinity constant
03314 const TheoryArithNew::EpsRational TheoryArithNew::EpsRational::MinusInfinity(MINUS_INFINITY);
03315 // The negative infinity constant
03316 const TheoryArithNew::EpsRational TheoryArithNew::EpsRational::Zero;
03317 
03318 
03319 Theorem TheoryArithNew::substAndCanonizeModTableaux(const Theorem& eq) {
03320 
03321   // If subst is empty, just return
03322     if(tableaux.empty()) return eq;
03323 
03324   // Get the expression of the equality
03325   const Expr& eqExpr = eq.getExpr();
03326 
03327   // Check if the equality if in the canonic form
03328   DebugAssert(eqExpr.getKind() == EQ, "TheoryArithNew::substAndCanonize, expected equality " + eqExpr.toString());
03329   
03330   // Get the left side of the equality
03331   const Expr& rightSide = eqExpr[1];
03332 
03333     // Do the actual substitution in the eqExpr
03334     Theorem thm = substAndCanonizeModTableaux(rightSide);
03335   
03336     // If the substitution had no result just return the original equation
03337     if(thm.getRHS() == rightSide) return eq;
03338     
03339     // Return the theorem 
03340     return iffMP(eq, substitutivityRule(eq.getExpr(), 1, thm));
03341 }
03342 
03343 Theorem TheoryArithNew::substAndCanonizeModTableaux(const Expr& sum) {
03344     
03345     Theorem res;                      // The resulting theorem    
03346     vector<Theorem> thms;             // The theorems of the changed terms for the substitution
03347     vector<unsigned> changed;         // The indices of the changed terms for the substitution
03348   
03349   // Trace the canonisation
03350   TRACE("simplex", "TheoryArithNew::substAndCanonizeModTableaux(", sum, ")");
03351   
03352   // Check if the equality if in the canonic form
03353   DebugAssert(sum.getKind() == PLUS, "TheoryArithNew::substAndCanonize, expected sum " + sum.toString());
03354    
03355     // Go throught the sum and try to substitute the variables
03356     int term_index_end = sum.arity();
03357     for(int term_index = 0; term_index < term_index_end; ++ term_index) {
03358 
03359     const Expr& term = sum[term_index]; // The current term expression (a*x)
03360     const Expr& var  = term[1];         // The variable of the current term
03361   
03362       // Find the variable in the map
03363       TebleauxMap::iterator find = tableaux.find(var);
03364       
03365       // If found, substitute it
03366       if (find != tableaux.end()) {
03367         
03368         // Substitute the variable
03369         Theorem termTheorem = canonThm(substitutivityRule(term, 1, (*find).second));
03370         
03371         // Check that the result is not trivial 
03372         DebugAssert(termTheorem.getExpr() != term, "substAndCanonizeModTableaux: got the same term in substitution");
03373           
03374       // Push it to the theorems vector
03375           thms.push_back(termTheorem);
03376           
03377           // Add the index to the changed vector
03378           changed.push_back(term_index);      
03379       }     
03380     }
03381   
03382     // Do the actual substitution and canonize the result
03383     if(thms.size() > 0) {
03384       // Substitute the changed subterms into the term
03385       res = substitutivityRule(sum, changed, thms);
03386       // Canonise the result
03387       res = canonThm(res);
03388     } else
03389       // Nothing happened, just return the reflexivity
03390       res = reflexivityRule(sum);     
03391   
03392     // Return the result
03393     return res;
03394 }
03395 
03396 void TheoryArithNew::substAndCanonizeTableaux(const Theorem& eq) {
03397 
03398   Theorem result; // The explaining theorem
03399     
03400   // Trace 
03401   TRACE("simplex", "TheoryArithNew::substAndCanonizeTableaux(", eq.getExpr(), ")");
03402     
03403   // Get the expression of the equality
03404   const Expr& eqExpr = eq.getExpr();
03405 
03406   // Check if the equality if in the canonic form
03407   DebugAssert(eqExpr.getKind() == EQ, "TheoryArithNew::substAndCanonize, expected equality " + eqExpr.toString());
03408   
03409   // Get the variable of the substitution
03410   const Expr& var = eqExpr[0];
03411 
03412   // Check if there are variables that depend on x_j
03413   DependenciesMap::iterator find = dependenciesMap.find(var);
03414   if (find != dependenciesMap.end()) {
03415       
03416     // Go through all the variables that depend on x_j, and update their value
03417     set<Expr>& dependent = (*find).second;
03418     set<Expr>::iterator it     = dependent.begin();
03419     set<Expr>::iterator it_end = dependent.end();     
03420     for(; it != it_end; ++ it) {
03421       
03422       // Get the expression and the right side of the row from the tableaux
03423       const Expr& leftSide      = *it;
03424       TebleauxMap::iterator row = tableaux.find(leftSide);
03425       const Expr& rowExpr       = (*row).second.getExpr();
03426       const Expr& rowRightSide  = rowExpr[1];
03427       
03428       // Go through the sum and try to substitute
03429       int right = rowRightSide.arity() - 1;
03430       int left  = 0;
03431       int term_i;
03432       while (left <= right) { 
03433         
03434         // Get the middle one
03435         term_i = (left + right) / 2;
03436         
03437         // Get the comparison of the variables
03438         int cmp = compare(rowRightSide[term_i][1], var);
03439               
03440         // If the variable is found
03441         if (cmp == 0) {
03442           
03443           // Do the substitution and canonise
03444           result = canonThm(substitutivityRule(rowRightSide[term_i], 1, eq));
03445           // Do the substitution and canonise in the sum
03446           result = canonThm(substitutivityRule(rowRightSide, term_i, result));
03447           // Do the substitution
03448           result = substitutivityRule(rowExpr, 1, result);
03449           
03450           // Get the new expression
03451           const Expr& newRowRightSide = result.getRHS()[1];       
03452           // Update the dependencies of the varriables in the expression
03453           updateDependencies(rowRightSide, newRowRightSide, leftSide, var); 
03454           
03455           // That's it, remember the new row
03456           (*row).second = iffMP((*row).second, result);       
03457           
03458           // Variables don't repeat, we can just break out        
03459           break;
03460         } else if (cmp > 0) 
03461           left = term_i + 1;
03462         else 
03463           right = term_i - 1;
03464       }
03465     }
03466     
03467     // Now nobody depends on var anymore, just erase it
03468       dependent.clear();
03469   } 
03470 }
03471 
03472 Theorem TheoryArithNew::pivotRule(const Theorem& eq, const Expr& var) {
03473 
03474   Theorem result;  // The theorem explaining the result
03475 
03476   // Get the expression
03477   const Expr& eqExpr = eq.getExpr();
03478   const Expr& right_side = eqExpr[1];
03479   const Expr& left_side = eqExpr[0];
03480 
03481   // Trace if askedTheorem canonLeft  = d_rules->canonMult(thm.getExpr()[0]);
03482       
03483   TRACE("simplex", "TheoryArithNew::pivotRule(", eqExpr.toString() + ", ", var.toString() + ")"); 
03484 
03485   // Eq must be an equation with the sum on the left side and a leaf on the right side (variable)
03486   DebugAssert(eqExpr.getKind() == EQ, "TheoryArithNew::pivotRule, expected an equality : " + eqExpr.toString());
03487   DebugAssert(right_side.getKind() == PLUS, "TheoryArithNew::pivotRule, expected a sum on the left-hand side : " + eqExpr.toString());
03488   DebugAssert(isLeaf(left_side), "TheoryArithNew::pivotRule, expected a leaf (variable) on the right-hand side : " + eqExpr.toString());
03489 
03490   // Find the term of var in the left-hand side of eq
03491   int term_i_end = right_side.arity();
03492   for(int term_i = 0; term_i < term_i_end; ++ term_i)
03493     // If found do the stuff
03494     if (right_side[term_i][1] == var) {
03495   
03496       // This is the term we need and the rational we need
03497       const Expr& termExpr = right_side[term_i];
03498       const Expr& termVar = termExpr[1]; 
03499       const Rational& termRat = termExpr[0].getRational();
03500       
03501       // Construct the expression we will add to the equality
03502       const Expr& minusTermExpr = multExpr(rat(-termRat), termVar);
03503       const Expr& minusVarExpr  = multExpr(rat(-1), left_side); 
03504     
03505       // Add the above expressions to the to the equality
03506       result = iffMP(eq, d_rules->plusPredicate(left_side, right_side, plusExpr(minusTermExpr, minusVarExpr), EQ));
03507       // Canonise the right-hand side
03508       result = transitivityRule(result, canon(result.getExpr()[1]));
03509       // Canonise the left-hand side 
03510       result = transitivityRule(symmetryRule(canon(result.getExpr()[0])), result);
03511       // Multiply by the inverse of the rational constant (negated and ^-1)
03512       result = iffMP(result, d_rules->multEqn(result.getExpr()[0], result.getExpr()[1], rat(-1/termRat)));
03513       // Canonise the left-hand side
03514       result = transitivityRule(result, canon(result.getExpr()[1]));
03515       // Canonise the right=hand side 
03516       result = transitivityRule(symmetryRule(canon(result.getExpr()[0])), result);
03517       // Rewrite 1*x => x in the left-hand side
03518       result = transitivityRule(symmetryRule(d_rules->oneElimination(result.getExpr()[0])), result);
03519       
03520       // Trace the result
03521       TRACE("simplex", "TheoryArithNew::pivotRule ==> ", result.getExpr().toString(), "");
03522       
03523       // We are done, there is just one variable there
03524       return result;
03525     }
03526   
03527   // If not found, there is something wrong
03528   DebugAssert(false, "TheoryArithNew::pivotRule, " + var.toString() + " does not occur in " + eqExpr.toString());
03529 
03530   // Dummy return 
03531   return result;
03532 }
03533 
03534 Theorem TheoryArithNew::getLowerBoundExplanation(const TebleauxMap::iterator& var_it) {
03535 
03536   vector<Theorem> upperBounds; // Vector of the upper-bound theorems 
03537   
03538   // Get the tableaux theorem explaining the leftside = var
03539   Theorem tableauxTheorem = (*var_it).second;
03540   
03541   // Get the variable on the right side
03542   const Expr& var = (*var_it).first;
03543 
03544   // Get the expression involved (leftSide = var)
03545   const Expr& rightSide = tableauxTheorem.getExpr()[1];
03546   
03547   // Go through the left side and pick up the apropriate lower and upper bounds
03548   int leftSide_i_end = rightSide.arity();
03549   for(int leftSide_i = 0; leftSide_i < leftSide_i_end; ++ leftSide_i) {
03550 
03551     // Get the rational
03552     const Expr& a = rightSide[leftSide_i][0];
03553     
03554     // Get the variable
03555     const Expr& x = rightSide[leftSide_i][1];
03556     
03557     // The positive ones make up the upper bounds (x < u => a * x < a * u)
03558     if (a.getRational() > 0) {
03559       // Get the upper bound x < u
03560       Theorem thm = getUpperBoundThm(x);
03561       
03562       // Multiply if by a ==> u_i * a > x * a
03563       thm = iffMP(thm, d_rules->multIneqn(thm.getExpr(), a));
03564     
03565       // Canonise the left and the right side
03566       Theorem canonRight  = d_rules->canonMultTermConst(thm.getExpr()[1][1], thm.getExpr()[1][0]);
03567       Theorem canonLeft = d_rules->canonMultConstConst(thm.getExpr()[0][0], thm.getExpr()[0][1]); 
03568       
03569       // Substitute the canonised to the inequality
03570       thm = iffMP(thm, substitutivityRule(thm.getExpr(), canonLeft, canonRight));
03571       
03572       // Add the bound to the vector of upper bounds (|- c_x > a * x) 
03573       upperBounds.push_back(thm);
03574     } 
03575     // The negative ones make up the lower bounds (x > l => a * x < a * l)
03576     else {
03577       // Get the lower bound l < x
03578       Theorem thm = getLowerBoundThm(x);
03579       
03580       // Multiply if by a |- l * a < x * a  
03581       thm = iffMP(thm, d_rules->multIneqn(thm.getExpr(), a));
03582     
03583       // Canonise the left and the right side
03584       Theorem canonRight  = d_rules->canonMultTermConst(thm.getExpr()[1][1], thm.getExpr()[1][0]);
03585       Theorem canonLeft = d_rules->canonMultConstConst(thm.getExpr()[0][0], thm.getExpr()[0][1]); 
03586       
03587       // Substitute the canonised to the inequality
03588       thm = iffMP(thm, substitutivityRule(thm.getExpr(), canonLeft, canonRight));
03589       
03590       // Add the bound to the vector of upper bounds (|- c_x > a * x)
03591       upperBounds.push_back(thm);
03592     }
03593   } 
03594   
03595   // Add up all the inequalities to get a C = \sum c_i > rightSide 
03596   Theorem sumInequalities = upperBounds[0];
03597   for(unsigned int i = 1; i < upperBounds.size(); i ++) { 
03598     // Add the inequalities
03599     sumInequalities = d_rules->addInequalities(sumInequalities, upperBounds[i]);
03600     
03601     // Canonise the left and the right side
03602     Theorem canonLeft  = d_rules->canonPlus(sumInequalities.getExpr()[0]);
03603     Theorem canonRight = d_rules->canonPlus(sumInequalities.getExpr()[1]);
03604     
03605     // Substitute the canonised to the inequality
03606     sumInequalities = iffMP(sumInequalities, substitutivityRule(sumInequalities.getExpr(), canonLeft, canonRight));   
03607   }
03608   
03609   // Substitute to get C > rightSide ==> C > var
03610   Theorem varUpperBound = substitutivityRule(sumInequalities.getExpr(), 1, symmetryRule(tableauxTheorem));  
03611   // MP to get C > var
03612   varUpperBound = iffMP(sumInequalities, varUpperBound);
03613   
03614   // Get the lower bound on the rigth side variable (l_var < var)
03615   Theorem varLowerBound = getLowerBoundThm(var);
03616 
03617   // Return the clashing bound theorem 
03618   return d_rules->clashingBounds(varLowerBound, varUpperBound);
03619 }
03620 
03621 Theorem TheoryArithNew::getUpperBoundExplanation(const TebleauxMap::iterator& var_it) {
03622 
03623   vector<Theorem> lowerBounds; // Vector of the upper-bound theorems 
03624   
03625   // Get the tableaux theorem explaining the leftside = var
03626   Theorem tableauxTheorem = (*var_it).second;
03627   
03628   // Get the variable on the right side
03629   const Expr& var = (*var_it).first;
03630 
03631   // Get the expression involved (var = rightSide)
03632   const Expr& rightSide = tableauxTheorem.getExpr()[1];
03633   
03634   // Trace if requested
03635   TRACE("explanations", "Generating explanation for the tableaux row ", tableauxTheorem.getExpr(), "");
03636   
03637   // Go through the right side and pick up the apropriate lower and upper bounds
03638   int rightSide_i_end = rightSide.arity();
03639   for(int rightSide_i = 0; rightSide_i < rightSide_i_end; ++ rightSide_i) {
03640 
03641     // Get the rational
03642     const Expr& a = rightSide[rightSide_i][0];
03643     
03644     // Get the variable
03645     const Expr& x = rightSide[rightSide_i][1];
03646     
03647     // The positive ones make up the lower bounds (x > l => a * x > a * l)
03648     if (a.getRational() > 0) {
03649       // Get the lower bound l < x
03650       Theorem thm = getLowerBoundThm(x);
03651       TRACE("explanations", "Got ", thm.getExpr(), "");
03652       
03653       // Multiply if by a ==> l * a < x * a
03654       thm = iffMP(thm, d_rules->multIneqn(thm.getExpr(), a));
03655       TRACE("explanations", "Got ", thm.getExpr(), "");
03656     
03657       // Canonise the left and the right side
03658       Theorem canonRight  = d_rules->canonMultTermConst(thm.getExpr()[1][1], thm.getExpr()[1][0]);
03659       Theorem canonLeft = d_rules->canonMultConstConst(thm.getExpr()[0][0], thm.getExpr()[0][1]); 
03660       
03661       // Substitute the canonised to the inequality
03662       thm = iffMP(thm, substitutivityRule(thm.getExpr(), canonLeft, canonRight));
03663       TRACE("explanations", "Got ", thm.getExpr(), "");
03664       
03665       // Add the bound to the vector of upper bounds (|- c_x < a * x) 
03666       lowerBounds.push_back(thm);
03667     } 
03668     // The negative ones make up the upper bounds (x < u => a * x > a * u)
03669     else {
03670       // Get the upper bound u > x
03671       Theorem thm = getUpperBoundThm(x);
03672       TRACE("explanations", "Got ", thm.getExpr(), "");
03673     
03674       // Multiply it by a |- u * a > x * a 
03675       thm = iffMP(thm, d_rules->multIneqn(thm.getExpr(), a));
03676       TRACE("explanations", "Got ", thm.getExpr(), "");
03677     
03678       // Canonise the left and the right side
03679       Theorem canonRight = d_rules->canonMultTermConst(thm.getExpr()[1][1], thm.getExpr()[1][0]);
03680       Theorem canonLeft  = d_rules->canonMultConstConst(thm.getExpr()[0][0], thm.getExpr()[0][1]); 
03681       
03682       // Substitute the canonised to the inequality
03683       thm = iffMP(thm, substitutivityRule(thm.getExpr(), canonLeft, canonRight));
03684       TRACE("explanations", "Got ", thm.getExpr(), "");
03685       
03686       // Add the bound to the vector of upper bounds (|- c_x < a * x)
03687       lowerBounds.push_back(thm);
03688     }
03689   } 
03690   
03691   // Add up all the inequalities to get a \sum c_i = C > rightSide
03692   Theorem sumInequalities = lowerBounds[0];
03693   for(unsigned int i = 1; i < lowerBounds.size(); i ++) { 
03694     // Add the inequalities
03695     sumInequalities = d_rules->addInequalities(sumInequalities, lowerBounds[i]);
03696     
03697     TRACE("explanations", "Got sum ", sumInequalities.getExpr(), "");
03698     
03699     // Canonise the left and the right side
03700     Theorem canonLeft  = d_rules->canonPlus(sumInequalities.getExpr()[0]);
03701     Theorem canonRight = d_rules->canonPlus(sumInequalities.getExpr()[1]);
03702     
03703     // Substitute the canonised to the inequality
03704     sumInequalities = iffMP(sumInequalities, substitutivityRule(sumInequalities.getExpr(), canonLeft, canonRight));   
03705   }
03706   
03707   // Trace if requested
03708   TRACE("explanations", "Got sum ", sumInequalities.getExpr(), "");
03709 
03710   // Substitute to get C < leftSide ==> C < var
03711   Theorem varLowerBound = substitutivityRule(sumInequalities.getExpr(), 1, symmetryRule(tableauxTheorem));  
03712   
03713   // MP to get C < var
03714   varLowerBound = iffMP(sumInequalities, varLowerBound);
03715 
03716   // Trace if requested
03717   TRACE("explanations", "Got lower bound ", varLowerBound.getExpr(), "");
03718   
03719   // Get the lower bound on the rigth side variable (var > l_var)
03720   Theorem varUpperBound = getUpperBoundThm(var);
03721 
03722   // Trace if requested
03723   TRACE("explanations", "Got upper bound ", varUpperBound.getExpr(), "");
03724 
03725   TRACE("explanations", "The var value (", var, ")" + getBeta(var).toString());
03726 
03727   // Return the clashing bound theorem 
03728   return d_rules->clashingBounds(varLowerBound, varUpperBound);
03729 }
03730 
03731 void TheoryArithNew::updateFreshVariables() {
03732 
03733   unsigned int size = assertedExpr.size();
03734   unsigned int i;
03735 
03736   for (i = assertedExprCount; i < size; ++ i)
03737     //Update the value
03738     updateValue(assertedExpr[i][0], assertedExpr[i][1]);
03739 
03740   // Update the asserted count to be the size of the vector
03741   assertedExprCount = i;
03742 
03743 }
03744 
03745 void TheoryArithNew::updateDependenciesAdd(const Expr& var, const Expr& sum) {
03746   
03747   // Trace it
03748   TRACE("tableaux_dep", "updateDependenciesAdd(", var.toString() + ", ", sum.toString() + ")");
03749   
03750   // Go through the sum and add var to the dependencies of that term variable
03751   Expr::iterator term = sum.begin();
03752   Expr::iterator term_end = sum.end();
03753   for(; term != term_end; term ++)
03754     dependenciesMap[(*term)[1]].insert(var);
03755 
03756 }
03757 
03758 void TheoryArithNew::updateDependenciesRemove(const Expr& var, const Expr& sum) {
03759   
03760   // Trace it
03761   TRACE("tableaux_dep", "updateDependenciesRemove(", var.toString() + ", ", sum.toString() + ")");
03762     
03763   // Go through the sum and remove var to the dependencies of that term variable
03764   Expr::iterator term = sum.begin();
03765   Expr::iterator term_end = sum.end();
03766   for(; term != term_end; term ++)
03767     dependenciesMap[(*term)[1]].erase(var);
03768 
03769 }
03770 
03771 void TheoryArithNew::updateDependencies(const Expr& oldSum, const Expr& newSum, const Expr& dependentVar, const Expr& skipVar) {
03772 
03773   // Trace it
03774   TRACE("tableaux_dep", "updateDependencies(", oldSum.toString() + ", " + newSum.toString() + ", " + dependentVar.toString(), ")");
03775   
03776   // Since the sums are ordered by variables, we can just to a "merge sort"
03777   int oldSum_i = 0, newSum_i = 0;
03778   int oldSum_end = oldSum.arity(), newSum_end = newSum.arity();
03779   // Get the first variables
03780   while (oldSum_i < oldSum_end && newSum_i < newSum_end) {
03781     
03782     // Get the variable references
03783     const Expr oldVar = oldSum[oldSum_i][1];
03784     const Expr newVar = newSum[newSum_i][1];
03785     
03786     // If variables equal, just skip, everything is ok
03787     if (oldVar == newVar) {
03788       ++ oldSum_i; ++ newSum_i; continue; 
03789     } 
03790     
03791     // Otherwise do the work with the smaller one
03792     if (oldVar > newVar) {
03793       // Old variable has dissapeared, remove dependent from its list 
03794       if (oldVar != skipVar)
03795         dependenciesMap[oldVar].erase(dependentVar);
03796       // Move the old variable forward
03797       ++ oldSum_i; 
03798     } else { 
03799       // New variable has appeared, insert dependent to its list
03800       if (newVar != skipVar)
03801         dependenciesMap[newVar].insert(dependentVar);
03802       // Move the new variable forward
03803       ++ newSum_i;
03804     }
03805   } 
03806 
03807   // If there is leftovers in the old sum, just do the removals
03808   while (oldSum_i < oldSum_end) {
03809     // Get the var, and increase the index
03810     const Expr& var = oldSum[oldSum_i ++][1];
03811     // Update the dependency
03812     if (var != skipVar) 
03813       dependenciesMap[var].erase(dependentVar);
03814   }
03815   
03816   while (newSum_i < newSum_end) {
03817     // Get the var, and increase the index
03818     const Expr& var = newSum[newSum_i ++][1];
03819     // Update the dependency
03820     if (var != skipVar) 
03821       dependenciesMap[var].insert(dependentVar);
03822   }
03823 }
03824 
03825 void TheoryArithNew::registerAtom(const Expr& e) {
03826   
03827   // Trace it
03828   TRACE("propagate_arith", "registerAtom(", e.toString(), ")");
03829   TRACE("arith_atoms", "", e.toString(), "");
03830   
03831   // If it is a atomic formula, add it to the map 
03832   if (e.isAbsAtomicFormula()) {
03833     Expr rightSide    = e[1];
03834     int kind          = e.getKind();
03835     Rational leftSide = e[0].getRational();
03836     
03837     // The eps bound we'll be using
03838     EpsRational bound; 
03839     
03840     // Depending on the type of the inequality define the bound
03841     switch (kind) {
03842       case LT: 
03843         bound = EpsRational(leftSide, +1);
03844         break;          
03845       case LE:
03846         bound = leftSide;
03847         break;
03848       case GT: 
03849         bound = EpsRational(leftSide, -1);
03850         break;
03851       case GE:
03852         bound = leftSide;
03853         break;      
03854       default:
03855         // How did we get here
03856             FatalAssert(false, "TheoryArithNew::registerAtom: control should not reach here" + e.toString()); 
03857     }
03858     
03859     // Bound has been set, now add the
03860     allBounds.insert(ExprBoundInfo(bound, e));    
03861   } 
03862 
03863 //  // Printout the current set of atoms in the set
03864 //  cout << "ALL BOUNDS:" << endl;
03865 //  BoundInfoSet::iterator it = allBounds.begin();
03866 //  while (it != allBounds.end()) {
03867 //    cout << (*it).e << endl;
03868 //    ++ it;
03869 //  }
03870 }
03871 
03872 void TheoryArithNew::propagateTheory(const Expr& assertExpr, const EpsRational& bound, const EpsRational& oldBound) {
03873 
03874   // Trace it
03875   TRACE("propagate_arith", "propagateTheory(", assertExpr.toString() + ", " + bound.toString(), ")");
03876 
03877   // Make the bound info object, so that we can search for it
03878   ExprBoundInfo boundInfo(bound, assertExpr);
03879   
03880     // Get the exression on the right side hand
03881   Expr rightSide = assertExpr[1];
03882   
03883   // Get the kid of the disequality
03884   int kind = assertExpr.getKind();
03885   
03886   // Check wheather the kind is one of LT, LE, GT, GE
03887   DebugAssert(kind == LT || kind == LE || kind == GT || kind == GE , "TheoryArithNew::propagateTheory : expected an inequality");
03888   
03889   // If the bound is of the type LT or LE we go up
03890   if (kind == LT || kind == LE) {
03891     // Find the asserted fact in the set (it must be there)
03892     BoundInfoSet::iterator find      = allBounds.lower_bound(boundInfo);
03893     BoundInfoSet::iterator find_end  = allBounds.lower_bound(ExprBoundInfo(oldBound, assertExpr));
03894     
03895     // If we are at the begining, or not found, just exit
03896     if (find == find_end) return;
03897     
03898     // Now, go up until reached wrong right side (skip the first one, its the same as given)
03899     while (find != find_end) {
03900       
03901       // Go up;
03902       -- find;
03903     
03904       // Get the theorem of the find
03905       const Expr& findExpr = (*find).e;
03906       
03907       // Get the bound of the find
03908       const EpsRational findRat = (*find).bound; 
03909     
03910       // Get the kind of the expression in the theorem
03911       int findExprKind = findExpr.getKind();
03912       
03913       // Check if the right sides are the same
03914       if (rightSide != findExpr[1]) break;
03915     
03916       // Construct the theorem object
03917       Theorem lemma;
03918     
03919       // Check the type and equeue the lemma 
03920       if (findExprKind == LT || findExprKind == LE)
03921         // Imply the weaker inequality
03922         lemma = d_rules->implyWeakerInequality(assertExpr, findExpr);
03923       else
03924         // Imply the negation of the impossible inequality
03925         lemma = d_rules->implyNegatedInequality(assertExpr, findExpr);
03926         
03927     
03928       TRACE("propagate_arith", "lemma ==>", lemma.toString(), "");
03929       TRACE("arith_atoms", "Propagate: ", lemma.getExpr().toString(), "");
03930     
03931       // Put it across
03932       enqueueFact(lemma);
03933     }
03934   }  
03935   // If the bound is of the type GT or GE we go down
03936   else {
03937     // Find the asserted fact in the set (it must be there)
03938     BoundInfoSet::iterator find          = allBounds.upper_bound(boundInfo);
03939     BoundInfoSet::iterator find_end      = allBounds.upper_bound(ExprBoundInfo(oldBound, assertExpr));
03940     
03941     // Now, go down until reached wrong right side (skip the first one, its the same as given)
03942     while (find != find_end) {
03943               
03944       // Get the bound of the find
03945       const EpsRational findRat = (*find).bound; 
03946     
03947       // Get the expression in the theorem
03948       const Expr& findExpr = (*find).e;
03949       int findExprKind = findExpr.getKind();
03950       
03951       // Check if the right sides are the same
03952       if (rightSide != findExpr[1]) break;
03953     
03954       // Construct the theorem object
03955       Theorem lemma;
03956     
03957       // Check the type and equeue the lemma 
03958       if (findExprKind == GT || findExprKind == GE)
03959         // Imply the weaker inequality
03960         lemma = d_rules->implyWeakerInequality(assertExpr, findExpr);
03961       else
03962         // Imply the negation of the impossible inequality (use oposite than above)
03963         lemma = d_rules->implyNegatedInequality(assertExpr, findExpr);
03964             
03965       TRACE("propagate_arith", "lemma ==>", lemma.toString(), "");
03966       TRACE("arith_atoms", "Propagate: ", lemma.getExpr().toString(), "");
03967         
03968       // Put it across
03969       enqueueFact(lemma);
03970       
03971       // Go to the next one
03972       ++ find;
03973     }
03974   } 
03975 }
03976 
03977 Theorem TheoryArithNew::deriveGomoryCut(const Expr& x_i) {
03978 
03979   Theorem res;
03980 
03981   // CHECK IF APPLICABLE
03982   DebugAssert(isBasic(x_i), "TheoryArithNew::deriveGomoryCut variable must be a basic variable : " + x_i.toString()); 
03983   DebugAssert(intVariables.count(x_i) > 0, "TheoryArithNew::deriveGomoryCut variable must be a basic variable : " + x_i.toString());
03984 
03985   // Get the rational value of x_i
03986   Rational x_i_Value = getBeta(x_i).getRational();
03987   
03988   // Compute f_0
03989   Rational f_0 = x_i_Value - floor(x_i_Value);
03990     
03991   return res;
03992 }

Generated on Tue Jul 3 14:33:39 2007 for CVC3 by  doxygen 1.5.1