1 |
|
/****************************************************************************** |
2 |
|
* Top contributors (to current version): |
3 |
|
* Gereon Kremer, Aina Niemetz, Mathias Preiner |
4 |
|
* |
5 |
|
* This file is part of the cvc5 project. |
6 |
|
* |
7 |
|
* Copyright (c) 2009-2021 by the authors listed in the file AUTHORS |
8 |
|
* in the top-level source directory and their institutional affiliations. |
9 |
|
* All rights reserved. See the file COPYING in the top-level source |
10 |
|
* directory for licensing information. |
11 |
|
* **************************************************************************** |
12 |
|
* |
13 |
|
* Utilities for converting to and from LibPoly objects. |
14 |
|
*/ |
15 |
|
|
16 |
|
#include "poly_conversion.h" |
17 |
|
|
18 |
|
#ifdef CVC5_POLY_IMP |
19 |
|
|
20 |
|
#include "expr/node.h" |
21 |
|
#include "expr/node_manager_attributes.h" |
22 |
|
#include "theory/arith/bound_inference.h" |
23 |
|
#include "util/poly_util.h" |
24 |
|
|
25 |
|
namespace cvc5 { |
26 |
|
namespace theory { |
27 |
|
namespace arith { |
28 |
|
namespace nl { |
29 |
|
|
30 |
2339 |
poly::Variable VariableMapper::operator()(const cvc5::Node& n) |
31 |
|
{ |
32 |
2339 |
auto it = mVarCVCpoly.find(n); |
33 |
2339 |
if (it == mVarCVCpoly.end()) |
34 |
|
{ |
35 |
174 |
std::string name; |
36 |
87 |
if (n.isVar()) |
37 |
|
{ |
38 |
87 |
if (!n.getAttribute(expr::VarNameAttr(), name)) |
39 |
|
{ |
40 |
|
Trace("poly::conversion") |
41 |
|
<< "Variable " << n << " has no name, using ID instead." |
42 |
|
<< std::endl; |
43 |
|
name = "v_" + std::to_string(n.getId()); |
44 |
|
} |
45 |
|
} |
46 |
|
else |
47 |
|
{ |
48 |
|
name = "v_" + std::to_string(n.getId()); |
49 |
|
} |
50 |
87 |
it = mVarCVCpoly.emplace(n, poly::Variable(name.c_str())).first; |
51 |
87 |
mVarpolyCVC.emplace(it->second, n); |
52 |
|
} |
53 |
2339 |
return it->second; |
54 |
|
} |
55 |
|
|
56 |
248 |
cvc5::Node VariableMapper::operator()(const poly::Variable& n) |
57 |
|
{ |
58 |
248 |
auto it = mVarpolyCVC.find(n); |
59 |
248 |
Assert(it != mVarpolyCVC.end()) |
60 |
|
<< "Expect variable " << n << " to be added already."; |
61 |
248 |
return it->second; |
62 |
|
} |
63 |
|
|
64 |
9 |
cvc5::Node as_cvc_upolynomial(const poly::UPolynomial& p, const cvc5::Node& var) |
65 |
|
{ |
66 |
18 |
Trace("poly::conversion") |
67 |
9 |
<< "Converting " << p << " over " << var << std::endl; |
68 |
|
|
69 |
18 |
std::vector<poly::Integer> coeffs = coefficients(p); |
70 |
|
|
71 |
9 |
auto* nm = NodeManager::currentNM(); |
72 |
|
|
73 |
9 |
Node res = nm->mkConst(Rational(0)); |
74 |
18 |
Node monomial = nm->mkConst(Rational(1)); |
75 |
36 |
for (std::size_t i = 0, n = coeffs.size(); i < n; ++i) |
76 |
|
{ |
77 |
27 |
if (!is_zero(coeffs[i])) |
78 |
|
{ |
79 |
36 |
Node coeff = nm->mkConst(poly_utils::toRational(coeffs[i])); |
80 |
36 |
Node term = nm->mkNode(Kind::MULT, coeff, monomial); |
81 |
18 |
res = nm->mkNode(Kind::PLUS, res, term); |
82 |
|
} |
83 |
27 |
monomial = nm->mkNode(Kind::NONLINEAR_MULT, monomial, var); |
84 |
|
} |
85 |
9 |
Trace("poly::conversion") << "-> " << res << std::endl; |
86 |
18 |
return res; |
87 |
|
} |
88 |
|
|
89 |
26 |
poly::UPolynomial as_poly_upolynomial_impl(const cvc5::Node& n, |
90 |
|
poly::Integer& denominator, |
91 |
|
const cvc5::Node& var) |
92 |
|
{ |
93 |
26 |
denominator = poly::Integer(1); |
94 |
26 |
if (n.isVar()) |
95 |
|
{ |
96 |
4 |
Assert(n == var) << "Unexpected variable: should be " << var << " but is " |
97 |
|
<< n; |
98 |
4 |
return poly::UPolynomial({0, 1}); |
99 |
|
} |
100 |
22 |
switch (n.getKind()) |
101 |
|
{ |
102 |
10 |
case Kind::CONST_RATIONAL: |
103 |
|
{ |
104 |
20 |
Rational r = n.getConst<Rational>(); |
105 |
10 |
denominator = poly_utils::toInteger(r.getDenominator()); |
106 |
10 |
return poly::UPolynomial(poly_utils::toInteger(r.getNumerator())); |
107 |
|
} |
108 |
4 |
case Kind::PLUS: |
109 |
|
{ |
110 |
8 |
poly::UPolynomial res; |
111 |
8 |
poly::Integer denom; |
112 |
12 |
for (const auto& child : n) |
113 |
|
{ |
114 |
16 |
poly::UPolynomial tmp = as_poly_upolynomial_impl(child, denom, var); |
115 |
|
/** Normalize denominators |
116 |
|
*/ |
117 |
16 |
poly::Integer g = gcd(denom, denominator); |
118 |
8 |
res = res * (denom / g) + tmp * (denominator / g); |
119 |
8 |
denominator *= (denom / g); |
120 |
|
} |
121 |
4 |
return res; |
122 |
|
} |
123 |
8 |
case Kind::MULT: |
124 |
|
case Kind::NONLINEAR_MULT: |
125 |
|
{ |
126 |
16 |
poly::UPolynomial res(denominator); |
127 |
16 |
poly::Integer denom; |
128 |
24 |
for (const auto& child : n) |
129 |
|
{ |
130 |
16 |
res = res * as_poly_upolynomial_impl(child, denom, var); |
131 |
16 |
denominator *= denom; |
132 |
|
} |
133 |
8 |
return res; |
134 |
|
} |
135 |
|
default: |
136 |
|
Warning() << "Unhandled node " << n << " with kind " << n.getKind() |
137 |
|
<< std::endl; |
138 |
|
} |
139 |
|
return poly::UPolynomial(); |
140 |
|
} |
141 |
|
|
142 |
2 |
poly::UPolynomial as_poly_upolynomial(const cvc5::Node& n, |
143 |
|
const cvc5::Node& var) |
144 |
|
{ |
145 |
4 |
poly::Integer denom; |
146 |
4 |
return as_poly_upolynomial_impl(n, denom, var); |
147 |
|
} |
148 |
|
|
149 |
4603 |
poly::Polynomial as_poly_polynomial_impl(const cvc5::Node& n, |
150 |
|
poly::Integer& denominator, |
151 |
|
VariableMapper& vm) |
152 |
|
{ |
153 |
4603 |
denominator = poly::Integer(1); |
154 |
4603 |
if (n.isVar()) |
155 |
|
{ |
156 |
2197 |
return poly::Polynomial(vm(n)); |
157 |
|
} |
158 |
2406 |
switch (n.getKind()) |
159 |
|
{ |
160 |
1224 |
case Kind::CONST_RATIONAL: |
161 |
|
{ |
162 |
2448 |
Rational r = n.getConst<Rational>(); |
163 |
1224 |
denominator = poly_utils::toInteger(r.getDenominator()); |
164 |
1224 |
return poly::Polynomial(poly_utils::toInteger(r.getNumerator())); |
165 |
|
} |
166 |
256 |
case Kind::PLUS: |
167 |
|
{ |
168 |
512 |
poly::Polynomial res; |
169 |
512 |
poly::Integer denom; |
170 |
813 |
for (const auto& child : n) |
171 |
|
{ |
172 |
1114 |
poly::Polynomial tmp = as_poly_polynomial_impl(child, denom, vm); |
173 |
|
/** Normalize denominators |
174 |
|
*/ |
175 |
1114 |
poly::Integer g = gcd(denom, denominator); |
176 |
557 |
res = res * (denom / g) + tmp * (denominator / g); |
177 |
557 |
denominator *= (denom / g); |
178 |
|
} |
179 |
256 |
return res; |
180 |
|
} |
181 |
926 |
case Kind::MULT: |
182 |
|
case Kind::NONLINEAR_MULT: |
183 |
|
{ |
184 |
1852 |
poly::Polynomial res(denominator); |
185 |
1852 |
poly::Integer denom; |
186 |
2954 |
for (const auto& child : n) |
187 |
|
{ |
188 |
2028 |
res *= as_poly_polynomial_impl(child, denom, vm); |
189 |
2028 |
denominator *= denom; |
190 |
|
} |
191 |
926 |
return res; |
192 |
|
} |
193 |
|
default: return poly::Polynomial(vm(n)); |
194 |
|
} |
195 |
|
return poly::Polynomial(); |
196 |
|
} |
197 |
|
poly::Polynomial as_poly_polynomial(const cvc5::Node& n, VariableMapper& vm) |
198 |
|
{ |
199 |
|
poly::Integer denom; |
200 |
|
return as_poly_polynomial_impl(n, denom, vm); |
201 |
|
} |
202 |
30 |
poly::Polynomial as_poly_polynomial(const cvc5::Node& n, |
203 |
|
VariableMapper& vm, |
204 |
|
poly::Rational& denominator) |
205 |
|
{ |
206 |
60 |
poly::Integer denom; |
207 |
30 |
auto res = as_poly_polynomial_impl(n, denom, vm); |
208 |
30 |
denominator = poly::Rational(denom); |
209 |
60 |
return res; |
210 |
|
} |
211 |
|
|
212 |
|
namespace { |
213 |
|
/** |
214 |
|
* Utility class that collects the monomial terms (as nodes) from the polynomial |
215 |
|
* we are converting. |
216 |
|
*/ |
217 |
5 |
struct CollectMonomialData |
218 |
|
{ |
219 |
10 |
CollectMonomialData(VariableMapper& v) : d_vm(v) {} |
220 |
|
|
221 |
|
/** Mapper from poly variables to cvc5 variables */ |
222 |
|
VariableMapper& d_vm; |
223 |
|
/** Collections of the monomial terms */ |
224 |
|
std::vector<Node> d_terms; |
225 |
|
/** Caches the current node manager */ |
226 |
5 |
NodeManager* d_nm = NodeManager::currentNM(); |
227 |
|
}; |
228 |
|
/** |
229 |
|
* Callback for lp_polynomial_traverse. Assumes data is actually a |
230 |
|
* CollectMonomialData object and puts the polynomial into it. |
231 |
|
*/ |
232 |
5 |
void collect_monomials(const lp_polynomial_context_t* ctx, |
233 |
|
lp_monomial_t* m, |
234 |
|
void* data) |
235 |
|
{ |
236 |
5 |
CollectMonomialData* d = static_cast<CollectMonomialData*>(data); |
237 |
|
// constant |
238 |
|
Node term = |
239 |
10 |
d->d_nm->mkConst<Rational>(poly_utils::toRational(poly::Integer(&m->a))); |
240 |
12 |
for (std::size_t i = 0; i < m->n; ++i) |
241 |
|
{ |
242 |
|
// variable exponent pair |
243 |
14 |
Node var = d->d_vm(m->p[i].x); |
244 |
7 |
if (m->p[i].d > 1) |
245 |
|
{ |
246 |
|
Node exp = d->d_nm->mkConst<Rational>(m->p[i].d); |
247 |
|
term = d->d_nm->mkNode( |
248 |
|
Kind::NONLINEAR_MULT, term, d->d_nm->mkNode(Kind::POW, var, exp)); |
249 |
|
} |
250 |
|
else |
251 |
|
{ |
252 |
7 |
term = d->d_nm->mkNode(Kind::NONLINEAR_MULT, term, var); |
253 |
|
} |
254 |
|
} |
255 |
5 |
d->d_terms.emplace_back(term); |
256 |
5 |
} |
257 |
|
} // namespace |
258 |
|
|
259 |
5 |
cvc5::Node as_cvc_polynomial(const poly::Polynomial& p, VariableMapper& vm) |
260 |
|
{ |
261 |
10 |
CollectMonomialData cmd(vm); |
262 |
|
// Do the actual conversion |
263 |
5 |
lp_polynomial_traverse(p.get_internal(), collect_monomials, &cmd); |
264 |
|
|
265 |
5 |
if (cmd.d_terms.empty()) |
266 |
|
{ |
267 |
|
return cmd.d_nm->mkConst<Rational>(0); |
268 |
|
} |
269 |
5 |
if (cmd.d_terms.size() == 1) |
270 |
|
{ |
271 |
5 |
return cmd.d_terms.front(); |
272 |
|
} |
273 |
|
return cmd.d_nm->mkNode(Kind::PLUS, cmd.d_terms); |
274 |
|
} |
275 |
|
|
276 |
994 |
poly::SignCondition normalize_kind(cvc5::Kind kind, |
277 |
|
bool negated, |
278 |
|
poly::Polynomial& lhs) |
279 |
|
{ |
280 |
994 |
switch (kind) |
281 |
|
{ |
282 |
750 |
case Kind::EQUAL: |
283 |
|
{ |
284 |
750 |
return negated ? poly::SignCondition::NE : poly::SignCondition::EQ; |
285 |
|
} |
286 |
|
case Kind::LT: |
287 |
|
{ |
288 |
|
if (negated) |
289 |
|
{ |
290 |
|
lhs = -lhs; |
291 |
|
return poly::SignCondition::LE; |
292 |
|
} |
293 |
|
return poly::SignCondition::LT; |
294 |
|
} |
295 |
|
case Kind::LEQ: |
296 |
|
{ |
297 |
|
if (negated) |
298 |
|
{ |
299 |
|
lhs = -lhs; |
300 |
|
return poly::SignCondition::LT; |
301 |
|
} |
302 |
|
return poly::SignCondition::LE; |
303 |
|
} |
304 |
|
case Kind::GT: |
305 |
|
{ |
306 |
|
if (negated) |
307 |
|
{ |
308 |
|
return poly::SignCondition::LE; |
309 |
|
} |
310 |
|
lhs = -lhs; |
311 |
|
return poly::SignCondition::LT; |
312 |
|
} |
313 |
244 |
case Kind::GEQ: |
314 |
|
{ |
315 |
244 |
if (negated) |
316 |
|
{ |
317 |
224 |
return poly::SignCondition::LT; |
318 |
|
} |
319 |
20 |
lhs = -lhs; |
320 |
20 |
return poly::SignCondition::LE; |
321 |
|
} |
322 |
|
default: |
323 |
|
Assert(false) << "This function only deals with arithmetic relations."; |
324 |
|
return poly::SignCondition::EQ; |
325 |
|
} |
326 |
|
} |
327 |
|
|
328 |
994 |
std::pair<poly::Polynomial, poly::SignCondition> as_poly_constraint( |
329 |
|
Node n, VariableMapper& vm) |
330 |
|
{ |
331 |
994 |
bool negated = false; |
332 |
1988 |
Node origin = n; |
333 |
994 |
if (n.getKind() == Kind::NOT) |
334 |
|
{ |
335 |
850 |
Assert(n.getNumChildren() == 1) |
336 |
|
<< "Expect negations to have a single child."; |
337 |
850 |
negated = true; |
338 |
850 |
n = *n.begin(); |
339 |
|
} |
340 |
994 |
Assert((n.getKind() == Kind::EQUAL) || (n.getKind() == Kind::GT) |
341 |
|
|| (n.getKind() == Kind::GEQ) || (n.getKind() == Kind::LT) |
342 |
|
|| (n.getKind() == Kind::LEQ)) |
343 |
|
<< "Found a constraint with unsupported relation " << n.getKind(); |
344 |
|
|
345 |
994 |
Assert(n.getNumChildren() == 2) |
346 |
|
<< "Supported relations only have two children."; |
347 |
994 |
auto childit = n.begin(); |
348 |
1988 |
poly::Integer ldenom; |
349 |
1988 |
poly::Polynomial left = as_poly_polynomial_impl(*childit++, ldenom, vm); |
350 |
1988 |
poly::Integer rdenom; |
351 |
1988 |
poly::Polynomial right = as_poly_polynomial_impl(*childit++, rdenom, vm); |
352 |
994 |
Assert(childit == n.end()) << "Screwed up iterator handling."; |
353 |
994 |
Assert(ldenom > poly::Integer(0) && rdenom > poly::Integer(0)) |
354 |
|
<< "Expected denominators to be always positive."; |
355 |
|
|
356 |
1988 |
poly::Integer g = gcd(ldenom, rdenom); |
357 |
1988 |
poly::Polynomial lhs = left * (rdenom / g) - right * (ldenom / g); |
358 |
994 |
poly::SignCondition sc = normalize_kind(n.getKind(), negated, lhs); |
359 |
1988 |
return {lhs, sc}; |
360 |
|
} |
361 |
|
|
362 |
2 |
Node ran_to_node(const RealAlgebraicNumber& ran, const Node& ran_variable) |
363 |
|
{ |
364 |
2 |
return ran_to_node(ran.getValue(), ran_variable); |
365 |
|
} |
366 |
|
|
367 |
12 |
Node ran_to_node(const poly::AlgebraicNumber& an, const Node& ran_variable) |
368 |
|
{ |
369 |
12 |
auto* nm = NodeManager::currentNM(); |
370 |
|
|
371 |
12 |
const poly::DyadicInterval& di = get_isolating_interval(an); |
372 |
12 |
if (is_point(di)) |
373 |
|
{ |
374 |
3 |
return nm->mkConst(poly_utils::toRational(get_point(di))); |
375 |
|
} |
376 |
9 |
Assert(di.get_internal()->a_open && di.get_internal()->b_open) |
377 |
|
<< "We assume an open interval here."; |
378 |
|
|
379 |
18 |
Node poly = as_cvc_upolynomial(get_defining_polynomial(an), ran_variable); |
380 |
18 |
Node lower = nm->mkConst(poly_utils::toRational(get_lower(di))); |
381 |
18 |
Node upper = nm->mkConst(poly_utils::toRational(get_upper(di))); |
382 |
|
|
383 |
|
// Construct witness: |
384 |
|
return nm->mkNode(Kind::AND, |
385 |
|
// poly(var) == 0 |
386 |
18 |
nm->mkNode(Kind::EQUAL, poly, nm->mkConst(Rational(0))), |
387 |
|
// lower_bound < var |
388 |
18 |
nm->mkNode(Kind::LT, lower, ran_variable), |
389 |
|
// var < upper_bound |
390 |
45 |
nm->mkNode(Kind::LT, ran_variable, upper)); |
391 |
|
} |
392 |
|
|
393 |
53 |
Node value_to_node(const poly::Value& v, const Node& ran_variable) |
394 |
|
{ |
395 |
53 |
Assert(!is_minus_infinity(v)) << "Can not convert minus infinity."; |
396 |
53 |
Assert(!is_none(v)) << "Can not convert none."; |
397 |
53 |
Assert(!is_plus_infinity(v)) << "Can not convert plus infinity."; |
398 |
|
|
399 |
53 |
if (is_algebraic_number(v)) |
400 |
|
{ |
401 |
10 |
return ran_to_node(as_algebraic_number(v), ran_variable); |
402 |
|
} |
403 |
43 |
auto* nm = NodeManager::currentNM(); |
404 |
43 |
if (is_dyadic_rational(v)) |
405 |
|
{ |
406 |
|
return nm->mkConst(poly_utils::toRational(as_dyadic_rational(v))); |
407 |
|
} |
408 |
43 |
if (is_integer(v)) |
409 |
|
{ |
410 |
7 |
return nm->mkConst(poly_utils::toRational(as_integer(v))); |
411 |
|
} |
412 |
36 |
if (is_rational(v)) |
413 |
|
{ |
414 |
36 |
return nm->mkConst(poly_utils::toRational(as_rational(v))); |
415 |
|
} |
416 |
|
Assert(false) << "All cases should be covered."; |
417 |
|
return nm->mkConst(Rational(0)); |
418 |
|
} |
419 |
|
|
420 |
|
Node lower_bound_as_node(const Node& var, |
421 |
|
const poly::Value& lower, |
422 |
|
bool open, |
423 |
|
bool allowNonlinearLemma) |
424 |
|
{ |
425 |
|
auto* nm = NodeManager::currentNM(); |
426 |
|
if (!poly::is_algebraic_number(lower)) |
427 |
|
{ |
428 |
|
return nm->mkNode(open ? Kind::LEQ : Kind::LT, |
429 |
|
var, |
430 |
|
nm->mkConst(poly_utils::toRationalAbove(lower))); |
431 |
|
} |
432 |
|
if (poly::represents_rational(lower)) |
433 |
|
{ |
434 |
|
return nm->mkNode( |
435 |
|
open ? Kind::LEQ : Kind::LT, |
436 |
|
var, |
437 |
|
nm->mkConst(poly_utils::toRationalAbove(poly::get_rational(lower)))); |
438 |
|
} |
439 |
|
if (!allowNonlinearLemma) |
440 |
|
{ |
441 |
|
return Node(); |
442 |
|
} |
443 |
|
|
444 |
|
const poly::AlgebraicNumber& alg = as_algebraic_number(lower); |
445 |
|
|
446 |
|
Node poly = as_cvc_upolynomial(get_defining_polynomial(alg), var); |
447 |
|
Rational l = poly_utils::toRational( |
448 |
|
poly::get_lower(poly::get_isolating_interval(alg))); |
449 |
|
Rational u = poly_utils::toRational( |
450 |
|
poly::get_upper(poly::get_isolating_interval(alg))); |
451 |
|
int sl = poly::sign_at(get_defining_polynomial(alg), |
452 |
|
poly::get_lower(poly::get_isolating_interval(alg))); |
453 |
|
#ifdef CVC5_ASSERTIONS |
454 |
|
int su = poly::sign_at(get_defining_polynomial(alg), |
455 |
|
poly::get_upper(poly::get_isolating_interval(alg))); |
456 |
|
Assert(sl != 0 && su != 0 && sl != su); |
457 |
|
#endif |
458 |
|
|
459 |
|
// open: var <= l or (var < u and sgn(poly(var)) == sl) |
460 |
|
// !open: var <= l or (var < u and sgn(poly(var)) == sl/0) |
461 |
|
Kind relation; |
462 |
|
if (open) |
463 |
|
{ |
464 |
|
relation = (sl < 0) ? Kind::LEQ : Kind::GEQ; |
465 |
|
} |
466 |
|
else |
467 |
|
{ |
468 |
|
relation = (sl < 0) ? Kind::LT : Kind::GT; |
469 |
|
} |
470 |
|
return nm->mkNode( |
471 |
|
Kind::OR, |
472 |
|
nm->mkNode(Kind::LEQ, var, nm->mkConst(l)), |
473 |
|
nm->mkNode(Kind::AND, |
474 |
|
nm->mkNode(Kind::LT, var, nm->mkConst(u)), |
475 |
|
nm->mkNode(relation, poly, nm->mkConst(Rational(0))))); |
476 |
|
} |
477 |
|
|
478 |
|
Node upper_bound_as_node(const Node& var, |
479 |
|
const poly::Value& upper, |
480 |
|
bool open, |
481 |
|
bool allowNonlinearLemma) |
482 |
|
{ |
483 |
|
auto* nm = NodeManager::currentNM(); |
484 |
|
if (!poly::is_algebraic_number(upper)) |
485 |
|
{ |
486 |
|
return nm->mkNode(open ? Kind::GEQ : Kind::GT, |
487 |
|
var, |
488 |
|
nm->mkConst(poly_utils::toRationalAbove(upper))); |
489 |
|
} |
490 |
|
if (poly::represents_rational(upper)) |
491 |
|
{ |
492 |
|
return nm->mkNode( |
493 |
|
open ? Kind::GEQ : Kind::GT, |
494 |
|
var, |
495 |
|
nm->mkConst(poly_utils::toRationalAbove(poly::get_rational(upper)))); |
496 |
|
} |
497 |
|
if (!allowNonlinearLemma) |
498 |
|
{ |
499 |
|
return Node(); |
500 |
|
} |
501 |
|
|
502 |
|
const poly::AlgebraicNumber& alg = as_algebraic_number(upper); |
503 |
|
|
504 |
|
Node poly = as_cvc_upolynomial(get_defining_polynomial(alg), var); |
505 |
|
Rational l = poly_utils::toRational( |
506 |
|
poly::get_lower(poly::get_isolating_interval(alg))); |
507 |
|
Rational u = poly_utils::toRational( |
508 |
|
poly::get_upper(poly::get_isolating_interval(alg))); |
509 |
|
#ifdef CVC5_ASSERTIONS |
510 |
|
int sl = poly::sign_at(get_defining_polynomial(alg), |
511 |
|
poly::get_lower(poly::get_isolating_interval(alg))); |
512 |
|
#endif |
513 |
|
int su = poly::sign_at(get_defining_polynomial(alg), |
514 |
|
poly::get_upper(poly::get_isolating_interval(alg))); |
515 |
|
Assert(sl != 0 && su != 0 && sl != su); |
516 |
|
|
517 |
|
// open: var >= u or (var > l and sgn(poly(var)) == su) |
518 |
|
// !open: var >= u or (var > l and sgn(poly(var)) == su/0) |
519 |
|
Kind relation; |
520 |
|
if (open) |
521 |
|
{ |
522 |
|
relation = (su < 0) ? Kind::LEQ : Kind::GEQ; |
523 |
|
} |
524 |
|
else |
525 |
|
{ |
526 |
|
relation = (su < 0) ? Kind::LT : Kind::GT; |
527 |
|
} |
528 |
|
return nm->mkNode( |
529 |
|
Kind::OR, |
530 |
|
nm->mkNode(Kind::GEQ, var, nm->mkConst(u)), |
531 |
|
nm->mkNode(Kind::AND, |
532 |
|
nm->mkNode(Kind::GT, var, nm->mkConst(l)), |
533 |
|
nm->mkNode(relation, poly, nm->mkConst(Rational(0))))); |
534 |
|
} |
535 |
|
|
536 |
|
Node excluding_interval_to_lemma(const Node& variable, |
537 |
|
const poly::Interval& interval, |
538 |
|
bool allowNonlinearLemma) |
539 |
|
{ |
540 |
|
auto* nm = NodeManager::currentNM(); |
541 |
|
const auto& lv = poly::get_lower(interval); |
542 |
|
const auto& uv = poly::get_upper(interval); |
543 |
|
if (bitsize(lv) > 100 || bitsize(uv) > 100) return Node(); |
544 |
|
bool li = poly::is_minus_infinity(lv); |
545 |
|
bool ui = poly::is_plus_infinity(uv); |
546 |
|
if (li && ui) return nm->mkConst(true); |
547 |
|
if (poly::is_point(interval)) |
548 |
|
{ |
549 |
|
if (is_algebraic_number(lv)) |
550 |
|
{ |
551 |
|
const poly::AlgebraicNumber& alg = as_algebraic_number(lv); |
552 |
|
if (poly::is_rational(alg)) |
553 |
|
{ |
554 |
|
Trace("nl-cad") << "Rational point interval: " << interval << std::endl; |
555 |
|
return nm->mkNode(Kind::DISTINCT, |
556 |
|
variable, |
557 |
|
nm->mkConst(poly_utils::toRational( |
558 |
|
poly::to_rational_approximation(alg)))); |
559 |
|
} |
560 |
|
Trace("nl-cad") << "Algebraic point interval: " << interval << std::endl; |
561 |
|
// p(x) != 0 or x <= lb or ub <= x |
562 |
|
if (allowNonlinearLemma) |
563 |
|
{ |
564 |
|
Node poly = as_cvc_upolynomial(get_defining_polynomial(alg), variable); |
565 |
|
return nm->mkNode( |
566 |
|
Kind::OR, |
567 |
|
nm->mkNode(Kind::DISTINCT, poly, nm->mkConst(Rational(0))), |
568 |
|
nm->mkNode(Kind::LT, |
569 |
|
variable, |
570 |
|
nm->mkConst(poly_utils::toRationalBelow(lv))), |
571 |
|
nm->mkNode(Kind::GT, |
572 |
|
variable, |
573 |
|
nm->mkConst(poly_utils::toRationalAbove(lv)))); |
574 |
|
} |
575 |
|
return Node(); |
576 |
|
} |
577 |
|
else |
578 |
|
{ |
579 |
|
Trace("nl-cad") << "Rational point interval: " << interval << std::endl; |
580 |
|
return nm->mkNode(Kind::DISTINCT, |
581 |
|
variable, |
582 |
|
nm->mkConst(poly_utils::toRationalBelow(lv))); |
583 |
|
} |
584 |
|
} |
585 |
|
if (li) |
586 |
|
{ |
587 |
|
Trace("nl-cad") << "Only upper bound: " << interval << std::endl; |
588 |
|
return upper_bound_as_node( |
589 |
|
variable, uv, poly::get_upper_open(interval), allowNonlinearLemma); |
590 |
|
} |
591 |
|
if (ui) |
592 |
|
{ |
593 |
|
Trace("nl-cad") << "Only lower bound: " << interval << std::endl; |
594 |
|
return lower_bound_as_node( |
595 |
|
variable, lv, poly::get_lower_open(interval), allowNonlinearLemma); |
596 |
|
} |
597 |
|
Trace("nl-cad") << "Proper interval: " << interval << std::endl; |
598 |
|
Node lb = lower_bound_as_node( |
599 |
|
variable, lv, poly::get_lower_open(interval), allowNonlinearLemma); |
600 |
|
Node ub = upper_bound_as_node( |
601 |
|
variable, uv, poly::get_upper_open(interval), allowNonlinearLemma); |
602 |
|
if (lb.isNull() || ub.isNull()) return Node(); |
603 |
|
return nm->mkNode(Kind::OR, lb, ub); |
604 |
|
} |
605 |
|
|
606 |
4 |
std::optional<Rational> get_lower_bound(const Node& n) |
607 |
|
{ |
608 |
4 |
if (n.getNumChildren() != 2) return std::optional<Rational>(); |
609 |
4 |
if (n.getKind() == Kind::LT) |
610 |
|
{ |
611 |
2 |
if (!n[0].isConst()) return std::optional<Rational>(); |
612 |
2 |
if (!n[1].isVar()) return std::optional<Rational>(); |
613 |
2 |
return n[0].getConst<Rational>(); |
614 |
|
} |
615 |
2 |
else if (n.getKind() == Kind::GT) |
616 |
|
{ |
617 |
|
if (!n[0].isVar()) return std::optional<Rational>(); |
618 |
|
if (!n[1].isConst()) return std::optional<Rational>(); |
619 |
|
return n[1].getConst<Rational>(); |
620 |
|
} |
621 |
2 |
return std::optional<Rational>(); |
622 |
|
} |
623 |
6 |
std::optional<Rational> get_upper_bound(const Node& n) |
624 |
|
{ |
625 |
6 |
if (n.getNumChildren() != 2) return std::optional<Rational>(); |
626 |
6 |
if (n.getKind() == Kind::LT) |
627 |
|
{ |
628 |
4 |
if (!n[0].isVar()) return std::optional<Rational>(); |
629 |
2 |
if (!n[1].isConst()) return std::optional<Rational>(); |
630 |
2 |
return n[1].getConst<Rational>(); |
631 |
|
} |
632 |
2 |
else if (n.getKind() == Kind::GT) |
633 |
|
{ |
634 |
|
if (!n[0].isConst()) return std::optional<Rational>(); |
635 |
|
if (!n[1].isVar()) return std::optional<Rational>(); |
636 |
|
return n[0].getConst<Rational>(); |
637 |
|
} |
638 |
2 |
return std::optional<Rational>(); |
639 |
|
} |
640 |
|
|
641 |
|
/** Returns indices of appropriate parts of ran encoding. |
642 |
|
* Returns (poly equation ; lower bound ; upper bound) |
643 |
|
*/ |
644 |
2 |
std::tuple<Node, Rational, Rational> detect_ran_encoding(const Node& n) |
645 |
|
{ |
646 |
2 |
Assert(n.getKind() == Kind::AND) << "Invalid node structure."; |
647 |
2 |
Assert(n.getNumChildren() == 3) << "Invalid node structure."; |
648 |
|
|
649 |
4 |
Node poly_eq; |
650 |
2 |
if (n[0].getKind() == Kind::EQUAL) |
651 |
2 |
poly_eq = n[0]; |
652 |
|
else if (n[1].getKind() == Kind::EQUAL) |
653 |
|
poly_eq = n[1]; |
654 |
|
else if (n[2].getKind() == Kind::EQUAL) |
655 |
|
poly_eq = n[2]; |
656 |
|
else |
657 |
|
Assert(false) << "Could not identify polynomial equation."; |
658 |
|
|
659 |
4 |
Node poly; |
660 |
2 |
Assert(poly_eq.getNumChildren() == 2) << "Invalid polynomial equation."; |
661 |
2 |
if (poly_eq[0].isConst()) |
662 |
|
{ |
663 |
|
Assert(poly_eq[0].getConst<Rational>() == Rational(0)) |
664 |
|
<< "Invalid polynomial equation."; |
665 |
|
poly = poly_eq[1]; |
666 |
|
} |
667 |
2 |
else if (poly_eq[1].isConst()) |
668 |
|
{ |
669 |
2 |
Assert(poly_eq[1].getConst<Rational>() == Rational(0)) |
670 |
|
<< "Invalid polynomial equation."; |
671 |
2 |
poly = poly_eq[0]; |
672 |
|
} |
673 |
|
else |
674 |
|
{ |
675 |
|
Assert(false) << "Invalid polynomial equation."; |
676 |
|
} |
677 |
|
|
678 |
4 |
std::optional<Rational> lower = get_lower_bound(n[0]); |
679 |
2 |
if (!lower) lower = get_lower_bound(n[1]); |
680 |
2 |
if (!lower) lower = get_lower_bound(n[2]); |
681 |
2 |
Assert(lower) << "Could not identify lower bound."; |
682 |
|
|
683 |
4 |
std::optional<Rational> upper = get_upper_bound(n[0]); |
684 |
2 |
if (!upper) upper = get_upper_bound(n[1]); |
685 |
2 |
if (!upper) upper = get_upper_bound(n[2]); |
686 |
2 |
Assert(upper) << "Could not identify upper bound."; |
687 |
|
|
688 |
|
return std::tuple<Node, Rational, Rational>( |
689 |
4 |
poly, lower.value(), upper.value()); |
690 |
|
} |
691 |
|
|
692 |
2 |
poly::AlgebraicNumber node_to_poly_ran(const Node& n, const Node& ran_variable) |
693 |
|
{ |
694 |
|
// Identify poly, lower and upper |
695 |
4 |
auto encoding = detect_ran_encoding(n); |
696 |
|
// Construct polynomial |
697 |
|
poly::UPolynomial pol = |
698 |
4 |
as_poly_upolynomial(std::get<0>(encoding), ran_variable); |
699 |
|
// Construct algebraic number |
700 |
|
return poly_utils::toPolyRanWithRefinement( |
701 |
4 |
std::move(pol), std::get<1>(encoding), std::get<2>(encoding)); |
702 |
|
} |
703 |
2 |
RealAlgebraicNumber node_to_ran(const Node& n, const Node& ran_variable) |
704 |
|
{ |
705 |
2 |
return RealAlgebraicNumber(node_to_poly_ran(n, ran_variable)); |
706 |
|
} |
707 |
|
|
708 |
152 |
poly::Value node_to_value(const Node& n, const Node& ran_variable) |
709 |
|
{ |
710 |
152 |
if (n.isConst()) |
711 |
|
{ |
712 |
152 |
return poly_utils::toRational(n.getConst<Rational>()); |
713 |
|
} |
714 |
|
return node_to_poly_ran(n, ran_variable); |
715 |
|
} |
716 |
|
|
717 |
|
/** Bitsize of a dyadic rational. */ |
718 |
|
std::size_t bitsize(const poly::DyadicRational& v) |
719 |
|
{ |
720 |
|
return bit_size(numerator(v)) + bit_size(denominator(v)); |
721 |
|
} |
722 |
|
/** Bitsize of an integer. */ |
723 |
|
std::size_t bitsize(const poly::Integer& v) { return bit_size(v); } |
724 |
|
/** Bitsize of a rational. */ |
725 |
168 |
std::size_t bitsize(const poly::Rational& v) |
726 |
|
{ |
727 |
168 |
return bit_size(numerator(v)) + bit_size(denominator(v)); |
728 |
|
} |
729 |
|
/** Bitsize of a univariate polynomial. */ |
730 |
|
std::size_t bitsize(const poly::UPolynomial& v) |
731 |
|
{ |
732 |
|
std::size_t sum = 0; |
733 |
|
for (const auto& c : coefficients(v)) |
734 |
|
{ |
735 |
|
sum += bitsize(c); |
736 |
|
} |
737 |
|
return sum; |
738 |
|
} |
739 |
|
/** Bitsize of an algebraic number. */ |
740 |
|
std::size_t bitsize(const poly::AlgebraicNumber& v) |
741 |
|
{ |
742 |
|
if (is_rational(v)) |
743 |
|
{ |
744 |
|
return bitsize(to_rational_approximation(v)); |
745 |
|
} |
746 |
|
return bitsize(get_lower_bound(v)) + bitsize(get_upper_bound(v)) |
747 |
|
+ bitsize(get_defining_polynomial(v)); |
748 |
|
} |
749 |
|
|
750 |
208 |
std::size_t bitsize(const poly::Value& v) |
751 |
|
{ |
752 |
208 |
if (is_algebraic_number(v)) |
753 |
|
{ |
754 |
|
return bitsize(as_algebraic_number(v)); |
755 |
|
} |
756 |
208 |
else if (is_dyadic_rational(v)) |
757 |
|
{ |
758 |
|
return bitsize(as_dyadic_rational(v)); |
759 |
|
} |
760 |
208 |
else if (is_integer(v)) |
761 |
|
{ |
762 |
|
return bitsize(as_integer(v)); |
763 |
|
} |
764 |
208 |
else if (is_minus_infinity(v)) |
765 |
|
{ |
766 |
22 |
return 1; |
767 |
|
} |
768 |
186 |
else if (is_none(v)) |
769 |
|
{ |
770 |
|
return 0; |
771 |
|
} |
772 |
186 |
else if (is_plus_infinity(v)) |
773 |
|
{ |
774 |
18 |
return 1; |
775 |
|
} |
776 |
168 |
else if (is_rational(v)) |
777 |
|
{ |
778 |
168 |
return bitsize(as_rational(v)); |
779 |
|
} |
780 |
|
Assert(false); |
781 |
|
return 0; |
782 |
|
} |
783 |
|
|
784 |
48 |
poly::IntervalAssignment getBounds(VariableMapper& vm, const BoundInference& bi) |
785 |
|
{ |
786 |
48 |
poly::IntervalAssignment res; |
787 |
140 |
for (const auto& vb : bi.get()) |
788 |
|
{ |
789 |
92 |
poly::Variable v = vm(vb.first); |
790 |
92 |
poly::Value l = vb.second.lower_value.isNull() |
791 |
|
? poly::Value::minus_infty() |
792 |
184 |
: node_to_value(vb.second.lower_value, vb.first); |
793 |
92 |
poly::Value u = vb.second.upper_value.isNull() |
794 |
|
? poly::Value::plus_infty() |
795 |
184 |
: node_to_value(vb.second.upper_value, vb.first); |
796 |
184 |
poly::Interval i(l, vb.second.lower_strict, u, vb.second.upper_strict); |
797 |
92 |
res.set(v, i); |
798 |
|
} |
799 |
48 |
return res; |
800 |
|
} |
801 |
|
|
802 |
|
} // namespace nl |
803 |
|
} // namespace arith |
804 |
|
} // namespace theory |
805 |
29508 |
} // namespace cvc5 |
806 |
|
|
807 |
|
#endif |