CVC3 2.2
|
00001 /*****************************************************************************/ 00002 /*! 00003 * \file rational-native.cpp 00004 * 00005 * \brief Implementation of class Rational using native (bounded 00006 * precision) computer arithmetic 00007 * 00008 * Author: Sergey Berezin 00009 * 00010 * Created: Mon Jul 28 12:18:03 2003 00011 * 00012 * <hr> 00013 * License to use, copy, modify, sell and/or distribute this software 00014 * and its documentation for any purpose is hereby granted without 00015 * royalty, subject to the terms and conditions defined in the \ref 00016 * LICENSE file provided with this distribution. 00017 * 00018 * <hr> 00019 * 00020 */ 00021 /*****************************************************************************/ 00022 00023 #ifdef RATIONAL_NATIVE 00024 00025 #include "compat_hash_set.h" 00026 #include "rational.h" 00027 // For atol() (ASCII to long) 00028 #include <stdlib.h> 00029 #include <limits.h> 00030 #include <errno.h> 00031 #include <sstream> 00032 #include <math.h> 00033 00034 #define OVERFLOW_MSG "\nThis is NOT a bug, but an explicit feature to preserve soundness\nwhen CVC3 uses native computer arithmetic (finite precision). To\navoid these types of errors, please recompile CVC3 with GMP library." 00035 00036 namespace CVC3 { 00037 00038 using namespace std; 00039 00040 //! Add two integers and check for overflows 00041 static long int plus(long int x, long int y) { 00042 long int res = x+y; 00043 FatalAssert(((x > 0) != (y > 0)) || ((x > 0) == (res > 0)), 00044 "plus(x,y): arithmetic overflow" OVERFLOW_MSG); 00045 return res; 00046 } 00047 00048 //! Unary minus which checks for overflows 00049 static long int uminus(long int x) { 00050 FatalAssert(x == 0 || x != -x, "uminus(x): arithmetic overflow" 00051 OVERFLOW_MSG); 00052 return -x; 00053 } 00054 00055 //! Multiply two integers and check for overflows 00056 static long int mult(long int x, long int y) { 00057 long int res = x*y; 00058 FatalAssert(x==0 || res/x == y, "mult(x,y): arithmetic overflow" 00059 OVERFLOW_MSG); 00060 return res; 00061 } 00062 00063 //! Compute GCD using Euclid's algorithm (from Aaron Stump's code) 00064 static long int gcd(long int n1, long int n2) { 00065 DebugAssert(n1!=0 && n2!=0, 00066 "gcd("+int2string(n1)+", "+int2string(n2)+"): bad args"); 00067 // First, make the arguments positive 00068 if(n1 < 0) n1 = uminus(n1); 00069 if(n2 < 0) n2 = uminus(n2); 00070 00071 if (n1 < n2) { 00072 long int tmp = n1; 00073 n1 = n2; 00074 n2 = tmp; 00075 } 00076 00077 // at this point, n1 >= n2 00078 long int r = n1 % n2; 00079 if (!r) 00080 return n2; 00081 00082 return gcd(n2, r); 00083 } 00084 00085 //! Compute LCM 00086 static long int lcm(long int n1, long int n2) { 00087 long int g = gcd(n1,n2); 00088 return mult(n1/g, n2); 00089 } 00090 00091 // Implementation of the forward-declared internal class 00092 class Rational::Impl { 00093 long int d_num; //!< Numerator 00094 long int d_denom; //!< Denominator 00095 //! Make the rational number canonical 00096 void canonicalize(); 00097 public: 00098 //! Default Constructor 00099 Impl(): d_num(0), d_denom(1) { } 00100 //! Copy constructor 00101 Impl(const Impl &x) : d_num(x.d_num), d_denom(x.d_denom) { } 00102 //! Constructor from a pair of integers 00103 Impl(long int n, long int d): d_num(n), d_denom(d) { canonicalize(); } 00104 //! Constructor from a string 00105 Impl(const string &n, int base); 00106 //! Constructor from a pair of strings 00107 Impl(const string &n, const string& d, int base); 00108 // Destructor 00109 virtual ~Impl() { } 00110 //! Get numerator 00111 long int getNum() const { return d_num; } 00112 //! Get denominator 00113 long int getDen() const { return d_denom; } 00114 00115 //! Unary minus 00116 Impl operator-() const; 00117 00118 //! Equality 00119 friend bool operator==(const Impl& x, const Impl& y) { 00120 return (x.d_num == y.d_num && x.d_denom == y.d_denom); 00121 } 00122 //! Dis-equality 00123 friend bool operator!=(const Impl& x, const Impl& y) { 00124 return (x.d_num != y.d_num || x.d_denom != y.d_denom); 00125 } 00126 /*! 00127 * Compare x=n1/d1 and y=n2/d2 as n1*f2 < n2*f1, where f1=d1/f, 00128 * f2=d2/f, and f=lcm(d1,d2) 00129 */ 00130 friend bool operator<(const Impl& x, const Impl& y) { 00131 Impl diff(x-y); 00132 return diff.d_num < 0; 00133 } 00134 00135 friend bool operator<=(const Impl& x, const Impl& y) { 00136 return (x == y || x < y); 00137 } 00138 00139 friend bool operator>(const Impl& x, const Impl& y) { 00140 Impl diff(x-y); 00141 return diff.d_num > 0; 00142 } 00143 00144 friend bool operator>=(const Impl& x, const Impl& y) { 00145 return (x == y || x > y); 00146 } 00147 00148 /*! Addition of x=n1/d1 and y=n2/d2: n1*g2 + n2*g1, where g1=d1/g, 00149 * g2=d2/g, and g=lcm(d1,d2) 00150 */ 00151 friend Impl operator+(const Impl& x, const Impl& y) { 00152 long int d1(x.getDen()), d2(y.getDen()); 00153 long int f(lcm(d1,d2)), f1(f/d1), f2(f/d2); 00154 long int n = plus(mult(x.getNum(), f1), mult(y.getNum(), f2)); 00155 return Impl(n, f); 00156 } 00157 00158 friend Impl operator-(const Impl& x, const Impl& y) { 00159 TRACE("rational", "operator-(", x, ", "+y.toString()+")"); 00160 long int d1(x.getDen()), d2(y.getDen()); 00161 long int f(lcm(d1,d2)), f1(f/d1), f2(f/d2); 00162 long int n = plus(mult(x.getNum(), f1), uminus(mult(y.getNum(), f2))); 00163 Impl res(n, f); 00164 TRACE("rational", " => ", res, ""); 00165 return res; 00166 } 00167 00168 /*! 00169 * Multiplication of x=n1/d1, y=n2/d2: 00170 * (n1/g1)*(n2/g2)/(d1/g2)*(d2/g1), where g1=gcd(n1,d2) and 00171 * g2=gcd(n2,d1) 00172 */ 00173 friend Impl operator*(const Impl& x, const Impl& y) { 00174 long int n1(x.getNum()), d1(x.getDen()), n2(y.getNum()), d2(y.getDen()); 00175 long int g1(n1? gcd(n1,d2) : 1), g2(n2? gcd(n2,d1) : 1); 00176 long int n(mult(n1/g1, n2/g2)), d(mult(d1/g2, d2/g1)); 00177 return Impl(n,d); 00178 } 00179 00180 /*! 00181 * Division of x=n1/d1, y=n2/d2: 00182 * (n1/g1)*(d2/g2)/(d1/g2)*(n2/g1), where g1=gcd(n1,n2) and 00183 * g2=gcd(d1,d2) 00184 */ 00185 friend Impl operator/(const Impl& x, const Impl& y) { 00186 long int n1(x.getNum()), d1(x.getDen()), n2(y.getNum()), d2(y.getDen()); 00187 DebugAssert(n2 != 0, "Impl::operator/: divisor is 0"); 00188 long int g1(n1? gcd(n1,n2) : 1), g2(gcd(d1,d2)); 00189 long int n(n1? mult(n1/g1, d2/g2) : 0), d(n1? mult(d1/g2, n2/g1) : 1); 00190 return Impl(n,d); 00191 } 00192 00193 friend Impl operator%(const Impl& x, const Impl& y) { 00194 DebugAssert(x.getDen() == 1 && y.getDen() == 1, 00195 "Impl % Impl: x and y must be integers"); 00196 return Impl(x.getNum() % y.getNum(), 1); 00197 } 00198 00199 //! Print to string 00200 string toString(int base = 10) const { 00201 ostringstream ss; 00202 if (d_num == 0) ss << "0"; 00203 else if (base == 10) { 00204 ss << d_num; 00205 if (d_denom != 1) 00206 ss << "/" << d_denom; 00207 } 00208 else { 00209 vector<int> vec; 00210 long num = d_num; 00211 while (num) { 00212 vec.push_back(num % base); 00213 num = num / base; 00214 } 00215 while (!vec.empty()) { 00216 if (base > 10 && vec.back() > 10) { 00217 ss << (char)('A' + (vec.back()-10)); 00218 } 00219 else ss << vec.back(); 00220 vec.pop_back(); 00221 } 00222 if(d_denom != 1) { 00223 ss << "/"; 00224 if (d_denom == 0) ss << "0"; 00225 else { 00226 num = d_denom; 00227 while (num) { 00228 vec.push_back(num % base); 00229 num = num / base; 00230 } 00231 while (!vec.empty()) { 00232 if (base > 10 && vec.back() > 10) { 00233 ss << (char)('A' + (vec.back()-10)); 00234 } 00235 else ss << vec.back(); 00236 vec.pop_back(); 00237 } 00238 } 00239 } 00240 } 00241 return(ss.str()); 00242 } 00243 00244 //! Printing to ostream 00245 friend ostream& operator<<(ostream& os, const Rational::Impl& n) { 00246 return os << n.toString(); 00247 } 00248 00249 }; 00250 00251 // Make the rational number canonical 00252 void Rational::Impl::canonicalize() { 00253 DebugAssert(d_denom != 0, 00254 "Rational::Impl::canonicalize: bad denominator: " 00255 +int2string(d_denom)); 00256 if(d_num == 0) { 00257 d_denom = 1; 00258 } else { 00259 if(d_denom < 0) { 00260 d_num = uminus(d_num); 00261 d_denom = uminus(d_denom); 00262 } 00263 long int d = gcd(d_num, d_denom); 00264 if(d != 1) { 00265 d_num /= d; 00266 d_denom /= d; 00267 } 00268 } 00269 } 00270 00271 // Constructor from a string 00272 Rational::Impl::Impl(const string &n, int base) { 00273 size_t i, iend; 00274 for(i=0,iend=n.size(); i<iend && n[i] != '/'; ++i); 00275 if(i<iend) 00276 // Found slash at i 00277 *this = Impl(n.substr(0,i-1), n.substr(i+1,iend-1), base); 00278 else 00279 *this = Impl(n, "1", base); 00280 canonicalize(); 00281 } 00282 // Constructor from a pair of strings 00283 Rational::Impl::Impl(const string &n, const string& d, int base) { 00284 d_num = strtol(n.c_str(), NULL, base); 00285 FatalAssert(d_num != LONG_MIN && d_num != LONG_MAX, 00286 "Rational::Impl(string, string): arithmetic overflow:" 00287 "n = "+n+", d="+d+", base="+int2string(base) 00288 +OVERFLOW_MSG); 00289 d_denom = strtol(d.c_str(), NULL, base); 00290 FatalAssert(d_denom != LONG_MIN && d_denom != LONG_MAX, 00291 "Rational::Impl(string, string): arithmetic overflow:" 00292 "n = "+n+", d="+d+", base="+int2string(base) 00293 +OVERFLOW_MSG); 00294 canonicalize(); 00295 } 00296 00297 Rational::Impl Rational::Impl::operator-() const { 00298 return Impl(uminus(d_num), d_denom); 00299 } 00300 00301 00302 //! Default constructor 00303 Rational::Rational() : d_n(new Impl) { 00304 #ifdef _DEBUG_RATIONAL_ 00305 int &num_created = getCreated(); 00306 num_created++; 00307 printStats(); 00308 #endif 00309 } 00310 //! Copy constructor 00311 Rational::Rational(const Rational &n) : d_n(new Impl(*n.d_n)) { 00312 #ifdef _DEBUG_RATIONAL_ 00313 int &num_created = getCreated(); 00314 num_created++; 00315 printStats(); 00316 #endif 00317 } 00318 00319 //! Private constructor 00320 Rational::Rational(const Impl& t): d_n(new Impl(t)) { 00321 #ifdef _DEBUG_RATIONAL_ 00322 int &num_created = getCreated(); 00323 num_created++; 00324 printStats(); 00325 #endif 00326 } 00327 Rational::Rational(int n, int d): d_n(new Impl(n, d)) { 00328 #ifdef _DEBUG_RATIONAL_ 00329 int &num_created = getCreated(); 00330 num_created++; 00331 printStats(); 00332 #endif 00333 } 00334 // Constructors from strings 00335 Rational::Rational(const char* n, int base) 00336 : d_n(new Impl(string(n), base)) { 00337 #ifdef _DEBUG_RATIONAL_ 00338 int &num_created = getCreated(); 00339 num_created++; 00340 printStats(); 00341 #endif 00342 } 00343 Rational::Rational(const string& n, int base) 00344 : d_n(new Impl(n, base)) { 00345 #ifdef _DEBUG_RATIONAL_ 00346 int &num_created = getCreated(); 00347 num_created++; 00348 printStats(); 00349 #endif 00350 } 00351 Rational::Rational(const char* n, const char* d, int base) 00352 : d_n(new Impl(string(n), string(d), base)) { 00353 #ifdef _DEBUG_RATIONAL_ 00354 int &num_created = getCreated(); 00355 num_created++; 00356 printStats(); 00357 #endif 00358 } 00359 Rational::Rational(const string& n, const string& d, int base) 00360 : d_n(new Impl(n, d, base)) { 00361 #ifdef _DEBUG_RATIONAL_ 00362 int &num_created = getCreated(); 00363 num_created++; 00364 printStats(); 00365 #endif 00366 } 00367 // Destructor 00368 Rational::~Rational() { 00369 delete d_n; 00370 #ifdef _DEBUG_RATIONAL_ 00371 int &num_deleted = getDeleted(); 00372 num_deleted++; 00373 printStats(); 00374 #endif 00375 } 00376 00377 // Get components 00378 Rational Rational::getNumerator() const 00379 { return Rational(Impl(d_n->getNum(), 1)); } 00380 Rational Rational::getDenominator() const 00381 { return Rational(Impl(d_n->getDen(), 1)); } 00382 00383 bool Rational::isInteger() const { return 1 == d_n->getDen(); } 00384 00385 // Assignment 00386 Rational& Rational::operator=(const Rational& n) { 00387 if(this == &n) return *this; 00388 delete d_n; 00389 d_n = new Impl(*n.d_n); 00390 return *this; 00391 } 00392 00393 ostream &operator<<(ostream &os, const Rational &n) { 00394 return(os << n.toString()); 00395 } 00396 00397 00398 // Check that argument is an int and print an error message otherwise 00399 00400 static void checkInt(const Rational& n, const string& funName) { 00401 DebugAssert(n.isInteger(), 00402 ("CVC3::Rational::" + funName 00403 + ": argument is not an integer: " + n.toString()).c_str()); 00404 } 00405 00406 /* Computes gcd and lcm on *integer* values. Result is always a 00407 positive integer. In this implementation, it is guaranteed by 00408 GMP. */ 00409 00410 Rational gcd(const Rational &x, const Rational &y) { 00411 checkInt(x, "gcd(*x*,y)"); 00412 checkInt(y, "gcd(x,*y*)"); 00413 return Rational(Rational::Impl(gcd(x.d_n->getNum(), y.d_n->getNum()), 1)); 00414 } 00415 00416 Rational gcd(const vector<Rational> &v) { 00417 long int g(1); 00418 if(v.size() > 0) { 00419 checkInt(v[0], "gcd(vector<Rational>[0])"); 00420 g = v[0].d_n->getNum(); 00421 } 00422 for(size_t i=1; i<v.size(); i++) { 00423 checkInt(v[i], "gcd(vector<Rational>)"); 00424 if(g == 0) g = v[i].d_n->getNum(); 00425 else if(v[i].d_n->getNum() != 0) 00426 g = gcd(g, v[i].d_n->getNum()); 00427 } 00428 return Rational(Rational::Impl(g,1)); 00429 } 00430 00431 Rational lcm(const Rational &x, const Rational &y) { 00432 long int g; 00433 checkInt(x, "lcm(*x*,y)"); 00434 checkInt(y, "lcm(x,*y*)"); 00435 g = lcm(x.d_n->getNum(), y.d_n->getNum()); 00436 return Rational(Rational::Impl(g, 1)); 00437 } 00438 00439 Rational lcm(const vector<Rational> &v) { 00440 long int g(1); 00441 for(size_t i=0; i<v.size(); i++) { 00442 checkInt(v[i], "lcm(vector<Rational>)"); 00443 if(v[i].d_n->getNum() != 0) 00444 g = lcm(g, v[i].d_n->getNum()); 00445 } 00446 return Rational(Rational::Impl(g,1)); 00447 } 00448 00449 Rational abs(const Rational &x) { 00450 long int n(x.d_n->getNum()); 00451 if(n>=0) return x; 00452 return Rational(Rational::Impl(-n, x.d_n->getDen())); 00453 } 00454 00455 Rational floor(const Rational &x) { 00456 if(x.d_n->getDen() == 1) return x; 00457 long int n = x.d_n->getNum(); 00458 long int nAbs = (n<0)? uminus(n) : n; 00459 long int q = nAbs / x.d_n->getDen(); 00460 if(n < 0) q = plus(uminus(q), -1); 00461 return Rational(Rational::Impl(q,1)); 00462 } 00463 00464 Rational ceil(const Rational &x) { 00465 if(x.d_n->getDen() == 1) return x; 00466 long int n = x.d_n->getNum(); 00467 long int nAbs = (n<0)? -n : n; 00468 long int q = nAbs / x.d_n->getDen(); 00469 if(n > 0) q = plus(q, 1); 00470 else q = uminus(q); 00471 return Rational(Rational::Impl(q,1)); 00472 } 00473 00474 Rational mod(const Rational &x, const Rational &y) { 00475 checkInt(x, "mod(*x*,y)"); 00476 checkInt(y, "mod(x,*y*)"); 00477 long int r = x.d_n->getNum() % y.d_n->getNum(); 00478 return(Rational(Rational::Impl(r,1))); 00479 } 00480 00481 Rational intRoot(const Rational& base, unsigned long int n) { 00482 checkInt(base, "intRoot(*x*,y)"); 00483 checkInt(n, "intRoot(x,*y*)"); 00484 double b = base.d_n->getNum(); 00485 double root = n; 00486 root = 1/root; 00487 b = ::pow(b, root); 00488 long res = (long) ::floor(b); 00489 if (::pow((long double)res, (int)n) == base.d_n->getNum()) { 00490 return Rational(Rational::Impl(res,1)); 00491 } 00492 return Rational(Rational::Impl((long)0,(long)1)); 00493 } 00494 00495 string Rational::toString(int base) const { 00496 return(d_n->toString(base)); 00497 } 00498 00499 size_t Rational::hash() const { 00500 std::hash<const char *> h; 00501 return h(toString().c_str()); 00502 } 00503 00504 void Rational::print() const { 00505 cout << (*this) << endl; 00506 } 00507 00508 // Unary minus 00509 Rational Rational::operator-() const { 00510 return Rational(Rational::Impl(-(d_n->getNum()), d_n->getDen())); 00511 } 00512 00513 Rational &Rational::operator+=(const Rational &n2) { 00514 *this = (*this) + n2; 00515 return *this; 00516 } 00517 00518 Rational &Rational::operator-=(const Rational &n2) { 00519 *this = (*this) - n2; 00520 return *this; 00521 } 00522 00523 Rational &Rational::operator*=(const Rational &n2) { 00524 *this = (*this) * n2; 00525 return *this; 00526 } 00527 00528 Rational &Rational::operator/=(const Rational &n2) { 00529 *this = (*this) / n2; 00530 return *this; 00531 } 00532 00533 int Rational::getInt() const { 00534 checkInt(*this, "getInt()"); 00535 long int res = d_n->getNum(); 00536 FatalAssert(res >= INT_MIN && res <= INT_MAX, 00537 "Rational::getInt(): arithmetic overflow on "+toString() + 00538 OVERFLOW_MSG); 00539 return (int)res; 00540 } 00541 00542 unsigned int Rational::getUnsigned() const { 00543 checkInt(*this, "getUnsigned()"); 00544 long int res = d_n->getNum(); 00545 FatalAssert(res >= 0 && res <= (long int)UINT_MAX, 00546 "Rational::getUnsigned(): arithmetic overflow on " + toString() + 00547 OVERFLOW_MSG); 00548 return (unsigned int)res; 00549 } 00550 00551 #ifdef _DEBUG_RATIONAL_ 00552 void Rational::printStats() { 00553 int &num_created = getCreated(); 00554 int &num_deleted = getDeleted(); 00555 if(num_created % 1000 == 0 || num_deleted % 1000 == 0) { 00556 std::cerr << "Rational(" << *d_n << "): created " << num_created 00557 << ", deleted " << num_deleted 00558 << ", currently alive " << num_created-num_deleted 00559 << std::endl; 00560 } 00561 } 00562 #endif 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 bool operator<=(const Rational &n1, const Rational &n2) { 00573 return(*n1.d_n <= *n2.d_n); 00574 } 00575 00576 bool operator>(const Rational &n1, const Rational &n2) { 00577 return(*n1.d_n > *n2.d_n); 00578 } 00579 00580 bool operator>=(const Rational &n1, const Rational &n2) { 00581 return(*n1.d_n >= *n2.d_n); 00582 } 00583 00584 bool operator!=(const Rational &n1, const Rational &n2) { 00585 return(*n1.d_n != *n2.d_n); 00586 } 00587 00588 Rational operator+(const Rational &n1, const Rational &n2) { 00589 return Rational(Rational::Impl(*n1.d_n + *n2.d_n)); 00590 } 00591 00592 Rational operator-(const Rational &n1, const Rational &n2) { 00593 return Rational(Rational::Impl((*n1.d_n) - (*n2.d_n))); 00594 } 00595 00596 Rational operator*(const Rational &n1, const Rational &n2) { 00597 return Rational(Rational::Impl((*n1.d_n) * (*n2.d_n))); 00598 } 00599 00600 Rational operator/(const Rational &n1, const Rational &n2) { 00601 return Rational(Rational::Impl(*n1.d_n / *n2.d_n)); 00602 } 00603 00604 Rational operator%(const Rational &n1, const Rational &n2) { 00605 return Rational(Rational::Impl(*n1.d_n % *n2.d_n)); 00606 } 00607 00608 } /* close namespace */ 00609 00610 #endif