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