Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	src/util/poly_util.cpp	Lines:	102	163	62.6 %
Date:	2021-09-18	Branches:	115	360	31.9 %


/******************************************************************************
 * Top contributors (to current version):
 *   Gereon Kremer
 *
 * This file is part of the cvc5 project.
 *
 * Copyright (c) 2009-2021 by the authors listed in the file AUTHORS
 * in the top-level source directory and their institutional affiliations.
 * All rights reserved.  See the file COPYING in the top-level source
 * directory for licensing information.
 * ****************************************************************************
 *
 * Utilities for working with LibPoly.
 */

#include "poly_util.h"

#ifdef CVC5_POLY_IMP

#include <poly/polyxx.h>

#include <map>
#include <sstream>

#include "base/check.h"
#include "util/integer.h"
#include "util/rational.h"
#include "util/real_algebraic_number.h"

namespace cvc5 {
namespace poly_utils {

namespace {
/**
 * Convert arbitrary data using a string as intermediary.
 * Assumes the existence of operator<<(std::ostream&, const From&) and To(const
 * std::string&); should be the last resort for type conversions: it may not
 * only yield bad performance, but is also dependent on compatible string
 * representations. Use with care!
 */
template <typename To, typename From>
To cast_by_string(const From& f)
{
  std::stringstream s;
  s << f;
  return To(s.str());
}
}  // namespace

Integer toInteger(const poly::Integer& i)
{
  const mpz_class& gi = *poly::detail::cast_to_gmp(&i);
#ifdef CVC5_GMP_IMP
  return Integer(gi);
#endif
#ifdef CVC5_CLN_IMP
  if (std::numeric_limits<long>::min() <= gi
      && gi <= std::numeric_limits<long>::max())
  {
    return Integer(gi.get_si());
  }
  else
  {
    return cast_by_string<Integer, poly::Integer>(i);
  }
#endif
}
Rational toRational(const poly::Integer& i) { return Rational(toInteger(i)); }
Rational toRational(const poly::Rational& r)
{
#ifdef CVC5_GMP_IMP
  return Rational(*poly::detail::cast_to_gmp(&r));
#endif
#ifdef CVC5_CLN_IMP
  return Rational(toInteger(numerator(r)), toInteger(denominator(r)));
#endif
}
Rational toRational(const poly::DyadicRational& dr)
{
  return Rational(toInteger(numerator(dr)), toInteger(denominator(dr)));
}
Rational toRationalAbove(const poly::Value& v)
{
  if (is_algebraic_number(v))
  {
    return toRational(get_upper_bound(as_algebraic_number(v)));
  }
  else if (is_dyadic_rational(v))
  {
    return toRational(as_dyadic_rational(v));
  }
  else if (is_integer(v))
  {
    return toRational(as_integer(v));
  }
  else if (is_rational(v))
  {
    return toRational(as_rational(v));
  }
  Assert(false) << "Can not convert " << v << " to rational.";
  return Rational();
}
Rational toRationalBelow(const poly::Value& v)
{
  if (is_algebraic_number(v))
  {
    return toRational(get_lower_bound(as_algebraic_number(v)));
  }
  else if (is_dyadic_rational(v))
  {
    return toRational(as_dyadic_rational(v));
  }
  else if (is_integer(v))
  {
    return toRational(as_integer(v));
  }
  else if (is_rational(v))
  {
    return toRational(as_rational(v));
  }
  Assert(false) << "Can not convert " << v << " to rational.";
  return Rational();
}

poly::Integer toInteger(const Integer& i)
{
#ifdef CVC5_GMP_IMP
  return poly::Integer(i.getValue());
#endif
#ifdef CVC5_CLN_IMP
  if (std::numeric_limits<long>::min() <= i.getValue()
      && i.getValue() <= std::numeric_limits<long>::max())
  {
    return poly::Integer(cln::cl_I_to_long(i.getValue()));
  }
  else
  {
    return poly::Integer(cast_by_string<mpz_class, Integer>(i));
  }
#endif
}
std::vector<poly::Integer> toInteger(const std::vector<Integer>& vi)
{
  std::vector<poly::Integer> res;
  for (const auto& i : vi) res.emplace_back(toInteger(i));
  return res;
}
poly::Rational toRational(const Rational& r)
{
#ifdef CVC5_GMP_IMP
  return poly::Rational(r.getValue());
#endif
#ifdef CVC5_CLN_IMP
  return poly::Rational(toInteger(r.getNumerator()),
                        toInteger(r.getDenominator()));
#endif
}

std::optional<poly::DyadicRational> toDyadicRational(const Rational& r)
{
  Integer den = r.getDenominator();
  if (den.isOne())
  {  // It's an integer anyway.
    return poly::DyadicRational(toInteger(r.getNumerator()));
  }
  unsigned long exp = den.isPow2();
  if (exp > 0)
  {
    // It's a dyadic rational.
    return div_2exp(poly::DyadicRational(toInteger(r.getNumerator())), exp - 1);
  }
  return std::optional<poly::DyadicRational>();
}

std::optional<poly::DyadicRational> toDyadicRational(const poly::Rational& r)
{
  poly::Integer den = denominator(r);
  if (den == poly::Integer(1))
  {  // It's an integer anyway.
    return poly::DyadicRational(numerator(r));
  }
  // Use bit_size as an estimate for the dyadic exponent.
  unsigned long size = bit_size(den) - 1;
  if (mul_pow2(poly::Integer(1), size) == den)
  {
    // It's a dyadic rational.
    return div_2exp(poly::DyadicRational(numerator(r)), size);
  }
  return std::optional<poly::DyadicRational>();
}

poly::Rational approximateToDyadic(const poly::Rational& r,
                                   const poly::Rational& original)
{
  // Multiply both numerator and denominator by two.
  // Increase or decrease the numerator, depending on whether r is too small or
  // too large.
  poly::Integer n = mul_pow2(numerator(r), 1);
  if (r < original)
  {
    ++n;
  }
  else if (r > original)
  {
    --n;
  }
  return poly::Rational(n, mul_pow2(denominator(r), 1));
}

poly::AlgebraicNumber toPolyRanWithRefinement(poly::UPolynomial&& p,
                                              const Rational& lower,
                                              const Rational& upper)
{
  std::optional<poly::DyadicRational> ml = toDyadicRational(lower);
  std::optional<poly::DyadicRational> mu = toDyadicRational(upper);
  if (ml && mu)
  {
    return poly::AlgebraicNumber(std::move(p),
                                 poly::DyadicInterval(ml.value(), mu.value()));
  }
  // The encoded real algebraic number did not have dyadic rational endpoints.
  poly::Rational origl = toRational(lower);
  poly::Rational origu = toRational(upper);
  poly::Rational l(floor(origl));
  poly::Rational u(ceil(origu));
  poly::RationalInterval ri(l, u);
  while (count_real_roots(p, ri) != 1)
  {
    l = approximateToDyadic(l, origl);
    u = approximateToDyadic(u, origu);
    ri = poly::RationalInterval(l, u);
  }
  Assert(count_real_roots(p, poly::RationalInterval(l, u)) == 1);
  ml = toDyadicRational(l);
  mu = toDyadicRational(u);
  Assert(ml && mu) << "Both bounds should be dyadic by now.";
  return poly::AlgebraicNumber(std::move(p),
                               poly::DyadicInterval(ml.value(), mu.value()));
}

RealAlgebraicNumber toRanWithRefinement(poly::UPolynomial&& p,
                                        const Rational& lower,
                                        const Rational& upper)
{
  return RealAlgebraicNumber(
      toPolyRanWithRefinement(std::move(p), lower, upper));
}

std::size_t totalDegree(const poly::Polynomial& p)
{
  std::size_t tdeg = 0;

  lp_polynomial_traverse_f f =
      [](const lp_polynomial_context_t* ctx, lp_monomial_t* m, void* data) {
        std::size_t sum = 0;
        for (std::size_t i = 0; i < m->n; ++i)
        {
          sum += m->p[i].d;
        }

        std::size_t* td = static_cast<std::size_t*>(data);
        *td = std::max(*td, sum);
      };

  lp_polynomial_traverse(p.get_internal(), f, &tdeg);

  return tdeg;
}

std::ostream& operator<<(std::ostream& os, const VariableInformation& vi)
{
  if (vi.var == poly::Variable())
  {
    os << "Totals: ";
    os << "max deg " << vi.max_degree;
    os << ", sum term deg " << vi.sum_term_degree;
    os << ", sum poly deg " << vi.sum_poly_degree;
    os << ", num polys " << vi.num_polynomials;
    os << ", num terms " << vi.num_terms;
  }
  else
  {
    os << "Info for " << vi.var << ": ";
    os << "max deg " << vi.max_degree;
    os << ", max lc deg: " << vi.max_lc_degree;
    os << ", max term tdeg: " << vi.max_terms_tdegree;
    os << ", sum term deg " << vi.sum_term_degree;
    os << ", sum poly deg " << vi.sum_poly_degree;
    os << ", num polys " << vi.num_polynomials;
    os << ", num terms " << vi.num_terms;
  }
  return os;
}

struct GetVarInfo
{
  VariableInformation* info;
  std::size_t cur_var_degree = 0;
  std::size_t cur_lc_degree = 0;
};
void getVariableInformation(VariableInformation& vi,
                            const poly::Polynomial& poly)
{
  GetVarInfo varinfo;
  varinfo.info = &vi;
  lp_polynomial_traverse_f f =
      [](const lp_polynomial_context_t* ctx, lp_monomial_t* m, void* data) {
        GetVarInfo* gvi = static_cast<GetVarInfo*>(data);
        VariableInformation* info = gvi->info;
        // Total degree of this term
        std::size_t tdeg = 0;
        // Degree of this variable within this term
        std::size_t vardeg = 0;
        for (std::size_t i = 0; i < m->n; ++i)
        {
          tdeg += m->p[i].d;
          if (m->p[i].x == info->var)
          {
            info->max_degree = std::max(info->max_degree, m->p[i].d);
            info->sum_term_degree += m->p[i].d;
            vardeg = m->p[i].d;
          }
        }
        if (info->var == poly::Variable())
        {
          ++info->num_terms;
          info->max_degree = std::max(info->max_degree, tdeg);
          info->sum_term_degree += tdeg;
        }
        else if (vardeg > 0)
        {
          ++info->num_terms;
          if (gvi->cur_var_degree < vardeg)
          {
            gvi->cur_lc_degree = tdeg - vardeg;
          }
          info->max_terms_tdegree = std::max(info->max_terms_tdegree, tdeg);
        }
      };
  std::size_t tmp_max_degree = vi.max_degree;
  std::size_t tmp_num_terms = vi.num_terms;
  vi.max_degree = 0;
  vi.num_terms = 0;
  lp_polynomial_traverse(poly.get_internal(), f, &varinfo);
  vi.max_lc_degree = std::max(vi.max_lc_degree, varinfo.cur_lc_degree);
  if (vi.num_terms > 0)
  {
    ++vi.num_polynomials;
  }
  vi.sum_poly_degree += vi.max_degree;
  vi.max_degree = std::max(vi.max_degree, tmp_max_degree);
  vi.num_terms += tmp_num_terms;
}

}  // namespace poly_utils
}  // namespace cvc5

#endif


Generated by: GCOVR (Version 3.2)

Line	Exec	Source
1		/******************************************************************************
2		* Top contributors (to current version):
3		* Gereon Kremer
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 working with LibPoly.
14		*/
15
16		#include "poly_util.h"
17
18		#ifdef CVC5_POLY_IMP
19
20		#include <poly/polyxx.h>
21
22		#include <map>
23		#include <sstream>
24
25		#include "base/check.h"
26		#include "util/integer.h"
27		#include "util/rational.h"
28		#include "util/real_algebraic_number.h"
29
30		namespace cvc5 {
31		namespace poly_utils {
32
33		namespace {
34		/**
35		* Convert arbitrary data using a string as intermediary.
36		* Assumes the existence of operator<<(std::ostream&, const From&) and To(const
37		* std::string&); should be the last resort for type conversions: it may not
38		* only yield bad performance, but is also dependent on compatible string
39		* representations. Use with care!
40		*/
41		template <typename To, typename From>
42	332	To cast_by_string(const From& f)
43		{
44	664	std::stringstream s;
45	332	s << f;
46	664	return To(s.str());
47		}
48		} // namespace
49
50	144	Integer toInteger(const poly::Integer& i)
51		{
52	144	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	288	if (std::numeric_limits<long>::min() <= gi
58	144	&& gi <= std::numeric_limits<long>::max())
59		{
60	144	return Integer(gi.get_si());
61		}
62		else
63		{
64		return cast_by_string<Integer, poly::Integer>(i);
65		}
66		#endif
67		}
68	30	Rational toRational(const poly::Integer& i) { return Rational(toInteger(i)); }
69	36	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	36	return Rational(toInteger(numerator(r)), toInteger(denominator(r)));
76		#endif
77		}
78	21	Rational toRational(const poly::DyadicRational& dr)
79		{
80	21	return Rational(toInteger(numerator(dr)), toInteger(denominator(dr)));
81		}
82		Rational toRationalAbove(const poly::Value& v)
83		{
84		if (is_algebraic_number(v))
85		{
86		return toRational(get_upper_bound(as_algebraic_number(v)));
87		}
88		else if (is_dyadic_rational(v))
89		{
90		return toRational(as_dyadic_rational(v));
91		}
92		else if (is_integer(v))
93		{
94		return toRational(as_integer(v));
95		}
96		else if (is_rational(v))
97		{
98		return toRational(as_rational(v));
99		}
100		Assert(false) << "Can not convert " << v << " to rational.";
101		return Rational();
102		}
103		Rational toRationalBelow(const poly::Value& v)
104		{
105		if (is_algebraic_number(v))
106		{
107		return toRational(get_lower_bound(as_algebraic_number(v)));
108		}
109		else if (is_dyadic_rational(v))
110		{
111		return toRational(as_dyadic_rational(v));
112		}
113		else if (is_integer(v))
114		{
115		return toRational(as_integer(v));
116		}
117		else if (is_rational(v))
118		{
119		return toRational(as_rational(v));
120		}
121		Assert(false) << "Can not convert " << v << " to rational.";
122		return Rational();
123		}
124
125	2804	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	8412	if (std::numeric_limits<long>::min() <= i.getValue()
132	8412	&& i.getValue() <= std::numeric_limits<long>::max())
133		{
134	2472	return poly::Integer(cln::cl_I_to_long(i.getValue()));
135		}
136		else
137		{
138	332	return poly::Integer(cast_by_string<mpz_class, Integer>(i));
139		}
140		#endif
141		}
142		std::vector<poly::Integer> toInteger(const std::vector<Integer>& vi)
143		{
144		std::vector<poly::Integer> res;
145		for (const auto& i : vi) res.emplace_back(toInteger(i));
146		return res;
147		}
148	158	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	316	return poly::Rational(toInteger(r.getNumerator()),
155	474	toInteger(r.getDenominator()));
156		#endif
157		}
158
159	10	std::optional<poly::DyadicRational> toDyadicRational(const Rational& r)
160		{
161	20	Integer den = r.getDenominator();
162	10	if (den.isOne())
163		{ // It's an integer anyway.
164	2	return poly::DyadicRational(toInteger(r.getNumerator()));
165		}
166	8	unsigned long exp = den.isPow2();
167	8	if (exp > 0)
168		{
169		// It's a dyadic rational.
170	4	return div_2exp(poly::DyadicRational(toInteger(r.getNumerator())), exp - 1);
171		}
172	4	return std::optional<poly::DyadicRational>();
173		}
174
175	4	std::optional<poly::DyadicRational> toDyadicRational(const poly::Rational& r)
176		{
177	8	poly::Integer den = denominator(r);
178	4	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		unsigned long size = bit_size(den) - 1;
184		if (mul_pow2(poly::Integer(1), size) == den)
185		{
186		// It's a dyadic rational.
187		return div_2exp(poly::DyadicRational(numerator(r)), size);
188		}
189		return std::optional<poly::DyadicRational>();
190		}
191
192		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		poly::Integer n = mul_pow2(numerator(r), 1);
199		if (r < original)
200		{
201		++n;
202		}
203		else if (r > original)
204		{
205		--n;
206		}
207		return poly::Rational(n, mul_pow2(denominator(r), 1));
208		}
209
210	4	poly::AlgebraicNumber toPolyRanWithRefinement(poly::UPolynomial&& p,
211		const Rational& lower,
212		const Rational& upper)
213		{
214	8	std::optional<poly::DyadicRational> ml = toDyadicRational(lower);
215	8	std::optional<poly::DyadicRational> mu = toDyadicRational(upper);
216	4	if (ml && mu)
217		{
218	2	return poly::AlgebraicNumber(std::move(p),
219	4	poly::DyadicInterval(ml.value(), mu.value()));
220		}
221		// The encoded real algebraic number did not have dyadic rational endpoints.
222	4	poly::Rational origl = toRational(lower);
223	4	poly::Rational origu = toRational(upper);
224	4	poly::Rational l(floor(origl));
225	4	poly::Rational u(ceil(origu));
226	4	poly::RationalInterval ri(l, u);
227	2	while (count_real_roots(p, ri) != 1)
228		{
229		l = approximateToDyadic(l, origl);
230		u = approximateToDyadic(u, origu);
231		ri = poly::RationalInterval(l, u);
232		}
233	2	Assert(count_real_roots(p, poly::RationalInterval(l, u)) == 1);
234	2	ml = toDyadicRational(l);
235	2	mu = toDyadicRational(u);
236	2	Assert(ml && mu) << "Both bounds should be dyadic by now.";
237	2	return poly::AlgebraicNumber(std::move(p),
238	4	poly::DyadicInterval(ml.value(), mu.value()));
239		}
240
241	2	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	98150	std::size_t totalDegree(const poly::Polynomial& p)
250		{
251	98150	std::size_t tdeg = 0;
252
253	98150	lp_polynomial_traverse_f f =
254	392560	[](const lp_polynomial_context_t* ctx, lp_monomial_t* m, void* data) {
255	196280	std::size_t sum = 0;
256	317116	for (std::size_t i = 0; i < m->n; ++i)
257		{
258	120836	sum += m->p[i].d;
259		}
260
261	196280	std::size_t* td = static_cast<std::size_t*>(data);
262	196280	td = std::max(td, sum);
263	392560	};
264
265	98150	lp_polynomial_traverse(p.get_internal(), f, &tdeg);
266
267	98150	return tdeg;
268		}
269
270		std::ostream& operator<<(std::ostream& os, const VariableInformation& vi)
271		{
272		if (vi.var == poly::Variable())
273		{
274		os << "Totals: ";
275		os << "max deg " << vi.max_degree;
276		os << ", sum term deg " << vi.sum_term_degree;
277		os << ", sum poly deg " << vi.sum_poly_degree;
278		os << ", num polys " << vi.num_polynomials;
279		os << ", num terms " << vi.num_terms;
280		}
281		else
282		{
283		os << "Info for " << vi.var << ": ";
284		os << "max deg " << vi.max_degree;
285		os << ", max lc deg: " << vi.max_lc_degree;
286		os << ", max term tdeg: " << vi.max_terms_tdegree;
287		os << ", sum term deg " << vi.sum_term_degree;
288		os << ", sum poly deg " << vi.sum_poly_degree;
289		os << ", num polys " << vi.num_polynomials;
290		os << ", num terms " << vi.num_terms;
291		}
292		return os;
293		}
294
295	2248	struct GetVarInfo
296		{
297		VariableInformation* info;
298		std::size_t cur_var_degree = 0;
299		std::size_t cur_lc_degree = 0;
300		};
301	2248	void getVariableInformation(VariableInformation& vi,
302		const poly::Polynomial& poly)
303		{
304	2248	GetVarInfo varinfo;
305	2248	varinfo.info = &vi;
306	2248	lp_polynomial_traverse_f f =
307	9626	[](const lp_polynomial_context_t* ctx, lp_monomial_t* m, void* data) {
308	4813	GetVarInfo* gvi = static_cast<GetVarInfo*>(data);
309	4813	VariableInformation* info = gvi->info;
310		// Total degree of this term
311	4813	std::size_t tdeg = 0;
312		// Degree of this variable within this term
313	4813	std::size_t vardeg = 0;
314	8577	for (std::size_t i = 0; i < m->n; ++i)
315		{
316	3764	tdeg += m->p[i].d;
317	3764	if (m->p[i].x == info->var)
318		{
319	1597	info->max_degree = std::max(info->max_degree, m->p[i].d);
320	1597	info->sum_term_degree += m->p[i].d;
321	1597	vardeg = m->p[i].d;
322		}
323		}
324	4813	if (info->var == poly::Variable())
325		{
326		++info->num_terms;
327		info->max_degree = std::max(info->max_degree, tdeg);
328		info->sum_term_degree += tdeg;
329		}
330	4813	else if (vardeg > 0)
331		{
332	1597	++info->num_terms;
333	1597	if (gvi->cur_var_degree < vardeg)
334		{
335	1597	gvi->cur_lc_degree = tdeg - vardeg;
336		}
337	1597	info->max_terms_tdegree = std::max(info->max_terms_tdegree, tdeg);
338		}
339	9626	};
340	2248	std::size_t tmp_max_degree = vi.max_degree;
341	2248	std::size_t tmp_num_terms = vi.num_terms;
342	2248	vi.max_degree = 0;
343	2248	vi.num_terms = 0;
344	2248	lp_polynomial_traverse(poly.get_internal(), f, &varinfo);
345	2248	vi.max_lc_degree = std::max(vi.max_lc_degree, varinfo.cur_lc_degree);
346	2248	if (vi.num_terms > 0)
347		{
348	1405	++vi.num_polynomials;
349		}
350	2248	vi.sum_poly_degree += vi.max_degree;
351	2248	vi.max_degree = std::max(vi.max_degree, tmp_max_degree);
352	2248	vi.num_terms += tmp_num_terms;
353	2248	}
354
355		} // namespace poly_utils
356	29574	} // namespace cvc5
357
358		#endif