src/util/poly_util.cpp

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