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