carl  24.04
Computer ARithmetic Library
substitute.tpp
Go to the documentation of this file.
1 #include <cmath>
2 #include <limits>
3 
4 #include <carl-arith/poly/umvpoly/functions/Derivative.h>
5 #include <carl-arith/poly/umvpoly/functions/SoSDecomposition.h>
6 #include <carl-arith/constraint/IntervalEvaluation.h>
7 
8 #include <carl-arith/vs/Substitution.h>
9 
10 //#define VS_DEBUG_SUBSTITUTION
11 const unsigned MAX_NUM_OF_TERMS = 512;
12 
13 
14 namespace carl::vs::detail
15 {
16  template <class combineType>
17  bool combine( const std::vector< std::vector< std::vector< combineType > > >& _toCombine, std::vector< std::vector< combineType > >& _combination )
18  {
19  // Initialize the current combination. If there is nothing to combine or an empty vector to combine with, return the empty vector.
20  if( _toCombine.empty() ) return true;
21  std::vector<typename std::vector< std::vector< combineType > >::const_iterator > combination
22  = std::vector<typename std::vector< std::vector< combineType > >::const_iterator >();
23  unsigned estimatedResultSize = 1;
24  for( auto iter = _toCombine.begin(); iter != _toCombine.end(); ++iter )
25  {
26  if( iter->empty() )
27  return true;
28  estimatedResultSize *= (unsigned)iter->size();
29  if( estimatedResultSize > MAX_NUM_OF_COMBINATION_RESULT )
30  return false;
31  else
32  combination.push_back( iter->begin() );
33  }
34  // As long as no more combination for the last vector in first vector of vectors exists.
35  bool lastCombinationReached = false;
36  while( !lastCombinationReached )
37  {
38  // Create a new combination of vectors.
39  _combination.push_back( std::vector< combineType >() );
40  bool previousCounterIncreased = false;
41  // For each vector in the vector of vectors, choose a vector. Finally we combine
42  // those vectors by merging them to a new vector and add this to the result.
43  auto currentVector = _toCombine.begin();
44  for( auto combineElement = combination.begin(); combineElement != combination.end(); ++combineElement )
45  {
46  // Add the current vectors elements to the combination.
47  _combination.back().insert( _combination.back().end(), (*combineElement)->begin(), (*combineElement)->end() );
48  // Set the counter.
49  if( !previousCounterIncreased )
50  {
51  ++(*combineElement);
52  if( *combineElement != currentVector->end() )
53  previousCounterIncreased = true;
54  else
55  {
56  if( combineElement != --combination.end() )
57  *combineElement = currentVector->begin();
58  else
59  lastCombinationReached = true;
60  }
61  }
62  ++currentVector;
63  }
64  }
65  return true;
66  }
67 
68  template<typename Poly>
69  void simplify( CaseDistinction<Poly>& _toSimplify )
70  {
71  bool containsEmptyDisjunction = false;
72  auto conj = _toSimplify.begin();
73  while( conj != _toSimplify.end() )
74  {
75  bool conjInconsistent = false;
76  auto cons = (*conj).begin();
77  while( cons != (*conj).end() )
78  {
79  if( *cons == Constraint<Poly>( false ) )
80  {
81  conjInconsistent = true;
82  break;
83  }
84  else if( *cons == Constraint<Poly>( true ) )
85  cons = (*conj).erase( cons );
86  else
87  cons++;
88  }
89  bool conjEmpty = (*conj).empty();
90  if( conjInconsistent || (containsEmptyDisjunction && conjEmpty) )
91  {
92  // Delete the conjunction.
93  conj->clear();
94  conj = _toSimplify.erase( conj );
95  }
96  else
97  conj++;
98  if( !containsEmptyDisjunction && conjEmpty )
99  containsEmptyDisjunction = true;
100  }
101  }
102 
103  template<typename Poly>
104  void simplify( CaseDistinction<Poly>& _toSimplify, Variables& _conflictingVars, const detail::EvalDoubleIntervalMap& _solutionSpace )
105  {
106  bool containsEmptyDisjunction = false;
107  auto conj = _toSimplify.begin();
108  while( conj != _toSimplify.end() )
109  {
110  bool conjInconsistent = false;
111  auto cons = (*conj).begin();
112  while( cons != (*conj).end() )
113  {
114  if( *cons == Constraint<Poly>( false ) )
115  {
116  conjInconsistent = true;
117  break;
118  }
119  else if( *cons == Constraint<Poly>( true ) )
120  cons = (*conj).erase( cons );
121  else
122  {
123  unsigned conflictingWithSolutionSpace = consistent_with(cons->constr(), _solutionSpace );
124 
125 // std::cout << "Is " << cons << std::endl;
126 // std::cout << std::endl;
127 // std::cout << "consistent with " << std::endl;
128 // for( auto iter = _solutionSpace.begin(); iter != _solutionSpace.end(); ++iter )
129 // std::cout << iter->first << " in " << iter->second << std::endl;
130 // std::cout << " -> " << conflictingWithSolutionSpace << std::endl;
131 
132  if( conflictingWithSolutionSpace == 0 )
133  {
134  auto vars = variables(*cons);
135  _conflictingVars.insert( vars.begin(), vars.end() );
136  conjInconsistent = true;
137  break;
138  }
139  else if( conflictingWithSolutionSpace == 1 )
140  cons = (*conj).erase( cons );
141  else
142  ++cons;
143  }
144  }
145  bool conjEmpty = (*conj).empty();
146  if( conjInconsistent || (containsEmptyDisjunction && conjEmpty) )
147  {
148  // Delete the conjunction.
149  conj->clear();
150  conj = _toSimplify.erase( conj );
151  }
152  else
153  ++conj;
154  if( !containsEmptyDisjunction && conjEmpty )
155  containsEmptyDisjunction = true;
156  }
157  }
158 
159  template<typename Poly>
160  bool splitProducts( CaseDistinction<Poly>& _toSimplify, bool _onlyNeq )
161  {
162  std::map<const Constraint<Poly>, CaseDistinction<Poly>> result_cache;
163  bool result = true;
164  size_t toSimpSize = _toSimplify.size();
165  for( size_t pos = 0; pos < toSimpSize; )
166  {
167  if( !_toSimplify.begin()->empty() )
168  {
169  CaseDistinction<Poly> temp;
170  if( !splitProducts( _toSimplify[pos], temp, result_cache, _onlyNeq ) )
171  result = false;
172  _toSimplify.erase( _toSimplify.begin() );
173  _toSimplify.insert( _toSimplify.end(), temp.begin(), temp.end() );
174  --toSimpSize;
175  }
176  else
177  ++pos;
178  }
179  return result;
180  }
181 
182  template<typename Poly>
183  bool splitProducts( const ConstraintConjunction<Poly>& _toSimplify, CaseDistinction<Poly>& _result, std::map<const Constraint<Poly>, CaseDistinction<Poly>>& result_cache, bool _onlyNeq )
184  {
185  std::vector<CaseDistinction<Poly>> toCombine;
186  for( auto constraint = _toSimplify.begin(); constraint != _toSimplify.end(); ++constraint )
187  {
188  auto i = result_cache.find(*constraint);
189  if (i == result_cache.end()) {
190  auto res = splitProducts(*constraint, _onlyNeq);
191  i = result_cache.emplace(*constraint, res).first;
192  }
193  toCombine.push_back(i->second);
194  }
195  bool result = true;
196  if( !combine( toCombine, _result ) )
197  result = false;
198  simplify( _result );
199  return result;
200  }
201 
202  template<typename Poly>
203  CaseDistinction<Poly> splitProducts( const Constraint<Poly>& _constraint, bool _onlyNeq )
204  {
205  CaseDistinction<Poly> result;
206  auto& factorization = _constraint.lhs_factorization();
207  if( !carl::is_trivial(factorization) )
208  {
209  switch( _constraint.relation() )
210  {
211  case Relation::EQ:
212  {
213  if( !_onlyNeq )
214  {
215  for( auto factor = factorization.begin(); factor != factorization.end(); ++factor )
216  {
217  result.emplace_back();
218  result.back().push_back( Constraint<Poly>( factor->first, Relation::EQ ) );
219  }
220  }
221  else
222  {
223  result.emplace_back();
224  result.back().push_back( _constraint );
225  }
226  simplify( result );
227  break;
228  }
229  case Relation::NEQ:
230  {
231  result.emplace_back();
232  for( auto factor = factorization.begin(); factor != factorization.end(); ++factor )
233  result.back().push_back( Constraint<Poly>( factor->first, Relation::NEQ ) );
234  simplify( result );
235  break;
236  }
237  default:
238  {
239  if( !_onlyNeq )
240  {
241  result = getSignCombinations( _constraint );
242  simplify( result );
243  }
244  else
245  {
246  result.emplace_back();
247  result.back().push_back( _constraint );
248  }
249  }
250 
251  }
252  }
253  else
254  {
255  result.emplace_back();
256  result.back().push_back( _constraint );
257  }
258  return result;
259  }
260 
261  template<typename Poly>
262  void splitSosDecompositions( CaseDistinction<Poly>& _toSimplify )
263  {
264  // TODO this needs to be reviewed and fixed
265  // It seems that is is assumed that lcoeffNeg * constraint.lhs() is positive if
266  // a SOS decomposition exists. However, we can only follow that it's non-negative
267  // (in the univariate case), for the multivariate case more requirements need to be made
268  // (therefor constraint.lhs().has_constant_term() ???).
269  // This lead to wrong simplifications in very rare cases, for example
270  // -100 + 140*z + -49*y^2 + -49*z^2 is wrongly simplified to true (see McsatVSBug).
271 
272  // Temporarily disabled (what can be followed in the multivariate case?).
273  return;
274 
275  for( size_t i = 0; i < _toSimplify.size(); )
276  {
277  auto& cc = _toSimplify[i];
278  bool foundNoInvalidConstraint = true;
279  size_t pos = 0;
280  while( foundNoInvalidConstraint && pos < cc.size() )
281  {
282  const Constraint<Poly>& constraint = cc[pos];
283  std::vector<std::pair<typename Poly::NumberType,Poly>> sosDec;
284  bool lcoeffNeg = carl::is_negative(constraint.lhs().lcoeff());
285  if (lcoeffNeg)
286  sosDec = carl::sos_decomposition(-constraint.lhs());
287  else
288  sosDec = carl::sos_decomposition(constraint.lhs());
289  if( sosDec.size() > 1 )
290  {
291 // std::cout << "Sum-of-squares decomposition of " << constraint.lhs() << " = " << sosDec << std::endl;
292  bool addSquares = true;
293  bool constraintValid = false;
294  switch( constraint.relation() )
295  {
296  case carl::Relation::EQ:
297  {
298  if( constraint.lhs().has_constant_term() )
299  {
300  foundNoInvalidConstraint = false;
301  addSquares = false;
302  }
303  break;
304  }
305  case carl::Relation::NEQ:
306  {
307  addSquares = false;
308  if( constraint.lhs().has_constant_term() )
309  {
310  constraintValid = true;
311  }
312  break;
313  }
314  case carl::Relation::LEQ:
315  {
316  if( lcoeffNeg )
317  {
318  addSquares = false;
319  constraintValid = true;
320  }
321  else if( constraint.lhs().has_constant_term() )
322  {
323  addSquares = false;
324  foundNoInvalidConstraint = false;
325  }
326  break;
327  }
328  case carl::Relation::LESS:
329  {
330  addSquares = false;
331  if( lcoeffNeg )
332  {
333  if( constraint.lhs().has_constant_term() )
334  constraintValid = true;
335  }
336  else
337  foundNoInvalidConstraint = false;
338  break;
339  }
340  case carl::Relation::GEQ:
341  {
342  if( !lcoeffNeg )
343  {
344  addSquares = false;
345  constraintValid = true;
346  }
347  else if( constraint.lhs().has_constant_term() )
348  {
349  addSquares = false;
350  foundNoInvalidConstraint = false;
351  }
352  break;
353  }
354  default:
355  {
356  assert( constraint.relation() == carl::Relation::GREATER );
357  addSquares = false;
358  if( lcoeffNeg )
359  foundNoInvalidConstraint = false;
360  else
361  {
362  if( constraint.lhs().has_constant_term() )
363  constraintValid = true;
364  }
365  }
366  }
367  assert( !(!foundNoInvalidConstraint && constraintValid) );
368  assert( !(!foundNoInvalidConstraint && addSquares) );
369  if( constraintValid || addSquares )
370  {
371  cc[pos] = cc.back();
372  cc.pop_back();
373  }
374  else
375  ++pos;
376  if( addSquares )
377  {
378  for( auto it = sosDec.begin(); it != sosDec.end(); ++it )
379  {
380  cc.emplace_back( it->second, carl::Relation::EQ );
381  }
382  }
383  }
384  else
385  ++pos;
386  }
387  if( foundNoInvalidConstraint )
388  ++i;
389  else
390  {
391  cc = _toSimplify.back();
392  _toSimplify.pop_back();
393  }
394  }
395  }
396 
397  template<typename Poly>
398  CaseDistinction<Poly> getSignCombinations( const Constraint<Poly>& _constraint )
399  {
400  CaseDistinction<Poly> combinations;
401  auto& factorization = _constraint.lhs_factorization();
402  if( !carl::is_trivial(factorization) && factorization.size() <= MAX_PRODUCT_SPLIT_NUMBER )
403  {
404  assert( _constraint.relation() == Relation::GREATER || _constraint.relation() == Relation::LESS
405  || _constraint.relation() == Relation::GEQ || _constraint.relation() == Relation::LEQ );
406  Relation relPos = Relation::GREATER;
407  Relation relNeg = Relation::LESS;
408  if( _constraint.relation() == Relation::GEQ || _constraint.relation() == Relation::LEQ )
409  {
410  relPos = Relation::GEQ;
411  relNeg = Relation::LEQ;
412  }
413  bool positive = (_constraint.relation() == Relation::GEQ || _constraint.relation() == Relation::GREATER);
414  ConstraintConjunction<Poly> positives;
415  ConstraintConjunction<Poly> alwayspositives;
416  ConstraintConjunction<Poly> negatives;
417  ConstraintConjunction<Poly> alwaysnegatives;
418  unsigned numOfAlwaysNegatives = 0;
419  for( auto factor = factorization.begin(); factor != factorization.end(); ++factor )
420  {
421  Constraint<Poly> consPos = Constraint<Poly>( factor->first, relPos );
422  unsigned posConsistent = consPos.is_consistent();
423  if( posConsistent != 0 )
424  positives.push_back( consPos );
425  Constraint<Poly> consNeg = Constraint<Poly>( factor->first, relNeg );
426  unsigned negConsistent = consNeg.is_consistent();
427  if( negConsistent == 0 )
428  {
429  if( posConsistent == 0 )
430  {
431  combinations.emplace_back();
432  combinations.back().push_back( consNeg );
433  return combinations;
434  }
435  if( posConsistent != 1 )
436  alwayspositives.push_back( positives.back() );
437  positives.pop_back();
438  }
439  else
440  {
441  if( posConsistent == 0 )
442  {
443  ++numOfAlwaysNegatives;
444  if( negConsistent != 1 )
445  alwaysnegatives.push_back( consNeg );
446  }
447  else negatives.push_back( consNeg );
448  }
449  }
450  assert( positives.size() == negatives.size() );
451  if( positives.size() > 0 )
452  {
453  std::vector< std::bitset<MAX_PRODUCT_SPLIT_NUMBER> > combSelector = std::vector< std::bitset<MAX_PRODUCT_SPLIT_NUMBER> >();
454  if( fmod( numOfAlwaysNegatives, 2.0 ) != 0.0 )
455  {
456  if( positive )
457  getOddBitStrings( positives.size(), combSelector );
458  else
459  getEvenBitStrings( positives.size(), combSelector );
460  }
461  else
462  {
463  if( positive )
464  getEvenBitStrings( positives.size(), combSelector );
465  else
466  getOddBitStrings( positives.size(), combSelector );
467  }
468  for( auto comb = combSelector.begin(); comb != combSelector.end(); ++comb )
469  {
470  combinations.emplace_back( alwaysnegatives );
471  combinations.back().insert( combinations.back().end(), alwayspositives.begin(), alwayspositives.end() );
472  for( unsigned pos = 0; pos < positives.size(); ++pos )
473  {
474  if( (*comb)[pos] )
475  combinations.back().push_back( negatives[pos] );
476  else
477  combinations.back().push_back( positives[pos] );
478  }
479  }
480  }
481  else
482  {
483  combinations.emplace_back( alwaysnegatives );
484  combinations.back().insert( combinations.back().end(), alwayspositives.begin(), alwayspositives.end() );
485  }
486  }
487  else
488  {
489  combinations.emplace_back();
490  combinations.back().push_back( _constraint );
491  }
492  return combinations;
493  }
494 
495  void getOddBitStrings( std::size_t _length, std::vector< std::bitset<MAX_PRODUCT_SPLIT_NUMBER> >& _strings )
496  {
497  assert( _length > 0 );
498  if( _length == 1 ) _strings.push_back( std::bitset<MAX_PRODUCT_SPLIT_NUMBER>( 1 ) );
499  else
500  {
501  std::size_t pos = _strings.size();
502  getEvenBitStrings( _length - 1, _strings );
503  for( ; pos < _strings.size(); ++pos )
504  {
505  _strings[pos] <<= 1;
506  _strings[pos].flip(0);
507  }
508  getOddBitStrings( _length - 1, _strings );
509  for( ; pos < _strings.size(); ++pos ) _strings[pos] <<= 1;
510  }
511  }
512 
513  void getEvenBitStrings( std::size_t _length, std::vector< std::bitset<MAX_PRODUCT_SPLIT_NUMBER> >& _strings )
514  {
515  assert( _length > 0 );
516  if( _length == 1 ) _strings.push_back( std::bitset<MAX_PRODUCT_SPLIT_NUMBER>( 0 ) );
517  else
518  {
519  std::size_t pos = _strings.size();
520  getEvenBitStrings( _length - 1, _strings );
521  for( ; pos < _strings.size(); ++pos ) _strings[pos] <<= 1;
522  getOddBitStrings( _length - 1, _strings );
523  for( ; pos < _strings.size(); ++pos )
524  {
525  _strings[pos] <<= 1;
526  _strings[pos].flip(0);
527  }
528  }
529  }
530 
531  template<typename Poly>
532  void print( CaseDistinction<Poly>& _substitutionResults )
533  {
534  auto conj = _substitutionResults.begin();
535  while( conj != _substitutionResults.end() )
536  {
537  if( conj != _substitutionResults.begin() )
538  std::cout << " or (";
539  else
540  std::cout << " (";
541  auto cons = (*conj).begin();
542  while( cons != (*conj).end() )
543  {
544  if( cons != (*conj).begin() )
545  std::cout << " and ";
546  std::cout << *cons;
547  cons++;
548  }
549  std::cout << ")" << std::endl;
550  conj++;
551  }
552  std::cout << std::endl;
553  }
554 
555  template<typename Poly>
556  bool substitute( const Constraint<Poly>& _cons,
557  const Substitution<Poly>& _subs,
558  CaseDistinction<Poly>& _result,
559  bool _accordingPaper,
560  Variables& _conflictingVariables,
561  const detail::EvalDoubleIntervalMap& _solutionSpace)
562  {
563  #ifdef VS_DEBUG_SUBSTITUTION
564  std::cout << "substitute: ( " << _cons << " )" << _subs << std::endl;
565  #endif
566  bool result = true;
567  // Apply the substitution according to its type.
568  switch( _subs.term().type() )
569  {
570  case TermType::NORMAL:
571  {
572  result = substituteNormal( _cons, _subs, _result, _accordingPaper, _conflictingVariables, _solutionSpace );
573  break;
574  }
575  case TermType::PLUS_EPSILON:
576  {
577  result = substitutePlusEps( _cons, _subs, _result, _accordingPaper, _conflictingVariables, _solutionSpace );
578  break;
579  }
580  case TermType::MINUS_INFINITY:
581  {
582  substituteInf( _cons, _subs, _result, _conflictingVariables, _solutionSpace );
583  break;
584  }
585  case TermType::PLUS_INFINITY:
586  {
587  substituteInf( _cons, _subs, _result, _conflictingVariables, _solutionSpace );
588  break;
589  }
590  default:
591  {
592  std::cout << "Error in substitute: unexpected type of substitution." << std::endl;
593  }
594  }
595  #ifdef VS_DEBUG_SUBSTITUTION
596  print( _result );
597  #endif
598  return result;
599  }
600 
601  template<typename Poly>
602  bool substituteNormal( const Constraint<Poly>& _cons,
603  const Substitution<Poly>& _subs,
604  CaseDistinction<Poly>& _result,
605  bool _accordingPaper,
606  Variables& _conflictingVariables,
607  const detail::EvalDoubleIntervalMap& _solutionSpace )
608  {
609 
610  bool result = true;
611  if( _cons.variables().has( _subs.variable() ) )
612  {
613  using Rational = typename Poly::NumberType;
614  // Collect all necessary left hand sides to create the new conditions of all cases referring to the virtual substitution.
615  if( carl::pow( Rational(Rational(_subs.term().sqrt_ex().constant_part().size()) + Rational(_subs.term().sqrt_ex().factor().size()) * Rational(_subs.term().sqrt_ex().radicand().size())), _cons.maxDegree( _subs.variable() )) > (MAX_NUM_OF_TERMS*MAX_NUM_OF_TERMS) )
616  {
617  return false;
618  }
619  carl::SqrtEx sub = carl::substitute( _cons.lhs(), _subs.variable(), _subs.term().sqrt_ex() );
620  #ifdef VS_DEBUG_SUBSTITUTION
621  std::cout << "Result of common substitution:" << sub << std::endl;
622  #endif
623  // The term then looks like: q/s
624  if( !sub.has_sqrt() )
625  {
626  // Create the new decision tuples.
627  if( _cons.relation() == Relation::EQ || _cons.relation() == Relation::NEQ )
628  {
629  // Add conjunction (sub.constant_part() = 0) to the substitution result.
630  _result.emplace_back();
631  _result.back().push_back( Constraint<Poly>( sub.constant_part(), _cons.relation() ) );
632  }
633  else
634  {
635  if( !_subs.term().sqrt_ex().denominator().is_constant() )
636  {
637  // Add conjunction (sub.denominator()>0 and sub.constant_part() </>/<=/>= 0) to the substitution result.
638  _result.emplace_back();
639  _result.back().push_back( Constraint<Poly>( sub.denominator(), Relation::GREATER ) );
640  _result.back().push_back( Constraint<Poly>( sub.constant_part(), _cons.relation() ) );
641  // Add conjunction (sub.denominator()<0 and sub.constant_part() >/</>=/<= 0) to the substitution result.
642  Relation inverseRelation;
643  switch( _cons.relation() )
644  {
645  case Relation::LESS:
646  inverseRelation = Relation::GREATER;
647  break;
648  case Relation::GREATER:
649  inverseRelation = Relation::LESS;
650  break;
651  case Relation::LEQ:
652  inverseRelation = Relation::GEQ;
653  break;
654  case Relation::GEQ:
655  inverseRelation = Relation::LEQ;
656  break;
657  default:
658  assert( false );
659  inverseRelation = Relation::EQ;
660  }
661  _result.emplace_back();
662  _result.back().push_back( Constraint<Poly>( sub.denominator(), Relation::LESS ) );
663  _result.back().push_back( Constraint<Poly>( sub.constant_part(), inverseRelation ) );
664  }
665  else
666  {
667  // Add conjunction (f(-c/b)*b^k </>/<=/>= 0) to the substitution result.
668  _result.emplace_back();
669  _result.back().push_back( Constraint<Poly>( sub.constant_part(), _cons.relation() ) );
670  }
671  }
672  }
673  // The term then looks like: (q+r*sqrt(b^2-4ac))/s
674  else
675  {
676  Poly s = Poly(1);
677  if( !_subs.term().sqrt_ex().denominator().is_constant() )
678  s = sub.denominator();
679  switch( _cons.relation() )
680  {
681  case Relation::EQ:
682  {
683  result = substituteNormalSqrtEq( sub.radicand(), sub.constant_part(), sub.factor(), _result, _accordingPaper );
684  break;
685  }
686  case Relation::NEQ:
687  {
688  result = substituteNormalSqrtNeq( sub.radicand(), sub.constant_part(), sub.factor(), _result, _accordingPaper );
689  break;
690  }
691  case Relation::LESS:
692  {
693  result = substituteNormalSqrtLess( sub.radicand(), sub.constant_part(), sub.factor(), s, _result, _accordingPaper );
694  break;
695  }
696  case Relation::GREATER:
697  {
698  result = substituteNormalSqrtLess( sub.radicand(), sub.constant_part(), sub.factor(), -s, _result, _accordingPaper );
699  break;
700  }
701  case Relation::LEQ:
702  {
703  result = substituteNormalSqrtLeq( sub.radicand(), sub.constant_part(), sub.factor(), s, _result, _accordingPaper );
704  break;
705  }
706  case Relation::GEQ:
707  {
708  result = substituteNormalSqrtLeq( sub.radicand(), sub.constant_part(), sub.factor(), -s, _result, _accordingPaper );
709  break;
710  }
711  default:
712  std::cout << "Error in substituteNormal: Unexpected relation symbol" << std::endl;
713  assert( false );
714  }
715  }
716  }
717  else
718  {
719  _result.emplace_back();
720  _result.back().push_back( _cons );
721  }
722  simplify( _result, _conflictingVariables, _solutionSpace );
723  return result;
724  }
725 
726  template<typename Poly>
727  bool substituteNormalSqrtEq( const Poly& _radicand,
728  const Poly& _q,
729  const Poly& _r,
730  CaseDistinction<Poly>& _result,
731  bool _accordingPaper )
732  {
733  if( _q.size() > MAX_NUM_OF_TERMS || _r.size() > MAX_NUM_OF_TERMS || _radicand.size() > MAX_NUM_OF_TERMS )
734  return false;
735  Poly lhs = carl::pow(_q, 2) - carl::pow(_r, 2) * _radicand;
736  if( _accordingPaper )
737  {
738  Poly qr = _q * _r;
739  // Add conjunction (q*r<=0 and q^2-r^2*radicand=0) to the substitution result.
740  _result.emplace_back();
741  _result.back().push_back( Constraint<Poly>( qr, Relation::LEQ ) );
742  _result.back().push_back( Constraint<Poly>( lhs, Relation::EQ ) );
743  }
744  else
745  {
746  // Add conjunction (q=0 and r=0) to the substitution result.
747  _result.emplace_back();
748  _result.back().push_back( Constraint<Poly>( _q, Relation::EQ ) );
749  _result.back().push_back( Constraint<Poly>( _r, Relation::EQ ) );
750  // Add conjunction (q=0 and radicand=0) to the substitution result.
751  _result.emplace_back();
752  _result.back().push_back( Constraint<Poly>( _q, Relation::EQ ) );
753  _result.back().push_back( Constraint<Poly>( _radicand, Relation::EQ ) );
754  // Add conjunction (q<0 and r>0 and q^2-r^2*radicand=0) to the substitution result.
755  _result.emplace_back();
756  _result.back().push_back( Constraint<Poly>( _q, Relation::LESS ) );
757  _result.back().push_back( Constraint<Poly>( _r, Relation::GREATER ) );
758  _result.back().push_back( Constraint<Poly>( lhs, Relation::EQ ) );
759  // Add conjunction (q>0 and r<0 and q^2-r^2*radicand=0) to the substitution result.
760  _result.emplace_back();
761  _result.back().push_back( Constraint<Poly>( _q, Relation::GREATER ) );
762  _result.back().push_back( Constraint<Poly>( _r, Relation::LESS ) );
763  _result.back().push_back( Constraint<Poly>( lhs, Relation::EQ ) );
764  }
765  return true;
766  }
767 
768  template<typename Poly>
769  bool substituteNormalSqrtNeq( const Poly& _radicand,
770  const Poly& _q,
771  const Poly& _r,
772  CaseDistinction<Poly>& _result,
773  bool _accordingPaper )
774  {
775  if( _q.size() > MAX_NUM_OF_TERMS || _r.size() > MAX_NUM_OF_TERMS || _radicand.size() > MAX_NUM_OF_TERMS )
776  return false;
777  Poly lhs = carl::pow(_q, 2) - carl::pow(_r, 2) * _radicand;
778  if( _accordingPaper )
779  {
780  Poly qr = _q * _r;
781  // Add conjunction (q*r>0 and q^2-r^2*radicand!=0) to the substitution result.
782  _result.emplace_back();
783  _result.back().push_back( Constraint<Poly>( qr, Relation::GREATER ) );
784  _result.back().push_back( Constraint<Poly>( lhs, Relation::NEQ ) );
785  }
786  else
787  {
788  // Add conjunction (q>0 and r>0) to the substitution result.
789  _result.emplace_back();
790  _result.back().push_back( Constraint<Poly>( _q, Relation::GREATER ) );
791  _result.back().push_back( Constraint<Poly>( _r, Relation::GREATER ) );
792  // Add conjunction (q<0 and r<0) to the substitution result.
793  _result.emplace_back();
794  _result.back().push_back( Constraint<Poly>( _q, Relation::LESS ) );
795  _result.back().push_back( Constraint<Poly>( _r, Relation::LESS ) );
796  // Add conjunction (q^2-r^2*radicand!=0) to the substitution result.
797  _result.emplace_back();
798  _result.back().push_back( Constraint<Poly>( lhs, Relation::NEQ ) );
799  }
800  return true;
801  }
802 
803  template<typename Poly>
804  bool substituteNormalSqrtLess( const Poly& _radicand,
805  const Poly& _q,
806  const Poly& _r,
807  const Poly& _s,
808  CaseDistinction<Poly>& _result,
809  bool _accordingPaper )
810  {
811  if( _q.size() > MAX_NUM_OF_TERMS || _r.size() > MAX_NUM_OF_TERMS || _radicand.size() > MAX_NUM_OF_TERMS )
812  return false;
813  Poly lhs = carl::pow(_q, 2) - carl::pow(_r, 2) * _radicand;
814  if( _accordingPaper )
815  {
816  Poly qs = _q * _s;
817  Poly rs = _r * _s;
818  // Add conjunction (q*s<0 and q^2-r^2*radicand>0) to the substitution result.
819  _result.emplace_back();
820  _result.back().push_back( Constraint<Poly>( qs, Relation::LESS ) );
821  _result.back().push_back( Constraint<Poly>( lhs, Relation::GREATER ) );
822  // Add conjunction (r*s<=0 and q*s<0) to the substitution result.
823  _result.emplace_back();
824  _result.back().push_back( Constraint<Poly>( rs, Relation::LEQ ) );
825  _result.back().push_back( Constraint<Poly>( qs, Relation::LESS ) );
826  // Add conjunction (r*s<=0 and q^2-r^2*radicand<0) to the substitution result.
827  _result.emplace_back();
828  _result.back().push_back( Constraint<Poly>( rs, Relation::LEQ ) );
829  _result.back().push_back( Constraint<Poly>( lhs, Relation::LESS ) );
830  }
831  else
832  {
833  // Add conjunction (q<0 and s>0 and q^2-r^2*radicand>0) to the substitution result.
834  _result.emplace_back();
835  _result.back().push_back( Constraint<Poly>( _q, Relation::LESS ) );
836  _result.back().push_back( Constraint<Poly>( _s, Relation::GREATER ) );
837  _result.back().push_back( Constraint<Poly>( lhs, Relation::GREATER ) );
838  // Add conjunction (q>0 and s<0 and q^2-r^2*radicand>0) to the substitution result.
839  _result.emplace_back();
840  _result.back().push_back( Constraint<Poly>( _q, Relation::GREATER ) );
841  _result.back().push_back( Constraint<Poly>( _s, Relation::LESS ) );
842  _result.back().push_back( Constraint<Poly>( lhs, Relation::GREATER ) );
843  // Add conjunction (r>0 and s<0 and q^2-r^2*radicand<0) to the substitution result.
844  _result.emplace_back();
845  _result.back().push_back( Constraint<Poly>( _r, Relation::GREATER ) );
846  _result.back().push_back( Constraint<Poly>( _s, Relation::LESS ) );
847  _result.back().push_back( Constraint<Poly>( lhs, Relation::LESS ) );
848  // Add conjunction (r<0 and s>0 and q^2-r^2*radicand<0) to the substitution result.
849  _result.emplace_back();
850  _result.back().push_back( Constraint<Poly>( _r, Relation::LESS ) );
851  _result.back().push_back( Constraint<Poly>( _s, Relation::GREATER ) );
852  _result.back().push_back( Constraint<Poly>( lhs, Relation::LESS ) );
853  // Add conjunction (r>=0 and q<0 and s>0) to the substitution result.
854  _result.emplace_back();
855  _result.back().push_back( Constraint<Poly>( _r, Relation::GEQ ) );
856  _result.back().push_back( Constraint<Poly>( _q, Relation::GREATER ) );
857  _result.back().push_back( Constraint<Poly>( _s, Relation::LESS ) );
858  // Add conjunction (r<=0 and q>0 and s<0) to the substitution result.
859  _result.emplace_back();
860  _result.back().push_back( Constraint<Poly>( _r, Relation::LEQ ) );
861  _result.back().push_back( Constraint<Poly>( _q, Relation::LESS ) );
862  _result.back().push_back( Constraint<Poly>( _s, Relation::GREATER ) );
863  }
864  return true;
865  }
866 
867  template<typename Poly>
868  bool substituteNormalSqrtLeq( const Poly& _radicand,
869  const Poly& _q,
870  const Poly& _r,
871  const Poly& _s,
872  CaseDistinction<Poly>& _result,
873  bool _accordingPaper )
874  {
875  if( _q.size() > MAX_NUM_OF_TERMS || _r.size() > MAX_NUM_OF_TERMS || _radicand.size() > MAX_NUM_OF_TERMS )
876  return false;
877  Poly lhs = carl::pow(_q, 2) - carl::pow(_r, 2) * _radicand;
878  if( _accordingPaper )
879  {
880  Poly qs = _q * _s;
881  Poly rs = _r * _s;
882  // Add conjunction (q*s<=0 and q^2-r^2*radicand>=0) to the substitution result.
883  _result.emplace_back();
884  _result.back().push_back( Constraint<Poly>( qs, Relation::LEQ ) );
885  _result.back().push_back( Constraint<Poly>( lhs, Relation::GEQ ) );
886  // Add conjunction (r*s<=0 and q^2-r^2*radicand<=0) to the substitution result.
887  _result.emplace_back();
888  _result.back().push_back( Constraint<Poly>( rs, Relation::LEQ ) );
889  _result.back().push_back( Constraint<Poly>( lhs, Relation::LEQ ) );
890  }
891  else
892  {
893  // Add conjunction (q<0 and s>0 and q^2-r^2*radicand>=0) to the substitution result.
894  _result.emplace_back();
895  _result.back().push_back( Constraint<Poly>( _q, Relation::LESS ) );
896  _result.back().push_back( Constraint<Poly>( _s, Relation::GREATER ) );
897  _result.back().push_back( Constraint<Poly>( lhs, Relation::GEQ ) );
898  // Add conjunction (q>0 and s<0 and q^2-r^2*radicand>=0) to the substitution result.
899  _result.emplace_back();
900  _result.back().push_back( Constraint<Poly>( _q, Relation::GREATER ) );
901  _result.back().push_back( Constraint<Poly>( _s, Relation::LESS ) );
902  _result.back().push_back( Constraint<Poly>( lhs, Relation::GEQ ) );
903  // Add conjunction (r>0 and s<0 and q^2-r^2*radicand<=0) to the substitution result.
904  _result.emplace_back();
905  _result.back().push_back( Constraint<Poly>( _r, Relation::GREATER ) );
906  _result.back().push_back( Constraint<Poly>( _s, Relation::LESS ) );
907  _result.back().push_back( Constraint<Poly>( lhs, Relation::LEQ ) );
908  // Add conjunction (r<0 and s>0 and q^2-r^2*radicand<=0) to the substitution result.
909  _result.emplace_back();
910  _result.back().push_back( Constraint<Poly>( _r, Relation::LESS ) );
911  _result.back().push_back( Constraint<Poly>( _s, Relation::GREATER ) );
912  _result.back().push_back( Constraint<Poly>( lhs, Relation::LEQ ) );
913  // Add conjunction (r=0 and q=0) to the substitution result.
914  _result.emplace_back();
915  _result.back().push_back( Constraint<Poly>( _r, Relation::EQ ) );
916  _result.back().push_back( Constraint<Poly>( _q, Relation::EQ ) );
917  // Add conjunction (radicand=0 and q=0) to the substitution result.
918  _result.emplace_back();
919  _result.back().push_back( Constraint<Poly>( _radicand, Relation::EQ ) );
920  _result.back().push_back( Constraint<Poly>( _q, Relation::EQ ) );
921  }
922  return true;
923  }
924 
925  template<typename Poly>
926  bool substitutePlusEps( const Constraint<Poly>& _cons,
927  const Substitution<Poly>& _subs,
928  CaseDistinction<Poly>& _result,
929  bool _accordingPaper,
930  Variables& _conflictingVariables,
931  const detail::EvalDoubleIntervalMap& _solutionSpace )
932  {
933  bool result = true;
934  auto vars = variables(_cons);
935  if( !vars.empty() )
936  {
937  if( vars.has( _subs.variable() ) )
938  {
939  switch( _cons.relation() )
940  {
941  case Relation::EQ:
942  {
943  substituteTrivialCase( _cons, _subs, _result );
944  break;
945  }
946  case Relation::NEQ:
947  {
948  substituteNotTrivialCase( _cons, _subs, _result );
949  break;
950  }
951  case Relation::LESS:
952  {
953  result = substituteEpsGradients( _cons, _subs, Relation::LESS, _result, _accordingPaper, _conflictingVariables, _solutionSpace );
954  break;
955  }
956  case Relation::GREATER:
957  {
958  result = substituteEpsGradients( _cons, _subs, Relation::GREATER, _result, _accordingPaper, _conflictingVariables, _solutionSpace );
959  break;
960  }
961  case Relation::LEQ:
962  {
963  substituteTrivialCase( _cons, _subs, _result );
964  result = substituteEpsGradients( _cons, _subs, Relation::LESS, _result, _accordingPaper, _conflictingVariables, _solutionSpace );
965  break;
966  }
967  case Relation::GEQ:
968  {
969  substituteTrivialCase( _cons, _subs, _result );
970  result = substituteEpsGradients( _cons, _subs, Relation::GREATER, _result, _accordingPaper, _conflictingVariables, _solutionSpace );
971  break;
972  }
973  default:
974  assert( false );
975  }
976  simplify( _result, _conflictingVariables, _solutionSpace );
977  }
978  else
979  {
980  _result.emplace_back();
981  _result.back().push_back( _cons );
982  }
983  }
984  else
985  {
986  assert( false );
987  std::cerr << "Warning in substitutePlusEps: The chosen constraint has no variable" << std::endl;
988  }
989  return result;
990  }
991 
992  template<typename Poly>
993  bool substituteEpsGradients( const Constraint<Poly>& _cons,
994  const Substitution<Poly>& _subs,
995  const Relation _relation,
996  CaseDistinction<Poly>& _result,
997  bool _accordingPaper,
998  Variables& _conflictingVariables,
999  const detail::EvalDoubleIntervalMap& _solutionSpace )
1000  {
1001  assert( _cons.variables().has( _subs.variable() ) );
1002  // Create a substitution formed by the given one without an addition of epsilon.
1003  auto term = Term<Poly>::normal(_subs.term().sqrt_ex());
1004  // Call the method substituteNormal with the constraint f(x)~0 and the substitution [x -> t], where the parameter relation is ~.
1005  Constraint<Poly> firstCaseInequality = Constraint<Poly>( _cons.lhs(), _relation );
1006  if( !substituteNormal( firstCaseInequality, {_subs.variable(), term}, _result, _accordingPaper, _conflictingVariables, _solutionSpace ) )
1007  return false;
1008  // Create a vector to store the results of each single substitution.
1009  std::vector<CaseDistinction<Poly>> substitutionResultsVector;
1010  /*
1011  * Let k be the maximum degree of x in f, then call for every 1<=i<=k substituteNormal with:
1012  *
1013  * f^0(x)=0 and ... and f^i-1(x)=0 and f^i(x)~0 as constraints and
1014  * [x -> t] as substitution,
1015  *
1016  * where the relation is ~.
1017  */
1018  Poly deriv = Poly( _cons.lhs() );
1019  while( deriv.has( _subs.variable() ) )
1020  {
1021  // Change the relation symbol of the last added constraint to "=".
1022  Constraint<Poly> equation = Constraint<Poly>( deriv, Relation::EQ );
1023  // Form the derivate of the left hand side of the last added constraint.
1024  deriv = carl::derivative(deriv, _subs.variable(), 1);
1025  // Add the constraint f^i(x)~0.
1026  Constraint<Poly> inequality = Constraint<Poly>( deriv, _relation );
1027  // Apply the substitution (without epsilon) to the new constraints.
1028  substitutionResultsVector.emplace_back();
1029  if( !substituteNormal( equation, {_subs.variable(), term}, substitutionResultsVector.back(), _accordingPaper, _conflictingVariables, _solutionSpace ) )
1030  return false;
1031  substitutionResultsVector.emplace_back();
1032  if( !substituteNormal( inequality, {_subs.variable(), term}, substitutionResultsVector.back(), _accordingPaper, _conflictingVariables, _solutionSpace ) )
1033  return false;
1034  if( !combine( substitutionResultsVector, _result ) )
1035  return false;
1036  simplify( _result, _conflictingVariables, _solutionSpace );
1037  // Remove the last substitution result.
1038  substitutionResultsVector.pop_back();
1039  }
1040  return true;
1041  }
1042 
1043  template<typename Poly>
1044  void substituteInf( const Constraint<Poly>& _cons, const Substitution<Poly>& _subs, CaseDistinction<Poly>& _result, Variables& _conflictingVariables, const detail::EvalDoubleIntervalMap& _solutionSpace )
1045  {
1046  auto vars = variables(_cons);
1047  if( !vars.empty() )
1048  {
1049  if( vars.has( _subs.variable() ) )
1050  {
1051  if( _cons.relation() == Relation::EQ )
1052  substituteTrivialCase( _cons, _subs, _result );
1053  else if( _cons.relation() == Relation::NEQ )
1054  substituteNotTrivialCase( _cons, _subs, _result );
1055  else
1056 
1057  substituteInfLessGreater( _cons, _subs, _result );
1058  simplify( _result, _conflictingVariables, _solutionSpace );
1059  }
1060  else
1061  {
1062  _result.emplace_back();
1063  _result.back().push_back( _cons );
1064  }
1065  }
1066  else
1067  std::cout << "Warning in substituteInf: The chosen constraint has no variable" << std::endl;
1068  }
1069 
1070  template<typename Poly>
1071  void substituteInfLessGreater( const Constraint<Poly>& _cons, const Substitution<Poly>& _subs, CaseDistinction<Poly>& _result )
1072  {
1073  assert( _cons.relation() != Relation::EQ );
1074  assert( _cons.relation() != Relation::NEQ );
1075  // Determine the relation for the coefficients of the odd and even degrees.
1076  Relation oddRelationType = Relation::GREATER;
1077  Relation evenRelationType = Relation::LESS;
1078  if( _subs.term().is_minus_infty())
1079  {
1080  if( _cons.relation() == Relation::GREATER || _cons.relation() == Relation::GEQ )
1081  {
1082  oddRelationType = Relation::LESS;
1083  evenRelationType = Relation::GREATER;
1084  }
1085  }
1086  else
1087  {
1088  assert( _subs.term().is_plus_infty() );
1089  if( _cons.relation() == Relation::LESS || _cons.relation() == Relation::LEQ )
1090  {
1091  oddRelationType = Relation::LESS;
1092  evenRelationType = Relation::GREATER;
1093  }
1094  }
1095  // Check all cases according to the substitution rules.
1096  carl::uint varDegree = _cons.maxDegree( _subs.variable() );
1097  assert( varDegree > 0 );
1098  for( carl::uint i = varDegree + 1; i > 0; --i )
1099  {
1100  // Add conjunction (a_n=0 and ... and a_i~0) to the substitution result.
1101  _result.emplace_back();
1102  for( carl::uint j = varDegree; j > i - 1; --j )
1103  _result.back().push_back( Constraint<Poly>( _cons.coefficient( _subs.variable(), j ), Relation::EQ ) );
1104  if( i > 1 )
1105  {
1106  if( fmod( i - 1, 2.0 ) != 0.0 )
1107  _result.back().push_back( Constraint<Poly>( _cons.coefficient( _subs.variable(), i - 1 ), oddRelationType ) );
1108  else
1109  _result.back().push_back( Constraint<Poly>( _cons.coefficient( _subs.variable(), i - 1 ), evenRelationType ) );
1110  }
1111  else
1112  _result.back().push_back( Constraint<Poly>( _cons.coefficient( _subs.variable(), i - 1 ), _cons.relation() ) );
1113  }
1114  }
1115 
1116  template<typename Poly>
1117  void substituteTrivialCase( const Constraint<Poly>& _cons, const Substitution<Poly>& _subs, CaseDistinction<Poly>& _result )
1118  {
1119  assert( _cons.relation() == Relation::EQ || _cons.relation() == Relation::LEQ || _cons.relation() == Relation::GEQ );
1120  carl::uint varDegree = _cons.maxDegree( _subs.variable() );
1121  // Check the cases (a_0=0 and ... and a_n=0)
1122  _result.emplace_back();
1123  for( carl::uint i = 0; i <= varDegree; ++i )
1124  _result.back().push_back( Constraint<Poly>( _cons.coefficient( _subs.variable(), i ), Relation::EQ ) );
1125  }
1126 
1127  template<typename Poly>
1128  void substituteNotTrivialCase( const Constraint<Poly>& _cons, const Substitution<Poly>& _subs, CaseDistinction<Poly>& _result )
1129  {
1130  assert( _cons.relation() == Relation::NEQ );
1131  carl::uint varDegree = _cons.maxDegree( _subs.variable() );
1132  for( carl::uint i = 0; i <= varDegree; ++i )
1133  {
1134  // Add conjunction (a_i!=0) to the substitution result.
1135  _result.emplace_back();
1136  _result.back().push_back( Constraint<Poly>( _cons.coefficient( _subs.variable(), i ), Relation::NEQ ) );
1137  }
1138  }
1139 } // end namspace vs