1 /*++
2 Copyright (c) 2006 Microsoft Corporation
3
4 Module Name:
5
6 theory_arith_nl.h
7
8 Abstract:
9
10 <abstract>
11
12 Author:
13
14 Leonardo de Moura (leonardo) 2008-12-08.
15
16 Revision History:
17
18 --*/
19 #pragma once
20
21 #include "ast/ast_smt2_pp.h"
22
23 namespace smt {
24
25 template<typename Ext>
mk_nary_mul(unsigned sz,expr * const * args,bool is_int)26 expr * theory_arith<Ext>::mk_nary_mul(unsigned sz, expr * const * args, bool is_int) {
27 if (sz == 0)
28 return m_util.mk_numeral(rational(1), is_int);
29 if (sz == 1)
30 return args[0];
31 if (sz == 2)
32 return m_util.mk_mul(args[0], args[1]);
33 if (m_util.is_numeral(args[0]))
34 return m_util.mk_mul(args[0], m_util.mk_mul(sz - 1, args + 1));
35 return m_util.mk_mul(sz, args);
36 }
37
38 template<typename Ext>
mk_nary_add(unsigned sz,expr * const * args,bool is_int)39 expr * theory_arith<Ext>::mk_nary_add(unsigned sz, expr * const * args, bool is_int) {
40 if (sz == 0)
41 return m_util.mk_numeral(rational(0), is_int);
42 if (sz == 1)
43 return args[0];
44 return m_util.mk_add(sz, args);
45 }
46
47 template<typename Ext>
mk_nary_add(unsigned sz,expr * const * args)48 expr * theory_arith<Ext>::mk_nary_add(unsigned sz, expr * const * args) {
49 SASSERT(sz != 0);
50 return mk_nary_add(sz, args, false);
51 }
52
53 /**
54 \brief Insert v into vars and already_found if v is not already in already_found.
55 */
56 template<typename Ext>
mark_var(theory_var v,svector<theory_var> & vars,var_set & already_found)57 void theory_arith<Ext>::mark_var(theory_var v, svector<theory_var> & vars, var_set & already_found) {
58 if (already_found.contains(v))
59 return;
60 already_found.insert(v);
61 vars.push_back(v);
62 }
63
64 /**
65 \brief Invoke mark_var for all variables in rows that contain v.
66 */
67 template<typename Ext>
mark_dependents(theory_var v,svector<theory_var> & vars,var_set & already_found,row_set & already_visited_rows)68 void theory_arith<Ext>::mark_dependents(theory_var v, svector<theory_var> & vars, var_set & already_found, row_set & already_visited_rows) {
69 if (is_pure_monomial(v)) {
70 expr * n = var2expr(v);
71 SASSERT(m_util.is_mul(n));
72 for (expr * curr : *to_app(n)) {
73 if (ctx.e_internalized(curr)) {
74 theory_var v = expr2var(curr);
75 SASSERT(v != null_theory_var);
76 mark_var(v, vars, already_found);
77 }
78 }
79 }
80 if (is_fixed(v))
81 return;
82 column & c = m_columns[v];
83 typename svector<col_entry>::iterator it = c.begin_entries();
84 typename svector<col_entry>::iterator end = c.end_entries();
85 for (; it != end; ++it) {
86 if (it->is_dead() || already_visited_rows.contains(it->m_row_id))
87 continue;
88 TRACE("non_linear_bug", tout << "visiting row: " << it->m_row_id << "\n";);
89 already_visited_rows.insert(it->m_row_id);
90 row & r = m_rows[it->m_row_id];
91 theory_var s = r.get_base_var();
92 // ignore quasi base vars... actually they should not be used if the problem is non linear...
93 if (is_quasi_base(s))
94 continue;
95 // If s is a base variable different from v and it is free, then this row can be ignored.
96 // It doesn't need to be part of the non linear cluster. For all purposes, this variable
97 // was eliminated by substitution.
98 if (s != null_theory_var && is_free(s) && s != v)
99 continue;
100 typename vector<row_entry>::const_iterator it2 = r.begin_entries();
101 typename vector<row_entry>::const_iterator end2 = r.end_entries();
102 for (; it2 != end2; ++it2) {
103 if (!it2->is_dead() && !is_fixed(it2->m_var))
104 mark_var(it2->m_var, vars, already_found);
105 if (!it2->is_dead() && is_fixed(it2->m_var)) {
106 TRACE("non_linear", tout << "skipped fixed\n";);
107 }
108 }
109 }
110 }
111
112 /**
113 \brief Store in vars the variables that are in the non linear cluster of constraints,
114 and are not satisfied by the current assignment.
115 */
116 template<typename Ext>
get_non_linear_cluster(svector<theory_var> & vars)117 void theory_arith<Ext>::get_non_linear_cluster(svector<theory_var> & vars) {
118 if (m_nl_monomials.empty())
119 return;
120 var_set already_found;
121 row_set already_visited_rows;
122 for (theory_var v : m_nl_monomials) {
123 expr * n = var2expr(v);
124 if (ctx.is_relevant(n))
125 mark_var(v, vars, already_found);
126 }
127 // NB: vars changes inside of loop
128 for (unsigned idx = 0; idx < vars.size(); ++idx) {
129 theory_var v = vars[idx];
130 TRACE("non_linear", tout << "marking dependents of: v" << v << "\n";);
131 mark_dependents(v, vars, already_found, already_visited_rows);
132 }
133 TRACE("non_linear", tout << "variables in non linear cluster:\n";
134 for (theory_var v : vars) tout << "v" << v << " ";
135 tout << "\n";);
136 }
137
138
139 /**
140 \brief Return the number of variables that
141 do not have bounds associated with it.
142 The result is 0, 1, or 2. The value 2 means "2 or more".
143 The second value is the idx of the variable that does not
144 have bounds associated with it. It is only useful when the first value is 1.
145 The second value is -1 if such variable does not exist, that is, the first
146 value is 0.
147
148 \remark if a variables has an even number of occurrences, then
149 I consider that it has a bound associated with it.
150
151 Examples:
152 1) Assume x1, x4 have bounds:
153 analyze_monomial(x1 * x2 * x2 * x3 * x3 * x3 * x4)
154 -->
155 (1,2)
156 Explanation: x2 doesn't have bounds, but x2 has an even power.
157 So x2*x2 has bound [0, oo). So, there is one variable without bounds x3.
158 It is the third variable in the monomial, then its idx is 2.
159 */
160
161 template<typename Ext>
analyze_monomial(expr * m)162 typename theory_arith<Ext>::n_var_power_pair theory_arith<Ext>::analyze_monomial(expr * m) const {
163
164 buffer<var_power_pair> vp;
165 decompose_monomial(m, vp);
166 unsigned c = 0;
167 var_power_pair q(nullptr, 0);
168 for (auto const& p : vp) {
169 if (p.second % 2 == 1 && is_free(p.first)) {
170 c++;
171 q = p;
172 if (c > 1)
173 break;
174 }
175 }
176 return std::make_pair(c, q);
177 }
178
179 template<typename Ext>
decompose_monomial(expr * m,buffer<var_power_pair> & vp)180 rational theory_arith<Ext>::decompose_monomial(expr* m, buffer<var_power_pair>& vp) const {
181 rational coeff(1);
182 vp.reset();
183 expr_fast_mark1 mark;
184 auto insert = [&](expr* arg) {
185 rational r;
186 if (m_util.is_numeral(arg, r))
187 coeff *= r;
188 else if (mark.is_marked(arg)) {
189 for (unsigned i = vp.size(); i-- > 0; ) {
190 if (vp[i].first == arg) {
191 vp[i].second++;
192 break;
193 }
194 }
195 }
196 else {
197 mark.mark(arg);
198 vp.push_back(var_power_pair(arg, 1));
199 }
200 };
201 while (m_util.is_mul(m)) {
202 unsigned sz = to_app(m)->get_num_args();
203 for (unsigned i = 0; i + 1 < sz; ++i) {
204 insert(to_app(m)->get_arg(i));
205 }
206 m = to_app(m)->get_arg(sz - 1);
207 }
208 insert(m);
209 return coeff;
210 }
211
212
213 /**
214 \brief Return an interval using the bounds for v.
215 */
216 template<typename Ext>
mk_interval_for(theory_var v)217 interval theory_arith<Ext>::mk_interval_for(theory_var v) {
218 bound * l = lower(v);
219 bound * u = upper(v);
220 if (l && u) {
221 // optimization may introduce non-standard bounds.
222 if (l->get_value() == u->get_value() && !l->get_value().get_infinitesimal().to_rational().is_zero()) {
223 return interval(m_dep_manager);
224 }
225 return interval(m_dep_manager,
226 l->get_value().get_rational().to_rational(),
227 l->get_value().get_infinitesimal().to_rational().is_pos(),
228 m_dep_manager.mk_leaf(l),
229 u->get_value().get_rational().to_rational(),
230 u->get_value().get_infinitesimal().to_rational().is_neg(),
231 m_dep_manager.mk_leaf(u));
232 }
233 else if (l) {
234 return interval(m_dep_manager,
235 l->get_value().get_rational().to_rational(),
236 l->get_value().get_infinitesimal().to_rational().is_pos(),
237 true,
238 m_dep_manager.mk_leaf(l));
239 }
240 else if (u) {
241 return interval(m_dep_manager,
242 u->get_value().get_rational().to_rational(),
243 u->get_value().get_infinitesimal().to_rational().is_neg(),
244 false,
245 m_dep_manager.mk_leaf(u));
246 }
247 else {
248 return interval(m_dep_manager);
249 }
250 }
251
252 /**
253 \brief Return an interval for the given expression using its bounds.
254 */
255 template<typename Ext>
mk_interval_for(expr * n)256 interval theory_arith<Ext>::mk_interval_for(expr * n) {
257 if (!has_var(n))
258 return interval(m_dep_manager);
259 return mk_interval_for(expr2var(n));
260 }
261
262 /**
263 \brief target *= [lower(var), upper(var)]^power
264 */
265 template<typename Ext>
mul_bound_of(expr * var,unsigned power,interval & target)266 void theory_arith<Ext>::mul_bound_of(expr * var, unsigned power, interval & target) {
267 theory_var v = expr2var(var);
268 interval i = mk_interval_for(v);
269
270 TRACE("non_linear",
271 display_interval(tout << "bound: ",i); tout << i << "\n";
272 tout << mk_pp(var, get_manager()) << "\n";
273 tout << "power " << power << ": " << expt(i, power) << "\n";
274 display_interval(tout << "target before: ", target); tout << "\n";);
275
276 i.expt(power);
277 target *= i;
278
279 get_manager().limit().inc((target.is_lower_open() || target.minus_infinity()) ? 1 : target.get_lower_value().bitsize());
280 get_manager().limit().inc((target.is_upper_open() || target.plus_infinity()) ? 1 : target.get_upper_value().bitsize());
281
282 TRACE("non_linear", display_interval(tout << "target after: ", target); tout << "\n";);
283 }
284
285 /**
286 \brief Evaluate the given expression using interval arithmetic.
287
288 - If a subexpression is internalized, then mk_interval_for is used to
289 compute its interval.
290
291 - Only +, *, and numerals are handled.
292 */
293 template<typename Ext>
evaluate_as_interval(expr * n)294 interval theory_arith<Ext>::evaluate_as_interval(expr * n) {
295 expr* arg;
296 rational val;
297 #define TR() TRACE("nl_evaluate", tout << "eval: " << mk_bounded_pp(n, get_manager(), 10) << "\n";\
298 display_nested_form(tout, n); tout << "\ninterval: " << r << "\n";);
299
300 if (has_var(n)) {
301 interval r = mk_interval_for(n);
302 TR();
303 return r;
304 }
305 else if (m_util.is_add(n)) {
306 interval r(m_dep_manager, rational(0));
307 for (expr* arg : *to_app(n)) {
308 r += evaluate_as_interval(arg);
309 }
310 TR();
311 return r;
312 }
313 else if (m_util.is_mul(n)) {
314 buffer<var_power_pair> vp;
315 rational coeff = decompose_monomial(n, vp);
316 interval r(m_dep_manager, coeff);
317 for (var_power_pair const& p : vp) {
318 expr * var = p.first;
319 unsigned power = p.second;
320 interval it = evaluate_as_interval(var);
321 it.expt(power);
322 r *= it;
323 }
324 TR();
325 return r;
326 }
327 else if (m_util.is_to_real(n, arg)) {
328 return evaluate_as_interval(arg);
329 }
330 else if (m_util.is_numeral(n, val)) {
331 interval r = interval(m_dep_manager, val);
332 TR();
333 return r;
334 }
335 else {
336 TRACE("nl_evaluate", tout << "is unknown\n";);
337 return interval(m_dep_manager);
338 }
339 }
340
341 template<typename Ext>
display_monomial(std::ostream & out,expr * n)342 void theory_arith<Ext>::display_monomial(std::ostream & out, expr * n) const {
343 bool first = true;
344 buffer<var_power_pair> vp;
345 rational coeff = decompose_monomial(n, vp);
346 if (!coeff.is_one()) {
347 out << coeff;
348 first = false;
349 }
350 for (auto const& p : vp) {
351 SASSERT(p.first != 0);
352 if (first) first = false; else out << " * ";
353 out << mk_bounded_pp(p.first, get_manager()) << "^" << p.second;
354 }
355 }
356
357 template<typename Ext>
dependency2new_bound(v_dependency * dep,derived_bound & new_bound)358 void theory_arith<Ext>::dependency2new_bound(v_dependency * dep, derived_bound& new_bound) {
359 ptr_vector<void> bounds;
360 m_dep_manager.linearize(dep, bounds);
361 m_tmp_lit_set.reset();
362 m_tmp_eq_set.reset();
363 for (void* _b : bounds) {
364 bound * b = static_cast<bound*>(_b);
365 accumulate_justification(*b, new_bound, numeral::zero(), m_tmp_lit_set, m_tmp_eq_set);
366 }
367 }
368
369 /**
370 \brief Create a new derived bound. The justification is stored in the object dep.
371 */
372 template<typename Ext>
mk_derived_nl_bound(theory_var v,inf_numeral const & coeff,bound_kind k,v_dependency * dep)373 void theory_arith<Ext>::mk_derived_nl_bound(theory_var v, inf_numeral const & coeff, bound_kind k, v_dependency * dep) {
374 inf_numeral coeff_norm = normalize_bound(v, coeff, k);
375 derived_bound * new_bound = alloc(derived_bound, v, coeff_norm, k);
376 m_bounds_to_delete.push_back(new_bound);
377 m_asserted_bounds.push_back(new_bound);
378 // copy justification to new bound
379 dependency2new_bound(dep, *new_bound);
380 TRACE("non_linear", new_bound->display(*this, tout); tout << "\n";);
381 }
382
383 /**
384 \brief Update the bounds of v, using the interval i.
385 Return true if i improves the bounds of v.
386 */
387 template<typename Ext>
update_bounds_using_interval(theory_var v,interval const & i)388 bool theory_arith<Ext>::update_bounds_using_interval(theory_var v, interval const & i) {
389 SASSERT(v != null_theory_var);
390 bool r = false;
391 if (!i.minus_infinity()) {
392 inf_numeral new_lower(i.get_lower_value());
393 if (i.is_lower_open()) {
394 if (is_int(v)) {
395 if (new_lower.is_int()) {
396 new_lower += rational::one();
397 }
398 else {
399 new_lower = ceil(new_lower.get_rational());
400 }
401 }
402 else {
403 new_lower += get_epsilon(v);
404 }
405 }
406 bound * old_lower = lower(v);
407 if (old_lower == nullptr || new_lower > old_lower->get_value()) {
408 TRACE("non_linear", tout << "NEW lower bound for v" << v << " " << mk_pp(var2expr(v), get_manager())
409 << " " << new_lower << "\n";
410 display_interval(tout, i); tout << "\n";);
411 mk_derived_nl_bound(v, new_lower, B_LOWER, i.get_lower_dependencies());
412 r = true;
413 }
414 }
415 if (!i.plus_infinity()) {
416 inf_numeral new_upper(i.get_upper_value());
417 if (i.is_upper_open()) {
418 if (is_int(v)) {
419 if (new_upper.is_int()) {
420 new_upper -= rational::one();
421 }
422 else {
423 new_upper = floor(new_upper.get_rational());
424 }
425 }
426 else {
427 new_upper -= get_epsilon(v);
428 }
429 }
430 bound * old_upper = upper(v);
431 if (old_upper == nullptr || new_upper < old_upper->get_value()) {
432 TRACE("non_linear", tout << "NEW upper bound for v" << v << " " << new_upper << "\n";
433 display_interval(tout, i); tout << "\n";);
434 mk_derived_nl_bound(v, new_upper, B_UPPER, i.get_upper_dependencies());
435 r = true;
436 }
437 }
438 return r;
439 }
440
441 template<typename Ext>
update_bounds_using_interval(expr * n,interval const & i)442 bool theory_arith<Ext>::update_bounds_using_interval(expr * n, interval const & i) {
443 SASSERT(expr2var(n) != null_theory_var);
444 TRACE("non_linear", tout << "NL bounds for m: " << i << "\n" << mk_pp(n, get_manager()) << "\n";);
445 return update_bounds_using_interval(expr2var(n), i);
446 }
447
448 /**
449 \brief Use the bounds of the variables to build a bound for m.
450 */
451 template<typename Ext>
propagate_nl_upward(expr * m)452 bool theory_arith<Ext>::propagate_nl_upward(expr * m) {
453 SASSERT(is_pure_monomial(m));
454 TRACE("nl_arith_bug", tout << "processing upward:\n" << mk_pp(m, get_manager()) << "\n";);
455 buffer<var_power_pair> vp;
456 rational coeff = decompose_monomial(m, vp);
457 interval new_bounds(m_dep_manager, coeff);
458 for (var_power_pair const& p : vp) {
459 expr * var = p.first;
460 unsigned power = p.second;
461 TRACE("nl_arith_bug", tout << "interval before: " << new_bounds << "\n";
462 theory_var v = expr2var(var);
463 interval i = mk_interval_for(v);
464 display_var(tout, v);
465 tout << "interval for v" << i << " " << mk_pp(var, get_manager()) << "\npower: " << power << " " << expt(i, power) << "\n";);
466 mul_bound_of(var, power, new_bounds);
467 TRACE("nl_arith_bug", tout << "interval after: " << new_bounds << "\n";);
468 }
469 return update_bounds_using_interval(m, new_bounds);
470 }
471
472 /**
473 \brief Propagate a bound to the i-th variable of the given monomial
474 using the bounds of m and other variables in m.
475
476 \remark We do not support roots in interval... so, if the i-th var has power != 1
477 the method returns without doing anything.
478 */
479 template<typename Ext>
propagate_nl_downward(expr * n,var_power_pair const & p)480 bool theory_arith<Ext>::propagate_nl_downward(expr * n, var_power_pair const& p) {
481 SASSERT(is_pure_monomial(n));
482 expr * v = p.first;
483 unsigned power = p.second;
484 if (power != 1)
485 return false; // TODO: remove, when the n-th root is implemented in interval.
486
487 buffer<var_power_pair> vp;
488 rational coeff = decompose_monomial(n, vp);
489
490 interval other_bounds(m_dep_manager, coeff);
491 // TODO: the following code can be improved it is quadratic on the degree of the monomial.
492 for (var_power_pair const& p : vp) {
493 if (p.first == v)
494 continue;
495 expr * var = p.first;
496 unsigned power = p.second;
497 mul_bound_of(var, power, other_bounds);
498 }
499 if (other_bounds.contains_zero())
500 return false; // interval division requires that divisor doesn't contain 0.
501 interval r = mk_interval_for(n);
502 TRACE("nl_arith_bug", tout << "m: " << mk_ismt2_pp(n, get_manager()) << "\nv: " << mk_ismt2_pp(v, get_manager()) <<
503 "\npower: " << power << "\n";
504 display_interval(tout << "monomial bounds\n", r);
505 display_interval(tout << "other bounds\n", other_bounds);
506 );
507 r /= other_bounds;
508 return update_bounds_using_interval(v, r);
509 }
510
511
512 /**
513 \brief The given monomial and its elements have bounds.
514 Propagate bounds to all of them.
515 Return true if some bound was propagated.
516 */
517 template<typename Ext>
propagate_nl_bounds(expr * m)518 bool theory_arith<Ext>::propagate_nl_bounds(expr * m) {
519 TRACE("non_linear", tout << "propagate several bounds using:\n"; display_monomial(tout, m); tout << "\n";);
520 bool result = propagate_nl_upward(m);
521 buffer<var_power_pair> vp;
522 rational coeff = decompose_monomial(m, vp);
523 for (auto const& p : vp) {
524 if (propagate_nl_downward(m, p)) {
525 m_stats.m_nl_bounds++;
526 result = true;
527 }
528 }
529 return result;
530 }
531
532 /**
533 \brief Try to propagate bounds using non linear monomials.
534 Return true if some bound was propagated.
535 */
536 template<typename Ext>
propagate_nl_bounds()537 bool theory_arith<Ext>::propagate_nl_bounds() {
538 m_dep_manager.reset();
539 bool propagated = false;
540 for (unsigned i = 0; i < m_nl_monomials.size(); i++) {
541 theory_var v = m_nl_monomials[i];
542 expr * m = var2expr(v);
543 if (!ctx.is_relevant(m))
544 continue;
545 auto p = analyze_monomial(m);
546 TRACE("propagate_nl_bound", tout << "m: " << mk_ismt2_pp(m, get_manager()) << "\n" << "p: " << p.first << "\n";);
547 unsigned num_bad_vars = p.first;
548 var_power_pair q = p.second;
549 SASSERT(num_bad_vars != 1 || q.first != nullptr);
550 if (num_bad_vars >= 2)
551 continue;
552 bool is_free_m = is_free(m);
553 TRACE("propagate_nl_bound", tout << "is_free_m: " << is_free_m << "\n";);
554 if (num_bad_vars == 1 && is_free_m)
555 continue;
556 if (num_bad_vars == 0) {
557 if (!is_free_m) {
558 if (propagate_nl_bounds(m))
559 propagated = true;
560 }
561 else {
562 if (propagate_nl_upward(m)) {
563 m_stats.m_nl_bounds++;
564 propagated = true;
565 }
566 }
567 }
568 else {
569 SASSERT (!is_free_m);
570 if (propagate_nl_downward(m, q)) {
571 m_stats.m_nl_bounds++;
572 propagated = true;
573 }
574 }
575 }
576 return propagated;
577 }
578
579 /**
580 \brief Return the value of v as a rational. If computed_epsilon = false and v has an infinitesimal, then
581 compute_epsilon() is invoked.
582 */
583 template<typename Ext>
get_value(theory_var v,bool & computed_epsilon)584 rational theory_arith<Ext>::get_value(theory_var v, bool & computed_epsilon) {
585 inf_numeral const & val = get_value(v);
586 if (!val.get_infinitesimal().is_zero() && !computed_epsilon) {
587 compute_epsilon();
588 refine_epsilon();
589 computed_epsilon = true;
590 m_model_depends_on_computed_epsilon = true;
591 }
592 return val.get_rational().to_rational() + m_epsilon.to_rational() * val.get_infinitesimal().to_rational();
593 }
594
595 /**
596 \brief Return true if for the monomial x_1 * ... * x_n associated with v,
597 the following holds:
598
599 get_value(x_1) * ... * get_value(x_n) = get_value(v)
600 */
601 template<typename Ext>
check_monomial_assignment(theory_var v,bool & computed_epsilon)602 bool theory_arith<Ext>::check_monomial_assignment(theory_var v, bool & computed_epsilon) {
603 SASSERT(is_pure_monomial(var2expr(v)));
604 expr * m = var2expr(v);
605 rational val(1), v_val;
606 for (expr* arg : *to_app(m)) {
607 theory_var curr = expr2var(arg);
608 SASSERT(curr != null_theory_var);
609 v_val = get_value(curr, computed_epsilon);
610 TRACE("non_linear", tout << mk_pp(arg, get_manager()) << " = " << v_val << "\n";);
611 val *= v_val;
612 }
613 v_val = get_value(v, computed_epsilon);
614 TRACE("non_linear", tout << "v" << v << " := " << v_val << " == " << val << "\n";);
615 return v_val == val;
616 }
617
618
619 /**
620 \brief Return true if for every monomial x_1 * ... * x_n,
621 get_value(x_1) * ... * get_value(x_n) = get_value(x_1 * ... * x_n)
622 */
623 template<typename Ext>
check_monomial_assignments()624 bool theory_arith<Ext>::check_monomial_assignments() {
625 bool computed_epsilon = false;
626 for (theory_var v : m_nl_monomials) {
627 TRACE("non_linear", tout << "v" << v << " is relevant: " << ctx.is_relevant(get_enode(v)) << "\n";
628 tout << "check_monomial_assignments result: " << check_monomial_assignment(v, computed_epsilon) << "\n";
629 tout << "computed_epsilon: " << computed_epsilon << "\n";);
630 if (ctx.is_relevant(get_enode(v)) && !check_monomial_assignment(v, computed_epsilon)) {
631 TRACE("non_linear_failed", tout << "check_monomial_assignment failed for:\n" << mk_ismt2_pp(var2expr(v), get_manager()) << "\n";
632 display_var(tout, v););
633 return false;
634 }
635 }
636 return true;
637 }
638
639 /**
640 \brief Try to find an integer variable for performing branching
641 in the non linear cluster.
642
643 The idea is select a variable in a monomial with an invalid
644 assignment. I give preference to variables with small ranges.
645 If no variable is bounded, then select a random one.
646
647 Free variables are not considered.
648 */
649 template<typename Ext>
find_nl_var_for_branching()650 theory_var theory_arith<Ext>::find_nl_var_for_branching() {
651 TRACE("nl_branching", tout << "looking for variable to branch...\n"; display(tout););
652 theory_var target = null_theory_var;
653 bool bounded = false;
654 unsigned n = 0;
655 numeral range;
656 for (unsigned j = 0; j < m_nl_monomials.size(); ++j) {
657 theory_var v = m_nl_monomials[j];
658 if (is_real(v))
659 continue;
660 bool computed_epsilon = false;
661 bool r = check_monomial_assignment(v, computed_epsilon);
662 if (!r) {
663 expr * m = get_enode(v)->get_expr();
664 SASSERT(is_pure_monomial(m));
665 for (expr * arg : *to_app(m)) {
666 theory_var curr = ctx.get_enode(arg)->get_th_var(get_id());
667 TRACE("nl_branching", tout << "target: v" << target << ", curr: v" << curr << "\n";);
668 if (!is_fixed(curr) && is_int(curr)) {
669 if (is_bounded(curr)) {
670 numeral new_range;
671 new_range = upper_bound(curr).get_rational();
672 new_range -= lower_bound(curr).get_rational();
673 if (!bounded || new_range < range) {
674 target = curr;
675 range = new_range;
676 bounded = true;
677 }
678 }
679 else if (!bounded) {
680 n++;
681 TRACE("nl_branching", tout << "n: " << n << "\n";);
682 if (m_random()%n == 0)
683 target = curr;
684 SASSERT(target != null_theory_var);
685 }
686 SASSERT(target != null_theory_var);
687 }
688 TRACE("nl_branching", tout << "after target: v" << target << "\n";);
689 }
690 }
691 }
692 return target;
693 }
694
695 /**
696 \brief Branch on an integer variable. This method is invoked when v is part
697 of a non linear monomial that is not satisfied by the current assignment.
698 if v >= l, then create the case split v >= l+1
699 else v <= u, then create the case split v <= u-1
700 else create the bound v = 0 and case split on it.
701 */
702 template<typename Ext>
branch_nl_int_var(theory_var v)703 bool theory_arith<Ext>::branch_nl_int_var(theory_var v) {
704 TRACE("non_linear", tout << "BRANCHING on v" << v << "\n";);
705 m_stats.m_nl_branching++;
706 SASSERT(is_int(v));
707 expr_ref bound(get_manager());
708 if (lower(v))
709 bound = m_util.mk_le(var2expr(v), m_util.mk_numeral(lower_bound(v).get_rational().to_rational(), true));
710 else if (upper(v))
711 bound = m_util.mk_ge(var2expr(v), m_util.mk_numeral(upper_bound(v).get_rational().to_rational(), true));
712 else
713 bound = m_util.mk_eq(var2expr(v), m_util.mk_numeral(rational(0), true));
714 TRACE("non_linear", tout << "new bound:\n" << mk_pp(bound, get_manager()) << "\n";);
715 ast_manager & m = get_manager();
716 {
717 std::function<expr*(void)> fn = [&]() { return m.mk_or(bound, m.mk_not(bound)); };
718 scoped_trace_stream _sts(*this, fn);
719 ctx.internalize(bound, true);
720 }
721 ctx.mark_as_relevant(bound.get());
722 literal l = ctx.get_literal(bound);
723 SASSERT(!l.sign());
724 ctx.set_true_first_flag(l.var()); // force the context to case split to true first, independently of the phase selection strategy.
725 return true;
726 }
727
728 /**
729 \brief Return true if the given monomial is linear.
730 */
731 template<typename Ext>
is_monomial_linear(expr * m)732 bool theory_arith<Ext>::is_monomial_linear(expr * m) const {
733 SASSERT(is_pure_monomial(m));
734 unsigned num_nl_vars = 0;
735 for (expr* arg : *to_app(m)) {
736 if (!ctx.e_internalized(arg))
737 return false;
738 theory_var _var = expr2var(arg);
739 CTRACE("non_linear", is_fixed(_var),
740 tout << mk_pp(arg, get_manager()) << " is fixed: " << lower_bound(_var) << "\n";);
741 if (!is_fixed(_var)) {
742 num_nl_vars++;
743 }
744 else if (lower_bound(_var).is_zero()) {
745 return true;
746 }
747 }
748 return num_nl_vars <= 1;
749 }
750
751 /**
752 \brief Return the product of the value of the fixed variables in the
753 monomial m.
754 */
755 template<typename Ext>
get_monomial_fixed_var_product(expr * m)756 typename theory_arith<Ext>::numeral theory_arith<Ext>::get_monomial_fixed_var_product(expr * m) const {
757 SASSERT(is_pure_monomial(m));
758 numeral r(1);
759 for (expr * arg : *to_app(m)) {
760 theory_var _var = expr2var(arg);
761 if (is_fixed(_var))
762 r *= lower_bound(_var).get_rational();
763 }
764 TRACE("arith", tout << mk_pp(m, get_manager()) << " " << r << "\n";);
765 return r;
766 }
767
768 /**
769 \brief Return the first non fixed variable in the given monomial.
770 Return 0, if the monomial does not have a non fixed variable.
771 */
772 template<typename Ext>
get_monomial_non_fixed_var(expr * m)773 expr * theory_arith<Ext>::get_monomial_non_fixed_var(expr * m) const {
774 SASSERT(is_pure_monomial(m));
775 for (unsigned i = 0; i < to_app(m)->get_num_args(); i++) {
776 expr * arg = to_app(m)->get_arg(i);
777 theory_var _var = expr2var(arg);
778 if (!is_fixed(_var))
779 return arg;
780 }
781 return nullptr;
782 }
783
784 /**
785 \brief Propagate linear monomial. Check whether the give
786 monomial became linear and propagate.
787 */
788 template<typename Ext>
propagate_linear_monomial(theory_var v)789 bool theory_arith<Ext>::propagate_linear_monomial(theory_var v) {
790 TRACE("non_linear", tout << "checking whether v" << v << " became linear...\n";);
791 if (m_data[v].m_nl_propagated)
792 return false; // already propagated this monomial.
793 expr * m = var2expr(v);
794 if (!is_monomial_linear(m))
795 return false; // monomial is not linear.
796
797 m_stats.m_nl_linear++;
798
799 m_data[v].m_nl_propagated = true;
800 m_nl_propagated.push_back(v);
801 TRACE("non_linear", tout << "v" << v << " is linear " << mk_pp(m, get_manager()) << "\n";);
802
803
804 numeral k = get_monomial_fixed_var_product(m);
805 TRACE("non_linear", tout << "new linear monomial... k: " << k << "\n";);
806 expr * x_n = k.is_zero() ? nullptr : get_monomial_non_fixed_var(m);
807 TRACE("non_linear_bug", if (x_n != 0) { tout << "x_n: " << mk_bounded_pp(x_n, get_manager()) << "\nx_n: #" << x_n->get_id() << "\n"; });
808 derived_bound * new_lower = nullptr;
809 derived_bound * new_upper = nullptr;
810 if (x_n != nullptr) {
811 // All but one of the x_i variables are assigned.
812 // Let x_n be the unassigned variable.
813 // Then, we know that x_1*...*x_n = k*x_n, where k is the product of beta(x_1)*...*beta(x_{n-1})
814 // beta(x_i) == lower(x_i)
815
816 // Let m be (* x_1 ... x_n), then assert equality
817 // (= (+ (* x_1 ... x_n) (* -k x_n)) 0) when x_1 ... x_{n-1} are fixed variables.
818 // where k = lower(x_1)*...*lower(x_{n-1})
819 TRACE("non_linear", tout << "x_n: " << mk_pp(x_n, get_manager()) << "\n";);
820 k.neg();
821 expr * k_x_n = k.is_one() ? x_n : m_util.mk_mul(m_util.mk_numeral(k.to_rational(), is_int(v)), x_n);
822 expr * rhs = m_util.mk_add(m, k_x_n);
823 TRACE("non_linear_bug", tout << "rhs: " << mk_bounded_pp(rhs, get_manager(),5) << "\ninternalized: " << ctx.e_internalized(rhs) << "\n";);
824 if (!has_var(rhs)) {
825 ctx.internalize(rhs, false);
826 ctx.mark_as_relevant(rhs);
827 }
828 TRACE("non_linear_bug", tout << "enode: " << ctx.get_enode(rhs) << " enode_id: " << ctx.get_enode(rhs)->get_owner_id() << "\n";);
829 theory_var new_v = expr2var(rhs);
830 SASSERT(new_v != null_theory_var);
831 new_lower = alloc(derived_bound, new_v, inf_numeral(0), B_LOWER);
832 new_upper = alloc(derived_bound, new_v, inf_numeral(0), B_UPPER);
833 }
834 else {
835 // One of the x_i variables is zero,
836 // or all of them are assigned.
837
838 // Assert the equality
839 // (= (* x_1 ... x_n) k)
840 TRACE("non_linear", tout << "all variables are fixed, and bound is: " << k << "\n";);
841 new_lower = alloc(derived_bound, v, inf_numeral(k), B_LOWER);
842 new_upper = alloc(derived_bound, v, inf_numeral(k), B_UPPER);
843 }
844 SASSERT(new_lower);
845 SASSERT(new_upper);
846 m_bounds_to_delete.push_back(new_lower);
847 m_asserted_bounds.push_back(new_lower);
848 m_bounds_to_delete.push_back(new_upper);
849 m_asserted_bounds.push_back(new_upper);
850
851 // Add the justification for new_lower and new_upper.
852 // The justification is the lower and upper bounds of all fixed variables.
853 m_tmp_lit_set.reset();
854 m_tmp_eq_set.reset();
855
856 SASSERT(is_pure_monomial(m));
857 bool found_zero = false;
858 for (unsigned i = 0; !found_zero && i < to_app(m)->get_num_args(); i++) {
859 expr * arg = to_app(m)->get_arg(i);
860 theory_var _var = expr2var(arg);
861 if (is_fixed(_var)) {
862 bound * l = lower(_var);
863 bound * u = upper(_var);
864 if (l->get_value().is_zero()) {
865 /* if zero was found, then it is the explanation */
866 SASSERT(k.is_zero());
867 found_zero = true;
868 m_tmp_lit_set.reset();
869 m_tmp_eq_set.reset();
870 new_lower->m_lits.reset();
871 new_lower->m_eqs.reset();
872 }
873 accumulate_justification(*l, *new_lower, numeral::zero(), m_tmp_lit_set, m_tmp_eq_set);
874
875 TRACE("non_linear",
876 for (literal l : new_lower->m_lits) {
877 ctx.display_detailed_literal(tout, l) << " ";
878 }
879 tout << "\n";);
880
881 accumulate_justification(*u, *new_lower, numeral::zero(), m_tmp_lit_set, m_tmp_eq_set);
882
883 TRACE("non_linear",
884 for (literal l : new_lower->m_lits) {
885 ctx.display_detailed_literal(tout, l) << " ";
886 }
887 tout << "\n";);
888
889 }
890 }
891 new_upper->m_lits.append(new_lower->m_lits);
892 new_upper->m_eqs.append(new_lower->m_eqs);
893
894 TRACE("non_linear",
895 new_lower->display(*this, tout << "lower: "); tout << "\n";
896 new_upper->display(*this, tout << "upper: "); tout << "\n";
897 for (literal lit : new_upper->m_lits) {
898 ctx.display_detailed_literal(tout, lit) << " ";
899 }
900 tout << "\n";);
901
902 return true;
903 }
904
905 /**
906 \brief Traverse all non linear monomials, and check the ones that became
907 linear and propagate. Return true if propagated.
908 */
909 template<typename Ext>
propagate_linear_monomials()910 bool theory_arith<Ext>::propagate_linear_monomials() {
911 if (!reflection_enabled())
912 return false;
913 TRACE("non_linear", tout << "propagating linear monomials...\n";);
914 bool p = false;
915 // CMW: m_nl_monomials can grow during this loop, so
916 // don't use iterators.
917 for (unsigned i = 0; i < m_nl_monomials.size(); i++) {
918 if (propagate_linear_monomial(m_nl_monomials[i]))
919 p = true;
920 }
921 CTRACE("non_linear", p, display(tout););
922 return p;
923 }
924
925 /*
926 Interval arithmetic does not satisfy distributivity.
927 Actually, it satisfies the sub-distributivity property:
928
929 x*(y + z) \subseteq x*y + x*z
930
931 The sub-distributivity property only holds if condensation
932 is not used. For example:
933
934 x * (x^3 + 1) \subseteq x*x^3 + x,
935
936 but it is not the case that
937
938 x * (x^3 + 1) \subseteq x^4 + x
939
940 for example, for x = [-2,1]
941
942 x*(x^3+1) = [-7, 14]
943 x^4 + x = [-2, 17]
944
945 This weakness of AI is known as the "dependency problem",
946 which comes from the decorrelation of the multiple occurrences
947 of one variable during interval evaluation.
948
949 Given a polynomial:
950 p(x) = a_0 + a_1 * x + ... + a_n * x^n
951 The horner extension is:
952 h_p(x) = a_0 + x*(a_1 + ... + x*(a_{n-1} + a_n * x) + ...)
953
954 The horner extension of p(x) = x^4 + x^3 + 2*x is:
955 h_p(x) = x(2 + x^3(1 + x))
956
957 The horner extension evaluates tighter intervals when
958 condensation is not used.
959
960 Remark: there is no guarantee that horner extension will
961 provide a tighter interval than a sum of monomials when
962 condensation is used.
963
964 For multivariate polynomials nested (or cross nested) forms
965 are used. The idea is to select one variable, and pretend the
966 other are parameters. The horner form is computed for the selected
967 variable, and the computation continues for the polynomials on the
968 parameters.
969
970 As described above, the horner form is not optimal with respect to
971 to condensation. I use the following two properties to deal with
972 monovariate polynomials with two monomials:
973
974 p(x) = a*x^n + b*x^{n+m} for n >= m
975 is equivalent to
976
977 b*x^{n-m}*[(x^{m} + a/(2b))^2 - (a/2b)^2]
978
979 This polynomial provides tight bound when n and m have the same parity and:
980 1) a*b > 0 and (lower(x) >= 0 or upper(x)^m <= -a/b)
981 2) a*b < 0 and (upper(x) <= 0 or lower(x)^m >= a/b)
982
983 This polynomial also provides tight bounds when n = m,
984 and the polynomial is simplified to, and n and m may have arbitrary parities:
985
986 b*[(x^{n} + a/(2b))^2 - (a/2b)^2]
987
988 Example:
989 x > 1
990 x^2 - x <= 0
991 is unsatisfiable
992
993 If we compute the bounds for x^2 - x we obtain
994 (-oo, oo).
995
996 On the other hand, if we compute the bounds for
997 (x - 1/2)^2 - 1/4
998 we obtain the bounds (0, oo), and the inconsistency
999 is detected.
1000
1001 Remark: In Z3, I condensate multiple occurrences of a variable
1002 when evaluating monomials. So, the interval for a monomial is
1003 always tight.
1004
1005 Remark: M1*(M2 + M3) is more precise than M1 * M2 + M1 * M3,
1006 if intersection(Vars(M1), union(Vars(M2), Vars(M3))) = empty-set,
1007
1008 Remark: A trivial consequence of Moore's theorem for interval
1009 arithmetic. If two monomials M1 and M2 do not share variables,
1010 then the interval for M1 + M2 is tight.
1011 */
1012
1013 /**
1014 \brief Check whether the same variable occurs in two different monomials.
1015
1016 \remark Fixed variables are ignored.
1017
1018 \remark A trivial consequence of Moore's theorem for interval
1019 arithmetic. If two monomials M1 and M2 do not share variables,
1020 then the interval for M1 + M2 is tight.
1021 */
1022 template<typename Ext>
is_problematic_non_linear_row(row const & r)1023 bool theory_arith<Ext>::is_problematic_non_linear_row(row const & r) {
1024 m_tmp_var_set.reset();
1025 typename vector<row_entry>::const_iterator it = r.begin_entries();
1026 typename vector<row_entry>::const_iterator end = r.end_entries();
1027 for (; it != end; ++it) {
1028 if (!it->is_dead()) {
1029 theory_var v = it->m_var;
1030 if (is_fixed(v))
1031 continue;
1032 if (is_pure_monomial(v)) {
1033 expr * m = var2expr(v);
1034 for (expr* arg : *to_app(m)) {
1035 theory_var curr = expr2var(arg);
1036 if (m_tmp_var_set.contains(curr))
1037 return true;
1038 }
1039 SASSERT(m == var2expr(v));
1040 for (expr* arg : *to_app(m)) {
1041 theory_var curr = expr2var(arg);
1042 if (!is_fixed(curr))
1043 m_tmp_var_set.insert(curr);
1044 }
1045 }
1046 else {
1047 if (m_tmp_var_set.contains(v))
1048 return true;
1049 SASSERT(!is_fixed(v));
1050 m_tmp_var_set.insert(v);
1051 }
1052 }
1053 }
1054 return false;
1055 }
1056
1057 /**
1058 \brief Return true if the row mixes real and integer variables.
1059 This kind of row cannot be converted back to an expression, since
1060 expressions in Z3 cannot have mixed sorts.
1061 */
1062 template<typename Ext>
is_mixed_real_integer(row const & r)1063 bool theory_arith<Ext>::is_mixed_real_integer(row const & r) const {
1064 bool found_int = false;
1065 bool found_real = false;
1066 typename vector<row_entry>::const_iterator it = r.begin_entries();
1067 typename vector<row_entry>::const_iterator end = r.end_entries();
1068 for (; it != end; ++it) {
1069 if (it->is_dead())
1070 continue;
1071 theory_var v = it->m_var;
1072 // TODO: possible improvement... ignore fixed variables.
1073 // If we implement this improvement, we are actually changing the contract of this function
1074 // and we will also have to fix the affected functions.
1075 if (is_int(v))
1076 found_int = true;
1077 if (is_real(v))
1078 found_real = true;
1079 if (found_int && found_real)
1080 return true;
1081 }
1082 return false;
1083 }
1084
1085 /**
1086 \brief Return true if the row contains only integer variables.
1087 */
1088 template<typename Ext>
is_integer(row const & r)1089 bool theory_arith<Ext>::is_integer(row const & r) const {
1090 typename vector<row_entry>::const_iterator it = r.begin_entries();
1091 typename vector<row_entry>::const_iterator end = r.end_entries();
1092 for (; it != end; ++it) {
1093 if (it->is_dead())
1094 continue;
1095 theory_var v = it->m_var;
1096 // TODO: possible improvement... ignore fixed variables.
1097 if (!is_int(v))
1098 return false;
1099 }
1100 return true;
1101 }
1102
1103 template<typename Ext>
display_coeff_exprs(std::ostream & out,buffer<coeff_expr> const & p)1104 void theory_arith<Ext>::display_coeff_exprs(std::ostream & out, buffer<coeff_expr> const & p) const {
1105 bool first = true;
1106 for (coeff_expr const& ce : p) {
1107 if (first)
1108 first = false;
1109 else
1110 out << "+\n";
1111 out << ce.first << " * " << mk_pp(ce.second, get_manager()) << "\n";
1112 }
1113 }
1114
1115 /**
1116 \brief Traverse p and store in vars the (non-fixed) variables
1117 that occur in more than one monomial. The number of
1118 occurrences is also stored.
1119 */
1120 template<typename Ext>
get_polynomial_info(buffer<coeff_expr> const & p,sbuffer<var_num_occs> & varinfo)1121 bool theory_arith<Ext>::get_polynomial_info(buffer<coeff_expr> const & p, sbuffer<var_num_occs> & varinfo) {
1122 varinfo.reset();
1123 m_var2num_occs.reset();
1124
1125 auto add_occ = [&](expr* v) {
1126 if (has_var(v) && !is_fixed(expr2var(v))) {
1127 TRACE("nl_info", tout << "adding occ: " << mk_bounded_pp(v, get_manager()) << "\n";);
1128 unsigned occs = 0;
1129 m_var2num_occs.find(v, occs);
1130 m_var2num_occs.insert(v, 1 + occs);
1131 }
1132 };
1133
1134 for (auto const& ce : p) {
1135 expr * m = ce.second;
1136 if (m_util.is_numeral(m)) {
1137 continue;
1138 }
1139 else if (false && m_util.is_add(m)) // introduced by #4532, disabled for #4765
1140 return false;
1141 else if (ctx.e_internalized(m) && !is_pure_monomial(m))
1142 add_occ(m);
1143 else {
1144 buffer<var_power_pair> vp;
1145 decompose_monomial(m, vp);
1146 for (auto const& p : vp) {
1147 add_occ(p.first);
1148 }
1149 }
1150 }
1151
1152 // Update the number of occurrences in the result vector.
1153 for (auto const& vn : m_var2num_occs) {
1154 if (vn.m_value > 1)
1155 varinfo.push_back(var_num_occs(vn.m_key, vn.m_value));
1156 }
1157 return true;
1158 }
1159
1160 /**
1161 \brief Convert p into an expression.
1162 */
1163 template<typename Ext>
p2expr(buffer<coeff_expr> & p)1164 expr_ref theory_arith<Ext>::p2expr(buffer<coeff_expr> & p) {
1165 SASSERT(!p.empty());
1166 TRACE("p2expr_bug", display_coeff_exprs(tout, p););
1167 ptr_buffer<expr> args;
1168 rational c2;
1169 for (coeff_expr const& ce : p) {
1170 rational const & c = ce.first;
1171 expr * var = ce.second;
1172 if (m_util.is_numeral(var, c2)) {
1173 expr* m = m_util.mk_numeral(c * c2, c.is_int() && m_util.is_int(var));
1174 m_nl_new_exprs.push_back(m);
1175 args.push_back(m);
1176 }
1177 else if (!c.is_one()) {
1178 expr * m = m_util.mk_mul(m_util.mk_numeral(c, c.is_int() && m_util.is_int(var)), var);
1179 m_nl_new_exprs.push_back(m);
1180 args.push_back(m);
1181 }
1182 else {
1183 args.push_back(var);
1184 }
1185 }
1186 SASSERT(!args.empty());
1187 expr_ref r(mk_nary_add(args.size(), args.data()), get_manager());
1188 m_nl_new_exprs.push_back(r);
1189 return r;
1190 }
1191
1192 /**
1193 \brief Return expression representing: var^power
1194 */
1195 template<typename Ext>
power(expr * var,unsigned power)1196 expr * theory_arith<Ext>::power(expr * var, unsigned power) {
1197 SASSERT(power > 0);
1198 expr * r = var;
1199 for (unsigned i = 1; i < power; i++)
1200 r = m_util.mk_mul(var, r);
1201 m_nl_new_exprs.push_back(r);
1202 return r;
1203 }
1204
1205 /**
1206 \brief Return true if var only occurs in two monovariate monomials,
1207 and return its power and coefficients and these monomials.
1208 The arguments i1 and i2 contain the position in p of the two monomials.
1209 */
1210 template<typename Ext>
in_monovariate_monomials(buffer<coeff_expr> & p,expr * var,unsigned & i1,rational & c1,unsigned & n1,unsigned & i2,rational & c2,unsigned & n2)1211 bool theory_arith<Ext>::in_monovariate_monomials(buffer<coeff_expr> & p, expr * var,
1212 unsigned & i1, rational & c1, unsigned & n1, unsigned & i2, rational & c2, unsigned & n2) {
1213 int idx = 0;
1214
1215 auto set_result = [&](unsigned i, unsigned power, coeff_expr const& ce) {
1216 if (idx == 0) {
1217 c1 = ce.first;
1218 n1 = power;
1219 idx = 1;
1220 i1 = i;
1221 }
1222 else if (idx == 1) {
1223 c2 = ce.first;
1224 n2 = power;
1225 idx = 2;
1226 i2 = i;
1227 }
1228 else {
1229 idx = 3;
1230 }
1231 };
1232
1233
1234 for (unsigned i = 0; i < p.size() && idx != 3; ++i) {
1235 auto const& ce = p[i];
1236 expr * m = ce.second;
1237 if (is_pure_monomial(m)) {
1238 buffer<var_power_pair> vp;
1239 decompose_monomial(m, vp);
1240 for (auto const& p : vp) {
1241 if (p.first == var) {
1242 if (vp.size() > 1)
1243 return false;
1244 set_result(i, p.second, ce); }
1245 }
1246 }
1247 else if (m == var) {
1248 set_result(i, 1, ce);
1249 }
1250 }
1251 if (idx != 2)
1252 return false;
1253 return true;
1254 }
1255
1256
1257 /**
1258 \brief Display a nested form expression
1259 */
1260 template<typename Ext>
display_nested_form(std::ostream & out,expr * p)1261 void theory_arith<Ext>::display_nested_form(std::ostream & out, expr * p) {
1262 if (has_var(p)) {
1263 out << "#" << p->get_id();
1264 }
1265 else if (m_util.is_add(p)) {
1266 SASSERT(!has_var(p));
1267 out << "(";
1268 for (unsigned i = 0; i < to_app(p)->get_num_args(); i++) {
1269 if (i > 0) out << " + ";
1270 display_nested_form(out, to_app(p)->get_arg(i));
1271 }
1272 out << ")";
1273 }
1274 else if (m_util.is_mul(p)) {
1275 buffer<var_power_pair> vp;
1276 rational c = decompose_monomial(p, vp);
1277 bool first = true;
1278 if (!c.is_one()) {
1279 out << c;
1280 first = false;
1281 }
1282 for (auto const& pair : vp) {
1283 if (first) first = false; else out << "*";
1284 expr * var = pair.first;
1285 unsigned power = pair.second;
1286 display_nested_form(out, var);
1287 if (power != 1)
1288 out << "^" << power;
1289 }
1290 }
1291 else {
1292 rational val;
1293 if (m_util.is_numeral(p, val))
1294 out << val;
1295 else
1296 out << "[unknown #" << p->get_id() << "]";
1297 }
1298 }
1299
1300 /**
1301 \brief Return the degree of var in m.
1302 */
1303 template<typename Ext>
get_degree_of(expr * m,expr * var)1304 unsigned theory_arith<Ext>::get_degree_of(expr * m, expr * var) {
1305 if (m == var)
1306 return 1;
1307 if (is_pure_monomial(m)) {
1308 buffer<var_power_pair> vp;
1309 decompose_monomial(m, vp);
1310 for (auto const& p : vp) {
1311 if (p.first == var)
1312 return p.second;
1313 }
1314 }
1315 return 0;
1316 }
1317
1318 /**
1319 \brief Return the minimal degree of var in the polynomial p.
1320 */
1321 template<typename Ext>
get_min_degree(buffer<coeff_expr> & p,expr * var)1322 unsigned theory_arith<Ext>::get_min_degree(buffer<coeff_expr> & p, expr * var) {
1323 SASSERT(!p.empty());
1324 SASSERT(var != 0);
1325 // get monomial where the degree of var is min.
1326 unsigned d = UINT_MAX; // min. degree of var
1327 for (coeff_expr const& ce : p) {
1328 expr * m = ce.second;
1329 d = std::min(d, get_degree_of(m, var));
1330 if (d == 0)
1331 return d;
1332 }
1333 SASSERT(d != UINT_MAX);
1334 return d;
1335 }
1336
1337 /**
1338 \brief Divide m by var^d.
1339 */
1340 template<typename Ext>
factor(expr * m,expr * var,unsigned d)1341 expr * theory_arith<Ext>::factor(expr * m, expr * var, unsigned d) {
1342 TRACE("factor", tout << "m: " << mk_pp(m, get_manager()) << "\nvar: " << mk_pp(var, get_manager()) << "\nd: " << d << "\n";);
1343 if (d == 0)
1344 return m;
1345 if (m == var) {
1346 SASSERT(d == 1);
1347 expr * result = m_util.mk_numeral(rational(1), m_util.is_int(var));
1348 m_nl_new_exprs.push_back(result);
1349 return result;
1350 }
1351
1352 ptr_buffer<expr> new_args;
1353 unsigned idx = 0;
1354 auto insert = [&](expr* arg) {
1355 if (idx < d && var == arg)
1356 ++idx;
1357 else
1358 new_args.push_back(arg);
1359 };
1360 while (m_util.is_mul(m) && idx < d) {
1361 unsigned sz = to_app(m)->get_num_args();
1362 for (unsigned i = 0; i + 1 < sz; ++i) {
1363 insert(to_app(m)->get_arg(i));
1364 }
1365 m = to_app(m)->get_arg(sz - 1);
1366 }
1367 insert(m);
1368 SASSERT(idx == d);
1369 TRACE("factor_bug", tout << "new_args:\n"; for(unsigned i = 0; i < new_args.size(); i++) tout << mk_pp(new_args[i], get_manager()) << "\n";);
1370 expr * result = mk_nary_mul(new_args.size(), new_args.data(), m_util.is_int(var));
1371 m_nl_new_exprs.push_back(result);
1372 TRACE("factor", tout << "result: " << mk_pp(result, get_manager()) << "\n";);
1373 return result;
1374 }
1375
1376 /**
1377 \brief Return the horner extension of p with respect to var.
1378 */
1379 template<typename Ext>
horner(unsigned depth,buffer<coeff_expr> & p,expr * var)1380 expr_ref theory_arith<Ext>::horner(unsigned depth, buffer<coeff_expr> & p, expr * var) {
1381 SASSERT(!p.empty());
1382 SASSERT(var != 0);
1383 unsigned d = get_min_degree(p, var);
1384 TRACE("horner_bug", tout << "poly:\n";
1385 for (unsigned i = 0; i < p.size(); i++) { if (i > 0) tout << " + "; tout << p[i].first << "*" << mk_pp(p[i].second, get_manager()); } tout << "\n";
1386 tout << "var: " << mk_pp(var, get_manager()) << "\n";
1387 tout << "min_degree: " << d << "\n";);
1388 buffer<coeff_expr> e; // monomials/x^d where var occurs with degree d
1389 buffer<coeff_expr> r; // rest
1390 for (auto const& kv : p) {
1391 expr * m = kv.second;
1392 expr * f = factor(m, var, d);
1393 if (get_degree_of(m, var) == d) {
1394 e.push_back(coeff_expr(kv.first, f));
1395 }
1396 else {
1397 SASSERT(get_degree_of(m, var) > d);
1398 r.push_back(coeff_expr(kv.first, f));
1399 }
1400 }
1401 expr_ref s = cross_nested(depth + 1, e, nullptr);
1402 if (!r.empty()) {
1403 expr_ref q = horner(depth + 1, r, var);
1404 // TODO: improve here
1405 s = m_util.mk_add(q, s);
1406 }
1407
1408 expr_ref result = s;
1409 if (d != 0) {
1410 expr * xd = power(var, d);
1411 result = m_util.mk_mul(xd, s);
1412 }
1413 m_nl_new_exprs.push_back(result);
1414 return result;
1415 }
1416
1417
1418 /**
1419 \brief Convert the polynomial p into an equivalent cross nested
1420 expression. The idea is to obtain an expression e where
1421 evaluate_as_interval(e) is more precise than
1422 evaluate_as_interval(p).
1423
1424 If var != 0, then it is used for performing the horner extension
1425 */
1426 template<typename Ext>
cross_nested(unsigned depth,buffer<coeff_expr> & p,expr * var)1427 expr_ref theory_arith<Ext>::cross_nested(unsigned depth, buffer<coeff_expr> & p, expr * var) {
1428 TRACE("non_linear", tout << "p.size: " << p.size() << "\n";);
1429 if (var == nullptr) {
1430 sbuffer<var_num_occs> varinfo;
1431 if (!get_polynomial_info(p, varinfo))
1432 return p2expr(p);
1433 if (varinfo.empty())
1434 return p2expr(p);
1435 unsigned max = 0;
1436 for (auto const& kv : varinfo) {
1437 if (kv.second >= max) {
1438 max = kv.second;
1439 var = kv.first;
1440 }
1441 }
1442 }
1443 if (depth > 20)
1444 return p2expr(p);
1445 SASSERT(var != nullptr);
1446 unsigned i1 = UINT_MAX;
1447 unsigned i2 = UINT_MAX;
1448 rational a, b;
1449 unsigned n = UINT_MAX;
1450 unsigned nm = UINT_MAX;
1451 if (in_monovariate_monomials(p, var, i1, a, n, i2, b, nm)) {
1452 CTRACE("in_monovariate_monomials", n == nm,
1453 for (unsigned i = 0; i < p.size(); i++) {
1454 if (i > 0) tout << " + "; tout << p[i].first << "*" << mk_pp(p[i].second, get_manager());
1455 }
1456 tout << "\n";
1457 tout << "var: " << mk_pp(var, get_manager()) << "\n";
1458 tout << "i1: " << i1 << "\n";
1459 tout << "a: " << a << "\n";
1460 tout << "n: " << n << "\n";
1461 tout << "i2: " << i2 << "\n";
1462 tout << "b: " << b << "\n";
1463 tout << "nm: " << nm << "\n";);
1464 if (n == nm) return horner(depth, p, var);
1465 SASSERT(n != nm);
1466 expr * new_expr = nullptr;
1467 if (nm < n) {
1468 std::swap(n, nm);
1469 std::swap(a, b);
1470 }
1471 SASSERT(nm > n);
1472 unsigned m = nm - n;
1473 if (n % 2 == m % 2 && n >= m) {
1474 // b*x^{n-m}*[(x^{m} + a/(2b))^2 - (a/2b)^2]
1475 // b*[(x^{m} + a/(2b))^2 - (a/2b)^2] for n == m
1476 rational a2b = a;
1477 expr_ref xm(power(var, m), get_manager());
1478 a2b /= (rational(2) * b);
1479 // we cannot create a numeral that has sort int, but it is a rational.
1480 if (!m_util.is_int(var) || a2b.is_int()) {
1481 rational ma2b2 = a2b * a2b;
1482 ma2b2.neg();
1483 expr * xm_a2b = m_util.mk_add(m_util.mk_numeral(a2b, m_util.is_int(var)), xm);
1484 expr * xm_a2b2 = m_util.mk_mul(xm_a2b, xm_a2b);
1485 expr * rhs = m_util.mk_add(xm_a2b2, m_util.mk_numeral(ma2b2, m_util.is_int(var)));
1486 expr * rhs2 = nullptr;
1487 if (n > m)
1488 rhs2 = m_util.mk_mul(power(var, n - m), rhs);
1489 else
1490 rhs2 = rhs;
1491 new_expr = b.is_one() ? rhs2 : m_util.mk_mul(m_util.mk_numeral(b, m_util.is_int(var)), rhs2);
1492 m_nl_new_exprs.push_back(new_expr);
1493 TRACE("non_linear", tout << "new_expr:\n"; display_nested_form(tout, new_expr); tout << "\n";);
1494 buffer<coeff_expr> rest;
1495 unsigned sz = p.size();
1496 for (unsigned i = 0; i < sz; i++) {
1497 if (i != i1 && i != i2)
1498 rest.push_back(p[i]);
1499 }
1500 if (rest.empty())
1501 return expr_ref(new_expr, get_manager());
1502 TRACE("non_linear", tout << "rest size: " << rest.size() << ", i1: " << i1 << ", i2: " << i2 << "\n";);
1503 expr_ref h = cross_nested(depth + 1, rest, nullptr);
1504 expr * r = m_util.mk_add(new_expr, h);
1505 m_nl_new_exprs.push_back(r);
1506 return expr_ref(r, get_manager());
1507 }
1508 }
1509 }
1510 return horner(depth, p, var);
1511 }
1512
1513 /**
1514 \brief Check whether the given polynomial is consistent with respect to the known
1515 bounds. The polynomial is converted into an equivalent cross nested form.
1516 */
1517 template<typename Ext>
is_cross_nested_consistent(buffer<coeff_expr> & p)1518 bool theory_arith<Ext>::is_cross_nested_consistent(buffer<coeff_expr> & p) {
1519 sbuffer<var_num_occs> varinfo;
1520 if (!get_polynomial_info(p, varinfo))
1521 return true;
1522 if (varinfo.empty())
1523 return true;
1524 std::stable_sort(varinfo.begin(), varinfo.end(), var_num_occs_lt());
1525 TRACE("cross_nested", tout << "var num occs:\n";
1526 for (auto const& kv : varinfo) {
1527 tout << mk_bounded_pp(kv.first, get_manager()) << " -> " << kv.second << "\n";
1528 });
1529 for (auto const& kv : varinfo) {
1530 m_nl_new_exprs.reset();
1531 expr * var = kv.first;
1532 expr_ref cn = cross_nested(0, p, var);
1533 // Remark: cn may not be well-sorted because, since a row may contain mixed integer/real monomials.
1534 // This is not really a problem, since evaluate_as_interval will work even if cn is not well-sorted.
1535 if (!cn)
1536 continue;
1537 TRACE("cross_nested", tout << "nested form for var:\n" << mk_ismt2_pp(var, get_manager()) << "\n";
1538 display_nested_form(tout, cn); tout << "\n";
1539 tout << "c:\n" << mk_ismt2_pp(cn, get_manager()) << "\n";);
1540 interval i = evaluate_as_interval(cn);
1541 TRACE("cross_nested", tout << "interval: " << i << "\n";);
1542 v_dependency * d = nullptr;
1543 if (!i.minus_infinity() && (i.get_lower_value().is_pos() || (i.get_lower_value().is_zero() && i.is_lower_open())))
1544 d = i.get_lower_dependencies();
1545 else if (!i.plus_infinity() && (i.get_upper_value().is_neg() || (i.get_upper_value().is_zero() && i.is_upper_open())))
1546 d = i.get_upper_dependencies();
1547 if (d) {
1548 TRACE("cross_nested", tout << "nested form conflict: " << i << "\n";);
1549 set_conflict(d);
1550 return false;
1551 }
1552 }
1553 return true;
1554 }
1555
1556 /**
1557 \brief Check whether the polynomial represented by the current row is
1558 consistent with respect to the known bound when converted into a
1559 equivalent cross nested form.
1560 */
1561 template<typename Ext>
is_cross_nested_consistent(row const & r)1562 bool theory_arith<Ext>::is_cross_nested_consistent(row const & r) {
1563 if (!is_problematic_non_linear_row(r))
1564 return true;
1565
1566 TRACE("cross_nested", tout << "is_cross_nested_consistent:\n"; display_row(tout, r, false);
1567 display(tout);
1568 );
1569
1570 /*
1571 The method is_cross_nested converts rows back to expressions.
1572 The conversion back to expressions may create sort incorrect expressions.
1573 This is in some sense ok, since these expressions are temporary, but
1574 the sort incorrect expressions may generate assertion violations.
1575
1576 Sort incorrect expressions may be created in the following cases:
1577
1578 1) mixed real int rows.
1579
1580 2) int rows that contain non integer coefficients.
1581
1582 3) int rows that when converted to cross nested form use non integer coefficients.
1583
1584 There are several ways to deal with this problem:
1585
1586 a) Ignore the assertion violations. Disadvantage: it will prevent us from running Z3 in debug mode on some benchmarks.
1587
1588 b) Remove the assertions. Disadvantage: these assertions helped us to find many important bugs in Z3
1589
1590 c) Disable the assertions temporally. This sounds like a big HACK.
1591
1592 d) Use a different data-structure to represent polynomials in cross-nested form. Disadvantage: code duplication, the data-structure
1593 is essentially identical to the ASTs we are using right now.
1594
1595 e) Disable the test when we cannot create a well-sorted expression.
1596 I'm temporally using this solution.
1597 I implemented the following logic:
1598 1) (mixed real int) Disable the test. Most benchmarks do not contain mixed int real variables.
1599 2) (int coeffs) I multiply the row by a constant to force it to have only integer coefficients.
1600 3) (new non-int coeffs) This only happens in an optional step in the conversion. Now, for int rows, I only apply this optional step only if non-int coeffs are not created.
1601 */
1602
1603 if (!get_manager().int_real_coercions() && is_mixed_real_integer(r))
1604 return true; // giving up... see comment above
1605
1606 TRACE("cross_nested", tout << "checking problematic row...\n";);
1607
1608 rational c = rational::one();
1609 if (is_integer(r))
1610 c = r.get_denominators_lcm().to_rational();
1611
1612 TRACE("non_linear", tout << "check problematic row:\n"; display_row(tout, r); display_row(tout, r, false););
1613 buffer<coeff_expr> p;
1614 for (auto & col : r) {
1615 if (!col.is_dead()) {
1616 p.push_back(coeff_expr(col.m_coeff.to_rational() * c, var2expr(col.m_var)));
1617 }
1618 }
1619 SASSERT(!p.empty());
1620 CTRACE("cross_nested_bug", !c.is_one(), tout << "c: " << c << "\n"; display_row(tout, r); tout << "---> p (coeffs, exprs):\n"; display_coeff_exprs(tout, p););
1621 return is_cross_nested_consistent(p);
1622 }
1623
1624 /**
1625 \brief Check whether an inconsistency can be found using cross nested
1626 form in the non linear cluster.
1627 */
1628 template<typename Ext>
is_cross_nested_consistent(svector<theory_var> const & nl_cluster)1629 bool theory_arith<Ext>::is_cross_nested_consistent(svector<theory_var> const & nl_cluster) {
1630 for (theory_var v : nl_cluster) {
1631 if (!is_base(v))
1632 continue;
1633 m_stats.m_nl_cross_nested++;
1634 row const & r = m_rows[get_var_row(v)];
1635 if (!is_cross_nested_consistent(r))
1636 return false;
1637 }
1638 return true;
1639 }
1640
1641 #define FIXED 0
1642 #define QUOTED_FIXED 1
1643 #define BOUNDED 2
1644 #define QUOTED_BOUNDED 3
1645 #define NOT_FREE 4
1646 #define QUOTED_NOT_FREE 5
1647 #define FREE 6
1648 #define QUOTED_FREE 7
1649 #define MAX_DEFAULT_WEIGHT 7
1650
1651 /**
1652 \brief Initialize variable order for grobner basis computation.
1653 Make:
1654 "quoted free vars" > "free vars" > "quoted variables with lower or upper bounds" >
1655 "variables with lower or upper bounds" > "quoted bounded variables" >
1656 "bounded variables" > "quoted fixed variables" > "fixed variables"
1657 */
1658 template<typename Ext>
init_grobner_var_order(svector<theory_var> const & nl_cluster,grobner & gb)1659 void theory_arith<Ext>::init_grobner_var_order(svector<theory_var> const & nl_cluster, grobner & gb) {
1660 // Initialize variable order
1661 for (theory_var v : nl_cluster) {
1662 expr * var = var2expr(v);
1663
1664 if (is_fixed(v)) {
1665 gb.set_weight(var, is_pure_monomial(var) ? QUOTED_FIXED : FIXED);
1666 }
1667 else if (is_bounded(v)) {
1668 gb.set_weight(var, is_pure_monomial(var) ? QUOTED_BOUNDED : BOUNDED);
1669 }
1670 else if (lower(v) || upper(v)) {
1671 gb.set_weight(var, is_pure_monomial(var) ? QUOTED_NOT_FREE : NOT_FREE);
1672 }
1673 else {
1674 SASSERT(is_free(v));
1675 gb.set_weight(var, is_pure_monomial(var) ? QUOTED_FREE : FREE);
1676 }
1677 }
1678 }
1679
1680
1681 /**
1682 \brief Create a new monomial using the given coeff and m. Fixed
1683 variables in m are substituted by their values. The arg dep is
1684 updated to store these dependencies. The set already_found is
1685 updated with the fixed variables in m. A variable is only
1686 added to dep if it is not already in already_found.
1687
1688 Return null if the monomial was simplified to 0.
1689 */
1690 template<typename Ext>
mk_gb_monomial(rational const & _coeff,expr * m,grobner & gb,v_dependency * & dep,var_set & already_found)1691 grobner::monomial * theory_arith<Ext>::mk_gb_monomial(rational const & _coeff, expr * m, grobner & gb, v_dependency * & dep, var_set & already_found) {
1692 ptr_buffer<expr> vars;
1693 rational coeff = _coeff;
1694 rational r;
1695 auto proc_var = [&](expr* v) {
1696 if (m_util.is_numeral(v, r)) {
1697 coeff *= r;
1698 return;
1699 }
1700 theory_var _var = expr2var(v);
1701 if (is_fixed(_var)) {
1702 if (!already_found.contains(_var)) {
1703 already_found.insert(_var);
1704 dep = m_dep_manager.mk_join(dep, m_dep_manager.mk_join(m_dep_manager.mk_leaf(lower(_var)), m_dep_manager.mk_leaf(upper(_var))));
1705 }
1706 coeff *= lower_bound(_var).get_rational().to_rational();
1707 }
1708 else {
1709 vars.push_back(v);
1710 }
1711 };
1712
1713 while (m_util.is_mul(m)) {
1714 unsigned sz = to_app(m)->get_num_args();
1715 for (unsigned i = 0; i + 1 < sz; ++i) {
1716 proc_var(to_app(m)->get_arg(i));
1717 }
1718 m = to_app(m)->get_arg(sz-1);
1719 }
1720 proc_var(m);
1721
1722 if (!coeff.is_zero())
1723 return gb.mk_monomial(coeff, vars.size(), vars.data());
1724 else
1725 return nullptr;
1726 }
1727
1728 /**
1729 \brief Send the given row to the grobner basis object.
1730 All fixed variables are substituted before sending the row to gb.
1731 */
1732 template<typename Ext>
add_row_to_gb(row const & r,grobner & gb)1733 void theory_arith<Ext>::add_row_to_gb(row const & r, grobner & gb) {
1734 TRACE("non_linear", tout << "adding row to gb\n"; display_row(tout, r););
1735 ptr_buffer<grobner::monomial> monomials;
1736 v_dependency * dep = nullptr;
1737 m_tmp_var_set.reset();
1738 typename vector<row_entry>::const_iterator it = r.begin_entries();
1739 typename vector<row_entry>::const_iterator end = r.end_entries();
1740 for (; it != end; ++it) {
1741 if (!it->is_dead()) {
1742 rational coeff = it->m_coeff.to_rational();
1743 expr * m = var2expr(it->m_var);
1744 TRACE("non_linear", tout << "monomial: " << mk_pp(m, get_manager()) << "\n";);
1745 grobner::monomial * new_m = mk_gb_monomial(coeff, m, gb, dep, m_tmp_var_set);
1746 TRACE("non_linear", tout << "new monomial:\n"; if (new_m) gb.display_monomial(tout, *new_m); else tout << "null"; tout << "\n";);
1747 if (new_m)
1748 monomials.push_back(new_m);
1749 }
1750 }
1751 gb.assert_eq_0(monomials.size(), monomials.data(), dep);
1752 }
1753
1754 /**
1755 \brief v must be a pure monomial. That is, v = (quote (* x_1 ... x_n))
1756 Add the monomial (quote (* x_1 ... x_n)) = x_1 * ... * x_n.
1757 Fixed variables are substituted.
1758 */
1759 template<typename Ext>
add_monomial_def_to_gb(theory_var v,grobner & gb)1760 void theory_arith<Ext>::add_monomial_def_to_gb(theory_var v, grobner & gb) {
1761 ptr_buffer<grobner::monomial> monomials;
1762 v_dependency * dep = nullptr;
1763 m_tmp_var_set.reset();
1764 expr * m = var2expr(v);
1765 SASSERT(is_pure_monomial(m));
1766 grobner::monomial * new_m = mk_gb_monomial(rational(1), m, gb, dep, m_tmp_var_set);
1767 if (new_m)
1768 monomials.push_back(new_m);
1769 rational coeff(-1);
1770 if (is_fixed(v)) {
1771 dep = m_dep_manager.mk_join(dep, m_dep_manager.mk_join(m_dep_manager.mk_leaf(lower(v)), m_dep_manager.mk_leaf(upper(v))));
1772 coeff *= lower_bound(v).get_rational().to_rational();
1773 if (!coeff.is_zero())
1774 monomials.push_back(gb.mk_monomial(coeff, 0, nullptr));
1775 }
1776 else {
1777 monomials.push_back(gb.mk_monomial(coeff, 1, &m));
1778 }
1779 gb.assert_eq_0(monomials.size(), monomials.data(), dep);
1780 }
1781
1782 /**
1783 Initialize grobner basis data structure using the non linear cluster.
1784 The GB is initialized using rows and non linear monomials.
1785 */
1786 template<typename Ext>
init_grobner(svector<theory_var> const & nl_cluster,grobner & gb)1787 void theory_arith<Ext>::init_grobner(svector<theory_var> const & nl_cluster, grobner & gb) {
1788 init_grobner_var_order(nl_cluster, gb);
1789 for (theory_var v : nl_cluster) {
1790 if (is_base(v)) {
1791 row const & r = m_rows[get_var_row(v)];
1792 add_row_to_gb(r, gb);
1793 }
1794 if (is_pure_monomial(v) && !m_data[v].m_nl_propagated && is_fixed(v)) {
1795 add_monomial_def_to_gb(v, gb);
1796 }
1797 }
1798 }
1799
1800 /**
1801 \brief Return the interval for the given monomial
1802 */
1803 template<typename Ext>
mk_interval_for(grobner::monomial const * m)1804 interval theory_arith<Ext>::mk_interval_for(grobner::monomial const * m) {
1805 interval r(m_dep_manager, rational(m->get_coeff()));
1806 expr * var = nullptr;
1807 unsigned power = 0;
1808 unsigned num_vars = m->get_degree();
1809 for (unsigned i = 0; i < num_vars; i++) {
1810 expr * curr = m->get_var(i);
1811 if (var == nullptr) {
1812 var = curr;
1813 power = 1;
1814 }
1815 else if (curr == var) {
1816 power++;
1817 }
1818 else {
1819 mul_bound_of(var, power, r);
1820 var = curr;
1821 power = 1;
1822 }
1823 }
1824 if (var != nullptr)
1825 mul_bound_of(var, power, r);
1826 return r;
1827 }
1828
1829 /**
1830 \brief Set a conflict using a dependency object.
1831 */
1832 template<typename Ext>
set_conflict(v_dependency * d)1833 void theory_arith<Ext>::set_conflict(v_dependency * d) {
1834 antecedents ante(*this);
1835 derived_bound b(null_theory_var, inf_numeral(0), B_LOWER);
1836 dependency2new_bound(d, b);
1837 set_conflict(b, ante, "arith_nl");
1838 TRACE("non_linear", for (literal lit : b.m_lits) ctx.display_literal_verbose(tout, lit) << "\n"; tout << "\n";);
1839 }
1840
1841 /**
1842 \brief Return true if I.get_lower() <= - M_1 - ... - M_n <= I.get_upper() is inconsistent.
1843 Where M_i is monomials[i] and n = num_monomials.
1844 A conflict will also be set using the bounds of the variables occurring in the monomials M_i's.
1845 */
1846 template<typename Ext>
is_inconsistent(interval const & I,unsigned num_monomials,grobner::monomial * const * monomials,v_dependency * dep)1847 bool theory_arith<Ext>::is_inconsistent(interval const & I, unsigned num_monomials, grobner::monomial * const * monomials, v_dependency * dep) {
1848 interval r(I);
1849 for (unsigned i = 0; i < num_monomials; i++) {
1850 grobner::monomial const * m = monomials[i];
1851 r += mk_interval_for(m);
1852 if (r.minus_infinity() && r.plus_infinity())
1853 return false;
1854 }
1855 TRACE("non_linear_bug", tout << "is_inconsistent, r: " << r << "\n";);
1856 v_dependency * interval_deps = nullptr;
1857 bool conflict = false;
1858 if (!r.minus_infinity() && (r.get_lower_value().is_pos() || (r.get_lower_value().is_zero() && r.is_lower_open()))) {
1859 interval_deps = r.get_lower_dependencies();
1860 conflict = true;
1861 TRACE("non_linear_bug", tout << "is inconsistent, interval_deps: " << interval_deps << "\n";);
1862 }
1863 else if (!r.plus_infinity() && (r.get_upper_value().is_neg() || (r.get_upper_value().is_zero() && r.is_upper_open()))) {
1864 interval_deps = r.get_upper_dependencies();
1865 conflict = true;
1866 TRACE("non_linear_bug", tout << "is inconsistent, interval_deps: " << interval_deps << "\n";);
1867 }
1868 // interval_deps cannot be used to check if a conflict was detected, since interval_deps may be 0 even when r does not contain 0
1869 if (conflict) {
1870 TRACE("non_linear", tout << "conflicting interval for = 0 equation: " << r << "\n";);
1871 set_conflict(m_dep_manager.mk_join(interval_deps, dep));
1872 return true;
1873 }
1874 return false;
1875 }
1876
1877 /**
1878 \brief Return true if the equation is inconsistent,
1879 and sign a conflict.
1880 */
1881 template<typename Ext>
is_inconsistent(grobner::equation const * eq,grobner & gb)1882 bool theory_arith<Ext>::is_inconsistent(grobner::equation const * eq, grobner & gb) {
1883 interval zero(m_dep_manager, rational(0));
1884 if (is_inconsistent(zero, eq->get_num_monomials(), eq->get_monomials(), eq->get_dependency())) {
1885 TRACE("non_linear", tout << "found conflict\n"; gb.display_equation(tout, *eq););
1886 return true;
1887 }
1888 return false;
1889 }
1890
1891 /**
1892 \brief Return true if the given monomial c*M is squared.
1893 The square root of the c is stored in r.
1894 */
is_perfect_square(grobner::monomial const * m,rational & r)1895 bool is_perfect_square(grobner::monomial const * m, rational & r) {
1896 unsigned num_vars = m->get_degree();
1897 if (num_vars % 2 == 1)
1898 return false;
1899 if (!m->get_coeff().is_perfect_square(r))
1900 return false;
1901 expr * var = nullptr;
1902 unsigned power = 0;
1903 for (unsigned i = 0; i < num_vars; i++) {
1904 expr * curr = m->get_var(i);
1905 if (var == nullptr) {
1906 var = curr;
1907 power = 1;
1908 }
1909 else if (var == curr) {
1910 power++;
1911 }
1912 else {
1913 if (power % 2 == 1)
1914 return false;
1915 var = curr;
1916 power = 1;
1917 }
1918 }
1919 return power % 2 == 0;
1920 }
1921
1922 /**
1923 \brief Return m1m2 is of the form (-2ab)*M1*M2
1924 assuming that
1925 m1_sq = a^2*M1*M1
1926 m2_sq = b^2*M2*M2
1927 */
is_perfect_square(grobner::monomial const * m1_sq,rational const & a,grobner::monomial const * m2_sq,rational const & b,grobner::monomial const * m1m2)1928 bool is_perfect_square(grobner::monomial const * m1_sq, rational const & a, grobner::monomial const * m2_sq, rational const & b, grobner::monomial const * m1m2) {
1929 DEBUG_CODE({
1930 rational a1;
1931 rational b1;
1932 SASSERT(is_perfect_square(m1_sq, a1) && a == a1 && is_perfect_square(m2_sq, b1) && b == b1);
1933 });
1934 if (m1m2->get_coeff().is_nonneg())
1935 return false;
1936 rational c(-2);
1937 c *= a;
1938 c *= b;
1939 if (m1m2->get_coeff() != c)
1940 return false;
1941 unsigned num1 = m1_sq->get_degree();
1942 unsigned num2 = m2_sq->get_degree();
1943 unsigned num12 = m1m2->get_degree();
1944 if (num1 + num2 != num12 * 2)
1945 return false;
1946 unsigned i1, i2, i12;
1947 i1 = i2 = i12 = 0;
1948 while (true) {
1949 expr * v1 = nullptr;
1950 expr * v2 = nullptr;
1951 expr * v12 = nullptr;
1952 if (i1 < num1)
1953 v1 = m1_sq->get_var(i1);
1954 if (i2 < num2)
1955 v2 = m2_sq->get_var(i2);
1956 if (i12 < num12)
1957 v12 = m1m2->get_var(i12);
1958 if (v1 == nullptr && v2 == nullptr && v12 == nullptr)
1959 return true;
1960 if (v12 == nullptr)
1961 return false;
1962 if (v1 == v12) {
1963 SASSERT(m1_sq->get_var(i1+1) == v1);
1964 i1 += 2;
1965 i12 ++;
1966 }
1967 else if (v2 == v12) {
1968 SASSERT(m2_sq->get_var(i2+1) == v2);
1969 i2 += 2;
1970 i12 ++;
1971 }
1972 else {
1973 return false;
1974 }
1975 }
1976 }
1977
1978 /**
1979 \brief Return true if the equation is inconsistent. In this
1980 version, perfect squares are eliminated, and replaced with the
1981 interval [0, oo), if the interval associated with them is less
1982 precise than [0, oo).
1983
1984 \remark I track only simple perfect squares of the form (M1 - M2)^2,
1985 where M1 and M2 are arbitrary monomials.
1986 */
1987 template<typename Ext>
is_inconsistent2(grobner::equation const * eq,grobner & gb)1988 bool theory_arith<Ext>::is_inconsistent2(grobner::equation const * eq, grobner & gb) {
1989 // TODO: a possible improvement: create a quotation for (M1 - M2)^2
1990 // instead of trying to find it in a specific equation.
1991 // This approach is more precise, but more expensive
1992 // since a new row must be created.
1993 buffer<interval> intervals;
1994 unsigned num = eq->get_num_monomials();
1995 for (unsigned i = 0; i < num; i++) {
1996 grobner::monomial const * m = eq->get_monomial(i);
1997 intervals.push_back(mk_interval_for(m));
1998 }
1999 sbuffer<bool> deleted;
2000 deleted.resize(num, false);
2001 ptr_buffer<grobner::monomial> monomials;
2002 // try to eliminate monomials that form perfect squares of the form (M1 - M2)^2
2003 for (unsigned i = 0; i < num; i++) {
2004 grobner::monomial const * m1 = eq->get_monomial(i);
2005 rational a;
2006 if (deleted[i])
2007 continue;
2008 if (!is_perfect_square(m1, a)) {
2009 monomials.push_back(const_cast<grobner::monomial*>(m1));
2010 continue;
2011 }
2012 TRACE("non_linear", tout << "found perfect square monomial m1: "; gb.display_monomial(tout, *m1); tout << "\n";);
2013 // try to find another perfect square
2014 unsigned j = i + 1;
2015 for (; j < num; j++) {
2016 if (deleted[j])
2017 continue;
2018 grobner::monomial const * m2 = eq->get_monomial(j);
2019 rational b;
2020 if (!is_perfect_square(m2, b))
2021 continue;
2022 TRACE("non_linear", tout << "found perfect square monomial m2: "; gb.display_monomial(tout, *m2); tout << "\n";);
2023 // try to find -2*root(m1)*root(m2)
2024 // This monomial must be smaller than m1, since m2 is smaller than m1.
2025 unsigned k = i + 1;
2026 for (; k < num; k++) {
2027 if (deleted[k])
2028 continue;
2029 grobner::monomial const * m1m2 = eq->get_monomial(k);
2030 if (!is_perfect_square(m1, a, m2, b, m1m2))
2031 continue;
2032 // m1, m2, and m1m2 form a perfect square.
2033 // check if [0, oo) provides a better lowerbound than adding the intervals of m1, m2 and m1m2;
2034 TRACE("non_linear", tout << "found perfect square (M1-M2)^2:\n";
2035 gb.display_monomial(tout, *m1); tout << "\n";
2036 gb.display_monomial(tout, *m2); tout << "\n";
2037 gb.display_monomial(tout, *m1m2); tout << "\n";);
2038 interval I = intervals[i];
2039 I += intervals[j];
2040 I += intervals[k];
2041 if (I.minus_infinity() || I.get_lower_value().is_neg()) {
2042 TRACE("non_linear", tout << "the lower bound improved when perfect square is eliminated.\n";);
2043 // Found improvement...
2044 // mark these monomials as deleted
2045 deleted[i] = true;
2046 deleted[j] = true;
2047 deleted[k] = true;
2048 break;
2049 }
2050 }
2051 if (k < num)
2052 break; // found perfect square
2053 }
2054 if (j == num) {
2055 // didn't find perfect square of the form (M1-M2)^2
2056 monomials.push_back(const_cast<grobner::monomial*>(m1));
2057 }
2058 }
2059 if (monomials.size() == num)
2060 return false; // didn't find any perfect square.
2061 interval ge_zero(m_dep_manager, rational(0), false, true, nullptr);
2062 if (is_inconsistent(ge_zero, monomials.size(), monomials.data(), eq->get_dependency())) {
2063 TRACE("non_linear", tout << "found conflict\n"; gb.display_equation(tout, *eq););
2064 return true;
2065 }
2066 return false;
2067 }
2068
2069 template<typename Ext>
monomial2expr(grobner::monomial const * m,bool is_int)2070 expr * theory_arith<Ext>::monomial2expr(grobner::monomial const * m, bool is_int) {
2071 unsigned num_vars = m->get_degree();
2072 ptr_buffer<expr> args;
2073 if (!m->get_coeff().is_one())
2074 args.push_back(m_util.mk_numeral(m->get_coeff(), is_int));
2075 for (unsigned j = 0; j < num_vars; j++)
2076 args.push_back(m->get_var(j));
2077 return mk_nary_mul(args.size(), args.data(), is_int);
2078 }
2079
2080 /**
2081 \brief Assert the new equation in the simplex tableau.
2082 */
2083 template<typename Ext>
internalize_gb_eq(grobner::equation const * eq)2084 bool theory_arith<Ext>::internalize_gb_eq(grobner::equation const * eq) {
2085 bool is_int = false;
2086 unsigned num_monomials = eq->get_num_monomials();
2087 for (unsigned i = 0; i < num_monomials; i++) {
2088 grobner::monomial const * m = eq->get_monomial(i);
2089 unsigned degree = m->get_degree();
2090 if (degree > m_params.m_nl_arith_max_degree)
2091 return false;
2092 if (degree > 0)
2093 is_int = m_util.is_int(m->get_var(0));
2094 }
2095 rational k;
2096 ptr_buffer<expr> args;
2097 for (unsigned i = 0; i < num_monomials; i++) {
2098 grobner::monomial const * m = eq->get_monomial(i);
2099 if (m->get_degree() == 0)
2100 k -= m->get_coeff();
2101 else
2102 args.push_back(monomial2expr(eq->get_monomial(i), is_int));
2103 }
2104 th_rewriter& s = ctx.get_rewriter();
2105 expr_ref pol(get_manager());
2106 SASSERT(!args.empty());
2107 pol = mk_nary_add(args.size(), args.data());
2108 expr_ref s_pol(get_manager());
2109 proof_ref pr(get_manager());
2110 TRACE("gb_bug", tout << mk_ll_pp(pol, get_manager()) << "\n";);
2111 s(pol, s_pol, pr);
2112 if (!has_var(s_pol)) {
2113 TRACE("spol_bug", tout << "internalizing...\n" << mk_ll_pp(s_pol, get_manager()) << "\n";);
2114 ctx.internalize(s_pol, false);
2115 ctx.mark_as_relevant(s_pol.get());
2116 }
2117 SASSERT(has_var(s_pol.get()));
2118 // s_pol = k
2119 theory_var v = expr2var(s_pol);
2120 // v = k
2121 CTRACE("spol_bug", v == null_theory_var, tout << mk_ll_pp(s_pol, get_manager()) << "\n"; display(tout););
2122 SASSERT(v != null_theory_var);
2123 // assert bounds for s_pol
2124 mk_derived_nl_bound(v, inf_numeral(k), B_LOWER, eq->get_dependency());
2125 mk_derived_nl_bound(v, inf_numeral(k), B_UPPER, eq->get_dependency());
2126 TRACE("non_linear", tout << "inserted new equation into the tableau\n"; display_var(tout, v););
2127 return true;
2128 }
2129
2130 template<typename Ext>
update_statistics(grobner & gb)2131 void theory_arith<Ext>::update_statistics(grobner & gb) {
2132 m_stats.m_gb_simplify += gb.m_stats.m_simplify;
2133 m_stats.m_gb_superpose += gb.m_stats.m_superpose;
2134 m_stats.m_gb_num_processed += gb.m_stats.m_num_processed;
2135 m_stats.m_gb_compute_basis++;
2136 }
2137
2138 template<typename Ext>
set_gb_exhausted()2139 void theory_arith<Ext>::set_gb_exhausted() {
2140 IF_VERBOSE(3, verbose_stream() << "Grobner basis computation interrupted. Increase threshold using NL_ARITH_GB_THRESHOLD=<limit>\n";);
2141 ctx.push_trail(value_trail<bool>(m_nl_gb_exhausted));
2142 m_nl_gb_exhausted = true;
2143 }
2144
2145 // Scan the grobner basis eqs, and look for inconsistencies.
2146 template<typename Ext>
get_gb_eqs_and_look_for_conflict(ptr_vector<grobner::equation> & eqs,grobner & gb)2147 bool theory_arith<Ext>::get_gb_eqs_and_look_for_conflict(ptr_vector<grobner::equation>& eqs, grobner& gb) {
2148 TRACE("grobner", );
2149
2150 eqs.reset();
2151 gb.get_equations(eqs);
2152 TRACE("grobner_bug", tout << "after gb\n";);
2153 for (grobner::equation* eq : eqs) {
2154 TRACE("grobner_bug", gb.display_equation(tout, *eq););
2155 if (is_inconsistent(eq, gb) || is_inconsistent2(eq, gb)) {
2156 TRACE("grobner", tout << "inconsistent: "; gb.display_equation(tout, *eq););
2157 return true;
2158 }
2159 }
2160 TRACE("grobner", tout << "not found\n"; );
2161 return false;
2162 }
2163
2164 // Try to change the variable order... in such a way the leading term is modified.
2165 // I only consider linear equations... (HACK)
2166 // Moreover, I do not change the weight of a variable more than once in this loop.
2167 template<typename Ext>
try_to_modify_eqs(ptr_vector<grobner::equation> & eqs,grobner & gb,unsigned & next_weight)2168 bool theory_arith<Ext>::try_to_modify_eqs(ptr_vector<grobner::equation>& eqs, grobner& gb, unsigned & next_weight) {
2169 bool modified = false;
2170 for (grobner::equation const* eq : eqs) {
2171 unsigned num_monomials = eq->get_num_monomials();
2172 CTRACE("grobner_bug", num_monomials <= 0, gb.display_equation(tout, *eq););
2173 if (num_monomials == 0)
2174 continue; // HACK: the equation 0 = 0, should have been discarded by the GB module.
2175 if (eq->get_monomial(0)->get_degree() != 1)
2176 continue;
2177 for (unsigned j = 1; j < num_monomials; j++) {
2178 grobner::monomial const * m = eq->get_monomial(j);
2179 if (m->get_degree() != 1)
2180 continue;
2181 expr * var = m->get_var(0);
2182 if (gb.get_weight(var) > MAX_DEFAULT_WEIGHT)
2183 continue; // variable was already updated
2184 TRACE("non_linear", tout << "increased weight of: " << mk_bounded_pp(var, get_manager()) << "\n";);
2185 gb.set_weight(var, next_weight);
2186 next_weight++;
2187 gb.update_order();
2188 TRACE("non_linear", tout << "after updating order\n"; gb.display(tout););
2189 modified = true;
2190 break;
2191 }
2192 if (modified)
2193 break;
2194 }
2195 return modified;
2196 }
2197
2198 // Scan the grobner basis eqs for equations of the form x - k = 0 or x = 0 is found, and x is not fixed,
2199 // then assert bounds for x, and continue
2200 template<typename Ext>
scan_for_linear(ptr_vector<grobner::equation> & eqs,grobner & gb)2201 bool theory_arith<Ext>::scan_for_linear(ptr_vector<grobner::equation>& eqs, grobner& gb) {
2202 bool result = false;
2203 if (m_params.m_nl_arith_gb_eqs) { // m_nl_arith_gb_eqs is false by default
2204 SASSERT(false);
2205 for (grobner::equation* eq : eqs) {
2206 if (!eq->is_linear_combination()) {
2207 TRACE("non_linear", tout << "processing new equality:\n"; gb.display_equation(tout, *eq););
2208 TRACE("non_linear_bug", tout << "processing new equality:\n"; gb.display_equation(tout, *eq););
2209 if (internalize_gb_eq(eq))
2210 result = true;
2211 }
2212 }
2213 }
2214 return result;
2215 }
2216
2217
2218 template<typename Ext>
compute_basis_loop(grobner & gb)2219 bool theory_arith<Ext>::compute_basis_loop(grobner & gb) {
2220 while (gb.get_num_new_equations() < m_params.m_nl_arith_gb_threshold) {
2221 if (ctx.get_cancel_flag()) {
2222 return false;
2223 }
2224 if (gb.compute_basis_step())
2225 return true;
2226 }
2227 return false;
2228 }
2229 template<typename Ext>
compute_basis(grobner & gb,bool & warn)2230 void theory_arith<Ext>::compute_basis(grobner& gb, bool& warn) {
2231 gb.compute_basis_init();
2232 if (!compute_basis_loop(gb) && !warn) {
2233 set_gb_exhausted();
2234 warn = true;
2235 }
2236 }
2237 /**
2238 \brief Compute Grobner basis, return true if a conflict or new fixed variables were detected.
2239 */
2240 template<typename Ext>
compute_grobner(svector<theory_var> const & nl_cluster)2241 typename theory_arith<Ext>::gb_result theory_arith<Ext>::compute_grobner(svector<theory_var> const & nl_cluster) {
2242 if (m_nl_gb_exhausted)
2243 return GB_FAIL;
2244 grobner gb(get_manager(), m_dep_manager);
2245 init_grobner(nl_cluster, gb);
2246 TRACE("non_linear", display(tout););
2247 bool warn = false;
2248 unsigned next_weight = MAX_DEFAULT_WEIGHT + 1; // next weight using during perturbation phase.
2249 ptr_vector<grobner::equation> eqs;
2250
2251 do {
2252 TRACE("non_linear_gb", tout << "before:\n"; gb.display(tout););
2253 compute_basis(gb, warn);
2254 update_statistics(gb);
2255 TRACE("non_linear_gb", tout << "after:\n"; gb.display(tout););
2256 if (ctx.get_cancel_flag())
2257 return GB_FAIL;
2258 if (get_gb_eqs_and_look_for_conflict(eqs, gb))
2259 return GB_PROGRESS;
2260 }
2261 while(scan_for_linear(eqs, gb) && m_params.m_nl_arith_gb_perturbate &&
2262 (!m_nl_gb_exhausted) && try_to_modify_eqs(eqs, gb, next_weight));
2263 return GB_FAIL;
2264 }
2265
2266 /**
2267 \brief Maximize/Minimize variables in non linear monomials.
2268 */
2269 template<typename Ext>
max_min_nl_vars()2270 bool theory_arith<Ext>::max_min_nl_vars() {
2271 var_set already_found;
2272 svector<theory_var> vars;
2273 for (theory_var v : m_nl_monomials) {
2274 mark_var(v, vars, already_found);
2275 expr * n = var2expr(v);
2276 SASSERT(is_pure_monomial(n));
2277 for (expr * curr : *to_app(n)) {
2278 if (ctx.e_internalized(curr)) {
2279 theory_var v = expr2var(curr);
2280 SASSERT(v != null_theory_var);
2281 mark_var(v, vars, already_found);
2282 }
2283 }
2284 }
2285 return max_min(vars);
2286 }
2287
2288 /**
2289 \brief Process non linear constraints.
2290 */
2291 template<typename Ext>
process_non_linear()2292 final_check_status theory_arith<Ext>::process_non_linear() {
2293 m_model_depends_on_computed_epsilon = false;
2294 if (m_nl_monomials.empty())
2295 return FC_DONE;
2296
2297 if (!reflection_enabled())
2298 return FC_GIVEUP;
2299
2300 if (check_monomial_assignments()) {
2301 return FC_DONE;
2302 }
2303
2304 if (!m_params.m_nl_arith) {
2305 TRACE("non_linear", tout << "Non-linear is not enabled\n";);
2306 return FC_GIVEUP;
2307 }
2308
2309 TRACE("process_non_linear", display(tout););
2310
2311 if (m_nl_rounds > m_params.m_nl_arith_rounds) {
2312 TRACE("non_linear", tout << "GIVEUP non linear problem...\n";);
2313 IF_VERBOSE(3, verbose_stream() << "Max. non linear arithmetic rounds. Increase threshold using NL_ARITH_ROUNDS=<limit>\n";);
2314 return FC_GIVEUP;
2315 }
2316
2317 ctx.push_trail(value_trail<unsigned>(m_nl_rounds));
2318 m_nl_rounds++;
2319
2320 elim_quasi_base_rows();
2321 move_non_base_vars_to_bounds();
2322 TRACE("non_linear_verbose", tout << "processing non linear constraints...\n"; ctx.display(tout););
2323 if (!make_feasible()) {
2324 TRACE("non_linear", tout << "failed to move variables to bounds.\n";);
2325 failed();
2326 return FC_CONTINUE;
2327 }
2328
2329 if (!max_min_nl_vars())
2330 return FC_CONTINUE;
2331
2332 if (check_monomial_assignments()) {
2333 return m_liberal_final_check || !m_changed_assignment ? FC_DONE : FC_CONTINUE;
2334 }
2335
2336 svector<theory_var> vars;
2337 get_non_linear_cluster(vars);
2338
2339 bool progress;
2340 unsigned old_idx = m_nl_strategy_idx;
2341 ctx.push_trail(value_trail<unsigned>(m_nl_strategy_idx));
2342
2343 do {
2344 progress = false;
2345 switch (m_nl_strategy_idx) {
2346 case 0:
2347 if (propagate_nl_bounds()) {
2348 propagate_core();
2349 progress = true;
2350 }
2351 break;
2352 case 1:
2353 if (!is_cross_nested_consistent(vars))
2354 progress = true;
2355 break;
2356 case 2:
2357 if (m_params.m_nl_arith_gb) {
2358 switch(compute_grobner(vars)) {
2359 case GB_PROGRESS:
2360 progress = true;
2361 break;
2362 case GB_NEW_EQ:
2363 progress = true;
2364 propagate_core();
2365 break;
2366 case GB_FAIL:
2367 break;
2368 }
2369 }
2370 break;
2371 case 3:
2372 if (m_params.m_nl_arith_branching) {
2373 theory_var target = find_nl_var_for_branching();
2374 if (target != null_theory_var && branch_nl_int_var(target))
2375 progress = true;
2376 }
2377 break;
2378 }
2379
2380 m_nl_strategy_idx = (m_nl_strategy_idx + 1) % 4;
2381 if (progress)
2382 return FC_CONTINUE;
2383 }
2384 while (m_nl_strategy_idx != old_idx);
2385
2386 if (check_monomial_assignments()) {
2387 return m_liberal_final_check || !m_changed_assignment ? FC_DONE : FC_CONTINUE;
2388 }
2389
2390 TRACE("non_linear", display(tout););
2391
2392 return FC_GIVEUP;
2393 }
2394
2395
2396 };
2397
2398
2399
2400