Branch data Line data Source code
1 : : /******************************************************************************
2 : : * This file is part of the cvc5 project.
3 : : *
4 : : * Copyright (c) 2009-2026 by the authors listed in the file AUTHORS
5 : : * in the top-level source directory and their institutional affiliations.
6 : : * All rights reserved. See the file COPYING in the top-level source
7 : : * directory for licensing information.
8 : : * ****************************************************************************
9 : : *
10 : : * Utilities for working with LibPoly.
11 : : */
12 : :
13 : : #include "poly_util.h"
14 : :
15 : : #ifdef CVC5_POLY_IMP
16 : :
17 : : #include <poly/polyxx.h>
18 : :
19 : : #include <map>
20 : : #include <sstream>
21 : :
22 : : #include "base/check.h"
23 : : #include "util/integer.h"
24 : : #include "util/rational.h"
25 : : #include "util/real_algebraic_number.h"
26 : :
27 : : namespace cvc5::internal {
28 : : namespace poly_utils {
29 : :
30 : : namespace {
31 : : /**
32 : : * Convert arbitrary data using a string as intermediary.
33 : : * Assumes the existence of operator<<(std::ostream&, const From&) and To(const
34 : : * std::string&); should be the last resort for type conversions: it may not
35 : : * only yield bad performance, but is also dependent on compatible string
36 : : * representations. Use with care!
37 : : *
38 : : * Only instantiated by the CLN branches below; GMP builds never use it,
39 : : * which newer clang reports under -Wunused-template.
40 : : */
41 : : template <typename To, typename From>
42 : 23 : [[maybe_unused]] To cast_by_string(const From& f)
43 : : {
44 : 23 : std::stringstream s;
45 : 23 : s << f;
46 : 46 : return To(s.str());
47 : 23 : }
48 : : } // namespace
49 : :
50 : 2529 : Integer toInteger(const poly::Integer& i)
51 : : {
52 : 2529 : const mpz_class& gi = *poly::detail::cast_to_gmp(&i);
53 : : #ifdef CVC5_GMP_IMP
54 : : return Integer(gi);
55 : : #endif
56 : : #ifdef CVC5_CLN_IMP
57 : 2529 : if (std::numeric_limits<long>::min() <= gi
58 [ + + ][ + + ]: 2529 : && gi <= std::numeric_limits<long>::max())
[ + + ]
59 : : {
60 : 2523 : return Integer(gi.get_si());
61 : : }
62 : : else
63 : : {
64 : 6 : return cast_by_string<Integer, poly::Integer>(i);
65 : : }
66 : : #endif
67 : : }
68 : 4278 : Rational toRational(const poly::Integer& i) { return Rational(toInteger(i)); }
69 : 171 : Rational toRational(const poly::Rational& r)
70 : : {
71 : : #ifdef CVC5_GMP_IMP
72 : : return Rational(*poly::detail::cast_to_gmp(&r));
73 : : #endif
74 : : #ifdef CVC5_CLN_IMP
75 : 342 : return Rational(toInteger(numerator(r)), toInteger(denominator(r)));
76 : : #endif
77 : : }
78 : 24 : Rational toRational(const poly::DyadicRational& dr)
79 : : {
80 : 48 : return Rational(toInteger(numerator(dr)), toInteger(denominator(dr)));
81 : : }
82 : 0 : Rational toRationalAbove(const poly::Value& v)
83 : : {
84 [ - - ]: 0 : if (is_algebraic_number(v))
85 : : {
86 : 0 : return toRational(get_upper_bound(as_algebraic_number(v)));
87 : : }
88 [ - - ]: 0 : else if (is_dyadic_rational(v))
89 : : {
90 : 0 : return toRational(as_dyadic_rational(v));
91 : : }
92 [ - - ]: 0 : else if (is_integer(v))
93 : : {
94 : 0 : return toRational(as_integer(v));
95 : : }
96 [ - - ]: 0 : else if (is_rational(v))
97 : : {
98 : 0 : return toRational(as_rational(v));
99 : : }
100 : 0 : DebugUnhandled() << "Can not convert " << v << " to rational.";
101 : : return Rational();
102 : : }
103 : 0 : Rational toRationalBelow(const poly::Value& v)
104 : : {
105 [ - - ]: 0 : if (is_algebraic_number(v))
106 : : {
107 : 0 : return toRational(get_lower_bound(as_algebraic_number(v)));
108 : : }
109 [ - - ]: 0 : else if (is_dyadic_rational(v))
110 : : {
111 : 0 : return toRational(as_dyadic_rational(v));
112 : : }
113 [ - - ]: 0 : else if (is_integer(v))
114 : : {
115 : 0 : return toRational(as_integer(v));
116 : : }
117 [ - - ]: 0 : else if (is_rational(v))
118 : : {
119 : 0 : return toRational(as_rational(v));
120 : : }
121 : 0 : DebugUnhandled() << "Can not convert " << v << " to rational.";
122 : : return Rational();
123 : : }
124 : :
125 : 20098 : poly::Integer toInteger(const Integer& i)
126 : : {
127 : : #ifdef CVC5_GMP_IMP
128 : : return poly::Integer(i.getValue());
129 : : #endif
130 : : #ifdef CVC5_CLN_IMP
131 [ + + ][ - - ]: 40196 : if (std::numeric_limits<long>::min() <= i.getValue()
132 [ + + ][ + + ]: 40196 : && i.getValue() <= std::numeric_limits<long>::max())
[ + + ][ + - ]
[ - - ]
133 : : {
134 : 20081 : return poly::Integer(cln::cl_I_to_long(i.getValue()));
135 : : }
136 : : else
137 : : {
138 : 34 : return poly::Integer(cast_by_string<mpz_class, Integer>(i));
139 : : }
140 : : #endif
141 : : }
142 : 0 : std::vector<poly::Integer> toInteger(const std::vector<Integer>& vi)
143 : : {
144 : 0 : std::vector<poly::Integer> res;
145 [ - - ]: 0 : for (const auto& i : vi) res.emplace_back(toInteger(i));
146 : 0 : return res;
147 : 0 : }
148 : 770 : poly::Rational toRational(const Rational& r)
149 : : {
150 : : #ifdef CVC5_GMP_IMP
151 : : return poly::Rational(r.getValue());
152 : : #endif
153 : : #ifdef CVC5_CLN_IMP
154 : 1540 : return poly::Rational(toInteger(r.getNumerator()),
155 : 2310 : toInteger(r.getDenominator()));
156 : : #endif
157 : : }
158 : :
159 : 758 : std::optional<poly::DyadicRational> toDyadicRational(const Rational& r)
160 : : {
161 : 758 : Integer den = r.getDenominator();
162 [ + + ]: 758 : if (den.isOne())
163 : : { // It's an integer anyway.
164 : 1466 : return poly::DyadicRational(toInteger(r.getNumerator()));
165 : : }
166 : 25 : unsigned long exp = den.isPow2();
167 [ + + ]: 25 : if (exp > 0)
168 : : {
169 : : // It's a dyadic rational.
170 : 8 : return div_2exp(poly::DyadicRational(toInteger(r.getNumerator())), exp - 1);
171 : : }
172 : 21 : return std::optional<poly::DyadicRational>();
173 : 758 : }
174 : :
175 : 2 : std::optional<poly::DyadicRational> toDyadicRational(const poly::Rational& r)
176 : : {
177 : 2 : poly::Integer den = denominator(r);
178 [ + - ]: 2 : if (den == poly::Integer(1))
179 : : { // It's an integer anyway.
180 : 4 : return poly::DyadicRational(numerator(r));
181 : : }
182 : : // Use bit_size as an estimate for the dyadic exponent.
183 : 0 : unsigned long size = bit_size(den) - 1;
184 [ - - ]: 0 : if (mul_pow2(poly::Integer(1), size) == den)
185 : : {
186 : : // It's a dyadic rational.
187 : 0 : return div_2exp(poly::DyadicRational(numerator(r)), size);
188 : : }
189 : 0 : return std::optional<poly::DyadicRational>();
190 : 2 : }
191 : :
192 : 0 : poly::Rational approximateToDyadic(const poly::Rational& r,
193 : : const poly::Rational& original)
194 : : {
195 : : // Multiply both numerator and denominator by two.
196 : : // Increase or decrease the numerator, depending on whether r is too small or
197 : : // too large.
198 : 0 : poly::Integer n = mul_pow2(numerator(r), 1);
199 [ - - ]: 0 : if (r < original)
200 : : {
201 : 0 : ++n;
202 : : }
203 [ - - ]: 0 : else if (r > original)
204 : : {
205 : 0 : --n;
206 : : }
207 : 0 : return poly::Rational(n, mul_pow2(denominator(r), 1));
208 : 0 : }
209 : :
210 : 2 : poly::AlgebraicNumber toPolyRanWithRefinement(poly::UPolynomial&& p,
211 : : const Rational& lower,
212 : : const Rational& upper)
213 : : {
214 : 2 : std::optional<poly::DyadicRational> ml = toDyadicRational(lower);
215 : 2 : std::optional<poly::DyadicRational> mu = toDyadicRational(upper);
216 [ + + ][ + - ]: 2 : if (ml && mu)
[ + + ]
217 : : {
218 : 1 : return poly::AlgebraicNumber(std::move(p),
219 : 3 : poly::DyadicInterval(ml.value(), mu.value()));
220 : : }
221 : : // The encoded real algebraic number did not have dyadic rational endpoints.
222 : 1 : poly::Rational origl = toRational(lower);
223 : 1 : poly::Rational origu = toRational(upper);
224 : 1 : poly::Rational l(floor(origl));
225 : 1 : poly::Rational u(ceil(origu));
226 : 1 : poly::RationalInterval ri(l, u);
227 [ - + ]: 1 : while (count_real_roots(p, ri) != 1)
228 : : {
229 : 0 : l = approximateToDyadic(l, origl);
230 : 0 : u = approximateToDyadic(u, origu);
231 : 0 : ri = poly::RationalInterval(l, u);
232 : : }
233 [ - + ][ - + ]: 1 : Assert(count_real_roots(p, poly::RationalInterval(l, u)) == 1);
[ - - ]
234 : 1 : ml = toDyadicRational(l);
235 : 1 : mu = toDyadicRational(u);
236 [ + - ][ + - ]: 1 : Assert(ml && mu) << "Both bounds should be dyadic by now.";
[ - + ][ - + ]
[ - - ]
237 : 1 : return poly::AlgebraicNumber(std::move(p),
238 : 3 : poly::DyadicInterval(ml.value(), mu.value()));
239 : 2 : }
240 : :
241 : 1 : RealAlgebraicNumber toRanWithRefinement(poly::UPolynomial&& p,
242 : : const Rational& lower,
243 : : const Rational& upper)
244 : : {
245 : : return RealAlgebraicNumber(
246 : 2 : toPolyRanWithRefinement(std::move(p), lower, upper));
247 : : }
248 : :
249 : 1762922 : std::size_t totalDegree(const poly::Polynomial& p)
250 : : {
251 : 1762922 : std::size_t tdeg = 0;
252 : :
253 : 1762922 : lp_polynomial_traverse_f f =
254 : 2732597 : [](const lp_polynomial_context_t*, lp_monomial_t* m, void* data) {
255 : 2732597 : std::size_t sum = 0;
256 [ + + ]: 6425029 : for (std::size_t i = 0; i < m->n; ++i)
257 : : {
258 : 3692432 : sum += m->p[i].d;
259 : : }
260 : :
261 : 2732597 : std::size_t* td = static_cast<std::size_t*>(data);
262 : 2732597 : *td = std::max(*td, sum);
263 : 2732597 : };
264 : :
265 : 1762922 : lp_polynomial_traverse(p.get_internal(), f, &tdeg);
266 : :
267 : 1762922 : return tdeg;
268 : : }
269 : :
270 : 0 : std::ostream& operator<<(std::ostream& os, const VariableInformation& vi)
271 : : {
272 [ - - ]: 0 : if (vi.var == poly::Variable())
273 : : {
274 : 0 : os << "Totals: ";
275 : 0 : os << "max deg " << vi.max_degree;
276 : 0 : os << ", sum term deg " << vi.sum_term_degree;
277 : 0 : os << ", sum poly deg " << vi.sum_poly_degree;
278 : 0 : os << ", num polys " << vi.num_polynomials;
279 : 0 : os << ", num terms " << vi.num_terms;
280 : : }
281 : : else
282 : : {
283 : 0 : os << "Info for " << stream_variable(*(vi.polyCtx), vi.var) << ": ";
284 : 0 : os << "max deg " << vi.max_degree;
285 : 0 : os << ", max lc deg: " << vi.max_lc_degree;
286 : 0 : os << ", max term tdeg: " << vi.max_terms_tdegree;
287 : 0 : os << ", sum term deg " << vi.sum_term_degree;
288 : 0 : os << ", sum poly deg " << vi.sum_poly_degree;
289 : 0 : os << ", num polys " << vi.num_polynomials;
290 : 0 : os << ", num terms " << vi.num_terms;
291 : : }
292 : 0 : return os;
293 : : }
294 : :
295 : : struct GetVarInfo
296 : : {
297 : : VariableInformation* info;
298 : : std::size_t cur_var_degree = 0;
299 : : std::size_t cur_lc_degree = 0;
300 : : };
301 : 32591 : void getVariableInformation(VariableInformation& vi,
302 : : const poly::Polynomial& poly)
303 : : {
304 : 32591 : GetVarInfo varinfo;
305 : 32591 : varinfo.info = &vi;
306 : 32591 : lp_polynomial_traverse_f f =
307 : 53370 : [](const lp_polynomial_context_t*, lp_monomial_t* m, void* data) {
308 : 53370 : GetVarInfo* gvi = static_cast<GetVarInfo*>(data);
309 : 53370 : VariableInformation* info = gvi->info;
310 : : // Total degree of this term
311 : 53370 : std::size_t tdeg = 0;
312 : : // Degree of this variable within this term
313 : 53370 : std::size_t vardeg = 0;
314 [ + + ]: 118495 : for (std::size_t i = 0; i < m->n; ++i)
315 : : {
316 : 65125 : tdeg += m->p[i].d;
317 [ + + ]: 65125 : if (poly::Variable(m->p[i].x) == info->var)
318 : : {
319 : 11831 : info->max_degree = std::max(info->max_degree, m->p[i].d);
320 : 11831 : info->sum_term_degree += m->p[i].d;
321 : 11831 : vardeg = m->p[i].d;
322 : : }
323 : : }
324 [ - + ]: 53370 : if (info->var == poly::Variable())
325 : : {
326 : 0 : ++info->num_terms;
327 : 0 : info->max_degree = std::max(info->max_degree, tdeg);
328 : 0 : info->sum_term_degree += tdeg;
329 : : }
330 [ + + ]: 53370 : else if (vardeg > 0)
331 : : {
332 : 11831 : ++info->num_terms;
333 [ + - ]: 11831 : if (gvi->cur_var_degree < vardeg)
334 : : {
335 : 11831 : gvi->cur_lc_degree = tdeg - vardeg;
336 : : }
337 : 11831 : info->max_terms_tdegree = std::max(info->max_terms_tdegree, tdeg);
338 : : }
339 : 53370 : };
340 : 32591 : std::size_t tmp_max_degree = vi.max_degree;
341 : 32591 : std::size_t tmp_num_terms = vi.num_terms;
342 : 32591 : vi.max_degree = 0;
343 : 32591 : vi.num_terms = 0;
344 : 32591 : lp_polynomial_traverse(poly.get_internal(), f, &varinfo);
345 : 32591 : vi.max_lc_degree = std::max(vi.max_lc_degree, varinfo.cur_lc_degree);
346 [ + + ]: 32591 : if (vi.num_terms > 0)
347 : : {
348 : 10727 : ++vi.num_polynomials;
349 : : }
350 : 32591 : vi.sum_poly_degree += vi.max_degree;
351 : 32591 : vi.max_degree = std::max(vi.max_degree, tmp_max_degree);
352 : 32591 : vi.num_terms += tmp_num_terms;
353 : 32591 : }
354 : :
355 : : } // namespace poly_utils
356 : : } // namespace cvc5::internal
357 : :
358 : : #endif
|
