Branch data     Line data    Source code

       1                 :            : /******************************************************************************
       2                 :            :  * Top contributors (to current version):
       3                 :            :  *   Andrew Reynolds, Gereon Kremer, Daniel Larraz
       4                 :            :  *
       5                 :            :  * This file is part of the cvc5 project.
       6                 :            :  *
       7                 :            :  * Copyright (c) 2009-2025 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                 :            :  * Implementation of RealAlgebraicNumber based on libpoly.
      14                 :            :  */
      15                 :            : 
      16                 :            : #include "base/cvc5config.h"
      17                 :            : #include "util/real_algebraic_number.h"
      18                 :            : 
      19                 :            : #ifdef CVC5_POLY_IMP
      20                 :            : #include <poly/polyxx.h>
      21                 :            : #endif
      22                 :            : 
      23                 :            : #include <limits>
      24                 :            : #include <sstream>
      25                 :            : 
      26                 :            : #include "base/check.h"
      27                 :            : #include "util/poly_util.h"
      28                 :            : 
      29                 :            : #define RAN_UNREACHABLE \
      30                 :            :   Unreachable() << "RealAlgebraicNumber is not available without libpoly."
      31                 :            : 
      32                 :            : namespace cvc5::internal {
      33                 :            : 
      34                 :    2413120 : RealAlgebraicNumber::RealAlgebraicNumber()
      35                 :            :     :
      36                 :            : #ifdef CVC5_POLY_IMP
      37                 :            :       d_isRational(true),
      38                 :            : #endif
      39                 :    2413120 :       d_rat()
      40                 :            : {
      41                 :    2413120 : }
      42                 :            : 
      43                 :            : #ifdef CVC5_POLY_IMP
      44                 :        168 : RealAlgebraicNumber::RealAlgebraicNumber(poly::AlgebraicNumber&& an)
      45                 :        168 :     : d_isRational(false), d_value(std::move(an))
      46                 :            : {
      47                 :        168 : }
      48                 :            : #endif
      49                 :            : 
      50                 :   23408800 : RealAlgebraicNumber::RealAlgebraicNumber(const Integer& i)
      51                 :            :     :
      52                 :            : #ifdef CVC5_POLY_IMP
      53                 :            :       d_isRational(true),
      54                 :            : #endif
      55                 :   23408800 :       d_rat(i)
      56                 :            : {
      57                 :   23408800 : }
      58                 :            : 
      59                 :   28240300 : RealAlgebraicNumber::RealAlgebraicNumber(const Rational& r)
      60                 :            :     :
      61                 :            : #ifdef CVC5_POLY_IMP
      62                 :            :       d_isRational(true),
      63                 :            : #endif
      64                 :   28240300 :       d_rat(r)
      65                 :            : {
      66                 :   28240300 : }
      67                 :            : 
      68                 :         21 : RealAlgebraicNumber::RealAlgebraicNumber(const std::vector<long>& coefficients,
      69                 :            :                                          long lower,
      70                 :         21 :                                          long upper)
      71                 :            : {
      72                 :            : #ifdef CVC5_ASSERTIONS
      73         [ +  + ]:         84 :   for (long c : coefficients)
      74                 :            :   {
      75 [ +  - ][ +  - ]:        126 :     Assert(std::numeric_limits<int32_t>::min() <= c
         [ -  + ][ -  - ]
      76                 :            :            && c <= std::numeric_limits<int32_t>::max())
      77                 :            :         << "Coefficients need to fit within 32 bit integers. Please use the "
      78                 :          0 :            "constructor based on Integer instead.";
      79                 :            :   }
      80                 :            : #endif
      81                 :            : #ifdef CVC5_POLY_IMP
      82                 :         21 :   d_isRational = false;
      83                 :         42 :   d_value = poly::AlgebraicNumber(poly::UPolynomial(coefficients),
      84                 :         63 :                                   poly::DyadicInterval(lower, upper));
      85                 :            : #else
      86                 :            :   RAN_UNREACHABLE;
      87                 :            : #endif
      88                 :         21 : }
      89                 :            : 
      90                 :          0 : RealAlgebraicNumber::RealAlgebraicNumber(
      91                 :            :     const std::vector<Integer>& coefficients,
      92                 :            :     const Rational& lower,
      93                 :          0 :     const Rational& upper)
      94                 :            : {
      95                 :            : #ifdef CVC5_POLY_IMP
      96                 :          0 :   d_isRational = false;
      97                 :          0 :   *this = poly_utils::toRanWithRefinement(
      98                 :          0 :       poly::UPolynomial(poly_utils::toInteger(coefficients)), lower, upper);
      99                 :            : #else
     100                 :            :   RAN_UNREACHABLE;
     101                 :            : #endif
     102                 :          0 : }
     103                 :          1 : RealAlgebraicNumber::RealAlgebraicNumber(
     104                 :            :     const std::vector<Rational>& coefficients,
     105                 :            :     const Rational& lower,
     106                 :          1 :     const Rational& upper)
     107                 :            : {
     108                 :            : #ifdef CVC5_POLY_IMP
     109                 :          1 :   d_isRational = false;
     110                 :          2 :   Integer factor = Integer(1);
     111         [ +  + ]:          4 :   for (const auto& c : coefficients)
     112                 :            :   {
     113                 :          3 :     factor = factor.lcm(c.getDenominator());
     114                 :            :   }
     115                 :          1 :   std::vector<poly::Integer> coeffs;
     116         [ +  + ]:          4 :   for (const auto& c : coefficients)
     117                 :            :   {
     118 [ -  + ][ -  + ]:          3 :     Assert((c * factor).getDenominator() == Integer(1));
                 [ -  - ]
     119                 :          3 :     coeffs.emplace_back(poly_utils::toInteger((c * factor).getNumerator()));
     120                 :            :   }
     121                 :          2 :   *this = poly_utils::toRanWithRefinement(
     122                 :          3 :       poly::UPolynomial(std::move(coeffs)), lower, upper);
     123                 :            : #else
     124                 :            :   RAN_UNREACHABLE;
     125                 :            : #endif
     126                 :          1 : }
     127                 :            : 
     128                 :   36745200 : bool RealAlgebraicNumber::isRational() const
     129                 :            : {
     130                 :            : #ifdef CVC5_POLY_IMP
     131         [ +  + ]:   36745200 :   if (!d_isRational)
     132                 :            :   {
     133                 :       1404 :     return poly::is_rational(getValue());
     134                 :            :   }
     135                 :            : #endif
     136                 :   36743800 :   return true;
     137                 :            : }
     138                 :   35027900 : Rational RealAlgebraicNumber::toRational() const
     139                 :            : {
     140                 :            : #ifdef CVC5_POLY_IMP
     141         [ +  + ]:   35027900 :   if (!d_isRational)
     142                 :            :   {
     143                 :         88 :     return poly_utils::toRational(poly::to_rational_approximation(getValue()));
     144                 :            :   }
     145                 :            : #endif
     146                 :   35027900 :   return getRationalValue();
     147                 :            : }
     148                 :            : 
     149                 :          4 : std::string RealAlgebraicNumber::toString() const
     150                 :            : {
     151                 :          8 :   std::stringstream ss;
     152                 :            : #ifdef CVC5_POLY_IMP
     153         [ +  - ]:          4 :   if (!d_isRational)
     154                 :            :   {
     155                 :          4 :     ss << getValue();
     156                 :          4 :     return ss.str();
     157                 :            :   }
     158                 :            : #endif
     159                 :          0 :   ss << getRationalValue();
     160                 :          0 :   return ss.str();
     161                 :            : }
     162                 :            : 
     163                 :            : #ifdef CVC5_POLY_IMP
     164                 :       4670 : poly::AlgebraicNumber RealAlgebraicNumber::convertToPoly(
     165                 :            :     const RealAlgebraicNumber& r)
     166                 :            : {
     167                 :            :   // if we are already poly, just return the value
     168         [ +  + ]:       4670 :   if (!r.d_isRational)
     169                 :            :   {
     170                 :       4031 :     return r.d_value;
     171                 :            :   }
     172                 :            :   // otherwise, this converts the rational value of r to poly
     173                 :        639 :   const Rational& rr = r.getRationalValue();
     174                 :       1278 :   poly::Rational pr = poly_utils::toRational(rr);
     175                 :       1278 :   auto dr = poly_utils::toDyadicRational(rr);
     176         [ +  + ]:        639 :   if (dr)
     177                 :            :   {
     178                 :        611 :     return poly::AlgebraicNumber(dr.value());
     179                 :            :   }
     180                 :            :   return poly::AlgebraicNumber(
     181 [ +  + ][ -  - ]:        112 :       poly::UPolynomial({-numerator(pr), denominator(pr)}),
     182                 :         84 :       poly::DyadicInterval(floor(pr), ceil(pr)));
     183                 :            : }
     184                 :            : #endif
     185                 :            : 
     186                 :       1707 : bool RealAlgebraicNumber::operator==(const RealAlgebraicNumber& rhs) const
     187                 :            : {
     188                 :            : #ifdef CVC5_POLY_IMP
     189 [ +  + ][ +  - ]:       1707 :   if (!d_isRational || !rhs.d_isRational)
     190                 :            :   {
     191                 :       1707 :     return convertToPoly(*this) == convertToPoly(rhs);
     192                 :            :   }
     193                 :            : #endif
     194                 :          0 :   return getRationalValue() == rhs.getRationalValue();
     195                 :            : }
     196                 :          0 : bool RealAlgebraicNumber::operator!=(const RealAlgebraicNumber& rhs) const
     197                 :            : {
     198                 :          0 :   return !(*this == rhs);
     199                 :            : }
     200                 :         12 : bool RealAlgebraicNumber::operator<(const RealAlgebraicNumber& rhs) const
     201                 :            : {
     202                 :            : #ifdef CVC5_POLY_IMP
     203 [ +  + ][ +  - ]:         12 :   if (!d_isRational || !rhs.d_isRational)
     204                 :            :   {
     205                 :         12 :     return convertToPoly(*this) < convertToPoly(rhs);
     206                 :            :   }
     207                 :            : #endif
     208                 :          0 :   return getRationalValue() < rhs.getRationalValue();
     209                 :            : }
     210                 :        152 : bool RealAlgebraicNumber::operator<=(const RealAlgebraicNumber& rhs) const
     211                 :            : {
     212                 :            : #ifdef CVC5_POLY_IMP
     213 [ +  - ][ +  - ]:        152 :   if (!d_isRational || !rhs.d_isRational)
     214                 :            :   {
     215                 :        152 :     return convertToPoly(*this) <= convertToPoly(rhs);
     216                 :            :   }
     217                 :            : #endif
     218                 :          0 :   return getRationalValue() <= rhs.getRationalValue();
     219                 :            : }
     220                 :          0 : bool RealAlgebraicNumber::operator>(const RealAlgebraicNumber& rhs) const
     221                 :            : {
     222                 :          0 :   return rhs < *this;
     223                 :            : }
     224                 :        152 : bool RealAlgebraicNumber::operator>=(const RealAlgebraicNumber& rhs) const
     225                 :            : {
     226                 :        152 :   return rhs <= *this;
     227                 :            :   ;
     228                 :            : }
     229                 :            : 
     230                 :          3 : RealAlgebraicNumber RealAlgebraicNumber::operator+(
     231                 :            :     const RealAlgebraicNumber& rhs) const
     232                 :            : {
     233                 :            : #ifdef CVC5_POLY_IMP
     234                 :            :   // if either is poly, we convert both and return the result
     235 [ -  + ][ -  - ]:          3 :   if (!d_isRational || !rhs.d_isRational)
     236                 :            :   {
     237                 :          6 :     return convertToPoly(*this) + convertToPoly(rhs);
     238                 :            :   }
     239                 :            : #endif
     240                 :          0 :   return getRationalValue() + rhs.getRationalValue();
     241                 :            : }
     242                 :          0 : RealAlgebraicNumber RealAlgebraicNumber::operator-(
     243                 :            :     const RealAlgebraicNumber& rhs) const
     244                 :            : {
     245                 :            : #ifdef CVC5_POLY_IMP
     246 [ -  - ][ -  - ]:          0 :   if (!d_isRational || !rhs.d_isRational)
     247                 :            :   {
     248                 :          0 :     return convertToPoly(*this) - convertToPoly(rhs);
     249                 :            :   }
     250                 :            : #endif
     251                 :          0 :   return getRationalValue() - rhs.getRationalValue();
     252                 :            : }
     253                 :    3515440 : RealAlgebraicNumber RealAlgebraicNumber::operator-() const
     254                 :            : {
     255                 :            : #ifdef CVC5_POLY_IMP
     256         [ +  + ]:    3515440 :   if (!d_isRational)
     257                 :            :   {
     258                 :         20 :     return -getValue();
     259                 :            :   }
     260                 :            : #endif
     261                 :    7030860 :   return -getRationalValue();
     262                 :            : }
     263                 :     145211 : RealAlgebraicNumber RealAlgebraicNumber::operator*(
     264                 :            :     const RealAlgebraicNumber& rhs) const
     265                 :            : {
     266                 :            : #ifdef CVC5_POLY_IMP
     267 [ +  + ][ -  + ]:     145211 :   if (!d_isRational || !rhs.d_isRational)
     268                 :            :   {
     269                 :        142 :     return convertToPoly(*this) * convertToPoly(rhs);
     270                 :            :   }
     271                 :            : #endif
     272                 :     290280 :   return getRationalValue() * rhs.getRationalValue();
     273                 :            : }
     274                 :     414247 : RealAlgebraicNumber RealAlgebraicNumber::operator/(
     275                 :            :     const RealAlgebraicNumber& rhs) const
     276                 :            : {
     277 [ -  + ][ -  + ]:     414247 :   Assert(!rhs.isZero()) << "Can not divide by zero";
                 [ -  - ]
     278                 :            : #ifdef CVC5_POLY_IMP
     279 [ +  + ][ +  + ]:     414247 :   if (!d_isRational || !rhs.d_isRational)
     280                 :            :   {
     281                 :         22 :     return convertToPoly(*this) / convertToPoly(rhs);
     282                 :            :   }
     283                 :            : #endif
     284                 :     828472 :   return getRationalValue() / rhs.getRationalValue();
     285                 :            : }
     286                 :            : 
     287                 :     726234 : RealAlgebraicNumber& RealAlgebraicNumber::operator+=(
     288                 :            :     const RealAlgebraicNumber& rhs)
     289                 :            : {
     290                 :            : #ifdef CVC5_POLY_IMP
     291 [ +  + ][ +  + ]:     726234 :   if (!d_isRational || !rhs.d_isRational)
     292                 :            :   {
     293                 :         78 :     getValue() = convertToPoly(*this) + convertToPoly(rhs);
     294                 :            :     // ensure it is no longer marked as rational
     295                 :         78 :     d_isRational = false;
     296                 :         78 :     return *this;
     297                 :            :   }
     298                 :            : #endif
     299                 :     726156 :   getRationalValue() = getRationalValue() + rhs.getRationalValue();
     300                 :     726156 :   return *this;
     301                 :            : }
     302                 :          0 : RealAlgebraicNumber& RealAlgebraicNumber::operator-=(
     303                 :            :     const RealAlgebraicNumber& rhs)
     304                 :            : {
     305                 :            : #ifdef CVC5_POLY_IMP
     306 [ -  - ][ -  - ]:          0 :   if (!d_isRational || !rhs.d_isRational)
     307                 :            :   {
     308                 :          0 :     getValue() = convertToPoly(*this) - convertToPoly(rhs);
     309                 :            :     // ensure it is no longer marked as rational
     310                 :          0 :     d_isRational = false;
     311                 :          0 :     return *this;
     312                 :            :   }
     313                 :            : #endif
     314                 :          0 :   getRationalValue() = getRationalValue() - rhs.getRationalValue();
     315                 :          0 :   return *this;
     316                 :            : }
     317                 :   17379100 : RealAlgebraicNumber& RealAlgebraicNumber::operator*=(
     318                 :            :     const RealAlgebraicNumber& rhs)
     319                 :            : {
     320                 :            : #ifdef CVC5_POLY_IMP
     321 [ +  + ][ +  + ]:   17379100 :   if (!d_isRational || !rhs.d_isRational)
     322                 :            :   {
     323                 :        301 :     getValue() = convertToPoly(*this) * convertToPoly(rhs);
     324                 :            :     // ensure it is no longer marked as rational
     325                 :        301 :     d_isRational = false;
     326                 :        301 :     return *this;
     327                 :            :   }
     328                 :            : #endif
     329                 :   17378800 :   getRationalValue() = getRationalValue() * rhs.getRationalValue();
     330                 :   17378800 :   return *this;
     331                 :            : }
     332                 :            : 
     333                 :    4185130 : int RealAlgebraicNumber::sgn() const
     334                 :            : {
     335                 :            : #ifdef CVC5_POLY_IMP
     336         [ +  + ]:    4185130 :   if (!d_isRational)
     337                 :            :   {
     338                 :          2 :     return poly::sgn(getValue());
     339                 :            :   }
     340                 :            : #endif
     341                 :    4185130 :   return getRationalValue().sgn();
     342                 :            : }
     343                 :   18642000 : bool RealAlgebraicNumber::isZero() const
     344                 :            : {
     345                 :            : #ifdef CVC5_POLY_IMP
     346         [ +  + ]:   18642000 :   if (!d_isRational)
     347                 :            :   {
     348                 :        210 :     return poly::is_zero(getValue());
     349                 :            :   }
     350                 :            : #endif
     351                 :   18641800 :   return getRationalValue().isZero();
     352                 :            : }
     353                 :     742997 : bool RealAlgebraicNumber::isOne() const
     354                 :            : {
     355                 :            : #ifdef CVC5_POLY_IMP
     356         [ +  + ]:     742997 :   if (!d_isRational)
     357                 :            :   {
     358                 :          1 :     return poly::is_one(getValue());
     359                 :            :   }
     360                 :            : #endif
     361                 :     742996 :   return getRationalValue().isOne();
     362                 :            : }
     363                 :          0 : RealAlgebraicNumber RealAlgebraicNumber::inverse() const
     364                 :            : {
     365                 :          0 :   Assert(!isZero()) << "Can not invert zero";
     366                 :            : #ifdef CVC5_POLY_IMP
     367         [ -  - ]:          0 :   if (!d_isRational)
     368                 :            :   {
     369                 :          0 :     return poly::inverse(getValue());
     370                 :            :   }
     371                 :            : #endif
     372                 :          0 :   return getRationalValue().inverse();
     373                 :            : }
     374                 :            : 
     375                 :       1032 : size_t RealAlgebraicNumber::hash() const
     376                 :            : {
     377                 :            : #ifdef CVC5_POLY_IMP
     378         [ +  - ]:       1032 :   if (!d_isRational)
     379                 :            :   {
     380                 :       1032 :     return lp_algebraic_number_hash_approx(getValue().get_internal(), 2);
     381                 :            :   }
     382                 :            : #endif
     383                 :          0 :   return getRationalValue().hash();
     384                 :            : }
     385                 :            : 
     386                 :          4 : std::ostream& operator<<(std::ostream& os, const RealAlgebraicNumber& ran)
     387                 :            : {
     388                 :          4 :   return os << ran.toString();
     389                 :            : }
     390                 :            : 
     391                 :            : }  // namespace cvc5::internal
     392                 :            : 
     393                 :            : namespace std {
     394                 :       1032 : size_t hash<cvc5::internal::RealAlgebraicNumber>::operator()(
     395                 :            :     const cvc5::internal::RealAlgebraicNumber& ran) const
     396                 :            : {
     397                 :       1032 :   return ran.hash();
     398                 :            : }
     399                 :            : 
     400                 :            : }  // namespace std