CVC3

rational-gmp.cpp

Go to the documentation of this file.
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     Unsigned getUnsignedMP() const;
00110 
00111     //! Unary minus
00112     Impl operator-() const;
00113     //! Get the floor
00114     Impl floor() const;
00115     //! Get the ceiling
00116     Impl ceil() const;
00117 
00118     //! Testing whether the denominator is 1
00119     bool isInteger() const;
00120 
00121     //! Equality
00122     friend bool operator==(const Impl& x, const Impl& y) {
00123       return mpq_equal(x.d_n, y.d_n);
00124     }
00125 
00126     //! Dis-equality
00127     friend bool operator!=(const Impl& x, const Impl& y) {
00128       return !mpq_equal(x.d_n, y.d_n);
00129     }
00130     //! Less than
00131     friend bool operator<(const Impl& x, const Impl& y) {
00132       return mpq_cmp(x.d_n, y.d_n) < 0;
00133     }
00134 
00135     friend bool operator<=(const Impl& x, const Impl& y) {
00136       return mpq_cmp(x.d_n, y.d_n) <= 0;
00137     }
00138 
00139     friend bool operator>(const Impl& x, const Impl& y) {
00140       return mpq_cmp(x.d_n, y.d_n) > 0;
00141     }
00142 
00143     friend bool operator>=(const Impl& x, const Impl& y) {
00144       return mpq_cmp(x.d_n, y.d_n) >= 0;
00145     }
00146 
00147     //! Addition
00148     friend Impl operator+(const Impl& x, const Impl& y) {
00149       Impl res;
00150       mpq_add(res.d_n, x.d_n, y.d_n);
00151       return res;
00152     }
00153 
00154     //! Subtraction
00155     friend Impl operator-(const Impl& x, const Impl& y) {
00156       Impl res;
00157       mpq_sub(res.d_n, x.d_n, y.d_n);
00158       return res;
00159     }
00160 
00161     //! Multiplication
00162     friend Impl operator*(const Impl& x, const Impl& y) {
00163       Impl res;
00164       mpq_mul(res.d_n, x.d_n, y.d_n);
00165       return res;
00166     }
00167 
00168     //! Division
00169     friend Impl operator/(const Impl& x, const Impl& y) {
00170       Impl res;
00171       mpq_div(res.d_n, x.d_n, y.d_n);
00172       return res;
00173     }
00174 
00175     friend Impl operator%(const Impl& x, const Impl& y) {
00176       DebugAssert(x.isInteger() && y.isInteger(),
00177       "Impl % Impl: x and y must be integers");
00178       mpz_t res;
00179       mpz_init(res);
00180       // Put the remainder directly into 'res'
00181       mpz_fdiv_r(res, mpq_numref(x.d_n), mpq_numref(y.d_n));
00182       Impl r(res);
00183       mpz_clear(res);
00184       return r;
00185     }
00186 
00187     //! Compute the remainder of x/y
00188     friend Impl mod(const Impl& x, const Impl& y) {
00189       DebugAssert(x.isInteger() && y.isInteger(),
00190       "Rational::Impl::mod(): x="+x.toString()
00191       +", y="+y.toString());
00192       mpz_t res;
00193       mpz_init(res);
00194       mpz_mod(res, mpq_numref(x.d_n), mpq_numref(y.d_n));
00195       Impl r(res);
00196       mpz_clear(res);
00197       return r;
00198     }
00199 
00200     friend Impl intRoot(const Impl& x, unsigned long int y) {
00201       DebugAssert(x.isInteger(),
00202       "Rational::Impl::intRoot(): x="+x.toString());
00203       mpz_t res;
00204       mpz_init(res);
00205       int exact = mpz_root(res, mpq_numref(x.d_n), y);
00206       if (!exact) {
00207         mpz_set_ui(res, 0);
00208       }
00209       Impl r(res);
00210       mpz_clear(res);
00211       return r;
00212     }
00213 
00214     friend Impl gcd(const Impl& x, const Impl& y) {
00215       DebugAssert(x.isInteger() && y.isInteger(),
00216       "Rational::Impl::gcd(): x="+x.toString()
00217       +", y="+y.toString());
00218       TRACE("rational", "Impl::gcd(", x, ", "+y.toString()+") {");
00219       mpz_t res;
00220       mpz_init(res);
00221       mpz_gcd(res, mpq_numref(x.d_n), mpq_numref(y.d_n));
00222       Impl r(res);
00223       mpz_clear(res);
00224       TRACE("rational", "Impl::gcd => ", r, "}");
00225       return r;
00226     }
00227 
00228     friend Impl lcm(const Impl& x, const Impl& y) {
00229       DebugAssert(x.isInteger() && y.isInteger(),
00230       "Rational::Impl::lcm(): x="+x.toString()
00231       +", y="+y.toString());
00232       mpz_t res;
00233       mpz_init(res);
00234       mpz_lcm(res, mpq_numref(x.d_n), mpq_numref(y.d_n));
00235       Impl r(res);
00236       mpz_clear(res);
00237       return r;
00238     }
00239 
00240     //! Print to string
00241     string toString(int base = 10) const {
00242       char* str = (char*)malloc(mpz_sizeinbase(mpq_numref(d_n), base)
00243         +mpz_sizeinbase(mpq_denref(d_n), base)+3);
00244       mpq_get_str(str, base, d_n);
00245       string res(str);
00246       free(str);
00247       return res;
00248     }
00249   };
00250 
00251   // Constructor from a pair of integers
00252   Rational::Impl::Impl(long int n, long int d) {
00253     mpq_init(d_n);
00254     DebugAssert(d > 0, "Rational::Impl(long n, long d): d = "+int2string(d));
00255     mpq_set_si(d_n, n, (unsigned long int)d);
00256     canonicalize();
00257   }
00258 
00259   // Constructor from a pair of unsigned integers
00260   Rational::Impl::Impl(unsigned int n, unsigned int d,
00261            unsigned int /* dummy arg, to disambiguate */) {
00262     mpq_init(d_n);
00263     mpq_set_ui(d_n, n, (unsigned long int)d);
00264     canonicalize();
00265   }
00266 
00267   // Constructor from a string
00268   Rational::Impl::Impl(const string &n, int base) {
00269     mpq_init(d_n);
00270     mpq_set_str(d_n, n.c_str(), base);
00271     canonicalize();
00272   }
00273 
00274   // Constructor from a pair of strings
00275   Rational::Impl::Impl(const string &n, const string& d, int base) {
00276     mpq_init(d_n);
00277     mpq_set_str(d_n, (n+"/"+d).c_str(), base);
00278     canonicalize();
00279   }
00280 
00281   Rational::Impl Rational::Impl::operator-() const {
00282     Impl res;
00283     mpq_neg(res.d_n, d_n);
00284     return res;
00285   }
00286 
00287   Rational::Impl Rational::Impl::floor() const {
00288     mpz_t res;
00289     mpz_init(res);
00290     mpz_fdiv_q(res, mpq_numref(d_n), mpq_denref(d_n));
00291     Impl r(res);
00292     mpz_clear(res);
00293     return r;
00294   }
00295 
00296   Rational::Impl Rational::Impl::ceil() const {
00297     mpz_t res;
00298     mpz_init(res);
00299     mpz_cdiv_q(res, mpq_numref(d_n), mpq_denref(d_n));
00300     Impl r(res);
00301     mpz_clear(res);
00302     return r;
00303   }
00304 
00305   bool Rational::Impl::isInteger() const {
00306     bool res(mpz_cmp_ui(mpq_denref(d_n), 1) == 0);
00307     TRACE("rational", "Impl::isInteger(", *this,
00308     ") => "+string(res? "true" : "false"));
00309     return res;
00310   }
00311 
00312   //! Default constructor
00313   Rational::Rational() : d_n(new Impl) {
00314 #ifdef _DEBUG_RATIONAL_
00315     int &num_created = getCreated();
00316     num_created++;
00317     printStats();
00318 #endif
00319   }
00320   //! Copy constructor
00321   Rational::Rational(const Rational &n) : d_n(new Impl(*n.d_n)) {
00322 #ifdef _DEBUG_RATIONAL_
00323     int &num_created = getCreated();
00324     num_created++;
00325     printStats();
00326 #endif
00327   }
00328 
00329   Rational::Rational(const Unsigned& n): d_n(new Impl(n.toString(), 10)) {
00330 #ifdef _DEBUG_RATIONAL_
00331     int &num_created = getCreated();
00332     num_created++;
00333     printStats();
00334 #endif
00335   }
00336   //! Private constructor
00337   Rational::Rational(const Impl& t): d_n(new Impl(t)) {
00338 #ifdef _DEBUG_RATIONAL_
00339     int &num_created = getCreated();
00340     num_created++;
00341     printStats();
00342 #endif
00343   }
00344   Rational::Rational(int n, int d): d_n(new Impl(n, d)) {
00345 #ifdef _DEBUG_RATIONAL_
00346     int &num_created = getCreated();
00347     num_created++;
00348     printStats();
00349 #endif
00350   }
00351   // Constructors from strings
00352   Rational::Rational(const char* n, int base)
00353     : d_n(new Impl(string(n), base)) {
00354 #ifdef _DEBUG_RATIONAL_
00355     int &num_created = getCreated();
00356     num_created++;
00357     printStats();
00358 #endif
00359   }
00360   Rational::Rational(const string& n, int base)
00361     : d_n(new Impl(n, base)) {
00362 #ifdef _DEBUG_RATIONAL_
00363     int &num_created = getCreated();
00364     num_created++;
00365     printStats();
00366 #endif
00367   }
00368   Rational::Rational(const char* n, const char* d, int base)
00369     : d_n(new Impl(string(n), string(d), base)) {
00370 #ifdef _DEBUG_RATIONAL_
00371     int &num_created = getCreated();
00372     num_created++;
00373     printStats();
00374 #endif
00375   }
00376   Rational::Rational(const string& n, const string& d, int base)
00377     : d_n(new Impl(n, d, base)) {
00378 #ifdef _DEBUG_RATIONAL_
00379     int &num_created = getCreated();
00380     num_created++;
00381     printStats();
00382 #endif
00383   }
00384   // Destructor
00385   Rational::~Rational() {
00386 #ifdef _DEBUG_RATIONAL_
00387     int &num_deleted = getDeleted();
00388     num_deleted++;
00389     printStats();
00390 #endif
00391     delete d_n;
00392   }
00393 
00394   // Get components
00395   Rational Rational::getNumerator() const
00396     { return Rational(d_n->getNum()); }
00397   Rational Rational::getDenominator() const
00398     { return Rational(d_n->getDen()); }
00399 
00400   bool Rational::isInteger() const { return d_n->isInteger(); }
00401 
00402   // Assignment
00403   Rational& Rational::operator=(const Rational& n) {
00404     if(this == &n) return *this;
00405     delete d_n;
00406     d_n = new Impl(*n.d_n);
00407     return *this;
00408   }
00409 
00410   ostream &operator<<(ostream &os, const Rational &n) {
00411     return(os << n.toString());
00412   }
00413 
00414   //! Printing to ostream
00415   std::ostream& operator<<(std::ostream& os, const Rational::Impl& n) {
00416     return os << n.toString();
00417   }
00418 
00419 
00420   // Check that argument is an int and print an error message otherwise
00421 
00422   static void checkInt(const Rational& n, const string& funName) {
00423     TRACE("rational", "checkInt(", n, ")");
00424     DebugAssert(n.isInteger(),
00425     "CVC3::Rational::" + funName
00426     + ": argument is not an integer: " + n.toString());
00427   }
00428 
00429     /* Computes gcd and lcm on *integer* values. Result is always a
00430        positive integer.  In this implementation, it is guaranteed by
00431        GMP. */
00432 
00433   Rational gcd(const Rational &x, const Rational &y) {
00434     checkInt(x, "gcd(*x*,y)");
00435     checkInt(y, "gcd(x,*y*)");
00436     return Rational(gcd(*x.d_n, *y.d_n));
00437   }
00438 
00439   Rational gcd(const vector<Rational> &v) {
00440     Rational::Impl g(1,1), zero;
00441     if(v.size() > 0) {
00442       checkInt(v[0], "gcd(vector<Rational>[0])");
00443       g = *v[0].d_n;
00444     }
00445     for(size_t i=1; i<v.size(); i++) {
00446       checkInt(v[i], "gcd(vector<Rational>)");
00447       if(g == zero)
00448   g = *(v[i].d_n);
00449       else if(*(v[i].d_n) != zero)
00450   g = gcd(g, *(v[i].d_n));
00451     }
00452     return Rational(g);
00453   }
00454 
00455   Rational lcm(const Rational &x, const Rational &y) {
00456     checkInt(x, "lcm(*x*,y)");
00457     checkInt(y, "lcm(x,*y*)");
00458     return Rational(lcm(*x.d_n, *y.d_n));
00459   }
00460 
00461   Rational lcm(const vector<Rational> &v) {
00462     Rational::Impl g(1,1), zero;
00463     for(size_t i=0; i<v.size(); i++) {
00464       checkInt(v[i], "lcm(vector<Rational>)");
00465       if(*v[i].d_n != zero)
00466   g = lcm(g, *v[i].d_n);
00467     }
00468     return Rational(g);
00469   }
00470 
00471   Rational abs(const Rational &x) {
00472     if(x<0) return -x;
00473     return x;
00474   }
00475 
00476   Rational floor(const Rational &x) {
00477     return Rational(x.d_n->floor());
00478   }
00479 
00480   Rational ceil(const Rational &x) {
00481     return Rational(x.d_n->ceil());
00482   }
00483 
00484   Rational mod(const Rational &x, const Rational &y) {
00485     checkInt(x, "mod(*x*,y)");
00486     checkInt(y, "mod(x,*y*)");
00487     return(Rational(mod(*x.d_n, *y.d_n)));
00488   }
00489 
00490   Rational intRoot(const Rational& base, unsigned long int n) {
00491     checkInt(base, "intRoot(*x*,y)");
00492     return Rational(intRoot(*base.d_n, n));
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(-(*d_n));
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     return d_n->getInt();
00536   }
00537 
00538   unsigned int Rational::getUnsigned() const {
00539     checkInt(*this, "getUnsigned()");
00540     return d_n->getUnsigned();
00541   }
00542 
00543   Unsigned Rational::getUnsignedMP() const {
00544     checkInt(*this, "getUnsignedMP()");
00545     return d_n->getUnsignedMP();
00546   }
00547 
00548 #ifdef _DEBUG_RATIONAL_
00549   void Rational::printStats() {
00550       int &num_created = getCreated();
00551       int &num_deleted = getDeleted();
00552       if(num_created % 1000 == 0 || num_deleted % 1000 == 0) {
00553   std::cerr << "Rational(" << *d_n << "): created " << num_created
00554       << ", deleted " << num_deleted
00555       << ", currently alive " << num_created-num_deleted
00556       << std::endl;
00557       }
00558     }
00559 #endif
00560 
00561     bool operator==(const Rational &n1, const Rational &n2) {
00562       return(*n1.d_n == *n2.d_n);
00563     }
00564 
00565     bool operator<(const Rational &n1, const Rational &n2) {
00566       return(*n1.d_n < *n2.d_n);
00567     }
00568 
00569     bool operator<=(const Rational &n1, const Rational &n2) {
00570       return(*n1.d_n <= *n2.d_n);
00571     }
00572 
00573     bool operator>(const Rational &n1, const Rational &n2) {
00574       return(*n1.d_n > *n2.d_n);
00575     }
00576 
00577     bool operator>=(const Rational &n1, const Rational &n2) {
00578       return(*n1.d_n >= *n2.d_n);
00579     }
00580 
00581     bool operator!=(const Rational &n1, const Rational &n2) {
00582       return(*n1.d_n != *n2.d_n);
00583     }
00584 
00585     Rational operator+(const Rational &n1, const Rational &n2) {
00586       return Rational((*n1.d_n) + (*n2.d_n));
00587     }
00588 
00589     Rational operator-(const Rational &n1, const Rational &n2) {
00590       return Rational((*n1.d_n) - (*n2.d_n));
00591     }
00592 
00593     Rational operator*(const Rational &n1, const Rational &n2) {
00594       return Rational((*n1.d_n) * (*n2.d_n));
00595     }
00596 
00597     Rational operator/(const Rational &n1, const Rational &n2) {
00598       return Rational((*n1.d_n) / (*n2.d_n));
00599     }
00600 
00601     Rational operator%(const Rational &n1, const Rational &n2) {
00602       return Rational((*n1.d_n) % (*n2.d_n));
00603     }
00604 
00605   //! Implementation of the forward-declared internal class
00606   class Unsigned::Impl {
00607     mpz_t d_n;
00608   public:
00609     //! Default Constructor
00610     Impl() { mpz_init(d_n); }
00611     //! Copy constructor (assumes x is canonicalized)
00612     Impl(const Impl &x) { mpz_init(d_n); mpz_set(d_n, x.d_n); }
00613     //! Constructor from mpz_t
00614     Impl(const mpz_t n) {
00615       mpz_init(d_n);
00616       mpz_set(d_n, n);
00617     }
00618     //! Constructor from an unsigned integer
00619     Impl(unsigned int n);
00620     //! Constructor from a string
00621     Impl(const string &n, int base);
00622     // Destructor
00623     virtual ~Impl() { mpz_clear(d_n); }
00624 
00625     //! Assignment
00626     Impl& operator=(const Impl& x) {
00627       if(this == &x) return *this;
00628       mpz_set(d_n, x.d_n);
00629       return *this;
00630     }
00631 
00632     int getInt() const {
00633       // Check for overflow
00634       static Impl max((int)INT_MAX);
00635       FatalAssert(*this <= max,
00636                   "Unsigned::getInt(): Arithmetic overflow for "
00637                   +toString());
00638       return mpz_get_si(d_n);
00639     }
00640 
00641     unsigned long getUnsigned() const {
00642       // Check for overflow
00643       static Impl max(UINT_MAX);
00644       FatalAssert(*this <= max,
00645                   "Unsigned::getUnsigned(): Arithmetic overflow for "
00646                   +toString());
00647       return mpz_get_ui(d_n);
00648     }
00649 
00650     //! Unary minus
00651     Impl operator-() const;
00652 
00653     //! Equality
00654     friend bool operator==(const Impl& x, const Impl& y) {
00655       return mpz_cmp(x.d_n, y.d_n) == 0;
00656     }
00657 
00658     //! Dis-equality
00659     friend bool operator!=(const Impl& x, const Impl& y) {
00660       return mpz_cmp(x.d_n, y.d_n) != 0;
00661     }
00662     //! Less than
00663     friend bool operator<(const Impl& x, const Impl& y) {
00664       return mpz_cmp(x.d_n, y.d_n) < 0;
00665     }
00666 
00667     friend bool operator<=(const Impl& x, const Impl& y) {
00668       return mpz_cmp(x.d_n, y.d_n) <= 0;
00669     }
00670 
00671     friend bool operator>(const Impl& x, const Impl& y) {
00672       return mpz_cmp(x.d_n, y.d_n) > 0;
00673     }
00674 
00675     friend bool operator>=(const Impl& x, const Impl& y) {
00676       return mpz_cmp(x.d_n, y.d_n) >= 0;
00677     }
00678 
00679     //! Addition
00680     friend Impl operator+(const Impl& x, const Impl& y) {
00681       Impl res;
00682       mpz_add(res.d_n, x.d_n, y.d_n);
00683       return res;
00684     }
00685 
00686     //! Subtraction
00687     friend Impl operator-(const Impl& x, const Impl& y) {
00688       Impl res;
00689       mpz_sub(res.d_n, x.d_n, y.d_n);
00690       return res;
00691     }
00692 
00693     //! Multiplication
00694     friend Impl operator*(const Impl& x, const Impl& y) {
00695       Impl res;
00696       mpz_mul(res.d_n, x.d_n, y.d_n);
00697       return res;
00698     }
00699 
00700     //! Division
00701     friend Impl operator/(const Impl& x, const Impl& y) {
00702       Impl res;
00703       mpz_div(res.d_n, x.d_n, y.d_n);
00704       return res;
00705     }
00706 
00707     friend Impl operator%(const Impl& x, const Impl& y) {
00708       mpz_t res;
00709       mpz_init(res);
00710       // Put the remainder directly into 'res'
00711       mpz_fdiv_r(res, x.d_n, y.d_n);
00712       Impl r(res);
00713       mpz_clear(res);
00714       return r;
00715     }
00716 
00717     friend Impl operator<<(const Impl& x, unsigned y) {
00718       mpz_t res;
00719       mpz_init(res);
00720       mpz_mul_2exp (res, x.d_n, y);
00721       Impl r(res);
00722       mpz_clear(res);
00723       return r;
00724     }
00725 
00726     friend Impl operator&(const Impl& x, const Impl& y) {
00727       mpz_t res;
00728       mpz_init(res);
00729       mpz_and (res, x.d_n, y.d_n);
00730       Impl r(res);
00731       mpz_clear(res);
00732       return r;
00733     }
00734 
00735     //! Compute the remainder of x/y
00736     friend Impl mod(const Impl& x, const Impl& y) {
00737       mpz_t res;
00738       mpz_init(res);
00739       mpz_mod(res, x.d_n, y.d_n);
00740       Impl r(res);
00741       mpz_clear(res);
00742       return r;
00743     }
00744 
00745     friend Impl intRoot(const Impl& x, unsigned long int y) {
00746       mpz_t res;
00747       mpz_init(res);
00748       int exact = mpz_root(res, x.d_n, y);
00749       if (!exact) {
00750         mpz_set_ui(res, 0);
00751       }
00752       Impl r(res);
00753       mpz_clear(res);
00754       return r;
00755     }
00756 
00757     friend Impl gcd(const Impl& x, const Impl& y) {
00758       mpz_t res;
00759       mpz_init(res);
00760       mpz_gcd(res, x.d_n, y.d_n);
00761       Impl r(res);
00762       mpz_clear(res);
00763       return r;
00764     }
00765 
00766     friend Impl lcm(const Impl& x, const Impl& y) {
00767       mpz_t res;
00768       mpz_init(res);
00769       mpz_lcm(res, x.d_n, y.d_n);
00770       Impl r(res);
00771       mpz_clear(res);
00772       return r;
00773     }
00774 
00775     //! Print to string
00776     string toString(int base = 10) const {
00777       char* str = (char*)malloc(mpz_sizeinbase(d_n, base)+2);
00778       mpz_get_str(str, base, d_n);
00779       string res(str);
00780       free(str);
00781       return res;
00782     }
00783 
00784     //! Printing to ostream
00785     friend ostream& operator<<(ostream& os, const Unsigned::Impl& n) {
00786       return os << n.toString();
00787     }
00788 
00789   };
00790 
00791   // Constructor from a pair of unsigned integers
00792   Unsigned::Impl::Impl(unsigned int n) {
00793     mpz_init(d_n);
00794     mpz_set_ui(d_n, n);
00795   }
00796 
00797   // Constructor from a string
00798   Unsigned::Impl::Impl(const string &n, int base) {
00799     mpz_init(d_n);
00800     mpz_set_str(d_n, n.c_str(), base);
00801   }
00802 
00803   Unsigned::Impl Unsigned::Impl::operator-() const {
00804     Impl res;
00805     mpz_neg(res.d_n, d_n);
00806     return res;
00807   }
00808 
00809   //! Default constructor
00810   Unsigned::Unsigned() : d_n(new Impl) { }
00811   //! Copy constructor
00812   Unsigned::Unsigned(const Unsigned &n) : d_n(new Impl(*n.d_n)) { }
00813 
00814   Unsigned::Unsigned(int n) : d_n(new Impl((unsigned)n)) { }
00815 
00816   Unsigned::Unsigned(unsigned n) : d_n(new Impl(n)) { }
00817   //! Private constructor
00818   Unsigned::Unsigned(const Impl& t): d_n(new Impl(t)) { }
00819   // Constructors from strings
00820   Unsigned::Unsigned(const char* n, int base)
00821     : d_n(new Impl(string(n), base)) { }
00822   Unsigned::Unsigned(const string& n, int base)
00823     : d_n(new Impl(n, base)) { }
00824   // Destructor
00825   Unsigned::~Unsigned() {
00826     delete d_n;
00827   }
00828 
00829   // Assignment
00830   Unsigned& Unsigned::operator=(const Unsigned& n) {
00831     if(this == &n) return *this;
00832     delete d_n;
00833     d_n = new Impl(*n.d_n);
00834     return *this;
00835   }
00836 
00837   ostream &operator<<(ostream &os, const Unsigned &n) {
00838     return(os << n.toString());
00839   }
00840 
00841 
00842   // Check that argument is an int and print an error message otherwise
00843 
00844     /* Computes gcd and lcm on *integer* values. Result is always a
00845        positive integer.  In this implementation, it is guaranteed by
00846        GMP. */
00847 
00848   Unsigned gcd(const Unsigned &x, const Unsigned &y) {
00849     return Unsigned(gcd(*x.d_n, *y.d_n));
00850   }
00851 
00852   Unsigned gcd(const vector<Unsigned> &v) {
00853     Unsigned::Impl g(1), zero;
00854     if(v.size() > 0) {
00855       g = *v[0].d_n;
00856     }
00857     for(size_t i=1; i<v.size(); i++) {
00858       if(g == zero)
00859   g = *(v[i].d_n);
00860       else if(*(v[i].d_n) != zero)
00861   g = gcd(g, *(v[i].d_n));
00862     }
00863     return Unsigned(g);
00864   }
00865 
00866   Unsigned lcm(const Unsigned &x, const Unsigned &y) {
00867     return Unsigned(lcm(*x.d_n, *y.d_n));
00868   }
00869 
00870   Unsigned lcm(const vector<Unsigned> &v) {
00871     Unsigned::Impl g(1), zero;
00872     for(size_t i=0; i<v.size(); i++) {
00873       if(*v[i].d_n != zero)
00874   g = lcm(g, *v[i].d_n);
00875     }
00876     return Unsigned(g);
00877   }
00878 
00879   Unsigned mod(const Unsigned &x, const Unsigned &y) {
00880     return(Unsigned(mod(*x.d_n, *y.d_n)));
00881   }
00882 
00883   Unsigned intRoot(const Unsigned& base, unsigned long int n) {
00884     return Unsigned(intRoot(*base.d_n, n));
00885   }
00886 
00887   string Unsigned::toString(int base) const {
00888     return(d_n->toString(base));
00889   }
00890 
00891   size_t Unsigned::hash() const {
00892     std::hash<const char *> h;
00893     return h(toString().c_str());
00894   }
00895 
00896   void Unsigned::print() const {
00897     cout << (*this) << endl;
00898   }
00899 
00900   Unsigned &Unsigned::operator+=(const Unsigned &n2) {
00901     *this = (*this) + n2;
00902     return *this;
00903   }
00904 
00905   Unsigned &Unsigned::operator-=(const Unsigned &n2) {
00906     *this = (*this) - n2;
00907     return *this;
00908   }
00909 
00910   Unsigned &Unsigned::operator*=(const Unsigned &n2) {
00911     *this = (*this) * n2;
00912     return *this;
00913   }
00914 
00915   Unsigned &Unsigned::operator/=(const Unsigned &n2) {
00916     *this = (*this) / n2;
00917     return *this;
00918   }
00919 
00920   unsigned long Unsigned::getUnsigned() const {
00921     return d_n->getUnsigned();
00922   }
00923 
00924   bool operator==(const Unsigned &n1, const Unsigned &n2) {
00925     return(*n1.d_n == *n2.d_n);
00926   }
00927 
00928   bool operator<(const Unsigned &n1, const Unsigned &n2) {
00929     return(*n1.d_n < *n2.d_n);
00930   }
00931 
00932   bool operator<=(const Unsigned &n1, const Unsigned &n2) {
00933     return(*n1.d_n <= *n2.d_n);
00934   }
00935 
00936   bool operator>(const Unsigned &n1, const Unsigned &n2) {
00937     return(*n1.d_n > *n2.d_n);
00938   }
00939 
00940   bool operator>=(const Unsigned &n1, const Unsigned &n2) {
00941     return(*n1.d_n >= *n2.d_n);
00942   }
00943 
00944   bool operator!=(const Unsigned &n1, const Unsigned &n2) {
00945     return(*n1.d_n != *n2.d_n);
00946   }
00947 
00948   Unsigned operator+(const Unsigned &n1, const Unsigned &n2) {
00949     return Unsigned((*n1.d_n) + (*n2.d_n));
00950   }
00951 
00952   Unsigned operator-(const Unsigned &n1, const Unsigned &n2) {
00953     return Unsigned((*n1.d_n) - (*n2.d_n));
00954   }
00955 
00956   Unsigned operator*(const Unsigned &n1, const Unsigned &n2) {
00957     return Unsigned((*n1.d_n) * (*n2.d_n));
00958   }
00959 
00960   Unsigned operator/(const Unsigned &n1, const Unsigned &n2) {
00961     return Unsigned((*n1.d_n) / (*n2.d_n));
00962   }
00963 
00964   Unsigned operator%(const Unsigned &n1, const Unsigned &n2) {
00965     return Unsigned((*n1.d_n) % (*n2.d_n));
00966   }
00967 
00968   Unsigned operator<<(const Unsigned &n1, unsigned n2) {
00969     return Unsigned((*n1.d_n) << n2);
00970   }
00971 
00972   Unsigned operator&(const Unsigned &n1, const Unsigned &n2) {
00973     return Unsigned((*n1.d_n) & (*n2.d_n));
00974   }
00975 
00976   Unsigned Rational::Impl::getUnsignedMP() const {
00977     return Unsigned(Unsigned::Impl(mpq_numref(d_n)));
00978   }
00979 
00980 
00981 } /* close namespace */
00982 
00983 #endif