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