CVC3 2.2
|
00001 /*****************************************************************************/ 00002 /*! 00003 * \file rational-gmp.cpp 00004 * 00005 * \brief Implementation of class Rational using GMP library (C interface) 00006 * 00007 * Author: Sergey Berezin 00008 * 00009 * Created: Mon Aug 4 10:06:04 2003 00010 * 00011 * <hr> 00012 * License to use, copy, modify, sell and/or distribute this software 00013 * and its documentation for any purpose is hereby granted without 00014 * royalty, subject to the terms and conditions defined in the \ref 00015 * LICENSE file provided with this distribution. 00016 * 00017 * <hr> 00018 * 00019 */ 00020 /*****************************************************************************/ 00021 00022 #ifdef RATIONAL_GMP 00023 00024 #include "compat_hash_set.h" 00025 #include "rational.h" 00026 00027 #include <climits> 00028 #include <sstream> 00029 #include <gmp.h> 00030 #include <limits.h> 00031 00032 namespace CVC3 { 00033 00034 using namespace std; 00035 00036 //! Implementation of the forward-declared internal class 00037 class Rational::Impl { 00038 mpq_t d_n; 00039 //! Make the rational number canonical 00040 void canonicalize() { mpq_canonicalize(d_n); } 00041 public: 00042 //! Default Constructor 00043 Impl() { mpq_init(d_n); } 00044 //! Copy constructor (assumes x is canonicalized) 00045 Impl(const Impl &x) { mpq_init(d_n); mpq_set(d_n, x.d_n); } 00046 //! Constructor from mpq_t 00047 Impl(mpq_t n) { 00048 mpq_init(d_n); 00049 mpq_set(d_n, n); 00050 canonicalize(); 00051 } 00052 //! Constructor from a pair of mpz_t (integers) 00053 Impl(mpz_t n, mpz_t d) { 00054 mpq_init(d_n); 00055 mpq_set_num(d_n, n); mpq_set_den(d_n, d); 00056 canonicalize(); 00057 } 00058 //! Constructor from a single mpz_t (integer) 00059 Impl(mpz_t n) { 00060 mpq_init(d_n); 00061 mpq_set_num(d_n, n); 00062 canonicalize(); 00063 } 00064 //! Constructor from a pair of integers 00065 Impl(long int n, long int d); 00066 //! Constructor from a pair of unsigned integers 00067 Impl(unsigned int n, unsigned int d, unsigned int /* dummy arg */); 00068 //! Constructor from a string 00069 Impl(const string &n, int base); 00070 //! Constructor from a pair of strings 00071 Impl(const string &n, const string& d, int base); 00072 // Destructor 00073 virtual ~Impl() { mpq_clear(d_n); } 00074 00075 //! Assignment 00076 Impl& operator=(const Impl& x) { 00077 if(this == &x) return *this; 00078 mpq_set(d_n, x.d_n); 00079 return *this; 00080 } 00081 00082 //! Get numerator 00083 Impl getNum() const { 00084 return Impl(mpq_numref(const_cast<Impl*>(this)->d_n)); 00085 } 00086 //! Get denominator 00087 Impl getDen() const { 00088 return Impl(mpq_denref(const_cast<Impl*>(this)->d_n)); 00089 } 00090 00091 int getInt() const { 00092 // Check for overflow 00093 static Impl min((int)INT_MIN, 1), max((int)INT_MAX, 1); 00094 FatalAssert(isInteger() && min <= *this && *this <= max, 00095 "Rational::getInt(): Arithmetic overflow for " 00096 +toString()); 00097 return mpz_get_si(mpq_numref(d_n)); 00098 } 00099 00100 unsigned int getUnsigned() const { 00101 // Check for overflow 00102 static Impl min(0, 1, 0), max(UINT_MAX, 1, 0); 00103 FatalAssert(min <= *this && *this <= max, 00104 "Rational::getUnsigned(): Arithmetic overflow for " 00105 +toString()); 00106 return mpz_get_ui(mpq_numref(d_n)); 00107 } 00108 00109 //! Unary minus 00110 Impl operator-() const; 00111 //! Get the floor 00112 Impl floor() const; 00113 //! Get the ceiling 00114 Impl ceil() const; 00115 00116 //! Testing whether the denominator is 1 00117 bool isInteger() const; 00118 00119 //! Equality 00120 friend bool operator==(const Impl& x, const Impl& y) { 00121 return mpq_equal(x.d_n, y.d_n); 00122 } 00123 00124 //! Dis-equality 00125 friend bool operator!=(const Impl& x, const Impl& y) { 00126 return !mpq_equal(x.d_n, y.d_n); 00127 } 00128 //! Less than 00129 friend bool operator<(const Impl& x, const Impl& y) { 00130 return mpq_cmp(x.d_n, y.d_n) < 0; 00131 } 00132 00133 friend bool operator<=(const Impl& x, const Impl& y) { 00134 return mpq_cmp(x.d_n, y.d_n) <= 0; 00135 } 00136 00137 friend bool operator>(const Impl& x, const Impl& y) { 00138 return mpq_cmp(x.d_n, y.d_n) > 0; 00139 } 00140 00141 friend bool operator>=(const Impl& x, const Impl& y) { 00142 return mpq_cmp(x.d_n, y.d_n) >= 0; 00143 } 00144 00145 //! Addition 00146 friend Impl operator+(const Impl& x, const Impl& y) { 00147 Impl res; 00148 mpq_add(res.d_n, x.d_n, y.d_n); 00149 return res; 00150 } 00151 00152 //! Subtraction 00153 friend Impl operator-(const Impl& x, const Impl& y) { 00154 Impl res; 00155 mpq_sub(res.d_n, x.d_n, y.d_n); 00156 return res; 00157 } 00158 00159 //! Multiplication 00160 friend Impl operator*(const Impl& x, const Impl& y) { 00161 Impl res; 00162 mpq_mul(res.d_n, x.d_n, y.d_n); 00163 return res; 00164 } 00165 00166 //! Division 00167 friend Impl operator/(const Impl& x, const Impl& y) { 00168 Impl res; 00169 mpq_div(res.d_n, x.d_n, y.d_n); 00170 return res; 00171 } 00172 00173 friend Impl operator%(const Impl& x, const Impl& y) { 00174 DebugAssert(x.isInteger() && y.isInteger(), 00175 "Impl % Impl: x and y must be integers"); 00176 mpz_t res; 00177 mpz_init(res); 00178 // Put the remainder directly into 'res' 00179 mpz_fdiv_r(res, mpq_numref(x.d_n), mpq_numref(y.d_n)); 00180 Impl r(res); 00181 mpz_clear(res); 00182 return r; 00183 } 00184 00185 //! Compute the remainder of x/y 00186 friend Impl mod(const Impl& x, const Impl& y) { 00187 DebugAssert(x.isInteger() && y.isInteger(), 00188 "Rational::Impl::mod(): x="+x.toString() 00189 +", y="+y.toString()); 00190 mpz_t res; 00191 mpz_init(res); 00192 mpz_mod(res, mpq_numref(x.d_n), mpq_numref(y.d_n)); 00193 Impl r(res); 00194 mpz_clear(res); 00195 return r; 00196 } 00197 00198 friend Impl intRoot(const Impl& x, unsigned long int y) { 00199 DebugAssert(x.isInteger(), 00200 "Rational::Impl::intRoot(): x="+x.toString()); 00201 mpz_t res; 00202 mpz_init(res); 00203 int exact = mpz_root(res, mpq_numref(x.d_n), y); 00204 if (!exact) { 00205 mpz_set_ui(res, 0); 00206 } 00207 Impl r(res); 00208 mpz_clear(res); 00209 return r; 00210 } 00211 00212 friend Impl gcd(const Impl& x, const Impl& y) { 00213 DebugAssert(x.isInteger() && y.isInteger(), 00214 "Rational::Impl::gcd(): x="+x.toString() 00215 +", y="+y.toString()); 00216 TRACE("rational", "Impl::gcd(", x, ", "+y.toString()+") {"); 00217 mpz_t res; 00218 mpz_init(res); 00219 mpz_gcd(res, mpq_numref(x.d_n), mpq_numref(y.d_n)); 00220 Impl r(res); 00221 mpz_clear(res); 00222 TRACE("rational", "Impl::gcd => ", r, "}"); 00223 return r; 00224 } 00225 00226 friend Impl lcm(const Impl& x, const Impl& y) { 00227 DebugAssert(x.isInteger() && y.isInteger(), 00228 "Rational::Impl::lcm(): x="+x.toString() 00229 +", y="+y.toString()); 00230 mpz_t res; 00231 mpz_init(res); 00232 mpz_lcm(res, mpq_numref(x.d_n), mpq_numref(y.d_n)); 00233 Impl r(res); 00234 mpz_clear(res); 00235 return r; 00236 } 00237 00238 //! Print to string 00239 string toString(int base = 10) const { 00240 char* str = (char*)malloc(mpz_sizeinbase(mpq_numref(d_n), base) 00241 +mpz_sizeinbase(mpq_denref(d_n), base)+3); 00242 mpq_get_str(str, base, d_n); 00243 string res(str); 00244 free(str); 00245 return res; 00246 } 00247 00248 //! Printing to ostream 00249 friend ostream& operator<<(ostream& os, const Rational::Impl& n) { 00250 return os << n.toString(); 00251 } 00252 00253 }; 00254 00255 // Constructor from a pair of integers 00256 Rational::Impl::Impl(long int n, long int d) { 00257 mpq_init(d_n); 00258 DebugAssert(d > 0, "Rational::Impl(long n, long d): d = "+int2string(d)); 00259 mpq_set_si(d_n, n, (unsigned long int)d); 00260 canonicalize(); 00261 } 00262 00263 // Constructor from a pair of unsigned integers 00264 Rational::Impl::Impl(unsigned int n, unsigned int d, 00265 unsigned int /* dummy arg, to disambiguate */) { 00266 mpq_init(d_n); 00267 mpq_set_ui(d_n, n, (unsigned long int)d); 00268 canonicalize(); 00269 } 00270 00271 // Constructor from a string 00272 Rational::Impl::Impl(const string &n, int base) { 00273 mpq_init(d_n); 00274 mpq_set_str(d_n, n.c_str(), base); 00275 canonicalize(); 00276 } 00277 00278 // Constructor from a pair of strings 00279 Rational::Impl::Impl(const string &n, const string& d, int base) { 00280 mpq_init(d_n); 00281 mpq_set_str(d_n, (n+"/"+d).c_str(), base); 00282 canonicalize(); 00283 } 00284 00285 Rational::Impl Rational::Impl::operator-() const { 00286 Impl res; 00287 mpq_neg(res.d_n, d_n); 00288 return res; 00289 } 00290 00291 Rational::Impl Rational::Impl::floor() const { 00292 mpz_t res; 00293 mpz_init(res); 00294 mpz_fdiv_q(res, mpq_numref(d_n), mpq_denref(d_n)); 00295 Impl r(res); 00296 mpz_clear(res); 00297 return r; 00298 } 00299 00300 Rational::Impl Rational::Impl::ceil() const { 00301 mpz_t res; 00302 mpz_init(res); 00303 mpz_cdiv_q(res, mpq_numref(d_n), mpq_denref(d_n)); 00304 Impl r(res); 00305 mpz_clear(res); 00306 return r; 00307 } 00308 00309 bool Rational::Impl::isInteger() const { 00310 bool res(mpz_cmp_ui(mpq_denref(d_n), 1) == 0); 00311 TRACE("rational", "Impl::isInteger(", *this, 00312 ") => "+string(res? "true" : "false")); 00313 return res; 00314 } 00315 00316 //! Default constructor 00317 Rational::Rational() : d_n(new Impl) { 00318 #ifdef _DEBUG_RATIONAL_ 00319 int &num_created = getCreated(); 00320 num_created++; 00321 printStats(); 00322 #endif 00323 } 00324 //! Copy constructor 00325 Rational::Rational(const Rational &n) : d_n(new Impl(*n.d_n)) { 00326 #ifdef _DEBUG_RATIONAL_ 00327 int &num_created = getCreated(); 00328 num_created++; 00329 printStats(); 00330 #endif 00331 } 00332 00333 //! Private constructor 00334 Rational::Rational(const Impl& t): d_n(new Impl(t)) { 00335 #ifdef _DEBUG_RATIONAL_ 00336 int &num_created = getCreated(); 00337 num_created++; 00338 printStats(); 00339 #endif 00340 } 00341 Rational::Rational(int n, int d): d_n(new Impl(n, d)) { 00342 #ifdef _DEBUG_RATIONAL_ 00343 int &num_created = getCreated(); 00344 num_created++; 00345 printStats(); 00346 #endif 00347 } 00348 // Constructors from strings 00349 Rational::Rational(const char* n, int base) 00350 : d_n(new Impl(string(n), base)) { 00351 #ifdef _DEBUG_RATIONAL_ 00352 int &num_created = getCreated(); 00353 num_created++; 00354 printStats(); 00355 #endif 00356 } 00357 Rational::Rational(const string& n, int base) 00358 : d_n(new Impl(n, base)) { 00359 #ifdef _DEBUG_RATIONAL_ 00360 int &num_created = getCreated(); 00361 num_created++; 00362 printStats(); 00363 #endif 00364 } 00365 Rational::Rational(const char* n, const char* d, int base) 00366 : d_n(new Impl(string(n), string(d), base)) { 00367 #ifdef _DEBUG_RATIONAL_ 00368 int &num_created = getCreated(); 00369 num_created++; 00370 printStats(); 00371 #endif 00372 } 00373 Rational::Rational(const string& n, const string& d, int base) 00374 : d_n(new Impl(n, d, base)) { 00375 #ifdef _DEBUG_RATIONAL_ 00376 int &num_created = getCreated(); 00377 num_created++; 00378 printStats(); 00379 #endif 00380 } 00381 // Destructor 00382 Rational::~Rational() { 00383 #ifdef _DEBUG_RATIONAL_ 00384 int &num_deleted = getDeleted(); 00385 num_deleted++; 00386 printStats(); 00387 #endif 00388 delete d_n; 00389 } 00390 00391 // Get components 00392 Rational Rational::getNumerator() const 00393 { return Rational(d_n->getNum()); } 00394 Rational Rational::getDenominator() const 00395 { return Rational(d_n->getDen()); } 00396 00397 bool Rational::isInteger() const { return d_n->isInteger(); } 00398 00399 // Assignment 00400 Rational& Rational::operator=(const Rational& n) { 00401 if(this == &n) return *this; 00402 delete d_n; 00403 d_n = new Impl(*n.d_n); 00404 return *this; 00405 } 00406 00407 ostream &operator<<(ostream &os, const Rational &n) { 00408 return(os << n.toString()); 00409 } 00410 00411 00412 // Check that argument is an int and print an error message otherwise 00413 00414 static void checkInt(const Rational& n, const string& funName) { 00415 TRACE("rational", "checkInt(", n, ")"); 00416 DebugAssert(n.isInteger(), 00417 "CVC3::Rational::" + funName 00418 + ": argument is not an integer: " + n.toString()); 00419 } 00420 00421 /* Computes gcd and lcm on *integer* values. Result is always a 00422 positive integer. In this implementation, it is guaranteed by 00423 GMP. */ 00424 00425 Rational gcd(const Rational &x, const Rational &y) { 00426 checkInt(x, "gcd(*x*,y)"); 00427 checkInt(y, "gcd(x,*y*)"); 00428 return Rational(gcd(*x.d_n, *y.d_n)); 00429 } 00430 00431 Rational gcd(const vector<Rational> &v) { 00432 Rational::Impl g(1,1), zero; 00433 if(v.size() > 0) { 00434 checkInt(v[0], "gcd(vector<Rational>[0])"); 00435 g = *v[0].d_n; 00436 } 00437 for(size_t i=1; i<v.size(); i++) { 00438 checkInt(v[i], "gcd(vector<Rational>)"); 00439 if(g == zero) 00440 g = *(v[i].d_n); 00441 else if(*(v[i].d_n) != zero) 00442 g = gcd(g, *(v[i].d_n)); 00443 } 00444 return Rational(g); 00445 } 00446 00447 Rational lcm(const Rational &x, const Rational &y) { 00448 checkInt(x, "lcm(*x*,y)"); 00449 checkInt(y, "lcm(x,*y*)"); 00450 return Rational(lcm(*x.d_n, *y.d_n)); 00451 } 00452 00453 Rational lcm(const vector<Rational> &v) { 00454 Rational::Impl g(1,1), zero; 00455 for(size_t i=0; i<v.size(); i++) { 00456 checkInt(v[i], "lcm(vector<Rational>)"); 00457 if(*v[i].d_n != zero) 00458 g = lcm(g, *v[i].d_n); 00459 } 00460 return Rational(g); 00461 } 00462 00463 Rational abs(const Rational &x) { 00464 if(x<0) return -x; 00465 return x; 00466 } 00467 00468 Rational floor(const Rational &x) { 00469 return Rational(x.d_n->floor()); 00470 } 00471 00472 Rational ceil(const Rational &x) { 00473 return Rational(x.d_n->ceil()); 00474 } 00475 00476 Rational mod(const Rational &x, const Rational &y) { 00477 checkInt(x, "mod(*x*,y)"); 00478 checkInt(y, "mod(x,*y*)"); 00479 return(Rational(mod(*x.d_n, *y.d_n))); 00480 } 00481 00482 Rational intRoot(const Rational& base, unsigned long int n) { 00483 checkInt(base, "intRoot(*x*,y)"); 00484 return Rational(intRoot(*base.d_n, n)); 00485 } 00486 00487 string Rational::toString(int base) const { 00488 return(d_n->toString(base)); 00489 } 00490 00491 size_t Rational::hash() const { 00492 std::hash<const char *> h; 00493 return h(toString().c_str()); 00494 } 00495 00496 void Rational::print() const { 00497 cout << (*this) << endl; 00498 } 00499 00500 // Unary minus 00501 Rational Rational::operator-() const { 00502 return Rational(-(*d_n)); 00503 } 00504 00505 Rational &Rational::operator+=(const Rational &n2) { 00506 *this = (*this) + n2; 00507 return *this; 00508 } 00509 00510 Rational &Rational::operator-=(const Rational &n2) { 00511 *this = (*this) - n2; 00512 return *this; 00513 } 00514 00515 Rational &Rational::operator*=(const Rational &n2) { 00516 *this = (*this) * n2; 00517 return *this; 00518 } 00519 00520 Rational &Rational::operator/=(const Rational &n2) { 00521 *this = (*this) / n2; 00522 return *this; 00523 } 00524 00525 int Rational::getInt() const { 00526 checkInt(*this, "getInt()"); 00527 return d_n->getInt(); 00528 } 00529 00530 unsigned int Rational::getUnsigned() const { 00531 checkInt(*this, "getUnsigned()"); 00532 return d_n->getUnsigned(); 00533 } 00534 00535 #ifdef _DEBUG_RATIONAL_ 00536 void Rational::printStats() { 00537 int &num_created = getCreated(); 00538 int &num_deleted = getDeleted(); 00539 if(num_created % 1000 == 0 || num_deleted % 1000 == 0) { 00540 std::cerr << "Rational(" << *d_n << "): created " << num_created 00541 << ", deleted " << num_deleted 00542 << ", currently alive " << num_created-num_deleted 00543 << std::endl; 00544 } 00545 } 00546 #endif 00547 00548 bool operator==(const Rational &n1, const Rational &n2) { 00549 return(*n1.d_n == *n2.d_n); 00550 } 00551 00552 bool operator<(const Rational &n1, const Rational &n2) { 00553 return(*n1.d_n < *n2.d_n); 00554 } 00555 00556 bool operator<=(const Rational &n1, const Rational &n2) { 00557 return(*n1.d_n <= *n2.d_n); 00558 } 00559 00560 bool operator>(const Rational &n1, const Rational &n2) { 00561 return(*n1.d_n > *n2.d_n); 00562 } 00563 00564 bool operator>=(const Rational &n1, const Rational &n2) { 00565 return(*n1.d_n >= *n2.d_n); 00566 } 00567 00568 bool operator!=(const Rational &n1, const Rational &n2) { 00569 return(*n1.d_n != *n2.d_n); 00570 } 00571 00572 Rational operator+(const Rational &n1, const Rational &n2) { 00573 return Rational((*n1.d_n) + (*n2.d_n)); 00574 } 00575 00576 Rational operator-(const Rational &n1, const Rational &n2) { 00577 return Rational((*n1.d_n) - (*n2.d_n)); 00578 } 00579 00580 Rational operator*(const Rational &n1, const Rational &n2) { 00581 return Rational((*n1.d_n) * (*n2.d_n)); 00582 } 00583 00584 Rational operator/(const Rational &n1, const Rational &n2) { 00585 return Rational((*n1.d_n) / (*n2.d_n)); 00586 } 00587 00588 Rational operator%(const Rational &n1, const Rational &n2) { 00589 return Rational((*n1.d_n) % (*n2.d_n)); 00590 } 00591 00592 } /* close namespace */ 00593 00594 #endif