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 : 12866 : inline Node makeIntegerVariable(NodeManager* nm)
31 : : {
32 : 12866 : return NodeManager::mkDummySkolem("intvar", nm->integerType());
33 : : }
34 : :
35 : 51213 : DioSolver::DioSolver(Env& env)
36 : : : EnvObj(env),
37 : 51213 : d_lastUsedProofVariable(context(), 0),
38 : 51213 : d_inputConstraints(context()),
39 : 51213 : d_nextInputConstraintToEnqueue(context(), 0),
40 : 51213 : d_trail(context()),
41 : 51213 : d_subs(context()),
42 : 51213 : d_currentF(),
43 : 51213 : d_savedQueue(context()),
44 : 51213 : d_savedQueueIndex(context(), 0),
45 : 51213 : d_conflictIndex(context()),
46 : 51213 : d_maxInputCoefficientLength(context(), 0),
47 : 51213 : d_usedDecomposeIndex(context(), false),
48 : 51213 : d_lastPureSubstitution(context(), 0),
49 : 51213 : d_pureSubstitionIter(context(), 0),
50 : 51213 : d_decompositionLemmaQueue(context()),
51 : 153639 : d_statistics(statisticsRegistry())
52 : : {
53 : 51213 : }
54 : :
55 : 51213 : DioSolver::Statistics::Statistics(StatisticsRegistry& sr)
56 : 51213 : : d_conflictCalls(sr.registerInt("theory::arith::dio::conflictCalls")),
57 : 51213 : d_cutCalls(sr.registerInt("theory::arith::dio::cutCalls")),
58 : 51213 : d_cuts(sr.registerInt("theory::arith::dio::cuts")),
59 : 51213 : d_conflicts(sr.registerInt("theory::arith::dio::conflicts")),
60 : 51213 : d_conflictTimer(sr.registerTimer("theory::arith::dio::conflictTimer")),
61 : 51213 : d_cutTimer(sr.registerTimer("theory::arith::dio::cutTimer"))
62 : : {
63 : 51213 : }
64 : :
65 : 89608 : bool DioSolver::queueConditions(TrailIndex t)
66 : : {
67 [ + - ]: 89608 : Trace("queueConditions") << !inConflict() << std::endl;
68 [ + - ]: 89608 : Trace("queueConditions") << gcdIsOne(t) << std::endl;
69 [ + - ]: 89608 : Trace("queueConditions") << !debugAnySubstitionApplies(t) << std::endl;
70 [ + - ]: 89608 : Trace("queueConditions") << !triviallySat(t) << std::endl;
71 [ + - ]: 89608 : Trace("queueConditions") << !triviallyUnsat(t) << std::endl;
72 : :
73 [ + - ][ + - ]: 179216 : return !inConflict() && gcdIsOne(t) && !debugAnySubstitionApplies(t)
74 [ + - ][ + - ]: 179216 : && !triviallySat(t) && !triviallyUnsat(t);
[ + - ]
75 : : }
76 : :
77 : 41978 : size_t DioSolver::allocateProofVariable()
78 : : {
79 [ - + ][ - + ]: 41978 : Assert(d_lastUsedProofVariable <= d_proofVariablePool.size());
[ - - ]
80 [ + + ]: 41978 : if (d_lastUsedProofVariable == d_proofVariablePool.size())
81 : : {
82 [ - + ][ - + ]: 11392 : Assert(d_lastUsedProofVariable == d_proofVariablePool.size());
[ - - ]
83 : 11392 : Node intVar = makeIntegerVariable(nodeManager());
84 : 11392 : d_proofVariablePool.push_back(Variable(intVar));
85 : 11392 : }
86 : 41978 : size_t res = d_lastUsedProofVariable;
87 : 41978 : d_lastUsedProofVariable = d_lastUsedProofVariable + 1;
88 : 41978 : return res;
89 : : }
90 : :
91 : 0 : Node DioSolver::nextPureSubstitution()
92 : : {
93 : 0 : Assert(hasMorePureSubstitutions());
94 : 0 : SubIndex curr = d_pureSubstitionIter;
95 : 0 : d_pureSubstitionIter = d_pureSubstitionIter + 1;
96 : :
97 : 0 : Assert(d_subs[curr].d_fresh.isNull());
98 : 0 : Variable v = d_subs[curr].d_eliminated;
99 : :
100 : 0 : SumPair sp = d_trail[d_subs[curr].d_constraint].d_eq;
101 : 0 : Polynomial p = sp.getPolynomial();
102 : 0 : Constant c = -sp.getConstant();
103 : 0 : Polynomial cancelV = p + Polynomial::mkPolynomial(v);
104 : 0 : Node eq = nodeManager()->mkNode(Kind::EQUAL, v.getNode(), cancelV.getNode());
105 : 0 : return eq;
106 : 0 : }
107 : :
108 : 41978 : bool DioSolver::debugEqualityInInputEquations(Node eq)
109 : : {
110 : : typedef context::CDList<InputConstraint>::const_iterator const_iterator;
111 : 41978 : const_iterator i = d_inputConstraints.begin(), end = d_inputConstraints.end();
112 [ + + ]: 1027991 : for (; i != end; ++i)
113 : : {
114 : 986013 : Node reason_i = (*i).d_reason;
115 [ - + ]: 986013 : if (eq == reason_i)
116 : : {
117 : 0 : return true;
118 : : }
119 [ + - ]: 986013 : }
120 : 41978 : return false;
121 : : }
122 : :
123 : 41978 : void DioSolver::pushInputConstraint(const Comparison& eq, Node reason)
124 : : {
125 [ - + ][ - + ]: 41978 : Assert(!debugEqualityInInputEquations(reason));
[ - - ]
126 [ - + ][ - + ]: 41978 : Assert(eq.debugIsIntegral());
[ - - ]
127 [ - + ][ - + ]: 41978 : Assert(eq.getNode().getKind() == Kind::EQUAL);
[ - - ]
128 : :
129 : 41978 : SumPair sp = eq.toSumPair();
130 [ - + ][ - + ]: 41978 : Assert(!sp.isNonlinear());
[ - - ]
131 : :
132 : 41978 : uint32_t length = sp.maxLength();
133 [ + + ]: 41978 : if (length > d_maxInputCoefficientLength)
134 : : {
135 : 4973 : d_maxInputCoefficientLength = length;
136 : : }
137 : :
138 : 41978 : size_t varIndex = allocateProofVariable();
139 : 41978 : Variable proofVariable(d_proofVariablePool[varIndex]);
140 : : // Variable proofVariable(makeIntegerVariable(nodeManager()));
141 : :
142 : 41978 : TrailIndex posInTrail = d_trail.size();
143 [ + - ]: 83956 : Trace("dio::pushInputConstraint")
144 : 0 : << "pushInputConstraint @ " << posInTrail << " " << eq.getNode() << " "
145 : 41978 : << reason << endl;
146 : 41978 : d_trail.push_back(Constraint(sp, Polynomial::mkPolynomial(proofVariable)));
147 : :
148 : 41978 : size_t posInConstraintList = d_inputConstraints.size();
149 : 41978 : d_inputConstraints.push_back(InputConstraint(reason, posInTrail));
150 : :
151 : 41978 : d_varToInputConstraintMap[proofVariable.getNode()] = posInConstraintList;
152 : 41978 : }
153 : :
154 : 26350 : DioSolver::TrailIndex DioSolver::scaleEqAtIndex(DioSolver::TrailIndex i,
155 : : const Integer& g)
156 : : {
157 [ - + ][ - + ]: 26350 : Assert(g != 0);
[ - - ]
158 : 52700 : Constant invg = Constant::mkConstant(nodeManager(), Rational(Integer(1), g));
159 : 26350 : const SumPair& sp = d_trail[i].d_eq;
160 : 26350 : const Polynomial& proof = d_trail[i].d_proof;
161 : :
162 : 26350 : SumPair newSP = sp * invg;
163 : 26350 : Polynomial newProof = proof * invg;
164 : :
165 [ - + ][ - + ]: 26350 : Assert(newSP.isIntegral());
[ - - ]
166 [ - + ][ - + ]: 26350 : Assert(newSP.gcd() == 1);
[ - - ]
167 : :
168 : 26350 : TrailIndex j = d_trail.size();
169 : :
170 : 26350 : d_trail.push_back(Constraint(newSP, newProof));
171 : :
172 [ + - ]: 26350 : Trace("arith::dio") << "scaleEqAtIndex(" << i << "," << g << ")" << endl;
173 [ + - ]: 52700 : Trace("arith::dio") << "derived " << newSP.getNode() << " with proof "
174 : 26350 : << newProof.getNode() << endl;
175 : 26350 : return j;
176 : 26350 : }
177 : :
178 : 255 : Node DioSolver::proveIndex(TrailIndex i)
179 : : {
180 [ - + ][ - + ]: 255 : Assert(inRange(i));
[ - - ]
181 : 255 : const Polynomial& proof = d_trail[i].d_proof;
182 [ - + ][ - + ]: 255 : Assert(!proof.isConstant());
[ - - ]
183 : :
184 : 255 : NodeBuilder nb(nodeManager(), Kind::AND);
185 : 1050 : for (Polynomial::iterator iter = proof.begin(), end = proof.end();
186 [ + + ]: 1050 : iter != end;
187 : : ++iter)
188 : : {
189 : 795 : Monomial m = (*iter);
190 [ - + ][ - + ]: 795 : Assert(!m.isConstant());
[ - - ]
191 : 795 : VarList vl = m.getVarList();
192 [ - + ][ - + ]: 795 : Assert(vl.singleton());
[ - - ]
193 : 795 : Variable v = vl.getHead();
194 : :
195 : 795 : Node input = proofVariableToReason(v);
196 [ + + ]: 795 : if (input.getKind() == Kind::AND)
197 : : {
198 : 79 : for (Node::iterator input_iter = input.begin(), input_end = input.end();
199 [ + + ]: 245 : input_iter != input_end;
200 : 166 : ++input_iter)
201 : : {
202 : 166 : Node inputChild = *input_iter;
203 : 166 : nb << inputChild;
204 : 166 : }
205 : : }
206 : : else
207 : : {
208 : 716 : nb << input;
209 : : }
210 : 1050 : }
211 : :
212 [ - + ]: 255 : Node result = (nb.getNumChildren() == 1) ? nb[0] : (Node)nb;
213 [ + - ]: 510 : Trace("arith::dio") << "Proof at " << i << " is " << d_trail[i].d_eq.getNode()
214 : 0 : << endl
215 : 0 : << d_trail[i].d_proof.getNode() << endl
216 : 255 : << " which becomes " << result << endl;
217 : 510 : return result;
218 : 255 : }
219 : :
220 : 60022 : bool DioSolver::anyCoefficientExceedsMaximum(TrailIndex j) const
221 : : {
222 : 60022 : uint32_t length = d_trail[j].d_eq.maxLength();
223 : 60022 : uint32_t nmonos = d_trail[j].d_eq.getPolynomial().numMonomials();
224 : :
225 : : bool result =
226 [ + + ][ + + ]: 60022 : nmonos >= 2 && length > d_maxInputCoefficientLength + MAX_GROWTH_RATE;
227 : 60022 : if (TraceIsOn("arith::dio::max") && result)
228 : : {
229 : 0 : const SumPair& eq = d_trail[j].d_eq;
230 : 0 : const Polynomial& proof = d_trail[j].d_proof;
231 : :
232 [ - - ]: 0 : Trace("arith::dio::max") << "about to drop:" << std::endl;
233 [ - - ]: 0 : Trace("arith::dio::max")
234 : 0 : << "d_trail[" << j << "].d_eq = " << eq.getNode() << std::endl;
235 [ - - ]: 0 : Trace("arith::dio::max")
236 : 0 : << "d_trail[" << j << "].d_proof = " << proof.getNode() << std::endl;
237 : : }
238 : 60022 : return result;
239 : : }
240 : :
241 : 7084 : void DioSolver::enqueueInputConstraints()
242 : : {
243 [ - + ][ - + ]: 7084 : Assert(d_currentF.empty());
[ - - ]
244 [ - + ]: 7084 : while (d_savedQueueIndex < d_savedQueue.size())
245 : : {
246 : 0 : d_currentF.push_back(d_savedQueue[d_savedQueueIndex]);
247 : 0 : d_savedQueueIndex = d_savedQueueIndex + 1;
248 : : }
249 : :
250 : 47631 : while (d_nextInputConstraintToEnqueue < d_inputConstraints.size()
251 [ + + ][ + + ]: 47631 : && !inConflict())
[ + + ]
252 : : {
253 : 40547 : size_t curr = d_nextInputConstraintToEnqueue;
254 : 40547 : d_nextInputConstraintToEnqueue = d_nextInputConstraintToEnqueue + 1;
255 : :
256 : 40547 : TrailIndex i = d_inputConstraints[curr].d_trailPos;
257 : 40547 : TrailIndex j = applyAllSubstitutionsToIndex(i);
258 : :
259 [ + + ]: 40547 : if (!triviallySat(j))
260 : : {
261 [ - + ]: 39616 : if (triviallyUnsat(j))
262 : : {
263 : 0 : raiseConflict(j);
264 : : }
265 : : else
266 : : {
267 : 39616 : TrailIndex k = reduceByGCD(j);
268 : :
269 [ + + ]: 39616 : if (!inConflict())
270 : : {
271 [ - + ]: 38726 : if (triviallyUnsat(k))
272 : : {
273 : 0 : raiseConflict(k);
274 : : }
275 [ + - ][ + + ]: 38726 : else if (!(triviallySat(k) || anyCoefficientExceedsMaximum(k)))
[ + + ]
276 : : {
277 : 38698 : pushToQueueBack(k);
278 : : }
279 : : }
280 : : }
281 : : }
282 : : }
283 : 7084 : }
284 : :
285 : : /*TODO Currently linear in the size of the Queue
286 : : *It is not clear if am O(log n) strategy would be better.
287 : : *Right before this in the algorithm is a substitution which could potentially
288 : : *effect the key values of everything in the queue.
289 : : */
290 : 28458 : void DioSolver::moveMinimumByAbsToQueueFront()
291 : : {
292 [ - + ][ - + ]: 28458 : Assert(!queueEmpty());
[ - - ]
293 : :
294 : : // Select the minimum element.
295 : 28458 : size_t indexInQueue = 0;
296 : 28458 : Monomial minMonomial = d_trail[d_currentF[indexInQueue]].d_minimalMonomial;
297 : :
298 : 28458 : size_t N = d_currentF.size();
299 [ + + ]: 470630 : for (size_t i = 1; i < N; ++i)
300 : : {
301 : 442172 : Monomial curr = d_trail[d_currentF[i]].d_minimalMonomial;
302 [ + + ]: 442172 : if (curr.absCmp(minMonomial) < 0)
303 : : {
304 : 1631 : indexInQueue = i;
305 : 1631 : minMonomial = curr;
306 : : }
307 : 442172 : }
308 : :
309 : 28458 : TrailIndex tmp = d_currentF[indexInQueue];
310 : 28458 : d_currentF[indexInQueue] = d_currentF.front();
311 : 28458 : d_currentF.front() = tmp;
312 : 28458 : }
313 : :
314 : 64000 : bool DioSolver::queueEmpty() const { return d_currentF.empty(); }
315 : :
316 : 1754 : Node DioSolver::columnGcdIsOne() const
317 : : {
318 : 1754 : std::unordered_map<Node, Integer> gcdMap;
319 : :
320 : 1754 : std::deque<TrailIndex>::const_iterator iter, end;
321 [ + + ]: 5181 : for (iter = d_currentF.begin(), end = d_currentF.end(); iter != end; ++iter)
322 : : {
323 : 3707 : TrailIndex curr = *iter;
324 : 3707 : Polynomial p = d_trail[curr].d_eq.getPolynomial();
325 : 3707 : Polynomial::iterator monoIter = p.begin(), monoEnd = p.end();
326 [ + + ]: 10882 : for (; monoIter != monoEnd; ++monoIter)
327 : : {
328 : 7455 : Monomial m = *monoIter;
329 : 7455 : VarList vl = m.getVarList();
330 : 7455 : Node vlNode = vl.getNode();
331 : :
332 : 7455 : Constant c = m.getConstant();
333 : 7455 : Integer zc = c.getValue().getNumerator();
334 [ + + ]: 7455 : if (gcdMap.find(vlNode) == gcdMap.end())
335 : : {
336 : : // have not seen vl yet.
337 : : // gcd is c
338 [ - + ][ - + ]: 5392 : Assert(c.isIntegral());
[ - - ]
339 [ - + ][ - + ]: 5392 : Assert(!m.absCoefficientIsOne());
[ - - ]
340 : 5392 : gcdMap.insert(make_pair(vlNode, zc.abs()));
341 : : }
342 : : else
343 : : {
344 : 2063 : const Integer& currentGcd = gcdMap[vlNode];
345 : 2063 : Integer newGcd = currentGcd.gcd(zc);
346 [ + + ]: 2063 : if (newGcd == 1)
347 : : {
348 : 280 : return vlNode;
349 : : }
350 : : else
351 : : {
352 : 1783 : gcdMap[vlNode] = newGcd;
353 : : }
354 [ + + ]: 2063 : }
355 [ + + ][ + + ]: 8575 : }
[ + + ][ + + ]
[ + + ]
356 [ + + ][ + + ]: 4267 : }
[ + + ]
357 : 1474 : return Node::null();
358 : 1754 : }
359 : :
360 : 0 : void DioSolver::saveQueue()
361 : : {
362 : 0 : std::deque<TrailIndex>::const_iterator iter, end;
363 [ - - ]: 0 : for (iter = d_currentF.begin(), end = d_currentF.end(); iter != end; ++iter)
364 : : {
365 : 0 : d_savedQueue.push_back(*iter);
366 : : }
367 : 0 : }
368 : :
369 : 1754 : DioSolver::TrailIndex DioSolver::impliedGcdOfOne()
370 : : {
371 : 1754 : Node canReduce = columnGcdIsOne();
372 [ + + ]: 1754 : if (canReduce.isNull())
373 : : {
374 : 1474 : return 0;
375 : : }
376 : : else
377 : : {
378 : 280 : VarList vl = VarList::parseVarList(canReduce);
379 : :
380 : : TrailIndex current;
381 : 280 : Integer currentCoeff, currentGcd;
382 : :
383 : : // step 1 find the first equation containing vl
384 : : // Set current and currentCoefficient
385 : 280 : std::deque<TrailIndex>::const_iterator iter, end;
386 : 588 : for (iter = d_currentF.begin(), end = d_currentF.end(); true; ++iter)
387 : : {
388 [ - + ][ - + ]: 588 : Assert(iter != end);
[ - - ]
389 : 588 : current = *iter;
390 : 588 : Constant coeff = d_trail[current].d_eq.getPolynomial().getCoefficient(vl);
391 [ + + ]: 588 : if (!coeff.isZero())
392 : : {
393 : 280 : currentCoeff = coeff.getValue().getNumerator();
394 : 280 : currentGcd = currentCoeff.abs();
395 : :
396 : 280 : ++iter;
397 : 280 : break;
398 : : }
399 [ + + ]: 896 : }
400 : :
401 : : // For the rest of the equations keep reducing until the coefficient is one
402 [ + - ]: 282 : for (; iter != end; ++iter)
403 : : {
404 [ + - ]: 564 : Trace("arith::dio") << "next round : " << currentCoeff << " "
405 : 282 : << currentGcd << endl;
406 : 282 : TrailIndex inQueue = *iter;
407 : 282 : Constant iqc = d_trail[inQueue].d_eq.getPolynomial().getCoefficient(vl);
408 [ + - ]: 282 : if (!iqc.isZero())
409 : : {
410 : 282 : Integer inQueueCoeff = iqc.getValue().getNumerator();
411 : :
412 : : // mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, mpz_t a, mpz_t b);
413 : 282 : Integer g, s, t;
414 : : // g = a*s + b*t
415 : 282 : Integer::extendedGcd(g, s, t, currentCoeff, inQueueCoeff);
416 : :
417 [ + - ]: 282 : Trace("arith::dio") << "extendedReduction : " << endl;
418 [ + - ]: 564 : Trace("arith::dio") << g << " = " << s << "*" << currentCoeff << " + "
419 : 282 : << t << "*" << inQueueCoeff << endl;
420 : :
421 [ - + ][ - + ]: 282 : Assert(g <= currentGcd);
[ - - ]
422 [ + + ]: 282 : if (g < currentGcd)
423 : : {
424 [ - + ]: 280 : if (s.sgn() == 0)
425 : : {
426 [ - - ]: 0 : Trace("arith::dio") << "extendedReduction drop" << endl;
427 : 0 : Assert(inQueueCoeff.divides(currentGcd));
428 : 0 : current = *iter;
429 : 0 : currentCoeff = inQueueCoeff;
430 : 0 : currentGcd = inQueueCoeff.abs();
431 : : }
432 : : else
433 : : {
434 [ + - ]: 280 : Trace("arith::dio") << "extendedReduction combine" << endl;
435 : 280 : TrailIndex next = combineEqAtIndexes(current, s, inQueue, t);
436 : :
437 [ - + ][ - + ]: 280 : Assert(d_trail[next]
[ - - ]
438 : : .d_eq.getPolynomial()
439 : : .getCoefficient(vl)
440 : : .getValue()
441 : : .getNumerator()
442 : : == g);
443 : :
444 : 280 : current = next;
445 : 280 : currentCoeff = g;
446 : 280 : currentGcd = g;
447 [ + - ]: 280 : if (currentGcd == 1)
448 : : {
449 : 280 : return current;
450 : : }
451 : : }
452 : : }
453 [ + + ][ + + ]: 1122 : }
[ + + ][ + + ]
454 [ + + ]: 282 : }
455 : : // This is not reachble as it is assured that the gcd of the column is 1
456 : 0 : Unreachable();
457 : 280 : }
458 : 1754 : }
459 : :
460 : 7084 : bool DioSolver::processEquations(bool allowDecomposition)
461 : : {
462 [ - + ][ - + ]: 7084 : Assert(!inConflict());
[ - - ]
463 : :
464 : 7084 : enqueueInputConstraints();
465 [ + + ][ + + ]: 35542 : while (!queueEmpty() && !inConflict())
[ + + ]
466 : : {
467 : 28458 : moveMinimumByAbsToQueueFront();
468 : :
469 : 28458 : TrailIndex minimum = d_currentF.front();
470 : : TrailIndex reduceIndex;
471 : :
472 [ - + ][ - + ]: 28458 : Assert(inRange(minimum));
[ - - ]
473 [ - + ][ - + ]: 28458 : Assert(!inConflict());
[ - - ]
474 : :
475 [ + - ]: 56916 : Trace("arith::dio") << "processEquations " << minimum << " : "
476 : 28458 : << d_trail[minimum].d_eq.getNode() << endl;
477 : :
478 [ - + ][ - + ]: 28458 : Assert(queueConditions(minimum));
[ - - ]
479 : :
480 : : bool canDirectlySolve =
481 : 28458 : d_trail[minimum].d_minimalMonomial.absCoefficientIsOne();
482 : :
483 : 28458 : std::pair<SubIndex, TrailIndex> p;
484 [ + + ]: 28458 : if (canDirectlySolve)
485 : : {
486 : 26704 : d_currentF.pop_front();
487 : 26704 : p = solveIndex(minimum);
488 : 26704 : reduceIndex = minimum;
489 : : }
490 : : else
491 : : {
492 : 1754 : TrailIndex implied = impliedGcdOfOne();
493 : :
494 [ + + ]: 1754 : if (implied != 0)
495 : : {
496 : 280 : p = solveIndex(implied);
497 : 280 : reduceIndex = implied;
498 : : }
499 [ + - ]: 1474 : else if (allowDecomposition)
500 : : {
501 : 1474 : d_currentF.pop_front();
502 : 1474 : p = decomposeIndex(minimum);
503 : 1474 : reduceIndex = minimum;
504 : : }
505 : : else
506 : : {
507 : : // Cannot make progress without decomposeIndex
508 : 0 : saveQueue();
509 : 0 : break;
510 : : }
511 : : }
512 : :
513 : 28458 : SubIndex subIndex = p.first;
514 : 28458 : TrailIndex next = p.second;
515 : 28458 : subAndReduceCurrentFByIndex(subIndex);
516 : :
517 [ + + ]: 28458 : if (next != reduceIndex)
518 : : {
519 [ - + ]: 1474 : if (triviallyUnsat(next))
520 : : {
521 : 0 : raiseConflict(next);
522 : : }
523 [ + - ]: 1474 : else if (!triviallySat(next))
524 : : {
525 : 1474 : pushToQueueBack(next);
526 : : }
527 : : }
528 : : }
529 : :
530 : 7084 : d_currentF.clear();
531 : 7084 : return inConflict();
532 : : }
533 : :
534 : 4181 : Node DioSolver::processEquationsForConflict()
535 : : {
536 : 4181 : TimerStat::CodeTimer codeTimer(d_statistics.d_conflictTimer);
537 : 4181 : ++(d_statistics.d_conflictCalls);
538 : :
539 [ - + ][ - + ]: 4181 : Assert(!inConflict());
[ - - ]
540 [ + + ]: 4181 : if (processEquations(true))
541 : : {
542 : 255 : ++(d_statistics.d_conflicts);
543 : 255 : return proveIndex(getConflictIndex());
544 : : }
545 : : else
546 : : {
547 : 3926 : return Node::null();
548 : : }
549 : 4181 : }
550 : :
551 : 2903 : SumPair DioSolver::processEquationsForCut()
552 : : {
553 : 2903 : TimerStat::CodeTimer codeTimer(d_statistics.d_cutTimer);
554 : 2903 : ++(d_statistics.d_cutCalls);
555 : :
556 [ - + ][ - + ]: 2903 : Assert(!inConflict());
[ - - ]
557 [ + + ]: 2903 : if (processEquations(true))
558 : : {
559 : 1944 : ++(d_statistics.d_cuts);
560 : 1944 : return purifyIndex(getConflictIndex());
561 : : }
562 : : else
563 : : {
564 : 959 : return SumPair::mkZero(nodeManager());
565 : : }
566 : 2903 : }
567 : :
568 : 1944 : SumPair DioSolver::purifyIndex(TrailIndex i)
569 : : {
570 : : // TODO: "This uses the substitution trail to reverse the substitutions from
571 : : // the sum term. Using the proof term should be more efficient."
572 : :
573 : 1944 : SumPair curr = d_trail[i].d_eq;
574 : :
575 : 1944 : Constant negOne = Constant::mkConstant(nodeManager(), -1);
576 : :
577 [ + + ]: 26700 : for (uint32_t revIter = d_subs.size(); revIter > 0; --revIter)
578 : : {
579 : 24756 : uint32_t i2 = revIter - 1;
580 : 24756 : Node freshNode = d_subs[i2].d_fresh;
581 [ + + ]: 24756 : if (freshNode.isNull())
582 : : {
583 : 23240 : continue;
584 : : }
585 : : else
586 : : {
587 : 1516 : Variable var(freshNode);
588 : 1516 : Polynomial vsum = curr.getPolynomial();
589 : :
590 : 3032 : Constant a = vsum.getCoefficient(VarList(var));
591 [ + + ]: 1516 : if (!a.isZero())
592 : : {
593 : 696 : const SumPair& sj = d_trail[d_subs[i2].d_constraint].d_eq;
594 [ - + ][ - + ]: 696 : Assert(sj.getPolynomial().getCoefficient(VarList(var)).isOne());
[ - - ]
595 : : // Use currTimesNegOne and sjTimesA to ensure deterministic node ID
596 : : // assignments
597 : 696 : SumPair currTimesNegOne = curr * negOne;
598 : 696 : SumPair sjTimesA = sj * a;
599 : 696 : SumPair newSi = currTimesNegOne + sjTimesA;
600 [ - + ][ - + ]: 696 : Assert(newSi.getPolynomial().getCoefficient(VarList(var)).isZero());
[ - - ]
601 : 696 : curr = newSi;
602 : 696 : }
603 : 1516 : }
604 [ + + ]: 24756 : }
605 : 3888 : return curr;
606 : 1944 : }
607 : :
608 : 44169 : DioSolver::TrailIndex DioSolver::combineEqAtIndexes(DioSolver::TrailIndex i,
609 : : const Integer& q,
610 : : DioSolver::TrailIndex j,
611 : : const Integer& r)
612 : : {
613 : 44169 : Constant cq = Constant::mkConstant(nodeManager(), q);
614 : 44169 : Constant cr = Constant::mkConstant(nodeManager(), r);
615 : :
616 : 44169 : const SumPair& si = d_trail[i].d_eq;
617 : 44169 : const SumPair& sj = d_trail[j].d_eq;
618 : :
619 [ + - ]: 88338 : Trace("arith::dio") << "combineEqAtIndexes(" << i << "," << q << "," << j
620 : 44169 : << "," << r << ")" << endl;
621 [ + - ]: 88338 : Trace("arith::dio") << "d_facts[i] = " << si.getNode() << endl
622 : 44169 : << "d_facts[j] = " << sj.getNode() << endl;
623 : :
624 : : // Use siTimesCq and sjTimesCr to ensure deterministic node ID assignments
625 : 44169 : SumPair siTimesCq = si * cq;
626 : 44169 : SumPair sjTimesCr = sj * cr;
627 : 44169 : SumPair newSi = siTimesCq + sjTimesCr;
628 : :
629 : 44169 : const Polynomial& pi = d_trail[i].d_proof;
630 : 44169 : const Polynomial& pj = d_trail[j].d_proof;
631 : : // Use piTimesCq and pjTimesCr to ensure deterministic node ID assignments
632 : 44169 : Polynomial piTimesCq = pi * cq;
633 : 44169 : Polynomial pjTimesCr = pj * cr;
634 : 44169 : Polynomial newPi = piTimesCq + pjTimesCr;
635 : :
636 : 44169 : TrailIndex k = d_trail.size();
637 : 44169 : d_trail.push_back(Constraint(newSi, newPi));
638 : :
639 [ + - ]: 88338 : Trace("arith::dio") << "derived " << newSi.getNode() << " with proof "
640 : 44169 : << newPi.getNode() << endl;
641 : :
642 : 44169 : return k;
643 : 44169 : }
644 : :
645 : 0 : void DioSolver::printQueue()
646 : : {
647 [ - - ]: 0 : Trace("arith::dio") << "DioSolver::printQueue()" << endl;
648 [ - - ]: 0 : for (TrailIndex i = 0, last = d_trail.size(); i < last; ++i)
649 : : {
650 [ - - ]: 0 : Trace("arith::dio") << "d_trail[i].d_eq = " << d_trail[i].d_eq.getNode()
651 : 0 : << endl;
652 [ - - ]: 0 : Trace("arith::dio") << "d_trail[i].d_proof = "
653 : 0 : << d_trail[i].d_proof.getNode() << endl;
654 : : }
655 : :
656 [ - - ]: 0 : Trace("arith::dio") << "DioSolver::printSubs()" << endl;
657 [ - - ]: 0 : for (SubIndex si = 0, sN = d_subs.size(); si < sN; ++si)
658 : : {
659 [ - - ]: 0 : Trace("arith::dio") << "d_subs[i] = {"
660 : 0 : << "d_fresh=" << d_subs[si].d_fresh << ","
661 : 0 : << "d_eliminated=" << d_subs[si].d_eliminated.getNode()
662 : 0 : << ","
663 : 0 : << "d_constraint=" << d_subs[si].d_constraint << "}"
664 : 0 : << endl;
665 [ - - ]: 0 : Trace("arith::dio") << "d_trail[d_subs[i].d_constraint].d_eq="
666 : 0 : << d_trail[d_subs[si].d_constraint].d_eq.getNode()
667 : 0 : << endl;
668 : : }
669 : 0 : }
670 : :
671 : 40547 : DioSolver::TrailIndex DioSolver::applyAllSubstitutionsToIndex(
672 : : DioSolver::TrailIndex trailIndex)
673 : : {
674 : 40547 : TrailIndex currentIndex = trailIndex;
675 [ + + ]: 329521 : for (SubIndex subIter = 0, siEnd = d_subs.size(); subIter < siEnd; ++subIter)
676 : : {
677 : 288974 : currentIndex = applySubstitution(subIter, currentIndex);
678 : : }
679 : 40547 : return currentIndex;
680 : : }
681 : :
682 : 1102020 : bool DioSolver::debugSubstitutionApplies(DioSolver::SubIndex si,
683 : : DioSolver::TrailIndex ti)
684 : : {
685 : 1102020 : Variable var = d_subs[si].d_eliminated;
686 : :
687 : 1102020 : const SumPair& curr = d_trail[ti].d_eq;
688 : 1102020 : Polynomial vsum = curr.getPolynomial();
689 : :
690 : 2204040 : Constant a = vsum.getCoefficient(VarList(var));
691 : 2204040 : return !a.isZero();
692 : 1102020 : }
693 : :
694 : 89608 : bool DioSolver::debugAnySubstitionApplies(DioSolver::TrailIndex i)
695 : : {
696 [ + + ]: 1191628 : for (SubIndex subIter = 0, siEnd = d_subs.size(); subIter < siEnd; ++subIter)
697 : : {
698 [ - + ]: 1102020 : if (debugSubstitutionApplies(subIter, i))
699 : : {
700 : 0 : return true;
701 : : }
702 : : }
703 : 89608 : return false;
704 : : }
705 : :
706 : 26984 : std::pair<DioSolver::SubIndex, DioSolver::TrailIndex> DioSolver::solveIndex(
707 : : DioSolver::TrailIndex i)
708 : : {
709 : 26984 : const SumPair& si = d_trail[i].d_eq;
710 : :
711 [ + - ]: 53968 : Trace("arith::dio") << "before solveIndex(" << i << ":" << si.getNode() << ")"
712 : 26984 : << endl;
713 : :
714 : : #ifdef CVC5_ASSERTIONS
715 : 26984 : const Polynomial& p = si.getPolynomial();
716 : : #endif
717 : :
718 [ - + ][ - + ]: 26984 : Assert(p.isIntegral());
[ - - ]
719 : :
720 [ - + ][ - + ]: 26984 : Assert(p.selectAbsMinimum() == d_trail[i].d_minimalMonomial);
[ - - ]
721 : 26984 : const Monomial av = d_trail[i].d_minimalMonomial;
722 : :
723 : 26984 : VarList vl = av.getVarList();
724 [ - + ][ - + ]: 26984 : Assert(vl.singleton());
[ - - ]
725 : 26984 : Variable var = vl.getHead();
726 : 26984 : Constant a = av.getConstant();
727 : 26984 : Integer a_abs = a.getValue().getNumerator().abs();
728 : :
729 [ - + ][ - + ]: 26984 : Assert(a_abs == 1);
[ - - ]
730 : :
731 [ + + ][ + + ]: 26984 : TrailIndex ci = !a.isNegative() ? scaleEqAtIndex(i, Integer(-1)) : i;
[ - - ]
732 : :
733 : 26984 : SubIndex subBy = d_subs.size();
734 : 26984 : d_subs.push_back(Substitution(Node::null(), var, ci));
735 : :
736 [ + - ]: 53968 : Trace("arith::dio") << "after solveIndex " << d_trail[ci].d_eq.getNode()
737 : 26984 : << " for " << av.getNode() << endl;
738 [ - + ][ - + ]: 134920 : AssertEqual(d_trail[ci].d_eq.getPolynomial().getCoefficient(vl),
[ - - ]
739 : : Constant::mkConstant(nodeManager(), -1));
740 : :
741 : 53968 : return make_pair(subBy, i);
742 : 26984 : }
743 : :
744 : 1474 : std::pair<DioSolver::SubIndex, DioSolver::TrailIndex> DioSolver::decomposeIndex(
745 : : DioSolver::TrailIndex i)
746 : : {
747 : 1474 : const SumPair& si = d_trail[i].d_eq;
748 : :
749 : 1474 : d_usedDecomposeIndex = true;
750 : :
751 [ + - ]: 2948 : Trace("arith::dio") << "before decomposeIndex(" << i << ":" << si.getNode()
752 : 1474 : << ")" << endl;
753 : :
754 : : #ifdef CVC5_ASSERTIONS
755 : 1474 : const Polynomial& p = si.getPolynomial();
756 : : #endif
757 : :
758 [ - + ][ - + ]: 1474 : Assert(p.isIntegral());
[ - - ]
759 : :
760 [ - + ][ - + ]: 1474 : Assert(p.selectAbsMinimum() == d_trail[i].d_minimalMonomial);
[ - - ]
761 : 1474 : const Monomial& av = d_trail[i].d_minimalMonomial;
762 : :
763 : 1474 : VarList vl = av.getVarList();
764 [ - + ][ - + ]: 1474 : Assert(vl.singleton());
[ - - ]
765 : 1474 : Variable var = vl.getHead();
766 : 1474 : Constant a = av.getConstant();
767 : 1474 : Integer a_abs = a.getValue().getNumerator().abs();
768 : :
769 [ - + ][ - + ]: 1474 : Assert(a_abs > 1);
[ - - ]
770 : :
771 : : // It is not sufficient to reduce the case where abs(a) == 1 to abs(a) > 1.
772 : : // We need to handle both cases seperately to ensure termination.
773 : 1474 : Node qr = SumPair::computeQR(si, a.getValue().getNumerator());
774 : :
775 [ - + ][ - + ]: 1474 : Assert(qr.getKind() == Kind::ADD);
[ - - ]
776 [ - + ][ - + ]: 1474 : Assert(qr.getNumChildren() == 2);
[ - - ]
777 : 2948 : SumPair q = SumPair::parseSumPair(qr[0]);
778 : 2948 : SumPair r = SumPair::parseSumPair(qr[1]);
779 : :
780 : 1474 : NodeManager* nm = nodeManager();
781 [ - + ][ - + ]: 7370 : AssertEqual(q.getPolynomial().getCoefficient(vl),
[ - - ]
782 : : Constant::mkConstant(nm, 1));
783 : :
784 [ - + ][ - + ]: 1474 : Assert(!r.isZero());
[ - - ]
785 : 1474 : Node freshNode = makeIntegerVariable(nodeManager());
786 : 1474 : Variable fresh(freshNode);
787 : 1474 : SumPair fresh_one = SumPair::mkSumPair(fresh);
788 : 1474 : SumPair fresh_a = fresh_one * a;
789 : :
790 : 1474 : SumPair newSI = SumPair(fresh_one) - q;
791 : : // this normalizes the coefficient of var to -1
792 : :
793 : 1474 : TrailIndex ci = d_trail.size();
794 : 1474 : d_trail.push_back(Constraint(newSI, Polynomial::mkZero(nm)));
795 : : // no longer reference av safely!
796 : 1474 : addTrailElementAsLemma(ci);
797 : :
798 [ + - ]: 2948 : Trace("arith::dio") << "Decompose ci(" << ci << ":"
799 : 0 : << d_trail[ci].d_eq.getNode() << ") for "
800 : 1474 : << d_trail[i].d_minimalMonomial.getNode() << endl;
801 [ - + ][ - + ]: 7370 : AssertEqual(d_trail[ci].d_eq.getPolynomial().getCoefficient(vl),
[ - - ]
802 : : Constant::mkConstant(nm, -1));
803 : :
804 : 1474 : SumPair newFact = r + fresh_a;
805 : :
806 : 1474 : TrailIndex nextIndex = d_trail.size();
807 : 1474 : d_trail.push_back(Constraint(newFact, d_trail[i].d_proof));
808 : :
809 : 1474 : SubIndex subBy = d_subs.size();
810 : 1474 : d_subs.push_back(Substitution(freshNode, var, ci));
811 : :
812 [ + - ]: 2948 : Trace("arith::dio") << "Decompose nextIndex "
813 : 1474 : << d_trail[nextIndex].d_eq.getNode() << endl;
814 : 2948 : return make_pair(subBy, nextIndex);
815 : 1474 : }
816 : :
817 : 728467 : DioSolver::TrailIndex DioSolver::applySubstitution(DioSolver::SubIndex si,
818 : : DioSolver::TrailIndex ti)
819 : : {
820 : 728467 : Variable var = d_subs[si].d_eliminated;
821 : 728467 : TrailIndex subIndex = d_subs[si].d_constraint;
822 : :
823 : 728467 : const SumPair& curr = d_trail[ti].d_eq;
824 : 728467 : Polynomial vsum = curr.getPolynomial();
825 : :
826 : 1456934 : Constant a = vsum.getCoefficient(VarList(var));
827 [ - + ][ - + ]: 728467 : Assert(a.isIntegral());
[ - - ]
828 [ + + ]: 728467 : if (!a.isZero())
829 : : {
830 : 43889 : Integer one(1);
831 : : TrailIndex afterSub =
832 : 43889 : combineEqAtIndexes(ti, one, subIndex, a.getValue().getNumerator());
833 [ - + ][ - + ]: 43889 : Assert(d_trail[afterSub]
[ - - ]
834 : : .d_eq.getPolynomial()
835 : : .getCoefficient(VarList(var))
836 : : .isZero());
837 : 43889 : return afterSub;
838 : 43889 : }
839 : : else
840 : : {
841 : 684578 : return ti;
842 : : }
843 : 728467 : }
844 : :
845 : 62221 : DioSolver::TrailIndex DioSolver::reduceByGCD(DioSolver::TrailIndex ti)
846 : : {
847 : 62221 : const SumPair& sp = d_trail[ti].d_eq;
848 : 62221 : Polynomial vsum = sp.getPolynomial();
849 : 62221 : Constant c = sp.getConstant();
850 : :
851 [ + - ]: 62221 : Trace("arith::dio") << "reduceByGCD " << vsum.getNode() << endl;
852 [ - + ][ - + ]: 62221 : Assert(!vsum.isConstant());
[ - - ]
853 : 62221 : Integer g = vsum.gcd();
854 [ - + ][ - + ]: 62221 : Assert(g >= 1);
[ - - ]
855 [ + - ]: 124442 : Trace("arith::dio") << "gcd(" << vsum.getNode() << ")=" << g << " "
856 : 62221 : << c.getValue() << endl;
857 [ + + ]: 62221 : if (g.divides(c.getValue().getNumerator()))
858 : : {
859 [ + + ]: 60022 : if (g > 1)
860 : : {
861 : 2369 : return scaleEqAtIndex(ti, g);
862 : : }
863 : : else
864 : : {
865 : 57653 : return ti;
866 : : }
867 : : }
868 : : else
869 : : {
870 : 2199 : raiseConflict(ti);
871 : 2199 : return ti;
872 : : }
873 : 62221 : }
874 : :
875 : 196430 : bool DioSolver::triviallySat(TrailIndex i)
876 : : {
877 : 196430 : const SumPair& eq = d_trail[i].d_eq;
878 [ + + ]: 196430 : if (eq.isConstant())
879 : : {
880 : 4401 : return eq.getConstant().isZero();
881 : : }
882 : : else
883 : : {
884 : 192029 : return false;
885 : : }
886 : : }
887 : :
888 : 195499 : bool DioSolver::triviallyUnsat(DioSolver::TrailIndex i)
889 : : {
890 : 195499 : const SumPair& eq = d_trail[i].d_eq;
891 [ + + ]: 195499 : if (eq.isConstant())
892 : : {
893 : 3470 : return !eq.getConstant().isZero();
894 : : }
895 : : else
896 : : {
897 : 192029 : return false;
898 : : }
899 : : }
900 : :
901 : 89608 : bool DioSolver::gcdIsOne(DioSolver::TrailIndex i)
902 : : {
903 : 89608 : const SumPair& eq = d_trail[i].d_eq;
904 : 89608 : return eq.gcd() == Integer(1);
905 : : }
906 : :
907 : 28458 : void DioSolver::subAndReduceCurrentFByIndex(DioSolver::SubIndex subIndex)
908 : : {
909 : 28458 : size_t N = d_currentF.size();
910 : :
911 : 28458 : size_t readIter = 0, writeIter = 0;
912 [ + + ][ + + ]: 467951 : for (; readIter < N && !inConflict(); ++readIter)
[ + + ]
913 : : {
914 : 439493 : TrailIndex curr = d_currentF[readIter];
915 : 439493 : TrailIndex nextTI = applySubstitution(subIndex, curr);
916 [ + + ]: 439493 : if (nextTI == curr)
917 : : {
918 : 413418 : d_currentF[writeIter] = curr;
919 : 413418 : ++writeIter;
920 : : }
921 : : else
922 : : {
923 [ - + ][ - + ]: 26075 : Assert(nextTI != curr);
[ - - ]
924 : :
925 [ - + ]: 26075 : if (triviallyUnsat(nextTI))
926 : : {
927 : 0 : raiseConflict(nextTI);
928 : : }
929 [ + + ]: 26075 : else if (!triviallySat(nextTI))
930 : : {
931 : 22605 : TrailIndex nextNextTI = reduceByGCD(nextTI);
932 : :
933 [ + + ][ + + ]: 22605 : if (!(inConflict() || anyCoefficientExceedsMaximum(nextNextTI)))
[ + + ]
934 : : {
935 [ - + ][ - + ]: 20978 : Assert(queueConditions(nextNextTI));
[ - - ]
936 : 20978 : d_currentF[writeIter] = nextNextTI;
937 : 20978 : ++writeIter;
938 : : }
939 : : }
940 : : }
941 : : }
942 [ + + ][ + + ]: 28458 : if (!inConflict() && writeIter < N)
[ + + ]
943 : : {
944 : 2459 : d_currentF.resize(writeIter);
945 : : }
946 : 28458 : }
947 : :
948 : 1474 : void DioSolver::addTrailElementAsLemma(TrailIndex i)
949 : : {
950 [ - + ]: 1474 : if (options().arith.exportDioDecompositions)
951 : : {
952 : 0 : d_decompositionLemmaQueue.push(i);
953 : : }
954 : 1474 : }
955 : :
956 : 0 : Node DioSolver::trailIndexToEquality(TrailIndex i) const
957 : : {
958 : 0 : const SumPair& sp = d_trail[i].d_eq;
959 : 0 : Node n = sp.getNode();
960 : 0 : Node zero = nodeManager()->mkConstRealOrInt(n.getType(), Rational(0));
961 : 0 : Node eq = n.eqNode(zero);
962 : 0 : return eq;
963 : 0 : }
964 : :
965 : : } // namespace arith::linear
966 : : } // namespace theory
967 : : } // namespace cvc5::internal
|
