Branch data Line data Source code
1 : : /******************************************************************************
2 : : * This file is part of the cvc5 project.
3 : : *
4 : : * Copyright (c) 2009-2026 by the authors listed in the file AUTHORS
5 : : * in the top-level source directory and their institutional affiliations.
6 : : * All rights reserved. See the file COPYING in the top-level source
7 : : * directory for licensing information.
8 : : * ****************************************************************************
9 : : *
10 : : * Diophantine equation solver
11 : : *
12 : : * A Diophantine equation solver for the theory of arithmetic.
13 : : */
14 : : #include "theory/arith/linear/dio_solver.h"
15 : :
16 : : #include <iostream>
17 : :
18 : : #include "base/output.h"
19 : : #include "expr/skolem_manager.h"
20 : : #include "options/arith_options.h"
21 : : #include "smt/env.h"
22 : : #include "theory/arith/linear/partial_model.h"
23 : :
24 : : using namespace std;
25 : :
26 : : namespace cvc5::internal {
27 : : namespace theory {
28 : : namespace arith::linear {
29 : :
30 : 12960 : inline Node makeIntegerVariable(NodeManager* nm)
31 : : {
32 : 12960 : return NodeManager::mkDummySkolem("intvar", nm->integerType());
33 : : }
34 : :
35 : 50263 : DioSolver::DioSolver(Env& env)
36 : : : EnvObj(env),
37 : 50263 : d_lastUsedProofVariable(context(), 0),
38 : 50263 : d_inputConstraints(context()),
39 : 50263 : d_nextInputConstraintToEnqueue(context(), 0),
40 : 50263 : d_trail(context()),
41 : 50263 : d_subs(context()),
42 : 50263 : d_currentF(),
43 : 50263 : d_savedQueue(context()),
44 : 50263 : d_savedQueueIndex(context(), 0),
45 : 50263 : d_conflictIndex(context()),
46 : 50263 : d_maxInputCoefficientLength(context(), 0),
47 : 50263 : d_usedDecomposeIndex(context(), false),
48 : 50263 : d_lastPureSubstitution(context(), 0),
49 : 50263 : d_pureSubstitionIter(context(), 0),
50 : 50263 : d_decompositionLemmaQueue(context()),
51 : 150789 : d_statistics(statisticsRegistry())
52 : : {
53 : 50263 : }
54 : :
55 : 50263 : DioSolver::Statistics::Statistics(StatisticsRegistry& sr)
56 : 50263 : : d_conflictCalls(sr.registerInt("theory::arith::dio::conflictCalls")),
57 : 50263 : d_cutCalls(sr.registerInt("theory::arith::dio::cutCalls")),
58 : 50263 : d_cuts(sr.registerInt("theory::arith::dio::cuts")),
59 : 50263 : d_conflicts(sr.registerInt("theory::arith::dio::conflicts")),
60 : 50263 : d_conflictTimer(sr.registerTimer("theory::arith::dio::conflictTimer")),
61 : 50263 : d_cutTimer(sr.registerTimer("theory::arith::dio::cutTimer"))
62 : : {
63 : 50263 : }
64 : :
65 : 90168 : bool DioSolver::queueConditions(TrailIndex t){
66 [ + - ]: 90168 : Trace("queueConditions") << !inConflict() << std::endl;
67 [ + - ]: 90168 : Trace("queueConditions") << gcdIsOne(t) << std::endl;
68 [ + - ]: 90168 : Trace("queueConditions") << !debugAnySubstitionApplies(t) << std::endl;
69 [ + - ]: 90168 : Trace("queueConditions") << !triviallySat(t) << std::endl;
70 [ + - ]: 90168 : Trace("queueConditions") << !triviallyUnsat(t) << std::endl;
71 : :
72 : : return
73 : 90168 : !inConflict() &&
74 [ + - ]: 90168 : gcdIsOne(t) &&
75 [ + - ]: 90168 : !debugAnySubstitionApplies(t) &&
76 [ + - ][ + - ]: 270504 : !triviallySat(t) &&
77 [ + - ]: 180336 : !triviallyUnsat(t);
78 : : }
79 : :
80 : 42326 : size_t DioSolver::allocateProofVariable() {
81 [ - + ][ - + ]: 42326 : Assert(d_lastUsedProofVariable <= d_proofVariablePool.size());
[ - - ]
82 [ + + ]: 42326 : if(d_lastUsedProofVariable == d_proofVariablePool.size()){
83 [ - + ][ - + ]: 11486 : Assert(d_lastUsedProofVariable == d_proofVariablePool.size());
[ - - ]
84 : 11486 : Node intVar = makeIntegerVariable(nodeManager());
85 : 11486 : d_proofVariablePool.push_back(Variable(intVar));
86 : 11486 : }
87 : 42326 : size_t res = d_lastUsedProofVariable;
88 : 42326 : d_lastUsedProofVariable = d_lastUsedProofVariable + 1;
89 : 42326 : return res;
90 : : }
91 : :
92 : :
93 : 0 : Node DioSolver::nextPureSubstitution(){
94 : 0 : Assert(hasMorePureSubstitutions());
95 : 0 : SubIndex curr = d_pureSubstitionIter;
96 : 0 : d_pureSubstitionIter = d_pureSubstitionIter + 1;
97 : :
98 : 0 : Assert(d_subs[curr].d_fresh.isNull());
99 : 0 : Variable v = d_subs[curr].d_eliminated;
100 : :
101 : 0 : SumPair sp = d_trail[d_subs[curr].d_constraint].d_eq;
102 : 0 : Polynomial p = sp.getPolynomial();
103 : 0 : Constant c = -sp.getConstant();
104 : 0 : Polynomial cancelV = p + Polynomial::mkPolynomial(v);
105 : 0 : Node eq = nodeManager()->mkNode(Kind::EQUAL, v.getNode(), cancelV.getNode());
106 : 0 : return eq;
107 : 0 : }
108 : :
109 : :
110 : 42326 : bool DioSolver::debugEqualityInInputEquations(Node eq){
111 : : typedef context::CDList<InputConstraint>::const_iterator const_iterator;
112 : 42326 : const_iterator i=d_inputConstraints.begin(), end = d_inputConstraints.end();
113 [ + + ]: 1224742 : for(; i != end; ++i){
114 : 1182416 : Node reason_i = (*i).d_reason;
115 [ - + ]: 1182416 : if(eq == reason_i){
116 : 0 : return true;
117 : : }
118 [ + - ]: 1182416 : }
119 : 42326 : return false;
120 : : }
121 : :
122 : :
123 : 42326 : void DioSolver::pushInputConstraint(const Comparison& eq, Node reason){
124 [ - + ][ - + ]: 42326 : Assert(!debugEqualityInInputEquations(reason));
[ - - ]
125 [ - + ][ - + ]: 42326 : Assert(eq.debugIsIntegral());
[ - - ]
126 [ - + ][ - + ]: 42326 : Assert(eq.getNode().getKind() == Kind::EQUAL);
[ - - ]
127 : :
128 : 42326 : SumPair sp = eq.toSumPair();
129 [ - + ][ - + ]: 42326 : Assert(!sp.isNonlinear());
[ - - ]
130 : :
131 : :
132 : :
133 : 42326 : uint32_t length = sp.maxLength();
134 [ + + ]: 42326 : if(length > d_maxInputCoefficientLength){
135 : 4726 : d_maxInputCoefficientLength = length;
136 : : }
137 : :
138 : 42326 : size_t varIndex = allocateProofVariable();
139 : 42326 : Variable proofVariable(d_proofVariablePool[varIndex]);
140 : : // Variable proofVariable(makeIntegerVariable(nodeManager()));
141 : :
142 : 42326 : TrailIndex posInTrail = d_trail.size();
143 [ + - ]: 84652 : Trace("dio::pushInputConstraint") << "pushInputConstraint @ " << posInTrail
144 : 0 : << " " << eq.getNode()
145 : 42326 : << " " << reason << endl;
146 : 42326 : d_trail.push_back(Constraint(sp,Polynomial::mkPolynomial(proofVariable)));
147 : :
148 : 42326 : size_t posInConstraintList = d_inputConstraints.size();
149 : 42326 : d_inputConstraints.push_back(InputConstraint(reason, posInTrail));
150 : :
151 : 42326 : d_varToInputConstraintMap[proofVariable.getNode()] = posInConstraintList;
152 : 42326 : }
153 : :
154 : :
155 : 26335 : DioSolver::TrailIndex DioSolver::scaleEqAtIndex(DioSolver::TrailIndex i, const Integer& g){
156 [ - + ][ - + ]: 26335 : Assert(g != 0);
[ - - ]
157 : 52670 : Constant invg = Constant::mkConstant(nodeManager(), Rational(Integer(1), g));
158 : 26335 : const SumPair& sp = d_trail[i].d_eq;
159 : 26335 : const Polynomial& proof = d_trail[i].d_proof;
160 : :
161 : 26335 : SumPair newSP = sp * invg;
162 : 26335 : Polynomial newProof = proof * invg;
163 : :
164 [ - + ][ - + ]: 26335 : Assert(newSP.isIntegral());
[ - - ]
165 [ - + ][ - + ]: 26335 : Assert(newSP.gcd() == 1);
[ - - ]
166 : :
167 : 26335 : TrailIndex j = d_trail.size();
168 : :
169 : 26335 : d_trail.push_back(Constraint(newSP, newProof));
170 : :
171 [ + - ]: 26335 : Trace("arith::dio") << "scaleEqAtIndex(" << i <<","<<g<<")"<<endl;
172 [ + - ]: 52670 : Trace("arith::dio") << "derived "<< newSP.getNode()
173 : 26335 : <<" with proof " << newProof.getNode() << endl;
174 : 26335 : return j;
175 : 26335 : }
176 : :
177 : 244 : Node DioSolver::proveIndex(TrailIndex i){
178 [ - + ][ - + ]: 244 : Assert(inRange(i));
[ - - ]
179 : 244 : const Polynomial& proof = d_trail[i].d_proof;
180 [ - + ][ - + ]: 244 : Assert(!proof.isConstant());
[ - - ]
181 : :
182 : 244 : NodeBuilder nb(nodeManager(), Kind::AND);
183 [ + + ]: 975 : for(Polynomial::iterator iter = proof.begin(), end = proof.end(); iter!= end; ++iter){
184 : 731 : Monomial m = (*iter);
185 [ - + ][ - + ]: 731 : Assert(!m.isConstant());
[ - - ]
186 : 731 : VarList vl = m.getVarList();
187 [ - + ][ - + ]: 731 : Assert(vl.singleton());
[ - - ]
188 : 731 : Variable v = vl.getHead();
189 : :
190 : 731 : Node input = proofVariableToReason(v);
191 [ + + ]: 731 : if (input.getKind() == Kind::AND)
192 : : {
193 [ + + ]: 233 : for(Node::iterator input_iter = input.begin(), input_end = input.end(); input_iter != input_end; ++input_iter){
194 : 158 : Node inputChild = *input_iter;
195 : 158 : nb << inputChild;
196 : 158 : }
197 : : }
198 : : else
199 : : {
200 : 656 : nb << input;
201 : : }
202 : 975 : }
203 : :
204 [ - + ]: 244 : Node result = (nb.getNumChildren() == 1) ? nb[0] : (Node)nb;
205 [ + - ]: 488 : Trace("arith::dio") << "Proof at " << i << " is "
206 : 0 : << d_trail[i].d_eq.getNode() << endl
207 : 0 : << d_trail[i].d_proof.getNode() << endl
208 : 244 : << " which becomes " << result << endl;
209 : 488 : return result;
210 : 244 : }
211 : :
212 : 60422 : bool DioSolver::anyCoefficientExceedsMaximum(TrailIndex j) const{
213 : 60422 : uint32_t length = d_trail[j].d_eq.maxLength();
214 : 60422 : uint32_t nmonos = d_trail[j].d_eq.getPolynomial().numMonomials();
215 : :
216 : : bool result =
217 [ + + ]: 98127 : nmonos >= 2 &&
218 [ + + ]: 37705 : length > d_maxInputCoefficientLength + MAX_GROWTH_RATE;
219 : 60422 : if(TraceIsOn("arith::dio::max") && result){
220 : :
221 : 0 : const SumPair& eq = d_trail[j].d_eq;
222 : 0 : const Polynomial& proof = d_trail[j].d_proof;
223 : :
224 [ - - ]: 0 : Trace("arith::dio::max") << "about to drop:" << std::endl;
225 [ - - ]: 0 : Trace("arith::dio::max") << "d_trail[" << j << "].d_eq = " << eq.getNode() << std::endl;
226 [ - - ]: 0 : Trace("arith::dio::max") << "d_trail[" << j << "].d_proof = " << proof.getNode() << std::endl;
227 : : }
228 : 60422 : return result;
229 : : }
230 : :
231 : 6785 : void DioSolver::enqueueInputConstraints(){
232 [ - + ][ - + ]: 6785 : Assert(d_currentF.empty());
[ - - ]
233 [ - + ]: 6785 : while(d_savedQueueIndex < d_savedQueue.size()){
234 : 0 : d_currentF.push_back(d_savedQueue[d_savedQueueIndex]);
235 : 0 : d_savedQueueIndex = d_savedQueueIndex + 1;
236 : : }
237 : :
238 [ + + ][ + + ]: 47795 : while(d_nextInputConstraintToEnqueue < d_inputConstraints.size() && !inConflict()){
[ + + ]
239 : 41010 : size_t curr = d_nextInputConstraintToEnqueue;
240 : 41010 : d_nextInputConstraintToEnqueue = d_nextInputConstraintToEnqueue + 1;
241 : :
242 : 41010 : TrailIndex i = d_inputConstraints[curr].d_trailPos;
243 : 41010 : TrailIndex j = applyAllSubstitutionsToIndex(i);
244 : :
245 [ + + ]: 41010 : if(!triviallySat(j)){
246 [ - + ]: 40104 : if(triviallyUnsat(j)){
247 : 0 : raiseConflict(j);
248 : : }else{
249 : 40104 : TrailIndex k = reduceByGCD(j);
250 : :
251 [ + + ]: 40104 : if(!inConflict()){
252 [ - + ]: 39238 : if(triviallyUnsat(k)){
253 : 0 : raiseConflict(k);
254 [ + - ][ + + ]: 39238 : }else if(!(triviallySat(k) || anyCoefficientExceedsMaximum(k))){
[ + + ]
255 : 39210 : pushToQueueBack(k);
256 : : }
257 : : }
258 : : }
259 : : }
260 : : }
261 : 6785 : }
262 : :
263 : :
264 : : /*TODO Currently linear in the size of the Queue
265 : : *It is not clear if am O(log n) strategy would be better.
266 : : *Right before this in the algorithm is a substitution which could potentially
267 : : *effect the key values of everything in the queue.
268 : : */
269 : 28608 : void DioSolver::moveMinimumByAbsToQueueFront(){
270 [ - + ][ - + ]: 28608 : Assert(!queueEmpty());
[ - - ]
271 : :
272 : : //Select the minimum element.
273 : 28608 : size_t indexInQueue = 0;
274 : 28608 : Monomial minMonomial = d_trail[d_currentF[indexInQueue]].d_minimalMonomial;
275 : :
276 : 28608 : size_t N = d_currentF.size();
277 [ + + ]: 545727 : for(size_t i=1; i < N; ++i){
278 : 517119 : Monomial curr = d_trail[d_currentF[i]].d_minimalMonomial;
279 [ + + ]: 517119 : if(curr.absCmp(minMonomial) < 0){
280 : 1563 : indexInQueue = i;
281 : 1563 : minMonomial = curr;
282 : : }
283 : 517119 : }
284 : :
285 : 28608 : TrailIndex tmp = d_currentF[indexInQueue];
286 : 28608 : d_currentF[indexInQueue] = d_currentF.front();
287 : 28608 : d_currentF.front() = tmp;
288 : 28608 : }
289 : :
290 : 64001 : bool DioSolver::queueEmpty() const{
291 : 64001 : return d_currentF.empty();
292 : : }
293 : :
294 : 1754 : Node DioSolver::columnGcdIsOne() const{
295 : 1754 : std::unordered_map<Node, Integer> gcdMap;
296 : :
297 : 1754 : std::deque<TrailIndex>::const_iterator iter, end;
298 [ + + ]: 5181 : for(iter = d_currentF.begin(), end = d_currentF.end(); iter != end; ++iter){
299 : 3707 : TrailIndex curr = *iter;
300 : 3707 : Polynomial p = d_trail[curr].d_eq.getPolynomial();
301 : 3707 : Polynomial::iterator monoIter = p.begin(), monoEnd = p.end();
302 [ + + ]: 10882 : for(; monoIter != monoEnd; ++monoIter){
303 : 7455 : Monomial m = *monoIter;
304 : 7455 : VarList vl = m.getVarList();
305 : 7455 : Node vlNode = vl.getNode();
306 : :
307 : 7455 : Constant c = m.getConstant();
308 : 7455 : Integer zc = c.getValue().getNumerator();
309 [ + + ]: 7455 : if(gcdMap.find(vlNode) == gcdMap.end()){
310 : : // have not seen vl yet.
311 : : // gcd is c
312 [ - + ][ - + ]: 5392 : Assert(c.isIntegral());
[ - - ]
313 [ - + ][ - + ]: 5392 : Assert(!m.absCoefficientIsOne());
[ - - ]
314 : 5392 : gcdMap.insert(make_pair(vlNode, zc.abs()));
315 : : }else{
316 : 2063 : const Integer& currentGcd = gcdMap[vlNode];
317 : 2063 : Integer newGcd = currentGcd.gcd(zc);
318 [ + + ]: 2063 : if(newGcd == 1){
319 : 280 : return vlNode;
320 : : }else{
321 : 1783 : gcdMap[vlNode] = newGcd;
322 : : }
323 [ + + ]: 2063 : }
324 [ + + ][ + + ]: 8575 : }
[ + + ][ + + ]
[ + + ]
325 [ + + ][ + + ]: 4267 : }
[ + + ]
326 : 1474 : return Node::null();
327 : 1754 : }
328 : :
329 : 0 : void DioSolver::saveQueue(){
330 : 0 : std::deque<TrailIndex>::const_iterator iter, end;
331 [ - - ]: 0 : for(iter = d_currentF.begin(), end = d_currentF.end(); iter != end; ++iter){
332 : 0 : d_savedQueue.push_back(*iter);
333 : : }
334 : 0 : }
335 : :
336 : 1754 : DioSolver::TrailIndex DioSolver::impliedGcdOfOne(){
337 : 1754 : Node canReduce = columnGcdIsOne();
338 [ + + ]: 1754 : if(canReduce.isNull()){
339 : 1474 : return 0;
340 : : }else{
341 : 280 : VarList vl = VarList::parseVarList(canReduce);
342 : :
343 : : TrailIndex current;
344 : 280 : Integer currentCoeff, currentGcd;
345 : :
346 : : //step 1 find the first equation containing vl
347 : : //Set current and currentCoefficient
348 : 280 : std::deque<TrailIndex>::const_iterator iter, end;
349 : 588 : for(iter = d_currentF.begin(), end = d_currentF.end(); true; ++iter){
350 [ - + ][ - + ]: 588 : Assert(iter != end);
[ - - ]
351 : 588 : current = *iter;
352 : 588 : Constant coeff = d_trail[current].d_eq.getPolynomial().getCoefficient(vl);
353 [ + + ]: 588 : if(!coeff.isZero()){
354 : 280 : currentCoeff = coeff.getValue().getNumerator();
355 : 280 : currentGcd = currentCoeff.abs();
356 : :
357 : 280 : ++iter;
358 : 280 : break;
359 : : }
360 [ + + ]: 896 : }
361 : :
362 : : //For the rest of the equations keep reducing until the coefficient is one
363 [ + - ]: 282 : for(; iter != end; ++iter){
364 [ + - ]: 282 : Trace("arith::dio") << "next round : " << currentCoeff << " " << currentGcd << endl;
365 : 282 : TrailIndex inQueue = *iter;
366 : 282 : Constant iqc = d_trail[inQueue].d_eq.getPolynomial().getCoefficient(vl);
367 [ + - ]: 282 : if(!iqc.isZero()){
368 : 282 : Integer inQueueCoeff = iqc.getValue().getNumerator();
369 : :
370 : : //mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, mpz_t a, mpz_t b);
371 : 282 : Integer g, s, t;
372 : : // g = a*s + b*t
373 : 282 : Integer::extendedGcd(g, s, t, currentCoeff, inQueueCoeff);
374 : :
375 [ + - ]: 282 : Trace("arith::dio") << "extendedReduction : " << endl;
376 [ + - ]: 282 : Trace("arith::dio") << g << " = " << s <<"*"<< currentCoeff << " + " << t <<"*"<< inQueueCoeff << endl;
377 : :
378 [ - + ][ - + ]: 282 : Assert(g <= currentGcd);
[ - - ]
379 [ + + ]: 282 : if(g < currentGcd){
380 [ - + ]: 280 : if(s.sgn() == 0){
381 [ - - ]: 0 : Trace("arith::dio") << "extendedReduction drop" << endl;
382 : 0 : Assert(inQueueCoeff.divides(currentGcd));
383 : 0 : current = *iter;
384 : 0 : currentCoeff = inQueueCoeff;
385 : 0 : currentGcd = inQueueCoeff.abs();
386 : : }else{
387 : :
388 [ + - ]: 280 : Trace("arith::dio") << "extendedReduction combine" << endl;
389 : 280 : TrailIndex next = combineEqAtIndexes(current, s, inQueue, t);
390 : :
391 [ - + ][ - + ]: 280 : Assert(d_trail[next]
[ - - ]
392 : : .d_eq.getPolynomial()
393 : : .getCoefficient(vl)
394 : : .getValue()
395 : : .getNumerator()
396 : : == g);
397 : :
398 : 280 : current = next;
399 : 280 : currentCoeff = g;
400 : 280 : currentGcd = g;
401 [ + - ]: 280 : if(currentGcd == 1){
402 : 280 : return current;
403 : : }
404 : : }
405 : : }
406 [ + + ][ + + ]: 1122 : }
[ + + ][ + + ]
407 [ + + ]: 282 : }
408 : : //This is not reachble as it is assured that the gcd of the column is 1
409 : 0 : Unreachable();
410 : 280 : }
411 : 1754 : }
412 : :
413 : 6785 : bool DioSolver::processEquations(bool allowDecomposition){
414 [ - + ][ - + ]: 6785 : Assert(!inConflict());
[ - - ]
415 : :
416 : 6785 : enqueueInputConstraints();
417 [ + + ][ + + ]: 35393 : while(! queueEmpty() && !inConflict()){
[ + + ]
418 : 28608 : moveMinimumByAbsToQueueFront();
419 : :
420 : 28608 : TrailIndex minimum = d_currentF.front();
421 : : TrailIndex reduceIndex;
422 : :
423 [ - + ][ - + ]: 28608 : Assert(inRange(minimum));
[ - - ]
424 [ - + ][ - + ]: 28608 : Assert(!inConflict());
[ - - ]
425 : :
426 [ + - ]: 28608 : Trace("arith::dio") << "processEquations " << minimum << " : " << d_trail[minimum].d_eq.getNode() << endl;
427 : :
428 [ - + ][ - + ]: 28608 : Assert(queueConditions(minimum));
[ - - ]
429 : :
430 : 28608 : bool canDirectlySolve = d_trail[minimum].d_minimalMonomial.absCoefficientIsOne();
431 : :
432 : 28608 : std::pair<SubIndex, TrailIndex> p;
433 [ + + ]: 28608 : if(canDirectlySolve){
434 : 26854 : d_currentF.pop_front();
435 : 26854 : p = solveIndex(minimum);
436 : 26854 : reduceIndex = minimum;
437 : : }else{
438 : 1754 : TrailIndex implied = impliedGcdOfOne();
439 : :
440 [ + + ]: 1754 : if(implied != 0){
441 : 280 : p = solveIndex(implied);
442 : 280 : reduceIndex = implied;
443 [ + - ]: 1474 : }else if(allowDecomposition){
444 : 1474 : d_currentF.pop_front();
445 : 1474 : p = decomposeIndex(minimum);
446 : 1474 : reduceIndex = minimum;
447 : : }else {
448 : : // Cannot make progress without decomposeIndex
449 : 0 : saveQueue();
450 : 0 : break;
451 : : }
452 : : }
453 : :
454 : 28608 : SubIndex subIndex = p.first;
455 : 28608 : TrailIndex next = p.second;
456 : 28608 : subAndReduceCurrentFByIndex(subIndex);
457 : :
458 [ + + ]: 28608 : if(next != reduceIndex){
459 [ - + ]: 1474 : if(triviallyUnsat(next)){
460 : 0 : raiseConflict(next);
461 [ + - ]: 1474 : }else if(! triviallySat(next) ){
462 : 1474 : pushToQueueBack(next);
463 : : }
464 : : }
465 : : }
466 : :
467 : 6785 : d_currentF.clear();
468 : 6785 : return inConflict();
469 : : }
470 : :
471 : 4016 : Node DioSolver::processEquationsForConflict(){
472 : 4016 : TimerStat::CodeTimer codeTimer(d_statistics.d_conflictTimer);
473 : 4016 : ++(d_statistics.d_conflictCalls);
474 : :
475 [ - + ][ - + ]: 4016 : Assert(!inConflict());
[ - - ]
476 [ + + ]: 4016 : if(processEquations(true)){
477 : 244 : ++(d_statistics.d_conflicts);
478 : 244 : return proveIndex(getConflictIndex());
479 : : }else{
480 : 3772 : return Node::null();
481 : : }
482 : 4016 : }
483 : :
484 : 2769 : SumPair DioSolver::processEquationsForCut(){
485 : 2769 : TimerStat::CodeTimer codeTimer(d_statistics.d_cutTimer);
486 : 2769 : ++(d_statistics.d_cutCalls);
487 : :
488 [ - + ][ - + ]: 2769 : Assert(!inConflict());
[ - - ]
489 [ + + ]: 2769 : if(processEquations(true)){
490 : 1815 : ++(d_statistics.d_cuts);
491 : 1815 : return purifyIndex(getConflictIndex());
492 : : }else{
493 : 954 : return SumPair::mkZero(nodeManager());
494 : : }
495 : 2769 : }
496 : :
497 : :
498 : 1815 : SumPair DioSolver::purifyIndex(TrailIndex i){
499 : : // TODO: "This uses the substitution trail to reverse the substitutions from the sum term. Using the proof term should be more efficient."
500 : :
501 : 1815 : SumPair curr = d_trail[i].d_eq;
502 : :
503 : 1815 : Constant negOne = Constant::mkConstant(nodeManager(), -1);
504 : :
505 [ + + ]: 26755 : for(uint32_t revIter = d_subs.size(); revIter > 0; --revIter){
506 : 24940 : uint32_t i2 = revIter - 1;
507 : 24940 : Node freshNode = d_subs[i2].d_fresh;
508 [ + + ]: 24940 : if(freshNode.isNull()){
509 : 23424 : continue;
510 : : }else{
511 : 1516 : Variable var(freshNode);
512 : 1516 : Polynomial vsum = curr.getPolynomial();
513 : :
514 : 3032 : Constant a = vsum.getCoefficient(VarList(var));
515 [ + + ]: 1516 : if(!a.isZero()){
516 : 696 : const SumPair& sj = d_trail[d_subs[i2].d_constraint].d_eq;
517 [ - + ][ - + ]: 696 : Assert(sj.getPolynomial().getCoefficient(VarList(var)).isOne());
[ - - ]
518 : 1392 : SumPair newSi = (curr * negOne) + (sj * a);
519 [ - + ][ - + ]: 696 : Assert(newSi.getPolynomial().getCoefficient(VarList(var)).isZero());
[ - - ]
520 : 696 : curr = newSi;
521 : 696 : }
522 : 1516 : }
523 [ + + ]: 24940 : }
524 : 3630 : return curr;
525 : 1815 : }
526 : :
527 : 44629 : DioSolver::TrailIndex DioSolver::combineEqAtIndexes(DioSolver::TrailIndex i, const Integer& q, DioSolver::TrailIndex j, const Integer& r){
528 : 44629 : Constant cq = Constant::mkConstant(nodeManager(), q);
529 : 44629 : Constant cr = Constant::mkConstant(nodeManager(), r);
530 : :
531 : 44629 : const SumPair& si = d_trail[i].d_eq;
532 : 44629 : const SumPair& sj = d_trail[j].d_eq;
533 : :
534 [ + - ]: 44629 : Trace("arith::dio") << "combineEqAtIndexes(" << i <<","<<q<<","<<j<<","<<r<<")"<<endl;
535 [ + - ]: 89258 : Trace("arith::dio") << "d_facts[i] = " << si.getNode() << endl
536 : 44629 : << "d_facts[j] = " << sj.getNode() << endl;
537 : :
538 : :
539 : 89258 : SumPair newSi = (si * cq) + (sj * cr);
540 : :
541 : :
542 : 44629 : const Polynomial& pi = d_trail[i].d_proof;
543 : 44629 : const Polynomial& pj = d_trail[j].d_proof;
544 : 89258 : Polynomial newPi = (pi * cq) + (pj * cr);
545 : :
546 : 44629 : TrailIndex k = d_trail.size();
547 : 44629 : d_trail.push_back(Constraint(newSi, newPi));
548 : :
549 : :
550 [ + - ]: 89258 : Trace("arith::dio") << "derived "<< newSi.getNode()
551 : 44629 : <<" with proof " << newPi.getNode() << endl;
552 : :
553 : 44629 : return k;
554 : :
555 : 44629 : }
556 : :
557 : 0 : void DioSolver::printQueue(){
558 [ - - ]: 0 : Trace("arith::dio") << "DioSolver::printQueue()" << endl;
559 [ - - ]: 0 : for(TrailIndex i = 0, last = d_trail.size(); i < last; ++i){
560 [ - - ]: 0 : Trace("arith::dio") << "d_trail[i].d_eq = " << d_trail[i].d_eq.getNode() << endl;
561 [ - - ]: 0 : Trace("arith::dio") << "d_trail[i].d_proof = " << d_trail[i].d_proof.getNode() << endl;
562 : : }
563 : :
564 [ - - ]: 0 : Trace("arith::dio") << "DioSolver::printSubs()" << endl;
565 [ - - ]: 0 : for(SubIndex si=0, sN=d_subs.size(); si < sN; ++si){
566 [ - - ]: 0 : Trace("arith::dio") << "d_subs[i] = {"
567 : 0 : << "d_fresh="<< d_subs[si].d_fresh <<","
568 : 0 : << "d_eliminated="<< d_subs[si].d_eliminated.getNode() <<","
569 : 0 : << "d_constraint="<< d_subs[si].d_constraint <<"}" << endl;
570 [ - - ]: 0 : Trace("arith::dio") << "d_trail[d_subs[i].d_constraint].d_eq="
571 : 0 : << d_trail[d_subs[si].d_constraint].d_eq.getNode() << endl;
572 : : }
573 : 0 : }
574 : :
575 : 41010 : DioSolver::TrailIndex DioSolver::applyAllSubstitutionsToIndex(DioSolver::TrailIndex trailIndex){
576 : 41010 : TrailIndex currentIndex = trailIndex;
577 [ + + ]: 401594 : for(SubIndex subIter = 0, siEnd = d_subs.size(); subIter < siEnd; ++subIter){
578 : 360584 : currentIndex = applySubstitution(subIter, currentIndex);
579 : : }
580 : 41010 : return currentIndex;
581 : : }
582 : :
583 : 1299332 : bool DioSolver::debugSubstitutionApplies(DioSolver::SubIndex si, DioSolver::TrailIndex ti){
584 : 1299332 : Variable var = d_subs[si].d_eliminated;
585 : :
586 : 1299332 : const SumPair& curr = d_trail[ti].d_eq;
587 : 1299332 : Polynomial vsum = curr.getPolynomial();
588 : :
589 : 2598664 : Constant a = vsum.getCoefficient(VarList(var));
590 : 2598664 : return !a.isZero();
591 : 1299332 : }
592 : :
593 : 90168 : bool DioSolver::debugAnySubstitionApplies(DioSolver::TrailIndex i){
594 [ + + ]: 1389500 : for(SubIndex subIter = 0, siEnd = d_subs.size(); subIter < siEnd; ++subIter){
595 [ - + ]: 1299332 : if(debugSubstitutionApplies(subIter, i)){
596 : 0 : return true;
597 : : }
598 : : }
599 : 90168 : return false;
600 : : }
601 : :
602 : 27134 : std::pair<DioSolver::SubIndex, DioSolver::TrailIndex> DioSolver::solveIndex(DioSolver::TrailIndex i){
603 : 27134 : const SumPair& si = d_trail[i].d_eq;
604 : :
605 [ + - ]: 27134 : Trace("arith::dio") << "before solveIndex("<<i<<":"<<si.getNode()<< ")" << endl;
606 : :
607 : : #ifdef CVC5_ASSERTIONS
608 : 27134 : const Polynomial& p = si.getPolynomial();
609 : : #endif
610 : :
611 [ - + ][ - + ]: 27134 : Assert(p.isIntegral());
[ - - ]
612 : :
613 [ - + ][ - + ]: 27134 : Assert(p.selectAbsMinimum() == d_trail[i].d_minimalMonomial);
[ - - ]
614 : 27134 : const Monomial av = d_trail[i].d_minimalMonomial;
615 : :
616 : 27134 : VarList vl = av.getVarList();
617 [ - + ][ - + ]: 27134 : Assert(vl.singleton());
[ - - ]
618 : 27134 : Variable var = vl.getHead();
619 : 27134 : Constant a = av.getConstant();
620 : 27134 : Integer a_abs = a.getValue().getNumerator().abs();
621 : :
622 [ - + ][ - + ]: 27134 : Assert(a_abs == 1);
[ - - ]
623 : :
624 [ + + ][ + + ]: 27134 : TrailIndex ci = !a.isNegative() ? scaleEqAtIndex(i, Integer(-1)) : i;
[ - - ]
625 : :
626 : 27134 : SubIndex subBy = d_subs.size();
627 : 27134 : d_subs.push_back(Substitution(Node::null(), var, ci));
628 : :
629 [ + - ]: 27134 : Trace("arith::dio") << "after solveIndex " << d_trail[ci].d_eq.getNode() << " for " << av.getNode() << endl;
630 [ - + ][ - + ]: 27134 : Assert(d_trail[ci].d_eq.getPolynomial().getCoefficient(vl)
[ - - ]
631 : : == Constant::mkConstant(nodeManager(), -1));
632 : :
633 : 54268 : return make_pair(subBy, i);
634 : 27134 : }
635 : :
636 : 1474 : std::pair<DioSolver::SubIndex, DioSolver::TrailIndex> DioSolver::decomposeIndex(DioSolver::TrailIndex i){
637 : 1474 : const SumPair& si = d_trail[i].d_eq;
638 : :
639 : 1474 : d_usedDecomposeIndex = true;
640 : :
641 [ + - ]: 1474 : Trace("arith::dio") << "before decomposeIndex("<<i<<":"<<si.getNode()<< ")" << endl;
642 : :
643 : : #ifdef CVC5_ASSERTIONS
644 : 1474 : const Polynomial& p = si.getPolynomial();
645 : : #endif
646 : :
647 [ - + ][ - + ]: 1474 : Assert(p.isIntegral());
[ - - ]
648 : :
649 [ - + ][ - + ]: 1474 : Assert(p.selectAbsMinimum() == d_trail[i].d_minimalMonomial);
[ - - ]
650 : 1474 : const Monomial& av = d_trail[i].d_minimalMonomial;
651 : :
652 : 1474 : VarList vl = av.getVarList();
653 [ - + ][ - + ]: 1474 : Assert(vl.singleton());
[ - - ]
654 : 1474 : Variable var = vl.getHead();
655 : 1474 : Constant a = av.getConstant();
656 : 1474 : Integer a_abs = a.getValue().getNumerator().abs();
657 : :
658 [ - + ][ - + ]: 1474 : Assert(a_abs > 1);
[ - - ]
659 : :
660 : : //It is not sufficient to reduce the case where abs(a) == 1 to abs(a) > 1.
661 : : //We need to handle both cases seperately to ensure termination.
662 : 1474 : Node qr = SumPair::computeQR(si, a.getValue().getNumerator());
663 : :
664 [ - + ][ - + ]: 1474 : Assert(qr.getKind() == Kind::ADD);
[ - - ]
665 [ - + ][ - + ]: 1474 : Assert(qr.getNumChildren() == 2);
[ - - ]
666 : 2948 : SumPair q = SumPair::parseSumPair(qr[0]);
667 : 2948 : SumPair r = SumPair::parseSumPair(qr[1]);
668 : :
669 : 1474 : NodeManager* nm = nodeManager();
670 [ - + ][ - + ]: 1474 : Assert(q.getPolynomial().getCoefficient(vl) == Constant::mkConstant(nm, 1));
[ - - ]
671 : :
672 [ - + ][ - + ]: 1474 : Assert(!r.isZero());
[ - - ]
673 : 1474 : Node freshNode = makeIntegerVariable(nodeManager());
674 : 1474 : Variable fresh(freshNode);
675 : 1474 : SumPair fresh_one=SumPair::mkSumPair(fresh);
676 : 1474 : SumPair fresh_a = fresh_one * a;
677 : :
678 : 1474 : SumPair newSI = SumPair(fresh_one) - q;
679 : : // this normalizes the coefficient of var to -1
680 : :
681 : :
682 : 1474 : TrailIndex ci = d_trail.size();
683 : 1474 : d_trail.push_back(Constraint(newSI, Polynomial::mkZero(nm)));
684 : : // no longer reference av safely!
685 : 1474 : addTrailElementAsLemma(ci);
686 : :
687 [ + - ]: 2948 : Trace("arith::dio") << "Decompose ci(" << ci <<":" << d_trail[ci].d_eq.getNode()
688 : 1474 : << ") for " << d_trail[i].d_minimalMonomial.getNode() << endl;
689 [ - + ][ - + ]: 1474 : Assert(d_trail[ci].d_eq.getPolynomial().getCoefficient(vl)
[ - - ]
690 : : == Constant::mkConstant(nm, -1));
691 : :
692 : 1474 : SumPair newFact = r + fresh_a;
693 : :
694 : 1474 : TrailIndex nextIndex = d_trail.size();
695 : 1474 : d_trail.push_back(Constraint(newFact, d_trail[i].d_proof));
696 : :
697 : 1474 : SubIndex subBy = d_subs.size();
698 : 1474 : d_subs.push_back(Substitution(freshNode, var, ci));
699 : :
700 [ + - ]: 1474 : Trace("arith::dio") << "Decompose nextIndex " << d_trail[nextIndex].d_eq.getNode() << endl;
701 : 2948 : return make_pair(subBy, nextIndex);
702 : 1474 : }
703 : :
704 : :
705 : 874812 : DioSolver::TrailIndex DioSolver::applySubstitution(DioSolver::SubIndex si, DioSolver::TrailIndex ti){
706 : 874812 : Variable var = d_subs[si].d_eliminated;
707 : 874812 : TrailIndex subIndex = d_subs[si].d_constraint;
708 : :
709 : 874812 : const SumPair& curr = d_trail[ti].d_eq;
710 : 874812 : Polynomial vsum = curr.getPolynomial();
711 : :
712 : 1749624 : Constant a = vsum.getCoefficient(VarList(var));
713 [ - + ][ - + ]: 874812 : Assert(a.isIntegral());
[ - - ]
714 [ + + ]: 874812 : if(!a.isZero()){
715 : 44349 : Integer one(1);
716 : 44349 : TrailIndex afterSub = combineEqAtIndexes(ti, one, subIndex, a.getValue().getNumerator());
717 [ - + ][ - + ]: 44349 : Assert(d_trail[afterSub]
[ - - ]
718 : : .d_eq.getPolynomial()
719 : : .getCoefficient(VarList(var))
720 : : .isZero());
721 : 44349 : return afterSub;
722 : 44349 : }else{
723 : 830463 : return ti;
724 : : }
725 : 874812 : }
726 : :
727 : :
728 : 62481 : DioSolver::TrailIndex DioSolver::reduceByGCD(DioSolver::TrailIndex ti){
729 : 62481 : const SumPair& sp = d_trail[ti].d_eq;
730 : 62481 : Polynomial vsum = sp.getPolynomial();
731 : 62481 : Constant c = sp.getConstant();
732 : :
733 [ + - ]: 62481 : Trace("arith::dio") << "reduceByGCD " << vsum.getNode() << endl;
734 [ - + ][ - + ]: 62481 : Assert(!vsum.isConstant());
[ - - ]
735 : 62481 : Integer g = vsum.gcd();
736 [ - + ][ - + ]: 62481 : Assert(g >= 1);
[ - - ]
737 [ + - ]: 62481 : Trace("arith::dio") << "gcd("<< vsum.getNode() <<")=" << g << " " << c.getValue() << endl;
738 [ + + ]: 62481 : if(g.divides(c.getValue().getNumerator())){
739 [ + + ]: 60422 : if(g > 1){
740 : 2296 : return scaleEqAtIndex(ti, g);
741 : : }else{
742 : 58126 : return ti;
743 : : }
744 : : }else{
745 : 2059 : raiseConflict(ti);
746 : 2059 : return ti;
747 : : }
748 : 62481 : }
749 : :
750 : 198156 : bool DioSolver::triviallySat(TrailIndex i){
751 : 198156 : const SumPair& eq = d_trail[i].d_eq;
752 [ + + ]: 198156 : if(eq.isConstant()){
753 : 4795 : return eq.getConstant().isZero();
754 : : }else{
755 : 193361 : return false;
756 : : }
757 : : }
758 : :
759 : 197250 : bool DioSolver::triviallyUnsat(DioSolver::TrailIndex i){
760 : 197250 : const SumPair& eq = d_trail[i].d_eq;
761 [ + + ]: 197250 : if(eq.isConstant()){
762 : 3889 : return !eq.getConstant().isZero();
763 : : }else{
764 : 193361 : return false;
765 : : }
766 : : }
767 : :
768 : :
769 : 90168 : bool DioSolver::gcdIsOne(DioSolver::TrailIndex i){
770 : 90168 : const SumPair& eq = d_trail[i].d_eq;
771 : 90168 : return eq.gcd() == Integer(1);
772 : : }
773 : :
774 : 28608 : void DioSolver::subAndReduceCurrentFByIndex(DioSolver::SubIndex subIndex){
775 : 28608 : size_t N = d_currentF.size();
776 : :
777 : 28608 : size_t readIter = 0, writeIter = 0;
778 [ + + ][ + + ]: 542836 : for(; readIter < N && !inConflict(); ++readIter){
[ + + ]
779 : 514228 : TrailIndex curr = d_currentF[readIter];
780 : 514228 : TrailIndex nextTI = applySubstitution(subIndex, curr);
781 [ + + ]: 514228 : if(nextTI == curr){
782 : 487962 : d_currentF[writeIter] = curr;
783 : 487962 : ++writeIter;
784 : : }else{
785 [ - + ][ - + ]: 26266 : Assert(nextTI != curr);
[ - - ]
786 : :
787 [ - + ]: 26266 : if(triviallyUnsat(nextTI)){
788 : 0 : raiseConflict(nextTI);
789 [ + + ]: 26266 : }else if(!triviallySat(nextTI)){
790 : 22377 : TrailIndex nextNextTI = reduceByGCD(nextTI);
791 : :
792 [ + + ][ + + ]: 22377 : if(!(inConflict() || anyCoefficientExceedsMaximum(nextNextTI))){
[ + + ]
793 [ - + ][ - + ]: 20876 : Assert(queueConditions(nextNextTI));
[ - - ]
794 : 20876 : d_currentF[writeIter] = nextNextTI;
795 : 20876 : ++writeIter;
796 : : }
797 : : }
798 : : }
799 : : }
800 [ + + ][ + + ]: 28608 : if(!inConflict() && writeIter < N){
[ + + ]
801 : 2800 : d_currentF.resize(writeIter);
802 : : }
803 : 28608 : }
804 : :
805 : 1474 : void DioSolver::addTrailElementAsLemma(TrailIndex i) {
806 [ - + ]: 1474 : if (options().arith.exportDioDecompositions)
807 : : {
808 : 0 : d_decompositionLemmaQueue.push(i);
809 : : }
810 : 1474 : }
811 : :
812 : 0 : Node DioSolver::trailIndexToEquality(TrailIndex i) const {
813 : 0 : const SumPair& sp = d_trail[i].d_eq;
814 : 0 : Node n = sp.getNode();
815 : 0 : Node zero = nodeManager()->mkConstRealOrInt(n.getType(), Rational(0));
816 : 0 : Node eq = n.eqNode(zero);
817 : 0 : return eq;
818 : 0 : }
819 : :
820 : : } // namespace arith
821 : : } // namespace theory
822 : : } // namespace cvc5::internal
|
