GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/dio_solver.cpp Lines: 427 489 87.3 %
Date: 2021-03-23 Branches: 844 2372 35.6 %

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